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"

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* prePos; // pre constrained positions 

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

  const double kB = 8.31451e-7;     // boltzmann constant in amu*Ang^2*fs^-2/K
  const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
  const int maxIteration = 300;
  const double tol = 1.0e-6;
  
  double* pos;
  double* vel;
  double* frc;
  double* trq;
  double* Amat;

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

class NVE : public Integrator{

  NVE ( void ):
    Integrator( theInfo, the_ff ){}
  virtual ~NVE(){}

  
};

class NVT : public Integrator{

public:

  NVT ( SimInfo &theInfo, ForceFields* the_ff) :
    Integrator( theInfo, the_ff );
  virtual ~NVT();

  void setQmass(double q) {qmass = q; have_qmass = 1;}
  void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
  void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}

protected:

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

  int readyCheck();

  Atom** atoms;

  // zeta is a propagated degree of freedom.

  double zeta;

  // targetTemp must be set.  One of qmass or tauThermostat must be set;

  double qmass;
  double targetTemp;
  double tauThermostat;

  short int have_tau_thermostat, have_target_temp, have_qmass;

};


class NPT : public Integrator{

public:

  NPT ( SimInfo &theInfo, ForceFields* the_ff) :
    Integrator( theInfo, the_ff );
  virtual ~NPT();

  void setQmass(double q) {qmass = q; have_qmass = 1;}
  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 moveA( void );
  virtual moveB( void );

  int readyCheck();

  Atom** atoms;

  // zeta and epsilonDot are the propagated degrees of freedom.

  double zeta;
  double epsilonDot;

  // targetTemp, targetPressure, and tauBarostat must be set.  
  // One of qmass or tauThermostat must be set;

  double qmass;
  double targetTemp;
  double targetPressure;
  double tauThermostat;
  double tauBarostat;

  short int have_tau_thermostat, have_tau_barostat, have_target_temp;
  short int have_target_pressure, have_qmass;

};

#endif
Revision:	560
Committed:	Fri Jun 20 16:49:33 2003 UTC (21 years ago) by gezelter
File size:	3845 byte(s)
Log Message:	NVT additions
#	User	Rev	Content
1	mmeineke	377	#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	mmeineke	540	#include "Thermo.hpp"
10			#include "ReadWrite.hpp"
11	mmeineke	377
12	mmeineke	540	class Integrator : public BaseIntegrator {
13	mmeineke	377
14			public:
15	mmeineke	548	Integrator( SimInfo &theInfo, ForceFields* the_ff );
16			virtual ~Integrator();
17	mmeineke	377	void integrate( void );
18
19
20	mmeineke	540	protected:
21	mmeineke	377
22	mmeineke	540	virtual void integrateStep( int calcPot, int calcStress );
23	mmeineke	548	virtual void preMove( void );
24	mmeineke	540	virtual void moveA( void );
25			virtual void moveB( void );
26			virtual void constrainA( void );
27			virtual void constrainB( void );
28	mmeineke	559	virtual int readyCheck( void ) { return 1; }
29	mmeineke	377
30	mmeineke	540	void checkConstraints( void );
31	mmeineke	377	void rotate( int axes1, int axes2, double angle, double j[3],
32			double A[3][3] );
33
34	mmeineke	540
35	mmeineke	377	ForceFields* myFF;
36
37	mmeineke	540	SimInfo *info; // all the info we'll ever need
38			int nAtoms; /* the number of atoms */
39	mmeineke	548	int oldAtoms;
40	mmeineke	540	Atom *atoms; / array of atom pointers */
41	mmeineke	423	Molecule* molecules;
42			int nMols;
43
44	mmeineke	548	int isConstrained; // boolean to know whether the systems contains
45			// constraints.
46			int nConstrained; // counter for number of constraints
47			int *constrainedA; // the i of a constraint pair
48			int *constrainedB; // the j of a constraint pair
49			double *constrainedDsqr; // the square of the constraint distance
50
51			int* moving; // tells whether we are moving atom i
52			int* moved; // tells whether we have moved atom i
53			double* prePos; // pre constrained positions
54
55	mmeineke	540	short isFirst; /boolean for the first time integrate is called /
56
57			double dt;
58	mmeineke	541	double dt2;
59	gezelter	560
60			const double kB = 8.31451e-7; // boltzmann constant in amuAng^2fs^-2/K
61	mmeineke	542	const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
62	mmeineke	548	const int maxIteration = 300;
63			const double tol = 1.0e-6;
64	mmeineke	541
65			double* pos;
66			double* vel;
67			double* frc;
68			double* trq;
69			double* Amat;
70	mmeineke	377
71	mmeineke	540	Thermo *tStats;
72			StatWriter* statOut;
73			DumpWriter* dumpOut;
74	mmeineke	377
75			};
76
77	mmeineke	555	class NVE : public Integrator{
78	mmeineke	540
79	mmeineke	555	NVE ( void ):
80			Integrator( theInfo, the_ff ){}
81			virtual ~NVE(){}
82
83
84
85			};
86
87	mmeineke	540	class NVT : public Integrator{
88
89	gezelter	560	public:
90
91			NVT ( SimInfo &theInfo, ForceFields* the_ff) :
92	mmeineke	555	Integrator( theInfo, the_ff );
93	mmeineke	540	virtual ~NVT();
94
95	gezelter	560	void setQmass(double q) {qmass = q; have_qmass = 1;}
96			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	540	virtual moveA( void );
102			virtual moveB( void );
103
104	gezelter	560	int readyCheck();
105
106			Atom** atoms;
107
108			// zeta is a propagated degree of freedom.
109
110			double zeta;
111
112			// targetTemp must be set. One of qmass or tauThermostat must be set;
113
114			double qmass;
115			double targetTemp;
116			double tauThermostat;
117
118			short int have_tau_thermostat, have_target_temp, have_qmass;
119
120	mmeineke	540	};
121
122
123	gezelter	560	class NPT : public Integrator{
124	mmeineke	540
125	gezelter	560	public:
126
127			NPT ( SimInfo &theInfo, ForceFields* the_ff) :
128			Integrator( theInfo, the_ff );
129			virtual ~NPT();
130
131			void setQmass(double q) {qmass = q; have_qmass = 1;}
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 moveA( void );
140			virtual moveB( void );
141
142			int readyCheck();
143
144			Atom** atoms;
145
146			// zeta and epsilonDot are the propagated degrees of freedom.
147
148			double zeta;
149			double epsilonDot;
150
151			// targetTemp, targetPressure, and tauBarostat must be set.
152			// One of qmass or tauThermostat must be set;
153
154			double qmass;
155			double targetTemp;
156			double targetPressure;
157			double tauThermostat;
158			double tauBarostat;
159
160			short int have_tau_thermostat, have_tau_barostat, have_target_temp;
161			short int have_target_pressure, have_qmass;
162
163			};
164
165	mmeineke	377	#endif