--- trunk/OOPSE/libmdtools/Integrator.cpp	2003/06/25 21:11:42	566
+++ trunk/OOPSE/libmdtools/Integrator.cpp	2003/06/25 21:12:14	567
@@ -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";
     }
 
     
@@ -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];
@@ -431,8 +429,7 @@ void Integrator::constrainA(){
     moving[i] = 0;
     moved[i]  = 1;
   }
-  
-  
+
   iteration = 0;
   done = 0;
   while( !done && (iteration < maxIteration )){
@@ -451,31 +448,31 @@ 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];
 
-	//periodic boundary condition
+ 	//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;
+
+ 	pabsq = pxab * pxab + pyab * pyab + pzab * pzab;
+
 	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) 
@@ -484,14 +481,12 @@ void Integrator::constrainA(){
 	    * (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();
@@ -505,8 +500,9 @@ 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;
@@ -545,7 +541,6 @@ void Integrator::constrainA(){
     }
 
     iteration++;
-    cerr << "iterainA = " << iteration << "\n";
   }
 
   if( !done ){
@@ -582,18 +577,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] ){
 	
@@ -605,13 +602,14 @@ void Integrator::constrainB( void ){
 	ryab = pos[ay] - pos[by];
 	rzab = 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);
-
+	
 	rma = 1.0 / atoms[a]->getMass();
 	rmb = 1.0 / atoms[b]->getMass();