--- trunk/src/integrators/Velocitizer.cpp	2005/03/01 14:45:45	381
+++ trunk/src/integrators/Velocitizer.cpp	2005/03/01 23:02:33	388
@@ -43,9 +43,35 @@
 #include "math/SquareMatrix3.hpp"
 #include "primitives/Molecule.hpp"
 #include "primitives/StuntDouble.hpp"
-#include "math/MersenneTwister.hpp"
+
 namespace oopse {
 
+Velocitizer::Velocitizer(SimInfo* info) {
+
+    int seedValue;
+    Globals * simParams = info->getSimParams();
+
+#ifndef IS_MPI
+    if (simParams->haveSeed()) {
+        seedValue = simParams->getSeed();
+        randNumGen_ = new MTRand(seedValue);
+    }else {
+        randNumGen_ = new MTRand();
+    }    
+#else
+    if (simParams->haveSeed()) {
+        seedValue = simParams->getSeed();
+        randNumGen_ = new ParallelRandNumGen(seedValue);
+    }else {
+        randNumGen_ = new ParallelRandNumGen();
+    }    
+#endif 
+}
+
+Velocitizer::~Velocitizer() {
+    delete randNumGen_;
+}
+
 void Velocitizer::velocitize(double temperature) {
     Vector3d aVel;
     Vector3d aJ;
@@ -65,15 +91,8 @@ void Velocitizer::velocitize(double temperature) {
     Molecule * mol;
     StuntDouble * integrableObject;
 
-    
-#ifndef IS_MPI
-    MTRand randNumGen(info_->getSeed());
-#else
-    int nProcessors;
-    MPI_Comm_size(MPI_COMM_WORLD, &nProcessors);
-    MTRand randNumGen(info_->getSeed(), nProcessors, worldRank);
-#endif
 
+
     kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf());
 
     for( mol = info_->beginMolecule(i); mol != NULL;
@@ -91,7 +110,7 @@ void Velocitizer::velocitize(double temperature) {
             // centered on vbar
 
             for( int k = 0; k < 3; k++ ) {
-                aVel[k] = vbar * randNumGen.randNorm(0.0, 1.0);
+                aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0);
             }
 
             integrableObject->setVel(aVel);
@@ -106,13 +125,13 @@ void Velocitizer::velocitize(double temperature) {
 
                     aJ[l] = 0.0;
                     vbar = sqrt(2.0 * kebar * I(m, m));
-                    aJ[m] = vbar * randNumGen.randNorm(0.0, 1.0);
+                    aJ[m] = vbar * randNumGen_->randNorm(0.0, 1.0);
                     vbar = sqrt(2.0 * kebar * I(n, n));
-                    aJ[n] = vbar * randNumGen.randNorm(0.0, 1.0);
+                    aJ[n] = vbar * randNumGen_->randNorm(0.0, 1.0);
                 } else {
                     for( int k = 0; k < 3; k++ ) {
                         vbar = sqrt(2.0 * kebar * I(k, k));
-                        aJ[k] = vbar * randNumGen.randNorm(0.0, 1.0);
+                        aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0);
                     }
                 } // else isLinear