--- trunk/OOPSE/libmdtools/NPTf.cpp	2003/09/15 16:52:02	763
+++ trunk/OOPSE/libmdtools/NPTf.cpp	2003/09/16 20:02:11	767
@@ -35,8 +35,27 @@ template<typename T> NPTf<T>::NPTf ( SimInfo *theInfo,
   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> NPTf<T>::~NPTf() {
+  delete[] oldPos;
+  delete[] oldVel;
+  delete[] oldJi;
+}
+
 template<typename T> void NPTf<T>::moveA() {
   
   int i, j, k;
@@ -53,6 +72,7 @@ template<typename T> void NPTf<T>::moveA() {
   double eta2ij;
   double press[3][3], vScale[3][3], hm[3][3], hmnew[3][3], scaleMat[3][3];
   double bigScale, smallScale, offDiagMax;
+  double COM[3];
 
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
@@ -60,64 +80,42 @@ template<typename T> void NPTf<T>::moveA() {
   instaTemp = tStats->getTemperature();
   tStats->getPressureTensor(press);
   instaVol = tStats->getVolume();
-   
-  // first evolve chi a half step
   
-  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
+  tStats->getCOM(COM);
 
+  //calculate scale factor of veloity
   for (i = 0; i < 3; i++ ) {
     for (j = 0; j < 3; j++ ) {
+      vScale[i][j] = eta[i][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];
-        
-      }
+        vScale[i][j] += chi;          
+      }               
     }
   }
-
+  
+  //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();
     
-    // velocity half step
-        
     info->matVecMul3( vScale, vel, sc );
-    
-    for (j = 0; j < 3; j++) {
+
+    for (j=0; j < 3; j++) {
+      // velocity half step  (use chi from previous step here):
       vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
-      rj[j] = pos[j];
+  
     }
 
     atoms[i]->setVel( vel );
-
-    // position whole step    
-
-    info->wrapVector(rj);
-
-    info->matVecMul3( eta, rj, sc );
-
-    for (j = 0; j < 3; j++ ) 
-      pos[j] += dt * (vel[j] + sc[j]);
-
-    atoms[i]->setPos( pos );
    
     if( atoms[i]->isDirectional() ){
 
       dAtom = (DirectionalAtom *)atoms[i];
-	  
+
       // get and convert the torque to body frame
       
       dAtom->getTrq( Tb );
@@ -158,9 +156,56 @@ template<typename T> void NPTf<T>::moveA() {
       
       dAtom->setJ( ji );
       dAtom->setA( A  );    
-    }                    
+    }    
   }
+
+  // advance chi half step
+  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
+
+  //calculate the integral of chidt
+  integralOfChidt += dt2*chi;
+
+  //advance eta half step
+  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] += dt2 * instaVol * press[i][j] / ( NkBT*tb2);
+    }
+    
+  //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]; 
+      
+      info->matVecMul3( eta, rj, sc );
+      
+      for(j = 0; j < 3; j++)
+        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
+
+      atoms[i]->setPos( pos );
+
+    }
+
+  }  
+
+ 
   // Scale the box after all the positions have been moved:
   
   // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
@@ -230,7 +275,7 @@ template<typename T> void NPTf<T>::moveB( void ){
 
 template<typename T> void NPTf<T>::moveB( void ){
 
-  int i, j;
+  int i, j, k;
   DirectionalAtom* dAtom;
   double Tb[3], ji[3];
   double vel[3], frc[3];
@@ -240,74 +285,119 @@ template<typename T> void NPTf<T>::moveB( void ){
   double tt2, tb2;
   double sc[3];
   double press[3][3], vScale[3][3];
+  double oldChi, prevChi;
+  double oldEta[3][3], preEta[3][3], diffEta;
   
   tt2 = tauThermostat * tauThermostat;
   tb2 = tauBarostat * tauBarostat;
 
-  instaTemp = tStats->getTemperature();
-  tStats->getPressureTensor(press);
-  instaVol = tStats->getVolume();
-   
-  // first evolve chi a half step
-  
-  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
-  
-  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);
+  // Set things up for the iteration:
 
-        vScale[i][j] = eta[i][j] + chi;
-        
-      } else {
-        
-        eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
+  oldChi = chi;
+  
+  for(i = 0; i < 3; i++)
+    for(j = 0; j < 3; j++)
+      oldEta[i][j] = eta[i][j];
 
-        vScale[i][j] = eta[i][j];
-        
-      }
-    }
-  }
-
   for( i=0; i<nAtoms; i++ ){
 
     atoms[i]->getVel( vel );
-    atoms[i]->getFrc( frc );
 
-    mass = atoms[i]->getMass();
-    
-    // velocity half step
-        
-    info->matVecMul3( vScale, vel, sc );
-    
-    for (j = 0; j < 3; j++) {
-      vel[j] += dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
-    }
+    for (j=0; j < 3; j++)
+      oldVel[3*i + j]  = vel[j];
 
-    atoms[i]->setVel( vel );
-    
     if( atoms[i]->isDirectional() ){
 
       dAtom = (DirectionalAtom *)atoms[i];
-	  
-      // get and convert the torque to body frame
+
+      dAtom->getJ( ji );
+
+      for (j=0; j < 3; j++)
+        oldJi[3*i + j] = ji[j];
+
+    }
+  }
+
+  // do the iteration:
+
+  instaVol = tStats->getVolume();
+  
+  for (k=0; k < 4; k++) {
+    
+    instaTemp = tStats->getTemperature();
+    tStats->getPressureTensor(press);
+
+    // evolve chi another half step using the temperature at t + dt/2
+
+    prevChi = chi;
+    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
+    
+    for(i = 0; i < 3; i++)
+      for(j = 0; j < 3; j++)
+        preEta[i][j] = eta[i][j];
+
+    //advance eta half step and calculate scale factor for velocity
+    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);
+          vScale[i][j] = eta[i][j] + chi;
+        }
+        else
+        {
+          eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
+          vScale[i][j] = eta[i][j];
+        }
+    }      
+
+    //advance velocity half step
+    for( i=0; i<nAtoms; i++ ){
+
+      atoms[i]->getFrc( frc );
+      atoms[i]->getVel(vel);
       
-      dAtom->getTrq( Tb );
-      dAtom->lab2Body( Tb );
+      mass = atoms[i]->getMass();
       
-      // get the angular momentum, and propagate a half step
+      info->matVecMul3( vScale, vel, sc );
+
+      for (j=0; j < 3; j++) {
+        // velocity half step  (use chi from previous step here):
+        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j]  / mass) * eConvert - sc[j]);
+      }
       
-      dAtom->getJ( ji );
+      atoms[i]->setVel( vel );
       
-      for (j=0; j < 3; j++) 
-        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
+      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 );
+          dAtom->setJ( ji );
+      }
+    }
 
-    }                    
+    
+    diffEta = 0;
+    for(i = 0; i < 3; i++)
+      diffEta += pow(preEta[i][i] - eta[i][i], 2);    
+    
+    if (fabs(prevChi - chi) <= chiTolerance && sqrt(diffEta / 3) <= etaTolerance) 
+      break;
   }
+
+  //calculate integral of chida
+  integralOfChidt += dt2*chi;
+
+  
 }
 
 template<typename T> void NPTf<T>::resetIntegrator() {
@@ -374,7 +464,8 @@ template<typename T> int NPTf<T>::readyCheck() {
 
   // 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;
 }
@@ -383,23 +474,39 @@ template<typename T> double NPTf<T>::getConservedQuant
 
   double conservedQuantity;
   double tb2;
-  double eta2[3][3];  
-  double trEta;
+  double trEta;  
+  double U;
+  double thermo;
+  double integral;
+  double baro;
+  double PV;
 
-  //HNVE
-  conservedQuantity = tStats->getTotalE();
+  U = tStats->getTotalE();
+  thermo = (fkBT * tauThermostat * tauThermostat * chi * chi / 2.0) / eConvert;
 
-  //HNVT
-  conservedQuantity += (info->getNDF() * kB * targetTemp * 
-    (integralOfChidt + tauThermostat * tauThermostat * chi * chi /2)) / eConvert;
+  tb2 = tauBarostat * tauBarostat;
+  trEta = info->matTrace3(eta);
+  baro = ((double)info->ndfTrans * kB * targetTemp * tb2 * trEta * trEta / 2.0) / eConvert;
 
-  //HNPT
-  tb2 = tauBarostat *tauBarostat;
+  integral = ((double)(info->ndf + 1) * kB * targetTemp * integralOfChidt) /eConvert;
 
-  trEta = info->matTrace3(eta);
+  PV = (targetPressure * tStats->getVolume() / p_convert) / eConvert;
+
+
+  cout.width(8);
+  cout.precision(8);
   
-  conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
-			3*NkBT/2 * tb2 * trEta * trEta) / eConvert;
+  cout << info->getTime() << "\t" 
+       << chi << "\t"
+       << trEta << "\t"
+       << U << "\t" 
+       << thermo << "\t"
+       << baro << "\t"
+       << integral << "\t"
+       << PV << "\t"
+       << U+thermo+integral+PV+baro << endl;
+
+  conservedQuantity = U+thermo+integral+PV+baro;
+  return conservedQuantity;
   
-  return conservedQuantity; 
 }