| 1 | #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 | #include "Thermo.hpp" | 
| 10 | #include "ReadWrite.hpp" | 
| 11 |  | 
| 12 | 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 | const double p_convert = 1.63882576e8; //converts amu*fs^-2*Ang^-1 -> atm | 
| 15 | const int maxIteration = 300; | 
| 16 | const double tol = 1.0e-6; | 
| 17 |  | 
| 18 | class Integrator : public BaseIntegrator { | 
| 19 |  | 
| 20 | public: | 
| 21 | Integrator( SimInfo *theInfo, ForceFields* the_ff ); | 
| 22 | virtual ~Integrator(); | 
| 23 | void integrate( void ); | 
| 24 |  | 
| 25 |  | 
| 26 | protected: | 
| 27 |  | 
| 28 | virtual void integrateStep( int calcPot, int calcStress ); | 
| 29 | virtual void preMove( void ); | 
| 30 | virtual void moveA( void ); | 
| 31 | virtual void moveB( void ); | 
| 32 | virtual void constrainA( void ); | 
| 33 | virtual void constrainB( void ); | 
| 34 | virtual int  readyCheck( void ) { return 1; } | 
| 35 |  | 
| 36 | void checkConstraints( void ); | 
| 37 | void rotate( int axes1, int axes2, double angle, double j[3], | 
| 38 | double A[3][3] ); | 
| 39 |  | 
| 40 |  | 
| 41 | ForceFields* myFF; | 
| 42 |  | 
| 43 | SimInfo *info; // all the info we'll ever need | 
| 44 | int nAtoms;  /* the number of atoms */ | 
| 45 | int oldAtoms; | 
| 46 | Atom **atoms; /* array of atom pointers */ | 
| 47 | Molecule* molecules; | 
| 48 | int nMols; | 
| 49 |  | 
| 50 | 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 | double* oldPos; // pre constrained positions | 
| 60 |  | 
| 61 | short isFirst; /*boolean for the first time integrate is called */ | 
| 62 |  | 
| 63 | double dt; | 
| 64 | double dt2; | 
| 65 |  | 
| 66 | double* pos; | 
| 67 | double* vel; | 
| 68 | double* frc; | 
| 69 | double* trq; | 
| 70 | double* Amat; | 
| 71 |  | 
| 72 | Thermo *tStats; | 
| 73 | StatWriter*  statOut; | 
| 74 | DumpWriter*  dumpOut; | 
| 75 |  | 
| 76 | }; | 
| 77 |  | 
| 78 | class NVE : public Integrator{ | 
| 79 |  | 
| 80 | public: | 
| 81 | NVE ( SimInfo *theInfo, ForceFields* the_ff ): | 
| 82 | Integrator( theInfo, the_ff ){} | 
| 83 | virtual ~NVE(){} | 
| 84 |  | 
| 85 |  | 
| 86 |  | 
| 87 | }; | 
| 88 |  | 
| 89 | class NVT : public Integrator{ | 
| 90 |  | 
| 91 | public: | 
| 92 |  | 
| 93 | NVT ( SimInfo *theInfo, ForceFields* the_ff); | 
| 94 | virtual ~NVT() {} | 
| 95 |  | 
| 96 | void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} | 
| 97 | void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} | 
| 98 |  | 
| 99 | protected: | 
| 100 |  | 
| 101 | virtual void moveA( void ); | 
| 102 | virtual void moveB( void ); | 
| 103 |  | 
| 104 | virtual int readyCheck(); | 
| 105 |  | 
| 106 | // chi is a propagated degree of freedom. | 
| 107 |  | 
| 108 | double chi; | 
| 109 |  | 
| 110 | // targetTemp must be set.  tauThermostat must also be set; | 
| 111 |  | 
| 112 | double targetTemp; | 
| 113 | double tauThermostat; | 
| 114 |  | 
| 115 | short int have_tau_thermostat, have_target_temp; | 
| 116 |  | 
| 117 | }; | 
| 118 |  | 
| 119 |  | 
| 120 | class NPTi : public Integrator{ | 
| 121 |  | 
| 122 | public: | 
| 123 |  | 
| 124 | NPTi ( SimInfo *theInfo, ForceFields* the_ff); | 
| 125 | virtual ~NPTi() {}; | 
| 126 |  | 
| 127 | virtual void integrateStep( int calcPot, int calcStress ){ | 
| 128 | calcStress = 1; | 
| 129 | Integrator::integrateStep( calcPot, calcStress ); | 
| 130 | } | 
| 131 |  | 
| 132 | 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 | virtual void  moveA( void ); | 
| 140 | virtual void moveB( void ); | 
| 141 |  | 
| 142 | virtual int readyCheck(); | 
| 143 |  | 
| 144 | // chi and eta are the propagated degrees of freedom | 
| 145 |  | 
| 146 | double chi; | 
| 147 | double eta; | 
| 148 | double NkBT; | 
| 149 |  | 
| 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 | short int have_target_pressure; | 
| 160 |  | 
| 161 | }; | 
| 162 |  | 
| 163 | 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 | class NPTf : public Integrator{ | 
| 207 |  | 
| 208 | public: | 
| 209 |  | 
| 210 | NPTf ( SimInfo *theInfo, ForceFields* the_ff); | 
| 211 | virtual ~NPTf() {}; | 
| 212 |  | 
| 213 | virtual void integrateStep( int calcPot, int calcStress ){ | 
| 214 | calcStress = 1; | 
| 215 | Integrator::integrateStep( calcPot, calcStress ); | 
| 216 | } | 
| 217 |  | 
| 218 | 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 | double eta[3][3]; | 
| 234 | 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 | 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 | #endif |