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;

  std::cerr << "integ nAtoms = "  << nAtoms << "\n";

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


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

    std::cerr << "calcPot = " << calcPot << "; calcStress = "
              << calcStress << "\n";

    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;
  double A[3][3];


  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;

    std::cerr<< "MoveA vel[" << i << "] = "
             << vel[atomIndex] << "\t"
             << vel[atomIndex+1]<< "\t"
             << vel[atomIndex+2]<< "\n";

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

    std::cerr<< "MoveA pos[" << i << "] = " 
             << pos[atomIndex] << "\t"
             << pos[atomIndex+1]<< "\t"
             << pos[atomIndex+2]<< "\n";

    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
      
          // get the atom's rotation matrix
          
      A[0][0] = dAtom->getAxx();
      A[0][1] = dAtom->getAxy();
      A[0][2] = dAtom->getAxz();
      
      A[1][0] = dAtom->getAyx();
      A[1][1] = dAtom->getAyy();
      A[1][2] = dAtom->getAyz();
      
      A[2][0] = dAtom->getAzx();
      A[2][1] = dAtom->getAzy();
      A[2][2] = dAtom->getAzz();
      
      // rotate about the x-axis      
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, A ); 
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, A );
      
      // rotate about the z-axis
      angle = dt * ji[2] / dAtom->getIzz();
      this->rotate( 0, 1, angle, ji, A );
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, A );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, A );
      
      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;

    std::cerr<< "MoveB vel[" << i << "] = "
             << vel[atomIndex] << "\t"
             << vel[atomIndex+1]<< "\t"
             << vel[atomIndex+2]<< "\n";


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