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"


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

template<typename T> Integrator<T>::~Integrator() {
  
  if( nConstrained ){
    delete[] constrainedA;
    delete[] constrainedB;
    delete[] constrainedDsqr;
    delete[] moving;
    delete[] moved;
    delete[] oldPos;
  }
  
}

template<typename T> void Integrator<T>::checkConstraints( void ){


  isConstrained = 0;

  Constraint *temp_con;
  Constraint *dummy_plug;
  temp_con = new Constraint[info->n_SRI];
  nConstrained = 0;
  int constrained = 0;
  
  SRI** theArray;
  for(int i = 0; i < nMols; i++){
    
    theArray = (SRI**) molecules[i].getMyBonds();
    for(int j=0; j<molecules[i].getNBonds(); j++){
      
      constrained = theArray[j]->is_constrained();

      if(constrained){

        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a( dummy_plug->get_a() );
        temp_con[nConstrained].set_b( dummy_plug->get_b() );
        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
        
        nConstrained++;
        constrained = 0;
      } 
    }

    theArray = (SRI**) molecules[i].getMyBends();
    for(int j=0; j<molecules[i].getNBends(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a( dummy_plug->get_a() );
        temp_con[nConstrained].set_b( dummy_plug->get_b() );
        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
        
        nConstrained++;
        constrained = 0;
      }
    }

    theArray = (SRI**) molecules[i].getMyTorsions();
    for(int j=0; j<molecules[i].getNTorsions(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a( dummy_plug->get_a() );
        temp_con[nConstrained].set_b( dummy_plug->get_b() );
        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
        
        nConstrained++;
        constrained = 0;
      }
    }
  }

  if(nConstrained > 0){
    
    isConstrained = 1;

    if(constrainedA != NULL )    delete[] constrainedA;
    if(constrainedB != NULL )    delete[] constrainedB;
    if(constrainedDsqr != NULL ) delete[] constrainedDsqr;

    constrainedA =    new int[nConstrained];
    constrainedB =    new int[nConstrained];
    constrainedDsqr = new double[nConstrained];
    
    for( int i = 0; i < nConstrained; i++){
      
      constrainedA[i] = temp_con[i].get_a();
      constrainedB[i] = temp_con[i].get_b();
      constrainedDsqr[i] = temp_con[i].get_dsqr();

    }

    
    // save oldAtoms to check for lode balanceing later on.
    
    oldAtoms = nAtoms;
    
    moving = new int[nAtoms];
    moved  = new int[nAtoms];

    oldPos = new double[nAtoms*3];
  }
  
  delete[] temp_con;
}


template<typename T> void Integrator<T>::integrate( void ){

  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;

  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

  calcForce(1, 1);
  // myFF->doForces(1,1);
        
  if( info->setTemp ){
    
    thermalize();
  }
  
  calcPot     = 0;
  calcStress  = 0;
  currSample  = sampleTime;
  currThermal = thermalTime;
  currStatus  = statusTime;

  dumpOut->writeDump( info->getTime() );
  statOut->writeStat( info->getTime() );

  readyCheck();

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

  while( info->getTime() < runTime ){

    if( (info->getTime()+dt) >= currStatus ){
      calcPot = 1;
      calcStress = 1;
    }

    integrateStep( calcPot, calcStress );
      
    info->incrTime(dt);

    if( info->setTemp ){
      if( info->getTime() >= currThermal ){
        thermalize();
        currThermal += thermalTime;
      }
    }

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

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

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

  }

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

  delete dumpOut;
  delete statOut;
}

template<typename T> void Integrator<T>::integrateStep( int calcPot, int calcStress ){


  // Position full step, and velocity half step

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

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

  // calc forces

  calcForce(calcPot,calcStress);

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

  // finish the velocity  half step
  
  moveB();
  if( nConstrained ) constrainB();
  
#ifdef IS_MPI
  strcpy( checkPointMsg, "Succesful moveB\n" );
  MPIcheckPoint();
#endif // is_mpi
  
 
}


template<typename T> void Integrator<T>::moveA( void ){
  
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double A[3][3], I[3][3];
  double angle;
  double vel[3], pos[3], frc[3];
  double mass;

  for( i=0; i<nAtoms; i++ ){

    atoms[i]->getVel( vel );
    atoms[i]->getPos( pos );
    atoms[i]->getFrc( frc );

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

    for (j=0; j < 3; j++) {
      // velocity half step
      vel[j] += ( dt2 * frc[j] / mass ) * eConvert;
      // position whole step
      pos[j] += dt * vel[j];
    }

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

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

      dAtom = (DirectionalAtom *)atoms[i];
          
      // get and convert the torque to body frame
      
      dAtom->getTrq( Tb );
      dAtom->lab2Body( Tb );

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

      dAtom->getJ( ji );

      for (j=0; j < 3; j++) 
        ji[j] += (dt2 * Tb[j]) * eConvert;
      
      // use the angular velocities to propagate the rotation matrix a
      // full time step

      dAtom->getA(A);
      dAtom->getI(I);
    
      // rotate about the x-axis      
      angle = dt2 * ji[0] / I[0][0];
      this->rotate( 1, 2, angle, ji, A ); 

      // rotate about the y-axis
      angle = dt2 * ji[1] / I[1][1];
      this->rotate( 2, 0, angle, ji, A );
      
      // rotate about the z-axis
      angle = dt * ji[2] / I[2][2];
      this->rotate( 0, 1, angle, ji, A);
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / I[1][1];
      this->rotate( 2, 0, angle, ji, A );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / I[0][0];
      this->rotate( 1, 2, angle, ji, A );
      

      dAtom->setJ( ji );
      dAtom->setA( A  );
          
    }    
  }
}


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

  for( i=0; i<nAtoms; i++ ){
 
    atoms[i]->getVel( vel );
    atoms[i]->getFrc( frc );

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

    // velocity half step
    for (j=0; j < 3; j++) 
      vel[j] += ( dt2 * frc[j] / mass ) * eConvert;
    
    atoms[i]->setVel( vel );
 
    if( atoms[i]->isDirectional() ){

      dAtom = (DirectionalAtom *)atoms[i];

      // get and convert the torque to body frame      

      dAtom->getTrq( Tb );
      dAtom->lab2Body( Tb );

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

      dAtom->getJ( ji );

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

      dAtom->setJ( ji );
    }
  }
}

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,k;
  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,k;
  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 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] ){

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