--- trunk/OOPSE/libmdtools/Integrator.cpp	2003/06/23 21:24:03	563
+++ trunk/OOPSE/libmdtools/Integrator.cpp	2003/07/02 21:26:55	572
@@ -138,8 +138,6 @@ void Integrator::checkConstraints( void ){
       constrainedB[i] = temp_con[i].get_b();
       constrainedDsqr[i] = temp_con[i].get_dsqr();
 
-      cerr << "constraint " << constrainedA[i] << " <-> " << constrainedB[i]
-	   << " => " << constrainedDsqr[i] << "\n";
     }
 
     
@@ -259,7 +257,7 @@ void Integrator::integrate( void ){
 
   }
 
-  dumpOut->writeFinal();
+  dumpOut->writeFinal(currTime);
 
   delete dumpOut;
   delete statOut;
@@ -296,19 +294,19 @@ void Integrator::moveA( void ){
   double ji[3];
   double angle;
 
+
+
   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];
-
-   
+    for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j];
+    
     if( atoms[i]->isDirectional() ){
 
       dAtom = (DirectionalAtom *)atoms[i];
@@ -413,8 +411,8 @@ void Integrator::constrainA(){
 
   int i,j,k;
   int done;
-  double pxab, pyab, pzab;
-  double rxab, ryab, rzab;
+  double pab[3];
+  double rab[3];
   int a, b, ax, ay, az, bx, by, bz;
   double rma, rmb;
   double dx, dy, dz;
@@ -431,8 +429,7 @@ void Integrator::constrainA(){
     moving[i] = 0;
     moved[i]  = 1;
   }
-  
-  
+
   iteration = 0;
   done = 0;
   while( !done && (iteration < maxIteration )){
@@ -451,47 +448,36 @@ void Integrator::constrainA(){
       by = (b*3) + 1;
       bz = (b*3) + 2;
 
-
       if( moved[a] || moved[b] ){
 	
-	pxab = pos[ax] - pos[bx];
-	pyab = pos[ay] - pos[by];
-	pzab = pos[az] - pos[bz];
+	pab[0] = pos[ax] - pos[bx];
+	pab[1] = pos[ay] - pos[by];
+	pab[2] = pos[az] - pos[bz];
 
-	//periodic boundary condition
-	pxab = pxab - info->box_x * copysign(1, pxab) 
-	  * (int)( fabs(pxab / info->box_x) + 0.5);
-	pyab = pyab - info->box_y * copysign(1, pyab) 
-	  * (int)( fabs(pyab / info->box_y) + 0.5);
-	pzab = pzab - info->box_z * copysign(1, pzab) 
-	  * (int)( fabs(pzab / info->box_z) + 0.5);
-      
-	pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
+ 	//periodic boundary condition
+
+	info->wrapVector( pab );
+
+ 	pabsq = pab[0] * pab[0] + pab[1] * pab[1] + pab[2] * pab[2];
+
 	rabsq = constrainedDsqr[i];
-	diffsq = pabsq - rabsq;
+	diffsq = rabsq - pabsq;
 
 	// the original rattle code from alan tidesley
 	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)( fabs(rxab / info->box_x) + 0.5);
- 	  ryab = ryab - info->box_y * copysign(1, ryab) 
-	    * (int)( fabs(ryab / info->box_y) + 0.5);
-	  rzab = rzab - info->box_z * copysign(1, rzab) 
-	    * (int)( fabs(rzab / info->box_z) + 0.5);
+	  rab[0] = oldPos[ax] - oldPos[bx];
+	  rab[1] = oldPos[ay] - oldPos[by];
+	  rab[2] = oldPos[az] - oldPos[bz];
 
-	  rpab = rxab * pxab + ryab * pyab + rzab * pzab;
+	  info->wrapVector( rab );
+
+	  rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
+
 	  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();
@@ -505,12 +491,13 @@ void Integrator::constrainA(){
 
 	  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;
 
+          dx = rab[0] * gab;
+          dy = rab[1] * gab;
+          dz = rab[2] * gab;
+
 	  pos[ax] += rma * dx;
 	  pos[ay] += rma * dy;
 	  pos[az] += rma * dz;
@@ -545,7 +532,6 @@ void Integrator::constrainA(){
     }
 
     iteration++;
-    cerr << "iterainA = " << iteration << "\n";
   }
 
   if( !done ){
@@ -564,7 +550,7 @@ void Integrator::constrainB( void ){
   int i,j,k;
   int done;
   double vxab, vyab, vzab;
-  double rxab, ryab, rzab;
+  double rab[3];
   int a, b, ax, ay, az, bx, by, bz;
   double rma, rmb;
   double dx, dy, dz;
@@ -582,18 +568,20 @@ void Integrator::constrainB( void ){
   iteration = 0;
   while( !done && (iteration < maxIteration ) ){
 
+    done = 1;
+
     for(i=0; i<nConstrained; i++){
       
       a = constrainedA[i];
       b = constrainedB[i];
 
-      ax = 3*a +0;
-      ay = 3*a +1;
-      az = 3*a +2;
+      ax = (a*3) + 0;
+      ay = (a*3) + 1;
+      az = (a*3) + 2;
 
-      bx = 3*b +0;
-      by = 3*b +1;
-      bz = 3*b +2;
+      bx = (b*3) + 0;
+      by = (b*3) + 1;
+      bz = (b*3) + 2;
 
       if( moved[a] || moved[b] ){
 	
@@ -601,29 +589,24 @@ void Integrator::constrainB( void ){
 	vyab = vel[ay] - vel[by];
 	vzab = vel[az] - vel[bz];
 
-	rxab = pos[ax] - pos[bx];
-	ryab = pos[ay] - pos[by];
-	rzab = pos[az] - pos[bz];
+	rab[0] = pos[ax] - pos[bx];
+	rab[1] = pos[ay] - pos[by];
+	rab[2] = pos[az] - pos[bz];
 	
-	rxab = rxab - info->box_x * copysign(1, rxab) 
-	  * (int)( fabs(rxab / info->box_x) + 0.5);
-	ryab = ryab - info->box_y * copysign(1, ryab) 
-	  * (int)( fabs(ryab / info->box_y) + 0.5);
-	rzab = rzab - info->box_z * copysign(1, rzab) 
-	  * (int)( fabs(rzab / info->box_z) + 0.5);
-
+	info->wrapVector( rab );
+	
 	rma = 1.0 / atoms[a]->getMass();
 	rmb = 1.0 / atoms[b]->getMass();
 
-	rvab = rxab * vxab + ryab * vyab + rzab * vzab;
+	rvab = rab[0] * vxab + rab[1] * vyab + rab[2] * vzab;
 	  
 	gab = -rvab / ( ( rma + rmb ) * constrainedDsqr[i] );
 
 	if (fabs(gab) > tol) {
 	  
-	  dx = rxab * gab;
-	  dy = ryab * gab;
-	  dz = rzab * gab;
+	  dx = rab[0] * gab;
+	  dy = rab[1] * gab;
+	  dz = rab[2] * gab;
 	  
 	  vel[ax] += rma * dx;
 	  vel[ay] += rma * dy;