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

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

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

void Integrator::constrainA(){

  int i,j,k;
  int done;
  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] ){
        
        pab[0] = pos[ax] - pos[bx];
        pab[1] = pos[ay] - pos[by];
        pab[2] = pos[az] - pos[bz];

        //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;

          pos[ax] += rma * dx;
          pos[ay] += rma * dy;
          pos[az] += rma * dz;

          pos[bx] -= rmb * dx;
          pos[by] -= rmb * dy;
          pos[bz] -= rmb * dz;

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

          vel[ax] += rma * dx;
          vel[ay] += rma * dy;
          vel[az] += rma * dz;

          vel[bx] -= rmb * dx;
          vel[by] -= rmb * dy;
          vel[bz] -= 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 rab[3];
  int a, b, ax, ay, az, bx, by, bz;
  double rma, rmb;
  double dx, dy, dz;
  double rabsq, pabsq, rvab;
  double diffsq;
  double gab;
  int iteration;

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

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

    done = 1;

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

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

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

      if( moved[a] || moved[b] ){
        
        vxab = vel[ax] - vel[bx];
        vyab = vel[ay] - vel[by];
        vzab = vel[az] - vel[bz];

        rab[0] = pos[ax] - pos[bx];
        rab[1] = pos[ay] - pos[by];
        rab[2] = pos[az] - pos[bz];
        
        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;
          
          vel[ax] += rma * dx;
          vel[ay] += rma * dy;
          vel[az] += rma * dz;

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