--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/06/20 20:29:36 561 +++ trunk/OOPSE/libmdtools/Integrator.hpp 2003/08/14 16:16:39 696 @@ -1,21 +1,28 @@ #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" +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; -class Integrator : public BaseIntegrator { +template class Integrator : public T { + public: Integrator( SimInfo *theInfo, ForceFields* the_ff ); virtual ~Integrator(); @@ -31,12 +38,14 @@ class Integrator : public BaseIntegrator { (protected) 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[9] ); - - + double A[3][3] ); + ForceFields* myFF; SimInfo *info; // all the info we'll ever need @@ -62,37 +71,30 @@ class Integrator : public BaseIntegrator { (protected) double dt; double dt2; - double* pos; - double* vel; - double* frc; - double* trq; - double* Amat; - Thermo *tStats; StatWriter* statOut; DumpWriter* dumpOut; }; -class NVE : public Integrator{ +typedef Integrator RealIntegrator; +template class NVE : public T { + public: NVE ( SimInfo *theInfo, ForceFields* the_ff ): - Integrator( theInfo, the_ff ){} - virtual ~NVE(){} - - - + T( theInfo, the_ff ){} + virtual ~NVE(){} }; -class NVT : public Integrator{ +template class NVT : public T { + public: NVT ( SimInfo *theInfo, ForceFields* the_ff); virtual ~NVT() {} - void setQmass(double q) {qmass = q; have_qmass = 1;} void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} @@ -103,31 +105,33 @@ class NVT : public Integrator{ (protected) virtual int readyCheck(); - // zeta is a propagated degree of freedom. + // chi is a propagated degree of freedom. - double zeta; + double chi; - // targetTemp must be set. One of qmass or tauThermostat must be set; + // targetTemp must be set. tauThermostat must also be set; - double qmass; double targetTemp; double tauThermostat; - - double NkBT; - short int have_tau_thermostat, have_target_temp, have_qmass; + short int have_tau_thermostat, have_target_temp; }; -class NPT : public Integrator{ +template class NPTi : public T{ + public: - NPT ( SimInfo *theInfo, ForceFields* the_ff); - virtual ~NPT(); + NPTi ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTi() {}; - void setQmass(double q) {qmass = q; have_qmass = 1;} + 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;} @@ -140,25 +144,240 @@ class NPT : public Integrator{ (protected) virtual int readyCheck(); - Atom** atoms; + // chi and eta are the propagated degrees of freedom - // zeta and epsilonDot are the propagated degrees of freedom. + double chi; + double eta; + double NkBT; - double zeta; - double epsilonDot; + // 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 qmass; double targetTemp; double targetPressure; double tauThermostat; double tauBarostat; short int have_tau_thermostat, have_tau_barostat, have_target_temp; - short int have_target_pressure, have_qmass; + 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 calcMovingMolsCOMVel(); + double calcSysCOMVel(); + double calcTotalForce(); + +}; + +#endif