--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/10/22 21:17:32 812 +++ trunk/OOPSE/libmdtools/Integrator.hpp 2004/06/04 20:29:36 1249 @@ -4,6 +4,7 @@ #include #include #include "Atom.hpp" +#include "StuntDouble.hpp" #include "Molecule.hpp" #include "SRI.hpp" #include "AbstractClasses.hpp" @@ -12,6 +13,7 @@ #include "Thermo.hpp" #include "ReadWrite.hpp" #include "ZConsWriter.hpp" +#include "Restraints.hpp" using namespace std; const double kB = 8.31451e-7;// boltzmann constant amu*Ang^2*fs^-2/K @@ -20,7 +22,7 @@ const double tol = 1.0e-6; const int maxIteration = 300; const double tol = 1.0e-6; - +class RattleFramework; template class Integrator : public T { public: @@ -28,67 +30,75 @@ template class Integrator virtual ~Integrator(); void integrate( void ); virtual double getConservedQuantity(void); + virtual string getAdditionalParameters(void); protected: virtual void integrateStep( int calcPot, int calcStress ); - virtual void preMove( void ); + //virtual void preMove( void ); virtual void moveA( void ); virtual void moveB( void ); - virtual void constrainA( void ); - virtual void constrainB( void ); + //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 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], + virtual bool stopIntegrator() {return false;} + + virtual void rotationPropagation( StuntDouble* sd, double ji[3] ); + + //void checkConstraints( void ); + void rotate( int axes1, int axes2, double angle, double j[3], double A[3][3] ); - + ForceFields* myFF; SimInfo *info; // all the info we'll ever need + vector integrableObjects; int nAtoms; /* the number of atoms */ int oldAtoms; Atom **atoms; /* array of atom pointers */ Molecule* molecules; int nMols; - int isConstrained; // boolean to know whether the systems contains - // constraints. - int nConstrained; // counter for number of constraints - int *constrainedA; // the i of a constraint pair - int *constrainedB; // the j of a constraint pair - double *constrainedDsqr; // the square of the constraint distance - - int* moving; // tells whether we are moving atom i - int* moved; // tells whether we have moved atom i - double* oldPos; // pre constrained positions + RattleFramework* rattle; + //int isConstrained; // boolean to know whether the systems contains constraints. + //int nConstrained; // counter for number of constraints + //int *constrainedA; // the i of a constraint pair + //int *constrainedB; // the j of a constraint pair + //double *constrainedDsqr; // the square of the constraint distance + + //int* moving; // tells whether we are moving atom i + //int* moved; // tells whether we have moved atom i + //double* oldPos; // pre constrained positions + short isFirst; /*boolean for the first time integrate is called */ - + double dt; double dt2; Thermo *tStats; StatWriter* statOut; DumpWriter* dumpOut; - + }; typedef Integrator RealIntegrator; +// ansi instantiation +template class Integrator; + template class NVE : public T { public: NVE ( SimInfo *theInfo, ForceFields* the_ff ): T( theInfo, the_ff ){} - virtual ~NVE(){} + virtual ~NVE(){} }; @@ -103,6 +113,7 @@ template class NVT : public T { (public) void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} void setChiTolerance(double tol) {chiTolerance = tol;} virtual double getConservedQuantity(void); + virtual string getAdditionalParameters(void); protected: @@ -124,7 +135,7 @@ template class NVT : public T { (public) double targetTemp; double tauThermostat; - + short int have_tau_thermostat, have_target_temp; double *oldVel; @@ -143,13 +154,17 @@ template class NPT : public T{ (public) 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; + virtual string getAdditionalParameters(void) = 0; + + double myTauThermo( void ) { return tauThermostat; } + double myTauBaro( void ) { return tauBarostat; } void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;} void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;} @@ -170,12 +185,14 @@ template class NPT : public T{ (public) 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], + virtual void getPosScale(double pos[3], double COM[3], int index, double sc[3]) = 0; + virtual void calcVelScale( void ) = 0; + virtual bool chiConverged( void ); virtual bool etaConverged( void ) = 0; - + virtual void evolveChiA( void ); virtual void evolveEtaA( void ) = 0; virtual void evolveChiB( void ); @@ -201,7 +218,7 @@ template class NPT : public T{ (public) double integralOfChidt; - // targetTemp, targetPressure, and tauBarostat must be set. + // targetTemp, targetPressure, and tauBarostat must be set. // One of qmass or tauThermostat must be set; double targetTemp; @@ -226,18 +243,18 @@ template class NPTi : public T{ }; template class NPTi : public T{ - + public: NPTi( SimInfo *theInfo, ForceFields* the_ff); ~NPTi(); virtual double getConservedQuantity(void); virtual void resetIntegrator(void); - + virtual string getAdditionalParameters(void); protected: - + virtual void evolveEtaA(void); virtual void evolveEtaB(void); @@ -247,10 +264,13 @@ template class NPTi : public T{ 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], + virtual void getPosScale(double pos[3], double COM[3], int index, double sc[3]); + virtual void calcVelScale( void ); + double eta, oldEta, prevEta; + double vScale; }; template class NPTf : public T{ @@ -261,6 +281,7 @@ template class NPTf : public T{ (public) virtual ~NPTf(); virtual double getConservedQuantity(void); + virtual string getAdditionalParameters(void); virtual void resetIntegrator(void); protected: @@ -274,12 +295,15 @@ template class NPTf : public T{ (public) 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], + virtual void getPosScale(double pos[3], double COM[3], int index, double sc[3]); + virtual void calcVelScale( void ); + double eta[3][3]; double oldEta[3][3]; double prevEta[3][3]; + double vScale[3][3]; }; template class NPTxyz : public T{ @@ -290,6 +314,7 @@ template class NPTxyz : public T{ (public) virtual ~NPTxyz(); virtual double getConservedQuantity(void); + virtual string getAdditionalParameters(void); virtual void resetIntegrator(void); protected: @@ -303,28 +328,31 @@ template class NPTxyz : public T{ (public) 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], + virtual void getPosScale(double pos[3], double COM[3], int index, double sc[3]); + virtual void calcVelScale( void ); + double eta[3][3]; double oldEta[3][3]; double prevEta[3][3]; + double vScale[3][3]; }; template class ZConstraint : public T { - - public: + + public: class ForceSubtractionPolicy{ public: ForceSubtractionPolicy(ZConstraint* integrator) {zconsIntegrator = integrator;} - virtual void update() = 0; + 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; }; @@ -332,13 +360,13 @@ template class ZConstraint : public T { class PolicyByNumber : public ForceSubtractionPolicy{ public: - PolicyByNumber(ZConstraint* integrator) :ForceSubtractionPolicy(integrator) {} - virtual void update(); + 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; }; @@ -346,9 +374,9 @@ template class ZConstraint : public T { class PolicyByMass : public ForceSubtractionPolicy{ public: - PolicyByMass(ZConstraint* integrator) :ForceSubtractionPolicy(integrator) {} - - virtual void update(); + 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); @@ -362,30 +390,30 @@ template class ZConstraint : public T { 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 - 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 @@ -393,13 +421,13 @@ template class ZConstraint : public T { protected: - - virtual void calcForce( int calcPot, int calcStress ); + + virtual void calcForce( int calcPot, int calcStress ); virtual void thermalize(void); - + void zeroOutVel(); void doZconstraintForce(); - void doHarmonic(); + void doHarmonic(vector& resPos); bool checkZConsState(); bool haveFixedZMols(); @@ -414,40 +442,50 @@ template class ZConstraint : public T { double zconsTol; //tolerance of z-contratint double zForceConst; //base force constant term //which is estimate by OOPSE - - vector massOfZConsMols; //mass of z-constraint molecule + + vector massOfZConsMols; //mass of z-constraint molecule vector kz; //force constant array vector zPos; // - - + + vector unconsMols; //unconstraint molecules array vector massOfUnconsMols; //mass array 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; - + vector indexOfZConsMols; //index of local Z-Constraint Molecules + vector fz; + vector curZPos; + bool usingSMD; + vector prevCantPos; + vector cantPos; + vector cantVel; - int whichDirection; //constraint direction + double zconsFixTime; + double zconsGap; + bool hasZConsGap; + vector endFixTime; + int whichDirection; //constraint direction + private: - + string zconsOutput; //filename of zconstraint output ZConsWriter* fzOut; //z-constraint writer - double curZconsTime; + double curZconsTime; double calcMovingMolsCOMVel(); double calcSysCOMVel(); double calcTotalForce(); + void updateZPos(); + void updateCantPos(); ForceSubtractionPolicy* forcePolicy; //force subtraction policy friend class ForceSubtractionPolicy;