OOPSE/libmdtools/Integrator.cpp

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

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

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


Integrator::Integrator( SimInfo *theInfo, ForceFields* the_ff ){
  
  info = theInfo;
  myFF = the_ff;
  isFirst = 1;

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

  // give a little love back to the SimInfo object
  
  if( info->the_integrator != NULL ) delete info->the_integrator;
  info->the_integrator = this;

  nAtoms = info->n_atoms;

  // check for constraints
  
  constrainedA    = NULL;
  constrainedB    = NULL;
  constrainedDsqr = NULL;
  moving          = NULL;
  moved           = NULL;
  oldPos          = NULL;
  
  nConstrained = 0;

  checkConstraints();
}

Integrator::~Integrator() {
  
  if( nConstrained ){
    delete[] constrainedA;
    delete[] constrainedB;
    delete[] constrainedDsqr;
    delete[] moving;
    delete[] moved;
    delete[] oldPos;
  }
  
}

void Integrator::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;
}


void Integrator::integrate( void ){

  int i, j;                         // loop counters

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

  double currSample;
  double currThermal;
  double currStatus;
  double currTime;

  int calcPot, calcStress;
  int isError;


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

  atoms = info->atoms;
  DirectionalAtom* dAtom;

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

  // initialize the forces before the first step

  myFF->doForces(1,1);
  
  if( info->setTemp ){
    
    tStats->velocitize();
  }
  
  dumpOut->writeDump( 0.0 );
  statOut->writeStat( 0.0 );
  
  calcPot     = 0;
  calcStress  = 0;
  currSample  = sampleTime;
  currThermal = thermalTime;
  currStatus  = statusTime;
  currTime    = 0.0;;


  readyCheck();

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


  pos  = Atom::getPosArray();
  vel  = Atom::getVelArray();
  frc  = Atom::getFrcArray();
  trq  = Atom::getTrqArray();
  Amat = Atom::getAmatArray();

  while( currTime < runTime ){

    if( (currTime+dt) >= currStatus ){
      calcPot = 1;
      calcStress = 1;
    }

    integrateStep( calcPot, calcStress );
      
    currTime += dt;

    if( info->setTemp ){
      if( currTime >= currThermal ){
        tStats->velocitize();
        currThermal += thermalTime;
      }
    }

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

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

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

  }

  dumpOut->writeFinal();

  delete dumpOut;
  delete statOut;
}

void Integrator::integrateStep( int calcPot, int calcStress ){


  // Position full step, and velocity half step

  preMove();
  moveA();
  if( nConstrained ) constrainA();

  // calc forces

  myFF->doForces(calcPot,calcStress);

  // finish the velocity  half step
  
  moveB();
  if( nConstrained ) constrainB();
   
}


void Integrator::moveA( void ){
  
  int i,j,k;
  int atomIndex, aMatIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];
  double angle;

  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;
    aMatIndex = i * 9;
    
    // velocity half step
    for( j=atomIndex; j<(atomIndex+3); j++ )
      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;

    // position whole step    
    for( j=atomIndex; j<(atomIndex+3); j++ )
      pos[j] += dt * vel[j];

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

      dAtom = (DirectionalAtom *)atoms[i];
          
      // get and convert the torque to body frame
      
      Tb[0] = dAtom->getTx();
      Tb[1] = dAtom->getTy();
      Tb[2] = dAtom->getTz();
      
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and propagate a half step
      
      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
      
      // use the angular velocities to propagate the rotation matrix a
      // full time step
      
      // rotate about the x-axis      
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); 
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
      
      // rotate about the z-axis
      angle = dt * ji[2] / dAtom->getIzz();
      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
    
  }
}


void Integrator::moveB( void ){
  int i,j,k;
  int atomIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];

  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;

    // velocity half step
    for( j=atomIndex; j<(atomIndex+3); j++ )
      vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;

    if( atoms[i]->isDirectional() ){
      
      dAtom = (DirectionalAtom *)atoms[i];
      
      // get and convert the torque to body frame
      
      Tb[0] = dAtom->getTx();
      Tb[1] = dAtom->getTy();
      Tb[2] = dAtom->getTz();
      
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and complete the angular momentum
      // half step
      
      ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
      ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * eConvert;
      ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * eConvert;
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
  }

}

void Integrator::preMove( void ){
  int i;

  if( nConstrained ){

//    if( oldAtoms != nAtoms ){
      
//       // save oldAtoms to check for lode balanceing later on.
      
//       oldAtoms = nAtoms;
      
//       delete[] moving;
//       delete[] moved;
//       delete[] oldPos;
      
//       moving = new int[nAtoms];
//       moved  = new int[nAtoms];
      
//       oldPos = new double[nAtoms*3];
//     }
  
    for(i=0; i<(nAtoms*3); i++) oldPos[i] = pos[i];
  }
}  

void Integrator::constrainA(){

  int i,j,k;
  int done;
  double pxab, pyab, pzab;
  double rxab, ryab, rzab;
  int a, b;
  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];
    
      if( moved[a] || moved[b] ){
        
        pxab = pos[3*a+0] - pos[3*b+0];
        pyab = pos[3*a+1] - pos[3*b+1];
        pzab = pos[3*a+2] - pos[3*b+2];

        //periodic boundary condition
        pxab = pxab - info->box_x * copysign(1, pxab) 
          * int( fabs(pxab) / info->box_x + 0.5);
        pyab = pyab - info->box_y * copysign(1, pyab) 
          * int( fabs(pyab) / info->box_y + 0.5);
        pzab = pzab - info->box_z * copysign(1, pzab) 
          * int( fabs(pzab) / info->box_z + 0.5);
      
        pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
        rabsq = constrainedDsqr[i];
        diffsq = pabsq - rabsq;

        // the original rattle code from alan tidesley
        if (fabs(diffsq) > tol*rabsq*2) {
          rxab = oldPos[3*a+0] - oldPos[3*b+0];
          ryab = oldPos[3*a+1] - oldPos[3*b+1];
          rzab = oldPos[3*a+2] - oldPos[3*b+2];
 
          rxab = rxab - info->box_x * copysign(1, rxab) 
            * int( fabs(rxab) / info->box_x + 0.5);
          ryab = ryab - info->box_y * copysign(1, ryab) 
            * int( fabs(ryab) / info->box_y + 0.5);
          rzab = rzab - info->box_z * copysign(1, rzab) 
            * int( fabs(rzab) / info->box_z + 0.5);

          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
          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 = rxab * gab;
          dy = ryab * gab;
          dz = rzab * gab;

          pos[3*a+0] += rma * dx;
          pos[3*a+1] += rma * dy;
          pos[3*a+2] += rma * dz;

          pos[3*b+0] -= rmb * dx;
          pos[3*b+1] -= rmb * dy;
          pos[3*b+2] -= rmb * dz;

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

          vel[3*a+0] += rma * dx;
          vel[3*a+1] += rma * dy;
          vel[3*a+2] += rma * dz;

          vel[3*b+0] -= rmb * dx;
          vel[3*b+1] -= rmb * dy;
          vel[3*b+2] -= rmb * dz;

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

}

void Integrator::constrainB( void ){
  
  int i,j,k;
  int done;
  double vxab, vyab, vzab;
  double rxab, ryab, rzab;
  int a, b;
  double rma, rmb;
  double dx, dy, dz;
  double rabsq, pabsq, rvab;
  double diffsq;
  double gab;
  int iteration;

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

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

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

      if( moved[a] || moved[b] ){
        
        vxab = vel[3*a+0] - vel[3*b+0];
        vyab = vel[3*a+1] - vel[3*b+1];
        vzab = vel[3*a+2] - vel[3*b+2];

        rxab = pos[3*a+0] - pos[3*b+0];
        ryab = pos[3*a+1] - pos[3*b+1];
        rzab = pos[3*a+2] - pos[3*b+2];
        
        rxab = rxab - info->box_x * copysign(1, rxab) 
          * int( fabs(rxab) / info->box_x + 0.5);
        ryab = ryab - info->box_y * copysign(1, ryab) 
          * int( fabs(ryab) / info->box_y + 0.5);
        rzab = rzab - info->box_z * copysign(1, rzab) 
          * int( fabs(rzab) / info->box_z + 0.5);

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

        rvab = rxab * vxab + ryab * vyab + rzab * vzab;
          
        gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );

        if (fabs(gab) > tol) {
          
          dx = rxab * gab;
          dy = ryab * gab;
          dz = rzab * gab;
          
          vel[3*a+0] += rma * dx;
          vel[3*a+1] += rma * dy;
          vel[3*a+2] += rma * dz;

          vel[3*b+0] -= rmb * dx;
          vel[3*b+1] -= rmb * dy;
          vel[3*b+2] -= rmb * dz;
          
          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();
  } 

}


void Integrator::rotate( int axes1, int axes2, double angle, double ji[3], 
                         double A[9] ){

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