--- trunk/OOPSE/libmdtools/Integrator.cpp	2003/07/11 22:34:48	594
+++ trunk/OOPSE/libmdtools/Integrator.cpp	2003/07/14 21:28:54	597
@@ -224,8 +224,6 @@ void Integrator::integrate( void ){
   pos  = Atom::getPosArray();
   vel  = Atom::getVelArray();
   frc  = Atom::getFrcArray();
-  trq  = Atom::getTrqArray();
-  Amat = Atom::getAmatArray();
 
   while( currTime < runTime ){
 
@@ -234,8 +232,7 @@ void Integrator::integrate( void ){
       calcStress = 1;
     }
 
-    std::cerr << "calcPot = " << calcPot << "; calcStress = "
-	      << calcStress << "\n";
+    std::cerr << currTime << "\n";
 
     integrateStep( calcPot, calcStress );
       
@@ -282,7 +279,7 @@ void Integrator::integrateStep( int calcPot, int calcS
 
   preMove();
   moveA();
-  if( nConstrained ) constrainA();
+  //if( nConstrained ) constrainA();
 
   // calc forces
 
@@ -304,7 +301,7 @@ void Integrator::moveA( void ){
   double Tb[3];
   double ji[3];
   double angle;
-  double A[3][3];
+  double A[3][3], At[3][3];
 
 
   for( i=0; i<nAtoms; i++ ){
@@ -315,20 +312,11 @@ void Integrator::moveA( void ){
     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];
@@ -338,9 +326,9 @@ void Integrator::moveA( void ){
       Tb[0] = dAtom->getTx();
       Tb[1] = dAtom->getTy();
       Tb[2] = dAtom->getTz();
-      
+
       dAtom->lab2Body( Tb );
-      
+
       // get the angular momentum, and propagate a half step
       
       ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
@@ -349,44 +337,39 @@ 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, A ); 
-      
+      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, A );
+      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, A );
+      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, A );
+      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, A );
+      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
       
       dAtom->setJx( ji[0] );
       dAtom->setJy( ji[1] );
       dAtom->setJz( ji[2] );
+
+      std::cerr << "Amat[" << i << "]\n";
+      info->printMat9( &Amat[aMatIndex] );
+	  
+      std::cerr << "ji[" << i << "]\t"
+		<< ji[0] << "\t"
+		<< ji[1] << "\t"
+		<< ji[2] << "\n";
+	  
     }
     
   }
@@ -395,24 +378,20 @@ void Integrator::moveB( void ){
 
 void Integrator::moveB( void ){
   int i,j,k;
-  int atomIndex;
+  int atomIndex, aMatIndex;
   DirectionalAtom* dAtom;
   double Tb[3];
   double ji[3];
 
   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;
 
-    std::cerr<< "MoveB vel[" << i << "] = "
-	     << vel[atomIndex] << "\t"
-	     << vel[atomIndex+1]<< "\t"
-	     << vel[atomIndex+2]<< "\n";
-
-
+ 
     if( atoms[i]->isDirectional() ){
       
       dAtom = (DirectionalAtom *)atoms[i];
@@ -423,6 +402,11 @@ void Integrator::moveB( void ){
       Tb[1] = dAtom->getTy();
       Tb[2] = dAtom->getTz();
       
+      std::cerr << "TrqB[" << i << "]\t"
+		<< Tb[0] << "\t"
+		<< Tb[1] << "\t"
+		<< Tb[2] << "\n";
+
       dAtom->lab2Body( Tb );
       
       // get the angular momentum, and complete the angular momentum
@@ -435,6 +419,15 @@ void Integrator::moveB( void ){
       dAtom->setJx( ji[0] );
       dAtom->setJy( ji[1] );
       dAtom->setJz( ji[2] );
+
+
+      std::cerr << "Amat[" << i << "]\n";
+      info->printMat9( &Amat[aMatIndex] );
+	  
+      std::cerr << "ji[" << i << "]\t"
+		<< ji[0] << "\t"
+		<< ji[1] << "\t"
+		<< ji[2] << "\n";
     }
   }
 
@@ -690,7 +683,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;
@@ -702,11 +695,12 @@ void Integrator::rotate( int axes1, int axes2, double 
   double tempA[3][3];
   double tempJ[3];
 
+
   // initialize the tempA
 
   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];
     }
   }
 
@@ -763,9 +757,9 @@ void Integrator::rotate( int axes1, int axes2, double 
 
   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];
       }
     }
   }