| 1 | mmeineke | 377 | #ifndef _INTEGRATOR_H_ | 
| 2 |  |  | #define _INTEGRATOR_H_ | 
| 3 |  |  |  | 
| 4 |  |  | #include "Atom.hpp" | 
| 5 |  |  | #include "SRI.hpp" | 
| 6 |  |  | #include "AbstractClasses.hpp" | 
| 7 |  |  | #include "SimInfo.hpp" | 
| 8 |  |  | #include "ForceFields.hpp" | 
| 9 | mmeineke | 540 | #include "Thermo.hpp" | 
| 10 |  |  | #include "ReadWrite.hpp" | 
| 11 | mmeineke | 377 |  | 
| 12 | mmeineke | 561 | const double kB = 8.31451e-7;// boltzmann constant amu*Ang^2*fs^-2/K | 
| 13 |  |  | const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2 | 
| 14 | mmeineke | 586 | const double p_convert = 1.63882576e8; //converts amu*fs^-2*Ang^-1 -> atm | 
| 15 | mmeineke | 561 | const int maxIteration = 300; | 
| 16 |  |  | const double tol = 1.0e-6; | 
| 17 |  |  |  | 
| 18 | mmeineke | 540 | class Integrator : public BaseIntegrator { | 
| 19 | mmeineke | 377 |  | 
| 20 |  |  | public: | 
| 21 | mmeineke | 561 | Integrator( SimInfo *theInfo, ForceFields* the_ff ); | 
| 22 | mmeineke | 548 | virtual ~Integrator(); | 
| 23 | mmeineke | 377 | void integrate( void ); | 
| 24 |  |  |  | 
| 25 |  |  |  | 
| 26 | mmeineke | 540 | protected: | 
| 27 | mmeineke | 377 |  | 
| 28 | mmeineke | 540 | virtual void integrateStep( int calcPot, int calcStress ); | 
| 29 | mmeineke | 548 | virtual void preMove( void ); | 
| 30 | mmeineke | 540 | virtual void moveA( void ); | 
| 31 |  |  | virtual void moveB( void ); | 
| 32 |  |  | virtual void constrainA( void ); | 
| 33 |  |  | virtual void constrainB( void ); | 
| 34 | mmeineke | 559 | virtual int  readyCheck( void ) { return 1; } | 
| 35 | mmeineke | 377 |  | 
| 36 | mmeineke | 540 | void checkConstraints( void ); | 
| 37 | mmeineke | 377 | void rotate( int axes1, int axes2, double angle, double j[3], | 
| 38 | mmeineke | 597 | double A[9] ); | 
| 39 | mmeineke | 377 |  | 
| 40 | mmeineke | 540 |  | 
| 41 | mmeineke | 377 | ForceFields* myFF; | 
| 42 |  |  |  | 
| 43 | mmeineke | 540 | SimInfo *info; // all the info we'll ever need | 
| 44 |  |  | int nAtoms;  /* the number of atoms */ | 
| 45 | mmeineke | 548 | int oldAtoms; | 
| 46 | mmeineke | 540 | Atom **atoms; /* array of atom pointers */ | 
| 47 | mmeineke | 423 | Molecule* molecules; | 
| 48 |  |  | int nMols; | 
| 49 |  |  |  | 
| 50 | mmeineke | 548 | int isConstrained; // boolean to know whether the systems contains | 
| 51 |  |  | // constraints. | 
| 52 |  |  | int nConstrained;  // counter for number of constraints | 
| 53 |  |  | int *constrainedA; // the i of a constraint pair | 
| 54 |  |  | int *constrainedB; // the j of a constraint pair | 
| 55 |  |  | double *constrainedDsqr; // the square of the constraint distance | 
| 56 |  |  |  | 
| 57 |  |  | int* moving; // tells whether we are moving atom i | 
| 58 |  |  | int* moved;  // tells whether we have moved atom i | 
| 59 | mmeineke | 561 | double* oldPos; // pre constrained positions | 
| 60 | mmeineke | 548 |  | 
| 61 | mmeineke | 540 | short isFirst; /*boolean for the first time integrate is called */ | 
| 62 |  |  |  | 
| 63 |  |  | double dt; | 
| 64 | mmeineke | 541 | double dt2; | 
| 65 | gezelter | 560 |  | 
| 66 | mmeineke | 541 | double* pos; | 
| 67 |  |  | double* vel; | 
| 68 |  |  | double* frc; | 
| 69 |  |  | double* trq; | 
| 70 |  |  | double* Amat; | 
| 71 | mmeineke | 377 |  | 
| 72 | mmeineke | 540 | Thermo *tStats; | 
| 73 |  |  | StatWriter*  statOut; | 
| 74 |  |  | DumpWriter*  dumpOut; | 
| 75 | mmeineke | 377 |  | 
| 76 |  |  | }; | 
| 77 |  |  |  | 
| 78 | mmeineke | 555 | class NVE : public Integrator{ | 
| 79 | mmeineke | 540 |  | 
| 80 | mmeineke | 561 | public: | 
| 81 |  |  | NVE ( SimInfo *theInfo, ForceFields* the_ff ): | 
| 82 | mmeineke | 555 | Integrator( theInfo, the_ff ){} | 
| 83 |  |  | virtual ~NVE(){} | 
| 84 |  |  |  | 
| 85 |  |  |  | 
| 86 |  |  |  | 
| 87 |  |  | }; | 
| 88 |  |  |  | 
| 89 | mmeineke | 540 | class NVT : public Integrator{ | 
| 90 |  |  |  | 
| 91 | gezelter | 560 | public: | 
| 92 |  |  |  | 
| 93 | mmeineke | 561 | NVT ( SimInfo *theInfo, ForceFields* the_ff); | 
| 94 |  |  | virtual ~NVT() {} | 
| 95 | mmeineke | 540 |  | 
| 96 | gezelter | 560 | void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} | 
| 97 |  |  | void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} | 
| 98 |  |  |  | 
| 99 | mmeineke | 540 | protected: | 
| 100 | gezelter | 560 |  | 
| 101 | mmeineke | 561 | virtual void moveA( void ); | 
| 102 |  |  | virtual void moveB( void ); | 
| 103 | mmeineke | 540 |  | 
| 104 | mmeineke | 561 | virtual int readyCheck(); | 
| 105 | gezelter | 560 |  | 
| 106 | gezelter | 565 | // chi is a propagated degree of freedom. | 
| 107 | gezelter | 560 |  | 
| 108 | gezelter | 565 | double chi; | 
| 109 | gezelter | 560 |  | 
| 110 | gezelter | 565 | // targetTemp must be set.  tauThermostat must also be set; | 
| 111 | gezelter | 560 |  | 
| 112 |  |  | double targetTemp; | 
| 113 |  |  | double tauThermostat; | 
| 114 | mmeineke | 561 |  | 
| 115 | gezelter | 565 | short int have_tau_thermostat, have_target_temp; | 
| 116 | gezelter | 560 |  | 
| 117 | mmeineke | 540 | }; | 
| 118 |  |  |  | 
| 119 |  |  |  | 
| 120 | gezelter | 574 | class NPTi : public Integrator{ | 
| 121 | mmeineke | 540 |  | 
| 122 | gezelter | 560 | public: | 
| 123 |  |  |  | 
| 124 | gezelter | 574 | NPTi ( SimInfo *theInfo, ForceFields* the_ff); | 
| 125 |  |  | virtual ~NPTi() {}; | 
| 126 | gezelter | 560 |  | 
| 127 | mmeineke | 594 | virtual void integrateStep( int calcPot, int calcStress ){ | 
| 128 |  |  | calcStress = 1; | 
| 129 |  |  | Integrator::integrateStep( calcPot, calcStress ); | 
| 130 |  |  | } | 
| 131 |  |  |  | 
| 132 | gezelter | 560 | void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} | 
| 133 |  |  | void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} | 
| 134 |  |  | void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} | 
| 135 |  |  | void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} | 
| 136 |  |  |  | 
| 137 |  |  | protected: | 
| 138 |  |  |  | 
| 139 | mmeineke | 561 | virtual void  moveA( void ); | 
| 140 |  |  | virtual void moveB( void ); | 
| 141 | gezelter | 560 |  | 
| 142 | mmeineke | 561 | virtual int readyCheck(); | 
| 143 | gezelter | 560 |  | 
| 144 | gezelter | 565 | // chi and eta are the propagated degrees of freedom | 
| 145 | gezelter | 560 |  | 
| 146 | gezelter | 565 | double chi; | 
| 147 |  |  | double eta; | 
| 148 | gezelter | 574 | double NkBT; | 
| 149 | gezelter | 560 |  | 
| 150 |  |  | // targetTemp, targetPressure, and tauBarostat must be set. | 
| 151 |  |  | // One of qmass or tauThermostat must be set; | 
| 152 |  |  |  | 
| 153 |  |  | double targetTemp; | 
| 154 |  |  | double targetPressure; | 
| 155 |  |  | double tauThermostat; | 
| 156 |  |  | double tauBarostat; | 
| 157 |  |  |  | 
| 158 |  |  | short int have_tau_thermostat, have_tau_barostat, have_target_temp; | 
| 159 | gezelter | 565 | short int have_target_pressure; | 
| 160 | gezelter | 560 |  | 
| 161 |  |  | }; | 
| 162 |  |  |  | 
| 163 | gezelter | 596 | class NPTim : public Integrator{ | 
| 164 |  |  |  | 
| 165 |  |  | public: | 
| 166 |  |  |  | 
| 167 |  |  | NPTim ( SimInfo *theInfo, ForceFields* the_ff); | 
| 168 |  |  | virtual ~NPTim() {}; | 
| 169 |  |  |  | 
| 170 |  |  | virtual void integrateStep( int calcPot, int calcStress ){ | 
| 171 |  |  | calcStress = 1; | 
| 172 |  |  | Integrator::integrateStep( calcPot, calcStress ); | 
| 173 |  |  | } | 
| 174 |  |  |  | 
| 175 |  |  | void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} | 
| 176 |  |  | void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} | 
| 177 |  |  | void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} | 
| 178 |  |  | void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} | 
| 179 |  |  |  | 
| 180 |  |  | protected: | 
| 181 |  |  |  | 
| 182 |  |  | virtual void  moveA( void ); | 
| 183 |  |  | virtual void moveB( void ); | 
| 184 |  |  |  | 
| 185 |  |  | virtual int readyCheck(); | 
| 186 |  |  |  | 
| 187 |  |  | // chi and eta are the propagated degrees of freedom | 
| 188 |  |  |  | 
| 189 |  |  | double chi; | 
| 190 |  |  | double eta; | 
| 191 |  |  | double NkBT; | 
| 192 |  |  |  | 
| 193 |  |  | // targetTemp, targetPressure, and tauBarostat must be set. | 
| 194 |  |  | // One of qmass or tauThermostat must be set; | 
| 195 |  |  |  | 
| 196 |  |  | double targetTemp; | 
| 197 |  |  | double targetPressure; | 
| 198 |  |  | double tauThermostat; | 
| 199 |  |  | double tauBarostat; | 
| 200 |  |  |  | 
| 201 |  |  | short int have_tau_thermostat, have_tau_barostat, have_target_temp; | 
| 202 |  |  | short int have_target_pressure; | 
| 203 |  |  |  | 
| 204 |  |  | }; | 
| 205 |  |  |  | 
| 206 | gezelter | 576 | class NPTf : public Integrator{ | 
| 207 |  |  |  | 
| 208 |  |  | public: | 
| 209 |  |  |  | 
| 210 |  |  | NPTf ( SimInfo *theInfo, ForceFields* the_ff); | 
| 211 |  |  | virtual ~NPTf() {}; | 
| 212 |  |  |  | 
| 213 | mmeineke | 594 | virtual void integrateStep( int calcPot, int calcStress ){ | 
| 214 |  |  | calcStress = 1; | 
| 215 |  |  | Integrator::integrateStep( calcPot, calcStress ); | 
| 216 |  |  | } | 
| 217 |  |  |  | 
| 218 | gezelter | 576 | void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} | 
| 219 |  |  | void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} | 
| 220 |  |  | void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} | 
| 221 |  |  | void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} | 
| 222 |  |  |  | 
| 223 |  |  | protected: | 
| 224 |  |  |  | 
| 225 |  |  | virtual void  moveA( void ); | 
| 226 |  |  | virtual void moveB( void ); | 
| 227 |  |  |  | 
| 228 |  |  | virtual int readyCheck(); | 
| 229 |  |  |  | 
| 230 |  |  | // chi and eta are the propagated degrees of freedom | 
| 231 |  |  |  | 
| 232 |  |  | double chi; | 
| 233 | gezelter | 588 | double eta[3][3]; | 
| 234 | gezelter | 576 | double NkBT; | 
| 235 |  |  |  | 
| 236 |  |  | // targetTemp, targetPressure, and tauBarostat must be set. | 
| 237 |  |  | // One of qmass or tauThermostat must be set; | 
| 238 |  |  |  | 
| 239 |  |  | double targetTemp; | 
| 240 |  |  | double targetPressure; | 
| 241 |  |  | double tauThermostat; | 
| 242 |  |  | double tauBarostat; | 
| 243 |  |  |  | 
| 244 |  |  | short int have_tau_thermostat, have_tau_barostat, have_target_temp; | 
| 245 |  |  | short int have_target_pressure; | 
| 246 |  |  |  | 
| 247 |  |  | }; | 
| 248 |  |  |  | 
| 249 | gezelter | 596 | class NPTfm : public Integrator{ | 
| 250 |  |  |  | 
| 251 |  |  | public: | 
| 252 |  |  |  | 
| 253 |  |  | NPTfm ( SimInfo *theInfo, ForceFields* the_ff); | 
| 254 |  |  | virtual ~NPTfm() {}; | 
| 255 |  |  |  | 
| 256 |  |  | virtual void integrateStep( int calcPot, int calcStress ){ | 
| 257 |  |  | calcStress = 1; | 
| 258 |  |  | Integrator::integrateStep( calcPot, calcStress ); | 
| 259 |  |  | } | 
| 260 |  |  |  | 
| 261 |  |  | void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} | 
| 262 |  |  | void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} | 
| 263 |  |  | void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} | 
| 264 |  |  | void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} | 
| 265 |  |  |  | 
| 266 |  |  | protected: | 
| 267 |  |  |  | 
| 268 |  |  | virtual void  moveA( void ); | 
| 269 |  |  | virtual void moveB( void ); | 
| 270 |  |  |  | 
| 271 |  |  | virtual int readyCheck(); | 
| 272 |  |  |  | 
| 273 |  |  | // chi and eta are the propagated degrees of freedom | 
| 274 |  |  |  | 
| 275 |  |  | double chi; | 
| 276 |  |  | double eta[3][3]; | 
| 277 |  |  | double NkBT; | 
| 278 |  |  |  | 
| 279 |  |  | // targetTemp, targetPressure, and tauBarostat must be set. | 
| 280 |  |  | // One of qmass or tauThermostat must be set; | 
| 281 |  |  |  | 
| 282 |  |  | double targetTemp; | 
| 283 |  |  | double targetPressure; | 
| 284 |  |  | double tauThermostat; | 
| 285 |  |  | double tauBarostat; | 
| 286 |  |  |  | 
| 287 |  |  | short int have_tau_thermostat, have_tau_barostat, have_target_temp; | 
| 288 |  |  | short int have_target_pressure; | 
| 289 |  |  |  | 
| 290 |  |  | }; | 
| 291 |  |  |  | 
| 292 | mmeineke | 377 | #endif |