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

Comparing trunk/OOPSE/libmdtools/Integrator.hpp (file contents):
Revision 1010 by tim, Tue Feb 3 20:43:08 2004 UTC vs.
Revision 1452 by tim, Mon Aug 23 15:11:36 2004 UTC

# Line 4 | Line 4
4   #include <string>
5   #include <vector>
6   #include "Atom.hpp"
7 + #include "StuntDouble.hpp"
8   #include "Molecule.hpp"
9   #include "SRI.hpp"
10   #include "AbstractClasses.hpp"
# Line 12 | Line 13
13   #include "Thermo.hpp"
14   #include "ReadWrite.hpp"
15   #include "ZConsWriter.hpp"
16 + #include "Restraints.hpp"
17 + #include "Quaternion.hpp"
18 + #include "Mat4x4d.hpp"
19  
20   using namespace std;
21   const double kB = 8.31451e-7;// boltzmann constant amu*Ang^2*fs^-2/K
# Line 20 | Line 24 | const double tol = 1.0e-6;
24   const int maxIteration = 300;
25   const double tol = 1.0e-6;
26  
27 <
27 > class VelVerletConsFramework;
28   template<typename T = BaseIntegrator> class Integrator : public T {
29  
30   public:
# Line 33 | Line 37 | template<typename T = BaseIntegrator> class Integrator
37   protected:
38  
39    virtual void integrateStep( int calcPot, int calcStress );
40 <  virtual void preMove( void );
40 >  //virtual void preMove( void );
41    virtual void moveA( void );
42    virtual void moveB( void );
43 <  virtual void constrainA( void );
44 <  virtual void constrainB( void );
43 >  //virtual void constrainA( void );
44 >  //virtual void constrainB( void );
45    virtual int  readyCheck( void ) { return 1; }
46  
47    virtual void resetIntegrator( void ) { }
# Line 45 | Line 49 | template<typename T = BaseIntegrator> class Integrator
49    virtual void calcForce( int calcPot, int calcStress );
50    virtual void thermalize();
51  
52 <  virtual void rotationPropagation( DirectionalAtom* dAtom, double ji[3] );
52 >  virtual bool stopIntegrator() {return false;}
53  
54 <  void checkConstraints( void );
54 >  virtual void rotationPropagation( StuntDouble* sd, double ji[3] );
55 >
56 >  //void checkConstraints( void );
57    void rotate( int axes1, int axes2, double angle, double j[3],
58           double A[3][3] );
59  
60 +  void printQuaternion(StuntDouble* sd);
61 +  Mat4x4d getS(const Quaternion& q);
62 +
63    ForceFields* myFF;
64  
65    SimInfo *info; // all the info we'll ever need
66 +  vector<StuntDouble*> integrableObjects;
67    int nAtoms;  /* the number of atoms */
68    int oldAtoms;
69    Atom **atoms; /* array of atom pointers */
70    Molecule* molecules;
71    int nMols;
72  
73 <  int isConstrained; // boolean to know whether the systems contains
64 <         // constraints.
65 <  int nConstrained;  // counter for number of constraints
66 <  int *constrainedA; // the i of a constraint pair
67 <  int *constrainedB; // the j of a constraint pair
68 <  double *constrainedDsqr; // the square of the constraint distance
73 >  VelVerletConsFramework* consFramework;
74  
75 <  int* moving; // tells whether we are moving atom i
76 <  int* moved;  // tells whether we have moved atom i
77 <  double* oldPos; // pre constrained positions
75 >  //int isConstrained; // boolean to know whether the systems contains constraints.
76 >  //int nConstrained;  // counter for number of constraints
77 >  //int *constrainedA; // the i of a constraint pair
78 >  //int *constrainedB; // the j of a constraint pair
79 >  //double *constrainedDsqr; // the square of the constraint distance
80 >
81 >  //int* moving; // tells whether we are moving atom i
82 >  //int* moved;  // tells whether we have moved atom i
83 >  //double* oldPos; // pre constrained positions
84  
85    short isFirst; /*boolean for the first time integrate is called */
86  
# Line 84 | Line 95 | typedef Integrator<BaseIntegrator> RealIntegrator;
95  
96   typedef Integrator<BaseIntegrator> RealIntegrator;
97  
98 + // ansi instantiation
99 + // template class Integrator<BaseIntegrator>;
100 +
101 +
102   template<typename T> class NVE : public T {
103  
104   public:
# Line 418 | Line 433 | template<typename T> class ZConstraint : public T { (p
433  
434    void zeroOutVel();
435    void doZconstraintForce();
436 <  void doHarmonic();
436 >  void doHarmonic(vector<double>& resPos);
437    bool checkZConsState();
438  
439    bool haveFixedZMols();
# Line 449 | Line 464 | template<typename T> class ZConstraint : public T { (p
464  
465    vector<int> indexOfAllZConsMols;     //index of All Z-Constraint Molecuels
466  
467 <  int* indexOfZConsMols;                   //index of local Z-Constraint Molecules
468 <  double* fz;
469 <  double* curZPos;
455 <
467 >  vector<int> indexOfZConsMols;                   //index of local Z-Constraint Molecules
468 >  vector<double> fz;
469 >  vector<double> curZPos;
470  
471 +  bool usingSMD;
472 +  vector<double> prevCantPos;
473 +  vector<double> cantPos;
474 +  vector<double> cantVel;
475  
476 +  double zconsFixTime;  
477 +  double zconsGap;
478 +  bool hasZConsGap;
479 +  vector<double> endFixTime;
480 +  
481    int whichDirection;                           //constraint direction
482  
483   private:
# Line 467 | Line 490 | template<typename T> class ZConstraint : public T { (p
490    double calcMovingMolsCOMVel();
491    double calcSysCOMVel();
492    double calcTotalForce();
493 <
493 >  void updateZPos();
494 >  void updateCantPos();
495 >  
496    ForceSubtractionPolicy* forcePolicy; //force subtraction policy
497    friend class ForceSubtractionPolicy;
498  
499   };
500  
501 < class OOPSEMinimizerBase : public RealIntegrator {
501 >
502 > //Sympletic quaternion Scheme Integrator
503 > //Reference:
504 > // T.F. Miller, M. Eleftheriou, P. Pattnaik, A. Ndirango, D. Newns and G.J. Martyna
505 > //Symplectic quaternion Scheme for biophysical molecular dynamics
506 > //116(20), 8649, J. Chem. Phys. (2002)
507 > template<typename T> class SQSIntegrator : public T{
508    public:
509 +    SQSIntegrator( SimInfo *theInfo, ForceFields* the_ff );
510 +    virtual void moveA();
511 +    virtual void moveB();
512 +  protected:
513 +    void freeRotor(double dt, Quaternion& q, Vector4d& qdot, Vector3d& I);
514 +    void rotate(int k, double dt,  Quaternion& q,  Vector4d& qdot, double Ik);
515 +  private:
516 +    Quaternion getPk(int i, const Vector4d& q);
517 +    Mat4x4d getS(const Quaternion& q);
518 +    vector<Quaternion> q;
519 +    vector<Vector4d> p_qua;
520 +    vector<double> I0;
521  
522 <    OOPSEMinimizerBase ( SimInfo *theInfo, ForceFields* the_ff );
523 <    virtual ~OOPSEMinimizerBase();
481 <    
482 <    double calcGradient(vector<double>& x, vector<double>& grad);
483 <    void setOptCoor(vector<double>& x);
484 <    void getOptCoor(vector<double>& x);
485 <    void output();
486 <    
522 >                Vector4d p2j(const Vector4d& p, Mat4x4d& m);    
523 >                Vector4d j2p(const Vector4d& j, Mat4x4d& m);
524   };
488 –
489 – template<typename TMinimizer> class OOPSEMinimizer : public OOPSEMinimizerBase, TMinimizer{
490 –  public:
491 –    void writeOutput();
492 – };
493 –  
525   #endif

Diff Legend

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