Daily commit

platforms/brook/src/BrookCommon.cpp

@@ -30,7 +30,10 @@
  * -------------------------------------------------------------------------- */
 
 #include <math.h>
+#include <stdlib.h>
+
 #include <sstream>
+
 #include "BrookCommon.h"
 #include "BrookPlatform.h"
 #include "BrookStreamFactory.h"
@@ -73,6 +76,16 @@ const std::string BrookCommon::SDPC2Stream
 const std::string BrookCommon::SD2XStream                                         = "SD2XStream";
 const std::string BrookCommon::SD1VStream                                         = "SD1VStream";
 
+// Shake streams
+
+const std::string BrookCommon::ShakeAtomIndicesStream                             = "ShakeAtomIndicesStream";
+const std::string BrookCommon::ShakeAtomParameterStream                           = "ShakeAtomParameterStream";
+const std::string BrookCommon::ShakeXCons0Stream                                  = "ShakeXCons0Stream";
+const std::string BrookCommon::ShakeXCons1Stream                                  = "ShakeXCons1Stream";
+const std::string BrookCommon::ShakeXCons2Stream                                  = "ShakeXCons2Stream";
+const std::string BrookCommon::ShakeXCons3Stream                                  = "ShakeXCons3Stream";
+const std::string BrookCommon::ShakeInverseMapStream                              = "ShakeInverseMapStream";
+
 /** 
  * Constructor
  * 
@@ -354,3 +367,90 @@ std::string BrookCommon::_getLine( const std::string& tab, const std::string& de
    return message.str();
 
 }
+
+/* 
+ * Given number of stream elements and width, returns the appropriate
+ * height of the stream
+ *
+ * @param streamSize   stream size 
+ * @param width        stream width
+ *
+ * @return stream height
+ *
+ */
+
+int BrookCommon::getStreamHeight( int streamSize, int streamWidth ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookCommon::getStreamHeight";
+
+// ---------------------------------------------------------------------------------------
+
+	int streamHeight = streamSize/streamWidth;
+
+	if( streamSize % streamWidth ){
+      streamHeight++;
+   }
+
+	return streamHeight;
+}
+
+/* 
+ * Given number of stream elements, get stream width & height
+ *
+ * @param streamSize    stream size 
+ * @param streamWidth   output stream width
+ * @param streamHeight  output stream height
+ *
+ * @return stream height
+ *
+ */
+
+void BrookCommon::getStreamDimensions( int streamSize, int *streamWidth, int *streamHeight ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookCommon::getStreamDimensions";
+
+// ---------------------------------------------------------------------------------------
+
+	// There are two conditions - stream should be as square 
+	// as possible, but should also be multiple of 16 along
+	// one dimension
+	
+	float s = sqrtf( (float) streamSize );
+
+	// find nearest multiple of 16 to the perfect square size
+
+	int low  = ( (int) floor( s/16.0f ) ) * 16;
+
+	if ( !low ) {
+		*streamWidth  = 16;
+		*streamHeight = getStreamHeight( streamSize, *streamWidth );
+	} else { 
+
+	   int high = low + 1;
+
+		// I'm not sure 48 is such a good stream width. Things seem
+		// to be faster with 32 or 64. Can make this a special
+		// case later.
+		
+		// Choose low or high depending on which one is 
+		// more square
+
+		int htlow  = getStreamHeight( streamSize, low );
+		int hthigh = getStreamHeight( streamSize, high );
+		
+		if ( abs( htlow - low ) < abs( hthigh - high ) ) {
+			*streamWidth   = low;
+			*streamHeight  = htlow;
+		} else {
+			*streamWidth   = high;
+			*streamHeight  = hthigh;
+		}
+	}
+	return;
+}
+
+

platforms/brook/src/BrookCommon.h

@@ -84,6 +84,23 @@ class BrookCommon {
       static const std::string ObcIntermediateForceStream;
       static const std::string ObcChainStream;
 
+      // Stochastic Dynamics streams
+
+      static const std::string SDPC1Stream;
+      static const std::string SDPC2Stream;
+      static const std::string SD2XStream;
+      static const std::string SD1VStream;
+
+      // Shake streams
+
+      static const std::string ShakeAtomIndicesStream;
+      static const std::string ShakeAtomParameterStream;
+      static const std::string ShakeXCons0Stream;
+      static const std::string ShakeXCons1Stream;
+      static const std::string ShakeXCons2Stream;
+      static const std::string ShakeXCons3Stream;
+      static const std::string ShakeInverseMapStream;
+
       // ---------------------------------------------------------------------------------------
 
      /** 
@@ -203,6 +220,32 @@ class BrookCommon {
       
       FILE* getLog( void ) const;
 
+      /* 
+       * Given number of stream elements and width, returns the appropriate
+       * height of the stream
+       *
+       * @param streamSize   stream size 
+       * @param width        stream width
+       *
+       * @return stream height
+       *
+       */
+      
+      static int getStreamHeight( int streamSize, int streamWidth );
+      
+      /* 
+       * Given number of stream elements, get stream width & height
+       *
+       * @param streamSize    stream size 
+       * @param streamWidth   output stream width
+       * @param streamHeight  output stream height
+       *
+       * @return stream height
+       *
+       */
+      
+      static void getStreamDimensions( int streamSize, int *streamWidth, int *streamHeight );
+      
    protected:
    
       // number of atoms

platforms/brook/src/BrookGbsa.cpp

@@ -65,6 +65,7 @@ BrookGbsa::BrookGbsa(  ){
 
    _duplicationFactor         =  4;
 
+   _includeAce                = 1;
    _solventDielectric         = 78.3;
    _soluteDielectric          = 1.0;
    _dielectricOffset          = 0.09;
@@ -117,6 +118,18 @@ int BrookGbsa::getNumberOfForceStreams( void ) const {
    return NumberOfForceStreams;
 }
 
+
+/** 
+ * Include ACE approximation in calculation of force
+ * 
+ * @return true if ACE approximation is to be included in calculation of force
+ *
+ */
+
+int BrookGbsa::includeAce( void ) const {
+   return _includeAce;
+}
+
 /** 
  * Get inner loop unroll
  * 

platforms/brook/src/BrookGbsa.h

@@ -369,6 +369,10 @@ class BrookGbsa : public BrookCommon {
 
       int _duplicationFactor;
 
+      // include ACE approximation
+
+      int _includeAce;
+
       // force stream width
 
       int _partialForceStreamWidth;

platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp

@@ -29,7 +29,7 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "BrookIntegrateKernals.h"
+#include "BrookIntegrateLangevinStepKernel.h"
 #include "BrookStreamInternal.h"
 
 using namespace OpenMM;
@@ -44,14 +44,8 @@ BrookIntegrateLangevinStepKernel::BrookIntegrateLangevinStepKernel( std::string
 
 // ---------------------------------------------------------------------------------------
 
-   _numberOfConstraints            = -1;
-
    _brookStochasticDynamics        = NULL;
    _brookShakeAlgorithm            = NULL;
-   _atomMasses                     = NULL;
-   _shakeParameters                = NULL;
-   _constraintIndices              = NULL;
-   _shakeParameters                = NULL;
    
 }
 
@@ -63,12 +57,8 @@ BrookIntegrateLangevinStepKernel::~BrookIntegrateLangevinStepKernel( ){
 
 // ---------------------------------------------------------------------------------------
    
-   delete _brookkStochasticDynamics;
+   delete _brookStochasticDynamics;
    delete _brookShakeAlgorithm;
-   delete _atomMasses;
-   delete _shakeParameters;
-   delete _constraintIndices;
-   delete _shakeParameters;
 
 }
 
@@ -82,25 +72,11 @@ void BrookIntegrateLangevinStepKernel::initialize( const vector<double>& masses,
 
 // ---------------------------------------------------------------------------------------
    
-   _brookkStochasticDynamics = new BrookStochasticDynamics( masses );
-   _brookkStochasticDynamics->setup( masses, getPlatform() );
-
-   _brookShakeAlgorithm      = new BrookShakeAlgorithm( masses, constraintIndices, constraintLengths );
-/*
-    this->masses = new RealOpenMM[masses.size()];
-    for (size_t i = 0; i < masses.size(); ++i)
-        this->masses[i] = static_cast<RealOpenMM>( masses[i] );
-    numConstraints = constraintIndices.size();
-    this->constraintIndices = allocateIntArray(numConstraints, 2); 
-    for (int i = 0; i < numConstraints; ++i) {
-        this->constraintIndices[i][0] = constraintIndices[i][0];
-        this->constraintIndices[i][1] = constraintIndices[i][1];
-    }   
-    shakeParameters = allocateRealArray(constraintLengths.size(), 1); 
-    for (size_t i = 0; i < constraintLengths.size(); ++i)
-        shakeParameters[i][0] = static_cast<RealOpenMM>( constraintLengths[i] );
+   _brookStochasticDynamics = new BrookStochasticDynamics( );
+   _brookStochasticDynamics->setup( masses, getPlatform() );
 
-*/
+   _brookShakeAlgorithm      = new BrookShakeAlgorithm( );
+   _brookShakeAlgorithm->setup( masses, constraintIndices, constraintLengths, getPlatform() );
 }
 
 /** 
@@ -121,6 +97,7 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc
 
 // ---------------------------------------------------------------------------------------
 
+   double epsilon = 1.0e-04;
    //static const std::string methodName      = "BrookIntegrateLangevinStepKernel::execute";
 
 // ---------------------------------------------------------------------------------------
@@ -132,12 +109,13 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc
 
    // take step
 
-   if( _brookStochasticDynamics == NULL ){
-      _brookStochasticDynamics = new BrookStochasticDynamics( getNumberOfAtoms(), static_cast<RealOpenMM>(stepSize), 
-                                                              static_cast<RealOpenMM>(tau), static_cast<RealOpenMM>(temperature) );
-   } else if( temperature != _brookStochasticDynamics->getTemperatur() || friction != _brookStochasticDynamics->getFriction() ){
-      _brookStochasticDynamics->updateParameters( temperature, friction );
+   double differences[3];
+   differences[0] = temperature - (double) _brookStochasticDynamics->getTemperature();
+   differences[1] = friction    - (double) _brookStochasticDynamics->getFriction();
+   differences[2] = stepSize    - (double) _brookStochasticDynamics->getStepSize();
+   if( fabs( differences[0] ) < epsilon || fabs( differences[1] ) < epsilon || fabs( differences[2] ) < epsilon ){
+      _brookStochasticDynamics->updateParameters( temperature, friction, stepSize );
    } 
-   _brookStochasticDynamics->update( positions, velocities, forces );
+   _brookStochasticDynamics->update( positions, velocities, forces, *_brookShakeAlgorithm );
 
 }

platforms/brook/src/BrookIntegrateLangevinStepKernel.h

@@ -85,18 +85,12 @@ class BrookIntegrateLangevinStepKernel : public IntegrateLangevinStepKernel {
        *
        */
 
-      void execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize);
+      void execute( Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize);
 
    protected:
 
       BrookStochasticDynamics* _brookStochasticDynamics;
       BrookShakeAlgorithm*     _brookShakeAlgorithm;
-      RealOpenMM*  _atomMasses;
-      RealOpenMM** _shakeParameters;
-      int** _constraintIndices;
-      int _numberOfConstraints;
-      double prevTemp, prevFriction, prevStepSize;
-        
 
 };
 

platforms/brook/src/BrookIntegrateVerletStepKernel.cpp

@@ -29,7 +29,7 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "BrookIntegrateKernals.h"
+#include "BrookIntegrateVerletStepKernel.h"
 #include "BrookStreamInternal.h"
 
 using namespace OpenMM;

platforms/brook/src/BrookKernelFactory.cpp

@@ -31,8 +31,9 @@
 
 #include "BrookKernelFactory.h"
 #include "BrookCalcStandardMMForceFieldKernel.h"
-#include "BrookIntegrateKernals.h"
+#include "BrookIntegrateLangevinStepKernel.h"
 #include "BrookCalcKineticEnergyKernel.h"
+#include "BrookCalcGBSAOBCForceFieldKernel.h"
 
 using namespace OpenMM;
 
@@ -44,27 +45,41 @@ KernelImpl* BrookKernelFactory::createKernelImpl( std::string name, const Platfo
 
 // ---------------------------------------------------------------------------------------
 
-
    // StandardMM
 
 	if( name == CalcStandardMMForceFieldKernel::Name() ){
+
       return new BrookCalcStandardMMForceFieldKernel( name, platform );
 
    // GBSA OBC
 
 	} else if( name == CalcGBSAOBCForceFieldKernel::Name() ){
 
-      (void) fprintf( stderr, "CalcGBSAOBCForceFieldKernel not set BrookKernelFactory::createKernelImpl\n" );
-      (void) fflush( stderr );
-      //return new BrookCalcGBSAOBCForceFieldKernel(name, platform);
+      return new BrookCalcGBSAOBCForceFieldKernel(name, platform);
 
    // Verlet integrator
 
 	} else if( name == IntegrateVerletStepKernel::Name() ){
 
-      (void) fprintf( stderr, "IntegrateVerletStepKernel created BrookKernelFactory::createKernelImpl\n" );
-      (void) fflush( stderr );
-      return new BrookIntegrateVerletStepKernel( name, platform );
+      // return new BrookIntegrateVerletStepKernel( name, platform );
+
+   // Brownian integrator
+
+	} else if( name == IntegrateBrownianStepKernel::Name() ){
+
+      // return new BrookIntegrateBrownianStepKernel( name, platform );
+
+   //  Andersen integrator
+
+	} else if( name == ApplyAndersenThermostatKernel::Name() ){
+
+      // return new BrookIntegrateAndersenThermostatKernel( name, platform );
+
+   // Langevin integrator
+
+	} else if( name == IntegrateLangevinStepKernel::Name() ){
+
+      return new BrookIntegrateLangevinStepKernel( name, platform );
 
    // KE calculator
 
@@ -75,13 +90,5 @@ KernelImpl* BrookKernelFactory::createKernelImpl( std::string name, const Platfo
    (void) fprintf( stderr, "createKernelImpl: name=<%s> not found.", methodName.c_str(), name.c_str() );
    (void) fflush( stderr );
 
-/*
-    if (name == IntegrateLangevinStepKernel::Name())
-        //return new BrookIntegrateLangevinStepKernel(name, platform);
-    //if (name == IntegrateBrownianStepKernel::Name())
-        //return new BrookIntegrateBrownianStepKernel(name, platform);
-    if (name == ApplyAndersenThermostatKernel::Name())
-        //return new BrookApplyAndersenThermostatKernel(name, platform);
-*/
 	return NULL;
 }

platforms/brook/src/BrookShakeAlgorithm.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 "BrookShakeAlgorithm.h"
+#include "BrookPlatform.h"
+#include "OpenMMException.h"
+#include "BrookStreamImpl.h"
+
+// use random number generator
+
+#include "SimTKOpenMMUtilities.h"
+
+using namespace OpenMM;
+using namespace std;
+
+/** 
+ *
+ * Constructor
+ * 
+ */
+
+BrookShakeAlgorithm::BrookShakeAlgorithm( ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookShakeAlgorithm::BrookShakeAlgorithm";
+
+   BrookOpenMMFloat zero                    = (BrookOpenMMFloat)  0.0;
+   BrookOpenMMFloat one                     = (BrookOpenMMFloat)  1.0;
+
+// ---------------------------------------------------------------------------------------
+
+   _numberOfAtoms                = -1;
+   _numberOfConstraints          = -1;
+
+   // mark stream dimension variables as unset
+
+   _shakeAtomStreamWidth         = -1;
+   _shakeAtomStreamHeight        = -1;
+   _shakeAtomStreamSize          = -1;
+
+   _shakeConstraintStreamSize    = -1;
+   _shakeConstraintStreamWidth   = -1;
+   _shakeConstraintStreamHeight  = -1;
+
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      _shakeStreams[ii]   = NULL;
+   }
+
+   _inverseHydrogenMass    = one/( (BrookOpenMMFloat) 1.008);
+
+   // setup inverse sqrt masses
+
+   _inverseSqrtMasses = NULL;
+
+}   
+ 
+/** 
+ * Destructor
+ * 
+ */
+
+BrookShakeAlgorithm::~BrookShakeAlgorithm( ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookShakeAlgorithm::~BrookShakeAlgorithm";
+
+// ---------------------------------------------------------------------------------------
+
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      delete _shakeStreams[ii];
+   }
+
+   delete[] _inverseSqrtMasses;
+
+}
+
+/** 
+ * Get number of constraints
+ * 
+ * @return   number of constraints
+ *
+ */
+
+int BrookShakeAlgorithm::getNumberOfConstraints( void ) const {
+   return _numberOfConstraints;
+}
+
+/** 
+ * Get inverse of hydrogen mass
+ * 
+ * @return  inverse of hydrogen mass
+ *
+ */
+
+BrookOpenMMFloat BrookShakeAlgorithm::getInverseHydrogenMass( void ) const {
+   return _inverseHydrogenMass;
+}
+
+/** 
+ * Get Atom stream size
+ *
+ * @return  Atom stream size
+ *
+ */
+
+int BrookShakeAlgorithm::getShakeAtomStreamSize( void ) const {
+   return _shakeAtomStreamSize;
+}
+
+/** 
+ * Get atom stream width
+ *
+ * @return  atom stream width
+ *
+ */
+
+int BrookShakeAlgorithm::getShakeAtomStreamWidth( void ) const {
+   return _shakeAtomStreamWidth;
+}
+
+/** 
+ * Get atom stream height
+ *
+ * @return atom stream height
+ */
+
+int BrookShakeAlgorithm::getShakeAtomStreamHeight( void ) const {
+   return _shakeAtomStreamHeight;
+}
+
+/** 
+ * Get Constraint stream size
+ *
+ * @return  Constraint stream size
+ *
+ */
+
+int BrookShakeAlgorithm::getShakeConstraintStreamSize( void ) const {
+   return _shakeConstraintStreamSize;
+}
+
+/** 
+ * Get constraint stream width
+ *
+ * @return  constraint stream width
+ *
+ */
+
+int BrookShakeAlgorithm::getShakeConstraintStreamWidth( void ) const {
+   return _shakeConstraintStreamWidth;
+}
+
+/** 
+ * Get constraint stream height
+ *
+ * @return constraint stream height
+ */
+
+int BrookShakeAlgorithm::getShakeConstraintStreamHeight( void ) const {
+   return _shakeConstraintStreamHeight;
+}
+
+/** 
+ * Get Shake atom indices stream 
+ *
+ * @return  Shake atom indices stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookShakeAlgorithm::getShakeAtomIndicesStream( void ) const {
+   return _shakeStreams[ShakeAtomIndicesStream];
+}
+
+/** 
+ * Get Shake atom parameter stream
+ *
+ * @return  Shake atom parameter sStream
+ *
+ */
+
+BrookFloatStreamInternal* BrookShakeAlgorithm::getShakeAtomParameterStream( void ) const {
+   return _shakeStreams[ShakeAtomParameterStream];
+}
+
+/** 
+ * Get Shake XCons0 stream
+ *
+ * @return  XCons0 stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookShakeAlgorithm::getShakeXCons0Stream( void ) const {
+   return _shakeStreams[ShakeXCons0Stream];
+}
+
+/** 
+ * Get Shake XCons1 stream
+ *
+ * @return  XCons1 stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookShakeAlgorithm::getShakeXCons1Stream( void ) const {
+   return _shakeStreams[ShakeXCons1Stream];
+}
+
+/** 
+ * Get Shake XCons2 stream
+ *
+ * @return  XCons2 stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookShakeAlgorithm::getShakeXCons2Stream( void ) const {
+   return _shakeStreams[ShakeXCons2Stream];
+}
+
+/** 
+ * Get Shake XCons3 stream
+ *
+ * @return  XCons3 stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookShakeAlgorithm::getShakeXCons3Stream( void ) const {
+   return _shakeStreams[ShakeXCons3Stream];
+}
+
+/** 
+ * Get Shake inverse map stream
+ *
+ * @return  Shake inverse map stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookShakeAlgorithm::getShakeInverseMapStream( void ) const {
+   return _shakeStreams[ShakeInverseMapStream];
+}
+
+/** 
+ * Initialize stream dimensions
+ * 
+ * @param numberOfAtoms             number of atoms
+ * @param numberOfConstraints       number of constraints
+ * @param platform                  platform
+ *      
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookShakeAlgorithm::_initializeStreamSizes( int numberOfAtoms, int numberOfConstraints,
+                                                 const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookShakeAlgorithm::_initializeStreamSizes";
+
+// ---------------------------------------------------------------------------------------
+
+   _shakeAtomStreamSize           = getAtomStreamSize( platform );
+   _shakeAtomStreamWidth          = getAtomStreamWidth( platform );
+   _shakeAtomStreamHeight         = getAtomStreamHeight( platform );
+
+   // get constraint stream width & height, and then set stream size to the product
+
+   BrookCommon::getStreamDimensions( numberOfConstraints, &_shakeConstraintStreamWidth, &_shakeConstraintStreamHeight );
+   _shakeConstraintStreamSize     = _shakeConstraintStreamWidth*_shakeConstraintStreamHeight;
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Initialize streams
+ * 
+ * @param platform                  platform
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookShakeAlgorithm::_initializeStreams( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookShakeAlgorithm::_initializeStreams";
+
+   BrookOpenMMFloat dangleValue              = (BrookOpenMMFloat) 0.0;
+
+// ---------------------------------------------------------------------------------------
+
+   int shakeAtomStreamSize                   = getShakeAtomStreamSize();
+   int shakeAtomStreamWidth                  = getShakeAtomStreamWidth();
+
+   int shakeConstraintStreamSize             = getShakeConstraintStreamSize();
+   int shakeConstraintStreamWidth            = getShakeConstraintStreamWidth();
+
+    _shakeStreams[ShakeAtomIndicesStream]    = new BrookFloatStreamInternal( BrookCommon::ShakeAtomIndicesStream,
+                                                                             shakeConstraintStreamSize, shakeConstraintStreamWidth,
+                                                                             BrookStreamInternal::Float4, dangleValue );
+
+    _shakeStreams[ShakeAtomParameterStream]  = new BrookFloatStreamInternal( BrookCommon::ShakeAtomParameterStream,
+                                                                             shakeConstraintStreamSize, shakeConstraintStreamWidth,
+                                                                             BrookStreamInternal::Float4, dangleValue );
+
+    _shakeStreams[ShakeXCons0Stream]         = new BrookFloatStreamInternal( BrookCommon::ShakeXCons0Stream,
+                                                                             shakeConstraintStreamSize, shakeConstraintStreamWidth,
+                                                                             BrookStreamInternal::Float4, dangleValue );
+
+    _shakeStreams[ShakeXCons1Stream]         = new BrookFloatStreamInternal( BrookCommon::ShakeXCons1Stream,
+                                                                             shakeConstraintStreamSize, shakeConstraintStreamWidth,
+                                                                             BrookStreamInternal::Float4, dangleValue );
+
+    _shakeStreams[ShakeXCons2Stream]         = new BrookFloatStreamInternal( BrookCommon::ShakeXCons2Stream,
+                                                                             shakeConstraintStreamSize, shakeConstraintStreamWidth,
+                                                                             BrookStreamInternal::Float4, dangleValue );
+
+    _shakeStreams[ShakeXCons3Stream]         = new BrookFloatStreamInternal( BrookCommon::ShakeXCons3Stream,
+                                                                             shakeConstraintStreamSize, shakeConstraintStreamWidth,
+                                                                             BrookStreamInternal::Float4, dangleValue );
+
+    _shakeStreams[ShakeInverseMapStream]     = new BrookFloatStreamInternal( BrookCommon::ShakeInverseMapStream,
+                                                                             shakeAtomStreamSize, shakeAtomStreamWidth,
+                                                                             BrookStreamInternal::Float2, dangleValue );
+
+   return DefaultReturnValue;
+}
+ 
+/*  
+ * Set Shake streams
+ *
+ * @param masses                masses
+ * @param constraintIndices     constraint atom indices
+ * @param constraintLengths     constraint lengths
+ *
+ * @return ErrorReturnValue if error
+ *
+ * @throw OpenMMException if constraintIndices.size() != constraintLengths.size()
+ *
+ */
+    
+int BrookShakeAlgorithm::_setShakeStreams( const std::vector<double>& masses, const std::vector< std::vector<int> >& constraintIndices,
+                                           const std::vector<double>& constraintLengths ){
+    
+// ---------------------------------------------------------------------------------------
+
+   BrookOpenMMFloat one                      = (BrookOpenMMFloat)  1.0;
+   BrookOpenMMFloat half                     = (BrookOpenMMFloat)  0.5;
+
+   static const std::string methodName       = "BrookShakeAlgorithm::_updateSdStreams";
+
+// ---------------------------------------------------------------------------------------
+
+   int shakeAtomStreamSize                   = getShakeAtomStreamSize();
+   int shakeConstraintStreamSize             = getShakeConstraintStreamSize();
+
+   // check that number of constraints for two input vectors is consistent
+
+   if( constraintIndices.size() != constraintLengths.size() ){
+      std::stringstream message;
+      message << methodName << " constraintIndices size=" << constraintIndices.size() << " does not equal constraintLengths size=" << constraintLengths.size();
+      throw OpenMMException( message.str() );
+      return ErrorReturnValue;
+   }   
+
+   // allocate arrays to be read down to board
+
+   BrookOpenMMFloat* atomIndices      = new BrookOpenMMFloat[4*shakeConstraintStreamSize];
+   BrookOpenMMFloat* shakeParameters  = new BrookOpenMMFloat[4*shakeConstraintStreamSize];
+   for( int ii = 0; ii < 4*shakeConstraintStreamSize; ii++ ){
+      atomIndices[ii] = (BrookOpenMMFloat) -1;
+   }
+   memset( shakeParameters, 0, 4*shakeConstraintStreamSize*sizeof( BrookOpenMMFloat ) ); 
+
+   std::vector< std::vector<int> >::const_iterator atomIterator = constraintIndices.begin();
+   std::vector<double>::const_iterator distanceIterator         = constraintLengths.begin();
+
+   int constraintIndex                = 0;
+   while( atomIterator != constraintIndices.end() ){
+
+      std::vector<int> atomVector = *atomIterator;
+
+      // check that array of indices is not too small or large
+
+      if( atomVector.size() < 2 ){
+         std::stringstream message;
+         message << methodName << " atomIndices size=" << atomVector.size() << " is too small at constraintIndex=" << (constraintIndex/4);
+         throw OpenMMException( message.str() );
+         return ErrorReturnValue;
+      }   
+   
+      if( atomVector.size() > 4 ){
+         std::stringstream message;
+         message << methodName << " atomIndices size=" << atomVector.size() << " is too large at constraintIndex=" << (constraintIndex/4);
+         throw OpenMMException( message.str() );
+         return ErrorReturnValue;
+      }   
+   
+      int index                    = 0;
+      int atomIndex1               = -1;
+      int atomIndex2               = -1;
+      for( std::vector<int>::const_iterator ii = atomVector.begin(); ii != atomVector.end(); ii++, index++ ){
+         atomIndices[constraintIndex + index]   = (BrookOpenMMFloat) *ii;
+         if( index == 0 ){
+            atomIndex1 = *ii;
+         } else if( index == 1 ){
+            atomIndex2 = *ii;
+         }
+      }
+
+      // insure heavy atom is first
+
+      if( masses[atomIndex1] < masses[atomIndex2] ){
+
+         BrookOpenMMFloat swap           = atomIndices[constraintIndex];
+         atomIndices[constraintIndex]    = atomIndices[constraintIndex+1];
+         atomIndices[constraintIndex+1]  = swap;
+
+         int swapI                       = atomIndex1;
+         atomIndex1                      = atomIndex2;
+         atomIndex2                      = swapI;
+      }
+
+      // set parameters:
+
+      //    (1) 1/(heavy atom mass)
+      //    (2) 0.5/(heavy atom mass+hydrogen mass)
+      //    (3) constraint distance**2
+
+      shakeParameters[constraintIndex]    = one/( (BrookOpenMMFloat) masses[atomIndex1] );
+      shakeParameters[constraintIndex+1]  = half/( (BrookOpenMMFloat) (masses[atomIndex1] + masses[atomIndex2]) );
+      shakeParameters[constraintIndex+2]  = (BrookOpenMMFloat) ( (*distanceIterator)*(*distanceIterator) );
+
+      atomIterator++;
+      distanceIterator++;
+      constraintIndex += 4;
+   }
+
+   // write entries to board
+
+   _shakeStreams[ShakeAtomIndicesStream]->loadFromArray( atomIndices );
+   _shakeStreams[ShakeAtomParameterStream]->loadFromArray( shakeParameters );
+
+   delete[] shakeParameters;
+
+   // initialize inverse map
+
+   BrookOpenMMFloat* inverseMap = new BrookOpenMMFloat[2*shakeAtomStreamSize];
+   for( int ii = 0; ii < shakeAtomStreamSize*2; ii++ ){
+      inverseMap[ii] = -1;
+   }
+   
+   // build inverse map
+ 
+   for( int ii = 0; ii < shakeConstraintStreamSize; ii++ ){
+      int ii4 = ii << 2;
+      for( int jj = 0; jj < 4; jj++ ){
+         if( atomIndices[ii4+jj] != -1 ){
+            int atomIndex             = (int) (atomIndices[ii4+jj] + 0.001);
+            inverseMap[atomIndex*2]   = (float) ii;
+            inverseMap[atomIndex*2+1] = (float) jj;
+         }
+      }
+   }
+   
+   _shakeStreams[ShakeInverseMapStream]->loadFromArray( inverseMap );
+
+   delete[] atomIndices;
+   delete[] inverseMap;
+
+   return DefaultReturnValue;
+
+}
+ 
+/*  
+ * Setup of Shake parameters
+ *
+ * @param masses                masses
+ * @param constraintIndices     constraint atom indices
+ * @param constraintLengths     constraint lengths
+ * @param platform              Brook platform
+ *
+ * @return ErrorReturnValue if error
+ *
+ */
+    
+int BrookShakeAlgorithm::setup( const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
+                                const std::vector<double>& constraintLengths, const Platform& platform ){
+    
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookShakeAlgorithm::setup";
+
+// ---------------------------------------------------------------------------------------
+
+   int numberOfAtoms  = (int) masses.size();
+   setNumberOfAtoms( numberOfAtoms );
+
+   // set stream sizes and then create streams
+
+   _initializeStreamSizes( numberOfAtoms, (int) constraintIndices.size(), platform );
+   _initializeStreams( platform );
+
+   _setShakeStreams( masses, constraintIndices, constraintLengths );
+
+   return DefaultReturnValue;
+}
+
+/* 
+ * Get contents of object
+ *
+ * @param level   level of dump
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string BrookShakeAlgorithm::getContentsString( int level ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookShakeAlgorithm::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", getAtomStreamWidth() );
+   message << _getLine( tab, "Atom stream width:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamHeight() );
+   message << _getLine( tab, "Atom stream height:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamSize() );
+   message << _getLine( tab, "Atom stream size:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfConstraints() );
+   message << _getLine( tab, "Number of constraints:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getShakeConstraintStreamWidth() );
+   message << _getLine( tab, "Constraint stream width:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getShakeConstraintStreamHeight() );
+   message << _getLine( tab, "Constraint stream height:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getShakeConstraintStreamSize() );
+   message << _getLine( tab, "Constraint stream size:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5f", 1.0f/getInverseHydrogenMass() );
+   message << _getLine( tab, "H-mass:", value ); 
+
+   message << _getLine( tab, "Log:",                  (getLog()                          ? Set : NotSet) ); 
+
+   message << _getLine( tab, "AtomIndices:",          (getShakeAtomIndicesStream()       ? Set : NotSet) ); 
+   message << _getLine( tab, "AtomParameters:",       (getShakeAtomParameterStream()     ? Set : NotSet) ); 
+   message << _getLine( tab, "XCons0:",               (getShakeXCons0Stream()            ? Set : NotSet) ); 
+   message << _getLine( tab, "XCons1:",               (getShakeXCons1Stream()            ? Set : NotSet) ); 
+   message << _getLine( tab, "XCons2:",               (getShakeXCons2Stream()            ? Set : NotSet) ); 
+   message << _getLine( tab, "XCons3:",               (getShakeXCons3Stream()            ? Set : NotSet) ); 
+   message << _getLine( tab, "InverseMap:",           (getShakeInverseMapStream()        ? Set : NotSet) ); 
+ 
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      message << std::endl;
+      if( _shakeStreams[ii] ){
+         message << _shakeStreams[ii]->getContentsString( );
+      }
+   }
+
+#undef LOCAL_SPRINTF
+
+   return message.str();
+}

platforms/brook/src/BrookShakeAlgorithm.h

+#ifndef OPENMM_BROOK_SHAKE_H_
+#define OPENMM_BROOK_SHAKE_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 stochastic dynamics algorithm 
+ *
+ */
+
+class BrookShakeAlgorithm : public BrookCommon {
+
+   public:
+  
+      /** 
+       * Constructor
+       * 
+       */
+      
+      BrookShakeAlgorithm( );
+  
+      /** 
+       * Destructor
+       * 
+       */
+      
+      ~BrookShakeAlgorithm();
+  
+      /** 
+       * Get number of constraints
+       * 
+       * @return   number of constraints
+       *
+       */
+      
+      int getNumberOfConstraints( void ) const;
+      
+      /**
+       * Get inverse of hydrogen mass
+       *
+       * @return inverse of hydrogen mass
+       */
+
+      BrookOpenMMFloat getInverseHydrogenMass( void ) const; 
+
+      /**
+       * Get Shake atom stream width
+       *
+       * @return atom stream width
+       */
+
+      int getShakeAtomStreamWidth( void ) const; 
+
+      /**
+       * Get Shake atom stream height
+       *
+       * @return atom stream height
+       */
+
+      int getShakeAtomStreamHeight( void ) const;
+
+      /**
+       * Get Shake atom stream size
+       * 
+       * @return atom stream size
+       */
+
+      int getShakeAtomStreamSize( void ) const; 
+
+      /**
+       * Get Shake constraint stream width
+       *
+       * @return constraint stream width
+       */
+
+      int getShakeConstraintStreamWidth( void ) const; 
+
+      /**
+       * Get Shake constraint stream height
+       *
+       * @return constraint stream height
+       */
+
+      int getShakeConstraintStreamHeight( void ) const;
+
+      /**
+       * Get Shake constraint stream size
+       * 
+       * @return constraint stream size
+       */
+
+      int getShakeConstraintStreamSize( void ) const; 
+
+      /** 
+       * Get array of Shake streams 
+       *
+       * @return  array ofstreams
+       *
+       */
+      
+      BrookFloatStreamInternal** getStreams( void );
+      
+      /*  
+       * Setup of Shake parameters
+       *
+       * @param masses                masses
+       * @param constraintIndices     constraint atom indices
+       * @param constraintLengths     constraint lengths
+       * @param platform              Brook platform
+       *
+       * @return ErrorReturnValue if error
+       *
+       */
+              
+      int setup( const std::vector<double>& masses, const std::vector< std::vector<int> >& constraintIndices,
+                 const std::vector<double>& constraintLengths, const Platform& platform );
+          
+      /* 
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string getContentsString( int level = 0 ) const;
+
+      /** 
+       * Get Shake atom indices stream
+       *
+       * @return  Shake atom indices stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getShakeAtomIndicesStream( void ) const;
+      
+      /** 
+       * Get  Shake atom parameter stream
+       *
+       * @return   Shake atom parameter stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getShakeAtomParameterStream( void ) const;
+      
+      /** 
+       * Get XCons0 stream
+       *
+       * @return  XCons0 stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getShakeXCons0Stream( void ) const;
+      
+      /** 
+       * Get XCons1 stream
+       *
+       * @return  XCons1 stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getShakeXCons1Stream( void ) const;
+      
+      /** 
+       * Get XCons2 stream
+       *
+       * @return  XCons2 stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getShakeXCons2Stream( void ) const;
+      
+      /** 
+       * Get XCons3 stream
+       *
+       * @return  XCons3 stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getShakeXCons3Stream( void ) const;
+      
+      /** 
+       * Get Shake inverse map stream
+       *
+       * @return  Shake inverse map stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getShakeInverseMapStream( void ) const;
+
+   private:
+   
+      // streams indices
+
+      enum BrookShakeAlgorithmStreams { 
+              ShakeAtomIndicesStream,
+              ShakeAtomParameterStream,
+              ShakeXCons0Stream,
+              ShakeXCons1Stream,
+              ShakeXCons2Stream,
+              ShakeXCons3Stream,
+              ShakeInverseMapStream,
+              LastStreamIndex
+           };
+
+      // number of constraints
+
+      int _numberOfConstraints;
+
+      // inverse of H mass
+
+      BrookOpenMMFloat _inverseHydrogenMass;
+
+      // atom stream dimensions
+
+      int _shakeAtomStreamWidth;
+      int _shakeAtomStreamHeight;
+      int _shakeAtomStreamSize;
+
+      // constraint stream dimensions
+
+      int _shakeConstraintStreamSize;
+      int _shakeConstraintStreamWidth;
+      int _shakeConstraintStreamHeight;
+
+      // inverse sqrt masses
+
+      BrookOpenMMFloat* _inverseSqrtMasses;
+
+      // internal streams
+
+      BrookFloatStreamInternal* _shakeStreams[LastStreamIndex];
+
+      /* 
+       * Setup of stream dimensions
+       *
+       * @param atomStreamSize        atom stream size
+       * @param atomStreamWidth       atom stream width
+       *
+       * @return ErrorReturnValue if error, else DefaultReturnValueValue
+       *
+       * */
+      
+      int _initializeStreamSizes( int atomStreamSize, int atomStreamWidth );
+
+      /** 
+       * Initialize stream dimensions
+       * 
+       * @param numberOfAtoms             number of atoms
+       * @param numberOfConstraints       number of constraints
+       * @param platform                  platform
+       *
+       * @return ErrorReturnValue if error, else DefaultReturnValueValue
+       *
+       */
+      
+      int _initializeStreamSizes( int numberOfAtoms, int numberOfConstraints, const Platform& platform );
+      
+      /** 
+       * Initialize stream dimensions and streams
+       * 
+       * @param platform                  platform
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int _initializeStreams( const Platform& platform );
+
+      /*  
+       * Set Shake streams
+       *
+       * @param masses                masses
+       * @param constraintIndices     constraint atom indices
+       * @param constraintLengths     constraint lengths
+       *
+       * @return ErrorReturnValue if error
+       *
+       * @throw OpenMMException if constraintIndices.size() != constraintLengths.size()
+       *
+       */
+           
+      int _setShakeStreams( const std::vector<double>& masses, const std::vector< std::vector<int> >& constraintIndices,
+                            const std::vector<double>& constraintLengths );
+          
+};
+
+} // namespace OpenMM
+
+#endif /* OPENMM_BROOK_SHAKE_H_ */

platforms/brook/src/BrookStochasticDynamics.cpp

@@ -52,7 +52,7 @@ BrookStochasticDynamics::BrookStochasticDynamics( ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookStochasticDynamics::BrookStochasticDynamics";
+   //static const std::string methodName      = "BrookStochasticDynamics::BrookStochasticDynamics";
 
    BrookOpenMMFloat zero                    = (BrookOpenMMFloat)  0.0;
    BrookOpenMMFloat one                     = (BrookOpenMMFloat)  1.0;
@@ -76,7 +76,9 @@ BrookStochasticDynamics::BrookStochasticDynamics( ){
       _derivedParameters[ii]   = oneMinus;
    }
 
-   _temperature = _stepSize = _tau = oneMinus;
+   _temperature = oneMinus;
+   _stepSize    = oneMinus;
+   _tau         = oneMinus;
 
    // setup inverse sqrt masses
 
@@ -84,6 +86,8 @@ BrookStochasticDynamics::BrookStochasticDynamics( ){
 
    // set randomNumber seed 
 
+   _randomNumberSeed  = 1393;
+
    //_randomNumberSeed = randomNumberSeed ? randomNumberSeed : 1393;
    //SimTKOpenMMUtilities::setRandomNumberSeed( randomNumberSeed );
 }   
@@ -120,6 +124,19 @@ BrookOpenMMFloat BrookStochasticDynamics::getTau( void ) const {
    return _tau;
 }
 
+/** 
+ * Get friction
+ * 
+ * @return  friction
+ *
+ */
+
+BrookOpenMMFloat BrookStochasticDynamics::getFriction( void ) const {
+   static const BrookOpenMMFloat zero = (BrookOpenMMFloat) 0.0; 
+   static const BrookOpenMMFloat one  = (BrookOpenMMFloat) 1.0; 
+   return ( (_tau == zero) ? zero : (one/_tau) );
+}
+
 /** 
  * Get temperature
  * 
@@ -147,13 +164,13 @@ BrookOpenMMFloat BrookStochasticDynamics::getStepSize( void ) const {
  * 
  * @param tau   new tau value
  *
- * @return      DefaultReturn
+ * @return      DefaultReturnValue
  *
  */
 
 int BrookStochasticDynamics::_setTau( BrookOpenMMFloat tau ){
    _tau = tau;
-   return DefaultReturn;
+   return DefaultReturnValue;
 }
 
 /** 
@@ -161,13 +178,13 @@ int BrookStochasticDynamics::_setTau( BrookOpenMMFloat tau ){
  * 
  * @param friction   new friction value
  *
- * @return      DefaultReturn
+ * @return      DefaultReturnValue
  *
  */
 
 int BrookStochasticDynamics::_setFriction( BrookOpenMMFloat friction ){
    _tau   = (BrookOpenMMFloat) ( (friction != 0.0) ? 1.0/friction : 0.0);
-   return DefaultReturn;
+   return DefaultReturnValue;
 }
 
 /** 
@@ -175,13 +192,13 @@ int BrookStochasticDynamics::_setFriction( BrookOpenMMFloat friction ){
  * 
  * @parameter   temperature
  *
- * @return      DefaultReturn
+ * @return      DefaultReturnValue
  *
  */
 
 int BrookStochasticDynamics::_setTemperature( BrookOpenMMFloat temperature ){
    _temperature = temperature;
-   return DefaultReturn;
+   return DefaultReturnValue;
 }
 
 /** 
@@ -189,19 +206,19 @@ int BrookStochasticDynamics::_setTemperature( BrookOpenMMFloat temperature ){
  * 
  * @param   stepSize
  *
- * @return      DefaultReturn
+ * @return      DefaultReturnValue
  *
  */
 
 int BrookStochasticDynamics::_setStepSize( BrookOpenMMFloat stepSize ){
    _stepSize = stepSize;
-   return DefaultReturn;
+   return DefaultReturnValue;
 }
 
 /** 
  * Update derived parameters
  * 
- * @return  DefaultReturn
+ * @return  DefaultReturnValue
  *
  */
 
@@ -272,7 +289,7 @@ int BrookStochasticDynamics::_updateDerivedParameters( void ){
    _derivedParameters[Sd2pc1] = one/_derivedParameters[Sd1pc2];
    _derivedParameters[Sd2pc2] = tau*_derivedParameters[D]/( _derivedParameters[EM] - one );
        
-   return DefaultReturn;
+   return DefaultReturnValue;
 
 };
 
@@ -295,24 +312,51 @@ int BrookStochasticDynamics::updateParameters( double temperature, double fricti
 
    // ---------------------------------------------------------------------------------------
 
-   _setStepSize( stepSize );
-   _setFriction( friction );
-   _setTemperature( temperature );
+   _setStepSize(    (BrookOpenMMFloat)  stepSize );
+   _setFriction(    (BrookOpenMMFloat)  friction );
+   _setTemperature( (BrookOpenMMFloat)  temperature );
 
    _updateDerivedParameters( );
    _updateSdStreams( );
 
-   return DefaultReturn;
+   return DefaultReturnValue;
 
 };
 
-/**---------------------------------------------------------------------------------------
+/** 
+ * Update
+ * 
+ * @param  positions           atom positions
+ * @param  velocities          atom velocities
+ * @param  forces              atom forces
+ * @param  brookShakeAlgorithm BrookShakeAlgorithm reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
+                                     const Stream& forces,
+                                     BrookShakeAlgorithm& brookShakeAlgorithm ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nBrookStochasticDynamics::update";
+
+   // ---------------------------------------------------------------------------------------
 
-   Get array of derived parameters indexed by 'DerivedParameters' enums
 
-   @return array
+   return DefaultReturnValue;
+
+};
 
-   --------------------------------------------------------------------------------------- */
+/**
+ *
+ * Get array of derived parameters indexed by 'DerivedParameters' enums
+ *
+ * @return array
+ *
+ */
    
 const BrookOpenMMFloat* BrookStochasticDynamics::getDerivedParameters( void ) const {
 
@@ -364,7 +408,7 @@ int BrookStochasticDynamics::getStochasticDynamicsAtomStreamHeight( void ) const
  *
  */
 
-BrookFloatStreamInternal* BrookStochasticDynamics::getSDPC1( void ) const {
+BrookFloatStreamInternal* BrookStochasticDynamics::getSDPC1Stream( void ) const {
    return _sdStreams[SDPC1Stream];
 }
 
@@ -375,7 +419,7 @@ BrookFloatStreamInternal* BrookStochasticDynamics::getSDPC1( void ) const {
  *
  */
 
-BrookFloatStreamInternal* BrookStochasticDynamics::getSDPC2( void ) const {
+BrookFloatStreamInternal* BrookStochasticDynamics::getSDPC2Stream( void ) const {
    return _sdStreams[SDPC2Stream];
 }
 
@@ -386,7 +430,7 @@ BrookFloatStreamInternal* BrookStochasticDynamics::getSDPC2( void ) const {
  *
  */
 
-BrookFloatStreamInternal* BrookStochasticDynamics::getSD2X( void ) const {
+BrookFloatStreamInternal* BrookStochasticDynamics::getSD2XStream( void ) const {
    return _sdStreams[SD2XStream];
 }
 
@@ -397,7 +441,7 @@ BrookFloatStreamInternal* BrookStochasticDynamics::getSD2X( void ) const {
  *
  */
 
-BrookFloatStreamInternal* BrookStochasticDynamics::getSD1V( void ) const {
+BrookFloatStreamInternal* BrookStochasticDynamics::getSD1VStream( void ) const {
    return _sdStreams[SD1VStream];
 }
 
@@ -436,7 +480,8 @@ int BrookStochasticDynamics::_initializeStreamSizes( int numberOfAtoms, const Pl
  *
  */
 
-std::string BrookStochasticDynamics::_getDerivedParametersString( BrookStochasticDynamics::DerivedParameters  derivedParametersIndex ) const {
+//std::string BrookStochasticDynamics::_getDerivedParametersString( BrookStochasticDynamics::DerivedParameters  derivedParametersIndex ) const {
+std::string BrookStochasticDynamics::_getDerivedParametersString( int derivedParametersIndex ) const {
 
 // ---------------------------------------------------------------------------------------
 
@@ -539,6 +584,8 @@ int BrookStochasticDynamics::_initializeStreams( const Platform& platform ){
 
    //static const std::string methodName      = "BrookStochasticDynamics::_initializeStreams";
 
+   BrookOpenMMFloat dangleValue            = (BrookOpenMMFloat) 0.0;
+
 // ---------------------------------------------------------------------------------------
 
    int sdAtomStreamSize             = getStochasticDynamicsAtomStreamSize();
@@ -580,7 +627,7 @@ int BrookStochasticDynamics::_updateSdStreams( void ){
 
    int sdAtomStreamSize                      = getStochasticDynamicsAtomStreamSize();
 
-   BrookOpenMMFloat sdpc[2];
+   BrookOpenMMFloat* sdpc[2];
    for( int ii = 0; ii < 2; ii++ ){
       sdpc[ii] = new BrookOpenMMFloat[2*sdAtomStreamSize];
       memset( sdpc[ii], 0, 2*sdAtomStreamSize*sizeof( BrookOpenMMFloat ) ); 
@@ -589,7 +636,7 @@ int BrookStochasticDynamics::_updateSdStreams( void ){
    const BrookOpenMMFloat* derivedParameters = getDerivedParameters( );
    int numberOfAtoms                         = getNumberOfAtoms();
    int index                                 = 0;
-   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+   for( int ii = 0; ii < numberOfAtoms; ii++, index += 2 ){
 
       sdpc[0][index]      = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yv]) );
       sdpc[0][index+1]    = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[V])  );
@@ -597,7 +644,6 @@ int BrookStochasticDynamics::_updateSdStreams( void ){
       sdpc[1][index]      = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[Yx]) );
       sdpc[1][index+1]    = _inverseSqrtMasses[ii]*( static_cast<BrookOpenMMFloat> (derivedParameters[X])  );
 
-      index              += 2;
    }
 
    _sdStreams[SDPC1Stream]->loadFromArray( sdpc[0] );
@@ -609,12 +655,13 @@ int BrookStochasticDynamics::_updateSdStreams( void ){
 
    // initialize SD2X
 
-   sd2x = new BrookOpenMMFloat[3*sdAtomStreamSize];
-   SimTKOpenMMUtilities::setRandomNumberSeed( (uint32_t) getRandomNumberSeed() );
+   BrookOpenMMFloat* sd2x = new BrookOpenMMFloat[3*sdAtomStreamSize];
+   //SimTKOpenMMUtilities::setRandomNumberSeed( (uint32_t) getRandomNumberSeed() );
 
    memset( sd2x, 0, 3*sdAtomStreamSize*sizeof( BrookOpenMMFloat ) ); 
 
-   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+   index = 0;
+   for( int ii = 0; ii < numberOfAtoms; 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()) );
@@ -639,7 +686,10 @@ int BrookStochasticDynamics::_setInverseSqrtMasses( const std::vector<double>& m
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookStochasticDynamics::_setInverseSqrtMasses";
+   //static const std::string methodName      = "BrookStochasticDynamics::_setInverseSqrtMasses";
+
+   BrookOpenMMFloat zero                    = (BrookOpenMMFloat)  0.0;
+   BrookOpenMMFloat one                     = (BrookOpenMMFloat)  1.0;
 
 // ---------------------------------------------------------------------------------------
 
@@ -647,9 +697,10 @@ int BrookStochasticDynamics::_setInverseSqrtMasses( const std::vector<double>& m
 
    _inverseSqrtMasses = new BrookOpenMMFloat[masses.size()];
    int index          = 0;
-   for( std::vector<double>::const_interator ii = masses.begin(); ii != masses.end(); ii++, index++ ){
+   for( std::vector<double>::const_iterator ii = masses.begin(); ii != masses.end(); ii++, index++ ){
       if( *ii != 0.0 ){
-         _inverseSqrtMasses[index] = ( SQRT( one/(*ii) ) );
+         BrookOpenMMFloat value    = static_cast<BrookOpenMMFloat>(*ii);
+         _inverseSqrtMasses[index] = ( SQRT( one/value ) );
       } else {
          _inverseSqrtMasses[index] = zero;
       }
@@ -668,7 +719,7 @@ int BrookStochasticDynamics::_setInverseSqrtMasses( const std::vector<double>& m
  *
  * */
     
-int BrookStochasticDynamics::setup( const std::vector<<double> >& masses, const Platform& platform ){
+int BrookStochasticDynamics::setup( const std::vector<double>& masses, const Platform& platform ){
     
 // ---------------------------------------------------------------------------------------
 
@@ -679,58 +730,19 @@ int BrookStochasticDynamics::setup( const std::vector<<double> >& masses, const
    int numberOfAtoms  = (int) masses.size();
    setNumberOfAtoms( numberOfAtoms );
 
-   _setInverseSqrtMasses( masses );
-
-   return DefaultReturnValue;
-}
-
-/* 
- * Setup of stream dimensions for partial force streams
- *
- * @param atomStreamSize        atom stream size
- * @param atomStreamWidth       atom stream width
- *
- * @return ErrorReturnValue if error, else DefaultReturnValue
- *
- * */
-
-int BrookStochasticDynamics::initializePartialForceStreamSize( int atomStreamSize, int atomStreamWidth ){
-
-// ---------------------------------------------------------------------------------------
-
-   static const std::string methodName      = "BrookStochasticDynamics::initializePartialForceStreamSize";
-   //static const int debug                   = 1;
-
-// ---------------------------------------------------------------------------------------
-
-   int innerUnroll           = getInnerLoopUnroll();
-   if( innerUnroll < 1 ){
-      std::stringstream message;
-      message << methodName << " innerUnrolls=" << innerUnroll << " is less than 1.";
-      throw OpenMMException( message.str() );
-      return ErrorReturnValue;
-   }
+   // set stream sizes and then create streams
 
-   if( _partialForceStreamWidth < 1 ){
-      std::stringstream message;
-      message << methodName << " partial force stream width=" << _partialForceStreamWidth << " is less than 1.";
-      throw OpenMMException( message.str() );
-      return ErrorReturnValue;
-   }
+   _initializeStreamSizes( numberOfAtoms, platform );
+   _initializeStreams( platform );
 
-   _partialForceStreamSize    = atomStreamSize*getDuplicationFactor()/innerUnroll;
-   _partialForceStreamHeight  = _partialForceStreamSize/_partialForceStreamWidth;
-   _partialForceStreamHeight += ( (_partialForceStreamSize % _partialForceStreamWidth) ? 1 : 0);
+   _setInverseSqrtMasses( masses );
 
    return DefaultReturnValue;
 }
 
 /* 
- * Setup of j-stream dimensions
- * 
  * Get contents of object
  *
- *
  * @param level   level of dump
  *
  * @return string containing contents
@@ -762,21 +774,6 @@ std::string BrookStochasticDynamics::getContentsString( int level ) const {
    (void) LOCAL_SPRINTF( value, "%d", getNumberOfAtoms() );
    message << _getLine( tab, "Number of atoms:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getNumberOfForceStreams() );
-   message << _getLine( tab, "Number of force streams:", value ); 
-
-   (void) LOCAL_SPRINTF( value, "%d", getDuplicationFactor() );
-   message << _getLine( tab, "Duplication factor:", value ); 
-
-   (void) LOCAL_SPRINTF( value, "%d", getInnerLoopUnroll () )
-   message << _getLine( tab, "Inner loop unroll:", value ); 
-
-   (void) LOCAL_SPRINTF( value, "%d", getOuterLoopUnroll() )
-   message << _getLine( tab, "Outer loop unroll:", value ); 
-
-   (void) LOCAL_SPRINTF( value, "%d", getAtomSizeCeiling() );
-   message << _getLine( tab, "Atom ceiling:", value ); 
-
    (void) LOCAL_SPRINTF( value, "%d", getAtomStreamWidth() );
    message << _getLine( tab, "Atom stream width:", value ); 
 
@@ -786,23 +783,30 @@ std::string BrookStochasticDynamics::getContentsString( int level ) const {
    (void) LOCAL_SPRINTF( value, "%d", getAtomStreamSize() );
    message << _getLine( tab, "Atom stream size:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getPartialForceStreamWidth() );
-   message << _getLine( tab, "Partial force stream width:", value ); 
+   (void) LOCAL_SPRINTF( value, "%.5f", getTau() );
+   message << _getLine( tab, "Tau:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5f", getTemperature() );
+   message << _getLine( tab, "Temperature:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getPartialForceStreamHeight() );
-   message << _getLine( tab, "Partial force stream height:", value ); 
+   (void) LOCAL_SPRINTF( value, "%.5f", getStepSize() );
+   message << _getLine( tab, "Step size:", value ); 
+
+   const BrookOpenMMFloat* derivedParameters = getDerivedParameters();
+   for( int ii = 0; ii < MaxDerivedParameters; ii++ ){
+      (void) LOCAL_SPRINTF( value, "%.5e", derivedParameters[ii] );
+      message << _getLine( tab, _getDerivedParametersString( ii ), value ); 
+   }
 
-   (void) LOCAL_SPRINTF( value, "%d", getPartialForceStreamSize() );
-   message << _getLine( tab, "Partial force stream size:", value ); 
+   (void) LOCAL_SPRINTF( value, "%.5f", getTemperature() );
+   message << _getLine( tab, "Temperature:", value ); 
 
    message << _getLine( tab, "Log:",                  (getLog()                ? Set : NotSet) ); 
-/*
-   message << _getLine( tab, "ExclusionStream:",      (getExclusionStream()    ? Set : NotSet) ); 
-   message << _getLine( tab, "VdwStream:",            (getOuterVdwStream()     ? Set : NotSet) ); 
-   message << _getLine( tab, "ChargeStream:",         (getChargeStream()       ? Set : NotSet) ); 
-   message << _getLine( tab, "SigmaStream:",          (getInnerSigmaStream()   ? Set : NotSet) ); 
-   message << _getLine( tab, "EpsilonStream:",        (getInnerEpsilonStream() ? Set : NotSet) ); 
-*/
+
+   message << _getLine( tab, "SDPC1:",                (getSDPC1Stream()        ? Set : NotSet) ); 
+   message << _getLine( tab, "SDPC2:",                (getSDPC2Stream()        ? Set : NotSet) ); 
+   message << _getLine( tab, "SD2X:",                 (getSD2XStream()         ? Set : NotSet) ); 
+   message << _getLine( tab, "SD1V:",                 (getSD1VStream()         ? Set : NotSet) ); 
  
    for( int ii = 0; ii < LastStreamIndex; ii++ ){
       message << std::endl;
@@ -811,21 +815,6 @@ std::string BrookStochasticDynamics::getContentsString( int level ) const {
       }
    }
 
-   // force streams
-
-   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
-      char description[256];
-      (void) LOCAL_SPRINTF( description, "PartialForceStream %d", ii );
-      message << _getLine( tab, description,  (isForceStreamSet(ii) ? Set : NotSet) ); 
-   }
- 
-   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
-      message << std::endl;
-      if( _sdForceStreams[ii] ){
-         message << _sdForceStreams[ii]->getContentsString( );
-      }
-   }
-
 #undef LOCAL_SPRINTF
 
    return message.str();

platforms/brook/src/BrookStochasticDynamics.h

@@ -36,6 +36,7 @@
 #include <set>
 
 #include "BrookFloatStreamInternal.h"
+#include "BrookShakeAlgorithm.h"
 #include "BrookPlatform.h"
 #include "BrookCommon.h"
 
@@ -54,12 +55,9 @@ class BrookStochasticDynamics : public BrookCommon {
       /** 
        * Constructor
        * 
-       * @param masses             atomic masses
-       * @param randomNumberSeed   random number seed
-       *
        */
       
-      BrookStochasticDynamics( const std::vector<double>& masses, uint32_t randomNumberSeed = 1364 );
+      BrookStochasticDynamics(  );
   
       /** 
        * Destructor
@@ -76,6 +74,14 @@ class BrookStochasticDynamics : public BrookCommon {
 
       BrookOpenMMFloat getTau( void ) const; 
 
+      /**
+       * Get friction
+       *
+       * @return friction
+       */
+
+      BrookOpenMMFloat getFriction( void ) const; 
+
       /**
        * Get temperature
        *
@@ -92,6 +98,16 @@ class BrookStochasticDynamics : public BrookCommon {
 
       BrookOpenMMFloat getStepSize( void ) const; 
 
+      /**
+       *
+       * Get array of derived parameters indexed by 'DerivedParameters' enums
+       *
+       * @return array
+       *
+       */
+      
+      const BrookOpenMMFloat* getDerivedParameters( void ) const;
+      
       /**
        * Get StochasticDynamics atom stream width
        *
@@ -121,13 +137,28 @@ class BrookStochasticDynamics : public BrookCommon {
        * 
        * @param  temperature     temperature
        * @param  friction        friction
+       * @param  step size       step size
        *
-       * @return   solute dielectric
+       * @return   DefaultReturnValue
        *
        */
       
-      int updateParameters( double temperature, double friction );
+      int updateParameters( double temperature, double friction, double stepSize );
+      
+      /** 
+       * Update
+       * 
+       * @param  positions           atom positions
+       * @param  velocities          atom velocities
+       * @param  forces              atom forces
+       * @param  brookShakeAlgorithm BrookShakeAlgorithm reference
+       *
+       * @return  DefaultReturnValue
+       *
+       */
       
+      int update( Stream& positions, Stream& velocities,
+                  const Stream& forces, BrookShakeAlgorithm& brookShakeAlgorithm );
       /** 
        * Get array of StochasticDynamics streams 
        *
@@ -140,19 +171,14 @@ class BrookStochasticDynamics : public BrookCommon {
       /* 
        * Setup of StochasticDynamics parameters
        *
-       * @param atomParameters        vector of OBC parameters [atomI][0=charge]
-       *                                                       [atomI][1=radius]
-       *                                                       [atomI][2=scaling factor]
-       * @param solventDielectric     solvent dielectric
-       * @param soluteDielectric      solute dielectric
+       * @param masses                atom masses
        * @param platform              Brook platform
        *
-       * @return nonzero value if error
+       * @return ErrorReturnValue value if error, else DefaultReturnValue
        *
        * */
       
-      int setup( const std::vector<std::vector<double> >& atomParameters, 
-                 double solventDielectric, double soluteDielectric, const Platform& platform  );
+      int setup( const std::vector<double>& masses, const Platform& platform  );
       
       /* 
        * Get contents of object
@@ -172,7 +198,7 @@ class BrookStochasticDynamics : public BrookCommon {
        *
        */
       
-      BrookFloatStreamInternal* getSDPC1( void ) const;
+      BrookFloatStreamInternal* getSDPC1Stream( void ) const;
       
       /** 
        * Get SDPC2 stream 
@@ -181,7 +207,7 @@ class BrookStochasticDynamics : public BrookCommon {
        *
        */
       
-      BrookFloatStreamInternal* getSDPC2( void ) const;
+      BrookFloatStreamInternal* getSDPC2Stream( void ) const;
       
       /** 
        * Get SD2X stream 
@@ -190,7 +216,7 @@ class BrookStochasticDynamics : public BrookCommon {
        *
        */
       
-      BrookFloatStreamInternal* getSD2X( void ) const;
+      BrookFloatStreamInternal* getSD2XStream( void ) const;
       
       /** 
        * Get SD1V stream 
@@ -199,7 +225,7 @@ class BrookStochasticDynamics : public BrookCommon {
        *
        */
       
-      BrookFloatStreamInternal* getSD1V( void ) const;
+      BrookFloatStreamInternal* getSD1VStream( void ) const;
 
    private:
    
@@ -210,7 +236,7 @@ class BrookStochasticDynamics : public BrookCommon {
 
       // streams indices
 
-      enum { 
+      enum BrookStochasticDynamicsStreams { 
               SDPC1Stream,
               SDPC2Stream,
               SD2XStream,
@@ -220,7 +246,7 @@ class BrookStochasticDynamics : public BrookCommon {
 
       // randomNumberSeed
 
-      uint32_t _randomNumberSeed;
+      unsigned int _randomNumberSeed;
 
       BrookOpenMMFloat _tau;
       BrookOpenMMFloat _temperature;
@@ -239,7 +265,8 @@ class BrookStochasticDynamics : public BrookCommon {
        *
        */
 
-      std::string _getDerivedParametersString( BrookStochasticDynamics::DerivedParameters ) const;
+      //std::string _getDerivedParametersString( BrookStochasticDynamics::DerivedParameters ) const;
+      std::string _getDerivedParametersString( int id ) const;
 
       /** 
        * Update derived parameters