Implemented CMAPTorsionForce, including reference implementation

olla/include/openmm/kernels.h

@@ -34,6 +34,7 @@
 
 #include "openmm/AndersenThermostat.h"
 #include "openmm/BrownianIntegrator.h"
+#include "openmm/CMAPTorsionForce.h"
 #include "openmm/CMMotionRemover.h"
 #include "openmm/CustomAngleForce.h"
 #include "openmm/CustomBondForce.h"
@@ -408,6 +409,38 @@ public:
     virtual double executeEnergy(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by CMAPTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCMAPTorsionForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcCMAPTorsionForce";
+    }
+    CalcCMAPTorsionForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CMAPTorsionForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const CMAPTorsionForce& force) = 0;
+    /**
+     * Execute the kernel to calculate the forces.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    virtual void executeForces(ContextImpl& context) = 0;
+    /**
+     * Execute the kernel to calculate the energy.
+     *
+     * @param context    the context in which to execute this kernel
+     * @return the potential energy due to the CMAPTorsionForce
+     */
+    virtual double executeEnergy(ContextImpl& context) = 0;
+};
+
 /**
  * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
  */

openmmapi/include/OpenMM.h

@@ -34,6 +34,7 @@
 
 #include "openmm/AndersenThermostat.h"
 #include "openmm/BrownianIntegrator.h"
+#include "openmm/CMAPTorsionForce.h"
 #include "openmm/CMMotionRemover.h"
 #include "openmm/CustomBondForce.h"
 #include "openmm/CustomAngleForce.h"

openmmapi/include/openmm/CMAPTorsionForce.h

+#ifndef OPENMM_CMAPTORSIONFORCE_H_
+#define OPENMM_CMAPTORSIONFORCE_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) 2010 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "Force.h"
+#include "Vec3.h"
+#include <vector>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class implements an interaction between pairs of dihedral angles.  The interaction energy is
+ * defined by an "energy correction map" (CMAP), which is simply a set of tabulated energy values
+ * on a regular grid of (phi, psi) angles.  Natural cubic spline interpolation is used to compute
+ * forces and energies at arbitrary values of the two angles.
+ *
+ * To use this class, first create one or more energy correction maps by calling addMap().  For each
+ * one, you provide an array of energies at uniformly spaced values of the two angles.  Next,
+ * add interactions by calling addTorsion().  For each one, you specify the sequence of particles used
+ * to calculate each of the two dihedral angles, and the index of the map used to calculate their
+ * interaction energy.
+ */
+
+class OPENMM_EXPORT CMAPTorsionForce : public Force {
+public:
+    /**
+     * Create a CMAPTorsionForce.
+     */
+    CMAPTorsionForce();
+    /**
+     * Get the number of maps that have been defined.
+     */
+    int getNumMaps() const {
+        return maps.size();
+    }
+    /**
+     * Get the number of CMAP torsion terms in the potential function
+     */
+    int getNumTorsions() const {
+        return torsions.size();
+    }
+    /**
+     * Create a new map that can be used for torsion pairs.
+     *
+     * @param size    the size of the map along each dimension
+     * @param energy  the energy values for the map.  This must be of length size*size.
+     *                The element energy[i+size*j] contains the energy when the first
+     *                torsion angle equals i*2*PI/size and the second torsion angle
+     *                equals j*2*PI/size.
+     * @return the index of the map that was added
+     */
+    int addMap(int size, const std::vector<double>& energy);
+    /**
+     * Get the energy values of a map.
+     *
+     * @param index   the index of the map for which to get energy values
+     * @param size    the size of the map along each dimension
+     * @param energy  the energy values for the map.  This must be of length size*size.
+     *                The element energy[i+size*j] contains the energy when the first
+     *                torsion angle equals i*2*PI/size and the second torsion angle
+     *                equals j*2*PI/size.
+     */
+    void getMapParameters(int index, int& size, std::vector<double>& energy) const;
+    /**
+     * Set the energy values of a map.
+     *
+     * @param index   the index of the map for which to set energy values
+     * @param size    the size of the map along each dimension
+     * @param energy  the energy values for the map.  This must be of length size*size.
+     *                The element energy[i+size*j] contains the energy when the first
+     *                torsion angle equals i*2*PI/size and the second torsion angle
+     *                equals j*2*PI/size.
+     */
+    void setMapParameters(int index, int size, const std::vector<double>& energy);
+    /**
+     * Add a CMAP torsion term to the force field.
+     *
+     * @param map   the index of the map to use for this term
+     * @param a1    the index of the first particle forming the first torsion
+     * @param a2    the index of the second particle forming the first torsion
+     * @param a3    the index of the third particle forming the first torsion
+     * @param a4    the index of the fourth particle forming the first torsion
+     * @param b1    the index of the first particle forming the second torsion
+     * @param b2    the index of the second particle forming the second torsion
+     * @param b3    the index of the third particle forming the second torsion
+     * @param b4    the index of the fourth particle forming the second torsion
+     * @return the index of the torsion that was added
+     */
+    int addTorsion(int map, int a1, int a2, int a3, int a4, int b1, int b2, int b3, int b4);
+    /**
+     * Get the force field parameters for a CMAP torsion term.
+     *
+     * @param index the index of the torsion for which to get parameters
+     * @param map   the index of the map to use for this term
+     * @param a1    the index of the first particle forming the first torsion
+     * @param a2    the index of the second particle forming the first torsion
+     * @param a3    the index of the third particle forming the first torsion
+     * @param a4    the index of the fourth particle forming the first torsion
+     * @param b1    the index of the first particle forming the second torsion
+     * @param b2    the index of the second particle forming the second torsion
+     * @param b3    the index of the third particle forming the second torsion
+     * @param b4    the index of the fourth particle forming the second torsion
+     */
+    void getTorsionParameters(int index, int& map, int& a1, int& a2, int& a3, int& a4, int& b1, int& b2, int& b3, int& b4) const;
+    /**
+     * Set the force field parameters for a CMAP torsion term.
+     *
+     * @param index the index of the torsion for which to set parameters
+     * @param map   the index of the map to use for this term
+     * @param a1    the index of the first particle forming the first torsion
+     * @param a2    the index of the second particle forming the first torsion
+     * @param a3    the index of the third particle forming the first torsion
+     * @param a4    the index of the fourth particle forming the first torsion
+     * @param b1    the index of the first particle forming the second torsion
+     * @param b2    the index of the second particle forming the second torsion
+     * @param b3    the index of the third particle forming the second torsion
+     * @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);
+protected:
+    ForceImpl* createImpl();
+private:
+    class MapInfo;
+    class CMAPTorsionInfo;
+    std::vector<MapInfo> maps;
+    std::vector<CMAPTorsionInfo> torsions;
+};
+
+/**
+ * This is an internal class used to record information about a map.
+ * @private
+ */
+class CMAPTorsionForce::MapInfo {
+public:
+    int size;
+    std::vector<double> energy;
+    MapInfo() {
+        size = -1;
+    }
+    MapInfo(int size, const std::vector<double>& energy) :
+        size(size), energy(energy) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a torsion.
+ * @private
+ */
+class CMAPTorsionForce::CMAPTorsionInfo {
+public:
+    int map, a1, a2, a3, a4, b1, b2, b3, b4;
+    CMAPTorsionInfo() {
+        map = a1 = a2 = a3 = a4 = b1 = b2 = b3 = b4 = -1;
+    }
+    CMAPTorsionInfo(int map, int a1, int a2, int a3, int a4, int b1, int b2, int b3, int b4) :
+        map(map), a1(a1), a2(a2), a3(a3), a4(a4), b1(b1), b2(b2), b3(b3), b4(b4) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CMAPTORSIONFORCE_H_*/

openmmapi/include/openmm/internal/CMAPTorsionForceImpl.h

+#ifndef OPENMM_CMAPTORSIONFORCEIMPL_H_
+#define OPENMM_CMAPTORSIONFORCEIMPL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "ForceImpl.h"
+#include "openmm/CMAPTorsionForce.h"
+#include "openmm/Kernel.h"
+#include <map>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of CMAPTorsionForce.
+ */
+
+class CMAPTorsionForceImpl : public ForceImpl {
+public:
+    CMAPTorsionForceImpl(CMAPTorsionForce& owner);
+    ~CMAPTorsionForceImpl();
+    void initialize(ContextImpl& context);
+    CMAPTorsionForce& getOwner() {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context) {
+        // This force field doesn't update the state directly.
+    }
+    void calcForces(ContextImpl& context);
+    double calcEnergy(ContextImpl& context);
+    std::map<std::string, double> getDefaultParameters() {
+        return std::map<std::string, double>(); // This force field doesn't define any parameters.
+    }
+    std::vector<std::string> getKernelNames();
+    /**
+     * Given the energy values for a map, compute the spline coefficients at each point of the map.
+     */
+    static void calcMapDerivatives(int size, const std::vector<double>& energy, std::vector<std::vector<double> >& c);
+private:
+    CMAPTorsionForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CMAPTORSIONFORCEIMPL_H_*/

openmmapi/src/CMAPTorsionForce.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) 2010 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/Force.h"
+#include "openmm/OpenMMException.h"
+#include "openmm/CMAPTorsionForce.h"
+#include "openmm/internal/CMAPTorsionForceImpl.h"
+
+using namespace OpenMM;
+
+CMAPTorsionForce::CMAPTorsionForce() {
+}
+
+int CMAPTorsionForce::addMap(int size, const std::vector<double>& energy) {
+    if (energy.size() != size*size)
+        throw OpenMMException("CMAPTorsionForce: incorrect number of energy values");
+    maps.push_back(MapInfo(size, energy));
+    return maps.size()-1;
+}
+
+void CMAPTorsionForce::getMapParameters(int index, int& size, std::vector<double>& energy) const {
+    size = maps[index].size;
+    energy = maps[index].energy;
+}
+
+void CMAPTorsionForce::setMapParameters(int index, int size, const std::vector<double>& energy) {
+    if (energy.size() != size*size)
+        throw OpenMMException("CMAPTorsionForce: incorrect number of energy values");
+    maps[index].size = size;
+    maps[index].energy = energy;
+}
+
+int CMAPTorsionForce::addTorsion(int map, int a1, int a2, int a3, int a4, int b1, int b2, int b3, int b4) {
+    torsions.push_back(CMAPTorsionInfo(map, a1, a2, a3, a4, b1, b2, b3, b4));
+    return torsions.size()-1;
+}
+
+void CMAPTorsionForce::getTorsionParameters(int index, int& map, int& a1, int& a2, int& a3, int& a4, int& b1, int& b2, int& b3, int& b4) const {
+    map = torsions[index].map;
+    a1 = torsions[index].a1;
+    a2 = torsions[index].a2;
+    a3 = torsions[index].a3;
+    a4 = torsions[index].a4;
+    b1 = torsions[index].b1;
+    b2 = torsions[index].b2;
+    b3 = torsions[index].b3;
+    b4 = torsions[index].b4;
+}
+
+void CMAPTorsionForce::setTorsionParameters(int index, int map, int a1, int a2, int a3, int a4, int b1, int b2, int b3, int b4) {
+    torsions[index].map = map;
+    torsions[index].a1 = a1;
+    torsions[index].a2 = a2;
+    torsions[index].a3 = a3;
+    torsions[index].a4 = a4;
+    torsions[index].b1 = b1;
+    torsions[index].b2 = b2;
+    torsions[index].b3 = b3;
+    torsions[index].b4 = b4;
+}
+
+ForceImpl* CMAPTorsionForce::createImpl() {
+    return new CMAPTorsionForceImpl(*this);
+}
+

openmmapi/src/CMAPTorsionForceImpl.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) 2010 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#ifdef WIN32
+  #define _USE_MATH_DEFINES // Needed to get M_PI
+#endif
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/CMAPTorsionForceImpl.h"
+#include "openmm/internal/SplineFitter.h"
+#include "openmm/kernels.h"
+#include <cmath>
+
+using namespace OpenMM;
+using namespace std;
+
+CMAPTorsionForceImpl::CMAPTorsionForceImpl(CMAPTorsionForce& owner) : owner(owner) {
+}
+
+CMAPTorsionForceImpl::~CMAPTorsionForceImpl() {
+}
+
+void CMAPTorsionForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcCMAPTorsionForceKernel::Name(), context);
+    dynamic_cast<CalcCMAPTorsionForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
+}
+
+void CMAPTorsionForceImpl::calcForces(ContextImpl& context) {
+    dynamic_cast<CalcCMAPTorsionForceKernel&>(kernel.getImpl()).executeForces(context);
+}
+
+double CMAPTorsionForceImpl::calcEnergy(ContextImpl& context) {
+    return dynamic_cast<CalcCMAPTorsionForceKernel&>(kernel.getImpl()).executeEnergy(context);
+}
+
+vector<string> CMAPTorsionForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCMAPTorsionForceKernel::Name());
+    return names;
+}
+
+void CMAPTorsionForceImpl::calcMapDerivatives(int size, const vector<double>& energy, vector<vector<double> >& c) {
+    vector<double> d1(size*size), d2(size*size), d12(size*size);
+    vector<double> x(size+1), y(size+1), deriv(size+1);
+    for (int i = 0; i < size+1; i++)
+        x[i] = i*2*M_PI/size;
+
+    // Compute derivatives with respect to the first angle.
+
+    for (int i = 0; i < size; i++) {
+        for (int j = 0; j < size; j++)
+            y[j] = energy[j+size*i];
+        y[size] = energy[size*i];
+        SplineFitter::createPeriodicSpline(x, y, deriv);
+        for (int j = 0; j < size; j++) {
+            d1[j+size*i] = SplineFitter::evaluateSplineDerivative(x, y, deriv, x[j]);
+        }
+    }
+
+    // Compute derivatives with respect to the second angle.
+
+    for (int i = 0; i < size; i++) {
+        for (int j = 0; j < size; j++)
+            y[j] = energy[i+size*j];
+        y[size] = energy[i];
+        SplineFitter::createPeriodicSpline(x, y, deriv);
+        for (int j = 0; j < size; j++)
+            d2[i+size*j] = SplineFitter::evaluateSplineDerivative(x, y, deriv, x[j]);
+    }
+
+    // Compute cross derivatives.
+
+    for (int i = 0; i < size; i++) {
+        for (int j = 0; j < size; j++)
+            y[j] = d2[j+size*i];
+        y[size] = d2[size*i];
+        SplineFitter::createPeriodicSpline(x, y, deriv);
+        for (int j = 0; j < size; j++)
+            d12[j+size*i] = SplineFitter::evaluateSplineDerivative(x, y, deriv, x[j]);
+    }
+
+    // Now compute the coefficients.
+
+    const int wt[] = {
+        1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4,
+        0, 0, 0, 0, 0, 0, 0, 0, 3, 0, -9, 6, -2, 0, 6, -4,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4,
+        0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4,
+        0, 0, 0, 0, 1, 0, -3, 2, -2, 0, 6, -4, 1, 0, -3, 2,
+        0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2,
+        0, 0, 0, 0, 0, 0, 3, -2, 0, 0, -6, 4, 0, 0, 3, -2,
+        0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 2, -2,
+        0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2,
+        0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 1, -2, 1,
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
+        0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
+    };
+    vector<double> rhs(16);
+    double delta = 2*M_PI/size;
+    c.resize(size*size);
+    for (int i = 0; i < size; i++) {
+        for (int j = 0; j < size; j++) {
+            // Compute the 16 coefficients for patch (i, j).
+
+            int nexti = (i+1)%size;
+            int nextj = (j+1)%size;
+            double e[] = {energy[i+j*size], energy[nexti+j*size], energy[nexti+nextj*size], energy[i+nextj*size]};
+            double e1[] = {d1[i+j*size], d1[nexti+j*size], d1[nexti+nextj*size], d1[i+nextj*size]};
+            double e2[] = {d2[i+j*size], d2[nexti+j*size], d2[nexti+nextj*size], d2[i+nextj*size]};
+            double e12[] = {d12[i+j*size], d12[nexti+j*size], d12[nexti+nextj*size], d12[i+nextj*size]};
+
+            for (int k = 0; k < 4; k++) {
+                rhs[k] = e[k];
+                rhs[k+4] = e1[k]*delta;
+                rhs[k+8] = e2[k]*delta;
+                rhs[k+12] = e12[k]*delta*delta;
+            }
+            vector<double>& coeff = c[i+j*size];
+            coeff.resize(16);
+            for (int k = 0; k < 16; k++) {
+                double sum = 0.0;
+                for (int m = 0; m < 16; m++)
+                    sum += wt[k+16*m]*rhs[m];
+                coeff[k] = sum;
+            }
+        }
+    }
+}
\ No newline at end of file

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -62,6 +62,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcPeriodicTorsionForceKernel(name, platform);
     if (name == CalcRBTorsionForceKernel::Name())
         return new ReferenceCalcRBTorsionForceKernel(name, platform);
+    if (name == CalcCMAPTorsionForceKernel::Name())
+        return new ReferenceCalcCMAPTorsionForceKernel(name, platform);
     if (name == CalcCustomTorsionForceKernel::Name())
         return new ReferenceCalcCustomTorsionForceKernel(name, platform);
     if (name == CalcGBSAOBCForceKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -37,6 +37,7 @@
 #include "SimTKReference/ReferenceBondForce.h"
 #include "SimTKReference/ReferenceBrownianDynamics.h"
 #include "SimTKReference/ReferenceCCMAAlgorithm.h"
+#include "SimTKReference/ReferenceCMAPTorsionIxn.h"
 #include "SimTKReference/ReferenceCustomAngleIxn.h"
 #include "SimTKReference/ReferenceCustomBondIxn.h"
 #include "SimTKReference/ReferenceCustomExternalIxn.h"
@@ -59,6 +60,7 @@
 #include "openmm/System.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/CustomHbondForceImpl.h"
+#include "openmm/internal/CMAPTorsionForceImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "openmm/Integrator.h"
 #include "openmm/OpenMMException.h"
@@ -602,6 +604,48 @@ double ReferenceCalcRBTorsionForceKernel::executeEnergy(ContextImpl& context) {
     return energy;
 }
 
+void ReferenceCalcCMAPTorsionForceKernel::initialize(const System& system, const CMAPTorsionForce& force) {
+    int numMaps = force.getNumMaps();
+    int numTorsions = force.getNumTorsions();
+    coeff.resize(numMaps);
+    vector<double> energy;
+    vector<vector<double> > c;
+    for (int i = 0; i < numMaps; i++) {
+        int size;
+        force.getMapParameters(i, size, energy);
+        CMAPTorsionForceImpl::calcMapDerivatives(size, energy, c);
+        coeff[i].resize(size*size);
+        for (int j = 0; j < size*size; j++) {
+            coeff[i][j].resize(16);
+            for (int k = 0; k < 16; k++)
+                coeff[i][j][k] = c[j][k];
+        }
+    }
+    torsionMaps.resize(numTorsions);
+    torsionIndices.resize(numTorsions);
+    for (int i = 0; i < numTorsions; i++) {
+        torsionIndices[i].resize(8);
+        force.getTorsionParameters(i, torsionMaps[i], torsionIndices[i][0], torsionIndices[i][1], torsionIndices[i][2],
+            torsionIndices[i][3], torsionIndices[i][4], torsionIndices[i][5], torsionIndices[i][6], torsionIndices[i][7]);
+    }
+}
+
+void ReferenceCalcCMAPTorsionForceKernel::executeForces(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = extractForces(context);
+    ReferenceCMAPTorsionIxn torsion(coeff, torsionMaps, torsionIndices);
+    torsion.calculateIxn(posData, forceData, 0);
+}
+
+double ReferenceCalcCMAPTorsionForceKernel::executeEnergy(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
+    RealOpenMM totalEnergy = 0;
+    ReferenceCMAPTorsionIxn torsion(coeff, torsionMaps, torsionIndices);
+    torsion.calculateIxn(posData, forceData, &totalEnergy);
+    return totalEnergy;
+}
+
 ReferenceCalcCustomTorsionForceKernel::~ReferenceCalcCustomTorsionForceKernel() {
     disposeIntArray(torsionIndexArray, numTorsions);
     disposeRealArray(torsionParamArray, numTorsions);

platforms/reference/src/ReferenceKernels.h

@@ -416,6 +416,39 @@ private:
     RealOpenMM **torsionParamArray;
 };
 
+/**
+ * This kernel is invoked by CMAPTorsionForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcCMAPTorsionForceKernel : public CalcCMAPTorsionForceKernel {
+public:
+    ReferenceCalcCMAPTorsionForceKernel(std::string name, const Platform& platform) : CalcCMAPTorsionForceKernel(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CMAPTorsionForce this kernel will be used for
+     */
+    void initialize(const System& system, const CMAPTorsionForce& force);
+    /**
+     * Execute the kernel to calculate the forces.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    void executeForces(ContextImpl& context);
+    /**
+     * Execute the kernel to calculate the energy.
+     *
+     * @param context    the context in which to execute this kernel
+     * @return the potential energy due to the CMAPTorsionForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    std::vector<std::vector<std::vector<RealOpenMM> > > coeff;
+    std::vector<int> torsionMaps;
+    std::vector<std::vector<int> > torsionIndices;
+};
+
 /**
  * This kernel is invoked by CustomTorsionForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/reference/src/ReferencePlatform.cpp

@@ -49,6 +49,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcCustomAngleForceKernel::Name(), factory);
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
+    registerKernelFactory(CalcCMAPTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomNonbondedForceKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceCMAPTorsionIxn.cpp

+
+/* Portions copyright (c) 2010 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include "ReferenceCMAPTorsionIxn.h"
+#include "ReferenceForce.h"
+
+using std::vector;
+
+/**---------------------------------------------------------------------------------------
+
+   Constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCMAPTorsionIxn::ReferenceCMAPTorsionIxn(const vector<vector<vector<RealOpenMM> > >& coeff,
+        const vector<int>& torsionMaps, const vector<vector<int> >& torsionIndices) :
+        coeff(coeff), torsionMaps(torsionMaps), torsionIndices(torsionIndices) {
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate torsion interaction
+
+   @param atomIndices      bond indices
+   @param atomCoordinates  atom coordinates
+   @param parameters       parameters
+   @param forces           force array (forces added)
+   @param energyByBond     bond energy
+   @param energy           atom energy
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceCMAPTorsionIxn::calculateIxn(RealOpenMM** atomCoordinates, RealOpenMM** forces, RealOpenMM* totalEnergy) const {
+    for (int i = 0; i < torsionMaps.size(); i++)
+        calculateOneIxn(i, atomCoordinates, forces, totalEnergy);
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate the interaction due to a single torsion pair
+
+   @param index            the index of the torsion
+   @param atomCoordinates  atom coordinates
+   @param forces           force array (forces added)
+   @param totalEnergy      total energy
+
+     --------------------------------------------------------------------------------------- */
+
+void ReferenceCMAPTorsionIxn::calculateOneIxn(int index, RealOpenMM** atomCoordinates, RealOpenMM** forces,
+                     RealOpenMM* totalEnergy) const {
+    int map = torsionMaps[index];
+    int a1 = torsionIndices[index][0];
+    int a2 = torsionIndices[index][1];
+    int a3 = torsionIndices[index][2];
+    int a4 = torsionIndices[index][3];
+    int b1 = torsionIndices[index][4];
+    int b2 = torsionIndices[index][5];
+    int b3 = torsionIndices[index][6];
+    int b4 = torsionIndices[index][7];
+
+    // Compute deltas between the various atoms involved.
+
+    RealOpenMM deltaA[3][ReferenceForce::LastDeltaRIndex];
+    ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a1], deltaA[0]);
+    ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a3], deltaA[1]);
+    ReferenceForce::getDeltaR(atomCoordinates[a4], atomCoordinates[a3], deltaA[2]);
+    RealOpenMM deltaB[3][ReferenceForce::LastDeltaRIndex];
+    ReferenceForce::getDeltaR(atomCoordinates[b2], atomCoordinates[b1], deltaB[0]);
+    ReferenceForce::getDeltaR(atomCoordinates[b2], atomCoordinates[b3], deltaB[1]);
+    ReferenceForce::getDeltaR(atomCoordinates[b4], atomCoordinates[b3], deltaB[2]);
+
+    // Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
+
+    RealOpenMM crossProductMemory[12];
+    RealOpenMM* cpA[2];
+    cpA[0] = crossProductMemory;
+    cpA[1] = crossProductMemory + 3;
+    RealOpenMM* cpB[2];
+    cpB[0] = crossProductMemory + 6;
+    cpB[1] = crossProductMemory + 9;
+
+   // Compute the dihedral angles.
+
+    RealOpenMM dotDihedral;
+    RealOpenMM signOfAngle;
+    RealOpenMM angleA =  getDihedralAngleBetweenThreeVectors(deltaA[0], deltaA[1], deltaA[2],
+         cpA, &dotDihedral, deltaA[0], &signOfAngle, 1);
+    RealOpenMM angleB =  getDihedralAngleBetweenThreeVectors(deltaB[0], deltaB[1], deltaB[2],
+         cpB, &dotDihedral, deltaB[0], &signOfAngle, 1);
+    angleA = fmod(angleA+2.0*M_PI, 2.0*M_PI);
+    angleB = fmod(angleB+2.0*M_PI, 2.0*M_PI);
+
+    // Identify which patch this is in.
+
+    int size = (int) SQRT(coeff[map].size());
+    RealOpenMM delta = 2*M_PI/size;
+    int s = (int) (angleA/delta);
+    int t = (int) (angleB/delta);
+    const vector<RealOpenMM>& c = coeff[map][s+size*t];
+    RealOpenMM da = angleA/delta-s;
+    RealOpenMM db = angleB/delta-t;
+
+    // Evaluate the spline to determine the energy and gradients.
+
+    RealOpenMM energy = 0;
+    RealOpenMM dEdA = 0;
+    RealOpenMM dEdB = 0;
+    for (int i = 3; i >= 0; i--) {
+        energy = da*energy + ((c[i*4+3]*db + c[i*4+2])*db + c[i*4+1])*db + c[i*4+0];
+        dEdA = db*dEdA + (3.0*c[i+3*4]*da + 2.0*c[i+2*4])*da + c[i+1*4];
+        dEdB = da*dEdB + (3.0*c[i*4+3]*db + 2.0*c[i*4+2])*db + c[i*4+1];
+    }
+    dEdA /= delta;
+    dEdB /= delta;
+    if (totalEnergy != NULL)
+        *totalEnergy += energy;
+
+    // Apply the force to the first torsion.
+
+    RealOpenMM forceFactors[4];
+    RealOpenMM normCross1 = DOT3(cpA[0], cpA[0]);
+    RealOpenMM normBC = deltaA[1][ReferenceForce::RIndex];
+    forceFactors[0] = (-dEdA*normBC)/normCross1;
+    RealOpenMM normCross2 = DOT3(cpA[1], cpA[1]);
+    forceFactors[3] = (dEdA*normBC)/normCross2;
+    forceFactors[1] = DOT3(deltaA[0], deltaA[1]);
+    forceFactors[1] /= deltaA[1][ReferenceForce::R2Index];
+    forceFactors[2] = DOT3(deltaA[2], deltaA[1]);
+    forceFactors[2] /= deltaA[1][ReferenceForce::R2Index];
+    for (int i = 0; i < 3; i++) {
+        RealOpenMM f0 = forceFactors[0]*cpA[0][i];
+        RealOpenMM f3 = forceFactors[3]*cpA[1][i];
+        RealOpenMM s = forceFactors[1]*f0 - forceFactors[2]*f3;
+        forces[a1][i] += f0;
+        forces[a2][i] -= f0-s;
+        forces[a3][i] -= f3+s;
+        forces[a4][i] += f3;
+    }
+
+    // Apply the force to the second torsion.
+
+    normCross1 = DOT3(cpB[0], cpB[0]);
+    normBC = deltaB[1][ReferenceForce::RIndex];
+    forceFactors[0] = (-dEdB*normBC)/normCross1;
+    normCross2 = DOT3(cpB[1], cpB[1]);
+    forceFactors[3] = (dEdB*normBC)/normCross2;
+    forceFactors[1] = DOT3(deltaB[0], deltaB[1]);
+    forceFactors[1] /= deltaB[1][ReferenceForce::R2Index];
+    forceFactors[2] = DOT3(deltaB[2], deltaB[1]);
+    forceFactors[2] /= deltaB[1][ReferenceForce::R2Index];
+    for (int i = 0; i < 3; i++) {
+        RealOpenMM f0 = forceFactors[0]*cpB[0][i];
+        RealOpenMM f3 = forceFactors[3]*cpB[1][i];
+        RealOpenMM s = forceFactors[1]*f0 - forceFactors[2]*f3;
+        forces[b1][i] += f0;
+        forces[b2][i] -= f0-s;
+        forces[b3][i] -= f3+s;
+        forces[b4][i] += f3;
+    }
+}
+
+// ---------------------------------------------------------------------------------------
+
+/**---------------------------------------------------------------------------------------
+
+   This is present only because we must define it to subclass ReferenceBondIxn.  It is never called.
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceCMAPTorsionIxn::calculateBondIxn(int* atomIndices, RealOpenMM** atomCoordinates,
+     RealOpenMM* parameters, RealOpenMM** forces, RealOpenMM* energyByBond, RealOpenMM* energyByAtom) const {
+}

platforms/reference/src/SimTKReference/ReferenceCMAPTorsionIxn.h

+
+/* Portions copyright (c) 2010 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef __ReferenceCMAPTorsionIxn_H__
+#define __ReferenceCMAPTorsionIxn_H__
+
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "ReferenceBondIxn.h"
+#include <vector>
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceCMAPTorsionIxn : public ReferenceBondIxn {
+
+private:
+
+    std::vector<std::vector<std::vector<RealOpenMM> > > coeff;
+    std::vector<int> torsionMaps;
+    std::vector<std::vector<int> > torsionIndices;
+
+    /**---------------------------------------------------------------------------------------
+
+       Calculate the interaction due to a single torsion pair
+
+       @param index            the index of the torsion
+       @param atomCoordinates  atom coordinates
+       @param forces           force array (forces added)
+       @param totalEnergy      total energy
+
+         --------------------------------------------------------------------------------------- */
+
+    void calculateOneIxn(int index, RealOpenMM** atomCoordinates, RealOpenMM** forces,
+                         RealOpenMM* totalEnergy) const;
+
+public:
+
+    /**---------------------------------------------------------------------------------------
+
+       Constructor
+
+       --------------------------------------------------------------------------------------- */
+
+    ReferenceCMAPTorsionIxn(const std::vector<std::vector<std::vector<RealOpenMM> > >& coeff,
+                            const std::vector<int>& torsionMaps,
+                            const std::vector<std::vector<int> >& torsionIndices);
+
+    /**---------------------------------------------------------------------------------------
+
+       Calculate torsion interaction
+
+       @param atomCoordinates    atom coordinates
+       @param forces             force array (forces added)
+       @param totalEnergy        total energy
+
+         --------------------------------------------------------------------------------------- */
+
+    void calculateIxn(RealOpenMM** atomCoordinates, RealOpenMM** forces, RealOpenMM* totalEnergy) const;
+
+    /**---------------------------------------------------------------------------------------
+
+       This is present only because we must define it to subclass ReferenceBondIxn.  It is never called.
+
+       --------------------------------------------------------------------------------------- */
+
+    int calculateBondIxn(int* atomIndices, RealOpenMM** atomCoordinates,
+                         RealOpenMM* parameters, RealOpenMM** forces,
+                         RealOpenMM* energyByBond, RealOpenMM* energyByAtom) const;
+
+// ---------------------------------------------------------------------------------------
+
+};
+
+#endif // __ReferenceCMAPTorsionIxn_H__

platforms/reference/tests/TestReferenceCMAPTorsionForce.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) 2010 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests all the different force terms in the reference implementation of CMAPTorsionForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/CMAPTorsionForce.h"
+#include "openmm/PeriodicTorsionForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testCMAPTorsions() {
+    const int mapSize = 36;
+
+    // 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);
+    PeriodicTorsionForce* periodic = new PeriodicTorsionForce();
+    periodic->addTorsion(0, 1, 2, 3, 2, M_PI/4, 1.5);
+    periodic->addTorsion(1, 2, 3, 4, 3, M_PI/3, 2.0);
+    system1.addForce(periodic);
+    System system2;
+    for (int i = 0; i < 5; i++)
+        system2.addParticle(1.0);
+    CMAPTorsionForce* cmap = new CMAPTorsionForce();
+    vector<double> mapEnergy(mapSize*mapSize);
+    for (int i = 0; i < mapSize; i++) {
+        double angle1 = i*2*M_PI/mapSize;
+        double energy1 = 1.5*(1+cos(2*angle1-M_PI/4));
+        for (int j = 0; j < mapSize; j++) {
+            double angle2 = j*2*M_PI/mapSize;
+            double energy2 = 2.0*(1+cos(3*angle2-M_PI/3));
+            mapEnergy[i+j*mapSize] = energy1+energy2;
+        }
+    }
+    cmap->addMap(mapSize, mapEnergy);
+    cmap->addTorsion(0, 0, 1, 2, 3, 1, 2, 3, 4);
+    system2.addForce(cmap);
+
+    // Set the atoms in various positions, and verify that both systems give equal forces and energy.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(5);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    Context c1(system1, integrator1, platform);
+    Context c2(system2, integrator2, platform);
+    for (int i = 0; i < 50; i++) {
+        for (int j = 0; j < (int) positions.size(); j++)
+            positions[j] = Vec3(5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt), 5.0*genrand_real2(sfmt));
+        c1.setPositions(positions);
+        c2.setPositions(positions);
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        for (int i = 0; i < system1.getNumParticles(); i++)
+            ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 0.05);
+        ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-3);
+    }
+}
+
+int main() {
+    try {
+        testCMAPTorsions();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+