Mods for OBC Gbsa

platforms/brook/src/BrookBonded.cpp

@@ -1054,7 +1054,7 @@ int BrookBonded::loadInvMaps( int nbondeds, int natoms, int *atoms, const BrookP
 
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
       for( int jj = 0; jj < getMaxInverseMapStreamCount(ii); jj++ ){
-         _inverseStreamMaps[ii][jj] = new BrookFloatStreamInternal( BrookStreamFactory::BondedInverseMapStreams, atomStreamSize,
+         _inverseStreamMaps[ii][jj] = new BrookFloatStreamInternal( BrookCommon::BondedInverseMapStreams, atomStreamSize,
                                                                     atomStreamWidth, BrookStreamInternal::Float4, dangleValue );
       }
    }
@@ -1239,16 +1239,16 @@ int BrookBonded::setup( int numberOfAtoms,
 
    // build streams
 
-   _atomIndicesStream    = new BrookFloatStreamInternal( BrookStreamFactory::BondedAtomIndicesStream, nbondeds, atomStreamWidth,
+   _atomIndicesStream    = new BrookFloatStreamInternal( BrookCommon::BondedAtomIndicesStream, nbondeds, atomStreamWidth,
                                                          BrookStreamInternal::Float4, dangleValue );
    _atomIndicesStream->loadFromArray( atoms, BrookStreamInternal::Integer ); 
 
-   _chargeStream         = new BrookFloatStreamInternal( BrookStreamFactory::BondedChargeStream, numberOfAtoms, atomStreamWidth,
+   _chargeStream         = new BrookFloatStreamInternal( BrookCommon::BondedChargeStream, numberOfAtoms, atomStreamWidth,
                                                          BrookStreamInternal::Float, dangleValue );
    _chargeStream->loadFromArray( charges ); 
 
    for( int ii = 0; ii < getNumberOfParameterStreams(); ii++ ){
-      _bondedParameters[ii]  = new BrookFloatStreamInternal( BrookStreamFactory::BondedParametersStream, nbondeds, atomStreamWidth,
+      _bondedParameters[ii]  = new BrookFloatStreamInternal( BrookCommon::BondedParametersStream, nbondeds, atomStreamWidth,
                                                              BrookStreamInternal::Float4, dangleValue );
       _bondedParameters[ii]->loadFromArray( params[ii] );
    }
@@ -1316,7 +1316,7 @@ int BrookBonded::setup( int numberOfAtoms,
    // initialize output streams
 
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
-      _bondedForceStreams[ii] = new BrookFloatStreamInternal( BrookStreamFactory::UnrolledForceStream, nbondeds, atomStreamWidth, 
+      _bondedForceStreams[ii] = new BrookFloatStreamInternal( BrookCommon::UnrolledForceStream, nbondeds, atomStreamWidth, 
                                                               BrookStreamInternal::Float3, dangleValue );
    }
 

platforms/brook/src/BrookCalcGBSAOBCForceFieldKernel.cpp

@@ -182,6 +182,12 @@ void BrookCalcGBSAOBCForceFieldKernel::executeForces( const Stream& positions, S
    float includeAce                                    = (float) (_brookGbsa->includeAce());
    BrookFloatStreamInternal**  gbsaForceStreams        = _brookGbsa->getForceStreams();
 
+   // calculate Born radii first time thru and initialize on board
+
+   if( _brookGbsa->haveBornRadiiBeenInitialized() ){
+      _brookGbsa->calculateBornRadii( positions );  
+   }
+
    // first major loop
 
    kObcLoop1( (float) _brookGbsa->getNumberOfAtoms(),

platforms/brook/src/BrookCommon.cpp

@@ -39,6 +39,33 @@
 using namespace OpenMM;
 using namespace std;
 
+// bonded streams
+    
+const std::string BrookCommon::BondedAtomIndicesStream                            = "BondedAtomIndicesStream";
+const std::string BrookCommon::BondedParametersStream                             = "BondedParametersStream";
+const std::string BrookCommon::UnrolledForceStream                                = "UnrolledForceStream";
+const std::string BrookCommon::BondedChargeStream                                 = "BondedChargeStream";
+const std::string BrookCommon::BondedInverseMapStreams                            = "BondedInverseMapStreams";
+
+// non-bonded streams
+
+const std::string BrookCommon::NonBondedExclusionStream                           = "NonBondedExclusionStream";
+const std::string BrookCommon::OuterVdwStream                                     = "OuterVdwStream";
+const std::string BrookCommon::InnerSigmaStream                                   = "InnerSigmaStream";
+const std::string BrookCommon::InnerEpsilonStream                                 = "InnerEpsilonStream";
+const std::string BrookCommon::NonBondedChargeStream                              = "NonBondedChargeStream";
+const std::string BrookCommon::PartialForceStream                                 = "PartialForceStream";
+
+// OBC Gbsa streams
+
+const std::string BrookCommon::ObcAtomicRadiiStream                               = "ObcAtomicRadiiStream";
+const std::string BrookCommon::ObcScaledAtomicRadiiStream                         = "ObcScaledAtomicRadiiStream";
+const std::string BrookCommon::ObcAtomicRadiiWithDielectricOffsetStream           = "ObcAtomicRadiiWithDielectricOffsetStream";
+const std::string BrookCommon::ObcBornRadiiStream                                 = "ObcBornRadiiStream";
+const std::string BrookCommon::ObcBornRadii2Stream                                = "ObcBornRadii2Stream";
+const std::string BrookCommon::ObcIntermediateForceStream                         = "ObcIntermediateForceStream";
+const std::string BrookCommon::ObcChainStream                                     = "ObcChainStream";
+
 /** 
  * Constructor
  * 

platforms/brook/src/BrookCommon.h

@@ -53,8 +53,51 @@ class BrookCommon {
       static const int DefaultReturnValue = 0;
       static const int ErrorReturnValue   = -1;
 
+      // ---------------------------------------------------------------------------------------
+
+      // Stream names
+     
+      // bonded stream names
+
+      static const std::string BondedAtomIndicesStream;
+      static const std::string BondedParametersStream;
+      static const std::string UnrolledForceStream;
+      static const std::string BondedChargeStream;
+      static const std::string BondedInverseMapStreams;
+
+      // nonbonded stream names
+
+      static const std::string NonBondedExclusionStream;
+      static const std::string OuterVdwStream;
+      static const std::string InnerSigmaStream;
+      static const std::string InnerEpsilonStream;
+      static const std::string NonBondedChargeStream;
+      static const std::string PartialForceStream;
+
+      // OBC Gbsa streams
+
+      static const std::string ObcAtomicRadiiStream;
+      static const std::string ObcScaledAtomicRadiiStream;
+      static const std::string ObcAtomicRadiiWithDielectricOffsetStream;
+      static const std::string ObcBornRadiiStream;
+      static const std::string ObcBornRadii2Stream;
+      static const std::string ObcIntermediateForceStream;
+      static const std::string ObcChainStream;
+
+      // ---------------------------------------------------------------------------------------
+
+     /** 
+      * Constructor
+      * 
+      */
+     
       BrookCommon( );
   
+     /** 
+      * Destructor
+      * 
+      */
+     
       ~BrookCommon();
   
       /**

platforms/brook/src/BrookFloatStreamInternal.cpp

@@ -357,34 +357,29 @@ void BrookFloatStreamInternal::fillWithValue( void* value ){
 
 }
 
-const RealOpenMM* const * BrookFloatStreamInternal::getData() const {
-   return NULL;
-}
-
-// problem w/ const here -- _data is modified (cast away?)
+/** 
+ * Get data
+ * 
+ * @param readFromBoard if set, read values on board 
+ *
+ * @return data array
+ *
+ */
 
-/*
-const RealOpenMM* const * BrookFloatStreamInternal::getData() const {
+void* BrookFloatStreamInternal::getData( int readFromBoard ){
 
-   // retrieve _data from GPU
+// ---------------------------------------------------------------------------------------
 
-   _aStream.write( _data );
+   // static const std::string methodName      = "BrookFloatStreamInternal::getData";
 
-   // check if RealOpenMM is float; if not, then 
-   // copy into realOpenMMData[][] array
+// ---------------------------------------------------------------------------------------
 
-   if( realOpenMMData ){
-      for( int i = 0; i < getSize(); i++ ){
-         for( int j = 0; j < getWidth(); j++ ){
-            realOpenMMData[i][j] = (RealOpenMM) _data[i][j];
-         }
-      }
-      return realOpenMMData;
-   } else {
-      return _data;
+   if( readFromBoard ){
+      _aStream.write( _data );
    }
+
+   return (void*) _data;
 }
-*/
 
 /** 
  * Get data
@@ -393,28 +388,15 @@ const RealOpenMM* const * BrookFloatStreamInternal::getData() const {
  *
  */
 
-RealOpenMM** BrookFloatStreamInternal::getData( void ){
+void* BrookFloatStreamInternal::getData( void ){
 
 // ---------------------------------------------------------------------------------------
 
    // static const std::string methodName      = "BrookFloatStreamInternal::getData";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
-/*
-   _aStream.write( _data );
-   if( _realOpenMMData ){
-      int totalSize                             = getSize()*getWidth();
-      for( int ii = 0; ii < totalSize; ii++ ){
-         _realOpenMMData[ii] = (RealOpenMM) _data[ii];
-      }
-      return _realOpenMMData;
-   } else {
-      return _data;
-   }
-*/
-   return NULL;
+   return getData( 0 );
 }
 
 /** 

platforms/brook/src/BrookFloatStreamInternal.h

@@ -104,10 +104,33 @@ class BrookFloatStreamInternal : public BrookStreamInternal {
 
       void saveToArray( void* array );
 
+      /** 
+       * Fill data w/ input value
+       * 
+       * @param  value to set array to
+       *
+       *
+       */
+
       void fillWithValue( void* value );
 
-      const RealOpenMM* const * getData( ) const;
-      RealOpenMM** getData( void );
+      /** 
+       * Get data array -- no read from board!
+       * 
+       * @return  data array
+       */
+      
+      void* getData( void );
+
+      /** 
+       * Get data array
+       * 
+       * @param readFromBoard if set, read values on board 
+       *
+       * @return  data array
+       */
+      
+      void* getData( int readFromBoard );
 
       /** 
        * Get dangle value

platforms/brook/src/BrookGbsa.cpp

@@ -32,6 +32,8 @@
 #include <sstream>
 #include "BrookGbsa.h"
 #include "BrookPlatform.h"
+#include "OpenMMException.h"
+#include "BrookStreamImpl.h"
 
 using namespace OpenMM;
 using namespace std;
@@ -57,14 +59,15 @@ BrookGbsa::BrookGbsa(  ){
    _partialForceStreamHeight  = -1;
    _partialForceStreamSize    = -1;
 
-   _gbsaStreamWidth           = -1;
-   _gbsaStreamHeight          = -1;
-   _gbsaStreamSize            = -1;
+   _gbsaAtomStreamWidth       = -1;
+   _gbsaAtomStreamHeight      = -1;
+   _gbsaAtomStreamSize        = -1;
 
    _duplicationFactor         =  4;
 
    _solventDielectric         = 78.3;
    _soluteDielectric          = 1.0;
+   _dielectricOffset          = 0.09;
 
    for( int ii = 0; ii < LastStreamIndex; ii++ ){
       _gbsaStreams[ii] = NULL;
@@ -74,6 +77,8 @@ BrookGbsa::BrookGbsa(  ){
       _gbsaForceStreams[ii] = NULL;
    }
 
+   _bornRadiiInitialized      = 0;
+   _cpuObc                    = NULL;
 }   
  
 /** 
@@ -97,6 +102,8 @@ BrookGbsa::~BrookGbsa( ){
       delete _gbsaForceStreams[ii];
    }
 
+   delete _cpuObc;
+
 }
 
 /** 
@@ -139,8 +146,8 @@ int BrookGbsa::getOuterLoopUnroll( void ) const {
  *
  */
 
-double BrookGbsa::getSoluteDielectric( void ) const {
-   return _soluteDielectric;
+float BrookGbsa::getSoluteDielectric( void ) const {
+   return (float) _soluteDielectric;
 }
 
 /** 
@@ -150,8 +157,19 @@ double BrookGbsa::getSoluteDielectric( void ) const {
  *
  */
 
-double BrookGbsa::getSolventDielectric( void ) const {
-   return _solventDielectric;
+float BrookGbsa::getSolventDielectric( void ) const {
+   return (float) _solventDielectric;
+}
+
+/** 
+ * Get OBC dielectric offset
+ * 
+ * @return  OBC dielectric offset
+ *
+ */
+
+float BrookGbsa::getDielectricOffset( void ) const {
+   return (float) _dielectricOffset;
 }
 
 /** 
@@ -208,31 +226,6 @@ int BrookGbsa::getDuplicationFactor( void ) const {
    return _duplicationFactor;
 }
 
-/** 
- * Get partial force stream width
- * 
- * @param platform  platform
- *
- * @return  partial force stream width
- *
- */
-
-int BrookGbsa::getPartialForceStreamWidth( const Platform& platform ){
-
-// ---------------------------------------------------------------------------------------
-
-   // static const std::string methodName      = "BrookGbsa::getPartialForceStreamWidth";
-
-// ---------------------------------------------------------------------------------------
-
-   // get partial force stream width
-
-   if( _partialForceStreamWidth < 0 ){
-      //_getPartialForceStreamDimensions( platform );
-   }
-   return _partialForceStreamWidth;
-}
-
 /** 
  * Get partial force stream width
  * 
@@ -252,31 +245,6 @@ int BrookGbsa::getPartialForceStreamWidth( void ) const {
    return _partialForceStreamWidth;
 }
 
-/** 
- * Get partial force stream height
- * 
- * @param platform platform
- *
- * @return  partial force stream height 
- *
- */
-
-int BrookGbsa::getPartialForceStreamHeight( const Platform& platform ){
-
-// ---------------------------------------------------------------------------------------
-
-   // static const std::string methodName      = "BrookGbsa::getPartialForceStreamHeight";
-
-// ---------------------------------------------------------------------------------------
-
-   // get partial force stream height
-
-   if( _partialForceStreamHeight < 0 ){
-      //_getPartialForceStreamDimensions( platform );
-   }
-   return _partialForceStreamHeight;
-}
-
 /** 
  * Get partial force stream height
  * 
@@ -298,59 +266,44 @@ int BrookGbsa::getPartialForceStreamHeight( void ) const {
 /** 
  * Get partial force stream size
  * 
- * @param platform  platform
- *
  * @return  partial force stream size
  *
  */
 
-int BrookGbsa::getPartialForceStreamSize( const Platform& platform ){
-
-// ---------------------------------------------------------------------------------------
-
-   // static const std::string methodName      = "BrookGbsa::getPartialForceStreamSize";
-
-// ---------------------------------------------------------------------------------------
-
-   // get partial force stream size
-
-   if( _partialForceStreamSize < 0 ){
-      //_getPartialForceStreamDimensions( platform );
-   }
+int BrookGbsa::getPartialForceStreamSize( void ) const {
    return _partialForceStreamSize;
 }
 
 /** 
- * Get partial force stream size
+ * Get Atom stream size
  *
- * @return  partial force stream size
+ * @return  Atom stream size
  *
  */
 
-int BrookGbsa::getPartialForceStreamSize( void ) const {
-   return _partialForceStreamSize;
+int BrookGbsa::getGbsaAtomStreamSize( void ) const {
+   return _gbsaAtomStreamSize;
 }
 
 /** 
- * Get Gbsa stream size
+ * Get atom stream width
  *
- * @return  Gbsa stream size
+ * @return  atom stream width
  *
  */
 
-int BrookGbsa::getGbsaStreamSize( void ) const {
-   return _gbsaStreamSize;
+int BrookGbsa::getGbsaAtomStreamWidth( void ) const {
+   return _gbsaAtomStreamWidth;
 }
 
 /** 
- * Get gbsa stream width
- *
- * @return  gbsa stream width
+ * Get atom stream height
  *
+ * @return atom stream height
  */
 
-int BrookGbsa::getGbsaStreamWidth( void ) const {
-   return _gbsaStreamWidth;
+int BrookGbsa::getGbsaAtomStreamHeight( void ) const {
+   return _gbsaAtomStreamHeight;
 }
 
 /** 
@@ -453,6 +406,80 @@ int BrookGbsa::isForceStreamSet( int index ) const {
    return (index >= 0 && index < getNumberOfForceStreams() && _gbsaForceStreams[index]) ? 1 : 0;
 }
 
+/** 
+ * Return true if Born radii have been initialized
+ *
+ * @return   true if Born radii have been initialized
+ *
+ */
+
+int BrookGbsa::haveBornRadiiBeenInitialized( void ) const {
+   return  _bornRadiiInitialized;
+}
+
+/** 
+ * Calculate Born radii
+ *
+ * @return  calculate Born radii
+ *
+ */
+
+int BrookGbsa::calculateBornRadii( const Stream& positions ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName                 = "BrookGbsa::calculateBornRadii";
+
+// ---------------------------------------------------------------------------------------
+
+   const BrookStreamImpl& positionStreamC              = dynamic_cast<const BrookStreamImpl&> (positions.getImpl());
+   BrookStreamImpl& positionStream                     = const_cast<BrookStreamImpl&>         (positionStreamC);
+   const RealOpenMM* coordinates                       = (RealOpenMM*) positionStream.getData( );
+
+   // load coordinates into RealOpenMM 2d array
+
+   int numberOfAtoms                                   = getNumberOfAtoms();
+
+   RealOpenMM** atomCoordinates                        = new RealOpenMM*[numberOfAtoms];
+   RealOpenMM* atomCoordinatesBlk                      = new RealOpenMM[3*numberOfAtoms];
+
+   // Born radii array size needs to match stream size since it will
+   // be written down to board
+
+   int streamSize                                      = getGbsaAtomStreamSize();
+   RealOpenMM* bornRadii                               = new RealOpenMM[streamSize];
+   memset( bornRadii, 0, sizeof( RealOpenMM )*streamSize );
+
+   int index = 0;
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+      atomCoordinates[ii]    = atomCoordinatesBlk;
+      atomCoordinatesBlk    += 3;
+
+      atomCoordinates[ii][0] = coordinates[index++];
+      atomCoordinates[ii][1] = coordinates[index++];
+      atomCoordinates[ii][2] = coordinates[index++];
+
+   }
+
+   // calculate Born radii
+
+   _cpuObc->computeBornRadii( atomCoordinates, bornRadii );
+
+   // write radii to board and set flag to indicate radii calculated once
+
+   _gbsaStreams[ObcBornRadiiStream]->loadFromArray( bornRadii );
+   _bornRadiiInitialized   = 1;
+
+   // free memory
+
+   delete[] atomCoordinatesBlk;
+   delete[] atomCoordinates;
+   delete[] bornRadii;
+
+   return DefaultReturnValue;
+}
+
 /** 
  * Initialize stream dimensions
  * 
@@ -468,17 +495,12 @@ int BrookGbsa::initializeStreamSizes( int numberOfAtoms, const Platform& platfor
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName      = "BrookGbsa::initializeStreamSizes";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
-   int atomStreamSize                  = getAtomStreamSize( platform );
-   int atomStreamWidth                 = getAtomStreamWidth( platform );
-
-   initializeExclusionStreamSize( atomStreamSize, atomStreamWidth );
-   initializeJStreamSize( atomStreamSize, atomStreamWidth );
-   initializeOuterVdwStreamSize( atomStreamSize, atomStreamWidth );
-   initializePartialForceStreamSize( atomStreamSize, atomStreamWidth );
+   _gbsaAtomStreamSize     = getAtomStreamSize( platform );
+   _gbsaAtomStreamWidth    = getAtomStreamWidth( platform );
+   _gbsaAtomStreamHeight   = getAtomStreamHeight( platform );
 
    return DefaultReturnValue;
 }
@@ -498,39 +520,57 @@ int BrookGbsa::initializeStreams( const Platform& platform ){
 
    static const std::string methodName      = "BrookGbsa::initializeStreams";
    static const double dangleValue          = 0.0;
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
-   const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (platform.getDefaultStreamFactory());
+   int gbsaAtomStreamSize   = getGbsaAtomStreamSize();
+   int gbsaAtomStreamWidth  = getGbsaAtomStreamWidth();
+
+   // atomic radii & charge
 
-   // exclusion
+    _gbsaStreams[ObcAtomicRadiiStream]                          = new BrookFloatStreamInternal( BrookCommon::ObcAtomicRadiiStream,
+                                                                                                gbsaAtomStreamSize, gbsaAtomStreamWidth,
+                                                                                                BrookStreamInternal::Float2, dangleValue );
 
-    _gbsaStreams[ExclusionStream]          = new BrookFloatStreamInternal( BrookStreamFactory::NonBondedExclusionStream,
-                                                                                getExclusionStreamSize(), getExclusionStreamWidth(),
+   // scaled atomic radii
+
+    _gbsaStreams[ObcScaledAtomicRadiiStream]                    = new BrookFloatStreamInternal( BrookCommon::ObcScaledAtomicRadiiStream,
+                                                                                                gbsaAtomStreamSize, gbsaAtomStreamWidth,
+                                                                                                BrookStreamInternal::Float2, dangleValue );
+
+   // atomic radii w/ DielectricOffset
+
+    _gbsaStreams[ObcAtomicRadiiWithDielectricOffsetStream]      = new BrookFloatStreamInternal( BrookCommon::ObcAtomicRadiiWithDielectricOffsetStream,
+                                                                                                gbsaAtomStreamSize, gbsaAtomStreamWidth,
                                                                                                 BrookStreamInternal::Float, dangleValue );
 
-   // outer vdw
+   // Born radii
+
+    _gbsaStreams[ObcBornRadiiStream]                            = new BrookFloatStreamInternal( BrookCommon::ObcBornRadiiStream,
+                                                                                                gbsaAtomStreamSize, gbsaAtomStreamWidth,
+                                                                                                BrookStreamInternal::Float, dangleValue );
 
-   _gbsaStreams[OuterVdwStream]            = new BrookFloatStreamInternal( BrookStreamFactory::OuterVdwStream, getAtomStreamSize(), 
-                                                                                getAtomStreamWidth(), BrookStreamInternal::Float2, dangleValue );
+   // Born2 radii
 
-   // inner sigma & epsilon
+    _gbsaStreams[ObcBornRadii2Stream]                           = new BrookFloatStreamInternal( BrookCommon::ObcBornRadii2Stream,
+                                                                                                gbsaAtomStreamSize, gbsaAtomStreamWidth,
+                                                                                                BrookStreamInternal::Float, dangleValue );
 
-   _gbsaStreams[InnerSigmaStream]         = new BrookFloatStreamInternal( BrookStreamFactory::InnerSigmaStream, getJStreamSize(), 
-                                                                               getJStreamWidth(), BrookStreamInternal::Float4, dangleValue );
+   // IntermediateForce
 
-   _gbsaStreams[InnerEpsilonStream]       = new BrookFloatStreamInternal( BrookStreamFactory::InnerEpsilonStream, getJStreamSize(),
-                                                                               getJStreamWidth(), BrookStreamInternal::Float4, dangleValue );
+    _gbsaStreams[ObcIntermediateForceStream]                    = new BrookFloatStreamInternal( BrookCommon::ObcIntermediateForceStream,
+                                                                                                gbsaAtomStreamSize, gbsaAtomStreamWidth,
+                                                                                                BrookStreamInternal::Float, dangleValue );
 
-   // charge stream
+   // Born2 radii
 
-   _gbsaStreams[ChargeStream]              = new BrookFloatStreamInternal( BrookStreamFactory::NonBondedChargeStream, getAtomStreamSize(),
-                                                                                getAtomStreamWidth(), BrookStreamInternal::Float, dangleValue );
+    _gbsaStreams[ObcChainStream]                                = new BrookFloatStreamInternal( BrookCommon::ObcChainStream,
+                                                                                                gbsaAtomStreamSize, gbsaAtomStreamWidth,
+                                                                                                BrookStreamInternal::Float, dangleValue );
 
+   // partial force streams
 
-   // partial force stream
-   std::string partialForceStream = BrookStreamFactory::PartialForceStream;
+   std::string partialForceStream = BrookCommon::PartialForceStream;
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
       std::stringstream name;
       name << partialForceStream << ii;
@@ -544,34 +584,116 @@ int BrookGbsa::initializeStreams( const Platform& platform ){
 /*  
  * Setup of Gbsa parameters
  *
- * @param atomParameters        vector of OBC parameters [atomI][0=index]
- *                                                       [atomI][1=charge]
- *                                                       [atomI][2=radius]
+ * @param atomParameters        vector of OBC parameters [atomI][0=charge]
+ *                                                       [atomI][1=radius]
  *                                                       [atomI][2=scaling factor]
  * @param solventDielectric     solvent dielectric
  * @param soluteDielectric      solute dielectric
+ * @param platform              Brook platform
  *
  * @return nonzero value if error
  *
  * */
     
-int BrookGbsa::setup( const std::vector<std::vector<double> >& atomParameters, 
-                      double solventDielectric, double soluteDielectric );
+int BrookGbsa::setup( const std::vector<std::vector<double> >& vectorOfAtomParameters, 
+                      double solventDielectric, double soluteDielectric, const Platform& platform ){
     
 // ---------------------------------------------------------------------------------------
 
+   static const int atomParametersSize      = 4; 
+   static const int maxErrors               = 20; 
    static const std::string methodName      = "BrookGbsa::setup";
 
 // ---------------------------------------------------------------------------------------
 
-   setNumberOfAtoms( (int) atomParameters.size() );
+   int numberOfAtoms  = (int) vectorOfAtomParameters.size();
+   setNumberOfAtoms( numberOfAtoms );
+
+   _solventDielectric = solventDielectric;
+   _soluteDielectric  = soluteDielectric;
 
-   solventDielectric = _solventDielectric;
-   soluteDielectric  = _soluteDielectric;
+   // initialize stream sizes and then Brook streams
 
-   initializeStreamSizes( getNumberOfAtoms(), platform );
+   initializeStreamSizes( numberOfAtoms, platform );
    initializeStreams( platform );
 
+   BrookOpenMMFloat* radiiAndCharge         = new BrookOpenMMFloat[getNumberOfAtoms()*2];
+   BrookOpenMMFloat* scaledRadiiAndOffset   = new BrookOpenMMFloat[getNumberOfAtoms()*2];
+
+   // used by CpuObc to calculate initial Born radii
+
+   vector<RealOpenMM> atomicRadii(numberOfAtoms);
+   vector<RealOpenMM> scaleFactors(numberOfAtoms);
+
+   float dielectricOffset                  = getDielectricOffset();
+
+   // loop over atom parameters
+   // track any errors and then throw exception
+   //    check parameter vector is right size
+   // set parameter entries or board and arrays used by CpuObc
+
+   int vectorIndex  = 0;
+   int errors       = 0;
+   std::stringstream message;
+
+   typedef std::vector< std::vector<double> > VectorOfDoubleVectors;
+   typedef VectorOfDoubleVectors::const_iterator VectorOfDoubleVectorsCI;
+
+   for( VectorOfDoubleVectorsCI ii = vectorOfAtomParameters.begin(); ii != vectorOfAtomParameters.end(); ii++ ){
+
+      std::vector<double> atomParameters = *ii;
+
+      if( atomParameters.size() != atomParametersSize && errors < maxErrors ){
+         message << methodName << " parameter size=" << atomParameters.size() << " for parameter vector index=" << vectorIndex << " is less than expected.\n";
+         errors++;
+      } else {
+
+         double charge                            = atomParameters[0];     
+         double radius                            = atomParameters[1];     
+         double scalingFactor                     = atomParameters[2];     
+
+         int streamIndex                          = 2*vectorIndex;
+
+         atomicRadii[vectorIndex]                 = static_cast<RealOpenMM> (radius);
+         scaleFactors[vectorIndex]                = static_cast<RealOpenMM> (scalingFactor);
+
+         radiiAndCharge[streamIndex]              = static_cast<BrookOpenMMFloat> (radius);
+         radiiAndCharge[streamIndex+1]            = static_cast<BrookOpenMMFloat> (charge);
+
+         scaledRadiiAndOffset[streamIndex]        = static_cast<BrookOpenMMFloat> (radius*scalingFactor);
+         scaledRadiiAndOffset[streamIndex+1]      = static_cast<BrookOpenMMFloat> (radius - dielectricOffset);
+
+      }
+
+      vectorIndex++;
+   }
+
+   // throw exception if errors detected
+
+   if( errors ){
+      throw OpenMMException( message.str() );
+   }
+
+   // load streams
+
+   _gbsaStreams[ObcAtomicRadiiStream]->loadFromArray( radiiAndCharge );
+   _gbsaStreams[ObcScaledAtomicRadiiStream]->loadFromArray( scaledRadiiAndOffset );
+
+   delete[] radiiAndCharge;
+   delete[] scaledRadiiAndOffset;
+
+   // setup for Born radii
+
+   ObcParameters* obcParameters  = new ObcParameters( numberOfAtoms, ObcParameters::ObcTypeII );
+   obcParameters->setAtomicRadii( atomicRadii, SimTKOpenMMCommon::MdUnits);
+
+   obcParameters->setScaledRadiusFactors( scaleFactors );
+   obcParameters->setSolventDielectric( static_cast<RealOpenMM>(solventDielectric) );
+   obcParameters->setSoluteDielectric(  static_cast<RealOpenMM>(soluteDielectric)  );
+
+   _cpuObc  = new CpuObc(obcParameters);
+   _cpuObc->setIncludeAceApproximation( true );
+
    return DefaultReturnValue;
 }
 
@@ -686,9 +808,6 @@ std::string BrookGbsa::getContentsString( int level ) const {
    (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 ); 
-
    message << _getLine( tab, "Log:",                  (getLog()                ? Set : NotSet) ); 
 /*
    message << _getLine( tab, "ExclusionStream:",      (getExclusionStream()    ? Set : NotSet) ); 

platforms/brook/src/BrookGbsa.h

@@ -38,6 +38,7 @@
 #include "BrookFloatStreamInternal.h"
 #include "BrookPlatform.h"
 #include "BrookCommon.h"
+#include "../../../platforms/reference/src/gbsa/CpuObc.h"
 
 namespace OpenMM {
 
@@ -53,8 +54,18 @@ class BrookGbsa : public BrookCommon {
       static const int DefaultReturnValue = 0;
       static const int ErrorReturnValue   = -1;
 
+      /** 
+       * Constructor
+       * 
+       */
+      
       BrookGbsa( );
   
+      /** 
+       * Destructor
+       * 
+       */
+      
       ~BrookGbsa();
   
       /**
@@ -144,6 +155,30 @@ class BrookGbsa : public BrookCommon {
 
       int getPartialForceStreamSize( void ) const; 
 
+      /**
+       * Get Gbsa atom stream width
+       *
+       * @return atom stream width
+       */
+
+      int getGbsaAtomStreamWidth( void ) const; 
+
+      /**
+       * Get Gbsa atom stream height
+       *
+       * @return atom stream height
+       */
+
+      int getGbsaAtomStreamHeight( void ) const;
+
+      /**
+       * Get Gbsa atom stream size
+       * 
+       * @return atom stream size
+       */
+
+      int getGbsaAtomStreamSize( void ) const; 
+
       /**
        * Get solute dielectric
        * 
@@ -160,6 +195,14 @@ class BrookGbsa : public BrookCommon {
 
       float getSolventDielectric( void ) const; 
 
+      /**
+       * Get OBC dielectric offset
+       * 
+       * @return dielectric offset
+       */
+
+      float getDielectricOffset( void ) const; 
+
       /** 
        * Get atomic radii
        *
@@ -232,6 +275,15 @@ class BrookGbsa : public BrookCommon {
       
       BrookFloatStreamInternal** getForceStreams( void );
       
+      /** 
+       * Get array of Gbsa streams 
+       *
+       * @return  array ofstreams
+       *
+       */
+      
+      BrookFloatStreamInternal** getStreams( void );
+      
       /** 
        * Return true if force[index] stream is set 
        *
@@ -241,22 +293,40 @@ class BrookGbsa : public BrookCommon {
       
       int isForceStreamSet( int index ) const;
       
+      /** 
+       * Return true if Born radii have been initialized
+       *
+       * @return  true  if Born radii have been initialized
+       *
+       */
+      
+      int haveBornRadiiBeenInitialized( void ) const;
+      
+      /** 
+       * Calculate Born radii
+       *
+       * @return  calculate Born radii
+       *
+       */
+      
+      int BrookGbsa::calculateBornRadii( const Stream& positions );
+      
       /* 
        * Setup of Gbsa parameters
        *
-       * @param atomParameters        vector of OBC parameters [atomI][0=index]
-       *                                                       [atomI][1=charge]
-       *                                                       [atomI][2=radius]
+       * @param atomParameters        vector of OBC parameters [atomI][0=charge]
+       *                                                       [atomI][1=radius]
        *                                                       [atomI][2=scaling factor]
        * @param solventDielectric     solvent dielectric
        * @param soluteDielectric      solute dielectric
+       * @param platform              Brook platform
        *
        * @return nonzero value if error
        *
        * */
       
       int setup( const std::vector<std::vector<double> >& atomParameters, 
-                 double solventDielectric, double soluteDielectric );
+                 double solventDielectric, double soluteDielectric, const Platform& platform  );
       
       /* 
        * Get contents of object
@@ -279,7 +349,6 @@ class BrookGbsa : public BrookCommon {
 
       enum { 
               ObcAtomicRadiiStream,
-              ObcRadiiStreamStream,
               ObcScaledAtomicRadiiStream,
               ObcAtomicRadiiWithDielectricOffsetStream,
               ObcBornRadiiStream,
@@ -308,22 +377,32 @@ class BrookGbsa : public BrookCommon {
       int _partialForceStreamHeight;
       int _partialForceStreamSize;
 
-      // Gbsa stream dimensions
+      // Atom stream dimensions
 
-      int _gbsaStreamWidth;
-      int _gbsaStreamHeight;
-      int _gbsaStreamSize;
+      int _gbsaAtomStreamWidth;
+      int _gbsaAtomStreamHeight;
+      int _gbsaAtomStreamSize;
 
       // dielectrics
 
-      double _solventDielectric
-      double _soluteDielectric
+      double _solventDielectric;
+      double _soluteDielectric;
+
+      // dielectric offset
+
+      double _dielectricOffset;
 
       // internal streams
 
       BrookFloatStreamInternal* _gbsaStreams[LastStreamIndex];
       BrookFloatStreamInternal* _gbsaForceStreams[NumberOfForceStreams];
 
+      int _bornRadiiInitialized;
+
+      // CpuObc reference
+
+      CpuObc* _cpuObc;
+
       /* 
        * Setup of stream dimensions
        *

platforms/brook/src/BrookIntStreamInternal.cpp

@@ -39,18 +39,21 @@ using namespace OpenMM;
  * BrookIntStreamInternal constructor
  * 
  * @param name                      stream name
- * @param size                      stream size
- * @param platform                  platform
+ * @param size                      array size
+ * @param streamWidth               stream width
+ * @param type                      stream type (Integer, Integer2, ...)
+ * @param dangleValue               fill value for tail of stream beyond array size
  *
  */
 
-BrookIntStreamInternal::BrookIntStreamInternal( std::string name, int size, int width, BrookStreamInternal::DataType type,
-                                                int dangleValue ) : BrookStreamInternal( name, size, width, type ){
+BrookIntStreamInternal::BrookIntStreamInternal( std::string name, int size, int streamWidth,
+                                                BrookStreamInternal::DataType type,
+                                                int dangleValue ) : 
+                        BrookStreamInternal( name, size, streamWidth, type ){
 
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName      = "BrookIntStreamInternal::BrookIntStreamInternal";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
@@ -59,32 +62,34 @@ BrookIntStreamInternal::BrookIntStreamInternal( std::string name, int size, int
    switch( type ){
 
       case BrookStreamInternal::Integer:
-          width = 1;
+
+          _width = 1;
           break;
 
       case BrookStreamInternal::Integer2:
-          width = 2;
+
+          _width = 2;
           break;
 
       case BrookStreamInternal::Integer3:
-          width = 3;
+
+          _width = 3;
           break;
 
       case BrookStreamInternal::Integer4:
-          width = 4;
+
+          _width = 4;
           break;
 
       default:
+
          std::stringstream message;
          message << methodName << " type=" << type << " not recognized.";
          throw OpenMMException( message.str() );
    }
 
-   data = new int*[size];
+   _data = new int[size*_width];
 
-   for( int ii = 0; ii < size; ii++ ){
-       data[ii] = new int[width];
-   }
 }
 
 /** 
@@ -93,7 +98,7 @@ BrookIntStreamInternal::BrookIntStreamInternal( std::string name, int size, int
  */
 
 BrookIntStreamInternal::~BrookIntStreamInternal() {
-   delete[] data;
+   delete[] _data;
 }
 
 /** 
@@ -127,7 +132,6 @@ void BrookIntStreamInternal::loadFromArray( const void* array, BrookStreamIntern
 // ---------------------------------------------------------------------------------------
 
    static const std::string methodName      = "BrookIntStreamInternal::loadFromArray(1)";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
@@ -138,11 +142,10 @@ void BrookIntStreamInternal::loadFromArray( const void* array, BrookStreamIntern
    }
 
    int* arrayData = (int*) array;
-   int index      = 0;
-   for( int ii = 0; ii < getSize(); ii++ ){
-      for( int jj = 0; jj < getWidth(); jj++ ){
-         data[ii][jj] = arrayData[index++];
-      }
+   int totalSize  = getSize()*getWidth();
+
+   for( int ii = 0; ii < totalSize; ii++ ){
+      _data[ii] = arrayData[ii];
    }
 }
 
@@ -158,16 +161,13 @@ void BrookIntStreamInternal::saveToArray( void* array ){
 // ---------------------------------------------------------------------------------------
 
    // static const std::string methodName      = "BrookIntStreamInternal::saveToArray";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
    int* arrayData = (int*) array;
-   int index      = 0;
-   for( int ii = 0; ii < getSize(); ii++ ){
-      for( int jj = 0; jj < getWidth(); jj++ ){
-         arrayData[index++] = data[ii][jj];
-      }
+   int totalSize  = getSize()*getWidth();
+   for( int ii = 0; ii < totalSize; ii++ ){
+      arrayData[ii] = _data[ii];
    }
 }
 
@@ -183,27 +183,15 @@ void BrookIntStreamInternal::fillWithValue( void* value ){
 // ---------------------------------------------------------------------------------------
 
    // static const std::string methodName      = "BrookIntStreamInternal::fillWithValue";
-   // static const int debug                   = 1;
 
 // ---------------------------------------------------------------------------------------
 
    int valueData = *((int*) value);
-   for( int ii = 0; ii < getSize(); ii++ ){
-      for (int jj = 0; jj < getWidth(); jj++ ){
-         data[ii][jj] = valueData;
-      }
-   }
-}
-
-/** 
- * Get data
- * 
- * @return data ptr
- *
- */
+   int totalSize  = getSize()*getWidth();
 
-const int* const * BrookIntStreamInternal::getData( void ) const {
-   return data;
+   for( int ii = 0; ii < totalSize; ii++ ){
+         _data[ii] = valueData;
+   }
 }
 
 /** 
@@ -213,7 +201,7 @@ const int* const * BrookIntStreamInternal::getData( void ) const {
  *
  */
 
-int** BrookIntStreamInternal::getData( void ){
-   return data;
+void* BrookIntStreamInternal::getData( void ){
+   return _data;
 }
 

platforms/brook/src/BrookIntStreamInternal.h

@@ -45,10 +45,36 @@ class BrookIntStreamInternal : public BrookStreamInternal {
 
 public:
 
+      /** 
+       * BrookIntStreamInternal constructor
+       * 
+       * @param name                      stream name
+       * @param size                      size of array
+       * @param streamWidth               stream width
+       * @param type                      stream type (float, float2, ...)
+       * @param inputDefaultDangleValue   default dangle value
+       *
+       */
+
       BrookIntStreamInternal( std::string name, int size, int streamWidth, BrookStreamInternal::DataType type, int dangleValue );
 
+      /** 
+       * BrookIntStreamInternal destructor
+       * 
+       */
+
       ~BrookIntStreamInternal( );
 
+      /** 
+       * Copy the contents of an array into this stream.
+       * 
+       * @param  array a pointer to the start of the array.  The array is assumed to have the same length as this stream,
+       * and to contain elements of the correct data type for this stream.  If the stream has a compound data type, all
+       * the values should be packed into a single array: all the values for the first element, followed by all the values
+       * for the next element, etc.
+       *
+       */
+
       void loadFromArray( const void* array );
 
       /** 
@@ -58,19 +84,50 @@ public:
        * and to contain elements of the correct data type for this stream.  If the stream has a compound data type, all
        * the values should be packed into a single array: all the values for the first element, followed by all the values
        * for the next element, etc.
+       *
+       * @param baseType data type of input array (float, double, int)
+       *
        */
+
       void loadFromArray( const void* array, BrookStreamInternal::DataType baseType );
 
+      /** 
+       * Save data to input array
+       * 
+       * @param  array a pointer to the start of the array.  The array is assumed to have the same length as this stream,
+       * and to contain elements of the correct _data type for this stream.  If the stream has a compound _data type, all
+       * the values should be packed into a single array: all the values for the first element, followed by all the values
+       * for the next element, etc.
+       *
+       * @throw exception if baseType not float or double
+       *
+       */
+
       void saveToArray( void* array );
+
+      /** 
+       * Fill data w/ input value
+       * 
+       * @param  value to set array to
+       *
+       *
+       */
+
       void fillWithValue( void* value );
-    const int* const * getData( void  ) const;
-    int** getData( void );
+
+      /** 
+       * Get data array
+       * 
+       * @return  data array
+       */
+
+      void* getData( void );
 
 private:
 
     int _dangleValue;
 
-    int** data;
+    int* _data;
 };
 
 } // namespace OpenMM

platforms/brook/src/BrookNonBonded.cpp

@@ -516,31 +516,31 @@ int BrookNonBonded::initializeStreams( const Platform& platform ){
 
    // exclusion
 
-    _nonbondedStreams[ExclusionStream]          = new BrookFloatStreamInternal( BrookStreamFactory::NonBondedExclusionStream,
+    _nonbondedStreams[ExclusionStream]          = new BrookFloatStreamInternal( BrookCommon::NonBondedExclusionStream,
                                                                                 getExclusionStreamSize(), getExclusionStreamWidth(),
                                                                                 BrookStreamInternal::Float, dangleValue );
 
    // outer vdw
 
-   _nonbondedStreams[OuterVdwStream]            = new BrookFloatStreamInternal( BrookStreamFactory::OuterVdwStream, getAtomStreamSize(), 
+   _nonbondedStreams[OuterVdwStream]            = new BrookFloatStreamInternal( BrookCommon::OuterVdwStream, getAtomStreamSize(), 
                                                                                 getAtomStreamWidth(), BrookStreamInternal::Float2, dangleValue );
 
    // inner sigma & epsilon
 
-   _nonbondedStreams[InnerSigmaStream]         = new BrookFloatStreamInternal( BrookStreamFactory::InnerSigmaStream, getJStreamSize(), 
+   _nonbondedStreams[InnerSigmaStream]         = new BrookFloatStreamInternal( BrookCommon::InnerSigmaStream, getJStreamSize(), 
                                                                                getJStreamWidth(), BrookStreamInternal::Float4, dangleValue );
 
-   _nonbondedStreams[InnerEpsilonStream]       = new BrookFloatStreamInternal( BrookStreamFactory::InnerEpsilonStream, getJStreamSize(),
+   _nonbondedStreams[InnerEpsilonStream]       = new BrookFloatStreamInternal( BrookCommon::InnerEpsilonStream, getJStreamSize(),
                                                                                getJStreamWidth(), BrookStreamInternal::Float4, dangleValue );
 
    // charge stream
 
-   _nonbondedStreams[ChargeStream]              = new BrookFloatStreamInternal( BrookStreamFactory::NonBondedChargeStream, getAtomStreamSize(),
+   _nonbondedStreams[ChargeStream]              = new BrookFloatStreamInternal( BrookCommon::NonBondedChargeStream, getAtomStreamSize(),
                                                                                 getAtomStreamWidth(), BrookStreamInternal::Float, dangleValue );
 
    
    // partial force stream
-   std::string partialForceStream = BrookStreamFactory::PartialForceStream;
+   std::string partialForceStream = BrookCommon::PartialForceStream;
    for( int ii = 0; ii < getNumberOfForceStreams(); ii++ ){
       std::stringstream name;
       name << partialForceStream << ii;

platforms/brook/src/BrookStreamFactory.cpp

@@ -40,23 +40,6 @@ const std::string BrookStreamFactory::AtomPositions              = "atomPosition
 const std::string BrookStreamFactory::AtomVelocities             = "atomVelocities";
 const std::string BrookStreamFactory::AtomForces                 = "atomForces";
 
-// bonded streams
-                                    
-const std::string BrookStreamFactory::BondedAtomIndicesStream    = "BondedAtomIndicesStream";
-const std::string BrookStreamFactory::BondedParametersStream     = "BondedParametersStream";
-const std::string BrookStreamFactory::UnrolledForceStream        = "UnrolledForceStream";
-const std::string BrookStreamFactory::BondedChargeStream         = "BondedChargeStream";
-const std::string BrookStreamFactory::BondedInverseMapStreams    = "BondedInverseMapStreams";
-
-// non-bonded streams
-
-const std::string BrookStreamFactory::NonBondedExclusionStream   = "NonBondedExclusionStream";
-const std::string BrookStreamFactory::OuterVdwStream             = "OuterVdwStream";
-const std::string BrookStreamFactory::InnerSigmaStream           = "InnerSigmaStream";
-const std::string BrookStreamFactory::InnerEpsilonStream         = "InnerEpsilonStream";
-const std::string BrookStreamFactory::NonBondedChargeStream      = "NonBondedChargeStream";
-const std::string BrookStreamFactory::PartialForceStream         = "PartialForceStream";
-
 const double DefaultDangleValue                                  = 1.0e+38;
 /** 
  * BrookStreamFactory constructor

platforms/brook/src/BrookStreamImpl.cpp

@@ -262,6 +262,15 @@ void BrookStreamImpl::fillWithValue( void* value ){
    return _brookStreamInternal->fillWithValue( value );
 }
 
+/** 
+ * Set every element of this stream to the same value.
+ * 
+ * @return data array
+ */
+void* BrookStreamImpl::getData( void* value ){
+   return _brookStreamInternal->getData( );
+}
+
 /** 
  * Get Brook stream
  * 

platforms/brook/src/BrookStreamImpl.h

@@ -108,6 +108,13 @@ class BrookStreamImpl : public StreamImpl {
        */
       void fillWithValue( void* value );
   
+      /**
+       * Get data array
+       * 
+       * @return data array
+       */
+      void* getData( void );
+  
       /**
        * Get Brook stream
        * 

platforms/brook/src/BrookStreamInternal.h

@@ -130,6 +130,13 @@ class BrookStreamInternal {
        */
       virtual void fillWithValue(void* value) = 0;
   
+      /**
+       * Get data
+       * 
+       * @return data array
+       */
+      virtual void* getData( void ) = 0;
+  
       /** 
        * Get type string
        *