--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/03 20:21:54	447
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/04/04 01:57:11	454
@@ -125,10 +125,7 @@ double Thermo::getTemperature(){
   return total;
 }
 
-double Thermo::getTemperature(){
-
-  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
-  double temperature;
+int Thermo::getNDF(){
   int ndf_local, ndf;
   
   ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
@@ -141,8 +138,32 @@ double Thermo::getTemperature(){
 #endif
 
   ndf = ndf - 3;
+
+  return ndf;
+}
+
+int Thermo::getNDFraw() {
+  int ndfRaw_local, ndfRaw;
+
+  // Raw degrees of freedom that we have to set
+  ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented;
   
-  temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
+#ifdef IS_MPI
+  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+#else
+  ndfRaw = ndfRaw_local;
+#endif
+
+  return ndfRaw;
+}
+
+
+double Thermo::getTemperature(){
+
+  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
+  double temperature;
+  
+  temperature = ( 2.0 * this->getKinetic() ) / ( (double)this->getNDF() * kb );
   return temperature;
 }
 
@@ -180,24 +201,9 @@ void Thermo::velocitize() {
   n_oriented    = entry_plug->n_oriented;
   n_constraints = entry_plug->n_constraints;
   
-  // Raw degrees of freedom that we have to set
-  ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented;
-
-  // Degrees of freedom that can contain kinetic energy
-  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
-    - entry_plug->n_constraints;
+  kebar = kb * temperature * (double)this->getNDF() / 
+    ( 2.0 * (double)this->getNDFraw() );
   
-#ifdef IS_MPI
-  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
-  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
-#else
-  ndfRaw = ndfRaw_local;
-  ndf = ndf_local;
-#endif
-  ndf = ndf - 3;
-
-  kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
-  
   for(vr = 0; vr < n_atoms; vr++){
     
     // uses equipartition theory to solve for vbar in angstrom/fs
@@ -253,7 +259,7 @@ void Thermo::velocitize() {
 
 	vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
 	jy = vbar * gaussStream->getGaussian();
-
+        
 	vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
 	jz = vbar * gaussStream->getGaussian();