--- trunk/OOPSE/libmdtools/NPTxyz.cpp	2003/10/31 18:28:52	847
+++ trunk/OOPSE/libmdtools/NPTxyz.cpp	2004/04/12 20:32:20	1097
@@ -1,4 +1,5 @@
 #include <math.h>
+#include "MatVec3.h"
 #include "Atom.hpp"
 #include "SRI.hpp"
 #include "AbstractClasses.hpp"
@@ -39,23 +40,26 @@ template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theI
     }
   }
 
-  // retrieve eta array from simInfo if it exists
-  data = info->getProperty(ETAVALUE_ID);
-  if(data){
-    etaValue = dynamic_cast<DoubleArrayData*>(data);
 
-    if(etaValue){
-      etaArray = etaValue->getData();
+  if( theInfo->useInitXSstate ){
 
-      for(i = 0; i < 3; i++){
-        for (j = 0; j < 3; j++){
-          eta[i][j] = etaArray[3*i+j];
-          oldEta[i][j] = eta[i][j];
-        }
+    // retrieve eta array from simInfo if it exists
+    data = info->getProperty(ETAVALUE_ID);
+    if(data){
+      etaValue = dynamic_cast<DoubleArrayData*>(data);
+      
+      if(etaValue){
+	etaArray = etaValue->getData();
+	
+	for(i = 0; i < 3; i++){
+	  for (j = 0; j < 3; j++){
+	    eta[i][j] = etaArray[3*i+j];
+	    oldEta[i][j] = eta[i][j];
+	  }
+	}
       }
     }
   }
-
 }
 
 template<typename T> NPTxyz<T>::~NPTxyz() {
@@ -113,9 +117,8 @@ template<typename T> void NPTxyz<T>::evolveEtaB() {
   }
 }
 
-template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
+template<typename T> void NPTxyz<T>::calcVelScale(void) {
   int i,j;
-  double vScale[3][3];
 
   for (i = 0; i < 3; i++ ) {
     for (j = 0; j < 3; j++ ) {
@@ -126,29 +129,20 @@ template<typename T> void NPTxyz<T>::getVelScaleA(doub
       }
     }
   }
+}
 
-  info->matVecMul3( vScale, vel, sc );
+template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
+  matVecMul3( vScale, vel, sc );
 }
 
 template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
-  int i,j;
+  int j;
   double myVel[3];
-  double vScale[3][3];
 
-  for (i = 0; i < 3; i++ ) {
-    for (j = 0; j < 3; j++ ) {
-      vScale[i][j] = eta[i][j];
-
-      if (i == j) {
-        vScale[i][j] += chi;
-      }
-    }
-  }
-
   for (j = 0; j < 3; j++)
     myVel[j] = oldVel[3*index + j];
 
-  info->matVecMul3( vScale, myVel, sc );
+  matVecMul3( vScale, myVel, sc );
 }
 
 template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
@@ -159,7 +153,7 @@ template<typename T> void NPTxyz<T>::getPosScale(doubl
   for(j=0; j<3; j++)
     rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
 
-  info->matVecMul3( eta, rj, sc );
+  matVecMul3( eta, rj, sc );
 }
 
 template<typename T> void NPTxyz<T>::scaleSimBox( void ){
@@ -240,7 +234,7 @@ template<typename T> void NPTxyz<T>::scaleSimBox( void
     simError();
   } else {
     info->getBoxM(hm);
-    info->matMul3(hm, scaleMat, hmnew);
+    matMul3(hm, scaleMat, hmnew);
     info->setBoxM(hmnew);
   }
 }
@@ -276,9 +270,9 @@ template<typename T> double NPTxyz<T>::getConservedQua
 
   thermostat_potential = fkBT* integralOfChidt / eConvert;
 
-  info->transposeMat3(eta, a);
-  info->matMul3(a, eta, b);
-  trEta = info->matTrace3(b);
+  transposeMat3(eta, a);
+  matMul3(a, eta, b);
+  trEta = matTrace3(b);
 
   barostat_kinetic = NkBT * tb2 * trEta /
     (2.0 * eConvert);
@@ -305,11 +299,11 @@ template<typename T> string NPTxyz<T>::getAdditionalPa
   const int BUFFERSIZE = 2000; // size of the read buffer
   char buffer[BUFFERSIZE];
 
-  sprintf(buffer,"\t%lf\t%lf;", chi, integralOfChidt);
+  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
   parameters += buffer;
 
   for(int i = 0; i < 3; i++){
-    sprintf(buffer,"\t%lf\t%lf\t%lf;", eta[3*i], eta[3*i+1], eta[3*i+2]);
+    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
     parameters += buffer;
   }