--- trunk/OOPSE/libmdtools/NPTf.cpp	2003/07/09 01:41:11	577
+++ trunk/OOPSE/libmdtools/NPTf.cpp	2003/07/09 22:14:06	586
@@ -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,35 +39,38 @@ 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;
   double press[9];
-  const double p_convert = 1.63882576e8;
 
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
 
   instaTemp = tStats->getTemperature();
   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[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2);
+  eta[0] += dt2 * instaVol * (press[0] - targetPressure/p_convert) / 
+    (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[4] += dt2 * instaVol * (press[4] - targetPressure/p_convert) / 
+    (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);
+  eta[8] += dt2 * instaVol * (press[8] - targetPressure/p_convert) / 
+    (NkBT*tb2);
   
   for( i=0; i<nAtoms; i++ ){
     atomIndex = i * 3;
@@ -160,51 +163,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 + i] = 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();
+  tStats->getPressureTensor(press);
   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/p_convert) / 
+    (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/p_convert) / 
+    (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/p_convert) / 
+    (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];
@@ -235,14 +294,14 @@ void NPTi::moveB( void ){
   }
 }
 
-int NPTi::readyCheck() {
+int NPTf::readyCheck() {
  
   // First check to see if we have a target temperature. 
   // Not having one is fatal. 
   
   if (!have_target_temp) {
     sprintf( painCave.errMsg,
-             "NPTi error: You can't use the NPTi integrator\n"
+             "NPTf error: You can't use the NPTf integrator\n"
              "   without a targetTemp!\n"
              );
     painCave.isFatal = 1;
@@ -252,7 +311,7 @@ int NPTi::readyCheck() {
 
   if (!have_target_pressure) {
     sprintf( painCave.errMsg,
-             "NPTi error: You can't use the NPTi integrator\n"
+             "NPTf error: You can't use the NPTf integrator\n"
              "   without a targetPressure!\n"
              );
     painCave.isFatal = 1;
@@ -264,7 +323,7 @@ int NPTi::readyCheck() {
    
   if (!have_tau_thermostat) {
     sprintf( painCave.errMsg,
-             "NPTi error: If you use the NPTi\n"
+             "NPTf error: If you use the NPTf\n"
              "   integrator, you must set tauThermostat.\n");
     painCave.isFatal = 1;
     simError();
@@ -275,7 +334,7 @@ int NPTi::readyCheck() {
    
   if (!have_tau_barostat) {
     sprintf( painCave.errMsg,
-             "NPTi error: If you use the NPTi\n"
+             "NPTf error: If you use the NPTf\n"
              "   integrator, you must set tauBarostat.\n");
     painCave.isFatal = 1;
     simError();