Implemented updateParametersInContext() for five more Force classes

openmmapi/src/CustomNonbondedForceImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -123,3 +123,7 @@ map<string, double> CustomNonbondedForceImpl::getDefaultParameters() {
         parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
     return parameters;
 }
+
+void CustomNonbondedForceImpl::updateParametersInContext(ContextImpl& context) {
+    kernel.getAs<CalcCustomNonbondedForceKernel>().copyParametersToContext(context, owner);
+}

platforms/cuda/src/CudaKernels.cpp

@@ -1033,6 +1033,10 @@ void CudaCalcCustomNonbondedForceKernel::updateGlobalParams(ContextImpl& context
         SetCustomNonbondedGlobalParams(globalParamValues);
 }
 
+void CudaCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force) {
+    throw OpenMMException("CudaPlatform does not support copyParametersToContext");
+}
+
 class CudaCalcGBSAOBCForceKernel::ForceInfo : public CudaForceInfo {
 public:
     ForceInfo(const GBSAOBCForce& force) : force(force) {
@@ -1204,6 +1208,10 @@ void CudaCalcCustomExternalForceKernel::updateGlobalParams(ContextImpl& context)
         SetCustomExternalGlobalParams(globalParamValues);
 }
 
+void CudaCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl& context, const CustomExternalForce& force) {
+    throw OpenMMException("CudaPlatform does not support copyParametersToContext");
+}
+
 void OPENMMCUDA_EXPORT OpenMM::cudaOpenMMInitializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
 
     // Initialize any terms that haven't already been handled by a Force.

platforms/cuda/src/CudaKernels.h

@@ -596,6 +596,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomNonbondedForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
 private:
     class ForceInfo;
     void updateGlobalParams(ContextImpl& context);
@@ -697,6 +704,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomNonbondedForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
 private:
     class ForceInfo;
     void updateGlobalParams(ContextImpl& context);

platforms/opencl/src/OpenCLKernels.cpp

@@ -98,11 +98,11 @@ void OpenCLCalcForcesAndEnergyKernel::beginComputation(ContextImpl& context, boo
     OpenCLNonbondedUtilities& nb = cl.getNonbondedUtilities();
     bool includeNonbonded = ((groups&(1<<nb.getForceGroup())) != 0);
     cl.setAtomsWereReordered(false);
-    if (cl.getMoleculesAreInvalid() || (nb.getUseCutoff() && includeNonbonded && cl.getComputeForceCount()%100 == 0)) {
+    if (nb.getUseCutoff() && includeNonbonded && (cl.getMoleculesAreInvalid() || cl.getComputeForceCount()%100 == 0)) {
         cl.reorderAtoms(!cl.getMoleculesAreInvalid());
         nb.updateNeighborListSize();
-        cl.setComputeForceCount(cl.getComputeForceCount()+1);
     }
+    cl.setComputeForceCount(cl.getComputeForceCount()+1);
     cl.clearAutoclearBuffers();
     if (includeNonbonded)
         nb.prepareInteractions();
@@ -1769,6 +1769,28 @@ double OpenCLCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bool
     return 0.0;
 }
 
+void OpenCLCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force) {
+    int numParticles = force.getNumParticles();
+    if (numParticles != cl.getNumAtoms())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    
+    // Record the per-particle parameters.
+    
+    vector<vector<cl_float> > paramVector(numParticles);
+    vector<double> parameters;
+    for (int i = 0; i < numParticles; i++) {
+        force.getParticleParameters(i, parameters);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (cl_float) parameters[j];
+    }
+    params->setParameterValues(paramVector);
+    
+    // Mark that the current reordering may be invalid.
+    
+    cl.invalidateMolecules();
+}
+
 class OpenCLGBSAOBCForceInfo : public OpenCLForceInfo {
 public:
     OpenCLGBSAOBCForceInfo(int requiredBuffers, const GBSAOBCForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
@@ -1946,7 +1968,7 @@ double OpenCLCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeF
         force1Kernel.setArg<mm_float4>(8, cl.getInvPeriodicBoxSize());
         if (maxTiles < nb.getInteractingTiles().getSize()) {
             maxTiles = nb.getInteractingTiles().getSize();
-            computeBornSumKernel.setArg<cl::Buffer>(4, nb.getInteractingTiles().getDeviceBuffer());
+            computeBornSumKernel.setArg<cl::Buffer>(3, nb.getInteractingTiles().getDeviceBuffer());
             computeBornSumKernel.setArg<cl_uint>(7, maxTiles);
             force1Kernel.setArg<cl::Buffer>(5, nb.getInteractingTiles().getDeviceBuffer());
             force1Kernel.setArg<cl_uint>(9, maxTiles);
@@ -2901,6 +2923,28 @@ double OpenCLCalcCustomGBForceKernel::execute(ContextImpl& context, bool include
     return 0.0;
 }
 
+void OpenCLCalcCustomGBForceKernel::copyParametersToContext(ContextImpl& context, const CustomGBForce& force) {
+    int numParticles = force.getNumParticles();
+    if (numParticles != cl.getNumAtoms())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    
+    // Record the per-particle parameters.
+    
+    vector<vector<cl_float> > paramVector(numParticles);
+    vector<double> parameters;
+    for (int i = 0; i < numParticles; i++) {
+        force.getParticleParameters(i, parameters);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (cl_float) parameters[j];
+    }
+    params->setParameterValues(paramVector);
+    
+    // Mark that the current reordering may be invalid.
+    
+    cl.invalidateMolecules();
+}
+
 class OpenCLCustomExternalForceInfo : public OpenCLForceInfo {
 public:
     OpenCLCustomExternalForceInfo(const CustomExternalForce& force, int numParticles) : OpenCLForceInfo(0), force(force), indices(numParticles, -1) {
@@ -3026,6 +3070,31 @@ double OpenCLCalcCustomExternalForceKernel::execute(ContextImpl& context, bool i
     return 0.0;
 }
 
+void OpenCLCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl& context, const CustomExternalForce& force) {
+    int numContexts = cl.getPlatformData().contexts.size();
+    int startIndex = cl.getContextIndex()*force.getNumParticles()/numContexts;
+    int endIndex = (cl.getContextIndex()+1)*force.getNumParticles()/numContexts;
+    if (numParticles != endIndex-startIndex)
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+    
+    // Record the per-particle parameters.
+    
+    vector<vector<cl_float> > paramVector(numParticles);
+    vector<double> parameters;
+    for (int i = 0; i < numParticles; i++) {
+        int particle;
+        force.getParticleParameters(startIndex+i, particle, parameters);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (cl_float) parameters[j];
+    }
+    params->setParameterValues(paramVector);
+    
+    // Mark that the current reordering may be invalid.
+    
+    cl.invalidateMolecules();
+}
+
 class OpenCLCustomHbondForceInfo : public OpenCLForceInfo {
 public:
     OpenCLCustomHbondForceInfo(int requiredBuffers, const CustomHbondForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {
@@ -3536,6 +3605,45 @@ double OpenCLCalcCustomHbondForceKernel::execute(ContextImpl& context, bool incl
     return 0.0;
 }
 
+void OpenCLCalcCustomHbondForceKernel::copyParametersToContext(ContextImpl& context, const CustomHbondForce& force) {
+    int numContexts = cl.getPlatformData().contexts.size();
+    int startIndex = cl.getContextIndex()*force.getNumDonors()/numContexts;
+    int endIndex = (cl.getContextIndex()+1)*force.getNumDonors()/numContexts;
+    if (numDonors != endIndex-startIndex)
+        throw OpenMMException("updateParametersInContext: The number of donors has changed");
+    if (numAcceptors != force.getNumAcceptors())
+        throw OpenMMException("updateParametersInContext: The number of acceptors has changed");
+    
+    // Record the per-donor parameters.
+    
+    vector<vector<cl_float> > donorParamVector(numDonors);
+    vector<double> parameters;
+    for (int i = 0; i < numDonors; i++) {
+        int d1, d2, d3;
+        force.getDonorParameters(startIndex+i, d1, d2, d3, parameters);
+        donorParamVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            donorParamVector[i][j] = (cl_float) parameters[j];
+    }
+    donorParams->setParameterValues(donorParamVector);
+    
+    // Record the per-acceptor parameters.
+    
+    vector<vector<cl_float> > acceptorParamVector(numAcceptors);
+    for (int i = 0; i < numAcceptors; i++) {
+        int a1, a2, a3;
+        force.getAcceptorParameters(i, a1, a2, a3, parameters);
+        acceptorParamVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            acceptorParamVector[i][j] = (cl_float) parameters[j];
+    }
+    acceptorParams->setParameterValues(acceptorParamVector);
+    
+    // Mark that the current reordering may be invalid.
+    
+    cl.invalidateMolecules();
+}
+
 class OpenCLCustomCompoundBondForceInfo : public OpenCLForceInfo {
 public:
     OpenCLCustomCompoundBondForceInfo(const CustomCompoundBondForce& force) : OpenCLForceInfo(0), force(force) {
@@ -3828,6 +3936,31 @@ double OpenCLCalcCustomCompoundBondForceKernel::execute(ContextImpl& context, bo
     return 0.0;
 }
 
+void OpenCLCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force) {
+    int numContexts = cl.getPlatformData().contexts.size();
+    int startIndex = cl.getContextIndex()*force.getNumBonds()/numContexts;
+    int endIndex = (cl.getContextIndex()+1)*force.getNumBonds()/numContexts;
+    if (numBonds != endIndex-startIndex)
+        throw OpenMMException("updateParametersInContext: The number of bonds has changed");
+    
+    // Record the per-bond parameters.
+    
+    vector<vector<cl_float> > paramVector(numBonds);
+    vector<int> particles;
+    vector<double> parameters;
+    for (int i = 0; i < numBonds; i++) {
+        force.getBondParameters(startIndex+i, particles, parameters);
+        paramVector[i].resize(parameters.size());
+        for (int j = 0; j < (int) parameters.size(); j++)
+            paramVector[i][j] = (cl_float) parameters[j];
+    }
+    params->setParameterValues(paramVector);
+    
+    // Mark that the current reordering may be invalid.
+    
+    cl.invalidateMolecules();
+}
+
 OpenCLIntegrateVerletStepKernel::~OpenCLIntegrateVerletStepKernel() {
 }
 

platforms/opencl/src/OpenCLKernels.h

@@ -647,6 +647,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomNonbondedForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
 private:
     bool hasInitializedKernel;
     OpenCLContext& cl;
@@ -736,6 +743,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomGBForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomGBForce& force);
 private:
     bool hasInitializedKernels, needParameterGradient;
     int maxTiles, numComputedValues;
@@ -781,6 +795,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomExternalForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
 private:
     int numParticles;
     bool hasInitializedKernel;
@@ -819,6 +840,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomHbondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomHbondForce& force);
 private:
     int numDonors, numAcceptors;
     bool hasInitializedKernel;
@@ -865,6 +893,14 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomCompoundBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
+
 private:
     int numBonds;
     OpenCLContext& cl;

platforms/opencl/src/OpenCLParallelKernels.cpp

@@ -598,6 +598,11 @@ double OpenCLParallelCalcCustomNonbondedForceKernel::execute(ContextImpl& contex
     return 0.0;
 }
 
+void OpenCLParallelCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force) {
+    for (int i = 0; i < (int) kernels.size(); i++)
+        getKernel(i).copyParametersToContext(context, force);
+}
+
 class OpenCLParallelCalcCustomExternalForceKernel::Task : public OpenCLContext::WorkTask {
 public:
     Task(ContextImpl& context, OpenCLCalcCustomExternalForceKernel& kernel, bool includeForce,
@@ -634,6 +639,11 @@ double OpenCLParallelCalcCustomExternalForceKernel::execute(ContextImpl& context
     return 0.0;
 }
 
+void OpenCLParallelCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl& context, const CustomExternalForce& force) {
+    for (int i = 0; i < (int) kernels.size(); i++)
+        getKernel(i).copyParametersToContext(context, force);
+}
+
 class OpenCLParallelCalcCustomHbondForceKernel::Task : public OpenCLContext::WorkTask {
 public:
     Task(ContextImpl& context, OpenCLCalcCustomHbondForceKernel& kernel, bool includeForce,
@@ -670,6 +680,11 @@ double OpenCLParallelCalcCustomHbondForceKernel::execute(ContextImpl& context, b
     return 0.0;
 }
 
+void OpenCLParallelCalcCustomHbondForceKernel::copyParametersToContext(ContextImpl& context, const CustomHbondForce& force) {
+    for (int i = 0; i < (int) kernels.size(); i++)
+        getKernel(i).copyParametersToContext(context, force);
+}
+
 class OpenCLParallelCalcCustomCompoundBondForceKernel::Task : public OpenCLContext::WorkTask {
 public:
     Task(ContextImpl& context, OpenCLCalcCustomCompoundBondForceKernel& kernel, bool includeForce,
@@ -705,3 +720,8 @@ double OpenCLParallelCalcCustomCompoundBondForceKernel::execute(ContextImpl& con
     }
     return 0.0;
 }
+
+void OpenCLParallelCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force) {
+    for (int i = 0; i < (int) kernels.size(); i++)
+        getKernel(i).copyParametersToContext(context, force);
+}

platforms/opencl/src/OpenCLParallelKernels.h

@@ -450,6 +450,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomNonbondedForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
 private:
     class Task;
     OpenCLPlatform::PlatformData& data;
@@ -481,6 +488,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomExternalForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
 private:
     class Task;
     OpenCLPlatform::PlatformData& data;
@@ -512,6 +526,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomHbondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomHbondForce& force);
 private:
     class Task;
     OpenCLPlatform::PlatformData& data;
@@ -543,6 +564,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomCompoundBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
 private:
     class Task;
     OpenCLPlatform::PlatformData& data;

platforms/opencl/tests/TestOpenCLCustomCompoundBondForce.cpp

@@ -121,6 +121,24 @@ void testBond() {
             ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
         ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
     }
+    
+    // Try changing the bond parameters and make sure it's still correct.
+    
+    parameters[0] = 1.6;
+    parameters[3] = 1.3;
+    custom->setBondParameters(0, particles, parameters);
+    custom->updateParametersInContext(c1);
+    bonds->setBondParameters(0, 0, 1, 1.3, 1.6);
+    bonds->setBondParameters(1, 1, 3, 1.3, 1.6);
+    bonds->updateParametersInContext(c2);
+    {
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = s1.getForces();
+        for (int i = 0; i < customSystem.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
+        ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
+    }
 }
 
 void testPositionDependence() {

platforms/opencl/tests/TestOpenCLCustomExternalForce.cpp

@@ -75,11 +75,28 @@ void testForce() {
     positions[2] = Vec3(1, 0, 1);
     context.setPositions(positions);
     State state = context.getState(State::Forces | State::Energy);
+    {
         const vector<Vec3>& forces = state.getForces();
         ASSERT_EQUAL_VEC(Vec3(-0.5, -0.5*2.0*2.0*1.5, 0), forces[0], TOL);
         ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
         ASSERT_EQUAL_VEC(Vec3(-0.5, 0.5*3.0*2.0*1.5, 0), forces[2], TOL);
         ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.0*1.5*1.5), state.getPotentialEnergy(), TOL);
+    }
+    
+    // Try changing the parameters and make sure it's still correct.
+    
+    parameters[0] = 1.4;
+    parameters[1] = 3.5;
+    forceField->setParticleParameters(1, 2, parameters);
+    forceField->updateParametersInContext(context);
+    state = context.getState(State::Forces | State::Energy);
+    {
+        const vector<Vec3>& forces = state.getForces();
+        ASSERT_EQUAL_VEC(Vec3(-0.5, -0.5*2.0*2.0*1.5, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
+        ASSERT_EQUAL_VEC(Vec3(-0.5, 0.5*3.5*2.0*1.4, 0), forces[2], TOL);
+        ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.5*1.4*1.4), state.getPotentialEnergy(), TOL);
+    }
 }
 
 void testManyParameters() {

platforms/opencl/tests/TestOpenCLCustomGBForce.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -139,6 +139,29 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
     for (int i = 0; i < numParticles; i++) {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
     }
+    
+    // Try changing the particle parameters and make sure it's still correct.
+    
+    for (int i = 0; i < numMolecules/2; i++) {
+        obc->setParticleParameters(2*i, 1.1, 0.3, 0.6);
+        params[0] = 1.1;
+        params[1] = 0.3;
+        params[2] = 0.6;
+        custom->setParticleParameters(2*i, params);
+        obc->setParticleParameters(2*i+1, -1.1, 0.2, 0.4);
+        params[0] = -1.1;
+        params[1] = 0.2;
+        params[2] = 0.4;
+        custom->setParticleParameters(2*i+1, params);
+    }
+    obc->updateParametersInContext(context1);
+    custom->updateParametersInContext(context2);
+    state1 = context1.getState(State::Forces | State::Energy);
+    state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
+    }
 }
 
 void testMembrane() {

platforms/opencl/tests/TestOpenCLCustomHbondForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -100,7 +100,7 @@ void testHbond() {
     angle->addAngle(1, 2, 3, 1.9, 0.6);
     standardSystem.addForce(angle);
     PeriodicTorsionForce* torsion = new PeriodicTorsionForce();
-    torsion->addTorsion(1, 2, 3, 4, 2, 2.1, 0.7);;
+    torsion->addTorsion(1, 2, 3, 4, 2, 2.1, 0.7);
     standardSystem.addForce(torsion);
 
     // Set the atoms in various positions, and verify that both systems give identical forces and energy.
@@ -109,11 +109,11 @@ void testHbond() {
     init_gen_rand(0, sfmt);
 
     vector<Vec3> positions(5);
-    for (int i = 0; i < 10; i++) {
     VerletIntegrator integrator1(0.01);
     VerletIntegrator integrator2(0.01);
     Context c1(customSystem, integrator1, platform);
     Context c2(standardSystem, integrator2, platform);
+    for (int i = 0; i < 10; i++) {
         for (int j = 0; j < (int) positions.size(); j++)
             positions[j] = Vec3(2.0*genrand_real2(sfmt), 2.0*genrand_real2(sfmt), 2.0*genrand_real2(sfmt));
         c1.setPositions(positions);
@@ -124,6 +124,28 @@ void testHbond() {
             ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], TOL);
         ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), TOL);
     }
+    
+    // Try changing the parameters and make sure it's still correct.
+    
+    parameters.resize(3);
+    parameters[0] = 1.4;
+    parameters[1] = 1.7;
+    parameters[2] = 1.9;
+    custom->setDonorParameters(0, 1, 0, -1, parameters);
+    parameters.resize(2);
+    parameters[0] = 2.2;
+    parameters[1] = 2;
+    custom->setAcceptorParameters(0, 2, 3, 4, parameters);
+    bond->setBondParameters(0, 1, 2, 1.4, 0.4);
+    torsion->setTorsionParameters(0, 1, 2, 3, 4, 2, 2.2, 0.7);
+    custom->updateParametersInContext(c1);
+    bond->updateParametersInContext(c2);
+    torsion->updateParametersInContext(c2);
+    State s1 = c1.getState(State::Forces | State::Energy);
+    State s2 = c2.getState(State::Forces | State::Energy);
+    for (int i = 0; i < customSystem.getNumParticles(); i++)
+        ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], TOL);
+    ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), TOL);
 }
 
 void testExclusions() {

platforms/opencl/tests/TestOpenCLCustomNonbondedForce.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -105,6 +105,9 @@ void testParameters() {
     ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
     ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
     ASSERT_EQUAL_TOL(3.0*(10*10*10), state.getPotentialEnergy(), TOL);
+    
+    // Try changing the global parameters and make sure it's still correct.
+    
     context.setParameter("scale", 1.5);
     context.setParameter("c", 1.0);
     state = context.getState(State::Forces | State::Energy);
@@ -113,6 +116,22 @@ void testParameters() {
     ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
     ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
     ASSERT_EQUAL_TOL(1.5*3.0*(12*12*12), state.getPotentialEnergy(), TOL);
+    
+    // Try changing the per-particle parameters and make sure it's still correct.
+    
+    params[0] = 1.6;
+    params[1] = 2.1;
+    forceField->setParticleParameters(0, params);
+    params[0] = 1.9;
+    params[1] = 2.8;
+    forceField->setParticleParameters(1, params);
+    forceField->updateParametersInContext(context);
+    state = context.getState(State::Forces | State::Energy);
+    forces = state.getForces();
+    force = -1.5*1.6*1.9*3*5.9*(11.8*11.8);
+    ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+    ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[1], TOL);
+    ASSERT_EQUAL_TOL(1.5*1.6*1.9*(11.8*11.8*11.8), state.getPotentialEnergy(), TOL);
 }
 
 void testManyParameters() {

platforms/reference/src/ReferenceKernels.cpp

@@ -1041,6 +1041,22 @@ double ReferenceCalcCustomNonbondedForceKernel::execute(ContextImpl& context, bo
     return energy;
 }
 
+void ReferenceCalcCustomNonbondedForceKernel::copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force) {
+    if (numParticles != force.getNumParticles())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerParticleParameters();
+    vector<double> params;
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+}
+
 ReferenceCalcGBSAOBCForceKernel::~ReferenceCalcGBSAOBCForceKernel() {
     if (obc) {
         // delete obc->getObcParameters();
@@ -1303,6 +1319,22 @@ double ReferenceCalcCustomGBForceKernel::execute(ContextImpl& context, bool incl
     return energy;
 }
 
+void ReferenceCalcCustomGBForceKernel::copyParametersToContext(ContextImpl& context, const CustomGBForce& force) {
+    if (numParticles != force.getNumParticles())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerParticleParameters();
+    vector<double> params;
+    for (int i = 0; i < numParticles; ++i) {
+        vector<double> parameters;
+        force.getParticleParameters(i, parameters);
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+}
+
 ReferenceCalcCustomExternalForceKernel::~ReferenceCalcCustomExternalForceKernel() {
     disposeRealArray(particleParamArray, numParticles);
 }
@@ -1348,6 +1380,25 @@ double ReferenceCalcCustomExternalForceKernel::execute(ContextImpl& context, boo
     return energy;
 }
 
+void ReferenceCalcCustomExternalForceKernel::copyParametersToContext(ContextImpl& context, const CustomExternalForce& force) {
+    if (numParticles != force.getNumParticles())
+        throw OpenMMException("updateParametersInContext: The number of particles has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerParticleParameters();
+    vector<double> params;
+    for (int i = 0; i < numParticles; ++i) {
+        int particle;
+        vector<double> parameters;
+        force.getParticleParameters(i, particle, parameters);
+        if (particle != particles[i])
+            throw OpenMMException("updateParametersInContext: A particle index has changed");
+        for (int j = 0; j < numParameters; j++)
+            particleParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+}
+
 ReferenceCalcCustomHbondForceKernel::~ReferenceCalcCustomHbondForceKernel() {
     disposeRealArray(donorParamArray, numDonors);
     disposeRealArray(acceptorParamArray, numAcceptors);
@@ -1458,6 +1509,37 @@ double ReferenceCalcCustomHbondForceKernel::execute(ContextImpl& context, bool i
     return energy;
 }
 
+void ReferenceCalcCustomHbondForceKernel::copyParametersToContext(ContextImpl& context, const CustomHbondForce& force) {
+    if (numDonors != force.getNumDonors())
+        throw OpenMMException("updateParametersInContext: The number of donors has changed");
+    if (numAcceptors != force.getNumAcceptors())
+        throw OpenMMException("updateParametersInContext: The number of acceptors has changed");
+
+    // Record the values.
+
+    vector<double> parameters;
+    int numDonorParameters = force.getNumPerDonorParameters();
+    const vector<vector<int> >& donorAtoms = ixn->getDonorAtoms();
+    for (int i = 0; i < numDonors; ++i) {
+        int d1, d2, d3;
+        force.getDonorParameters(i, d1, d2, d3, parameters);
+        if (d1 != donorAtoms[i][0] || d2 != donorAtoms[i][1] || d3 != donorAtoms[i][2])
+            throw OpenMMException("updateParametersInContext: The set of particles in a donor group has changed");
+        for (int j = 0; j < numDonorParameters; j++)
+            donorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+    int numAcceptorParameters = force.getNumPerAcceptorParameters();
+    const vector<vector<int> >& acceptorAtoms = ixn->getAcceptorAtoms();
+    for (int i = 0; i < numAcceptors; ++i) {
+        int a1, a2, a3;
+        force.getAcceptorParameters(i, a1, a2, a3, parameters);
+        if (a1 != acceptorAtoms[i][0] || a2 != acceptorAtoms[i][1] || a3 != acceptorAtoms[i][2])
+            throw OpenMMException("updateParametersInContext: The set of particles in an acceptor group has changed");
+        for (int j = 0; j < numAcceptorParameters; j++)
+            acceptorParamArray[i][j] = static_cast<RealOpenMM>(parameters[j]);
+    }
+}
+
 ReferenceCalcCustomCompoundBondForceKernel::~ReferenceCalcCustomCompoundBondForceKernel() {
     disposeRealArray(bondParamArray, numBonds);
     if (ixn != NULL)
@@ -1521,6 +1603,26 @@ double ReferenceCalcCustomCompoundBondForceKernel::execute(ContextImpl& context,
     return energy;
 }
 
+void ReferenceCalcCustomCompoundBondForceKernel::copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force) {
+    if (numBonds != force.getNumBonds())
+        throw OpenMMException("updateParametersInContext: The number of bonds has changed");
+
+    // Record the values.
+
+    int numParameters = force.getNumPerBondParameters();
+    const vector<vector<int> >& bondAtoms = ixn->getBondAtoms();
+    vector<int> particles;
+    vector<double> params;
+    for (int i = 0; i < numBonds; ++i) {
+        force.getBondParameters(i, particles, params);
+        for (int j = 0; j < numParticles; j++)
+            if (particles[j] != bondAtoms[i][j])
+                throw OpenMMException("updateParametersInContext: The set of particles in a bond has changed");
+        for (int j = 0; j < numParameters; j++)
+            bondParamArray[i][j] = (RealOpenMM) params[j];
+    }
+}
+
 ReferenceIntegrateVerletStepKernel::~ReferenceIntegrateVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferenceKernels.h

@@ -598,6 +598,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomNonbondedForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomNonbondedForce& force);
 private:
     int numParticles;
     int **exclusionArray;
@@ -702,6 +709,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomGBForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomGBForce& force);
 private:
     int numParticles;
     bool isPeriodic;
@@ -745,6 +759,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomExternalForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomExternalForce& force);
 private:
     int numParticles;
     std::vector<int> particles;
@@ -777,6 +798,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomHbondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomHbondForce& force);
 private:
     int numDonors, numAcceptors, numParticles;
     bool isPeriodic;
@@ -812,6 +840,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CustomCompoundBondForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CustomCompoundBondForce& force);
 private:
     int numBonds, numParticles;
     RealOpenMM **bondParamArray;

platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.h

@@ -91,6 +91,15 @@ class ReferenceCustomCompoundBondIxn : public ReferenceBondIxn {
 
        ~ReferenceCustomCompoundBondIxn();
 
+      /**---------------------------------------------------------------------------------------
+
+         Get the list atoms in each bond.
+
+         --------------------------------------------------------------------------------------- */
+
+       const std::vector<std::vector<int> >& getBondAtoms() const {
+           return bondAtoms;
+       }
 
       /**---------------------------------------------------------------------------------------
 

platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.h

@@ -116,6 +116,26 @@ class ReferenceCustomHbondIxn : public ReferenceBondIxn {
 
       void setPeriodic(OpenMM::RealVec& boxSize);
 
+      /**---------------------------------------------------------------------------------------
+
+         Get the list of atoms for each donor group.
+
+         --------------------------------------------------------------------------------------- */
+
+      const std::vector<std::vector<int> >& getDonorAtoms() const {
+          return donorAtoms;
+      }
+
+      /**---------------------------------------------------------------------------------------
+
+         Get the list of atoms for each acceptor group.
+
+         --------------------------------------------------------------------------------------- */
+
+      const std::vector<std::vector<int> >& getAcceptorAtoms() const {
+          return acceptorAtoms;
+      }
+
       /**---------------------------------------------------------------------------------------
 
          Calculate custom hbond interaction

platforms/reference/tests/TestReferenceCustomCompoundBondForce.cpp

@@ -121,6 +121,24 @@ void testBond() {
             ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
         ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
     }
+    
+    // Try changing the bond parameters and make sure it's still correct.
+    
+    parameters[0] = 1.6;
+    parameters[3] = 1.3;
+    custom->setBondParameters(0, particles, parameters);
+    custom->updateParametersInContext(c1);
+    bonds->setBondParameters(0, 0, 1, 1.3, 1.6);
+    bonds->setBondParameters(1, 1, 3, 1.3, 1.6);
+    bonds->updateParametersInContext(c2);
+    {
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = s1.getForces();
+        for (int i = 0; i < customSystem.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], TOL);
+        ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
+    }
 }
 
 void testPositionDependence() {

platforms/reference/tests/TestReferenceCustomExternalForce.cpp

@@ -74,11 +74,28 @@ void testForce() {
     positions[2] = Vec3(1, 0, 1);
     context.setPositions(positions);
     State state = context.getState(State::Forces | State::Energy);
+    {
         const vector<Vec3>& forces = state.getForces();
         ASSERT_EQUAL_VEC(Vec3(-0.5, -0.5*2.0*2.0*1.5, 0), forces[0], TOL);
         ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
         ASSERT_EQUAL_VEC(Vec3(-0.5, 0.5*3.0*2.0*1.5, 0), forces[2], TOL);
         ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.0*1.5*1.5), state.getPotentialEnergy(), TOL);
+    }
+    
+    // Try changing the parameters and make sure it's still correct.
+    
+    parameters[0] = 1.4;
+    parameters[1] = 3.5;
+    forceField->setParticleParameters(1, 2, parameters);
+    forceField->updateParametersInContext(context);
+    state = context.getState(State::Forces | State::Energy);
+    {
+        const vector<Vec3>& forces = state.getForces();
+        ASSERT_EQUAL_VEC(Vec3(-0.5, -0.5*2.0*2.0*1.5, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(0, 0, 0), forces[1], TOL);
+        ASSERT_EQUAL_VEC(Vec3(-0.5, 0.5*3.5*2.0*1.4, 0), forces[2], TOL);
+        ASSERT_EQUAL_TOL(0.5*(1.0 + 2.0*1.5*1.5 + 3.5*1.4*1.4), state.getPotentialEnergy(), TOL);
+    }
 }
 
 int main() {

platforms/reference/tests/TestReferenceCustomGBForce.cpp

@@ -7,7 +7,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -142,6 +142,29 @@ void testOBC(GBSAOBCForce::NonbondedMethod obcMethod, CustomGBForce::NonbondedMe
     for (int i = 0; i < numParticles; i++) {
         ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
     }
+    
+    // Try changing the particle parameters and make sure it's still correct.
+    
+    for (int i = 0; i < numMolecules/2; i++) {
+        obc->setParticleParameters(2*i, 1.1, 0.3, 0.6);
+        params[0] = 1.1;
+        params[1] = 0.3;
+        params[2] = 0.6;
+        custom->setParticleParameters(2*i, params);
+        obc->setParticleParameters(2*i+1, -1.1, 0.2, 0.4);
+        params[0] = -1.1;
+        params[1] = 0.2;
+        params[2] = 0.4;
+        custom->setParticleParameters(2*i+1, params);
+    }
+    obc->updateParametersInContext(context1);
+    custom->updateParametersInContext(context2);
+    state1 = context1.getState(State::Forces | State::Energy);
+    state2 = context2.getState(State::Forces | State::Energy);
+    ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
+    for (int i = 0; i < numParticles; i++) {
+        ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
+    }
 }
 
 void testMembrane() {

platforms/reference/tests/TestReferenceCustomHbondForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -100,7 +100,7 @@ void testHbond() {
     angle->addAngle(1, 2, 3, 1.9, 0.6);
     standardSystem.addForce(angle);
     PeriodicTorsionForce* torsion = new PeriodicTorsionForce();
-    torsion->addTorsion(1, 2, 3, 4, 2, 2.1, 0.7);;
+    torsion->addTorsion(1, 2, 3, 4, 2, 2.1, 0.7);
     standardSystem.addForce(torsion);
 
     // Set the atoms in various positions, and verify that both systems give identical forces and energy.
@@ -109,11 +109,11 @@ void testHbond() {
     init_gen_rand(0, sfmt);
 
     vector<Vec3> positions(5);
-    for (int i = 0; i < 10; i++) {
     VerletIntegrator integrator1(0.01);
     VerletIntegrator integrator2(0.01);
     Context c1(customSystem, integrator1, platform);
     Context c2(standardSystem, integrator2, platform);
+    for (int i = 0; i < 10; i++) {
         for (int j = 0; j < (int) positions.size(); j++)
             positions[j] = Vec3(2.0*genrand_real2(sfmt), 2.0*genrand_real2(sfmt), 2.0*genrand_real2(sfmt));
         c1.setPositions(positions);
@@ -124,6 +124,28 @@ void testHbond() {
             ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], TOL);
         ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), TOL);
     }
+    
+    // Try changing the parameters and make sure it's still correct.
+    
+    parameters.resize(3);
+    parameters[0] = 1.4;
+    parameters[1] = 1.7;
+    parameters[2] = 1.9;
+    custom->setDonorParameters(0, 1, 0, -1, parameters);
+    parameters.resize(2);
+    parameters[0] = 2.2;
+    parameters[1] = 2;
+    custom->setAcceptorParameters(0, 2, 3, 4, parameters);
+    bond->setBondParameters(0, 1, 2, 1.4, 0.4);
+    torsion->setTorsionParameters(0, 1, 2, 3, 4, 2, 2.2, 0.7);
+    custom->updateParametersInContext(c1);
+    bond->updateParametersInContext(c2);
+    torsion->updateParametersInContext(c2);
+    State s1 = c1.getState(State::Forces | State::Energy);
+    State s2 = c2.getState(State::Forces | State::Energy);
+    for (int i = 0; i < customSystem.getNumParticles(); i++)
+        ASSERT_EQUAL_VEC(s2.getForces()[i], s1.getForces()[i], TOL);
+    ASSERT_EQUAL_TOL(s2.getPotentialEnergy(), s1.getPotentialEnergy(), TOL);
 }
 
 void testExclusions() {