--- trunk/OOPSE/libmdtools/Integrator.cpp	2003/06/19 22:02:44	559
+++ trunk/OOPSE/libmdtools/Integrator.cpp	2003/06/23 21:24:03	563
@@ -1,5 +1,6 @@
 #include <iostream>
 #include <cstdlib>
+#include <cmath>
 
 #ifdef IS_MPI
 #include "mpiSimulation.hpp"
@@ -10,7 +11,7 @@
 #include "simError.h"
 
 
-Integrator::Integrator( SimInfo* theInfo, ForceFields* the_ff ){
+Integrator::Integrator( SimInfo *theInfo, ForceFields* the_ff ){
   
   info = theInfo;
   myFF = the_ff;
@@ -33,7 +34,7 @@ Integrator::Integrator( SimInfo* theInfo, ForceFields*
   constrainedDsqr = NULL;
   moving          = NULL;
   moved           = NULL;
-  prePos          = NULL;
+  oldPos          = NULL;
   
   nConstrained = 0;
 
@@ -48,7 +49,7 @@ Integrator::~Integrator() {
     delete[] constrainedDsqr;
     delete[] moving;
     delete[] moved;
-    delete[] prePos;
+    delete[] oldPos;
   }
   
 }
@@ -136,6 +137,9 @@ void Integrator::checkConstraints( void ){
       constrainedA[i] = temp_con[i].get_a();
       constrainedB[i] = temp_con[i].get_b();
       constrainedDsqr[i] = temp_con[i].get_dsqr();
+
+      cerr << "constraint " << constrainedA[i] << " <-> " << constrainedB[i]
+	   << " => " << constrainedDsqr[i] << "\n";
     }
 
     
@@ -146,7 +150,7 @@ void Integrator::checkConstraints( void ){
     moving = new int[nAtoms];
     moved  = new int[nAtoms];
 
-    prePos = new double[nAtoms*3];
+    oldPos = new double[nAtoms*3];
   }
   
   delete[] temp_con;
@@ -156,24 +160,7 @@ void Integrator::integrate( void ){
 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];
-
-  double dt          = info->dt;
   double runTime     = info->run_time;
   double sampleTime  = info->sampleTime;
   double statusTime  = info->statusTime;
@@ -187,12 +174,16 @@ void Integrator::integrate( void ){
   int calcPot, calcStress;
   int isError;
 
-  tStats   = new Thermo( info );
-  e_out    = new StatWriter( info );
-  dump_out = new DumpWriter( info );
 
-  Atom** atoms = info->atoms;
+
+  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
@@ -204,8 +195,8 @@ void Integrator::integrate( void ){
     tStats->velocitize();
   }
   
-  dump_out->writeDump( 0.0 );
-  e_out->writeStat( 0.0 );
+  dumpOut->writeDump( 0.0 );
+  statOut->writeStat( 0.0 );
   
   calcPot     = 0;
   calcStress  = 0;
@@ -223,13 +214,20 @@ void Integrator::integrate( void ){
   MPIcheckPoint();
 #endif // is_mpi
 
+
+  pos  = Atom::getPosArray();
+  vel  = Atom::getVelArray();
+  frc  = Atom::getFrcArray();
+  trq  = Atom::getTrqArray();
+  Amat = Atom::getAmatArray();
+
   while( currTime < runTime ){
 
     if( (currTime+dt) >= currStatus ){
       calcPot = 1;
       calcStress = 1;
     }
-    
+
     integrateStep( calcPot, calcStress );
       
     currTime += dt;
@@ -242,12 +240,12 @@ void Integrator::integrate( void ){
     }
 
     if( currTime >= currSample ){
-      dump_out->writeDump( currTime );
+      dumpOut->writeDump( currTime );
       currSample += sampleTime;
     }
 
     if( currTime >= currStatus ){ 
-      e_out->writeStat( time * dt ); 
+      statOut->writeStat( currTime ); 
       calcPot = 0; 
       calcStress = 0;
       currStatus += statusTime;
@@ -261,17 +259,19 @@ void Integrator::integrate( void ){
 
   }
 
-  dump_out->writeFinal();
+  dumpOut->writeFinal();
 
-  delete dump_out;
-  delete e_out;
+  delete dumpOut;
+  delete statOut;
 }
 
 void Integrator::integrateStep( int calcPot, int calcStress ){
 
+
+      
   // Position full step, and velocity half step
 
-  //preMove();
+  preMove();
   moveA();
   if( nConstrained ) constrainA();
 
@@ -294,6 +294,7 @@ void Integrator::moveA( void ){
   DirectionalAtom* dAtom;
   double Tb[3];
   double ji[3];
+  double angle;
 
   for( i=0; i<nAtoms; i++ ){
     atomIndex = i * 3;
@@ -331,23 +332,23 @@ void Integrator::moveA( void ){
       
       // rotate about the x-axis      
       angle = dt2 * ji[0] / dAtom->getIxx();
-      this->rotate( 1, 2, angle, ji, &aMat[aMatIndex] ); 
+      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] );
+      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] );
+      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] );
+      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] );
+      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
       
       dAtom->setJx( ji[0] );
       dAtom->setJy( ji[1] );
@@ -391,10 +392,6 @@ void Integrator::moveB( void ){
       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] );
@@ -407,22 +404,7 @@ 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];
   }
 }  
@@ -433,9 +415,10 @@ void Integrator::constrainA(){
   int done;
   double pxab, pyab, pzab;
   double rxab, ryab, rzab;
-  int a, b;
+  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;
@@ -459,49 +442,62 @@ void Integrator::constrainA(){
 
       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] ){
 	
-	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];
+	pxab = pos[ax] - pos[bx];
+	pyab = pos[ay] - pos[by];
+	pzab = pos[az] - pos[bz];
 
 	//periodic boundary condition
 	pxab = pxab - info->box_x * copysign(1, pxab) 
-	  * int(pxab / info->box_x + 0.5);
+	  * (int)( fabs(pxab / info->box_x) + 0.5);
 	pyab = pyab - info->box_y * copysign(1, pyab) 
-	  * int(pyab / info->box_y + 0.5);
+	  * (int)( fabs(pyab / info->box_y) + 0.5);
 	pzab = pzab - info->box_z * copysign(1, pzab) 
-	  * int(pzab / info->box_z + 0.5);
+	  * (int)( fabs(pzab / info->box_z) + 0.5);
       
 	pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
-	rabsq = constraintedDsqr[i];
+	rabsq = constrainedDsqr[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];
+	if (fabs(diffsq) > (tol*rabsq*2)) {
+	  rxab = oldPos[ax] - oldPos[bx];
+	  ryab = oldPos[ay] - oldPos[by];
+	  rzab = oldPos[az] - oldPos[bz];
  
 	  rxab = rxab - info->box_x * copysign(1, rxab) 
-	    * int(rxab / info->box_x + 0.5);
+	    * (int)( fabs(rxab / info->box_x) + 0.5);
  	  ryab = ryab - info->box_y * copysign(1, ryab) 
-	    * int(ryab / info->box_y + 0.5);
+	    * (int)( fabs(ryab / info->box_y) + 0.5);
 	  rzab = rzab - info->box_z * copysign(1, rzab) 
-	    * int(rzab / info->box_z + 0.5);
+	    * (int)( fabs(rzab / info->box_z) + 0.5);
 
 	  rpab = rxab * pxab + ryab * pyab + rzab * pzab;
 	  rpabsq = rpab * rpab;
 
 
 	  if (rpabsq < (rabsq * -diffsq)){
+
+	    cerr << "rpabsq = " << rpabsq << ", rabsq = " << rabsq 
+		 << ", -diffsq = " << -diffsq << "\n";
+
 #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.",
+		     "Constraint failure in constrainA at atom %d and %d.\n",
 		     a, b );
 	    painCave.isFatal = 1;
 	    simError();
@@ -515,25 +511,25 @@ void Integrator::constrainA(){
           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[ax] += rma * dx;
+	  pos[ay] += rma * dy;
+	  pos[az] += rma * dz;
 
-	  pos[3*b+0] -= rmb * dx;
-	  pos[3*b+1] -= rmb * dy;
-	  pos[3*b+2] -= rmb * dz;
+	  pos[bx] -= rmb * dx;
+	  pos[by] -= rmb * dy;
+	  pos[bz] -= 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[ax] += rma * dx;
+	  vel[ay] += rma * dy;
+	  vel[az] += rma * dz;
 
-	  vel[3*b+0] -= rmb * dx;
-	  vel[3*b+1] -= rmb * dy;
-	  vel[3*b+2] -= rmb * dz;
+	  vel[bx] -= rmb * dx;
+	  vel[by] -= rmb * dy;
+	  vel[bz] -= rmb * dz;
 
 	  moving[a] = 1;
 	  moving[b] = 1;
@@ -549,13 +545,14 @@ void Integrator::constrainA(){
     }
 
     iteration++;
+    cerr << "iterainA = " << iteration << "\n";
   }
 
   if( !done ){
 
-    sprintf( painCae.errMsg,
+    sprintf( painCave.errMsg,
 	     "Constraint failure in constrainA, too many iterations: %d\n",
-	     iterations );
+	     iteration );
     painCave.isFatal = 1;
     simError();
   }
@@ -568,7 +565,7 @@ void Integrator::constrainB( void ){
   int done;
   double vxab, vyab, vzab;
   double rxab, ryab, rzab;
-  int a, b;
+  int a, b, ax, ay, az, bx, by, bz;
   double rma, rmb;
   double dx, dy, dz;
   double rabsq, pabsq, rvab;
@@ -576,12 +573,13 @@ void Integrator::constrainB( void ){
   double gab;
   int iteration;
 
-  for(i=0; i<nAtom; i++){
+  for(i=0; i<nAtoms; i++){
     moving[i] = 0;
     moved[i] = 1;
   }
 
   done = 0;
+  iteration = 0;
   while( !done && (iteration < maxIteration ) ){
 
     for(i=0; i<nConstrained; i++){
@@ -589,29 +587,37 @@ void Integrator::constrainB( void ){
       a = constrainedA[i];
       b = constrainedB[i];
 
+      ax = 3*a +0;
+      ay = 3*a +1;
+      az = 3*a +2;
+
+      bx = 3*b +0;
+      by = 3*b +1;
+      bz = 3*b +2;
+
       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];
+	vxab = vel[ax] - vel[bx];
+	vyab = vel[ay] - vel[by];
+	vzab = vel[az] - vel[bz];
 
-	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 = pos[ax] - pos[bx];
+	ryab = pos[ay] - pos[by];
+	rzab = pos[az] - pos[bz];
 	
 	rxab = rxab - info->box_x * copysign(1, rxab) 
-	  * int(rxab / info->box_x + 0.5);
+	  * (int)( fabs(rxab / info->box_x) + 0.5);
 	ryab = ryab - info->box_y * copysign(1, ryab) 
-	  * int(ryab / info->box_y + 0.5);
+	  * (int)( fabs(ryab / info->box_y) + 0.5);
 	rzab = rzab - info->box_z * copysign(1, rzab) 
-	  * int(rzab / info->box_z + 0.5);
+	  * (int)( fabs(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] );
+	gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
 
 	if (fabs(gab) > tol) {
 	  
@@ -619,13 +625,13 @@ void Integrator::constrainB( void ){
 	  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[ax] += rma * dx;
+	  vel[ay] += rma * dy;
+	  vel[az] += rma * dz;
 
-	  vel[3*b+0] -= rmb * dx;
-	  vel[3*b+1] -= rmb * dy;
-	  vel[3*b+2] -= rmb * dz;
+	  vel[bx] -= rmb * dx;
+	  vel[by] -= rmb * dy;
+	  vel[bz] -= rmb * dz;
 	  
 	  moving[a] = 1;
 	  moving[b] = 1;
@@ -645,9 +651,9 @@ void Integrator::constrainB( void ){
   if( !done ){
 
    
-    sprintf( painCae.errMsg,
+    sprintf( painCave.errMsg,
 	     "Constraint failure in constrainB, too many iterations: %d\n",
-	     iterations );
+	     iteration );
     painCave.isFatal = 1;
     simError();
   } 
@@ -661,7 +667,7 @@ void Integrator::rotate( int axes1, int axes2, double 
 
 
 void Integrator::rotate( int axes1, int axes2, double angle, double ji[3], 
-			 double A[3][3] ){
+			 double A[9] ){
 
   int i,j,k;
   double sinAngle;
@@ -677,7 +683,7 @@ void Integrator::rotate( int axes1, int axes2, double 
 
   for(i=0; i<3; i++){
     for(j=0; j<3; j++){
-      tempA[j][i] = A[i][j];
+      tempA[j][i] = A[3*i + j];
     }
   }
 
@@ -728,15 +734,15 @@ void Integrator::rotate( int axes1, int axes2, double 
   //            A[][] = A[][] * transpose(rot[][])
 
 
-  // NOte for as yet unknown reason, we are setting the performing the
+  // 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;
+      A[3*j + i] = 0.0;
       for(k=0; k<3; k++){
-	A[j][i] += tempA[i][k] * rot[j][k];
+	A[3*j + i] += tempA[i][k] * rot[j][k];
       }
     }
   }