Mods

platforms/brook/src/BrookCalcKineticEnergyKernel.cpp

@@ -107,6 +107,13 @@ void BrookCalcKineticEnergyKernel::initialize( const System& system ){
       _masses[ii]  =  static_cast<BrookOpenMMFloat>(system.getParticleMass(ii));
    }
 
+std::vector<double> masses;
+for( unsigned int ii = 0; ii < (unsigned int) _numberOfParticles; ii++ ){
+   masses.push_back( system.getParticleMass(ii) );
+}
+_brookVelocityCenterOfMassRemoval = new BrookVelocityCenterOfMassRemoval();
+_brookVelocityCenterOfMassRemoval->setup( masses, getPlatform() );
+
    return;
 }
 
@@ -127,19 +134,17 @@ double BrookCalcKineticEnergyKernel::execute( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   if( _masses == NULL ){
-      std::stringstream message;
-      message << methodName << " masses not set.";
-      throw OpenMMException( message.str() );
-   }    
+//BrookStreamImpl* velocities = _openMMBrookInterface.getParticleVelocities();
+//_brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( *velocities );
 
    void* dataV                            = _openMMBrookInterface.getParticleVelocities()->getData( 1 );
    float* velocity                        = (float*) dataV;
    double energy                          = 0.0;
    int index                              = 0;
 
+/*
 printf( "BrookCalcKineticEnergyKernel:\n" );
-double com[3] = { 0.0, 0.0, 0.0 };
+float com[3] = { 0.0, 0.0, 0.0 };
 for ( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){
    com[0] += velocity[index];
    com[1] += velocity[index+1];
@@ -147,13 +152,24 @@ for ( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){
    printf( "  %d %.3f [%12.5e %12.5e %12.5e]\n", ii, _masses[ii], velocity[index], velocity[index+1], velocity[index+2] );
 }
 printf( "Com [%12.5e %12.5e %12.5e]\n", com[0], com[1], com[2] );
-
 index = 0;
+float newcom[3] = { 0.0, 0.0, 0.0 };
+for ( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){
+   velocity[index]   -= com[0];
+   velocity[index+1] -= com[1];
+   velocity[index+2] -= com[2];
+   newcom[0]         += velocity[index];
+   newcom[1]         += velocity[index+1];
+   newcom[2]         += velocity[index+2];
+}
+printf( "NewCom [%12.5e %12.5e %12.5e]\n", newcom[0], newcom[1], newcom[2] );
+index = 0;
+*/
 
-   for ( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){
+   for( int ii = 0; ii < _numberOfParticles; ii++, index += 3 ){
       energy += _masses[ii]*(velocity[index]*velocity[index] + velocity[index + 1]*velocity[index + 1] + velocity[index + 2]*velocity[index + 2]);
    }
 
-printf( " Ke=%12.5e\n", 0.5*energy );
+//printf( " KQ=%12.5e\n", 0.5*energy );
    return 0.5*energy;
 }

platforms/brook/src/BrookCalcKineticEnergyKernel.h

@@ -36,6 +36,8 @@
 #include "BrookFloatStreamInternal.h"
 #include "OpenMMBrookInterface.h"
 
+#include "BrookVelocityCenterOfMassRemoval.h"
+
 namespace OpenMM {
 
 /**
@@ -97,6 +99,7 @@ class BrookCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
 
       System& _system;
 
+BrookVelocityCenterOfMassRemoval* _brookVelocityCenterOfMassRemoval;
 
 
 };

platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp

@@ -191,16 +191,18 @@ void BrookIntegrateLangevinStepKernel::initialize( const System& system, const L
    _brookShakeAlgorithm->setMaxIterations( 30 );
    _brookShakeAlgorithm->setLog( log );
 
-   // assert( (_brookShakeAlgorithm->getNumberOfConstraints() > 0) );
+   // random number generator
+
+   _brookRandomNumberGenerator   = new BrookRandomNumberGenerator( );
+   _brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() );
 
    if( printOn && log ){
-      (void) fprintf( log, "%s done shake setup const=%d\n", methodName.c_str(), numberOfConstraints );
+      (void) fprintf( log, "%s done setup:\nBrookShakeAlgorithm:\n%s\nBrookRandomNumberGenerator:\n%s\n\n", methodName.c_str(),
+                      _brookShakeAlgorithm->getContentsString().c_str(), 
+                      _brookRandomNumberGenerator->getContentsString().c_str() );
       (void) fflush( log );
    }
 
-   _brookRandomNumberGenerator   = new BrookRandomNumberGenerator( );
-   _brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() );
-
 }
 
 /** 

platforms/brook/src/BrookIntegrateVerletStepKernel.cpp

@@ -38,13 +38,16 @@ using namespace std;
 /** 
  * BrookIntegrateVerletStepKernel constructor
  * 
- * @param name        name of the stream to create
+ * @param name                  name of the kernel
  * @param platform              platform
+ * @param openMMBrookInterface  OpenMMBrookInterface reference
+ * @param system                System reference  
  *
  */
 
-BrookIntegrateVerletStepKernel::BrookIntegrateVerletStepKernel( std::string name, const Platform& platform ) :
-                                IntegrateVerletStepKernel( name, platform ){
+BrookIntegrateVerletStepKernel::BrookIntegrateVerletStepKernel( std::string name, const Platform& platform,
+                                  OpenMMBrookInterface& openMMBrookInterface, System& system ) : 
+                                  IntegrateVerletStepKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -54,6 +57,14 @@ BrookIntegrateVerletStepKernel::BrookIntegrateVerletStepKernel( std::string name
    
    _brookVerletDynamics   = NULL;
    _brookShakeAlgorithm   = NULL;
+
+   const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
+   if( brookPlatform.getLog() != NULL ){
+      setLog( brookPlatform.getLog() );
+   } else {
+      _log                = NULL;
+   }
+
 }
 
 /** 
@@ -74,45 +85,107 @@ BrookIntegrateVerletStepKernel::~BrookIntegrateVerletStepKernel( ){
    
 }
 
+/** 
+ * Get log file reference
+ * 
+ * @return  log file reference
+ *
+ */
+
+FILE* BrookIntegrateVerletStepKernel::getLog( void ) const {
+   return _log;
+}
+
+/** 
+ * Set log file reference
+ * 
+ * @param  log file reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int BrookIntegrateVerletStepKernel::setLog( FILE* log ){
+   _log = log;
+   return DefaultReturnValue;
+}
+
 /** 
  * Initialize the kernel, setting up all parameters related to integrator.
  * 
- * @param masses             the mass of each particle
- * @param constraintIndices  each element contains the indices of two particles whose distance should be constrained
- * @param constraintLengths  the required distance between each pair of constrained particles
+ * @param system                System reference  
+ * @param integrator            VerletIntegrator reference
  *
  */
 
-void BrookIntegrateVerletStepKernel::initialize( const vector<double>& masses,
-                                                 const vector<vector<int> >& constraintIndices,
-                                                 const vector<double>& constraintLengths ){
+void BrookIntegrateVerletStepKernel::initialize(  const System& system, const VerletIntegrator& integrator  ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookIntegrateVerletStepKernel::initialize";
+   int printOn                              = 1;
+   static const std::string methodName      = "BrookIntegrateVerletStepKernel::initialize";
 
 // ---------------------------------------------------------------------------------------
    
+   FILE* log             = getLog();
+   int numberOfParticles = system.getNumParticles();
+
+   // masses
+
+   std::vector<double> masses;
+   masses.resize( numberOfParticles );
+
+   for( int ii = 0; ii < numberOfParticles; ii++ ){
+      masses[ii] = static_cast<double>(system.getParticleMass(ii));
+   }
+
+   // constraints
+
+   int numberOfConstraints = system.getNumConstraints();
+
+   std::vector<std::vector<int> > constraintIndicesVector;
+   constraintIndicesVector.resize( numberOfConstraints );
+   std::vector<double> constraintLengths;
+   for( int ii = 0; ii < numberOfConstraints; ii++ ){
+
+      int particle1, particle2;
+      double distance;
+
+      system.getConstraintParameters( ii, particle1, particle2, distance );
+
+      constraintIndicesVector[ii].push_back( particle1 );
+      constraintIndicesVector[ii].push_back( particle2 );
+      constraintLengths.push_back( distance );
+   }
+
    _brookVerletDynamics          = new BrookVerletDynamics( );
    _brookVerletDynamics->setup( masses, getPlatform() );
+   _brookVerletDynamics->setLog( log );
 
    _brookShakeAlgorithm          = new BrookShakeAlgorithm( );
-   _brookShakeAlgorithm->setup( masses, constraintIndices, constraintLengths, getPlatform() );
+   _brookShakeAlgorithm->setup( masses, constraintIndicesVector, constraintLengths, getPlatform() );
+
+   BrookOpenMMFloat tolerance = static_cast<BrookOpenMMFloat>( integrator.getConstraintTolerance() );
+   _brookShakeAlgorithm->setShakeTolerance( tolerance );
+   _brookShakeAlgorithm->setMaxIterations( 30 );
+   _brookShakeAlgorithm->setLog( log );
+
+   if( printOn && log ){
+      (void) fprintf( log, "%s done w/ setup: particles=%d const=%d\n", methodName.c_str(), numberOfParticles, numberOfConstraints );
+      (void) fflush( log );
+   }
 
 }
 
 /** 
  * Execute kernel
  * 
- * @param positions          particle coordinates
- * @param velocities         particle velocities
- * @param forces             particle forces
- * @param stepSize           integration step size
+ * @param context            OpenMMContextImpl reference
+ * @param integrator         VerletIntegrator reference
  *
  */
 
-void BrookIntegrateVerletStepKernel::execute( Stream& positions, Stream& velocities,
-                                              const Stream& forces, double stepSize ){
+void BrookIntegrateVerletStepKernel::execute( OpenMMContextImpl& context, const VerletIntegrator& integrator ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -121,19 +194,19 @@ void BrookIntegrateVerletStepKernel::execute( Stream& positions, Stream& velocit
 
 // ---------------------------------------------------------------------------------------
 
-   // first time through initialize _brookVerletDynamics
-
    // for each subsequent call, check if parameters need to be updated due to a change
    // in the step size
 
    // take step
 
+   double stepSize   = integrator.getStepSize();
    double difference = stepSize - (double) _brookVerletDynamics->getStepSize();
    if( fabs( difference ) > epsilon ){
-//printf( "%s calling updateParameters\n", methodName.c_str() );
       _brookVerletDynamics->updateParameters( stepSize );
    }
-   _brookVerletDynamics->update( positions, velocities, forces, *_brookShakeAlgorithm );
 
-}
+   _brookVerletDynamics->update( *(_openMMBrookInterface.getParticlePositions()), *(_openMMBrookInterface.getParticleVelocities()),
+                                 *(_openMMBrookInterface.getParticleForces()), *_brookShakeAlgorithm );
 
+
+}

platforms/brook/src/BrookIntegrateVerletStepKernel.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "kernels.h"
+#include "OpenMMBrookInterface.h"
 #include "BrookVerletDynamics.h"
 #include "BrookShakeAlgorithm.h"
 
@@ -46,6 +47,11 @@ class BrookIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
 
    public:
 
+      // return values
+
+      static const int DefaultReturnValue = 0;
+      static const int ErrorReturnValue   = -1; 
+
       /**
        * BrookIntegrateVerletStepKernel constructor
        * 
@@ -54,7 +60,7 @@ class BrookIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
        *
        */
   
-      BrookIntegrateVerletStepKernel( std::string name, const Platform& platform );
+      BrookIntegrateVerletStepKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system  );
 
       /**
        * BrookIntegrateVerletStepKernel destructor
@@ -66,30 +72,68 @@ class BrookIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
       /** 
        * Initialize the kernel, setting up all parameters related to integrator.
        * 
-       * @param masses             the mass of each particle
-       * @param constraintIndices  each element contains the indices of two particles whose distance should be constrained
-       * @param constraintLengths  the required distance between each pair of constrained particles
-       *
+       * @param system             System reference
+       * @param integrator         VerletIntegrator reference
        */
-      void initialize( const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
-                       const std::vector<double>& constraintLengths );
+
+      void initialize( const System& system, const VerletIntegrator& integrator );
+
       /** 
        * Execute kernel
        * 
-       * @param positions          particle coordinates
-       * @param velocities         particle velocities
-       * @param forces             particle forces
-       * @param stepSize           integration step size
+       * @param context            OpenMMContextImpl reference
+       * @param integrator         VerletIntegrator reference
+       *
+       */
+
+      void execute( OpenMMContextImpl& context, const VerletIntegrator& integrator );
+
+      /** 
+       * Set log file reference
+       * 
+       * @param  log file reference
+       *
+       * @return DefaultReturnValue
+       *
+       */
+    
+      int setLog( FILE* log );
+
+      /*  
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+    
+      std::string getContents( int level ) const;
+
+      /** 
+       * Get log file reference
+       * 
+       * @return  log file reference
        *
        */
     
-      void execute( Stream& positions, Stream& velocities, const Stream& forces, double stepSize );
+      FILE* getLog( void ) const;
 
-   protected:
+   private:
+
+      FILE* _log;
 
       BrookVerletDynamics*      _brookVerletDynamics;
       BrookShakeAlgorithm*      _brookShakeAlgorithm;
  
+      // interface
+
+      OpenMMBrookInterface& _openMMBrookInterface;
+
+      // System reference
+
+      System& _system;
+
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookKernelFactory.cpp

@@ -102,7 +102,7 @@ KernelImpl* BrookKernelFactory::createKernelImpl( std::string name, const Platfo
 
 	} else if( name == IntegrateVerletStepKernel::Name() ){
 
-//      return new BrookIntegrateVerletStepKernel( name, platform, openMMBrookInterface );
+      return new BrookIntegrateVerletStepKernel( name, platform, openMMBrookInterface, context.getSystem() );
 
    // Brownian integrator
 

platforms/brook/src/BrookRandomNumberGenerator.cpp

@@ -31,6 +31,7 @@
 
 #include <sstream>
 #include "BrookRandomNumberGenerator.h"
+#include "../../reference/src/SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "OpenMMException.h"
 #include "gpu/kupdatesd.h"
 
@@ -65,7 +66,7 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
    _rvStreamIndex                   = 0;
    _rvStreamOffset                  = 0;
    _numberOfShuffles                = 0;
-   _maxShuffles                     = 3;
+   _maxShuffles                     = 0;
    //_maxShuffles                     = 100;
 
    _loadBuffer                      = NULL;
@@ -78,6 +79,8 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
    // set randomNumber seed 
 
    _randomNumberSeed      = 1393;
+   _randomNumberGenerator = Mersenne;
+   _randomNumberGenerator = Kiss;
 
    //_randomNumberSeed = randomNumberSeed ? randomNumberSeed : 1393;
    //SimTKOpenMMUtilities::setRandomNumberSeed( randomNumberSeed );
@@ -360,7 +363,7 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsKiss( void ){
 
    static unsigned int state[4];
    static int stateInitialized    = 0;
-   static const int reseed        = 5;
+   static const int reseed        = 10000;
  
    // static const char* methodName = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsKiss";
 
@@ -376,8 +379,8 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsKiss( void ){
       state[3] = rand();
 
       if( getLog() ){
-         (void) fprintf( getLog(), "LoadGVStreamsKiss: reset state seeds stateInitialized=%d reseed=%d\n",
-                         stateInitialized, reseed );
+         (void) fprintf( getLog(), "LoadGVStreamsKiss: reset state seeds stateInitialized=%d reseed=%d [%u %u %u %u]\n",
+                         stateInitialized, reseed,  state[0], state[1], state[2], state[3] );
          (void) fflush( getLog() );
       }
 
@@ -415,6 +418,48 @@ state[3] = 27587;
 	   getRandomNumberStream( jj )->loadFromArray( loadBuffer );
    }
 
+   if( getLog() ){
+      (void) fprintf( getLog(), "LoadGVStreamsKiss: stats\n%s\n", getStatisticsString().c_str() );
+      (void) fflush( getLog() );
+   }
+
+   return DefaultReturnValue;
+}
+
+
+/** 
+ * Load random number streams using Mersenne algorithm
+ * 
+ *
+ * @return DefaultReturnValue;
+ */
+
+int BrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne( void ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne";
+
+   // ---------------------------------------------------------------------------------------
+   
+   // allocate memory once for download of random nos.
+
+   float* loadBuffer = _getLoadBuffer();
+
+   for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
+      for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii += 3 ){
+         loadBuffer[ii]   = (float) SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+         loadBuffer[ii+1] = (float) SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+         loadBuffer[ii+2] = (float) SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+   	}
+	   getRandomNumberStream( jj )->loadFromArray( loadBuffer );
+   }
+
+   if( getLog() ){
+      (void) fprintf( getLog(), "_loadRandomNumberStreamsMersenne: stats\n%s\n", getStatisticsString().c_str() );
+      (void) fflush( getLog() );
+   }
+
    return DefaultReturnValue;
 }
 
@@ -631,9 +676,8 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
 
    // ---------------------------------------------------------------------------------------
 
-   static const char* methodName = "\nBrookRandomNumberGenerator::advanceGVCursor";
-
-   static const int PrintOn      = 1;
+   static const char* methodName = "BrookRandomNumberGenerator::advanceGVCursor";
+   static const int PrintOn      = 0;
 
    // ---------------------------------------------------------------------------------------
 	
@@ -668,12 +712,15 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
 
 			} else { //Need to refresh random numbers from cpu
 
-            if( UseOriginalRng ){
-               action = "loaded new values from GPU using original rng";
-               _loadGVStreamsOriginal( );
-            } else {
-               _loadRandomNumberStreamsKiss( );
+            if( _randomNumberGenerator == Mersenne ){
+               action = "loaded new values from GPU using Mersenne rng";
+               _loadRandomNumberStreamsMersenne( );
+            } else if( _randomNumberGenerator == Kiss ){
                action = "loaded new values from GPU using KISS rng";
+               _loadRandomNumberStreamsKiss( );
+            } else if( _randomNumberGenerator == Original ){
+               action = "loaded new values from GPU using original Gromac's rng";
+               _loadGVStreamsOriginal( );
             }
             _numberOfShuffles  = 0;
 			}
@@ -681,7 +728,7 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
 		}
 
       if( PrintOn && getLog() ){
-         (void) fprintf( getLog(), "%s StrmIdx=%d action=%s", methodName, _rvStreamIndex, action );
+         (void) fprintf( getLog(), "%s StrmIdx=%d action=%s\n", methodName, _rvStreamIndex, action );
          (void) fflush( getLog() );
       }
 	}
@@ -898,10 +945,13 @@ int BrookRandomNumberGenerator::setup( int numberOfParticles,  const Platform& p
 
    _initializeStreamSizes( numberOfParticles, platform );
    _initializeStreams( platform );
-   if( UseOriginalRng ){
-      _loadGVStreamsOriginal( );
-   } else {
+
+   if( _randomNumberGenerator == Mersenne ){
+      _loadRandomNumberStreamsMersenne( );
+   } else if( _randomNumberGenerator == Kiss ){
       _loadRandomNumberStreamsKiss( );
+   } else if( _randomNumberGenerator == Original ){
+      _loadGVStreamsOriginal( );
    }
 
    return DefaultReturnValue;
@@ -938,7 +988,14 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
 #define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
 #endif
 
-   (void) LOCAL_SPRINTF( value, "%s", (UseOriginalRng ? "Original Rng" : "Kiss Rng") );
+   
+   if( _randomNumberGenerator == Mersenne ){
+      (void) LOCAL_SPRINTF( value, "%s", "Mersenne Rng");
+   } else if( _randomNumberGenerator == Kiss ){
+      (void) LOCAL_SPRINTF( value, "%s", "Kiss Rng");
+   } else if( _randomNumberGenerator == Original ){
+      (void) LOCAL_SPRINTF( value, "%s", "Original Gromacs Rng");
+   }
    message << _getLine( tab, "Random number generator:", value ); 
 
    (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamSize() );
@@ -972,10 +1029,54 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
 
    message << _getLine( tab, "Shuffle:",              (_getShuffleStream()      ? Set : NotSet) ); 
  
+   message << getStatisticsString( );
+
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      message << std::endl;
+      if( _auxiliaryStreams[ii] ){
+         message << _auxiliaryStreams[ii]->getContentsString( );
+      }
+   }
+
+#undef LOCAL_SPRINTF
+
+   return message.str();
+}
+
+/* 
+ * Get statistics
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string BrookRandomNumberGenerator::getStatisticsString( void ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookRandomNumberGenerator::getStatisticsString";
+
+   static const unsigned int MAX_LINE_CHARS = 256;
+   char value[MAX_LINE_CHARS];
+   static const char* Set                   = "Set";
+   static const char* NotSet                = "Not set";
+   static double cumulativeStatistics[7]    = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0e+99, -1.0e+99 };
+
+// ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   std::string tab   = "   ";
+
+#ifdef WIN32
+#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
+#else
+#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
+#endif
+
    message << "\n   Statistics:\n";
 
    double statistics[7];
-   getStatistics( statistics, -1 );
+   getStatistics( statistics, -1, cumulativeStatistics );
 
    (void) LOCAL_SPRINTF( value, "%.5e", statistics[4] );
    message << _getLine( tab, "Count:",  value ); 
@@ -995,11 +1096,35 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
    (void) LOCAL_SPRINTF( value, "%.6e", statistics[6] );
    message << _getLine( tab, "Max:",  value ); 
 
-   for( int ii = 0; ii < LastStreamIndex; ii++ ){
-      message << std::endl;
-      if( _auxiliaryStreams[ii] ){
-         message << _auxiliaryStreams[ii]->getContentsString( );
+   if( cumulativeStatistics[4] > 1000.0 ){ 
+
+      for( int ii = 0; ii < 7; ii++ ){
+         statistics[ii] = cumulativeStatistics[ii];
+      }
+
+      statistics[0] /= statistics[4];
+      statistics[1]  = statistics[1] - statistics[4]*statistics[0]*statistics[0];
+      if( statistics[4] > 1.0 ){
+         statistics[1] = sqrt( statistics[1]/( statistics[4] - 1.0 ) );
       }
+      statistics[3]  = (statistics[3]/(statistics[4]*statistics[1]*statistics[1]) ) - 3.0;
+      (void) LOCAL_SPRINTF( value, "%.5e", statistics[4] );
+      message << _getLine( tab, "Cumulative Count:",  value ); 
+   
+      (void) LOCAL_SPRINTF( value, "%.5e", statistics[0] );
+      message << _getLine( tab, "Cumulative Average:",  value ); 
+   
+      (void) LOCAL_SPRINTF( value, "%.5e", statistics[1] );
+      message << _getLine( tab, "Cumulative StdDev:",  value ); 
+   
+      (void) LOCAL_SPRINTF( value, "%.5e", statistics[3] );
+      message << _getLine( tab, "Cumulative Kurtosis:",  value ); 
+   
+      (void) LOCAL_SPRINTF( value, "%.5e", statistics[5] );
+      message << _getLine( tab, "Cumulative Min:",  value ); 
+   
+      (void) LOCAL_SPRINTF( value, "%.6e", statistics[6] );
+      message << _getLine( tab, "Cumulative Max:",  value ); 
    }
    
 #undef LOCAL_SPRINTF
@@ -1025,7 +1150,7 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
  *
  * */
 
-int BrookRandomNumberGenerator::getStatistics( double statistics[7], int streamIndex ) const {
+int BrookRandomNumberGenerator::getStatistics( double statistics[7], int streamIndex, double cumulativeStatistics[7] ) const {
 
 // ---------------------------------------------------------------------------------------
 
@@ -1067,6 +1192,16 @@ int BrookRandomNumberGenerator::getStatistics( double statistics[7], int streamI
       }   
    }   
 
+   for( int ii = 0; ii < 5; ii++ ){
+      cumulativeStatistics[ii] += statistics[ii];
+   }
+   if( statistics[5] < cumulativeStatistics[5] ){
+      cumulativeStatistics[5] = statistics[5];
+   }
+   if( statistics[6] > cumulativeStatistics[6] ){
+      cumulativeStatistics[6] = statistics[6];
+   }
+
    if( statistics[4] > 0.0 ){ 
       statistics[0] /= statistics[4];
       statistics[1]  = statistics[1] - statistics[4]*statistics[0]*statistics[0];

platforms/brook/src/BrookRandomNumberGenerator.h

@@ -46,9 +46,9 @@ class BrookRandomNumberGenerator : public BrookCommon {
 
    public:
   
-      // toggle between original rng & Kiss (Nvidia) code
+      // toggle between original, Mersenne, & Kiss (Nvidia) random generators
 
-      static const int UseOriginalRng = 0;
+      enum Rngs { Original, Kiss, Mersenne };
  
       /** 
        * Constructor
@@ -129,6 +129,15 @@ class BrookRandomNumberGenerator : public BrookCommon {
       
       std::string getContentsString( int level = 0 ) const;
 
+      /* 
+       * Get stats
+       *
+       * @return string containing stats
+       *
+       * */
+      
+      std::string getStatisticsString( void ) const;
+
       /** 
        * Get random number stream 
        *
@@ -218,12 +227,13 @@ class BrookRandomNumberGenerator : public BrookCommon {
        *                       6: max
        *
        * @param streamIndex  stream index to analyze
+       * @param cumulativeStatistics accumulate stats array entries  same as statistics
        *
        * @return DefaultReturnValue
        *
-       * */
+       **/
       
-      int getStatistics( double statistics[7], int streamIndex ) const;
+      int getStatistics( double statistics[7], int streamIndex, double cumulativeStatistics[7] ) const;
       
       // ---------------------------------------------------------------------------------------
       
@@ -263,6 +273,8 @@ class BrookRandomNumberGenerator : public BrookCommon {
       float* _loadBuffer;
       int*   _shuffleIndices;
 
+      Rngs _randomNumberGenerator;
+
       /* 
        * Setup of stream dimensions
        *
@@ -351,6 +363,15 @@ class BrookRandomNumberGenerator : public BrookCommon {
       
       int _loadRandomNumberStreamsKiss( void );
 
+      /** 
+       * Load random number streams using Mersenne algorithm
+       * 
+       *
+       * @return DefaultReturnValue;
+       */
+      
+      int _loadRandomNumberStreamsMersenne( void );
+      
       /** 
        * Load random number streams using original gpu algorithm
        * 

platforms/brook/src/BrookVelocityCenterOfMassRemoval.cpp

@@ -148,7 +148,7 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImp
    // ---------------------------------------------------------------------------------------
 
    static const char* methodName  = "BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
-   static const int debug         = 1;
+   static const int debug         = 0;
 
    // ---------------------------------------------------------------------------------------
 
@@ -162,9 +162,10 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImp
       const float* vArray              = (float*) velV;
 
       int index                        = 0;
+      if( debug > 1 ){
          for( int ii = 0; ii < getNumberOfParticles(); ii++, index += 3 ){
             (void) fprintf( getLog(), "V %d [%12.5e %12.5e %12.5e]\n", ii, vArray[index], vArray[index+1], vArray[index+2] );
-
+         }
       }
 
       (void) fflush( getLog() );
@@ -197,12 +198,14 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImp
       void* w1                        = getWorkStream()->getData( 1 );
       const float* w2                 = (float*) w1;
 
+      if( debug > 1 ){
          int index                        = 0;
          for( int ii = 0; ii < getNumberOfParticles(); ii++, index += 3 ){
             (void) fprintf( getLog(), "V %d [%12.5e %12.5e %12.5e] [%12.5e %12.5e %12.5e]\n", ii,
                             vArray[index], vArray[index+1], vArray[index+2], w2[index], w2[index+1], w2[index+2] );
    
          }
+      }
 
       (void) fflush( getLog() );
 

platforms/brook/src/BrookVerletDynamics.cpp

@@ -171,31 +171,22 @@ int BrookVerletDynamics::updateParameters( double stepSize ){
  * @param  velocities                  particle velocities
  * @param  forces                      particle forces
  * @param  brookShakeAlgorithm         BrookShakeAlgorithm reference
- * @param  brookRandomNumberGenerator  BrookRandomNumberGenerator reference
  *
  * @return  DefaultReturnValue
  *
  */
 
-int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
-                                 const Stream& forces,
+int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImpl& velocityStream,
+                                 const BrookStreamImpl& forceStreamC,
                                  BrookShakeAlgorithm& brookShakeAlgorithm ){
 
 // ---------------------------------------------------------------------------------------
 
-   // unused Shake parameter
-
-   float numberOfIterations       = 25.0f;
-
    static const char* methodName  = "\nBrookVerletDynamics::update";
-
-   static const int PrintOn       = 0;
+   static const int PrintOn       = 1;
 
 // ---------------------------------------------------------------------------------------
 
-   BrookStreamImpl& positionStream                     = dynamic_cast<BrookStreamImpl&>       (positions.getImpl());
-   BrookStreamImpl& velocityStream                     = dynamic_cast<BrookStreamImpl&>       (velocities.getImpl());
-   const BrookStreamImpl& forceStreamC                 = dynamic_cast<const BrookStreamImpl&> (forces.getImpl());
    BrookStreamImpl& forceStream = const_cast<BrookStreamImpl&> (forceStreamC);
 
    if( (1 || PrintOn) && getLog() ){
@@ -220,14 +211,17 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
 
    }
 
+   // Shake step
+
+   if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
+
       // integration step
    
-   kupdate_md_verlet( (float) getStepSize(),
+      kupdate_md1( (float) getStepSize(),
                    positionStream.getBrookStream(),
                    velocityStream.getBrookStream(),
                    forceStream.getBrookStream(),
                    getInverseMassStream()->getBrookStream(),
-                       velocityStream.getBrookStream(),
                    getXPrimeStream()->getBrookStream() 
                  );
     
@@ -256,11 +250,6 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
          getXPrimeStream()->printToFile( getLog() ); 
       }   
 
-   // second Shake step
-
-   if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
-
-/*
       kshakeh_fix1( 
                     (float) brookShakeAlgorithm.getMaxIterations(),
                     (float) getVerletDynamicsParticleStreamWidth(),
@@ -274,8 +263,6 @@ int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
                     brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
-*/
-fprintf( stderr, "\nVerlet shake off!!\n" );
    
       if( (1|| PrintOn) && getLog() ){
 
@@ -351,9 +338,22 @@ fprintf( stderr, "\nVerlet shake off!!\n" );
 
       }   
 
+      // second integration step
+      kupdate_md2( (float) getStepSize(),
+                   getXPrimeStream()->getBrookStream(), 
+                   positionStream.getBrookStream(),
+                   velocityStream.getBrookStream(),
+                   positionStream.getBrookStream() );
    } else {
-      //kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
-      ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
+      //ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
+      float inverseStepSize = 1.0f/getStepSize();
+      kupdateMdNoShake( inverseStepSize,
+                        positionStream.getBrookStream(),
+                        velocityStream.getBrookStream(),
+                        forceStream.getBrookStream(),
+                        getInverseMassStream()->getBrookStream(),
+                        velocityStream.getBrookStream(),
+                        positionStream.getBrookStream() );
    }
 
 //_brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );

platforms/brook/src/BrookVerletDynamics.h

@@ -35,7 +35,7 @@
 #include <vector>
 #include <set>
 
-#include "BrookFloatStreamInternal.h"
+#include "BrookStreamImpl.h"
 #include "BrookShakeAlgorithm.h"
 #include "BrookPlatform.h"
 #include "BrookCommon.h"
@@ -131,8 +131,8 @@ class BrookVerletDynamics : public BrookCommon {
        *
        */
       
-      int update( Stream& positions, Stream& velocities,
-                  const Stream& forces, BrookShakeAlgorithm& brookShakeAlgorithm );
+      int update( BrookStreamImpl& positions, BrookStreamImpl& velocities,
+                  const BrookStreamImpl& forces, BrookShakeAlgorithm& brookShakeAlgorithm );
 
       /** 
        * Get array of VerletDynamics streams 

platforms/brook/src/gpu/kupdatemd.br

@@ -29,150 +29,39 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-kernel void kupdate_md_verlet( 
+kernel void kupdate_md1( 
 		float dt,
 		float3 posq<>, 
 		float3 v<>, 
 		float3 f<>,
 		float invmass<>,
-		out float3 vnew<>, 
 		out float3 posqp<> ){
 
-	vnew   = v + f*invmass*dt;
-	posqp  = posq + vnew*dt;
+   posqp  = v + dt*invmass*f;
+	posqp *= dt;
 }
 
-kernel void kupdate_md1_berendsen( 
-		float dt,
-		float3 lg, //Berendsen coupling, assuming only one group
-		float3 uold, //Mean velocity from previous step
-		float4 posq<>, 
-		float3 v<>, 
-		float3 f<>,
-		float invmass<>,
-		out float3 vnew<>, 
-		out float4 posqp<> )
-{
-	float3 vb;
-	float3 one;
-	posqp = posq;
-
-	one = float3( 1.0f, 1.0f, 1.0f );
-
-	vb = ( one - lg ) * uold;
-	
-	vnew = lg* ( v + f * invmass * dt ) + vb;
- 	
-	posqp.xyz += vnew * dt;
-}
-
-//Nose-Hoover / Parinello Rahman
-kernel void kupdate_md1_extended (
-		float dt,
-		float3 lg, 
-		float xi,
-		float3 M0, 
-		float3 M1,
-		float3 M2,
-		float3 uold,
-		float4 posq<>, 
-		float3 v<>, 
-		float3 f<>,
-		float invmass<>,
-		out float3 vnew<>, 
-		out float4 posqp<> )
-{
-	float3 vrel, vnrel;
-	float3 vtrans;
-
-	vrel  = v - uold;
-	vtrans = float3( dot(M0, vrel), dot(M1, vrel), dot(M2, vrel) );
-	vrel += dt * ( invmass * f - xi * vrel - vtrans );
-	vnew = uold + lg * vrel;
-	posqp = posq;
-	posqp.xyz += vnew * dt;
-}
-
-		
 kernel void kupdate_md2(
 		float dtinv, //1/dt
-		float4 posqp<>, //positions after constraints
-		float4 posq<>,  //positions before update
+		float3 posqp<>, //positions after constraints
+		float3 posq<>,  //positions before update
 		out float3 vnew<>, //Corrected velocities
-		out float4 posqnew<> //equal to posqp, avoids an extra call to copy
-		)
-{
-	posqnew = posqp;
-	vnew = ( posqp - posq ) * dtinv;
+		out float3 posqnew<> //equal to posqp, avoids an extra call to copy
+		){
+   posqnew     = posq  + posqp;
+   vnew        = posqp * dtinv;
 }
 
-/* Version that handles terms of order dt more carefully
- * Update1 computes deltax's rather than x + deltax 
- * Corresponding shake implementation modifies the delta's
- * Update2 adds the constrained deltas to x.
- * */
-
-kernel void kupdate_md1_berendsen_fix1( 
+kernel void kupdateMdNoShake(
       float dt,
-		float3 lg, //Berendsen coupling, assuming only one group
-		float3 uold, //Mean velocity from previous step
-		float4 posq<>, 
-		float3 v<>, 
-		float3 f<>,
-		float invmass<>,
-		out float3 vnew<>, 
-		out float4 posqp<> )
-{
-	float3 vb;
-	float3 one;
-
-	one = float3( 1.0f, 1.0f, 1.0f );
-
-	vb = ( one - lg ) * uold;
-	
-	vnew = lg* ( v + f * invmass * dt ) + vb;
- 	
-	posqp.xyz = vnew * dt;
-	posqp.w = posq.w;
-}
-
-//Nose-Hoover / Parinello Rahman
-kernel void kupdate_md1_extended_fix1 (
-		float dt,
-		float3 lg, 
-		float xi,
-		float3 M0, 
-		float3 M1,
-		float3 M2,
-		float3 uold,
-		float4 posq<>, 
+      float3 posq<>,
       float3 v<>, 
       float3 f<>,
       float invmass<>,
-		out float3 vnew<>, 
-		out float4 posqp<> )
-{
-	float3 vrel, vnrel;
-	float3 vtrans;
-
-	vrel  = v - uold;
-	vtrans = float3( dot(M0, vrel), dot(M1, vrel), dot(M2, vrel) );
-	vrel += dt * ( invmass * f - xi * vrel - vtrans );
-	vnew = uold + lg * vrel;
-	posqp.w = posq.w;
-	posqp.xyz = vnew * dt;
-}
-
-kernel void kupdate_md2_fix1(
-		float dtinv, //1/dt
-		float4 posqp<>, //positions after constraints
-		float4 posq<>,  //positions before update
-		out float3 vnew<>, //Corrected velocities
-		out float4 posqnew<> //equal to posqp, avoids an extra call to copy
-		)
-{
-	posqnew.xyz = posq.xyz + posqp.xyz;
-	posqnew.w   = posq.w;
-	vnew = posqp.xyz * dtinv;
+      out float3 outv<>,
+      out float3 posqp<> ){
+   posqp     = posq;
+   outv      = v + dt*invmass*f;
+   posqp    += dt*outv;
 }  
 

platforms/brook/src/gpu/kupdatemd.h

@@ -29,75 +29,27 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-void  kupdate_md2 (const float  dtinv,
-		::brook::stream posqp,
-		::brook::stream posq,
-		::brook::stream vnew,
-		::brook::stream posqnew
-	    ); 
-
-void  kupdate_md1_extended (
-		const float  dt,
-		const float3  lg,
-		const float  xi,
-		const float3  M0,
-		const float3  M1,
-		const float3  M2,
-		const float3  uold,
-		::brook::stream posq,
-		::brook::stream v,
-		::brook::stream f,
-		::brook::stream invmass,
-		::brook::stream vnew,
-		::brook::stream posqp);
-
-void  kupdate_md_verlet(
-		const float  dt,
-		::brook::stream posq,
-		::brook::stream v,
-		::brook::stream f,
-		::brook::stream invmass,
-		::brook::stream vnew,
-		::brook::stream posqp);
-
-
-void  kupdate_md1_berendsen (
-		const float  dt,
-		const float3  lg,
-		const float3  uold,
+void  kupdate_md1(
+         const float  dtinv,
    		::brook::stream posq,
    		::brook::stream v,
    		::brook::stream f, 
    		::brook::stream invmass, 
-		::brook::stream vnew,
-		::brook::stream posqp); 
+   		::brook::stream posqnew ); 
 
-void  kupdate_md2_fix1 (const float  dtinv,
+void  kupdate_md2(
+         const float  dtinv,
    		::brook::stream posqp,
    		::brook::stream posq,
 	   	::brook::stream vnew,
-		::brook::stream posqnew); 
+	   	::brook::stream posqnew ); 
 
-void  kupdate_md1_extended_fix1 (const float  dt,
-		const float3  lg,
-		const float  xi,
-		const float3  M0,
-		const float3  M1,
-		const float3  M2,
-		const float3  uold,
+void  kupdateMdNoShake(
+         const float  dtinv,
    		::brook::stream posq,
 	   	::brook::stream v,
 	   	::brook::stream f,
    		::brook::stream invmass, 
-		::brook::stream vnew,
-		::brook::stream posqp);
+	   	::brook::stream outv,
+	   	::brook::stream posqp ); 
 
-void  kupdate_md1_berendsen_fix1 (const float  dt,
-		const float3  lg,
-		const float3  uold,
-		::brook::stream posq,
-		::brook::stream v,
-		::brook::stream f,
-		::brook::stream invmass,
-		::brook::stream vnew,
-		::brook::stream posqp);