OOPSE/libmdtools/Symplectic.cpp

#include <iostream>
#include <cstdlib>

#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;
  prePos          = NULL;
  
  nConstrained = 0;

  checkConstraints();
}

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

void Integrator::checkConstraints( void ){


  isConstrained = 0;

  Constraint *temp_con;
  Constraint *dummy_plug;
  temp_con = new Constraint[info->n_SRI];
  nConstrained = 0;
  int constrained = 0;
  
  SRI** theArray;
  for(int i = 0; i < nMols; i++){
    
    theArray = (SRI**) molecules[i].getMyBonds();
    for(int j=0; j<molecules[i].getNBonds(); j++){
      
      constrained = theArray[j]->is_constrained();
      
      if(constrained){
        
        dummy_plug = theArray[j]->get_constraint();
        temp_con[nConstrained].set_a( dummy_plug->get_a() );
        temp_con[nConstrained].set_b( dummy_plug->get_b() );
        temp_con[nConstrained].set_dsqr( dummy_plug->get_dsqr() );
        
        nConstrained++;
        constrained = 0;
      }
    }

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

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

  if(nConstrained > 0){
    
    isConstrained = 1;

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

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

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

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


void Integrator::integrate( void ){

  int i, j;                         // loop counters
  double kE = 0.0;                  // the kinetic energy  
  double rot_kE;
  double trans_kE;
  int tl;                        // the time loop conter
  double dt2;                       // half the dt

  double vx, vy, vz;    // the velocities
  double vx2, vy2, vz2; // the square of the velocities
  double rx, ry, rz;    // the postitions
  
  double ji[3];   // the body frame angular momentum
  double jx2, jy2, jz2; // the square of the angular momentums
  double Tb[3];   // torque in the body frame
  double angle;   // the angle through which to rotate the rotation matrix
  double A[3][3]; // the rotation matrix
  double press[9];

  int time;

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

  int n_loops  = (int)( runTime / dt );
  int sample_n = (int)( sampleTime / dt );
  int status_n = (int)( statusTime / dt );
  int vel_n    = (int)( thermalTime / dt );

  int calcPot, calcStress;
  int isError;

  tStats   = new Thermo( info );
  e_out    = new StatWriter( info );
  dump_out = new DumpWriter( info );

  Atom** atoms = info->atoms;
  DirectionalAtom* dAtom;
  dt2 = 0.5 * dt;

  // initialize the forces before the first step

  myFF->doForces(1,1);
  
  if( info->setTemp ){
    
    tStats->velocitize();
  }
  
  dump_out->writeDump( 0.0 );
  e_out->writeStat( 0.0 );
  
  calcPot = 0;

  for( tl=0; tl<nLoops; tl++){

    integrateStep( calcPot, calcStress );
      
    time = tl + 1;
    
    if( info->setTemp ){
      if( !(time % vel_n) ) tStats->velocitize();
    }
    if( !(time % sample_n) ) dump_out->writeDump( time * dt );
    if( !((time+1) % status_n) ) {
      calcPot = 1;
      calcStress = 1;
    }
    if( !(time % status_n) ){ 
      e_out->writeStat( time * dt ); 
      calcPot = 0; 
      if (!strcasecmp(info->ensemble, "NPT")) calcStress = 1;
      else calcStress = 0;
    }    

   
  }

  dump_out->writeFinal();

  delete dump_out;
  delete e_out;
}


void Integrator::moveA( void ){
  
  int i,j,k;
  int atomIndex, aMatIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[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;

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

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

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


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

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

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

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

}

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

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

void Integrator::constrainA(){

  int i,j,k;
  int done;
  double pxab, pyab, pzab;
  double rxab, ryab, rzab;
  int a, b;
  double rma, rmb;
  double dx, dy, dz;
  double rabsq, pabsq, rpabsq;
  double diffsq;
  double gab;
  int iteration;


  for( i=0; i<nAtoms; i++){
    
    moving[i] = 0;
    moved[i]  = 1;
  }
  
  
  iteration = 0;
  done = 0;
  while( !done && (iteration < maxIteration )){

    done = 1;
    for(i=0; i<nConstrained; i++){

      a = constrainedA[i];
      b = constrainedB[i];
    
      if( moved[a] || moved[b] ){
        
        pxab = pos[3*a+0] - pos[3*b+0];
        pyab = pos[3*a+1] - pos[3*b+1];
        pzab = pos[3*a+2] - pos[3*b+2];

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

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

          rpab = rxab * pxab + ryab * pyab + rzab * pzab;
          rpabsq = rpab * rpab;


          if (rpabsq < (rabsq * -diffsq)){
#ifdef IS_MPI
            a = atoms[a]->getGlobalIndex();
            b = atoms[b]->getGlobalIndex();
#endif //is_mpi
            sprintf( painCave.errMsg,
                     "Constraint failure in constrainA at atom %d and %d\n.",
                     a, b );
            painCave.isFatal = 1;
            simError();
          }

          rma = 1.0 / atoms[a]->getMass();
          rmb = 1.0 / atoms[b]->getMass();
          
          gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
          dx = rxab * gab;
          dy = ryab * gab;
          dz = rzab * gab;

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

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

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

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

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

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

    iteration++;
  }

  if( !done ){

    sprintf( painCae.errMsg,
             "Constraint failure in constrainA, too many iterations: %d\n",
             iterations );
    painCave.isFatal = 1;
    simError();
  }

}

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

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

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

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

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

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

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

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

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

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

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

  if( !done ){

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