--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/05/30 21:31:48 542 +++ trunk/OOPSE/libmdtools/Integrator.hpp 2003/08/13 19:21:53 693 @@ -1,85 +1,382 @@ #ifndef _INTEGRATOR_H_ #define _INTEGRATOR_H_ +#include +#include #include "Atom.hpp" +#include "Molecule.hpp" #include "SRI.hpp" #include "AbstractClasses.hpp" #include "SimInfo.hpp" #include "ForceFields.hpp" #include "Thermo.hpp" #include "ReadWrite.hpp" +#include "ZConsWriter.hpp" -class Integrator : public BaseIntegrator { +using namespace std; +const double kB = 8.31451e-7;// boltzmann constant amu*Ang^2*fs^-2/K +const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2 +const double p_convert = 1.63882576e8; //converts amu*fs^-2*Ang^-1 -> atm +const int maxIteration = 300; +const double tol = 1.0e-6; + +template class Integrator : public T { + public: - Symplectic( SimInfo &theInfo, ForceFields* the_ff ); - virtual ~Symplectic(); + Integrator( SimInfo *theInfo, ForceFields* the_ff ); + virtual ~Integrator(); void integrate( void ); protected: - virtual void integrateStep( int calcPot, int calcStress ); + virtual void preMove( void ); virtual void moveA( void ); virtual void moveB( void ); virtual void constrainA( void ); virtual void constrainB( void ); + virtual int readyCheck( void ) { return 1; } + + virtual void calcForce( int calcPot, int calcStress ); + virtual void thermalize(); - void checkConstraints( void ); void rotate( int axes1, int axes2, double angle, double j[3], double A[3][3] ); - - + ForceFields* myFF; SimInfo *info; // all the info we'll ever need int nAtoms; /* the number of atoms */ + int oldAtoms; Atom **atoms; /* array of atom pointers */ Molecule* molecules; int nMols; - int isConstrained; /*boolean to know whether the systems contains - constraints. */ - int nConstrained; /*counter for number of constraints */ - int *constrainedI; /* the i of a constraint pair */ - int *constrainedJ; /* the j of a constraint pair */ - double *constrainedDsqr; /* the square of the constraint distance */ + int isConstrained; // boolean to know whether the systems contains + // constraints. + int nConstrained; // counter for number of constraints + int *constrainedA; // the i of a constraint pair + int *constrainedB; // the j of a constraint pair + double *constrainedDsqr; // the square of the constraint distance + + int* moving; // tells whether we are moving atom i + int* moved; // tells whether we have moved atom i + double* oldPos; // pre constrained positions + short isFirst; /*boolean for the first time integrate is called */ double dt; double dt2; - const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2 - - double* pos; - double* vel; - double* frc; - double* trq; - double* Amat; - - - Thermo *tStats; StatWriter* statOut; DumpWriter* dumpOut; }; +typedef Integrator RealIntegrator; -class NVT : public Integrator{ +template class NVE : public T { - NVT ( void ); - virtual ~NVT(); +public: + NVE ( SimInfo *theInfo, ForceFields* the_ff ): + T( theInfo, the_ff ){} + virtual ~NVE(){} +}; + +template class NVT : public T { + +public: + + NVT ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NVT() {} + + void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} + void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} + protected: - virtual moveA( void ); - virtual moveB( void ); + virtual void moveA( void ); + virtual void moveB( void ); + + virtual int readyCheck(); + + // chi is a propagated degree of freedom. + + double chi; + + // targetTemp must be set. tauThermostat must also be set; + + double targetTemp; + double tauThermostat; + + short int have_tau_thermostat, have_target_temp; + }; + + +template class NPTi : public T{ + +public: + + NPTi ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTi() {}; + + virtual void integrateStep( int calcPot, int calcStress ){ + calcStress = 1; + T::integrateStep( calcPot, calcStress ); + } + + void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} + void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} + void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} + void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} + +protected: + + virtual void moveA( void ); + virtual void moveB( void ); + + virtual int readyCheck(); + + // chi and eta are the propagated degrees of freedom + + double chi; + double eta; + double NkBT; + + // targetTemp, targetPressure, and tauBarostat must be set. + // One of qmass or tauThermostat must be set; + + double targetTemp; + double targetPressure; + double tauThermostat; + double tauBarostat; + + short int have_tau_thermostat, have_tau_barostat, have_target_temp; + short int have_target_pressure; + +}; + +template class NPTim : public T{ + +public: + + NPTim ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTim() {}; + + virtual void integrateStep( int calcPot, int calcStress ){ + calcStress = 1; + T::integrateStep( calcPot, calcStress ); + } + + void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} + void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} + void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} + void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} + +protected: + + virtual void moveA( void ); + virtual void moveB( void ); + + virtual int readyCheck(); + + Molecule* myMolecules; + Atom** myAtoms; + + // chi and eta are the propagated degrees of freedom + + double chi; + double eta; + double NkBT; + + // targetTemp, targetPressure, and tauBarostat must be set. + // One of qmass or tauThermostat must be set; + + double targetTemp; + double targetPressure; + double tauThermostat; + double tauBarostat; + + short int have_tau_thermostat, have_tau_barostat, have_target_temp; + short int have_target_pressure; + +}; + +template class NPTf : public T{ + +public: + + NPTf ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTf() {}; + + virtual void integrateStep( int calcPot, int calcStress ){ + calcStress = 1; + T::integrateStep( calcPot, calcStress ); + } + + void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} + void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} + void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} + void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} + +protected: + + virtual void moveA( void ); + virtual void moveB( void ); + + virtual int readyCheck(); + + // chi and eta are the propagated degrees of freedom + + double chi; + double eta[3][3]; + double NkBT; + + // targetTemp, targetPressure, and tauBarostat must be set. + // One of qmass or tauThermostat must be set; + + double targetTemp; + double targetPressure; + double tauThermostat; + double tauBarostat; + + short int have_tau_thermostat, have_tau_barostat, have_target_temp; + short int have_target_pressure; + +}; + +template class NPTfm : public T{ + +public: + + NPTfm ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTfm() {}; + + virtual void integrateStep( int calcPot, int calcStress ){ + calcStress = 1; + T::integrateStep( calcPot, calcStress ); + } + + void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} + void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} + void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} + void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} + +protected: + + virtual void moveA( void ); + virtual void moveB( void ); + + virtual int readyCheck(); + + Molecule* myMolecules; + Atom** myAtoms; + + // chi and eta are the propagated degrees of freedom + + double chi; + double eta[3][3]; + double NkBT; + + // targetTemp, targetPressure, and tauBarostat must be set. + // One of qmass or tauThermostat must be set; + + double targetTemp; + double targetPressure; + double tauThermostat; + double tauBarostat; + + short int have_tau_thermostat, have_tau_barostat, have_target_temp; + short int have_target_pressure; + +}; + +template class ZConstraint : public T { + +public: + + ZConstraint( SimInfo *theInfo, ForceFields* the_ff); + ~ZConstraint(); + + void setZConsTime(double time) {this->zconsTime = time;} + void getZConsTime() {return zconsTime;} + void setIndexOfAllZConsMols(vector index) {indexOfAllZConsMols = index;} + void getIndexOfAllZConsMols() {return indexOfAllZConsMols;} + + void setZConsOutput(const char * fileName) {zconsOutput = fileName;} + string getZConsOutput() {return zconsOutput;} + + virtual void integrate(); + +#ifdef IS_MPI + virtual void update(); //which is called to indicate the molecules' migration +#endif +protected: + + enum ZConsState {zcsMoving, zcsFixed}; + + + + virtual void calcForce( int calcPot, int calcStress ); + virtual void thermalize(void); + + void zeroOutVel(); + void doZconstraintForce(); + void doHarmonic(); + bool checkZConsState(); + + bool haveFixedZMols(); + bool haveMovingZMols(); + + double calcZSys(); + + int isZConstraintMol(Molecule* mol); + + + double zconsTime; + double zconsTol; + double zForceConst; + + vector zconsMols; + vector massOfZConsMols; + vector kz; + vector states; + vector zPos; + + + vector unconsMols; + vector massOfUnconsMols; + double totalMassOfUncons; + + vector* parameters; + + vector indexOfAllZConsMols; //index of All Z-Constraint Molecuels + + int* indexOfZConsMols; //index of local Z-Constraint Molecules + double* fz; + + int totNumOfUnconsAtoms; + + int whichDirection; //constraint direction + +private: + + string zconsOutput; + ZConsWriter* fzOut; + + double calcCOMVel(); + double calcCOMVel2(); + +}; + #endif