OOPSE/libmdtools/Integrator.cpp

#include <iostream>
#include <stdlib.h>
#include <math.h>

#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;
  }

  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){

  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;

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

  atoms = info->atoms;

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

  readyCheck();

  // initialize the forces before the first step

  std::cerr << "Before initial Force calc\n";

  calcForce(1, 1);

  if (nConstrained){
    preMove();
    constrainA();
    calcForce(1, 1);
    constrainB();
    std::cerr << "premove done\n";
  }


  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());


#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;
      }
    }

    std::cerr << "done with time = " << info->getTime() << "\n";

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


  // write the last frame
  dumpOut->writeDump(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 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;

      this->rotationPropagation( dAtom, ji );

      dAtom->setJ(ji);
    }
  }

  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;
  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;
  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 rvab;
  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>::rotationPropagation
( DirectionalAtom* dAtom, double ji[3] ){

  double angle;
  double A[3][3], I[3][3];

  // 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->setA( A  );
}

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