Name changes and minor mods to sync w/ new interface

platforms/brook/src/BrookCalcPeriodicTorsionForceKernel.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "OpenMMException.h"
+#include <sstream>
+#include "BrookCalcPeriodicTorsionForceKernel.h"
+
+using namespace OpenMM;
+using namespace std;
+
+const std::string BrookCalcPeriodicTorsionForceKernel::BondName = "ProperDihedral";
+
+/** 
+ * BrookCalcPeriodicTorsionForceKernel constructor
+ * 
+ * @param name                      kernel name
+ * @param platform                  platform
+ * @param OpenMMBrookInterface      OpenMM-Brook interface
+ * @param System                    System reference
+ *
+ */
+
+BrookCalcPeriodicTorsionForceKernel::BrookCalcPeriodicTorsionForceKernel( std::string name, const Platform& platform,
+                                                                          OpenMMBrookInterface& openMMBrookInterface, System& system ) :
+                     CalcPeriodicTorsionForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookCalcPeriodicTorsionForceKernel::BrookCalcPeriodicTorsionForceKernel";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   _brookBondParameters              = NULL;
+   _log                              = NULL;
+
+   const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
+   if( brookPlatform.getLog() != NULL ){
+      setLog( brookPlatform.getLog() );
+   }
+      
+}   
+
+/** 
+ * BrookCalcPeriodicTorsionForceKernel destructor
+ * 
+ */
+
+BrookCalcPeriodicTorsionForceKernel::~BrookCalcPeriodicTorsionForceKernel( ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookCalcPeriodicTorsionForceKernel::BrookCalcPeriodicTorsionForceKernel";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   delete _brookBondParameters;
+}
+
+/** 
+ * Get log file reference
+ * 
+ * @return  log file reference
+ *
+ */
+
+FILE* BrookCalcPeriodicTorsionForceKernel::getLog( void ) const {
+   return _log;
+}
+
+/** 
+ * Set log file reference
+ * 
+ * @param  log file reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int BrookCalcPeriodicTorsionForceKernel::setLog( FILE* log ){
+   _log = log;
+   return BrookCommon::DefaultReturnValue;
+}
+
+/** 
+ * Initialize the kernel, setting up the values of all the force field parameters.
+ * 
+ * @param system                    System reference
+ * @param force                     ProperDihedralForce reference
+ *
+ */
+
+void BrookCalcPeriodicTorsionForceKernel::initialize( const System& system, const PeriodicTorsionForce& force ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookCalcPeriodicTorsionForceKernel::initialize";
+
+// ---------------------------------------------------------------------------------------
+
+   FILE* log                 = getLog();
+
+   // ---------------------------------------------------------------------------------------
+
+   // create _brookBondParameters object containing atom indices/parameters
+
+   int numberOfBonds         = force.getNumTorsions();
+
+   if( _brookBondParameters ){
+      delete _brookBondParameters;
+   }
+   _brookBondParameters = new BrookBondParameters( BondName, NumberOfAtomsInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
+
+   for( int ii = 0; ii < numberOfBonds; ii++ ){
+
+      int particle1, particle2, particle3, particle4, periodicity;
+      double phase, k;
+
+      int particles[NumberOfAtomsInBond];
+      double parameters[NumberOfParametersInBond];
+
+      force.getTorsionParameters( ii, particle1, particle2, particle3, particle4, periodicity, phase, k ); 
+      particles[0]    = particle1;
+      particles[1]    = particle2;
+      particles[2]    = particle3;
+      particles[3]    = particle4;
+ 
+      parameters[0]   = (double) periodicity;
+      parameters[1]   = phase;
+      parameters[2]   = k;
+
+      _brookBondParameters->setBond( ii, particles, parameters );
+   }   
+   _openMMBrookInterface.setPeriodicTorsionForceParameters( _brookBondParameters );
+   _openMMBrookInterface.setTriggerForceKernel( this );
+   _openMMBrookInterface.setTriggerEnergyKernel( this );
+
+   if( log ){
+      std::string contents = _brookBondParameters->getContentsString( ); 
+      (void) fprintf( log, "%s brookGbsa::contents\n%s", methodName.c_str(), contents.c_str() );
+      (void) fflush( log );
+   }
+
+   // ---------------------------------------------------------------------------------------
+    
+}
+
+/** 
+ * Compute forces given atom coordinates
+ * 
+ * @param context OpenMMContextImpl context
+ *
+ */
+
+void BrookCalcPeriodicTorsionForceKernel::executeForces( OpenMMContextImpl& context ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName   = "BrookCalcPeriodicTorsionForceKernel::executeForces";
+
+// ---------------------------------------------------------------------------------------
+
+   if( _openMMBrookInterface.getTriggerForceKernel() == this ){
+      _openMMBrookInterface.computeForces( context );
+   }
+
+   return;
+
+   // ---------------------------------------------------------------------------------------
+}
+
+/**
+ * Execute the kernel to calculate the energy
+ * 
+ * @param context OpenMMContextImpl context
+ *
+ * @return  potential energy
+ *
+ */
+
+double BrookCalcPeriodicTorsionForceKernel::executeEnergy( OpenMMContextImpl& context ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookCalcPeriodicTorsionForceKernel::executeEnergy";
+
+// ---------------------------------------------------------------------------------------
+
+   if( _openMMBrookInterface.getTriggerEnergyKernel() == this ){
+      return (double) _openMMBrookInterface.computeEnergy( context );
+   } else {
+      return 0.0;
+   }
+
+}

platforms/brook/src/BrookCalcPeriodicTorsionForceKernel.h

+#ifndef OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_
+#define OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "kernels.h"
+#include "BrookPlatform.h"
+#include "BrookBondParameters.h"
+#include "OpenMMBrookInterface.h"
+
+namespace OpenMM {
+
+/**
+ * This kernel is invoked to calculate the periodic torsion forces acting on the system.
+ */
+
+class BrookCalcPeriodicTorsionForceKernel : public CalcPeriodicTorsionForceKernel {
+
+   public:
+  
+      /**
+       * BrookCalcPeriodicTorsionForceKernel constructor
+       */
+
+      BrookCalcPeriodicTorsionForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
+  
+      /**
+       * BrookCalcPeriodicTorsionForceKernel destructor
+       */
+
+      ~BrookCalcPeriodicTorsionForceKernel();
+  
+      /**
+       * Initialize the kernel, setting up the values to calculate harmonic bond force & energy
+       * 
+       * @param system                    System reference
+       * @param force                     ProperDihedralForce reference
+       *
+       */
+
+      void initialize( const System& system, const PeriodicTorsionForce& force );
+  
+      /**
+       * Execute the kernel to calculate the forces.
+       * 
+       * @param positions   atom coordiantes
+       * @param forces      output forces
+       *
+       */
+
+      void executeForces( OpenMMContextImpl& context );
+  
+      /**
+       * Execute the kernel to calculate the energy.
+       * 
+       * @param context    the context in which to execute this kernel
+       *
+       * @return  potential energy associated with the harmonic angle force
+       *
+       */
+
+      double executeEnergy( OpenMMContextImpl& context );
+
+      /** 
+       * Set log file reference
+       * 
+       * @param  log file reference
+       *
+       * @return DefaultReturnValue
+       *
+       */
+      
+      int setLog( FILE* log );
+
+      /* 
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string getContents( int level ) const;
+
+      /** 
+       * Get log file reference
+       * 
+       * @return  log file reference
+       *
+       */
+      
+      FILE* getLog( void ) const;
+      
+      /** 
+       * Get number of bonds
+       * 
+       * @return  number of bonds
+       *
+       */
+      
+      int getNumberOfBonds( void ) const;
+      
+      /** 
+       * Get indices/parameters
+       * 
+       * @return  BrookBondParameters containing atom indices/parameters
+       *
+       */
+      
+      BrookBondParameters* getBrookBondParameters( void ) const;
+      
+   private:
+   
+      static const int NumberOfAtomsInBond      = 4;
+      static const int NumberOfParametersInBond = 3;
+
+      // bond name
+
+      static const std::string BondName;
+
+      // log file reference
+
+      FILE* _log;
+
+
+      // Brook bond parameters
+
+      BrookBondParameters* _brookBondParameters;
+
+      // interface
+
+      OpenMMBrookInterface& _openMMBrookInterface;
+
+      // System reference
+
+      System& _system;
+
+};
+
+} // namespace OpenMM
+
+#endif /* OPENMM_BROOK_CALC_PROPER_DIHEDRAL_FORCE_KERNEL_H_ */

platforms/brook/src/BrookCalcRBTorsionForceKernel.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "OpenMMException.h"
+#include <sstream>
+#include "BrookCalcRBTorsionForceKernel.h"
+
+using namespace OpenMM;
+using namespace std;
+
+const std::string BrookCalcRBTorsionForceKernel::BondName = "RbDihedral";
+
+/** 
+ * BrookCalcRBTorsionForceKernel constructor
+ * 
+ * @param name                      kernel name
+ * @param platform                  platform
+ * @param OpenMMBrookInterface      OpenMM-Brook interface
+ * @param System                    System reference
+ *
+ */
+
+BrookCalcRBTorsionForceKernel::BrookCalcRBTorsionForceKernel( std::string name, const Platform& platform,
+                                                                      OpenMMBrookInterface& openMMBrookInterface, System& system ) :
+                     CalcRBTorsionForceKernel( name, platform ), _openMMBrookInterface( openMMBrookInterface ), _system( system ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookCalcRBTorsionForceKernel::BrookCalcRBTorsionForceKernel";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   _brookBondParameters              = NULL;
+   _log                              = NULL;
+
+   const BrookPlatform brookPlatform = dynamic_cast<const BrookPlatform&> (platform);
+   if( brookPlatform.getLog() != NULL ){
+      setLog( brookPlatform.getLog() );
+   }
+      
+}   
+
+/** 
+ * BrookCalcRBTorsionForceKernel destructor
+ * 
+ */
+
+BrookCalcRBTorsionForceKernel::~BrookCalcRBTorsionForceKernel( ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "BrookCalcRBTorsionForceKernel::BrookCalcRBTorsionForceKernel";
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+   delete _brookBondParameters;
+}
+
+/** 
+ * Get log file reference
+ * 
+ * @return  log file reference
+ *
+ */
+
+FILE* BrookCalcRBTorsionForceKernel::getLog( void ) const {
+   return _log;
+}
+
+/** 
+ * Set log file reference
+ * 
+ * @param  log file reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int BrookCalcRBTorsionForceKernel::setLog( FILE* log ){
+   _log = log;
+   return BrookCommon::DefaultReturnValue;
+}
+
+/** 
+ * Initialize the kernel, setting up the values of all the force field parameters.
+ * 
+ * @param system                    System reference
+ * @param force                     RbDihedralForce reference
+ *
+ */
+
+void BrookCalcRBTorsionForceKernel::initialize( const System& system, const RBTorsionForce& force ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "BrookCalcRBTorsionForceKernel::initialize";
+
+// ---------------------------------------------------------------------------------------
+
+   FILE* log                 = getLog();
+
+   // ---------------------------------------------------------------------------------------
+
+   // create _brookBondParameters object containing atom indices/parameters
+
+   int numberOfBonds         = force.getNumTorsions();
+
+   if( _brookBondParameters ){
+      delete _brookBondParameters;
+   }
+   _brookBondParameters = new BrookBondParameters( BondName, NumberOfAtomsInBond, NumberOfParametersInBond, numberOfBonds, getLog() );
+
+   for( int ii = 0; ii < numberOfBonds; ii++ ){
+
+      int particle1, particle2, particle3, particle4;
+      double c0, c1, c2, c3, c4, c5;
+
+      int particles[NumberOfAtomsInBond];
+      double parameters[NumberOfParametersInBond];
+
+      force.getTorsionParameters( ii, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5 ); 
+      particles[0]    = particle1;
+      particles[1]    = particle2;
+      particles[2]    = particle3;
+      particles[3]    = particle4;
+ 
+      parameters[0]   = c0;
+      parameters[1]   = c1;
+      parameters[2]   = c2;
+      parameters[3]   = c3;
+      parameters[4]   = c4;
+      parameters[5]   = c5;
+
+      _brookBondParameters->setBond( ii, particles, parameters );
+   }   
+   _openMMBrookInterface.setRBTorsionForceParameters( _brookBondParameters );
+   _openMMBrookInterface.setTriggerForceKernel( this );
+   _openMMBrookInterface.setTriggerEnergyKernel( this );
+
+   if( log ){
+      std::string contents = _brookBondParameters->getContentsString( ); 
+      (void) fprintf( log, "%s brookGbsa::contents\n%s", methodName.c_str(), contents.c_str() );
+      (void) fflush( log );
+   }
+
+   // ---------------------------------------------------------------------------------------
+    
+}
+
+/** 
+ * Compute forces given atom coordinates
+ * 
+ * @param context OpenMMContextImpl context
+ *
+ */
+
+void BrookCalcRBTorsionForceKernel::executeForces( OpenMMContextImpl& context ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName   = "BrookCalcRBTorsionForceKernel::executeForces";
+
+// ---------------------------------------------------------------------------------------
+
+   if( _openMMBrookInterface.getTriggerForceKernel() == this ){
+      _openMMBrookInterface.computeForces( context );
+   }
+
+   return;
+
+   // ---------------------------------------------------------------------------------------
+}
+
+/**
+ * Execute the kernel to calculate the energy
+ * 
+ * @param context OpenMMContextImpl context
+ *
+ * @return  potential energy
+ *
+ */
+
+double BrookCalcRBTorsionForceKernel::executeEnergy( OpenMMContextImpl& context ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "BrookCalcRBTorsionForceKernel::executeEnergy";
+
+// ---------------------------------------------------------------------------------------
+
+   if( _openMMBrookInterface.getTriggerEnergyKernel() == this ){
+      return (double) _openMMBrookInterface.computeEnergy( context );
+   } else {
+      return 0.0;
+   }
+
+}

platforms/brook/src/BrookCalcRBTorsionForceKernel.h

+#ifndef OPENMM_BROOK_CALC_RB_TORSION_FORCE_KERNEL_H_
+#define OPENMM_BROOK_CALC_RB_TORSION_FORCE_KERNEL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "kernels.h"
+#include "BrookPlatform.h"
+#include "BrookBondParameters.h"
+#include "OpenMMBrookInterface.h"
+
+namespace OpenMM {
+
+/**
+ * This kernel is invoked to calculate the RB torsion forces acting on the system.
+ */
+
+class BrookCalcRBTorsionForceKernel : public CalcRBTorsionForceKernel {
+
+   public:
+  
+      /**
+       * BrookCalcRBTorsionForceKernel constructor
+       */
+
+      BrookCalcRBTorsionForceKernel( std::string name, const Platform& platform, OpenMMBrookInterface& openMMBrookInterface, System& system );
+  
+      /**
+       * BrookCalcRBTorsionForceKernel destructor
+       */
+
+      ~BrookCalcRBTorsionForceKernel();
+  
+      /**
+       * Initialize the kernel, setting up the values to calculate harmonic bond force & energy
+       * 
+       * @param system                    System reference
+       * @param force                     RbDihedralForce reference
+       *
+       */
+
+      void initialize( const System& system, const RBTorsionForce& force );
+  
+      /**
+       * Execute the kernel to calculate the forces.
+       * 
+       * @param positions   atom coordiantes
+       * @param forces      output forces
+       *
+       */
+
+      void executeForces( OpenMMContextImpl& context );
+  
+      /**
+       * Execute the kernel to calculate the energy.
+       * 
+       * @param context    the context in which to execute this kernel
+       *
+       * @return  potential energy associated with the harmonic angle force
+       *
+       */
+
+      double executeEnergy( OpenMMContextImpl& context );
+
+      /** 
+       * Set log file reference
+       * 
+       * @param  log file reference
+       *
+       * @return DefaultReturnValue
+       *
+       */
+      
+      int setLog( FILE* log );
+
+      /* 
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string getContents( int level ) const;
+
+      /** 
+       * Get log file reference
+       * 
+       * @return  log file reference
+       *
+       */
+      
+      FILE* getLog( void ) const;
+      
+      /** 
+       * Get number of bonds
+       * 
+       * @return  number of bonds
+       *
+       */
+      
+      int getNumberOfBonds( void ) const;
+      
+      /** 
+       * Get indices/parameters
+       * 
+       * @return  BrookBondParameters containing atom indices/parameters
+       *
+       */
+      
+      BrookBondParameters* getBrookBondParameters( void ) const;
+      
+   private:
+   
+      static const int NumberOfAtomsInBond      = 4;
+      static const int NumberOfParametersInBond = 6;
+
+      // bond name
+
+      static const std::string BondName;
+
+      // log file reference
+
+      FILE* _log;
+
+
+      // Brook bond parameters
+
+      BrookBondParameters* _brookBondParameters;
+
+      // interface
+
+      OpenMMBrookInterface& _openMMBrookInterface;
+
+      // System reference
+
+      System& _system;
+
+};
+
+} // namespace OpenMM
+
+#endif /* OPENMM_BROOK_CALC_RB_TORSION_FORCE_KERNEL_H_ */

platforms/brook/src/OpenMMBrookInterface.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <limits>
+#include "OpenMMException.h"
+#include <sstream>
+
+#include "BrookStreamImpl.h"
+#include "OpenMMBrookInterface.h"
+#include "gpu/kforce.h"
+#include "gpu/kinvmap_gather.h"
+#include "NonbondedForce.h"
+
+using namespace OpenMM;
+using namespace std;
+
+/** 
+ * OpenMMBrookInterface constructor
+ * 
+ */
+
+OpenMMBrookInterface::OpenMMBrookInterface( void ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "OpenMMBrookInterface::OpenMMBrookInterface";
+
+// ---------------------------------------------------------------------------------------
+
+   _numberOfAtoms                           = 0;
+   _brookBonded                             = NULL;
+   _brookNonBonded                          = NULL;
+   _refForceField                           = NULL;
+   _log                                     = NULL;
+
+   _refForceField                           = NULL;
+   _refSystem                               = NULL;
+   _refOpenMMContext                        = NULL;
+   _referencePlatform                       = NULL;
+   _refVerletIntegrator                     = NULL;
+
+}   
+
+/** 
+ * OpenMMBrookInterface destructor
+ * 
+ */
+
+OpenMMBrookInterface::~OpenMMBrookInterface( ){
+
+// ---------------------------------------------------------------------------------------
+
+   // static const std::string methodName      = "OpenMMBrookInterface::OpenMMBrookInterface";
+
+// ---------------------------------------------------------------------------------------
+
+   delete _brookBonded;
+   delete _brookNonBonded;
+
+   // deleted w/ kernel delete? If activated, program crashes
+
+   //delete _refForceField;
+/*
+   delete _refSystem;
+   delete _refOpenMMContext;
+   delete _referencePlatform;
+   delete _refVerletIntegrator;
+*/
+}
+
+/** 
+ * Get log file reference
+ * 
+ * @return  log file reference
+ *
+ */
+
+FILE* OpenMMBrookInterface::getLog( void ) const {
+   return _log;
+}
+
+/** 
+ * Set log file reference
+ * 
+ * @param  log file reference
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int OpenMMBrookInterface::setLog( FILE* log ){
+   _log = log;
+   return BrookCommon::DefaultReturnValue;
+}
+
+/** 
+ * Set BrookBondParameters at specified index
+ * 
+ * @param index
+ * @param brookBondParameters brookBondParameters for BrookBondParameters
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int OpenMMBrookInterface::_setBondParameters( BondParameterIndices index, BrookBondParameters* brookBondParameters ){
+   _bondParameters[index] = brookBondParameters;
+   return BrookCommon::DefaultReturnValue;
+}
+
+/** 
+ * Set BrookBondParameters for harmonic bond force
+ * 
+ * @param brookBondParameters brookBondParameters for harmonic bond force
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int OpenMMBrookInterface::setHarmonicBondForceParameters( BrookBondParameters* brookBondParameters ){
+   return _setBondParameters( HarmonicAngleIndex, brookBondParameters );
+}
+
+/** 
+ * Set BrookBondParameters for harmonic angle force
+ * 
+ * @param brookBondParameters brookBondParameters for harmonic angle force
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int OpenMMBrookInterface::setHarmonicAngleForceParameters( BrookBondParameters* brookBondParameters ){
+   return _setBondParameters( HarmonicAngleIndex, brookBondParameters );
+}
+
+/** 
+ * Set BrookBondParameters for proper dihedral force
+ * 
+ * @param brookBondParameters brookBondParameters for proper dihedral force
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int OpenMMBrookInterface::setPeriodicTorsionForceParameters( BrookBondParameters* brookBondParameters ){
+   return _setBondParameters( PeriodicTorsionForceIndex, brookBondParameters );
+}
+
+/** 
+ * Set BrookBondParameters for RB dihedral force
+ * 
+ * @param brookBondParameters brookBondParameters for RB force
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int OpenMMBrookInterface::setRBTorsionForceParameters( BrookBondParameters* brookBondParameters ){
+   return _setBondParameters( RbTorsionForceIndex, brookBondParameters );
+}
+
+/** 
+ * Set BrookBondParameters for LJ 14 force
+ * 
+ * @param brookBondParameters brookBondParameters for LJ 14 force
+ *
+ * @return  DefaultReturnValue
+ *
+ */
+
+int OpenMMBrookInterface::setLJ14( BrookBondParameters* brookBondParameters ){
+   return _setBondParameters( LJ14Index, brookBondParameters );
+}
+
+/** 
+ * 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 lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
+ * @param coulomb14Scale            the factor by which Coulomb interactions should be reduced for bonded 1-4 pairs
+ * @param exclusions                the i'th element lists the indices of all atoms with which the i'th atom should not interact through
+ *                                  nonbonded forces.  Bonded 1-4 pairs are also included in this list, since they should be omitted from
+ *                                  the standard nonbonded calculation.
+ * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+ * @param nonbondedMethod           the method to use for handling long range nonbonded interactions
+ * @param nonbondedCutoff           the cutoff distance for nonbonded interactions (if nonbondedMethod involves a cutoff)
+ * @param periodicBoxSize           the size of the periodic box (if nonbondedMethod involves a periodic boundary conditions)
+ */
+ 
+void OpenMMBrookInterface::initialize( 
+        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,        double lj14Scale, double coulomb14Scale,
+        const vector<set<int> >& exclusions, const vector<vector<double> >& nonbondedParameters,
+         //NonbondedMethod nonbondedMethod,
+        double nonbondedCutoff, double periodicBoxSize[3] ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "OpenMMBrookInterface::initialize";
+
+// ---------------------------------------------------------------------------------------
+
+   FILE* log                 = getLog();
+/*
+   _numberOfAtoms            = nonbondedParameters.size();
+
+   // ---------------------------------------------------------------------------------------
+
+   // bonded
+
+   if( _brookBonded ){
+      delete _brookBonded;
+   }
+   _brookBonded              = new BrookBonded();
+   _brookBonded->setLog( log );
+    
+   _brookBonded->setup( _numberOfAtoms, 
+                        bondIndices,            bondParameters, 
+                        angleIndices,           angleParameters,
+                        periodicTorsionIndices, periodicTorsionParameters,
+                        rbTorsionIndices,       rbTorsionParameters,
+                        bonded14Indices,        nonbondedParameters,
+                        lj14Scale, coulomb14Scale, getPlatform() );
+
+   // echo contents
+
+   if( log ){
+      std::string contents = _brookBonded->getContentsString( ); 
+      (void) fprintf( log, "%s brookBonded::contents\n%s", methodName.c_str(), contents.c_str() );
+      (void) fflush( log );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // nonbonded
+
+   if( _brookNonBonded ){
+      delete _brookNonBonded;
+   }
+   _brookNonBonded           = new BrookNonBonded();
+   _brookNonBonded->setLog( log );
+
+   _brookNonBonded->setup( _numberOfAtoms, nonbondedParameters, exclusions, getPlatform() );
+
+   // echo contents
+
+   if( log ){
+      std::string contents = _brookNonBonded->getContentsString( ); 
+      (void) fprintf( log, "%s brookNonBonded::contents\n%s", methodName.c_str(), contents.c_str() );
+      (void) fflush( log );
+   }
+
+   // ---------------------------------------------------------------------------------------
+
+   // used for calculating energy
+
+
+   _referenceCalcNonbondedForceKernel = new ReferenceCalcNonbondedForceKernel( getName(), getPlatform() );
+   _referenceCalcNonbondedForceKernel->initialize( bondIndices,            bondParameters, 
+                                                         angleIndices,           angleParameters,
+                                                         periodicTorsionIndices, periodicTorsionParameters,
+                                                         rbTorsionIndices,       rbTorsionParameters,
+                                                         bonded14Indices,        lj14Scale, coulomb14Scale,
+                                                         exclusions, nonbondedParameters,
+                                                         nonbondedMethod, nonbondedCutoff, periodicBoxSize );
+
+   _refForceField = new NonbondedForce( _numberOfAtoms, bondIndices.size(), angleIndices.size(), 
+                                              periodicTorsionIndices.size(),  rbTorsionIndices.size() );
+
+   // bonds
+
+   vector<vector<int> >::const_iterator bondI_it    = bondIndices.begin();
+   vector<vector<double> >::const_iterator bondP_it = bondParameters.begin();
+   for( unsigned int ii = 0; ii < bondIndices.size(); ii++, bondI_it++, bondP_it++ ){
+      const vector<int>    bondInd = *bondI_it;
+      const vector<double> bondPrm = *bondP_it;
+      _refForceField->setBondParameters( ii, bondInd[0], bondInd[1], bondPrm[0], bondPrm[1] );
+
+printf( "%s Bond: [%d %d] [%.5e %.5e ]\n", methodName.c_str(),
+        bondInd[0], bondInd[1], bondPrm[0], bondPrm[1] ); fflush( stdout );
+   }
+
+   // angles
+
+   vector<vector<int> >::const_iterator angleI_it    = angleIndices.begin();
+   vector<vector<double> >::const_iterator angleP_it = angleParameters.begin();
+   for( unsigned int ii = 0; ii < angleIndices.size(); ii++, angleI_it++, angleP_it++ ){
+      const vector<int>    angleInd = *angleI_it;
+      const vector<double> anglePrm = *angleP_it;
+      _refForceField->setAngleParameters( ii, angleInd[0], angleInd[1], angleInd[2], anglePrm[0], anglePrm[1] );
+   }
+
+   // periodicTorsions
+
+   vector<vector<int> >::const_iterator periodicTorsionI_it    = periodicTorsionIndices.begin();
+   vector<vector<double> >::const_iterator periodicTorsionP_it = periodicTorsionParameters.begin();
+   for( unsigned int ii = 0; ii < periodicTorsionIndices.size(); ii++, periodicTorsionI_it++, periodicTorsionP_it++ ){
+      const vector<int>    periodicTorsionInd = *periodicTorsionI_it;
+      const vector<double> periodicTorsionPrm = *periodicTorsionP_it;
+      _refForceField->setPeriodicTorsionParameters( ii, periodicTorsionInd[0], periodicTorsionInd[1], periodicTorsionInd[2], periodicTorsionInd[3],
+                                                    (int) (periodicTorsionPrm[2] + 0.001), periodicTorsionPrm[1], periodicTorsionPrm[0] );
+
+printf( "PeriodicTor: [%d %d %d %d] [%.5e %.5e %.5e]\n", periodicTorsionInd[0], periodicTorsionInd[1], periodicTorsionInd[2], periodicTorsionInd[3],
+                                                         periodicTorsionPrm[2], periodicTorsionPrm[1], periodicTorsionPrm[0] ); fflush( stdout );
+   }
+
+   // rbTorsions
+
+   vector<vector<int> >::const_iterator rbTorsionI_it    = rbTorsionIndices.begin();
+   vector<vector<double> >::const_iterator rbTorsionP_it = rbTorsionParameters.begin();
+   for( unsigned int ii = 0; ii < rbTorsionIndices.size(); ii++, rbTorsionI_it++, rbTorsionP_it++ ){
+      const vector<int>    rbTorsionInd = *rbTorsionI_it;
+      const vector<double> rbTorsionPrm = *rbTorsionP_it;
+      _refForceField->setRBTorsionParameters( ii, rbTorsionInd[0], rbTorsionInd[1], rbTorsionInd[2], rbTorsionInd[3],
+                                              rbTorsionPrm[0], rbTorsionPrm[1], rbTorsionPrm[2], rbTorsionPrm[3], 
+                                              rbTorsionPrm[4], rbTorsionPrm[5] );
+*/
+
+/*
+printf( "RbTor: [%d %d %d %d] [%.5e %.5e %.5e %.5e %.5e %.5e]\n", rbTorsionInd[0], rbTorsionInd[1], rbTorsionInd[2], rbTorsionInd[3],
+                                                         rbTorsionPrm[0], rbTorsionPrm[1], rbTorsionPrm[2], rbTorsionPrm[3], rbTorsionPrm[4], rbTorsionPrm[5] ); fflush( stdout );
+   }
+
+*/
+   // nonbonded
+
+/*
+   for( unsigned int ii = 0; ii < nonbondedParameters.size(); ii++ ){
+
+      vector<double> nonbondedParameterVector 
+                                     = nonbondedParameters[ii];
+
+      double c6                      = nonbondedParameterVector[1];
+      double c12                     = nonbondedParameterVector[2];
+      double charge                  = nonbondedParameterVector[0];
+
+      // int index, double charge, double radius, double depth
+
+      _refForceField->setAtomParameters( ii, charge, c6, c12 );
+
+   }
+*/
+
+   // ---------------------------------------------------------------------------------------
+
+   return;
+}
+
+/** 
+ * Execute the kernel to calculate the bonded & nonbonded forces
+ * 
+ * @param positions   stream of type Double3 containing the position (x, y, z) of each atom
+ * @param forces      stream of type Double3 containing the force (x, y, z) on each atom.  On entry, this contains the forces that
+ *                    have been calculated so far.  The kernel should add its own forces to the values already in the stream.
+ */
+
+void OpenMMBrookInterface::executeForces( const Stream& positions, Stream& forces ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName      = "OpenMMBrookInterface::executeForces";
+
+   static const int I_Stream                = 0;
+   static const int J_Stream                = 1;
+   static const int K_Stream                = 2;
+   static const int L_Stream                = 3;
+
+   static const int PrintOn                 = 0;
+   static const int MaxErrorMessages        = 2;
+   static       int ErrorMessages           = 0;
+
+   static const float4 dummyParameters( 0.0, 0.0, 0.0, 0.0 );
+
+   // static const int debug                   = 1;
+
+// ---------------------------------------------------------------------------------------
+
+
+   const BrookStreamImpl& positionStreamC              = dynamic_cast<const BrookStreamImpl&> (positions.getImpl());
+   BrookStreamImpl& positionStream                     = const_cast<BrookStreamImpl&>         (positionStreamC);
+   BrookStreamImpl& forceStream                        = dynamic_cast<BrookStreamImpl&>       (forces.getImpl());
+   
+   // nonbonded forces
+
+   // added charge stream to knbforce_CDLJ4
+   // libs generated from ~/src/gmxgpu-nsqOpenMM
+
+   BrookFloatStreamInternal**  nonbondedForceStreams = _brookNonBonded->getForceStreams();
+
+   float epsfac                                      = 138.935485f;
+
+   knbforce_CDLJ4(
+             (float) _brookNonBonded->getNumberOfAtoms(),
+             (float) _brookNonBonded->getAtomSizeCeiling(),
+             (float) _brookNonBonded->getDuplicationFactor(),
+             (float) _brookNonBonded->getAtomStreamHeight( ),
+             (float) _brookNonBonded->getAtomStreamWidth( ),
+             (float) _brookNonBonded->getJStreamWidth( ),
+             (float) _brookNonBonded->getPartialForceStreamWidth( ),
+             epsfac,
+             dummyParameters,
+             positionStream.getBrookStream(),
+             _brookNonBonded->getChargeStream()->getBrookStream(),
+             _brookNonBonded->getOuterVdwStream()->getBrookStream(),
+             _brookNonBonded->getInnerSigmaStream()->getBrookStream(),
+             _brookNonBonded->getInnerEpsilonStream()->getBrookStream(),
+             _brookNonBonded->getExclusionStream()->getBrookStream(),
+             nonbondedForceStreams[0]->getBrookStream(),
+             nonbondedForceStreams[1]->getBrookStream(),
+             nonbondedForceStreams[2]->getBrookStream(),
+             nonbondedForceStreams[3]->getBrookStream()
+           );
+
+/*
+float zerof = 0.0f;
+nonbondedForceStreams[0]->fillWithValue( &zerof );
+nonbondedForceStreams[1]->fillWithValue( &zerof );
+nonbondedForceStreams[2]->fillWithValue( &zerof );
+nonbondedForceStreams[3]->fillWithValue( &zerof );
+*/
+
+   // diagnostics
+
+   if( 1 && PrintOn ){
+      (void) fprintf( getLog(), "\nPost knbforce_CDLJ4: atoms=%6d ceiling=%3d dupFac=%3d", _brookNonBonded->getNumberOfAtoms(),  
+                                                                                           _brookNonBonded->getAtomSizeCeiling(),
+                                                                                           _brookNonBonded->getDuplicationFactor()  );
+
+      (void) fprintf( getLog(), "\n                      hght=%6d   width=%3d   jWid=%3d", _brookNonBonded->getAtomStreamHeight( ),
+                                                                                           _brookNonBonded->getAtomStreamWidth( ),
+                                                                                           _brookNonBonded->getJStreamWidth( ) );
+      (void) fprintf( getLog(), "\n                      pFrc=%6d     eps=%12.5e\n",       _brookNonBonded->getPartialForceStreamWidth( ), epsfac );
+
+      (void) fprintf( getLog(), "\nOuterVdwStreamd\n" );
+      _brookNonBonded->getOuterVdwStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nInnerSigmaStream\n" );
+      _brookNonBonded->getInnerSigmaStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nInnerEpsilonStream\n" );
+      _brookNonBonded->getInnerEpsilonStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nExclusionStream\n" );
+      _brookNonBonded->getExclusionStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nChargeStream\n" );
+      _brookNonBonded->getChargeStream()->printToFile( getLog() );
+
+      for( int ii = 0; ii < 4; ii++ ){
+         (void) fprintf( getLog(), "\nForce stream %d\n", ii );
+         nonbondedForceStreams[ii]->printToFile( getLog() );
+      }
+   }
+
+// ---------------------------------------------------------------------------------------
+
+   // gather forces
+
+   kMergeFloat3_4_nobranch( (float) _brookNonBonded->getDuplicationFactor(),
+                            (float) _brookNonBonded->getAtomStreamWidth(),
+                            (float) _brookNonBonded->getPartialForceStreamWidth(),
+                            (float) _brookNonBonded->getNumberOfAtoms(),
+                            (float) _brookNonBonded->getAtomSizeCeiling(),
+                            (float) _brookNonBonded->getOuterLoopUnroll(),
+                            nonbondedForceStreams[0]->getBrookStream(),
+                            nonbondedForceStreams[1]->getBrookStream(),
+                            nonbondedForceStreams[2]->getBrookStream(),
+                            nonbondedForceStreams[3]->getBrookStream(),
+                            forceStream.getBrookStream() );
+
+   // diagnostics
+
+   if( 0 && PrintOn ){
+
+      (void) fprintf( getLog(), "\nNB forces" );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      brookStreamInternalF->printToFile( getLog() );
+
+   }
+
+// ---------------------------------------------------------------------------------------
+
+   // bonded
+
+         epsfac                                        = (float) (_brookBonded->getLJ_14Scale()*_brookBonded->getCoulombFactor());
+   float width                                         = (float) (_brookBonded->getInverseMapStreamWidth());
+
+   // bonded forces
+
+   BrookFloatStreamInternal**  bondedParameters        = _brookBonded->getBondedParameterStreams();
+   BrookFloatStreamInternal**  bondedForceStreams      = _brookBonded->getBondedForceStreams();
+
+   BrookFloatStreamInternal**  inverseStreamMaps[4];
+   inverseStreamMaps[0]                                = _brookBonded->getInverseStreamMapsStreams( 0 );
+   inverseStreamMaps[1]                                = _brookBonded->getInverseStreamMapsStreams( 1 );
+   inverseStreamMaps[2]                                = _brookBonded->getInverseStreamMapsStreams( 2 );
+   inverseStreamMaps[3]                                = _brookBonded->getInverseStreamMapsStreams( 3 );
+
+   kbonded_CDLJ( epsfac, 
+                 (float) bondedForceStreams[0]->getStreamWidth(),
+                 dummyParameters,
+                 positionStream.getBrookStream(),
+                 _brookBonded->getChargeStream()->getBrookStream(),
+                 _brookBonded->getAtomIndicesStream()->getBrookStream(),
+                 bondedParameters[0]->getBrookStream(),
+                 bondedParameters[1]->getBrookStream(),
+                 bondedParameters[2]->getBrookStream(),
+                 bondedParameters[3]->getBrookStream(),
+                 bondedParameters[4]->getBrookStream(),
+                 bondedForceStreams[0]->getBrookStream(),
+                 bondedForceStreams[1]->getBrookStream(),
+                 bondedForceStreams[2]->getBrookStream(),
+                 bondedForceStreams[3]->getBrookStream() );
+
+
+   // diagnostics
+
+   if( 1 && PrintOn ){
+
+      int countPrintInvMap[4] = { 3, 5, 2, 4 }; 
+
+      (void) fprintf( getLog(), "\nPost kbonded_CDLJ: epsFac=%.6f %.6f %.6f", epsfac, _brookBonded->getLJ_14Scale(), _brookBonded->getCoulombFactor());
+      (void) fprintf( getLog(), "\nAtom indices stream\n" );
+      _brookBonded->getAtomIndicesStream()->printToFile( getLog() );
+
+      (void) fprintf( getLog(), "\nCharge stream\n" );
+      _brookBonded->getChargeStream()->printToFile( getLog() );
+
+      for( int ii = 0; ii < 5; ii++ ){
+         (void) fprintf( getLog(), "\nParam stream %d\n", ii );
+         bondedParameters[ii]->printToFile( getLog() );
+      }
+      for( int ii = 0; ii < 4; ii++ ){
+         (void) fprintf( getLog(), "\nForce stream %d\n", ii );
+         bondedForceStreams[ii]->printToFile( getLog() );
+      }
+
+/*
+      (void) fprintf( getLog(), "\nNB1 forces" );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      brookStreamInternalF->printToFile( getLog() );
+*/
+
+      (void) fprintf( getLog(), "\nInverse map streams -- K_Stream cnt=%d\n", _brookBonded->getInverseMapStreamCount( K_Stream ) );
+      for( int ii = 0; ii < 4; ii++ ){
+         for( int jj = 0; jj < countPrintInvMap[ii]; jj++ ){
+            (void) fprintf( getLog(), "\n   Inverse map streams index=%d %d\n", ii, jj );
+            inverseStreamMaps[ii][jj]->printToFile( getLog() );
+         }
+      }
+   }
+
+   // gather forces
+
+   if( _brookBonded->getInverseMapStreamCount( I_Stream ) == 3 && _brookBonded->getInverseMapStreamCount( K_Stream ) == 3 ){
+
+      kinvmap_gather3_3( width,
+
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][2]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][2]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+
+   } else if( _brookBonded->getInverseMapStreamCount( I_Stream ) == 3 && _brookBonded->getInverseMapStreamCount( K_Stream ) == 4 ){
+
+      kinvmap_gather3_4( width,
+
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][2]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][2]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][3]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+
+   } else if( _brookBonded->getInverseMapStreamCount( I_Stream ) == 3 && _brookBonded->getInverseMapStreamCount( K_Stream ) == 5 ){
+
+      kinvmap_gather3_5( width,
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][2]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][2]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][3]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][4]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+
+   } else {
+
+      // case not handled -- throw an exception
+
+      if( _brookBonded->getLog() && ErrorMessages++ < MaxErrorMessages ){
+         (void) fprintf( _brookBonded->getLog(), "%s case: I-map=%d K-map=%d -- not handled.\n",
+                          methodName.c_str(), _brookBonded->getInverseMapStreamCount( I_Stream ),
+                                              _brookBonded->getInverseMapStreamCount( K_Stream ) );
+         (void) fflush(  _brookBonded->getLog() );
+      }
+
+      kinvmap_gather3_3( width,
+
+                         inverseStreamMaps[I_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[I_Stream][2]->getBrookStream(),
+                         bondedForceStreams[I_Stream]->getBrookStream(),
+
+                         inverseStreamMaps[K_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[K_Stream][2]->getBrookStream(),
+                         bondedForceStreams[K_Stream]->getBrookStream(),
+
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+/*
+      std::stringstream message;
+      message << methodName << "I-maps=" << _brookBonded->getInverseMapStreamCount( I_Stream ) << " and " << 
+                               "K-maps=" << _brookBonded->getInverseMapStreamCount( K_Stream ) << " not handled.";
+      throw OpenMMException( message.str() );
+*/
+
+   }
+
+   // diagnostics
+
+   if( 0 && PrintOn ){
+
+      (void) fprintf( getLog(), "\nPost 3_4/3_5 && NB forces" );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      brookStreamInternalF->printToFile( getLog() );
+
+   }
+
+   if( _brookBonded->getInverseMapStreamCount( J_Stream ) == 5 && _brookBonded->getInverseMapStreamCount( L_Stream ) == 2 ){
+      kinvmap_gather5_2( width,
+                         inverseStreamMaps[J_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][2]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][3]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][4]->getBrookStream(),
+                         bondedForceStreams[J_Stream]->getBrookStream(),
+                         inverseStreamMaps[L_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[L_Stream][1]->getBrookStream(),
+                         bondedForceStreams[L_Stream]->getBrookStream(),
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+   
+   } else {
+
+      // case not handled -- throw an exception
+
+      if( _brookBonded->getLog() && ErrorMessages++ < MaxErrorMessages ){
+         (void) fprintf( _brookBonded->getLog(), "%s case: J-map=%d L-map=%d -- not handled.\n",
+                          methodName.c_str(), _brookBonded->getInverseMapStreamCount( J_Stream ),
+                                              _brookBonded->getInverseMapStreamCount( L_Stream ) );
+         (void) fflush(  _brookBonded->getLog() );
+      }
+
+      kinvmap_gather5_2( width,
+                         inverseStreamMaps[J_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][1]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][2]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][3]->getBrookStream(),
+                         inverseStreamMaps[J_Stream][4]->getBrookStream(),
+                         bondedForceStreams[J_Stream]->getBrookStream(),
+                         inverseStreamMaps[L_Stream][0]->getBrookStream(),
+                         inverseStreamMaps[L_Stream][1]->getBrookStream(),
+                         bondedForceStreams[L_Stream]->getBrookStream(),
+                         forceStream.getBrookStream(), forceStream.getBrookStream() );
+/*
+      std::stringstream message;
+      message << methodName << "J-maps=" << _brookBonded->getInverseMapStreamCount( J_Stream ) << " and " << 
+                               "L-maps=" << _brookBonded->getInverseMapStreamCount( L_Stream ) << " not handled.";
+      throw OpenMMException( message.str() );
+*/
+
+   }
+
+   // diagnostics
+
+   if( 1 && PrintOn ){
+
+      (void) fprintf( getLog(), "\nFinal NB & bonded forces" );
+      BrookStreamInternal* brookStreamInternalF   = forceStream.getBrookStreamImpl();
+      brookStreamInternalF->printToFile( getLog() );
+/*
+      void* dataV = brookStreamInternalF->getData(1);
+      float* data = (float*) dataV;
+      (void) fprintf( getLog(), "\nFinal NB & bonded forces RAW\n" );
+      for( int ii = 0; ii < _brookNonBonded->getNumberOfAtoms()*3; ii += 3 ){
+         (void) fprintf( getLog(), "%d [%.6e %.6e %.6e]\n", ii, data[ii], data[ii+1], data[ii+2] );
+      }
+*/
+
+   }
+
+   // ---------------------------------------------------------------------------------------
+}
+
+/**
+ * Execute the kernel to calculate the energy.
+ * 
+ * @param positions   atom positions
+ *
+ * @return  potential energy due to the NonbondedForce
+ * Currently always return 0.0 since energies not calculated on gpu
+ *
+ */
+
+double OpenMMBrookInterface::executeEnergy( const Stream& atomPositions ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "OpenMMBrookInterface::executeEnergy";
+
+// ---------------------------------------------------------------------------------------
+
+   const BrookStreamImpl& positionStreamC              = dynamic_cast<const BrookStreamImpl&> (atomPositions.getImpl());
+   BrookStreamImpl& positionStream                     = const_cast<BrookStreamImpl&>         (positionStreamC);
+   BrookOpenMMFloat* positionsF                        = (BrookOpenMMFloat*) positionStream.getData();
+
+   OpenMMContext* context                              = getReferenceOpenMMContext( atomPositions.getSize() );
+
+   vector<Vec3> positions( positionStream.getSize() );
+   int index = 0;
+   for( int ii = 0; ii < positionStream.getSize(); ii++, index += 3 ){
+      positions[ii] = Vec3( positionsF[index], positionsF[index+1], positionsF[index+2] );
+   }
+   context->setPositions(positions);
+
+   State state     = context->getState( State::Energy );
+   double energy   = state.getPotentialEnergy();
+
+(void) fprintf( stdout, "BrookCalcStandardMMForceFieldKernel::executeEnergy E=%.5e\n", energy );
+for( int ii = 0; ii < positionStream.getSize(); ii++, index += 3 ){
+   (void) fprintf( stdout, "   %5d Pos[%.5f %.5f %.5f]\n", ii, positions[ii][0], positions[ii][1], positions[ii][2] );
+}
+fflush( stdout );
+
+   return energy;
+
+}
+
+/**
+ * Get reference Context
+ * 
+ * @param numberOfAtoms  number of atoms
+ *
+ * @return  OpenMMContext
+ *
+ */
+
+OpenMMContext* OpenMMBrookInterface::getReferenceOpenMMContext( int numberOfAtoms ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "OpenMMBrookInterface::getReferenceOpenMMContext";
+
+// ---------------------------------------------------------------------------------------
+
+   if( _refOpenMMContext == NULL ){
+      
+      _referencePlatform       = new ReferencePlatform();
+      _refSystem               = new System( numberOfAtoms, 0 ); 
+		_refVerletIntegrator     = new VerletIntegrator( 0.01 );
+
+      _refSystem->addForce( _refForceField );
+
+      _refOpenMMContext        = new OpenMMContext( *_refSystem, *_refVerletIntegrator, *_referencePlatform );
+   }
+
+   return _refOpenMMContext;
+
+}
+
+/**
+ * Execute the kernel to calculate the energy.
+ * 
+ * @param positions   atom positions
+ *
+ * @return  potential energy due to the NonbondedForce
+ * Currently always return 0.0 since energies not calculated on gpu
+ *
+ */
+
+double OpenMMBrookInterface::executeEnergyOld( const Stream& atomPositions ){
+
+// ---------------------------------------------------------------------------------------
+
+   //static const std::string methodName      = "OpenMMBrookInterface::executeEnergy";
+
+// ---------------------------------------------------------------------------------------
+
+   const BrookStreamImpl& positionStreamC              = dynamic_cast<const BrookStreamImpl&> (atomPositions.getImpl());
+   BrookStreamImpl& positionStream                     = const_cast<BrookStreamImpl&>         (positionStreamC);
+   BrookOpenMMFloat* positionsF                        = (BrookOpenMMFloat*) positionStream.getData();
+
+   ReferencePlatform refPlatform;
+   System system( atomPositions.getSize(), 0); 
+   VerletIntegrator integrator( 0.01 );
+
+   system.addForce( _refForceField );
+
+   OpenMMContext context( system, integrator, refPlatform );
+
+   vector<Vec3> positions( positionStream.getSize() );
+   int index = 0;
+   for( int ii = 0; ii < positionStream.getSize(); ii++, index += 3 ){
+      positions[ii] = Vec3( positionsF[index], positionsF[index+1], positionsF[index+2] );
+   }
+   context.setPositions(positions);
+
+   State state     = context.getState( State::Energy );
+   double energy   = state.getPotentialEnergy();
+
+   //delete forceField;
+
+// (void) fprintf( stdout, "OpenMMBrookInterface::executeEnergy E=%.5e\n", energy ); fflush( stdout );
+
+   return energy;
+
+}

platforms/brook/src/OpenMMBrookInterface.h

+#ifndef OPENMM_BROOK_INTERFACE_H_
+#define OPENMM_BROOK_INTERFACE_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Authors: Peter Eastman, Mark Friedrichs                                    *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "kernels.h"
+#include "../../reference/src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "BrookBonded.h"
+#include "BrookNonBonded.h"
+#include "NonbondedForce.h"
+#include "OpenMMContext.h"
+#include "System.h"
+#include "ReferencePlatform.h"
+#include "VerletIntegrator.h"
+
+class BrookBondParameters;
+
+namespace OpenMM {
+
+/**
+ * OpenMM-Brook interface methods
+ */
+class OpenMMBrookInterface {
+
+   public:
+  
+      OpenMMBrookInterface( void );
+  
+      ~OpenMMBrookInterface();
+  
+      /**
+       * 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 lj14Scale                 the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
+       * @param coulomb14Scale            the factor by which Coulomb interactions should be reduced for bonded 1-4 pairs
+       * @param exclusions                the i'th element lists the indices of all atoms with which the i'th atom should not interact through
+       *                                  nonbonded forces.  Bonded 1-4 pairs are also included in this list, since they should be omitted from
+       *                                  the standard nonbonded calculation.
+       * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
+       * @param nonbondedMethod           the method to use for handling long range nonbonded interactions
+       * @param nonbondedCutoff           the cutoff distance for nonbonded interactions (if nonbondedMethod involves a cutoff)
+       * @param periodicBoxSize           the size of the periodic box (if nonbondedMethod involves a periodic boundary conditions)
+       *
+       */
+
+
+      void initialize( const std::vector<std::vector<int> >& bondIndices,      const std::vector<std::vector<double> >& bondParameters,
+                       const std::vector<std::vector<int> >& angleIndices,     const std::vector<std::vector<double> >& angleParameters,
+                       const std::vector<std::vector<int> >& periodicTorsionIndices, const std::vector<std::vector<double> >& periodicTorsionParameters,
+                       const std::vector<std::vector<int> >& rbTorsionIndices, const std::vector<std::vector<double> >& rbTorsionParameters,
+                       const std::vector<std::vector<int> >& bonded14Indices,  double lj14Scale, double coulomb14Scale,
+                       const std::vector<std::set   <int> >& exclusions,       const std::vector<std::vector<double> >& nonbondedParameters,
+                       //NonbondedMethod nonbondedMethod,
+                       double nonbondedCutoff, double periodicBoxSize[3] );
+
+  
+      /**
+       * Execute the kernel to calculate the forces.
+       * 
+       * @param positions   a Stream of type Double3 containing the position (x, y, z) of each atom
+       * @param forces      a Stream of type Double3 containing the force (x, y, z) on each atom.  On entry, this contains the forces that
+       *                    have been calculated so far.  The kernel should add its own forces to the values already in the stream.
+       */
+
+      void executeForces( const Stream& positions, Stream& forces );
+  
+      /**
+       * Execute the kernel to calculate the energy.
+       * 
+       * @param positions   a Stream of type Double3 containing the position (x, y, z) of each atom
+       *
+       * @return the potential energy due to the NonbondedForce
+       *
+       * Currently always return 0.0 since energies not calculated on gpu
+       */
+
+      double executeEnergy( const Stream& positions );
+      double executeEnergyOld( const Stream& positions );
+
+      /**
+       * Get reference Context
+       * 
+       * @param numberOfAtoms  number of atoms
+       *
+       * @return  OpenMMContext
+       *
+       */
+      
+      OpenMMContext* getReferenceOpenMMContext( int numberOfAtoms );
+      
+      
+      /** 
+       * Set log file reference
+       * 
+       * @param  log file reference
+       *
+       * @return DefaultReturnValue
+       *
+       */
+      
+      int setLog( FILE* log );
+
+      /* 
+       * Get contents of object
+       *
+       * @param level of dump
+       *
+       * @return string containing contents
+       *
+       * */
+      
+      std::string getContents( int level ) const;
+
+      /** 
+       * Get log file reference
+       * 
+       * @return  log file reference
+       *
+       */
+      
+      FILE* getLog( void ) const;
+      
+      /** 
+       * Compute forces
+       *
+       * @param context OpenMMContextImpl context 
+       *
+       */
+      
+      void computeForces( OpenMMContextImpl& context );
+      
+      /** 
+       * Compute energy
+       *
+       * @param context OpenMMContextImpl context 
+       *
+       */
+      
+      float computeEnergy( OpenMMContextImpl& context );
+      
+      /** 
+       * Set trigger Force Kernel
+       *
+       * @param triggerForceKernel kernel to calculate force
+       *
+       */
+      
+      void setTriggerForceKernel( KernelImpl* triggerForceKernel );
+      
+      /** 
+       * Set trigger Energy Kernel
+       *
+       * @param triggerEnergyKernel kernel to calculate energy
+       *
+       */
+      
+      void setTriggerEnergyKernel( KernelImpl* triggerForceKernel );
+      
+      /** 
+       * Get trigger Force Kernel
+       *
+       * @return triggerForceKernel kernel to calculate force
+       *
+       */
+      
+      KernelImpl* getTriggerForceKernel( void ) const;
+      
+      /** 
+       * Get trigger Energy Kernel
+       *
+       * @return triggerEnergyKernel kernel to calculate energy
+       *
+       */
+      
+      KernelImpl* getTriggerEnergyKernel( void ) const;
+      
+      /** 
+       * Set BrookBondParameters for harmonic angle force
+       * 
+       * @param brookBondParameters brookBondParameters for BrookBondParameters for harmonic angle force
+       *
+       * @return DefaultReturnValue
+       *
+       */
+      
+      int setHarmonicBondForceParameters( BrookBondParameters* brookBondParameters );
+      
+      /** 
+       * Set BrookBondParameters for harmonic angle force
+       * 
+       * @param brookBondParameters brookBondParameters for BrookBondParameters for harmonic angle force
+       *
+       * @return DefaultReturnValue
+       *
+       */
+      
+      int setHarmonicAngleForceParameters( BrookBondParameters* brookBondParameters );
+      
+      /** 
+       * Set BrookBondParameters for proper dihedral force
+       * 
+       * @param brookBondParameters brookBondParameters for proper dihedral force
+       *
+       * @return  DefaultReturnValue
+       *
+       */
+      
+      int setPeriodicTorsionForceParameters( BrookBondParameters* brookBondParameters );
+      
+      /** 
+       * Set BrookBondParameters for RB dihedral force
+       * 
+       * @param brookBondParameters brookBondParameters for RB force
+       *
+       * @return  DefaultReturnValue
+       *
+       */
+      
+      int setRBTorsionForceParameters( BrookBondParameters* brookBondParameters );
+      
+      /** 
+       * Set BrookBondParameters for LJ 14 force
+       * 
+       * @param brookBondParameters brookBondParameters for LJ 14 force
+       *
+       * @return  DefaultReturnValue
+       *
+       */
+      
+      int setLJ14( BrookBondParameters* brookBondParameters );
+      
+   private:
+   
+      enum BondParameterIndices { HarmonicBondIndex, HarmonicAngleIndex, PeriodicTorsionForceIndex, RbTorsionForceIndex, LJ14Index, LastBondForce };
+
+      // log file reference
+
+      FILE* _log;
+
+       // number of atoms
+
+       int _numberOfAtoms;
+   
+       // Brook bonded & nonbonded
+
+       BrookBonded* _brookBonded;
+       BrookNonBonded* _brookNonBonded;
+
+       BrookBondParameters* _bondParameters[LastBondForce];
+
+       // used to calculate energy
+
+       NonbondedForce*       _refForceField;
+       System*               _refSystem;
+       OpenMMContext*        _refOpenMMContext;
+       ReferencePlatform*    _referencePlatform;
+       VerletIntegrator*     _refVerletIntegrator;
+
+      /** 
+       * Set BrookBondParameters 
+       * 
+       * @param index
+       * @param brookBondParameters brookBondParameters for BrookBondParameters
+       *
+       * @return DefaultReturnValue
+       *
+       */
+      
+      int _setBondParameters( BondParameterIndices index, BrookBondParameters* brookBondParameters );
+      
+};
+
+} // namespace OpenMM
+
+#endif /* OPENMM_BROOK_INTERFACE_H_ */