OOPSE/libmdtools/Integrator.cpp

#include <iostream>
#include <cstdlib>
#include <cmath>

#ifdef IS_MPI
#include "mpiSimulation.hpp"
#include <unistd.h>
#endif //is_mpi

#include "Integrator.hpp"
#include "simError.h"


template<typename T> Integrator<T>::Integrator(SimInfo* theInfo,
                                               ForceFields* the_ff){
  info = theInfo;
  myFF = the_ff;
  isFirst = 1;

  molecules = info->molecules;
  nMols = info->n_mol;

  // give a little love back to the SimInfo object

  if (info->the_integrator != NULL){
    delete info->the_integrator;
  }
  info->the_integrator = this;

  nAtoms = info->n_atoms;

  // check for constraints

  constrainedA = NULL;
  constrainedB = NULL;
  constrainedDsqr = NULL;
  moving = NULL;
  moved = NULL;
  oldPos = NULL;

  nConstrained = 0;

  checkConstraints();
}

template<typename T> Integrator<T>::~Integrator(){
  if (nConstrained){
    delete[] constrainedA;
    delete[] constrainedB;
    delete[] constrainedDsqr;
    delete[] moving;
    delete[] moved;
    delete[] oldPos;
  }
}

template<typename T> void Integrator<T>::checkConstraints(void){
  isConstrained = 0;

  Constraint* temp_con;
  Constraint* dummy_plug;
  temp_con = new Constraint[info->n_SRI];
  nConstrained = 0;
  int constrained = 0;

  SRI** theArray;
  for (int i = 0; i < nMols; i++){
    theArray = (SRI * *) molecules[i].getMyBonds();
    for (int j = 0; j < molecules[i].getNBonds(); j++){
      constrained = theArray[j]->is_constrained();

      if (constrained){
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a(dummy_plug->get_a());
        temp_con[nConstrained].set_b(dummy_plug->get_b());
        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());

        nConstrained++;
        constrained = 0;
      }
    }

    theArray = (SRI * *) molecules[i].getMyBends();
    for (int j = 0; j < molecules[i].getNBends(); j++){
      constrained = theArray[j]->is_constrained();

      if (constrained){
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a(dummy_plug->get_a());
        temp_con[nConstrained].set_b(dummy_plug->get_b());
        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());

        nConstrained++;
        constrained = 0;
      }
    }

    theArray = (SRI * *) molecules[i].getMyTorsions();
    for (int j = 0; j < molecules[i].getNTorsions(); j++){
      constrained = theArray[j]->is_constrained();

      if (constrained){
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a(dummy_plug->get_a());
        temp_con[nConstrained].set_b(dummy_plug->get_b());
        temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());

        nConstrained++;
        constrained = 0;
      }
    }
  }

  if (nConstrained > 0){
    isConstrained = 1;

    if (constrainedA != NULL)
      delete[] constrainedA;
    if (constrainedB != NULL)
      delete[] constrainedB;
    if (constrainedDsqr != NULL)
      delete[] constrainedDsqr;

    constrainedA = new int[nConstrained];
    constrainedB = new int[nConstrained];
    constrainedDsqr = new double[nConstrained];

    for (int i = 0; i < nConstrained; i++){
      constrainedA[i] = temp_con[i].get_a();
      constrainedB[i] = temp_con[i].get_b();
      constrainedDsqr[i] = temp_con[i].get_dsqr();
    }


    // save oldAtoms to check for lode balanceing later on.

    oldAtoms = nAtoms;

    moving = new int[nAtoms];
    moved = new int[nAtoms];

    oldPos = new double[nAtoms * 3];
  }

  delete[] temp_con;
}


template<typename T> void Integrator<T>::integrate(void){
  int i, j;                         // loop counters

  double runTime = info->run_time;
  double sampleTime = info->sampleTime;
  double statusTime = info->statusTime;
  double thermalTime = info->thermalTime;
  double resetTime = info->resetTime;


  double currSample;
  double currThermal;
  double currStatus;
  double currReset;
  
  int calcPot, calcStress;
  int isError;

  tStats = new Thermo(info);
  statOut = new StatWriter(info);
  dumpOut = new DumpWriter(info);

  atoms = info->atoms;
  DirectionalAtom* dAtom;

  dt = info->dt;
  dt2 = 0.5 * dt;

  // initialize the forces before the first step

  calcForce(1, 1);
  
  if (info->setTemp){
    thermalize();
  }

  calcPot     = 0;
  calcStress  = 0;
  currSample  = sampleTime + info->getTime();
  currThermal = thermalTime+ info->getTime();
  currStatus  = statusTime + info->getTime();
  currReset   = resetTime  + info->getTime();

  dumpOut->writeDump(info->getTime());
  statOut->writeStat(info->getTime());

  readyCheck();

#ifdef IS_MPI
  strcpy(checkPointMsg, "The integrator is ready to go.");
  MPIcheckPoint();
#endif // is_mpi

  while (info->getTime() < runTime){
    if ((info->getTime() + dt) >= currStatus){
      calcPot = 1;
      calcStress = 1;
    }

    integrateStep(calcPot, calcStress);

    info->incrTime(dt);

    if (info->setTemp){
      if (info->getTime() >= currThermal){
        thermalize();
        currThermal += thermalTime;
      }
    }

    if (info->getTime() >= currSample){
      dumpOut->writeDump(info->getTime());
      currSample += sampleTime;
    }

    if (info->getTime() >= currStatus){
      statOut->writeStat(info->getTime()); 
      calcPot = 0; 
      calcStress = 0;
      currStatus += statusTime;
    } 

    if (info->resetIntegrator){
      if (info->getTime() >= currReset){
        this->resetIntegrator();
        currReset += resetTime;
      }
    }

#ifdef IS_MPI
    strcpy(checkPointMsg, "successfully took a time step.");
    MPIcheckPoint();
#endif // is_mpi
  }

  dumpOut->writeFinal(info->getTime());

  delete dumpOut;
  delete statOut;
}

template<typename T> void Integrator<T>::integrateStep(int calcPot,
                                                       int calcStress){
  // Position full step, and velocity half step
  preMove();

  moveA();


#ifdef IS_MPI
  strcpy(checkPointMsg, "Succesful moveA\n");
  MPIcheckPoint();
#endif // is_mpi


  // calc forces

  calcForce(calcPot, calcStress);

#ifdef IS_MPI
  strcpy(checkPointMsg, "Succesful doForces\n");
  MPIcheckPoint();
#endif // is_mpi


  // finish the velocity  half step

  moveB();


#ifdef IS_MPI
  strcpy(checkPointMsg, "Succesful moveB\n");
  MPIcheckPoint();
#endif // is_mpi
}


template<typename T> void Integrator<T>::moveA(void){
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double A[3][3], I[3][3];
  double angle;
  double vel[3], pos[3], frc[3];
  double mass;

  for (i = 0; i < nAtoms; i++){
    atoms[i]->getVel(vel);
    atoms[i]->getPos(pos);
    atoms[i]->getFrc(frc);

    mass = atoms[i]->getMass();

    for (j = 0; j < 3; j++){
      // velocity half step
      vel[j] += (dt2 * frc[j] / mass) * eConvert;
      // position whole step
      pos[j] += dt * vel[j];
    }

    atoms[i]->setVel(vel);
    atoms[i]->setPos(pos);

    if (atoms[i]->isDirectional()){
      dAtom = (DirectionalAtom *) atoms[i];

      // get and convert the torque to body frame

      dAtom->getTrq(Tb);
      dAtom->lab2Body(Tb);

      // get the angular momentum, and propagate a half step

      dAtom->getJ(ji);

      for (j = 0; j < 3; j++)
        ji[j] += (dt2 * Tb[j]) * eConvert;

      // use the angular velocities to propagate the rotation matrix a
      // full time step

      dAtom->getA(A);
      dAtom->getI(I);

      // rotate about the x-axis      
      angle = dt2 * ji[0] / I[0][0];
      this->rotate(1, 2, angle, ji, A); 

      // rotate about the y-axis
      angle = dt2 * ji[1] / I[1][1];
      this->rotate(2, 0, angle, ji, A);

      // rotate about the z-axis
      angle = dt * ji[2] / I[2][2];
      this->rotate(0, 1, angle, ji, A);

      // rotate about the y-axis
      angle = dt2 * ji[1] / I[1][1];
      this->rotate(2, 0, angle, ji, A);

      // rotate about the x-axis
      angle = dt2 * ji[0] / I[0][0];
      this->rotate(1, 2, angle, ji, A);

      dAtom->setJ(ji);
      dAtom->setA(A);
    }
  }

  if (nConstrained){
    constrainA();
  }
}


template<typename T> void Integrator<T>::moveB(void){
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double vel[3], frc[3];
  double mass;

  for (i = 0; i < nAtoms; i++){
    atoms[i]->getVel(vel);
    atoms[i]->getFrc(frc);

    mass = atoms[i]->getMass();

    // velocity half step
    for (j = 0; j < 3; j++)
      vel[j] += (dt2 * frc[j] / mass) * eConvert;

    atoms[i]->setVel(vel);

    if (atoms[i]->isDirectional()){
      dAtom = (DirectionalAtom *) atoms[i];

      // get and convert the torque to body frame      

      dAtom->getTrq(Tb);
      dAtom->lab2Body(Tb);

      // get the angular momentum, and propagate a half step

      dAtom->getJ(ji);

      for (j = 0; j < 3; j++)
        ji[j] += (dt2 * Tb[j]) * eConvert;


      dAtom->setJ(ji);
    }
  }

  if (nConstrained){
    constrainB();
  }
}

template<typename T> void Integrator<T>::preMove(void){
  int i, j;
  double pos[3];

  if (nConstrained){
    for (i = 0; i < nAtoms; i++){
      atoms[i]->getPos(pos);

      for (j = 0; j < 3; j++){
        oldPos[3 * i + j] = pos[j];
      }
    }
  }
}

template<typename T> void Integrator<T>::constrainA(){
  int i, j, k;
  int done;
  double posA[3], posB[3];
  double velA[3], velB[3];
  double pab[3];
  double rab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double dx, dy, dz;
  double rpab;
  double rabsq, pabsq, rpabsq;
  double diffsq;
  double gab;
  int iteration;

  for (i = 0; i < nAtoms; i++){
    moving[i] = 0;
    moved[i] = 1;
  }

  iteration = 0;
  done = 0;
  while (!done && (iteration < maxIteration)){
    done = 1;
    for (i = 0; i < nConstrained; i++){
      a = constrainedA[i];
      b = constrainedB[i];

      ax = (a * 3) + 0;
      ay = (a * 3) + 1;
      az = (a * 3) + 2;

      bx = (b * 3) + 0;
      by = (b * 3) + 1;
      bz = (b * 3) + 2;

      if (moved[a] || moved[b]){
        atoms[a]->getPos(posA);
        atoms[b]->getPos(posB);

        for (j = 0; j < 3; j++)
          pab[j] = posA[j] - posB[j];

        //periodic boundary condition

        info->wrapVector(pab);

        pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];

        rabsq = constrainedDsqr[i];
        diffsq = rabsq - pabsq;

        // the original rattle code from alan tidesley
        if (fabs(diffsq) > (tol * rabsq * 2)){
          rab[0] = oldPos[ax] - oldPos[bx];
          rab[1] = oldPos[ay] - oldPos[by];
          rab[2] = oldPos[az] - oldPos[bz];

          info->wrapVector(rab);

          rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];

          rpabsq = rpab * rpab;


          if (rpabsq < (rabsq * -diffsq)){
#ifdef IS_MPI
            a = atoms[a]->getGlobalIndex();
            b = atoms[b]->getGlobalIndex();
#endif //is_mpi
            sprintf(painCave.errMsg,
                    "Constraint failure in constrainA at atom %d and %d.\n", a,
                    b);
            painCave.isFatal = 1;
            simError();
          }

          rma = 1.0 / atoms[a]->getMass();
          rmb = 1.0 / atoms[b]->getMass();

          gab = diffsq / (2.0 * (rma + rmb) * rpab);

          dx = rab[0] * gab;
          dy = rab[1] * gab;
          dz = rab[2] * gab;

          posA[0] += rma * dx;
          posA[1] += rma * dy;
          posA[2] += rma * dz;

          atoms[a]->setPos(posA);

          posB[0] -= rmb * dx;
          posB[1] -= rmb * dy;
          posB[2] -= rmb * dz;

          atoms[b]->setPos(posB);

          dx = dx / dt;
          dy = dy / dt;
          dz = dz / dt;

          atoms[a]->getVel(velA);

          velA[0] += rma * dx;
          velA[1] += rma * dy;
          velA[2] += rma * dz;

          atoms[a]->setVel(velA);

          atoms[b]->getVel(velB);

          velB[0] -= rmb * dx;
          velB[1] -= rmb * dy;
          velB[2] -= rmb * dz;

          atoms[b]->setVel(velB);

          moving[a] = 1;
          moving[b] = 1;
          done = 0;
        }
      }
    }

    for (i = 0; i < nAtoms; i++){
      moved[i] = moving[i];
      moving[i] = 0;
    }

    iteration++;
  }

  if (!done){
    sprintf(painCave.errMsg,
            "Constraint failure in constrainA, too many iterations: %d\n",
            iteration);
    painCave.isFatal = 1;
    simError();
  }

}

template<typename T> void Integrator<T>::constrainB(void){
  int i, j, k;
  int done;
  double posA[3], posB[3];
  double velA[3], velB[3];
  double vxab, vyab, vzab;
  double rab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double dx, dy, dz;
  double rabsq, pabsq, rvab;
  double diffsq;
  double gab;
  int iteration;

  for (i = 0; i < nAtoms; i++){
    moving[i] = 0;
    moved[i] = 1;
  }

  done = 0;
  iteration = 0;
  while (!done && (iteration < maxIteration)){
    done = 1;

    for (i = 0; i < nConstrained; i++){
      a = constrainedA[i];
      b = constrainedB[i];

      ax = (a * 3) + 0;
      ay = (a * 3) + 1;
      az = (a * 3) + 2;

      bx = (b * 3) + 0;
      by = (b * 3) + 1;
      bz = (b * 3) + 2;

      if (moved[a] || moved[b]){
        atoms[a]->getVel(velA);
        atoms[b]->getVel(velB);

        vxab = velA[0] - velB[0];
        vyab = velA[1] - velB[1];
        vzab = velA[2] - velB[2];

        atoms[a]->getPos(posA);
        atoms[b]->getPos(posB);

        for (j = 0; j < 3; j++)
          rab[j] = posA[j] - posB[j];

        info->wrapVector(rab);

        rma = 1.0 / atoms[a]->getMass();
        rmb = 1.0 / atoms[b]->getMass();

        rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;

        gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);

        if (fabs(gab) > tol){
          dx = rab[0] * gab;
          dy = rab[1] * gab;
          dz = rab[2] * gab;

          velA[0] += rma * dx;
          velA[1] += rma * dy;
          velA[2] += rma * dz;

          atoms[a]->setVel(velA);

          velB[0] -= rmb * dx;
          velB[1] -= rmb * dy;
          velB[2] -= rmb * dz;

          atoms[b]->setVel(velB);

          moving[a] = 1;
          moving[b] = 1;
          done = 0;
        }
      }
    }

    for (i = 0; i < nAtoms; i++){
      moved[i] = moving[i];
      moving[i] = 0;
    }

    iteration++;
  }

  if (!done){
    sprintf(painCave.errMsg,
            "Constraint failure in constrainB, too many iterations: %d\n",
            iteration);
    painCave.isFatal = 1;
    simError();
  }
}

template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
                                                double angle, double ji[3],
                                                double A[3][3]){
  int i, j, k;
  double sinAngle;
  double cosAngle;
  double angleSqr;
  double angleSqrOver4;
  double top, bottom;
  double rot[3][3];
  double tempA[3][3];
  double tempJ[3];

  // initialize the tempA

  for (i = 0; i < 3; i++){
    for (j = 0; j < 3; j++){
      tempA[j][i] = A[i][j];
    }
  }

  // initialize the tempJ

  for (i = 0; i < 3; i++)
    tempJ[i] = ji[i];

  // initalize rot as a unit matrix

  rot[0][0] = 1.0;
  rot[0][1] = 0.0;
  rot[0][2] = 0.0;

  rot[1][0] = 0.0;
  rot[1][1] = 1.0;
  rot[1][2] = 0.0;

  rot[2][0] = 0.0;
  rot[2][1] = 0.0;
  rot[2][2] = 1.0;

  // use a small angle aproximation for sin and cosine

  angleSqr = angle * angle;
  angleSqrOver4 = angleSqr / 4.0;
  top = 1.0 - angleSqrOver4;
  bottom = 1.0 + angleSqrOver4;

  cosAngle = top / bottom;
  sinAngle = angle / bottom;

  rot[axes1][axes1] = cosAngle;
  rot[axes2][axes2] = cosAngle;

  rot[axes1][axes2] = sinAngle;
  rot[axes2][axes1] = -sinAngle;

  // rotate the momentum acoording to: ji[] = rot[][] * ji[]

  for (i = 0; i < 3; i++){
    ji[i] = 0.0;
    for (k = 0; k < 3; k++){
      ji[i] += rot[i][k] * tempJ[k];
    }
  }

  // rotate the Rotation matrix acording to: 
  //            A[][] = A[][] * transpose(rot[][])


  // NOte for as yet unknown reason, we are performing the
  // calculation as:
  //                transpose(A[][]) = transpose(A[][]) * transpose(rot[][])

  for (i = 0; i < 3; i++){
    for (j = 0; j < 3; j++){
      A[j][i] = 0.0;
      for (k = 0; k < 3; k++){
        A[j][i] += tempA[i][k] * rot[j][k];
      }
    }
  }
}

template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
  myFF->doForces(calcPot, calcStress);
}

template<typename T> void Integrator<T>::thermalize(){
  tStats->velocitize();
}

template<typename T> double Integrator<T>::getConservedQuantity(void){
  return tStats->getTotalE();
}
Revision:	768
Committed:	Wed Sep 17 14:22:15 2003 UTC (20 years, 9 months ago) by mmeineke
File size:	16481 byte(s)
Log Message:	fixed NPTi to now work with constraints.
#	User	Rev	Content
1	mmeineke	558	#include <iostream>
2			#include <cstdlib>
3	mmeineke	561	#include <cmath>
4	mmeineke	558
5			#ifdef IS_MPI
6			#include "mpiSimulation.hpp"
7			#include <unistd.h>
8			#endif //is_mpi
9
10			#include "Integrator.hpp"
11			#include "simError.h"
12
13
14	tim	725	template<typename T> Integrator<T>::Integrator(SimInfo* theInfo,
15			ForceFields* the_ff){
16	mmeineke	558	info = theInfo;
17			myFF = the_ff;
18			isFirst = 1;
19
20			molecules = info->molecules;
21			nMols = info->n_mol;
22
23			// give a little love back to the SimInfo object
24	tim	725
25			if (info->the_integrator != NULL){
26			delete info->the_integrator;
27			}
28	mmeineke	558	info->the_integrator = this;
29
30			nAtoms = info->n_atoms;
31
32			// check for constraints
33	tim	725
34			constrainedA = NULL;
35			constrainedB = NULL;
36	mmeineke	558	constrainedDsqr = NULL;
37	tim	725	moving = NULL;
38			moved = NULL;
39			oldPos = NULL;
40
41	mmeineke	558	nConstrained = 0;
42
43			checkConstraints();
44			}
45
46	tim	725	template<typename T> Integrator<T>::~Integrator(){
47			if (nConstrained){
48	mmeineke	558	delete[] constrainedA;
49			delete[] constrainedB;
50			delete[] constrainedDsqr;
51			delete[] moving;
52			delete[] moved;
53	mmeineke	561	delete[] oldPos;
54	mmeineke	558	}
55			}
56
57	tim	725	template<typename T> void Integrator<T>::checkConstraints(void){
58	mmeineke	558	isConstrained = 0;
59
60	tim	725	Constraint* temp_con;
61			Constraint* dummy_plug;
62	mmeineke	558	temp_con = new Constraint[info->n_SRI];
63			nConstrained = 0;
64			int constrained = 0;
65	tim	725
66	mmeineke	558	SRI** theArray;
67	tim	725	for (int i = 0; i < nMols; i++){
68			theArray = (SRI * *) molecules[i].getMyBonds();
69			for (int j = 0; j < molecules[i].getNBonds(); j++){
70	mmeineke	558	constrained = theArray[j]->is_constrained();
71	mmeineke	594
72	tim	725	if (constrained){
73			dummy_plug = theArray[j]->get_constraint();
74			temp_con[nConstrained].set_a(dummy_plug->get_a());
75			temp_con[nConstrained].set_b(dummy_plug->get_b());
76			temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
77	mmeineke	594
78	tim	725	nConstrained++;
79			constrained = 0;
80			}
81	mmeineke	558	}
82
83	tim	725	theArray = (SRI * *) molecules[i].getMyBends();
84			for (int j = 0; j < molecules[i].getNBends(); j++){
85	mmeineke	558	constrained = theArray[j]->is_constrained();
86	tim	725
87			if (constrained){
88			dummy_plug = theArray[j]->get_constraint();
89			temp_con[nConstrained].set_a(dummy_plug->get_a());
90			temp_con[nConstrained].set_b(dummy_plug->get_b());
91			temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
92
93			nConstrained++;
94			constrained = 0;
95	mmeineke	558	}
96			}
97
98	tim	725	theArray = (SRI * *) molecules[i].getMyTorsions();
99			for (int j = 0; j < molecules[i].getNTorsions(); j++){
100	mmeineke	558	constrained = theArray[j]->is_constrained();
101	tim	725
102			if (constrained){
103			dummy_plug = theArray[j]->get_constraint();
104			temp_con[nConstrained].set_a(dummy_plug->get_a());
105			temp_con[nConstrained].set_b(dummy_plug->get_b());
106			temp_con[nConstrained].set_dsqr(dummy_plug->get_dsqr());
107
108			nConstrained++;
109			constrained = 0;
110	mmeineke	558	}
111			}
112			}
113
114	tim	725	if (nConstrained > 0){
115	mmeineke	558	isConstrained = 1;
116
117	tim	725	if (constrainedA != NULL)
118			delete[] constrainedA;
119			if (constrainedB != NULL)
120			delete[] constrainedB;
121			if (constrainedDsqr != NULL)
122			delete[] constrainedDsqr;
123	mmeineke	558
124	tim	725	constrainedA = new int[nConstrained];
125			constrainedB = new int[nConstrained];
126	mmeineke	558	constrainedDsqr = new double[nConstrained];
127	tim	725
128			for (int i = 0; i < nConstrained; i++){
129	mmeineke	558	constrainedA[i] = temp_con[i].get_a();
130			constrainedB[i] = temp_con[i].get_b();
131			constrainedDsqr[i] = temp_con[i].get_dsqr();
132			}
133
134	tim	725
135	mmeineke	558	// save oldAtoms to check for lode balanceing later on.
136	tim	725
137	mmeineke	558	oldAtoms = nAtoms;
138	tim	725
139	mmeineke	558	moving = new int[nAtoms];
140	tim	725	moved = new int[nAtoms];
141	mmeineke	558
142	tim	725	oldPos = new double[nAtoms * 3];
143	mmeineke	558	}
144	tim	725
145	mmeineke	558	delete[] temp_con;
146			}
147
148
149	tim	725	template<typename T> void Integrator<T>::integrate(void){
150	mmeineke	558	int i, j; // loop counters
151
152	tim	725	double runTime = info->run_time;
153			double sampleTime = info->sampleTime;
154			double statusTime = info->statusTime;
155	mmeineke	558	double thermalTime = info->thermalTime;
156	mmeineke	746	double resetTime = info->resetTime;
157	mmeineke	558
158	mmeineke	746
159	mmeineke	558	double currSample;
160			double currThermal;
161			double currStatus;
162	mmeineke	746	double currReset;
163
164	mmeineke	558	int calcPot, calcStress;
165			int isError;
166
167	tim	725	tStats = new Thermo(info);
168			statOut = new StatWriter(info);
169			dumpOut = new DumpWriter(info);
170	mmeineke	558
171	mmeineke	561	atoms = info->atoms;
172	mmeineke	558	DirectionalAtom* dAtom;
173	mmeineke	561
174			dt = info->dt;
175	mmeineke	558	dt2 = 0.5 * dt;
176
177			// initialize the forces before the first step
178
179	tim	677	calcForce(1, 1);
180	tim	733
181	tim	725	if (info->setTemp){
182	tim	677	thermalize();
183	mmeineke	558	}
184	tim	725
185	mmeineke	558	calcPot = 0;
186			calcStress = 0;
187	mmeineke	711	currSample = sampleTime + info->getTime();
188			currThermal = thermalTime+ info->getTime();
189			currStatus = statusTime + info->getTime();
190	mmeineke	746	currReset = resetTime + info->getTime();
191	mmeineke	558
192	tim	725	dumpOut->writeDump(info->getTime());
193			statOut->writeStat(info->getTime());
194	mmeineke	559
195			readyCheck();
196
197			#ifdef IS_MPI
198	tim	725	strcpy(checkPointMsg, "The integrator is ready to go.");
199	mmeineke	559	MPIcheckPoint();
200			#endif // is_mpi
201
202	tim	725	while (info->getTime() < runTime){
203			if ((info->getTime() + dt) >= currStatus){
204	mmeineke	558	calcPot = 1;
205			calcStress = 1;
206			}
207	mmeineke	561
208	tim	725	integrateStep(calcPot, calcStress);
209
210	mmeineke	643	info->incrTime(dt);
211	mmeineke	558
212	tim	725	if (info->setTemp){
213			if (info->getTime() >= currThermal){
214			thermalize();
215			currThermal += thermalTime;
216	mmeineke	558	}
217			}
218
219	tim	725	if (info->getTime() >= currSample){
220			dumpOut->writeDump(info->getTime());
221	mmeineke	558	currSample += sampleTime;
222			}
223
224	tim	725	if (info->getTime() >= currStatus){
225			statOut->writeStat(info->getTime());
226	mmeineke	558	calcPot = 0;
227			calcStress = 0;
228			currStatus += statusTime;
229			}
230	mmeineke	559
231	mmeineke	746	if (info->resetIntegrator){
232			if (info->getTime() >= currReset){
233			this->resetIntegrator();
234			currReset += resetTime;
235			}
236			}
237
238	mmeineke	559	#ifdef IS_MPI
239	tim	725	strcpy(checkPointMsg, "successfully took a time step.");
240	mmeineke	559	MPIcheckPoint();
241			#endif // is_mpi
242	mmeineke	558	}
243
244	mmeineke	643	dumpOut->writeFinal(info->getTime());
245	mmeineke	558
246	mmeineke	561	delete dumpOut;
247			delete statOut;
248	mmeineke	558	}
249
250	tim	725	template<typename T> void Integrator<T>::integrateStep(int calcPot,
251			int calcStress){
252	mmeineke	558	// Position full step, and velocity half step
253	tim	725	preMove();
254	mmeineke	558
255			moveA();
256
257	tim	725
258
259	mmeineke	768
260	mmeineke	614	#ifdef IS_MPI
261	tim	725	strcpy(checkPointMsg, "Succesful moveA\n");
262	mmeineke	614	MPIcheckPoint();
263			#endif // is_mpi
264
265	tim	725
266	mmeineke	558	// calc forces
267
268	tim	725	calcForce(calcPot, calcStress);
269	mmeineke	558
270	mmeineke	614	#ifdef IS_MPI
271	tim	725	strcpy(checkPointMsg, "Succesful doForces\n");
272	mmeineke	614	MPIcheckPoint();
273			#endif // is_mpi
274
275	tim	725
276	mmeineke	558	// finish the velocity half step
277	tim	725
278	mmeineke	558	moveB();
279	tim	725
280
281	mmeineke	768
282	mmeineke	614	#ifdef IS_MPI
283	tim	725	strcpy(checkPointMsg, "Succesful moveB\n");
284	mmeineke	614	MPIcheckPoint();
285			#endif // is_mpi
286	mmeineke	558	}
287
288
289	tim	725	template<typename T> void Integrator<T>::moveA(void){
290	gezelter	600	int i, j;
291	mmeineke	558	DirectionalAtom* dAtom;
292	gezelter	600	double Tb[3], ji[3];
293			double A[3][3], I[3][3];
294	mmeineke	561	double angle;
295	gezelter	600	double vel[3], pos[3], frc[3];
296			double mass;
297	mmeineke	558
298	tim	725	for (i = 0; i < nAtoms; i++){
299			atoms[i]->getVel(vel);
300			atoms[i]->getPos(pos);
301			atoms[i]->getFrc(frc);
302	mmeineke	567
303	gezelter	600	mass = atoms[i]->getMass();
304	mmeineke	594
305	tim	725	for (j = 0; j < 3; j++){
306	gezelter	600	// velocity half step
307	tim	725	vel[j] += (dt2 * frc[j] / mass) * eConvert;
308	gezelter	600	// position whole step
309			pos[j] += dt * vel[j];
310			}
311	mmeineke	594
312	tim	725	atoms[i]->setVel(vel);
313			atoms[i]->setPos(pos);
314	gezelter	600
315	tim	725	if (atoms[i]->isDirectional()){
316			dAtom = (DirectionalAtom *) atoms[i];
317	mmeineke	558
318			// get and convert the torque to body frame
319	mmeineke	597
320	tim	725	dAtom->getTrq(Tb);
321			dAtom->lab2Body(Tb);
322
323	mmeineke	558	// get the angular momentum, and propagate a half step
324	gezelter	600
325	tim	725	dAtom->getJ(ji);
326	gezelter	600
327	tim	725	for (j = 0; j < 3; j++)
328	gezelter	600	ji[j] += (dt2 * Tb[j]) * eConvert;
329	tim	725
330	mmeineke	558	// use the angular velocities to propagate the rotation matrix a
331			// full time step
332	gezelter	600
333			dAtom->getA(A);
334			dAtom->getI(I);
335	tim	725
336	mmeineke	558	// rotate about the x-axis
337	gezelter	600	angle = dt2 * ji[0] / I[0][0];
338	tim	725	this->rotate(1, 2, angle, ji, A);
339	mmeineke	597
340	mmeineke	558	// rotate about the y-axis
341	gezelter	600	angle = dt2 * ji[1] / I[1][1];
342	tim	725	this->rotate(2, 0, angle, ji, A);
343
344	mmeineke	558	// rotate about the z-axis
345	gezelter	600	angle = dt * ji[2] / I[2][2];
346	tim	725	this->rotate(0, 1, angle, ji, A);
347
348	mmeineke	558	// rotate about the y-axis
349	gezelter	600	angle = dt2 * ji[1] / I[1][1];
350	tim	725	this->rotate(2, 0, angle, ji, A);
351
352			// rotate about the x-axis
353	gezelter	600	angle = dt2 * ji[0] / I[0][0];
354	tim	725	this->rotate(1, 2, angle, ji, A);
355	mmeineke	597
356	tim	725	dAtom->setJ(ji);
357			dAtom->setA(A);
358			}
359	mmeineke	558	}
360	mmeineke	768
361			if (nConstrained){
362			constrainA();
363			}
364	mmeineke	558	}
365
366
367	tim	725	template<typename T> void Integrator<T>::moveB(void){
368	gezelter	600	int i, j;
369	mmeineke	558	DirectionalAtom* dAtom;
370	gezelter	600	double Tb[3], ji[3];
371			double vel[3], frc[3];
372			double mass;
373	mmeineke	558
374	tim	725	for (i = 0; i < nAtoms; i++){
375			atoms[i]->getVel(vel);
376			atoms[i]->getFrc(frc);
377	mmeineke	558
378	gezelter	600	mass = atoms[i]->getMass();
379
380	mmeineke	558	// velocity half step
381	tim	725	for (j = 0; j < 3; j++)
382			vel[j] += (dt2 * frc[j] / mass) * eConvert;
383	gezelter	600
384	tim	725	atoms[i]->setVel(vel);
385	mmeineke	597
386	tim	725	if (atoms[i]->isDirectional()){
387			dAtom = (DirectionalAtom *) atoms[i];
388
389	gezelter	600	// get and convert the torque to body frame
390
391	tim	725	dAtom->getTrq(Tb);
392			dAtom->lab2Body(Tb);
393	gezelter	600
394			// get the angular momentum, and propagate a half step
395
396	tim	725	dAtom->getJ(ji);
397	gezelter	600
398	tim	725	for (j = 0; j < 3; j++)
399	gezelter	600	ji[j] += (dt2 * Tb[j]) * eConvert;
400	mmeineke	597
401	tim	725
402			dAtom->setJ(ji);
403	mmeineke	558	}
404			}
405	mmeineke	768
406			if (nConstrained){
407			constrainB();
408			}
409	mmeineke	558	}
410
411	tim	725	template<typename T> void Integrator<T>::preMove(void){
412	gezelter	600	int i, j;
413			double pos[3];
414	mmeineke	558
415	tim	725	if (nConstrained){
416			for (i = 0; i < nAtoms; i++){
417			atoms[i]->getPos(pos);
418	mmeineke	561
419	tim	725	for (j = 0; j < 3; j++){
420			oldPos[3 * i + j] = pos[j];
421	gezelter	600	}
422			}
423	tim	725	}
424	gezelter	600	}
425
426	tim	645	template<typename T> void Integrator<T>::constrainA(){
427	tim	725	int i, j, k;
428	mmeineke	558	int done;
429	gezelter	600	double posA[3], posB[3];
430			double velA[3], velB[3];
431	mmeineke	572	double pab[3];
432			double rab[3];
433	mmeineke	563	int a, b, ax, ay, az, bx, by, bz;
434	mmeineke	558	double rma, rmb;
435			double dx, dy, dz;
436	mmeineke	561	double rpab;
437	mmeineke	558	double rabsq, pabsq, rpabsq;
438			double diffsq;
439			double gab;
440			int iteration;
441
442	tim	725	for (i = 0; i < nAtoms; i++){
443	mmeineke	558	moving[i] = 0;
444	tim	725	moved[i] = 1;
445	mmeineke	558	}
446	mmeineke	567
447	mmeineke	558	iteration = 0;
448			done = 0;
449	tim	725	while (!done && (iteration < maxIteration)){
450	mmeineke	558	done = 1;
451	tim	725	for (i = 0; i < nConstrained; i++){
452	mmeineke	558	a = constrainedA[i];
453			b = constrainedB[i];
454	mmeineke	563
455	tim	725	ax = (a * 3) + 0;
456			ay = (a * 3) + 1;
457			az = (a * 3) + 2;
458	mmeineke	563
459	tim	725	bx = (b * 3) + 0;
460			by = (b * 3) + 1;
461			bz = (b * 3) + 2;
462
463			if (moved[a] \|\| moved[b]){
464			atoms[a]->getPos(posA);
465			atoms[b]->getPos(posB);
466
467			for (j = 0; j < 3; j++)
468	gezelter	600	pab[j] = posA[j] - posB[j];
469	mmeineke	567
470	tim	725	//periodic boundary condition
471	mmeineke	567
472	tim	725	info->wrapVector(pab);
473	mmeineke	572
474	tim	725	pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
475	mmeineke	558
476	tim	725	rabsq = constrainedDsqr[i];
477			diffsq = rabsq - pabsq;
478	mmeineke	567
479	tim	725	// the original rattle code from alan tidesley
480			if (fabs(diffsq) > (tol * rabsq * 2)){
481			rab[0] = oldPos[ax] - oldPos[bx];
482			rab[1] = oldPos[ay] - oldPos[by];
483			rab[2] = oldPos[az] - oldPos[bz];
484	mmeineke	558
485	tim	725	info->wrapVector(rab);
486	mmeineke	567
487	tim	725	rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
488	mmeineke	558
489	tim	725	rpabsq = rpab * rpab;
490	mmeineke	558
491	mmeineke	563
492	tim	725	if (rpabsq < (rabsq * -diffsq)){
493	mmeineke	558	#ifdef IS_MPI
494	tim	725	a = atoms[a]->getGlobalIndex();
495			b = atoms[b]->getGlobalIndex();
496	mmeineke	558	#endif //is_mpi
497	tim	725	sprintf(painCave.errMsg,
498			"Constraint failure in constrainA at atom %d and %d.\n", a,
499			b);
500			painCave.isFatal = 1;
501			simError();
502			}
503	mmeineke	558
504	tim	725	rma = 1.0 / atoms[a]->getMass();
505			rmb = 1.0 / atoms[b]->getMass();
506	mmeineke	567
507	tim	725	gab = diffsq / (2.0 * (rma + rmb) * rpab);
508	mmeineke	567
509	mmeineke	572	dx = rab[0] * gab;
510			dy = rab[1] * gab;
511			dz = rab[2] * gab;
512	mmeineke	558
513	tim	725	posA[0] += rma * dx;
514			posA[1] += rma * dy;
515			posA[2] += rma * dz;
516	mmeineke	558
517	tim	725	atoms[a]->setPos(posA);
518	mmeineke	558
519	tim	725	posB[0] -= rmb * dx;
520			posB[1] -= rmb * dy;
521			posB[2] -= rmb * dz;
522	gezelter	600
523	tim	725	atoms[b]->setPos(posB);
524	gezelter	600
525	mmeineke	558	dx = dx / dt;
526			dy = dy / dt;
527			dz = dz / dt;
528
529	tim	725	atoms[a]->getVel(velA);
530	mmeineke	558
531	tim	725	velA[0] += rma * dx;
532			velA[1] += rma * dy;
533			velA[2] += rma * dz;
534	mmeineke	558
535	tim	725	atoms[a]->setVel(velA);
536	gezelter	600
537	tim	725	atoms[b]->getVel(velB);
538	gezelter	600
539	tim	725	velB[0] -= rmb * dx;
540			velB[1] -= rmb * dy;
541			velB[2] -= rmb * dz;
542	gezelter	600
543	tim	725	atoms[b]->setVel(velB);
544	gezelter	600
545	tim	725	moving[a] = 1;
546			moving[b] = 1;
547			done = 0;
548			}
549	mmeineke	558	}
550			}
551	tim	725
552			for (i = 0; i < nAtoms; i++){
553	mmeineke	558	moved[i] = moving[i];
554			moving[i] = 0;
555			}
556
557			iteration++;
558			}
559
560	tim	725	if (!done){
561			sprintf(painCave.errMsg,
562			"Constraint failure in constrainA, too many iterations: %d\n",
563			iteration);
564	mmeineke	558	painCave.isFatal = 1;
565			simError();
566			}
567	mmeineke	768
568	mmeineke	558	}
569
570	tim	725	template<typename T> void Integrator<T>::constrainB(void){
571			int i, j, k;
572	mmeineke	558	int done;
573	gezelter	600	double posA[3], posB[3];
574			double velA[3], velB[3];
575	mmeineke	558	double vxab, vyab, vzab;
576	mmeineke	572	double rab[3];
577	mmeineke	563	int a, b, ax, ay, az, bx, by, bz;
578	mmeineke	558	double rma, rmb;
579			double dx, dy, dz;
580			double rabsq, pabsq, rvab;
581			double diffsq;
582			double gab;
583			int iteration;
584
585	tim	725	for (i = 0; i < nAtoms; i++){
586	mmeineke	558	moving[i] = 0;
587			moved[i] = 1;
588			}
589
590			done = 0;
591	mmeineke	561	iteration = 0;
592	tim	725	while (!done && (iteration < maxIteration)){
593	mmeineke	567	done = 1;
594
595	tim	725	for (i = 0; i < nConstrained; i++){
596	mmeineke	558	a = constrainedA[i];
597			b = constrainedB[i];
598
599	tim	725	ax = (a * 3) + 0;
600			ay = (a * 3) + 1;
601			az = (a * 3) + 2;
602	mmeineke	563
603	tim	725	bx = (b * 3) + 0;
604			by = (b * 3) + 1;
605			bz = (b * 3) + 2;
606	mmeineke	563
607	tim	725	if (moved[a] \|\| moved[b]){
608			atoms[a]->getVel(velA);
609			atoms[b]->getVel(velB);
610	mmeineke	558
611	tim	725	vxab = velA[0] - velB[0];
612			vyab = velA[1] - velB[1];
613			vzab = velA[2] - velB[2];
614	gezelter	600
615	tim	725	atoms[a]->getPos(posA);
616			atoms[b]->getPos(posB);
617	gezelter	600
618	tim	725	for (j = 0; j < 3; j++)
619	gezelter	600	rab[j] = posA[j] - posB[j];
620	mmeineke	558
621	tim	725	info->wrapVector(rab);
622	mmeineke	558
623	tim	725	rma = 1.0 / atoms[a]->getMass();
624			rmb = 1.0 / atoms[b]->getMass();
625	mmeineke	558
626	tim	725	rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
627	gezelter	600
628	tim	725	gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);
629	gezelter	600
630	tim	725	if (fabs(gab) > tol){
631			dx = rab[0] * gab;
632			dy = rab[1] * gab;
633			dz = rab[2] * gab;
634
635			velA[0] += rma * dx;
636			velA[1] += rma * dy;
637			velA[2] += rma * dz;
638
639			atoms[a]->setVel(velA);
640
641			velB[0] -= rmb * dx;
642			velB[1] -= rmb * dy;
643			velB[2] -= rmb * dz;
644
645			atoms[b]->setVel(velB);
646
647			moving[a] = 1;
648			moving[b] = 1;
649			done = 0;
650			}
651	mmeineke	558	}
652			}
653
654	tim	725	for (i = 0; i < nAtoms; i++){
655	mmeineke	558	moved[i] = moving[i];
656			moving[i] = 0;
657			}
658	tim	725
659	mmeineke	558	iteration++;
660			}
661
662	tim	725	if (!done){
663			sprintf(painCave.errMsg,
664			"Constraint failure in constrainB, too many iterations: %d\n",
665			iteration);
666	mmeineke	558	painCave.isFatal = 1;
667			simError();
668	tim	725	}
669	mmeineke	558	}
670
671	tim	725	template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
672			double angle, double ji[3],
673			double A[3][3]){
674			int i, j, k;
675	mmeineke	558	double sinAngle;
676			double cosAngle;
677			double angleSqr;
678			double angleSqrOver4;
679			double top, bottom;
680			double rot[3][3];
681			double tempA[3][3];
682			double tempJ[3];
683
684			// initialize the tempA
685
686	tim	725	for (i = 0; i < 3; i++){
687			for (j = 0; j < 3; j++){
688	gezelter	600	tempA[j][i] = A[i][j];
689	mmeineke	558	}
690			}
691
692			// initialize the tempJ
693
694	tim	725	for (i = 0; i < 3; i++)
695			tempJ[i] = ji[i];
696
697	mmeineke	558	// initalize rot as a unit matrix
698
699			rot[0][0] = 1.0;
700			rot[0][1] = 0.0;
701			rot[0][2] = 0.0;
702
703			rot[1][0] = 0.0;
704			rot[1][1] = 1.0;
705			rot[1][2] = 0.0;
706	tim	725
707	mmeineke	558	rot[2][0] = 0.0;
708			rot[2][1] = 0.0;
709			rot[2][2] = 1.0;
710	tim	725
711	mmeineke	558	// use a small angle aproximation for sin and cosine
712
713	tim	725	angleSqr = angle * angle;
714	mmeineke	558	angleSqrOver4 = angleSqr / 4.0;
715			top = 1.0 - angleSqrOver4;
716			bottom = 1.0 + angleSqrOver4;
717
718			cosAngle = top / bottom;
719			sinAngle = angle / bottom;
720
721			rot[axes1][axes1] = cosAngle;
722			rot[axes2][axes2] = cosAngle;
723
724			rot[axes1][axes2] = sinAngle;
725			rot[axes2][axes1] = -sinAngle;
726	tim	725
727	mmeineke	558	// rotate the momentum acoording to: ji[] = rot[][] * ji[]
728	tim	725
729			for (i = 0; i < 3; i++){
730	mmeineke	558	ji[i] = 0.0;
731	tim	725	for (k = 0; k < 3; k++){
732	mmeineke	558	ji[i] += rot[i][k] * tempJ[k];
733			}
734			}
735
736			// rotate the Rotation matrix acording to:
737			// A[][] = A[][] * transpose(rot[][])
738
739
740	mmeineke	561	// NOte for as yet unknown reason, we are performing the
741	mmeineke	558	// calculation as:
742			// transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
743
744	tim	725	for (i = 0; i < 3; i++){
745			for (j = 0; j < 3; j++){
746	gezelter	600	A[j][i] = 0.0;
747	tim	725	for (k = 0; k < 3; k++){
748			A[j][i] += tempA[i][k] * rot[j][k];
749	mmeineke	558	}
750			}
751			}
752			}
753	tim	677
754	tim	725	template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
755			myFF->doForces(calcPot, calcStress);
756	tim	677	}
757
758			template<typename T> void Integrator<T>::thermalize(){
759	tim	725	tStats->velocitize();
760	tim	677	}
761	tim	763
762			template<typename T> double Integrator<T>::getConservedQuantity(void){
763			return tStats->getTotalE();
764	mmeineke	768	}