OOPSE/libmdtools/Integrator.hpp

#ifndef _INTEGRATOR_H_
#define _INTEGRATOR_H_

#include <string>
#include <vector>
#include "Atom.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"

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


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

typedef Integrator<BaseIntegrator> RealIntegrator;

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

protected:

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

  virtual int readyCheck();

  // chi is a propagated degree of freedom.

  double chi;

  // targetTemp must be set.  tauThermostat must also be set;

  double targetTemp;
  double tauThermostat;
  
  short int have_tau_thermostat, have_target_temp;

};


template<typename T> class NPTi : public T{

public:

  NPTi ( SimInfo *theInfo, ForceFields* the_ff);
  virtual ~NPTi() {};

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

  // chi and eta are the propagated degrees of freedom

  double chi;
  double eta;
  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<typename T> class NPTim : public T{

public:

  NPTim ( SimInfo *theInfo, ForceFields* the_ff);
  virtual ~NPTim() {};

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

  Molecule* myMolecules;
  Atom** myAtoms;

  // chi and eta are the propagated degrees of freedom

  double chi;
  double eta;
  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<typename T> class NPTf : public T{

public:

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

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

  // 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<typename T> class NPTfm : public T{

public:

  NPTfm ( SimInfo *theInfo, ForceFields* the_ff);
  virtual ~NPTfm() {};

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

  Molecule* myMolecules;
  Atom** myAtoms;

  // 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<typename T> class ZConstraint : public T {

public:

  ZConstraint( SimInfo *theInfo, ForceFields* the_ff);
  ~ZConstraint();

  virtual void integrateStep( int calcPot, int calcStress );
  
  
  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;}

#ifdef IS_MPI
  virtual void update(); //which is called to indicate the molecules' migration
#endif

protected:

  double zconsTime;
  
  void resetZ(void);
  
  vector<Molecule*> zconsMols;
  vector<double> massOfZConsMols;
  
  vector<Molecule*> unconsMols;
  vector<double> massOfUnconsMols;
  double totalMassOfUncons;

  vector<double> allRefZ;    
  vector<double> refZ;
  
  vector<int> indexOfAllZConsMols;     //index of All Z-Constraint Molecuels
  int* indexOfZConsMols;               //index of local Z-Constraint Molecules
    
  double* fz;
  
private:

  int isZConstraintMol(Molecule* mol);
  string zconsOutput;
  ZConsWriter* fzOut;
};

#endif
Revision:	658
Committed:	Thu Jul 31 15:35:07 2003 UTC (20 years, 11 months ago) by tim
File size:	8660 byte(s)
Log Message:	Added Z constraint.
#	User	Rev	Content
1	mmeineke	377	#ifndef _INTEGRATOR_H_
2			#define _INTEGRATOR_H_
3
4	tim	658	#include <string>
5			#include <vector>
6	mmeineke	377	#include "Atom.hpp"
7	gezelter	604	#include "Molecule.hpp"
8	mmeineke	377	#include "SRI.hpp"
9			#include "AbstractClasses.hpp"
10			#include "SimInfo.hpp"
11			#include "ForceFields.hpp"
12	mmeineke	540	#include "Thermo.hpp"
13			#include "ReadWrite.hpp"
14	tim	658	#include "ZConsWriter.hpp"
15	mmeineke	377
16	tim	658	using namespace std;
17	mmeineke	561	const double kB = 8.31451e-7;// boltzmann constant amuAng^2fs^-2/K
18			const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
19	mmeineke	586	const double p_convert = 1.63882576e8; //converts amufs^-2Ang^-1 -> atm
20	mmeineke	561	const int maxIteration = 300;
21			const double tol = 1.0e-6;
22
23	mmeineke	377
24	tim	645	template<typename T = BaseIntegrator> class Integrator : public T {
25
26	mmeineke	377	public:
27	mmeineke	561	Integrator( SimInfo theInfo, ForceFields the_ff );
28	mmeineke	548	virtual ~Integrator();
29	mmeineke	377	void integrate( void );
30
31
32	mmeineke	540	protected:
33	mmeineke	377
34	mmeineke	540	virtual void integrateStep( int calcPot, int calcStress );
35	mmeineke	548	virtual void preMove( void );
36	mmeineke	540	virtual void moveA( void );
37			virtual void moveB( void );
38			virtual void constrainA( void );
39			virtual void constrainB( void );
40	mmeineke	559	virtual int readyCheck( void ) { return 1; }
41	mmeineke	377
42	mmeineke	540	void checkConstraints( void );
43	mmeineke	377	void rotate( int axes1, int axes2, double angle, double j[3],
44	gezelter	600	double A[3][3] );
45	mmeineke	377
46	mmeineke	540
47	mmeineke	377	ForceFields* myFF;
48
49	mmeineke	540	SimInfo *info; // all the info we'll ever need
50			int nAtoms; /* the number of atoms */
51	mmeineke	548	int oldAtoms;
52	mmeineke	540	Atom *atoms; / array of atom pointers */
53	mmeineke	423	Molecule* molecules;
54			int nMols;
55
56	mmeineke	548	int isConstrained; // boolean to know whether the systems contains
57			// constraints.
58			int nConstrained; // counter for number of constraints
59			int *constrainedA; // the i of a constraint pair
60			int *constrainedB; // the j of a constraint pair
61			double *constrainedDsqr; // the square of the constraint distance
62
63			int* moving; // tells whether we are moving atom i
64			int* moved; // tells whether we have moved atom i
65	mmeineke	561	double* oldPos; // pre constrained positions
66	mmeineke	548
67	mmeineke	540	short isFirst; /boolean for the first time integrate is called /
68
69			double dt;
70	mmeineke	541	double dt2;
71	gezelter	560
72	mmeineke	540	Thermo *tStats;
73			StatWriter* statOut;
74			DumpWriter* dumpOut;
75	mmeineke	377
76			};
77
78	tim	645	typedef Integrator<BaseIntegrator> RealIntegrator;
79	mmeineke	540
80	tim	645	template<typename T> class NVE : public T {
81
82	mmeineke	561	public:
83			NVE ( SimInfo theInfo, ForceFields the_ff ):
84	tim	645	T( theInfo, the_ff ){}
85			virtual ~NVE(){}
86			};
87	mmeineke	555
88
89	tim	645	template<typename T> class NVT : public T {
90	mmeineke	555
91	gezelter	560	public:
92
93	mmeineke	561	NVT ( SimInfo theInfo, ForceFields the_ff);
94			virtual ~NVT() {}
95	mmeineke	540
96	gezelter	560	void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
97			void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
98
99	mmeineke	540	protected:
100	gezelter	560
101	mmeineke	561	virtual void moveA( void );
102			virtual void moveB( void );
103	mmeineke	540
104	mmeineke	561	virtual int readyCheck();
105	gezelter	560
106	gezelter	565	// chi is a propagated degree of freedom.
107	gezelter	560
108	gezelter	565	double chi;
109	gezelter	560
110	gezelter	565	// targetTemp must be set. tauThermostat must also be set;
111	gezelter	560
112			double targetTemp;
113			double tauThermostat;
114	mmeineke	561
115	gezelter	565	short int have_tau_thermostat, have_target_temp;
116	gezelter	560
117	mmeineke	540	};
118
119
120
121	tim	645	template<typename T> class NPTi : public T{
122
123	gezelter	560	public:
124
125	gezelter	574	NPTi ( SimInfo theInfo, ForceFields the_ff);
126			virtual ~NPTi() {};
127	gezelter	560
128	mmeineke	594	virtual void integrateStep( int calcPot, int calcStress ){
129			calcStress = 1;
130	tim	645	T::integrateStep( calcPot, calcStress );
131	mmeineke	594	}
132
133	gezelter	560	void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
134			void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
135			void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
136			void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
137
138			protected:
139
140	mmeineke	561	virtual void moveA( void );
141			virtual void moveB( void );
142	gezelter	560
143	mmeineke	561	virtual int readyCheck();
144	gezelter	560
145	gezelter	565	// chi and eta are the propagated degrees of freedom
146	gezelter	560
147	gezelter	565	double chi;
148			double eta;
149	gezelter	574	double NkBT;
150	gezelter	560
151			// targetTemp, targetPressure, and tauBarostat must be set.
152			// One of qmass or tauThermostat must be set;
153
154			double targetTemp;
155			double targetPressure;
156			double tauThermostat;
157			double tauBarostat;
158
159			short int have_tau_thermostat, have_tau_barostat, have_target_temp;
160	gezelter	565	short int have_target_pressure;
161	gezelter	560
162			};
163
164	tim	645	template<typename T> class NPTim : public T{
165	gezelter	596
166			public:
167
168			NPTim ( SimInfo theInfo, ForceFields the_ff);
169			virtual ~NPTim() {};
170
171			virtual void integrateStep( int calcPot, int calcStress ){
172			calcStress = 1;
173	tim	645	T::integrateStep( calcPot, calcStress );
174	gezelter	596	}
175
176			void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
177			void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
178			void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
179			void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
180
181			protected:
182
183	gezelter	604	virtual void moveA( void );
184	gezelter	596	virtual void moveB( void );
185
186			virtual int readyCheck();
187
188	gezelter	604	Molecule* myMolecules;
189			Atom** myAtoms;
190
191	gezelter	596	// chi and eta are the propagated degrees of freedom
192
193			double chi;
194			double eta;
195			double NkBT;
196
197			// targetTemp, targetPressure, and tauBarostat must be set.
198			// One of qmass or tauThermostat must be set;
199
200			double targetTemp;
201			double targetPressure;
202			double tauThermostat;
203			double tauBarostat;
204
205			short int have_tau_thermostat, have_tau_barostat, have_target_temp;
206			short int have_target_pressure;
207
208			};
209
210	tim	645	template<typename T> class NPTf : public T{
211	gezelter	576
212			public:
213
214			NPTf ( SimInfo theInfo, ForceFields the_ff);
215			virtual ~NPTf() {};
216
217	mmeineke	594	virtual void integrateStep( int calcPot, int calcStress ){
218			calcStress = 1;
219	tim	645	T::integrateStep( calcPot, calcStress );
220	mmeineke	594	}
221
222	gezelter	576	void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
223			void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
224			void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
225			void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
226
227			protected:
228
229			virtual void moveA( void );
230			virtual void moveB( void );
231
232			virtual int readyCheck();
233
234			// chi and eta are the propagated degrees of freedom
235
236			double chi;
237	gezelter	588	double eta[3][3];
238	gezelter	576	double NkBT;
239
240			// targetTemp, targetPressure, and tauBarostat must be set.
241			// One of qmass or tauThermostat must be set;
242
243			double targetTemp;
244			double targetPressure;
245			double tauThermostat;
246			double tauBarostat;
247
248			short int have_tau_thermostat, have_tau_barostat, have_target_temp;
249			short int have_target_pressure;
250
251			};
252
253	tim	645	template<typename T> class NPTfm : public T{
254	gezelter	596
255			public:
256
257			NPTfm ( SimInfo theInfo, ForceFields the_ff);
258			virtual ~NPTfm() {};
259
260			virtual void integrateStep( int calcPot, int calcStress ){
261			calcStress = 1;
262	tim	645	T::integrateStep( calcPot, calcStress );
263	gezelter	596	}
264
265			void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
266			void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
267			void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
268			void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
269
270			protected:
271
272			virtual void moveA( void );
273			virtual void moveB( void );
274
275			virtual int readyCheck();
276
277	gezelter	605	Molecule* myMolecules;
278			Atom** myAtoms;
279
280	gezelter	596	// chi and eta are the propagated degrees of freedom
281
282			double chi;
283			double eta[3][3];
284			double NkBT;
285
286			// targetTemp, targetPressure, and tauBarostat must be set.
287			// One of qmass or tauThermostat must be set;
288
289			double targetTemp;
290			double targetPressure;
291			double tauThermostat;
292			double tauBarostat;
293
294			short int have_tau_thermostat, have_tau_barostat, have_target_temp;
295			short int have_target_pressure;
296
297			};
298
299	tim	658	template<typename T> class ZConstraint : public T {
300
301			public:
302
303			ZConstraint( SimInfo theInfo, ForceFields the_ff);
304			~ZConstraint();
305
306			virtual void integrateStep( int calcPot, int calcStress );
307
308
309			void setZConsTime(double time) {this->zconsTime = time;}
310			void getZConsTime() {return zconsTime;}
311
312			void setIndexOfAllZConsMols(vector<int> index) {indexOfAllZConsMols = index;}
313			void getIndexOfAllZConsMols() {return indexOfAllZConsMols;}
314
315			void setZConsOutput(const char * fileName) {zconsOutput = fileName;}
316			string getZConsOutput() {return zconsOutput;}
317
318			#ifdef IS_MPI
319			virtual void update(); //which is called to indicate the molecules' migration
320	mmeineke	377	#endif
321	tim	658
322			protected:
323
324			double zconsTime;
325
326			void resetZ(void);
327
328			vector<Molecule*> zconsMols;
329			vector<double> massOfZConsMols;
330
331			vector<Molecule*> unconsMols;
332			vector<double> massOfUnconsMols;
333			double totalMassOfUncons;
334
335			vector<double> allRefZ;
336			vector<double> refZ;
337
338			vector<int> indexOfAllZConsMols; //index of All Z-Constraint Molecuels
339			int* indexOfZConsMols; //index of local Z-Constraint Molecules
340
341			double* fz;
342
343			private:
344
345			int isZConstraintMol(Molecule* mol);
346			string zconsOutput;
347			ZConsWriter* fzOut;
348			};
349
350			#endif