Merge pull request #854 from peastman/cmap

CMAPTorsionForce implements updateParametersInContext()

olla/include/openmm/kernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -494,6 +494,13 @@ public:
      * @return the potential energy due to the force
      */
     virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
+    /**
+     * Copy changed parameters over to a context.
+     *
+     * @param context    the context to copy parameters to
+     * @param force      the CMAPTorsionForce to copy the parameters from
+     */
+    virtual void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) = 0;
 };
 
 /**

openmmapi/include/openmm/CMAPTorsionForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -148,6 +148,17 @@ public:
      * @param b4    the index of the fourth particle forming the second torsion
      */
     void setTorsionParameters(int index, int map, int a1, int a2, int a3, int a4, int b1, int b2, int b3, int b4);
+    /**
+     * Update the map and torsion parameters in a Context to match those stored in this Force object.  This method provides
+     * an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
+     * Simply call setMapParameters() and setTorsionParameters() to modify this object's parameters, then call updateParametersInContext()
+     * to copy them over to the Context.
+     *
+     * The only information that can be updated with this method is the energy values for a map, and the map index
+     * for a torsion.  The size of a map and the set of particles involved in a torsion cannot be changed.  Also,
+     * new bonds and torsions cannot be added.
+     */
+    void updateParametersInContext(Context& context);
     /**
      * Returns whether or not this force makes use of periodic boundary
      * conditions.

openmmapi/include/openmm/internal/CMAPTorsionForceImpl.h

@@ -9,7 +9,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-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -60,6 +60,7 @@ public:
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }
     std::vector<std::string> getKernelNames();
+    void updateParametersInContext(ContextImpl& context);
     /**
      * Given the energy values for a map, compute the spline coefficients at each point of the map.
      */

openmmapi/src/CMAPTorsionForce.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) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -96,3 +96,6 @@ ForceImpl* CMAPTorsionForce::createImpl() const {
     return new CMAPTorsionForceImpl(*this);
 }
 
+void CMAPTorsionForce::updateParametersInContext(Context& context) {
+    dynamic_cast<CMAPTorsionForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
+}

openmmapi/src/CMAPTorsionForceImpl.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) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -154,4 +154,8 @@ void CMAPTorsionForceImpl::calcMapDerivatives(int size, const vector<double>& en
             }
         }
     }
-}
\ No newline at end of file
+}
+
+void CMAPTorsionForceImpl::updateParametersInContext(ContextImpl& context) {
+    kernel.getAs<CalcCMAPTorsionForceKernel>().copyParametersToContext(context, owner);
+}

platforms/cuda/include/CudaKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -520,11 +520,19 @@ 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 CMAPTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
 private:
     int numTorsions;
     bool hasInitializedKernel;
     CudaContext& cu;
     const System& system;
+    std::vector<int2> mapPositionsVec;
     CudaArray* coefficients;
     CudaArray* mapPositions;
     CudaArray* torsionMaps;

platforms/cuda/include/CudaParallelKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -345,6 +345,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 CMAPTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
 private:
     class Task;
     CudaPlatform::PlatformData& data;

platforms/cuda/src/CudaKernels.cpp

@@ -1127,7 +1127,7 @@ void CudaCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAP
         return;
     int numMaps = force.getNumMaps();
     vector<float4> coeffVec;
-    vector<int2> mapPositionsVec(numMaps);
+    mapPositionsVec.resize(numMaps);
     vector<double> energy;
     vector<vector<double> > c;
     int currentPosition = 0;
@@ -1166,6 +1166,49 @@ double CudaCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool includ
     return 0.0;
 }
 
+void CudaCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) {
+    int numMaps = force.getNumMaps();
+    int numContexts = cu.getPlatformData().contexts.size();
+    int startIndex = cu.getContextIndex()*force.getNumTorsions()/numContexts;
+    int endIndex = (cu.getContextIndex()+1)*force.getNumTorsions()/numContexts;
+    numTorsions = endIndex-startIndex;
+    if (mapPositions->getSize() != numMaps)
+        throw OpenMMException("updateParametersInContext: The number of maps has changed");
+    if (torsionMaps->getSize() != numTorsions)
+        throw OpenMMException("updateParametersInContext: The number of CMAP torsions has changed");
+
+    // Update the maps.
+
+    vector<float4> coeffVec;
+    vector<double> energy;
+    vector<vector<double> > c;
+    int currentPosition = 0;
+    for (int i = 0; i < numMaps; i++) {
+        int size;
+        force.getMapParameters(i, size, energy);
+        if (size != mapPositionsVec[i].y)
+            throw OpenMMException("updateParametersInContext: The size of a map has changed");
+        CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
+        currentPosition += 4*size*size;
+        for (int j = 0; j < size*size; j++) {
+            coeffVec.push_back(make_float4((float) c[j][0], (float) c[j][1], (float) c[j][2], (float) c[j][3]));
+            coeffVec.push_back(make_float4((float) c[j][4], (float) c[j][5], (float) c[j][6], (float) c[j][7]));
+            coeffVec.push_back(make_float4((float) c[j][8], (float) c[j][9], (float) c[j][10], (float) c[j][11]));
+            coeffVec.push_back(make_float4((float) c[j][12], (float) c[j][13], (float) c[j][14], (float) c[j][15]));
+        }
+    }
+    coefficients->upload(coeffVec);
+
+    // Update the indices.
+
+    vector<int> torsionMapsVec(numTorsions);
+    for (int i = 0; i < numTorsions; i++) {
+        int index[8];
+        force.getTorsionParameters(i, torsionMapsVec[i], index[0], index[1], index[2], index[3], index[4], index[5], index[6], index[7]);
+    }
+    torsionMaps->upload(torsionMapsVec);
+}
+
 class CudaCustomTorsionForceInfo : public CudaForceInfo {
 public:
     CudaCustomTorsionForceInfo(const CustomTorsionForce& force) : force(force) {

platforms/cuda/src/CudaParallelKernels.cpp

@@ -536,6 +536,11 @@ double CudaParallelCalcCMAPTorsionForceKernel::execute(ContextImpl& context, boo
     return 0.0;
 }
 
+void CudaParallelCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) {
+    for (int i = 0; i < (int) kernels.size(); i++)
+        getKernel(i).copyParametersToContext(context, force);
+}
+
 class CudaParallelCalcCustomTorsionForceKernel::Task : public CudaContext::WorkTask {
 public:
     Task(ContextImpl& context, CudaCalcCustomTorsionForceKernel& kernel, bool includeForce,

platforms/cuda/tests/TestCudaCMAPTorsionForce.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) 2010-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -105,11 +105,68 @@ void testCMAPTorsions() {
     }
 }
 
+void testChangingParameters() {
+    // Create a system with two maps and one torsion.
+
+    const int mapSize = 8;
+    System system;
+    for (int i = 0; i < 5; i++)
+        system.addParticle(1.0);
+    CMAPTorsionForce* cmap = new CMAPTorsionForce();
+    vector<double> mapEnergy1(mapSize*mapSize);
+    vector<double> mapEnergy2(mapSize*mapSize);
+    for (int i = 0; i < mapSize; i++) {
+        double angle1 = i*2*M_PI/mapSize;
+        double energy1 = cos(angle1);
+        for (int j = 0; j < mapSize; j++) {
+            double angle2 = j*2*M_PI/mapSize;
+            double energy2 = 10*sin(angle2);
+            mapEnergy1[i+j*mapSize] = energy1+energy2;
+            mapEnergy2[i+j*mapSize] = energy1-energy2;
+        }
+    }
+    cmap->addMap(mapSize, mapEnergy1);
+    cmap->addMap(mapSize, mapEnergy2);
+    cmap->addTorsion(0, 0, 1, 2, 3, 1, 2, 3, 4);
+    system.addForce(cmap);
+
+    // Set particle positions so angle1=0 and angle2=PI/4.
+
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(0, 0, 1);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 0, 1);
+    positions[4] = Vec3(0.5, -0.5, 1);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+
+    // Check that the energy is correct.
+
+    double energy = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(1+10*sin(M_PI/4), energy, 1e-5);
+
+    // Modify the parameters.
+
+    cmap->setTorsionParameters(0, 1, 0, 1, 2, 3, 1, 2, 3, 4);
+    for (int i = 0; i < mapSize*mapSize; i++)
+        mapEnergy2[i] *= 2.0;
+    cmap->setMapParameters(1, mapSize, mapEnergy2);
+    cmap->updateParametersInContext(context);
+
+    // See if the results are correct.
+
+    energy = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(2-20*sin(M_PI/4), energy, 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
             platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
         testCMAPTorsions();
+        testChangingParameters();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/include/OpenCLKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -498,11 +498,19 @@ 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 CMAPTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
 private:
     int numTorsions;
     bool hasInitializedKernel;
     OpenCLContext& cl;
     const System& system;
+    std::vector<mm_int2> mapPositionsVec;
     OpenCLArray* coefficients;
     OpenCLArray* mapPositions;
     OpenCLArray* torsionMaps;

platforms/opencl/include/OpenCLParallelKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -344,6 +344,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 CMAPTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
 private:
     class Task;
     OpenCLPlatform::PlatformData& data;

platforms/opencl/src/OpenCLKernels.cpp

@@ -1125,7 +1125,7 @@ void OpenCLCalcCMAPTorsionForceKernel::initialize(const System& system, const CM
         return;
     int numMaps = force.getNumMaps();
     vector<mm_float4> coeffVec;
-    vector<mm_int2> mapPositionsVec(numMaps);
+    mapPositionsVec.resize(numMaps);
     vector<double> energy;
     vector<vector<double> > c;
     int currentPosition = 0;
@@ -1164,6 +1164,49 @@ double OpenCLCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool incl
     return 0.0;
 }
 
+void OpenCLCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) {
+    int numMaps = force.getNumMaps();
+    int numContexts = cl.getPlatformData().contexts.size();
+    int startIndex = cl.getContextIndex()*force.getNumTorsions()/numContexts;
+    int endIndex = (cl.getContextIndex()+1)*force.getNumTorsions()/numContexts;
+    numTorsions = endIndex-startIndex;
+    if (mapPositions->getSize() != numMaps)
+        throw OpenMMException("updateParametersInContext: The number of maps has changed");
+    if (torsionMaps->getSize() != numTorsions)
+        throw OpenMMException("updateParametersInContext: The number of CMAP torsions has changed");
+
+    // Update the maps.
+
+    vector<mm_float4> coeffVec;
+    vector<double> energy;
+    vector<vector<double> > c;
+    int currentPosition = 0;
+    for (int i = 0; i < numMaps; i++) {
+        int size;
+        force.getMapParameters(i, size, energy);
+        if (size != mapPositionsVec[i].y)
+            throw OpenMMException("updateParametersInContext: The size of a map has changed");
+        CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
+        currentPosition += 4*size*size;
+        for (int j = 0; j < size*size; j++) {
+            coeffVec.push_back(mm_float4((float) c[j][0], (float) c[j][1], (float) c[j][2], (float) c[j][3]));
+            coeffVec.push_back(mm_float4((float) c[j][4], (float) c[j][5], (float) c[j][6], (float) c[j][7]));
+            coeffVec.push_back(mm_float4((float) c[j][8], (float) c[j][9], (float) c[j][10], (float) c[j][11]));
+            coeffVec.push_back(mm_float4((float) c[j][12], (float) c[j][13], (float) c[j][14], (float) c[j][15]));
+        }
+    }
+    coefficients->upload(coeffVec);
+
+    // Update the indices.
+
+    vector<int> torsionMapsVec(numTorsions);
+    for (int i = 0; i < numTorsions; i++) {
+        int index[8];
+        force.getTorsionParameters(i, torsionMapsVec[i], index[0], index[1], index[2], index[3], index[4], index[5], index[6], index[7]);
+    }
+    torsionMaps->upload(torsionMapsVec);
+}
+
 class OpenCLCustomTorsionForceInfo : public OpenCLForceInfo {
 public:
     OpenCLCustomTorsionForceInfo(const CustomTorsionForce& force) : OpenCLForceInfo(0), force(force) {

platforms/opencl/src/OpenCLParallelKernels.cpp

@@ -492,6 +492,11 @@ double OpenCLParallelCalcCMAPTorsionForceKernel::execute(ContextImpl& context, b
     return 0.0;
 }
 
+void OpenCLParallelCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) {
+    for (int i = 0; i < (int) kernels.size(); i++)
+        getKernel(i).copyParametersToContext(context, force);
+}
+
 class OpenCLParallelCalcCustomTorsionForceKernel::Task : public OpenCLContext::WorkTask {
 public:
     Task(ContextImpl& context, OpenCLCalcCustomTorsionForceKernel& kernel, bool includeForce,

platforms/opencl/tests/TestOpenCLCMAPTorsionForce.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) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -105,11 +105,68 @@ void testCMAPTorsions() {
     }
 }
 
+void testChangingParameters() {
+    // Create a system with two maps and one torsion.
+
+    const int mapSize = 8;
+    System system;
+    for (int i = 0; i < 5; i++)
+        system.addParticle(1.0);
+    CMAPTorsionForce* cmap = new CMAPTorsionForce();
+    vector<double> mapEnergy1(mapSize*mapSize);
+    vector<double> mapEnergy2(mapSize*mapSize);
+    for (int i = 0; i < mapSize; i++) {
+        double angle1 = i*2*M_PI/mapSize;
+        double energy1 = cos(angle1);
+        for (int j = 0; j < mapSize; j++) {
+            double angle2 = j*2*M_PI/mapSize;
+            double energy2 = 10*sin(angle2);
+            mapEnergy1[i+j*mapSize] = energy1+energy2;
+            mapEnergy2[i+j*mapSize] = energy1-energy2;
+        }
+    }
+    cmap->addMap(mapSize, mapEnergy1);
+    cmap->addMap(mapSize, mapEnergy2);
+    cmap->addTorsion(0, 0, 1, 2, 3, 1, 2, 3, 4);
+    system.addForce(cmap);
+
+    // Set particle positions so angle1=0 and angle2=PI/4.
+
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(0, 0, 1);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 0, 1);
+    positions[4] = Vec3(0.5, -0.5, 1);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+
+    // Check that the energy is correct.
+
+    double energy = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(1+10*sin(M_PI/4), energy, 1e-5);
+
+    // Modify the parameters.
+
+    cmap->setTorsionParameters(0, 1, 0, 1, 2, 3, 1, 2, 3, 4);
+    for (int i = 0; i < mapSize*mapSize; i++)
+        mapEnergy2[i] *= 2.0;
+    cmap->setMapParameters(1, mapSize, mapEnergy2);
+    cmap->updateParametersInContext(context);
+
+    // See if the results are correct.
+
+    energy = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(2-20*sin(M_PI/4), energy, 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
             platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
         testCMAPTorsions();
+        testChangingParameters();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/include/ReferenceKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -493,6 +493,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 CMAPTorsionForce to copy the parameters from
+     */
+    void copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force);
 private:
     std::vector<std::vector<std::vector<RealOpenMM> > > coeff;
     std::vector<int> torsionMaps;

platforms/reference/src/ReferenceKernels.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-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -710,6 +710,40 @@ double ReferenceCalcCMAPTorsionForceKernel::execute(ContextImpl& context, bool i
     return totalEnergy;
 }
 
+void ReferenceCalcCMAPTorsionForceKernel::copyParametersToContext(ContextImpl& context, const CMAPTorsionForce& force) {
+    int numMaps = force.getNumMaps();
+    int numTorsions = force.getNumTorsions();
+    if (coeff.size() != numMaps)
+        throw OpenMMException("updateParametersInContext: The number of maps has changed");
+    if (torsionMaps.size() != numTorsions)
+        throw OpenMMException("updateParametersInContext: The number of CMAP torsions has changed");
+
+    // Update the maps.
+
+    vector<double> energy;
+    vector<vector<double> > c;
+    for (int i = 0; i < numMaps; i++) {
+        int size;
+        force.getMapParameters(i, size, energy);
+        if (coeff[i].size() != size*size)
+            throw OpenMMException("updateParametersInContext: The size of a map has changed");
+        CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
+        for (int j = 0; j < size*size; j++)
+            for (int k = 0; k < 16; k++)
+                coeff[i][j][k] = c[j][k];
+    }
+
+    // Update the indices.
+
+    for (int i = 0; i < numTorsions; i++) {
+        int index[8];
+        force.getTorsionParameters(i, torsionMaps[i], index[0], index[1], index[2], index[3], index[4], index[5], index[6], index[7]);
+        for (int j = 0; j < 8; j++)
+            if (index[j] != torsionIndices[i][j])
+                throw OpenMMException("updateParametersInContext: The set of particles in a CMAP torsion has changed");
+    }
+}
+
 ReferenceCalcCustomTorsionForceKernel::~ReferenceCalcCustomTorsionForceKernel() {
     disposeIntArray(torsionIndexArray, numTorsions);
     disposeRealArray(torsionParamArray, numTorsions);

platforms/reference/tests/TestReferenceCMAPTorsionForce.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) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,6 +48,8 @@
 using namespace OpenMM;
 using namespace std;
 
+ReferencePlatform platform;
+
 const double TOL = 1e-5;
 
 void testCMAPTorsions() {
@@ -56,7 +58,6 @@ void testCMAPTorsions() {
     // Create two systems: one with a pair of periodic torsions, and one with a CMAP torsion
     // that approximates the same force.
 
-    ReferencePlatform platform;
     System system1;
     for (int i = 0; i < 5; i++)
         system1.addParticle(1.0);
@@ -108,9 +109,66 @@ void testCMAPTorsions() {
     }
 }
 
+void testChangingParameters() {
+    // Create a system with two maps and one torsion.
+
+    const int mapSize = 8;
+    System system;
+    for (int i = 0; i < 5; i++)
+        system.addParticle(1.0);
+    CMAPTorsionForce* cmap = new CMAPTorsionForce();
+    vector<double> mapEnergy1(mapSize*mapSize);
+    vector<double> mapEnergy2(mapSize*mapSize);
+    for (int i = 0; i < mapSize; i++) {
+        double angle1 = i*2*M_PI/mapSize;
+        double energy1 = cos(angle1);
+        for (int j = 0; j < mapSize; j++) {
+            double angle2 = j*2*M_PI/mapSize;
+            double energy2 = 10*sin(angle2);
+            mapEnergy1[i+j*mapSize] = energy1+energy2;
+            mapEnergy2[i+j*mapSize] = energy1-energy2;
+        }
+    }
+    cmap->addMap(mapSize, mapEnergy1);
+    cmap->addMap(mapSize, mapEnergy2);
+    cmap->addTorsion(0, 0, 1, 2, 3, 1, 2, 3, 4);
+    system.addForce(cmap);
+
+    // Set particle positions so angle1=0 and angle2=PI/4.
+
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(0, 0, 1);
+    positions[1] = Vec3(0, 0, 0);
+    positions[2] = Vec3(1, 0, 0);
+    positions[3] = Vec3(1, 0, 1);
+    positions[4] = Vec3(0.5, -0.5, 1);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+
+    // Check that the energy is correct.
+
+    double energy = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(1+10*sin(M_PI/4), energy, 1e-5);
+
+    // Modify the parameters.
+
+    cmap->setTorsionParameters(0, 1, 0, 1, 2, 3, 1, 2, 3, 4);
+    for (int i = 0; i < mapSize*mapSize; i++)
+        mapEnergy2[i] *= 2.0;
+    cmap->setMapParameters(1, mapSize, mapEnergy2);
+    cmap->updateParametersInContext(context);
+
+    // See if the results are correct.
+
+    energy = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(2-20*sin(M_PI/4), energy, 1e-5);
+}
+
 int main() {
     try {
         testCMAPTorsions();
+        testChangingParameters();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;