Mods

platforms/brook/src/BrookBondParameters.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: 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 <sstream>
+#include "BrookBondParameters.h"
+#include "OpenMMException.h"
+
+using namespace OpenMM;
+using namespace std;
+
+/** 
+ * BrookBondParameters constructor
+ * 
+ * @param numberOfParticlesInBond       no. of particles in each bond
+ * @param numberOfParametersInBond  no. of parameters in each bond
+ * @param numberOfBonds             no. of bonds
+ * @param log                       optional log reference
+ *
+ */
+
+BrookBondParameters::BrookBondParameters( std::string bondName, int numberOfParticlesInBond,
+                                          int numberOfParametersInBond, int numberOfBonds, FILE* log = NULL ) :
+ _bondName( bondName ), _numberOfParticlesInBond( numberOfParticlesInBond ),
+ _numberOfParametersInBond( numberOfParametersInBond ), _numberOfBonds( numberOfBonds ), _log( log ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookBondParameters::BrookBondParameters";
+
+// ---------------------------------------------------------------------------------------
+
+   // allocate memory for particle indices/parameters
+
+   int* particleIndicesBlock                    = new int[numberOfParticlesInBond*numberOfBonds];
+   double* parametersBlock                  = new double[numberOfParametersInBond*numberOfBonds];
+   _particleIndices                             = new int*[numberOfBonds];
+   _bondParameters                          = new double*[numberOfBonds];
+
+   for( int ii = 0; ii < numberOfBonds; ii++ ){
+
+      _particleIndices[ii]                      = particleIndicesBlock; 
+      particleIndicesBlock                     += numberOfParticlesInBond;
+
+      _bondParameters[ii]                   = parametersBlock;
+      parametersBlock                      += numberOfParametersInBond;
+   }
+
+}   
+
+/** 
+ * BrookBondParameters destructor
+ * 
+ */
+
+BrookBondParameters::~BrookBondParameters( ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookBondParameters::BrookBondParameters";
+
+// ---------------------------------------------------------------------------------------
+
+   delete[] _particleIndices[0];
+   delete[] _particleIndices;
+
+   delete[] _bondParameters[0];
+   delete[] _bondParameters;
+
+}
+
+/** 
+ * Get log file reference
+ * 
+ * @return  log file reference
+ *
+ */
+
+FILE* BrookBondParameters::getLog( void ) const {
+   return _log;
+}
+
+/** 
+ * Set log file reference
+ * 
+ * @param  log file reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int BrookBondParameters::setLog( FILE* log ){
+   _log = log;
+   return DefaultReturnValue;
+}
+
+/** 
+ * Get number of bonds
+ * 
+ * @return  numberOfBonds
+ *
+ */
+
+int BrookBondParameters::getNumberOfBonds( void ) const {
+   return _numberOfBonds;
+}
+
+/** 
+ * Get number of particles in bond
+ * 
+ * @return  numberOfParticlesInBond
+ *
+ */
+
+int BrookBondParameters::getNumberOfParticlesInBond( void ) const {
+   return _numberOfParticlesInBond;
+}
+
+/** 
+ * Get number of parameters in bond
+ * 
+ * @return numberOfParametersInBond 
+ *
+ */
+
+int BrookBondParameters::getNumberOfParametersInBond( void ) const {
+   return _numberOfParametersInBond;
+}
+
+/** 
+ * Initialize the kernel, setting up the values of all the force field parameters.
+ * 
+ * @param bondIndex                 bond index to be set
+ * @param particleIndices           indices of particles (dimension=numberOfParticlesInBond)
+ * @param bondParameters            bond parameters (dimension=numberOfParametersInBond)
+ *
+ * @return DefaultReturnValue
+ *
+ * @throw OpenMMException exeception if bond index is invalid
+ *
+ */
+ 
+int BrookBondParameters::setBond( int bondIndex, int* particleIndices, double* bondParameters ){ 
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookBondParameters::setBond";
+
+// ---------------------------------------------------------------------------------------
+
+   FILE* log                 = getLog();
+
+   // ---------------------------------------------------------------------------------------
+
+   // validate bond index
+
+   if( bondIndex < 0 ){
+      std::stringstream message;
+      message << methodName << "BondIndex=" << bondIndex << " is < 0.";
+      throw OpenMMException( message.str() );
+   }
+
+   if( bondIndex >= getNumberOfBonds() ){
+      std::stringstream message;
+      message << methodName << "BondIndex=" << bondIndex << " is >= " << getNumberOfBonds() << ".";
+      throw OpenMMException( message.str() );
+   }
+
+   // load'em up
+
+   int numberOfParticlesInBond = getNumberOfParticlesInBond();
+   for( int ii = 0; ii < numberOfParticlesInBond; ii++ ){
+      _particleIndices[bondIndex][ii] = particleIndices[ii];
+   }
+
+   int numberOfParametersInBond = getNumberOfParametersInBond();
+   for( int ii = 0; ii < numberOfParametersInBond; ii++ ){
+      _bondParameters[bondIndex][ii] = bondParameters[ii];
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   return DefaultReturnValue;
+}
+
+/* 
+ * Get contents of object
+ *
+ * @param tab         tab
+ * @param description description
+ * @param value       value
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string BrookBondParameters::_getLine( const std::string& tab,
+                                           const std::string& description,
+                                           const std::string& value ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookStreamInternal::_getLine";
+
+   static const unsigned int MAX_LINE_CHARS = 256;
+   char line[MAX_LINE_CHARS];
+
+// ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   memset( line, ' ', MAX_LINE_CHARS );  
+#ifdef WIN32
+   (void) sprintf_s( line, MAX_LINE_CHARS, "%s %-40s %s", tab.c_str(), description.c_str(), value.c_str() );
+#else
+   (void) sprintf( line, "%s %-40s %s", tab.c_str(), description.c_str(), value.c_str() );
+#endif
+   message << std::string( line ) << std::endl;
+
+   return message.str();
+
+}
+
+/* 
+ * Get contents of object
+ *
+ *
+ * @param level   level of dump
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string BrookBondParameters::getContentsString( int level ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookBondParameters::getContentsString";
+
+   static const unsigned int MAX_LINE_CHARS = 1024;
+   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) );   
+#define LOCAL_2_SPRINTF(a,b,c,d) sprintf_s( (a), MAX_LINE_CHARS, (b), (c), (d) );   
+#else
+#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
+#define LOCAL_2_SPRINTF(a,b,c,d) sprintf( (a), (b), (c), (d) );   
+#endif
+
+   (void) LOCAL_SPRINTF( value, "%s", getBondName() );
+   message << _getLine( tab, "Bond name:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfBonds() );
+   message << _getLine( tab, "Number of bonds:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfParticlesInBond() );
+   message << _getLine( tab, "Particles/bond:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfParametersInBond() );
+   message << _getLine( tab, "Parameters/bond:", value ); 
+
+   message << "Bonds:"  << std::endl; 
+   for( int ii = 0; ii < getNumberOfBonds(); ii++ ){
+      char description[256];
+      char buffer[256];
+      (void) LOCAL_SPRINTF( description, "%6d [", ii );
+
+      // particle indices
+
+      for( int jj = 0; jj < getNumberOfParticlesInBond(); jj++ ){
+         (void) LOCAL_SPRINTF( buffer, "%6d ", _particleIndices[ii][jj] );
+         (void) strcat( description, buffer );
+      }
+      (void) strcat( description, "] [" );
+
+      // parameters
+
+      for( int jj = 0; jj < getNumberOfParametersInBond(); jj++ ){
+         (void) LOCAL_SPRINTF( buffer, "%18.10e ", _bondParameters[ii][jj] );
+         (void) strcat( description, buffer );
+      }
+      (void) strcat( description, "]" );
+      message << _getLine( tab, "", description );
+   }
+ 
+#undef LOCAL_SPRINTF
+#undef LOCAL_2_SPRINTF
+
+   return message.str();
+}
+

platforms/brook/src/BrookCalcProperDihedralBondForceKernel.h → platforms/brook/src/BrookBondParameters.h

-#ifndef OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_
-#define OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_
+#ifndef OPENMM_BROOK_BOND_PARAMETERS_H_
+#define OPENMM_BROOK_BOND_PARAMETERS_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -32,64 +32,24 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "kernels.h"
-#include "BrookPlatform.h"
-#include "BrookBondParameters.h"
-#include "OpenMMBrookInterface.h"
-
 namespace OpenMM {
 
 /**
- * This kernel is invoked to calculate the harmonic angle forces acting on the system.
+ * Container for bond parameters
  */
-
-class BrookCalcProperDihedralForceKernel : public CalcProperDihedralForceKernel {
+class BrookBondParameters {
 
    public:
   
-      /**
-       * BrookCalcProperDihedralForceKernel constructor
-       */
+     // return values
 
-      BrookCalcProperDihedralForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
-  
-      /**
-       * BrookCalcProperDihedralForceKernel destructor
-       */
+      static const int DefaultReturnValue = 0;
+      static const int ErrorReturnValue   = -1; 
 
-      ~BrookCalcProperDihedralForceKernel();
+      BrookBondParameters( std::string bondName, int numberOfParticlesInBond, int numberOfParametersInBond, int numberOfBonds, FILE* log );
   
-      /**
-       * Initialize the kernel, setting up the values to calculate harmonic bond force & energy
-       * 
-       * @param system                    System reference
-       * @param force                     ProperDihedralForce reference
-       *
-       */
-
-      void initialize( const System& system, const ProperDihedralForce& force );
+      ~BrookBondParameters();
   
-      /**
-       * Execute the kernel to calculate the forces.
-       * 
-       * @param positions   atom coordiantes
-       * @param forces      output forces
-       *
-       */
-
-      void executeForces( OpenMMContextImpl& context );
-  
-      /**
-       * Execute the kernel to calculate the energy.
-       * 
-       * @param context    the context in which to execute this kernel
-       *
-       * @return  potential energy associated with the harmonic angle force
-       *
-       */
-
-      double executeEnergy( OpenMMContextImpl& context );
-
       /** 
        * Set log file reference
        * 
@@ -122,51 +82,116 @@ class BrookCalcProperDihedralForceKernel : public CalcProperDihedralForceKernel
       FILE* getLog( void ) const;
       
       /** 
-       * Get number of bonds
+       * Set bond info
        * 
-       * @return  number of bonds
+       * @param bondIndex       index of bond
+       * @param particleIndices array of particle indices
+       * @param bondParameters  array of bond parameters
+       *
+       * @return  DefaultReturnValue
+       *
+       * @throw OpenMMException exeception if bond index is invalid
        *
        */
       
-      int getNumberOfBonds( void ) const;
+      int setBond( int bondIndex, int* particleIndices, double* bondParameters );
+      
+      /** 
+       * Get bond name
+       * 
+       * @return  bond name
+       *
+       */
+      
+      std::string getBondName( void ) const;
+      
+      /** 
+       * Set bond name
+       * 
+       * @param bondName       bond name
+       *
+       * @return  DefaultReturnValue
+       *
+       */
+      
+      //int setBondName( std::string bondName );
+      
+      /** 
+       * Get NumberOfParticlesInBond
+       * 
+       * @return NumberOfParticlesInBond
+       *
+       */
+      
+      int getNumberOfParticlesInBond( void ) const;
       
       /** 
-       * Get indices/parameters
+       * Get NumberOfParametersInBond
        * 
-       * @return  BrookBondParameters containing atom indices/parameters
+       * @return NumberOfParametersInBond
        *
        */
       
-      BrookBondParameters* getBrookBondParameters( void ) const;
+      int getNumberOfParametersInBond( void ) const;
+      
+      /** 
+       * Get NumberOfBonds
+       * 
+       * @return NumberOfBonds
+       *
+       */
       
+      int getNumberOfBonds( void ) const;
+      
+      /*  
+       * Get contents of object
+       *
+       *
+       * @param level   level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+    
+      std::string getContentsString( int level = 0 ) const;
+
    private:
    
-      static const int NumberOfAtomsInBond      = 3;
-      static const int NumberOfParametersInBond = 2;
-
-      // bond name
-
-      static const std::string BondName;
-
       // log file reference
 
       FILE* _log;
 
+      // bond name
 
-      // Brook bond parameters
-
-      BrookBondParameters* _brookBondParameters;
-
-      // interface
+      std::string _bondName;
 
-      OpenMMBrookInterface& _openMMBrookInterface;
+      // number of bonds
 
-      // System reference
+      int _numberOfBonds;
+      int _numberOfParticlesInBond;
+      int _numberOfParametersInBond;
+   
+      // particle indices and parameters
 
-      System& _system;
+      int** _particleIndices;
+      double** _bondParameters;
 
+      /* 
+       * Get contents of object
+       *
+       * @param tab         tab
+       * @param description description
+       * @param value       value
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string _getLine( const std::string& tab, const std::string& description,
+                            const std::string& value ) const;
+      
 };
 
 } // namespace OpenMM
 
-#endif /* OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_ */
+#endif /* OPENMM_BROOK_BOND_PARAMETERS_H_ */

platforms/brook/src/BrookBonded.cpp

@@ -35,13 +35,14 @@
 #include "BrookPlatform.h"
 #include "BrookStreamFactory.h"
 #include "OpenMMException.h"
+#include "gpu/kinvmap_gather.h"
 #include "gpu/invmap.h"
 #include "gpu/kforce.h"
 
 using namespace OpenMM;
 using namespace std;
 
-#define ATOMS(X,Y) (atoms[ 5*(X) + (Y) + 1 ])
+#define ATOMS(X,Y) (particles[ 5*(X) + (Y) + 1 ])
 #define PARAMS(X,Y,Z) (params[(Y)][4*(X) + Z])
 
 /** 
@@ -58,12 +59,12 @@ BrookBonded::BrookBonded( ){
 
 // ---------------------------------------------------------------------------------------
 
-   _numberOfAtoms             = 0;
+   _numberOfParticles         = 0;
    _ljScale                   = (BrookOpenMMFloat) 0.83333333;
    //_coulombFactor             = 332.0;
    _coulombFactor             = (BrookOpenMMFloat) 138.935485;
 
-   _atomIndicesStream         = NULL;
+   _particleIndicesStream     = NULL;
    _chargeStream              = NULL;
 
    // parameter streams
@@ -118,7 +119,7 @@ BrookBonded::~BrookBonded( ){
 
 // ---------------------------------------------------------------------------------------
 
-   delete _atomIndicesStream;
+   delete _particleIndicesStream;
    delete _chargeStream;
 
    for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
@@ -338,21 +339,21 @@ int BrookBonded::getInverseMapStreamCount( int index ) const {
 }
 
 /** 
- * Get bonded atom indices stream
+ * Get bonded particle indices stream
  * 
- * @return  atom indices stream
+ * @return  particle indices stream
  *
  */
 
-BrookFloatStreamInternal* BrookBonded::getAtomIndicesStream( void ) const {
+BrookFloatStreamInternal* BrookBonded::getParticleIndicesStream( void ) const {
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName = "BrookBonded::getAtomIndicesStream";
+   // static const std::string methodName = "BrookBonded::getParticleIndicesStream";
 
 // ---------------------------------------------------------------------------------------
 
-   return _atomIndicesStream;
+   return _particleIndicesStream;
 }
 
 /** 
@@ -450,7 +451,7 @@ BrookFloatStreamInternal** BrookBonded::getInverseStreamMapsStreams( int index )
 
 /*Flips i,j,k,l to l,k,j,i while correctly shuffling the params */
 
-void BrookBonded::flipQuartet( int ibonded, int *atoms ){
+void BrookBonded::flipQuartet( int ibonded, int *particles ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -473,7 +474,7 @@ void BrookBonded::flipQuartet( int ibonded, int *atoms ){
 
 }
 
-int BrookBonded::matchTorsion( int i, int j, int k, int l, int nbondeds, int *atoms ){
+int BrookBonded::matchTorsion( int i, int j, int k, int l, int nbondeds, int *particles ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -512,12 +513,12 @@ int BrookBonded::matchTorsion( int i, int j, int k, int l, int nbondeds, int *at
  * all dihedrals. But I think there are force fields that don't 
  * use all dihedrals.
  *
- * @return ErrorReturnValue if error; else atom index
+ * @return ErrorReturnValue if error; else particle index
  *
  **/
 
 int BrookBonded::matchAngle( int i, int j, int k, int nbondeds, 
-                             int *atoms, int *flag ){
+                             int *particles, int *flag ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -572,7 +573,7 @@ int BrookBonded::matchAngle( int i, int j, int k, int nbondeds,
  *
  * */
 
-int BrookBonded::matchBond( int i, int j, int nbondeds, int *atoms, int *flag ){
+int BrookBonded::matchBond( int i, int j, int nbondeds, int *particles, int *flag ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -638,7 +639,7 @@ int BrookBonded::matchBond( int i, int j, int nbondeds, int *atoms, int *flag ){
    return ErrorReturnValue;
 }
 
-int BrookBonded::matchPair( int i, int j, int nbondeds, int *atoms ){
+int BrookBonded::matchPair( int i, int j, int nbondeds, int *particles ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -652,7 +653,7 @@ int BrookBonded::matchPair( int i, int j, int nbondeds, int *atoms ){
 
       if( ATOMS(n, 0) == -1 ){
 
-         //If one of i,j matches the l atom
+         //If one of i,j matches the l particle
 
          if( ATOMS(n, 3) == i ){
             ATOMS(n, 0) = j;
@@ -664,7 +665,7 @@ int BrookBonded::matchPair( int i, int j, int nbondeds, int *atoms ){
          }
       }
 
-      //If the l-atom is available
+      //If the l-particle is available
       if( ATOMS(n, 3) == -1 ){
          if( ATOMS(n, 0) == i ){
             ATOMS(n, 3) = j;
@@ -686,12 +687,12 @@ int BrookBonded::matchPair( int i, int j, int nbondeds, int *atoms ){
 }
 
 /** 
- * Setup Ryckaert-Bellemans parameters/atom indices
+ * Setup Ryckaert-Bellemans parameters/particle indices
  * 
  * @param nbondeds                  number of bonded entries
- * @param atoms                     array of atom indices
+ * @param particles                     array of particle indices
  * @param params                    arrays of bond parameters
- * @param rbTorsionIndices          the four atoms connected by each Ryckaert-Bellemans torsion term
+ * @param rbTorsionIndices          the four particles connected by each Ryckaert-Bellemans torsion term
  * @param rbTorsionParameters       the coefficients (in order of increasing powers) for each Ryckaert-Bellemans torsion term
  * @param log                       log reference
  *
@@ -699,7 +700,7 @@ int BrookBonded::matchPair( int i, int j, int nbondeds, int *atoms ){
  *
  */
 
-int BrookBonded::addRBDihedrals( int *nbondeds, int *atoms, float *params[], 
+int BrookBonded::addRBDihedrals( int *nbondeds, int *particles, float *params[], 
                                  const vector<vector<int> >& rbTorsionIndices, 
                                  const vector<vector<double> >& rbTorsionParameters ){
 
@@ -716,16 +717,16 @@ int BrookBonded::addRBDihedrals( int *nbondeds, int *atoms, float *params[],
 
    for( unsigned int ii = 0; ii < rbTorsionIndices.size(); ii++ ){
 
-      vector<int> atomsIndices      = rbTorsionIndices[ii];
+      vector<int> particlesIndices      = rbTorsionIndices[ii];
       vector<double> rbParameters   = rbTorsionParameters[ii];
 
       int index = 0;
-      int i     = atomsIndices[index++];
-      int j     = atomsIndices[index++];
-      int k     = atomsIndices[index++];
-      int l     = atomsIndices[index++];
+      int i     = particlesIndices[index++];
+      int j     = particlesIndices[index++];
+      int k     = particlesIndices[index++];
+      int l     = particlesIndices[index++];
       
-      int ibonded = matchTorsion( i, j, k, l, *nbondeds, atoms );
+      int ibonded = matchTorsion( i, j, k, l, *nbondeds, particles );
       
       if( ibonded < 0 ){
          ibonded             = *nbondeds;
@@ -754,12 +755,12 @@ int BrookBonded::addRBDihedrals( int *nbondeds, int *atoms, float *params[],
 }
 
 /**
- * Setup periodic torsion parameters/atom indices
+ * Setup periodic torsion parameters/particle indices
  * 
  * @param nbondeds                  number of bonded entries
- * @param atoms                     array of atom indices
+ * @param particles                     array of particle indices
  * @param params                    arrays of bond parameters
- * @param periodicTorsionIndices    the four atoms connected by each periodic torsion term
+ * @param periodicTorsionIndices    the four particles connected by each periodic torsion term
  * @param periodicTorsionParameters the force parameters (k, phase, periodicity) for each periodic torsion term
  * @param log                       log reference
  *
@@ -767,7 +768,7 @@ int BrookBonded::addRBDihedrals( int *nbondeds, int *atoms, float *params[],
  *
  */
 
-int BrookBonded::addPDihedrals( int *nbondeds, int *atoms, BrookOpenMMFloat* params[],
+int BrookBonded::addPDihedrals( int *nbondeds, int *particles, BrookOpenMMFloat* params[],
                                 const vector<vector<int> >& periodicTorsionIndices, 
                                 const vector<vector<double> >& periodicTorsionParameters ){
 
@@ -784,16 +785,16 @@ int BrookBonded::addPDihedrals( int *nbondeds, int *atoms, BrookOpenMMFloat* par
 
    for( unsigned int ii = 0; ii < periodicTorsionIndices.size(); ii++ ){
 
-      vector<int> atomsIndices      = periodicTorsionIndices[ii];
+      vector<int> particlesIndices      = periodicTorsionIndices[ii];
       vector<double> pTParameters   = periodicTorsionParameters[ii];
 
       int index = 0;
-      int i     = atomsIndices[index++];
-      int j     = atomsIndices[index++];
-      int k     = atomsIndices[index++];
-      int l     = atomsIndices[index++];
+      int i     = particlesIndices[index++];
+      int j     = particlesIndices[index++];
+      int k     = particlesIndices[index++];
+      int l     = particlesIndices[index++];
 
-      int ibonded = matchTorsion( i, j, k, l, *nbondeds, atoms );
+      int ibonded = matchTorsion( i, j, k, l, *nbondeds, particles );
       
       if( ibonded < 0 ){
          ibonded = *nbondeds;
@@ -817,12 +818,12 @@ int BrookBonded::addPDihedrals( int *nbondeds, int *atoms, BrookOpenMMFloat* par
 }
 
 /**
- * Setup angle bond parameters/atom indices
+ * Setup angle bond parameters/particle indices
  * 
  * @param nbondeds                  number of bonded entries
- * @param atoms                     array of atom indices
+ * @param particles                     array of particle indices
  * @param params                    arrays of bond parameters
- * @param angleIndices              the angle bond atom indices
+ * @param angleIndices              the angle bond particle indices
  * @param angleParameters           the angle parameters (angle in radians, force constant)
  * @param log                       log reference
  *
@@ -830,7 +831,7 @@ int BrookBonded::addPDihedrals( int *nbondeds, int *atoms, BrookOpenMMFloat* par
  *
  */
 
-int BrookBonded::addAngles( int *nbondeds, int *atoms, float *params[], const std::vector<std::vector<int> >& angleIndices,
+int BrookBonded::addAngles( int *nbondeds, int *particles, float *params[], const std::vector<std::vector<int> >& angleIndices,
                             const std::vector<std::vector<double> >& angleParameters ){
 
 // ---------------------------------------------------------------------------------------
@@ -848,16 +849,16 @@ int BrookBonded::addAngles( int *nbondeds, int *atoms, float *params[], const st
 
    for( unsigned int ii = 0; ii < angleIndices.size(); ii++ ){
 
-      vector<int> atomsIndices       = angleIndices[ii];
+      vector<int> particlesIndices       = angleIndices[ii];
       vector<double> angParameters   = angleParameters[ii];
 
       int index = 0;
-      int i     = atomsIndices[index++];
-      int j     = atomsIndices[index++];
-      int k     = atomsIndices[index++];
+      int i     = particlesIndices[index++];
+      int j     = particlesIndices[index++];
+      int k     = particlesIndices[index++];
  
       int flag;
-      int ibonded = matchAngle( i, j, k, *nbondeds, atoms, &flag );
+      int ibonded = matchAngle( i, j, k, *nbondeds, particles, &flag );
       if( ibonded < 0 ){
          ibonded = *nbondeds;
          ATOMS( ibonded, 0 ) = i;
@@ -879,12 +880,12 @@ int BrookBonded::addAngles( int *nbondeds, int *atoms, float *params[], const st
 }
 
 /**
- * Setup harmonic bond parameters/atom indices
+ * Setup harmonic bond parameters/particle indices
  * 
  * @param nbondeds                  number of bonded entries
- * @param atoms                     array of atom indices
+ * @param particles                     array of particle indices
  * @param params                    arrays of bond parameters
- * @param bondIndices               two harmonic bond atom indices
+ * @param bondIndices               two harmonic bond particle indices
  * @param bondParameters            the force parameters (distance, k)
  * @param log                       log reference
  *
@@ -892,7 +893,7 @@ int BrookBonded::addAngles( int *nbondeds, int *atoms, float *params[], const st
  *
  */
 
-int BrookBonded::addBonds( int *nbondeds, int *atoms, float *params[], const vector<vector<int> >& bondIndices,
+int BrookBonded::addBonds( int *nbondeds, int *particles, float *params[], const vector<vector<int> >& bondIndices,
                            const vector<vector<double> >& bondParameters ){
 
 // ---------------------------------------------------------------------------------------
@@ -910,12 +911,12 @@ int BrookBonded::addBonds( int *nbondeds, int *atoms, float *params[], const vec
 
    for( unsigned int ii = 0; ii < bondIndices.size(); ii++ ){
 
-      vector<int> atomsIndices       = bondIndices[ii];
+      vector<int> particlesIndices       = bondIndices[ii];
       vector<double> bndParameters   = bondParameters[ii];
 
       int index = 0;
-      int i     = atomsIndices[index++];
-      int j     = atomsIndices[index++];
+      int i     = particlesIndices[index++];
+      int j     = particlesIndices[index++];
 
       // insure i < j
 
@@ -926,7 +927,7 @@ int BrookBonded::addBonds( int *nbondeds, int *atoms, float *params[], const vec
       }
 
       int flag;
-      int ibonded = matchBond( i, j, *nbondeds, atoms, &flag );
+      int ibonded = matchBond( i, j, *nbondeds, particles, &flag );
 int saveIbond = ibonded;
       if( ibonded < 0 ){
          ibonded = *nbondeds;
@@ -954,15 +955,15 @@ int saveIbond = ibonded;
 }
 
 /**
- * Setup LJ/Coulomb 1-4 parameters/atom indices
+ * Setup LJ/Coulomb 1-4 parameters/particle indices
  * 
  * @param nbondeds                  number of bonded entries
- * @param atoms                     array of atom indices
+ * @param particles                     array of particle indices
  * @param params                    arrays of bond parameters
  * @param charges                   array of charges
- * @param bonded14Indices           each element contains the indices of two atoms whose nonbonded interactions should be reduced since
+ * @param bonded14Indices           each element contains the indices of two particles whose nonbonded interactions should be reduced since
  *                                  they form a bonded 1-4 pair
- * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+ * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each particle
  * @param lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
  * @param log                       log reference
  *
@@ -970,7 +971,7 @@ int saveIbond = ibonded;
  *
  */
 
-int BrookBonded::addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[],
+int BrookBonded::addPairs( int *nbondeds, int *particles, BrookOpenMMFloat* params[],
                            BrookOpenMMFloat* charges,
                            const std::vector<std::vector<int> >& bonded14Indices,
                            const std::vector<std::vector<double> >& nonbondedParameters,
@@ -990,13 +991,13 @@ int BrookBonded::addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[]
 
    for( unsigned int ii = 0; ii < bonded14Indices.size(); ii++ ){
 
-      std::vector<int> atomsIndices        = bonded14Indices[ii];
+      std::vector<int> particlesIndices        = bonded14Indices[ii];
 
       int index = 0;
-      int i     = atomsIndices[index++];
-      int j     = atomsIndices[index++];
+      int i     = particlesIndices[index++];
+      int j     = particlesIndices[index++];
 
-      int ibonded = matchPair( i, j, *nbondeds, atoms );
+      int ibonded = matchPair( i, j, *nbondeds, particles );
       if( ibonded < 0 ){
          ibonded = *nbondeds;
          ATOMS(ibonded, 0) = i;
@@ -1042,8 +1043,8 @@ int BrookBonded::addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[]
  * Create and load inverse maps for bonded ixns
  * 
  * @param nbondeds                  number of bonded entries
- * @param natoms                    number of atoms
- * @param atoms                     arrays of atom indices (atoms[numberOfBonds][4])
+ * @param nparticles                    number of particles
+ * @param particles                     arrays of particle indices (particles[numberOfBonds][4])
  * @param platform                  BrookPlatform reference
  * @param log                       log file reference (optional)
  *
@@ -1051,7 +1052,7 @@ int BrookBonded::addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[]
  *
  */
 
-int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookPlatform& brookPlatform ){
+int BrookBonded::loadInvMaps( int nbondeds, int nparticles, int *particles, const BrookPlatform& brookPlatform ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -1061,28 +1062,28 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
 
 // ---------------------------------------------------------------------------------------
 
-   // get atom stream size
+   // get particle stream size
 
    const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (brookPlatform.getDefaultStreamFactory() );
-   int atomStreamWidth                          = brookStreamFactory.getDefaultAtomStreamWidth();
-   int atomStreamSize                           = brookPlatform.getStreamSize( getNumberOfAtoms(), atomStreamWidth, NULL );
-   _inverseMapStreamWidth                       = atomStreamWidth;
+   int particleStreamWidth                          = brookStreamFactory.getDefaultParticleStreamWidth();
+   int particleStreamSize                           = brookPlatform.getStreamSize( getNumberOfParticles(), particleStreamWidth, NULL );
+   _inverseMapStreamWidth                       = particleStreamWidth;
    
 // ---------------------------------------------------------------------------------------
 
    // allocate temp memory
 
    float4** invmaps = new float4*[getMaxInverseMapStreamCount()];
-   float* block     = new float[4*getMaxInverseMapStreamCount()*atomStreamSize];
+   float* block     = new float[4*getMaxInverseMapStreamCount()*particleStreamSize];
 
-   //memset( block, 0, 4*getMaxInverseMapStreamCount()*atomStreamSize*sizeof( float ) );
+   //memset( block, 0, 4*getMaxInverseMapStreamCount()*particleStreamSize*sizeof( float ) );
 
    float* blockPtr = block;
    for( int ii = 0; ii < getMaxInverseMapStreamCount(); ii++ ){
       invmaps[ii]  = (float4*) blockPtr;
-      blockPtr    += 4*atomStreamSize;
+      blockPtr    += 4*particleStreamSize;
    }
-   int* counts = new int[atomStreamSize];
+   int* counts = new int[particleStreamSize];
 
 // ---------------------------------------------------------------------------------------
 
@@ -1096,26 +1097,26 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
 
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
       for( int jj = 0; jj < getMaxInverseMapStreamCount(ii); jj++ ){
-         _inverseStreamMaps[ii][jj] = new BrookFloatStreamInternal( BrookCommon::BondedInverseMapStreams, atomStreamSize,
-                                                                    atomStreamWidth, BrookStreamInternal::Float4, dangleValue );
+         _inverseStreamMaps[ii][jj] = new BrookFloatStreamInternal( BrookCommon::BondedInverseMapStreams, particleStreamSize,
+                                                                    particleStreamWidth, BrookStreamInternal::Float4, dangleValue );
       }
    }
 
    if( PrintOn && getLog() ){
       (void) fprintf( getLog(), "%s force stream strms=%d nbondeds=%d max counts=[%d %d %d %d] strSz&Wd=%d %d\n", methodName.c_str(), getNumberOfForceStreams(),
                       nbondeds, getMaxInverseMapStreamCount(0), getMaxInverseMapStreamCount(1), getMaxInverseMapStreamCount(2), getMaxInverseMapStreamCount(3),
-                      atomStreamSize, atomStreamWidth );
+                      particleStreamSize, particleStreamWidth );
       (void) fflush( getLog() );
    }
 
    // load data
 
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
-      for( int jj = 0; jj < 4*getMaxInverseMapStreamCount()*atomStreamSize; jj++ ){
+      for( int jj = 0; jj < 4*getMaxInverseMapStreamCount()*particleStreamSize; jj++ ){
          block[jj] = -1.0f;
       }
-      gpuCalcInvMap( ii, 4, nbondeds, natoms, atoms, getInverseMapStreamCount( ii ), counts, invmaps, &(_inverseMapStreamCount[ii]) );
-//gpuPrintInvMaps( _inverseMapStreamCount[ii], natoms, counts, invmaps, getLog() );
+      gpuCalcInvMap( ii, 4, nbondeds, nparticles, particles, getInverseMapStreamCount( ii ), counts, invmaps, &(_inverseMapStreamCount[ii]) );
+//gpuPrintInvMaps( _inverseMapStreamCount[ii], nparticles, counts, invmaps, getLog() );
       validateInverseMapStreamCount( ii, _inverseMapStreamCount[ii] ); 
       for( int jj = 0; jj < _inverseMapStreamCount[ii]; jj++ ){
          _inverseStreamMaps[ii][jj]->loadFromArray( invmaps[jj] );
@@ -1125,7 +1126,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
                             methodName.c_str(), getNumberOfForceStreams(), _inverseMapStreamCount[ii], ii, jj,
                             getInverseMapStreamCount( ii ), getMaxInverseMapStreamCount( ii ) );
 
-            for( int kk = 0; kk < atomStreamSize; kk++ ){
+            for( int kk = 0; kk < particleStreamSize; kk++ ){
                (void) fprintf( getLog(), "%8d [ %.1f %.1f %.1f %.1f]\n", kk, invmaps[jj][kk].x, invmaps[jj][kk].y, invmaps[jj][kk].z, invmaps[jj][kk].w  );
             }
          }
@@ -1136,7 +1137,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
       // keep invalid entries from being included in forces
 
       if( _inverseMapStreamCount[ii] < getMaxInverseMapStreamCount( ii ) ){
-         for( int jj = 0; jj < 4*atomStreamSize; jj++ ){
+         for( int jj = 0; jj < 4*particleStreamSize; jj++ ){
             block[jj] = -1.0f;
          }
          for( int jj = _inverseMapStreamCount[ii]; jj < getMaxInverseMapStreamCount( ii ); jj++ ){
@@ -1165,16 +1166,16 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
 /* 
  * Setup for bonded ixns
  *
- * @param numberOfAtoms                number of atoms
- * @param bondIndices                  vector of vector of harmonic                   bond indices    -- one entry each bond (2 atoms     )
+ * @param numberOfParticles            number of particles
+ * @param bondIndices                  vector of vector of harmonic                   bond indices    -- one entry each bond (2 particles     )
  * @param bondParameters               vector of vector of harmonic                   bond parameters -- one entry each bond (2 parameters)
- * @param angleIndices                 vector of vector of angle                      bond indices    -- one entry each bond (3 atoms     )
+ * @param angleIndices                 vector of vector of angle                      bond indices    -- one entry each bond (3 particles     )
  * @param angleParameters              vector of vector of angle                      bond parameters -- one entry each bond (2 parameters)
- * @param periodicTorsionIndices       vector of vector of periodicTorsionIndices     bond indices    -- one entry each bond (4 atoms     )
+ * @param periodicTorsionIndices       vector of vector of periodicTorsionIndices     bond indices    -- one entry each bond (4 particles     )
  * @param periodicTorsionParameters    vector of vector of periodicTorsionParameters  bond parameters -- one entry each bond (3 parameters)
- * @param rbTorsionIndices             vector of vector of rb torsion                 bond indices    -- one entry each bond (4 atoms     )
+ * @param rbTorsionIndices             vector of vector of rb torsion                 bond indices    -- one entry each bond (4 particles     )
  * @param rbTorsionParameters          vector of vector of rb torsion                 bond parameters -- one entry each bond (5 parameters)
- * @param bonded14Indices              vector of vector of Lennard-Jones 14           atom indices    -- one entry each bond (2 atoms     )
+ * @param bonded14Indices              vector of vector of Lennard-Jones 14           particle indices    -- one entry each bond (2 particles     )
  * @param nonbondedParameters          vector of vector of Lennard-Jones 14           parameters      -- one entry each bond (3 parameters)
  * @param lj14Scale                    scaling factor for 1-4 ixns
  * @param coulombScale                 Coulomb scaling factor for 1-4 ixns
@@ -1192,7 +1193,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
  * the optimal fit, but should not be too bad. 
  * */
 
-int BrookBonded::setup( int numberOfAtoms,
+int BrookBonded::setup( int numberOfParticles,
                         const vector<vector<int> >& bondIndices,            const vector<vector<double> >& bondParameters,
                         const vector<vector<int> >& angleIndices,           const vector<vector<double> >& angleParameters,
                         const vector<vector<int> >& periodicTorsionIndices, const vector<vector<double> >& periodicTorsionParameters,
@@ -1208,15 +1209,15 @@ int BrookBonded::setup( int numberOfAtoms,
 
 // ---------------------------------------------------------------------------------------
 
-   _numberOfAtoms = numberOfAtoms;
+   _numberOfParticles = numberOfParticles;
 
    const BrookPlatform& brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
 
-   // check that atom indices & parameters agree
+   // check that particle indices & parameters agree
 
    if( bondIndices.size() != bondParameters.size() ){
       std::stringstream message;
-      message << methodName << " number of harmonic bond atom indices=" << bondIndices.size() << " does not equal number of harmonic bond parameter entries=" << bondParameters.size();
+      message << methodName << " number of harmonic bond particle indices=" << bondIndices.size() << " does not equal number of harmonic bond parameter entries=" << bondParameters.size();
       throw OpenMMException( message.str() );
    } else if( PrintOn && getLog() ){
       (void) fprintf( getLog(), "%s harmonic bonds=%d\n", methodName.c_str(), bondIndices.size() );
@@ -1225,7 +1226,7 @@ int BrookBonded::setup( int numberOfAtoms,
 
    if( angleIndices.size() != angleParameters.size() ){
       std::stringstream message;
-      message << methodName << " number of angle atom indices=" << angleIndices.size() << " does not equal number of angle parameter entries=" << angleParameters.size();
+      message << methodName << " number of angle particle indices=" << angleIndices.size() << " does not equal number of angle parameter entries=" << angleParameters.size();
       throw OpenMMException( message.str() );
    } else if( PrintOn && getLog() ){
       (void) fprintf( getLog(), "%s angle bonds=%d\n", methodName.c_str(), angleIndices.size() );
@@ -1234,7 +1235,7 @@ int BrookBonded::setup( int numberOfAtoms,
 
    if( periodicTorsionIndices.size() != periodicTorsionParameters.size() ){
       std::stringstream message;
-      message << methodName << " number of periodicTorsion atom indices=" << periodicTorsionIndices.size() << " does not equal number of periodicTorsion parameter entries=" << periodicTorsionParameters.size();
+      message << methodName << " number of periodicTorsion particle indices=" << periodicTorsionIndices.size() << " does not equal number of periodicTorsion parameter entries=" << periodicTorsionParameters.size();
       throw OpenMMException( message.str() );
    } else if( PrintOn && getLog() ){
       (void) fprintf( getLog(), "%s periodicTorsion bonds=%d\n", methodName.c_str(), periodicTorsionIndices.size() );
@@ -1243,16 +1244,16 @@ int BrookBonded::setup( int numberOfAtoms,
 
    if( rbTorsionIndices.size() != rbTorsionParameters.size() ){
       std::stringstream message;
-      message << methodName << " number of rbTorsion atom indices=" << rbTorsionIndices.size() << " does not equal number of rbTorsion parameter entries=" << rbTorsionParameters.size();
+      message << methodName << " number of rbTorsion particle indices=" << rbTorsionIndices.size() << " does not equal number of rbTorsion parameter entries=" << rbTorsionParameters.size();
       throw OpenMMException( message.str() );
    } else if( PrintOn && getLog() ){
       (void) fprintf( getLog(), "%s rbTorsion bonds=%d\n", methodName.c_str(), rbTorsionIndices.size() );
       (void) fflush( getLog() );
    }
 
-   if( (numberOfAtoms != (int) nonbondedParameters.size()) && bonded14Indices.size() > 0 ){
+   if( (numberOfParticles != (int) nonbondedParameters.size()) && bonded14Indices.size() > 0 ){
       std::stringstream message;
-      message << methodName << " number atoms=" << numberOfAtoms << " does not equal number of nb parameter entries=" << nonbondedParameters.size();
+      message << methodName << " number particles=" << numberOfParticles << " does not equal number of nb parameter entries=" << nonbondedParameters.size();
       throw OpenMMException( message.str() );
    } else if( PrintOn && getLog() ){
       (void) fprintf( getLog(), "%s LJ 14 ixns=%d\n", methodName.c_str(), bonded14Indices.size() );
@@ -1261,8 +1262,8 @@ int BrookBonded::setup( int numberOfAtoms,
 
    // allocate temp memory
 
-   int maxBonds = 10*numberOfAtoms;
-   int* atoms   = new int[5*maxBonds];
+   int maxBonds = 10*numberOfParticles;
+   int* particles   = new int[5*maxBonds];
 
    BrookOpenMMFloat** params = new BrookOpenMMFloat*[getNumberOfParameterStreams()];
    for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
@@ -1274,9 +1275,9 @@ int BrookBonded::setup( int numberOfAtoms,
    // build streams
 
    const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (brookPlatform.getDefaultStreamFactory() );
-   int atomStreamWidth                          = brookStreamFactory.getDefaultAtomStreamWidth();
+   int particleStreamWidth                          = brookStreamFactory.getDefaultParticleStreamWidth();
 
-   // Initialize all atom indices to -1 to indicate empty slots
+   // Initialize all particle indices to -1 to indicate empty slots
    // All parameters must be initialized to values that will 
    // produce zero for the corresponding force. 
 
@@ -1300,25 +1301,25 @@ int BrookBonded::setup( int numberOfAtoms,
 
    // nbondeds tracks number of ixn
    int nbondeds = 0;
-   addRBDihedrals( &nbondeds, atoms, params, rbTorsionIndices,       rbTorsionParameters       );
-   addPDihedrals ( &nbondeds, atoms, params, periodicTorsionIndices, periodicTorsionParameters );
-   addAngles(      &nbondeds, atoms, params, angleIndices,           angleParameters           );
-   addBonds(       &nbondeds, atoms, params, bondIndices,            bondParameters            );
+   addRBDihedrals( &nbondeds, particles, params, rbTorsionIndices,       rbTorsionParameters       );
+   addPDihedrals ( &nbondeds, particles, params, periodicTorsionIndices, periodicTorsionParameters );
+   addAngles(      &nbondeds, particles, params, angleIndices,           angleParameters           );
+   addBonds(       &nbondeds, particles, params, bondIndices,            bondParameters            );
 
 // ---------------------------------------------------------------------------------------
 
    // charge stream
 
-   _chargeStream         = new BrookFloatStreamInternal( BrookCommon::BondedChargeStream, numberOfAtoms, atomStreamWidth,
+   _chargeStream         = new BrookFloatStreamInternal( BrookCommon::BondedChargeStream, numberOfParticles, particleStreamWidth,
                                                          BrookStreamInternal::Float, dangleValue );
    BrookOpenMMFloat* charges = new BrookOpenMMFloat[_chargeStream->getStreamSize()];
    memset( charges, 0, _chargeStream->getStreamSize()*sizeof( BrookOpenMMFloat ) );
 
 // ---------------------------------------------------------------------------------------
 
-//(void) fprintf( getLog(), "%s Post addBonds atoms=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfAtoms, nbondeds, maxBonds );
+//(void) fprintf( getLog(), "%s Post addBonds particles=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfParticles, nbondeds, maxBonds );
 
-   addPairs( &nbondeds, atoms, params, charges, bonded14Indices, nonbondedParameters, lj14Scale, coulombScale );
+   addPairs( &nbondeds, particles, params, charges, bonded14Indices, nonbondedParameters, lj14Scale, coulombScale );
 
    // check that number of bonds not too large for memory allocated
 
@@ -1327,28 +1328,28 @@ int BrookBonded::setup( int numberOfAtoms,
       message << methodName << " number of bonds=" << nbondeds << " is greater than maxBonds=" << maxBonds;
       throw OpenMMException( message.str() );
    } else if( PrintOn && getLog() ){
-      (void) fprintf( getLog(), "%s atoms=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfAtoms, nbondeds, maxBonds );
+      (void) fprintf( getLog(), "%s particles=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfParticles, nbondeds, maxBonds );
       (void) fflush( getLog() );
    }
 
 // ---------------------------------------------------------------------------------------
 
-   // atom indices stream
+   // particle indices stream
 
-   _atomIndicesStream    = new BrookFloatStreamInternal( BrookCommon::BondedAtomIndicesStream, nbondeds, atomStreamWidth,
-                                                         BrookStreamInternal::Float4, dangleValue );
+   _particleIndicesStream    = new BrookFloatStreamInternal( BrookCommon::BondedParticleIndicesStream, nbondeds, particleStreamWidth,
+                                                             BrookStreamInternal::Float4, dangleValue );
 
-   int* buffer           = new int[4*_atomIndicesStream->getStreamSize()];
-   memset( buffer, 0, sizeof( int )*4*_atomIndicesStream->getStreamSize() );
+   int* buffer           = new int[4*_particleIndicesStream->getStreamSize()];
+   memset( buffer, 0, sizeof( int )*4*_particleIndicesStream->getStreamSize() );
 
    int index             = 0;
    for( int ii = 0; ii < nbondeds; ii++ ){
       for( int jj = 0; jj < 4; jj++ ){
          buffer[index++] = ATOMS( ii, jj );
-//(void) fprintf( getLog(), "%s atomIndices %d %d  %d buffer=%d atoms=%d\n", methodName.c_str(), ii, jj, index, buffer[index-1], ATOMS( ii, jj ) );
+//(void) fprintf( getLog(), "%s particleIndices %d %d  %d buffer=%d particles=%d\n", methodName.c_str(), ii, jj, index, buffer[index-1], ATOMS( ii, jj ) );
       }
    }
-   _atomIndicesStream->loadFromArray( buffer, BrookStreamInternal::Integer ); 
+   _particleIndicesStream->loadFromArray( buffer, BrookStreamInternal::Integer ); 
    delete[] buffer;
 
 // ---------------------------------------------------------------------------------------
@@ -1360,7 +1361,7 @@ int BrookBonded::setup( int numberOfAtoms,
    // bonded parameters
 
    for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
-      _bondedParameters[ii]  = new BrookFloatStreamInternal( BrookCommon::BondedParametersStream, nbondeds, atomStreamWidth,
+      _bondedParameters[ii]  = new BrookFloatStreamInternal( BrookCommon::BondedParametersStream, nbondeds, particleStreamWidth,
                                                              BrookStreamInternal::Float4, dangleValue );
       _bondedParameters[ii]->loadFromArray( params[ii] );
    }
@@ -1372,9 +1373,9 @@ int BrookBonded::setup( int numberOfAtoms,
    if( PrintOn && getLog() ){
 
       (void) fprintf( getLog(), "%s nbondeds=%d strDim [%d %d ] sz=%d\n", methodName.c_str(), nbondeds,
-                      _atomIndicesStream->getStreamWidth(),
-                      _atomIndicesStream->getStreamHeight(), 
-                      _atomIndicesStream->getStreamSize() );
+                      _particleIndicesStream->getStreamWidth(),
+                      _particleIndicesStream->getStreamHeight(), 
+                      _particleIndicesStream->getStreamSize() );
 
       int kIndex = 0;
       int jIndex = 1;
@@ -1411,7 +1412,7 @@ int BrookBonded::setup( int numberOfAtoms,
 
    // load inverse maps to streams
 
-   loadInvMaps( nbondeds, getNumberOfAtoms(), atoms, brookPlatform );
+   loadInvMaps( nbondeds, getNumberOfParticles(), particles, brookPlatform );
    
 // ---------------------------------------------------------------------------------------
 
@@ -1422,7 +1423,7 @@ int BrookBonded::setup( int numberOfAtoms,
    }
 
    delete[] params;
-   delete[] atoms;
+   delete[] particles;
    delete[] charges;
 
    // set the fudge factors
@@ -1433,7 +1434,7 @@ int BrookBonded::setup( int numberOfAtoms,
    // initialize output streams
 
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
-      _bondedForceStreams[ii] = new BrookFloatStreamInternal( BrookCommon::UnrolledForceStream, nbondeds, atomStreamWidth, 
+      _bondedForceStreams[ii] = new BrookFloatStreamInternal( BrookCommon::UnrolledForceStream, nbondeds, particleStreamWidth, 
                                                               BrookStreamInternal::Float3, dangleValue );
    }
 
@@ -1474,8 +1475,8 @@ std::string BrookBonded::getContentsString( int level ) const {
 #define LOCAL_2_SPRINTF(a,b,c,d) sprintf( (a), (b), (c), (d) );   
 #endif
 
-   (void) LOCAL_SPRINTF( value, "%d", getNumberOfAtoms() );
-   message << _getLine( tab, "Number of atoms:", value ); 
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfParticles() );
+   message << _getLine( tab, "Number of particles:", value ); 
 
    (void) LOCAL_SPRINTF( value, "%.5f", getLJ_14Scale() );
    message << _getLine( tab, "LJ 14 scaling:", value ); 
@@ -1487,14 +1488,14 @@ std::string BrookBonded::getContentsString( int level ) const {
    message << _getLine( tab, "Inverse map stream width:", value ); 
 
 /*
-   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamWidth() );
-   message << _getLine( tab, "Atom stream width:", value ); 
+   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamWidth() );
+   message << _getLine( tab, "Particle stream width:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamHeight() );
-   message << _getLine( tab, "Atom stream height:", value ); 
+   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamHeight() );
+   message << _getLine( tab, "Particle stream height:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamSize() );
-   message << _getLine( tab, "Atom stream size:", value ); 
+   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamSize() );
+   message << _getLine( tab, "Particle stream size:", value ); 
 */
 
    (void) LOCAL_SPRINTF( value, "%d", getNumberOfParameterStreams() );
@@ -1518,7 +1519,7 @@ std::string BrookBonded::getContentsString( int level ) const {
    }
 
    message << _getLine( tab, "Log:",                 (getLog()               ? Set : NotSet) ); 
-   message << _getLine( tab, "Atom indices stream:", (getAtomIndicesStream() ? Set : NotSet) ); 
+   message << _getLine( tab, "Particle indices stream:", (getParticleIndicesStream() ? Set : NotSet) ); 
    //message << _getLine( tab, "Charge stream:",        (getChargeStream()      ? Set : NotSet) ); 
  
    (void) LOCAL_SPRINTF( value, "%d", getNumberOfForceStreams() );
@@ -1540,23 +1541,23 @@ std::string BrookBonded::getContentsString( int level ) const {
  * Helper functions for building inverse maps for 
  * torsions, impropers and angles.
  * 
- * For each atom, calculates the positions at which it's
+ * For each particle, calculates the positions at which it's
  * forces are to be picked up from and stores the position
  * in the appropriate index.
  *
- * Input: number of dihedrals, the atom indices, and a flag indicating
+ * Input: number of dihedrals, the particle indices, and a flag indicating
  *        whether we're doing i(0), j(1), k(2) or l(3)
- * Output: an array of counts per atom
+ * Output: an array of counts per particle
  *         arrays of inversemaps
  *         nimaps - the number of invmaps actually used.
  *
- * @param posflag       0-niatoms-1
- * @param niatoms       3 for angles, 4 for torsions, impropers
+ * @param posflag       0-niparticles-1
+ * @param niparticles       3 for angles, 4 for torsions, impropers
  * @param nints         number of interactions
- * @param natoms        number of atoms
- * @param *atoms        gromacs interaction list
+ * @param nparticles        number of particles
+ * @param *particles        gromacs interaction list
  * @param nmaps         maximum number of inverse maps
- * @param   counts[]    output counts of how many places each atom occurs
+ * @param   counts[]    output counts of how many places each particle occurs
  * @param *invmaps[]    output array of nmaps inverse maps
  * @param *nimaps,      output max number of inverse maps actually used
  *
@@ -1564,14 +1565,14 @@ std::string BrookBonded::getContentsString( int level ) const {
  *
  **/
 
-int BrookBonded::gpuCalcInvMap( int posflag, int niatoms, int nints, int natoms,
-                                int *atoms, int nmaps, int counts[], float4 *invmaps[],
+int BrookBonded::gpuCalcInvMap( int posflag, int niparticles, int nints, int nparticles,
+                                int *particles, int nmaps, int counts[], float4 *invmaps[],
                                 int *nimaps ){
 
 // ---------------------------------------------------------------------------------------
 
    int i, j;
-   int atom;
+   int particle;
    int mapnum, mapcomp;
 
    static const std::string methodName      = "BrookBonded::gpuCalcInvMap";
@@ -1584,10 +1585,10 @@ int BrookBonded::gpuCalcInvMap( int posflag, int niatoms, int nints, int natoms,
 
 // ---------------------------------------------------------------------------------------
 
-   memset( counts, 0, sizeof( int )*natoms );
+   memset( counts, 0, sizeof( int )*nparticles );
 
    for( i = 0; i < nmaps; i++ ){
-      for( j = 0; j < natoms; j++ ){
+      for( j = 0; j < nparticles; j++ ){
          invmaps[i][j] = float4( -1.0, -1.0, -1.0, -1.0 );
       }
    }
@@ -1596,61 +1597,61 @@ int BrookBonded::gpuCalcInvMap( int posflag, int niatoms, int nints, int natoms,
 
    *nimaps = -1;
 
-   //Now note down the positions where each atom occurs
+   //Now note down the positions where each particle occurs
 
    if( PrintOn && getLog() ){
-      (void) fprintf( getLog(), "%s: pos=%d ni=%d nints=%d natoms=%d nmaps=<%d>\n", methodName.c_str(), posflag, niatoms, nints, natoms, nmaps ); 
+      (void) fprintf( getLog(), "%s: pos=%d ni=%d nints=%d nparticles=%d nmaps=<%d>\n", methodName.c_str(), posflag, niparticles, nints, nparticles, nmaps ); 
       (void) fflush( getLog() );
    }
 
-int atomRange[2]   = { 90000000, -90000000 };
+int particleRange[2]   = { 90000000, -90000000 };
 int mapnumRange[2] = { 90000000, -90000000 };
 
    for(  i = 0; i < nints; i++ ){
-      //This is our atom
-      atom = atoms[ (niatoms + 1) * i + posflag + 1 ];
+      //This is our particle
+      particle = particles[ (niparticles + 1) * i + posflag + 1 ];
 
       //Special for merged bondeds
-      if ( atom == -1 ){
+      if ( particle == -1 ){
          continue;
       }
 
-if( atom < atomRange[0] ){
-   atomRange[0] = atom;
+if( particle < particleRange[0] ){
+   particleRange[0] = particle;
 }
-if( atom > atomRange[1] ){
-   atomRange[1] = atom;
+if( particle > particleRange[1] ){
+   particleRange[1] = particle;
 }
       //Check to make sure we're inside the limits
-      if ( counts[atom] > nmaps * 4 ){
+      if ( counts[particle] > nmaps * 4 ){
          if( PrintOn && getLog() ){
-            (void) fprintf( getLog(), "%s Atom %d has too many proper dihedrals(%d, max %d)\n",
-                            methodName.c_str(), atom, counts[atom], nmaps*4 );
+            (void) fprintf( getLog(), "%s Particle %d has too many proper dihedrals(%d, max %d)\n",
+                            methodName.c_str(), particle, counts[particle], nmaps*4 );
             (void) fflush( getLog() );
          }
          std::stringstream message;
-         message << methodName << " Atom " << atom << " has too many proper dihedrals; valid range:(" << counts[atom] << ", " << nmaps*4 << ")";
+         message << methodName << " Particle " << particle << " has too many proper dihedrals; valid range:(" << counts[particle] << ", " << nmaps*4 << ")";
          throw OpenMMException( message.str() );
       }
       
       //Which invmap will this go into
 
-      mapnum = counts[atom] / 4;
+      mapnum = counts[particle] / 4;
 
       if ( mapnum > *nimaps )
          *nimaps = mapnum;
 
       //Which component will it be
-      mapcomp = counts[atom] % 4;
+      mapcomp = counts[particle] % 4;
 
       //Set it
       //This is silly, but otherwise I have to declare it as float*
       //and things get even more confusing. :)
       switch (mapcomp){
-         case 0: invmaps[mapnum][atom].x = (float) i; break;
-         case 1: invmaps[mapnum][atom].y = (float) i; break;
-         case 2: invmaps[mapnum][atom].z = (float) i; break;
-         case 3: invmaps[mapnum][atom].w = (float) i; break;
+         case 0: invmaps[mapnum][particle].x = (float) i; break;
+         case 1: invmaps[mapnum][particle].y = (float) i; break;
+         case 2: invmaps[mapnum][particle].z = (float) i; break;
+         case 3: invmaps[mapnum][particle].w = (float) i; break;
          default:
             if( PrintOn && getLog() ){
                (void) fprintf( getLog(), "mapcomp %d invalid -- impossible!\n", mapcomp );
@@ -1662,7 +1663,7 @@ if( atom > atomRange[1] ){
             break;
       }
       
-      counts[atom]++;
+      counts[particle]++;
 
 if( mapnum < mapnumRange[0] ){
    mapnumRange[0] = mapnum;
@@ -1671,15 +1672,15 @@ if( mapnum > mapnumRange[1] ){
    mapnumRange[1] = mapnum;
 }
 
-//fprintf( gpu->log, "%d atom=%d  mapcomp=%d counts[]=%d mapnum=%d\n", i, atom, mapcomp, counts[atom], mapnum );
+//fprintf( gpu->log, "%d particle=%d  mapcomp=%d counts[]=%d mapnum=%d\n", i, particle, mapcomp, counts[particle], mapnum );
 
    }
 
    (*nimaps)++;
 
 if( PrintOn && getLog() ){
-   (void) fprintf( getLog(), "%s mnmaps=%d Ranges: atom [%d %d] mapnum [%d %d]\n",
-                   methodName.c_str(), *nimaps, atomRange[0], atomRange[1], mapnumRange[0], mapnumRange[1] );
+   (void) fprintf( getLog(), "%s mnmaps=%d Ranges: particle [%d %d] mapnum [%d %d]\n",
+                   methodName.c_str(), *nimaps, particleRange[0], particleRange[1], mapnumRange[0], mapnumRange[1] );
    (void) fflush( getLog() );
 }
 
@@ -1687,10 +1688,10 @@ if( PrintOn && getLog() ){
 }
 
 
-void BrookBonded::gpuPrintInvMaps( int nmaps, int natoms, int counts[], float4 *invmap[], FILE* logFile ){
+void BrookBonded::gpuPrintInvMaps( int nmaps, int nparticles, int counts[], float4 *invmap[], FILE* logFile ){
    int i;
    int j;
-   for(  i = 0; i < natoms; i++ ){
+   for(  i = 0; i < nparticles; i++ ){
       fprintf( logFile, "%d %d ", i, counts[i] );
       for(  j = 0; j < nmaps; j++ ){
          fprintf( logFile, "%6.0f %6.0f %6.0f %6.0f", invmap[j][i].x, invmap[j][i].y, 
@@ -1708,25 +1709,26 @@ void BrookBonded::gpuPrintInvMaps( int nmaps, int natoms, int counts[], float4 *
  * lookup. This assumes that nints < 100000, preferably nints << 100000
  * which should always be true
  * */
+
 int BrookBonded::gpuCalcInvMap_merged( 
       int nints,    //number of interactions
-      int natoms,   //number of atoms
-      int *atoms,   //ijkl,ijkl,ijkl...
+      int nparticles,   //number of particles
+      int *particles,   //ijkl,ijkl,ijkl...
       int nmaps,      //maximum number of inverse maps
-        int counts[],   //output counts of how many places each atom occurs
+        int counts[],   //output counts of how many places each particle occurs
       float4 *invmaps[], //output array of nmaps inverse maps
       int *nimaps        //output max number of inverse maps actually used
       ){
    int i, j;
-   int atom;
+   int particle;
    int mapnum, mapcomp;
    int pos;
    
-   for(  i = 0; i < natoms; i++ )
+   for(  i = 0; i < nparticles; i++ )
       counts[i] = 0;
 
    for(  i = 0; i < nmaps; i++ ){
-      for(  j = 0; j < natoms; j++ ){
+      for(  j = 0; j < nparticles; j++ ){
          invmaps[i][j] = float4( -1.0, -1.0, -1.0, -1.0 );
       }
    }
@@ -1734,24 +1736,24 @@ int BrookBonded::gpuCalcInvMap_merged(
    //This will hold the number of imaps actually used
    *nimaps = -1;
 
-   //For each atom
+   //For each particle
    for(  i = 0; i < nints; i++ ){
       for(  j = 0; j < 4; j++ ){
          
-         atom = atoms[ i * 4 + j ];
+         particle = particles[ i * 4 + j ];
          
-         if ( atom == -1 ){
-         	//Nothing to be done for this atom, go to next
+         if ( particle == -1 ){
+         	//Nothing to be done for this particle, go to next
          	continue;
          }
          
          //Which map
-         mapnum = counts[ atom ] / 4;
+         mapnum = counts[ particle ] / 4;
          
          //Make sure we have space
          if ( mapnum >= nmaps ){
-         	printf( "Atom %d has too many bondeds(%d, max %d)\n",
-         			 atom, counts[atom], nmaps * 4 );
+         	printf( "Particle %d has too many bondeds(%d, max %d)\n",
+         			 particle, counts[particle], nmaps * 4 );
          	return 0;
          }
          	
@@ -1760,19 +1762,19 @@ int BrookBonded::gpuCalcInvMap_merged(
          }
 
          //Which component
-         mapcomp = counts[ atom ] % 4;
+         mapcomp = counts[ particle ] % 4;
          
          //Encode target stream and position
          pos = 100000 * j + i;
 
          switch ( mapcomp ){
-         	case 0: invmaps[mapnum][atom].x = (float) pos; break;
-         	case 1: invmaps[mapnum][atom].y = (float) pos; break;
-         	case 2: invmaps[mapnum][atom].z = (float) pos; break;
-         	case 3: invmaps[mapnum][atom].w = (float) pos; break;
+         	case 0: invmaps[mapnum][particle].x = (float) pos; break;
+         	case 1: invmaps[mapnum][particle].y = (float) pos; break;
+         	case 2: invmaps[mapnum][particle].z = (float) pos; break;
+         	case 3: invmaps[mapnum][particle].w = (float) pos; break;
          }
 
-         counts[ atom ]++;
+         counts[ particle ]++;
 
       }
    }
@@ -1784,28 +1786,282 @@ int BrookBonded::gpuCalcInvMap_merged(
 /* Repacks the invmap streams for more efficient access in the
  * merged inverse gather kernel
  *
- * buf should be nimaps * natoms large.
+ * buf should be nimaps * nparticles large.
  * */
-int BrookBonded::gpuRepackInvMap_merged( int natoms, int nmaps, int *counts, 
+int BrookBonded::gpuRepackInvMap_merged( int nparticles, int nmaps, int *counts, 
                                          float4 *invmaps[], float4 *buf ){
    int i, j;
    int nmaps_i;
 
-   for(  i = 0; i < natoms; i++ ){
+   for(  i = 0; i < nparticles; i++ ){
       for(  j = 0; j < nmaps; j++ ){
-         buf[ i + j*natoms ] = float4( -1.0f, -1.0f, -1.0f, -1.0f );
+         buf[ i + j*nparticles ] = float4( -1.0f, -1.0f, -1.0f, -1.0f );
       }
    }
    
-   for(  i = 0; i < natoms; i++ ){
+   for(  i = 0; i < nparticles; i++ ){
       
       nmaps_i = counts[i] / 4;
       if ( counts[i] % 4 ) 
          nmaps_i++;
       
       for(  j = 0; j < nmaps_i; j++ ){
-         buf[ i + j * natoms ] = invmaps[j][i];
+         buf[ i + j * nparticles ] = invmaps[j][i];
       }
    }
    return 1;
 }
+
+/** 
+ * Compute forces
+ * 
+ */
+
+void BrookBonded::computeForces( BrookStreamImpl& positionStream, BrookStreamImpl& forceStream ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookBonded::computeForces";
+
+   static const int I_Stream                = 0;
+   static const int J_Stream                = 1;
+   static const int K_Stream                = 2;
+   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 );
+
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   // bonded
+
+   float epsfac                                        = (float) (getLJ_14Scale()*getCoulombFactor());
+   float width                                         = (float) (getInverseMapStreamWidth());
+
+   // bonded forces
+
+   BrookFloatStreamInternal**  bondedParameters        = getBondedParameterStreams();
+   BrookFloatStreamInternal**  bondedForceStreams      = getBondedForceStreams();
+
+   BrookFloatStreamInternal**  inverseStreamMaps[4];
+   inverseStreamMaps[0]                                = getInverseStreamMapsStreams( 0 );
+   inverseStreamMaps[1]                                = getInverseStreamMapsStreams( 1 );
+   inverseStreamMaps[2]                                = getInverseStreamMapsStreams( 2 );
+   inverseStreamMaps[3]                                = getInverseStreamMapsStreams( 3 );
+
+   kbonded_CDLJ( epsfac, 
+                 (float) bondedForceStreams[0]->getStreamWidth(),
+                 dummyParameters,
+                 positionStream.getBrookStream(),
+                 getChargeStream()->getBrookStream(),
+                 getParticleIndicesStream()->getBrookStream(),
+                 bondedParameters[0]->getBrookStream(),
+                 bondedParameters[1]->getBrookStream(),
+                 bondedParameters[2]->getBrookStream(),
+                 bondedParameters[3]->getBrookStream(),
+                 bondedParameters[4]->getBrookStream(),
+                 bondedForceStreams[0]->getBrookStream(),
+                 bondedForceStreams[1]->getBrookStream(),
+                 bondedForceStreams[2]->getBrookStream(),
+                 bondedForceStreams[3]->getBrookStream() );
+
+
+   // diagnostics
+
+   if( 1 && PrintOn ){
+
+      int countPrintInvMap[4] = { 3, 5, 2, 4 }; 
+
+      (void) fprintf( getLog(), "\nPost kbonded_CDLJ: epsFac=%.6f %.6f %.6f", epsfac, getLJ_14Scale(), getCoulombFactor());
+      (void) fprintf( getLog(), "\nParticle indices stream\n" );
+      getParticleIndicesStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nCharge stream\n" );
+      getChargeStream()->printToFile( getLog() );
+
+      for( int ii = 0; ii < 5; ii++ ){
+         (void) fprintf( getLog(), "\nParam stream %d\n", ii );
+         bondedParameters[ii]->printToFile( getLog() );
+      }
+      for( int ii = 0; ii < 4; ii++ ){
+         (void) fprintf( getLog(), "\nForce stream %d\n", ii );
+         bondedForceStreams[ii]->printToFile( getLog() );
+      }
+
+/*
+      (void) fprintf( getLog(), "\nNB1 forces" );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      brookStreamInternalF->printToFile( getLog() );
+*/
+
+      (void) fprintf( getLog(), "\nInverse map streams -- K_Stream cnt=%d\n", getInverseMapStreamCount( K_Stream ) );
+      for( int ii = 0; ii < 4; ii++ ){
+         for( int jj = 0; jj < countPrintInvMap[ii]; jj++ ){
+            (void) fprintf( getLog(), "\n   Inverse map streams index=%d %d\n", ii, jj );
+            inverseStreamMaps[ii][jj]->printToFile( getLog() );
+         }
+      }
+   }
+
+   // gather forces
+
+   if( getInverseMapStreamCount( I_Stream ) == 3 && getInverseMapStreamCount( K_Stream ) == 3 ){
+
+      kinvmap_gather3_3( width,
+
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][2]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][2]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+
+   } else if( getInverseMapStreamCount( I_Stream ) == 3 && getInverseMapStreamCount( K_Stream ) == 4 ){
+
+      kinvmap_gather3_4( width,
+
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][2]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][2]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][3]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+
+   } else if( getInverseMapStreamCount( I_Stream ) == 3 && getInverseMapStreamCount( K_Stream ) == 5 ){
+
+      kinvmap_gather3_5( width,
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][2]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][2]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][3]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][4]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+
+   } else {
+
+      // case not handled -- throw an exception
+
+      if( getLog() && ErrorMessages++ < MaxErrorMessages ){
+         (void) fprintf( getLog(), "%s case: I-map=%d K-map=%d -- not handled.\n",
+                          methodName.c_str(), getInverseMapStreamCount( I_Stream ),
+                                              getInverseMapStreamCount( K_Stream ) );
+         (void) fflush(  getLog() );
+      }
+
+      kinvmap_gather3_3( width,
+
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][2]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][2]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+/*
+      std::stringstream message;
+      message << methodName << "I-maps=" << getInverseMapStreamCount( I_Stream ) << " and " << 
+                               "K-maps=" << getInverseMapStreamCount( K_Stream ) << " not handled.";
+      throw OpenMMException( message.str() );
+*/
+
+   }
+
+   // diagnostics
+
+   if( 0 && PrintOn ){
+
+      (void) fprintf( getLog(), "\nPost 3_4/3_5 && NB forces" );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      brookStreamInternalF->printToFile( getLog() );
+
+   }
+
+   if( getInverseMapStreamCount( J_Stream ) == 5 && getInverseMapStreamCount( L_Stream ) == 2 ){
+      kinvmap_gather5_2( width,
+                         inverseStreamMaps[J_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][2]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][3]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][4]->getBrookStream(),
+                         bondedForceStreams[J_Stream]->getBrookStream(),
+                         inverseStreamMaps[L_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[L_Stream][1]->getBrookStream(),
+                         bondedForceStreams[L_Stream]->getBrookStream(),
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+   
+   } else {
+
+      // case not handled -- throw an exception
+
+      if( getLog() && ErrorMessages++ < MaxErrorMessages ){
+         (void) fprintf( getLog(), "%s case: J-map=%d L-map=%d -- not handled.\n",
+                          methodName.c_str(), getInverseMapStreamCount( J_Stream ),
+                                              getInverseMapStreamCount( L_Stream ) );
+         (void) fflush(  getLog() );
+      }
+
+      kinvmap_gather5_2( width,
+                         inverseStreamMaps[J_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][2]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][3]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][4]->getBrookStream(),
+                         bondedForceStreams[J_Stream]->getBrookStream(),
+                         inverseStreamMaps[L_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[L_Stream][1]->getBrookStream(),
+                         bondedForceStreams[L_Stream]->getBrookStream(),
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+/*
+      std::stringstream message;
+      message << methodName << "J-maps=" << getInverseMapStreamCount( J_Stream ) << " and " << 
+                               "L-maps=" << getInverseMapStreamCount( L_Stream ) << " not handled.";
+      throw OpenMMException( message.str() );
+*/
+
+   }
+
+   // diagnostics
+
+   if( 1 && PrintOn ){
+
+      (void) fprintf( getLog(), "\nFinal NB & bonded forces" );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      brookStreamInternalF->printToFile( getLog() );
+/*
+      void* dataV = brookStreamInternalF->getData(1);
+      float* data = (float*) dataV;
+      (void) fprintf( getLog(), "\nFinal NB & bonded forces RAW\n" );
+      for( int ii = 0; ii < _brookNonBonded->getNumberOfParticles()*3; ii += 3 ){
+         (void) fprintf( getLog(), "%d [%.6e %.6e %.6e]\n", ii, data[ii], data[ii+1], data[ii+2] );
+      }
+*/
+
+   }
+
+   // ---------------------------------------------------------------------------------------
+}

platforms/brook/src/BrookBonded.h

@@ -34,8 +34,7 @@
 
 #include <vector>
 
-#include "BrookFloatStreamInternal.h"
-#include "BrookIntStreamInternal.h"
+#include "BrookStreamImpl.h"
 #include "BrookPlatform.h"
 #include "BrookCommon.h"
 #include "OpenMMContext.h"
@@ -68,17 +67,17 @@ class BrookBonded : public BrookCommon {
       /**
        * Initialize the kernel, setting up the values of all the force field parameters.
        * 
-       * @param bondIndices               the two atoms connected by each bond term
+       * @param bondIndices               the two particles connected by each bond term
        * @param bondParameters            the force parameters (length, k) for each bond term
-       * @param angleIndices              the three atoms connected by each angle term
+       * @param angleIndices              the three particles connected by each angle term
        * @param angleParameters           the force parameters (angle, k) for each angle term
-       * @param periodicTorsionIndices    the four atoms connected by each periodic torsion term
+       * @param periodicTorsionIndices    the four particles connected by each periodic torsion term
        * @param periodicTorsionParameters the force parameters (k, phase, periodicity) for each periodic torsion term
-       * @param rbTorsionIndices          the four atoms connected by each Ryckaert-Bellemans torsion term
+       * @param rbTorsionIndices          the four particles connected by each Ryckaert-Bellemans torsion term
        * @param rbTorsionParameters       the coefficients (in order of increasing powers) for each Ryckaert-Bellemans torsion term
-       * @param bonded14Indices           each element contains the indices of two atoms whose nonbonded interactions should be reduced since
+       * @param bonded14Indices           each element contains the indices of two particles whose nonbonded interactions should be reduced since
        *                                  they form a bonded 1-4 pair
-       * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+       * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each particle
        * @param lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
        * @param coulomb14Scale            the factor by which Coulomb interactions should be reduced for bonded 1-4 pairs
        * @param log                       log reference
@@ -87,7 +86,7 @@ class BrookBonded : public BrookCommon {
        *
        */
 
-      int setup( int numberOfAtoms, 
+      int setup( int numberOfParticles, 
                  const std::vector<std::vector<int> >& bondIndices,            const std::vector<std::vector<double> >& bondParameters,
                  const std::vector<std::vector<int> >& angleIndices,           const std::vector<std::vector<double> >& angleParameters,
                  const std::vector<std::vector<int> >& periodicTorsionIndices, const std::vector<std::vector<double> >& periodicTorsionParameters,
@@ -157,7 +156,7 @@ class BrookBonded : public BrookCommon {
        * @return 
        */
 
-      BrookFloatStreamInternal* getBrookAtomIndices( void ) const; 
+      BrookFloatStreamInternal* getBrookParticleIndices( void ) const; 
 
       /**
        * Get LJ 14 scale factor
@@ -178,13 +177,13 @@ class BrookBonded : public BrookCommon {
       BrookOpenMMFloat getCoulombFactor( void ) const;
 
       /** 
-       * Get bonded atom indices stream
+       * Get bonded particle indices stream
        * 
-       * @return  atom indices stream
+       * @return  particle indices stream
        *
        */
       
-      BrookFloatStreamInternal* getAtomIndicesStream( void ) const;
+      BrookFloatStreamInternal* getParticleIndicesStream( void ) const;
       
       /** 
        * Get bonded charge stream
@@ -268,6 +267,14 @@ class BrookBonded : public BrookCommon {
       
       std::string getContentsString( int level = 0 ) const;
 
+      /** 
+       * Compute forces
+       * 
+       */
+      
+      void computeForces( BrookStreamImpl& positionStream, BrookStreamImpl& forceStream );
+   
+   
    private:
    
       static const int NumberOfParameterStreams = 5;
@@ -291,7 +298,7 @@ class BrookBonded : public BrookCommon {
 
       // streams
 
-      BrookFloatStreamInternal* _atomIndicesStream;
+      BrookFloatStreamInternal* _particleIndicesStream;
       BrookFloatStreamInternal* _bondedParameters[NumberOfParameterStreams];
       BrookFloatStreamInternal* _bondedForceStreams[NumberOfForceStreams];
       BrookFloatStreamInternal* _chargeStream;
@@ -302,93 +309,93 @@ class BrookBonded : public BrookCommon {
 
       // helper methods in setup of parameters
 
-      void flipQuartet( int ibonded, int *atoms );
-      int matchTorsion( int i, int j, int k, int l, int nbondeds, int *atoms );
-      int matchAngle( int i, int j, int k, int nbondeds, int *atoms, int *flag );
-      int matchBond( int i, int j, int nbondeds, int *atoms, int *flag );
-      int matchPair( int i, int j, int nbondeds, int *atoms );
+      void flipQuartet( int ibonded, int *particles );
+      int matchTorsion( int i, int j, int k, int l, int nbondeds, int *particles );
+      int matchAngle( int i, int j, int k, int nbondeds, int *particles, int *flag );
+      int matchBond( int i, int j, int nbondeds, int *particles, int *flag );
+      int matchPair( int i, int j, int nbondeds, int *particles );
 
       /**
-       * Setup Ryckaert-Bellemans parameters/atom indices
+       * Setup Ryckaert-Bellemans parameters/particle indices
        * 
        * @param nbondeds                  number of bonded entries
-       * @param atoms                     array of atom indices
+       * @param particles                     array of particle indices
        * @param params                    arrays of bond parameters
-       * @param rbTorsionIndices          the four atoms connected by each Ryckaert-Bellemans torsion term
+       * @param rbTorsionIndices          the four particles connected by each Ryckaert-Bellemans torsion term
        * @param rbTorsionParameters       the coefficients (in order of increasing powers) for each Ryckaert-Bellemans torsion term
        *
        * @return nonzero value if error
        *
        */
 
-      int addRBDihedrals( int *nbondeds, int *atoms, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& rbTorsionIndices,  
+      int addRBDihedrals( int *nbondeds, int *particles, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& rbTorsionIndices,  
                           const std::vector<std::vector<double> >& rbTorsionParameters );
 
       /**
-       * Setup periodic torsion parameters/atom indices
+       * Setup periodic torsion parameters/particle indices
        * 
        * @param nbondeds                  number of bonded entries
-       * @param atoms                     array of atom indices
+       * @param particles                     array of particle indices
        * @param params                    arrays of bond parameters
-       * @param periodicTorsionIndices    the four atoms connected by each periodic torsion term
+       * @param periodicTorsionIndices    the four particles connected by each periodic torsion term
        * @param periodicTorsionParameters the force parameters (k, phase, periodicity) for each periodic torsion term
        *
        * @return nonzero value if error
        *
        */
 
-      int addPDihedrals( int *nbondeds, int *atoms, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& periodicTorsionIndices,  
+      int addPDihedrals( int *nbondeds, int *particles, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& periodicTorsionIndices,  
                           const std::vector<std::vector<double> >& periodicTorsionParameters );
 
       /**
-       * Setup angle bond parameters/atom indices
+       * Setup angle bond parameters/particle indices
        * 
        * @param nbondeds                  number of bonded entries
-       * @param atoms                     array of atom indices
+       * @param particles                     array of particle indices
        * @param params                    arrays of bond parameters
-       * @param angleIndices              the angle bond atom indices
+       * @param angleIndices              the angle bond particle indices
        * @param angleParameters           the angle parameters (angle in radians, force constant)
        *
        * @return nonzero value if error
        *
        */
 
-      int addAngles( int *nbondeds, int *atoms, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& angleIndices,
+      int addAngles( int *nbondeds, int *particles, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& angleIndices,
                      const std::vector<std::vector<double> >& angleParameters );
 
       /**
-       * Setup harmonic bond parameters/atom indices
+       * Setup harmonic bond parameters/particle indices
        * 
        * @param nbondeds                  number of bonded entries
-       * @param atoms                     array of atom indices
+       * @param particles                     array of particle indices
        * @param params                    arrays of bond parameters
-       * @param bondIndices               two harmonic bond atom indices
+       * @param bondIndices               two harmonic bond particle indices
        * @param bondParameters            the force parameters (distance, k)
        *
        * @return nonzero value if error
        *
        */
 
-      int addBonds( int *nbondeds, int *atoms, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& bondIndices,
+      int addBonds( int *nbondeds, int *particles, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& bondIndices,
                     const std::vector<std::vector<double> >& bondParameters );
 
       /**
-       * Setup LJ/Coulomb 1-4 parameters/atom indices
+       * Setup LJ/Coulomb 1-4 parameters/particle indices
        * 
        * @param nbondeds                  number of bonded entries
-       * @param atoms                     array of atom indices
+       * @param particles                     array of particle indices
        * @param params                    arrays of bond parameters
        * @param charges                   array of charges
-       * @param bonded14Indices           each element contains the indices of two atoms whose nonbonded interactions should be reduced since
+       * @param bonded14Indices           each element contains the indices of two particles whose nonbonded interactions should be reduced since
        *                                  they form a bonded 1-4 pair
-       * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+       * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each particle
        * @param lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
        *
        * @return nonzero value if error
        *
        */
 
-      int addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[], BrookOpenMMFloat* charges,
+      int addPairs( int *nbondeds, int *particles, BrookOpenMMFloat* params[], BrookOpenMMFloat* charges,
                     const std::vector<std::vector<int> >& bonded14Indices, const std::vector<std::vector<double> >& nonbondedParameters,
                     double lj14Scale, double coulombScale );
       
@@ -396,8 +403,8 @@ class BrookBonded : public BrookCommon {
        * Create and load inverse maps for bonded ixns
        * 
        * @param nbondeds                  number of bonded entries
-       * @param natoms                    number of atoms
-       * @param atoms                     arrays of atom indices (atoms[numberOfBonds][4])
+       * @param nparticles                    number of particles
+       * @param particles                     arrays of particle indices (particles[numberOfBonds][4])
        * @param platform                  BrookPlatform reference
        * @param log                       log file reference (optional)
        *
@@ -405,7 +412,7 @@ class BrookBonded : public BrookCommon {
        *
        */
 
-      int loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookPlatform& platform );
+      int loadInvMaps( int nbondeds, int nparticles, int *particles, const BrookPlatform& platform );
       
       /**
        * Validate inverse map count
@@ -425,23 +432,23 @@ class BrookBonded : public BrookCommon {
        * Helper functions for building inverse maps for 
        * torsions, impropers and angles.
        * 
-       * For each atom, calculates the positions at which it's
+       * For each particle, calculates the positions at which it's
        * forces are to be picked up from and stores the position
        * in the appropriate index.
        *
-       * Input: number of dihedrals, the atom indices, and a flag indicating
+       * Input: number of dihedrals, the particle indices, and a flag indicating
        *        whether we're doing i(0), j(1), k(2) or l(3)
-       * Output: an array of counts per atom
+       * Output: an array of counts per particle
        *         arrays of inversemaps
        *         nimaps - the number of invmaps actually used.
        *
-       * @param posflag       0-niatoms-1
-       * @param niatoms       3 for angles, 4 for torsions, impropers
+       * @param posflag       0-niparticles-1
+       * @param niparticles       3 for angles, 4 for torsions, impropers
        * @param nints         number of interactions
-       * @param natoms        number of atoms
-       * @param *atoms        gromacs interaction list
+       * @param nparticles        number of particles
+       * @param *particles        gromacs interaction list
        * @param nmaps         maximum number of inverse maps
-       * @param   counts[]    output counts of how many places each atom occurs
+       * @param   counts[]    output counts of how many places each particle occurs
        * @param *invmaps[]    output array of nmaps inverse maps
        * @param *nimaps,      output max number of inverse maps actually used
        *
@@ -449,11 +456,11 @@ class BrookBonded : public BrookCommon {
        *
        **/
       
-      int gpuCalcInvMap( int posflag, int niatoms, int nints, int natoms,
-                          int *atoms, int nmaps, int counts[], float4 *invmaps[],
+      int gpuCalcInvMap( int posflag, int niparticles, int nints, int nparticles,
+                          int *particles, int nmaps, int counts[], float4 *invmaps[],
                           int *nimaps );
       
-      void gpuPrintInvMaps( int nmaps, int natoms, int counts[], float4 *invmap[], FILE* logFile );
+      void gpuPrintInvMaps( int nmaps, int nparticles, int counts[], float4 *invmap[], FILE* logFile );
       
       /* We are still plagued by kernel call overheads. This is for a big fat
        * merged inverse gather kernel:
@@ -463,14 +470,14 @@ class BrookBonded : public BrookCommon {
        * lookup. This assumes that nints < 100000, preferably nints << 100000
        * which should always be true
        * */
-      int gpuCalcInvMap_merged( int nints, int natoms, int *atoms, int nmaps, int counts[], float4 *invmaps[], int *nimaps );
+      int gpuCalcInvMap_merged( int nints, int nparticles, int *particles, int nmaps, int counts[], float4 *invmaps[], int *nimaps );
       
       /* Repacks the invmap streams for more efficient access in the
        * merged inverse gather kernel
        *
-       * buf should be nimaps * natoms large.
+       * buf should be nimaps * nparticles large.
        * */
-      int gpuRepackInvMap_merged( int natoms, int nmaps, int *counts, float4 *invmaps[], float4 *buf );
+      int gpuRepackInvMap_merged( int nparticles, int nmaps, int *counts, float4 *invmaps[], float4 *buf );
       
         
 };

platforms/brook/src/BrookBrownianDynamics.cpp

@@ -64,13 +64,13 @@ BrookBrownianDynamics::BrookBrownianDynamics( ){
 
 // ---------------------------------------------------------------------------------------
 
-   _numberOfAtoms             = -1;
+   _numberOfParticles             = -1;
 
    // mark stream dimension variables as unset
 
-   _bdAtomStreamWidth         = -1;
-   _bdAtomStreamHeight        = -1;
-   _bdAtomStreamSize          = -1;
+   _bdParticleStreamWidth         = -1;
+   _bdParticleStreamHeight        = -1;
+   _bdParticleStreamSize          = -1;
 
    for( int ii = 0; ii < LastStreamIndex; ii++ ){
       _streams[ii]   = NULL;
@@ -376,9 +376,9 @@ int BrookBrownianDynamics::updateParameters( double temperature, double friction
 /** 
  * Update
  * 
- * @param  positions                   atom positions
- * @param  velocities                  atom velocities
- * @param  forces                      atom forces
+ * @param  positions                   particle positions
+ * @param  velocities                  particle velocities
+ * @param  forces                      particle forces
  * @param  brookShakeAlgorithm         BrookShakeAlgorithm reference
  * @param  brookRandomNumberGenerator  BrookRandomNumberGenerator reference
  *
@@ -413,7 +413,7 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
 
       static int showAux = 1;
 
-      if( PrintOn ){
+      if( showAux ){
          (void) fprintf( getLog(), "%s shake=%d\n", methodName, brookShakeAlgorithm.getNumberOfConstraints() );
          (void) fflush( getLog() );
       }
@@ -434,10 +434,51 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
       }
    }
 
+   // diagnostics
+
+   if( (1 || PrintOn) ){
+      (void) fprintf( getLog(), "\nPre kintegrate_bd: %d rngStrW=%3d rngOff=%5d "
+                                "ForceScale=%12.5e NoiseAmplitude=%12.5e\n",
+                                getBrownianDynamicsParticleStreamWidth(),
+                                brookRandomNumberGenerator.getRandomNumberStreamWidth(),
+                                brookRandomNumberGenerator.getRvStreamOffset(),
+                                getForceScale(),  getNoiseAmplitude() );
+
+      // (void) fprintf( getLog(), "\nInverseMassStream\n" );
+      //getInverseMassStream()->printToFile( getLog() );
+
+      //StreamImpl& positionStreamImpl               = positionStream.getImpl();
+      //const BrookStreamImpl brookPositions         = dynamic_cast<BrookStreamImpl&> (positionStreamImpl);
+/*
+      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
+      (void) fprintf( getLog(), "\nPositionStream\n" );
+      brookStreamInternalPos->printToFile( getLog() );
+*/
+
+      double forceSum[3];
+      BrookStreamInternal* brookStreamInternalFF       = forceStream.getBrookStreamImpl();
+      BrookFloatStreamInternal* brookStreamInternalF   = dynamic_cast<BrookFloatStreamInternal*> (brookStreamInternalFF);
+      brookStreamInternalF->sumByDimension( getNumberOfParticles(), forceSum );
+      (void) fprintf( getLog(), "\nForceStream [%18.10e %18.10e %18.10e]\n", forceSum[0], forceSum[1], forceSum[2] );
+      brookStreamInternalF->printToFile( getLog() );
+/*
+
+      (void) fprintf( getLog(), "\nXPrimeStream\n" );
+      getXPrimeStream()->printToFile( getLog() ); 
+
+      (void) fprintf( getLog(), "\nRvStreamIndex=%d\n", brookRandomNumberGenerator.getRvStreamIndex() );
+*/
+      // brookRandomNumberGenerator.getRandomNumberStream( brookRandomNumberGenerator.getRvStreamIndex() )->printToFile( getLog() ); 
+
+         BrookStreamInternal* brookStreamInternalVel  = velocityStream.getBrookStreamImpl();
+         (void) fprintf( getLog(), "\nVelocityStream\n" );
+         brookStreamInternalVel->printToFile( getLog() ); 
+   }   
+
    // integration step -- deltas returned in XPrime
 
    kintegrate_bd(
-         (float) getBrownianDynamicsAtomStreamWidth(),
+         (float) getBrownianDynamicsParticleStreamWidth(),
          (float) brookRandomNumberGenerator.getRandomNumberStreamWidth(),
          (float) brookRandomNumberGenerator.getRvStreamOffset(),
          getForceScale(), getNoiseAmplitude(),
@@ -450,7 +491,7 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
    if( PrintOn ){
       (void) fprintf( getLog(), "\nPost kintegrate_bd: %d rngStrW=%3d rngOff=%5d "
                                 "ForceScale=%12.5e NoiseAmplitude=%12.5e\n",
-                                getBrownianDynamicsAtomStreamWidth(),
+                                getBrownianDynamicsParticleStreamWidth(),
                                 brookRandomNumberGenerator.getRandomNumberStreamWidth(),
                                 brookRandomNumberGenerator.getRvStreamOffset(),
                                 getForceScale(),  getNoiseAmplitude() );
@@ -477,20 +518,20 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
 
    // advance random number cursor
 
-   brookRandomNumberGenerator.advanceGVCursor( getNumberOfAtoms() );
+   brookRandomNumberGenerator.advanceGVCursor( getNumberOfParticles() );
 
    // Shake
 
    if( brookShakeAlgorithm.getNumberOfConstraints() > 0 ){
       kshakeh_fix1( 
                     10.0f,
-                    (float) getBrownianDynamicsAtomStreamWidth(),
+                    (float) getBrownianDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getInverseHydrogenMass(),
                     omega, 
-                    brookShakeAlgorithm.getShakeAtomIndicesStream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeParticleIndicesStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),
-                    brookShakeAlgorithm.getShakeAtomParameterStream()->getBrookStream(),
+                    brookShakeAlgorithm.getShakeParticleParameterStream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons0Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons1Stream()->getBrookStream(),
                     brookShakeAlgorithm.getShakeXCons2Stream()->getBrookStream(),
@@ -499,7 +540,7 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
       // second Shake gather
    
       kshakeh_update2_fix1( 
-                    (float) getBrownianDynamicsAtomStreamWidth(),
+                    (float) getBrownianDynamicsParticleStreamWidth(),
                     brookShakeAlgorithm.getShakeInverseMapStream()->getBrookStream(),
                     positionStream.getBrookStream(),
                     getXPrimeStream()->getBrookStream(),
@@ -539,10 +580,10 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
 
       // diagnostics
    
-      if( PrintOn ){
+      if( (1 || PrintOn) ){
    
          (void) fprintf( getLog(), "\nPost kupdate_bd2: velocityScale=%12.5e\n", velocityScale );
-   
+
          (void) fprintf( getLog(), "\nXPrimeStream\n" );
          getXPrimeStream()->printToFile( getLog() ); 
 
@@ -563,35 +604,35 @@ int BrookBrownianDynamics::update( Stream& positions, Stream& velocities,
 };
 
 /** 
- * Get Atom stream size
+ * Get Particle stream size
  *
- * @return  Atom stream size
+ * @return  Particle stream size
  *
  */
 
-int BrookBrownianDynamics::getBrownianDynamicsAtomStreamSize( void ) const {
-   return _bdAtomStreamSize;
+int BrookBrownianDynamics::getBrownianDynamicsParticleStreamSize( void ) const {
+   return _bdParticleStreamSize;
 }
 
 /** 
- * Get atom stream width
+ * Get particle stream width
  *
- * @return  atom stream width
+ * @return  particle stream width
  *
  */
 
-int BrookBrownianDynamics::getBrownianDynamicsAtomStreamWidth( void ) const {
-   return _bdAtomStreamWidth;
+int BrookBrownianDynamics::getBrownianDynamicsParticleStreamWidth( void ) const {
+   return _bdParticleStreamWidth;
 }
 
 /** 
- * Get atom stream height
+ * Get particle stream height
  *
- * @return atom stream height
+ * @return particle stream height
  */
 
-int BrookBrownianDynamics::getBrownianDynamicsAtomStreamHeight( void ) const {
-   return _bdAtomStreamHeight;
+int BrookBrownianDynamics::getBrownianDynamicsParticleStreamHeight( void ) const {
+   return _bdParticleStreamHeight;
 }
 
 /** 
@@ -619,14 +660,14 @@ BrookFloatStreamInternal* BrookBrownianDynamics::getInverseMassStream( void ) co
 /** 
  * Initialize stream dimensions
  * 
- * @param numberOfAtoms             number of atoms
+ * @param numberOfParticles             number of particles
  * @param platform                  platform
  *      
  * @return ErrorReturnValue if error, else DefaultReturnValue
  *
  */
 
-int BrookBrownianDynamics::_initializeStreamSizes( int numberOfAtoms, const Platform& platform ){
+int BrookBrownianDynamics::_initializeStreamSizes( int numberOfParticles, const Platform& platform ){
 
 // ---------------------------------------------------------------------------------------
 
@@ -634,9 +675,9 @@ int BrookBrownianDynamics::_initializeStreamSizes( int numberOfAtoms, const Plat
 
 // ---------------------------------------------------------------------------------------
 
-   _bdAtomStreamSize     = getAtomStreamSize( platform );
-   _bdAtomStreamWidth    = getAtomStreamWidth( platform );
-   _bdAtomStreamHeight   = getAtomStreamHeight( platform );
+   _bdParticleStreamSize     = getParticleStreamSize( platform );
+   _bdParticleStreamWidth    = getParticleStreamWidth( platform );
+   _bdParticleStreamHeight   = getParticleStreamHeight( platform );
 
    return DefaultReturnValue;
 }
@@ -660,19 +701,19 @@ int BrookBrownianDynamics::_initializeStreams( const Platform& platform ){
 
 // ---------------------------------------------------------------------------------------
 
-   int sdAtomStreamSize             = getBrownianDynamicsAtomStreamSize();
-   int sdAtomStreamWidth            = getBrownianDynamicsAtomStreamWidth();
+   int sdParticleStreamSize             = getBrownianDynamicsParticleStreamSize();
+   int sdParticleStreamWidth            = getBrownianDynamicsParticleStreamWidth();
 
     _streams[VPrimeStream]          = new BrookFloatStreamInternal( BrookCommon::VPrimeStream,
-                                                                    sdAtomStreamSize, sdAtomStreamWidth,
+                                                                    sdParticleStreamSize, sdParticleStreamWidth,
                                                                     BrookStreamInternal::Float3, dangleValue );
 
     _streams[XPrimeStream]          = new BrookFloatStreamInternal( BrookCommon::XPrimeStream,
-                                                                    sdAtomStreamSize, sdAtomStreamWidth,
+                                                                    sdParticleStreamSize, sdParticleStreamWidth,
                                                                     BrookStreamInternal::Float3, dangleValue );
 
     _streams[InverseMassStream]     = new BrookFloatStreamInternal( BrookCommon::InverseMassStream,
-                                                                    sdAtomStreamSize, sdAtomStreamWidth,
+                                                                    sdParticleStreamSize, sdParticleStreamWidth,
                                                                     BrookStreamInternal::Float, dangleValue );
 
    return DefaultReturnValue;
@@ -693,7 +734,7 @@ int BrookBrownianDynamics::_updateStreams( void ){
 
 // ---------------------------------------------------------------------------------------
 
-   //int atomStreamSize                        = getBrownianDynamicsAtomStreamSize();
+   //int particleStreamSize                        = getBrownianDynamicsParticleStreamSize();
 
    return DefaultReturnValue;
 
@@ -702,7 +743,7 @@ int BrookBrownianDynamics::_updateStreams( void ){
 /** 
  * Set masses 
  * 
- * @param masses             atomic masses
+ * @param masses             particle masses
  *
  */
 
@@ -754,12 +795,12 @@ int BrookBrownianDynamics::setup( const std::vector<double>& masses, const Platf
    const BrookPlatform brookPlatform            = dynamic_cast<const BrookPlatform&> (platform);
    setLog( brookPlatform.getLog() );
 
-   int numberOfAtoms  = (int) masses.size();
-   setNumberOfAtoms( numberOfAtoms );
+   int numberOfParticles  = (int) masses.size();
+   setNumberOfParticles( numberOfParticles );
 
    // set stream sizes and then create streams
 
-   _initializeStreamSizes( numberOfAtoms, platform );
+   _initializeStreamSizes( numberOfParticles, platform );
    _initializeStreams( platform );
 
    _setInverseSqrtMasses( masses );
@@ -798,17 +839,17 @@ std::string BrookBrownianDynamics::getContentsString( int level ) const {
 #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", getNumberOfParticles() );
+   message << _getLine( tab, "Number of particles:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamWidth() );
-   message << _getLine( tab, "Atom stream width:", value ); 
+   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamWidth() );
+   message << _getLine( tab, "Particle stream width:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamHeight() );
-   message << _getLine( tab, "Atom stream height:", value ); 
+   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamHeight() );
+   message << _getLine( tab, "Particle stream height:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamSize() );
-   message << _getLine( tab, "Atom stream size:", value ); 
+   (void) LOCAL_SPRINTF( value, "%d", getParticleStreamSize() );
+   message << _getLine( tab, "Particle stream size:", value ); 
 
    (void) LOCAL_SPRINTF( value, "%.5f", getTau() );
    message << _getLine( tab, "Tau:", value ); 
@@ -819,10 +860,10 @@ std::string BrookBrownianDynamics::getContentsString( int level ) const {
    (void) LOCAL_SPRINTF( value, "%.5f", getStepSize() );
    message << _getLine( tab, "Step size:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%.5f", getForceScale() );
+   (void) LOCAL_SPRINTF( value, "%.5e", getForceScale() );
    message << _getLine( tab, "Force scale:", value ); 
 
-   (void) LOCAL_SPRINTF( value, "%.5f", getNoiseAmplitude() );
+   (void) LOCAL_SPRINTF( value, "%.5e", getNoiseAmplitude() );
    message << _getLine( tab, "Noise amplitude:", value ); 
 
    (void) LOCAL_SPRINTF( value, "%.5f", getTemperature() );

platforms/brook/src/BrookBrownianDynamics.h

@@ -116,28 +116,28 @@ class BrookBrownianDynamics : public BrookCommon {
       BrookOpenMMFloat getForceScale( void ) const; 
 
       /**
-       * Get BrownianDynamics atom stream width
+       * Get BrownianDynamics particle stream width
        *
-       * @return atom stream width
+       * @return particle stream width
        */
 
-      int getBrownianDynamicsAtomStreamWidth( void ) const; 
+      int getBrownianDynamicsParticleStreamWidth( void ) const; 
 
       /**
-       * Get BrownianDynamics atom stream height
+       * Get BrownianDynamics particle stream height
        *
-       * @return atom stream height
+       * @return particle stream height
        */
 
-      int getBrownianDynamicsAtomStreamHeight( void ) const;
+      int getBrownianDynamicsParticleStreamHeight( void ) const;
 
       /**
-       * Get BrownianDynamics atom stream size
+       * Get BrownianDynamics particle stream size
        * 
-       * @return atom stream size
+       * @return particle stream size
        */
 
-      int getBrownianDynamicsAtomStreamSize( void ) const; 
+      int getBrownianDynamicsParticleStreamSize( void ) const; 
 
       /** 
        * Update parameters
@@ -155,9 +155,9 @@ class BrookBrownianDynamics : public BrookCommon {
       /** 
        * Update
        * 
-       * @param  positions                   atom positions
-       * @param  velocities                  atom velocities
-       * @param  forces                      atom forces
+       * @param  positions                   particle positions
+       * @param  velocities                  particle velocities
+       * @param  forces                      particle forces
        * @param  brookShakeAlgorithm         BrookShakeAlgorithm reference
        * @param  brookRandomNumberGenerator  BrookRandomNumberGenerator reference
        *
@@ -181,7 +181,7 @@ class BrookBrownianDynamics : public BrookCommon {
       /* 
        * Setup of BrownianDynamics parameters
        *
-       * @param masses                atom masses
+       * @param masses                particle masses
        * @param platform              Brook platform
        *
        * @return ErrorReturnValue value if error, else DefaultReturnValue
@@ -240,11 +240,11 @@ class BrookBrownianDynamics : public BrookCommon {
       BrookOpenMMFloat _forceScale;
       BrookOpenMMFloat _noiseAmplitude;
 
-      // Atom stream dimensions
+      // Particle stream dimensions
 
-      int _bdAtomStreamWidth;
-      int _bdAtomStreamHeight;
-      int _bdAtomStreamSize;
+      int _bdParticleStreamWidth;
+      int _bdParticleStreamHeight;
+      int _bdParticleStreamSize;
 
       /*
        * Update streams
@@ -328,26 +328,26 @@ class BrookBrownianDynamics : public BrookCommon {
       /* 
        * Setup of stream dimensions
        *
-       * @param atomStreamSize        atom stream size
-       * @param atomStreamWidth       atom stream width
+       * @param particleStreamSize        particle stream size
+       * @param particleStreamWidth       particle stream width
        *
        * @return ErrorReturnValue if error, else DefaultReturnValue
        *
        * */
       
-      int _initializeStreamSizes( int atomStreamSize, int atomStreamWidth );
+      int _initializeStreamSizes( int particleStreamSize, int particleStreamWidth );
 
       /** 
        * Initialize stream dimensions
        * 
-       * @param numberOfAtoms             number of atoms
+       * @param numberOfParticles             number of particles
        * @param platform                  platform
        *
        * @return ErrorReturnValue if error, else DefaultReturnValue
        *
        */
       
-      int _initializeStreamSizes( int numberOfAtoms, const Platform& platform );
+      int _initializeStreamSizes( int numberOfParticles, const Platform& platform );
       
       /** 
        * Initialize stream dimensions and streams
@@ -363,7 +363,7 @@ class BrookBrownianDynamics : public BrookCommon {
       /** 
        * Set masses 
        * 
-       * @param masses             atomic masses
+       * @param masses             particleic masses
        *
        */
       

platforms/brook/src/BrookCalcGBSAOBCForceFieldKernel.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: 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 <cmath>
-#include <limits>
-#include "OpenMMException.h"
-#include <sstream>
-
-#include "BrookStreamImpl.h"
-#include "BrookCalcGBSAOBCForceKernel.h"
-#include "gpu/kgbsa.h"
-#include "gpu/kforce.h"
-#include "math.h"
-
-using namespace OpenMM;
-using namespace std;
-
-/** 
- * BrookCalcGBSAOBCForceKernel constructor
- * 
- * @param name                      kernel name
- * @param platform                  platform
- *
- */
-
-BrookCalcGBSAOBCForceKernel::BrookCalcGBSAOBCForceKernel( std::string name, const Platform& platform ) :
-                     CalcGBSAOBCForceKernel( name, platform ){
-
-// ---------------------------------------------------------------------------------------
-
-   // static const std::string methodName      = "BrookCalcGBSAOBCForceKernel::BrookCalcGBSAOBCForceKernel";
-   // static const int debug                   = 1;
-
-// ---------------------------------------------------------------------------------------
-
-   _numberOfAtoms                    = 0;
-   _brookGbsa                        = NULL;
-   _log                              = NULL;
-
-   const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
-   if( brookPlatform.getLog() != NULL ){
-      setLog( brookPlatform.getLog() );
-   }
-      
-}   
-
-/** 
- * BrookCalcGBSAOBCForceKernel destructor
- * 
- */
-
-BrookCalcGBSAOBCForceKernel::~BrookCalcGBSAOBCForceKernel( ){
-
-// ---------------------------------------------------------------------------------------
-
-   // static const std::string methodName      = "BrookCalcGBSAOBCForceKernel::BrookCalcGBSAOBCForceKernel";
-   // static const int debug                   = 1;
-
-// ---------------------------------------------------------------------------------------
-
-   delete _brookGbsa;
-}
-
-/** 
- * Get log file reference
- * 
- * @return  log file reference
- *
- */
-
-FILE* BrookCalcGBSAOBCForceKernel::getLog( void ) const {
-   return _log;
-}
-
-/** 
- * Set log file reference
- * 
- * @param  log file reference
- *
- * @return  DefaultReturnValue
- *
- */
-
-int BrookCalcGBSAOBCForceKernel::setLog( FILE* log ){
-   _log = log;
-   return BrookCommon::DefaultReturnValue;
-}
-
-/** 
- * Initialize the kernel, setting up the values of all the force field parameters.
- * 
- * @param atomParameters            vector containing atom index, charge, radius, scalingFactor
- * @param solventDielectric         solvent dielectric
- * @param soluteDielectric          solute dielectric
- *
- */
-
-void BrookCalcGBSAOBCForceKernel::initialize( const std::vector<std::vector<double> >& atomParameters, 
-                                                   double solventDielectric, double soluteDielectric ){
-
-// ---------------------------------------------------------------------------------------
-
-   static const std::string methodName      = "BrookCalcGBSAOBCForceKernel::initialize";
-
-// ---------------------------------------------------------------------------------------
-
-   FILE* log                 = getLog();
-   _numberOfAtoms            = (int) atomParameters.size();
-
-   // ---------------------------------------------------------------------------------------
-
-   // bonded
-
-   if( _brookGbsa ){
-      delete _brookGbsa;
-   }
-   _brookGbsa              = new BrookGbsa();
-   _brookGbsa->setLog( log );
-    
-   _brookGbsa->setup( atomParameters, solventDielectric, soluteDielectric, getPlatform() );
-
-   if( log ){
-      std::string contents = _brookGbsa->getContentsString( ); 
-      (void) fprintf( log, "%s brookGbsa::contents\n%s", methodName.c_str(), contents.c_str() );
-      (void) fflush( log );
-   }
-
-   // ---------------------------------------------------------------------------------------
-    
-}
-
-/** 
- * Compute forces given atom coordinates
- * 
- * @param positions                 atom coordinates
- * @param forces                    output forces
- *
- */
-
-void BrookCalcGBSAOBCForceKernel::executeForces( const Stream& positions, Stream& forces ){
-
-// ---------------------------------------------------------------------------------------
-
-   static const std::string methodName   = "BrookCalcGBSAOBCForceKernel::executeForces";
-   static const int PrintOn              = 0; 
-
-// ---------------------------------------------------------------------------------------
-
-   // OBC
-
-   const BrookStreamImpl& positionStreamC              = dynamic_cast<const BrookStreamImpl&> (positions.getImpl());
-   BrookStreamImpl& positionStream                     = const_cast<BrookStreamImpl&>         (positionStreamC);
-   BrookStreamImpl& forceStream                        = dynamic_cast<BrookStreamImpl&>       (forces.getImpl());
-
-   float includeAce                                    = (float) (_brookGbsa->includeAce());
-   BrookFloatStreamInternal**  gbsaForceStreams        = _brookGbsa->getForceStreams();
-
-   // calculate Born radii first time thru and initialize on board
-
-   if( !_brookGbsa->haveBornRadiiBeenInitialized() ){
-      _brookGbsa->calculateBornRadii( positions );  
-   }
-
-   // first major loop
-
-   kObcLoop1( (float) _brookGbsa->getNumberOfAtoms(),
-              (float) _brookGbsa->getAtomSizeCeiling(),
-              (float) _brookGbsa->getDuplicationFactor(),
-              (float) _brookGbsa->getAtomStreamWidth( ),
-              (float) _brookGbsa->getPartialForceStreamWidth( ),
-              _brookGbsa->getSoluteDielectric(),
-              _brookGbsa->getSolventDielectric(),
-              includeAce,
-
-              positionStream.getBrookStream(),
-
-              _brookGbsa->getObcBornRadii()->getBrookStream(),
-              _brookGbsa->getObcAtomicRadii()->getBrookStream(),
-
-              gbsaForceStreams[0]->getBrookStream(),
-              gbsaForceStreams[1]->getBrookStream(),
-              gbsaForceStreams[2]->getBrookStream(),
-              gbsaForceStreams[3]->getBrookStream()
-            );
-
-// ---------------------------------------------------------------------------------------
-
-   // diagnostics
-
-   if( 1 && PrintOn && getLog() ){
-
-      (void) fprintf( getLog(), "\nPost kObcLoop1: atms=%d ceil=%d dup=%d atomStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
-                      _brookGbsa->getNumberOfAtoms(),
-                      _brookGbsa->getAtomSizeCeiling(),
-                      _brookGbsa->getDuplicationFactor(),
-                      _brookGbsa->getAtomStreamWidth( ),
-                      _brookGbsa->getPartialForceStreamWidth( ),
-                      _brookGbsa->getSoluteDielectric(),
-                      _brookGbsa->getSolventDielectric(), includeAce );
-
-      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
-      (void) fprintf( getLog(), "\nPositionStream\n" );
-      brookStreamInternalPos->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\nBornR\n" );
-      _brookGbsa->getObcBornRadii()->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\nAtomR\n" );
-      _brookGbsa->getObcAtomicRadii()->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\nForceStreams output\n" );
-      for( int ii = 0; ii < 4; ii++ ){
-         gbsaForceStreams[ii]->printToFile( getLog() );
-      }
-
-   }
-
-// ---------------------------------------------------------------------------------------
-
-   // gather for first loop
-
-   kPostObcLoop1_nobranch(
-              (float) _brookGbsa->getDuplicationFactor(),
-              (float) _brookGbsa->getAtomStreamWidth( ),
-              (float) _brookGbsa->getPartialForceStreamWidth( ),
-              (float) _brookGbsa->getNumberOfAtoms(),
-              (float) _brookGbsa->getAtomSizeCeiling(),
-              (float) _brookGbsa->getInnerLoopUnroll(),
-              gbsaForceStreams[0]->getBrookStream(),
-              gbsaForceStreams[1]->getBrookStream(),
-              gbsaForceStreams[2]->getBrookStream(),
-              gbsaForceStreams[3]->getBrookStream(),
-              _brookGbsa->getObcChain()->getBrookStream(),
-              _brookGbsa->getObcBornRadii()->getBrookStream(),
-              _brookGbsa->getObcIntermediateForce()->getBrookStream(),
-              _brookGbsa->getObcBornRadii2()->getBrookStream() );
- 
-// ---------------------------------------------------------------------------------------
-
-   // diagnostics
-
-   if( PrintOn && getLog()){
-
-      (void) fprintf( getLog(), "\nPost kPostObcLoop1_nobranch: dup=%d aStrW=%d pStrW=%d no.atms=%3d ceil=%3d Unroll=%1d\n",
-                      _brookGbsa->getDuplicationFactor(),
-                      _brookGbsa->getAtomStreamWidth( ),
-                      _brookGbsa->getPartialForceStreamWidth( ),
-                      _brookGbsa->getNumberOfAtoms(),
-                      _brookGbsa->getAtomSizeCeiling(),
-                      _brookGbsa->getInnerLoopUnroll() );
-
-      (void) fprintf( getLog(), "\nForceStreams\n" );
-      for( int ii = 0; ii < 4; ii++ ){
-         gbsaForceStreams[ii]->printToFile( getLog() );
-      }
-
-      (void) fprintf( getLog(), "\nObcChain\n" );
-      _brookGbsa->getObcChain()->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\nBornR\n" );
-      _brookGbsa->getObcBornRadii()->printToFile( getLog() );
-
-      // output
-
-      (void) fprintf( getLog(), "\nObcIntermediateForce output\n" );
-      _brookGbsa->getObcIntermediateForce()->printToFile( getLog() );
-
-      // output
-
-      (void) fprintf( getLog(), "\nObcBornRadii2 output\n" );
-      _brookGbsa->getObcBornRadii2()->printToFile( getLog() );
-
-   }
-
-// ---------------------------------------------------------------------------------------
-
-   // second major loop
-
-   kObcLoop2( (float) _brookGbsa->getNumberOfAtoms(),
-              (float) _brookGbsa->getAtomSizeCeiling(),
-              (float) _brookGbsa->getDuplicationFactor(),
-              (float) _brookGbsa->getAtomStreamWidth( ),
-              (float) _brookGbsa->getPartialForceStreamWidth( ),
-              positionStream.getBrookStream(),
-              _brookGbsa->getObcScaledAtomicRadii()->getBrookStream(),
-              _brookGbsa->getObcBornRadii2()->getBrookStream(),
-              gbsaForceStreams[0]->getBrookStream(),
-              gbsaForceStreams[1]->getBrookStream(),
-              gbsaForceStreams[2]->getBrookStream(),
-              gbsaForceStreams[3]->getBrookStream()
-            );
-
-// ---------------------------------------------------------------------------------------
-
-   // diagnostics
-
-   if( PrintOn && getLog() ){
-
-      (void) fprintf( getLog(), "\nPost kObcLoop2: no.atms=%5d ceil=%3d dup=%3d strW=%3d pStrW=%3d\n",
-                      _brookGbsa->getNumberOfAtoms(),
-                      _brookGbsa->getAtomSizeCeiling(),
-                      _brookGbsa->getDuplicationFactor(),
-                      _brookGbsa->getAtomStreamWidth( ),
-                      _brookGbsa->getPartialForceStreamWidth( ) );
-
-      BrookStreamInternal* brookStreamInternalPos  = positionStream.getBrookStreamImpl();
-      (void) fprintf( getLog(), "\nPositionStream\n" );
-      brookStreamInternalPos->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\nObcScaledAtomicRadii\n" );
-      _brookGbsa->getObcScaledAtomicRadii()->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\ngetObcBornRadii2\n" );
-      _brookGbsa->getObcBornRadii2()->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\nForceStreams\n" );
-      for( int ii = 0; ii < 4; ii++ ){
-         gbsaForceStreams[ii]->printToFile( getLog() );
-      }
-
-   }
-
-// ---------------------------------------------------------------------------------------
-
-   // gather for second loop
-
-    float mergeNonObcForces = 1.0f;
-    float kcalMolTokJNM     = -0.4184f;
-    kPostObcLoop2_nobranch(
-              (float) _brookGbsa->getDuplicationFactor(),
-              (float) _brookGbsa->getAtomStreamWidth( ),
-              (float) _brookGbsa->getPartialForceStreamWidth( ),
-              (float) _brookGbsa->getNumberOfAtoms(),
-              (float) _brookGbsa->getAtomSizeCeiling(),
-              (float) _brookGbsa->getInnerLoopUnroll(),
-              kcalMolTokJNM,
-              mergeNonObcForces,
-              _brookGbsa->getObcIntermediateForce()->getBrookStream(),
-              forceStream.getBrookStream(),
-              gbsaForceStreams[0]->getBrookStream(),
-              gbsaForceStreams[1]->getBrookStream(),
-              gbsaForceStreams[2]->getBrookStream(),
-              gbsaForceStreams[3]->getBrookStream(),
-              _brookGbsa->getObcAtomicRadii()->getBrookStream(),
-              _brookGbsa->getObcBornRadii()->getBrookStream(),
-              _brookGbsa->getObcChain()->getBrookStream(),
-              forceStream.getBrookStream()
-           );
-
-// ---------------------------------------------------------------------------------------
-
-   // diagnostics
-
-   if( PrintOn && getLog() ){
-
-      (void) fprintf( getLog(), "\nPost kPostObcLoop2_nobranch: atms=%d ceil=%d dup=%d atomStrW=%3d prtlF=%3d diel=%.3f %.3f ACE=%.1f\n",
-                      _brookGbsa->getNumberOfAtoms(),
-                      _brookGbsa->getAtomSizeCeiling(),
-                      _brookGbsa->getDuplicationFactor(),
-                      _brookGbsa->getAtomStreamWidth( ),
-                      _brookGbsa->getPartialForceStreamWidth( ),
-                      _brookGbsa->getSoluteDielectric(),
-                      _brookGbsa->getSolventDielectric(), includeAce );
-
-      (void) fprintf( getLog(), "\nPartialForceStreams\n" );
-      for( int ii = 0; ii < 4; ii++ ){
-         gbsaForceStreams[ii]->printToFile( getLog() );
-      }
-
-      BrookStreamInternal* brookStreamInternalF  = forceStream.getBrookStreamImpl();
-      (void) fprintf( getLog(), "\nForceStream\n" );
-      brookStreamInternalF->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\nChain\n" );
-      _brookGbsa->getObcChain()->printToFile( getLog() );
-
-      (void) fprintf( getLog(), "\nBornR\n" );
-      _brookGbsa->getObcBornRadii()->printToFile( getLog() );
-
-   }
-
-   // ---------------------------------------------------------------------------------------
-}
-
-/**
- * Execute the kernel to calculate the OBC energy
- * 
- * @param positions   atom positions
- *
- * @return  potential energy due to the OBC forces
- *
- */
-
-double BrookCalcGBSAOBCForceKernel::executeEnergy( const Stream& positions ){
-
-// ---------------------------------------------------------------------------------------
-
-   //static const std::string methodName      = "BrookCalcGBSAOBCForceKernel::executeEnergy";
-
-// ---------------------------------------------------------------------------------------
-
-   return (double) _brookGbsa->getEnergy( positions );
-
-}

platforms/brook/src/BrookNonbonded14ForceKernel.cpp → platforms/brook/src/BrookCalcGBSAOBCForceKernel.cpp

@@ -29,17 +29,22 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include <cmath>
+#include <limits>
 #include "OpenMMException.h"
 #include <sstream>
-#include "BrookNonbonded14ForceKernel.h"
+
+#include "BrookStreamImpl.h"
+#include "BrookCalcGBSAOBCForceKernel.h"
+#include "gpu/kgbsa.h"
+#include "gpu/kforce.h"
+#include "math.h"
 
 using namespace OpenMM;
 using namespace std;
 
-const std::string BrookNonbonded14ForceKernel::BondName = "HarmonicLJ14";
-
 /** 
- * BrookNonbonded14ForceKernel constructor
+ * BrookCalcGBSAOBCForceKernel constructor
  * 
  * @param name                      kernel name
  * @param platform                  platform
@@ -48,18 +53,18 @@ const std::string BrookNonbonded14ForceKernel::BondName = "HarmonicLJ14";
  *
  */
 
-BrookNonbonded14ForceKernel::BrookNonbonded14ForceKernel( std::string name, const Platform& platform,
-                                                          OpenMMBrookInterface& openMMBrookInterface, System& system ) :
-                     CalcHarmonicLJ14ForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
+BrookCalcGBSAOBCForceKernel::BrookCalcGBSAOBCForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system ) :
+                             CalcGBSAOBCForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookNonbonded14ForceKernel::BrookNonbonded14ForceKernel";
+   // static const std::string methodName      = "BrookCalcGBSAOBCForceKernel::BrookCalcGBSAOBCForceKernel";
    // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
-   _brookBondParameters              = NULL;
+   _numberOfParticles                = 0;
+   _brookGbsa                        = NULL;
    _log                              = NULL;
 
    const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
@@ -70,20 +75,20 @@ BrookNonbonded14ForceKernel::BrookNonbonded14ForceKernel( std::string name, cons
 }   
 
 /** 
- * BrookNonbonded14ForceKernel destructor
+ * BrookCalcGBSAOBCForceKernel destructor
  * 
  */
 
-BrookNonbonded14ForceKernel::~BrookNonbonded14ForceKernel( ){
+BrookCalcGBSAOBCForceKernel::~BrookCalcGBSAOBCForceKernel( ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookNonbonded14ForceKernel::BrookNonbonded14ForceKernel";
+   // static const std::string methodName      = "BrookCalcGBSAOBCForceKernel::BrookCalcGBSAOBCForceKernel";
    // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
-   delete _brookBondParameters;
+   delete _brookGbsa;
 }
 
 /** 
@@ -93,7 +98,7 @@ BrookNonbonded14ForceKernel::~BrookNonbonded14ForceKernel( ){
  *
  */
 
-FILE* BrookNonbonded14ForceKernel::getLog( void ) const {
+FILE* BrookCalcGBSAOBCForceKernel::getLog( void ) const {
    return _log;
 }
 
@@ -106,7 +111,7 @@ FILE* BrookNonbonded14ForceKernel::getLog( void ) const {
  *
  */
 
-int BrookNonbonded14ForceKernel::setLog( FILE* log ){
+int BrookCalcGBSAOBCForceKernel::setLog( FILE* log ){
    _log = log;
    return BrookCommon::DefaultReturnValue;
 }
@@ -114,16 +119,16 @@ int BrookNonbonded14ForceKernel::setLog( FILE* log ){
 /** 
  * Initialize the kernel, setting up the values of all the force field parameters.
  * 
- * @param system                    System reference
- * @param force                     HarmonicLJ14Force reference
+ * @param system     system this kernel will be applied to
+ * @param force      GBSAOBCForce this kernel will be used for
  *
  */
 
-void BrookNonbonded14ForceKernel::initialize( const System& system, const HarmonicLJ14Force& force ){
+void BrookCalcGBSAOBCForceKernel::initialize( const System& system, const GBSAOBCForce& force ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookNonbonded14ForceKernel::initialize";
+   static const std::string methodName      = "BrookCalcGBSAOBCForceKernel::initialize";
 
 // ---------------------------------------------------------------------------------------
 
@@ -131,39 +136,36 @@ void BrookNonbonded14ForceKernel::initialize( const System& system, const Harmon
 
    // ---------------------------------------------------------------------------------------
 
-   // create _brookBondParameters object containing atom indices/parameters
-
-   int numberOfBonds         = force.getNumLJ14s();
-
-   if( _brookBondParameters ){
-      delete _brookBondParameters;
+   if( _brookGbsa ){
+      delete _brookGbsa;
    }
-   _brookBondParameters = new BrookBondParameters( BondName, NumberOfAtomsInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
-
-   for( int ii = 0; ii < numberOfBonds; ii++ ){
-
-      int particle1, particle2, particle3;
-      double angle, k;
-
-      int particles[NumberOfAtomsInBond];
-      double parameters[NumberOfParametersInBond];
-
-      force.getLJ14Parameters( ii, particle1, particle2, particle3, angle, k ); 
-      particles[0]    = particle1;
-      particles[1]    = particle2;
-      particles[2]    = particle3;
- 
-      parameters[0]   = angle;
-      parameters[1]   = k;
-
-      _brookBondParameters->setBond( ii, particles, parameters );
+   _brookGbsa              = new BrookGbsa();
+   _brookGbsa->setLog( log );
+    
+   // get parameters from force object
+   // and initialize brookGbsa
+
+   _numberOfParticles = system.getNumParticles();
+   std::vector<std::vector<double> > particleParameters;
+   for( int ii = 0; ii < _numberOfParticles; ii++ ){
+
+      double charge, radius, scalingFactor;
+      force.getParticleParameters( ii, charge, radius, scalingFactor );
+
+      std::vector<double> parameters;
+      particleParameters[ii] = parameters;
+       
+      parameters[0]      = charge;
+      parameters[1]      = radius;
+      parameters[2]      = scalingFactor;
    }   
-   _openMMBrookInterface.setHarmonicLJ14ForceParameters( _brookBondParameters );
+   _brookGbsa->setup( particleParameters, force.getSolventDielectric(), force.getSoluteDielectric(), getPlatform() );
+
    _openMMBrookInterface.setTriggerForceKernel( this );
    _openMMBrookInterface.setTriggerEnergyKernel( this );
 
    if( log ){
-      std::string contents = _brookBondParameters->getContentsString( ); 
+      std::string contents = _brookGbsa->getContentsString( ); 
       (void) fprintf( log, "%s brookGbsa::contents\n%s", methodName.c_str(), contents.c_str() );
       (void) fflush( log );
    }
@@ -173,43 +175,40 @@ void BrookNonbonded14ForceKernel::initialize( const System& system, const Harmon
 }
 
 /** 
- * Compute forces given atom coordinates
+ * Compute forces given particle coordinates
  * 
  * @param context OpenMMContextImpl context
  *
  */
 
-void BrookNonbonded14ForceKernel::executeForces( OpenMMContextImpl& context ){
+void BrookCalcGBSAOBCForceKernel::executeForces( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName   = "BrookNonbonded14ForceKernel::executeForces";
+   //static const std::string methodName   = "BrookCalcGBSAOBCForceKernel::executeForces";
 
 // ---------------------------------------------------------------------------------------
 
    if( _openMMBrookInterface.getTriggerForceKernel() == this ){
       _openMMBrookInterface.computeForces( context );
-   }
-
-   return;
+   }   
 
-   // ---------------------------------------------------------------------------------------
 }
 
 /**
- * Execute the kernel to calculate the energy
+ * Execute the kernel to calculate the OBC energy
  * 
  * @param context OpenMMContextImpl context
  *
- * @return  potential energy
+ * @return energy
  *
  */
 
-double BrookNonbonded14ForceKernel::executeEnergy( OpenMMContextImpl& context ){
+double BrookCalcGBSAOBCForceKernel::executeEnergy( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName      = "BrookNonbonded14ForceKernel::executeEnergy";
+   //static const std::string methodName      = "BrookCalcGBSAOBCForceKernel::executeEnergy";
 
 // ---------------------------------------------------------------------------------------
 

platforms/brook/src/BrookCalcGBSAOBCForceFieldKernel.h → platforms/brook/src/BrookCalcGBSAOBCForceKernel.h

-#ifndef OPENMM_BROOK__CALCL_GBSAOBC_FORCEFIELD_KERNEL_H_
-#define OPENMM_BROOK__CALCL_GBSAOBC_FORCEFIELD_KERNEL_H_
+#ifndef OPENMM_BROOK_CALC_GBSAOBC_FORCE_KERNEL_H_
+#define OPENMM_BROOK_CALC_GBSAOBC_FORCE_KERNEL_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -35,12 +35,14 @@
 #include "kernels.h"
 #include "../../reference/src/SimTKUtilities/SimTKOpenMMRealType.h"
 #include "BrookGbsa.h"
+#include "OpenMMBrookInterface.h"
 
 namespace OpenMM {
 
 /**
- * This kernel is invoked by NonbondedForce to calculate the forces acting on the system.
+ * This kernel is invoked to calculate the OBC forces acting on the system.
  */
+
 class BrookCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
 
    public:
@@ -49,7 +51,7 @@ class BrookCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
        * BrookCalcGBSAOBCForceKernel constructor
        */
 
-      BrookCalcGBSAOBCForceKernel( std::string name, const Platform& platform );
+      BrookCalcGBSAOBCForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
   
       /**
        * BrookCalcGBSAOBCForceKernel destructor
@@ -60,35 +62,34 @@ class BrookCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
       /**
        * Initialize the kernel, setting up the values of all the force field parameters.
        * 
-       * @param atomParameters            vector containing atom index, charge, radius, scalingFactor
-       * @param solventDielectric         solvent dielectric
-       * @param soluteDielectric          solute dielectric
+       * @param system     system this kernel will be applied to
+       * @param force      GBSAOBCForce this kernel will be used for
        *
        */
 
-      void initialize( const std::vector<std::vector<double> >& atomParameters, double solventDielectric, double soluteDielectric );
+      void initialize( const System& system, const GBSAOBCForce& force );
   
       /**
        * Execute the kernel to calculate the forces.
        * 
-       * @param positions   atom coordiantes
+       * @param positions   particle coordiantes
        * @param forces      output forces
        *
        */
 
-      void executeForces( const Stream& positions, Stream& forces );
+      void executeForces( OpenMMContextImpl& context );
   
       /**
        * Execute the kernel to calculate the energy.
        * 
-       * @param positions   atom positions
+       * @param positions   particle positions
        *
        * @return  potential energy due to the NonbondedForce
        * Currently always return 0.0 since energies not calculated on gpu
        *
        */
 
-      double executeEnergy( const Stream& positions );
+      double executeEnergy( OpenMMContextImpl& context );
 
       /** 
        * Set log file reference
@@ -127,16 +128,24 @@ class BrookCalcGBSAOBCForceKernel : public CalcGBSAOBCForceKernel {
 
       FILE* _log;
 
-      // number of atoms
+      // number of particles
 
-      int _numberOfAtoms;
+      int _numberOfParticles;
    
       // Brook Gbsa
 
       BrookGbsa* _brookGbsa;
 
+      // interface
+
+      OpenMMBrookInterface& _openMMBrookInterface;
+
+      // System reference
+
+      System& _system;
+
 };
 
 } // namespace OpenMM
 
-#endif /* OPENMM_BROOK__CALCL_GBSAOBC_FORCEFIELD_KERNEL_H_ */
+#endif /* OPENMM_BROOK_CALC_GBSAOBC_FORCE_KERNEL_H_ */

platforms/brook/src/BrookCalcLJ14ForceKernel.cpp → platforms/brook/src/BrookCalcHarmonicAngleForceKernel.cpp

@@ -31,15 +31,15 @@
 
 #include "OpenMMException.h"
 #include <sstream>
-#include "BrookCalcHarmonicLJ14ForceKernel.h"
+#include "BrookCalcHarmonicAngleForceKernel.h"
 
 using namespace OpenMM;
 using namespace std;
 
-const std::string BrookCalcHarmonicLJ14ForceKernel::BondName = "HarmonicLJ14";
+const std::string BrookCalcHarmonicAngleForceKernel::BondName = "HarmonicAngle";
 
 /** 
- * BrookCalcHarmonicLJ14ForceKernel constructor
+ * BrookCalcHarmonicAngleForceKernel constructor
  * 
  * @param name                      kernel name
  * @param platform                  platform
@@ -48,13 +48,13 @@ const std::string BrookCalcHarmonicLJ14ForceKernel::BondName = "HarmonicLJ14";
  *
  */
 
-BrookCalcHarmonicLJ14ForceKernel::BrookCalcHarmonicLJ14ForceKernel( std::string name, const Platform& platform,
+BrookCalcHarmonicAngleForceKernel::BrookCalcHarmonicAngleForceKernel( std::string name, const Platform& platform,
                                                                       OpenMMBrookInterface& openMMBrookInterface, System& system ) :
-                     CalcHarmonicLJ14ForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
+                     CalcHarmonicAngleForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookCalcHarmonicLJ14ForceKernel::BrookCalcHarmonicLJ14ForceKernel";
+   // static const std::string methodName      = "BrookCalcHarmonicAngleForceKernel::BrookCalcHarmonicAngleForceKernel";
    // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
@@ -70,15 +70,15 @@ BrookCalcHarmonicLJ14ForceKernel::BrookCalcHarmonicLJ14ForceKernel( std::string
 }   
 
 /** 
- * BrookCalcHarmonicLJ14ForceKernel destructor
+ * BrookCalcHarmonicAngleForceKernel destructor
  * 
  */
 
-BrookCalcHarmonicLJ14ForceKernel::~BrookCalcHarmonicLJ14ForceKernel( ){
+BrookCalcHarmonicAngleForceKernel::~BrookCalcHarmonicAngleForceKernel( ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookCalcHarmonicLJ14ForceKernel::BrookCalcHarmonicLJ14ForceKernel";
+   // static const std::string methodName      = "BrookCalcHarmonicAngleForceKernel::BrookCalcHarmonicAngleForceKernel";
    // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
@@ -93,7 +93,7 @@ BrookCalcHarmonicLJ14ForceKernel::~BrookCalcHarmonicLJ14ForceKernel( ){
  *
  */
 
-FILE* BrookCalcHarmonicLJ14ForceKernel::getLog( void ) const {
+FILE* BrookCalcHarmonicAngleForceKernel::getLog( void ) const {
    return _log;
 }
 
@@ -106,7 +106,7 @@ FILE* BrookCalcHarmonicLJ14ForceKernel::getLog( void ) const {
  *
  */
 
-int BrookCalcHarmonicLJ14ForceKernel::setLog( FILE* log ){
+int BrookCalcHarmonicAngleForceKernel::setLog( FILE* log ){
    _log = log;
    return BrookCommon::DefaultReturnValue;
 }
@@ -115,15 +115,15 @@ int BrookCalcHarmonicLJ14ForceKernel::setLog( FILE* log ){
  * Initialize the kernel, setting up the values of all the force field parameters.
  * 
  * @param system                    System reference
- * @param force                     HarmonicLJ14Force reference
+ * @param force                     HarmonicAngleForce reference
  *
  */
 
-void BrookCalcHarmonicLJ14ForceKernel::initialize( const System& system, const HarmonicLJ14Force& force ){
+void BrookCalcHarmonicAngleForceKernel::initialize( const System& system, const HarmonicAngleForce& force ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookCalcHarmonicLJ14ForceKernel::initialize";
+   static const std::string methodName      = "BrookCalcHarmonicAngleForceKernel::initialize";
 
 // ---------------------------------------------------------------------------------------
 
@@ -131,24 +131,24 @@ void BrookCalcHarmonicLJ14ForceKernel::initialize( const System& system, const H
 
    // ---------------------------------------------------------------------------------------
 
-   // create _brookBondParameters object containing atom indices/parameters
+   // create _brookBondParameters object containing particle indices/parameters
 
-   int numberOfBonds         = force.getNumLJ14s();
+   int numberOfBonds         = force.getNumAngles();
 
    if( _brookBondParameters ){
       delete _brookBondParameters;
    }
-   _brookBondParameters = new BrookBondParameters( BondName, NumberOfAtomsInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
+   _brookBondParameters = new BrookBondParameters( BondName, NumberOfParticlesInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
 
    for( int ii = 0; ii < numberOfBonds; ii++ ){
 
       int particle1, particle2, particle3;
       double angle, k;
 
-      int particles[NumberOfAtomsInBond];
+      int particles[NumberOfParticlesInBond];
       double parameters[NumberOfParametersInBond];
 
-      force.getLJ14Parameters( ii, particle1, particle2, particle3, angle, k ); 
+      force.getAngleParameters( ii, particle1, particle2, particle3, angle, k ); 
       particles[0]    = particle1;
       particles[1]    = particle2;
       particles[2]    = particle3;
@@ -158,7 +158,7 @@ void BrookCalcHarmonicLJ14ForceKernel::initialize( const System& system, const H
 
       _brookBondParameters->setBond( ii, particles, parameters );
    }   
-   _openMMBrookInterface.setHarmonicLJ14ForceParameters( _brookBondParameters );
+   _openMMBrookInterface.setHarmonicAngleForceParameters( _brookBondParameters );
    _openMMBrookInterface.setTriggerForceKernel( this );
    _openMMBrookInterface.setTriggerEnergyKernel( this );
 
@@ -173,17 +173,17 @@ void BrookCalcHarmonicLJ14ForceKernel::initialize( const System& system, const H
 }
 
 /** 
- * Compute forces given atom coordinates
+ * Compute forces given particle coordinates
  * 
  * @param context OpenMMContextImpl context
  *
  */
 
-void BrookCalcHarmonicLJ14ForceKernel::executeForces( OpenMMContextImpl& context ){
+void BrookCalcHarmonicAngleForceKernel::executeForces( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName   = "BrookCalcHarmonicLJ14ForceKernel::executeForces";
+   //static const std::string methodName   = "BrookCalcHarmonicAngleForceKernel::executeForces";
 
 // ---------------------------------------------------------------------------------------
 
@@ -205,11 +205,11 @@ void BrookCalcHarmonicLJ14ForceKernel::executeForces( OpenMMContextImpl& context
  *
  */
 
-double BrookCalcHarmonicLJ14ForceKernel::executeEnergy( OpenMMContextImpl& context ){
+double BrookCalcHarmonicAngleForceKernel::executeEnergy( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName      = "BrookCalcHarmonicLJ14ForceKernel::executeEnergy";
+   //static const std::string methodName      = "BrookCalcHarmonicAngleForceKernel::executeEnergy";
 
 // ---------------------------------------------------------------------------------------
 

platforms/brook/src/BrookCalcRBDihedralBondForceKernel.h → platforms/brook/src/BrookCalcHarmonicAngleForceKernel.h

-#ifndef OPENMM_BROOK_CALC_RB_DIHEDRAL_FORCE_KERNEL_H_
-#define OPENMM_BROOK_CALC_RB_DIHEDRAL_FORCE_KERNEL_H_
+#ifndef OPENMM_BROOK_CALC_HARMONIC_ANGLE_FORCE_KERNEL_H_
+#define OPENMM_BROOK_CALC_HARMONIC_ANGLE_FORCE_KERNEL_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -43,36 +43,36 @@ namespace OpenMM {
  * This kernel is invoked to calculate the harmonic angle forces acting on the system.
  */
 
-class BrookCalcRbDihedralForceKernel : public CalcRBTorsionForceKernel {
+class BrookCalcHarmonicAngleForceKernel : public CalcHarmonicAngleForceKernel {
 
    public:
   
       /**
-       * BrookCalcRbDihedralForceKernel constructor
+       * BrookCalcHarmonicAngleForceKernel constructor
        */
 
-      BrookCalcRbDihedralForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
+      BrookCalcHarmonicAngleForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
   
       /**
-       * BrookCalcRbDihedralForceKernel destructor
+       * BrookCalcHarmonicAngleForceKernel destructor
        */
 
-      ~BrookCalcRbDihedralForceKernel();
+      ~BrookCalcHarmonicAngleForceKernel();
   
       /**
        * Initialize the kernel, setting up the values to calculate harmonic bond force & energy
        * 
        * @param system                    System reference
-       * @param force                     RbDihedralForce reference
+       * @param force                     HarmonicAngleForce reference
        *
        */
 
-      void initialize( const System& system, const RBTorsionForce& force );
+      void initialize( const System& system, const HarmonicAngleForce& force );
   
       /**
        * Execute the kernel to calculate the forces.
        * 
-       * @param positions   atom coordiantes
+       * @param positions   particle coordiantes
        * @param forces      output forces
        *
        */
@@ -133,7 +133,7 @@ class BrookCalcRbDihedralForceKernel : public CalcRBTorsionForceKernel {
       /** 
        * Get indices/parameters
        * 
-       * @return  BrookBondParameters containing atom indices/parameters
+       * @return  BrookBondParameters containing particle indices/parameters
        *
        */
       
@@ -141,7 +141,7 @@ class BrookCalcRbDihedralForceKernel : public CalcRBTorsionForceKernel {
       
    private:
    
-      static const int NumberOfAtomsInBond      = 3;
+      static const int NumberOfParticlesInBond  = 3;
       static const int NumberOfParametersInBond = 2;
 
       // bond name
@@ -169,4 +169,4 @@ class BrookCalcRbDihedralForceKernel : public CalcRBTorsionForceKernel {
 
 } // namespace OpenMM
 
-#endif /* OPENMM_BROOK_CALC_RB_DIHEDRAL_FORCE_KERNEL_H_ */
+#endif /* OPENMM_BROOK_CALC_HARMONIC_ANGLE_FORCE_KERNEL_H_ */

platforms/brook/src/BrookCalcRBDihedralBondForceKernel.cpp → platforms/brook/src/BrookCalcHarmonicBondForceKernel.cpp

@@ -31,15 +31,15 @@
 
 #include "OpenMMException.h"
 #include <sstream>
-#include "BrookCalcRBDihedralBondForceKernel.h"
+#include "BrookCalcHarmonicBondForceKernel.h"
 
 using namespace OpenMM;
 using namespace std;
 
-const std::string BrookCalcRbDihedralForceKernel::BondName = "RbDihedral";
+const std::string BrookCalcHarmonicBondForceKernel::BondName = "HarmonicBond";
 
 /** 
- * BrookCalcRbDihedralForceKernel constructor
+ * BrookCalcHarmonicBondForceKernel constructor
  * 
  * @param name                      kernel name
  * @param platform                  platform
@@ -48,13 +48,13 @@ const std::string BrookCalcRbDihedralForceKernel::BondName = "RbDihedral";
  *
  */
 
-BrookCalcRbDihedralForceKernel::BrookCalcRbDihedralForceKernel( std::string name, const Platform& platform,
+BrookCalcHarmonicBondForceKernel::BrookCalcHarmonicBondForceKernel( std::string name, const Platform& platform,
                                                                       OpenMMBrookInterface& openMMBrookInterface, System& system ) :
-                     CalcRbDihedralForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
+                     CalcHarmonicBondForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookCalcRbDihedralForceKernel::BrookCalcRbDihedralForceKernel";
+   // static const std::string methodName      = "BrookCalcHarmonicBondForceKernel::BrookCalcHarmonicBondForceKernel";
    // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
@@ -70,15 +70,15 @@ BrookCalcRbDihedralForceKernel::BrookCalcRbDihedralForceKernel( std::string name
 }   
 
 /** 
- * BrookCalcRbDihedralForceKernel destructor
+ * BrookCalcHarmonicBondForceKernel destructor
  * 
  */
 
-BrookCalcRbDihedralForceKernel::~BrookCalcRbDihedralForceKernel( ){
+BrookCalcHarmonicBondForceKernel::~BrookCalcHarmonicBondForceKernel( ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookCalcRbDihedralForceKernel::BrookCalcRbDihedralForceKernel";
+   // static const std::string methodName      = "BrookCalcHarmonicBondForceKernel::BrookCalcHarmonicBondForceKernel";
    // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
@@ -93,7 +93,7 @@ BrookCalcRbDihedralForceKernel::~BrookCalcRbDihedralForceKernel( ){
  *
  */
 
-FILE* BrookCalcRbDihedralForceKernel::getLog( void ) const {
+FILE* BrookCalcHarmonicBondForceKernel::getLog( void ) const {
    return _log;
 }
 
@@ -106,7 +106,7 @@ FILE* BrookCalcRbDihedralForceKernel::getLog( void ) const {
  *
  */
 
-int BrookCalcRbDihedralForceKernel::setLog( FILE* log ){
+int BrookCalcHarmonicBondForceKernel::setLog( FILE* log ){
    _log = log;
    return BrookCommon::DefaultReturnValue;
 }
@@ -115,15 +115,15 @@ int BrookCalcRbDihedralForceKernel::setLog( FILE* log ){
  * Initialize the kernel, setting up the values of all the force field parameters.
  * 
  * @param system                    System reference
- * @param force                     RbDihedralForce reference
+ * @param force                     HarmonicBondForce reference
  *
  */
 
-void BrookCalcRbDihedralForceKernel::initialize( const System& system, const RbDihedralForce& force ){
+void BrookCalcHarmonicBondForceKernel::initialize( const System& system, const HarmonicBondForce& force ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookCalcRbDihedralForceKernel::initialize";
+   static const std::string methodName      = "BrookCalcHarmonicBondForceKernel::initialize";
 
 // ---------------------------------------------------------------------------------------
 
@@ -131,40 +131,39 @@ void BrookCalcRbDihedralForceKernel::initialize( const System& system, const RbD
 
    // ---------------------------------------------------------------------------------------
 
-   // create _brookBondParameters object containing atom indices/parameters
+   // create _brookBondParameters object containing particle indices/parameters
 
-   int numberOfBonds         = force.getNumAngles();
+   int numberOfBonds         = force.getNumBonds();
 
    if( _brookBondParameters ){
       delete _brookBondParameters;
    }
-   _brookBondParameters = new BrookBondParameters( BondName, NumberOfAtomsInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
+   _brookBondParameters = new BrookBondParameters( BondName, NumberOfParticlesInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
 
    for( int ii = 0; ii < numberOfBonds; ii++ ){
 
-      int particle1, particle2, particle3;
-      double angle, k;
+      int particle1, particle2;
+      double length, k;
 
-      int particles[NumberOfAtomsInBond];
+      int particles[NumberOfParticlesInBond];
       double parameters[NumberOfParametersInBond];
 
-      force.getAngleParameters( ii, particle1, particle2, particle3, angle, k ); 
+      force.getBondParameters( ii, particle1, particle2, length, k ); 
       particles[0]    = particle1;
       particles[1]    = particle2;
-      particles[2]    = particle3;
  
-      parameters[0]   = angle;
+      parameters[0]   = length;
       parameters[1]   = k;
 
       _brookBondParameters->setBond( ii, particles, parameters );
    }   
-   _openMMBrookInterface.setRbDihedralForceParameters( _brookBondParameters );
+   _openMMBrookInterface.setHarmonicBondForceParameters( _brookBondParameters );
    _openMMBrookInterface.setTriggerForceKernel( this );
    _openMMBrookInterface.setTriggerEnergyKernel( this );
 
    if( log ){
       std::string contents = _brookBondParameters->getContentsString( ); 
-      (void) fprintf( log, "%s brookGbsa::contents\n%s", methodName.c_str(), contents.c_str() );
+      (void) fprintf( log, "%s contents\n%s", methodName.c_str(), contents.c_str() );
       (void) fflush( log );
    }
 
@@ -173,17 +172,17 @@ void BrookCalcRbDihedralForceKernel::initialize( const System& system, const RbD
 }
 
 /** 
- * Compute forces given atom coordinates
+ * Compute forces given particle coordinates
  * 
  * @param context OpenMMContextImpl context
  *
  */
 
-void BrookCalcRbDihedralForceKernel::executeForces( OpenMMContextImpl& context ){
+void BrookCalcHarmonicBondForceKernel::executeForces( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName   = "BrookCalcRbDihedralForceKernel::executeForces";
+   //static const std::string methodName   = "BrookCalcHarmonicBondForceKernel::executeForces";
 
 // ---------------------------------------------------------------------------------------
 
@@ -205,11 +204,11 @@ void BrookCalcRbDihedralForceKernel::executeForces( OpenMMContextImpl& context )
  *
  */
 
-double BrookCalcRbDihedralForceKernel::executeEnergy( OpenMMContextImpl& context ){
+double BrookCalcHarmonicBondForceKernel::executeEnergy( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName      = "BrookCalcRbDihedralForceKernel::executeEnergy";
+   //static const std::string methodName      = "BrookCalcHarmonicBondForceKernel::executeEnergy";
 
 // ---------------------------------------------------------------------------------------
 

platforms/brook/src/BrookNonbonded14ForceKernel.h → platforms/brook/src/BrookCalcHarmonicBondForceKernel.h

-#ifndef OPENMM_BROOK_CALC_LJ14_FORCE_KERNEL_H_
-#define OPENMM_BROOK_CALC_LJ14_FORCE_KERNEL_H_
+#ifndef OPENMM_BROOK_CALC_HARMONIC_BOND_FORCE_KERNEL_H_
+#define OPENMM_BROOK_CALC_HARMONIC_BOND_FORCE_KERNEL_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -43,36 +43,36 @@ namespace OpenMM {
  * This kernel is invoked to calculate the harmonic angle forces acting on the system.
  */
 
-class BrookCalcHarmonicLJ14ForceKernel : public CalcHarmonicLJ14ForceKernel {
+class BrookCalcHarmonicBondForceKernel : public CalcHarmonicBondForceKernel {
 
    public:
   
       /**
-       * BrookCalcHarmonicLJ14ForceKernel constructor
+       * BrookCalcHarmonicBondForceKernel constructor
        */
 
-      BrookCalcHarmonicLJ14ForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
+      BrookCalcHarmonicBondForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
   
       /**
-       * BrookCalcHarmonicLJ14ForceKernel destructor
+       * BrookCalcHarmonicBondForceKernel destructor
        */
 
-      ~BrookCalcHarmonicLJ14ForceKernel();
+      ~BrookCalcHarmonicBondForceKernel();
   
       /**
        * Initialize the kernel, setting up the values to calculate harmonic bond force & energy
        * 
        * @param system                    System reference
-       * @param force                     HarmonicLJ14Force reference
+       * @param force                     HarmonicBondForce reference
        *
        */
 
-      void initialize( const System& system, const HarmonicLJ14Force& force );
+      void initialize( const System& system, const HarmonicBondForce& force );
   
       /**
        * Execute the kernel to calculate the forces.
        * 
-       * @param positions   atom coordiantes
+       * @param positions   particle coordiantes
        * @param forces      output forces
        *
        */
@@ -133,7 +133,7 @@ class BrookCalcHarmonicLJ14ForceKernel : public CalcHarmonicLJ14ForceKernel {
       /** 
        * Get indices/parameters
        * 
-       * @return  BrookBondParameters containing atom indices/parameters
+       * @return  BrookBondParameters containing particle indices/parameters
        *
        */
       
@@ -141,7 +141,7 @@ class BrookCalcHarmonicLJ14ForceKernel : public CalcHarmonicLJ14ForceKernel {
       
    private:
    
-      static const int NumberOfAtomsInBond      = 3;
+      static const int NumberOfParticlesInBond  = 2;
       static const int NumberOfParametersInBond = 2;
 
       // bond name
@@ -152,7 +152,6 @@ class BrookCalcHarmonicLJ14ForceKernel : public CalcHarmonicLJ14ForceKernel {
 
       FILE* _log;
 
-
       // Brook bond parameters
 
       BrookBondParameters* _brookBondParameters;
@@ -169,4 +168,4 @@ class BrookCalcHarmonicLJ14ForceKernel : public CalcHarmonicLJ14ForceKernel {
 
 } // namespace OpenMM
 
-#endif /* OPENMM_BROOK_CALC_LJ14_FORCE_KERNEL_H_ */
+#endif /* OPENMM_BROOK_CALC_HARMONIC_BOND_FORCE_KERNEL_H_ */

platforms/brook/src/BrookCalcKineticEnergyKernel.cpp

@@ -77,7 +77,7 @@ BrookCalcKineticEnergyKernel::~BrookCalcKineticEnergyKernel( ){
 /** 
  * Initialize the kernel
  * 
- * @param masses   mass of each atom
+ * @param masses   mass of each particle
  *
  */
 
@@ -107,7 +107,7 @@ void BrookCalcKineticEnergyKernel::initialize( const vector<double>& masses ){
 /** 
  * Execute kernel
  * 
- * @param velocities  stream of atom velocities
+ * @param velocities  stream of particle velocities
  *
  * @return kinetic energy of the system
  *

platforms/brook/src/BrookCalcKineticEnergyKernel.h

@@ -65,7 +65,7 @@ class BrookCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
       /** 
        * Initialize the kernel
        * 
-       * @param masses   mass of each atom
+       * @param masses   mass of each particle
        *
        */
       void initialize( const std::vector<double>& masses );
@@ -73,7 +73,7 @@ class BrookCalcKineticEnergyKernel : public CalcKineticEnergyKernel {
       /** 
        * Execute the kernel.
        * 
-       * @param velocities stream of atom velocities
+       * @param velocities stream of particle velocities
        *
        */
 

platforms/brook/src/BrookCalcLJ14ForceKernel.h

-#ifndef OPENMM_BROOK_CALC_LJ14_FORCE_KERNEL_H_
-#define OPENMM_BROOK_CALC_LJ14_FORCE_KERNEL_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 "kernels.h"
-#include "BrookPlatform.h"
-#include "BrookBondParameters.h"
-#include "OpenMMBrookInterface.h"
-
-namespace OpenMM {
-
-/**
- * This kernel is invoked to calculate the harmonic angle forces acting on the system.
- */
-
-class BrookCalcHarmonicLJ14ForceKernel : public CalcHarmonicLJ14ForceKernel {
-
-   public:
-  
-      /**
-       * BrookCalcHarmonicLJ14ForceKernel constructor
-       */
-
-      BrookCalcHarmonicLJ14ForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
-  
-      /**
-       * BrookCalcHarmonicLJ14ForceKernel destructor
-       */
-
-      ~BrookCalcHarmonicLJ14ForceKernel();
-  
-      /**
-       * Initialize the kernel, setting up the values to calculate harmonic bond force & energy
-       * 
-       * @param system                    System reference
-       * @param force                     HarmonicLJ14Force reference
-       *
-       */
-
-      void initialize( const System& system, const HarmonicLJ14Force& force );
-  
-      /**
-       * Execute the kernel to calculate the forces.
-       * 
-       * @param positions   atom coordiantes
-       * @param forces      output forces
-       *
-       */
-
-      void executeForces( OpenMMContextImpl& context );
-  
-      /**
-       * Execute the kernel to calculate the energy.
-       * 
-       * @param context    the context in which to execute this kernel
-       *
-       * @return  potential energy associated with the harmonic angle force
-       *
-       */
-
-      double executeEnergy( OpenMMContextImpl& context );
-
-      /** 
-       * Set log file reference
-       * 
-       * @param  log file reference
-       *
-       * @return DefaultReturnValue
-       *
-       */
-      
-      int setLog( FILE* log );
-
-      /* 
-       * Get contents of object
-       *
-       * @param level of dump
-       *
-       * @return string containing contents
-       *
-       * */
-      
-      std::string getContents( int level ) const;
-
-      /** 
-       * Get log file reference
-       * 
-       * @return  log file reference
-       *
-       */
-      
-      FILE* getLog( void ) const;
-      
-      /** 
-       * Get number of bonds
-       * 
-       * @return  number of bonds
-       *
-       */
-      
-      int getNumberOfBonds( void ) const;
-      
-      /** 
-       * Get indices/parameters
-       * 
-       * @return  BrookBondParameters containing atom indices/parameters
-       *
-       */
-      
-      BrookBondParameters* getBrookBondParameters( void ) const;
-      
-   private:
-   
-      static const int NumberOfAtomsInBond      = 3;
-      static const int NumberOfParametersInBond = 2;
-
-      // bond name
-
-      static const std::string BondName;
-
-      // log file reference
-
-      FILE* _log;
-
-
-      // Brook bond parameters
-
-      BrookBondParameters* _brookBondParameters;
-
-      // interface
-
-      OpenMMBrookInterface& _openMMBrookInterface;
-
-      // System reference
-
-      System& _system;
-
-};
-
-} // namespace OpenMM
-
-#endif /* OPENMM_BROOK_CALC_LJ14_FORCE_KERNEL_H_ */

platforms/brook/src/BrookCalcProperDihedralBondForceKernel.cpp → platforms/brook/src/BrookCalcNonbondedForceKernel.cpp

@@ -29,40 +29,51 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include <cmath>
+#include <limits>
 #include "OpenMMException.h"
 #include <sstream>
-#include "BrookCalcProperDihedralForceKernel.h"
+
+#include "BrookStreamImpl.h"
+#include "BrookCalcNonbondedForceKernel.h"
+#include "NonbondedForce.h"
 
 using namespace OpenMM;
 using namespace std;
 
-const std::string BrookCalcProperDihedralForceKernel::BondName = "ProperDihedral";
+const std::string BrookCalcNonbondedForceKernel::BondName = "LJ14";
 
 /** 
- * BrookCalcProperDihedralForceKernel constructor
+ * BrookCalcNonbondedForceKernel constructor
  * 
  * @param name                      kernel name
  * @param platform                  platform
- * @param OpenMMBrookInterface      OpenMM-Brook interface
- * @param System                    System reference
  *
  */
 
-BrookCalcProperDihedralForceKernel::BrookCalcProperDihedralForceKernel( std::string name, const Platform& platform,
-                                                                      OpenMMBrookInterface& openMMBrookInterface, System& system ) :
-                     CalcProperDihedralForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
+BrookCalcNonbondedForceKernel::BrookCalcNonbondedForceKernel( std::string name, const Platform& platform,
+                                                              OpenMMBrookInterface& openMMBrookInterface, System& system ) :
+                     CalcNonbondedForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookCalcProperDihedralForceKernel::BrookCalcProperDihedralForceKernel";
-   // static const int debug                   = 1;
+   // static const std::string methodName      = "BrookCalcNonbondedForceKernel::BrookCalcNonbondedForceKernel";
 
 // ---------------------------------------------------------------------------------------
 
-   _brookBondParameters              = NULL;
-   _log                              = NULL;
+   _numberOfParticles                       = 0;
+   _brookNonBonded                          = NULL;
+   _brookBondParameters                     = NULL;
+
+   _log                                     = NULL;
+
+   _refForceField                           = NULL;
+   _refSystem                               = NULL;
+   _refOpenMMContext                        = NULL;
+   _referencePlatform                       = NULL;
+   _refVerletIntegrator                     = NULL;
 
-   const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
+   const BrookPlatform brookPlatform        = dynamic_cast<const BrookPlatform&> (platform);
    if( brookPlatform.getLog() != NULL ){
       setLog( brookPlatform.getLog() );
    }
@@ -70,20 +81,30 @@ BrookCalcProperDihedralForceKernel::BrookCalcProperDihedralForceKernel( std::str
 }   
 
 /** 
- * BrookCalcProperDihedralForceKernel destructor
+ * BrookCalcNonbondedForceKernel destructor
  * 
  */
 
-BrookCalcProperDihedralForceKernel::~BrookCalcProperDihedralForceKernel( ){
+BrookCalcNonbondedForceKernel::~BrookCalcNonbondedForceKernel( ){
 
 // ---------------------------------------------------------------------------------------
 
-   // static const std::string methodName      = "BrookCalcProperDihedralForceKernel::BrookCalcProperDihedralForceKernel";
-   // static const int debug                   = 1;
+   // static const std::string methodName      = "BrookCalcNonbondedForceKernel::BrookCalcNonbondedForceKernel";
 
 // ---------------------------------------------------------------------------------------
 
-   delete _brookBondParameters;
+   //delete _brookBondParameters;
+   delete _brookNonBonded;
+
+   // deleted w/ kernel delete? If activated, program crashes
+
+   //delete _refForceField;
+/*
+   delete _refSystem;
+   delete _refOpenMMContext;
+   delete _referencePlatform;
+   delete _refVerletIntegrator;
+*/
 }
 
 /** 
@@ -93,7 +114,7 @@ BrookCalcProperDihedralForceKernel::~BrookCalcProperDihedralForceKernel( ){
  *
  */
 
-FILE* BrookCalcProperDihedralForceKernel::getLog( void ) const {
+FILE* BrookCalcNonbondedForceKernel::getLog( void ) const {
    return _log;
 }
 
@@ -106,24 +127,98 @@ FILE* BrookCalcProperDihedralForceKernel::getLog( void ) const {
  *
  */
 
-int BrookCalcProperDihedralForceKernel::setLog( FILE* log ){
+int BrookCalcNonbondedForceKernel::setLog( FILE* log ){
    _log = log;
    return BrookCommon::DefaultReturnValue;
 }
 
+/** 
+ * Set log file reference
+ * 
+ * @param  log file reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+void BrookCalcNonbondedForceKernel::initialize( const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookCalcNonbondedForceKernel::initialize";
+
+// ---------------------------------------------------------------------------------------
+
+    FILE* log                 = getLog();
+    _numberOfParticles        = force.getNumParticles();
+/*
+    nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
+    nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
+    Vec3 boxVectors[3];
+    force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
+    periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
+    periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
+    if (nonbondedMethod == NoCutoff)
+        neighborList = NULL;
+    else
+        neighborList = new NeighborList();
+*/
+   // ---------------------------------------------------------------------------------------
+
+   // nonbonded
+
+   if( _brookNonBonded ){
+      delete _brookNonBonded;
+   }
+   _brookNonBonded           = new BrookNonBonded();
+   _brookNonBonded->setLog( log );
+
+   // charge & LJ parameters
+
+   std::vector<std::vector<double> > nonbondedParameters;
+   for( int ii = 0; ii < _numberOfParticles; ii++ ){
+      double charge, radius, depth;
+      force.getParticleParameters( ii, charge, radius, depth );
+      std::vector<double> particleParamArray;
+      nonbondedParameters[ii] = particleParamArray;
+      particleParamArray[0]   = radius;
+      particleParamArray[1]   = depth;
+      particleParamArray[2]   = charge;
+   }   
+
+   _brookNonBonded->setup( _numberOfParticles, nonbondedParameters, exclusions, getPlatform() );
+   _openMMBrookInterface.setTriggerForceKernel( this );
+   _openMMBrookInterface.setTriggerEnergyKernel( this );
+
+
+   // echo contents
+
+   if( log ){
+      std::string contents = _brookNonBonded->getContentsString( );
+      (void) fprintf( log, "%s brookNonBonded::contents\n%s", methodName.c_str(), contents.c_str() );
+      (void) fflush( log );
+   }
+
+   // nonbonded 14 ixns
+
+   initialize14Interactions( system, force );
+
+}
+
 /** 
  * Initialize the kernel, setting up the values of all the force field parameters.
  * 
  * @param system                    System reference
- * @param force                     ProperDihedralForce reference
+ * @param force                     HarmonicLJ14Force reference
  *
  */
 
-void BrookCalcProperDihedralForceKernel::initialize( const System& system, const ProperDihedralForce& force ){
+void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& system, const NonbondedForce& force ){
 
 // ---------------------------------------------------------------------------------------
 
-   static const std::string methodName      = "BrookCalcProperDihedralForceKernel::initialize";
+   static const std::string methodName      = "BrookCalcNonbondedForceKernel::initialize14Interactions";
 
 // ---------------------------------------------------------------------------------------
 
@@ -131,85 +226,79 @@ void BrookCalcProperDihedralForceKernel::initialize( const System& system, const
 
    // ---------------------------------------------------------------------------------------
 
-   // create _brookBondParameters object containing atom indices/parameters
+   // create _brookBondParameters object containing particle indices/parameters
 
-   int numberOfBonds         = force.getNumAngles();
+   int numberOf14Forces         = force.getNumNonbonded14();
 
-   if( _brookBondParameters ){
-      delete _brookBondParameters;
-   }
-   _brookBondParameters = new BrookBondParameters( BondName, NumberOfAtomsInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
+   //delete _brookBondParameters;
+   _brookBondParameters = new BrookBondParameters( BondName, NumberOfParticlesInBond, NumberOfParametersInBond, numberOf14Forces, getLog() );
 
-   for( int ii = 0; ii < numberOfBonds; ii++ ){
+   for( int ii = 0; ii < numberOf14Forces; ii++ ){
 
-      int particle1, particle2, particle3;
-      double angle, k;
+      int particle1, particle2;
+      double  charge, radius, depth;
 
-      int particles[NumberOfAtomsInBond];
+      int particles[NumberOfParticlesInBond];
       double parameters[NumberOfParametersInBond];
 
-      force.getAngleParameters( ii, particle1, particle2, particle3, angle, k ); 
+      force.getNonbonded14Parameters( ii, particle1, particle2,  charge, radius, depth ); 
       particles[0]    = particle1;
       particles[1]    = particle2;
-      particles[2]    = particle3;
  
-      parameters[0]   = angle;
-      parameters[1]   = k;
+      parameters[0]   = charge;
+      parameters[1]   = radius;
+      parameters[2]   = depth;
 
       _brookBondParameters->setBond( ii, particles, parameters );
    }   
-   _openMMBrookInterface.setProperDihedralForceParameters( _brookBondParameters );
-   _openMMBrookInterface.setTriggerForceKernel( this );
-   _openMMBrookInterface.setTriggerEnergyKernel( this );
+   _openMMBrookInterface.setNonBonded14ForceParameters( _brookBondParameters );
 
    if( log ){
       std::string contents = _brookBondParameters->getContentsString( ); 
-      (void) fprintf( log, "%s brookGbsa::contents\n%s", methodName.c_str(), contents.c_str() );
+      (void) fprintf( log, "%s contents:\n%s", methodName.c_str(), contents.c_str() );
       (void) fflush( log );
    }
 
    // ---------------------------------------------------------------------------------------
     
 }
-
 /** 
- * Compute forces given atom coordinates
+ * Execute the kernel to calculate the nonbonded forces
  * 
  * @param context OpenMMContextImpl context
  *
  */
 
-void BrookCalcProperDihedralForceKernel::executeForces( OpenMMContextImpl& context ){
+void BrookCalcNonbondedForceKernel::executeForces( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName   = "BrookCalcProperDihedralForceKernel::executeForces";
+   static const std::string methodName      = "BrookCalcNonbondedForceKernel::executeForces";
 
 // ---------------------------------------------------------------------------------------
 
    if( _openMMBrookInterface.getTriggerForceKernel() == this ){
       _openMMBrookInterface.computeForces( context );
-   }
-
-   return;
+   }   
 
    // ---------------------------------------------------------------------------------------
 }
 
 /**
- * Execute the kernel to calculate the energy
+ * Execute the kernel to calculate the energy.
  * 
  * @param context OpenMMContextImpl context
  *
- * @return  potential energy
+ * @return  potential energy due to the NonbondedForce
+ * Currently always return 0.0 since energies not calculated on gpu
  *
  */
 
-double BrookCalcProperDihedralForceKernel::executeEnergy( OpenMMContextImpl& context ){
+double BrookCalcNonbondedForceKernel::executeEnergy( OpenMMContextImpl& context ){
 
 // ---------------------------------------------------------------------------------------
 
-   //static const std::string methodName      = "BrookCalcProperDihedralForceKernel::executeEnergy";
+   //static const std::string methodName      = "BrookCalcNonbondedForceKernel::executeEnergy";
 
 // ---------------------------------------------------------------------------------------
 
@@ -217,6 +306,38 @@ double BrookCalcProperDihedralForceKernel::executeEnergy( OpenMMContextImpl& con
       return (double) _openMMBrookInterface.computeEnergy( context );
    } else {
       return 0.0;
+   }   
+
+}
+
+/**
+ * Get reference Context
+ * 
+ * @param numberOfParticles  number of particles
+ *
+ * @return  OpenMMContext
+ *
+ */
+
+OpenMMContext* BrookCalcNonbondedForceKernel::getReferenceOpenMMContext( int numberOfParticles ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookCalcNonbondedForceKernel::getReferenceOpenMMContext";
+
+// ---------------------------------------------------------------------------------------
+
+   if( _refOpenMMContext == NULL ){
+      
+      _referencePlatform       = new ReferencePlatform();
+      _refSystem               = new System( numberOfParticles, 0 ); 
+		_refVerletIntegrator     = new VerletIntegrator( 0.01 );
+
+      _refSystem->addForce( _refForceField );
+
+      _refOpenMMContext        = new OpenMMContext( *_refSystem, *_refVerletIntegrator, *_referencePlatform );
    }
 
+   return _refOpenMMContext;
+
 }

platforms/brook/src/BrookCalcStandardMMForceFieldKernel.h → platforms/brook/src/BrookCalcNonbondedForceKernel.h

-#ifndef OPENMM_BROOK_CALC_STANDARD_MM_FORCEFIELD_KERNEL_H_
-#define OPENMM_BROOK_CALC_STANDARD_MM_FORCEFIELD_KERNEL_H_
+#ifndef OPENMM_BROOK_CALC_NONBONDED_FORCE_KERNEL_H_
+#define OPENMM_BROOK_CALC_NONBONDED_FORCE_KERNEL_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -33,14 +33,9 @@
  * -------------------------------------------------------------------------- */
 
 #include "kernels.h"
-#include "../../reference/src/SimTKUtilities/SimTKOpenMMRealType.h"
-#include "BrookBonded.h"
-#include "BrookNonBonded.h"
+//#include "../../reference/src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "OpenMMBrookInterface.h"
 #include "NonbondedForce.h"
-#include "OpenMMContext.h"
-#include "System.h"
-#include "ReferencePlatform.h"
-#include "VerletIntegrator.h"
 
 namespace OpenMM {
 
@@ -51,57 +46,53 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
 
    public:
   
-      BrookCalcNonbondedForceKernel( std::string name, const Platform& platform );
+      BrookCalcNonbondedForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
   
       ~BrookCalcNonbondedForceKernel();
   
-      /**
-       * Initialize the kernel, setting up the values of all the force field parameters.
+      /** 
+       * Initialize the kernel
        * 
-       * @param bondIndices               the two atoms connected by each bond term
-       * @param bondParameters            the force parameters (length, k) for each bond term
-       * @param angleIndices              the three atoms connected by each angle term
-       * @param angleParameters           the force parameters (angle, k) for each angle term
-       * @param periodicTorsionIndices    the four atoms connected by each periodic torsion term
-       * @param periodicTorsionParameters the force parameters (k, phase, periodicity) for each periodic torsion term
-       * @param rbTorsionIndices          the four atoms connected by each Ryckaert-Bellemans torsion term
-       * @param rbTorsionParameters       the coefficients (in order of increasing powers) for each Ryckaert-Bellemans torsion term
-       * @param bonded14Indices           each element contains the indices of two atoms whose nonbonded interactions should be reduced since
-       *                                  they form a bonded 1-4 pair
-       * @param lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
-       * @param coulomb14Scale            the factor by which Coulomb interactions should be reduced for bonded 1-4 pairs
-       * @param exclusions                the i'th element lists the indices of all atoms with which the i'th atom should not interact through
-       *                                  nonbonded forces.  Bonded 1-4 pairs are also included in this list, since they should be omitted from
-       *                                  the standard nonbonded calculation.
-       * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
-       * @param nonbondedMethod           the method to use for handling long range nonbonded interactions
-       * @param nonbondedCutoff           the cutoff distance for nonbonded interactions (if nonbondedMethod involves a cutoff)
-       * @param periodicBoxSize           the size of the periodic box (if nonbondedMethod involves a periodic boundary conditions)
-       *
+       * @param system     the System this kernel will be applied to
+       * @param force      the NonbondedForce this kernel will be used for
+       * @param exclusions the i'th element lists the indices of all particles with which the i'th particle should not interact through
+       *                   nonbonded forces.  Bonded 1-4 pairs are also included in this list, since they should be omitted from
+       *                   the standard nonbonded calculation.
+       */
+
+      void initialize( const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions );
+  
+      /** 
+       * Initialize the 14 ixns 
+       * 
+       * @param system     the System this kernel will be applied to
+       * @param force      the NonbondedForce this kernel will be used for
        */
 
-      void initialize( const std::vector<std::vector<int> >& bondIndices,      const std::vector<std::vector<double> >& bondParameters,
-                       const std::vector<std::vector<int> >& angleIndices,     const std::vector<std::vector<double> >& angleParameters,
-                       const std::vector<std::vector<int> >& periodicTorsionIndices, const std::vector<std::vector<double> >& periodicTorsionParameters,
-                       const std::vector<std::vector<int> >& rbTorsionIndices, const std::vector<std::vector<double> >& rbTorsionParameters,
-                       const std::vector<std::vector<int> >& bonded14Indices,  double lj14Scale, double coulomb14Scale,
-                       const std::vector<std::set   <int> >& exclusions,       const std::vector<std::vector<double> >& nonbondedParameters,
-                       NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3] );
+      void initialize14Interactions( const System& system, const NonbondedForce& force );
   
       /**
        * Execute the kernel to calculate the forces.
        * 
-       * @param positions   a Stream of type Double3 containing the position (x, y, z) of each atom
-       * @param forces      a Stream of type Double3 containing the force (x, y, z) on each atom.  On entry, this contains the forces that
+       * @param positions   a Stream of type Double3 containing the position (x, y, z) of each particle
+       * @param forces      a Stream of type Double3 containing the force (x, y, z) on each particle.  On entry, this contains the forces that
        *                    have been calculated so far.  The kernel should add its own forces to the values already in the stream.
        */
 
       void executeForces( const Stream& positions, Stream& forces );
+
+      /** 
+       * Execute the kernel to calculate the forces.
+       * 
+       * @param context    the context in which to execute this kernel
+       */
+
+      void executeForces( OpenMMContextImpl& context );
   
       /**
        * Execute the kernel to calculate the energy.
        * 
-       * @param positions   a Stream of type Double3 containing the position (x, y, z) of each atom
+       * @param positions   a Stream of type Double3 containing the position (x, y, z) of each particle
        *
        * @return the potential energy due to the NonbondedForce
        *
@@ -109,19 +100,26 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
        */
 
       double executeEnergy( const Stream& positions );
-      double executeEnergyOld( const Stream& positions );
+
+      /** 
+       * Execute the kernel to calculate the energy.
+       * 
+       * @param context    the context in which to execute this kernel
+       * @return the potential energy due to the NonbondedForce
+       */
+  
+      double executeEnergy( OpenMMContextImpl& context );
 
       /**
        * Get reference Context
        * 
-       * @param numberOfAtoms  number of atoms
+       * @param numberOfParticles  number of particles
        *
        * @return  OpenMMContext
        *
        */
       
-      OpenMMContext* getReferenceOpenMMContext( int numberOfAtoms );
-      
+      OpenMMContext* getReferenceOpenMMContext( int numberOfParticles );
       
       /** 
        * Set log file reference
@@ -156,19 +154,30 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
       
    private:
    
+      // LJ14 'bond' name
+
+      static const std::string BondName;
+
+      static const int NumberOfParticlesInBond  = 2;
+      static const int NumberOfParametersInBond = 3;
+
       // log file reference
 
       FILE* _log;
 
-       // number of atoms
+       // number of particles
 
-       int _numberOfAtoms;
+       int _numberOfParticles;
    
-       // Brook bonded & nonbonded
+       // Brook nonbonded
 
-       BrookBonded* _brookBonded;
        BrookNonBonded* _brookNonBonded;
 
+       OpenMMBrookInterface& _openMMBrookInterface;
+       System& _system;
+
+       BrookBondParameters* _brookBondParameters;
+
        // used to calculate energy
 
        NonbondedForce* _refForceField;
@@ -181,4 +190,4 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
 
 } // namespace OpenMM
 
-#endif /* OPENMM_BROOK_CALC_STANDARD_MM_FORCEFIELD_KERNEL_H_ */
+#endif /* OPENMM_BROOK_CALC_NONBONDED_FORCE_KERNEL_H_ */

platforms/brook/src/BrookCalcPeriodicTorsionForceKernel.cpp

@@ -131,21 +131,21 @@ void BrookCalcPeriodicTorsionForceKernel::initialize( const System& system, cons
 
    // ---------------------------------------------------------------------------------------
 
-   // create _brookBondParameters object containing atom indices/parameters
+   // create _brookBondParameters object containing particle indices/parameters
 
    int numberOfBonds         = force.getNumTorsions();
 
    if( _brookBondParameters ){
       delete _brookBondParameters;
    }
-   _brookBondParameters = new BrookBondParameters( BondName, NumberOfAtomsInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
+   _brookBondParameters = new BrookBondParameters( BondName, NumberOfParticlesInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
 
    for( int ii = 0; ii < numberOfBonds; ii++ ){
 
       int particle1, particle2, particle3, particle4, periodicity;
       double phase, k;
 
-      int particles[NumberOfAtomsInBond];
+      int particles[NumberOfParticlesInBond];
       double parameters[NumberOfParametersInBond];
 
       force.getTorsionParameters( ii, particle1, particle2, particle3, particle4, periodicity, phase, k ); 
@@ -175,7 +175,7 @@ void BrookCalcPeriodicTorsionForceKernel::initialize( const System& system, cons
 }
 
 /** 
- * Compute forces given atom coordinates
+ * Compute forces given particle coordinates
  * 
  * @param context OpenMMContextImpl context
  *

platforms/brook/src/BrookCalcPeriodicTorsionForceKernel.h

-#ifndef OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_
-#define OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_
+#ifndef OPENMM_BROOK_CALC_PERIODIC_TORSION_FORCE_KERNEL_H_
+#define OPENMM_BROOK_CALC_PERIODIC_TORSION_FORCE_KERNEL_H_
 
 /* -------------------------------------------------------------------------- *
  *                                   OpenMM                                   *
@@ -72,7 +72,7 @@ class BrookCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKerne
       /**
        * Execute the kernel to calculate the forces.
        * 
-       * @param positions   atom coordiantes
+       * @param positions   particle coordiantes
        * @param forces      output forces
        *
        */
@@ -133,7 +133,7 @@ class BrookCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKerne
       /** 
        * Get indices/parameters
        * 
-       * @return  BrookBondParameters containing atom indices/parameters
+       * @return  BrookBondParameters containing particle indices/parameters
        *
        */
       
@@ -141,7 +141,7 @@ class BrookCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKerne
       
    private:
    
-      static const int NumberOfAtomsInBond      = 4;
+      static const int NumberOfParticlesInBond  = 4;
       static const int NumberOfParametersInBond = 3;
 
       // bond name
@@ -169,4 +169,4 @@ class BrookCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKerne
 
 } // namespace OpenMM
 
-#endif /* OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_ */
+#endif /* OPENMM_BROOK_CALC_PERIODIC_TORSION_FORCE_KERNEL_H_ */