Continuing to implement ReferencePlatform (kernels for StandardMMForceField are done).

CMakeLists.txt

@@ -16,7 +16,7 @@
 
 PROJECT (OpenMM)
 
-SUBDIRS (tests)
+SUBDIRS (tests platforms/reference/tests)
 
 ADD_DEFINITIONS(-DOPENMM_BUILDING_SHARED_LIBRARY)
 

olla/include/Kernel.h

@@ -64,6 +64,10 @@ public:
      * Get the name of this Kernel.
      */
     std::string getName() const;
+    /**
+     * Get the object which implements this Kernel.
+     */
+    const KernelImpl& getImpl() const;
     /**
      * Get the object which implements this Kernel.
      */

olla/include/Stream.h

@@ -113,6 +113,10 @@ public:
      * @param a pointer to the value.  It is assumed to be of the correct data type for this stream.
      */
     void fillWithValue(void* value);
+    /**
+     * Get the object which implements this Stream.
+     */
+    const StreamImpl& getImpl() const;
     /**
      * Get the object which implements this Stream.
      */

olla/include/kernels.h

@@ -41,14 +41,14 @@
 namespace OpenMM {
 
 /**
- * This kernel is invoked by StandardMMForceField to calculate the forces acting on the system.
+ * This kernel is invoked by StandardMMForceField to calculate the forces acting on the system and the energy of the system.
  */
-class CalcStandardMMForcesKernel : public KernelImpl {
+class CalcStandardMMForceFieldKernel : public KernelImpl {
 public:
     static std::string Name() {
-        return "CalcStandardMMForces";
+        return "CalcStandardMMForceField";
     }
-    CalcStandardMMForcesKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    CalcStandardMMForceFieldKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
     }
     /**
      * Initialize the kernel, setting up the values of all the force field parameters.
@@ -58,84 +58,50 @@ public:
      * @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 (periodicity, phase, k) for 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, van der Waals radius, van der Waals depth) for each atom
+     * @param nonbondedParameters       the nonbonded force parameters (charge, sigma, epsilon) for each atom
      */
     virtual 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, const std::vector<std::set<int> >& exclusions,
-            const std::vector<std::vector<double> >& nonbondedParameters) = 0;
+            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) = 0;
     /**
-     * Execute the kernel.
+     * 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.
      */
-    virtual void execute(const Stream& positions, Stream& forces) = 0;
-};
-
-/**
- * This kernel is invoked by StandardMMForceField to calculate the energy of the system.
- */
-class CalcStandardMMEnergyKernel : public KernelImpl {
-public:
-    static std::string Name() {
-        return "CalcStandardMMEnergy";
-    }
-    CalcStandardMMEnergyKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
-    }
-    /**
-     * 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 (periodicity, phase, k) 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 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, van der Waals radius, van der Waals depth) for each atom
-     */
-    virtual 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, const std::vector<std::set<int> >& exclusions,
-            const std::vector<std::vector<double> >& nonbondedParameters) = 0;
+    virtual void executeForces(const Stream& positions, Stream& forces) = 0;
     /**
-     * Execute the kernel.
+     * 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 StandardMMForceField
      */
-    virtual double execute(const Stream& positions) = 0;
+    virtual double executeEnergy(const Stream& positions) = 0;
 };
 
 /**
- * This kernel is invoked by GBSAOBCForceField to calculate the forces acting on the system.
+ * This kernel is invoked by GBSAOBCForceField to calculate the forces acting on the system and the energy of the system.
  */
-class CalcGBSAOBCForcesKernel : public KernelImpl {
+class CalcGBSAOBCForceFieldKernel : public KernelImpl {
 public:
     static std::string Name() {
         return "CalcGBSAOBCForces";
     }
-    CalcGBSAOBCForcesKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    CalcGBSAOBCForceFieldKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
     }
     /**
      * Initialize the kernel, setting up the values of all the force field parameters.
@@ -148,42 +114,20 @@ public:
     virtual void initialize(const std::vector<double>& bornRadii, const std::vector<std::vector<double> >& atomParameters,
             double solventDielectric, double soluteDielectric) = 0;
     /**
-     * Execute the kernel.
+     * 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.
      */
-    virtual void execute(const Stream& positions, Stream& forces) = 0;
-};
-
-/**
- * This kernel is invoked by GBSAOBCForceField to calculate the energy of the system.
- */
-class CalcGBSAOBCEnergyKernel : public KernelImpl {
-public:
-    static std::string Name() {
-        return "CalcGBSAOBCEnergy";
-    }
-    CalcGBSAOBCEnergyKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
-    }
-    /**
-     * Initialize the kernel, setting up the values of all the force field parameters.
-     * 
-     * @param bornRadii           the initial value of the Born radius for each atom
-     * @param atomParameters      the force parameters (charge, atomic radius, scaling factor) for each atom
-     * @param solventDielectric   the dielectric constant of the solvent
-     * @param soluteDielectric    the dielectric constant of the solute
-     */
-    virtual void initialize(const std::vector<double>& bornRadii, const std::vector<std::vector<double> >& atomParameters,
-            double solventDielectric, double soluteDielectric) = 0;
+    virtual void executeForces(const Stream& positions, Stream& forces) = 0;
     /**
-     * Execute the kernel.
+     * 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 GBSAOBCForceField
      */
-    virtual double execute(const Stream& positions) = 0;
+    virtual double executeEnergy(const Stream& positions) = 0;
 };
 
 /**

olla/src/Kernel.cpp

@@ -62,6 +62,10 @@ string Kernel::Kernel::getName() const {
     return impl->getName();
 }
 
+const KernelImpl& Kernel::getImpl() const {
+    return *impl;
+}
+
 KernelImpl& Kernel::getImpl() {
     return *impl;
 }

olla/src/Stream.cpp

@@ -83,6 +83,10 @@ void Stream::fillWithValue(void* value) {
     impl->fillWithValue(value);
 }
 
+const StreamImpl& Stream::getImpl() const {
+    return *impl;
+}
+
 StreamImpl& Stream::getImpl() {
     return *impl;
 }

openmmapi/include/OpenMMContext.h

@@ -100,7 +100,7 @@ public:
      * @param types the set of data types which should be stored in the State object.  This
      * should be a union of DataType values, e.g. (State::Positions | State::Velocities).
      */
-    State getState(State::DataType types) const;
+    State getState(int types) const;
     /**
      * Get the current time of the simulation (in picoseconds).
      */

openmmapi/include/StandardMMForceField.h

@@ -104,8 +104,8 @@ public:
      * 
      * @param index     the index of the atom for which to set parameters
      * @param charge    the charge of the atom, measured in units of the proton charge
-     * @param radius    the van der Waals radius of the atom, measured in angstroms
-     * @param depth     the well depth of the van der Waals interaction, measured in kcal/mol
+     * @param radius    the van der Waals radius of the atom (sigma in the Lennard Jones potential), measured in angstroms
+     * @param depth     the well depth of the van der Waals interaction (epsilon in the Lennard Jones potential), measured in kcal/mol
      */
     void setAtomParameters(int index, double charge, double radius, double depth);
     /**

openmmapi/include/Vec3.h

@@ -33,6 +33,7 @@
  * -------------------------------------------------------------------------- */
 
 #include <cassert>
+#include <iosfwd>
 
 namespace OpenMM {
 
@@ -69,6 +70,12 @@ private:
     double data[3];
 };
 
+template <class CHAR, class TRAITS>
+std::basic_ostream<CHAR,TRAITS>& operator<<(std::basic_ostream<CHAR,TRAITS>& o, const Vec3& v) {
+    o<<'['<<v[0]<<", "<<v[1]<<", "<<v[2]<<']';
+    return o;
+}
+
 } // namespace OpenMM
 
 #endif /*OPENMM_VEC3_H_*/

openmmapi/include/internal/GBSAOBCForceFieldImpl.h

@@ -61,7 +61,7 @@ public:
 private:
     void initialize(OpenMMContextImpl& context);
     GBSAOBCForceField& owner;
-    Kernel forceKernel, energyKernel;
+    Kernel kernel;
     bool hasInitialized;
 };
 

openmmapi/include/internal/StandardMMForceFieldImpl.h

@@ -65,7 +65,7 @@ private:
     void findExclusions(const std::vector<std::vector<int> >& bondIndices, std::vector<std::set<int> >& exclusions, std::set<std::pair<int, int> >& bonded14Indices) const;
     void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseAtom, int fromAtom, int currentLevel) const;
     StandardMMForceField& owner;
-    Kernel forceKernel, energyKernel;
+    Kernel kernel;
 };
 
 } // namespace OpenMM

openmmapi/src/GBSAOBCForceFieldImpl.cpp

@@ -42,8 +42,7 @@ GBSAOBCForceFieldImpl::GBSAOBCForceFieldImpl(GBSAOBCForceField& owner, OpenMMCon
 
 void GBSAOBCForceFieldImpl::initialize(OpenMMContextImpl& context) {
     hasInitialized = true;
-    forceKernel = context.getPlatform().createKernel(CalcStandardMMForcesKernel::Name());
-    energyKernel = context.getPlatform().createKernel(CalcStandardMMEnergyKernel::Name());
+    kernel = context.getPlatform().createKernel(CalcGBSAOBCForceFieldKernel::Name());
     std::vector<double> bornRadii;
     // TODO calculate the initial Born radii.
     vector<vector<double> > atomParameters;
@@ -54,27 +53,24 @@ void GBSAOBCForceFieldImpl::initialize(OpenMMContextImpl& context) {
         atomParameters[i].push_back(radius);
         atomParameters[i].push_back(scalingFactor);
     }
-    dynamic_cast<CalcGBSAOBCForcesKernel&>(forceKernel.getImpl()).initialize(bornRadii, atomParameters,
-            owner.getSolventDielectric(), owner.getSoluteDielectric());
-    dynamic_cast<CalcGBSAOBCEnergyKernel&>(forceKernel.getImpl()).initialize(bornRadii, atomParameters,
+    dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).initialize(bornRadii, atomParameters,
             owner.getSolventDielectric(), owner.getSoluteDielectric());
 }
 
 void GBSAOBCForceFieldImpl::calcForces(OpenMMContextImpl& context, Stream& forces) {
     if (!hasInitialized)
         initialize(context);
-    dynamic_cast<CalcGBSAOBCForcesKernel&>(forceKernel.getImpl()).execute(context.getPositions(), forces);
+    dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).executeForces(context.getPositions(), forces);
 }
 
 double GBSAOBCForceFieldImpl::calcEnergy(OpenMMContextImpl& context) {
     if (!hasInitialized)
         initialize(context);
-    return dynamic_cast<CalcGBSAOBCEnergyKernel&>(energyKernel.getImpl()).execute(context.getPositions());
+    return dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).executeEnergy(context.getPositions());
 }
 
 std::vector<std::string> GBSAOBCForceFieldImpl::getKernelNames() {
     std::vector<std::string> names;
-    names.push_back(CalcGBSAOBCForcesKernel::Name());
-    names.push_back(CalcGBSAOBCEnergyKernel::Name());
+    names.push_back(CalcGBSAOBCForceFieldKernel::Name());
     return names;
 }

openmmapi/src/OpenMMContext.cpp

@@ -61,12 +61,14 @@ Integrator& OpenMMContext::getIntegrator() {
 
 }
 
-State OpenMMContext::getState(State::DataType types) const {
-    State state(impl->getTime(), impl->getSystem().getNumAtoms(), types);
+State OpenMMContext::getState(int types) const {
+    State state(impl->getTime(), impl->getSystem().getNumAtoms(), State::DataType(types));
     if (types&State::Energy)
         state.setEnergy(impl->calcKineticEnergy(), impl->calcPotentialEnergy());
-    if (types&State::Forces)
+    if (types&State::Forces) {
+        impl->calcForces();
         impl->getForces().saveToArray(&state.updForces()[0]);
+    }
     if (types&State::Parameters) {
         for (map<string, double>::const_iterator iter = impl->parameters.begin(); iter != impl->parameters.end(); iter++)
             state.updParameters()[iter->first] = iter->second;

openmmapi/src/StandardMMForceFieldImpl.cpp

@@ -39,8 +39,7 @@ using std::vector;
 using std::set;
 
 StandardMMForceFieldImpl::StandardMMForceFieldImpl(StandardMMForceField& owner, OpenMMContextImpl& context) : owner(owner) {
-    forceKernel = context.getPlatform().createKernel(CalcStandardMMForcesKernel::Name());
-    energyKernel = context.getPlatform().createKernel(CalcStandardMMEnergyKernel::Name());
+    kernel = context.getPlatform().createKernel(CalcStandardMMForceFieldKernel::Name());
     vector<vector<int> > bondIndices(owner.getNumBonds());
     vector<vector<double> > bondParameters(owner.getNumBonds());
     vector<vector<int> > angleIndices(owner.getNumAngles());
@@ -79,9 +78,9 @@ StandardMMForceFieldImpl::StandardMMForceFieldImpl(StandardMMForceField& owner,
         periodicTorsionIndices[i].push_back(atom2);
         periodicTorsionIndices[i].push_back(atom3);
         periodicTorsionIndices[i].push_back(atom4);
-        periodicTorsionParameters[i].push_back(periodicity);
-        periodicTorsionParameters[i].push_back(phase);
         periodicTorsionParameters[i].push_back(k);
+        periodicTorsionParameters[i].push_back(phase);
+        periodicTorsionParameters[i].push_back(periodicity);
     }
     for (int i = 0; i < owner.getNumRBTorsions(); ++i) {
         int atom1, atom2, atom3, atom4;
@@ -113,27 +112,24 @@ StandardMMForceFieldImpl::StandardMMForceFieldImpl(StandardMMForceField& owner,
         bonded14Indices[index].push_back(iter->first);
         bonded14Indices[index++].push_back(iter->second);
     }
-    dynamic_cast<CalcStandardMMForcesKernel&>(forceKernel.getImpl()).initialize(bondIndices, bondParameters, angleIndices, angleParameters,
-            periodicTorsionIndices, periodicTorsionParameters, rbTorsionIndices, rbTorsionParameters, bonded14Indices, exclusions, nonbondedParameters);
-    dynamic_cast<CalcStandardMMEnergyKernel&>(energyKernel.getImpl()).initialize(bondIndices, bondParameters, angleIndices, angleParameters,
-            periodicTorsionIndices, periodicTorsionParameters, rbTorsionIndices, rbTorsionParameters, bonded14Indices, exclusions, nonbondedParameters);
+    dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).initialize(bondIndices, bondParameters, angleIndices, angleParameters,
+            periodicTorsionIndices, periodicTorsionParameters, rbTorsionIndices, rbTorsionParameters, bonded14Indices, 0.5, 1.0/1.2, exclusions, nonbondedParameters);
 }
 
 StandardMMForceFieldImpl::~StandardMMForceFieldImpl() {
 }
 
 void StandardMMForceFieldImpl::calcForces(OpenMMContextImpl& context, Stream& forces) {
-    dynamic_cast<CalcStandardMMForcesKernel&>(forceKernel.getImpl()).execute(context.getPositions(), forces);
+    dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).executeForces(context.getPositions(), forces);
 }
 
 double StandardMMForceFieldImpl::calcEnergy(OpenMMContextImpl& context) {
-    return dynamic_cast<CalcStandardMMEnergyKernel&>(energyKernel.getImpl()).execute(context.getPositions());
+    return dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).executeEnergy(context.getPositions());
 }
 
 std::vector<std::string> StandardMMForceFieldImpl::getKernelNames() {
     std::vector<std::string> names;
-    names.push_back(CalcStandardMMForcesKernel::Name());
-    names.push_back(CalcStandardMMEnergyKernel::Name());
+    names.push_back(CalcStandardMMForceFieldKernel::Name());
     return names;
 }
 

platforms/reference/src/ReferenceFloatStreamImpl.cpp

@@ -121,6 +121,10 @@ void ReferenceFloatStreamImpl::fillWithValue(void* value) {
 
 }
 
+const RealOpenMM* const * ReferenceFloatStreamImpl::getData() const {
+    return data;
+}
+
 RealOpenMM** ReferenceFloatStreamImpl::getData() {
     return data;
 }

platforms/reference/src/ReferenceFloatStreamImpl.h

@@ -48,6 +48,7 @@ public:
     void loadFromArray(const void* array);
     void saveToArray(void* array);
     void fillWithValue(void* value);
+    const RealOpenMM* const * getData() const;
     RealOpenMM** getData();
 private:
     int width;

platforms/reference/src/ReferenceIntStreamImpl.cpp

@@ -78,6 +78,10 @@ void ReferenceIntStreamImpl::fillWithValue(void* value) {
             data[i][j] = valueData;
 }
 
+const int* const * ReferenceIntStreamImpl::getData() const {
+    return data;
+}
+
 int** ReferenceIntStreamImpl::getData() {
     return data;
 }

platforms/reference/src/ReferenceIntStreamImpl.h

@@ -47,6 +47,7 @@ public:
     void loadFromArray(const void* array);
     void saveToArray(void* array);
     void fillWithValue(void* value);
+    const int* const * getData() const;
     int** getData();
 private:
     int width;

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -35,14 +35,10 @@
 using namespace OpenMM;
 
 KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Platform& platform) const {
-    if (name == CalcStandardMMForcesKernel::Name())
-        return new ReferenceCalcStandardMMForcesKernel(name, platform);
-    if (name == CalcStandardMMEnergyKernel::Name())
-        return new ReferenceCalcStandardMMEnergyKernel(name, platform);
-    if (name == CalcGBSAOBCForcesKernel::Name())
-        return new ReferenceCalcGBSAOBCForcesKernel(name, platform);
-    if (name == CalcGBSAOBCEnergyKernel::Name())
-        return new ReferenceCalcGBSAOBCEnergyKernel(name, platform);
+    if (name == CalcStandardMMForceFieldKernel::Name())
+        return new ReferenceCalcStandardMMForceFieldKernel(name, platform);
+    if (name == CalcGBSAOBCForceFieldKernel::Name())
+        return new ReferenceCalcGBSAOBCForceFieldKernel(name, platform);
     if (name == IntegrateVerletStepKernel::Name())
         return new ReferenceIntegrateVerletStepKernel(name, platform);
     if (name == IntegrateLangevinStepKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -30,6 +30,15 @@
  * -------------------------------------------------------------------------- */
 
 #include "ReferenceKernels.h"
+#include "ReferenceFloatStreamImpl.h"
+#include "SimTKReference/ReferenceAngleBondIxn.h"
+#include "SimTKReference/ReferenceBondForce.h"
+#include "SimTKReference/ReferenceHarmonicBondIxn.h"
+#include "SimTKReference/ReferenceLJCoulomb14.h"
+#include "SimTKReference/ReferenceLJCoulombIxn.h"
+#include "SimTKReference/ReferenceProperDihedralBond.h"
+#include "SimTKReference/ReferenceRbDihedralBond.h"
+#include <cmath>
 
 using namespace OpenMM;
 using namespace std;
@@ -84,7 +93,7 @@ void disposeRealArray(RealOpenMM** array, int size) {
     }
 }
 
-ReferenceCalcStandardMMForcesKernel::~ReferenceCalcStandardMMForcesKernel() {
+ReferenceCalcStandardMMForceFieldKernel::~ReferenceCalcStandardMMForceFieldKernel() {
     disposeIntArray(bondIndexArray, numBonds);
     disposeRealArray(bondParamArray, numBonds);
     disposeIntArray(angleIndexArray, numAngles);
@@ -93,18 +102,24 @@ ReferenceCalcStandardMMForcesKernel::~ReferenceCalcStandardMMForcesKernel() {
     disposeRealArray(periodicTorsionParamArray, numPeriodicTorsions);
     disposeIntArray(rbTorsionIndexArray, numRBTorsions);
     disposeRealArray(rbTorsionParamArray, numRBTorsions);
+    disposeRealArray(atomParamArray, numAtoms);
+    disposeIntArray(exclusionArray, numAtoms);
+    disposeIntArray(bonded14IndexArray, num14);
+    disposeRealArray(bonded14ParamArray, num14);
 }
 
-void ReferenceCalcStandardMMForcesKernel::initialize(const vector<vector<int> >& bondIndices, const vector<vector<double> >& bondParameters,
+void ReferenceCalcStandardMMForceFieldKernel::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, const vector<set<int> >& exclusions,
-        const vector<vector<double> >& nonbondedParameters) {
+        const vector<vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
+        const vector<set<int> >& exclusions, const vector<vector<double> >& nonbondedParameters) {
+    numAtoms = nonbondedParameters.size();
     numBonds = bondIndices.size();
     numAngles = angleIndices.size();
     numPeriodicTorsions = periodicTorsionIndices.size();
     numRBTorsions = rbTorsionIndices.size();
+    num14 = bonded14Indices.size();
     bondIndexArray = copyToArray(bondIndices);
     bondParamArray = copyToArray(bondParameters);
     angleIndexArray = copyToArray(angleIndices);
@@ -113,40 +128,94 @@ void ReferenceCalcStandardMMForcesKernel::initialize(const vector<vector<int> >&
     periodicTorsionParamArray = copyToArray(periodicTorsionParameters);
     rbTorsionIndexArray = copyToArray(rbTorsionIndices);
     rbTorsionParamArray = copyToArray(rbTorsionParameters);
+    atomParamArray = allocateRealArray(numAtoms, 3);
+    RealOpenMM sqrtEps = std::sqrt(138.935485);
+    for (int i = 0; i < numAtoms; ++i) {
+        atomParamArray[i][0] = 0.5*nonbondedParameters[i][1];
+        atomParamArray[i][1] = 2.0*sqrt(nonbondedParameters[i][2]);
+        atomParamArray[i][2] = nonbondedParameters[i][0]*sqrtEps;
+    }
+    exclusionArray = new int*[numAtoms];
+    for (int i = 0; i < numAtoms; ++i) {
+        exclusionArray[i] = new int[exclusions[i].size()+1];
+        exclusionArray[i][0] = exclusions[i].size();
+        int index = 0;
+        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
+            exclusionArray[i][++index] = *iter;
+    }
+    bonded14IndexArray = copyToArray(bonded14Indices);
+    bonded14ParamArray = allocateRealArray(num14, 3);
+    for (int i = 0; i < num14; ++i) {
+        int atom1 = bonded14Indices[i][0];
+        int atom2 = bonded14Indices[i][1];
+        bonded14ParamArray[i][0] = atomParamArray[atom1][0]+atomParamArray[atom2][0];
+        bonded14ParamArray[i][1] = lj14Scale*(atomParamArray[atom1][1]*atomParamArray[atom2][1]);
+        bonded14ParamArray[i][2] = coulomb14Scale*(atomParamArray[atom1][2]*atomParamArray[atom2][2]);
+    }
 }
 
-void ReferenceCalcStandardMMForcesKernel::execute(const Stream& positions, Stream& forces) {
-    
+void ReferenceCalcStandardMMForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
+    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData());
+    RealOpenMM** forceData = ((ReferenceFloatStreamImpl&) forces.getImpl()).getData();
+    ReferenceBondForce refBondForce;
+    ReferenceHarmonicBondIxn harmonicBond;
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, 0, 0, 0, harmonicBond);
+    ReferenceAngleBondIxn angleBond;
+    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, angleBond);
+    ReferenceProperDihedralBond periodicTorsionBond;
+    refBondForce.calculateForce(numPeriodicTorsions, periodicTorsionIndexArray, posData, periodicTorsionParamArray, forceData, 0, 0, 0, periodicTorsionBond);
+    ReferenceRbDihedralBond rbTorsionBond;
+    refBondForce.calculateForce(numRBTorsions, rbTorsionIndexArray, posData, rbTorsionParamArray, forceData, 0, 0, 0, rbTorsionBond);
+    ReferenceLJCoulombIxn clj;
+    clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, 0);
+    ReferenceLJCoulomb14 nonbonded14;
+    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
 }
 
-void ReferenceCalcStandardMMEnergyKernel::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, const vector<set<int> >& exclusions,
-        const vector<vector<double> >& nonbondedParameters) {
-    
+double ReferenceCalcStandardMMForceFieldKernel::executeEnergy(const Stream& positions) {
+    RealOpenMM** posData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) positions.getImpl()).getData());
+    RealOpenMM** forceData = allocateRealArray(numAtoms, 3);
+    int arraySize = max(max(max(max(numAtoms, numBonds), numAngles), numPeriodicTorsions), numRBTorsions);
+    RealOpenMM* energyArray = new RealOpenMM[arraySize];
+    RealOpenMM energy = 0;
+    ReferenceBondForce refBondForce;
+    ReferenceHarmonicBondIxn harmonicBond;
+    for (int i = 0; i < arraySize; ++i)
+        energyArray[i] = 0;
+    refBondForce.calculateForce(numBonds, bondIndexArray, posData, bondParamArray, forceData, energyArray, 0, &energy, harmonicBond);
+    ReferenceAngleBondIxn angleBond;
+    for (int i = 0; i < arraySize; ++i)
+        energyArray[i] = 0;
+    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, angleBond);
+    ReferenceProperDihedralBond periodicTorsionBond;
+    for (int i = 0; i < arraySize; ++i)
+        energyArray[i] = 0;
+    refBondForce.calculateForce(numPeriodicTorsions, periodicTorsionIndexArray, posData, periodicTorsionParamArray, forceData, energyArray, 0, &energy, periodicTorsionBond);
+    ReferenceRbDihedralBond rbTorsionBond;
+    for (int i = 0; i < arraySize; ++i)
+        energyArray[i] = 0;
+    refBondForce.calculateForce(numRBTorsions, rbTorsionIndexArray, posData, rbTorsionParamArray, forceData, energyArray, 0, &energy, rbTorsionBond);
+    ReferenceLJCoulombIxn clj;
+    clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, &energy);
+    ReferenceLJCoulomb14 nonbonded14;
+    for (int i = 0; i < arraySize; ++i)
+        energyArray[i] = 0;
+    refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);
+    disposeRealArray(forceData, numAtoms);
+    delete[] energyArray;
+    return energy;
 }
 
-double ReferenceCalcStandardMMEnergyKernel::execute(const Stream& positions) {
-    return 0.0; // TODO implement correctly
-}
-
-void ReferenceCalcGBSAOBCForcesKernel::initialize(const vector<double>& bornRadii, const vector<vector<double> >& atomParameters,
+void ReferenceCalcGBSAOBCForceFieldKernel::initialize(const vector<double>& bornRadii, const vector<vector<double> >& atomParameters,
         double solventDielectric, double soluteDielectric) {
     
 }
 
-void ReferenceCalcGBSAOBCForcesKernel::execute(const Stream& positions, Stream& forces) {
-    
-}
-
-void ReferenceCalcGBSAOBCEnergyKernel::initialize(const vector<double>& bornRadii, const vector<vector<double> >& atomParameters,
-        double solventDielectric, double soluteDielectric) {
+void ReferenceCalcGBSAOBCForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
     
 }
 
-double ReferenceCalcGBSAOBCEnergyKernel::execute(const Stream& positions) {
+double ReferenceCalcGBSAOBCForceFieldKernel::executeEnergy(const Stream& positions) {
     return 0.0; // TODO implement correctly
 }