Added velocity-center-ofmass removal

platforms/brook/src/BrookCommon.cpp

@@ -93,6 +93,11 @@ const std::string BrookCommon::ShakeInverseMapStream
 
 const std::string BrookCommon::ShuffleStream                                      = "ShuffleStream";
 
+// Random number streams
+
+const std::string BrookCommon::BrookVelocityCenterOfMassRemovalWorkStream         = "VelocityCenterOfMassRemovalWorkStream";
+const std::string BrookCommon::BrookVelocityCenterOfMassRemovalMassStream         = "VelocityCenterOfMassRemovalMassStream";
+
 /** 
  * Constructor
  * 

platforms/brook/src/BrookCommon.h

@@ -108,6 +108,11 @@ class BrookCommon {
 
       static const std::string ShuffleStream;
 
+      // BrookVelocityCenterOfMassRemoval streams
+
+      static const std::string BrookVelocityCenterOfMassRemovalWorkStream;
+      static const std::string BrookVelocityCenterOfMassRemovalMassStream;
+
       // ---------------------------------------------------------------------------------------
 
      /** 

platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp

@@ -106,6 +106,7 @@ void BrookIntegrateLangevinStepKernel::initialize( const vector<double>& masses,
 
    _brookRandomNumberGenerator   = new BrookRandomNumberGenerator( );
    _brookRandomNumberGenerator->setup( (int) masses.size(), getPlatform() );
+   _brookRandomNumberGenerator->setVerbosity( 1 );
 }
 
 /** 

platforms/brook/src/BrookRandomNumberGenerator.cpp

@@ -723,13 +723,6 @@ BrookFloatStreamInternal* BrookRandomNumberGenerator::getRandomNumberStream( int
       return NULL;
    }
 
-   if( index < 0 || index >= _numberOfRandomNumberStreams ){
-      std::stringstream message;
-      message << methodName << " index=" << index << " is out of range: [0," << _numberOfRandomNumberStreams;
-      throw OpenMMException( message.str() );
-      return NULL;
-   }
-
    if( _randomNumberGeneratorStreams == NULL ){
       std::stringstream message;
       message << methodName << " rv streams not initialized; input index=" << index;
@@ -888,6 +881,11 @@ int BrookRandomNumberGenerator::setup( int numberOfAtoms,  const Platform& platf
 
    _initializeStreamSizes( numberOfAtoms, platform );
    _initializeStreams( platform );
+   if( UseOriginalRng ){
+      _loadGVStreamsOriginal( );
+   } else {
+      _loadRandomNumberStreamsKiss( );
+   }
 
    return DefaultReturnValue;
 }
@@ -926,8 +924,8 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
    (void) LOCAL_SPRINTF( value, "%s", (UseOriginalRng ? "Original Rng" : "Kiss Rng") );
    message << _getLine( tab, "Random number generator:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamWidth() );
-   message << _getLine( tab, "RandomNumber stream width:", value ); 
+   (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamSize() );
+   message << _getLine( tab, "RandomNumber stream size:", value ); 
 
    (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamWidth() );
    message << _getLine( tab, "RandomNumber stream width:", value ); 
@@ -957,6 +955,29 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
 
    message << _getLine( tab, "Shuffle:",              (_getShuffleStream()      ? Set : NotSet) ); 
  
+   message << "\n   Statistics:\n";
+
+   double statistics[7];
+   getStatistics( statistics, -1 );
+
+   (void) LOCAL_SPRINTF( value, "%.5e", statistics[4] );
+   message << _getLine( tab, "Count:",  value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5e", statistics[0] );
+   message << _getLine( tab, "Average:",  value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5e", statistics[1] );
+   message << _getLine( tab, "StdDev:",  value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5e", statistics[3] );
+   message << _getLine( tab, "Kurtosis:",  value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5e", statistics[5] );
+   message << _getLine( tab, "Min:",  value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.6e", statistics[6] );
+   message << _getLine( tab, "Max:",  value ); 
+
    for( int ii = 0; ii < LastStreamIndex; ii++ ){
       message << std::endl;
       if( _auxiliaryStreams[ii] ){
@@ -968,3 +989,75 @@ std::string BrookRandomNumberGenerator::getContentsString( int level ) const {
 
    return message.str();
 }
+
+/* 
+ * Get statistics
+ *
+ * @param statistics   array of size 5:
+ *                       0: mean
+ *                       1: std dev
+ *                       2: 3rd moment (not normalized)
+ *                       3: kurtosis
+ *                       4: count 
+ *                       5: min
+ *                       6: max
+ *
+ * @param streamIndex  stream index to analyze
+ *
+ * @return DefaultReturnValue
+ *
+ * */
+
+int BrookRandomNumberGenerator::getStatistics( double statistics[7], int streamIndex ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookRandomNumberGenerator::";
+
+// ---------------------------------------------------------------------------------------
+
+   statistics[0] = 0.0;
+   statistics[1] = 0.0;
+   statistics[2] = 0.0;
+   statistics[3] = 0.0;
+   statistics[4] = 0.0;
+   statistics[5] =  1.0e+99;
+   statistics[6] = -1.0e+99;
+
+   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;
+         const float* dataArray = (float*) dataArrayV;
+         for( int ii = 0; ii < streamSize; ii++, index++ ){
+
+             statistics[0] += dataArray[index];
+
+             double rv2     = dataArray[index]*dataArray[index];
+             statistics[1] += rv2;
+             statistics[2] += rv2*dataArray[index];
+             statistics[3] += rv2*rv2;
+
+             if( statistics[5] > dataArray[index] ){
+                statistics[5] = dataArray[index];
+             }
+             if( statistics[6] < dataArray[index] ){
+                statistics[6] = dataArray[index];
+             }
+         } 
+         statistics[4] += (double) index;
+      }   
+   }   
+
+   if( statistics[4] > 0.0 ){ 
+      statistics[0] /= statistics[4];
+      statistics[1]  = statistics[1] - statistics[4]*statistics[0]*statistics[0];
+      if( statistics[4] > 1.0 ){
+         statistics[1] = sqrt( statistics[1]/( statistics[4] - 1.0 ) );
+      }
+      statistics[3]  = (statistics[3]/(statistics[4]*statistics[1]*statistics[1]) ) - 3.0;
+   }
+
+   return DefaultReturnValue;
+}

platforms/brook/src/BrookRandomNumberGenerator.h

@@ -48,7 +48,7 @@ class BrookRandomNumberGenerator : public BrookCommon {
   
       // toggle between original rng & Kiss (Nvidia) code
 
-      static const int UseOriginalRng = 1;
+      static const int UseOriginalRng = 0;
 
       /** 
        * Constructor
@@ -205,6 +205,28 @@ class BrookRandomNumberGenerator : public BrookCommon {
       
       int getRvStreamOffset( void ) const;
       
+      /* 
+       * Get statistics
+       *
+       * @param statistics   array of size 7:
+       *                       0: mean
+       *                       1: std dev
+       *                       2: 3rd moment (not normalized)
+       *                       3: kurtosis
+       *                       4: count 
+       *                       5: min
+       *                       6: max
+       *
+       * @param streamIndex  stream index to analyze
+       *
+       * @return DefaultReturnValue
+       *
+       * */
+      
+      int getStatistics( double statistics[7], int streamIndex ) const;
+      
+      // ---------------------------------------------------------------------------------------
+      
    private:
    
       // streams indices

platforms/brook/src/BrookStochasticDynamics.cpp

@@ -91,6 +91,8 @@ BrookStochasticDynamics::BrookStochasticDynamics( ){
 
    _randomNumberSeed  = 1393;
 
+   _brookVelocityCenterOfMassRemoval = NULL;
+
    //_randomNumberSeed = randomNumberSeed ? randomNumberSeed : 1393;
    //SimTKOpenMMUtilities::setRandomNumberSeed( randomNumberSeed );
 }   
@@ -114,6 +116,8 @@ BrookStochasticDynamics::~BrookStochasticDynamics( ){
 
    delete[] _inverseSqrtMasses;
 
+   delete _brookVelocityCenterOfMassRemoval;
+
 }
 
 /** 
@@ -375,7 +379,7 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
 
    static const char* methodName = "\nBrookStochasticDynamics::update";
 
-   static const int PrintOn      = 1;
+   static const int PrintOn      = 0;
 
 // ---------------------------------------------------------------------------------------
 
@@ -386,10 +390,26 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
    const BrookStreamImpl& forceStreamC                 = dynamic_cast<const BrookStreamImpl&> (forces.getImpl());
    BrookStreamImpl& forceStream                        = const_cast<BrookStreamImpl&> (forceStreamC);
 
-   if( PrintOn && getLog() ){
-      (void) fprintf( getLog(), "%s shake=%d\n", methodName, brookShakeAlgorithm.getNumberOfConstraints() );
-      (void) fflush( getLog() );
+   if( (1 || PrintOn) && getLog() ){
+
+      static int showAux = 1;
+
+      if( PrintOn ){
+         (void) fprintf( getLog(), "%s shake=%d\n", methodName, brookShakeAlgorithm.getNumberOfConstraints() );
+         (void) fflush( getLog() );
+      }
 
+      // show update
+   
+      if( showAux ){
+         showAux = 0;
+         std::string contents = _brookVelocityCenterOfMassRemoval->getContentsString( );
+         (void) fprintf( getLog(), "%s VelocityCenterOfMassRemoval contents\n%s", methodName, contents.c_str() );
+         contents             = brookRandomNumberGenerator.getContentsString( );
+         (void) fprintf( getLog(), "%s RNG contents\n%s", methodName, contents.c_str() );
+         brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); 
+         (void) fflush( getLog() );
+      }
    }
 
    // first integration step
@@ -417,7 +437,7 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
    // diagnostics
 
    if( PrintOn ){
-      (void) fprintf( getLog(), "\nPost kupdate_sd1_fix1: atomStrW=%3d rngStrW=%3d rngOff=%3d"
+      (void) fprintf( getLog(), "\nPost kupdate_sd1_fix1: atomStrW=%3d rngStrW=%3d rngOff=%5d "
                                 "EM=%12.5e Sd1pc[]=[%12.5e %12.5e %12.5e]",
                                 getStochasticDynamicsAtomStreamWidth(),
                                 brookRandomNumberGenerator.getRandomNumberStreamWidth(),
@@ -455,6 +475,9 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
 
       (void) fprintf( getLog(), "\nXPrimeStream\n" );
       getXPrimeStream()->printToFile( getLog() ); 
+
+      (void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() );
+      // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); 
    }   
 
    // advance random number cursor
@@ -513,6 +536,43 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
          getXPrimeStream()->getBrookStream() 
          );
 
+   // diagnostics
+
+   if( PrintOn ){
+      (void) fprintf( getLog(), "\nPost kupdate_sd2_fix1: atomStrW=%3d rngStrW=%3d rngOff=%5d "
+                                "Sd2pc[]=[%12.5e %12.5e]",
+                                getStochasticDynamicsAtomStreamWidth(),
+                                brookRandomNumberGenerator.getRandomNumberStreamWidth(),
+                                brookRandomNumberGenerator.getRvStreamOffset(),
+                                derivedParameters[Sd2pc1], derivedParameters[Sd2pc2] );
+
+      (void) fprintf( getLog(), "\nSDPC2Stream\n" );
+      getSDPC2Stream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nSD2XStream\n" );
+      getSD1VStream()->printToFile( getLog() );
+
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+      (void) fprintf( getLog(), "\nPositionStream\n" );
+      brookStreamInternalPos->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nVPrimeStream\n" );
+      getVPrimeStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nXPrimeStream\n" );
+      getXPrimeStream()->printToFile( getLog() ); 
+
+      (void) fprintf( getLog(), "\ngetSD2XStream\n" );
+      getSD2XStream()->printToFile( getLog() );
+
+      BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamImpl();
+      (void) fprintf( getLog(), "\nVelocityStream\n" );
+      brookStreamInternalVel->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() );
+//      brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); 
+   }   
+
    // advance random number cursor
 
    brookRandomNumberGenerator.advanceGVCursor( 2*getNumberOfAtoms() );
@@ -552,6 +612,8 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
       ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
    }
 
+   _brookVelocityCenterOfMassRemoval->removeVelocityCenterOfMass( velocities );
+
    return DefaultReturnValue;
 
 };
@@ -997,6 +1059,9 @@ int BrookStochasticDynamics::setup( const std::vector<double>& masses, const Pla
 
    _setInverseSqrtMasses( masses );
 
+   _brookVelocityCenterOfMassRemoval = new BrookVelocityCenterOfMassRemoval();
+   _brookVelocityCenterOfMassRemoval->setup( masses, platform );
+
    return DefaultReturnValue;
 }
 

platforms/brook/src/BrookStochasticDynamics.h

@@ -38,6 +38,7 @@
 #include "BrookFloatStreamInternal.h"
 #include "BrookShakeAlgorithm.h"
 #include "BrookRandomNumberGenerator.h"
+#include "BrookVelocityCenterOfMassRemoval.h"
 #include "BrookPlatform.h"
 #include "BrookCommon.h"
 
@@ -323,6 +324,10 @@ class BrookStochasticDynamics : public BrookCommon {
 
       BrookOpenMMFloat* _inverseSqrtMasses;
 
+      // remove com
+
+      BrookVelocityCenterOfMassRemoval* _brookVelocityCenterOfMassRemoval;
+
       // internal streams
 
       BrookFloatStreamInternal* _sdStreams[LastStreamIndex];

platforms/brook/src/BrookVelocityCenterOfMassRemoval.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * 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        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Mark Friedrichs                                                   *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include <sstream>
+#include "BrookVelocityCenterOfMassRemoval.h"
+#include "BrookPlatform.h"
+#include "OpenMMException.h"
+#include "BrookStreamImpl.h"
+#include "gpu/kcom.h"
+
+using namespace OpenMM;
+using namespace std;
+
+/** 
+ *
+ * Constructor
+ * 
+ */
+
+BrookVelocityCenterOfMassRemoval::BrookVelocityCenterOfMassRemoval( ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVelocityCenterOfMassRemoval::BrookVelocityCenterOfMassRemoval";
+
+   BrookOpenMMFloat zero                    = (BrookOpenMMFloat)  0.0;
+
+// ---------------------------------------------------------------------------------------
+
+   _numberOfAtoms             = -1;
+
+   // mark stream dimension variables as unset
+
+   _atomStreamWidth         = -1;
+   _atomStreamHeight        = -1;
+   _atomStreamSize          = -1;
+
+   _totalInverseMass        = zero;
+
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      _streams[ii]   = NULL;
+   }
+
+}   
+ 
+/** 
+ * Destructor
+ * 
+ */
+
+BrookVelocityCenterOfMassRemoval::~BrookVelocityCenterOfMassRemoval( ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVelocityCenterOfMassRemoval::~BrookVelocityCenterOfMassRemoval";
+
+// ---------------------------------------------------------------------------------------
+
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      delete _streams[ii];
+   }
+
+}
+
+/** 
+ * Get inverse of total mass
+ * 
+ * @return  inverse of total mass
+ *
+ */
+
+BrookOpenMMFloat BrookVelocityCenterOfMassRemoval::getTotalInverseMass( void ) const {
+   return _totalInverseMass;
+}
+
+/** 
+ * Get inverse mass stream 
+ *
+ * @return  inverse mass stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookVelocityCenterOfMassRemoval::getMassStream( void ) const {
+   return _streams[MassStream];
+}
+
+/** 
+ * Get work stream 
+ *
+ * @return  work stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookVelocityCenterOfMassRemoval::getWorkStream( void ) const {
+   return _streams[WorkStream];
+}
+
+/** 
+ * Get linear momentum stream
+ *
+ * @return   linear momentum stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookVelocityCenterOfMassRemoval::getLinearMomentumStream( void ) const {
+   return _streams[LinearMomentumStream];
+}
+
+/**
+ * Remove velocity-COM
+ *
+ * @param velocities velocities
+ *
+ * @return DefaultReturnValue
+ *
+ */
+   
+int BrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass( Stream& velocities ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nBrookVelocityCenterOfMassRemoval::removeVelocityCenterOfMass";
+
+   static int debug         = 0;
+
+   // static char* testName[2] = { "kCalculateLinearMomentum", "kRemoveLinearMomentum" };
+   // char fileName[128];
+
+   float zero               = 0.0f;
+   float one                = 1.0f;
+
+   // ---------------------------------------------------------------------------------------
+
+   // calculate linear momentum via reduction
+   // subtract it (/totalMass) from velocities
+
+   BrookStreamImpl& velocityStream  = dynamic_cast<BrookStreamImpl&> (velocities.getImpl());
+
+   kCalculateLinearMomentum( getMassStream()->getBrookStream(), velocityStream.getBrookStream(), getWorkStream()->getBrookStream() );
+   kScale( zero, getLinearMomentumStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
+   kSumLinearMomentum( (float) getComAtomStreamWidth(), (float) getNumberOfAtoms(), getWorkStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
+   kScale( (float) getTotalInverseMass(), getLinearMomentumStream()->getBrookStream(), getLinearMomentumStream()->getBrookStream() );
+   kRemoveLinearMomentum( getLinearMomentumStream()->getBrookStream(), velocityStream.getBrookStream(), velocityStream.getBrookStream() );
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Get Atom stream size
+ *
+ * @return  Atom stream size
+ *
+ */
+
+int BrookVelocityCenterOfMassRemoval::getComAtomStreamSize( void ) const {
+   return _atomStreamSize;
+}
+
+/** 
+ * Get atom stream width
+ *
+ * @return  atom stream width
+ *
+ */
+
+int BrookVelocityCenterOfMassRemoval::getComAtomStreamWidth( void ) const {
+   return _atomStreamWidth;
+}
+
+/** 
+ * Get atom stream height
+ *
+ * @return atom stream height
+ */
+
+int BrookVelocityCenterOfMassRemoval::getComAtomStreamHeight( void ) const {
+   return _atomStreamHeight;
+}
+
+/** 
+ * Initialize stream dimensions
+ * 
+ * @param numberOfAtoms             number of atoms
+ * @param platform                  platform
+ *      
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookVelocityCenterOfMassRemoval::_initializeStreamSizes( int numberOfAtoms, const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVelocityCenterOfMassRemoval::_initializeStreamSizes";
+
+// ---------------------------------------------------------------------------------------
+
+   _atomStreamSize     = getAtomStreamSize( platform );
+   _atomStreamWidth    = getAtomStreamWidth( platform );
+   _atomStreamHeight   = getAtomStreamHeight( platform );
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Initialize streams
+ * 
+ * @param platform                  platform
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookVelocityCenterOfMassRemoval::_initializeStreams( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVelocityCenterOfMassRemoval::_initializeStreams";
+
+   BrookOpenMMFloat dangleValue            = (BrookOpenMMFloat) 0.0;
+
+// ---------------------------------------------------------------------------------------
+
+   int atomStreamSize               = getComAtomStreamSize();
+   int atomStreamWidth              = getComAtomStreamWidth();
+
+    _streams[WorkStream]            = new BrookFloatStreamInternal( BrookCommon::BrookVelocityCenterOfMassRemovalWorkStream,
+                                                                    atomStreamSize, atomStreamWidth,
+                                                                    BrookStreamInternal::Float, dangleValue );
+
+
+    _streams[MassStream]            = new BrookFloatStreamInternal( BrookCommon::BrookVelocityCenterOfMassRemovalMassStream,
+                                                                    atomStreamSize, atomStreamWidth,
+                                                                    BrookStreamInternal::Float, dangleValue );
+
+
+    _streams[LinearMomentumStream]  = new BrookFloatStreamInternal( "LinearMomentumStream",
+                                                                    1, 3, BrookStreamInternal::Float, dangleValue );
+
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Set inverse masses 
+ * 
+ * @param masses             atomic masses
+ *
+ */
+
+int BrookVelocityCenterOfMassRemoval::_setMasses( const std::vector<double>& masses ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookVelocityCenterOfMassRemoval::_setMasses";
+
+   BrookOpenMMFloat zero                    = (BrookOpenMMFloat)  0.0;
+   BrookOpenMMFloat one                     = (BrookOpenMMFloat)  1.0;
+
+// ---------------------------------------------------------------------------------------
+
+   if( (int) masses.size() > getComAtomStreamSize() ){
+      std::stringstream message;
+      message << methodName << " mass array size=" << masses.size() << " is larger than stream size=" << getComAtomStreamSize();
+      throw OpenMMException( message.str() );
+   }
+
+
+   // setup inverse sqrt masses
+
+   BrookOpenMMFloat* inverseMasses = new BrookOpenMMFloat[getComAtomStreamSize()];
+   memset( inverseMasses, 0, sizeof( BrookOpenMMFloat )*getComAtomStreamSize() );
+
+   int index          = 0;
+   _totalInverseMass  = (BrookOpenMMFloat) 0.0;
+   for( std::vector<double>::const_iterator ii = masses.begin(); ii != masses.end(); ii++, index++ ){
+      if( *ii != 0.0 ){
+         BrookOpenMMFloat value    = static_cast<BrookOpenMMFloat>(*ii);
+         inverseMasses[index]      = value;
+         _totalInverseMass        += value;
+      } else {
+         inverseMasses[index] = zero;
+      }
+   }
+   if( _totalInverseMass > zero ){
+      _totalInverseMass = one/_totalInverseMass;
+   }
+
+   _streams[MassStream]->loadFromArray( inverseMasses );
+
+   delete[] inverseMasses;
+
+   return DefaultReturnValue;
+}   
+ 
+/*  
+ * Setup of StochasticDynamics parameters
+ *
+ * @param masses                masses
+ * @param platform              Brook platform
+ *
+ * @return nonzero value if error
+ *
+ * */
+    
+int BrookVelocityCenterOfMassRemoval::setup( const std::vector<double>& masses, const Platform& platform ){
+    
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookVelocityCenterOfMassRemoval::setup";
+
+// ---------------------------------------------------------------------------------------
+
+   const BrookPlatform brookPlatform        = dynamic_cast<const BrookPlatform&> (platform);
+   setLog( brookPlatform.getLog() );
+
+   int numberOfAtoms  = (int) masses.size();
+   setNumberOfAtoms( numberOfAtoms );
+
+   // set stream sizes and then create streams
+
+   _initializeStreamSizes( numberOfAtoms, platform );
+   _initializeStreams( platform );
+
+   _setMasses( masses );
+
+   return DefaultReturnValue;
+}
+
+/* 
+ * Get contents of object
+ *
+ * @param level   level of dump
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string BrookVelocityCenterOfMassRemoval::getContentsString( int level ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookVelocityCenterOfMassRemoval::getContentsString";
+
+   static const unsigned int MAX_LINE_CHARS = 256;
+   char value[MAX_LINE_CHARS];
+   static const char* Set                   = "Set";
+   static const char* NotSet                = "Not set";
+
+// ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   std::string tab   = "   ";
+
+#ifdef WIN32
+#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
+#else
+#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
+#endif
+
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfAtoms() );
+   message << _getLine( tab, "Number of atoms:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getComAtomStreamWidth() );
+   message << _getLine( tab, "Atom stream width:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getComAtomStreamHeight() );
+   message << _getLine( tab, "Atom stream height:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getComAtomStreamSize() );
+   message << _getLine( tab, "Atom stream size:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5f", getTotalInverseMass() );
+   message << _getLine( tab, "TotalInverseMass:", value ); 
+
+   message << _getLine( tab, "Log:",                  (getLog()                  ? Set : NotSet) ); 
+
+   message << _getLine( tab, "Mass:",                 (getMassStream()           ? Set : NotSet) ); 
+   message << _getLine( tab, "Work:",                 (getWorkStream()           ? Set : NotSet) ); 
+   message << _getLine( tab, "LinearMomentum:",       (getLinearMomentumStream() ? Set : NotSet) ); 
+ 
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      message << std::endl;
+      if( _streams[ii] ){
+         message << _streams[ii]->getContentsString( );
+      }
+   }
+
+#undef LOCAL_SPRINTF
+
+   return message.str();
+}

platforms/brook/src/BrookVelocityCenterOfMassRemoval.h

+#ifndef OPENMM_BROOK_VELOCITY_CENTER_OF_MASS_REMOVAL_H_
+#define OPENMM_BROOK_VELOCITY_CENTER_OF_MASS_REMOVAL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * 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        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include <vector>
+#include <set>
+
+#include "BrookFloatStreamInternal.h"
+#include "BrookPlatform.h"
+#include "BrookCommon.h"
+
+namespace OpenMM {
+
+/**
+ *
+ * Encapsulates removal of center of mass
+ *
+ */
+
+class BrookVelocityCenterOfMassRemoval : public BrookCommon {
+
+   public:
+  
+      /** 
+       * Constructor
+       * 
+       */
+      
+      BrookVelocityCenterOfMassRemoval( );
+  
+      /** 
+       * Destructor
+       * 
+       */
+      
+      ~BrookVelocityCenterOfMassRemoval( );
+  
+      /**
+       * Get  atom stream width
+       *
+       * @return atom stream width
+       */
+
+      int getComAtomStreamWidth( void ) const; 
+
+      /**
+       * Get  atom stream height
+       *
+       * @return atom stream height
+       */
+
+      int getComAtomStreamHeight( void ) const;
+
+      /**
+       * Get  atom stream size
+       * 
+       * @return atom stream size
+       */
+
+      int getComAtomStreamSize( void ) const; 
+
+      /** 
+       * Update
+       * 
+       * @param  positions                   atom positions
+       * @param  velocities                  atom velocities
+       * @param  forces                      atom forces
+       * @param  brookShakeAlgorithm         BrookShakeAlgorithm reference
+       * @param  brookRandomNumberGenerator  BrookRandomNumberGenerator reference
+       *
+       * @return  DefaultReturnValue
+       *
+       */
+      
+      int removeVelocityCenterOfMass( Stream& velocities );
+
+      /* 
+       * Setup of  parameters
+       *
+       * @param masses                atom masses
+       * @param platform              Brook platform
+       *
+       * @return ErrorReturnValue value if error, else DefaultReturnValue
+       *
+       * */
+      
+      int setup( const std::vector<double>& masses, const Platform& platform  );
+      
+      /* 
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string getContentsString( int level = 0 ) const;
+
+      /**
+       * Get inverse of total mass
+       *
+       * @return inverse of total mass
+       *
+       */
+         
+      BrookOpenMMFloat getTotalInverseMass( void ) const;
+      
+      /** 
+       * Get inverse mass stream 
+       *
+       * @return  inverse mass stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getMassStream( void ) const;
+      
+      /** 
+       * Get work stream 
+       *
+       * @return  work stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getWorkStream( void ) const;
+      
+      /** 
+       * Get linear momentum stream 
+       *
+       * @return  linear momentum stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getLinearMomentumStream( void ) const;
+      
+   private:
+   
+      // streams indices
+
+      enum BrookVelocityCenterOfMassRemovalStreams { 
+              MassStream,
+              WorkStream,
+              LinearMomentumStream,
+              LastStreamIndex
+           };
+
+      BrookOpenMMFloat _totalInverseMass;
+
+      // Atom stream dimensions
+
+      int _atomStreamWidth;
+      int _atomStreamHeight;
+      int _atomStreamSize;
+
+      // internal streams
+
+      BrookFloatStreamInternal* _streams[LastStreamIndex];
+
+      /* 
+       * Setup of stream dimensions
+       *
+       * @param atomStreamSize        atom stream size
+       * @param atomStreamWidth       atom stream width
+       *
+       * @return ErrorReturnValue if error, else DefaultReturnValue
+       *
+       * */
+      
+      int _initializeStreamSizes( int atomStreamSize, int atomStreamWidth );
+
+      /** 
+       * Initialize stream dimensions
+       * 
+       * @param numberOfAtoms             number of atoms
+       * @param platform                  platform
+       *
+       * @return ErrorReturnValue if error, else DefaultReturnValue
+       *
+       */
+      
+      int _initializeStreamSizes( int numberOfAtoms, const Platform& platform );
+      
+      /** 
+       * Initialize stream dimensions and streams
+       * 
+       * @param platform                  platform
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int _initializeStreams( const Platform& platform );
+
+      /** 
+       * Set inverse masses 
+       * 
+       * @param masses             atomic masses
+       *
+       */
+      
+      int _setMasses( const std::vector<double>& masses );
+      
+      
+};
+
+} // namespace OpenMM
+
+#endif /* OPENMM_BROOK_VELOCITY_CENTER_OF_MASS_REMOVAL_H_ */

platforms/brook/src/gpu/kcom.h

+
+void  kCalculateLinearMomentum( ::brook::stream mass, ::brook::stream velocities, ::brook::stream linearMomentum );
+void  kReduceLinearMomentum( ::brook::stream momentum, ::brook::stream linearMomentum );
+//void  kSumLinearMomentum( ::brook::stream momentum, float3* linearMomentum );
+void  kScale( float scale, ::brook::stream linearMomentumIn, ::brook::stream linearMomentumOut );
+void  kRemoveLinearMomentum( ::brook::stream linearMomentum, ::brook::stream velocitiesIn, ::brook::stream velocitiesOut );
+void kSum( ::brook::stream array, ::brook::stream sum );
+
+/**---------------------------------------------------------------------------------------
+
+   This kernel calculates the total linear momentum via a reduction
+
+   @param atomStrWidth   atom stream width
+   @param numberOfAtoms  number of atoms
+   @param momentum       momentum
+   @param linearMomentum total momentum
+
+   --------------------------------------------------------------------------------------- */
+
+void kSumLinearMomentum( float atomStrWidth, float numberOfAtoms, ::brook::stream momentum, ::brook::stream linearMomentum );
+
+