--- trunk/OOPSE/libmdtools/NPTf.cpp	2003/07/08 21:10:16	576
+++ trunk/OOPSE/libmdtools/NPTf.cpp	2003/07/09 01:41:11	577
@@ -19,7 +19,7 @@ NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
 // 
 //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
 
-NPTi::NPTi ( SimInfo *theInfo, ForceFields* the_ff):
+NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
   Integrator( theInfo, the_ff )
 {
   int i;
@@ -31,7 +31,7 @@ void NPTi::moveA() {
   have_target_pressure = 0;
 }
 
-void NPTi::moveA() {
+void NPTf::moveA() {
   
   int i,j,k;
   int atomIndex, aMatIndex;
@@ -42,49 +42,70 @@ void NPTi::moveA() {
   double instaTemp, instaPress, instaVol;
   double tt2, tb2;
   double angle;
+  double press[9];
+  const double p_convert = 1.63882576e8;
 
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
 
   instaTemp = tStats->getTemperature();
-  instaPress = tStats->getPressure();
+  tStats->getPressureTensor(press);
+
+  for (i=0; i < 9; i++) press[i] *= p_convert;
+
   instaVol = tStats->getVolume();
    
   // first evolve chi a half step
   
   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
   
-  for (i = 0; i < 9; i++) {
-    eta[i] += dt2 * ( instaVol * (sigma[i] - targetPressure*identMat[i])) 
-      / (NkBT*tb2));
-}
-
+  eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
+  eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2);
+  eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2);
+  eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2);
+  eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2);
+  eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2);
+  eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2);
+  eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2);
+  eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2);
+  
   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 
-                       - vel[j]*(chi+eta));
+    
+    vx = vel[atomIndex];
+    vy = vel[atomIndex+1];
+    vz = vel[atomIndex+2];
+    
+    scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz;
+    scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz;
+    scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz;
+    
+    vx += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - scx);
+    vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy);
+    vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz);
+
+    vel[atomIndex] = vx;
+    vel[atomIndex+1] = vy;
+    vel[atomIndex+2] = vz;
 
     // position whole step    
 
-    for( j=atomIndex; j<(atomIndex+3); j=j+3 ) {
-      rj[0] = pos[j];
-      rj[1] = pos[j+1];
-      rj[2] = pos[j+2];
+    rj[0] = pos[atomIndex];
+    rj[1] = pos[atomIndex+1];
+    rj[2] = pos[atomIndex+2];
 
-      info->wrapVector(rj);
+    info->wrapVector(rj);
 
-      pos[j] += dt * (vel[j] + eta*rj[0]);
-      pos[j+1] += dt * (vel[j+1] + eta*rj[1]);
-      pos[j+2] += dt * (vel[j+2] + eta*rj[2]);
-    }
+    scx = eta[0]*rj[0] + eta[1]*rj[1] + eta[2]*rj[2];
+    scy = eta[3]*rj[0] + eta[4]*rj[1] + eta[5]*rj[2];
+    scz = eta[6]*rj[0] + eta[7]*rj[1] + eta[8]*rj[2];
 
-    // Scale the box after all the positions have been moved:
-
-    info->scaleBox(exp(dt*eta));
+    pos[atomIndex] += dt * (vel[atomIndex] + scx);
+    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
+    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
    
     if( atoms[i]->isDirectional() ){
 
@@ -137,6 +158,23 @@ void NPTi::moveA() {
     }
     
   }
+
+  // Scale the box after all the positions have been moved:
+  
+
+
+  // Use a taylor expansion for eta products
+   
+  info->getBoxM(hm);
+  
+
+
+
+
+
+   info->scaleBox(exp(dt*eta));
+
+
 }
 
 void NPTi::moveB( void ){