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

      std::cerr << "Is the folowing bond constrained \n";
      theArray[j]->printMe();
      
      if(constrained){
        
        std::cerr << "Yes\n";

        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;
      } 
      else std::cerr << "No.\n";
    }

    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

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


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

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

    }
  }  
}

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

}

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

}

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