--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/07/22 19:54:52 645 +++ trunk/OOPSE/libmdtools/Integrator.hpp 2003/09/22 21:23:25 780 @@ -1,6 +1,8 @@ #ifndef _INTEGRATOR_H_ #define _INTEGRATOR_H_ +#include +#include #include "Atom.hpp" #include "Molecule.hpp" #include "SRI.hpp" @@ -9,7 +11,9 @@ #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 @@ -23,10 +27,10 @@ template class Integrator Integrator( SimInfo *theInfo, ForceFields* the_ff ); virtual ~Integrator(); void integrate( void ); + virtual double getConservedQuantity(void); - protected: - + virtual void integrateStep( int calcPot, int calcStress ); virtual void preMove( void ); virtual void moveA( void ); @@ -34,12 +38,18 @@ template class Integrator virtual void constrainA( void ); virtual void constrainB( void ); virtual int readyCheck( void ) { return 1; } + + virtual void resetIntegrator( void ) { } + + virtual void calcForce( int calcPot, int calcStress ); + virtual void thermalize(); + virtual void rotationPropagation( DirectionalAtom* dAtom, double ji[3] ); + 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 @@ -50,7 +60,7 @@ template class Integrator 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 @@ -87,10 +97,12 @@ template class NVT : public T { (public) public: NVT ( SimInfo *theInfo, ForceFields* the_ff); - virtual ~NVT() {} + virtual ~NVT(); void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} + void setChiTolerance(double tol) {chiTolerance = tol;} + virtual double getConservedQuantity(void); protected: @@ -99,10 +111,15 @@ template class NVT : public T { (public) virtual int readyCheck(); + virtual void resetIntegrator( void ); + // chi is a propagated degree of freedom. double chi; + //integral of chi(t)dt + double integralOfChidt; + // targetTemp must be set. tauThermostat must also be set; double targetTemp; @@ -110,26 +127,37 @@ template class NVT : public T { (public) short int have_tau_thermostat, have_target_temp; + double *oldVel; + double *oldJi; + + double chiTolerance; + short int have_chi_tolerance; + }; -template class NPTi : public T{ +template class NPT : public T{ public: - NPTi ( SimInfo *theInfo, ForceFields* the_ff); - virtual ~NPTi() {}; - + NPT ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPT(); + virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; T::integrateStep( calcPot, calcStress ); } + virtual double getConservedQuantity(void) = 0; + 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;} + void setChiTolerance(double tol) {chiTolerance = tol; have_chi_tolerance = 1;} + void setPosIterTolerance(double tol) {posIterTolerance = tol; have_pos_iter_tolerance = 1;} + void setEtaTolerance(double tol) {etaTolerance = tol; have_eta_tolerance = 1;} protected: @@ -138,12 +166,41 @@ template class NPTi : public T{ (protected virtual int readyCheck(); + virtual void resetIntegrator( void ); + + virtual void getVelScaleA( double sc[3], double vel[3] ) = 0; + virtual void getVelScaleB( double sc[3], int index ) = 0; + virtual void getPosScale(double pos[3], double COM[3], + int index, double sc[3]) = 0; + + virtual bool chiConverged( void ); + virtual bool etaConverged( void ) = 0; + + virtual void evolveChiA( void ); + virtual void evolveEtaA( void ) = 0; + virtual void evolveChiB( void ); + virtual void evolveEtaB( void ) = 0; + + virtual void scaleSimBox( void ) = 0; + + void accIntegralOfChidt(void) { integralOfChidt += dt * chi;} + // chi and eta are the propagated degrees of freedom + double oldChi; + double prevChi; double chi; - double eta; double NkBT; + double fkBT; + double tt2, tb2; + double instaTemp, instaPress, instaVol; + double press[3][3]; + + int Nparticles; + + double integralOfChidt; + // targetTemp, targetPressure, and tauBarostat must be set. // One of qmass or tauThermostat must be set; @@ -155,14 +212,53 @@ template class NPTi : public T{ (protected short int have_tau_thermostat, have_tau_barostat, have_target_temp; short int have_target_pressure; + double *oldPos; + double *oldVel; + double *oldJi; + + double chiTolerance; + short int have_chi_tolerance; + double posIterTolerance; + short int have_pos_iter_tolerance; + double etaTolerance; + short int have_eta_tolerance; + }; -template class NPTim : public T{ +template class NPTi : public T{ + +public: + NPTi( SimInfo *theInfo, ForceFields* the_ff); + ~NPTi(); + virtual double getConservedQuantity(void); + virtual void resetIntegrator(void); + +protected: + + + + virtual void evolveEtaA(void); + virtual void evolveEtaB(void); + + virtual bool etaConverged( void ); + + virtual void scaleSimBox( void ); + + virtual void getVelScaleA( double sc[3], double vel[3] ); + virtual void getVelScaleB( double sc[3], int index ); + virtual void getPosScale(double pos[3], double COM[3], + int index, double sc[3]); + + double eta, oldEta, prevEta; +}; + +template class NPTzm : public T{ + public: - NPTim ( SimInfo *theInfo, ForceFields* the_ff); - virtual ~NPTim() {}; + NPTzm ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTzm() {}; virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; @@ -181,6 +277,8 @@ template class NPTim : public T{ (protecte virtual int readyCheck(); + virtual void resetIntegrator( void ); + Molecule* myMolecules; Atom** myAtoms; @@ -188,6 +286,7 @@ template class NPTim : public T{ (protecte double chi; double eta; + double etaZ; double NkBT; // targetTemp, targetPressure, and tauBarostat must be set. @@ -208,8 +307,37 @@ template class NPTf : public T{ (public) public: NPTf ( SimInfo *theInfo, ForceFields* the_ff); - virtual ~NPTf() {}; + virtual ~NPTf(); + virtual double getConservedQuantity(void); + virtual void resetIntegrator(void); + +protected: + + virtual void evolveEtaA(void); + virtual void evolveEtaB(void); + + virtual bool etaConverged( void ); + + virtual void scaleSimBox( void ); + + virtual void getVelScaleA( double sc[3], double vel[3] ); + virtual void getVelScaleB( double sc[3], int index ); + virtual void getPosScale(double pos[3], double COM[3], + int index, double sc[3]); + + double eta[3][3]; + double oldEta[3][3]; + double prevEta[3][3]; +}; + +template class NPTxym : public T{ + +public: + + NPTxym ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTxym() {}; + virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; T::integrateStep( calcPot, calcStress ); @@ -222,15 +350,22 @@ template class NPTf : public T{ (public) protected: - virtual void moveA( void ); + virtual void moveA( void ); virtual void moveB( void ); virtual int readyCheck(); + virtual void resetIntegrator( void ); + + Molecule* myMolecules; + Atom** myAtoms; + // chi and eta are the propagated degrees of freedom double chi; - double eta[3][3]; + double eta; + double etaX; + double etaY; double NkBT; // targetTemp, targetPressure, and tauBarostat must be set. @@ -246,6 +381,7 @@ template class NPTfm : public T{ }; + template class NPTfm : public T{ public: @@ -256,23 +392,85 @@ template class NPTfm : public T{ (public) virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; T::integrateStep( calcPot, calcStress ); + accIntegralOfChidt(); } + virtual double getConservedQuantity(void); + 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;} + void setChiTolerance(double tol) {chiTolerance = tol;} + void setPosIterTolerance(double tol) {posIterTolerance = tol;} protected: virtual void moveA( void ); virtual void moveB( void ); + virtual void resetIntegrator( void ); + virtual int readyCheck(); + void accIntegralOfChidt(void) { integralOfChidt += dt * chi;} + Molecule* myMolecules; Atom** myAtoms; + + // chi and eta are the propagated degrees of freedom + + double chi; + double eta[3][3]; + double NkBT; + double integralOfChidt; + // 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; + double chiTolerance; + short int have_chi_tolerance; + double posIterTolerance; + short int have_pos_iter_tolerance; + +}; + + +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;} + void setChiTolerance(double tol) {chiTolerance = tol;} + void setPosIterTolerance(double tol) {posIterTolerance = tol;} + +protected: + + virtual void moveA( void ); + virtual void moveB( void ); + + virtual int readyCheck(); + + virtual void resetIntegrator( void ); + // chi and eta are the propagated degrees of freedom double chi; @@ -289,7 +487,143 @@ template class NPTfm : public T{ (public) short int have_tau_thermostat, have_tau_barostat, have_target_temp; short int have_target_pressure; + double chiTolerance; + short int have_chi_tolerance; + double posIterTolerance; + short int have_pos_iter_tolerance; }; + +template class ZConstraint : public T { + + public: + class ForceSubtractionPolicy{ + public: + ForceSubtractionPolicy(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 : public ForceSubtractionPolicy{ + + public: + PolicyByNumber(ZConstraint* integrator) :ForceSubtractionPolicy(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 : public ForceSubtractionPolicy{ + + public: + PolicyByMass(ZConstraint* integrator) :ForceSubtractionPolicy(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(); + + ForceSubtractionPolicy* forcePolicy; //force subtraction policy + friend class ForceSubtractionPolicy; + +}; + +#endif