--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/07/31 15:35:07 658 +++ trunk/OOPSE/libmdtools/Integrator.hpp 2003/09/05 21:28:52 747 @@ -38,12 +38,16 @@ 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(); 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 @@ -54,7 +58,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 @@ -102,6 +106,8 @@ template class NVT : public T { (protected virtual void moveB( void ); virtual int readyCheck(); + + virtual void resetIntegrator( void ); // chi is a propagated degree of freedom. @@ -142,6 +148,8 @@ template class NPTi : public T{ (protected virtual int readyCheck(); + virtual void resetIntegrator( void ); + // chi and eta are the propagated degrees of freedom double chi; @@ -185,6 +193,8 @@ template class NPTim : public T{ (protecte virtual int readyCheck(); + virtual void resetIntegrator( void ); + Molecule* myMolecules; Atom** myAtoms; @@ -229,6 +239,8 @@ template class NPTf : public T{ (protected virtual void moveA( void ); virtual void moveB( void ); + virtual void resetIntegrator( void ); + virtual int readyCheck(); // chi and eta are the propagated degrees of freedom @@ -272,6 +284,8 @@ template class NPTfm : public T{ (protecte virtual void moveA( void ); virtual void moveB( void ); + virtual void resetIntegrator( void ); + virtual int readyCheck(); Molecule* myMolecules; @@ -296,55 +310,182 @@ template class ZConstraint : public T { }; -template class ZConstraint : public T { +template class NPTpr : public T{ + public: - ZConstraint( SimInfo *theInfo, ForceFields* the_ff); - ~ZConstraint(); + NPTpr ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTpr() {}; - virtual void integrateStep( int calcPot, int calcStress ); - + 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(); + + virtual void resetIntegrator( void ); + + // 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 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 setIndexOfAllZConsMols(vector index) {indexOfAllZConsMols = index;} + void getIndexOfAllZConsMols() {return indexOfAllZConsMols;} - void setZConsOutput(const char * fileName) {zconsOutput = fileName;} + 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 + virtual void update(); //which is called to indicate the molecules' migration #endif protected: - double zconsTime; + enum ZConsState {zcsMoving, zcsFixed}; + + virtual void calcForce( int calcPot, int calcStress ); + virtual void thermalize(void); - void resetZ(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; - vector massOfZConsMols; + 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; - vector massOfUnconsMols; - double totalMassOfUncons; + + vector unconsMols; //unconstraint molecules array + vector massOfUnconsMols; //mass array of unconstraint molecules + double totalMassOfUncons; //total mas of unconstraint molecules - vector allRefZ; - vector refZ; + vector* parameters; // vector indexOfAllZConsMols; //index of All Z-Constraint Molecuels - int* indexOfZConsMols; //index of local Z-Constraint Molecules - + + 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 - int isZConstraintMol(Molecule* mol); - string zconsOutput; - ZConsWriter* fzOut; + double curZconsTime; + + double calcMovingMolsCOMVel(); + double calcSysCOMVel(); + double calcTotalForce(); + + ForceSubtractionPolicy* forcePolicy; //force subtraction policy + friend class ForceSubtractionPolicy; + }; #endif