src/integrators/Integrator.cpp

#include <iostream>
#include <stdlib.h>
#include <math.h>
#ifdef IS_MPI
#include "brains/mpiSimulation.hpp"
#include <unistd.h>
#endif //is_mpi

#ifdef PROFILE
#include "profiling/mdProfile.hpp"
#endif // profile

#include "integrators/Integrator.hpp"
#include "utils/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;
  integrableObjects = info->integrableObjects;


  // 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 balancing 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 difference;
  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();

  // remove center of mass drift velocity (in case we passed in a configuration
  // that was drifting
  tStats->removeCOMdrift();

  // initialize the retraints if necessary
  if (info->useSolidThermInt && !info->useLiquidThermInt) {
    myFF->initRestraints();
  }

  // initialize the forces before the first step

  calcForce(1, 1);

  //execute constraint algorithm to make sure at the very beginning the system is constrained  
  if(nConstrained){
    preMove();
    constrainA();
    calcForce(1, 1);
    constrainB();
  }

  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 && !stopIntegrator()){
    difference = info->getTime() + dt - currStatus;
    if (difference > 0 || fabs(difference) < 1e-4 ){
      calcPot = 1;
      calcStress = 1;
    }

#ifdef PROFILE
    startProfile( pro1 );
#endif
    
    integrateStep(calcPot, calcStress);

#ifdef PROFILE
    endProfile( pro1 );

    startProfile( pro2 );
#endif // profile

    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 PROFILE
    endProfile( pro2 );
#endif //profile

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

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

  // dump out a file containing the omega values for the final configuration
  if (info->useSolidThermInt && !info->useLiquidThermInt)
    myFF->dumpzAngle();
  

  delete dumpOut;
  delete statOut;
}

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

#ifdef PROFILE
  startProfile(pro3);
#endif //profile

  //save old state (position, velocity etc)
  preMove();
#ifdef PROFILE
  endProfile(pro3);

  startProfile(pro4);
#endif // profile

  moveA();

#ifdef PROFILE
  endProfile(pro4);
  
  startProfile(pro5);
#endif//profile


#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

#ifdef PROFILE
  endProfile( pro5 );

  startProfile( pro6 );
#endif //profile

  // finish the velocity  half step

  moveB();

#ifdef PROFILE
  endProfile(pro6);
#endif // profile

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


template<typename T> void Integrator<T>::moveA(void){
  size_t i, j;
  DirectionalAtom* dAtom;
  Vector3d Tb;
  Vector3d ji;
  Vector3d vel;
  Vector3d pos;
  Vector3d frc;
  double mass;
  double omega;
 
  for (i = 0; i < integrableObjects.size() ; i++){
    vel = integrableObjects[i]->getVel();
    pos = integrableObjects[i]->getPos();
    integrableObjects[i]->getFrc(frc);
    
    mass = integrableObjects[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];
    }

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

    if (integrableObjects[i]->isDirectional()){

      // get and convert the torque to body frame

      Tb = integrableObjects[i]->getTrq();
      integrableObjects[i]->lab2Body(Tb);

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

      ji = integrableObjects[i]->getJ();

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

      this->rotationPropagation( integrableObjects[i], ji );

      integrableObjects[i]->setJ(ji);
    }
  }

  if(nConstrained)
    constrainA();
}


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

  for (i = 0; i < integrableObjects.size(); i++){
    vel = integrableObjects[i]->getVel();
    integrableObjects[i]->getFrc(frc);

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

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

    integrableObjects[i]->setVel(vel);

    if (integrableObjects[i]->isDirectional()){

      // get and convert the torque to body frame

      Tb = integrableObjects[i]->getTrq();
      integrableObjects[i]->lab2Body(Tb);

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

      ji = integrableObjects[i]->getJ();

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


      integrableObjects[i]->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++){
      pos = atoms[i]->getPos();

      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]){
        posA = atoms[a]->getPos();
        posB = atoms[b]->getPos();

        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;

          velA = atoms[a]->getVel();

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

          atoms[a]->setVel(velA);

          velB = atoms[b]->getVel();

          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]){
        velA = atoms[a]->getVel();
        velB = atoms[b]->getVel();

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

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

        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
( StuntDouble* sd, double ji[3] ){

  double angle;
  double A[3][3], I[3][3];
  int i, j, k;

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

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

  if (sd->isLinear()) {
    i = sd->linearAxis();
    j = (i+1)%3;
    k = (i+2)%3;
    
    angle = dt2 * ji[j] / I[j][j];
    this->rotate( k, i, angle, ji, A );

    angle = dt * ji[k] / I[k][k];
    this->rotate( i, j, angle, ji, A);

    angle = dt2 * ji[j] / I[j][j];
    this->rotate( k, i, angle, ji, A );

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