--- trunk/OOPSE/libmdtools/Integrator.cpp	2003/07/02 21:26:55	572
+++ trunk/OOPSE/libmdtools/Integrator.cpp	2003/07/11 22:34:48	594
@@ -27,6 +27,8 @@ Integrator::Integrator( SimInfo *theInfo, ForceFields*
 
   nAtoms = info->n_atoms;
 
+  std::cerr << "integ nAtoms = "  << nAtoms << "\n";
+
   // check for constraints
   
   constrainedA    = NULL;
@@ -72,9 +74,14 @@ void Integrator::checkConstraints( void ){
     for(int j=0; j<molecules[i].getNBonds(); j++){
       
       constrained = theArray[j]->is_constrained();
+
+      std::cerr << "Is the folowing bond constrained \n";
+      theArray[j]->printMe();
       
       if(constrained){
 	
+	std::cerr << "Yes\n";
+
 	dummy_plug = theArray[j]->get_constraint();
 	temp_con[nConstrained].set_a( dummy_plug->get_a() );
 	temp_con[nConstrained].set_b( dummy_plug->get_b() );
@@ -82,7 +89,8 @@ void Integrator::checkConstraints( void ){
 	
 	nConstrained++;
 	constrained = 0;
-      }
+      } 
+      else std::cerr << "No.\n";
     }
 
     theArray = (SRI**) molecules[i].getMyBends();
@@ -226,6 +234,9 @@ void Integrator::integrate( void ){
       calcStress = 1;
     }
 
+    std::cerr << "calcPot = " << calcPot << "; calcStress = "
+	      << calcStress << "\n";
+
     integrateStep( calcPot, calcStress );
       
     currTime += dt;
@@ -293,9 +304,9 @@ void Integrator::moveA( void ){
   double Tb[3];
   double ji[3];
   double angle;
+  double A[3][3];
 
 
-
   for( i=0; i<nAtoms; i++ ){
     atomIndex = i * 3;
     aMatIndex = i * 9;
@@ -303,10 +314,21 @@ void Integrator::moveA( void ){
     // velocity half step
     for( j=atomIndex; j<(atomIndex+3); j++ )
       vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
+
+    std::cerr<< "MoveA vel[" << i << "] = "
+	     << vel[atomIndex] << "\t"
+	     << vel[atomIndex+1]<< "\t"
+	     << vel[atomIndex+2]<< "\n";
 
     // position whole step    
     for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j];
     
+
+    std::cerr<< "MoveA pos[" << i << "] = " 
+	     << pos[atomIndex] << "\t"
+	     << pos[atomIndex+1]<< "\t"
+	     << pos[atomIndex+2]<< "\n";
+
     if( atoms[i]->isDirectional() ){
 
       dAtom = (DirectionalAtom *)atoms[i];
@@ -327,26 +349,40 @@ void Integrator::moveA( void ){
       
       // use the angular velocities to propagate the rotation matrix a
       // full time step
+      
+	  // get the atom's rotation matrix
+	  
+      A[0][0] = dAtom->getAxx();
+      A[0][1] = dAtom->getAxy();
+      A[0][2] = dAtom->getAxz();
+      
+      A[1][0] = dAtom->getAyx();
+      A[1][1] = dAtom->getAyy();
+      A[1][2] = dAtom->getAyz();
       
+      A[2][0] = dAtom->getAzx();
+      A[2][1] = dAtom->getAzy();
+      A[2][2] = dAtom->getAzz();
+      
       // 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, A ); 
       
       // 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, A );
       
       // 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, A );
       
       // 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, A );
       
        // 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, A );
       
       dAtom->setJx( ji[0] );
       dAtom->setJy( ji[1] );
@@ -371,6 +407,12 @@ void Integrator::moveB( void ){
     for( j=atomIndex; j<(atomIndex+3); j++ )
       vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
 
+    std::cerr<< "MoveB vel[" << i << "] = "
+	     << vel[atomIndex] << "\t"
+	     << vel[atomIndex+1]<< "\t"
+	     << vel[atomIndex+2]<< "\n";
+
+
     if( atoms[i]->isDirectional() ){
       
       dAtom = (DirectionalAtom *)atoms[i];
@@ -422,8 +464,6 @@ void Integrator::constrainA(){
   double gab;
   int iteration;
 
-
-  
   for( i=0; i<nAtoms; i++){
     
     moving[i] = 0;
@@ -650,7 +690,7 @@ void Integrator::rotate( int axes1, int axes2, double 
 
 
 void Integrator::rotate( int axes1, int axes2, double angle, double ji[3], 
-			 double A[9] ){
+			 double A[3][3] ){
 
   int i,j,k;
   double sinAngle;
@@ -666,7 +706,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[3*i + j];
+      tempA[j][i] = A[i][j];
     }
   }
 
@@ -723,9 +763,9 @@ void Integrator::rotate( int axes1, int axes2, double 
 
   for(i=0; i<3; i++){
     for(j=0; j<3; j++){
-      A[3*j + i] = 0.0;
+      A[j][i] = 0.0;
       for(k=0; k<3; k++){
-	A[3*j + i] += tempA[i][k] * rot[j][k];
+	A[j][i] += tempA[i][k] * rot[j][k];
       }
     }
   }