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

Comparing trunk/OOPSE-4/src/integrators/Integrator.hpp (file contents):
Revision 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 2008 by tim, Sun Feb 13 19:10:25 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 >        StatWriter* statWriter;
117 >        Thermo thermo;
118  
119 +        double runTime;
120 +        double sampleTime;
121 +        double statusTime;
122 +        double thermalTime;
123 +        double dt;
124  
125 <  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
125 >        Snapshot* currentSnapshot_; //During the integration, the address of currentSnapshot Will not change
126  
127 <  int* moving; // tells whether we are moving atom i
128 <  int* moved;  // tells whether we have moved atom i
129 <  double* oldPos; // pre constrained positions
127 >        
128 >    private:
129 >        
130 >        virtual double calcConservedQuantity() = 0;
131 >        
132 >        virtual DumpWriter* createDumpWriter() = 0;
133  
134 <  short isFirst; /*boolean for the first time integrate is called */
80 <
81 <  double dt;
82 <  double dt2;
83 <
84 <  Thermo *tStats;
85 <  StatWriter*  statOut;
86 <  DumpWriter*  dumpOut;
87 <
88 < };
89 <
90 < typedef Integrator<BaseIntegrator> RealIntegrator;
91 <
92 < // ansi instantiation
93 < // template class Integrator<BaseIntegrator>;
94 <
95 <
96 < template<typename T> class NVE : public T {
97 <
98 < public:
99 <  NVE ( SimInfo *theInfo, ForceFields* the_ff ):
100 <    T( theInfo, the_ff ){}
101 <  virtual ~NVE(){}
102 < };
103 <
104 <
105 < template<typename T> class NVT : public T {
106 <
107 < public:
108 <
109 <  NVT ( SimInfo *theInfo, ForceFields* the_ff);
110 <  virtual ~NVT();
111 <
112 <  void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
113 <  void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
114 <  void setChiTolerance(double tol) {chiTolerance = tol;}
115 <  virtual double  getConservedQuantity(void);
116 <  virtual string getAdditionalParameters(void);
117 <
118 < protected:
119 <
120 <  virtual void moveA( void );
121 <  virtual void moveB( void );
122 <
123 <  virtual int readyCheck();
124 <
125 <  virtual void resetIntegrator( void );
126 <
127 <  // chi is a propagated degree of freedom.
128 <
129 <  double chi;
130 <
131 <  //integral of chi(t)dt
132 <  double integralOfChidt;
133 <
134 <  // targetTemp must be set.  tauThermostat must also be set;
135 <
136 <  double targetTemp;
137 <  double tauThermostat;
138 <
139 <  short int have_tau_thermostat, have_target_temp;
140 <
141 <  double *oldVel;
142 <  double *oldJi;
143 <
144 <  double chiTolerance;
145 <  short int have_chi_tolerance;
146 <
147 < };
148 <
149 <
150 <
151 < template<typename T> class NPT : public T{
152 <
153 < public:
154 <
155 <  NPT ( SimInfo *theInfo, ForceFields* the_ff);
156 <  virtual ~NPT();
157 <
158 <  virtual void integrateStep( int calcPot, int calcStress ){
159 <    calcStress = 1;
160 <    T::integrateStep( calcPot, calcStress );
161 <  }
162 <
163 <  virtual double getConservedQuantity(void) = 0;
164 <  virtual string getAdditionalParameters(void) = 0;
165 <  
166 <  double myTauThermo( void ) { return tauThermostat; }
167 <  double myTauBaro( void ) { return tauBarostat; }
168 <
169 <  void setTauThermostat(double tt) {tauThermostat = tt; have_tau_thermostat=1;}
170 <  void setTauBarostat(double tb) {tauBarostat = tb; have_tau_barostat=1;}
171 <  void setTargetTemp(double tt) {targetTemp = tt; have_target_temp = 1;}
172 <  void setTargetPressure(double tp) {targetPressure = tp; have_target_pressure = 1;}
173 <  void setChiTolerance(double tol) {chiTolerance = tol; have_chi_tolerance = 1;}
174 <  void setPosIterTolerance(double tol) {posIterTolerance = tol; have_pos_iter_tolerance = 1;}
175 <  void setEtaTolerance(double tol) {etaTolerance = tol; have_eta_tolerance = 1;}
176 <
177 < protected:
178 <
179 <  virtual void  moveA( void );
180 <  virtual void moveB( void );
181 <
182 <  virtual int readyCheck();
183 <
184 <  virtual void resetIntegrator( void );
185 <
186 <  virtual void getVelScaleA( double sc[3], double vel[3] ) = 0;
187 <  virtual void getVelScaleB( double sc[3], int index ) = 0;
188 <  virtual void getPosScale(double pos[3], double COM[3],
189 <                           int index, double sc[3]) = 0;
190 <
191 <  virtual void calcVelScale( void ) = 0;
192 <
193 <  virtual bool chiConverged( void );
194 <  virtual bool etaConverged( void ) = 0;
195 <
196 <  virtual void evolveChiA( void );
197 <  virtual void evolveEtaA( void ) = 0;
198 <  virtual void evolveChiB( void );
199 <  virtual void evolveEtaB( void ) = 0;
200 <
201 <  virtual void scaleSimBox( void ) = 0;
202 <
203 <  void accIntegralOfChidt(void) { integralOfChidt += dt * chi;}
204 <
205 <  // chi and eta are the propagated degrees of freedom
206 <
207 <  double oldChi;
208 <  double prevChi;
209 <  double chi;
210 <  double NkBT;
211 <  double fkBT;
212 <
213 <  double tt2, tb2;
214 <  double instaTemp, instaPress, instaVol;
215 <  double press[3][3];
216 <
217 <  int Nparticles;
218 <
219 <  double integralOfChidt;
220 <
221 <  // targetTemp, targetPressure, and tauBarostat must be set.
222 <  // One of qmass or tauThermostat must be set;
223 <
224 <  double targetTemp;
225 <  double targetPressure;
226 <  double tauThermostat;
227 <  double tauBarostat;
228 <
229 <  short int have_tau_thermostat, have_tau_barostat, have_target_temp;
230 <  short int have_target_pressure;
231 <
232 <  double *oldPos;
233 <  double *oldVel;
234 <  double *oldJi;
235 <
236 <  double chiTolerance;
237 <  short int have_chi_tolerance;
238 <  double posIterTolerance;
239 <  short int have_pos_iter_tolerance;
240 <  double etaTolerance;
241 <  short int have_eta_tolerance;
242 <
134 >        virtual StatWriter* createStatWriter() = 0;
135   };
136  
245 – template<typename T> class NPTi : public T{
246 –
247 – public:
248 –  NPTi( SimInfo *theInfo, ForceFields* the_ff);
249 –  ~NPTi();
250 –
251 –  virtual double getConservedQuantity(void);
252 –  virtual void resetIntegrator(void);
253 –  virtual string getAdditionalParameters(void);
254 – protected:
255 –
256 –
257 –
258 –  virtual void evolveEtaA(void);
259 –  virtual void evolveEtaB(void);
260 –
261 –  virtual bool etaConverged( void );
262 –
263 –  virtual void scaleSimBox( void );
264 –
265 –  virtual void getVelScaleA( double sc[3], double vel[3] );
266 –  virtual void getVelScaleB( double sc[3], int index );
267 –  virtual void getPosScale(double pos[3], double COM[3],
268 –                           int index, double sc[3]);
269 –
270 –  virtual void calcVelScale( void );
271 –
272 –  double eta, oldEta, prevEta;
273 –  double vScale;
274 – };
275 –
276 – template<typename T> class NPTf : public T{
277 –
278 – public:
279 –
280 –  NPTf ( SimInfo *theInfo, ForceFields* the_ff);
281 –  virtual ~NPTf();
282 –
283 –  virtual double getConservedQuantity(void);
284 –  virtual string getAdditionalParameters(void);
285 –  virtual void resetIntegrator(void);
286 –
287 – protected:
288 –
289 –  virtual void evolveEtaA(void);
290 –  virtual void evolveEtaB(void);
291 –
292 –  virtual bool etaConverged( void );
293 –
294 –  virtual void scaleSimBox( void );
295 –
296 –  virtual void getVelScaleA( double sc[3], double vel[3] );
297 –  virtual void getVelScaleB( double sc[3], int index );
298 –  virtual void getPosScale(double pos[3], double COM[3],
299 –                           int index, double sc[3]);
300 –
301 –  virtual void calcVelScale( void );
302 –
303 –  double eta[3][3];
304 –  double oldEta[3][3];
305 –  double prevEta[3][3];
306 –  double vScale[3][3];
307 – };
308 –
309 – template<typename T> class NPTxyz : public T{
310 –
311 – public:
312 –
313 –  NPTxyz ( SimInfo *theInfo, ForceFields* the_ff);
314 –  virtual ~NPTxyz();
315 –
316 –  virtual double getConservedQuantity(void);
317 –  virtual string getAdditionalParameters(void);
318 –  virtual void resetIntegrator(void);
319 –
320 – protected:
321 –
322 –  virtual void evolveEtaA(void);
323 –  virtual void evolveEtaB(void);
324 –
325 –  virtual bool etaConverged( void );
326 –
327 –  virtual void scaleSimBox( void );
328 –
329 –  virtual void getVelScaleA( double sc[3], double vel[3] );
330 –  virtual void getVelScaleB( double sc[3], int index );
331 –  virtual void getPosScale(double pos[3], double COM[3],
332 –                           int index, double sc[3]);
333 –
334 –  virtual void calcVelScale( void );
335 –
336 –  double eta[3][3];
337 –  double oldEta[3][3];
338 –  double prevEta[3][3];
339 –  double vScale[3][3];
340 – };
341 –
342 –
343 – template<typename T> class ZConstraint : public T {
344 –
345 –  public:
346 –  class ForceSubtractionPolicy{
347 –    public:
348 –      ForceSubtractionPolicy(ZConstraint<T>* integrator) {zconsIntegrator = integrator;}
349 –
350 –      virtual void update() = 0;
351 –      virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0;
352 –      virtual double getZFOfMovingMols(Atom* atom, double totalForce) = 0;
353 –      virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) = 0;
354 –      virtual double getHFOfUnconsMols(Atom* atom, double totalForce) = 0;
355 –
356 –   protected:
357 –     ZConstraint<T>* zconsIntegrator;
358 –  };
359 –
360 –  class PolicyByNumber : public ForceSubtractionPolicy{
361 –
362 –    public:
363 –      PolicyByNumber(ZConstraint<T>* integrator) :ForceSubtractionPolicy(integrator) {}
364 –      virtual void update();
365 –      virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ;
366 –      virtual double getZFOfMovingMols(Atom* atom, double totalForce) ;
367 –      virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce);
368 –      virtual double getHFOfUnconsMols(Atom* atom, double totalForce);
369 –
370 –    private:
371 –      int totNumOfMovingAtoms;
372 –  };
373 –
374 –  class PolicyByMass : public ForceSubtractionPolicy{
375 –
376 –    public:
377 –      PolicyByMass(ZConstraint<T>* integrator) :ForceSubtractionPolicy(integrator) {}
378 –
379 –      virtual void update();
380 –      virtual double getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce) ;
381 –      virtual double getZFOfMovingMols(Atom* atom, double totalForce) ;
382 –      virtual double getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce);
383 –      virtual double getHFOfUnconsMols(Atom* atom, double totalForce);
384 –
385 –   private:
386 –     double totMassOfMovingAtoms;
387 –  };
388 –
389 – public:
390 –
391 –  ZConstraint( SimInfo *theInfo, ForceFields* the_ff);
392 –  ~ZConstraint();
393 –
394 –  void setZConsTime(double time)                  {this->zconsTime = time;}
395 –  void getZConsTime()                             {return zconsTime;}
396 –
397 –  void setIndexOfAllZConsMols(vector<int> index) {indexOfAllZConsMols = index;}
398 –  void getIndexOfAllZConsMols()                  {return indexOfAllZConsMols;}
399 –
400 –  void setZConsOutput(const char * fileName)          {zconsOutput = fileName;}
401 –  string getZConsOutput()                         {return zconsOutput;}
402 –
403 –  virtual void integrate();
404 –
405 –
406 – #ifdef IS_MPI
407 –  virtual void update();                      //which is called to indicate the molecules' migration
408 – #endif
409 –
410 –  enum ZConsState {zcsMoving, zcsFixed};
411 –
412 –  vector<Molecule*> zconsMols;              //z-constraint molecules array
413 –  vector<ZConsState> states;                 //state of z-constraint molecules
414 –
415 –
416 –
417 –  int totNumOfUnconsAtoms;              //total number of uncontraint atoms
418 –  double totalMassOfUncons;                //total mas of unconstraint molecules
419 –
420 –
421 – protected:
422 –
423 –
424 –
425 –  virtual void calcForce( int calcPot, int calcStress );
426 –  virtual void thermalize(void);
427 –
428 –  void zeroOutVel();
429 –  void doZconstraintForce();
430 –  void doHarmonic(vector<double>& resPos);
431 –  bool checkZConsState();
432 –
433 –  bool haveFixedZMols();
434 –  bool haveMovingZMols();
435 –
436 –  double calcZSys();
437 –
438 –  int isZConstraintMol(Molecule* mol);
439 –
440 –
441 –  double zconsTime;                              //sample time
442 –  double zconsTol;                                 //tolerance of z-contratint
443 –  double zForceConst;                           //base force constant term
444 –                                                          //which is estimate by OOPSE
445 –
446 –
447 –  vector<double> massOfZConsMols;       //mass of z-constraint molecule
448 –  vector<double> kz;                              //force constant array
449 –
450 –  vector<double> zPos;                          //
451 –
452 –
453 –  vector<Molecule*> unconsMols;           //unconstraint molecules array
454 –  vector<double> massOfUnconsMols;    //mass array of unconstraint molecules
455 –
456 –
457 –  vector<ZConsParaItem>* parameters; //
458 –
459 –  vector<int> indexOfAllZConsMols;     //index of All Z-Constraint Molecuels
460 –
461 –  vector<int> indexOfZConsMols;                   //index of local Z-Constraint Molecules
462 –  vector<double> fz;
463 –  vector<double> curZPos;
464 –
465 –  bool usingSMD;
466 –  vector<double> prevCantPos;
467 –  vector<double> cantPos;
468 –  vector<double> cantVel;
469 –
470 –  double zconsFixTime;  
471 –  double zconsGap;
472 –  bool hasZConsGap;
473 –  vector<double> endFixTime;
474 –  
475 –  int whichDirection;                           //constraint direction
476 –
477 – private:
478 –
479 –  string zconsOutput;                         //filename of zconstraint output
480 –  ZConsWriter* fzOut;                         //z-constraint writer
481 –
482 –  double curZconsTime;
483 –
484 –  double calcMovingMolsCOMVel();
485 –  double calcSysCOMVel();
486 –  double calcTotalForce();
487 –  void updateZPos();
488 –  void updateCantPos();
489 –  
490 –  ForceSubtractionPolicy* forcePolicy; //force subtraction policy
491 –  friend class ForceSubtractionPolicy;
492 –
493 – };
494 –
495 –
496 – //Sympletic quaternion Scheme Integrator
497 – //Reference:
498 – // T.F. Miller, M. Eleftheriou, P. Pattnaik, A. Ndirango, D. Newns and G.J. Martyna
499 – //Symplectic quaternion Scheme for biophysical molecular dynamics
500 – //116(20), 8649, J. Chem. Phys. (2002)
501 – template<typename T> class SQSIntegrator : public T{
502 –  public:
503 –    virtual void moveA();
504 –    virtual void moveB();
505 –  protected:
506 –    void freeRotor();
507 –    void rotate(int k, double dt);
137      
138 < };
139 < #endif
138 > }
139 > #endif //INTEGRATORS_INTEGRATOR_HPP

Diff Legend

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