OOPSE/libmdtools/Integrator.hpp

#ifndef _INTEGRATOR_H_
#define _INTEGRATOR_H_

#include "Atom.hpp"
#include "SRI.hpp"
#include "AbstractClasses.hpp"
#include "SimInfo.hpp"
#include "ForceFields.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"

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;

class Integrator : public BaseIntegrator {

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;

  double* pos;
  double* vel;
  double* frc;
  double* trq;
  double* Amat;

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

class NVE : public Integrator{

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

  
};

class NVT : public Integrator{

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;

};


class NPTi : public Integrator{

public:

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

  virtual void integrateStep( int calcPot, int calcStress ){
    calcStress = 1;
    Integrator::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;

};

class NPTf : public Integrator{

public:

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

  virtual void integrateStep( int calcPot, int calcStress ){
    calcStress = 1;
    Integrator::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;

};

#endif
Revision:	594
Committed:	Fri Jul 11 22:34:48 2003 UTC (20 years, 11 months ago) by mmeineke
File size:	4933 byte(s)
Log Message:	working on som integrator bugs
#	Content
1	#ifndef _INTEGRATOR_H_
2	#define _INTEGRATOR_H_
3
4	#include "Atom.hpp"
5	#include "SRI.hpp"
6	#include "AbstractClasses.hpp"
7	#include "SimInfo.hpp"
8	#include "ForceFields.hpp"
9	#include "Thermo.hpp"
10	#include "ReadWrite.hpp"
11
12	const double kB = 8.31451e-7;// boltzmann constant amuAng^2fs^-2/K
13	const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
14	const double p_convert = 1.63882576e8; //converts amufs^-2Ang^-1 -> atm
15	const int maxIteration = 300;
16	const double tol = 1.0e-6;
17
18	class Integrator : public BaseIntegrator {
19
20	public:
21	Integrator( SimInfo theInfo, ForceFields the_ff );
22	virtual ~Integrator();
23	void integrate( void );
24
25
26	protected:
27
28	virtual void integrateStep( int calcPot, int calcStress );
29	virtual void preMove( void );
30	virtual void moveA( void );
31	virtual void moveB( void );
32	virtual void constrainA( void );
33	virtual void constrainB( void );
34	virtual int readyCheck( void ) { return 1; }
35
36	void checkConstraints( void );
37	void rotate( int axes1, int axes2, double angle, double j[3],
38	double A[3][3] );
39
40
41	ForceFields* myFF;
42
43	SimInfo *info; // all the info we'll ever need
44	int nAtoms; /* the number of atoms */
45	int oldAtoms;
46	Atom *atoms; / array of atom pointers */
47	Molecule* molecules;
48	int nMols;
49
50	int isConstrained; // boolean to know whether the systems contains
51	// constraints.
52	int nConstrained; // counter for number of constraints
53	int *constrainedA; // the i of a constraint pair
54	int *constrainedB; // the j of a constraint pair
55	double *constrainedDsqr; // the square of the constraint distance
56
57	int* moving; // tells whether we are moving atom i
58	int* moved; // tells whether we have moved atom i
59	double* oldPos; // pre constrained positions
60
61	short isFirst; /boolean for the first time integrate is called /
62
63	double dt;
64	double dt2;
65
66	double* pos;
67	double* vel;
68	double* frc;
69	double* trq;
70	double* Amat;
71
72	Thermo *tStats;
73	StatWriter* statOut;
74	DumpWriter* dumpOut;
75
76	};
77
78	class NVE : public Integrator{
79
80	public:
81	NVE ( SimInfo theInfo, ForceFields the_ff ):
82	Integrator( theInfo, the_ff ){}
83	virtual ~NVE(){}
84
85
86
87	};
88
89	class NVT : public Integrator{
90
91	public:
92
93	NVT ( SimInfo theInfo, ForceFields the_ff);
94	virtual ~NVT() {}
95
96	void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
97	void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
98
99	protected:
100
101	virtual void moveA( void );
102	virtual void moveB( void );
103
104	virtual int readyCheck();
105
106	// chi is a propagated degree of freedom.
107
108	double chi;
109
110	// targetTemp must be set. tauThermostat must also be set;
111
112	double targetTemp;
113	double tauThermostat;
114
115	short int have_tau_thermostat, have_target_temp;
116
117	};
118
119
120	class NPTi : public Integrator{
121
122	public:
123
124	NPTi ( SimInfo theInfo, ForceFields the_ff);
125	virtual ~NPTi() {};
126
127	virtual void integrateStep( int calcPot, int calcStress ){
128	calcStress = 1;
129	Integrator::integrateStep( calcPot, calcStress );
130	}
131
132	void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
133	void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
134	void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
135	void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
136
137	protected:
138
139	virtual void moveA( void );
140	virtual void moveB( void );
141
142	virtual int readyCheck();
143
144	// chi and eta are the propagated degrees of freedom
145
146	double chi;
147	double eta;
148	double NkBT;
149
150	// targetTemp, targetPressure, and tauBarostat must be set.
151	// One of qmass or tauThermostat must be set;
152
153	double targetTemp;
154	double targetPressure;
155	double tauThermostat;
156	double tauBarostat;
157
158	short int have_tau_thermostat, have_tau_barostat, have_target_temp;
159	short int have_target_pressure;
160
161	};
162
163	class NPTf : public Integrator{
164
165	public:
166
167	NPTf ( SimInfo theInfo, ForceFields the_ff);
168	virtual ~NPTf() {};
169
170	virtual void integrateStep( int calcPot, int calcStress ){
171	calcStress = 1;
172	Integrator::integrateStep( calcPot, calcStress );
173	}
174
175	void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
176	void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
177	void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
178	void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
179
180	protected:
181
182	virtual void moveA( void );
183	virtual void moveB( void );
184
185	virtual int readyCheck();
186
187	// chi and eta are the propagated degrees of freedom
188
189	double chi;
190	double eta[3][3];
191	double NkBT;
192
193	// targetTemp, targetPressure, and tauBarostat must be set.
194	// One of qmass or tauThermostat must be set;
195
196	double targetTemp;
197	double targetPressure;
198	double tauThermostat;
199	double tauBarostat;
200
201	short int have_tau_thermostat, have_tau_barostat, have_target_temp;
202	short int have_target_pressure;
203
204	};
205
206	#endif