--- trunk/mdtools/md_code/Thermo.cpp	2002/07/09 21:46:32	25
+++ trunk/mdtools/md_code/Thermo.cpp	2003/01/02 20:14:08	221
@@ -1,4 +1,5 @@
 #include <cmath>
+#include <mpi++.h>
 
 #include "Thermo.hpp"
 #include "SRI.hpp"
@@ -17,12 +18,15 @@ double Thermo::getKinetic(){
   DirectionalAtom *dAtom;
 
   int n_atoms;
+  double kinetic_global;
   Atom** atoms;
+
   
   n_atoms = entry_plug->n_atoms;
   atoms = entry_plug->atoms;
 
   kinetic = 0.0;
+  kinetic_global = 0.0;
   for( kl=0; kl < n_atoms; kl++ ){
 
     vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
@@ -44,7 +48,11 @@ double Thermo::getKinetic(){
 	+ (jz2 / dAtom->getIzz());
     }
   }
-  
+#ifdef IS_MPI
+  MPI_COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
+  kinetic = kinetic_global;
+#endif
+
   kinetic = kinetic * 0.5 / e_convert;
 
   return kinetic;
@@ -53,6 +61,7 @@ double Thermo::getPotential(){
 double Thermo::getPotential(){
   
   double potential;
+  double potential_global;
   int el, nSRI;
   SRI** sris;
 
@@ -60,16 +69,21 @@ double Thermo::getPotential(){
   nSRI = entry_plug->n_SRI;
 
   potential = 0.0;
-
+  potential_global = 0.0;
   potential += entry_plug->longRange->get_potential();;
 
   // std::cerr << "long range potential: " << potential << "\n";
-
   for( el=0; el<nSRI; el++ ){
     
     potential += sris[el]->get_potential();
   }
 
+  // Get total potential for entire system from MPI.
+#ifdef IS_MPI
+  MPI_COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
+  potential = potential_global;
+#endif
+
   return potential;
 }
 
@@ -95,9 +109,9 @@ double Thermo::getPressure(){
 
 double Thermo::getPressure(){
 
-  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
-  const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
-  const double conv_A_m = 1.0E-10; //convert A -> m
+//  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
+// const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
+//  const double conv_A_m = 1.0E-10; //convert A -> m
 
   return 0.0;
 }
@@ -145,7 +159,19 @@ void Thermo::velocitize() {
     
     // picks random velocities from a gaussian distribution
     // centered on vbar
-    
+#ifndef USE_SPRNG
+    /* If we are using mpi, we need to use the SPRNG random
+       generator. The non drand48 generator will just repeat
+       the same numbers for every node creating a non-gaussian
+       distribution for the simulation. drand48 is fine for the
+       single processor version of the code, but SPRNG should
+       still be preferred for consistency.
+    */
+#ifdef IS_MPI
+#error "SPRNG random number generator must be used for MPI"
+#else
+#warning "Using drand48 for random number generation"
+#endif   
     x = drand48();
     y = drand48();
     vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
@@ -157,7 +183,14 @@ void Thermo::velocitize() {
     x = drand48();
     y = drand48();
     vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
-    
+#endif
+
+#ifdef USE_SPRNG
+    vx = vbar * entry_plug->gaussStream->getGaussian();
+    vy = vbar * entry_plug->gaussStream->getGaussian();
+    vz = vbar * entry_plug->gaussStream->getGaussian();
+#endif
+
     atoms[vr]->set_vx( vx ); 
     atoms[vr]->set_vy( vy );
     atoms[vr]->set_vz( vz );
@@ -205,6 +238,11 @@ void Thermo::velocitize() {
       if( atoms[i]->isDirectional() ){
 	
 	dAtom = (DirectionalAtom *)atoms[i];
+#ifdef IS_MPI
+#error "SPRNG random number generator must be used for MPI"
+#else
+#warning "Using drand48 for random number generation"
+#endif   
 	
 	vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
 	x = drand48();
@@ -220,6 +258,17 @@ void Thermo::velocitize() {
 	x = drand48();
 	y = drand48();
 	jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
+#endif
+#ifdef USE_SPRNG
+	vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
+	jx = vbar * entry_plug->gaussStream->getGaussian();
+
+	vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
+	jy = vbar * entry_plug->gaussStream->getGaussian();
+
+	vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
+	jz = vbar * entry_plug->gaussStream->getGaussian();
+#endif
 	
 	dAtom->setJx( jx );
 	dAtom->setJy( jy );