ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Integrator.hpp
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.

File Contents

# 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 1097 #include "StuntDouble.hpp"
8 gezelter 604 #include "Molecule.hpp"
9 mmeineke 377 #include "SRI.hpp"
10     #include "AbstractClasses.hpp"
11     #include "SimInfo.hpp"
12     #include "ForceFields.hpp"
13 mmeineke 540 #include "Thermo.hpp"
14     #include "ReadWrite.hpp"
15 tim 658 #include "ZConsWriter.hpp"
16 chrisfen 1180 #include "Restraints.hpp"
17 mmeineke 377
18 tim 658 using namespace std;
19 mmeineke 561 const double kB = 8.31451e-7;// boltzmann constant amu*Ang^2*fs^-2/K
20     const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
21 mmeineke 586 const double p_convert = 1.63882576e8; //converts amu*fs^-2*Ang^-1 -> atm
22 mmeineke 561 const int maxIteration = 300;
23     const double tol = 1.0e-6;
24    
25 tim 645 template<typename T = BaseIntegrator> class Integrator : public T {
26    
27 mmeineke 377 public:
28 mmeineke 561 Integrator( SimInfo *theInfo, ForceFields* the_ff );
29 mmeineke 548 virtual ~Integrator();
30 mmeineke 377 void integrate( void );
31 tim 763 virtual double getConservedQuantity(void);
32 tim 837 virtual string getAdditionalParameters(void);
33 mmeineke 377
34 mmeineke 540 protected:
35 mmeineke 778
36 mmeineke 540 virtual void integrateStep( int calcPot, int calcStress );
37 mmeineke 548 virtual void preMove( void );
38 mmeineke 540 virtual void moveA( void );
39     virtual void moveB( void );
40     virtual void constrainA( void );
41     virtual void constrainB( void );
42 mmeineke 559 virtual int readyCheck( void ) { return 1; }
43 tim 677
44 mmeineke 746 virtual void resetIntegrator( void ) { }
45 tim 837
46     virtual void calcForce( int calcPot, int calcStress );
47 tim 677 virtual void thermalize();
48 tim 837
49 tim 1108 virtual bool stopIntegrator() {return false;}
50    
51 gezelter 1097 virtual void rotationPropagation( StuntDouble* sd, double ji[3] );
52 mmeineke 778
53 mmeineke 540 void checkConstraints( void );
54 tim 837 void rotate( int axes1, int axes2, double angle, double j[3],
55 tim 701 double A[3][3] );
56 tim 837
57 mmeineke 377 ForceFields* myFF;
58    
59 mmeineke 540 SimInfo *info; // all the info we'll ever need
60 gezelter 1097 vector<StuntDouble*> integrableObjects;
61 mmeineke 540 int nAtoms; /* the number of atoms */
62 mmeineke 548 int oldAtoms;
63 mmeineke 540 Atom **atoms; /* array of atom pointers */
64 mmeineke 423 Molecule* molecules;
65     int nMols;
66    
67 mmeineke 548 int isConstrained; // boolean to know whether the systems contains
68 tim 701 // constraints.
69 tim 837 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 mmeineke 548 int* moving; // tells whether we are moving atom i
75     int* moved; // tells whether we have moved atom i
76 tim 837 double* oldPos; // pre constrained positions
77 mmeineke 548
78 mmeineke 540 short isFirst; /*boolean for the first time integrate is called */
79 tim 837
80 mmeineke 540 double dt;
81 mmeineke 541 double dt2;
82 gezelter 560
83 mmeineke 540 Thermo *tStats;
84     StatWriter* statOut;
85     DumpWriter* dumpOut;
86 tim 837
87 chrisfen 1180 int i; // just an int
88 mmeineke 377 };
89    
90 tim 645 typedef Integrator<BaseIntegrator> RealIntegrator;
91 mmeineke 540
92 tim 1161 // ansi instantiation
93     template class Integrator<BaseIntegrator>;
94    
95 tim 645 template<typename T> class NVE : public T {
96    
97 mmeineke 561 public:
98     NVE ( SimInfo *theInfo, ForceFields* the_ff ):
99 tim 645 T( theInfo, the_ff ){}
100 tim 837 virtual ~NVE(){}
101 tim 645 };
102 mmeineke 555
103    
104 tim 645 template<typename T> class NVT : public T {
105 mmeineke 555
106 gezelter 560 public:
107    
108 mmeineke 561 NVT ( SimInfo *theInfo, ForceFields* the_ff);
109 tim 763 virtual ~NVT();
110 mmeineke 540
111 gezelter 560 void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
112     void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
113 tim 763 void setChiTolerance(double tol) {chiTolerance = tol;}
114     virtual double getConservedQuantity(void);
115 tim 837 virtual string getAdditionalParameters(void);
116 gezelter 560
117 mmeineke 540 protected:
118 gezelter 560
119 mmeineke 561 virtual void moveA( void );
120     virtual void moveB( void );
121 mmeineke 540
122 mmeineke 561 virtual int readyCheck();
123 gezelter 560
124 mmeineke 746 virtual void resetIntegrator( void );
125    
126 gezelter 565 // chi is a propagated degree of freedom.
127 gezelter 560
128 gezelter 565 double chi;
129 gezelter 560
130 tim 763 //integral of chi(t)dt
131     double integralOfChidt;
132    
133 gezelter 565 // targetTemp must be set. tauThermostat must also be set;
134 gezelter 560
135     double targetTemp;
136     double tauThermostat;
137 tim 837
138 gezelter 565 short int have_tau_thermostat, have_target_temp;
139 gezelter 560
140 tim 763 double *oldVel;
141     double *oldJi;
142    
143     double chiTolerance;
144     short int have_chi_tolerance;
145    
146 mmeineke 540 };
147    
148    
149    
150 mmeineke 778 template<typename T> class NPT : public T{
151 tim 645
152 gezelter 560 public:
153    
154 mmeineke 778 NPT ( SimInfo *theInfo, ForceFields* the_ff);
155     virtual ~NPT();
156 tim 837
157 mmeineke 594 virtual void integrateStep( int calcPot, int calcStress ){
158     calcStress = 1;
159 tim 645 T::integrateStep( calcPot, calcStress );
160 mmeineke 594 }
161    
162 mmeineke 778 virtual double getConservedQuantity(void) = 0;
163 tim 837 virtual string getAdditionalParameters(void) = 0;
164 mmeineke 843
165     double myTauThermo( void ) { return tauThermostat; }
166     double myTauBaro( void ) { return tauBarostat; }
167    
168 gezelter 560 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 tim 763 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 gezelter 560
176     protected:
177    
178 mmeineke 561 virtual void moveA( void );
179     virtual void moveB( void );
180 gezelter 560
181 mmeineke 561 virtual int readyCheck();
182 gezelter 560
183 mmeineke 746 virtual void resetIntegrator( void );
184    
185 mmeineke 778 virtual void getVelScaleA( double sc[3], double vel[3] ) = 0;
186     virtual void getVelScaleB( double sc[3], int index ) = 0;
187 tim 837 virtual void getPosScale(double pos[3], double COM[3],
188 mmeineke 778 int index, double sc[3]) = 0;
189    
190 mmeineke 857 virtual void calcVelScale( void ) = 0;
191    
192 mmeineke 778 virtual bool chiConverged( void );
193     virtual bool etaConverged( void ) = 0;
194 tim 837
195 mmeineke 778 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 tim 763 void accIntegralOfChidt(void) { integralOfChidt += dt * chi;}
203    
204 gezelter 565 // chi and eta are the propagated degrees of freedom
205 gezelter 560
206 mmeineke 778 double oldChi;
207     double prevChi;
208 gezelter 565 double chi;
209 gezelter 574 double NkBT;
210 tim 763 double fkBT;
211 gezelter 560
212 mmeineke 778 double tt2, tb2;
213     double instaTemp, instaPress, instaVol;
214 mmeineke 780 double press[3][3];
215 mmeineke 778
216 tim 763 int Nparticles;
217    
218     double integralOfChidt;
219    
220 tim 837 // targetTemp, targetPressure, and tauBarostat must be set.
221 gezelter 560 // 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 gezelter 565 short int have_target_pressure;
230 gezelter 560
231 tim 763 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 gezelter 560 };
243    
244 mmeineke 778 template<typename T> class NPTi : public T{
245 tim 837
246 mmeineke 778 public:
247     NPTi( SimInfo *theInfo, ForceFields* the_ff);
248     ~NPTi();
249    
250     virtual double getConservedQuantity(void);
251     virtual void resetIntegrator(void);
252 tim 837 virtual string getAdditionalParameters(void);
253 mmeineke 778 protected:
254    
255    
256 tim 837
257 mmeineke 778 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 tim 837 virtual void getPosScale(double pos[3], double COM[3],
267 mmeineke 778 int index, double sc[3]);
268    
269 mmeineke 857 virtual void calcVelScale( void );
270    
271 mmeineke 778 double eta, oldEta, prevEta;
272 mmeineke 857 double vScale;
273 mmeineke 778 };
274    
275 tim 645 template<typename T> class NPTf : public T{
276 gezelter 576
277     public:
278    
279     NPTf ( SimInfo *theInfo, ForceFields* the_ff);
280 tim 767 virtual ~NPTf();
281 gezelter 576
282 tim 763 virtual double getConservedQuantity(void);
283 tim 837 virtual string getAdditionalParameters(void);
284 mmeineke 780 virtual void resetIntegrator(void);
285 mmeineke 594
286 gezelter 576 protected:
287    
288 mmeineke 780 virtual void evolveEtaA(void);
289     virtual void evolveEtaB(void);
290 gezelter 576
291 mmeineke 780 virtual bool etaConverged( void );
292 mmeineke 746
293 mmeineke 780 virtual void scaleSimBox( void );
294 gezelter 576
295 mmeineke 780 virtual void getVelScaleA( double sc[3], double vel[3] );
296     virtual void getVelScaleB( double sc[3], int index );
297 tim 837 virtual void getPosScale(double pos[3], double COM[3],
298 mmeineke 780 int index, double sc[3]);
299 tim 763
300 mmeineke 857 virtual void calcVelScale( void );
301    
302 gezelter 588 double eta[3][3];
303 mmeineke 780 double oldEta[3][3];
304     double prevEta[3][3];
305 mmeineke 857 double vScale[3][3];
306 gezelter 576 };
307    
308 mmeineke 812 template<typename T> class NPTxyz : public T{
309 mmeineke 755
310     public:
311    
312 mmeineke 812 NPTxyz ( SimInfo *theInfo, ForceFields* the_ff);
313     virtual ~NPTxyz();
314 mmeineke 755
315 mmeineke 812 virtual double getConservedQuantity(void);
316 tim 837 virtual string getAdditionalParameters(void);
317 mmeineke 812 virtual void resetIntegrator(void);
318 mmeineke 755
319     protected:
320    
321 mmeineke 812 virtual void evolveEtaA(void);
322     virtual void evolveEtaB(void);
323 mmeineke 755
324 mmeineke 812 virtual bool etaConverged( void );
325 mmeineke 755
326 mmeineke 812 virtual void scaleSimBox( void );
327 mmeineke 755
328 mmeineke 812 virtual void getVelScaleA( double sc[3], double vel[3] );
329     virtual void getVelScaleB( double sc[3], int index );
330 tim 837 virtual void getPosScale(double pos[3], double COM[3],
331 mmeineke 812 int index, double sc[3]);
332 mmeineke 755
333 mmeineke 857 virtual void calcVelScale( void );
334    
335 mmeineke 812 double eta[3][3];
336     double oldEta[3][3];
337     double prevEta[3][3];
338 mmeineke 857 double vScale[3][3];
339 mmeineke 812 };
340 mmeineke 755
341    
342 tim 658 template<typename T> class ZConstraint : public T {
343 tim 837
344     public:
345 tim 738 class ForceSubtractionPolicy{
346 tim 699 public:
347 tim 738 ForceSubtractionPolicy(ZConstraint<T>* integrator) {zconsIntegrator = integrator;}
348 tim 658
349 tim 837 virtual void update() = 0;
350 tim 699 virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0;
351     virtual double getZFOfMovingMols(Atom* atom, double totalForce) = 0;
352 tim 701 virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0;
353     virtual double getHFOfUnconsMols(Atom* atom, double totalForce) = 0;
354 tim 837
355 tim 701 protected:
356 mmeineke 788 ZConstraint<T>* zconsIntegrator;
357 tim 699 };
358    
359 tim 738 class PolicyByNumber : public ForceSubtractionPolicy{
360 gezelter 747
361 tim 699 public:
362 tim 837 PolicyByNumber(ZConstraint<T>* integrator) :ForceSubtractionPolicy(integrator) {}
363     virtual void update();
364 tim 699 virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ;
365     virtual double getZFOfMovingMols(Atom* atom, double totalForce) ;
366 tim 701 virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce);
367     virtual double getHFOfUnconsMols(Atom* atom, double totalForce);
368 tim 837
369 tim 699 private:
370 tim 763 int totNumOfMovingAtoms;
371 tim 699 };
372    
373 tim 738 class PolicyByMass : public ForceSubtractionPolicy{
374 gezelter 747
375 tim 699 public:
376 tim 837 PolicyByMass(ZConstraint<T>* integrator) :ForceSubtractionPolicy(integrator) {}
377    
378     virtual void update();
379 tim 699 virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ;
380     virtual double getZFOfMovingMols(Atom* atom, double totalForce) ;
381 tim 701 virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce);
382     virtual double getHFOfUnconsMols(Atom* atom, double totalForce);
383 tim 699
384 tim 701 private:
385     double totMassOfMovingAtoms;
386 tim 699 };
387    
388 tim 658 public:
389    
390     ZConstraint( SimInfo *theInfo, ForceFields* the_ff);
391     ~ZConstraint();
392 tim 837
393 tim 658 void setZConsTime(double time) {this->zconsTime = time;}
394     void getZConsTime() {return zconsTime;}
395 tim 837
396 tim 701 void setIndexOfAllZConsMols(vector<int> index) {indexOfAllZConsMols = index;}
397     void getIndexOfAllZConsMols() {return indexOfAllZConsMols;}
398 tim 837
399 tim 701 void setZConsOutput(const char * fileName) {zconsOutput = fileName;}
400 tim 658 string getZConsOutput() {return zconsOutput;}
401 tim 837
402 tim 677 virtual void integrate();
403 tim 658
404 tim 837
405 tim 658 #ifdef IS_MPI
406 tim 701 virtual void update(); //which is called to indicate the molecules' migration
407 mmeineke 377 #endif
408 tim 658
409 tim 837 enum ZConsState {zcsMoving, zcsFixed};
410 mmeineke 790
411     vector<Molecule*> zconsMols; //z-constraint molecules array
412     vector<ZConsState> states; //state of z-constraint molecules
413    
414    
415 tim 837
416 mmeineke 790 int totNumOfUnconsAtoms; //total number of uncontraint atoms
417     double totalMassOfUncons; //total mas of unconstraint molecules
418    
419    
420 tim 658 protected:
421    
422 mmeineke 790
423 tim 837
424     virtual void calcForce( int calcPot, int calcStress );
425 tim 677 virtual void thermalize(void);
426 tim 837
427 tim 677 void zeroOutVel();
428     void doZconstraintForce();
429 tim 1093 void doHarmonic(vector<double>& resPos);
430 tim 677 bool checkZConsState();
431    
432     bool haveFixedZMols();
433     bool haveMovingZMols();
434    
435     double calcZSys();
436    
437     int isZConstraintMol(Molecule* mol);
438    
439    
440 tim 701 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 mmeineke 790
445 tim 837
446     vector<double> massOfZConsMols; //mass of z-constraint molecule
447 tim 701 vector<double> kz; //force constant array
448 mmeineke 790
449 tim 701 vector<double> zPos; //
450 tim 837
451    
452 tim 701 vector<Molecule*> unconsMols; //unconstraint molecules array
453     vector<double> massOfUnconsMols; //mass array of unconstraint molecules
454 tim 682
455 mmeineke 790
456 tim 701 vector<ZConsParaItem>* parameters; //
457 tim 837
458 tim 658 vector<int> indexOfAllZConsMols; //index of All Z-Constraint Molecuels
459 tim 677
460 tim 1093 vector<int> indexOfZConsMols; //index of local Z-Constraint Molecules
461     vector<double> fz;
462     vector<double> curZPos;
463 tim 677
464 tim 1093 bool usingSMD;
465 tim 1141 vector<double> prevCantPos;
466 tim 1093 vector<double> cantPos;
467     vector<double> cantVel;
468    
469 tim 1091 double zconsFixTime;
470     double zconsGap;
471     bool hasZConsGap;
472     vector<double> endFixTime;
473    
474 tim 837 int whichDirection; //constraint direction
475    
476 tim 658 private:
477 tim 837
478 tim 701 string zconsOutput; //filename of zconstraint output
479     ZConsWriter* fzOut; //z-constraint writer
480 tim 677
481 tim 837 double curZconsTime;
482 tim 699
483 tim 696 double calcMovingMolsCOMVel();
484     double calcSysCOMVel();
485     double calcTotalForce();
486 tim 1091 void updateZPos();
487 tim 1093 void updateCantPos();
488    
489 gezelter 747 ForceSubtractionPolicy* forcePolicy; //force subtraction policy
490 tim 738 friend class ForceSubtractionPolicy;
491 tim 677
492 tim 658 };
493    
494 tim 1064 /*
495 tim 1108 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 tim 1091 template<typename T> class SMD : public T{
523 tim 969 public:
524 tim 1091 SMD( SimInfo *theInfo, ForceFields* the_ff);
525     ~SMD();
526    
527 tim 1108 virtual void integrate();
528     virtual void calcForce( int calcPot, int calcStress );
529     bool stopIntegrator();
530     private:
531    
532 tim 969 };
533 tim 1108
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 tim 658 #endif