--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/07/14 21:28:54 597 +++ trunk/OOPSE/libmdtools/Integrator.hpp 2003/08/15 19:24:13 699 @@ -1,22 +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(); @@ -32,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 @@ -63,31 +71,25 @@ 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); @@ -117,8 +119,9 @@ class NPTi : public Integrator{ }; -class NPTi : public Integrator{ +template class NPTi : public T{ + public: NPTi ( SimInfo *theInfo, ForceFields* the_ff); @@ -126,7 +129,7 @@ class NPTi : public Integrator{ (public) virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; - Integrator::integrateStep( calcPot, calcStress ); + T::integrateStep( calcPot, calcStress ); } void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} @@ -160,7 +163,7 @@ class NPTim : public Integrator{ }; -class NPTim : public Integrator{ +template class NPTim : public T{ public: @@ -169,7 +172,7 @@ class NPTim : public Integrator{ (public) virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; - Integrator::integrateStep( calcPot, calcStress ); + T::integrateStep( calcPot, calcStress ); } void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} @@ -179,11 +182,14 @@ class NPTim : public Integrator{ (public) protected: - virtual void moveA( void ); + 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; @@ -203,7 +209,7 @@ class NPTf : public Integrator{ }; -class NPTf : public Integrator{ +template class NPTf : public T{ public: @@ -212,7 +218,7 @@ class NPTf : public Integrator{ (public) virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; - Integrator::integrateStep( calcPot, calcStress ); + T::integrateStep( calcPot, calcStress ); } void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} @@ -246,7 +252,7 @@ class NPTfm : public Integrator{ }; -class NPTfm : public Integrator{ +template class NPTfm : public T{ public: @@ -255,7 +261,7 @@ class NPTfm : public Integrator{ (public) virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; - Integrator::integrateStep( calcPot, calcStress ); + T::integrateStep( calcPot, calcStress ); } void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} @@ -270,6 +276,9 @@ class NPTfm : public Integrator{ (public) virtual int readyCheck(); + Molecule* myMolecules; + Atom** myAtoms; + // chi and eta are the propagated degrees of freedom double chi; @@ -289,4 +298,134 @@ class NPTfm : public Integrator{ (public) }; +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; + 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; + double* curZPos; + + int totNumOfUnconsAtoms; + + int whichDirection; //constraint direction + +private: + + string zconsOutput; + ZConsWriter* fzOut; + + double curZconsTime; + + double calcMovingMolsCOMVel(); + double calcSysCOMVel(); + double calcTotalForce(); + + ForceSubstractionPolicy* forcePolicy; + friend class ForceSubstractionPolicy; + +}; + +#endif