--- trunk/OOPSE/libmdtools/Thermo.cpp	2003/03/26 14:55:50	402
+++ trunk/OOPSE/libmdtools/Thermo.cpp	2003/03/26 15:37:05	403
@@ -152,13 +152,13 @@ void Thermo::velocitize() {
   double vx, vy, vz;
   double jx, jy, jz;
   int i, vr, vd; // velocity randomizer loop counters
-  double *vdrift;
+  double vdrift[3];
   double vbar;
   const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
   double av2;
   double kebar;
-  int ndf; // number of degrees of freedom
-  int ndfRaw; // the raw number of degrees of freedom
+  int ndf, ndf_local; // number of degrees of freedom
+  int ndfRaw, ndfRaw_local; // the raw number of degrees of freedom
   int n_atoms;
   Atom** atoms;
   DirectionalAtom* dAtom;
@@ -172,9 +172,22 @@ 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;
 
-  ndfRaw = 3 * n_atoms + 3 * n_oriented;
-  ndf = ndfRaw - n_constraints - 3;
+  // Degrees of freedom that can contain kinetic energy
+  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
+    - entry_plug->n_constraints;
+  
+#ifdef IS_MPI
+  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
+  MPI::COMM_WORLD.Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM);
+#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++){
@@ -200,7 +213,7 @@ void Thermo::velocitize() {
 
   // Get the Center of Mass drift velocity.
 
-  vdrift = getCOMVel();
+  getCOMVel(vdrift);
   
   //  Corrects for the center of mass drift.
   // sums all the momentum and divides by total mass.
@@ -244,15 +257,13 @@ double* Thermo::getCOMVel(){
   }
 }
 
-double* Thermo::getCOMVel(){
+void Thermo::getCOMVel(double vdrift[3]){
 
   double mtot, mtot_local;
-  double* vdrift;
   double vdrift_local[3];
   int vd, n_atoms;
   Atom** atoms;
 
-  vdrift = new double[3];
   // We are very careless here with the distinction between n_atoms and n_local
   // We should really fix this before someone pokes an eye out.
 
@@ -275,7 +286,7 @@ double* Thermo::getCOMVel(){
 
 #ifdef IS_MPI
   MPI::COMM_WORLD.Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM);
-  MPI::COMM_WORLD.Allreduce(&vdrift_local,&vdrift,3,MPI_DOUBLE,MPI_SUM);
+  MPI::COMM_WORLD.Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM);
 #else
   mtot = mtot_local;
   for(vd = 0; vd < 3; vd++) {
@@ -287,6 +298,5 @@ double* Thermo::getCOMVel(){
     vdrift[vd] = vdrift[vd] / mtot;
   }
   
-  return vdrift;
 }