--- branches/new_design/OOPSE-2.0/src/integrators/NPTxyz.cpp	2004/11/03 16:08:43	1701
+++ branches/new_design/OOPSE-2.0/src/integrators/NPTxyz.cpp	2004/12/02 02:08:29	1822
@@ -1,19 +1,10 @@
-#include <math.h>
-#include "math/MatVec3.h"
-#include "primitives/Atom.hpp"
-#include "primitives/SRI.hpp"
-#include "primitives/AbstractClasses.hpp"
 #include "brains/SimInfo.hpp"
-#include "UseTheForce/ForceFields.hpp"
 #include "brains/Thermo.hpp"
-#include "io/ReadWrite.hpp"
-#include "integrators/Integrator.hpp"
+#include "integrators/NPTxyz.hpp"
+#include "primitives/Molecule.hpp"
+#include "utils/OOPSEConstant.hpp"
 #include "utils/simError.h"
 
-#ifdef IS_MPI
-#include "brains/mpiSimulation.hpp"
-#endif
-
 // Basic non-isotropic thermostating and barostating via the Melchionna
 // modification of the Hoover algorithm:
 //
@@ -24,9 +15,9 @@ template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theI
 //
 //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
 
-template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theInfo, ForceFields* the_ff):
-  T( theInfo, the_ff )
-{
+namespace oopse {
+/*
+NPTxyz::NPTxyz (SimInfo* info): NPT(info) {
   GenericData* data;
   DoubleVectorGenericData * etaValue;
   int i,j;
@@ -34,8 +25,8 @@ template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theI
   for(i = 0; i < 3; i++){
     for (j = 0; j < 3; j++){
 
-      eta[i][j] = 0.0;
-      oldEta[i][j] = 0.0;
+      eta(i, j) = 0.0;
+      oldEta(i, j) = 0.0;
     }
   }
 
@@ -51,8 +42,8 @@ template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theI
 	
 	for(i = 0; i < 3; i++){
 	  for (j = 0; j < 3; j++){
-	    eta[i][j] = (*etaValue)[3*i+j];
-	    oldEta[i][j] = eta[i][j];
+	    eta(i, j) = (*etaValue)[3*i+j];
+	    oldEta(i, j) = eta(i, j);
 	  }
 	}
       }
@@ -60,107 +51,135 @@ template<typename T> NPTxyz<T>::~NPTxyz() {
   }
 }
 
-template<typename T> NPTxyz<T>::~NPTxyz() {
 
-  // empty for now
-}
 
-template<typename T> void NPTxyz<T>::resetIntegrator() {
+void NPTxyz::evolveEtaA() {
 
   int i, j;
 
-  for(i = 0; i < 3; i++)
-    for (j = 0; j < 3; j++)
-      eta[i][j] = 0.0;
-
-  T::resetIntegrator();
-}
-
-template<typename T> void NPTxyz<T>::evolveEtaA() {
-
-  int i, j;
-
-  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);
-      else
-        eta[i][j] = 0.0;
+    for(i = 0; i < 3; i ++){
+        for(j = 0; j < 3; j++){
+            if( i == j) {
+                eta(i, j) += dt2 *  instaVol *(press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
+            } else {
+                eta(i, j) = 0.0;
+            }
+        }
     }
-  }
 
-  for(i = 0; i < 3; i++)
-    for (j = 0; j < 3; j++)
-      oldEta[i][j] = eta[i][j];
+    for(i = 0; i < 3; i++) {
+        for (j = 0; j < 3; j++) {
+            oldEta(i, j) = eta(i, j);
+        }
+    }
+    
 }
 
-template<typename T> void NPTxyz<T>::evolveEtaB() {
+void NPTxyz::evolveEtaB() {
 
   int i,j;
 
   for(i = 0; i < 3; i++)
     for (j = 0; j < 3; j++)
-      prevEta[i][j] = eta[i][j];
+      prevEta(i, j) = eta(i, j);
 
   for(i = 0; i < 3; i ++){
     for(j = 0; j < 3; j++){
       if( i == j) {
-	eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
-	  (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
+	eta(i, j) = oldEta(i, j) + dt2 *  instaVol *
+	  (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
       } else {
-	eta[i][j] = 0.0;
+	eta(i, j) = 0.0;
       }
     }
   }
 }
 
-template<typename T> void NPTxyz<T>::calcVelScale(void) {
+void NPTxyz::calcVelScale(void) {
   int i,j;
 
   for (i = 0; i < 3; i++ ) {
     for (j = 0; j < 3; j++ ) {
-      vScale[i][j] = eta[i][j];
+      vScale(i, j) = eta(i, j);
 
       if (i == j) {
-        vScale[i][j] += chi;
+        vScale(i, j) += chi;
       }
     }
   }
 }
 
-template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
-  matVecMul3( vScale, vel, sc );
+
+void NPTxyz::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
+    sc = vScale * vel;
 }
 
-template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
-  int j;
-  double myVel[3];
+void NPTxyz::getVelScaleB(Vector3d& sc, int index ) {
+  sc = vScale * oldVel[index];
+}
 
-  for (j = 0; j < 3; j++)
-    myVel[j] = oldVel[3*index + j];
+void NPTxyz::getPosScale(const Vector3d& pos, const Vector3d& COM,
+			   int index, Vector3d& sc) {
+			   
+    Vector3d rj = (oldPos[index] + pos[j])/2.0 -COM;
+    sc = eta * rj;
+}
 
-  matVecMul3( vScale, myVel, sc );
+bool NPTxyz::etaConverged() {
+  int i;
+  double diffEta, sumEta;
+
+  sumEta = 0;
+  for(i = 0; i < 3; i++)
+    sumEta += pow(prevEta(i, i) - eta(i, i), 2);
+
+  diffEta = sqrt( sumEta / 3.0 );
+
+  return ( diffEta <= etaTolerance );
 }
 
-template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
-					       int index, double sc[3]){
-  int j;
-  double rj[3];
+*/
+    
+double NPTxyz::calcConservedQuantity(){
 
-  for(j=0; j<3; j++)
-    rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
+  double conservedQuantity;
+  double totalEnergy;
+  double thermostat_kinetic;
+  double thermostat_potential;
+  double barostat_kinetic;
+  double barostat_potential;
+  double trEta;
 
-  matVecMul3( eta, rj, sc );
+  totalEnergy = thermo.getTotalE();
+
+  thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
+
+  thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
+
+    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
+    trEta = tmp.trace();
+
+  barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
+
+  barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
+
+  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
+    barostat_kinetic + barostat_potential;
+
+
+  return conservedQuantity;
+
 }
 
-template<typename T> void NPTxyz<T>::scaleSimBox( void ){
+    
+void NPTxyz::scaleSimBox(){
 
   int i,j,k;
-  double scaleMat[3][3];
+  Mat3x3d scaleMat;
   double eta2ij, scaleFactor;
   double bigScale, smallScale, offDiagMax;
-  double hm[3][3], hmnew[3][3];
+  Mat3x3d hm;
+  Mat3x3d hmnew;
 
 
 
@@ -175,8 +194,8 @@ template<typename T> void NPTxyz<T>::scaleSimBox( void
 
   for(i=0; i<3; i++){
     for(j=0; j<3; j++){
-      scaleMat[i][j] = 0.0;
-      if(i==j) scaleMat[i][j] = 1.0;
+      scaleMat(i, j) = 0.0;
+      if(i==j) scaleMat(i, j) = 1.0;
     }
   }
 
@@ -184,40 +203,14 @@ template<typename T> void NPTxyz<T>::scaleSimBox( void
 
     // calculate the scaleFactors
 
-    scaleFactor = exp(dt*eta[i][i]);
+    scaleFactor = exp(dt*eta(i, i));
 
-    scaleMat[i][i] = scaleFactor;
+    scaleMat(i, i) = scaleFactor;
 
-    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
-    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
+    if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
+    if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
   }
 
-//   for(i=0; i<3; i++){
-//     for(j=0; j<3; 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;
-
-//       if (i != j)
-//         if (fabs(scaleMat[i][j]) > offDiagMax)
-//           offDiagMax = fabs(scaleMat[i][j]);
-//     }
-
-//     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
-//     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
-//   }
-
   if ((bigScale > 1.1) || (smallScale < 0.9)) {
     sprintf( painCave.errMsg,
              "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
@@ -225,86 +218,17 @@ template<typename T> void NPTxyz<T>::scaleSimBox( void
              " scaleMat = [%lf\t%lf\t%lf]\n"
              "            [%lf\t%lf\t%lf]\n"
              "            [%lf\t%lf\t%lf]\n",
-             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
-             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
-             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
+             scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
+             scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
+             scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
     painCave.isFatal = 1;
     simError();
   } else {
-    info->getBoxM(hm);
-    matMul3(hm, scaleMat, hmnew);
-    info->setBoxM(hmnew);
+
+        Mat3x3d hmat = currentSnapshot_->getHmat();
+        hmat = hmat *scaleMat;
+        currentSnapshot_->setHmat(hmat);
   }
 }
 
-template<typename T> bool NPTxyz<T>::etaConverged() {
-  int i;
-  double diffEta, sumEta;
-
-  sumEta = 0;
-  for(i = 0; i < 3; i++)
-    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
-
-  diffEta = sqrt( sumEta / 3.0 );
-
-  return ( diffEta <= etaTolerance );
 }
-
-template<typename T> double NPTxyz<T>::getConservedQuantity(void){
-
-  double conservedQuantity;
-  double totalEnergy;
-  double thermostat_kinetic;
-  double thermostat_potential;
-  double barostat_kinetic;
-  double barostat_potential;
-  double trEta;
-  double a[3][3], b[3][3];
-
-  totalEnergy = tStats->getTotalE();
-
-  thermostat_kinetic = fkBT * tt2 * chi * chi /
-    (2.0 * eConvert);
-
-  thermostat_potential = fkBT* integralOfChidt / eConvert;
-
-  transposeMat3(eta, a);
-  matMul3(a, eta, b);
-  trEta = matTrace3(b);
-
-  barostat_kinetic = NkBT * tb2 * trEta /
-    (2.0 * eConvert);
-
-  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
-    eConvert;
-
-  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
-    barostat_kinetic + barostat_potential;
-
-//   cout.width(8);
-//   cout.precision(8);
-
-//   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
-//       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
-//       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
-
-  return conservedQuantity;
-
-}
-
-template<typename T> string NPTxyz<T>::getAdditionalParameters(void){
-  string parameters;
-  const int BUFFERSIZE = 2000; // size of the read buffer
-  char buffer[BUFFERSIZE];
-
-  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
-  parameters += buffer;
-
-  for(int i = 0; i < 3; i++){
-    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
-    parameters += buffer;
-  }
-
-  return parameters;
-
-}