Mods

platforms/brook/src/BrookBonded.cpp

@@ -1911,7 +1911,7 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
 /*
       (void) fprintf( log, "\nNB1 forces" );
-      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
       brookStreamInternalF->printToFile( log );
 */
 
@@ -2036,7 +2036,7 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
    if( printOn ){
       (void) fprintf( log, "%s Post 3_4/3_5 forces\n", methodName.c_str() );
-      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
       brookStreamInternalF->printToFile( log );
    }
 
@@ -2103,7 +2103,7 @@ void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImp
 
    if( printOn ){
       (void) fprintf( log, "%s Final forces", methodName.c_str() );
-      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
       brookStreamInternalF->printToFile( log );
    }
 

platforms/brook/src/BrookBrownianDynamics.cpp

@@ -448,13 +448,13 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
       //StreamImpl& positionStreamImpl               = positionStream.getImpl();
       //const BrookStreamImpl brookPositions         = dynamic_cast<BrookStreamImpl&> (positionStreamImpl);
 /*
-      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
       (void) fprintf( getLog(), "\nPositionStream\n" );
       brookStreamInternalPos->printToFile( getLog() );
 */
 
       double forceSum[3];
-      BrookStreamInternal* brookStreamInternalFF       = forceStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalFF       = forceStream.getBrookStreamInternal();
       BrookFloatStreamInternal* brookStreamInternalF   = dynamic_cast<BrookFloatStreamInternal*> (brookStreamInternalFF);
       brookStreamInternalF->sumByDimension( getNumberOfParticles(), forceSum );
       (void) fprintf( getLog(), "\nForceStream [%18.10e %18.10e %18.10e]\n", forceSum[0], forceSum[1], forceSum[2] );
@@ -468,7 +468,7 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
 */
       // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); 
 
-         BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamInternal();
          (void) fprintf( getLog(), "\nVelocityStream\n" );
          brookStreamInternalVel->printToFile( getLog() ); 
    }   
@@ -499,12 +499,12 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
 
       //StreamImpl& positionStreamImpl               = positionStream.getImpl();
       //const BrookStreamImpl brookPositions         = dynamic_cast<BrookStreamImpl&> (positionStreamImpl);
-      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
       (void) fprintf( getLog(), "\nPositionStream\n" );
       brookStreamInternalPos->printToFile( getLog() );
 
       (void) fprintf( getLog(), "\nForceStream\n" );
-      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
       brookStreamInternalF->printToFile( getLog() );
 
       (void) fprintf( getLog(), "\nXPrimeStream\n" );
@@ -557,14 +557,14 @@ fprintf( stderr, "\nBrownian shake off!!\n" );
    
          (void) fprintf( getLog(), "\nPost kupdate_bd:\n" );
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          (void) fprintf( getLog(), "\nPositionStream\n" );
          brookStreamInternalPos->printToFile( getLog() );
    
          (void) fprintf( getLog(), "\nXPrimeStream\n" );
          getXPrimeStream()->printToFile( getLog() ); 
    
-         BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamInternal();
          (void) fprintf( getLog(), "\nVelocityStream\n" );
          brookStreamInternalVel->printToFile( getLog() ); 
    
@@ -587,11 +587,11 @@ fprintf( stderr, "\nBrownian shake off!!\n" );
          (void) fprintf( getLog(), "\nXPrimeStream\n" );
          getXPrimeStream()->printToFile( getLog() ); 
 
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          (void) fprintf( getLog(), "\nPositionStream\n" );
          brookStreamInternalPos->printToFile( getLog() );
    
-         BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamInternal();
          (void) fprintf( getLog(), "\nVelocityStream\n" );
          brookStreamInternalVel->printToFile( getLog() ); 
    

platforms/brook/src/BrookCommon.cpp

@@ -92,6 +92,7 @@ const std::string BrookCommon::ShakeInverseMapStream
 // Random number streams
 
 const std::string BrookCommon::ShuffleStream                                      = "ShuffleStream";
+const std::string BrookCommon::RandomValuesStream                                 = "RandomValuesStream";
 
 // Random number streams
 

platforms/brook/src/BrookCommon.h

@@ -107,6 +107,7 @@ class BrookCommon {
       // Random number generator streams
 
       static const std::string ShuffleStream;
+      static const std::string RandomValuesStream;
 
       // BrookVelocityCenterOfMassRemoval streams
 

platforms/brook/src/BrookGbsa.cpp

@@ -965,7 +965,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
                       getParticleStreamWidth( ),
                       getPartialForceStreamWidth( ) );
 
-      BrookStreamInternal* brookStreamInternalF  = positionStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF  = positionStream.getBrookStreamInternal();
       (void) fprintf( log, "\nPositionStream\n" );
       brookStreamInternalF->printToFile( log );
 
@@ -1006,7 +1006,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
                       getParticleStreamWidth( ),
                       getPartialForceStreamWidth( ) );
 
-      BrookStreamInternal* brookStreamInternalF  = positionStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF  = positionStream.getBrookStreamInternal();
       (void) fprintf( log, "\nPositionStream\n" );
       brookStreamInternalF->printToFile( log );
 
@@ -1060,7 +1060,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
                       getSoluteDielectric(),
                       getSolventDielectric(), includeAceTerm );
 
-      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
       (void) fprintf( log, "\nPost kObcLoop1 PositionStream\n" );
       brookStreamInternalPos->printToFile( log );
 
@@ -1167,7 +1167,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
                       getParticleStreamWidth( ),
                       getPartialForceStreamWidth( ) );
 
-      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
       (void) fprintf( log, "\nPost kObcLoop2: PositionStream\n" );
       brookStreamInternalPos->printToFile( log );
 
@@ -1230,7 +1230,7 @@ void BrookGbsa::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl&
          gbsaForceStreams[ii]->printToFile( log );
       }
 
-      BrookStreamInternal* brookStreamInternalF  = forceStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF  = forceStream.getBrookStreamInternal();
       (void) fprintf( log, "\nPost kPostObcLoop2_nobranch:  ForceStream\n" );
       brookStreamInternalF->printToFile( log );
 

platforms/brook/src/BrookLangevinDynamics.cpp

 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
- * -------------------------------------------------------------------------- *
+ * -i------------------------------------------------------------------------- *
  * This is part of the OpenMM molecular simulation toolkit originating from   *
  * Simbios, the NIH National Center for Physics-Based Simulation of           *
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
@@ -727,11 +727,32 @@ int BrookLangevinDynamics::_updateSdStreams( void ){
    memset( sd2x, 0, 3*sdParticleStreamSize*sizeof( BrookOpenMMFloat ) ); 
 
    index = 0;
+   int useFixedRandomValue = 0;
+   if( useFixedRandomValue ){
+
+      // diagnostics only!
+
+      BrookOpenMMFloat fixedRandomValue = static_cast<BrookOpenMMFloat>( 0.1 );
+      for( int ii = 0; ii < numberOfParticles; ii++, index += 3 ){
+         BrookOpenMMFloat value = _inverseSqrtMasses[ii]*derivedParameters[X]*fixedRandomValue;
+         sd2x[index]            = value;
+         sd2x[index+1]          = value;
+         sd2x[index+2]          = value;
+      }
+
+      // print message letting user know non-random value being used
+
+      FILE* log = getLog() ? getLog() : stderr;
+      (void) fprintf( log, "%s using fixed 'random value'=%.3f to initialize sd2x\n", methodName.c_str(), fixedRandomValue );
+
+   } else {
+
       for( int ii = 0; ii < numberOfParticles; ii++, index += 3 ){
          sd2x[index]        = _inverseSqrtMasses[ii]*derivedParameters[X]*( static_cast<BrookOpenMMFloat> (SimTKOpenMMUtilities::getNormallyDistributedRandomNumber()) );
          sd2x[index+1]      = _inverseSqrtMasses[ii]*derivedParameters[X]*( static_cast<BrookOpenMMFloat> (SimTKOpenMMUtilities::getNormallyDistributedRandomNumber()) );
          sd2x[index+2]      = _inverseSqrtMasses[ii]*derivedParameters[X]*( static_cast<BrookOpenMMFloat> (SimTKOpenMMUtilities::getNormallyDistributedRandomNumber()) );
       }
+   }
    
    _sdStreams[SD2XStream]->loadFromArray( sd2x );
 
@@ -837,18 +858,24 @@ float BrookLangevinDynamics::getTemperature( BrookStreamInternal* velocities, Br
    float ke                   = 0.0f;
    int index                  = 0;
 
-   for( int ii = 0; ii < getNumberOfParticles(); ii++ ){
+   int numberOfParticles      = getNumberOfParticles();
+   for( int ii = 0; ii < numberOfParticles; ii++, index += 3 ){
       ke    += (velocitiesI[index]*velocitiesI[index] + velocitiesI[index+1]*velocitiesI[index+1] + velocitiesI[index+2]*velocitiesI[index+2] )/inverseMassStreamI[ii];
-//(void) fprintf( stderr, "%s %d Velocities[%14.5e %14.5e %14.5e]\n", methodName.c_str(), index/3, velocitiesI[index], velocitiesI[index+1], velocitiesI[index+2] ); 
-      index += 3;
+
+/*
+if( ii < 3 || ii >= (numberOfParticles - 3) ){
+(void) fprintf( stderr, "%s %6d m=%14.5e v[%14.5e %14.5e %14.5e]\n", methodName.c_str(), ii, 1.0/inverseMassStreamI[ii],
+                velocitiesI[index], velocitiesI[index+1], velocitiesI[index+2] ); 
+} */
+
    }
 
    int degreesOfFreedom = 3*getNumberOfParticles() - numberOfConstraints;
+   float denominator    = 1.0f/( ((float) BOLTZ)*((float) ( degreesOfFreedom )) );
+//(void) fprintf( stderr, "%s ke=%.5e T=%.3f dof=%d\n", methodName.c_str(), ke, (ke*denominator),  degreesOfFreedom ); 
+   ke *= denominator;
 
-(void) fprintf( stderr, "%s ke=%.5e dof=%d\n", methodName.c_str(), ke, degreesOfFreedom ); 
-   ke /= ((float) BOLTZ)*((float) ( degreesOfFreedom ));
-
-   return 0.5f*ke;
+   return ke;
 }
 
 /** 
@@ -930,6 +957,41 @@ int BrookLangevinDynamics::removeCom( BrookStreamInternal* velocities, BrookFloa
    return DefaultReturnValue;
 }
 
+/** 
+ * Reset velocities (diagnostics)
+ * 
+ * @param velocities             velocities
+ *
+ * @return DefaultReturnValue  
+ */
+
+int BrookLangevinDynamics::resetVelocities( BrookStreamInternal* velocities ) const {
+    
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookLangevinDynamics::resetVelocities";
+
+// ---------------------------------------------------------------------------------------
+
+(void) fprintf( stderr, "%s\n", methodName.c_str() ); fflush( stderr );
+
+   // resset v to determinisitic values
+
+   double* newVelocities = new double[velocities->getStreamSize()*velocities->getWidth()];
+   memset( newVelocities, 0, sizeof( double )*velocities->getStreamSize()*velocities->getWidth() );
+   for( int ii = 1; ii <= 3*getNumberOfParticles(); ii++ ){
+      int jj                = ii % 10;
+      double sign           = jj % 2 ? 0.1 : -0.1;
+      newVelocities[ii-1]   = sign*( (double) jj); 
+   }
+
+   velocities->loadFromArray( newVelocities );
+
+   delete[] newVelocities;
+
+   return DefaultReturnValue;
+}
+
 /* 
  * Get contents of object
  *
@@ -1031,14 +1093,14 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName = "\nBrookLangevinDynamics::update";
-   int printOn                         = 1;
+   int printOn                         = 0;
    FILE* log;
 
 // ---------------------------------------------------------------------------------------
 
    _internalStepCount++;
 
-// setLog( stderr );
+//setLog( stderr );
 
    if( printOn && getLog() ){
       log = getLog();
@@ -1075,13 +1137,34 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
       }
    }
 
-   if( 1 || (_internalStepCount % 10) == 0 ){
+   if( 0 && _internalStepCount == 1 ){ 
+      //resetVelocities( velocityStream.getBrookStreamInternal() );
+      std::string velocityFileName = "kupdate_sd1_strV.in";
+      std::string rvName0          = "kupdate_sd1_0.fgauss.in";
+      std::string rvName1          = "kupdate_sd1_1.fgauss.in";
+      std::string sd2xName         = "kupdate_sd1_strSD2X.in";
+      velocityStream.getBrookStreamInternal()->loadStreamGivenFileName( velocityFileName );
+      brookRandomNumberGenerator.getRandomNumberStream( 0 )->loadStreamGivenFileName( rvName0 );
+      brookRandomNumberGenerator.getRandomNumberStream( 1 )->loadStreamGivenFileName( rvName1 );
+      getSD2XStream()->loadStreamGivenFileName( sd2xName );
+   }
+
+   if( 0 && (_internalStepCount % 10) == 0 ){
       FILE*	log1 = stderr;
-      BrookStreamInternal* brookStreamInternalPos  = velocityStream.getBrookStreamImpl();
+      (void) fprintf( log1, "\nVelocityStream %d XX\n", _internalStepCount ); fflush( log1 );
+      BrookStreamInternal* brookStreamInternalPos  = velocityStream.getBrookStreamInternal();
       float temperature = getTemperature( brookStreamInternalPos, getInverseMassStream(), brookShakeAlgorithm.getNumberOfConstraints() );
 //      removeCom( brookStreamInternalPos, getInverseMassStream() );
       float temperaturePost = getTemperature( brookStreamInternalPos, getInverseMassStream(), brookShakeAlgorithm.getNumberOfConstraints() );
       (void) fprintf( log1, "\nVelocityStream %d Tp=%.3f %.3f\n", _internalStepCount, temperature, temperaturePost );
+      (void) fprintf( log1, "\n%s step=%d Post kupdate_sd1_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
+                                "EM=%12.5e Sd1pc[]=[%12.5e %12.5e %12.5e]", methodName.c_str(), _internalStepCount,
+                                getLangevinDynamicsParticleStreamWidth(),
+                                brookRandomNumberGenerator.getRandomNumberStreamWidth(),
+                                brookRandomNumberGenerator.getRvStreamOffset(),
+                                derivedParameters[EM], derivedParameters[Sd1pc1], derivedParameters[Sd1pc2], derivedParameters[Sd1pc3] );
+      (void) fprintf( log1, "\nRvStreamIndex=%d step=%d\n", brookRandomNumberGenerator.getRvStreamIndex(), _internalStepCount );
+      fflush( log1 );
    }
 
    // first integration step
@@ -1103,8 +1186,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
          getInverseMassStream()->getBrookStream(),
          getSD1VStream()->getBrookStream(),
          getVPrimeStream()->getBrookStream(),
-         getXPrimeStream()->getBrookStream() 
-         );
+         getXPrimeStream()->getBrookStream() );
  
    // diagnostics
 
@@ -1117,38 +1199,46 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
                                 derivedParameters[EM], derivedParameters[Sd1pc1], derivedParameters[Sd1pc2], derivedParameters[Sd1pc3] );
 
       if( _internalStepCount == 1 ){
-         (void) fprintf( log, "\nSDPC1Stream %d\n", _internalStepCount );
+         (void) fprintf( log, "\nSDPC1Stream fixed input sd1 step=%d\n", _internalStepCount );
          getSDPC1Stream()->printToFile( log );
    
-         (void) fprintf( log, "\nSD2XStream %d\n", _internalStepCount );
+         (void) fprintf( log, "\nSD2XStream fixed input sd1 step=%d\n", _internalStepCount );
          getSD2XStream()->printToFile( log );
    
-         (void) fprintf( log, "\nInverseMassStream %d\n", _internalStepCount );
+         (void) fprintf( log, "\nInverseMassStream fixed input sd1 step=%d\n", _internalStepCount );
          getInverseMassStream()->printToFile( log );
    
       }
 
-      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
-      (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
+      (void) fprintf( log, "\nPositionStream input sd1 step=%d\n", _internalStepCount );
       brookStreamInternalPos->printToFile( log );
 
-      (void) fprintf( log, "\nForceStream %d\n", _internalStepCount );
-      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      (void) fprintf( log, "\nForceStream input sd1 step=%d\n", _internalStepCount );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
       brookStreamInternalF->printToFile( log );
 
-      BrookStreamInternal* brookStreamInternalV   = velocityStream.getBrookStreamImpl();
-      (void) fprintf( log, "\nVelocityStream %d\n", _internalStepCount );
+      BrookStreamInternal* brookStreamInternalV   = velocityStream.getBrookStreamInternal();
+      (void) fprintf( log, "\nVelocityStream input sd1 step=%d\n", _internalStepCount );
       brookStreamInternalV->printToFile( log );
 
-      (void) fprintf( log, "\nSD1VStream %d\n", _internalStepCount );
+      (void) fprintf( log, "\nSD1VStream output sd1 step=%d\n", _internalStepCount );
       getSD1VStream()->printToFile( log );
 
-      (void) fprintf( log, "\nVPrimeStream %d\n", _internalStepCount );
+      (void) fprintf( log, "\nVPrimeStream output sd1 step=%d\n", _internalStepCount );
       getVPrimeStream()->printToFile( log );
 
-      (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
+      (void) fprintf( log, "\nXPrimeStream output sd1 step=%d\n", _internalStepCount );
       getXPrimeStream()->printToFile( log ); 
 
+      if( _internalStepCount == 1 ){
+            std::vector<BrookStreamInternal*> streams;
+            streams.push_back( brookStreamInternalPos );
+            streams.push_back( getXPrimeStream() );
+            std::stringstream fileNameBaseS;
+            fileNameBaseS << "Brook_Sd1PreShk_" << _internalStepCount << ".txt"; 
+            BrookStreamInternal::printStreamsToFile( fileNameBaseS.str(), streams );
+      }
       (void) fprintf( log, "\nRvStreamIndex=%d step=%d\n", brookRandomNumberGenerator.getRvStreamIndex(), _internalStepCount );
       // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log ); 
    }   
@@ -1186,35 +1276,45 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
                     brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
                     getXPrimeStream()->getBrookStream() );
 
-      if( 0 && printOn ){
+      if( printOn ){
 
-         (void) fprintf( log, "\n%s Post kshakeh_update2_fix1: ShakeConstraintStreamWidth=%3d",
-                                   methodName.c_str(), brookShakeAlgorithm.getShakeConstraintStreamWidth() );
+         (void) fprintf( log, "\n%s Post kshakeh_update1: sw=%d ShkCnstStrW=%3d tol=%.3f maxIt=%d",
+                         methodName.c_str(), getLangevinDynamicsParticleStreamWidth(),
+                         brookShakeAlgorithm.getShakeConstraintStreamWidth(),
+                         brookShakeAlgorithm.getShakeTolerance(), brookShakeAlgorithm.getMaxIterations() );
    
          if( _internalStepCount == 1 ){
-            (void) fprintf( log, "\nShakeInverseMapStream %d\n" );
+            (void) fprintf( log, "\nShakeInverseMapStream fixed input sd shake1 at step=%d\n" );
             brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( log );
          }
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
-         (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
+         (void) fprintf( log, "\nPositionStream input sd shake1 at step=%d\n", _internalStepCount );
          brookStreamInternalPos->printToFile( log );
    
-         (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
+         (void) fprintf( log, "\nXPrimeStream output sd shake1 at step=%d\n", _internalStepCount );
          getXPrimeStream()->printToFile( log );  
 
-         (void) fprintf( log, "\nShakeXCons0 %d\n", _internalStepCount );
+         (void) fprintf( log, "\nShakeXCons0 output sd shake1 at step=%d\n", _internalStepCount );
          brookShakeAlgorithm.getShakeXCons0Stream()->printToFile( log );  
 
-         (void) fprintf( log, "\nShakeXCons1 %d\n", _internalStepCount );
+         (void) fprintf( log, "\nShakeXCons1 output sd shake1 at step=%d\n", _internalStepCount );
          brookShakeAlgorithm.getShakeXCons1Stream()->printToFile( log );  
 
-         (void) fprintf( log, "\nShakeXCons2 %d\n", _internalStepCount );
+         (void) fprintf( log, "\nShakeXCons2 output sd shake1 at step=%d\n", _internalStepCount );
          brookShakeAlgorithm.getShakeXCons2Stream()->printToFile( log );  
 
-         (void) fprintf( log, "\nShakeXCons3 %d\n", _internalStepCount );
+         (void) fprintf( log, "\nShakeXCons3 output sd shake1 at step=%d\n", _internalStepCount );
          brookShakeAlgorithm.getShakeXCons3Stream()->printToFile( log );  
 
+         if( _internalStepCount < 2 ){
+            std::vector<BrookStreamInternal*> streams;
+            streams.push_back( brookStreamInternalPos );
+            streams.push_back( getXPrimeStream() );
+            std::stringstream fileNameBaseS;
+            fileNameBaseS << "Brook_Sd1PostShk_" << _internalStepCount << ".txt"; 
+            BrookStreamInternal::printStreamsToFile( fileNameBaseS.str(), streams );
+         }
       }   
 
    }
@@ -1248,31 +1348,40 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
                                 brookRandomNumberGenerator.getRvStreamOffset(),
                                 derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
 
-      (void) fprintf( log, "\nSDPC2Stream %d \n", _internalStepCount );
+      (void) fprintf( log, "\nSDPC2Stream input sd2 step=%d\n", _internalStepCount );
       getSDPC2Stream()->printToFile( log );
 
-      (void) fprintf( log, "\nSD2XStream %d \n", _internalStepCount );
+      (void) fprintf( log, "\ngetSD1VStream input sd2 step=%d\n", _internalStepCount );
       getSD1VStream()->printToFile( log );
 
-      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
-      (void) fprintf( log, "\nPositionStream %d \n", _internalStepCount );
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
+      (void) fprintf( log, "\nPositionStream input sd2 step=%d \n", _internalStepCount );
       brookStreamInternalPos->printToFile( log );
 
-      (void) fprintf( log, "\nVPrimeStream %d \n", _internalStepCount );
+      (void) fprintf( log, "\nVPrimeStream input sd2 step=%d\n", _internalStepCount );
       getVPrimeStream()->printToFile( log );
 
-      (void) fprintf( log, "\nXPrimeStream %d \n", _internalStepCount );
+      (void) fprintf( log, "\nXPrimeStream output sd2 step=%d\n", _internalStepCount );
       getXPrimeStream()->printToFile( log ); 
 
-      (void) fprintf( log, "\nSD2XStream %d \n", _internalStepCount );
+      (void) fprintf( log, "\nSD2XStream output sd2 step=%d\n", _internalStepCount );
       getSD2XStream()->printToFile( log );
 
-      BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamImpl();
-      (void) fprintf( log, "\nVelocityStream %d \n", _internalStepCount );
+      BrookStreamInternal* brookStreamInternalVel = velocityStream.getBrookStreamInternal();
+      (void) fprintf( log, "\nVelocityStream output sd2 step=%d\n", _internalStepCount );
       brookStreamInternalVel->printToFile( log );
 
       (void) fprintf( log, "\nRvStreamIndex=%d step=%d\n", brookRandomNumberGenerator.getRvStreamIndex(), _internalStepCount );
 //      brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log ); 
+ 
+      if( _internalStepCount < 3 ){
+         std::vector<BrookStreamInternal*> streams;
+         streams.push_back( brookStreamInternalVel );
+         streams.push_back( getXPrimeStream() );
+         std::stringstream fileNameBaseS;
+         fileNameBaseS << "Brook_Sd2Out_" << _internalStepCount << ".txt"; 
+         BrookStreamInternal::printStreamsToFile( fileNameBaseS.str(), streams );
+      }
    }   
 
    // advance random number cursor
@@ -1319,12 +1428,12 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
                               brookRandomNumberGenerator.getRvStreamOffset(),
                               derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, log, 0.0001f );
-         (void) fprintf( log, "\nPositionStream\n" );
+         (void) fprintf( log, "\nPositionStream output sd shake2 step=%d\n", _internalStepCount );
          brookStreamInternalPos->printToFile( log );
    
-         BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
          std::vector<std::vector<double> > forceStatistics;
          brookStreamInternalF->getStatistics( forceStatistics, getNumberOfParticles() );
 
@@ -1334,7 +1443,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
          (void) fprintf( log, "\nStep %d Force stats:\n%s", _internalStepCount, stats.c_str() );
          brookStreamInternalF->printToFile( log );
 
-         brookStreamInternalPos  = velocityStream.getBrookStreamImpl();
+         brookStreamInternalPos  = velocityStream.getBrookStreamInternal();
          std::vector<std::vector<double> > velocityStatistics;
          brookStreamInternalPos->getStatistics( velocityStatistics, getNumberOfParticles() );
          std::stringstream tagV;
@@ -1345,7 +1454,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
          (void) fprintf( log, "\nVelocityStream %d T=%.3f\n", _internalStepCount, temperature );
          brookStreamInternalPos->printToFile( log );
    
-         (void) fprintf( log, "\nXPrimeStream\n" );
+         (void) fprintf( log, "\nXPrimeStream input sd shake2 step=%d\n", _internalStepCount );
          getXPrimeStream()->printToFile( log ); 
    
 //      brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log ); 
@@ -1360,7 +1469,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 //         (void) fprintf( log, "\nXPrimeStream\n" );
 //         getXPrimeStream()->printToFile( log );
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          (void) fprintf( log, "\nPositionStream\n" );
          brookStreamInternalPos->printToFile( log );
       }
@@ -1376,7 +1485,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
          (void) fprintf( log, "\nXPrimeStream\n" );
          getXPrimeStream()->printToFile( log );
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          (void) fprintf( log, "\nPositionStream\n" );
          brookStreamInternalPos->printToFile( log );
    
@@ -1384,15 +1493,17 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
       }   
    }
 
-   if( (_internalStepCount % 10) == 0 ){
+   if( (_internalStepCount % 400) == 0 ){
       FILE*	log1 = stderr;
-      BrookStreamInternal* brookStreamInternalPos  = velocityStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalPos  = velocityStream.getBrookStreamInternal();
+/*
       std::vector<std::vector<double> > velocityStatistics;
       brookStreamInternalPos->getStatistics( velocityStatistics, getNumberOfParticles() );
       std::stringstream tagV;
       tagV << _internalStepCount << " Vxx "; 
       std::string stats = brookStreamInternalPos->printStatistics( tagV.str(), velocityStatistics );
       (void) fprintf( log1, "\nStep %d Velocity stats:\n%s", _internalStepCount, stats.c_str() );
+*/
       float temperature = getTemperature( brookStreamInternalPos, getInverseMassStream(), brookShakeAlgorithm.getNumberOfConstraints() );
       //removeCom( brookStreamInternalPos, getInverseMassStream() );
       (void) fprintf( log1, "\nVelocityStream %d T=%.3f\n", _internalStepCount, temperature );

platforms/brook/src/BrookLangevinDynamics.h

@@ -283,6 +283,16 @@ class BrookLangevinDynamics : public BrookCommon {
       
       int removeCom( BrookStreamInternal* velocities, BrookFloatStreamInternal* inverseMassStream ) const;
       
+      /** 
+       * Reset velocities (diagnostics)
+       * 
+       * @param velocities             velocities
+       *
+       * @return DefaultReturnValue  
+       */
+      
+      int resetVelocities( BrookStreamInternal* velocities ) const;
+      
    private:
    
       enum DerivedParameters { GDT, EPH, EMH, EP, EM, B, C, D, V, X, Yv, Yx,

platforms/brook/src/BrookNonBonded.cpp

@@ -1341,7 +1341,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
       (void) fprintf( log, "\n                      pFrc=%6d     eps=%12.5e\n",       getPartialForceStreamWidth( ), epsfac );
 
       (void) fprintf( log, "\nFinal NB & bonded forces" );
-      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
       brookStreamInternalF->printToFile( log );
  
       (void) fprintf( log, "\nOuterVdwStreamd\n" );
@@ -1387,7 +1387,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
    if( printOn ){
 
       (void) fprintf( log, "\n%s NB forces\n", methodName.c_str() );
-      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
       brookStreamInternalF->printToFile( log );
       (void) fprintf( log, "\n%s Done\n", methodName.c_str() ); (void) fflush( log );
 

platforms/brook/src/BrookPlatform.cpp

@@ -44,7 +44,6 @@
 
 using namespace OpenMM;
 
-
 /** 
  * BrookPlatformData constructor
  *
@@ -314,9 +313,12 @@ void BrookPlatform::_setBrookRuntime( const std::string& runtime ){
       throw OpenMMException( message.str() );
    }
 
-   if( getLog() ){
-      (void) fprintf( getLog(), "%s Brook initializing to runtime=<%s>\n", methodName.c_str(), _runtime.c_str() ); 
-      (void) fflush( getLog() );
+   // let user know runtime setting
+
+   if( 1 ){
+      FILE* log = getLog() ? getLog() : stderr;
+      (void) fprintf( log, "%s Brook initializing to runtime=<%s>\n", methodName.c_str(), _runtime.c_str() ); 
+      (void) fflush( log );
    }
 
    brook::initialize( _runtime.c_str(), NULL );

platforms/brook/src/BrookRandomNumberGenerator.cpp

@@ -81,6 +81,7 @@ BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){
    _randomNumberSeed      = 1393;
    //_randomNumberGenerator = Mersenne;
    _randomNumberGenerator = Kiss;
+   //_randomNumberGenerator = FixedValue;
 
    //SimTKOpenMMUtilities::setRandomNumberSeed( randomNumberSeed );
 }   
@@ -473,6 +474,38 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne( void ){
    return DefaultReturnValue;
 }
 
+/** 
+ * Load random number streams w/ fixed value
+ * 
+ *
+ * @return DefaultReturnValue;
+ */
+
+int BrookRandomNumberGenerator::_loadRandomNumberStreamsFixedValue( void ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const std::string methodName  = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsFixedValue";
+   int printOn                          = 1;
+
+   // ---------------------------------------------------------------------------------------
+   
+   // load fixed value
+
+   float fixedValue  = 0.1f;
+   for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
+	   getRandomNumberStream( jj )->fillWithValue( &fixedValue );
+   }
+
+   if( printOn ){
+      FILE* log = getLog() ? getLog() : stderr;
+      (void) fprintf( log, "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
+      (void) fflush( log );
+   }
+
+   return DefaultReturnValue;
+}
+
 /** 
  * Load random number streams using original gpu algorithm
  * 
@@ -595,7 +628,8 @@ int BrookRandomNumberGenerator::_loadGVShuffle( void ){
 	const int np  = sizeof(p) / sizeof(p[0]);
    const int pmax = p[np-1];
 
-   // static const std::string methodName = "\nBrookRandomNumberGenerator::loadGVShuffle";
+   static const std::string methodName  = "\nBrookRandomNumberGenerator::loadGVShuffle";
+   int printOn                          = 0;
 
    // ---------------------------------------------------------------------------------------
 	
@@ -636,6 +670,25 @@ int BrookRandomNumberGenerator::_loadGVShuffle( void ){
 	}
    _getShuffleStream()->loadFromArray( loadBuffer );
 	
+   if( printOn ){
+      FILE* log = getLog() ? getLog() : stderr;
+      (void) fprintf( log, "%s: sz=%d pmax=%d np=%d sample indices:\n", methodName.c_str(), rvSize, pmax, np );
+      float maxIdx = -1.0f;
+      float minIdx =  (float) (rvSize)*2.0f;
+      for( int ii = 0; ii < rvSize; ii++ ){
+         if( ii < 30 || ii > (rvSize-30) ){
+            (void) fprintf( log, "  %d %.8f\n", ii, loadBuffer[ii] );
+         }
+         if( loadBuffer[ii] < minIdx ){
+            minIdx = loadBuffer[ii];
+         }
+         if( loadBuffer[ii] > maxIdx ){
+            maxIdx = loadBuffer[ii];
+         }
+      }
+      (void) fprintf( log, "%s: min-max indices: %.8f %.8f\n", methodName.c_str(), minIdx, maxIdx );
+   }
+      
    return DefaultReturnValue;
 }
 
@@ -685,12 +738,12 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
    // ---------------------------------------------------------------------------------------
 
    static const std::string methodName = "BrookRandomNumberGenerator::advanceGVCursor";
-   int printOn                         = 0;
+   int printOn                         = 1;
    FILE* log;
 
    // ---------------------------------------------------------------------------------------
 	
-//setLog( stderr );
+setLog( stderr );
    if( printOn && getLog() ){
       log = getLog();
    } else {
@@ -728,14 +781,17 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
 			} else { //Need to refresh random numbers from cpu
 
             if( _randomNumberGenerator == Mersenne ){
-               action = "loaded new values from GPU using Mersenne rng";
+               action = "loaded new values to GPU using Mersenne rng";
                _loadRandomNumberStreamsMersenne( );
             } else if( _randomNumberGenerator == Kiss ){
-               action = "loaded new values from GPU using KISS rng";
+               action = "loaded new values to GPU using KISS rng";
                _loadRandomNumberStreamsKiss( );
             } else if( _randomNumberGenerator == Original ){
-               action = "loaded new values from GPU using original Gromac's rng";
+               action = "loaded new values to GPU using original Gromac's rng";
                _loadGVStreamsOriginal( );
+            } else if( _randomNumberGenerator == FixedValue ){
+               action = "loaded new fixed values to GPU";
+               _loadRandomNumberStreamsFixedValue( );
             }
             _numberOfShuffles  = 0;
 			}
@@ -743,7 +799,9 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
 		}
 
       if( printOn ){
-         (void) fprintf( log, "%s StrmIdx=%d action=%s\n", methodName.c_str(), _rvStreamIndex, action );
+         (void) fprintf( log, "%s offset=%d consume/itr=%d StrmSz=%d idx=%d shffle=%d action=%s\n",
+                         methodName.c_str(), _rvStreamOffset, numberOfRandomValuesConsumedPerIteration,
+                         rvStreamSize, _rvStreamIndex, _numberOfShuffles, action );
          (void) fflush( log );
       }
 
@@ -802,7 +860,7 @@ BrookFloatStreamInternal* BrookRandomNumberGenerator::_getShuffleStream( void )
  *
  */
 
-BrookFloatStreamInternal* BrookRandomNumberGenerator::getRandomNumberStream( int index ) const {
+BrookFloatStreamInternal* BrookRandomNumberGenerator::getRandomNumberStream( int index ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -944,16 +1002,20 @@ int BrookRandomNumberGenerator::_initializeStreams( const Platform& platform ){
    _randomNumberGeneratorStreams = new BrookFloatStreamInternal*[_numberOfRandomNumberStreams];
 
    for( int ii = 0; ii < _numberOfRandomNumberStreams; ii++ ){
-      _randomNumberGeneratorStreams[ii]    = new BrookFloatStreamInternal( BrookCommon::ShuffleStream,
+      _randomNumberGeneratorStreams[ii]    = new BrookFloatStreamInternal( BrookCommon::RandomValuesStream,
                                                                            randomNumberStreamSize, randomNumberStreamWidth,
                                                                            BrookStreamInternal::Float3, dangleValue );
    }
 
+   // create shuffle stream
+
+   _loadGVShuffle();
+
    return DefaultReturnValue;
 }
 
 /*  
- * Setup of StochasticDynamics parameters
+ * Setup of streams, ... associated w/ random number generator
  *
  * @param numberOfParticles     number of particles
  * @param platform              Brook platform
@@ -984,6 +1046,8 @@ int BrookRandomNumberGenerator::setup( int numberOfParticles,  const Platform& p
       _loadRandomNumberStreamsKiss( );
    } else if( _randomNumberGenerator == Original ){
       _loadGVStreamsOriginal( );
+   } else if( _randomNumberGenerator == FixedValue ){
+      _loadRandomNumberStreamsFixedValue( );
    }
 
    return DefaultReturnValue;
@@ -1027,6 +1091,8 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
       (void) LOCAL_SPRINTF( value, "%s", "Kiss Rng");
    } else if( _randomNumberGenerator == Original ){
       (void) LOCAL_SPRINTF( value, "%s", "Original Gromacs Rng");
+   } else if( _randomNumberGenerator == FixedValue ){
+      (void) LOCAL_SPRINTF( value, "%s", "Fixed value Rng");
    }
    message << _getLine( tab, "Random number generator:", value ); 
 
@@ -1061,6 +1127,8 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
 
    message << _getLine( tab, "Shuffle:",              (_getShuffleStream()      ? Set : NotSet) ); 
  
+   // show stats
+
    message << getStatisticsString( );
 
    for( int ii = 0; ii < LastStreamIndex; ii++ ){
@@ -1070,6 +1138,12 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
       }
    }
 
+   for( int ii = 0; ii < _numberOfRandomNumberStreams; ii++ ){
+      message << std::endl;
+      if( _randomNumberGeneratorStreams[ii] ){
+         message <<  _randomNumberGeneratorStreams[ii]->getContentsString();
+      }
+   }
 #undef LOCAL_SPRINTF
 
    return message.str();
@@ -1201,29 +1275,30 @@ int BrookRandomNumberGenerator::getStatistics( double statistics[7], int streamI
    for( int ii = 0; ii < _numberOfRandomNumberStreams; ii++ ){
       if( streamIndex < 0 || ii == streamIndex ){
          void* dataArrayV       = _randomNumberGeneratorStreams[ii]->getData( 1 );
-         int streamSize         = _randomNumberGeneratorStreams[ii]->getStreamSize();
-         int index              = 0;
+         int numberOfValues     = _randomNumberGeneratorStreams[ii]->getStreamSize()*_randomNumberGeneratorStreams[ii]->getWidth();
          const float* dataArray = (float*) dataArrayV;
-         for( int ii = 0; ii < streamSize; ii++, index++ ){
+         for( int ii = 0; ii < numberOfValues; ii++ ){
 
-             statistics[0] += dataArray[index];
+             statistics[0] += dataArray[ii];
 
-             double rv2     = dataArray[index]*dataArray[index];
+             double rv2     = dataArray[ii]*dataArray[ii];
              statistics[1] += rv2;
-             statistics[2] += rv2*dataArray[index];
+             statistics[2] += rv2*dataArray[ii];
              statistics[3] += rv2*rv2;
 
-             if( statistics[5] > dataArray[index] ){
-                statistics[5] = dataArray[index];
+             if( statistics[5] > dataArray[ii] ){
+                statistics[5] = dataArray[ii];
              }
-             if( statistics[6] < dataArray[index] ){
-                statistics[6] = dataArray[index];
+             if( statistics[6] < dataArray[ii] ){
+                statistics[6] = dataArray[ii];
              }
          } 
-         statistics[4] += (double) index;
+         statistics[4] += (double) numberOfValues;
       }   
    }   
 
+   // accumulate moments, ... in cumulativeStatistics array
+
    for( int ii = 0; ii < 5; ii++ ){
       cumulativeStatistics[ii] += statistics[ii];
    }

platforms/brook/src/BrookRandomNumberGenerator.h

@@ -46,9 +46,9 @@ class BrookRandomNumberGenerator : public BrookCommon {
 
    public:
   
-      // toggle between original, Mersenne, & Kiss (Nvidia) random generators
+      // toggle between original, Mersenne, & Kiss (Nvidia), fixed value random generators
 
-      enum Rngs { Original, Kiss, Mersenne };
+      enum Rngs { Original, Kiss, Mersenne, FixedValue };
  
       /** 
        * Constructor
@@ -147,7 +147,7 @@ class BrookRandomNumberGenerator : public BrookCommon {
        *
        */
       
-      BrookFloatStreamInternal* getRandomNumberStream( int index ) const;
+      BrookFloatStreamInternal* getRandomNumberStream( int index );
       
       /** 
        * Get random number seed
@@ -381,6 +381,16 @@ class BrookRandomNumberGenerator : public BrookCommon {
       
       int _loadGVStreamsOriginal( void );
       
+      /** 
+       * Load fixed value 'random number' streams using original gpu algorithm
+       * used for diagnostics
+       * 
+       *
+       * @return DefaultReturnValue;
+       */
+      
+      int _loadRandomNumberStreamsFixedValue( void );
+      
       /** 
        * Loads a permutation of indices from 0 to gvSize-1 in
        * sdp->strShuffle. To make sure that the order of the

platforms/brook/src/BrookStreamImpl.cpp

@@ -298,6 +298,6 @@ brook::stream& BrookStreamImpl::getBrookStream( void ){
  * @return Brook stream impl
  */
 
-BrookStreamInternal* BrookStreamImpl::getBrookStreamImpl( void ) const {
+BrookStreamInternal* BrookStreamImpl::getBrookStreamInternal( void ) const {
    return _brookStreamInternal;
 }

platforms/brook/src/BrookStreamImpl.h

@@ -173,7 +173,7 @@ class BrookStreamImpl : public StreamImpl {
        * @return Brook stream impl
        */
 
-      BrookStreamInternal* getBrookStreamImpl( void ) const;
+      BrookStreamInternal* getBrookStreamInternal( void ) const;
 
    protected:
 

platforms/brook/src/BrookStreamInternal.cpp

@@ -448,3 +448,222 @@ std::string BrookStreamInternal::printStatistics( std::string tag, std::vector<s
 
    return message.str();
 }
+
+/* 
+ * Print streams to file
+ *
+ * @param fileName     file name
+ * @param streams      streams to print
+ *
+ * @return DefaultReturnValue
+ *
+ * */
+
+int BrookStreamInternal::printStreamsToFile( std::string fileName, std::vector<BrookStreamInternal*>& streams ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookStreamInternal::printStreamsToFile";
+
+// ---------------------------------------------------------------------------------------
+
+   FILE* filePtr = fopen( fileName.c_str(), "w" );
+   if( !filePtr ){
+      (void) fprintf( stderr, "%s could not open file=<%s>\n", methodName.c_str(), fileName.c_str() );
+      (void) fflush( stderr );
+      return ErrorReturnValue;
+   }
+
+   // gather arrays, widths for eah stream, and set index for each stream
+   // also set minimum of stream sizes
+
+   int minIndex    = 10000000;
+   float** arrays  = new float*[streams.size()];
+   float** sums    = new float*[streams.size()];
+   int* widths     = new int[streams.size()];
+   int* indices    = new int[streams.size()];
+   for( unsigned int ii = 0; ii < streams.size(); ii++ ){
+      BrookStreamInternal* stream  = streams[ii];
+      void* dataArrayV         = stream->getData( 1 );
+      arrays[ii]               = (float*) dataArrayV;
+      widths[ii]               =  stream->getWidth();
+      indices[ii]              = 0;
+      sums[ii]                 = new float[4];
+      sums[ii][0] = sums[ii][1] = sums[ii][2] = sums[ii][3] = 0.0f;
+      if( minIndex > stream->getSize() ){
+         minIndex = stream->getSize();
+      }
+   }
+
+   // sum columns 
+
+   for( int ii = 0; ii < minIndex; ii++ ){
+      for( unsigned int kk = 0; kk < streams.size(); kk++ ){
+         for( int jj = 0; jj < widths[kk]; jj++ ){
+            sums[kk][jj] += arrays[kk][indices[kk]++];
+         }
+      }
+   }
+
+   // reinitialize indices
+
+   for( unsigned int kk = 0; kk < streams.size(); kk++ ){
+      indices[kk] = 0;
+   }
+
+   // show column sums
+
+   for( unsigned int kk = 0; kk < streams.size(); kk++ ){
+      (void) fprintf( filePtr, "Sms " );
+      for( int jj = 0; jj < widths[kk]; jj++ ){
+         (void) fprintf( filePtr, "%15.5e ", sums[kk][jj] );
+      }
+   }
+   (void) fprintf( filePtr, "\n" );
+
+   for( int ii = 0; ii < minIndex; ii++ ){
+      (void) fprintf( filePtr, "%6d ", ii );
+
+      // streams
+
+      for( unsigned int kk = 0; kk < streams.size(); kk++ ){
+
+//         (void) fprintf( filePtr, "[ " );
+ 
+         // ii elements of stream kk
+
+         for( int jj = 0; jj < widths[kk]; jj++ ){
+            (void) fprintf( filePtr, "%15.5e ", arrays[kk][indices[kk]++] );
+         }
+
+ //        (void) fprintf( filePtr, " ]", ii );
+      }
+      (void) fprintf( filePtr, "\n", ii );
+   }
+
+
+   // cleanup
+
+   (void) fclose( filePtr );
+
+   delete[] arrays;
+   delete[] widths;
+   delete[] indices;
+   for( int ii = 0; ii < 4; ii++ ){
+      delete[] sums[ii];
+   }
+   delete[] sums;
+
+   return DefaultReturnValue;
+
+}
+
+typedef struct {
+  unsigned int type;
+  unsigned int dimensions;
+  unsigned int dims[4];
+} STREAM_HEADER;
+
+#include <winsock.h>
+#include <stdarg.h>
+#include <limits>
+#include <set>
+
+int BrookStreamInternal::loadStreamGivenFileName( std::string& filename ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookStreamInternal::loadStreamGivenFileName";
+
+// ---------------------------------------------------------------------------------------
+
+   FILE* log = stderr;
+   
+   /* open file, read header */
+
+   FILE * filePtr = fopen( filename.c_str(), "rb" );
+   if(  filePtr == NULL ){
+      (void) fprintf( log, "%s Unable to open/read %s for stream creation\n", methodName.c_str(), filename.c_str() );
+      return ErrorReturnValue;
+   }
+
+   STREAM_HEADER header;
+   if( 1 != fread( &header, sizeof( header ), 1, filePtr ) ){
+      (void) fprintf( log, "%s Unable to read %s header for stream %s\n", methodName.c_str(), filename.c_str() );
+      (void) fclose( filePtr );
+      return ErrorReturnValue;
+   }
+ 
+   /* fix endian */
+
+//   header.type = ntohl( header.type );
+//  header.dimensions = ntohl( header.dimensions );
+/*
+   for( int ii = 0; ii < 4; ii++ ){
+   //   header.dims[ii] = ntohl( header.dims[ii] );
+      (void) fprintf( log, "%s header %d %d for stream %s from %s\n", methodName.c_str(), ii, header.dims[ii], getName().c_str(), filename.c_str() );
+   }
+
+   (void) fflush( log );
+*/
+/*
+   if( header.dimensions[0]
+     (unsigned int) header.dimensions, (const ::brook::StreamType *) &type,
+     (bool) false );
+*/
+ 
+   /* ok, load in the data */
+
+/*
+   if( getStreamSize()*getWidth() != header.dims[0]*header.dims[1] ){
+      (void) fprintf( log, "%s dimension inconsistency for stream %s from %s dim:[%d %d %d %d] sz=%d != sz=%d containerW=%d\n", 
+                      methodName.c_str(), getName().c_str(), filename.c_str(),
+                      header.dims[0], header.dims[1], header.dims[2], header.dims[3], header.dims[0]*header.dims[1],
+                      getStreamSize(), getWidth() );
+      (void) fclose( filePtr );
+      return ErrorReturnValue;
+   }
+*/
+
+   // always float
+
+   int bytesToRead   = getStreamSize()*getWidth()*sizeof( float );
+   void* dataBuffer  = (void*) malloc( bytesToRead );
+
+   if( dataBuffer == NULL ){
+      (void) fprintf( log, "%s Memory Error for file=%s\n", methodName.c_str(), filename.c_str() );
+      (void) fclose( filePtr );
+      return ErrorReturnValue;
+   }
+
+   size_t bytesRead = fread( dataBuffer, bytesToRead, 1, filePtr );
+   if( bytesRead != 1 ){
+      (void) fprintf( log, "%s Unable to read %d bytes=%d stream from %s\n", methodName.c_str(), bytesRead, bytesToRead, filename.c_str() );
+      free( dataBuffer );
+      (void) fclose( filePtr );
+      return ErrorReturnValue;
+   }
+
+   // float/integer case -- nothing more to do but load; if double, then convert float to double
+
+   if( getBaseDataType() == Float || getBaseDataType() == Integer ){
+      loadFromArray( dataBuffer );
+   } else {
+      double* loadBuffer  = (double*) malloc( bytesToRead*2 );
+      float* readBuffer   = (float*) dataBuffer;
+      for( int ii = 0; ii < getStreamSize()*getWidth(); ii++ ){
+         loadBuffer[ii] = (double) readBuffer[ii];
+      }
+      loadFromArray( loadBuffer );
+      free( loadBuffer );
+   }
+         
+   (void) fclose( filePtr );
+   free( dataBuffer );
+ 
+   (void) fprintf( log, "%s read %d bytes for stream %s from %s dim:[%d %d %d %d] container=%d %d\n", methodName.c_str(), bytesToRead, getName().c_str(), filename.c_str(),
+                   header.dims[0], header.dims[0], header.dims[0], header.dims[0], getStreamSize(), getWidth() );
+   (void) fflush( log );
+
+   return DefaultReturnValue;
+}

platforms/brook/src/BrookStreamInternal.h

@@ -263,10 +263,33 @@ class BrookStreamInternal {
        *
        * @return stat string
        *
-       * */
+       **/
       
       std::string printStatistics( std::string tag, std::vector<std::vector<double> >& statistics ) const;
       
+      /* 
+       * Read stream from file
+       *
+       * @param fileName     file name
+       *
+       * @return DefaultReturnValue or ErrorReturnValue if problems
+       *
+       **/
+      
+      int loadStreamGivenFileName( std::string& filename );
+      
+      /* 
+       * Print streams to file
+       *
+       * @param fileName     file name
+       * @param streams      streams to print
+       *
+       * @return DefaultReturnValue
+       *
+       **/
+       
+      static int printStreamsToFile( std::string fileName, std::vector<BrookStreamInternal*>& streams );
+      
    protected:
 
       std::string _name;

platforms/brook/src/BrookVelocityCenterOfMassRemoval.cpp

@@ -144,7 +144,7 @@ BrookFloatStreamInternal* BrookVelocityCenterOfMassRemoval::getLinearMomentumStr
  *
  */
    
-int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl& vStream, BrookOpenMMFloat velocityCom[3] ){
+int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl& vStream, BrookOpenMMFloat velocityCom[3], BrookOpenMMFloat* ke ){
 
    // ---------------------------------------------------------------------------------------
 
@@ -157,7 +157,8 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl&
    // calculate linear momentum via reduction
    // subtract it (/totalMass) from velocities
 
-   BrookStreamInternal* velocityStream  = dynamic_cast<BrookStreamInternal*> (vStream.getBrookStreamImpl());
+   BrookStreamInternal* velocityStream  = dynamic_cast<BrookStreamInternal*> (vStream.getBrookStreamInternal());
+   *ke                                  = static_cast<BrookOpenMMFloat>(0.0);
 
    void* velV                           = velocityStream->getData( 1 );
    const float* vArray                  = (float*) velV;
@@ -174,11 +175,14 @@ int BrookVelocityCenterOfMassRemoval::getVelocityCenterOfMass( BrookStreamImpl&
       velocityCom[0] += mArray[ii]*vArray[index];
       velocityCom[1] += mArray[ii]*vArray[index+1];
       velocityCom[2] += mArray[ii]*vArray[index+2];
-if( ii < 10 ){
-fprintf( stderr, "getVelocityCenterOfMass %d %d m=%15.6e v[%15.6e %15.6e %15.6e ]\n", ii, index, mArray[ii], vArray[index], vArray[index+1], vArray[index+2] );
+                 *ke += mArray[ii]*( vArray[index]*vArray[index] + vArray[index+1]*vArray[index+1] + vArray[index+2]*vArray[index+2] );
+if( ii < 3 ){
+fprintf( stderr, "getVelocityCenterOfMass %5d %5d m=%15.6e v[%15.6e %15.6e %15.6e ] %d\n", ii, index, mArray[ii], vArray[index], vArray[index+1], vArray[index+2], numberOfParticles );
 }
    }
 
+   *ke *= static_cast<BrookOpenMMFloat>(0.5);
+
    return DefaultReturnValue;
 }
 
@@ -316,7 +320,6 @@ int BrookVelocityCenterOfMassRemoval::_setMasses( const std::vector<double>& mas
          BrookOpenMMFloat value    = static_cast<BrookOpenMMFloat>(*ii);
          localMasses[index]        = value;
          _totalInverseMass        += value;
-fprintf( stderr, "%s %d mass=%.3f\n", methodName.c_str(), index, value );
       }
    }
 
@@ -455,15 +458,19 @@ int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImp
 
    static const std::string methodName  = "BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
    static const int debug               = 1;
+   static const double boltz            = 8.3145112119486e-03;
+   BrookOpenMMFloat ke;
 
    // ---------------------------------------------------------------------------------------
 
 setLog( stderr );
    if( debug && getLog() ){
       BrookOpenMMFloat com[3];
-      getVelocityCenterOfMass( velocityStream, com );
-      (void) fprintf( getLog(), "%s Pre removal com: [%12.5e %12.5e %12.5e]\n", methodName.c_str(), com[0], com[1], com[2] );
-      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamImpl());
+      getVelocityCenterOfMass( velocityStream, com, &ke );
+      BrookOpenMMFloat denomiator = (( (float) 3*getNumberOfParticles()))*boltz;
+      (void) fprintf( getLog(), "%s Pre removal com: [%12.5e %12.5e %12.5e] ke=%.3f ~T=%.3f\n", 
+                      methodName.c_str(), com[0], com[1], com[2], ke, ke/denomiator );
+      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamInternal());
 
       void* velV                       = outputVelocityStream->getData( 1 );
       const float* vArray              = (float*) velV;
@@ -483,20 +490,22 @@ setLog( stderr );
    kCalculateLinearMomentum( getMassStream()->getBrookStream(), velocityStream.getBrookStream(), getWorkStream()->getBrookStream() );
    kSumLinearMomentum( (float) getComParticleStreamWidth(), (float) getNumberOfParticles(), (float) getTotalInverseMass(),
                         getWorkStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
-//   kScale( (float) getTotalInverseMass(), getLinearMomentumStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
    kRemoveLinearMomentum( getLinearMomentumStream()->getBrookStream(), velocityStream.getBrookStream(), velocityStream.getBrookStream() );
 
    if( (0 || debug) && getLog() ){
       BrookOpenMMFloat com[3];
-      getVelocityCenterOfMass( velocityStream, com );
-      (void) fprintf( getLog(), "%s strW=%d iatm=%d invMass=%.4e Post removal com: [%12.5e %12.5e %12.5e]", methodName.c_str(),
-                      getComParticleStreamWidth(), getNumberOfParticles(), getTotalInverseMass(),  com[0], com[1], com[2] );
+      getVelocityCenterOfMass( velocityStream, com, &ke );
+      BrookOpenMMFloat denomiator = (( (float) 3*getNumberOfParticles()))*boltz;
+      (void) fprintf( getLog(), "%s strW=%d iatm=%d invMass=%.4e\n   Post removal com: [%12.5e %12.5e %12.5e]  ke=%.3f ~T=%.3f\n",
+                      methodName.c_str(),
+                      getComParticleStreamWidth(), getNumberOfParticles(), getTotalInverseMass(),  com[0], com[1], com[2],
+                      ke, ke/denomiator );
 
       void* linMoV = getLinearMomentumStream()->getData( 1 );
       float* linMo = (float*) linMoV;
-      (void) fprintf( getLog(), "LM [%12.5e %12.5e %12.5e]\n", linMo[0], linMo[1], linMo[2] );
+      (void) fprintf( getLog(), "   LM [%12.5e %12.5e %12.5e]\n", linMo[0], linMo[1], linMo[2] );
 
-      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamImpl());
+      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamInternal());
 
       void* velV                       = outputVelocityStream->getData( 1 );
       const float* vArray              = (float*) velV;

platforms/brook/src/BrookVelocityCenterOfMassRemoval.h

@@ -102,16 +102,17 @@ class BrookVelocityCenterOfMassRemoval : public BrookCommon {
       int removeVelocityCenterOfMass( BrookStreamImpl& velocityStream );
 
       /** 
-       * Get velocity center-of-mass
+       * Get velocity center-of-mass and kinetic energy (used for diagnostics)
        * 
        * @param  velocities                  particle velocities
        * @param  velocityCom                 output velocity com
+       * @param  ke                          output kinetic energy
        *
        * @return  DefaultReturnValue
        *
        */
       
-      int getVelocityCenterOfMass( BrookStreamImpl& vStream, BrookOpenMMFloat velocityCom[3] );
+      int getVelocityCenterOfMass( BrookStreamImpl& vStream, BrookOpenMMFloat velocityCom[3], BrookOpenMMFloat* ke );
 
       /* 
        * Setup of  parameters

platforms/brook/src/BrookVerletDynamics.cpp

@@ -519,12 +519,12 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
          (void) fprintf( log, "\nInverseMassStream %d\n", _internalStepCount );
          getInverseMassStream()->printToFile( log );
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
          brookStreamInternalPos->printToFile( log );
    
          (void) fprintf( log, "\nForceStream %d\n", _internalStepCount );
-         BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
          std::vector<std::vector<double> > forceStatistics;
          brookStreamInternalF->getStatistics( forceStatistics, getNumberOfParticles() );
 
@@ -535,7 +535,7 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
 
          brookStreamInternalF->printToFile( log );
    
-         BrookStreamInternal* brookStreamInternalV   = velocityStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalV   = velocityStream.getBrookStreamInternal();
          std::vector<std::vector<double> > velocityStatistics;
          brookStreamInternalV->getStatistics( velocityStatistics, getNumberOfParticles() );
          std::stringstream tagV;
@@ -571,7 +571,7 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
          (void) fprintf( log, "\nShakeParticleIndicesStream %d\n", _internalStepCount );
          brookShakeAlgorithm.getShakeParticleIndicesStream()->printToFile( log );
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
          brookStreamInternalPos->printToFile( log );
    
@@ -616,7 +616,7 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
          (void) fprintf( log, "\nShakeInverseMapStream %d\n", _internalStepCount );
          brookShakeAlgorithm.getShakeInverseMapStream()->printToFile( log );
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
          brookStreamInternalPos->printToFile( log );
    
@@ -651,7 +651,7 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
          (void) fprintf( log, "\n%s step=%d Post kupdate_md2: inverseStepSize=%3e",
                                    methodName.c_str(), _internalStepCount, inverseStepSize );
    
-         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamInternal();
          brookShakeAlgorithm.checkConstraints( brookStreamInternalPos, log, 0.0001f );
          (void) fprintf( log, "\nPositionStream %d\n", _internalStepCount );
          brookStreamInternalPos->printToFile( log );
@@ -659,7 +659,7 @@ int BrookVerletDynamics::update( BrookStreamImpl& positionStream, BrookStreamImp
          (void) fprintf( log, "\nXPrimeStream %d\n", _internalStepCount );
          getXPrimeStream()->printToFile( log ); 
 
-         brookStreamInternalPos = velocityStream.getBrookStreamImpl();
+         brookStreamInternalPos = velocityStream.getBrookStreamInternal();
          (void) fprintf( log, "\nVelocityStream %d\n", _internalStepCount );
          brookStreamInternalPos->printToFile( log ); 
       }   

platforms/brook/src/OpenMMBrookInterface.cpp

@@ -591,7 +591,7 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
    
    /*
          (void) fprintf( getLog(), "\nFinal NB & bonded forces" );
-         BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+         BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamInternal();
          brookStreamInternalF->printToFile( getLog() );
          void* dataV = brookStreamInternalF->getData(1);
          float* data = (float*) dataV;
@@ -636,7 +636,7 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
 
    // first step only?
 
-   if( step > 20 ){
+   if( step > 0 ){
       return;
    }
 
@@ -650,9 +650,9 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
    BrookStreamImpl* positions = getParticlePositions();
    BrookStreamImpl* forces    = getParticleForces();
 
-   std::vector<BrookStreamImpl*> streams;
-   streams.push_back( positions );
-   streams.push_back( forces );
+   std::vector<BrookStreamInternal*> streams;
+   streams.push_back( positions->getBrookStreamInternal() );
+   streams.push_back( forces->getBrookStreamInternal() );
 
 // ---------------------------------------------------------------------------------------
 
@@ -662,7 +662,7 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
       forces->fillWithValue( &zero );
       _brookNonBonded.computeForces( *positions, *forces );
       std::string fileName = fileNameBase + "NonBonded.txt";
-      printStreamsToFile( fileName, streams );
+      BrookStreamInternal::printStreamsToFile( fileName, streams );
    }
 
 // ---------------------------------------------------------------------------------------
@@ -684,7 +684,7 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
       forces->fillWithValue( &zero );
       _brookBonded.computeForces( *positions, *forces );
       std::string fileName = fileNameBase + "Bonded.txt";
-      printStreamsToFile( fileName, streams );
+      BrookStreamInternal::printStreamsToFile( fileName, streams );
    
    }
 
@@ -699,7 +699,7 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
       _brookGbsa.computeForces( *positions, *forces );
 
       std::string fileName = fileNameBase + "Obc.txt";
-      printStreamsToFile( fileName, streams );
+      BrookStreamInternal::printStreamsToFile( fileName, streams );
 
    }
    forces->fillWithValue( &zero );
@@ -723,7 +723,7 @@ void OpenMMBrookInterface::printForcesToFile( OpenMMContextImpl& context ){
       }
 
       std::string fileName = fileNameBase + "AllF.txt";
-      printStreamsToFile( fileName, streams );
+      BrookStreamInternal::printStreamsToFile( fileName, streams );
    }
 
 (void) fprintf( stderr, "%s done computeForces\n", methodName.c_str() );
@@ -770,114 +770,3 @@ double OpenMMBrookInterface::computeEnergy( OpenMMContextImpl& context, System&
    return refContext.getState(State::Energy).getPotentialEnergy();
 
 }
-
-/* 
- * Print contents of object to file
- *
- * @param fileName     file name
- * @param streams      streams to print
- *
- * @return DefaultReturnValue
- *
- * */
-
-int OpenMMBrookInterface::printStreamsToFile( std::string fileName, std::vector<BrookStreamImpl*>& streams ){
-
-// ---------------------------------------------------------------------------------------
-
-   static const std::string methodName      = "OpenMMBrookInterface::printStreamsToFile";
-
-// ---------------------------------------------------------------------------------------
-
-   FILE* filePtr = fopen( fileName.c_str(), "w" );
-   if( !filePtr ){
-      (void) fprintf( stderr, "%s coud not open file=<%s>\n", methodName.c_str(), fileName.c_str() );
-      (void) fflush( stderr );
-      return ErrorReturnValue;
-   }
-
-   // gather arrays, widths for eah stream, and set index for each stream
-   // also set minimum of stream sizes
-
-   int minIndex    = 10000000;
-   float** arrays  = new float*[streams.size()];
-   float** sums    = new float*[streams.size()];
-   int* widths     = new int[streams.size()];
-   int* indices    = new int[streams.size()];
-   for( unsigned int ii = 0; ii < streams.size(); ii++ ){
-      BrookStreamImpl* stream  = streams[ii];
-      void* dataArrayV         = stream->getData( 1 );
-      arrays[ii]               = (float*) dataArrayV;
-      widths[ii]               =  stream->getWidth();
-      indices[ii]              = 0;
-      sums[ii]                 = new float[4];
-      sums[ii][0] = sums[ii][1] = sums[ii][2] = sums[ii][3] = 0.0f;
-      if( minIndex > stream->getStreamSize() ){
-         minIndex = stream->getStreamSize();
-      }
-   }
-   minIndex = minIndex > getNumberOfParticles() ? getNumberOfParticles() : minIndex;
-
-   // sum columns 
-
-   for( int ii = 0; ii < minIndex; ii++ ){
-      for( unsigned int kk = 0; kk < streams.size(); kk++ ){
-         for( int jj = 0; jj < widths[kk]; jj++ ){
-            sums[kk][jj] += arrays[kk][indices[kk]++];
-         }
-      }
-   }
-
-   // reinitialize indices
-
-   for( unsigned int kk = 0; kk < streams.size(); kk++ ){
-      indices[kk] = 0;
-   }
-
-   // show column sums
-
-   for( unsigned int kk = 0; kk < streams.size(); kk++ ){
-      (void) fprintf( filePtr, "Sms " );
-      for( int jj = 0; jj < widths[kk]; jj++ ){
-         (void) fprintf( filePtr, "%12.5e ", sums[kk][jj] );
-      }
-   }
-   (void) fprintf( filePtr, "\n" );
-
-   for( int ii = 0; ii < minIndex; ii++ ){
-      (void) fprintf( filePtr, "%d ", ii );
-
-      // streams
-
-      for( unsigned int kk = 0; kk < streams.size(); kk++ ){
-
-//         (void) fprintf( filePtr, "[ " );
- 
-         // ii elements of stream kk
-
-         for( int jj = 0; jj < widths[kk]; jj++ ){
-            (void) fprintf( filePtr, "%12.5e ", arrays[kk][indices[kk]++] );
-         }
-
- //        (void) fprintf( filePtr, " ]", ii );
-      }
-      (void) fprintf( filePtr, "\n", ii );
-   }
-
-
-   // cleanup
-
-   (void) fclose( filePtr );
-
-   delete[] arrays;
-   delete[] widths;
-   delete[] indices;
-   for( int ii = 0; ii < 4; ii++ ){
-      delete[] sums[ii];
-   }
-   delete[] sums;
-
-   return DefaultReturnValue;
-
-}
-

platforms/brook/src/OpenMMBrookInterface.h

@@ -390,18 +390,6 @@ class OpenMMBrookInterface {
       
       void printForcesToFile( OpenMMContextImpl& context );
       
-      /* 
-       * Print contents of object to file
-       *
-       * @param fileName     file name
-       * @param streams      streams to print
-       *
-       * @return DefaultReturnValue
-       *
-       * */
-      
-      int printStreamsToFile( std::string fileName, std::vector<BrookStreamImpl*>& streams );
-      
    private:
    
       static const int DefaultReturnValue = 0;