--- trunk/OOPSE/libmdtools/Integrator.hpp 2003/10/29 00:19:10 837 +++ 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: @@ -33,11 +35,11 @@ template class Integrator 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 ) { } @@ -45,32 +47,36 @@ template class Integrator virtual void calcForce( int calcPot, int calcStress ); virtual void thermalize(); - virtual void rotationPropagation( DirectionalAtom* dAtom, double ji[3] ); + virtual bool stopIntegrator() {return false;} - void checkConstraints( void ); + 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 + RattleFramework* rattle; - int* moving; // tells whether we are moving atom i - int* moved; // tells whether we have moved atom i - double* oldPos; // pre constrained positions + //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; @@ -84,6 +90,9 @@ template class NVE : public T { typedef Integrator RealIntegrator; +// ansi instantiation +template class Integrator; + template class NVE : public T { public: @@ -153,6 +162,10 @@ template class NPT : public T{ (public) 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;} void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;} @@ -175,6 +188,8 @@ template class NPT : public T{ (public) 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; @@ -252,7 +267,10 @@ template class NPTi : public T{ (protected 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{ @@ -280,9 +298,12 @@ template class NPTf : public T{ (protected 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{ @@ -310,9 +331,12 @@ template class NPTxyz : public T{ (protect 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]; }; @@ -403,7 +427,7 @@ template class ZConstraint : public T { (p void zeroOutVel(); void doZconstraintForce(); - void doHarmonic(); + void doHarmonic(vector& resPos); bool checkZConsState(); bool haveFixedZMols(); @@ -434,12 +458,20 @@ template class ZConstraint : public T { (p 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; - + double zconsFixTime; + double zconsGap; + bool hasZConsGap; + vector endFixTime; + int whichDirection; //constraint direction private: @@ -452,7 +484,9 @@ template class ZConstraint : public T { (p double calcMovingMolsCOMVel(); double calcSysCOMVel(); double calcTotalForce(); - + void updateZPos(); + void updateCantPos(); + ForceSubtractionPolicy* forcePolicy; //force subtraction policy friend class ForceSubtractionPolicy;