Bonded & Nonbonded classes

platforms/brook/src/BrookBonded.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Mark Friedrichs                                                   *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <sstream>
+#include "BrookBonded.h"
+#include "BrookPlatform.h"
+#include "BrookStreamFactory.h"
+#include "OpenMMException.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 PARAMS(X,Y,Z) (params[(Y)][4*(X) + Z])
+
+/** 
+ * BrookBonded constructor
+ * 
+ */
+
+BrookBonded::BrookBonded( ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::BrookBonded";
+
+// ---------------------------------------------------------------------------------------
+
+   _numberOfAtoms             = 0;
+   _ljScale                   = 1.0;
+   _coulombFactor             = 332.0;
+
+   _atomIndicesStream         = NULL;
+   _chargeStream              = NULL;
+
+   // parameter streams
+
+   for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
+      _bondedParameters[ii] = NULL;
+   }
+
+   // inverse maps & force streams
+
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      _bondedForceStreams[ii]       = NULL;
+      _inverseMapStreamCount[ii]    = 0;
+      for( int jj = 0; jj < MaxNumberOfInverseMaps; jj++ ){
+         _inverseStreamMaps[ii][jj] = NULL;
+      }
+   }
+
+   _maxInverseMapStreamCount[StreamI] = 9;
+   _maxInverseMapStreamCount[StreamJ] = 6;
+   _maxInverseMapStreamCount[StreamK] = 6;
+   _maxInverseMapStreamCount[StreamL] = 9;
+
+   _maxNumberOfInverseMaps = _maxInverseMapStreamCount[0];
+   for( int ii = 1; ii < getNumberOfForceStreams(); ii++ ){
+      if( _maxNumberOfInverseMaps < _maxInverseMapStreamCount[ii] ){
+         _maxNumberOfInverseMaps  = _maxInverseMapStreamCount[ii];
+      }
+   }
+
+   // check that MaxNumberOfInverseMaps is big enough
+
+   if( _maxNumberOfInverseMaps > MaxNumberOfInverseMaps ){
+      std::stringstream message;
+      message << methodName << " max number of inverse maps=" << _maxNumberOfInverseMaps << " is greater than hardwired value=" << MaxNumberOfInverseMaps;
+      throw OpenMMException( message.str() );
+   }
+
+   _invMapStreamWidth = -1;
+}   
+ 
+/** 
+ * BrookBonded destructor
+ * 
+ */
+
+BrookBonded::~BrookBonded( ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName = "BrookBonded::~BrookBonded";
+
+// ---------------------------------------------------------------------------------------
+
+   delete _atomIndicesStream;
+   delete _chargeStream;
+
+   for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
+      delete _bondedParameters[ii];
+   }
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      delete _bondedForceStreams[ii];
+      for( int jj = 0; jj < MaxNumberOfInverseMaps; jj++ ){
+         delete _inverseStreamMaps[ii][jj];
+      }
+   }
+
+}
+
+/** 
+ * Get number of parameter streams
+ * 
+ * @return  number of parameter streams
+ *
+ */
+
+int BrookBonded::getNumberOfParameterStreams( void ) const {
+   return NumberOfParameterStreams;
+}
+
+/** 
+ * Get number of force streams
+ * 
+ * @return  number of force streams
+ *
+ */
+
+int BrookBonded::getNumberOfForceStreams( void ) const {
+   return NumberOfForceStreams;
+}
+
+/** 
+ * Get max number of inverse maps
+ * 
+ * @return  max number of inverse maps
+ *
+ */
+
+int BrookBonded::getMaxInverseMapStreamCount( void ) const {
+   return _maxNumberOfInverseMaps;
+}
+
+/** 
+ * Get width of inverse map streams
+ * 
+ * @return  width of inverse map streams
+ *
+ */
+
+int BrookBonded::getInvMapStreamWidth( void ) const {
+   return _invMapStreamWidth;
+}
+
+/** 
+ * Get LJ 14 scaling parameter
+ * 
+ * @return LJ 14 scaling parameter
+ *
+ */
+
+BrookOpenMMFloat BrookBonded::getLJ_14Scale( void ) const {
+   return _ljScale;
+}
+
+/** 
+ * Get Coulomb factor
+ * 
+ * @return Coulomb factor
+ *
+ */
+
+BrookOpenMMFloat BrookBonded::getCoulombFactor( void ) const {
+   return _coulombFactor;
+}
+
+/** 
+ * Return true if force[index] stream is set
+ *
+ * @param    index into force stream
+ * @return   true if index is valid && force[index] stream is set; else false
+ *
+ */
+
+int BrookBonded::isForceStreamSet( int index ) const {
+   return (index >= 0 && index < getNumberOfForceStreams() && _bondedForceStreams[index]) ? 1 : 0;
+}
+
+/** 
+ * Return string showing if all inverse map streams are set
+ *
+ * @param    index into inverse map stream array
+ *
+ * @return   informative string
+ *
+ */
+
+std::string BrookBonded::checkInverseMapStream( int index ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookBonded::getContents";
+   int ok                                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   if( index < 0 || index >= getNumberOfForceStreams() ){
+      message << "Inverse map index=" << index << " is out of range [0, "  << getNumberOfForceStreams() << ")";
+      return message.str();
+   }
+
+   message << "[ ";
+   for( int ii = 0; ii < getInverseMapStreamCount( index ); ii++ ){
+      message << (_inverseStreamMaps[index][ii] ? 1 : 0) << " ";
+      if( !_inverseStreamMaps[index][ii] ){
+         ok = 0;
+      }
+   }
+   message << " ]";
+   if( !ok ){
+      message << " ERROR!";
+   }
+
+   return message.str();
+   
+}
+
+/** 
+ * Return true if paramsterSet[index] stream is set
+ *
+ * @param    index into parameter stream
+ *
+ * @return   true if index is valid && paramsterSet[index] stream is set; else false
+ *
+ */
+
+int BrookBonded::isParameterStreamSet( int index ) const {
+   return (index >= 0 && index < getNumberOfParameterStreams() && _bondedParameters[index]) ? 1 : 0;
+}
+
+/** 
+ * Validate stream count
+ * 
+ * @return  DefaultReturnValue if count valid; else return ErrorReturnValue
+ *
+ * @exception OpenMMException is thrown if count is invalid
+ *
+ */
+
+int BrookBonded::validateInverseMapStreamCount( int index, int count ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::validateInverseMapStreamCount";
+
+// ---------------------------------------------------------------------------------------
+
+   int valid = DefaultReturnValue;
+
+   if( index == 2 ){
+      if( count > 5 ){
+         valid = ErrorReturnValue;
+      }
+   }
+
+   if( valid == ErrorReturnValue ){
+      std::stringstream message;
+      message << methodName << " input index=" << index << " has invalid count=" << count;
+      throw OpenMMException( message.str() );
+   }
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Get inverse map stream count
+ * 
+ * @return  count
+ *
+ * @exception OpenMMException is thrown if index is out of range [0, getNumberOfForceStreams() ]
+ *
+ */
+
+int BrookBonded::getInverseMapStreamCount( int index ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::getInverseMapStreamCount";
+
+// ---------------------------------------------------------------------------------------
+
+   if( index < 0 || index >= getNumberOfForceStreams() ){
+      std::stringstream message;
+      message << methodName << " input index=" << index << " is out of range [0, " << getNumberOfForceStreams() << " )";
+      throw OpenMMException( message.str() );
+   }
+
+   return _inverseMapStreamCount[index];
+}
+
+/** 
+ * Get bonded atom indices stream
+ * 
+ * @return  atom indices stream
+ *
+ */
+
+BrookFloatStreamImpl* BrookBonded::getAtomIndicesStream( void ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName = "BrookBonded::getAtomIndicesStream";
+
+// ---------------------------------------------------------------------------------------
+
+   return _atomIndicesStream;
+}
+
+/** 
+ * Get array of bonded parameter streams
+ * 
+ * @return  array of bonded parameter streams
+ *
+ */
+
+BrookFloatStreamImpl** BrookBonded::getBondedParameterStreams( void ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName = "BrookBonded::getBondedParameterStreams";
+
+// ---------------------------------------------------------------------------------------
+
+   return _bondedParameters;
+}
+
+/** 
+ * Get array of force streams
+ * 
+ * @return  array  force streams
+ *
+ */
+
+BrookFloatStreamImpl** BrookBonded::getBondedForceStreams( void ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName = "BrookBonded::getBondedForceStreams";
+
+// ---------------------------------------------------------------------------------------
+
+   return _bondedForceStreams;
+}
+
+/** 
+ * Get array of inverse map streams
+ * 
+ * @param index  array index 
+ *
+ * @return  array inverse map streams
+ *
+ */
+
+BrookFloatStreamImpl** BrookBonded::getInverseStreamMapsStreams( int index ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName = "BrookBonded::getInverseStreamMapsStreams";
+
+// ---------------------------------------------------------------------------------------
+
+   // no checking on index -- assume ok -- speed
+
+   return _inverseStreamMaps[index];
+}
+
+/* Gromacs sorts most bonded interactions
+ * We sort those that gromacs forgets
+ * The parameters are all symmetric with respect 
+ * to inversion of order.
+ * 
+ * Also we assume that the same set of indices are not
+ * used with different sets of parameters in gromacs
+ * This is just an assumption for convenience here,
+ * nothing stops us from evaluating such terms in the
+ * kernel.
+ * 
+ * To begin with all interactions have -1 -1 -1 -1
+ * After we fit in all the interactions, we have
+ * to convert the -1's to zeros to avoid indexing
+ * errors on the GPU.
+ * */
+
+/*Flips i,j,k,l to l,k,j,i while correctly shuffling the params */
+
+void BrookBonded::flipQuartet( int ibonded, int *atoms ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName = "BrookBonded::flipQuartet";
+
+// ---------------------------------------------------------------------------------------
+
+   int tmp;
+
+   //For now, simply flip the indices
+   //we're just studying the packing
+
+   tmp = ATOMS( ibonded, 0 );
+   ATOMS( ibonded, 0 ) = ATOMS( ibonded, 3 );
+   ATOMS( ibonded, 3 ) = ATOMS( ibonded, 0 );
+
+   tmp = ATOMS( ibonded, 1 );
+   ATOMS( ibonded, 1 ) = ATOMS( ibonded, 2 );
+   ATOMS( ibonded, 2 ) = ATOMS( ibonded, 1 );
+
+}
+
+int BrookBonded::matchTorsion( int i, int j, int k, int l, int nbondeds, int *atoms ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName = "BrookBonded::matchTorsion";
+
+// ---------------------------------------------------------------------------------------
+
+   int tmp;
+   if( j > k ){
+      //swap into order
+      tmp = i;
+      i = l;
+      l = tmp;
+
+      tmp = j;
+      j = k;
+      k = tmp;
+   }
+   
+   for( int n = 0; n < nbondeds; n++ ){
+      if( i == ATOMS(n, 0) && j == ATOMS(n, 1) && 
+          k == ATOMS(n, 2) && l == ATOMS(n, 3) ){
+         return n;
+      }
+   }
+
+   return ErrorReturnValue;
+}
+
+/* 
+ * We try to match i,j,k against the first 3 in a quartet, then
+ * against the next three. Then we check if there's a free slot
+ * in the fourth component and match the middle two components
+ *
+ * The last bit will not be needed as long as gromacs generates
+ * all dihedrals. But I think there are force fields that don't 
+ * use all dihedrals.
+ *
+ * @return ErrorReturnValue if error; else atom index
+ *
+ **/
+
+int BrookBonded::matchAngle( int i, int j, int k, int nbondeds, 
+                             int *atoms, int *flag ){
+
+// ---------------------------------------------------------------------------------------
+
+   int n;
+   static const std::string methodName = "BrookBonded::matchAngle";
+
+// ---------------------------------------------------------------------------------------
+
+   // validate i > k
+
+   if( i > k ){
+      if( getLog() ){
+         (void) fprintf( getLog(), "%s Invalid triplet %d-%d-%d\n", methodName.c_str(), i, j, k );
+         (void) fflush( getLog() );
+      }
+      std::stringstream message;
+      message << methodName << " invalid triplet: " << i << "-" << j << "-" << k  << std::endl;
+      throw OpenMMException( message.str() );
+      return ErrorReturnValue;
+   }
+   
+   for( n = 0; n < nbondeds; n++ ){
+      if( j == ATOMS(n, 1) ){
+         if(   (i == ATOMS(n, 0) && k == ATOMS(n, 2)) 
+             || (i == ATOMS(n, 2) && k == ATOMS(n, 0)) )  {
+            *flag = 0;
+            return n;
+         }
+      }
+      else if( j == ATOMS(n, 2) ){
+         if( i == ATOMS(n, 1) ){
+            if( ATOMS(n, 3) == -1 || k == ATOMS(n, 3) ){
+               ATOMS(n, 3) = k;
+               *flag = 1;
+               return n;
+            }
+         }
+         else if( k == ATOMS(n, 1) ){
+            if( ATOMS(n, 3) == -1 || i == ATOMS(n, 3) ){
+               ATOMS(n, 3) = i;
+               *flag = 1;
+               return n;
+            }
+         }
+      }
+   }
+   return ErrorReturnValue;
+}
+
+/*flag = 0 means match i-j, 1 means j-k and 2 means k-l
+ * Again, like above, if we have uninitialized slots, we take em
+ *
+ * */
+
+int BrookBonded::matchBond( int i, int j, int nbondeds, int *atoms, int *flag ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::matchBond";
+
+// ---------------------------------------------------------------------------------------
+
+   if( i > j ){
+      if( getLog() ){
+         (void) fprintf( getLog(), "%s invalid bond %d-%d\n", methodName.c_str(), i, j );
+         (void) fflush( getLog() );
+      }
+      std::stringstream message;
+      message << methodName << " invalid bond " << i << "-" << j << std::endl;
+      throw OpenMMException( message.str() );
+      return ErrorReturnValue;
+   }
+
+   for( int n = 0; n < nbondeds; n++ ){
+
+      if( ( i == ATOMS(n, 0) && j == ATOMS(n, 1) ) || ( i == ATOMS(n, 1) && j == ATOMS(n, 0) ) ){
+
+         *flag = 0;
+         return n;
+      }
+      
+      //Try second spot
+
+      if(  ATOMS(n, 2) == -1 ){ //One available spot
+         if( i == ATOMS(n, 1) ){
+            ATOMS(n, 2) = j;
+            *flag = 1;
+            return n;
+         }
+         if( j == ATOMS(n, 1) ){
+            ATOMS(n, 2) = i;
+            *flag = 1;
+            return n;
+         }
+      }
+      else if(  ( i == ATOMS(n, 1) && j == ATOMS(n, 2) ) ||( i == ATOMS(n, 2) && j == ATOMS(n, 1) ) ){
+         *flag = 1;
+         return n;
+      }
+      
+      //Try third spot
+      if( ATOMS(n, 3) == -1 ){
+         if( i == ATOMS(n, 2) ){
+            ATOMS(n, 3) = j;
+            *flag = 2;
+            return n;
+         }
+         if( j == ATOMS(n, 2) ){
+            ATOMS(n, 3) = i;
+            *flag = 2;
+            return n;
+         }
+      } else if(  ( i == ATOMS(n, 2) && j == ATOMS(n, 3) ) ||( i == ATOMS(n, 3) && j == ATOMS(n, 2) ) ){
+         *flag = 2;
+         return n;
+      }
+   }
+
+   return ErrorReturnValue;
+}
+
+int BrookBonded::matchPair( int i, int j, int nbondeds, int *atoms ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName = "BrookBonded::matchPair";
+
+// ---------------------------------------------------------------------------------------
+
+   for( int n = 0; n < nbondeds; n++ ){
+
+      //If there is an empty slot available:
+
+      if( ATOMS(n, 0) == -1 ){
+
+         //If one of i,j matches the l atom
+
+         if( ATOMS(n, 3) == i ){
+            ATOMS(n, 0) = j;
+            return n;
+         }
+         if( ATOMS(n, 3) == j ){
+            ATOMS(n, 0) = i;
+            return n;
+         }
+      }
+
+      //If the l-atom is available
+      if( ATOMS(n, 3) == -1 ){
+         if( ATOMS(n, 0) == i ){
+            ATOMS(n, 3) = j;
+            return n;
+         }
+         if( ATOMS(n, 0) == j ){
+            ATOMS(n, 3) = i;
+            return n;
+         }
+      }
+
+      //Both are unavailable, both much match
+      if(   ( i == ATOMS(n, 0) && j == ATOMS(n, 3) ) 
+          || ( i == ATOMS(n, 3) && j == ATOMS(n, 0) ) ){
+         return n;
+      }
+   }
+   return ErrorReturnValue;
+}
+
+/** 
+ * Setup Ryckaert-Bellemans parameters/atom indices
+ * 
+ * @param nbondeds                  number of bonded entries
+ * @param atoms                     array of atom indices
+ * @param params                    arrays of bond parameters
+ * @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 log                       log reference
+ *
+ * @return nonzero value if error
+ *
+ */
+
+int BrookBonded::addRBDihedrals( int *nbondeds, int *atoms, float *params[], 
+                                 const vector<vector<int> >& rbTorsionIndices, 
+                                 const vector<vector<double> >& rbTorsionParameters ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::BrookBonded";
+   static const int debug              = 0;
+
+// ---------------------------------------------------------------------------------------
+
+   if( debug && getLog() ){
+      (void) fprintf( getLog(), "%s number of bonds=%d\n", methodName.c_str(), rbTorsionIndices.size() );
+   }
+
+   for( unsigned int ii = 0; ii < rbTorsionIndices.size(); ii++ ){
+
+      vector<int> atomsIndices      = 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 ibonded = matchTorsion( i, j, k, l, *nbondeds, atoms );
+      
+      if( ibonded < 0 ){
+         ibonded             = *nbondeds;
+         ATOMS( ibonded, 0 ) = i;
+         ATOMS( ibonded, 1 ) = j;
+         ATOMS( ibonded, 2 ) = k;
+         ATOMS( ibonded, 3 ) = l;
+         (*nbondeds)++;
+      }
+
+      index                   = 0;
+      PARAMS( ibonded, 0, 0 ) = (BrookOpenMMFloat) rbParameters[index++];
+      PARAMS( ibonded, 0, 1 ) = (BrookOpenMMFloat) rbParameters[index++];
+      PARAMS( ibonded, 0, 2 ) = (BrookOpenMMFloat) rbParameters[index++];
+      PARAMS( ibonded, 0, 3 ) = (BrookOpenMMFloat) rbParameters[index++];
+      PARAMS( ibonded, 1, 0 ) = (BrookOpenMMFloat) rbParameters[index++];
+
+      if( debug && getLog() ){
+         (void) fprintf( getLog(), "   %d [%d %d %d %d] %.3e %.3e %.3e %.3e\n", ibonded, i, j, k, l,
+                         rbParameters[0], rbParameters[1],
+                         rbParameters[2], rbParameters[3], rbParameters[4] );
+      }
+   }
+
+   return DefaultReturnValue;
+}
+
+/**
+ * Setup periodic torsion parameters/atom indices
+ * 
+ * @param nbondeds                  number of bonded entries
+ * @param atoms                     array of atom indices
+ * @param params                    arrays of bond parameters
+ * @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 log                       log reference
+ *
+ * @return nonzero value if error
+ *
+ */
+
+int BrookBonded::addPDihedrals( int *nbondeds, int *atoms, BrookOpenMMFloat* params[],
+                                const vector<vector<int> >& periodicTorsionIndices, 
+                                const vector<vector<double> >& periodicTorsionParameters ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::addPDihedrals";
+   static const int debug              = 0;
+
+// ---------------------------------------------------------------------------------------
+
+   if( debug && getLog() ){
+      (void) fprintf( getLog(), "%s npdih=%d\n", methodName.c_str(), periodicTorsionIndices.size() );
+   }
+
+   for( unsigned int ii = 0; ii < periodicTorsionIndices.size(); ii++ ){
+
+      vector<int> atomsIndices      = 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 ibonded = matchTorsion( i, j, k, l, *nbondeds, atoms );
+      
+      if( ibonded < 0 ){
+         ibonded = *nbondeds;
+         ATOMS( ibonded, 0 ) = i;
+         ATOMS( ibonded, 1 ) = j;
+         ATOMS( ibonded, 2 ) = k;
+         ATOMS( ibonded, 3 ) = l;
+         (*nbondeds)++;
+      }
+      PARAMS( ibonded, 1, 1 ) = (BrookOpenMMFloat) pTParameters[0];
+      PARAMS( ibonded, 1, 2 ) = (BrookOpenMMFloat) pTParameters[1];
+      PARAMS( ibonded, 1, 3 ) = (BrookOpenMMFloat) pTParameters[2];
+
+      if( debug && getLog() ){
+         (void) fprintf( getLog(), "   %d [%d %d %d %d] %.3e %.3e %.3e\n", ibonded, i, j, k, l,
+                         pTParameters[0], pTParameters[1], pTParameters[2] );
+      }
+   }
+
+   return DefaultReturnValue;
+}
+
+/**
+ * Setup angle bond parameters/atom indices
+ * 
+ * @param nbondeds                  number of bonded entries
+ * @param atoms                     array of atom indices
+ * @param params                    arrays of bond parameters
+ * @param angleIndices              the angle bond atom indices
+ * @param angleParameters           the angle parameters (angle in radians, force constant)
+ * @param log                       log reference
+ *
+ * @return nonzero value if error
+ *
+ */
+
+int BrookBonded::addAngles( int *nbondeds, int *atoms, float *params[], const std::vector<std::vector<int> >& angleIndices,
+                            const std::vector<std::vector<double> >& angleParameters ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::addAngles";
+   static const int debug              = 0;
+
+// ---------------------------------------------------------------------------------------
+
+   if( debug && getLog() ){
+      (void) fprintf( getLog(), "%s nang=%d\n", methodName.c_str(), angleIndices.size() );
+   }
+
+   // loop over bonds
+
+   for( unsigned int ii = 0; ii < angleIndices.size(); ii++ ){
+
+      vector<int> atomsIndices       = angleIndices[ii];
+      vector<double> angParameters   = angleParameters[ii];
+
+      int index = 0;
+      int i     = atomsIndices[index++];
+      int j     = atomsIndices[index++];
+      int k     = atomsIndices[index++];
+ 
+      int flag;
+      int ibonded = matchAngle( i, j, k, *nbondeds, atoms, &flag );
+      if( ibonded < 0 ){
+         ibonded = *nbondeds;
+         ATOMS( ibonded, 0 ) = i;
+         ATOMS( ibonded, 1 ) = j;
+         ATOMS( ibonded, 2 ) = k;
+         flag = 0;
+         (*nbondeds)++;
+      }
+      PARAMS( ibonded, 2, flag*2 )     = (BrookOpenMMFloat) angParameters[0];
+      PARAMS( ibonded, 2, flag*2 + 1 ) = (BrookOpenMMFloat) angParameters[1];
+      
+      if( debug && getLog() ){
+         (void) fprintf( getLog(), "   %d [%d %d %d ] %.6e %.6e\n", ibonded, i, j, k,
+                         angParameters[0], angParameters[1] );
+      }
+   }
+
+   return DefaultReturnValue;
+}
+
+/**
+ * Setup harmonic bond parameters/atom indices
+ * 
+ * @param nbondeds                  number of bonded entries
+ * @param atoms                     array of atom indices
+ * @param params                    arrays of bond parameters
+ * @param bondIndices               two harmonic bond atom indices
+ * @param bondParameters            the force parameters (distance, k)
+ * @param log                       log reference
+ *
+ * @return nonzero value if error
+ *
+ */
+
+int BrookBonded::addBonds( int *nbondeds, int *atoms, float *params[], const vector<vector<int> >& bondIndices,
+                           const vector<vector<double> >& bondParameters ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::addBonds";
+   static const int debug              = 0;
+
+// ---------------------------------------------------------------------------------------
+
+   if( debug && getLog() ){
+      (void) fprintf( getLog(), "%s nbonds=%d\n", methodName.c_str(), bondIndices.size() );
+   }
+
+   // loop over bonds
+
+   for( unsigned int ii = 0; ii < bondIndices.size(); ii++ ){
+
+      vector<int> atomsIndices       = bondIndices[ii];
+      vector<double> bndParameters   = bondParameters[ii];
+
+      int index = 0;
+      int i     = atomsIndices[index++];
+      int j     = atomsIndices[index++];
+
+      int flag;
+      int ibonded = matchBond( i, j, *nbondeds, atoms, &flag );
+      if( ibonded < 0 ){
+         ibonded = *nbondeds;
+         ATOMS( ibonded, 0 ) = i;
+         ATOMS( ibonded, 1 ) = j;
+         flag = 0;
+         (*nbondeds)++;
+      }
+
+      if( flag < 2 ){
+         PARAMS( ibonded, 3, flag*2 )     = (BrookOpenMMFloat) bndParameters[0];
+         PARAMS( ibonded, 3, flag*2 + 1 ) = (BrookOpenMMFloat) bndParameters[1];
+      } else {
+         PARAMS( ibonded, 4, 0 )          = (BrookOpenMMFloat) bndParameters[0];
+         PARAMS( ibonded, 4, 1 )          = (BrookOpenMMFloat) bndParameters[1];
+      }
+
+      if( debug && getLog() ){
+         (void) fprintf( getLog(), "   %d [%d %d ] flag=%d %.3e %.3e\n", ibonded, i, j, flag,
+                         bndParameters[0], bndParameters[1] );
+      }
+   }
+
+   return DefaultReturnValue;
+}
+
+/**
+ * Setup LJ/Coulomb 1-4 parameters/atom indices
+ * 
+ * @param nbondeds                  number of bonded entries
+ * @param atoms                     array of atom 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
+ *                                  they form a bonded 1-4 pair
+ * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+ * @param lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
+ * @param log                       log reference
+ *
+ * @return nonzero value if error
+ *
+ */
+
+int BrookBonded::addPairs( int *nbondeds, int *atoms, BrookOpenMMFloat* params[],
+                           BrookOpenMMFloat* charges,
+                           const std::vector<std::vector<int> >& bonded14Indices,
+                           const std::vector<std::vector<double> >& nonbondedParameters,
+                           double lj14Scale, double coulombScale ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::BrookBonded";
+   static const double oneSixth        = 1.0/6.0;
+   static const int debug              = 0;
+
+// ---------------------------------------------------------------------------------------
+
+   if( debug && getLog() ){
+      (void) fprintf( getLog(), "%s npairs=%d\n", methodName.c_str(), bonded14Indices.size() );
+   }
+
+   for( unsigned int ii = 0; ii < bonded14Indices.size(); ii++ ){
+
+      std::vector<int> atomsIndices        = bonded14Indices[ii];
+
+      int index = 0;
+      int i     = atomsIndices[index++];
+      int j     = atomsIndices[index++];
+
+      int ibonded = matchPair( i, j, *nbondeds, atoms );
+      if( ibonded < 0 ){
+         ibonded = *nbondeds;
+         ATOMS(ibonded, 0) = i;
+         ATOMS(ibonded, 3) = j;
+         (*nbondeds)++;
+      }
+
+      vector<double> iParameters  = nonbondedParameters[i];
+      vector<double> jParameters  = nonbondedParameters[j];
+
+      double c6                   = iParameters[0] + jParameters[0];
+      double c12                  = lj14Scale*(iParameters[1] * jParameters[1]);
+
+      double sig, eps;
+      if( c12 != 0.0 ){
+         eps = c6*c6/c12;
+         sig = pow( c12/c6, oneSixth );
+      } else {
+         eps = 0.0;
+         sig = 1.0;
+      }
+
+      PARAMS( ibonded, 4, 2 ) = (BrookOpenMMFloat) sig;
+      PARAMS( ibonded, 4, 3 ) = (BrookOpenMMFloat) eps;
+
+      // a little wasteful, but ...
+
+      charges[i] = (BrookOpenMMFloat) iParameters[2];
+      charges[j] = (BrookOpenMMFloat) jParameters[2];
+
+      if( debug ){
+         (void) fprintf( getLog(), "\n %d [%d %d ] %.3e %.3e", ibonded, i, j, sig, eps );
+      }
+   }
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * 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 platform                  BrookPlatform reference
+ * @param log                       log file reference (optional)
+ *
+ * @return nonzero value if error
+ *
+ */
+
+int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookPlatform& brookPlatform ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::loadInvMaps";
+
+// ---------------------------------------------------------------------------------------
+
+   // get atom stream size
+
+   const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (brookPlatform.getDefaultStreamFactory() );
+   BrookStreamInfo* brookStreamInfo             = brookStreamFactory.getBrookStreamInfo( BrookStreamFactory::AtomPositions );
+   int atomStreamWidth                          = brookStreamInfo->getStreamWidth();
+   int atomStreamSize                           = brookPlatform.getStreamSize( getNumberOfAtoms(), atomStreamWidth, NULL );
+   _invMapStreamWidth                           = atomStreamWidth;
+   
+   // allocate temp memory
+
+   float4** invmaps = new float4*[getMaxInverseMapStreamCount()];
+   float* block     = new float[4*getMaxInverseMapStreamCount()*atomStreamSize];
+   for( int ii = 0; ii < getMaxInverseMapStreamCount(); ii++ ){
+      invmaps[ii]  = (float4*) block;
+      block       += 4*atomStreamSize;
+   }
+   int* counts = new int[atomStreamSize];
+
+   // get inverse maps and load into streams
+
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      gpuCalcInvMap( ii, 4, nbondeds, natoms, atoms, getInverseMapStreamCount( ii ), counts, invmaps, &(_inverseMapStreamCount[ii]) );
+      validateInverseMapStreamCount( ii, _inverseMapStreamCount[ii] ); 
+      for( int jj = 0; jj < _inverseMapStreamCount[ii]; jj++ ){
+         StreamImpl* inverseStreamMaps = brookStreamFactory.createStreamImpl( BrookStreamFactory::BondedInverseMapStreams, atomStreamSize, Stream::Float4, brookPlatform );
+         _inverseStreamMaps[ii][jj]    = dynamic_cast<BrookFloatStreamImpl*> ( inverseStreamMaps );
+         _inverseStreamMaps[ii][jj]->loadFromArray( invmaps[jj] );
+      }
+   }
+
+   if( getLog() ){
+      (void) fprintf( getLog(), "%s done\n", methodName.c_str() );
+      (void) fflush( getLog() );
+      //gpuPrintInvMaps( bp->nimaps, natoms, counts, invmaps, gpu->log );
+   }
+
+   // free memory
+
+   delete counts;
+   delete invmaps[0];
+   delete invmaps;
+
+   return DefaultReturnValue;
+}
+
+/* 
+ * 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 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 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 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 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 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
+ * @param platform                     Brook platform reference
+ *
+ * @return always 1
+ *
+ * We go through the interactions in the following order
+ * torsions
+ * angles
+ * 14 interactions
+ * bonds
+ * For each new interaction, we try to fit it in with 
+ * those already in. This may not necessarily result in
+ * the optimal fit, but should not be too bad. 
+ * */
+
+int BrookBonded::setup( int numberOfAtoms,
+                        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,
+                        const vector<vector<int> >& rbTorsionIndices,       const vector<vector<double> >& rbTorsionParameters,
+                        const vector<vector<int> >& bonded14Indices,        const vector<vector<double> >& nonbondedParameters,
+                        double lj14Scale, double coulombScale,  const BrookPlatform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "BrookBonded::setup";
+
+// ---------------------------------------------------------------------------------------
+
+   _numberOfAtoms = numberOfAtoms;
+
+   // check that atom 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();
+      throw OpenMMException( message.str() );
+   } else if( getLog() ){
+      (void) fprintf( getLog(), "%s harmonic bonds=%d\n", methodName.c_str(), bondIndices.size() );
+      (void) fflush( getLog() );
+   }
+
+   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();
+      throw OpenMMException( message.str() );
+   } else if( getLog() ){
+      (void) fprintf( getLog(), "%s angle bonds=%d\n", methodName.c_str(), angleIndices.size() );
+      (void) fflush( getLog() );
+   }
+
+   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();
+      throw OpenMMException( message.str() );
+   } else if( getLog() ){
+      (void) fprintf( getLog(), "%s periodicTorsion bonds=%d\n", methodName.c_str(), periodicTorsionIndices.size() );
+      (void) fflush( getLog() );
+   }
+
+   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();
+      throw OpenMMException( message.str() );
+   } else if( 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 ){
+      std::stringstream message;
+      message << methodName << " number atoms=" << numberOfAtoms << " does not equal number of nb parameter entries=" << nonbondedParameters.size();
+      throw OpenMMException( message.str() );
+   } else if( getLog() ){
+      (void) fprintf( getLog(), "%s LJ 14 ixns=%d\n", methodName.c_str(), bonded14Indices.size() );
+      (void) fflush( getLog() );
+   }
+
+   // allocate temp memory
+
+   int maxBonds = 10*numberOfAtoms;
+   int* atoms   = new int[5*maxBonds];
+
+   BrookOpenMMFloat** params = new BrookOpenMMFloat*[getNumberOfParameterStreams()];
+   for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
+      params[ii] = new BrookOpenMMFloat[4*maxBonds];
+   }
+   BrookOpenMMFloat* charges = new BrookOpenMMFloat[numberOfAtoms];
+
+   // Initialize all atom indices to -1 to indicate empty slots
+   // All parameters must be initialized to values that will 
+   // produce zero for the corresponding force. 
+
+   for( int ii = 0; ii < maxBonds; ii++ ){
+
+      ATOMS( ii, 0 ) = -1;
+      ATOMS( ii, 1 ) = -1;
+      ATOMS( ii, 2 ) = -1;
+      ATOMS( ii, 3 ) = -1;
+   
+      for( int jj = 0; jj < getNumberOfParameterStreams(); jj++ ){
+         PARAMS( ii, jj, 0 ) = 0.0; 
+         PARAMS( ii, jj, 1 ) = 0.0; 
+         PARAMS( ii, jj, 2 ) = 0.0; 
+         PARAMS( ii, jj, 3 ) = 0.0; 
+      }
+
+      //Set sigma negative to indicate no pair
+      PARAMS( ii, 4, 2 ) = -1.0;
+   }
+
+   // 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            );
+
+   addPairs( &nbondeds, atoms, params, charges, bonded14Indices, nonbondedParameters, lj14Scale, coulombScale );
+
+   // check that number of bonds not too large for memory allocated
+
+   if( nbondeds >= maxBonds ){
+      std::stringstream message;
+      message << methodName << " number of bonds=" << nbondeds << " is greater than maxBonds=" << maxBonds;
+      throw OpenMMException( message.str() );
+   } else if( getLog() ){
+      (void) fprintf( getLog(), "%s atoms=%d number of bonds=%d maxBonds=%d\n", methodName.c_str(), numberOfAtoms, nbondeds, maxBonds );
+      (void) fflush( getLog() );
+   }
+
+   // get factory
+
+   const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (platform.getDefaultStreamFactory());
+
+   // build streams
+
+   StreamImpl* atomIndicesStream    = brookStreamFactory.createStreamImpl( BrookStreamFactory::BondedAtomIndicesStream, nbondeds, Stream::Float4, platform );
+   _atomIndicesStream               = dynamic_cast<BrookFloatStreamImpl*> (atomIndicesStream); 
+   _atomIndicesStream->loadFromArray( atoms, Stream::Integer ); 
+
+   StreamImpl* chargeStream         = brookStreamFactory.createStreamImpl( BrookStreamFactory::BondedChargeStream, numberOfAtoms, Stream::Float, platform );
+   _chargeStream                    = dynamic_cast<BrookFloatStreamImpl*> (chargeStream); 
+   _chargeStream->loadFromArray( charges ); 
+
+   for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
+      StreamImpl* bondedParameters  = brookStreamFactory.createStreamImpl( BrookStreamFactory::BondedParametersStream, nbondeds, Stream::Float4, platform );
+      _bondedParameters[ii]         = dynamic_cast<BrookFloatStreamImpl*> (bondedParameters); 
+      _bondedParameters[ii]->loadFromArray( params[ii] );
+   }
+
+   // debug stuff
+
+   if( 1 && getLog() ){
+
+      BrookFloatStreamImpl* atomIndicesStream = dynamic_cast<BrookFloatStreamImpl*> (_atomIndicesStream); 
+      (void) fprintf( getLog(), "%s nbondeds=%d strDim [%d %d ] sz=%d\n", methodName.c_str(), nbondeds,
+                      atomIndicesStream->getStreamWidth(),
+                      atomIndicesStream->getStreamHeight(), 
+                      atomIndicesStream->getStreamSize() );
+
+      int kIndex = 0;
+      int jIndex = 1;
+      int iIndex = 2;
+      int lIndex = 3;
+      int mIndex = 4;
+      
+      /*
+       * float4 parm0   (rbc[1], rbc[2], rbc[3], rbc[4])
+       * float4 parm1   (rbc[5], cp, phi, mult) latter three are for pdih
+       * float4 parm2   (theta, k, theta, k ) angles i-j-k and j-k-l
+       * float4 parm3   (x0, k, x0, k) bonds i-j, j-k
+       * float4 parm4   (x0, k, sig14, eps14) bond k-l, i-l
+       */
+      
+      (void) fprintf( getLog(), "\nParams\n" );
+      int index = 0; 
+      for( int ii = 0; ii < 4*nbondeds; ii += 4, index++ ){
+         //  #define PARAMS(X,Y,Z) (params[(Y)][4*(X) + Z])
+         (void) fprintf( getLog(), "\n%4d   [%4d %4d %4d %4d]\n"
+                              "     rb[%10.6f %10.6f %10.6f %10.6f %10.6f]\n"
+                              "    phi[%10.6f %6.1f %5.1f]\n"
+                              "    ang[%6.3f %8.2f] [%6.3f %8.2f]\n"
+                              "      h[%8.5f %14.6e] [%8.5f %14.6e] [%8.5f %14.6e]\n"
+                              "     14[%15.6e %15.6e]\n", index,
+                         ATOMS(index, 0), ATOMS(index, 1), ATOMS(index, 2), ATOMS(index, 3),
+                         params[kIndex][ii], params[kIndex][ii+1], params[kIndex][ii+2], params[kIndex][ii+3], params[jIndex][ii],
+                         params[jIndex][ii+1], params[jIndex][ii+2], params[jIndex][ii+3],
+                         params[iIndex][ii], params[iIndex][ii+1], params[iIndex][ii+2], params[iIndex][ii+3], 
+                         params[lIndex][ii], params[lIndex][ii+1], params[lIndex][ii+2], params[lIndex][ii+3],
+                         params[mIndex][ii], params[mIndex][ii+1], params[mIndex][ii+2], params[mIndex][ii+3] );
+      }
+   }
+
+   // load inverse maps to streams
+
+   loadInvMaps( nbondeds, getNumberOfAtoms(), atoms, platform );
+   
+   // free memory
+
+   for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
+      delete[] params[ii];
+   }
+   delete[] params;
+   delete[] atoms;
+   delete[] charges;
+
+   // set the fudge factors
+
+   _ljScale        = (BrookOpenMMFloat) lj14Scale;
+   _coulombFactor  = (BrookOpenMMFloat) coulombScale;
+
+   // initialize output streams
+
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      StreamImpl* bondedForceStreams = brookStreamFactory.createStreamImpl( BrookStreamFactory::UnrolledForceStream, nbondeds, Stream::Float3, platform );
+      _bondedForceStreams[ii]        = dynamic_cast<BrookFloatStreamImpl*> (bondedForceStreams); 
+   }
+
+   return DefaultReturnValue;
+}
+
+/* 
+ * Get contents of object
+ *
+ *
+ * @param level   level of dump
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string BrookBonded::getContents( int level ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookBonded::getContents";
+
+   static const unsigned int MAX_LINE_CHARS = 256;
+   char value[MAX_LINE_CHARS];
+   static const char* Set                   = "Set";
+   static const char* NotSet                = "Not set";
+
+// ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   std::string tab   = "   ";
+
+#ifdef WIN32
+#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
+#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, "%d", getNumberOfAtoms() );
+   message << _getLine( tab, "Number of atoms:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5f", getLJ_14Scale() );
+   message << _getLine( tab, "LJ 14 scaling:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%.5f", getCoulombFactor() );
+   message << _getLine( tab, "Coulomb factor:", value ); 
+
+/*
+   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamWidth() );
+   message << _getLine( tab, "Atom stream width:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamHeight() );
+   message << _getLine( tab, "Atom stream height:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamSize() );
+   message << _getLine( tab, "Atom stream size:", value ); 
+*/
+
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfParameterStreams() );
+   message << _getLine( tab, "Number of parameter streams:", value ); 
+   for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
+      char description[256];
+      (void) LOCAL_SPRINTF( description, "Parameter stream %d", ii );
+      message << _getLine( tab, description,  (isParameterStreamSet(ii) ? Set : NotSet) ); 
+   }
+ 
+   (void) LOCAL_SPRINTF( value, "%d", getMaxInverseMapStreamCount() )
+   message << _getLine( tab, "Max inverseMap count:", value ); 
+
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      char description[256];
+      (void) LOCAL_SPRINTF( description, "Inverse map count %d", ii );
+      
+      std::string checkInverseMap = checkInverseMapStream( ii );
+      (void) LOCAL_2_SPRINTF( value, "%d %s", getInverseMapStreamCount( ii ), checkInverseMap.c_str() );
+      message << _getLine( tab, description, value ); 
+   }
+
+   message << _getLine( tab, "Log:",                 (getLog()               ? Set : NotSet) ); 
+   message << _getLine( tab, "Atom indices stream:", (getAtomIndicesStream() ? Set : NotSet) ); 
+   //message << _getLine( tab, "Charge stream:",        (getChargeStream()      ? Set : NotSet) ); 
+ 
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfForceStreams() );
+   message << _getLine( tab, "Number of force streams:", value ); 
+
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      char description[256];
+      (void) LOCAL_SPRINTF( description, "Force stream %d", ii );
+      message << _getLine( tab, description,  (isForceStreamSet(ii) ? Set : NotSet) ); 
+   }
+ 
+#undef LOCAL_SPRINTF
+#undef LOCAL_2_SPRINTF
+
+   return message.str();
+}
+

platforms/brook/src/BrookBonded.h

+#ifndef BrookBonded_H_
+#define BrookBonded_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include <vector>
+
+#include "BrookFloatStreamImpl.h"
+#include "BrookIntStreamImpl.h"
+#include "BrookPlatform.h"
+#include "BrookCommon.h"
+#include "OpenMMContext.h"
+
+namespace OpenMM {
+
+/**
+ * This kernel is invoked by StandardMMForceField to calculate the forces acting on the system.
+ */
+class BrookBonded : public BrookCommon {
+
+   public:
+  
+      static const int DefaultReturnValue = 0;
+      static const int ErrorReturnValue   = -1; 
+
+      BrookBonded( );
+  
+      ~BrookBonded();
+  
+      /**
+       * Initialize the kernel, setting up the values of all the force field parameters.
+       * 
+       * @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 nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+       * @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
+       *
+       * @return nonzero value if error
+       *
+       */
+
+      int setup( int numberOfAtoms, 
+                 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,        const std::vector<std::vector<double> >& nonbondedParameters,
+                 double lj14Scale, double coulomb14Scale, const BrookPlatform& platform );
+
+      /**
+       * Return number of parameter streams
+       * 
+       * @return number of parameter streams
+       *
+       */
+
+      int getNumberOfParameterStreams( void ) const; 
+
+      /**
+       * Return number of force streams
+       * 
+       * @return number of force streams
+       *
+       */
+
+      int getNumberOfForceStreams( void ) const; 
+
+      /**
+       * Return stream count for specifed index (i, j, k, l)
+       * 
+       * @return stream count for specifed index
+       *
+       * @throw throws OpenMMException if index out of range
+       */
+
+      int getInverseMapStreamCount( int index ) const; 
+
+      /**
+       * Return max stream count
+       * 
+       * @return max stream count
+       */
+
+      int getMaxInverseMapStreamCount( void ) const; 
+
+      /**
+       * Return Brook stream handle
+       * 
+       * @return 
+       */
+
+      BrookFloatStreamImpl* getBrookAtomIndices( void ) const; 
+
+      /**
+       * Get LJ 14 scale factor
+       *
+       * @return  LJ 14 scaling (fudge) factor
+       * 
+       */
+
+       BrookOpenMMFloat getLJ_14Scale( void ) const;
+
+      /** 
+       * Get Coulomb factor
+       * 
+       * @return Coulomb factor
+       *
+       */
+      
+      BrookOpenMMFloat getCoulombFactor( void ) const;
+
+      /** 
+       * Get bonded atom indices stream
+       * 
+       * @return  atom indices stream
+       *
+       */
+      
+      BrookFloatStreamImpl* getAtomIndicesStream( void ) const;
+      
+      /** 
+       * Get array of bonded parameter streams
+       * 
+       * @return  array of bonded parameter streams
+       *
+       */
+      
+      BrookFloatStreamImpl** getBondedParameterStreams( void );
+      
+      /** 
+       * Get array of force streams
+       * 
+       * @return  array
+       *
+       */
+      
+      BrookFloatStreamImpl** getBondedForceStreams( void );
+      
+      /** 
+       * Get array of inverse map streams
+       * 
+       * @param index  array index 
+       *
+       * @return  array inverse map streams
+       *
+       */
+      
+      BrookFloatStreamImpl** getInverseStreamMapsStreams( int index );
+      
+      /** 
+       * Return true if force[index] stream is set
+       *
+       * @param    index into force stream
+       * @return   true if index is valid && force[index] stream is set; else false
+       *
+       */
+      
+      int isForceStreamSet( int index ) const;
+       
+      /** 
+       * Return true if paramsterSet[index] stream is set
+       *
+       * @param    index into parameter stream
+       *
+       * @return   true if index is valid && paramsterSet[index] stream is set; else false
+       *
+       */
+      
+      int isParameterStreamSet( int index ) const;
+          
+      /** 
+       * Return string showing if all inverse map streams are set
+       *
+       * @param    index into inverse map stream array
+       *
+       * @return   informative string
+       *
+       */
+      
+      std::string checkInverseMapStream( int index ) const;
+
+      /* 
+       * Get contents of object
+       *
+       *
+       * @param level   level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string getContents( int level ) const;
+
+   private:
+   
+      static const int NumberOfParameterStreams = 5;
+      static const int NumberOfForceStreams     = 4;
+      static const int MaxNumberOfInverseMaps   = 9;
+
+      enum { StreamI, StreamJ, StreamK, StreamL, StreamMax };
+
+      // inverse map stream width
+
+      int _invMapStreamWidth;
+
+      // actual max number of inverse maps 
+
+      int _maxNumberOfInverseMaps;
+
+      // scale factors for 1-4 ixn
+
+      BrookOpenMMFloat _ljScale;
+      BrookOpenMMFloat _coulombFactor;
+
+      // streams
+
+      BrookFloatStreamImpl* _atomIndicesStream;
+      BrookFloatStreamImpl* _bondedParameters[NumberOfParameterStreams];
+      BrookFloatStreamImpl* _bondedForceStreams[NumberOfForceStreams];
+      BrookFloatStreamImpl* _chargeStream;
+      BrookFloatStreamImpl* _inverseStreamMaps[NumberOfForceStreams][MaxNumberOfInverseMaps];
+
+      int _maxInverseMapStreamCount[NumberOfForceStreams];
+      int _inverseMapStreamCount[NumberOfForceStreams];
+
+      // 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 );
+
+      /**
+       * Setup Ryckaert-Bellemans parameters/atom indices
+       * 
+       * @param nbondeds                  number of bonded entries
+       * @param atoms                     array of atom indices
+       * @param params                    arrays of bond parameters
+       * @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
+       *
+       * @return nonzero value if error
+       *
+       */
+
+      int addRBDihedrals( int *nbondeds, int *atoms, BrookOpenMMFloat* params[], const std::vector<std::vector<int> >& rbTorsionIndices,  
+                          const std::vector<std::vector<double> >& rbTorsionParameters );
+
+      /**
+       * Setup periodic torsion parameters/atom indices
+       * 
+       * @param nbondeds                  number of bonded entries
+       * @param atoms                     array of atom indices
+       * @param params                    arrays of bond parameters
+       * @param periodicTorsionIndices    the four atoms 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,  
+                          const std::vector<std::vector<double> >& periodicTorsionParameters );
+
+      /**
+       * Setup angle bond parameters/atom indices
+       * 
+       * @param nbondeds                  number of bonded entries
+       * @param atoms                     array of atom indices
+       * @param params                    arrays of bond parameters
+       * @param angleIndices              the angle bond atom 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,
+                     const std::vector<std::vector<double> >& angleParameters );
+
+      /**
+       * Setup harmonic bond parameters/atom indices
+       * 
+       * @param nbondeds                  number of bonded entries
+       * @param atoms                     array of atom indices
+       * @param params                    arrays of bond parameters
+       * @param bondIndices               two harmonic bond atom 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,
+                    const std::vector<std::vector<double> >& bondParameters );
+
+      /**
+       * Setup LJ/Coulomb 1-4 parameters/atom indices
+       * 
+       * @param nbondeds                  number of bonded entries
+       * @param atoms                     array of atom 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
+       *                                  they form a bonded 1-4 pair
+       * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+       * @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,
+                    const std::vector<std::vector<int> >& bonded14Indices, const std::vector<std::vector<double> >& nonbondedParameters,
+                    double lj14Scale, double coulombScale );
+      
+      /**
+       * 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 platform                  BrookPlatform reference
+       * @param log                       log file reference (optional)
+       *
+       * @return nonzero value if error
+       *
+       */
+
+      int loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookPlatform& platform );
+      
+      /**
+       * Get inverse amp stream width
+       * 
+       * @return stream width
+       *
+       */
+
+      int getInvMapStreamWidth( void ) const;
+      
+      /**
+       * Validate inverse map count
+       * 
+       * @param index index to check (1-4) 
+       * @param count count for index
+       *
+       * @return -1 if count invalid
+       *
+       * @tthrow OpenMMException exeception if count invalid
+       *
+       */
+
+      int validateInverseMapStreamCount( int index, int count ) const;
+      
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_BROOKKERNELS_H_*/

platforms/brook/src/BrookNonBonded.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Mark Friedrichs                                                   *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <sstream>
+#include "BrookNonBonded.h"
+#include "BrookPlatform.h"
+#include "BrookStreamFactory.h"
+#include "OpenMMException.h"
+#include "gpu/invmap.h"
+#include "gpu/kforce.h"
+
+using namespace OpenMM;
+using namespace std;
+
+/** 
+ * Constructor
+ * 
+ */
+
+BrookNonBonded::BrookNonBonded(  ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::BrookNonBonded";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   _atomSizeCeiling           = -1;
+   _outerUnroll               = 4;
+   _innerUnroll               = 4;
+
+   _partialForceStreamWidth   = 64;
+   _partialForceStreamHeight  = -1;
+   _partialForceStreamSize    = -1;
+
+   _duplicationFactor         = 4;
+
+   _exclusionStreamWidth      = -1;
+   _exclusionStreamHeight     = -1;
+   _exclusionStreamSize       = -1;
+
+   _jStreamWidth              = -1;
+   _jStreamHeight             = -1;
+   _jStreamSize               = -1;
+
+   _exclusionStream           = NULL;
+   _sigmaStream               = NULL;
+   _epsilonStream             = NULL;
+
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      _nonbondedForceStreams[ii] = NULL;
+   }
+
+}   
+ 
+/** 
+ * Destructor
+ * 
+ */
+
+BrookNonBonded::~BrookNonBonded( ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookNonBonded::~BrookNonBonded";
+   //static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   delete _exclusionStream;
+   delete _sigmaStream;
+   delete _epsilonStream;
+
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      delete _nonbondedForceStreams[ii];
+   }
+
+}
+
+/** 
+ * Get number of force streams
+ * 
+ * @return  number of force streams (fixed value)
+ *
+ */
+
+int BrookNonBonded::getNumberOfForceStreams( void ) const {
+   return NumberOfForceStreams;
+}
+
+/** 
+ * Get inner loop unroll
+ * 
+ * @return   inner loop unroll (fixed value)
+ *
+ */
+
+int BrookNonBonded::getInnerLoopUnroll( void ) const {
+   return _innerUnroll;
+}
+
+/** 
+ * Get outer loop unroll
+ * 
+ * @return   outer loop unroll (fixed value)
+ *
+ */
+
+int BrookNonBonded::getOuterLoopUnroll( void ) const {
+   return _outerUnroll;
+}
+
+/** 
+ * Set outer loop unroll
+ * 
+ * @param  outer loop unroll (fixed value)
+ *
+ * @return updated outer loop unroll (fixed value)
+ *
+ */
+
+int BrookNonBonded::setOuterLoopUnroll( int outerUnroll ){
+   if( outerUnroll != _outerUnroll ){
+      _atomSizeCeiling = -1;
+   }
+   _outerUnroll = _outerUnroll;
+   return _outerUnroll;
+}
+
+/** 
+ * Get atom ceiling parameter
+ * 
+ * @return atom ceiling parameter
+ *
+ */
+
+int BrookNonBonded::getAtomSizeCeiling( void ) const {
+
+   if( _atomSizeCeiling < 0 ){
+      BrookNonBonded* localThis = const_cast<BrookNonBonded* const>(this);
+      localThis->_atomSizeCeiling = localThis->getNumberOfAtoms() % localThis->getOuterLoopUnroll();
+      if( localThis->_atomSizeCeiling ){
+         localThis->_atomSizeCeiling = localThis->getOuterLoopUnroll() - localThis->_atomSizeCeiling;
+      }   
+   }
+
+   return _atomSizeCeiling;
+}
+
+/** 
+ * Get duplication factor
+ * 
+ * @return   duplication factor
+ *
+ */
+
+int BrookNonBonded::getDuplicationFactor( void ) const {
+   return _duplicationFactor;
+}
+
+/** 
+ * Get j-stream width
+ * 
+ * @param platform  platform
+ *
+ * @return  j-stream width
+ *
+ */
+
+int BrookNonBonded::getJStreamWidth( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getJStreamWidth";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   // get j-stream width
+
+   if( _jStreamWidth < 0 ){
+      //_getJStreamDimensions( platform );
+   }
+   return _jStreamWidth;
+}
+
+/** 
+ * Get j-stream width
+ * 
+ * @return  j-stream width
+ *
+ */
+
+int BrookNonBonded::getJStreamWidth( void ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getJStreamWidth";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   return _jStreamWidth;
+}
+
+/** 
+ * Get j-stream height
+ * 
+ * @param platform platform
+ *
+ * @return  j-stream height 
+ *
+ */
+
+int BrookNonBonded::getJStreamHeight( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getJStreamHeight";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   // get j-stream height
+
+   if( _jStreamHeight < 0 ){
+      //_getJStreamDimensions( platform );
+   }
+   return _jStreamHeight;
+}
+
+/** 
+ * Get j-stream height
+ * 
+ * @return  j-stream height 
+ *
+ */
+
+int BrookNonBonded::getJStreamHeight( void ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getJStreamHeight";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   return _jStreamHeight;
+}
+
+/** 
+ * Get j-stream size
+ * 
+ * @param platform  platform
+ *
+ * @return  j-stream size
+ *
+ */
+
+int BrookNonBonded::getJStreamSize( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getJStreamSize";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   // get j-stream size
+
+   if( _jStreamSize < 0 ){
+      //_getJStreamDimensions( platform );
+   }
+   return _jStreamSize;
+}
+
+/** 
+ * Get j-stream size
+ * 
+ * @return  j-stream size
+ *
+ */
+
+int BrookNonBonded::getJStreamSize( void ) const {
+   return _jStreamSize;
+}
+
+/** 
+ * Get partial force stream width
+ * 
+ * @param platform  platform
+ *
+ * @return  partial force stream width
+ *
+ */
+
+int BrookNonBonded::getPartialForceStreamWidth( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getPartialForceStreamWidth";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   // get partial force stream width
+
+   if( _partialForceStreamWidth < 0 ){
+      //_getPartialForceStreamDimensions( platform );
+   }
+   return _partialForceStreamWidth;
+}
+
+/** 
+ * Get partial force stream width
+ * 
+ * @return  partial force stream width
+ *
+ */
+
+int BrookNonBonded::getPartialForceStreamWidth( void ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getPartialForceStreamWidth";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   return _partialForceStreamWidth;
+}
+
+/** 
+ * Get partial force stream height
+ * 
+ * @param platform platform
+ *
+ * @return  partial force stream height 
+ *
+ */
+
+int BrookNonBonded::getPartialForceStreamHeight( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getPartialForceStreamHeight";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   // get partial force stream height
+
+   if( _partialForceStreamHeight < 0 ){
+      //_getPartialForceStreamDimensions( platform );
+   }
+   return _partialForceStreamHeight;
+}
+
+/** 
+ * Get partial force stream height
+ * 
+ * @return  partial force stream height 
+ *
+ */
+
+int BrookNonBonded::getPartialForceStreamHeight( void ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getPartialForceStreamHeight";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   return _partialForceStreamHeight;
+}
+
+/** 
+ * Get partial force stream size
+ * 
+ * @param platform  platform
+ *
+ * @return  partial force stream size
+ *
+ */
+
+int BrookNonBonded::getPartialForceStreamSize( const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::getPartialForceStreamSize";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   // get partial force stream size
+
+   if( _partialForceStreamSize < 0 ){
+      //_getPartialForceStreamDimensions( platform );
+   }
+   return _partialForceStreamSize;
+}
+
+/** 
+ * Get partial force stream size
+ * 
+ * @return  partial force stream size
+ *
+ */
+
+int BrookNonBonded::getPartialForceStreamSize( void ) const {
+   return _partialForceStreamSize;
+}
+
+/** 
+ * Get exclusion stream size
+ *
+ * @return  exclusion stream size
+ *
+ */
+
+int BrookNonBonded::getExclusionStreamSize( void ) const {
+   return _exclusionStreamSize;
+}
+
+/** 
+ * Get exclusion stream width
+ *
+ * @return  exclusion stream width
+ *
+ */
+
+int BrookNonBonded::getExclusionStreamWidth( void ) const {
+   return _exclusionStreamWidth;
+}
+
+/** 
+ * Get exclusion stream 
+ *
+ * @return  exclusion stream
+ *
+ */
+
+BrookFloatStreamImpl* BrookNonBonded::getExclusionStream( void ) const {
+   return _exclusionStream;
+}
+
+/** 
+ * Get vdw stream 
+ *
+ * @return  vdw stream
+ *
+ */
+
+BrookFloatStreamImpl* BrookNonBonded::getVdwStream( void ) const {
+   return _vdwStream;
+}
+
+/** 
+ * Get charge stream 
+ *
+ * @return  charge stream
+ *
+ */
+
+BrookFloatStreamImpl* BrookNonBonded::getChargeStream( void ) const {
+   return _chargeStream;
+}
+
+/** 
+ * Get sigma stream 
+ *
+ * @return  sigma stream
+ *
+ */
+
+BrookFloatStreamImpl* BrookNonBonded::getSigmaStream( void ) const {
+   return _sigmaStream;
+}
+
+/** 
+ * Get epsilon stream 
+ *
+ * @return  epsilon stream
+ *
+ */
+
+BrookFloatStreamImpl* BrookNonBonded::getEpsilonStream( void ) const {
+   return _epsilonStream;
+}
+
+/** 
+ * Get force streams 
+ *
+ * @return  force streams
+ *
+ */
+
+BrookFloatStreamImpl** BrookNonBonded::getForceStreams( void ){
+   return _nonbondedForceStreams;
+}
+
+/** 
+ * Return true if force[index] stream is set
+ *
+ * @param    index into force stream
+ * @return   true if index is valid && force[index] stream is set; else false
+ *
+ */
+
+int BrookNonBonded::isForceStreamSet( int index ) const {
+   return (index >= 0 && index < getNumberOfForceStreams() && _nonbondedForceStreams[index]) ? 1 : 0;
+}
+
+/** 
+ * Initialize stream dimensions and streams
+ * 
+ * @param numberOfAtoms             number of atoms
+ * @param platform                  platform
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookNonBonded::initializeStreams( int numberOfAtoms, const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookNonBonded::initializeStreams";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   int atomStreamSize        = getAtomStreamSize( platform );
+   _exclusionStreamWidth     = atomStreamSize;
+
+   int innerUnroll           = getInnerLoopUnroll();
+   if( innerUnroll < 1 ){
+      std::stringstream message;
+      message << methodName << " innerUnrolls=" << innerUnroll << " is less than 1.";
+      throw OpenMMException( message.str() );
+      return ErrorReturnValue;
+   }
+
+   const BrookPlatform brookPlatform            = dynamic_cast<const BrookPlatform&> (platform);
+   const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (platform.getDefaultStreamFactory() );
+
+   _exclusionStreamHeight    = _exclusionStreamWidth/innerUnroll; 
+   _exclusionStreamSize      = _exclusionStreamWidth*_exclusionStreamHeight;
+   StreamImpl* exStream      = brookStreamFactory.createStreamImpl( BrookStreamFactory::NonBondedExclusionStream, _exclusionStreamSize, Stream::Float, platform );
+   _exclusionStream          = dynamic_cast<BrookFloatStreamImpl*> ( exStream );
+
+
+   StreamImpl* vStream       = brookStreamFactory.createStreamImpl( BrookStreamFactory::NonBondedVdwStream, atomStreamSize, Stream::Float3, platform );
+   _vdwStream                = dynamic_cast<BrookFloatStreamImpl*> ( vStream );
+
+   StreamImpl* cStream       = brookStreamFactory.createStreamImpl( BrookStreamFactory::NonBondedVdwStream, atomStreamSize, Stream::Float3, platform );
+   _chargeStream             = dynamic_cast<BrookFloatStreamImpl*> ( cStream );
+
+   if( getLog() ){
+      int atomStreamWidth       = getAtomStreamWidth( platform );
+      int atomStreamHeight      = getAtomStreamHeight( platform );
+      (void) fprintf( getLog(), "%s initializeVdwStream: stream dimensions: [%d %d] size=%d\n", methodName.c_str(), atomStreamWidth, atomStreamHeight, atomStreamSize );
+   }
+
+   if( getLog() ){
+      (void) fprintf( getLog(), "%s initializeExclusionStream: stream dimensions: [%d %d] size=%d\n", methodName.c_str(), _exclusionStreamWidth, _exclusionStreamHeight, _exclusionStreamSize );
+   }
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Set exclusion (4x4)
+ * 
+ * @param i                         atom i index
+ * @param j                         atom j index
+ * @param exclusionStreamWidth      exclusion stream width
+ * @param exclusion                 array of packed exclusions
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookNonBonded::setExclusion( int i, int j, int exclusionStreamWidth, BrookOpenMMFloat* exclusion ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookNonBonded::setExclusion";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   int index    = ( (j/4)* exclusionStreamWidth + i ); //iindex is the x coordinate!
+   int joffset  = j % 4; 
+
+   BrookOpenMMFloat flag;
+   switch( joffset ) {
+
+      case 0: 
+         flag = (BrookOpenMMFloat) 2.0; 
+         break;
+
+      case 1:
+         flag = (BrookOpenMMFloat) 3.0; 
+         break;
+
+      case 2: 
+         flag = (BrookOpenMMFloat) 5.0; 
+         break;
+
+      case 3 :  
+         flag = (BrookOpenMMFloat) 7.0; 
+         break;
+
+   }
+
+   exclusion[index] /= flag;
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Initialize exclusions
+ * 
+ * @param exclusions            vector of sets containing exclusions (1 set entry for every atom)
+ * @param platform              Brook platform
+ * @param log                   optional Log file reference
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookNonBonded::initializeExclusions( const std::vector<std::set<int> >& exclusionsVector, const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookNonBonded::initializeExclusions";
+   static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   if( debug && getLog() ){
+      (void) fprintf( getLog(), "In initializeExclusions exclusions vector size=%d\n", exclusionsVector.size() );
+   }
+
+   int exclusionStreamSize = getExclusionStreamSize();
+   if( exclusionStreamSize < 1 ){
+      std::stringstream message;
+      message << "BrookNonBonded::initializeExclusions exclusionStreamSize=" << exclusionStreamSize << " is not set.";
+      throw OpenMMException( message.str() );
+      return ErrorReturnValue;
+   }
+
+   // allocate and initialize exclusions array
+
+   BrookOpenMMFloat* exclusions = new BrookOpenMMFloat[exclusionStreamSize];
+   for( unsigned int ii = 0; ii < (unsigned int) exclusionStreamSize; ii++ ){
+      exclusions[ii] = (BrookOpenMMFloat) 210.0;
+   }
+
+   // pack in values
+
+   int  exclusionStreamWidth = getExclusionStreamWidth();
+   int  numberOfAtoms        = getNumberOfAtoms();
+   int  atomStreamSize       = getAtomStreamSize();
+   for( unsigned int ii = 0; ii < exclusionsVector.size(); ii++ ){
+
+      set<int> exclusionIndices  = exclusionsVector[ii];
+      for( set<int>::const_iterator jj = exclusionIndices.begin(); jj != exclusionIndices.end(); jj++ ){
+         setExclusion( ii, *jj, exclusionStreamWidth, exclusions );
+      }
+
+      // explicitly exclude junk atoms from interacting w/ atom ii
+
+      for( int jj = numberOfAtoms; jj < atomStreamSize; jj++ ){
+         setExclusion( ii, jj, exclusionStreamWidth, exclusions );
+      }
+   }
+
+   // explicitly exclude junk atoms from interacting w/ all atoms
+
+   for( int ii = numberOfAtoms; ii < atomStreamSize; ii++ ){
+      for( int jj = 0; jj < atomStreamSize; jj++ ){
+         setExclusion( ii, jj, exclusionStreamWidth, exclusions );
+      }
+   }
+
+   // load stream
+
+   _exclusionStream->loadFromArray( exclusions ); 
+   delete[] exclusions;
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Set sigma & epsilon given c6 & c12 (geometric rule)
+ * 
+ * @param c6                        vdw c6
+ * @param c12                       vdw c12
+ * @param sigma                     massaged sigma
+ * @param epsilon                   massaged epsilon
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookNonBonded::setSigmaEpsilon( double c6, double c12, double* sigma , double* epsilon ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookNonBonded::setSigmaEpsilon";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   if( fabs( c12 ) < 1.0e-04 ){
+      *sigma    = 1.0;
+      *epsilon  = 0.0;
+   } else {
+      *sigma    = 0.5*pow( c12/c6, 1.0/6.0 );
+      *epsilon  = sqrt( c6*c6/c12 );
+   }
+
+   return DefaultReturnValue;
+}
+
+/** 
+ * Initialize vdw & charge
+ * 
+ * @param exclusions                vector of sets containing exclusions (1 set entry for every atom)
+ * @param platform                  platform
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ */
+
+int BrookNonBonded::initializeVdwAndCharge( const vector<vector<double> >& nonbondedParameters, const Platform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookNonBonded::initializeVdwAndCharge";
+   static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   if( debug && getLog() ){
+      (void) fprintf( getLog(), "In initializeVdwAndCharge nonbonded vector size=%d\n", nonbondedParameters.size() );
+   }
+
+   int atomStreamSize        = getAtomStreamSize();
+
+   // allocate and initialize vdw & charge array
+
+   unsigned int vdwParametersSize  = 2*atomStreamSize;
+   BrookOpenMMFloat* vdwParameters = new BrookOpenMMFloat[vdwParametersSize];
+   memset( vdwParameters, 0, sizeof( BrookOpenMMFloat )*vdwParametersSize );
+   
+   BrookOpenMMFloat* charges       = new BrookOpenMMFloat[atomStreamSize];
+   memset( charges, 0, sizeof( BrookOpenMMFloat )*atomStreamSize );
+   
+   // pack in values
+
+   int  numberOfAtoms        = getNumberOfAtoms();
+   double sigma, epsilon;
+   int trackingIndex         = 0;
+   for( unsigned int ii = 0; ii < nonbondedParameters.size(); ii++ ){
+
+      vector<double> nonbondedParameterVector = nonbondedParameters[ii];
+
+      double c6                      = nonbondedParameterVector[0];
+      double c12                     = nonbondedParameterVector[1];
+      double charge                  = nonbondedParameterVector[2];
+
+      // geometric combination rules
+
+      setSigmaEpsilon( c6, c12, &sigma, &epsilon );
+
+      vdwParameters[trackingIndex++] = (BrookOpenMMFloat) sigma;
+      vdwParameters[trackingIndex++] = (BrookOpenMMFloat) epsilon;
+      charges[ii]                    = (BrookOpenMMFloat) charge;
+         
+   }
+
+   // set outlier atoms
+
+   for( unsigned int ii = nonbondedParameters.size(); ii < (unsigned int) atomStreamSize; ii++ ){
+      vdwParameters[trackingIndex++] = (BrookOpenMMFloat) 1.0;
+      vdwParameters[trackingIndex++] = (BrookOpenMMFloat) 0.0;
+   }
+
+   // load stream
+
+   _vdwStream->loadFromArray( vdwParameters ); 
+   _chargeStream->loadFromArray( charges ); 
+
+   delete[] vdwParameters;
+   delete[] charges;
+
+   return DefaultReturnValue;
+}
+
+/* 
+ * Setup of nonbonded ixns
+ *
+ * @param numberOfAtoms         number of atoms
+ * @param nonbondedParameters   vector of nonbonded parameters [atomI][0=c6]
+ *                                                             [atomI][1=c12]
+ *                                                             [atomI][2=charge]
+ * @param platform              Brook platform
+ *
+ * @return ErrorReturnValue if error, else DefaultReturnValue
+ *
+ * */
+
+int BrookNonBonded::setup( int numberOfAtoms, const std::vector<std::vector<double> >& nonbondedParameters,
+                           const std::vector<std::set<int> >& exclusions, const BrookPlatform& platform ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookNonBonded::setup";
+   static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   setNumberOfAtoms( numberOfAtoms );
+
+   initializeStreams( numberOfAtoms, platform );
+
+   initializeExclusions( exclusions, platform );
+   initializeVdwAndCharge( nonbondedParameters, platform );
+   //initializeJStreamVdw( nonbondedParameters, platform );
+
+   return DefaultReturnValue;
+}
+
+/* 
+ * Get contents of object
+ *
+ *
+ * @param level   level of dump
+ *
+ * @return string containing contents
+ *
+ * */
+
+std::string BrookNonBonded::getContents( int level ) const {
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookNonBonded::getContents";
+
+   static const unsigned int MAX_LINE_CHARS = 256;
+   char value[MAX_LINE_CHARS];
+   static const char* Set                   = "Set";
+   static const char* NotSet                = "Not set";
+
+// ---------------------------------------------------------------------------------------
+
+   std::stringstream message;
+   std::string tab   = "   ";
+
+#ifdef WIN32
+#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
+#else
+#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
+#endif
+
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfAtoms() );
+   message << _getLine( tab, "Number of atoms:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getNumberOfForceStreams() );
+   message << _getLine( tab, "Number of force streams:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getDuplicationFactor() );
+   message << _getLine( tab, "Duplication factor:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getInnerLoopUnroll () )
+   message << _getLine( tab, "Inner loop unroll:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getOuterLoopUnroll() )
+   message << _getLine( tab, "Outer loop unroll:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getAtomSizeCeiling() );
+   message << _getLine( tab, "Atom ceiling:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamWidth() );
+   message << _getLine( tab, "Atom stream width:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamHeight() );
+   message << _getLine( tab, "Atom stream height:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getAtomStreamSize() );
+   message << _getLine( tab, "Atom stream size:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getJStreamWidth() );
+   message << _getLine( tab, "J-stream width:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getJStreamHeight() );
+   message << _getLine( tab, "J-stream height:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getJStreamSize() );
+   message << _getLine( tab, "J-stream size:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getPartialForceStreamWidth() );
+   message << _getLine( tab, "Partial force stream width:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getPartialForceStreamHeight() );
+   message << _getLine( tab, "Partial force stream height:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getPartialForceStreamSize() );
+   message << _getLine( tab, "Partial force stream size:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getExclusionStreamWidth() );
+   message << _getLine( tab, "Exclusion stream width:", value ); 
+
+   //(void) LOCAL_SPRINTF( value, "%d", getExclusionStreamHeight() );
+   //message << _getLine( tab, "Exclusion stream height:", value ); 
+
+   (void) LOCAL_SPRINTF( value, "%d", getExclusionStreamSize() );
+   message << _getLine( tab, "Exclusion stream size:", value ); 
+
+   message << _getLine( tab, "Log:",                  (getLog()             ? Set : NotSet) ); 
+   message << _getLine( tab, "ExclusionStream:",      (getExclusionStream() ? Set : NotSet) ); 
+   message << _getLine( tab, "VdwStream:",            (getVdwStream()       ? Set : NotSet) ); 
+   message << _getLine( tab, "ChargeStream:",         (getChargeStream()    ? Set : NotSet) ); 
+   message << _getLine( tab, "SigmaStream:",          (getSigmaStream()     ? Set : NotSet) ); 
+   message << _getLine( tab, "EpsilonStream:",        (getEpsilonStream()   ? Set : NotSet) ); 
+ 
+   for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
+      char description[256];
+      (void) LOCAL_SPRINTF( description, "PartialForceStream %d", ii );
+      message << _getLine( tab, description,  (isForceStreamSet(ii) ? Set : NotSet) ); 
+   }
+ 
+#undef LOCAL_SPRINTF
+
+   return message.str();
+}

platforms/brook/src/BrookNonBonded.h

+#ifndef BrookNonBonded_H_
+#define BrookNonBonded_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include <vector>
+#include <set>
+
+#include "BrookFloatStreamImpl.h"
+#include "BrookIntStreamImpl.h"
+#include "BrookPlatform.h"
+#include "BrookCommon.h"
+
+namespace OpenMM {
+
+/**
+ * This kernel is invoked by StandardMMForceField to calculate the forces acting on the system.
+ */
+class BrookNonBonded : public BrookCommon {
+
+   public:
+  
+      // return values
+
+      static const int DefaultReturnValue = 0;
+      static const int ErrorReturnValue   = -1;
+
+      BrookNonBonded( );
+  
+      ~BrookNonBonded();
+  
+      /**
+       * Return number of force streams
+       * 
+       * @return number of force streams
+       *
+       */
+
+      int getNumberOfForceStreams( void ) const; 
+
+      /** 
+       * Get duplication factor
+       * 
+       * @return   duplication factor
+       *
+       */
+      
+      int getDuplicationFactor( void ) const;
+      
+      /** 
+       * Get atom ceiling parameter
+       * 
+       * @return atom ceiling parameter
+       *
+       */
+         
+      int getAtomSizeCeiling( void ) const;
+ 
+      /** 
+       * Get outer loop unroll
+       * 
+       * @return   outer loop unroll (fixed value)
+       *
+       */
+      
+      int getOuterLoopUnroll( void ) const;
+      
+      /** 
+       * Set outer loop unroll
+       * 
+       * @param  outer loop unroll (fixed value)
+       *
+       * @return updated outer loop unroll (fixed value)
+       *
+       */
+      
+      int setOuterLoopUnroll( int outerUnroll );
+      
+      /**
+       * Return unrolling for inner loops
+       * 
+       * @return outer loop unrolling
+       */
+
+      int getInnerLoopUnroll( void ) const; 
+
+      /**
+       * Get j-stream width
+       *
+       * @param platform platform reference
+       *
+       * @return j-stream width
+       */
+
+      int getJStreamWidth( const Platform& platform ); 
+
+      /**
+       * Get j-stream width
+       *
+       * @return j-stream width
+       */
+
+      int getJStreamWidth( void ) const; 
+
+      /**
+       * Get j-stream height
+       *
+       * @param platform platform reference
+       *
+       * @return j-stream height
+       */
+
+      int getJStreamHeight( const Platform& platform ); 
+
+      /**
+       * Get j-stream height
+       *
+       * @return j-stream height
+       */
+
+      int getJStreamHeight( void ) const;
+
+      /**
+       * Get j-stream size
+       * 
+       * @param platform platform reference
+       *
+       * @return j-stream size
+       */
+
+      int getJStreamSize( const Platform& platform ); 
+
+      /**
+       * Get j-stream size
+       * 
+       * @return j-stream size
+       */
+
+      int getJStreamSize( void ) const; 
+
+      /**
+       * Get partial force stream width
+       *
+       * @param platform platform reference
+       *
+       * @return partial force stream width
+       */
+
+      int getPartialForceStreamWidth( const Platform& platform ); 
+
+      /**
+       * Get partial force stream width
+       *
+       * @return partial force stream width
+       */
+
+      int getPartialForceStreamWidth( void ) const; 
+
+      /**
+       * Get partial force stream height
+       *
+       * @param platform platform reference
+       *
+       * @return partial force stream height
+       */
+
+      int getPartialForceStreamHeight( const Platform& platform ); 
+
+      /**
+       * Get partial force stream height
+       *
+       * @return partial force stream height
+       */
+
+      int getPartialForceStreamHeight( void ) const;
+
+      /**
+       * Get partial force stream size
+       * 
+       * @param platform platform reference
+       *
+       * @return partial force stream size
+       */
+
+      int getPartialForceStreamSize( const Platform& platform ); 
+
+      /**
+       * Get partial force stream size
+       * 
+       * @return partial force stream size
+       */
+
+      int getPartialForceStreamSize( void ) const; 
+
+      /**
+       * Get partial force stream size
+       * 
+       * @return partial force stream size
+       */
+
+      /**
+       * Get exclusion stream width
+       * 
+       * @return exclusion stream width
+       */
+
+      int getExclusionStreamWidth( void ) const; 
+
+      /**
+       * Get exclusion stream size
+       * 
+       * @return exclusion stream size
+       */
+
+      int getExclusionStreamSize( void ) const; 
+
+      /** 
+       * Get exclusion stream 
+       *
+       * @return  exclusion stream
+       *
+       */
+      
+      BrookFloatStreamImpl* getExclusionStream( void ) const;
+      
+      /** 
+       * Get vdw stream 
+       *
+       * @return  vdw stream
+       *
+       */
+      
+      BrookFloatStreamImpl* getVdwStream( void ) const;
+      
+      /** 
+       * Get charge stream 
+       *
+       * @return  charge stream
+       *
+       */
+      
+      BrookFloatStreamImpl* getChargeStream( void ) const;
+      
+      /** 
+       * Get sigma-eps stream 
+       *
+       * @return  sigma-eps stream
+       *
+       */
+      
+      BrookFloatStreamImpl* getSigmaStream( void ) const;
+      
+      /** 
+       * Get epsilon stream 
+       *
+       * @return  epsilon stream
+       *
+       */
+      
+      BrookFloatStreamImpl* getEpsilonStream( void ) const;
+      
+      /** 
+       * Get force streams 
+       *
+       * @return  force streams
+       *
+       */
+      
+      BrookFloatStreamImpl** getForceStreams( void );
+      
+      /** 
+       * Return true if force[index] stream is set 
+       *
+       * @return  true  if index is valid && force[index] stream is set; else false
+       *
+       */
+      
+      int isForceStreamSet( int index ) const;
+      
+      /* 
+       * Setup of nonbonded ixns
+       *
+       * @param numberOfAtoms         number of atoms
+       * @param nonbondedParameters   vector of nonbonded parameters [atomI][0=c6]
+       *                                                             [atomI][1=c12]
+       *                                                             [atomI][2=charge]
+       * @param platform              Brook platform
+       * @param log                   optional Log file reference
+       *
+       * @return nonzero value if error
+       * */
+      
+      int setup( int numberOfAtoms, const std::vector<std::vector<double> >& nonbondedParameters,
+                 const std::vector<std::set<int> >& exclusions,  const BrookPlatform& platform );
+      
+      /* 
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string getContents( int level ) const;
+
+   private:
+   
+      // fixed number of force streams
+
+      static const int NumberOfForceStreams     = 4;
+
+      // atom ceiling
+
+      int _atomSizeCeiling;
+
+      // unroll in i/j dimensions
+
+      int _outerUnroll;
+      int _innerUnroll;
+
+      // duplication factor
+
+      int _duplicationFactor;
+
+      // force stream width
+
+      int _partialForceStreamWidth;
+      int _partialForceStreamHeight;
+      int _partialForceStreamSize;
+
+      // exclusions stream dimensions
+
+      int _exclusionStreamWidth;
+      int _exclusionStreamHeight;
+      int _exclusionStreamSize;
+
+      // j-stream dimensions
+
+      int _jStreamWidth;
+      int _jStreamHeight;
+      int _jStreamSize;
+
+      // streams
+
+      BrookFloatStreamImpl* _exclusionStream;
+      BrookFloatStreamImpl* _vdwStream;
+      BrookFloatStreamImpl* _chargeStream;
+      BrookFloatStreamImpl* _sigmaStream;
+      BrookFloatStreamImpl* _epsilonStream;
+      BrookFloatStreamImpl* _nonbondedForceStreams[NumberOfForceStreams];
+
+      /** 
+       * Initialize exclusion stream dimensions and stream
+       * 
+       * @param platform                  platform
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int initializeExclusionStream( const Platform& platform );
+            
+      /** 
+       * Set exclusion (4x4)
+       * 
+       * @param i                         atom i index
+       * @param j                         atom j index
+       * @param exclusionStreamWidth      exclusion stream width
+       * @param exclusion                 array of packed exclusions
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int setExclusion( int i, int j, int exclusionStreamWidth, BrookOpenMMFloat* exclusion );
+      
+      /** 
+       * Initialize exclusions
+       * 
+       * @param exclusions                vector of sets containing exclusions (1 set entry for every atom)
+       * @param platform                  platform
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int initializeExclusions( const std::vector<std::set<int> >& exclusionsVector, const Platform& platform );
+
+      /** 
+       * Initialize stream dimensions and streams
+       * 
+       * @param platform                  platform
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int initializeStreams( int numberOfAtoms, const Platform& platform );
+      
+      /** 
+       * Set sigma & epsilon given c6 & c12 (geometric rule)
+       * 
+       * @param c6                        vdw c6
+       * @param c12                       vdw c12
+       * @param sigma                     massaged sigma
+       * @param epsilon                   massaged epsilon
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int setSigmaEpsilon( double c6, double c12, double* sigma , double* epsilon );
+
+      /** 
+       * Initialize vdw & charge
+       * 
+       * @param exclusions                vector of sets containing exclusions (1 set entry for every atom)
+       * @param platform                  platform
+       *
+       * @return nonzero value if error
+       *
+       */
+      
+      int initializeVdwAndCharge( const std::vector<std::vector<double> >& nonbondedParameters, const Platform& platform );
+      
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_BROOKKERNELS_H_*/