--- trunk/OOPSE/libmdtools/NVT.cpp	2003/06/20 20:29:36	561
+++ trunk/OOPSE/libmdtools/NVT.cpp	2003/06/24 22:51:57	565
@@ -14,10 +14,9 @@ NVT::NVT ( SimInfo *theInfo, ForceFields* the_ff):
 NVT::NVT ( SimInfo *theInfo, ForceFields* the_ff):
   Integrator( theInfo, the_ff )
 {
-  zeta = 0.0;
+  chi = 0.0;
   have_tau_thermostat = 0;
   have_target_temp = 0;
-  have_qmass = 0;
 }
 
 void NVT::moveA() {
@@ -27,12 +26,14 @@ void NVT::moveA() {
   DirectionalAtom* dAtom;
   double Tb[3];
   double ji[3];
-  double ke;
+  double instTemp;
   double angle;
 
+  instTemp = tStats->getTemperature();
 
-  ke = tStats->getKinetic() * eConvert;
-  zeta += dt2 * ( (2.0 * ke  -  NkBT) / qmass );
+  // first evolve chi a half step
+  
+  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
 
   for( i=0; i<nAtoms; i++ ){
     atomIndex = i * 3;
@@ -40,7 +41,7 @@ void NVT::moveA() {
     
     // velocity half step
     for( j=atomIndex; j<(atomIndex+3); j++ )
-      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert - vel[j]*zeta);
+      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert - vel[j]*chi);
 
     // position whole step    
     for( j=atomIndex; j<(atomIndex+3); j++ )
@@ -65,9 +66,9 @@ void NVT::moveA() {
       ji[1] = dAtom->getJy();
       ji[2] = dAtom->getJz();
       
-      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*zeta);
-      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*zeta);
-      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*zeta);
+      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
+      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
+      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
       
       // use the angular velocities to propagate the rotation matrix a
       // full time step
@@ -106,18 +107,17 @@ void NVT::moveB( void ){
   DirectionalAtom* dAtom;
   double Tb[3];
   double ji[3];
-  double ke;
+  double instTemp;
   
-
-  ke = tStats->getKinetic() * eConvert;
-  zeta += dt2 * ( (2.0 * ke  -  NkBT) / qmass );
+  instTemp = tStats->getTemperature();
+  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
   
   for( i=0; i<nAtoms; i++ ){
     atomIndex = i * 3;
     
     // velocity half step
     for( j=atomIndex; j<(atomIndex+3); j++ )
-      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert - vel[j]*zeta);
+      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert - vel[j]*chi);
     
     if( atoms[i]->isDirectional() ){
       
@@ -138,9 +138,9 @@ void NVT::moveB( void ){
       ji[1] = dAtom->getJy();
       ji[2] = dAtom->getJz();
       
-      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*zeta);
-      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*zeta);
-      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*zeta);
+      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
+      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
+      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
       
       dAtom->setJx( ji[0] );
       dAtom->setJy( ji[1] );
@@ -162,42 +162,17 @@ int NVT::readyCheck() {
     simError();
     return -1;
   }
-     
-  // Next check to see that we have a reasonable number of degrees of freedom
-  // and then set NkBT if we do have it.   Unreasonable numbers of DOFs
-  // are also fatal.
-
-  if (info->ndf > 0) {
-    NkBT = (double)info->ndf * kB * targetTemp;
-  } else {
+  
+  // We must set tauThermostat.
+   
+  if (!have_tau_thermostat) {
     sprintf( painCave.errMsg,
-             "NVT error: We got a silly number of degrees of freedom!\n"
-             );
+             "NVT error: If you use the constant temperature\n"
+             "   integrator, you must set tauThermostat.\n");
     painCave.isFatal = 1;
     simError();
     return -1;
-  }
-    
-  // We have our choice on setting qmass or tauThermostat.  One of them
-  // must be set.
-
-  if (!have_qmass) {
-    if (have_tau_thermostat) {
-      sprintf( painCave.errMsg,
-               "NVT info: Setting qMass = %lf\n", tauThermostat * NkBT);
-      this->setQmass(tauThermostat * NkBT);      
-      painCave.isFatal = 0;
-      simError();
-    } else {
-      sprintf( painCave.errMsg,
-               "NVT error: If you use the constant temperature\n"
-               "   integrator, you must set either tauThermostat or qMass.\n");
-      painCave.isFatal = 1;
-      simError();
-      return -1;
-    }
-  }
-  
+  }    
   return 1;
 }