--- trunk/OOPSE/libmdtools/NPTi.cpp	2003/09/04 21:48:35	746
+++ trunk/OOPSE/libmdtools/NPTi.cpp	2003/09/15 16:52:02	763
@@ -9,7 +9,11 @@
 #include "Integrator.hpp"
 #include "simError.h" 
 
+#ifdef IS_MPI
+#include "mpiSimulation.hpp"
+#endif
 
+
 // Basic isotropic thermostating and barostating via the Melchionna
 // modification of the Hoover algorithm:
 //
@@ -25,15 +29,149 @@ template<typename T> NPTi<T>::NPTi ( SimInfo *theInfo,
 {
   chi = 0.0;
   eta = 0.0;
+  integralOfChidt = 0.0;
   have_tau_thermostat = 0;
   have_tau_barostat = 0;
   have_target_temp = 0;
   have_target_pressure = 0;
+  have_chi_tolerance = 0;
+  have_eta_tolerance = 0;
+  have_pos_iter_tolerance = 0;
+
+  oldPos = new double[3*nAtoms];
+  oldVel = new double[3*nAtoms];
+  oldJi = new double[3*nAtoms];
+#ifdef IS_MPI
+  Nparticles = mpiSim->getTotAtoms();
+#else
+  Nparticles = theInfo->n_atoms;
+#endif
+
 }
 
+template<typename T> NPTi<T>::~NPTi() {
+  delete[] oldPos;
+  delete[] oldVel;
+  delete[] oldJi;
+}
+
 template<typename T> void NPTi<T>::moveA() {
+
+ 
+//   int i, j;
+//   DirectionalAtom* dAtom;
+//   double Tb[3], ji[3];
+//   double A[3][3], I[3][3];
+//   double angle, mass;
+//   double vel[3], pos[3], frc[3];
+
+//   double rj[3];
+//   double instaTemp, instaPress, instaVol;
+//   double tt2, tb2, scaleFactor;
+
+//   tt2 = tauThermostat * tauThermostat;
+//   tb2 = tauBarostat * tauBarostat;
+
+//   instaTemp = tStats->getTemperature();
+//   instaPress = tStats->getPressure();
+//   instaVol = tStats->getVolume();
+   
+//    // first evolve chi a half step
   
-  int i, j;
+//   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
+//   eta += dt2 * ( instaVol * (instaPress - targetPressure) / 
+// 		 (p_convert*NkBT*tb2));
+
+//   integralOfChidt += dt2* chi;
+
+//   for( i=0; i<nAtoms; i++ ){
+//     atoms[i]->getVel( vel );
+//     atoms[i]->getPos( pos );
+//     atoms[i]->getFrc( frc );
+
+//     mass = atoms[i]->getMass();
+
+//     for (j=0; j < 3; j++) {
+//       vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
+//       rj[j] = pos[j];
+//     }
+
+//     atoms[i]->setVel( vel );
+
+//     info->wrapVector(rj);
+
+//     for (j = 0; j < 3; j++) 
+//       pos[j] += dt * (vel[j] + eta*rj[j]);
+
+//     atoms[i]->setPos( pos );
+
+//     if( atoms[i]->isDirectional() ){
+
+//       dAtom = (DirectionalAtom *)atoms[i];
+	  
+//       // get and convert the torque to body frame
+      
+//       dAtom->getTrq( Tb );
+//       dAtom->lab2Body( Tb );
+      
+//       // get the angular momentum, and propagate a half step
+
+//       dAtom->getJ( ji );
+
+//       for (j=0; j < 3; j++) 
+//         ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
+      
+//       // use the angular velocities to propagate the rotation matrix a
+//       // full time step
+
+//       dAtom->getA(A);
+//       dAtom->getI(I);
+    
+//       // rotate about the x-axis      
+//       angle = dt2 * ji[0] / I[0][0];
+//       this->rotate( 1, 2, angle, ji, A ); 
+
+//       // rotate about the y-axis
+//       angle = dt2 * ji[1] / I[1][1];
+//       this->rotate( 2, 0, angle, ji, A );
+      
+//       // rotate about the z-axis
+//       angle = dt * ji[2] / I[2][2];
+//       this->rotate( 0, 1, angle, ji, A);
+      
+//       // rotate about the y-axis
+//       angle = dt2 * ji[1] / I[1][1];
+//       this->rotate( 2, 0, angle, ji, A );
+      
+//        // rotate about the x-axis
+//       angle = dt2 * ji[0] / I[0][0];
+//       this->rotate( 1, 2, angle, ji, A );
+      
+//       dAtom->setJ( ji );
+//       dAtom->setA( A  );    
+//     }                
+
+//   }
+
+//   // Scale the box after all the positions have been moved:
+  
+//   scaleFactor = exp(dt*eta);
+
+//   if ((scaleFactor > 1.1) || (scaleFactor < 0.9)) {
+//     sprintf( painCave.errMsg,
+//              "NPTi error: Attempting a Box scaling of more than 10 percent"
+//              " check your tauBarostat, as it is probably too small!\n"
+//              " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
+//              );
+//     painCave.isFatal = 1;
+//     simError();
+//   } else {         
+//     info->scaleBox(exp(dt*eta));      
+//   }
+  
+ 
+  //new version of NPTi
+  int i, j, k;
   DirectionalAtom* dAtom;
   double Tb[3], ji[3];
   double A[3][3], I[3][3];
@@ -43,6 +181,7 @@ template<typename T> void NPTi<T>::moveA() {
   double rj[3];
   double instaTemp, instaPress, instaVol;
   double tt2, tb2, scaleFactor;
+  double COM[3];
 
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
@@ -50,38 +189,29 @@ template<typename T> void NPTi<T>::moveA() {
   instaTemp = tStats->getTemperature();
   instaPress = tStats->getPressure();
   instaVol = tStats->getVolume();
-   
-   // first evolve chi a half step
   
-  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
-  eta += dt2 * ( instaVol * (instaPress - targetPressure) / 
-		 (p_convert*NkBT*tb2));
-
+  tStats->getCOM(COM);
+  
+  //evolve velocity half step
   for( i=0; i<nAtoms; i++ ){
+
     atoms[i]->getVel( vel );
-    atoms[i]->getPos( pos );
     atoms[i]->getFrc( frc );
 
     mass = atoms[i]->getMass();
 
     for (j=0; j < 3; j++) {
-      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
-      rj[j] = pos[j];
+      // velocity half step  (use chi from previous step here):
+      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi + eta));
+  
     }
 
     atoms[i]->setVel( vel );
-
-    info->wrapVector(rj);
-
-    for (j = 0; j < 3; j++) 
-      pos[j] += dt * (vel[j] + eta*rj[j]);
-
-    atoms[i]->setPos( pos );
-
+   
     if( atoms[i]->isDirectional() ){
 
       dAtom = (DirectionalAtom *)atoms[i];
-	  
+
       // get and convert the torque to body frame
       
       dAtom->getTrq( Tb );
@@ -122,10 +252,50 @@ template<typename T> void NPTi<T>::moveA() {
       
       dAtom->setJ( ji );
       dAtom->setA( A  );    
-    }                
+    }    
+  }
 
+  // evolve chi and eta  half step
+  
+  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
+  eta += dt2 * ( instaVol * (instaPress - targetPressure) / (p_convert*NkBT*tb2));
+
+  //calculate the integral of chidt
+  integralOfChidt += dt2*chi;
+
+  //save the old positions
+  for(i = 0; i < nAtoms; i++){
+    atoms[i]->getPos(pos);
+    for(j = 0; j < 3; j++)
+      oldPos[i*3 + j] = pos[j];
   }
+  
+  //the first estimation of r(t+dt) is equal to  r(t) 
+    
+  for(k = 0; k < 4; k ++){
 
+    for(i =0 ; i < nAtoms; i++){
+
+      atoms[i]->getVel(vel);
+      atoms[i]->getPos(pos);
+
+      for(j = 0; j < 3; j++)
+        rj[j] = (oldPos[i*3 + j] + pos[j])/2 - COM[j];
+
+      
+      //wrapVector(r(t)) = r(t)-R0
+      //info->wrapVector(rj);
+      
+      for(j = 0; j < 3; j++)
+        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + eta*rj[j]);
+
+      atoms[i]->setPos( pos );
+
+    }
+
+  }
+    
+
   // Scale the box after all the positions have been moved:
   
   scaleFactor = exp(dt*eta);
@@ -139,13 +309,16 @@ template<typename T> void NPTi<T>::moveA() {
     painCave.isFatal = 1;
     simError();
   } else {         
-    info->scaleBox(exp(dt*eta));      
-  }
+    info->scaleBox(scaleFactor);      
+  }  
 
+  //advance volume;
+  volume = volume * exp(dt*eta);
 }
 
 template<typename T> void NPTi<T>::moveB( void ){
 
+/*
   int i, j;
   DirectionalAtom* dAtom;
   double Tb[3], ji[3];
@@ -165,6 +338,7 @@ template<typename T> void NPTi<T>::moveB( void ){
   chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
   eta += dt2 * ( instaVol * (instaPress - targetPressure) / 
 		 (p_convert*NkBT*tb2));
+  integralOfChidt += dt2*chi;
   
   for( i=0; i<nAtoms; i++ ){
 
@@ -198,6 +372,106 @@ template<typename T> void NPTi<T>::moveB( void ){
       dAtom->setJ( ji );
     }
   }
+ 
+*/
+  
+  //new version of NPTi
+  int i, j, k;
+  DirectionalAtom* dAtom;
+  double Tb[3], ji[3];
+  double vel[3], frc[3];
+  double mass;
+
+  double instTemp, instPress, instVol;
+  double tt2, tb2;
+  double oldChi, prevChi;
+  double oldEta, preEta;
+  
+  tt2 = tauThermostat * tauThermostat;
+  tb2 = tauBarostat * tauBarostat;
+
+
+  // Set things up for the iteration:
+
+  oldChi = chi;
+  oldEta = eta;
+
+  for( i=0; i<nAtoms; i++ ){
+
+    atoms[i]->getVel( vel );
+
+    for (j=0; j < 3; j++)
+      oldVel[3*i + j]  = vel[j];
+
+    if( atoms[i]->isDirectional() ){
+
+      dAtom = (DirectionalAtom *)atoms[i];
+
+      dAtom->getJ( ji );
+
+      for (j=0; j < 3; j++)
+        oldJi[3*i + j] = ji[j];
+
+    }
+  }
+
+  // do the iteration:
+
+  instVol = tStats->getVolume();
+  
+  for (k=0; k < 4; k++) {
+    
+    instTemp = tStats->getTemperature();
+    instPress = tStats->getPressure();
+
+    // evolve chi another half step using the temperature at t + dt/2
+
+    prevChi = chi;
+    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / 
+      (tauThermostat*tauThermostat);
+
+    preEta = eta;
+    eta = oldEta + dt2 * ( instVol * (instPress - targetPressure) / 
+       (p_convert*NkBT*tb2));
+
+  
+    for( i=0; i<nAtoms; i++ ){
+
+      atoms[i]->getFrc( frc );
+      atoms[i]->getVel(vel);
+      
+      mass = atoms[i]->getMass();
+      
+      // velocity half step
+      for (j=0; j < 3; j++) 
+        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*(chi + eta));
+      
+      atoms[i]->setVel( vel );
+      
+      if( atoms[i]->isDirectional() ){
+
+        dAtom = (DirectionalAtom *)atoms[i];
+  
+        // get and convert the torque to body frame      
+  
+        dAtom->getTrq( Tb );
+        dAtom->lab2Body( Tb );       
+             
+        for (j=0; j < 3; j++) 
+          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
+      
+          dAtom->setJ( ji );
+      }
+    }
+
+    if (fabs(prevChi - chi) <= chiTolerance && fabs(preEta -eta) <= etaTolerance) 
+      break;
+  }
+
+  //calculate integral of chida
+  integralOfChidt += dt2*chi;
+
+
 }
 
 template<typename T> void NPTi<T>::resetIntegrator() {
@@ -256,9 +530,104 @@ template<typename T> int NPTi<T>::readyCheck() {
     return -1;
   }    
 
+  if (!have_chi_tolerance) {
+    sprintf( painCave.errMsg,
+             "NPTi warning: setting chi tolerance to 1e-6\n");
+    chiTolerance = 1e-6;
+    have_chi_tolerance = 1;
+    painCave.isFatal = 0;
+    simError();
+  } 
+
+    if (!have_eta_tolerance) {
+    sprintf( painCave.errMsg,
+             "NPTi warning: setting eta tolerance to 1e-6\n");
+    etaTolerance = 1e-6;
+    have_eta_tolerance = 1;
+    painCave.isFatal = 0;
+    simError();
+  } 
   // We need NkBT a lot, so just set it here:
 
-  NkBT = (double)info->ndf * kB * targetTemp;
+  NkBT = (double)Nparticles * kB * targetTemp;
+  fkBT = (double)info->ndf * kB * targetTemp;
 
   return 1;
 }
+
+template<typename T> double NPTi<T>::getConservedQuantity(void){
+
+  double conservedQuantity;
+  double tb2;
+  double eta2;  
+  double E_NPT;
+  double U;
+  double TS;
+  double PV;
+  double extra;
+
+  static double pre_U;
+  static double pre_TS;
+  static double pre_PV;
+  static double pre_extra;
+  static int hackCount = 0;
+
+  double delta_U;
+  double delta_TS;
+  double delta_PV;
+  double delta_extra;
+
+  U = tStats->getTotalE();
+
+  TS = fkBT * 
+    (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
+
+  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
+
+  tb2 = tauBarostat * tauBarostat;
+  eta2 = eta * eta;
+
+  extra = (fkBT * tb2 * eta2 / 2.0 ) / eConvert;
+  /*
+  if(hackCount == 0){
+    pre_U = U;
+    pre_TS =TS;
+    pre_PV = PV;
+    pre_extra =extra;
+    hackCount ++;
+  }
+
+  delta_U = U - pre_U;
+  delta_TS = TS - pre_TS;
+  delta_PV = PV - pre_PV;
+  delta_extra = extra - pre_extra;
+*/
+  cout.width(8);
+  cout.precision(8);
+
+  
+  cout << info->getTime() << "\t" 
+       << chi << "\t"
+       << eta << "\t"
+       << U << "\t" 
+       << TS << "\t"
+       << PV << "\t"
+       << extra << "\t"
+       << U+TS+PV+extra << endl;
+
+/*
+    pre_U = U;
+    pre_TS =TS;
+    pre_PV = PV;
+    pre_extra =extra;
+
+
+  cout << info->getTime() << "\t" 
+       << U << "\t" 
+       << U+TS << "\t"
+       << U+TS+PV << "\t"
+       << U+TS+PV+extra << endl;
+*/
+  conservedQuantity = U+TS+PV+extra;
+  return conservedQuantity; 
+}