Revised molecule identification by Cuda platform to allow plugins to specify molecule structure

platforms/cuda/src/CudaForceInfo.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) 2009 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "CudaForceInfo.h"
+
+using namespace OpenMM;
+using namespace std;
+
+bool CudaForceInfo::areParticlesIdentical(int particle1, int particle2) {
+    return true;
+}
+
+int CudaForceInfo::getNumParticleGroups() {
+    return 0;
+}
+
+void CudaForceInfo::getParticlesInGroup(int index, vector<int>& particles) {
+    return;
+}
+
+bool CudaForceInfo::areGroupsIdentical(int group1, int group2) {
+    return true;
+}

platforms/cuda/src/CudaForceInfo.h

+#ifndef OPENMM_CUDAFORCEINFO_H_
+#define OPENMM_CUDAFORCEINFO_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) 2009 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * This program is free software: you can redistribute it and/or modify       *
+ * it under the terms of the GNU Lesser General Public License as published   *
+ * by the Free Software Foundation, either version 3 of the License, or       *
+ * (at your option) any later version.                                        *
+ *                                                                            *
+ * This program is distributed in the hope that it will be useful,            *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+ * GNU Lesser General Public License for more details.                        *
+ *                                                                            *
+ * You should have received a copy of the GNU Lesser General Public License   *
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/windowsExport.h"
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * This class is used by the Cuda implementation of a Force class to convey information
+ * about the behavior and requirements of that force.
+ */
+
+class OPENMM_EXPORT CudaForceInfo {
+public:
+    CudaForceInfo() {
+    }
+    /**
+     * Get whether or not two particles have identical force field parameters.
+     */
+    virtual bool areParticlesIdentical(int particle1, int particle2);
+    /**
+     * Get the number of particle groups defined by this force.
+     */
+    virtual int getNumParticleGroups();
+    /**
+     * Get the list of particles in a particular group.
+     */
+    virtual void getParticlesInGroup(int index, std::vector<int>& particles);
+    /**
+     * Get whether two particle groups are identical.
+     */
+    virtual bool areGroupsIdentical(int group1, int group2);
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUDAFORCEINFO_H_*/

platforms/cuda/src/CudaKernels.cpp

@@ -25,6 +25,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "CudaKernels.h"
+#include "CudaForceInfo.h"
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
@@ -184,6 +185,32 @@ void CudaApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
     kApplyConstraints(data.gpu);
 }
 
+class CudaCalcHarmonicBondForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const HarmonicBondForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumBonds();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        double length, k;
+        force.getBondParameters(index, particle1, particle2, length, k);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2;
+        double length1, length2, k1, k2;
+        force.getBondParameters(group1, particle1, particle2, length1, k1);
+        force.getBondParameters(group2, particle1, particle2, length2, k2);
+        return (length1 == length2 && k1 == k2);
+    }
+private:
+    const HarmonicBondForce& force;
+};
+
 CudaCalcHarmonicBondForceKernel::~CudaCalcHarmonicBondForceKernel() {
 }
 
@@ -201,12 +228,42 @@ void CudaCalcHarmonicBondForceKernel::initialize(const System& system, const Har
         k[i] = (float) kValue;
     }
     gpuSetBondParameters(data.gpu, particle1, particle2, length, k);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcHarmonicBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     return 0.0;
 }
 
+class CudaCalcCustomBondForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const CustomBondForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumBonds();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        vector<double> parameters;
+        force.getBondParameters(index, particle1, particle2, parameters);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2;
+        vector<double> parameters1, parameters2;
+        force.getBondParameters(group1, particle1, particle2, parameters1);
+        force.getBondParameters(group2, particle1, particle2, parameters2);
+        for (int i = 0; i < (int) parameters1.size(); i++)
+            if (parameters1[i] != parameters2[i])
+                return false;
+        return true;
+    }
+private:
+    const CustomBondForce& force;
+};
+
 CudaCalcCustomBondForceKernel::~CudaCalcCustomBondForceKernel() {
 }
 
@@ -229,6 +286,7 @@ void CudaCalcCustomBondForceKernel::initialize(const System& system, const Custo
     gpuSetCustomBondParameters(data.gpu, particle1, particle2, params, force.getEnergyFunction(), paramNames, globalParamNames);
     if (globalParamValues.size() > 0)
         SetCustomBondGlobalParams(globalParamValues);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcCustomBondForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -249,6 +307,33 @@ void CudaCalcCustomBondForceKernel::updateGlobalParams(ContextImpl& context) {
         SetCustomBondGlobalParams(globalParamValues);
 }
 
+class CudaCalcHarmonicAngleForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const HarmonicAngleForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumAngles();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2, particle3;
+        double angle, k;
+        force.getAngleParameters(index, particle1, particle2, particle3, angle, k);
+        particles.resize(3);
+        particles[0] = particle1;
+        particles[1] = particle2;
+        particles[2] = particle3;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2, particle3;
+        double angle1, angle2, k1, k2;
+        force.getAngleParameters(group1, particle1, particle2, particle3, angle1, k1);
+        force.getAngleParameters(group2, particle1, particle2, particle3, angle2, k2);
+        return (angle1 == angle2 && k1 == k2);
+    }
+private:
+    const HarmonicAngleForce& force;
+};
+
 CudaCalcHarmonicAngleForceKernel::~CudaCalcHarmonicAngleForceKernel() {
 }
 
@@ -268,12 +353,43 @@ void CudaCalcHarmonicAngleForceKernel::initialize(const System& system, const Ha
         k[i] = (float) kValue;
     }
     gpuSetBondAngleParameters(data.gpu, particle1, particle2, particle3, angle, k);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcHarmonicAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     return 0.0;
 }
 
+class CudaCalcCustomAngleForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const CustomAngleForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumAngles();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2, particle3;
+        vector<double> parameters;
+        force.getAngleParameters(index, particle1, particle2, particle3, parameters);
+        particles.resize(3);
+        particles[0] = particle1;
+        particles[1] = particle2;
+        particles[2] = particle3;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2, particle3;
+        vector<double> parameters1, parameters2;
+        force.getAngleParameters(group1, particle1, particle2, particle3, parameters1);
+        force.getAngleParameters(group2, particle1, particle2, particle3, parameters2);
+        for (int i = 0; i < (int) parameters1.size(); i++)
+            if (parameters1[i] != parameters2[i])
+                return false;
+        return true;
+    }
+private:
+    const CustomAngleForce& force;
+};
+
 CudaCalcCustomAngleForceKernel::~CudaCalcCustomAngleForceKernel() {
 }
 
@@ -297,6 +413,7 @@ void CudaCalcCustomAngleForceKernel::initialize(const System& system, const Cust
     gpuSetCustomAngleParameters(data.gpu, particle1, particle2, particle3, params, force.getEnergyFunction(), paramNames, globalParamNames);
     if (globalParamValues.size() > 0)
         SetCustomAngleGlobalParams(globalParamValues);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcCustomAngleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -317,6 +434,34 @@ void CudaCalcCustomAngleForceKernel::updateGlobalParams(ContextImpl& context) {
         SetCustomAngleGlobalParams(globalParamValues);
 }
 
+class CudaCalcPeriodicTorsionForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const PeriodicTorsionForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumTorsions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2, particle3, particle4, periodicity;
+        double phase, k;
+        force.getTorsionParameters(index, particle1, particle2, particle3, particle4, periodicity, phase, k);
+        particles.resize(4);
+        particles[0] = particle1;
+        particles[1] = particle2;
+        particles[2] = particle3;
+        particles[3] = particle4;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2, particle3, particle4, periodicity1, periodicity2;
+        double phase1, phase2, k1, k2;
+        force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, periodicity1, phase1, k1);
+        force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, periodicity2, phase2, k2);
+        return (periodicity1 == periodicity2 && phase1 == phase2 && k1 == k2);
+    }
+private:
+    const PeriodicTorsionForce& force;
+};
+
 CudaCalcPeriodicTorsionForceKernel::~CudaCalcPeriodicTorsionForceKernel() {
 }
 
@@ -338,12 +483,41 @@ void CudaCalcPeriodicTorsionForceKernel::initialize(const System& system, const
         phase[i] = (float) (phaseValue*RadiansToDegrees);
     }
     gpuSetDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, k, phase, periodicity);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcPeriodicTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     return 0.0;
 }
 
+class CudaCalcRBTorsionForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const RBTorsionForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumTorsions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2, particle3, particle4;
+        double c0, c1, c2, c3, c4, c5;
+        force.getTorsionParameters(index, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
+        particles.resize(4);
+        particles[0] = particle1;
+        particles[1] = particle2;
+        particles[2] = particle3;
+        particles[3] = particle4;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2, particle3, particle4;
+        double c0a, c0b, c1a, c1b, c2a, c2b, c3a, c3b, c4a, c4b, c5a, c5b;
+        force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, c0a, c1a, c2a, c3a, c4a, c5a);
+        force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, c0b, c1b, c2b, c3b, c4b, c5b);
+        return (c0a == c0b && c1a == c1b && c2a == c2b && c3a == c3b && c4a == c4b && c5a == c5b);
+    }
+private:
+    const RBTorsionForce& force;
+};
+
 CudaCalcRBTorsionForceKernel::~CudaCalcRBTorsionForceKernel() {
 }
 
@@ -371,12 +545,43 @@ void CudaCalcRBTorsionForceKernel::initialize(const System& system, const RBTors
         c5[i] = (float) c[5];
     }
     gpuSetRbDihedralParameters(data.gpu, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcRBTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     return 0.0;
 }
 
+class CudaCalcCMAPTorsionForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const CMAPTorsionForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumTorsions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int map, a1, a2, a3, a4, b1, b2, b3, b4;
+        force.getTorsionParameters(index, map, a1, a2, a3, a4, b1, b2, b3, b4);
+        particles.resize(8);
+        particles[0] = a1;
+        particles[1] = a2;
+        particles[2] = a3;
+        particles[3] = a4;
+        particles[4] = b1;
+        particles[5] = b2;
+        particles[6] = b3;
+        particles[7] = b4;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int map1, map2, a1, a2, a3, a4, b1, b2, b3, b4;
+        force.getTorsionParameters(group1, map1, a1, a2, a3, a4, b1, b2, b3, b4);
+        force.getTorsionParameters(group2, map2, a1, a2, a3, a4, b1, b2, b3, b4);
+        return (map1 == map2);
+    }
+private:
+    const CMAPTorsionForce& force;
+};
+
 CudaCalcCMAPTorsionForceKernel::~CudaCalcCMAPTorsionForceKernel() {
     if (coefficients != NULL)
         delete coefficients;
@@ -436,6 +641,7 @@ void CudaCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAP
         maxBuffers = max(maxBuffers, forceBufferCounter[i]);
     if (maxBuffers > data.gpu->sim.outputBuffers)
         data.gpu->sim.outputBuffers = maxBuffers;
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -443,6 +649,37 @@ double CudaCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includ
     return 0.0;
 }
 
+class CudaCalcCustomTorsionForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const CustomTorsionForce& force) : force(force) {
+    }
+    int getNumParticleGroups() {
+        return force.getNumTorsions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2, particle3, particle4;
+        vector<double> parameters;
+        force.getTorsionParameters(index, particle1, particle2, particle3, particle4, parameters);
+        particles.resize(4);
+        particles[0] = particle1;
+        particles[1] = particle2;
+        particles[2] = particle3;
+        particles[3] = particle4;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2, particle3, particle4;
+        vector<double> parameters1, parameters2;
+        force.getTorsionParameters(group1, particle1, particle2, particle3, particle4, parameters1);
+        force.getTorsionParameters(group2, particle1, particle2, particle3, particle4, parameters2);
+        for (int i = 0; i < (int) parameters1.size(); i++)
+            if (parameters1[i] != parameters2[i])
+                return false;
+        return true;
+    }
+private:
+    const CustomTorsionForce& force;
+};
+
 CudaCalcCustomTorsionForceKernel::~CudaCalcCustomTorsionForceKernel() {
 }
 
@@ -467,6 +704,7 @@ void CudaCalcCustomTorsionForceKernel::initialize(const System& system, const Cu
     gpuSetCustomTorsionParameters(data.gpu, particle1, particle2, particle3, particle4, params, force.getEnergyFunction(), paramNames, globalParamNames);
     if (globalParamValues.size() > 0)
         SetCustomTorsionGlobalParams(globalParamValues);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcCustomTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -487,6 +725,38 @@ void CudaCalcCustomTorsionForceKernel::updateGlobalParams(ContextImpl& context)
         SetCustomTorsionGlobalParams(globalParamValues);
 }
 
+class CudaCalcNonbondedForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const NonbondedForce& force) : force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        double charge1, charge2, sigma1, sigma2, epsilon1, epsilon2;
+        force.getParticleParameters(particle1, charge1, sigma1, epsilon1);
+        force.getParticleParameters(particle2, charge2, sigma2, epsilon2);
+        return (charge1 == charge2 && sigma1 == sigma2 && epsilon1 == epsilon2);
+    }
+    int getNumParticleGroups() {
+        return force.getNumExceptions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        double chargeProd, sigma, epsilon;
+        force.getExceptionParameters(index, particle1, particle2, chargeProd, sigma, epsilon);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        int particle1, particle2;
+        double chargeProd1, chargeProd2, sigma1, sigma2, epsilon1, epsilon2;
+        force.getExceptionParameters(group1, particle1, particle2, chargeProd1, sigma1, epsilon1);
+        force.getExceptionParameters(group2, particle1, particle2, chargeProd2, sigma2, epsilon2);
+        return (chargeProd1 == chargeProd2 && sigma1 == sigma2 && epsilon1 == epsilon2);
+    }
+private:
+    const NonbondedForce& force;
+};
+
 CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
 }
 
@@ -594,12 +864,44 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
         }
         gpuSetLJ14Parameters(gpu, (float) ONE_4PI_EPS0, 1.0f, particle1, particle2, c6, c12, q1, q2);
     }
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     return 0.0;
 }
 
+class CudaCalcCustomNonbondedForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const CustomNonbondedForce& force) : force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        vector<double> params1;
+        vector<double> params2;
+        force.getParticleParameters(particle1, params1);
+        force.getParticleParameters(particle2, params2);
+        for (int i = 0; i < (int) params1.size(); i++)
+            if (params1[i] != params2[i])
+                return false;
+        return true;
+    }
+    int getNumParticleGroups() {
+        return force.getNumExclusions();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        int particle1, particle2;
+        force.getExclusionParticles(index, particle1, particle2);
+        particles.resize(2);
+        particles[0] = particle1;
+        particles[1] = particle2;
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        return true;
+    }
+private:
+    const CustomNonbondedForce& force;
+};
+
 CudaCalcCustomNonbondedForceKernel::~CudaCalcCustomNonbondedForceKernel() {
 }
 
@@ -657,6 +959,7 @@ void CudaCalcCustomNonbondedForceKernel::initialize(const System& system, const
     gpuSetCustomNonbondedParameters(gpu, parameters, exclusionList, method, (float) force.getCutoffDistance(), force.getEnergyFunction(), paramNames, globalParamNames);
     if (globalParamValues.size() > 0)
         SetCustomNonbondedGlobalParams(globalParamValues);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -676,6 +979,20 @@ void CudaCalcCustomNonbondedForceKernel::updateGlobalParams(ContextImpl& context
         SetCustomNonbondedGlobalParams(globalParamValues);
 }
 
+class CudaCalcGBSAOBCForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const GBSAOBCForce& force) : force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        double charge1, charge2, radius1, radius2, scale1, scale2;
+        force.getParticleParameters(particle1, charge1, radius1, scale1);
+        force.getParticleParameters(particle2, charge2, radius2, scale2);
+        return (charge1 == charge2 && radius1 == radius2 && scale1 == scale2);
+    }
+private:
+    const GBSAOBCForce& force;
+};
+
 CudaCalcGBSAOBCForceKernel::~CudaCalcGBSAOBCForceKernel() {
 }
 
@@ -694,12 +1011,27 @@ void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCF
         charge[i] = (float) particleCharge;
     }
     gpuSetObcParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), radius, scale, charge);
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
 	return 0.0;
 }
 
+class CudaCalcGBVIForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const GBVIForce& force) : force(force) {
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        double charge1, charge2, radius1, radius2, gamma1, gamma2;
+        force.getParticleParameters(particle1, charge1, radius1, gamma1);
+        force.getParticleParameters(particle2, charge2, radius2, gamma2);
+        return (charge1 == charge2 && radius1 == radius2 && gamma1 == gamma2);
+    }
+private:
+    const GBVIForce& force;
+};
+
 CudaCalcGBVIForceKernel::~CudaCalcGBVIForceKernel() {
 }
 
@@ -723,12 +1055,45 @@ void CudaCalcGBVIForceKernel::initialize(const System& system, const GBVIForce&
     }
     gpuSetGBVIParameters(gpu, (float) force.getSoluteDielectric(), (float) force.getSolventDielectric(), particle,
                          radius, gammas, scaledRadii );
+    data.gpu->forces.push_back(new ForceInfo(force));
 }
 
 double CudaCalcGBVIForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     return 0.0;
 }
 
+class CudaCalcCustomExternalForceKernel::ForceInfo : public CudaForceInfo {
+public:
+    ForceInfo(const CustomExternalForce& force, int numParticles) : force(force), indices(numParticles, -1) {
+        vector<double> params;
+        for (int i = 0; i < force.getNumParticles(); i++) {
+            int particle;
+            force.getParticleParameters(i, particle, params);
+            indices[particle] = i;
+        }
+    }
+    bool areParticlesIdentical(int particle1, int particle2) {
+        particle1 = indices[particle1];
+        particle2 = indices[particle2];
+        if (particle1 == -1 && particle2 == -1)
+            return true;
+        if (particle1 == -1 || particle2 == -1)
+            return false;
+        int temp;
+        vector<double> params1;
+        vector<double> params2;
+        force.getParticleParameters(particle1, temp, params1);
+        force.getParticleParameters(particle2, temp, params2);
+        for (int i = 0; i < (int) params1.size(); i++)
+            if (params1[i] != params2[i])
+                return false;
+        return true;
+    }
+private:
+    const CustomExternalForce& force;
+    vector<int> indices;
+};
+
 CudaCalcCustomExternalForceKernel::~CudaCalcCustomExternalForceKernel() {
 }
 
@@ -750,6 +1115,7 @@ void CudaCalcCustomExternalForceKernel::initialize(const System& system, const C
     gpuSetCustomExternalParameters(data.gpu, particle, params, force.getEnergyFunction(), paramNames, globalParamNames);
     if (globalParamValues.size() > 0)
         SetCustomExternalGlobalParams(globalParamValues);
+    data.gpu->forces.push_back(new ForceInfo(force, system.getNumParticles()));
 }
 
 double CudaCalcCustomExternalForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {

platforms/cuda/src/CudaKernels.h

@@ -208,6 +208,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     int numBonds;
     CudaPlatform::PlatformData& data;
     System& system;
@@ -239,6 +240,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     void updateGlobalParams(ContextImpl& context);
     int numBonds;
     CudaPlatform::PlatformData& data;
@@ -272,6 +274,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     int numAngles;
     CudaPlatform::PlatformData& data;
     System& system;
@@ -303,6 +306,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     void updateGlobalParams(ContextImpl& context);
     int numAngles;
     CudaPlatform::PlatformData& data;
@@ -336,6 +340,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     int numTorsions;
     CudaPlatform::PlatformData& data;
     System& system;
@@ -366,6 +371,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     int numTorsions;
     CudaPlatform::PlatformData& data;
     System& system;
@@ -398,6 +404,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     CudaPlatform::PlatformData& data;
     System& system;
     int numTorsions;
@@ -433,6 +440,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     void updateGlobalParams(ContextImpl& context);
     int numTorsions;
     CudaPlatform::PlatformData& data;
@@ -466,6 +474,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     CudaPlatform::PlatformData& data;
     int numParticles;
     System& system;
@@ -496,6 +505,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     void updateGlobalParams(ContextImpl& context);
     CudaPlatform::PlatformData& data;
     int numParticles;
@@ -529,6 +539,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     CudaPlatform::PlatformData& data;
 };
 
@@ -558,6 +569,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     CudaPlatform::PlatformData& data;
 };
 
@@ -587,6 +599,7 @@ public:
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
 private:
+    class ForceInfo;
     void updateGlobalParams(ContextImpl& context);
     int numParticles;
     CudaPlatform::PlatformData& data;

platforms/cuda/src/kernels/gpu.cpp

@@ -53,6 +53,7 @@ using namespace std;
 #include "quern.h"
 #include "Lepton.h"
 #include "rng.h"
+#include "../CudaForceInfo.h"
 
 // In case we're using some primitive version of Visual Studio this will
 // make sure that erf() and erfc() are defined.
@@ -106,15 +107,8 @@ struct ConstraintOrderer : public binary_function<int, int, bool> {
 
 struct Molecule {
     vector<int> atoms;
-    vector<int> bonds;
-    vector<int> angles;
-    vector<int> periodicTorsions;
-    vector<int> rbTorsions;
     vector<int> constraints;
-    vector<int> lj14s;
-    vector<int> customBonds;
-    vector<int> customAngles;
-    vector<int> customTorsions;
+    vector<vector<int> > groups;
 };
 
 static const float dielectricOffset         =    0.009f;
@@ -2534,12 +2528,15 @@ static void findMoleculeGroups(gpuContext gpu)
 
     int numAtoms = gpu->natoms;
     vector<vector<int> > atomBonds(numAtoms);
-    for (int i = 0; i < (int)gpu->sim.bonds; i++)
-    {
-        int atom1 = (*gpu->psBondID)[i].x;
-        int atom2 = (*gpu->psBondID)[i].y;
-        atomBonds[atom1].push_back(atom2);
-        atomBonds[atom2].push_back(atom1);
+    for (int i = 0; i < (int) gpu->forces.size(); i++) {
+        for (int j = 0; j < gpu->forces[i]->getNumParticleGroups(); j++) {
+            vector<int> particles;
+            gpu->forces[i]->getParticlesInGroup(j, particles);
+            for (int k = 0; k < (int) particles.size(); k++)
+                for (int m = 0; m < (int) particles.size(); m++)
+                    if (k != m)
+                        atomBonds[particles[k]].push_back(particles[m]);
+        }
     }
     for (int i = 0; i < (int)constraints.size(); i++)
     {
@@ -2548,9 +2545,6 @@ static void findMoleculeGroups(gpuContext gpu)
         atomBonds[atom1].push_back(atom2);
         atomBonds[atom2].push_back(atom1);
     }
-    for (int i = 0; i < (int)gpu->exclusions.size(); i++)
-        for (int j = 0; j < (int)gpu->exclusions[i].size(); j++)
-            atomBonds[i].push_back(gpu->exclusions[i][j]);
 
     // Now tag atoms by which molecule they belong to.
 
@@ -2567,51 +2561,21 @@ static void findMoleculeGroups(gpuContext gpu)
 
     vector<Molecule> molecules(numMolecules);
     for (int i = 0; i < numMolecules; i++)
-        molecules[i].atoms = atomIndices[i];
-    for (int i = 0; i < (int)gpu->sim.bonds; i++)
-    {
-        int atom1 = (*gpu->psBondID)[i].x;
-        molecules[atomMolecule[atom1]].bonds.push_back(i);
-    }
-    for (int i = 0; i < (int)gpu->sim.bond_angles; i++)
-    {
-        int atom1 = (*gpu->psBondAngleID1)[i].x;
-        molecules[atomMolecule[atom1]].angles.push_back(i);
-    }
-    for (int i = 0; i < (int)gpu->sim.dihedrals; i++)
-    {
-        int atom1 = (*gpu->psDihedralID1)[i].x;
-        molecules[atomMolecule[atom1]].periodicTorsions.push_back(i);
-    }
-    for (int i = 0; i < (int)gpu->sim.rb_dihedrals; i++)
     {
-        int atom1 = (*gpu->psRbDihedralID1)[i].x;
-        molecules[atomMolecule[atom1]].rbTorsions.push_back(i);
+        molecules[i].atoms = atomIndices[i];
+        molecules[i].groups.resize(gpu->forces.size());
     }
+    for (int i = 0; i < (int) gpu->forces.size(); i++)
+        for (int j = 0; j < gpu->forces[i]->getNumParticleGroups(); j++)
+        {
+            vector<int> particles;
+            gpu->forces[i]->getParticlesInGroup(j, particles);
+            molecules[atomMolecule[particles[0]]].groups[i].push_back(j);
+        }
     for (int i = 0; i < (int)constraints.size(); i++)
     {
         molecules[atomMolecule[constraints[i].atom1]].constraints.push_back(i);
     }
-    for (int i = 0; i < (int)gpu->sim.LJ14s; i++)
-    {
-        int atom1 = (*gpu->psLJ14ID)[i].x;
-        molecules[atomMolecule[atom1]].lj14s.push_back(i);
-    }
-    for (int i = 0; i < (int)gpu->sim.customBonds; i++)
-    {
-        int atom1 = (*gpu->psCustomBondID)[i].x;
-        molecules[atomMolecule[atom1]].customBonds.push_back(i);
-    }
-    for (int i = 0; i < (int)gpu->sim.customAngles; i++)
-    {
-        int atom1 = (*gpu->psCustomAngleID1)[i].x;
-        molecules[atomMolecule[atom1]].customAngles.push_back(i);
-    }
-    for (int i = 0; i < (int)gpu->sim.customTorsions; i++)
-    {
-        int atom1 = (*gpu->psCustomTorsionID1)[i].x;
-        molecules[atomMolecule[atom1]].customTorsions.push_back(i);
-    }
 
     // Sort them into groups of identical molecules.
 
@@ -2627,112 +2591,36 @@ static void findMoleculeGroups(gpuContext gpu)
         for (int j = 0; j < (int)uniqueMolecules.size() && isNew; j++)
         {
             Molecule& mol2 = uniqueMolecules[j];
-            bool identical = true;
-            if (mol.atoms.size() != mol2.atoms.size() || mol.bonds.size() != mol2.bonds.size()
-                    || mol.angles.size() != mol2.angles.size() || mol.periodicTorsions.size() != mol2.periodicTorsions.size()
-                    || mol.rbTorsions.size() != mol2.rbTorsions.size() || mol.constraints.size() != mol2.constraints.size()
-                    || mol.lj14s.size() != mol2.lj14s.size() || mol.customBonds.size() != mol2.customBonds.size()
-                    || mol.customAngles.size() != mol2.customAngles.size() || mol.customTorsions.size() != mol2.customTorsions.size())
-                identical = false;
+            bool identical = (mol.atoms.size() == mol2.atoms.size() && mol.constraints.size() == mol2.constraints.size());
+
+            // See if the atoms are identical.
+
             int atomOffset = mol2.atoms[0]-mol.atoms[0];
-            float4* posq = gpu->psPosq4->_pSysData;
             float4* velm = gpu->psVelm4->_pSysData;
-            float2* sigeps = gpu->psSigEps2->_pSysData;
-            for (int i = 0; i < (int)mol.atoms.size() && identical; i++)
-                if (mol.atoms[i] != mol2.atoms[i]-atomOffset || posq[mol.atoms[i]].w != posq[mol2.atoms[i]].w ||
-                        velm[mol.atoms[i]].w != velm[mol2.atoms[i]].w || sigeps[mol.atoms[i]].x != sigeps[mol2.atoms[i]].x ||
-                        sigeps[mol.atoms[i]].y != sigeps[mol2.atoms[i]].y)
-                    identical = false;
-            int4* bondID = gpu->psBondID->_pSysData;
-            float2* bondParam = gpu->psBondParameter->_pSysData;
-            for (int i = 0; i < (int)mol.bonds.size() && identical; i++)
-                if (bondID[mol.bonds[i]].x != bondID[mol2.bonds[i]].x-atomOffset || bondID[mol.bonds[i]].y != bondID[mol2.bonds[i]].y-atomOffset ||
-                        bondParam[mol.bonds[i]].x != bondParam[mol2.bonds[i]].x || bondParam[mol.bonds[i]].y != bondParam[mol2.bonds[i]].y)
+            for (int i = 0; i < (int) mol.atoms.size() && identical; i++) {
+                if (mol.atoms[i] != mol2.atoms[i]-atomOffset || velm[mol.atoms[i]].w != velm[mol2.atoms[i]].w)
                     identical = false;
-            int4* angleID = gpu->psBondAngleID1->_pSysData;
-            float2* angleParam = gpu->psBondAngleParameter->_pSysData;
-            for (int i = 0; i < (int)mol.angles.size() && identical; i++)
-                if (angleID[mol.angles[i]].x != angleID[mol2.angles[i]].x-atomOffset ||
-                        angleID[mol.angles[i]].y != angleID[mol2.angles[i]].y-atomOffset ||
-                        angleID[mol.angles[i]].z != angleID[mol2.angles[i]].z-atomOffset ||
-                        angleParam[mol.angles[i]].x != angleParam[mol2.angles[i]].x ||
-                        angleParam[mol.angles[i]].y != angleParam[mol2.angles[i]].y)
-                    identical = false;
-            int4* periodicID = gpu->psDihedralID1->_pSysData;
-            float4* periodicParam = gpu->psDihedralParameter->_pSysData;
-            for (int i = 0; i < (int)mol.periodicTorsions.size() && identical; i++)
-                if (periodicID[mol.periodicTorsions[i]].x != periodicID[mol2.periodicTorsions[i]].x-atomOffset ||
-                        periodicID[mol.periodicTorsions[i]].y != periodicID[mol2.periodicTorsions[i]].y-atomOffset ||
-                        periodicID[mol.periodicTorsions[i]].z != periodicID[mol2.periodicTorsions[i]].z-atomOffset ||
-                        periodicID[mol.periodicTorsions[i]].w != periodicID[mol2.periodicTorsions[i]].w-atomOffset ||
-                        periodicParam[mol.periodicTorsions[i]].x != periodicParam[mol2.periodicTorsions[i]].x ||
-                        periodicParam[mol.periodicTorsions[i]].y != periodicParam[mol2.periodicTorsions[i]].y ||
-                        periodicParam[mol.periodicTorsions[i]].z != periodicParam[mol2.periodicTorsions[i]].z)
-                    identical = false;
-            int4* rbID = gpu->psRbDihedralID1->_pSysData;
-            float4* rbParam1 = gpu->psRbDihedralParameter1->_pSysData;
-            float2* rbParam2 = gpu->psRbDihedralParameter2->_pSysData;
-            for (int i = 0; i < (int)mol.rbTorsions.size() && identical; i++)
-                if (rbID[mol.rbTorsions[i]].x != rbID[mol2.rbTorsions[i]].x-atomOffset ||
-                        rbID[mol.rbTorsions[i]].y != rbID[mol2.rbTorsions[i]].y-atomOffset ||
-                        rbID[mol.rbTorsions[i]].z != rbID[mol2.rbTorsions[i]].z-atomOffset ||
-                        rbID[mol.rbTorsions[i]].w != rbID[mol2.rbTorsions[i]].w-atomOffset ||
-                        rbParam1[mol.rbTorsions[i]].x != rbParam1[mol2.rbTorsions[i]].x ||
-                        rbParam1[mol.rbTorsions[i]].y != rbParam1[mol2.rbTorsions[i]].y ||
-                        rbParam1[mol.rbTorsions[i]].z != rbParam1[mol2.rbTorsions[i]].z ||
-                        rbParam1[mol.rbTorsions[i]].w != rbParam1[mol2.rbTorsions[i]].w ||
-                        rbParam2[mol.rbTorsions[i]].x != rbParam2[mol2.rbTorsions[i]].x ||
-                        rbParam2[mol.rbTorsions[i]].y != rbParam2[mol2.rbTorsions[i]].y)
-                    identical = false;
-            for (int i = 0; i < (int)mol.constraints.size() && identical; i++)
+                for (int k = 0; k < (int) gpu->forces.size(); k++)
+                    if (!gpu->forces[k]->areParticlesIdentical(mol.atoms[i], mol2.atoms[i]))
+                        identical = false;
+            }
+
+            // See if the constraints are identical.
+
+            for (int i = 0; i < (int) mol.constraints.size() && identical; i++)
                 if (constraints[mol.constraints[i]].atom1 != constraints[mol2.constraints[i]].atom1-atomOffset ||
                         constraints[mol.constraints[i]].atom2 != constraints[mol2.constraints[i]].atom2-atomOffset ||
                         constraints[mol.constraints[i]].distance2 != constraints[mol2.constraints[i]].distance2)
                     identical = false;
-            int4* lj14ID = gpu->psLJ14ID->_pSysData;
-            float4* lj14Param = gpu->psLJ14Parameter->_pSysData;
-            for (int i = 0; i < (int)mol.lj14s.size() && identical; i++)
-                if (lj14ID[mol.lj14s[i]].x != lj14ID[mol2.lj14s[i]].x-atomOffset || lj14ID[mol.lj14s[i]].y != lj14ID[mol2.lj14s[i]].y-atomOffset ||
-                        lj14Param[mol.lj14s[i]].x != lj14Param[mol2.lj14s[i]].x || lj14Param[mol.lj14s[i]].y != lj14Param[mol2.lj14s[i]].y ||
-                        lj14Param[mol.lj14s[i]].z != lj14Param[mol2.lj14s[i]].z)
+
+            // See if the force groups are identical.
+
+            for (int i = 0; i < (int) gpu->forces.size() && identical; i++)
+            {
+                if (mol.groups[i].size() != mol2.groups[i].size())
                     identical = false;
-            if (mol.customBonds.size() > 0) {
-                int4* customBondID = gpu->psCustomBondID->_pSysData;
-                float4* customBondParam = gpu->psCustomBondParams->_pSysData;
-                for (int i = 0; i < (int)mol.customBonds.size() && identical; i++)
-                    if (customBondID[mol.customBonds[i]].x != customBondID[mol2.customBonds[i]].x-atomOffset ||
-                            customBondID[mol.customBonds[i]].y != customBondID[mol2.customBonds[i]].y-atomOffset ||
-                            (customBondParam[mol.customBonds[i]].x != customBondParam[mol2.customBonds[i]].x && gpu->sim.customBondParameters > 0) ||
-                            (customBondParam[mol.customBonds[i]].y != customBondParam[mol2.customBonds[i]].y && gpu->sim.customBondParameters > 1) ||
-                            (customBondParam[mol.customBonds[i]].z != customBondParam[mol2.customBonds[i]].z && gpu->sim.customBondParameters > 2) ||
-                            (customBondParam[mol.customBonds[i]].w != customBondParam[mol2.customBonds[i]].w && gpu->sim.customBondParameters > 3))
-                        identical = false;
-            }
-            if (mol.customAngles.size() > 0) {
-                int4* customAngleID = gpu->psCustomAngleID1->_pSysData;
-                float4* customAngleParam = gpu->psCustomAngleParams->_pSysData;
-                for (int i = 0; i < (int)mol.customAngles.size() && identical; i++)
-                    if (customAngleID[mol.customAngles[i]].x != customAngleID[mol2.customAngles[i]].x-atomOffset ||
-                            customAngleID[mol.customAngles[i]].y != customAngleID[mol2.customAngles[i]].y-atomOffset ||
-                            customAngleID[mol.customAngles[i]].z != customAngleID[mol2.customAngles[i]].z-atomOffset ||
-                            (customAngleParam[mol.customAngles[i]].x != customAngleParam[mol2.customAngles[i]].x && gpu->sim.customAngleParameters > 0) ||
-                            (customAngleParam[mol.customAngles[i]].y != customAngleParam[mol2.customAngles[i]].y && gpu->sim.customAngleParameters > 1) ||
-                            (customAngleParam[mol.customAngles[i]].z != customAngleParam[mol2.customAngles[i]].z && gpu->sim.customAngleParameters > 2) ||
-                            (customAngleParam[mol.customAngles[i]].w != customAngleParam[mol2.customAngles[i]].w && gpu->sim.customAngleParameters > 3))
-                        identical = false;
-            }
-            if (mol.customTorsions.size() > 0) {
-                int4* customTorsionID = gpu->psCustomTorsionID1->_pSysData;
-                float4* customTorsionParam = gpu->psCustomTorsionParams->_pSysData;
-                for (int i = 0; i < (int)mol.customTorsions.size() && identical; i++)
-                    if (customTorsionID[mol.customTorsions[i]].x != customTorsionID[mol2.customTorsions[i]].x-atomOffset ||
-                            customTorsionID[mol.customTorsions[i]].y != customTorsionID[mol2.customTorsions[i]].y-atomOffset ||
-                            customTorsionID[mol.customTorsions[i]].z != customTorsionID[mol2.customTorsions[i]].z-atomOffset ||
-                            customTorsionID[mol.customTorsions[i]].w != customTorsionID[mol2.customTorsions[i]].w-atomOffset ||
-                            (customTorsionParam[mol.customTorsions[i]].x != customTorsionParam[mol2.customTorsions[i]].x && gpu->sim.customTorsionParameters > 0) ||
-                            (customTorsionParam[mol.customTorsions[i]].y != customTorsionParam[mol2.customTorsions[i]].y && gpu->sim.customTorsionParameters > 1) ||
-                            (customTorsionParam[mol.customTorsions[i]].z != customTorsionParam[mol2.customTorsions[i]].z && gpu->sim.customTorsionParameters > 2) ||
-                            (customTorsionParam[mol.customTorsions[i]].w != customTorsionParam[mol2.customTorsions[i]].w && gpu->sim.customTorsionParameters > 3))
+                for (int k = 0; k < (int) mol.groups[i].size() && identical; k++)
+                    if (!gpu->forces[i]->areGroupsIdentical(mol.groups[i][k], mol2.groups[i][k]))
                         identical = false;
             }
             if (identical)

platforms/cuda/src/kernels/gputypes.h

@@ -32,6 +32,10 @@
 #include <vector>
 #include "windowsExportCuda.h"
 
+namespace OpenMM {
+    class CudaForceInfo;
+}
+
 struct gpuAtomType {
     std::string name;
     char symbol;
@@ -75,6 +79,7 @@ struct _gpuContext {
     unsigned int sharedMemoryPerBlock;
     cudaGmxSimulation sim;
     unsigned int* pOutputBufferCounter;
+    std::vector<OpenMM::CudaForceInfo*> forces;
     std::vector<std::vector<int> > exclusions;
     unsigned char* pAtomSymbol;
     std::vector<gpuMoleculeGroup> moleculeGroups;