ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/integrators/Integrator.hpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/integrators/Integrator.hpp (file contents):
Revision 1772 by chrisfen, Tue Nov 23 22:48:31 2004 UTC vs.
Revision 1960 by tim, Wed Jan 26 15:26:47 2005 UTC

# Line 1 | Line 1
1 < #ifndef _INTEGRATOR_H_
2 < #define _INTEGRATOR_H_
1 > /*
2 > * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 > *
4 > * The University of Notre Dame grants you ("Licensee") a
5 > * non-exclusive, royalty free, license to use, modify and
6 > * redistribute this software in source and binary code form, provided
7 > * that the following conditions are met:
8 > *
9 > * 1. Acknowledgement of the program authors must be made in any
10 > *    publication of scientific results based in part on use of the
11 > *    program.  An acceptable form of acknowledgement is citation of
12 > *    the article in which the program was described (Matthew
13 > *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 > *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 > *    Parallel Simulation Engine for Molecular Dynamics,"
16 > *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 > *
18 > * 2. Redistributions of source code must retain the above copyright
19 > *    notice, this list of conditions and the following disclaimer.
20 > *
21 > * 3. Redistributions in binary form must reproduce the above copyright
22 > *    notice, this list of conditions and the following disclaimer in the
23 > *    documentation and/or other materials provided with the
24 > *    distribution.
25 > *
26 > * This software is provided "AS IS," without a warranty of any
27 > * kind. All express or implied conditions, representations and
28 > * warranties, including any implied warranty of merchantability,
29 > * fitness for a particular purpose or non-infringement, are hereby
30 > * excluded.  The University of Notre Dame and its licensors shall not
31 > * be liable for any damages suffered by licensee as a result of
32 > * using, modifying or distributing the software or its
33 > * derivatives. In no event will the University of Notre Dame or its
34 > * licensors be liable for any lost revenue, profit or data, or for
35 > * direct, indirect, special, consequential, incidental or punitive
36 > * damages, however caused and regardless of the theory of liability,
37 > * arising out of the use of or inability to use software, even if the
38 > * University of Notre Dame has been advised of the possibility of
39 > * such damages.
40 > */
41 >
42 > /**
43 >  * @file Integrator.hpp
44 >  * @author tlin
45 >  * @date 11/08/2004
46 >  * @time 13:25am
47 >  * @version 1.0
48 >  */
49  
50 < #include <string>
51 < #include <vector>
6 < #include "primitives/Atom.hpp"
7 < #include "primitives/StuntDouble.hpp"
8 < #include "primitives/Molecule.hpp"
9 < #include "primitives/SRI.hpp"
10 < #include "primitives/AbstractClasses.hpp"
11 < #include "brains/SimInfo.hpp"
12 < #include "UseTheForce/ForceFields.hpp"
13 < #include "brains/Thermo.hpp"
14 < #include "io/ReadWrite.hpp"
15 < #include "io/ZConsWriter.hpp"
16 < #include "restraints/Restraints.hpp"
50 > #ifndef INTEGRATORS_INTEGRATOR_HPP
51 > #define INTEGRATORS_INTEGRATOR_HPP
52  
53 < using namespace std;
54 < const double kB = 8.31451e-7;// boltzmann constant amu*Ang^2*fs^-2/K
55 < const double eConvert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2
56 < const double p_convert = 1.63882576e8; //converts amu*fs^-2*Ang^-1 -> atm
22 < const int maxIteration = 300;
23 < const double tol = 1.0e-6;
53 > #include "brains/ForceManager.hpp"
54 > #include "io/DumpWriter.hpp"
55 > #include "io/StatWriter.hpp"
56 > #include "integrators/Velocitizer.hpp"
57  
58 < class VelVerletConsFramework;
26 < template<typename T = BaseIntegrator> class Integrator : public T {
58 > namespace oopse {
59  
28 – public:
29 –  Integrator( SimInfo *theInfo, ForceFields* the_ff );
30 –  virtual ~Integrator();
31 –  void integrate( void );
32 –  virtual double  getConservedQuantity(void);
33 –  virtual string getAdditionalParameters(void);
60  
61 < protected:
61 > /**
62 > * @class Integrator Integrator.hpp "integrators/Integrator.hpp"
63 > * @brief Base class of Integrator
64 > * @todo document
65 > */
66 > class Integrator {
67 >    public:
68  
69 <  virtual void integrateStep( int calcPot, int calcStress );
38 <  virtual void preMove( void );
39 <  virtual void moveA( void );
40 <  virtual void moveB( void );
41 <  virtual void constrainA( void );
42 <  virtual void constrainB( void );
43 <  virtual int  readyCheck( void ) { return 1; }
69 >        virtual ~Integrator();
70  
71 <  virtual void resetIntegrator( void ) { }
71 >        //avoid public virtual function        
72 >        void integrate() {
73 >            doIntegrate();
74 >        }
75  
76 <  virtual void calcForce( int calcPot, int calcStress );
77 <  virtual void thermalize();
76 >        void update() {
77 >            doUpdate();
78 >        }
79  
80 <  virtual bool stopIntegrator() {return false;}
80 >        void setForceManager(ForceManager* forceMan) {
81 >            if (forceMan_ != forceMan && forceMan_  != NULL) {
82 >                delete forceMan_;
83 >            }
84 >            forceMan_ = forceMan;
85 >        }
86  
87 <  virtual void rotationPropagation( StuntDouble* sd, double ji[3] );
87 >        void setVelocitizer(Velocitizer* velocitizer) {
88 >            if (velocitizer_ != velocitizer && velocitizer_  != NULL) {
89 >                delete velocitizer_;
90 >            }
91 >            velocitizer_  = velocitizer;
92 >        }
93 >        
94 >    protected:
95  
96 <  void checkConstraints( void );
55 <  void rotate( int axes1, int axes2, double angle, double j[3],
56 <         double A[3][3] );
96 >        Integrator(SimInfo* info);
97  
98 <  ForceFields* myFF;
98 >        virtual void doIntegrate() = 0;
99  
100 <  SimInfo *info; // all the info we'll ever need
101 <  vector<StuntDouble*> integrableObjects;
102 <  int nAtoms;  /* the number of atoms */
103 <  int oldAtoms;
104 <  Atom **atoms; /* array of atom pointers */
105 <  Molecule* molecules;
106 <  int nMols;
100 >        virtual void doUpdate() {}
101 >        
102 >        void saveConservedQuantity() {
103 >            currentSnapshot_->statData[Stats::CONSERVED_QUANTITY] = calcConservedQuantity();
104 >        }
105 >        
106 >        SimInfo* info_;
107 >        ForceManager* forceMan_;
108 >        bool needPotential;
109 >        bool needStress;
110 >        
111 >        Velocitizer* velocitizer_;
112 >        bool needVelocityScaling;
113 >        double targetScalingTemp;
114 >    
115 >        DumpWriter*dumpWriter;
116 >        DumpWriter* eorWriter;
117 >        StatWriter* statWriter;
118 >        Thermo thermo;
119  
120 +        double runTime;
121 +        double sampleTime;
122 +        double statusTime;
123 +        double thermalTime;
124 +        double dt;
125  
126 <  int isConstrained; // boolean to know whether the systems contains constraints.
70 <  int nConstrained;  // counter for number of constraints
71 <  int *constrainedA; // the i of a constraint pair
72 <  int *constrainedB; // the j of a constraint pair
73 <  double *constrainedDsqr; // the square of the constraint distance
126 >        Snapshot* currentSnapshot_; //During the integration, the address of currentSnapshot Will not change
127  
128 <  int* moving; // tells whether we are moving atom i
129 <  int* moved;  // tells whether we have moved atom i
130 <  double* oldPos; // pre constrained positions
131 <
132 <  short isFirst; /*boolean for the first time integrate is called */
133 <
134 <  double dt;
135 <  double dt2;
83 <
84 <  Thermo *tStats;
85 <  StatWriter*  statOut;
86 <  DumpWriter*  dumpOut;
87 <  RestraintWriter* restOut;
88 <  RestraintReader* initRestraints;
89 <
90 < };
91 <
92 < typedef Integrator<BaseIntegrator> RealIntegrator;
93 <
94 < // ansi instantiation
95 < // template class Integrator<BaseIntegrator>;
96 <
97 <
98 < template<typename T> class NVE : public T {
99 <
100 < public:
101 <  NVE ( SimInfo *theInfo, ForceFields* the_ff ):
102 <    T( theInfo, the_ff ){}
103 <  virtual ~NVE(){}
104 < };
105 <
106 <
107 < template<typename T> class NVT : public T {
108 <
109 < public:
110 <
111 <  NVT ( SimInfo *theInfo, ForceFields* the_ff);
112 <  virtual ~NVT();
113 <
114 <  void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
115 <  void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
116 <  void setChiTolerance(double tol) {chiTolerance = tol;}
117 <  virtual double  getConservedQuantity(void);
118 <  virtual string getAdditionalParameters(void);
119 <
120 < protected:
121 <
122 <  virtual void moveA( void );
123 <  virtual void moveB( void );
124 <
125 <  virtual int readyCheck();
126 <
127 <  virtual void resetIntegrator( void );
128 <
129 <  // chi is a propagated degree of freedom.
130 <
131 <  double chi;
132 <
133 <  //integral of chi(t)dt
134 <  double integralOfChidt;
135 <
136 <  // targetTemp must be set.  tauThermostat must also be set;
137 <
138 <  double targetTemp;
139 <  double tauThermostat;
140 <
141 <  short int have_tau_thermostat, have_target_temp;
142 <
143 <  double *oldVel;
144 <  double *oldJi;
145 <
146 <  double chiTolerance;
147 <  short int have_chi_tolerance;
148 <
149 < };
150 <
151 <
152 <
153 < template<typename T> class NPT : public T{
154 <
155 < public:
156 <
157 <  NPT ( SimInfo *theInfo, ForceFields* the_ff);
158 <  virtual ~NPT();
159 <
160 <  virtual void integrateStep( int calcPot, int calcStress ){
161 <    calcStress = 1;
162 <    T::integrateStep( calcPot, calcStress );
163 <  }
164 <
165 <  virtual double getConservedQuantity(void) = 0;
166 <  virtual string getAdditionalParameters(void) = 0;
167 <  
168 <  double myTauThermo( void ) { return tauThermostat; }
169 <  double myTauBaro( void ) { return tauBarostat; }
170 <
171 <  void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
172 <  void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
173 <  void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
174 <  void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
175 <  void setChiTolerance(double tol) {chiTolerance = tol; have_chi_tolerance = 1;}
176 <  void setPosIterTolerance(double tol) {posIterTolerance = tol; have_pos_iter_tolerance = 1;}
177 <  void setEtaTolerance(double tol) {etaTolerance = tol; have_eta_tolerance = 1;}
178 <
179 < protected:
180 <
181 <  virtual void  moveA( void );
182 <  virtual void moveB( void );
183 <
184 <  virtual int readyCheck();
185 <
186 <  virtual void resetIntegrator( void );
187 <
188 <  virtual void getVelScaleA( double sc[3], double vel[3] ) = 0;
189 <  virtual void getVelScaleB( double sc[3], int index ) = 0;
190 <  virtual void getPosScale(double pos[3], double COM[3],
191 <                           int index, double sc[3]) = 0;
192 <
193 <  virtual void calcVelScale( void ) = 0;
194 <
195 <  virtual bool chiConverged( void );
196 <  virtual bool etaConverged( void ) = 0;
197 <
198 <  virtual void evolveChiA( void );
199 <  virtual void evolveEtaA( void ) = 0;
200 <  virtual void evolveChiB( void );
201 <  virtual void evolveEtaB( void ) = 0;
202 <
203 <  virtual void scaleSimBox( void ) = 0;
204 <
205 <  void accIntegralOfChidt(void) { integralOfChidt += dt * chi;}
206 <
207 <  // chi and eta are the propagated degrees of freedom
128 >        
129 >    private:
130 >        
131 >        virtual double calcConservedQuantity() = 0;
132 >        
133 >        virtual DumpWriter* createDumpWriter() = 0;
134 >        
135 >        virtual DumpWriter* createEorWriter() = 0;
136  
137 <  double oldChi;
210 <  double prevChi;
211 <  double chi;
212 <  double NkBT;
213 <  double fkBT;
214 <
215 <  double tt2, tb2;
216 <  double instaTemp, instaPress, instaVol;
217 <  double press[3][3];
218 <
219 <  int Nparticles;
220 <
221 <  double integralOfChidt;
222 <
223 <  // targetTemp, targetPressure, and tauBarostat must be set.
224 <  // One of qmass or tauThermostat must be set;
225 <
226 <  double targetTemp;
227 <  double targetPressure;
228 <  double tauThermostat;
229 <  double tauBarostat;
230 <
231 <  short int have_tau_thermostat, have_tau_barostat, have_target_temp;
232 <  short int have_target_pressure;
233 <
234 <  double *oldPos;
235 <  double *oldVel;
236 <  double *oldJi;
237 <
238 <  double chiTolerance;
239 <  short int have_chi_tolerance;
240 <  double posIterTolerance;
241 <  short int have_pos_iter_tolerance;
242 <  double etaTolerance;
243 <  short int have_eta_tolerance;
244 <
137 >        virtual StatWriter* createStatWriter() = 0;
138   };
139  
247 – template<typename T> class NPTi : public T{
248 –
249 – public:
250 –  NPTi( SimInfo *theInfo, ForceFields* the_ff);
251 –  ~NPTi();
252 –
253 –  virtual double getConservedQuantity(void);
254 –  virtual void resetIntegrator(void);
255 –  virtual string getAdditionalParameters(void);
256 – protected:
257 –
258 –
259 –
260 –  virtual void evolveEtaA(void);
261 –  virtual void evolveEtaB(void);
262 –
263 –  virtual bool etaConverged( void );
264 –
265 –  virtual void scaleSimBox( void );
266 –
267 –  virtual void getVelScaleA( double sc[3], double vel[3] );
268 –  virtual void getVelScaleB( double sc[3], int index );
269 –  virtual void getPosScale(double pos[3], double COM[3],
270 –                           int index, double sc[3]);
271 –
272 –  virtual void calcVelScale( void );
273 –
274 –  double eta, oldEta, prevEta;
275 –  double vScale;
276 – };
277 –
278 – template<typename T> class NPTf : public T{
279 –
280 – public:
281 –
282 –  NPTf ( SimInfo *theInfo, ForceFields* the_ff);
283 –  virtual ~NPTf();
284 –
285 –  virtual double getConservedQuantity(void);
286 –  virtual string getAdditionalParameters(void);
287 –  virtual void resetIntegrator(void);
288 –
289 – protected:
290 –
291 –  virtual void evolveEtaA(void);
292 –  virtual void evolveEtaB(void);
293 –
294 –  virtual bool etaConverged( void );
295 –
296 –  virtual void scaleSimBox( void );
297 –
298 –  virtual void getVelScaleA( double sc[3], double vel[3] );
299 –  virtual void getVelScaleB( double sc[3], int index );
300 –  virtual void getPosScale(double pos[3], double COM[3],
301 –                           int index, double sc[3]);
302 –
303 –  virtual void calcVelScale( void );
304 –
305 –  double eta[3][3];
306 –  double oldEta[3][3];
307 –  double prevEta[3][3];
308 –  double vScale[3][3];
309 – };
310 –
311 – template<typename T> class NPTxyz : public T{
312 –
313 – public:
314 –
315 –  NPTxyz ( SimInfo *theInfo, ForceFields* the_ff);
316 –  virtual ~NPTxyz();
317 –
318 –  virtual double getConservedQuantity(void);
319 –  virtual string getAdditionalParameters(void);
320 –  virtual void resetIntegrator(void);
321 –
322 – protected:
323 –
324 –  virtual void evolveEtaA(void);
325 –  virtual void evolveEtaB(void);
326 –
327 –  virtual bool etaConverged( void );
328 –
329 –  virtual void scaleSimBox( void );
330 –
331 –  virtual void getVelScaleA( double sc[3], double vel[3] );
332 –  virtual void getVelScaleB( double sc[3], int index );
333 –  virtual void getPosScale(double pos[3], double COM[3],
334 –                           int index, double sc[3]);
335 –
336 –  virtual void calcVelScale( void );
337 –
338 –  double eta[3][3];
339 –  double oldEta[3][3];
340 –  double prevEta[3][3];
341 –  double vScale[3][3];
342 – };
343 –
344 –
345 – template<typename T> class ZConstraint : public T {
346 –
347 –  public:
348 –  class ForceSubtractionPolicy{
349 –    public:
350 –      ForceSubtractionPolicy(ZConstraint<T>* integrator) {zconsIntegrator = integrator;}
351 –
352 –      virtual void update() = 0;
353 –      virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0;
354 –      virtual double getZFOfMovingMols(Atom* atom, double totalForce) = 0;
355 –      virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0;
356 –      virtual double getHFOfUnconsMols(Atom* atom, double totalForce) = 0;
357 –
358 –   protected:
359 –     ZConstraint<T>* zconsIntegrator;
360 –  };
361 –
362 –  class PolicyByNumber : public ForceSubtractionPolicy{
363 –
364 –    public:
365 –      PolicyByNumber(ZConstraint<T>* integrator) :ForceSubtractionPolicy(integrator) {}
366 –      virtual void update();
367 –      virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ;
368 –      virtual double getZFOfMovingMols(Atom* atom, double totalForce) ;
369 –      virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce);
370 –      virtual double getHFOfUnconsMols(Atom* atom, double totalForce);
371 –
372 –    private:
373 –      int totNumOfMovingAtoms;
374 –  };
375 –
376 –  class PolicyByMass : public ForceSubtractionPolicy{
377 –
378 –    public:
379 –      PolicyByMass(ZConstraint<T>* integrator) :ForceSubtractionPolicy(integrator) {}
380 –
381 –      virtual void update();
382 –      virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ;
383 –      virtual double getZFOfMovingMols(Atom* atom, double totalForce) ;
384 –      virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce);
385 –      virtual double getHFOfUnconsMols(Atom* atom, double totalForce);
386 –
387 –   private:
388 –     double totMassOfMovingAtoms;
389 –  };
390 –
391 – public:
392 –
393 –  ZConstraint( SimInfo *theInfo, ForceFields* the_ff);
394 –  ~ZConstraint();
395 –
396 –  void setZConsTime(double time)                  {this->zconsTime = time;}
397 –  void getZConsTime()                             {return zconsTime;}
398 –
399 –  void setIndexOfAllZConsMols(vector<int> index) {indexOfAllZConsMols = index;}
400 –  void getIndexOfAllZConsMols()                  {return indexOfAllZConsMols;}
401 –
402 –  void setZConsOutput(const char * fileName)          {zconsOutput = fileName;}
403 –  string getZConsOutput()                         {return zconsOutput;}
404 –
405 –  virtual void integrate();
406 –
407 –
408 – #ifdef IS_MPI
409 –  virtual void update();                      //which is called to indicate the molecules' migration
410 – #endif
411 –
412 –  enum ZConsState {zcsMoving, zcsFixed};
413 –
414 –  vector<Molecule*> zconsMols;              //z-constraint molecules array
415 –  vector<ZConsState> states;                 //state of z-constraint molecules
416 –
417 –
418 –
419 –  int totNumOfUnconsAtoms;              //total number of uncontraint atoms
420 –  double totalMassOfUncons;                //total mas of unconstraint molecules
421 –
422 –
423 – protected:
424 –
425 –
426 –
427 –  virtual void calcForce( int calcPot, int calcStress );
428 –  virtual void thermalize(void);
429 –
430 –  void zeroOutVel();
431 –  void doZconstraintForce();
432 –  void doHarmonic(vector<double>& resPos);
433 –  bool checkZConsState();
434 –
435 –  bool haveFixedZMols();
436 –  bool haveMovingZMols();
437 –
438 –  double calcZSys();
439 –
440 –  int isZConstraintMol(Molecule* mol);
441 –
442 –
443 –  double zconsTime;                              //sample time
444 –  double zconsTol;                                 //tolerance of z-contratint
445 –  double zForceConst;                           //base force constant term
446 –                                                          //which is estimate by OOPSE
447 –
448 –
449 –  vector<double> massOfZConsMols;       //mass of z-constraint molecule
450 –  vector<double> kz;                              //force constant array
451 –
452 –  vector<double> zPos;                          //
453 –
454 –
455 –  vector<Molecule*> unconsMols;           //unconstraint molecules array
456 –  vector<double> massOfUnconsMols;    //mass array of unconstraint molecules
457 –
458 –
459 –  vector<ZConsParaItem>* parameters; //
460 –
461 –  vector<int> indexOfAllZConsMols;     //index of All Z-Constraint Molecuels
462 –
463 –  vector<int> indexOfZConsMols;                   //index of local Z-Constraint Molecules
464 –  vector<double> fz;
465 –  vector<double> curZPos;
466 –
467 –  bool usingSMD;
468 –  vector<double> prevCantPos;
469 –  vector<double> cantPos;
470 –  vector<double> cantVel;
471 –
472 –  double zconsFixTime;  
473 –  double zconsGap;
474 –  bool hasZConsGap;
475 –  vector<double> endFixTime;
476 –  
477 –  int whichDirection;                           //constraint direction
478 –
479 – private:
480 –
481 –  string zconsOutput;                         //filename of zconstraint output
482 –  ZConsWriter* fzOut;                         //z-constraint writer
483 –
484 –  double curZconsTime;
485 –
486 –  double calcMovingMolsCOMVel();
487 –  double calcSysCOMVel();
488 –  double calcTotalForce();
489 –  void updateZPos();
490 –  void updateCantPos();
491 –  
492 –  ForceSubtractionPolicy* forcePolicy; //force subtraction policy
493 –  friend class ForceSubtractionPolicy;
494 –
495 – };
496 –
497 –
498 – //Sympletic quaternion Scheme Integrator
499 – //Reference:
500 – // T.F. Miller, M. Eleftheriou, P. Pattnaik, A. Ndirango, D. Newns and G.J. Martyna
501 – //Symplectic quaternion Scheme for biophysical molecular dynamics
502 – //116(20), 8649, J. Chem. Phys. (2002)
503 – template<typename T> class SQSIntegrator : public T{
504 –  public:
505 –    virtual void moveA();
506 –    virtual void moveB();
507 –  protected:
508 –    void freeRotor();
509 –    void rotate(int k, double dt);
140      
141 < };
142 < #endif
141 > }
142 > #endif //INTEGRATORS_INTEGRATOR_HPP

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines