--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/06/20 16:49:33 560 +++ trunk/OOPSE/libmdtools/Integrator.hpp 2003/08/25 21:51:30 718 @@ -1,18 +1,30 @@ #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: - Integrator( SimInfo &theInfo, ForceFields* the_ff ); + Integrator( SimInfo *theInfo, ForceFields* the_ff ); virtual ~Integrator(); void integrate( void ); @@ -26,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[3][3] ); - - + double A[3][3] ); + ForceFields* myFF; SimInfo *info; // all the info we'll ever need @@ -42,7 +56,7 @@ class Integrator : public BaseIntegrator { (protected) int nMols; int isConstrained; // boolean to know whether the systems contains - // constraints. + // constraints. int nConstrained; // counter for number of constraints int *constrainedA; // the i of a constraint pair int *constrainedB; // the j of a constraint pair @@ -50,85 +64,163 @@ class Integrator : public BaseIntegrator { (protected) int* moving; // tells whether we are moving atom i int* moved; // tells whether we have moved atom i - double* prePos; // pre constrained positions + double* oldPos; // pre constrained positions short isFirst; /*boolean for the first time integrate is called */ double dt; double dt2; - const double kB = 8.31451e-7; // boltzmann constant in amu*Ang^2*fs^-2/K - const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2 - const int maxIteration = 300; - const double tol = 1.0e-6; - - double* pos; - double* vel; - double* frc; - double* trq; - double* Amat; - Thermo *tStats; StatWriter* statOut; DumpWriter* dumpOut; }; -class NVE : public Integrator{ +typedef Integrator RealIntegrator; - NVE ( void ): - Integrator( theInfo, the_ff ){} - virtual ~NVE(){} +template class NVE : public T { +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 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; }; -class NVT : public Integrator{ + +template class NPTi : public T{ + public: - NVT ( SimInfo &theInfo, ForceFields* the_ff) : - Integrator( theInfo, the_ff ); - virtual ~NVT(); + 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;} + void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;} protected: - virtual moveA( void ); - virtual moveB( void ); + virtual void moveA( void ); + virtual void moveB( void ); - int readyCheck(); + virtual int readyCheck(); - Atom** atoms; + // chi and eta are the propagated degrees of freedom - // zeta is a propagated degree of freedom. + double chi; + double eta; + double NkBT; - double zeta; + // targetTemp, targetPressure, and tauBarostat must be set. + // One of qmass or tauThermostat must be set; - // targetTemp 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_target_temp, have_qmass; + short int have_tau_thermostat, have_tau_barostat, have_target_temp; + short int have_target_pressure; }; +template class NPTim : public T{ -class NPT : public Integrator{ +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: - NPT ( SimInfo &theInfo, ForceFields* the_ff) : - Integrator( theInfo, the_ff ); - virtual ~NPT(); + NPTf ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTf() {}; - 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;} @@ -136,30 +228,248 @@ class NPT : public Integrator{ (protected) protected: - virtual moveA( void ); - virtual moveB( void ); + virtual void moveA( void ); + virtual void moveB( void ); - int readyCheck(); + 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[3][3]; + 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 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 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 NPTpr : public T{ + +public: + + NPTpr ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTpr() {}; + + 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 ZConstraint : public T { + + public: + class ForceSubstractionPolicy{ + public: + ForceSubstractionPolicy(ZConstraint* integrator) {zconsIntegrator = integrator;} + + virtual void update() = 0; + virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0; + virtual double getZFOfMovingMols(Atom* atom, double totalForce) = 0; + virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0; + virtual double getHFOfUnconsMols(Atom* atom, double totalForce) = 0; + + protected: + ZConstraint* zconsIntegrator;; + }; + + class PolicyByNumber : ForceSubstractionPolicy{ + public: + PolicyByNumber(ZConstraint* integrator) :ForceSubstractionPolicy(integrator) {} + virtual void update(); + virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ; + virtual double getZFOfMovingMols(Atom* atom, double totalForce) ; + virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce); + virtual double getHFOfUnconsMols(Atom* atom, double totalForce); + + private: + int totNumOfMovingAtoms; + }; + + class PolicyByMass :ForceSubstractionPolicy{ + public: + PolicyByMass(ZConstraint* integrator) :ForceSubstractionPolicy(integrator) {} + + virtual void update(); + virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ; + virtual double getZFOfMovingMols(Atom* atom, double totalForce) ; + virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce); + virtual double getHFOfUnconsMols(Atom* atom, double totalForce); + + private: + double totMassOfMovingAtoms; + }; + +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; //sample time + double zconsTol; //tolerance of z-contratint + double zForceConst; //base force constant term + //which is estimate by OOPSE + + vector zconsMols; //z-constraint molecules array + vector massOfZConsMols; //mass of z-constraint molecule + vector kz; //force constant array + vector states; //state of z-constraint molecules + vector zPos; // + + + vector unconsMols; //unconstraint molecules array + vector massOfUnconsMols; //mass array of unconstraint molecules + double totalMassOfUncons; //total mas of unconstraint molecules + + vector* parameters; // + + vector indexOfAllZConsMols; //index of All Z-Constraint Molecuels + + int* indexOfZConsMols; //index of local Z-Constraint Molecules + double* fz; + double* curZPos; + + int totNumOfUnconsAtoms; //total number of uncontraint atoms + + int whichDirection; //constraint direction + +private: + + string zconsOutput; //filename of zconstraint output + ZConsWriter* fzOut; //z-constraint writer + + double curZconsTime; + + double calcMovingMolsCOMVel(); + double calcSysCOMVel(); + double calcTotalForce(); + + ForceSubstractionPolicy* forcePolicy; //force substration policy + friend class ForceSubstractionPolicy; + +}; + +#endif