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;
  int i;
  int localIndex;

#ifdef IS_MPI
  int which_node;
#endif // is_mpi
  
  vector<StuntDouble*> particles;
  string inAngle;

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

  if (info->useSolidThermInt && !info->useLiquidThermInt) {
    restOut = new RestraintWriter(info);
    initRestraints = new RestraintReader(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) {
    initRestraints->zeroZangle();
    inAngle = info->zAngleName + "_in";
    initRestraints->readZangle(inAngle.c_str());
    initRestraints->readIdealCrystal();
  }

  // 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());
  if (info->useSolidThermInt && !info->useLiquidThermInt)
    restOut->writeZangle(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());
      // write a zAng file to coincide with each dump or eor file
      if (info->useSolidThermInt && !info->useLiquidThermInt)
        restOut->writeZangle(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());

  // write the file containing the omega values of the final configuration
  if (info->useSolidThermInt && !info->useLiquidThermInt){
    restOut->writeZangle(info->getTime());
    restOut->writeZangle(info->getTime(), inAngle.c_str());
  }

  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;
  double Tb[3], ji[3];
  double vel[3], pos[3], frc[3];
  double mass;
  double omega;
 
  for (i = 0; i < integrableObjects.size() ; i++){
    integrableObjects[i]->getVel(vel);
    integrableObjects[i]->getPos(pos);
    integrableObjects[i]->getFrc(frc);
    //    std::cerr << "f = " << frc[0] << "\t" << frc[1] << "\t" << frc[2] << "\n";
    
    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

      integrableObjects[i]->getTrq(Tb);

      //      std::cerr << "t = " << Tb[0] << "\t" << Tb[1] << "\t" << Tb[2] << "\n";
      integrableObjects[i]->lab2Body(Tb);

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

      integrableObjects[i]->getJ(ji);

      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++){
    integrableObjects[i]->getVel(vel);
    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

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

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

      integrableObjects[i]->getJ(ji);

      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++){
      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
( 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:	224
Committed:	Wed Nov 24 18:06:14 2004 UTC (20 years, 11 months ago) by chrisfen
File size:	19674 byte(s)
Log Message:	improved restraints
#	User	Rev	Content
1	gezelter	2	#include <iostream>
2			#include <stdlib.h>
3			#include <math.h>
4			#ifdef IS_MPI
5	tim	3	#include "brains/mpiSimulation.hpp"
6	gezelter	2	#include <unistd.h>
7			#endif //is_mpi
8
9			#ifdef PROFILE
10	tim	3	#include "profiling/mdProfile.hpp"
11	gezelter	2	#endif // profile
12
13	tim	3	#include "integrators/Integrator.hpp"
14			#include "utils/simError.h"
15	gezelter	2
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	chrisfen	221	int i;
175			int localIndex;
176	gezelter	2
177	chrisfen	221	#ifdef IS_MPI
178			int which_node;
179			#endif // is_mpi
180
181			vector<StuntDouble*> particles;
182			string inAngle;
183
184	gezelter	2	tStats = new Thermo(info);
185			statOut = new StatWriter(info);
186			dumpOut = new DumpWriter(info);
187
188	chrisfen	221	if (info->useSolidThermInt && !info->useLiquidThermInt) {
189			restOut = new RestraintWriter(info);
190			initRestraints = new RestraintReader(info);
191			}
192
193	gezelter	2	atoms = info->atoms;
194
195			dt = info->dt;
196			dt2 = 0.5 * dt;
197
198			readyCheck();
199
200			// remove center of mass drift velocity (in case we passed in a configuration
201			// that was drifting
202			tStats->removeCOMdrift();
203
204			// initialize the retraints if necessary
205			if (info->useSolidThermInt && !info->useLiquidThermInt) {
206	chrisfen	221	initRestraints->zeroZangle();
207			inAngle = info->zAngleName + "_in";
208			initRestraints->readZangle(inAngle.c_str());
209			initRestraints->readIdealCrystal();
210	gezelter	2	}
211
212			// initialize the forces before the first step
213			calcForce(1, 1);
214
215			//execute constraint algorithm to make sure at the very beginning the system is constrained
216			if(nConstrained){
217			preMove();
218			constrainA();
219			calcForce(1, 1);
220			constrainB();
221			}
222
223			if (info->setTemp){
224			thermalize();
225			}
226
227			calcPot = 0;
228			calcStress = 0;
229			currSample = sampleTime + info->getTime();
230			currThermal = thermalTime+ info->getTime();
231			currStatus = statusTime + info->getTime();
232			currReset = resetTime + info->getTime();
233
234			dumpOut->writeDump(info->getTime());
235			statOut->writeStat(info->getTime());
236	chrisfen	224	if (info->useSolidThermInt && !info->useLiquidThermInt)
237			restOut->writeZangle(info->getTime());
238	gezelter	2
239			#ifdef IS_MPI
240			strcpy(checkPointMsg, "The integrator is ready to go.");
241			MPIcheckPoint();
242			#endif // is_mpi
243
244			while (info->getTime() < runTime && !stopIntegrator()){
245			difference = info->getTime() + dt - currStatus;
246			if (difference > 0 \|\| fabs(difference) < 1e-4 ){
247			calcPot = 1;
248			calcStress = 1;
249			}
250
251			#ifdef PROFILE
252			startProfile( pro1 );
253			#endif
254
255			integrateStep(calcPot, calcStress);
256
257			#ifdef PROFILE
258			endProfile( pro1 );
259
260			startProfile( pro2 );
261			#endif // profile
262
263			info->incrTime(dt);
264
265			if (info->setTemp){
266			if (info->getTime() >= currThermal){
267			thermalize();
268			currThermal += thermalTime;
269			}
270			}
271
272			if (info->getTime() >= currSample){
273			dumpOut->writeDump(info->getTime());
274	chrisfen	221	// write a zAng file to coincide with each dump or eor file
275			if (info->useSolidThermInt && !info->useLiquidThermInt)
276			restOut->writeZangle(info->getTime());
277	gezelter	2	currSample += sampleTime;
278			}
279
280			if (info->getTime() >= currStatus){
281			statOut->writeStat(info->getTime());
282			calcPot = 0;
283			calcStress = 0;
284			currStatus += statusTime;
285			}
286
287			if (info->resetIntegrator){
288			if (info->getTime() >= currReset){
289			this->resetIntegrator();
290			currReset += resetTime;
291			}
292			}
293
294			#ifdef PROFILE
295			endProfile( pro2 );
296			#endif //profile
297
298			#ifdef IS_MPI
299			strcpy(checkPointMsg, "successfully took a time step.");
300			MPIcheckPoint();
301			#endif // is_mpi
302			}
303
304			dumpOut->writeFinal(info->getTime());
305
306	chrisfen	221	// write the file containing the omega values of the final configuration
307			if (info->useSolidThermInt && !info->useLiquidThermInt){
308			restOut->writeZangle(info->getTime());
309			restOut->writeZangle(info->getTime(), inAngle.c_str());
310			}
311	gezelter	2
312			delete dumpOut;
313			delete statOut;
314			}
315
316			template<typename T> void Integrator<T>::integrateStep(int calcPot,
317			int calcStress){
318			// Position full step, and velocity half step
319
320			#ifdef PROFILE
321			startProfile(pro3);
322			#endif //profile
323
324			//save old state (position, velocity etc)
325			preMove();
326			#ifdef PROFILE
327			endProfile(pro3);
328
329			startProfile(pro4);
330			#endif // profile
331
332			moveA();
333
334			#ifdef PROFILE
335			endProfile(pro4);
336
337			startProfile(pro5);
338			#endif//profile
339
340
341			#ifdef IS_MPI
342			strcpy(checkPointMsg, "Succesful moveA\n");
343			MPIcheckPoint();
344			#endif // is_mpi
345
346			// calc forces
347			calcForce(calcPot, calcStress);
348
349			#ifdef IS_MPI
350			strcpy(checkPointMsg, "Succesful doForces\n");
351			MPIcheckPoint();
352			#endif // is_mpi
353
354			#ifdef PROFILE
355			endProfile( pro5 );
356
357			startProfile( pro6 );
358			#endif //profile
359
360			// finish the velocity half step
361
362			moveB();
363
364			#ifdef PROFILE
365			endProfile(pro6);
366			#endif // profile
367
368			#ifdef IS_MPI
369			strcpy(checkPointMsg, "Succesful moveB\n");
370			MPIcheckPoint();
371			#endif // is_mpi
372			}
373
374
375			template<typename T> void Integrator<T>::moveA(void){
376			size_t i, j;
377			DirectionalAtom* dAtom;
378			double Tb[3], ji[3];
379			double vel[3], pos[3], frc[3];
380			double mass;
381			double omega;
382
383			for (i = 0; i < integrableObjects.size() ; i++){
384			integrableObjects[i]->getVel(vel);
385			integrableObjects[i]->getPos(pos);
386			integrableObjects[i]->getFrc(frc);
387	chrisfen	214	// std::cerr << "f = " << frc[0] << "\t" << frc[1] << "\t" << frc[2] << "\n";
388	gezelter	2
389			mass = integrableObjects[i]->getMass();
390
391			for (j = 0; j < 3; j++){
392			// velocity half step
393			vel[j] += (dt2 * frc[j] / mass) * eConvert;
394			// position whole step
395			pos[j] += dt * vel[j];
396			}
397
398			integrableObjects[i]->setVel(vel);
399			integrableObjects[i]->setPos(pos);
400
401	gezelter	204
402	gezelter	2	if (integrableObjects[i]->isDirectional()){
403
404			// get and convert the torque to body frame
405
406			integrableObjects[i]->getTrq(Tb);
407	gezelter	204
408	chrisfen	214	// std::cerr << "t = " << Tb[0] << "\t" << Tb[1] << "\t" << Tb[2] << "\n";
409	gezelter	2	integrableObjects[i]->lab2Body(Tb);
410
411			// get the angular momentum, and propagate a half step
412
413			integrableObjects[i]->getJ(ji);
414
415			for (j = 0; j < 3; j++)
416			ji[j] += (dt2 * Tb[j]) * eConvert;
417
418			this->rotationPropagation( integrableObjects[i], ji );
419
420			integrableObjects[i]->setJ(ji);
421			}
422			}
423
424			if(nConstrained)
425			constrainA();
426			}
427
428
429			template<typename T> void Integrator<T>::moveB(void){
430			int i, j;
431			double Tb[3], ji[3];
432			double vel[3], frc[3];
433			double mass;
434
435			for (i = 0; i < integrableObjects.size(); i++){
436			integrableObjects[i]->getVel(vel);
437			integrableObjects[i]->getFrc(frc);
438
439			mass = integrableObjects[i]->getMass();
440
441			// velocity half step
442			for (j = 0; j < 3; j++)
443			vel[j] += (dt2 * frc[j] / mass) * eConvert;
444
445			integrableObjects[i]->setVel(vel);
446
447			if (integrableObjects[i]->isDirectional()){
448
449			// get and convert the torque to body frame
450
451			integrableObjects[i]->getTrq(Tb);
452			integrableObjects[i]->lab2Body(Tb);
453
454			// get the angular momentum, and propagate a half step
455
456			integrableObjects[i]->getJ(ji);
457
458			for (j = 0; j < 3; j++)
459			ji[j] += (dt2 * Tb[j]) * eConvert;
460
461
462			integrableObjects[i]->setJ(ji);
463			}
464			}
465
466			if(nConstrained)
467			constrainB();
468			}
469
470
471			template<typename T> void Integrator<T>::preMove(void){
472			int i, j;
473			double pos[3];
474
475			if (nConstrained){
476			for (i = 0; i < nAtoms; i++){
477			atoms[i]->getPos(pos);
478
479			for (j = 0; j < 3; j++){
480			oldPos[3 * i + j] = pos[j];
481			}
482			}
483			}
484			}
485
486			template<typename T> void Integrator<T>::constrainA(){
487			int i, j;
488			int done;
489			double posA[3], posB[3];
490			double velA[3], velB[3];
491			double pab[3];
492			double rab[3];
493			int a, b, ax, ay, az, bx, by, bz;
494			double rma, rmb;
495			double dx, dy, dz;
496			double rpab;
497			double rabsq, pabsq, rpabsq;
498			double diffsq;
499			double gab;
500			int iteration;
501
502			for (i = 0; i < nAtoms; i++){
503			moving[i] = 0;
504			moved[i] = 1;
505			}
506
507			iteration = 0;
508			done = 0;
509			while (!done && (iteration < maxIteration)){
510			done = 1;
511			for (i = 0; i < nConstrained; i++){
512			a = constrainedA[i];
513			b = constrainedB[i];
514
515			ax = (a * 3) + 0;
516			ay = (a * 3) + 1;
517			az = (a * 3) + 2;
518
519			bx = (b * 3) + 0;
520			by = (b * 3) + 1;
521			bz = (b * 3) + 2;
522
523			if (moved[a] \|\| moved[b]){
524			atoms[a]->getPos(posA);
525			atoms[b]->getPos(posB);
526
527			for (j = 0; j < 3; j++)
528			pab[j] = posA[j] - posB[j];
529
530			//periodic boundary condition
531
532			info->wrapVector(pab);
533
534			pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
535
536			rabsq = constrainedDsqr[i];
537			diffsq = rabsq - pabsq;
538
539			// the original rattle code from alan tidesley
540			if (fabs(diffsq) > (tol * rabsq * 2)){
541			rab[0] = oldPos[ax] - oldPos[bx];
542			rab[1] = oldPos[ay] - oldPos[by];
543			rab[2] = oldPos[az] - oldPos[bz];
544
545			info->wrapVector(rab);
546
547			rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
548
549			rpabsq = rpab * rpab;
550
551
552			if (rpabsq < (rabsq * -diffsq)){
553			#ifdef IS_MPI
554			a = atoms[a]->getGlobalIndex();
555			b = atoms[b]->getGlobalIndex();
556			#endif //is_mpi
557			sprintf(painCave.errMsg,
558			"Constraint failure in constrainA at atom %d and %d.\n", a,
559			b);
560			painCave.isFatal = 1;
561			simError();
562			}
563
564			rma = 1.0 / atoms[a]->getMass();
565			rmb = 1.0 / atoms[b]->getMass();
566
567			gab = diffsq / (2.0 * (rma + rmb) * rpab);
568
569			dx = rab[0] * gab;
570			dy = rab[1] * gab;
571			dz = rab[2] * gab;
572
573			posA[0] += rma * dx;
574			posA[1] += rma * dy;
575			posA[2] += rma * dz;
576
577			atoms[a]->setPos(posA);
578
579			posB[0] -= rmb * dx;
580			posB[1] -= rmb * dy;
581			posB[2] -= rmb * dz;
582
583			atoms[b]->setPos(posB);
584
585			dx = dx / dt;
586			dy = dy / dt;
587			dz = dz / dt;
588
589			atoms[a]->getVel(velA);
590
591			velA[0] += rma * dx;
592			velA[1] += rma * dy;
593			velA[2] += rma * dz;
594
595			atoms[a]->setVel(velA);
596
597			atoms[b]->getVel(velB);
598
599			velB[0] -= rmb * dx;
600			velB[1] -= rmb * dy;
601			velB[2] -= rmb * dz;
602
603			atoms[b]->setVel(velB);
604
605			moving[a] = 1;
606			moving[b] = 1;
607			done = 0;
608			}
609			}
610			}
611
612			for (i = 0; i < nAtoms; i++){
613			moved[i] = moving[i];
614			moving[i] = 0;
615			}
616
617			iteration++;
618			}
619
620			if (!done){
621			sprintf(painCave.errMsg,
622			"Constraint failure in constrainA, too many iterations: %d\n",
623			iteration);
624			painCave.isFatal = 1;
625			simError();
626			}
627
628			}
629
630			template<typename T> void Integrator<T>::constrainB(void){
631			int i, j;
632			int done;
633			double posA[3], posB[3];
634			double velA[3], velB[3];
635			double vxab, vyab, vzab;
636			double rab[3];
637			int a, b, ax, ay, az, bx, by, bz;
638			double rma, rmb;
639			double dx, dy, dz;
640			double rvab;
641			double gab;
642			int iteration;
643
644			for (i = 0; i < nAtoms; i++){
645			moving[i] = 0;
646			moved[i] = 1;
647			}
648
649			done = 0;
650			iteration = 0;
651			while (!done && (iteration < maxIteration)){
652			done = 1;
653
654			for (i = 0; i < nConstrained; i++){
655			a = constrainedA[i];
656			b = constrainedB[i];
657
658			ax = (a * 3) + 0;
659			ay = (a * 3) + 1;
660			az = (a * 3) + 2;
661
662			bx = (b * 3) + 0;
663			by = (b * 3) + 1;
664			bz = (b * 3) + 2;
665
666			if (moved[a] \|\| moved[b]){
667			atoms[a]->getVel(velA);
668			atoms[b]->getVel(velB);
669
670			vxab = velA[0] - velB[0];
671			vyab = velA[1] - velB[1];
672			vzab = velA[2] - velB[2];
673
674			atoms[a]->getPos(posA);
675			atoms[b]->getPos(posB);
676
677			for (j = 0; j < 3; j++)
678			rab[j] = posA[j] - posB[j];
679
680			info->wrapVector(rab);
681
682			rma = 1.0 / atoms[a]->getMass();
683			rmb = 1.0 / atoms[b]->getMass();
684
685			rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
686
687			gab = -rvab / ((rma + rmb) * constrainedDsqr[i]);
688
689			if (fabs(gab) > tol){
690			dx = rab[0] * gab;
691			dy = rab[1] * gab;
692			dz = rab[2] * gab;
693
694			velA[0] += rma * dx;
695			velA[1] += rma * dy;
696			velA[2] += rma * dz;
697
698			atoms[a]->setVel(velA);
699
700			velB[0] -= rmb * dx;
701			velB[1] -= rmb * dy;
702			velB[2] -= rmb * dz;
703
704			atoms[b]->setVel(velB);
705
706			moving[a] = 1;
707			moving[b] = 1;
708			done = 0;
709			}
710			}
711			}
712
713			for (i = 0; i < nAtoms; i++){
714			moved[i] = moving[i];
715			moving[i] = 0;
716			}
717
718			iteration++;
719			}
720
721			if (!done){
722			sprintf(painCave.errMsg,
723			"Constraint failure in constrainB, too many iterations: %d\n",
724			iteration);
725			painCave.isFatal = 1;
726			simError();
727			}
728			}
729
730			template<typename T> void Integrator<T>::rotationPropagation
731			( StuntDouble* sd, double ji[3] ){
732
733			double angle;
734			double A[3][3], I[3][3];
735			int i, j, k;
736
737			// use the angular velocities to propagate the rotation matrix a
738			// full time step
739
740			sd->getA(A);
741			sd->getI(I);
742
743			if (sd->isLinear()) {
744	gezelter	204
745	gezelter	2	i = sd->linearAxis();
746			j = (i+1)%3;
747			k = (i+2)%3;
748	gezelter	204
749	gezelter	2	angle = dt2 * ji[j] / I[j][j];
750			this->rotate( k, i, angle, ji, A );
751
752			angle = dt * ji[k] / I[k][k];
753			this->rotate( i, j, angle, ji, A);
754
755			angle = dt2 * ji[j] / I[j][j];
756			this->rotate( k, i, angle, ji, A );
757
758			} else {
759			// rotate about the x-axis
760			angle = dt2 * ji[0] / I[0][0];
761			this->rotate( 1, 2, angle, ji, A );
762
763			// rotate about the y-axis
764			angle = dt2 * ji[1] / I[1][1];
765			this->rotate( 2, 0, angle, ji, A );
766
767			// rotate about the z-axis
768			angle = dt * ji[2] / I[2][2];
769			sd->addZangle(angle);
770			this->rotate( 0, 1, angle, ji, A);
771
772			// rotate about the y-axis
773			angle = dt2 * ji[1] / I[1][1];
774			this->rotate( 2, 0, angle, ji, A );
775
776			// rotate about the x-axis
777			angle = dt2 * ji[0] / I[0][0];
778			this->rotate( 1, 2, angle, ji, A );
779
780			}
781			sd->setA( A );
782			}
783
784			template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
785			double angle, double ji[3],
786			double A[3][3]){
787			int i, j, k;
788			double sinAngle;
789			double cosAngle;
790			double angleSqr;
791			double angleSqrOver4;
792			double top, bottom;
793			double rot[3][3];
794			double tempA[3][3];
795			double tempJ[3];
796
797			// initialize the tempA
798
799			for (i = 0; i < 3; i++){
800			for (j = 0; j < 3; j++){
801			tempA[j][i] = A[i][j];
802			}
803			}
804
805			// initialize the tempJ
806
807			for (i = 0; i < 3; i++)
808			tempJ[i] = ji[i];
809
810			// initalize rot as a unit matrix
811
812			rot[0][0] = 1.0;
813			rot[0][1] = 0.0;
814			rot[0][2] = 0.0;
815
816			rot[1][0] = 0.0;
817			rot[1][1] = 1.0;
818			rot[1][2] = 0.0;
819
820			rot[2][0] = 0.0;
821			rot[2][1] = 0.0;
822			rot[2][2] = 1.0;
823
824			// use a small angle aproximation for sin and cosine
825
826			angleSqr = angle * angle;
827			angleSqrOver4 = angleSqr / 4.0;
828			top = 1.0 - angleSqrOver4;
829			bottom = 1.0 + angleSqrOver4;
830
831			cosAngle = top / bottom;
832			sinAngle = angle / bottom;
833
834			rot[axes1][axes1] = cosAngle;
835			rot[axes2][axes2] = cosAngle;
836
837			rot[axes1][axes2] = sinAngle;
838			rot[axes2][axes1] = -sinAngle;
839
840			// rotate the momentum acoording to: ji[] = rot[][] * ji[]
841
842			for (i = 0; i < 3; i++){
843			ji[i] = 0.0;
844			for (k = 0; k < 3; k++){
845			ji[i] += rot[i][k] * tempJ[k];
846			}
847			}
848
849			// rotate the Rotation matrix acording to:
850			// A[][] = A[][] * transpose(rot[][])
851
852
853			// NOte for as yet unknown reason, we are performing the
854			// calculation as:
855			// transpose(A[][]) = transpose(A[][]) * transpose(rot[][])
856
857			for (i = 0; i < 3; i++){
858			for (j = 0; j < 3; j++){
859			A[j][i] = 0.0;
860			for (k = 0; k < 3; k++){
861			A[j][i] += tempA[i][k] * rot[j][k];
862			}
863			}
864			}
865			}
866
867			template<typename T> void Integrator<T>::calcForce(int calcPot, int calcStress){
868			myFF->doForces(calcPot, calcStress);
869			}
870
871			template<typename T> void Integrator<T>::thermalize(){
872			tStats->velocitize();
873			}
874
875			template<typename T> double Integrator<T>::getConservedQuantity(void){
876			return tStats->getTotalE();
877			}
878			template<typename T> string Integrator<T>::getAdditionalParameters(void){
879			//By default, return a null string
880			//The reason we use string instead of char* is that if we use char*, we will
881			//return a pointer point to local variable which might cause problem
882			return string();
883			}