Added Verlet & KE code

platforms/brook/src/BrookCalcKineticEnergyKernel.cpp

@@ -29,8 +29,10 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include <sstream>
+#include "OpenMMException.h"
 #include "BrookCalcKineticEnergyKernel.h"
-#include "BrookStreamInternal.h"
+#include "BrookStreamImpl.h"
 
 using namespace OpenMM;
 using namespace std;
@@ -52,6 +54,8 @@ BrookCalcKineticEnergyKernel::BrookCalcKineticEnergyKernel( std::string name, co
 
 // ---------------------------------------------------------------------------------------
 
+   _masses   = NULL;
+
 }
 
 /** 
@@ -67,19 +71,7 @@ BrookCalcKineticEnergyKernel::~BrookCalcKineticEnergyKernel( ){
 
 // ---------------------------------------------------------------------------------------
 
-/*
-    if (dynamics)
-        delete dynamics;
-    if (shake)
-        delete shake;
-    if (masses)
-        delete[] masses;
-    if (constraintIndices)
-        disposeIntArray(constraintIndices, numConstraints);
-    if (shakeParameters)
-        disposeRealArray(shakeParameters, numConstraints);
-*/
-
+   delete[] _masses;
 }
 
 /** 
@@ -97,7 +89,17 @@ void BrookCalcKineticEnergyKernel::initialize( const vector<double>& masses ){
 
 // ---------------------------------------------------------------------------------------
 
-    // this->masses = masses;
+   // masses
+
+   if( _masses ){
+      delete[] _masses;
+   }
+
+   _masses = new BrookOpenMMFloat[masses.size()];
+
+   for( unsigned int ii = 0; ii < masses.size(); ii++ ){
+      _masses[ii]  = static_cast<BrookOpenMMFloat> (masses[ii]);
+   }
 
    return;
 }
@@ -115,17 +117,35 @@ double BrookCalcKineticEnergyKernel::execute( const Stream& velocities ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookCalcKineticEnergyKernel::execute";
+   static const std::string methodName      = "BrookCalcKineticEnergyKernel::execute";
 
 // ---------------------------------------------------------------------------------------
 
+   const BrookStreamImpl& velocityStreamC = dynamic_cast<const BrookStreamImpl&> (velocities.getImpl());
+         BrookStreamImpl& velocityStream  = const_cast<BrookStreamImpl&> (velocityStreamC);
+   void* dataV                            = velocityStream.getData( );
+   float* velocity                        = (float*) dataV;
+
+   double energy                          = 0.0;
+   int index                              = 0;
+
+   if( _masses == NULL ){
+      std::stringstream message;
+      message << methodName << " masses not set.";
+      throw OpenMMException( message.str() );
+   }    
+
 /*
-    RealOpenMM** velData = const_cast<RealOpenMM**>(((BrookFloatStreamImpl&) velocities.getImpl()).getData()); // Brook code needs to be made const correct
-    double energy = 0.0;
-    for (size_t i = 0; i < masses.size(); ++i)
-        energy += masses[i]*(velData[i][0]*velData[i][0]+velData[i][1]*velData[i][1]+velData[i][2]*velData[i][2]);
-    return 0.5*energy;
+printf( "   BrookCalcKineticEnergyKernel Masses=%12.5e %12.5e", _masses[0], _masses[1] );
+printf( " [%12.5e %12.5e %12.5e]", velocity[index], velocity[index+1], velocity[index+2] );
+index += 3;
+printf( " [%12.5e %12.5e %12.5e]\n", velocity[index], velocity[index+1], velocity[index+2] );
+index = 0;
 */
 
-   return 0.0;
+   for ( int ii = 0; ii < velocityStream.getSize(); ii++, index += 3 ){
+      energy += _masses[ii]*(velocity[index]*velocity[index] + velocity[index + 1]*velocity[index + 1] + velocity[index + 2]*velocity[index + 2]);
+   }
+
+   return 0.5*energy;
 }

platforms/brook/src/BrookCalcKineticEnergyKernel.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "kernels.h"
+#include "BrookFloatStreamInternal.h"
 
 namespace OpenMM {
 
@@ -80,6 +81,10 @@ class BrookCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
 
    private:
 
+      // masses
+
+      BrookOpenMMFloat* _masses;
+
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookCalcStandardMMForceFieldKernel.cpp

@@ -38,8 +38,6 @@
 #include "BrookCalcStandardMMForceFieldKernel.h"
 #include "gpu/kforce.h"
 #include "gpu/kinvmap_gather.h"
-#include "ReferencePlatform.h"
-#include "VerletIntegrator.h"
 #include "StandardMMForceField.h"
 
 using namespace OpenMM;
@@ -66,6 +64,13 @@ BrookCalcStandardMMForceFieldKernel::BrookCalcStandardMMForceFieldKernel( std::s
    _brookBonded                             = NULL;
    _brookNonBonded                          = NULL;
    _refForceField                           = NULL;
+   _log                                     = NULL;
+
+   _refForceField                           = NULL;
+   _refSystem                               = NULL;
+   _refOpenMMContext                        = NULL;
+   _referencePlatform                       = NULL;
+   _refVerletIntegrator                     = NULL;
 
    const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
    if( brookPlatform.getLog() != NULL ){
@@ -93,6 +98,12 @@ BrookCalcStandardMMForceFieldKernel::~BrookCalcStandardMMForceFieldKernel( ){
    // deleted w/ kernel delete? If activated, program crashes
 
    //delete _refForceField;
+/*
+   delete _refSystem;
+   delete _refOpenMMContext;
+   delete _referencePlatform;
+   delete _refVerletIntegrator;
+*/
 }
 
 /** 
@@ -319,6 +330,8 @@ void BrookCalcStandardMMForceFieldKernel::executeForces( const Stream& positions
    static const int L_Stream                = 3;
 
    static const int PrintOn                 = 0;
+   static const int MaxErrorMessages        = 2;
+   static       int ErrorMessages           = 0;
 
    static const float4 dummyParameters( 0.0, 0.0, 0.0, 0.0 );
 
@@ -555,7 +568,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
 
       // case not handled -- throw an exception
 
-      if( _brookBonded->getLog() ){
+      if( _brookBonded->getLog() && ErrorMessages++ < MaxErrorMessages ){
          (void) fprintf( _brookBonded->getLog(), "%s case: I-map=%d K-map=%d -- not handled.\n",
                           methodName.c_str(), _brookBonded->getInverseMapStreamCount( I_Stream ),
                                               _brookBonded->getInverseMapStreamCount( K_Stream ) );
@@ -611,7 +624,7 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
 
       // case not handled -- throw an exception
 
-      if( _brookBonded->getLog() ){
+      if( _brookBonded->getLog() && ErrorMessages++ < MaxErrorMessages ){
          (void) fprintf( _brookBonded->getLog(), "%s case: J-map=%d L-map=%d -- not handled.\n",
                           methodName.c_str(), _brookBonded->getInverseMapStreamCount( J_Stream ),
                                               _brookBonded->getInverseMapStreamCount( L_Stream ) );
@@ -671,6 +684,78 @@ nonbondedForceStreams[3]->fillWithValue( &zerof );
 
 double BrookCalcStandardMMForceFieldKernel::executeEnergy( const Stream& atomPositions ){
 
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookCalcStandardMMForceFieldKernel::executeEnergy";
+
+// ---------------------------------------------------------------------------------------
+
+   const BrookStreamImpl& positionStreamC              = dynamic_cast<const BrookStreamImpl&> (atomPositions.getImpl());
+   BrookStreamImpl& positionStream                     = const_cast<BrookStreamImpl&>         (positionStreamC);
+   BrookOpenMMFloat* positionsF                        = (BrookOpenMMFloat*) positionStream.getData();
+
+   OpenMMContext* context                              = getReferenceOpenMMContext( atomPositions.getSize() );
+
+   vector<Vec3> positions( positionStream.getSize() );
+   int index = 0;
+   for( int ii = 0; ii < positionStream.getSize(); ii++, index += 3 ){
+      positions[ii] = Vec3( positionsF[index], positionsF[index+1], positionsF[index+2] );
+   }
+   context->setPositions(positions);
+
+   State state     = context->getState( State::Energy );
+   double energy   = state.getPotentialEnergy();
+
+// (void) fprintf( stdout, "BrookCalcStandardMMForceFieldKernel::executeEnergy E=%.5e\n", energy ); fflush( stdout );
+
+   return energy;
+
+}
+
+/**
+ * Get reference Context
+ * 
+ * @param numberOfAtoms  number of atoms
+ *
+ * @return  OpenMMContext
+ *
+ */
+
+OpenMMContext* BrookCalcStandardMMForceFieldKernel::getReferenceOpenMMContext( int numberOfAtoms ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookCalcStandardMMForceFieldKernel::getReferenceOpenMMContext";
+
+// ---------------------------------------------------------------------------------------
+
+   if( _refOpenMMContext == NULL ){
+      
+      _referencePlatform       = new ReferencePlatform();
+      _refSystem               = new System( numberOfAtoms, 0 ); 
+      _refVerletIntegrator     = new VerletIntegrator( 0.01 );
+
+      _refSystem->addForce( _refForceField );
+
+      _refOpenMMContext        = new OpenMMContext( *_refSystem, *_refVerletIntegrator, *_referencePlatform );
+   }
+
+   return _refOpenMMContext;
+
+}
+
+/**
+ * Execute the kernel to calculate the energy.
+ * 
+ * @param positions   atom positions
+ *
+ * @return  potential energy due to the StandardMMForceField
+ * Currently always return 0.0 since energies not calculated on gpu
+ *
+ */
+
+double BrookCalcStandardMMForceFieldKernel::executeEnergyOld( const Stream& atomPositions ){
+
 // ---------------------------------------------------------------------------------------
 
    //static const std::string methodName      = "BrookCalcStandardMMForceFieldKernel::executeEnergy";
@@ -687,7 +772,7 @@ double BrookCalcStandardMMForceFieldKernel::executeEnergy( const Stream& atomPos
 
    system.addForce( _refForceField );
 
-   OpenMMContext context(system, integrator, refPlatform);
+   OpenMMContext context( system, integrator, refPlatform );
 
    vector<Vec3> positions( positionStream.getSize() );
    int index = 0;

platforms/brook/src/BrookCalcStandardMMForceFieldKernel.h

@@ -37,6 +37,10 @@
 #include "BrookBonded.h"
 #include "BrookNonBonded.h"
 #include "StandardMMForceField.h"
+#include "OpenMMContext.h"
+#include "System.h"
+#include "ReferencePlatform.h"
+#include "VerletIntegrator.h"
 
 namespace OpenMM {
 
@@ -105,7 +109,20 @@ class BrookCalcStandardMMForceFieldKernel : public CalcStandardMMForceFieldKerne
        */
 
       double executeEnergy( const Stream& positions );
+      double executeEnergyOld( const Stream& positions );
 
+      /**
+       * Get reference Context
+       * 
+       * @param numberOfAtoms  number of atoms
+       *
+       * @return  OpenMMContext
+       *
+       */
+      
+      OpenMMContext* getReferenceOpenMMContext( int numberOfAtoms );
+      
+      
       /** 
        * Set log file reference
        * 
@@ -155,6 +172,10 @@ class BrookCalcStandardMMForceFieldKernel : public CalcStandardMMForceFieldKerne
        // used to calculate energy
 
        StandardMMForceField* _refForceField;
+       System*               _refSystem;
+       OpenMMContext*        _refOpenMMContext;
+       ReferencePlatform*    _referencePlatform;
+       VerletIntegrator*     _refVerletIntegrator;
 
 };
 

platforms/brook/src/BrookIntegrateLangevinStepKernel.cpp

@@ -35,6 +35,14 @@
 using namespace OpenMM;
 using namespace std;
 
+      /** 
+       * BrookIntegrateLangevinStepKernel constructor
+       * 
+       * @param name        name of the stream to create
+       * @param platform    platform
+       *
+       */
+
 BrookIntegrateLangevinStepKernel::BrookIntegrateLangevinStepKernel( std::string name, const Platform& platform ) :
                                   IntegrateLangevinStepKernel( name, platform ){
 
@@ -50,6 +58,11 @@ BrookIntegrateLangevinStepKernel::BrookIntegrateLangevinStepKernel( std::string
 
 }
 
+/** 
+ * BrookIntegrateVerletStepKernel destructor
+ * 
+ */
+  
 BrookIntegrateLangevinStepKernel::~BrookIntegrateLangevinStepKernel( ){
 
 // ---------------------------------------------------------------------------------------
@@ -64,6 +77,15 @@ BrookIntegrateLangevinStepKernel::~BrookIntegrateLangevinStepKernel( ){
 
 }
 
+/** 
+ * Initialize the kernel, setting up all parameters related to integrator.
+ * 
+ * @param masses             the mass of each atom
+ * @param constraintIndices  each element contains the indices of two atoms whose distance should be constrained
+ * @param constraintLengths  the required distance between each pair of constrained atoms
+ *
+ */
+
 void BrookIntegrateLangevinStepKernel::initialize( const vector<double>& masses,
                                                    const vector<vector<int> >& constraintIndices,
                                                    const vector<double>& constraintLengths ){
@@ -121,7 +143,7 @@ void BrookIntegrateLangevinStepKernel::execute( Stream& positions, Stream& veloc
    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 ){
-//printf( "%s calling updateParameters\n", methodName.c_str() );
+printf( "%s calling updateParameters\n", methodName.c_str() );
       _brookStochasticDynamics->updateParameters( temperature, friction, stepSize );
    } else {
 //printf( "%s NOT calling updateParameters\n", methodName.c_str() );

platforms/brook/src/BrookIntegrateVerletStepKernel.cpp

@@ -35,6 +35,14 @@
 using namespace OpenMM;
 using namespace std;
 
+/** 
+ * BrookIntegrateVerletStepKernel constructor
+ * 
+ * @param name        name of the stream to create
+ * @param platform    platform
+ *
+ */
+
 BrookIntegrateVerletStepKernel::BrookIntegrateVerletStepKernel( std::string name, const Platform& platform ) :
                                 IntegrateVerletStepKernel( name, platform ){
 
@@ -44,8 +52,15 @@ BrookIntegrateVerletStepKernel::BrookIntegrateVerletStepKernel( std::string name
 
 // ---------------------------------------------------------------------------------------
    
+   _brookVerletDynamics   = NULL;
+   _brookShakeAlgorithm   = NULL;
 }
 
+/** 
+ * BrookIntegrateVerletStepKernel destructor
+ * 
+ */
+  
 BrookIntegrateVerletStepKernel::~BrookIntegrateVerletStepKernel( ){
 
 // ---------------------------------------------------------------------------------------
@@ -53,21 +68,21 @@ BrookIntegrateVerletStepKernel::~BrookIntegrateVerletStepKernel( ){
    // static const std::string methodName      = "BrookIntegrateVerletStepKernel::~BrookIntegrateVerletStepKernel";
 
 // ---------------------------------------------------------------------------------------
+
+   delete _brookVerletDynamics;
+   delete _brookShakeAlgorithm;
    
-/*
-    if (dynamics)
-        delete dynamics;
-    if (shake)
-        delete shake;
-    if (masses)
-        delete[] masses;
-    if (constraintIndices)
-        disposeIntArray(constraintIndices, numConstraints);
-    if (shakeParameters)
-        disposeRealArray(shakeParameters, numConstraints);
-*/
 }
 
+/** 
+ * Initialize the kernel, setting up all parameters related to integrator.
+ * 
+ * @param masses             the mass of each atom
+ * @param constraintIndices  each element contains the indices of two atoms whose distance should be constrained
+ * @param constraintLengths  the required distance between each pair of constrained atoms
+ *
+ */
+
 void BrookIntegrateVerletStepKernel::initialize( const vector<double>& masses,
                                                  const vector<vector<int> >& constraintIndices,
                                                  const vector<double>& constraintLengths ){
@@ -78,47 +93,47 @@ void BrookIntegrateVerletStepKernel::initialize( const vector<double>& masses,
 
 // ---------------------------------------------------------------------------------------
    
-/*
-    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] );
-*/
+   _brookVerletDynamics          = new BrookVerletDynamics( );
+   _brookVerletDynamics->setup( masses, getPlatform() );
+
+   _brookShakeAlgorithm          = new BrookShakeAlgorithm( );
+   _brookShakeAlgorithm->setup( masses, constraintIndices, constraintLengths, getPlatform() );
+
 }
 
+/** 
+ * Execute kernel
+ * 
+ * @param positions          atom coordinates
+ * @param velocities         atom velocities
+ * @param forces             atom forces
+ * @param stepSize           integration step size
+ *
+ */
+
 void BrookIntegrateVerletStepKernel::execute( Stream& positions, Stream& velocities,
                                               const Stream& forces, double stepSize ){
+
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookIntegrateVerletStepKernel::execute";
+   double epsilon                           = 1.0e-04;
+   static const std::string methodName      = "BrookIntegrateVerletStepKernel::execute";
 
 // ---------------------------------------------------------------------------------------
-   
-/*
-    RealOpenMM** posData = ((BrookFloatStreamImpl&) positions.getImpl()).getData();
-    RealOpenMM** velData = ((BrookFloatStreamImpl&) velocities.getImpl()).getData();
-    RealOpenMM** forceData = const_cast<RealOpenMM**>(((BrookFloatStreamImpl&) forces.getImpl()).getData()); // Brook code needs to be made const correct
-    if (dynamics == 0 || stepSize != prevStepSize) {
-        // Recreate the computation objects with the new parameters.
-        
-        if (dynamics) {
-            delete dynamics;
-            delete shake;
-        }
-        dynamics = new BrookVerletDynamics(positions.getSize(), static_cast<RealOpenMM>(stepSize) );
-        shake = new BrookShakeAlgorithm(numConstraints, constraintIndices, shakeParameters);
-        dynamics->setBrookShakeAlgorithm(shake);
-        prevStepSize = stepSize;
-    }
-    dynamics->update(positions.getSize(), posData, velData, forceData, masses);
-*/
+
+   // first time through initialize _brookVerletDynamics
+
+   // for each subsequent call, check if parameters need to be updated due to a change
+   // in the step size
+
+   // take step
+
+   double difference = stepSize - (double) _brookVerletDynamics->getStepSize();
+   if( fabs( difference ) > epsilon ){
+//printf( "%s calling updateParameters\n", methodName.c_str() );
+      _brookVerletDynamics->updateParameters( stepSize );
+   }
+   _brookVerletDynamics->update( positions, velocities, forces, *_brookShakeAlgorithm );
+
 }
 

platforms/brook/src/BrookIntegrateVerletStepKernel.h

-#ifndef OPENMM_BROOK_INTEGRATE_KERNELS_H_
-#define OPENMM_BROOK_INTEGRATE_KERNELS_H_
+#ifndef OPENMM_BROOK_INTEGRATE_VERLET_STEP_KERNEL_H_
+#define OPENMM_BROOK_INTEGRATE_VERLET_STEP_KERNEL_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -33,6 +33,8 @@
  * -------------------------------------------------------------------------- */
 
 #include "kernels.h"
+#include "BrookVerletDynamics.h"
+#include "BrookShakeAlgorithm.h"
 
 namespace OpenMM {
 
@@ -72,23 +74,24 @@ class BrookIntegrateVerletStepKernel : public IntegrateVerletStepKernel {
       void initialize( const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
                        const std::vector<double>& constraintLengths );
       /** 
-       * Execute the kernel.
+       * Execute kernel
        * 
-       * @param positions          a Stream of type Double3 containing the position (x, y, z) of each atom
-       * @param velocities         a Stream of type Double3 containing the velocity (x, y, z) of each atom
-       * @param forces             a Stream of type Double3 containing the force (x, y, z) on each atom
-       * @param stepSize           the integration step size
+       * @param positions          atom coordinates
+       * @param velocities         atom velocities
+       * @param forces             atom forces
+       * @param stepSize           integration step size
        *
        */
-
+      
       void execute( Stream& positions, Stream& velocities, const Stream& forces, double stepSize );
 
    protected:
 
-  
-  
+      BrookVerletDynamics*      _brookVerletDynamics;
+      BrookShakeAlgorithm*      _brookShakeAlgorithm;
+ 
 };
 
 } // namespace OpenMM
 
-#endif /* OPENMM_BROOK_INTEGRATE_KERNELS_H_ */
+#endif /* OPENMM_BROOK_INTEGRATE_VERLET_STEP_KERNEL_H_ */

platforms/brook/src/BrookKernelFactory.cpp

@@ -32,6 +32,7 @@
 #include "BrookKernelFactory.h"
 #include "BrookCalcStandardMMForceFieldKernel.h"
 #include "BrookIntegrateLangevinStepKernel.h"
+#include "BrookIntegrateVerletStepKernel.h"
 #include "BrookCalcKineticEnergyKernel.h"
 #include "BrookCalcGBSAOBCForceFieldKernel.h"
 
@@ -61,7 +62,7 @@ KernelImpl* BrookKernelFactory::createKernelImpl( std::string name, const Platfo
 
 	} else if( name == IntegrateVerletStepKernel::Name() ){
 
-      // return new BrookIntegrateVerletStepKernel( name, platform );
+      return new BrookIntegrateVerletStepKernel( name, platform );
 
    // Brownian integrator
 
@@ -69,7 +70,7 @@ KernelImpl* BrookKernelFactory::createKernelImpl( std::string name, const Platfo
 
       // return new BrookIntegrateBrownianStepKernel( name, platform );
 
-   //  Andersen integrator
+   //  Andersen thermostat
 
 	} else if( name == ApplyAndersenThermostatKernel::Name() ){
 

platforms/brook/src/BrookPlatform.cpp

@@ -155,7 +155,7 @@ void BrookPlatform::_initializeKernelFactory( void ){
 
    registerKernelFactory( CalcStandardMMForceFieldKernel::Name(), factory);
    registerKernelFactory( CalcGBSAOBCForceFieldKernel::Name(),    factory);
-   //registerKernelFactory( IntegrateVerletStepKernel::Name(),      factory);
+   registerKernelFactory( IntegrateVerletStepKernel::Name(),      factory);
    registerKernelFactory( IntegrateLangevinStepKernel::Name(),    factory);
    //registerKernelFactory( IntegrateBrownianStepKernel::Name(),    factory);
    //registerKernelFactory( ApplyAndersenThermostatKernel::Name(),  factory);

platforms/brook/src/BrookStochasticDynamics.cpp

@@ -375,7 +375,7 @@ int BrookStochasticDynamics::update( Stream& positions, Stream& velocities,
 
    static const char* methodName = "\nBrookStochasticDynamics::update";
 
-   static const int PrintOn      = 0;
+   static const int PrintOn      = 1;
 
 // ---------------------------------------------------------------------------------------
 

platforms/brook/src/BrookVerletDynamics.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 "BrookVerletDynamics.h"
+#include "BrookPlatform.h"
+#include "OpenMMException.h"
+#include "BrookStreamImpl.h"
+#include "gpu/kshakeh.h"
+#include "gpu/kupdatemd.h"
+#include "gpu/kcommon.h"
+
+// use random number generator
+
+#include "../../reference/src/SimTKUtilities/SimTKOpenMMUtilities.h"
+
+using namespace OpenMM;
+using namespace std;
+
+/** 
+ *
+ * Constructor
+ * 
+ */
+
+BrookVerletDynamics::BrookVerletDynamics( ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVerletDynamics::BrookVerletDynamics";
+
+   BrookOpenMMFloat zero                    = (BrookOpenMMFloat)  0.0;
+   BrookOpenMMFloat one                     = (BrookOpenMMFloat)  1.0;
+   BrookOpenMMFloat oneMinus                = (BrookOpenMMFloat) -1.0;
+
+// ---------------------------------------------------------------------------------------
+
+   _numberOfAtoms                 = -1;
+
+   // mark stream dimension variables as unset
+
+   _verletAtomStreamWidth         = -1;
+   _verletAtomStreamHeight        = -1;
+   _verletAtomStreamSize          = -1;
+
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      _verletStreams[ii]   = NULL;
+   }
+
+   _stepSize                      = oneMinus;
+
+   // setup inverse sqrt masses
+
+   _inverseMasses             = NULL;
+
+}   
+ 
+/** 
+ * Destructor
+ * 
+ */
+
+BrookVerletDynamics::~BrookVerletDynamics( ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVerletDynamics::~BrookVerletDynamics";
+
+// ---------------------------------------------------------------------------------------
+
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      delete _verletStreams[ii];
+   }
+
+   delete[] _inverseMasses;
+
+}
+
+/** 
+ * Get stepSize
+ * 
+ * @return  stepSize
+ *
+ */
+
+BrookOpenMMFloat BrookVerletDynamics::getStepSize( void ) const {
+   return _stepSize;
+}
+
+/** 
+ * Set stepSize
+ * 
+ * @param   stepSize
+ *
+ * @return      DefaultReturnValue
+ *
+ */
+
+int BrookVerletDynamics::_setStepSize( BrookOpenMMFloat stepSize ){
+   _stepSize = stepSize;
+   return DefaultReturnValue;
+}
+
+/** 
+ * Update parameters -- only way parameters can be set
+ * 
+ * @param  step size       step size
+ *
+ * @return   DefaultReturnValue
+ *
+ */
+
+int BrookVerletDynamics::updateParameters( double stepSize ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static int showUpdate         = 1;
+   static int maxShowUpdate      = 3;
+   static const char* methodName = "\nBrookVerletDynamics::updateParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   _setStepSize(    (BrookOpenMMFloat)  stepSize );
+   _updateVerletStreams( );
+
+   // show update
+
+   if( showUpdate && getLog() && (showUpdate++ < maxShowUpdate) ){
+      std::string contents = getContentsString( );
+      (void) fprintf( getLog(), "%s contents %s", methodName, contents.c_str() );
+      (void) fflush( getLog() );
+
+   }
+
+   return DefaultReturnValue;
+
+};
+
+/** 
+ * Update
+ * 
+ * @param  positions                   atom positions
+ * @param  velocities                  atom velocities
+ * @param  forces                      atom forces
+ * @param  brookShakeAlgorithm         BrookShakeAlgorithm reference
+ * @param  brookRandomNumberGenerator  BrookRandomNumberGenerator reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int BrookVerletDynamics::update( Stream& positions, Stream& velocities,
+                                 const Stream& forces,
+                                 BrookShakeAlgorithm& brookShakeAlgorithm ){
+
+// ---------------------------------------------------------------------------------------
+
+   // unused Shake parameter
+
+   float omega                   = 1.0f;
+
+   static const char* methodName = "\nBrookVerletDynamics::update";
+
+   static const int PrintOn      = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   BrookStreamImpl& positionStream                     = dynamic_cast<BrookStreamImpl&>       (positions.getImpl());
+   BrookStreamImpl& velocityStream                     = dynamic_cast<BrookStreamImpl&>       (velocities.getImpl());
+   const BrookStreamImpl& forceStreamC                 = dynamic_cast<const BrookStreamImpl&> (forces.getImpl());
+   BrookStreamImpl& forceStream                        = const_cast<BrookStreamImpl&> (forceStreamC);
+
+   if( PrintOn && getLog() ){
+      (void) fprintf( getLog(), "%s shake=%d\n", methodName, brookShakeAlgorithm.getNumberOfConstraints() );
+      (void) fflush( getLog() );
+
+   }
+
+   // integration step
+
+   kupdate_md_verlet( (float) getStepSize(),
+                       positionStream.getBrookStream(),
+                       velocityStream.getBrookStream(),
+                       forceStream.getBrookStream(),
+                       getInverseMassStream()->getBrookStream(),
+                       velocityStream.getBrookStream(),
+                       getXPrimeStream()->getBrookStream() 
+                    );
+ 
+   // diagnostics
+
+   if( PrintOn && getLog() ){
+      (void) fprintf( getLog(), "\nPost kupdate_md_verlet: atomStrW=%3d step=%.5f",
+                                getVerletDynamicsAtomStreamWidth(), getStepSize() );
+
+      (void) fprintf( getLog(), "\nInverseMassStream\n" );
+      getInverseMassStream()->printToFile( getLog() );
+
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+      (void) fprintf( getLog(), "\nPositionStream\n" );
+      brookStreamInternalPos->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nForceStream\n" );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      brookStreamInternalF->printToFile( getLog() );
+
+      BrookStreamInternal* brookStreamInternalV   = velocityStream.getBrookStreamImpl();
+      (void) fprintf( getLog(), "\nVelocityStream\n" );
+      brookStreamInternalV->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nXPrimeStream\n" );
+      getXPrimeStream()->printToFile( getLog() ); 
+   }   
+
+   // second Shake step
+
+   if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
+      kshakeh_fix1( 
+                    10.0f,
+                    (float) getVerletDynamicsAtomStreamWidth(),
+                    brookShakeAlgorithm.getInverseHydrogenMass(),
+                    omega, 
+                    brookShakeAlgorithm.getShakeAtomIndicesStream()->getBrookStream(),
+                    positionStream.getBrookStream(),
+                    getXPrimeStream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeAtomParameterStream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream() );
+   
+      // second Shake gather
+   
+      kshakeh_update2_fix1( 
+                    (float) getVerletDynamicsAtomStreamWidth(),
+                    brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
+                    positionStream.getBrookStream(),
+                    getXPrimeStream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeXCons3Stream()->getBrookStream(),
+                    positionStream.getBrookStream() );
+   
+   } else {
+      //kadd3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
+      ksetStr3( getXPrimeStream()->getBrookStream(), positionStream.getBrookStream() );
+   }
+
+   return DefaultReturnValue;
+
+};
+
+/** 
+ * Get Atom stream size
+ *
+ * @return  Atom stream size
+ *
+ */
+
+int BrookVerletDynamics::getVerletDynamicsAtomStreamSize( void ) const {
+   return _verletAtomStreamSize;
+}
+
+/** 
+ * Get atom stream width
+ *
+ * @return  atom stream width
+ *
+ */
+
+int BrookVerletDynamics::getVerletDynamicsAtomStreamWidth( void ) const {
+   return _verletAtomStreamWidth;
+}
+
+/** 
+ * Get atom stream height
+ *
+ * @return atom stream height
+ */
+
+int BrookVerletDynamics::getVerletDynamicsAtomStreamHeight( void ) const {
+   return _verletAtomStreamHeight;
+}
+
+/** 
+ * Get VPrime stream 
+ *
+ * @return  Vprime stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookVerletDynamics::getVPrimeStream( void ) const {
+   return _verletStreams[VPrimeStream];
+}
+
+/** 
+ * Get XPrime stream 
+ *
+ * @return  Xprime stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookVerletDynamics::getXPrimeStream( void ) const {
+   return _verletStreams[XPrimeStream];
+}
+
+/** 
+ * Get InverseMass stream 
+ *
+ * @return  inverse mass stream
+ *
+ */
+
+BrookFloatStreamInternal* BrookVerletDynamics::getInverseMassStream( void ) const {
+   return _verletStreams[InverseMassStream];
+}
+
+/** 
+ * Initialize stream dimensions
+ * 
+ * @param numberOfAtoms             number of atoms
+ * @param platform                  platform
+ *      
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookVerletDynamics::_initializeStreamSizes( int numberOfAtoms, const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVerletDynamics::_initializeStreamSizes";
+
+// ---------------------------------------------------------------------------------------
+
+   _verletAtomStreamSize     = getAtomStreamSize( platform );
+   _verletAtomStreamWidth    = getAtomStreamWidth( platform );
+   _verletAtomStreamHeight   = getAtomStreamHeight( platform );
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Initialize streams
+ * 
+ * @param platform                  platform
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookVerletDynamics::_initializeStreams( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVerletDynamics::_initializeStreams";
+
+   BrookOpenMMFloat dangleValue            = (BrookOpenMMFloat) 0.0;
+
+// ---------------------------------------------------------------------------------------
+
+   int sdAtomStreamSize             = getVerletDynamicsAtomStreamSize();
+   int sdAtomStreamWidth            = getVerletDynamicsAtomStreamWidth();
+
+    _verletStreams[VPrimeStream]        = new BrookFloatStreamInternal( BrookCommon::VPrimeStream,
+                                                                    sdAtomStreamSize, sdAtomStreamWidth,
+                                                                    BrookStreamInternal::Float3, dangleValue );
+
+    _verletStreams[XPrimeStream]        = new BrookFloatStreamInternal( BrookCommon::XPrimeStream,
+                                                                    sdAtomStreamSize, sdAtomStreamWidth,
+                                                                    BrookStreamInternal::Float3, dangleValue );
+
+    _verletStreams[InverseMassStream]   = new BrookFloatStreamInternal( BrookCommon::InverseMassStream,
+                                                                    sdAtomStreamSize, sdAtomStreamWidth,
+                                                                    BrookStreamInternal::Float, dangleValue );
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Update sd streams -- called after parameters change
+ * 
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookVerletDynamics::_updateVerletStreams( void ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName       = "BrookVerletDynamics::_updateVerletStreams";
+
+// ---------------------------------------------------------------------------------------
+
+   int atomStreamSize                      = getVerletDynamicsAtomStreamSize();
+
+   BrookOpenMMFloat* inverseMass = new BrookOpenMMFloat[atomStreamSize];
+   memset( inverseMass, 0, atomStreamSize*sizeof( BrookOpenMMFloat ) ); 
+
+   int numberOfAtoms                         = getNumberOfAtoms();
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+      inverseMass[ii]  = _inverseMasses[ii];
+   }
+
+   _verletStreams[InverseMassStream]->loadFromArray( inverseMass );
+
+   delete[] inverseMass;
+
+   return DefaultReturnValue;
+
+}
+
+/** 
+ * Set masses 
+ * 
+ * @param masses             atomic masses
+ *
+ */
+
+int BrookVerletDynamics::_setInverseMasses( const std::vector<double>& masses ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookVerletDynamics::_setInverseSqrtMasses";
+
+   BrookOpenMMFloat zero                    = (BrookOpenMMFloat)  0.0;
+   BrookOpenMMFloat one                     = (BrookOpenMMFloat)  1.0;
+
+// ---------------------------------------------------------------------------------------
+
+   // setup inverse sqrt masses
+
+   _inverseMasses = new BrookOpenMMFloat[masses.size()];
+   int index          = 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] = ( one/value );
+      } else {
+         _inverseMasses[index] = zero;
+      }
+   }
+
+   return DefaultReturnValue;
+}   
+ 
+/*  
+ * Setup of VerletDynamics parameters
+ *
+ * @param masses                masses
+ * @param platform              Brook platform
+ *
+ * @return nonzero value if error
+ *
+ * */
+    
+int BrookVerletDynamics::setup( const std::vector<double>& masses, const Platform& platform ){
+    
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookVerletDynamics::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 );
+
+   _setInverseMasses( masses );
+
+   return DefaultReturnValue;
+}
+
+/* 
+ * Get contents of object
+ *
+ * @param level   level of dump
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string BrookVerletDynamics::getContentsString( int level ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookVerletDynamics::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, "%.5f", getStepSize() );
+   message << _getLine( tab, "Step size:", value ); 
+
+   message << _getLine( tab, "Log:",                  (getLog()                ? Set : NotSet) ); 
+
+   for( int ii = 0; ii < LastStreamIndex; ii++ ){
+      message << std::endl;
+      if( _verletStreams[ii] ){
+         message << _verletStreams[ii]->getContentsString( );
+      }
+   }
+
+#undef LOCAL_SPRINTF
+
+   return message.str();
+}

platforms/brook/src/BrookVerletDynamics.h

+#ifndef OPENMM_BROOK_VERLET_DYNAMCIS_H_
+#define OPENMM_BROOK_VERLET_DYNAMCIS_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 "BrookShakeAlgorithm.h"
+#include "BrookPlatform.h"
+#include "BrookCommon.h"
+
+namespace OpenMM {
+
+/**
+ *
+ * Encapsulates stochastic dynamics algorithm 
+ *
+ */
+
+class BrookVerletDynamics : public BrookCommon {
+
+   public:
+  
+      /** 
+       * Constructor
+       * 
+       */
+      
+      BrookVerletDynamics(  );
+  
+      /** 
+       * Destructor
+       * 
+       */
+      
+      ~BrookVerletDynamics();
+  
+      /**
+       * Get step size
+       *
+       * @return step size
+       */
+
+      BrookOpenMMFloat getStepSize( void ) const; 
+
+      /**
+       *
+       * Get array of derived parameters indexed by 'DerivedParameters' enums
+       *
+       * @return array
+       *
+       */
+      
+      const BrookOpenMMFloat* getDerivedParameters( void ) const;
+      
+      /**
+       * Get VerletDynamics atom stream width
+       *
+       * @return atom stream width
+       */
+
+      int getVerletDynamicsAtomStreamWidth( void ) const; 
+
+      /**
+       * Get VerletDynamics atom stream height
+       *
+       * @return atom stream height
+       */
+
+      int getVerletDynamicsAtomStreamHeight( void ) const;
+
+      /**
+       * Get VerletDynamics atom stream size
+       * 
+       * @return atom stream size
+       */
+
+      int getVerletDynamicsAtomStreamSize( void ) const; 
+
+      /** 
+       * Update parameters
+       * 
+       * @param  step size       step size
+       *
+       * @return   DefaultReturnValue
+       *
+       */
+      
+      int updateParameters( 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 VerletDynamics streams 
+       *
+       * @return  array ofstreams
+       *
+       */
+      
+      BrookFloatStreamInternal** getStreams( void );
+      
+      /* 
+       * Setup of VerletDynamics 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 V-prime stream 
+       *
+       * @return  V-prime stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getVPrimeStream( void ) const;
+
+      /** 
+       * Get X-prime stream 
+       *
+       * @return  X-prime stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getXPrimeStream( void ) const;
+
+      /** 
+       * Get inverse sqrt masses 
+       * 
+       * @return   inverse sqrt masses stream
+       *
+       */
+      
+      BrookFloatStreamInternal* getInverseMassStream( void ) const;
+
+   private:
+   
+      // streams indices
+
+      enum BrookVerletDynamicsStreams { 
+              VPrimeStream,
+              XPrimeStream,
+              InverseMassStream,
+              LastStreamIndex
+           };
+
+      BrookOpenMMFloat _stepSize;
+
+      // Atom stream dimensions
+
+      int _verletAtomStreamWidth;
+      int _verletAtomStreamHeight;
+      int _verletAtomStreamSize;
+
+      /*
+       * Update streams
+       * 
+       * @return  DefaultReturn
+       *
+       */
+      
+      int _updateVerletStreams( void );
+      
+      // inverse masses
+
+      BrookOpenMMFloat* _inverseMasses;
+
+      // internal streams
+
+      BrookFloatStreamInternal* _verletStreams[LastStreamIndex];
+
+      /** 
+       * Set stepSize
+       * 
+       * @param   stepSize
+       *
+       * @return      DefaultReturn
+       *
+       */
+      
+      int _setStepSize( BrookOpenMMFloat stepSize );
+      
+      /* 
+       * 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 masses 
+       * 
+       * @param masses             atomic masses
+       *
+       */
+      
+      int _setInverseMasses( const std::vector<double>& masses );
+      
+};
+
+} // namespace OpenMM
+
+#endif /* OPENMM_BROOK_VERLET_DYNAMCIS_H_ */

platforms/brook/src/gpu/kupdatemd.br

@@ -5,6 +5,19 @@
 * Copyright (C) Pande Group, Stanford, 2006
 *****************************************************************/
 
+kernel void kupdate_md_verlet( 
+		float dt,
+		float3 posq<>, 
+		float3 v<>, 
+		float3 f<>,
+		float invmass<>,
+		out float3 vnew<>, 
+		out float3 posqp<> ){
+
+	vnew   = v + f*invmass*dt;
+	posqp  = posq + vnew*dt;
+}
+
 kernel void kupdate_md1_berendsen( 
 		float dt,
 		float3 lg, //Berendsen coupling, assuming only one group

platforms/brook/src/gpu/kupdatemd.h

+void  kupdate_md2 (const float  dtinv,
+		::brook::stream posqp,
+		::brook::stream posq,
+		::brook::stream vnew,
+		::brook::stream posqnew
+	    ); 
+
+void  kupdate_md1_extended (
+		const float  dt,
+		const float3  lg,
+		const float  xi,
+		const float3  M0,
+		const float3  M1,
+		const float3  M2,
+		const float3  uold,
+		::brook::stream posq,
+		::brook::stream v,
+		::brook::stream f,
+		::brook::stream invmass,
+		::brook::stream vnew,
+		::brook::stream posqp);
+
+void  kupdate_md_verlet(
+		const float  dt,
+		::brook::stream posq,
+		::brook::stream v,
+		::brook::stream f,
+		::brook::stream invmass,
+		::brook::stream vnew,
+		::brook::stream posqp);
+
+
+void  kupdate_md1_berendsen (
+		const float  dt,
+		const float3  lg,
+		const float3  uold,
+		::brook::stream posq,
+		::brook::stream v,
+		::brook::stream f,
+		::brook::stream invmass,
+		::brook::stream vnew,
+		::brook::stream posqp); 
+
+void  kupdate_md2_fix1 (const float  dtinv,
+		::brook::stream posqp,
+		::brook::stream posq,
+		::brook::stream vnew,
+		::brook::stream posqnew); 
+
+void  kupdate_md1_extended_fix1 (const float  dt,
+		const float3  lg,
+		const float  xi,
+		const float3  M0,
+		const float3  M1,
+		const float3  M2,
+		const float3  uold,
+		::brook::stream posq,
+		::brook::stream v,
+		::brook::stream f,
+		::brook::stream invmass,
+		::brook::stream vnew,
+		::brook::stream posqp);
+
+void  kupdate_md1_berendsen_fix1 (const float  dt,
+		const float3  lg,
+		const float3  uold,
+		::brook::stream posq,
+		::brook::stream v,
+		::brook::stream f,
+		::brook::stream invmass,
+		::brook::stream vnew,
+		::brook::stream posqp);