--- trunk/OOPSE/libmdtools/NPTf.cpp	2003/07/09 22:14:06	586
+++ trunk/OOPSE/libmdtools/NPTf.cpp	2003/07/10 17:10:56	588
@@ -22,9 +22,13 @@ NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
 NPTf::NPTf ( SimInfo *theInfo, ForceFields* the_ff):
   Integrator( theInfo, the_ff )
 {
-  int i;
+  int i, j;
   chi = 0.0;
-  for(i = 0; i < 9; i++) eta[i] = 0.0;
+
+  for(i = 0; i < 3; i++) 
+    for (j = 0; j < 3; j_++) 
+      eta[i][j] = 0.0;
+
   have_tau_thermostat = 0;
   have_tau_barostat = 0;
   have_target_temp = 0;
@@ -38,15 +42,11 @@ void NPTf::moveA() {
   DirectionalAtom* dAtom;
   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 ri[3], vi[3], sc[3];
+  double instaTemp, instaVol;
+  double tt2, tb2, eta2ij;
   double angle;
-  double press[9];
+  double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
 
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
@@ -58,57 +58,59 @@ void NPTf::moveA() {
   // first evolve chi a half step
   
   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
-  
-  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 < 3; i++ ) {
+    for (j = 0; j < 3; j++ ) {
+      if (i == j) {
+
+        eta[i][j] += dt2 * instaVol * 
+          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
+
+        vScale[i][j] = eta[i][j] + chi;
+        
+      } else {
+        
+        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
+
+        vScale[i][j] = eta[i][j];
+        
+      }
+    }
+  }
+
   for( i=0; i<nAtoms; i++ ){
     atomIndex = i * 3;
     aMatIndex = i * 9;
     
     // velocity half step
     
-    vx = vel[atomIndex];
-    vy = vel[atomIndex+1];
-    vz = vel[atomIndex+2];
+    vi[0] = vel[atomIndex];
+    vi[1] = vel[atomIndex+1];
+    vi[2] = 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;
+    info->matVecMul3( vScale, vi, sc );
     
-    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);
+    vi[0] += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - sc[0]);
+    vi[1] += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - sc[1]);
+    vi[2] += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - sc[2]);
 
-    vel[atomIndex] = vx;
-    vel[atomIndex+1] = vy;
-    vel[atomIndex+2] = vz;
+    vel[atomIndex] = vi[0]
+    vel[atomIndex+1] = vi[1];
+    vel[atomIndex+2] = vi[2];
 
     // position whole step    
 
-    rj[0] = pos[atomIndex];
-    rj[1] = pos[atomIndex+1];
-    rj[2] = pos[atomIndex+2];
+    ri[0] = pos[atomIndex];
+    ri[1] = pos[atomIndex+1];
+    ri[2] = pos[atomIndex+2];
 
-    info->wrapVector(rj);
+    info->wrapVector(ri);
 
-    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];
+    info->matVecMul3( eta, ri, sc );
 
-    pos[atomIndex] += dt * (vel[atomIndex] + scx);
-    pos[atomIndex+1] += dt * (vel[atomIndex+1] + scy);
-    pos[atomIndex+2] += dt * (vel[atomIndex+2] + scz);
+    pos[atomIndex] += dt * (vel[atomIndex] + sc[0]);
+    pos[atomIndex+1] += dt * (vel[atomIndex+1] + sc[1]);
+    pos[atomIndex+2] += dt * (vel[atomIndex+2] + sc[2]);
    
     if( atoms[i]->isDirectional() ){
 
@@ -170,54 +172,40 @@ void NPTf::moveA() {
 
   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];
+
+      // Calculate the matrix Product of the eta array (we only need
+      // the ij element right now):
+
+      eta2ij = 0.0;
+      for(k=0; k<3; k++){
+        eta2ij += eta[i][k] * eta[k][j];
       }
+      
+      scaleMat[i][j] = 0.0;
+      // identity matrix (see above):
+      if (i == j) scaleMat[i][j] = 1.0;
+      // Taylor expansion for the exponential truncated at second order:
+      scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
+
     }
-    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]);
-      }
-    }
-  }
+  info->matMul3(hm, scaleMat, hmnew);
+  info->setBoxM(hmnew);
   
-  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->setBoxM(hm);
-  
 }
 
 void NPTf::moveB( void ){
-  int i,j,k;
+  int i,j, k;
   int atomIndex;
   DirectionalAtom* dAtom;
   double Tb[3];
   double ji[3];
-  double press[9];
+  double vi[3], sc[3];
   double instaTemp, instaVol;
   double tt2, tb2;
-  double vx, vy, vz;
-  double scx, scy, scz;
-  const double p_convert = 1.63882576e8;
+  double press[3][3], vScale[3][3];
   
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
@@ -230,39 +218,43 @@ void NPTf::moveB( void ){
   
   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
   
-  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 < 3; i++ ) {
+    for (j = 0; j < 3; j++ ) {
+      if (i == j) {
 
+        eta[i][j] += dt2 * instaVol * 
+          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
+
+        vScale[i][j] = eta[i][j] + chi;
+        
+      } else {
+        
+        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
+
+        vScale[i][j] = eta[i][j];
+        
+      }
+    }
+  }
+
   for( i=0; i<nAtoms; i++ ){
     atomIndex = i * 3;
 
     // velocity half step
     
-    vx = vel[atomIndex];
-    vy = vel[atomIndex+1];
-    vz = vel[atomIndex+2];
+    vi[0] = vel[atomIndex];
+    vi[1] = vel[atomIndex+1];
+    vi[2] = 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;
+    info->matVecMul3( vScale, vi, sc );
     
-    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);
+    vi[0] += dt2 * ((frc[atomIndex]  /atoms[i]->getMass())*eConvert - sc[0]);
+    vi[1] += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - sc[1]);
+    vi[2] += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - sc[2]);
 
-    vel[atomIndex] = vx;
-    vel[atomIndex+1] = vy;
-    vel[atomIndex+2] = vz;
+    vel[atomIndex] = vi[0]
+    vel[atomIndex+1] = vi[1];
+    vel[atomIndex+2] = vi[2];
     
     if( atoms[i]->isDirectional() ){