--- trunk/OOPSE/libmdtools/NVT.cpp	2003/09/04 21:48:35	746
+++ trunk/OOPSE/libmdtools/NVT.cpp	2003/09/15 16:52:02	763
@@ -17,8 +17,18 @@ template<typename T> NVT<T>::NVT ( SimInfo *theInfo, F
   chi = 0.0;
   have_tau_thermostat = 0;
   have_target_temp = 0;
+  have_chi_tolerance = 0;
+  integralOfChidt = 0.0;
+
+  oldVel = new double[3*nAtoms];
+  oldJi = new double[3*nAtoms];
 }
 
+template<typename T> NVT<T>::~NVT() {
+  delete[] oldVel;
+  delete[] oldJi;
+}
+
 template<typename T> void NVT<T>::moveA() {
   
   int i, j;
@@ -30,12 +40,10 @@ template<typename T> void NVT<T>::moveA() {
 
   double instTemp;
 
-  instTemp = tStats->getTemperature();
+  // We need the temperature at time = t for the chi update below:
 
-  // first evolve chi a half step
+  instTemp = tStats->getTemperature();
   
-  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
-
   for( i=0; i<nAtoms; i++ ){
 
     atoms[i]->getVel( vel );
@@ -45,7 +53,7 @@ template<typename T> void NVT<T>::moveA() {
     mass = atoms[i]->getMass();
 
     for (j=0; j < 3; j++) {
-      // velocity half step
+      // velocity half step  (use chi from previous step here):
       vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
       // position whole step
       pos[j] += dt * vel[j];
@@ -100,58 +108,98 @@ template<typename T> void NVT<T>::moveA() {
       dAtom->setA( A  );    
     }    
   }
+
+  // Finally, evolve chi a half step (just like a velocity) using 
+  // temperature at time t, not time t+dt/2
+
+  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
+  integralOfChidt += chi*dt2;
+
 }
 
 template<typename T> void NVT<T>::moveB( void ){
-  int i, j;
+  int i, j, k;
   DirectionalAtom* dAtom;
   double Tb[3], ji[3];
   double vel[3], frc[3];
   double mass;
-
   double instTemp;
-  
-  instTemp = tStats->getTemperature();
-  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
-  
+  double oldChi, prevChi;
+
+  // Set things up for the iteration:
+
+  oldChi = chi;
+
   for( i=0; i<nAtoms; i++ ){
 
     atoms[i]->getVel( vel );
-    atoms[i]->getFrc( frc );
 
-    mass = atoms[i]->getMass();
+    for (j=0; j < 3; j++)
+      oldVel[3*i + j]  = vel[j];
 
-    // velocity half step
-    for (j=0; j < 3; j++) 
-      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
-    
-    atoms[i]->setVel( vel );
-
     if( atoms[i]->isDirectional() ){
 
       dAtom = (DirectionalAtom *)atoms[i];
 
-      // get and convert the torque to body frame      
+      dAtom->getJ( ji );
 
-      dAtom->getTrq( Tb );
-      dAtom->lab2Body( Tb );
+      for (j=0; j < 3; j++)
+        oldJi[3*i + j] = ji[j];
 
-      // get the angular momentum, and propagate a half step
+    }
+  }
 
-      dAtom->getJ( ji );
+  // do the iteration:
 
+  for (k=0; k < 4; k++) {
+    
+    instTemp = tStats->getTemperature();
+
+    // evolve chi another half step using the temperature at t + dt/2
+
+    prevChi = chi;
+    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / 
+      (tauThermostat*tauThermostat);
+  
+    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++) 
-        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
+	vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*chi);
       
-
-      dAtom->setJ( ji );
+      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) break;
   }
+  
+  integralOfChidt += dt2*chi;
 }
 
 template<typename T> void NVT<T>::resetIntegrator( void ){
   
   chi = 0.0;
+  integralOfChidt = 0.0;
 }
 
 template<typename T> int NVT<T>::readyCheck() {
@@ -182,5 +230,35 @@ template<typename T> int NVT<T>::readyCheck() {
     simError();
     return -1;
   }    
-  return 1;
+
+  if (!have_chi_tolerance) {
+    sprintf( painCave.errMsg,
+             "NVT warning: setting chi tolerance to 1e-6\n");
+    chiTolerance = 1e-6;
+    have_chi_tolerance = 1;
+    painCave.isFatal = 0;
+    simError();
+  }    
+
+  return 1;    
+
 }
+
+template<typename T> double NVT<T>::getConservedQuantity(void){
+
+  double conservedQuantity;
+  double E_NVT;
+
+  //HNVE
+  conservedQuantity = tStats->getTotalE();
+  //HNVE
+  
+  E_NVT =  (info->getNDF() * kB * targetTemp * 
+                (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0 )) / eConvert;
+
+  conservedQuantity += E_NVT;
+
+  //cerr << info->getTime() << "\t" << chi << "\t" << integralOfChidt << "\t" << E_NVT << endl;
+
+  return conservedQuantity; 
+}