Mods

platforms/brook/src/BrookFloatStreamInternal.cpp

@@ -242,6 +242,12 @@ void BrookFloatStreamInternal::loadFromArray( const void* array, BrookStreamInte
    int totalSize                           = getSize()*getWidth();
    //int totalSize                           = getSize();
 
+/*
+fprintf( stderr, "%s %s data=%p stream=%d [%d %d] width=%d type=%d\n",
+          methodName.c_str(), getName().c_str(), _data, getStreamSize(), _streamHeight, _streamWidth, _width, baseType );
+fflush( stderr );
+*/
+
    if( baseType == BrookStreamInternal::Float ){
 
       memcpy( _data, array, sizeof( float )*totalSize );

platforms/brook/src/BrookLangevinDynamics.cpp

@@ -249,7 +249,7 @@ int BrookLangevinDynamics::_updateDerivedParameters( void ){
    BrookOpenMMFloat temperature = getTemperature();
    BrookOpenMMFloat stepSize    = getStepSize();
 
-   if(  fabsf( float( tau ) ) < epsilon ){
+   if( fabsf( (float) tau ) < epsilon ){
       std::stringstream message;
       message << methodName << " tau=" << tau << " too small.";
       throw OpenMMException( message.str() );
@@ -809,18 +809,18 @@ int BrookLangevinDynamics::setup( const std::vector<double>& masses, const Platf
    return DefaultReturnValue;
 }
 
-/*  
- * Setup of LangevinDynamics parameters
- *
- * @param masses                masses
- * @param platform              Brook platform
+/** 
+ * Get T
  * 
- * @return nonzero value if error
+ * @param velocities             velocities
+ * @param inverseMassStream      inverse masses
+ * @param numberOfConstraints    number of constraints
  *
- * */
+ * @return   temperature
+ */
 
 float BrookLangevinDynamics::getTemperature( BrookStreamInternal* velocities, BrookFloatStreamInternal* inverseMassStream,
-                                             BrookShakeAlgorithm& brookShakeAlgorithm ) const {
+                                             int numberOfConstraints ) const {
     
 // ---------------------------------------------------------------------------------------
 
@@ -839,14 +839,95 @@ float BrookLangevinDynamics::getTemperature( BrookStreamInternal* velocities, Br
 
    for( int ii = 0; ii < getNumberOfParticles(); ii++ ){
       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;
    }
 
-   int degreesOfFreedom = 3*getNumberOfParticles() - brookShakeAlgorithm.getNumberOfConstraints();
+   int degreesOfFreedom = 3*getNumberOfParticles() - numberOfConstraints;
 
+(void) fprintf( stderr, "%s ke=%.5e dof=%d\n", methodName.c_str(), ke, degreesOfFreedom ); 
    ke /= ((float) BOLTZ)*((float) ( degreesOfFreedom ));
 
-   return ke;
+   return 0.5f*ke;
+}
+
+/** 
+ * Remove velocity com
+ * 
+ * @param velocities             velocities
+ * @param inverseMassStream      inverse masses
+ *
+ * @return DefaultReturnValue  
+ */
+
+int BrookLangevinDynamics::removeCom( BrookStreamInternal* velocities, BrookFloatStreamInternal* inverseMassStream ) const {
+    
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookLangevinDynamics::removeCom";
+
+// ---------------------------------------------------------------------------------------
+
+ //  (void) fprintf( stderr, "%s\n", methodName.c_str() ); fflush( stderr );
+
+   void* dataArrayV           = velocities->getData( 1 );
+   float* velocitiesI         = (float*) dataArrayV;
+
+   void* inverseMassStreamV   = inverseMassStream->getData( 1 );
+   float* inverseMassStreamI  = (float*) inverseMassStreamV;
+
+   float totalMass            = 0.0f;
+   float com[3]               = { 0.0f, 0.0f, 0.0f };
+   int index                  = 0;
+
+//(void) fprintf( stderr, "%s strm %d %d\n", methodName.c_str(), velocities->getStreamSize(),velocities->getWidth() ); fflush( stderr );
+   for( int ii = 0; ii < getNumberOfParticles(); ii++ ){
+      float mass   = 1.0f/inverseMassStreamI[ii];
+      totalMass   += mass;
+      com[0]      += mass*velocitiesI[index];
+      com[1]      += mass*velocitiesI[index+1];
+      com[2]      += mass*velocitiesI[index+2];
+      index       += 3;
+   }
+   totalMass   = 1.0f/totalMass;
+   com[0]     *= totalMass;
+   com[1]     *= totalMass;
+   com[2]     *= totalMass;
+
+   index                = 0;
+   double* newVelocities = new double[velocities->getStreamSize()*velocities->getWidth()];
+   memset( newVelocities, 0, sizeof( double )*velocities->getStreamSize()*velocities->getWidth() );
+   for( int ii = 0; ii < getNumberOfParticles(); ii++ ){
+      newVelocities[index]    = (double) velocitiesI[index]   - com[0];
+      newVelocities[index+1]  = (double) velocitiesI[index+1] - com[1];
+      newVelocities[index+2]  = (double) velocitiesI[index+2] - com[2];
+      index                  += 3;
+   }
+
+/*
+   (void) fprintf( stderr, "%s com[%14.5e %14.5e %14.5e]\n", methodName.c_str(), com[0], com[1], com[2] ); 
+   for( int ii = 0; ii < velocities->getStreamSize()*3; ii += 3 ){
+      (void) fprintf( stderr, "%s %d newVelocities[%14.5e %14.5e %14.5e]\n", methodName.c_str(), ii/3, newVelocities[ii], newVelocities[ii+1], newVelocities[ii+2] ); 
+   }
+*/
+
+   velocities->loadFromArray( newVelocities );
+
+   dataArrayV           = velocities->getData( 1 );
+   velocitiesI         = (float*) dataArrayV;
+
+/*
+   (void) fprintf( stderr, "%s readback\n", methodName.c_str() );
+   for( int ii = 0; ii < velocities->getStreamSize()*3; ii += 3 ){
+      (void) fprintf( stderr, "%s %d velocitiesI[%14.5e %14.5e %14.5e]\n", methodName.c_str(), ii/3, velocitiesI[ii], velocitiesI[ii+1], velocitiesI[ii+2] ); 
+   }
+*/
+
+   velocities->loadFromArray( newVelocities );
+
+   delete[] newVelocities;
+
+   return DefaultReturnValue;
 }
 
 /* 
@@ -950,21 +1031,25 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName = "\nBrookLangevinDynamics::update";
-   int printOn                         = 0;
+   int printOn                         = 1;
    FILE* log;
 
 // ---------------------------------------------------------------------------------------
 
    _internalStepCount++;
 
-//setLog( stderr );
-   printOn = (printOn && getLog()) ? printOn : 0;
+// setLog( stderr );
+
+   if( printOn && getLog() ){
+      log = getLog();
+   } else { 
+      printOn = 0;
+   }
 
    const BrookOpenMMFloat* derivedParameters           = getDerivedParameters();
 
    if( printOn ){
 
-      log = getLog();
       static int showAux = 1;
 
       if( printOn ){
@@ -977,7 +1062,10 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
       if( showAux ){
          showAux = 0;
 
-         std::string contents             = brookRandomNumberGenerator.getContentsString( );
+         std::string contents             = getContentsString( );
+         (void) fprintf( log, "%s step=%d contents\n%s", methodName.c_str(), _internalStepCount, contents.c_str() );
+
+         contents             = brookRandomNumberGenerator.getContentsString( );
          (void) fprintf( log, "%s step=%d RNG contents\n%s", methodName.c_str(), _internalStepCount, contents.c_str() );
          // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( log ); 
 
@@ -987,6 +1075,15 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
       }
    }
 
+   if( 1 || (_internalStepCount % 10) == 0 ){
+      FILE*	log1 = stderr;
+      BrookStreamInternal* brookStreamInternalPos  = velocityStream.getBrookStreamImpl();
+      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 );
+   }
+
    // first integration step
 
    kupdate_sd1_fix1(
@@ -1011,7 +1108,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
  
    // diagnostics
 
-   if( 0 && printOn ){
+   if( printOn ){
       (void) fprintf( log, "\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(),
@@ -1143,7 +1240,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
 
    // diagnostics
 
-   if( 0 && printOn ){
+   if( printOn ){
       (void) fprintf( log, "\n%s step=%d Post kupdate_sd2_fix1: particleStrW=%3d rngStrW=%3d rngOff=%5d "
                                 "Sd2pc[]=[%12.5e %12.5e]", methodName.c_str(), _internalStepCount,
                                 getLangevinDynamicsParticleStreamWidth(),
@@ -1167,7 +1264,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
       (void) fprintf( log, "\nXPrimeStream %d \n", _internalStepCount );
       getXPrimeStream()->printToFile( log ); 
 
-      (void) fprintf( log, "\ngetSD2XStream %d \n", _internalStepCount );
+      (void) fprintf( log, "\nSD2XStream %d \n", _internalStepCount );
       getSD2XStream()->printToFile( log );
 
       BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamImpl();
@@ -1244,7 +1341,7 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
          tagV << _internalStepCount << " Vxx "; 
          stats = brookStreamInternalPos->printStatistics( tagV.str(), velocityStatistics );
          (void) fprintf( log, "\nStep %d Velocity stats:\n%s", _internalStepCount, stats.c_str() );
-         float temperature = getTemperature( brookStreamInternalPos, getInverseMassStream(), brookShakeAlgorithm );
+         float temperature = getTemperature( brookStreamInternalPos, getInverseMassStream(), brookShakeAlgorithm.getNumberOfConstraints() );
          (void) fprintf( log, "\nVelocityStream %d T=%.3f\n", _internalStepCount, temperature );
          brookStreamInternalPos->printToFile( log );
    
@@ -1287,6 +1384,20 @@ int BrookLangevinDynamics::update( BrookStreamImpl& positionStream, BrookStreamI
       }   
    }
 
+   if( (_internalStepCount % 10) == 0 ){
+      FILE*	log1 = stderr;
+      BrookStreamInternal* brookStreamInternalPos  = velocityStream.getBrookStreamImpl();
+      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 );
+   }
+
    return DefaultReturnValue;
 
 }

platforms/brook/src/BrookLangevinDynamics.h

@@ -258,16 +258,30 @@ class BrookLangevinDynamics : public BrookCommon {
       
       BrookFloatStreamInternal* getInverseMassStream( void ) const;
 
-
       /** 
-       * Get T
+       * Get Temperature
+       * 
+       * @param velocities             velocities
+       * @param inverseMassStream      inverse masses
+       * @param numberOfConstraints    number of constraints
        *
-       * @return   T
+       * @return temperature
        *
        */
       
       float getTemperature( BrookStreamInternal* velocities, BrookFloatStreamInternal* inverseMassStream, 
-                            BrookShakeAlgorithm& brookShakeAlgorithm ) const;
+                            int numberOfConstraints ) const;
+
+      /** 
+       * Remove velocity com
+       * 
+       * @param velocities             velocities
+       * @param inverseMassStream      inverse masses
+       *
+       * @return DefaultReturnValue  
+       */
+      
+      int removeCom( BrookStreamInternal* velocities, BrookFloatStreamInternal* inverseMassStream ) const;
       
    private:
    

platforms/brook/src/BrookRandomNumberGenerator.cpp

@@ -426,6 +426,7 @@ state[3] = 27587;
    }
 
    if( printOn ){
+      FILE* log = getLog() ? getLog() : stderr;
       (void) fprintf( log, "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
       (void) fflush( log );
    }
@@ -447,7 +448,6 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne( void ){
 
    static const std::string methodName  = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne";
    int printOn                          = 0;
-   FILE* log;
 
    // ---------------------------------------------------------------------------------------
    
@@ -465,6 +465,7 @@ int BrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne( void ){
    }
 
    if( printOn ){
+      FILE* log = getLog() ? getLog() : stderr;
       (void) fprintf( log, "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
       (void) fflush( log );
    }
@@ -689,6 +690,7 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
 
    // ---------------------------------------------------------------------------------------
 	
+//setLog( stderr );
    if( printOn && getLog() ){
       log = getLog();
    } else {
@@ -703,10 +705,9 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
 
 	//Check if we've used up this texture
 
+   char* action = "none";
 	if ( _rvStreamOffset > rvStreamSize - numberOfRandomValuesConsumedPerIteration ){
 
-      char* action;
-
 		// next one if available
 
       _rvStreamOffset = 0;
@@ -742,9 +743,27 @@ int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsume
 		}
 
       if( printOn ){
-         (void) fprintf( log, "%s StrmIdx=%d action=%s\n", methodName, _rvStreamIndex, action );
+         (void) fprintf( log, "%s StrmIdx=%d action=%s\n", methodName.c_str(), _rvStreamIndex, action );
          (void) fflush( log );
       }
+
+/*
+      // check rng distribution
+
+      static int count = 0;
+      if( count++ < 2 ){
+
+         // accumulate rng -- stats will be in cumulative fields in stat string
+
+         int testIterations = 1000;
+         for( int ii = 0; ii < testIterations; ii++ ){
+            //_loadRandomNumberStreamsKiss( );
+            _loadRandomNumberStreamsMersenne( );
+            getStatisticsString();
+         }
+         (void) fprintf( log, "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
+      }
+*/
 	}
 
    return DefaultReturnValue;

platforms/brook/src/BrookShakeAlgorithm.cpp

@@ -463,6 +463,7 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
    
    particleIterator                                     = constraintIndices.begin();
    std::vector<double>::const_iterator distanceIterator = constraintLengths.begin();
+   _numberOfConstraints                                 = static_cast<int>(constraintIndices.size());
    while( particleIterator != constraintIndices.end() ){
 
       float distance                  = static_cast<float>( *distanceIterator );
@@ -475,7 +476,7 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
       bool firstIsCentral;
       if( constraintCount[atomI] > 1 ){
          firstIsCentral = true;
-      } else if (constraintCount[atomJ] > 1 ){
+      } else if( constraintCount[atomJ] > 1 ){
          firstIsCentral = false;
       } else if( atomI < atomJ ){
          firstIsCentral = true;
@@ -531,7 +532,6 @@ int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, co
    // load indices & parameters
 
    int constraintIndex                   = 0;
-   _numberOfConstraints                  = static_cast<int>(_clusters.size());
    for( map<int, ShakeCluster>::const_iterator iter = _clusters.begin(); iter != _clusters.end(); ++iter ){
 
       const ShakeCluster& cluster        = iter->second;

platforms/brook/src/BrookVelocityCenterOfMassRemoval.cpp

@@ -134,87 +134,6 @@ BrookFloatStreamInternal* BrookVelocityCenterOfMassRemoval::getLinearMomentumStr
    return _streams[LinearMomentumStream];
 }
 
-/**
- * Remove velocity-COM
- *
- * @param velocities velocities
- *
- * @return DefaultReturnValue
- *
- */
-   
-int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImpl& velocityStream ){
-
-   // ---------------------------------------------------------------------------------------
-
-   static const std::string methodName  = "BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
-   static const int debug               = 0;
-
-   // ---------------------------------------------------------------------------------------
-
-   if( debug && getLog() ){
-      BrookOpenMMFloat com[3];
-      getVelocityCenterOfMass( velocityStream, com );
-      (void) fprintf( getLog(), "\n%s Pre removal com: [%12.5e %12.5e %12.5e]\n", methodName, com[0], com[1], com[2] );
-      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamImpl());
-
-      void* velV                       = outputVelocityStream->getData( 1 );
-      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() );
-   }
-
-   // calculate linear momentum via reduction
-   // subtract it (/totalMass) from velocities
-
-   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,
-                      getComParticleStreamWidth(), getNumberOfParticles(), getTotalInverseMass(),  com[0], com[1], com[2] );
-
-      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] );
-
-      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamImpl());
-
-      void* velV                       = outputVelocityStream->getData( 1 );
-      const float* vArray              = (float*) velV;
-
-      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() );
-
-//exit(0);
-   }
-
-   return DefaultReturnValue;
-}
-
 /**
  * Get velocity-COM
  *
@@ -255,6 +174,9 @@ 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] );
+}
    }
 
    return DefaultReturnValue;
@@ -394,6 +316,7 @@ 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 );
       }
    }
 
@@ -443,10 +366,12 @@ int BrookVelocityCenterOfMassRemoval::setup( const std::vector<double>& masses,
 
    _setMasses( masses );
 
-   if( 1 && getLog() ){
+   //if( 1 && getLog() ){
+   if( 1 ){
+FILE* log = stderr;
       std::string contents = getContentsString( 0 );
-      (void) fprintf( getLog(), "%s contents:\n%s\n", methodName.c_str(), contents.c_str() );
-      (void) fflush( getLog() );
+      (void) fprintf( log, "%s contents:\n%s\n", methodName.c_str(), contents.c_str() );
+      (void) fflush( log );
    }
    
    return DefaultReturnValue;
@@ -495,7 +420,7 @@ std::string BrookVelocityCenterOfMassRemoval::getContentsString( int level ) con
    (void) LOCAL_SPRINTF( value, "%d", getComParticleStreamSize() );
    message << _getLine( tab, "Particle stream size:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%.5f", getTotalInverseMass() );
+   (void) LOCAL_SPRINTF( value, "%.5e", getTotalInverseMass() );
    message << _getLine( tab, "TotalInverseMass:", value ); 
 
    message << _getLine( tab, "Log:",                  (getLog()                  ? Set : NotSet) ); 
@@ -515,3 +440,84 @@ std::string BrookVelocityCenterOfMassRemoval::getContentsString( int level ) con
 
    return message.str();
 }
+/**
+ * Remove velocity-COM
+ *
+ * @param velocities velocities
+ *
+ * @return DefaultReturnValue
+ *
+ */
+   
+int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( BrookStreamImpl& velocityStream ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const std::string methodName  = "BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
+   static const int debug               = 1;
+
+   // ---------------------------------------------------------------------------------------
+
+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());
+
+      void* velV                       = outputVelocityStream->getData( 1 );
+      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() );
+   }
+
+   // calculate linear momentum via reduction
+   // subtract it (/totalMass) from velocities
+
+   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] );
+
+      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] );
+
+      BrookStreamInternal* outputVelocityStream = dynamic_cast<BrookStreamInternal*> (velocityStream.getBrookStreamImpl());
+
+      void* velV                       = outputVelocityStream->getData( 1 );
+      const float* vArray              = (float*) velV;
+
+      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() );
+
+//exit(0);
+   }
+
+   return DefaultReturnValue;
+}
+

platforms/brook/src/OpenMMBrookInterface.cpp

@@ -546,7 +546,8 @@ void OpenMMBrookInterface::computeForces( OpenMMContextImpl& context ){
 
    // info
 
-   if( printOn > 1 ){
+   //if( printOn > 1 ){
+   if( 1 ){
       printForcesToFile( context );
    }
 

platforms/brook/src/gpu/kshakeh.h

@@ -52,8 +52,8 @@ void  kshakeh_update1_fix1 (
 		::brook::stream cposq3,
 		::brook::stream oposq); 
 
-
-void  kshakeh_update2_fix1 (const float  strwidth,
+void  kshakeh_update2_fix1 (
+      const float  strwidth,
 		::brook::stream invmap,
 		::brook::stream posq,
 		::brook::stream posqp,