--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/08/14 16:16:39 696 +++ trunk/OOPSE/libmdtools/Integrator.hpp 2003/09/29 21:16:11 790 @@ -27,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 ); @@ -39,13 +39,17 @@ template class Integrator 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 @@ -56,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 @@ -93,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: @@ -105,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; @@ -116,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: @@ -144,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; @@ -161,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; @@ -187,6 +277,8 @@ template class NPTim : public T{ (protecte virtual int readyCheck(); + virtual void resetIntegrator( void ); + Molecule* myMolecules; Atom** myAtoms; @@ -194,6 +286,7 @@ template class NPTim : public T{ (protecte double chi; double eta; + double etaZ; double NkBT; // targetTemp, targetPressure, and tauBarostat must be set. @@ -214,50 +307,36 @@ template class NPTf : public T{ (public) public: NPTf ( SimInfo *theInfo, ForceFields* the_ff); - virtual ~NPTf() {}; + virtual ~NPTf(); - virtual void integrateStep( int calcPot, int calcStress ){ - calcStress = 1; - T::integrateStep( calcPot, calcStress ); - } + virtual double getConservedQuantity(void); + virtual void resetIntegrator(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;} - protected: - virtual void moveA( void ); - virtual void moveB( void ); + virtual void evolveEtaA(void); + virtual void evolveEtaB(void); - virtual int readyCheck(); + virtual bool etaConverged( void ); - // chi and eta are the propagated degrees of freedom + virtual void scaleSimBox( void ); - double chi; - double eta[3][3]; - double NkBT; + 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]); - // 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 eta[3][3]; + double oldEta[3][3]; + double prevEta[3][3]; }; -template class NPTfm : public T{ +template class NPTxym : public T{ public: - NPTfm ( SimInfo *theInfo, ForceFields* the_ff); - virtual ~NPTfm() {}; + NPTxym ( SimInfo *theInfo, ForceFields* the_ff); + virtual ~NPTxym() {}; virtual void integrateStep( int calcPot, int calcStress ){ calcStress = 1; @@ -271,18 +350,22 @@ template class NPTfm : public T{ (protecte 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. @@ -299,7 +382,51 @@ template class ZConstraint : public T { }; 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); @@ -308,24 +435,33 @@ template class ZConstraint : public T { 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: + enum ZConsState {zcsMoving, zcsFixed}; - enum ZConsState {zcsMoving, zcsFixed}; + vector zconsMols; //z-constraint molecules array + vector states; //state of z-constraint molecules + + int totNumOfUnconsAtoms; //total number of uncontraint atoms + double totalMassOfUncons; //total mas of unconstraint molecules + + +protected: + + virtual void calcForce( int calcPot, int calcStress ); virtual void thermalize(void); @@ -343,40 +479,47 @@ template class ZConstraint : public T { int isZConstraintMol(Molecule* mol); - double zconsTime; - double zconsTol; - double zForceConst; + 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 kz; - vector states; - vector zPos; + + vector massOfZConsMols; //mass of z-constraint molecule + vector kz; //force constant array + + vector zPos; // - vector unconsMols; - vector massOfUnconsMols; - double totalMassOfUncons; + vector unconsMols; //unconstraint molecules array + vector massOfUnconsMols; //mass array of unconstraint molecules - vector* parameters; + + 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; + 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; };