--- trunk/OOPSE/libmdtools/NPTf.cpp	2003/07/09 01:41:11	577
+++ trunk/OOPSE/libmdtools/NPTf.cpp	2003/07/09 02:15:29	578
@@ -9,7 +9,7 @@
 #include "simError.h" 
 
 
-// Basic isotropic thermostating and barostating via the Melchionna
+// Basic non-isotropic thermostating and barostating via the Melchionna
 // modification of the Hoover algorithm:
 //
 //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
@@ -39,6 +39,10 @@ void NPTf::moveA() {
   double Tb[3];
   double ji[3];
   double rj[3];
+  double ident[3][3], eta1[3][3], eta2[3][3], hmnew[3][3];
+  double hm[9];
+  double vx, vy, vz;
+  double scx, scy, scz;
   double instaTemp, instaPress, instaVol;
   double tt2, tb2;
   double angle;
@@ -160,51 +164,107 @@ void NPTf::moveA() {
   }
 
   // Scale the box after all the positions have been moved:
-  
 
+  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
+  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
 
-  // Use a taylor expansion for eta products
+
+  for(i=0; i<3; i++){
+    for(j=0; j<3; j++){
+      ident[i][j] = 0.0;
+      eta1[i][j] = eta[3*i+j];
+      eta2[i][j] = 0.0;
+      for(k=0; k<3; k++){
+        eta2[i][j] += eta[3*i+k] * eta[3*k+j];
+      }
+    }
+    ident[i][i] = 1.0;
+  }
+
    
   info->getBoxM(hm);
+ 
+  for(i=0; i<3; i++){
+    for(j=0; j<3; j++){      
+      hmnew[i][j] = 0.0;
+      for(k=0; k<3; k++){
+        // remember that hmat has transpose ordering for Fortran compat:
+        hmnew[i][j] += hm[3*k+i] * (ident[k][j] 
+                                    + dt * eta1[k][j] 
+                                    + 0.5 * dt * dt * eta2[k][j]);
+      }
+    }
+  }
   
+  for (i = 0; i < 3; i++) {
+    for (j = 0; j < 3; j++) {
+      // remember that hmat has transpose ordering for Fortran compat:
+      hm[3*j + 1] = hmnew[i][j];
+    }
+  }
 
-
-
-
-
-   info->scaleBox(exp(dt*eta));
-
-
+  info->setBoxM(hm);
+  
 }
 
-void NPTi::moveB( void ){
+void NPTf::moveB( void ){
   int i,j,k;
   int atomIndex;
   DirectionalAtom* dAtom;
   double Tb[3];
   double ji[3];
-  double instaTemp, instaPress, instaVol;
+  double press[9];
+  double instaTemp, instaVol;
   double tt2, tb2;
+  double vx, vy, vz;
+  double scx, scy, scz;
+  const double p_convert = 1.63882576e8;
   
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
 
   instaTemp = tStats->getTemperature();
-  instaPress = tStats->getPressure();
-  instaVol = tStats->getVolume();
+  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;
-  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (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;
-    
+
     // velocity half step
-    for( j=atomIndex; j<(atomIndex+3); j++ )
-    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;
+    
     if( atoms[i]->isDirectional() ){
       
       dAtom = (DirectionalAtom *)atoms[i];