OOPSE/libmdtools/Integrator.hpp

#ifndef _INTEGRATOR_H_
#define _INTEGRATOR_H_

#include <string>
#include <vector>
#include "Atom.hpp"
#include "StuntDouble.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"
#include "Restraints.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;

template<typename T = BaseIntegrator> class Integrator : public T {

public:
  Integrator( SimInfo *theInfo, ForceFields* the_ff );
  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 moveA( void );
  virtual void moveB( 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 thermalize();

  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<StuntDouble*> 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

  short isFirst; /*boolean for the first time integrate is called */

  double dt;
  double dt2;

  Thermo *tStats;
  StatWriter*  statOut;
  DumpWriter*  dumpOut;

  int i; // just an int
};

typedef Integrator<BaseIntegrator> RealIntegrator;

// ansi instantiation
template class Integrator<BaseIntegrator>;

template<typename T> class NVE : public T {

public:
  NVE ( SimInfo *theInfo, ForceFields* the_ff ):
    T( theInfo, the_ff ){}
  virtual ~NVE(){}
};


template<typename T> class NVT : public T {

public:

  NVT ( SimInfo *theInfo, ForceFields* the_ff);
  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);
  virtual string getAdditionalParameters(void);

protected:

  virtual void moveA( void );
  virtual void moveB( void );

  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;
  double tauThermostat;

  short int have_tau_thermostat, have_target_temp;

  double *oldVel;
  double *oldJi;

  double chiTolerance;
  short int have_chi_tolerance;

};


template<typename T> 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;}
  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:

  virtual void  moveA( void );
  virtual void moveB( void );

  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 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 );
  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 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;

  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 *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<typename T> 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);

  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]);

  virtual void calcVelScale( void );

  double eta, oldEta, prevEta;
  double vScale;
};

template<typename T> class NPTf : public T{

public:

  NPTf ( SimInfo *theInfo, ForceFields* the_ff);
  virtual ~NPTf();

  virtual double getConservedQuantity(void);
  virtual string getAdditionalParameters(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]);

  virtual void calcVelScale( void );

  double eta[3][3];
  double oldEta[3][3];
  double prevEta[3][3];
  double vScale[3][3];
};

template<typename T> class NPTxyz : public T{

public:

  NPTxyz ( SimInfo *theInfo, ForceFields* the_ff);
  virtual ~NPTxyz();

  virtual double getConservedQuantity(void);
  virtual string getAdditionalParameters(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]);

  virtual void calcVelScale( void );

  double eta[3][3];
  double oldEta[3][3];
  double prevEta[3][3];
  double vScale[3][3];
};


template<typename T> class ZConstraint : public T {

  public:
  class ForceSubtractionPolicy{
    public:
      ForceSubtractionPolicy(ZConstraint<T>* 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<T>* zconsIntegrator;
  };

  class PolicyByNumber : public ForceSubtractionPolicy{

    public:
      PolicyByNumber(ZConstraint<T>* 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<T>* 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<int> 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};

  vector<Molecule*> zconsMols;              //z-constraint molecules array
  vector<ZConsState> 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);

  void zeroOutVel();
  void doZconstraintForce();
  void doHarmonic(vector<double>& resPos);
  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<double> massOfZConsMols;       //mass of z-constraint molecule
  vector<double> kz;                              //force constant array

  vector<double> zPos;                          //


  vector<Molecule*> unconsMols;           //unconstraint molecules array
  vector<double> massOfUnconsMols;    //mass array of unconstraint molecules


  vector<ZConsParaItem>* parameters; //

  vector<int> indexOfAllZConsMols;     //index of All Z-Constraint Molecuels

  vector<int> indexOfZConsMols;                   //index of local Z-Constraint Molecules
  vector<double> fz;
  vector<double> curZPos;

  bool usingSMD;
  vector<double> prevCantPos;
  vector<double> cantPos;
  vector<double> cantVel;

  double zconsFixTime;  
  double zconsGap;
  bool hasZConsGap;
  vector<double> endFixTime;
  
  int whichDirection;                           //constraint direction

private:

  string zconsOutput;                         //filename of zconstraint output
  ZConsWriter* fzOut;                         //z-constraint writer

  double curZconsTime;

  double calcMovingMolsCOMVel();
  double calcSysCOMVel();
  double calcTotalForce();
  void updateZPos();
  void updateCantPos();
  
  ForceSubtractionPolicy* forcePolicy; //force subtraction policy
  friend class ForceSubtractionPolicy;

};

/*
template<typename T> class SingleZConstrain : public T{


};
*/

template<typename T> class NonEquMD : public T {
  public:
    

};


//
template<typename T> class SingleZConstraint : public T{
  public:
    SingleZConstraint(SimInfo *theInfo, ForceFields* the_ff);
    ~SingleZConstraint();
    
    bool stopIntegrator();
    
  protected:
    
};

//Steered Molecular Dynamics, curret implement only support one steered molecule
 template<typename T> class SMD : public T{
  public:
    SMD( SimInfo *theInfo, ForceFields* the_ff);
    ~SMD();
  
    virtual void integrate();
    virtual void calcForce( int calcPot, int calcStress );  
    bool stopIntegrator();
  private:
    
};

//By using state pattern, Coordinate Drive is responsible for switching back and forth between
//Driven Molecular Dynamics and ZConstraint Method. 
template<typename T> class CoordinateDriver : public T {
  public:
    typedef T ParentIntegrator;

    CoordinateDriver(SimInfo*, ForceFields*, BaseIntegrator*, BaseIntegrator*);
    ~CoordinateDriver();
    
    virtual void integrate();

  private:    
    BaseIntegrator* zconsIntegrator;
    BaseIntegrator* drivenIntegrator;
    
};

#endif
Revision:	1180
Committed:	Thu May 20 20:24:07 2004 UTC (20 years, 2 months ago) by chrisfen
File size:	14406 byte(s)
Log Message:	Several additions... Restraints.cpp and .hpp were included for restraining particles in thermodynamic integration. By including these, changes were made in Integrator, SimInfo, ForceFeilds, SimSetup, StatWriter, and possibly some other files. Two bass keywords were also added for performing thermodynamic integration: a lambda value one and a k power one.
#	Content
1	#ifndef _INTEGRATOR_H_
2	#define _INTEGRATOR_H_
3
4	#include <string>
5	#include <vector>
6	#include "Atom.hpp"
7	#include "StuntDouble.hpp"
8	#include "Molecule.hpp"
9	#include "SRI.hpp"
10	#include "AbstractClasses.hpp"
11	#include "SimInfo.hpp"
12	#include "ForceFields.hpp"
13	#include "Thermo.hpp"
14	#include "ReadWrite.hpp"
15	#include "ZConsWriter.hpp"
16	#include "Restraints.hpp"
17
18	using namespace std;
19	const double kB = 8.31451e-7;// boltzmann constant amuAng^2fs^-2/K
20	const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
21	const double p_convert = 1.63882576e8; //converts amufs^-2Ang^-1 -> atm
22	const int maxIteration = 300;
23	const double tol = 1.0e-6;
24
25	template<typename T = BaseIntegrator> class Integrator : public T {
26
27	public:
28	Integrator( SimInfo theInfo, ForceFields the_ff );
29	virtual ~Integrator();
30	void integrate( void );
31	virtual double getConservedQuantity(void);
32	virtual string getAdditionalParameters(void);
33
34	protected:
35
36	virtual void integrateStep( int calcPot, int calcStress );
37	virtual void preMove( void );
38	virtual void moveA( void );
39	virtual void moveB( void );
40	virtual void constrainA( void );
41	virtual void constrainB( void );
42	virtual int readyCheck( void ) { return 1; }
43
44	virtual void resetIntegrator( void ) { }
45
46	virtual void calcForce( int calcPot, int calcStress );
47	virtual void thermalize();
48
49	virtual bool stopIntegrator() {return false;}
50
51	virtual void rotationPropagation( StuntDouble* sd, double ji[3] );
52
53	void checkConstraints( void );
54	void rotate( int axes1, int axes2, double angle, double j[3],
55	double A[3][3] );
56
57	ForceFields* myFF;
58
59	SimInfo *info; // all the info we'll ever need
60	vector<StuntDouble*> integrableObjects;
61	int nAtoms; /* the number of atoms */
62	int oldAtoms;
63	Atom *atoms; / array of atom pointers */
64	Molecule* molecules;
65	int nMols;
66
67	int isConstrained; // boolean to know whether the systems contains
68	// constraints.
69	int nConstrained; // counter for number of constraints
70	int *constrainedA; // the i of a constraint pair
71	int *constrainedB; // the j of a constraint pair
72	double *constrainedDsqr; // the square of the constraint distance
73
74	int* moving; // tells whether we are moving atom i
75	int* moved; // tells whether we have moved atom i
76	double* oldPos; // pre constrained positions
77
78	short isFirst; /boolean for the first time integrate is called /
79
80	double dt;
81	double dt2;
82
83	Thermo *tStats;
84	StatWriter* statOut;
85	DumpWriter* dumpOut;
86
87	int i; // just an int
88	};
89
90	typedef Integrator<BaseIntegrator> RealIntegrator;
91
92	// ansi instantiation
93	template class Integrator<BaseIntegrator>;
94
95	template<typename T> class NVE : public T {
96
97	public:
98	NVE ( SimInfo theInfo, ForceFields the_ff ):
99	T( theInfo, the_ff ){}
100	virtual ~NVE(){}
101	};
102
103
104	template<typename T> class NVT : public T {
105
106	public:
107
108	NVT ( SimInfo theInfo, ForceFields the_ff);
109	virtual ~NVT();
110
111	void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
112	void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
113	void setChiTolerance(double tol) {chiTolerance = tol;}
114	virtual double getConservedQuantity(void);
115	virtual string getAdditionalParameters(void);
116
117	protected:
118
119	virtual void moveA( void );
120	virtual void moveB( void );
121
122	virtual int readyCheck();
123
124	virtual void resetIntegrator( void );
125
126	// chi is a propagated degree of freedom.
127
128	double chi;
129
130	//integral of chi(t)dt
131	double integralOfChidt;
132
133	// targetTemp must be set. tauThermostat must also be set;
134
135	double targetTemp;
136	double tauThermostat;
137
138	short int have_tau_thermostat, have_target_temp;
139
140	double *oldVel;
141	double *oldJi;
142
143	double chiTolerance;
144	short int have_chi_tolerance;
145
146	};
147
148
149
150	template<typename T> class NPT : public T{
151
152	public:
153
154	NPT ( SimInfo theInfo, ForceFields the_ff);
155	virtual ~NPT();
156
157	virtual void integrateStep( int calcPot, int calcStress ){
158	calcStress = 1;
159	T::integrateStep( calcPot, calcStress );
160	}
161
162	virtual double getConservedQuantity(void) = 0;
163	virtual string getAdditionalParameters(void) = 0;
164
165	double myTauThermo( void ) { return tauThermostat; }
166	double myTauBaro( void ) { return tauBarostat; }
167
168	void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
169	void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
170	void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
171	void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
172	void setChiTolerance(double tol) {chiTolerance = tol; have_chi_tolerance = 1;}
173	void setPosIterTolerance(double tol) {posIterTolerance = tol; have_pos_iter_tolerance = 1;}
174	void setEtaTolerance(double tol) {etaTolerance = tol; have_eta_tolerance = 1;}
175
176	protected:
177
178	virtual void moveA( void );
179	virtual void moveB( void );
180
181	virtual int readyCheck();
182
183	virtual void resetIntegrator( void );
184
185	virtual void getVelScaleA( double sc[3], double vel[3] ) = 0;
186	virtual void getVelScaleB( double sc[3], int index ) = 0;
187	virtual void getPosScale(double pos[3], double COM[3],
188	int index, double sc[3]) = 0;
189
190	virtual void calcVelScale( void ) = 0;
191
192	virtual bool chiConverged( void );
193	virtual bool etaConverged( void ) = 0;
194
195	virtual void evolveChiA( void );
196	virtual void evolveEtaA( void ) = 0;
197	virtual void evolveChiB( void );
198	virtual void evolveEtaB( void ) = 0;
199
200	virtual void scaleSimBox( void ) = 0;
201
202	void accIntegralOfChidt(void) { integralOfChidt += dt * chi;}
203
204	// chi and eta are the propagated degrees of freedom
205
206	double oldChi;
207	double prevChi;
208	double chi;
209	double NkBT;
210	double fkBT;
211
212	double tt2, tb2;
213	double instaTemp, instaPress, instaVol;
214	double press[3][3];
215
216	int Nparticles;
217
218	double integralOfChidt;
219
220	// targetTemp, targetPressure, and tauBarostat must be set.
221	// One of qmass or tauThermostat must be set;
222
223	double targetTemp;
224	double targetPressure;
225	double tauThermostat;
226	double tauBarostat;
227
228	short int have_tau_thermostat, have_tau_barostat, have_target_temp;
229	short int have_target_pressure;
230
231	double *oldPos;
232	double *oldVel;
233	double *oldJi;
234
235	double chiTolerance;
236	short int have_chi_tolerance;
237	double posIterTolerance;
238	short int have_pos_iter_tolerance;
239	double etaTolerance;
240	short int have_eta_tolerance;
241
242	};
243
244	template<typename T> class NPTi : public T{
245
246	public:
247	NPTi( SimInfo theInfo, ForceFields the_ff);
248	~NPTi();
249
250	virtual double getConservedQuantity(void);
251	virtual void resetIntegrator(void);
252	virtual string getAdditionalParameters(void);
253	protected:
254
255
256
257	virtual void evolveEtaA(void);
258	virtual void evolveEtaB(void);
259
260	virtual bool etaConverged( void );
261
262	virtual void scaleSimBox( void );
263
264	virtual void getVelScaleA( double sc[3], double vel[3] );
265	virtual void getVelScaleB( double sc[3], int index );
266	virtual void getPosScale(double pos[3], double COM[3],
267	int index, double sc[3]);
268
269	virtual void calcVelScale( void );
270
271	double eta, oldEta, prevEta;
272	double vScale;
273	};
274
275	template<typename T> class NPTf : public T{
276
277	public:
278
279	NPTf ( SimInfo theInfo, ForceFields the_ff);
280	virtual ~NPTf();
281
282	virtual double getConservedQuantity(void);
283	virtual string getAdditionalParameters(void);
284	virtual void resetIntegrator(void);
285
286	protected:
287
288	virtual void evolveEtaA(void);
289	virtual void evolveEtaB(void);
290
291	virtual bool etaConverged( void );
292
293	virtual void scaleSimBox( void );
294
295	virtual void getVelScaleA( double sc[3], double vel[3] );
296	virtual void getVelScaleB( double sc[3], int index );
297	virtual void getPosScale(double pos[3], double COM[3],
298	int index, double sc[3]);
299
300	virtual void calcVelScale( void );
301
302	double eta[3][3];
303	double oldEta[3][3];
304	double prevEta[3][3];
305	double vScale[3][3];
306	};
307
308	template<typename T> class NPTxyz : public T{
309
310	public:
311
312	NPTxyz ( SimInfo theInfo, ForceFields the_ff);
313	virtual ~NPTxyz();
314
315	virtual double getConservedQuantity(void);
316	virtual string getAdditionalParameters(void);
317	virtual void resetIntegrator(void);
318
319	protected:
320
321	virtual void evolveEtaA(void);
322	virtual void evolveEtaB(void);
323
324	virtual bool etaConverged( void );
325
326	virtual void scaleSimBox( void );
327
328	virtual void getVelScaleA( double sc[3], double vel[3] );
329	virtual void getVelScaleB( double sc[3], int index );
330	virtual void getPosScale(double pos[3], double COM[3],
331	int index, double sc[3]);
332
333	virtual void calcVelScale( void );
334
335	double eta[3][3];
336	double oldEta[3][3];
337	double prevEta[3][3];
338	double vScale[3][3];
339	};
340
341
342	template<typename T> class ZConstraint : public T {
343
344	public:
345	class ForceSubtractionPolicy{
346	public:
347	ForceSubtractionPolicy(ZConstraint<T>* integrator) {zconsIntegrator = integrator;}
348
349	virtual void update() = 0;
350	virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0;
351	virtual double getZFOfMovingMols(Atom* atom, double totalForce) = 0;
352	virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0;
353	virtual double getHFOfUnconsMols(Atom* atom, double totalForce) = 0;
354
355	protected:
356	ZConstraint<T>* zconsIntegrator;
357	};
358
359	class PolicyByNumber : public ForceSubtractionPolicy{
360
361	public:
362	PolicyByNumber(ZConstraint<T>* integrator) :ForceSubtractionPolicy(integrator) {}
363	virtual void update();
364	virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ;
365	virtual double getZFOfMovingMols(Atom* atom, double totalForce) ;
366	virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce);
367	virtual double getHFOfUnconsMols(Atom* atom, double totalForce);
368
369	private:
370	int totNumOfMovingAtoms;
371	};
372
373	class PolicyByMass : public ForceSubtractionPolicy{
374
375	public:
376	PolicyByMass(ZConstraint<T>* integrator) :ForceSubtractionPolicy(integrator) {}
377
378	virtual void update();
379	virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ;
380	virtual double getZFOfMovingMols(Atom* atom, double totalForce) ;
381	virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce);
382	virtual double getHFOfUnconsMols(Atom* atom, double totalForce);
383
384	private:
385	double totMassOfMovingAtoms;
386	};
387
388	public:
389
390	ZConstraint( SimInfo theInfo, ForceFields the_ff);
391	~ZConstraint();
392
393	void setZConsTime(double time) {this->zconsTime = time;}
394	void getZConsTime() {return zconsTime;}
395
396	void setIndexOfAllZConsMols(vector<int> index) {indexOfAllZConsMols = index;}
397	void getIndexOfAllZConsMols() {return indexOfAllZConsMols;}
398
399	void setZConsOutput(const char * fileName) {zconsOutput = fileName;}
400	string getZConsOutput() {return zconsOutput;}
401
402	virtual void integrate();
403
404
405	#ifdef IS_MPI
406	virtual void update(); //which is called to indicate the molecules' migration
407	#endif
408
409	enum ZConsState {zcsMoving, zcsFixed};
410
411	vector<Molecule*> zconsMols; //z-constraint molecules array
412	vector<ZConsState> states; //state of z-constraint molecules
413
414
415
416	int totNumOfUnconsAtoms; //total number of uncontraint atoms
417	double totalMassOfUncons; //total mas of unconstraint molecules
418
419
420	protected:
421
422
423
424	virtual void calcForce( int calcPot, int calcStress );
425	virtual void thermalize(void);
426
427	void zeroOutVel();
428	void doZconstraintForce();
429	void doHarmonic(vector<double>& resPos);
430	bool checkZConsState();
431
432	bool haveFixedZMols();
433	bool haveMovingZMols();
434
435	double calcZSys();
436
437	int isZConstraintMol(Molecule* mol);
438
439
440	double zconsTime; //sample time
441	double zconsTol; //tolerance of z-contratint
442	double zForceConst; //base force constant term
443	//which is estimate by OOPSE
444
445
446	vector<double> massOfZConsMols; //mass of z-constraint molecule
447	vector<double> kz; //force constant array
448
449	vector<double> zPos; //
450
451
452	vector<Molecule*> unconsMols; //unconstraint molecules array
453	vector<double> massOfUnconsMols; //mass array of unconstraint molecules
454
455
456	vector<ZConsParaItem>* parameters; //
457
458	vector<int> indexOfAllZConsMols; //index of All Z-Constraint Molecuels
459
460	vector<int> indexOfZConsMols; //index of local Z-Constraint Molecules
461	vector<double> fz;
462	vector<double> curZPos;
463
464	bool usingSMD;
465	vector<double> prevCantPos;
466	vector<double> cantPos;
467	vector<double> cantVel;
468
469	double zconsFixTime;
470	double zconsGap;
471	bool hasZConsGap;
472	vector<double> endFixTime;
473
474	int whichDirection; //constraint direction
475
476	private:
477
478	string zconsOutput; //filename of zconstraint output
479	ZConsWriter* fzOut; //z-constraint writer
480
481	double curZconsTime;
482
483	double calcMovingMolsCOMVel();
484	double calcSysCOMVel();
485	double calcTotalForce();
486	void updateZPos();
487	void updateCantPos();
488
489	ForceSubtractionPolicy* forcePolicy; //force subtraction policy
490	friend class ForceSubtractionPolicy;
491
492	};
493
494	/*
495	template<typename T> class SingleZConstrain : public T{
496
497
498	};
499	*/
500
501	template<typename T> class NonEquMD : public T {
502	public:
503
504
505
506	};
507
508
509	//
510	template<typename T> class SingleZConstraint : public T{
511	public:
512	SingleZConstraint(SimInfo theInfo, ForceFields the_ff);
513	~SingleZConstraint();
514
515	bool stopIntegrator();
516
517	protected:
518
519	};
520
521	//Steered Molecular Dynamics, curret implement only support one steered molecule
522	template<typename T> class SMD : public T{
523	public:
524	SMD( SimInfo theInfo, ForceFields the_ff);
525	~SMD();
526
527	virtual void integrate();
528	virtual void calcForce( int calcPot, int calcStress );
529	bool stopIntegrator();
530	private:
531
532	};
533
534	//By using state pattern, Coordinate Drive is responsible for switching back and forth between
535	//Driven Molecular Dynamics and ZConstraint Method.
536	template<typename T> class CoordinateDriver : public T {
537	public:
538	typedef T ParentIntegrator;
539
540	CoordinateDriver(SimInfo, ForceFields, BaseIntegrator, BaseIntegrator);
541	~CoordinateDriver();
542
543	virtual void integrate();
544
545	private:
546	BaseIntegrator* zconsIntegrator;
547	BaseIntegrator* drivenIntegrator;
548
549	};
550
551	#endif