--- trunk/mdtools/md_code/Thermo.cpp	2002/12/29 19:11:05	218
+++ trunk/mdtools/md_code/Thermo.cpp	2003/01/03 22:04:50	223
@@ -6,7 +6,18 @@
 #include "LRI.hpp"
 #include "Integrator.hpp"
 
+#define BASE_SEED 123456789
 
+Thermo::Thermo( SimInfo* the_entry_plug ) { 
+  entry_plug = the_entry_plug;
+  baseSeed = BASE_SEED;
+  gaussStream = new gaussianSPRNG( baseSeed );
+}
+
+Thermo::~Thermo(){
+  delete gaussStream;
+}
+
 double Thermo::getKinetic(){
 
   const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
@@ -159,7 +170,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);
@@ -171,7 +194,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 * gaussStream->getGaussian();
+    vy = vbar * gaussStream->getGaussian();
+    vz = vbar * gaussStream->getGaussian();
+#endif
+
     atoms[vr]->set_vx( vx ); 
     atoms[vr]->set_vy( vy );
     atoms[vr]->set_vz( vz );
@@ -219,6 +249,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();
@@ -234,6 +269,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 * gaussStream->getGaussian();
+
+	vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
+	jy = vbar * gaussStream->getGaussian();
+
+	vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
+	jz = vbar * gaussStream->getGaussian();
+#endif
 	
 	dAtom->setJx( jx );
 	dAtom->setJy( jy );