| 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" | 
| 13 | 
  | 
#include "Thermo.hpp" | 
| 14 | 
  | 
#include "ReadWrite.hpp" | 
| 15 | 
  | 
#include "ZConsWriter.hpp" | 
| 16 | 
+ | 
#include "Restraints.hpp" | 
| 17 | 
  | 
 | 
| 18 | 
  | 
using namespace std; | 
| 19 | 
  | 
const double kB = 8.31451e-7;// boltzmann constant amu*Ang^2*fs^-2/K | 
| 22 | 
  | 
const int maxIteration = 300; | 
| 23 | 
  | 
const double tol = 1.0e-6; | 
| 24 | 
  | 
 | 
| 25 | 
< | 
 | 
| 25 | 
> | 
class VelVerletConsFramework; | 
| 26 | 
  | 
template<typename T = BaseIntegrator> class Integrator : public T { | 
| 27 | 
  | 
 | 
| 28 | 
  | 
public: | 
| 35 | 
  | 
protected: | 
| 36 | 
  | 
 | 
| 37 | 
  | 
  virtual void integrateStep( int calcPot, int calcStress ); | 
| 38 | 
< | 
  virtual void preMove( void ); | 
| 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 ); | 
| 41 | 
> | 
  //virtual void constrainA( void ); | 
| 42 | 
> | 
  //virtual void constrainB( void ); | 
| 43 | 
  | 
  virtual int  readyCheck( void ) { return 1; } | 
| 44 | 
  | 
 | 
| 45 | 
  | 
  virtual void resetIntegrator( void ) { } | 
| 47 | 
  | 
  virtual void calcForce( int calcPot, int calcStress ); | 
| 48 | 
  | 
  virtual void thermalize(); | 
| 49 | 
  | 
 | 
| 50 | 
< | 
  virtual void rotationPropagation( DirectionalAtom* dAtom, double ji[3] ); | 
| 50 | 
> | 
  virtual bool stopIntegrator() {return false;} | 
| 51 | 
  | 
 | 
| 52 | 
< | 
  void checkConstraints( void ); | 
| 52 | 
> | 
  virtual void rotationPropagation( StuntDouble* sd, double ji[3] ); | 
| 53 | 
> | 
 | 
| 54 | 
> | 
  //void checkConstraints( void ); | 
| 55 | 
  | 
  void rotate( int axes1, int axes2, double angle, double j[3], | 
| 56 | 
  | 
         double A[3][3] ); | 
| 57 | 
  | 
 | 
| 58 | 
  | 
  ForceFields* myFF; | 
| 59 | 
  | 
 | 
| 60 | 
  | 
  SimInfo *info; // all the info we'll ever need | 
| 61 | 
+ | 
  vector<StuntDouble*> integrableObjects; | 
| 62 | 
  | 
  int nAtoms;  /* the number of atoms */ | 
| 63 | 
  | 
  int oldAtoms; | 
| 64 | 
  | 
  Atom **atoms; /* array of atom pointers */ | 
| 65 | 
  | 
  Molecule* molecules; | 
| 66 | 
  | 
  int nMols; | 
| 67 | 
  | 
 | 
| 68 | 
< | 
  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 | 
| 68 | 
> | 
  VelVerletConsFramework* consFramework; | 
| 69 | 
  | 
 | 
| 70 | 
< | 
  int* moving; // tells whether we are moving atom i | 
| 71 | 
< | 
  int* moved;  // tells whether we have moved atom i | 
| 72 | 
< | 
  double* oldPos; // pre constrained positions | 
| 70 | 
> | 
  //int isConstrained; // boolean to know whether the systems contains constraints. | 
| 71 | 
> | 
  //int nConstrained;  // counter for number of constraints | 
| 72 | 
> | 
  //int *constrainedA; // the i of a constraint pair | 
| 73 | 
> | 
  //int *constrainedB; // the j of a constraint pair | 
| 74 | 
> | 
  //double *constrainedDsqr; // the square of the constraint distance | 
| 75 | 
> | 
 | 
| 76 | 
> | 
  //int* moving; // tells whether we are moving atom i | 
| 77 | 
> | 
  //int* moved;  // tells whether we have moved atom i | 
| 78 | 
> | 
  //double* oldPos; // pre constrained positions | 
| 79 | 
  | 
 | 
| 80 | 
  | 
  short isFirst; /*boolean for the first time integrate is called */ | 
| 81 | 
  | 
 | 
| 90 | 
  | 
 | 
| 91 | 
  | 
typedef Integrator<BaseIntegrator> RealIntegrator; | 
| 92 | 
  | 
 | 
| 93 | 
+ | 
// ansi instantiation | 
| 94 | 
+ | 
// template class Integrator<BaseIntegrator>; | 
| 95 | 
+ | 
 | 
| 96 | 
  | 
template<typename T> class NVE : public T { | 
| 97 | 
  | 
 | 
| 98 | 
  | 
public: | 
| 427 | 
  | 
 | 
| 428 | 
  | 
  void zeroOutVel(); | 
| 429 | 
  | 
  void doZconstraintForce(); | 
| 430 | 
< | 
  void doHarmonic(); | 
| 430 | 
> | 
  void doHarmonic(vector<double>& resPos); | 
| 431 | 
  | 
  bool checkZConsState(); | 
| 432 | 
  | 
 | 
| 433 | 
  | 
  bool haveFixedZMols(); | 
| 457 | 
  | 
  vector<ZConsParaItem>* parameters; // | 
| 458 | 
  | 
 | 
| 459 | 
  | 
  vector<int> indexOfAllZConsMols;     //index of All Z-Constraint Molecuels | 
| 451 | 
– | 
 | 
| 452 | 
– | 
  int* indexOfZConsMols;                   //index of local Z-Constraint Molecules | 
| 453 | 
– | 
  double* fz; | 
| 454 | 
– | 
  double* curZPos; | 
| 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: | 
| 484 | 
  | 
  double calcMovingMolsCOMVel(); | 
| 485 | 
  | 
  double calcSysCOMVel(); | 
| 486 | 
  | 
  double calcTotalForce(); | 
| 487 | 
< | 
 | 
| 487 | 
> | 
  void updateZPos(); | 
| 488 | 
> | 
  void updateCantPos(); | 
| 489 | 
> | 
   | 
| 490 | 
  | 
  ForceSubtractionPolicy* forcePolicy; //force subtraction policy | 
| 491 | 
  | 
  friend class ForceSubtractionPolicy; | 
| 492 | 
  | 
 | 
| 493 | 
  | 
}; | 
| 494 | 
  | 
 | 
| 476 | 
– | 
class OOPSEMinimizerBase : public RealIntegrator { | 
| 477 | 
– | 
  public: | 
| 495 | 
  | 
 | 
| 496 | 
< | 
    OOPSEMinimizerBase ( SimInfo *theInfo, ForceFields* the_ff ); | 
| 497 | 
< | 
    virtual ~OOPSEMinimizerBase(); | 
| 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); | 
| 508 | 
  | 
     | 
| 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 | 
– | 
     | 
| 509 | 
  | 
}; | 
| 488 | 
– | 
 | 
| 489 | 
– | 
template<typename TMinimizer> class OOPSEMinimizer : public OOPSEMinimizerBase, TMinimizer{ | 
| 490 | 
– | 
  public: | 
| 491 | 
– | 
    void writeOutput(); | 
| 492 | 
– | 
}; | 
| 493 | 
– | 
   | 
| 510 | 
  | 
#endif |