--- trunk/OOPSE/libmdtools/Integrator.cpp	2003/06/25 21:12:14	567
+++ trunk/OOPSE/libmdtools/Integrator.cpp	2003/07/02 21:26:55	572
@@ -257,7 +257,7 @@ void Integrator::integrate( void ){
 
   }
 
-  dumpOut->writeFinal();
+  dumpOut->writeFinal(currTime);
 
   delete dumpOut;
   delete statOut;
@@ -411,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;
@@ -450,37 +450,28 @@ void Integrator::constrainA(){
 
       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;
+	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)) {
-	  rxab = oldPos[ax] - oldPos[bx];
-	  ryab = oldPos[ay] - oldPos[by];
-	  rzab = oldPos[az] - oldPos[bz];
+	  rab[0] = oldPos[ax] - oldPos[bx];
+	  rab[1] = oldPos[ay] - oldPos[by];
+	  rab[2] = 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);
+	  info->wrapVector( rab );
 
-	  rpab = rxab * pxab + ryab * pyab + rzab * pzab;
+	  rpab = rab[0] * pab[0] + rab[1] * pab[1] + rab[2] * pab[2];
 
 	  rpabsq = rpab * rpab;
 
@@ -503,9 +494,9 @@ void Integrator::constrainA(){
 
 	  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;
@@ -559,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;
@@ -598,30 +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;