Created CustomAngleForce, including reference implementation

olla/include/openmm/kernels.h

@@ -35,6 +35,7 @@
 #include "openmm/AndersenThermostat.h"
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMMotionRemover.h"
+#include "openmm/CustomAngleForce.h"
 #include "openmm/CustomBondForce.h"
 #include "openmm/CustomExternalForce.h"
 #include "openmm/CustomGBForce.h"
@@ -267,6 +268,38 @@ public:
     virtual double executeEnergy(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class CalcCustomAngleForceKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "CalcCustomAngleForce";
+    }
+    CalcCustomAngleForceKernel(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 CustomAngleForce this kernel will be used for
+     */
+    virtual void initialize(const System& system, const CustomAngleForce& 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 CustomAngleForce
+     */
+    virtual double executeEnergy(ContextImpl& context) = 0;
+};
+
 /**
  * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
  */

openmmapi/include/openmm/CustomAngleForce.h

+#ifndef OPENMM_CUSTOMANGLEFORCE_H_
+#define OPENMM_CUSTOMANGLEFORCE_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 interactions between sets of three particles that depend on the angle between them.
+ * Unlike HarmonicAngleForce, the functional form of the interaction is completely customizable, and may
+ * involve arbitrary algebraic expressions.  In addition to the angle formed by the particles, it may depend
+ * on arbitrary global and per-angle parameters.
+ *
+ * To use this class, create a CustomAngleForce object, passing an algebraic expression to the constructor
+ * that defines the interaction energy between each set of particles.  The expression may depend on theta, the angle
+ * formed by the particles, as well as on any parameters you choose.  Then call addPerAngleParameter() to define per-angle
+ * parameters, and addGlobalParameter() to define global parameters.  The values of per-angle parameters are specified as
+ * part of the system definition, while values of global parameters may be modified during a simulation by calling Context::setParameter().
+ * Finally, call addAngle() once for each angle.  After an angle has been added, you can modify its parameters by calling setAngleParameters().
+ *
+ * As an example, the following code creates a CustomAngleForce that implements a harmonic potential:
+ *
+ * <tt>CustomAngleForce* force = new CustomAngleForce("0.5*k*(theta-theta0)^2");</tt>
+ *
+ * This force depends on two parameters: the spring constant k and equilibrium angle theta0.  The following code defines these parameters:
+ *
+ * <tt><pre>
+ * force->addPerAngleParameter("k");
+ * force->addPerAngleParameter("theta0");
+ * </pre></tt>
+ *
+ * Expressions may involve the operators + (add), - (subtract), * (multiply), / (divide), and ^ (power), and the following
+ * functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, step.  All trigonometric functions
+ * are defined in radians, and log is the natural logarithm.  step(x) = 0 if x is less than 0, 1 otherwise.
+ */
+
+class OPENMM_EXPORT CustomAngleForce : public Force {
+public:
+    /**
+     * Create a CustomAngleForce.
+     *
+     * @param energy    an algebraic expression giving the interaction energy between three particles as a function
+     *                  of theta, the angle between them
+     */
+    CustomAngleForce(const std::string& energy);
+    /**
+     * Get the number of angles for which force field parameters have been defined.
+     */
+    int getNumAngles() const {
+        return angles.size();
+    }
+    /**
+     * Get the number of per-angle parameters that the interaction depends on.
+     */
+    int getNumPerAngleParameters() const {
+        return parameters.size();
+    }
+    /**
+     * Get the number of global parameters that the interaction depends on.
+     */
+    int getNumGlobalParameters() const {
+        return globalParameters.size();
+    }
+    /**
+     * Get the algebraic expression that gives the interaction energy for each angle
+     */
+    const std::string& getEnergyFunction() const;
+    /**
+     * Set the algebraic expression that gives the interaction energy for each angle
+     */
+    void setEnergyFunction(const std::string& energy);
+    /**
+     * Add a new per-angle parameter that the interaction may depend on.
+     *
+     * @param name      the name of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addPerAngleParameter(const std::string& name);
+    /**
+     * Get the name of a per-angle parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getPerAngleParameterName(int index) const;
+    /**
+     * Set the name of a per-angle parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setPerAngleParameterName(int index, const std::string& name);
+    /**
+     * Add a new global parameter that the interaction may depend on.
+     *
+     * @param name             the name of the parameter
+     * @param defaultValue     the default value of the parameter
+     * @return the index of the parameter that was added
+     */
+    int addGlobalParameter(const std::string& name, double defaultValue);
+    /**
+     * Get the name of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the name
+     * @return the parameter name
+     */
+    const std::string& getGlobalParameterName(int index) const;
+    /**
+     * Set the name of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the name
+     * @param name           the name of the parameter
+     */
+    void setGlobalParameterName(int index, const std::string& name);
+    /**
+     * Get the default value of a global parameter.
+     *
+     * @param index     the index of the parameter for which to get the default value
+     * @return the parameter default value
+     */
+    double getGlobalParameterDefaultValue(int index) const;
+    /**
+     * Set the default value of a global parameter.
+     *
+     * @param index          the index of the parameter for which to set the default value
+     * @param name           the default value of the parameter
+     */
+    void setGlobalParameterDefaultValue(int index, double defaultValue);
+    /**
+     * Add an angle term to the force field.
+     *
+     * @param particle1     the index of the first particle connected by the angle
+     * @param particle2     the index of the second particle connected by the angle
+     * @param particle3     the index of the third particle connected by the angle
+     * @param parameters    the list of parameters for the new angle
+     * @return the index of the angle that was added
+     */
+    int addAngle(int particle1, int particle2, int particle3, const std::vector<double>& parameters);
+    /**
+     * Get the force field parameters for an angle term.
+     *
+     * @param index         the index of the angle for which to get parameters
+     * @param particle1     the index of the first particle connected by the angle
+     * @param particle2     the index of the second particle connected by the angle
+     * @param particle3     the index of the third particle connected by the angle
+     * @param parameters    the list of parameters for the angle
+     */
+    void getAngleParameters(int index, int& particle1, int& particle2, int& particle3, std::vector<double>& parameters) const;
+    /**
+     * Set the force field parameters for an angle term.
+     *
+     * @param index         the index of the angle for which to set parameters
+     * @param particle1     the index of the first particle connected by the angle
+     * @param particle2     the index of the second particle connected by the angle
+     * @param particle3     the index of the third particle connected by the angle
+     * @param parameters    the list of parameters for the angle
+     */
+    void setAngleParameters(int index, int particle1, int particle2, int particle3, const std::vector<double>& parameters);
+protected:
+    ForceImpl* createImpl();
+private:
+    class AngleInfo;
+    class AngleParameterInfo;
+    class GlobalParameterInfo;
+    std::string energyExpression;
+    std::vector<AngleParameterInfo> parameters;
+    std::vector<GlobalParameterInfo> globalParameters;
+    std::vector<AngleInfo> angles;
+};
+
+/**
+ * This is an internal class used to record information about an angle.
+ * @private
+ */
+class CustomAngleForce::AngleInfo {
+public:
+    int particle1, particle2, particle3;
+    std::vector<double> parameters;
+    AngleInfo() : particle1(-1), particle2(-1), particle3(-1) {
+    }
+    AngleInfo(int particle1, int particle2, int particle3, const std::vector<double>& parameters) :
+            particle1(particle1), particle2(particle2), particle3(particle3), parameters(parameters) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a per-angle parameter.
+ * @private
+ */
+class CustomAngleForce::AngleParameterInfo {
+public:
+    std::string name;
+    AngleParameterInfo() {
+    }
+    AngleParameterInfo(const std::string& name) : name(name) {
+    }
+};
+
+/**
+ * This is an internal class used to record information about a global parameter.
+ * @private
+ */
+class CustomAngleForce::GlobalParameterInfo {
+public:
+    std::string name;
+    double defaultValue;
+    GlobalParameterInfo() {
+    }
+    GlobalParameterInfo(const std::string& name, double defaultValue) : name(name), defaultValue(defaultValue) {
+    }
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMANGLEFORCE_H_*/

openmmapi/include/openmm/internal/CustomAngleForceImpl.h

+#ifndef OPENMM_CUSTOMANGLEFORCEIMPL_H_
+#define OPENMM_CUSTOMANGLEFORCEIMPL_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 "ForceImpl.h"
+#include "openmm/CustomAngleForce.h"
+#include "openmm/Kernel.h"
+#include <utility>
+#include <map>
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of CustomAngleForce.
+ */
+
+class CustomAngleForceImpl : public ForceImpl {
+public:
+    CustomAngleForceImpl(CustomAngleForce& owner);
+    ~CustomAngleForceImpl();
+    void initialize(ContextImpl& context);
+    CustomAngleForce& 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();
+    std::vector<std::string> getKernelNames();
+private:
+    CustomAngleForce& owner;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_CUSTOMANGLEFORCEIMPL_H_*/

openmmapi/src/CustomAngleForce.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/CustomAngleForce.h"
+#include "openmm/internal/CustomAngleForceImpl.h"
+#include <cmath>
+#include <map>
+#include <sstream>
+#include <utility>
+
+using namespace OpenMM;
+using std::string;
+using std::stringstream;
+using std::vector;
+
+CustomAngleForce::CustomAngleForce(const string& energy) : energyExpression(energy) {
+}
+
+const string& CustomAngleForce::getEnergyFunction() const {
+    return energyExpression;
+}
+
+void CustomAngleForce::setEnergyFunction(const std::string& energy) {
+    energyExpression = energy;
+}
+
+int CustomAngleForce::addPerAngleParameter(const string& name) {
+    parameters.push_back(AngleParameterInfo(name));
+    return parameters.size()-1;
+}
+
+const string& CustomAngleForce::getPerAngleParameterName(int index) const {
+    return parameters[index].name;
+}
+
+void CustomAngleForce::setPerAngleParameterName(int index, const string& name) {
+    parameters[index].name = name;
+}
+
+int CustomAngleForce::addGlobalParameter(const string& name, double defaultValue) {
+    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
+    return globalParameters.size()-1;
+}
+
+const string& CustomAngleForce::getGlobalParameterName(int index) const {
+    return globalParameters[index].name;
+}
+
+void CustomAngleForce::setGlobalParameterName(int index, const string& name) {
+    globalParameters[index].name = name;
+}
+
+double CustomAngleForce::getGlobalParameterDefaultValue(int index) const {
+    return globalParameters[index].defaultValue;
+}
+
+void CustomAngleForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
+    globalParameters[index].defaultValue = defaultValue;
+}
+
+int CustomAngleForce::addAngle(int particle1, int particle2, int particle3, const vector<double>& parameters) {
+    angles.push_back(AngleInfo(particle1, particle2, particle3, parameters));
+    return angles.size()-1;
+}
+
+void CustomAngleForce::getAngleParameters(int index, int& particle1, int& particle2, int& particle3, std::vector<double>& parameters) const {
+    particle1 = angles[index].particle1;
+    particle2 = angles[index].particle2;
+    particle3 = angles[index].particle3;
+    parameters = angles[index].parameters;
+}
+
+void CustomAngleForce::setAngleParameters(int index, int particle1, int particle2, int particle3, const vector<double>& parameters) {
+    angles[index].parameters = parameters;
+    angles[index].particle1 = particle1;
+    angles[index].particle2 = particle2;
+    angles[index].particle3 = particle3;
+}
+
+ForceImpl* CustomAngleForce::createImpl() {
+    return new CustomAngleForceImpl(*this);
+}
+

openmmapi/src/CustomAngleForceImpl.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/OpenMMException.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/CustomAngleForceImpl.h"
+#include "openmm/kernels.h"
+#include <sstream>
+
+using namespace OpenMM;
+using std::map;
+using std::pair;
+using std::vector;
+using std::set;
+using std::string;
+using std::stringstream;
+
+CustomAngleForceImpl::CustomAngleForceImpl(CustomAngleForce& owner) : owner(owner) {
+}
+
+CustomAngleForceImpl::~CustomAngleForceImpl() {
+}
+
+void CustomAngleForceImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(CalcCustomAngleForceKernel::Name(), context);
+
+    // Check for errors in the specification of bonds.
+
+    System& system = context.getSystem();
+    vector<double> parameters;
+    int numParameters = owner.getNumPerAngleParameters();
+    for (int i = 0; i < owner.getNumAngles(); i++) {
+        int particle1, particle2, particle3;
+        owner.getAngleParameters(i, particle1, particle2, particle3, parameters);
+        if (particle1 < 0 || particle1 >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomAngleForce: Illegal particle index for an angle: ";
+            msg << particle1;
+            throw OpenMMException(msg.str());
+        }
+        if (particle2 < 0 || particle2 >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomAngleForce: Illegal particle index for an angle: ";
+            msg << particle2;
+            throw OpenMMException(msg.str());
+        }
+        if (particle3 < 0 || particle3 >= system.getNumParticles()) {
+            stringstream msg;
+            msg << "CustomAngleForce: Illegal particle index for an angle: ";
+            msg << particle3;
+            throw OpenMMException(msg.str());
+        }
+        if (parameters.size() != numParameters) {
+            stringstream msg;
+            msg << "CustomAngleForce: Wrong number of parameters for angle ";
+            msg << i;
+            throw OpenMMException(msg.str());
+        }
+    }
+    dynamic_cast<CalcCustomAngleForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
+}
+
+void CustomAngleForceImpl::calcForces(ContextImpl& context) {
+    dynamic_cast<CalcCustomAngleForceKernel&>(kernel.getImpl()).executeForces(context);
+}
+
+double CustomAngleForceImpl::calcEnergy(ContextImpl& context) {
+    return dynamic_cast<CalcCustomAngleForceKernel&>(kernel.getImpl()).executeEnergy(context);
+}
+
+vector<string> CustomAngleForceImpl::getKernelNames() {
+    vector<string> names;
+    names.push_back(CalcCustomAngleForceKernel::Name());
+    return names;
+}
+
+map<string, double> CustomAngleForceImpl::getDefaultParameters() {
+    map<string, double> parameters;
+    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
+        parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
+    return parameters;
+}
+

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -54,6 +54,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceCalcHarmonicAngleForceKernel(name, platform);
     if (name == CalcHarmonicAngleForceKernel::Name())
         return new ReferenceCalcHarmonicAngleForceKernel(name, platform);
+    if (name == CalcCustomAngleForceKernel::Name())
+        return new ReferenceCalcCustomAngleForceKernel(name, platform);
     if (name == CalcPeriodicTorsionForceKernel::Name())
         return new ReferenceCalcPeriodicTorsionForceKernel(name, platform);
     if (name == CalcRBTorsionForceKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -37,6 +37,7 @@
 #include "SimTKReference/ReferenceBondForce.h"
 #include "SimTKReference/ReferenceBrownianDynamics.h"
 #include "SimTKReference/ReferenceCCMAAlgorithm.h"
+#include "SimTKReference/ReferenceCustomAngleIxn.h"
 #include "SimTKReference/ReferenceCustomBondIxn.h"
 #include "SimTKReference/ReferenceCustomExternalIxn.h"
 #include "SimTKReference/ReferenceCustomGBIxn.h"
@@ -277,7 +278,7 @@ void ReferenceCalcCustomBondForceKernel::initialize(const System& system, const
 
     // Build the arrays.
 
-    bondIndexArray = allocateIntArray(numBonds, numParameters);
+    bondIndexArray = allocateIntArray(numBonds, 2);
     bondParamArray = allocateRealArray(numBonds, numParameters);
     vector<double> params;
     for (int i = 0; i < force.getNumBonds(); ++i) {
@@ -373,6 +374,70 @@ double ReferenceCalcHarmonicAngleForceKernel::executeEnergy(ContextImpl& context
     return energy;
 }
 
+ReferenceCalcCustomAngleForceKernel::~ReferenceCalcCustomAngleForceKernel() {
+    disposeIntArray(angleIndexArray, numAngles);
+    disposeRealArray(angleParamArray, numAngles);
+}
+
+void ReferenceCalcCustomAngleForceKernel::initialize(const System& system, const CustomAngleForce& force) {
+    numAngles = force.getNumAngles();
+    int numParameters = force.getNumPerAngleParameters();
+
+    // Build the arrays.
+
+    angleIndexArray = allocateIntArray(numAngles, 3);
+    angleParamArray = allocateRealArray(numAngles, numParameters);
+    vector<double> params;
+    for (int i = 0; i < force.getNumAngles(); ++i) {
+        int particle1, particle2, particle3;
+        force.getAngleParameters(i, particle1, particle2, particle3, params);
+        angleIndexArray[i][0] = particle1;
+        angleIndexArray[i][1] = particle2;
+        angleIndexArray[i][2] = particle3;
+        for (int j = 0; j < numParameters; j++)
+            angleParamArray[i][j] = (RealOpenMM) params[j];
+    }
+
+    // Parse the expression used to calculate the force.
+
+    Lepton::ParsedExpression expression = Lepton::Parser::parse(force.getEnergyFunction()).optimize();
+    energyExpression = expression.createProgram();
+    forceExpression = expression.differentiate("theta").optimize().createProgram();
+    for (int i = 0; i < numParameters; i++)
+        parameterNames.push_back(force.getPerAngleParameterName(i));
+    for (int i = 0; i < force.getNumGlobalParameters(); i++)
+        globalParameterNames.push_back(force.getGlobalParameterName(i));
+}
+
+void ReferenceCalcCustomAngleForceKernel::executeForces(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = extractForces(context);
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ReferenceBondForce refBondForce;
+    ReferenceCustomAngleIxn harmonicAngle(energyExpression, forceExpression, parameterNames, globalParameters);
+    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, 0, 0, 0, harmonicAngle);
+}
+
+double ReferenceCalcCustomAngleForceKernel::executeEnergy(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    RealOpenMM** forceData = allocateRealArray(context.getSystem().getNumParticles(), 3);
+    RealOpenMM* energyArray = new RealOpenMM[numAngles];
+    RealOpenMM energy = 0;
+    map<string, double> globalParameters;
+    for (int i = 0; i < (int) globalParameterNames.size(); i++)
+        globalParameters[globalParameterNames[i]] = context.getParameter(globalParameterNames[i]);
+    ReferenceBondForce refBondForce;
+    ReferenceCustomAngleIxn harmonicAngle(energyExpression, forceExpression, parameterNames, globalParameters);
+    for (int i = 0; i < numAngles; ++i)
+        energyArray[i] = 0;
+    refBondForce.calculateForce(numAngles, angleIndexArray, posData, angleParamArray, forceData, energyArray, 0, &energy, harmonicAngle);
+    disposeRealArray(forceData, context.getSystem().getNumParticles());
+    delete[] energyArray;
+    return energy;
+}
+
 ReferenceCalcPeriodicTorsionForceKernel::~ReferenceCalcPeriodicTorsionForceKernel() {
     disposeIntArray(torsionIndexArray, numTorsions);
     disposeRealArray(torsionParamArray, numTorsions);

platforms/reference/src/ReferenceKernels.h

@@ -262,6 +262,42 @@ private:
     RealOpenMM **angleParamArray;
 };
 
+/**
+ * This kernel is invoked by CustomAngleForce to calculate the forces acting on the system and the energy of the system.
+ */
+class ReferenceCalcCustomAngleForceKernel : public CalcCustomAngleForceKernel {
+public:
+    ReferenceCalcCustomAngleForceKernel(std::string name, const Platform& platform) : CalcCustomAngleForceKernel(name, platform) {
+    }
+    ~ReferenceCalcCustomAngleForceKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param force      the CustomAngleForce this kernel will be used for
+     */
+    void initialize(const System& system, const CustomAngleForce& 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 CustomAngleForce
+     */
+    double executeEnergy(ContextImpl& context);
+private:
+    int numAngles;
+    int **angleIndexArray;
+    RealOpenMM **angleParamArray;
+    Lepton::ExpressionProgram energyExpression, forceExpression;
+    std::vector<std::string> parameterNames, globalParameterNames;
+};
+
 /**
  * This kernel is invoked by PeriodicTorsionForce to calculate the forces acting on the system and the energy of the system.
  */

platforms/reference/src/ReferencePlatform.cpp

@@ -45,6 +45,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(CalcHarmonicBondForceKernel::Name(), factory);
     registerKernelFactory(CalcCustomBondForceKernel::Name(), factory);
     registerKernelFactory(CalcHarmonicAngleForceKernel::Name(), factory);
+    registerKernelFactory(CalcCustomAngleForceKernel::Name(), factory);
     registerKernelFactory(CalcPeriodicTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcRBTorsionForceKernel::Name(), factory);
     registerKernelFactory(CalcNonbondedForceKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.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 <string.h>
+#include <sstream>
+
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "../SimTKUtilities/SimTKOpenMMLog.h"
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "ReferenceCustomAngleIxn.h"
+#include "ReferenceForce.h"
+
+using namespace std;
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomAngleIxn constructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomAngleIxn::ReferenceCustomAngleIxn(const Lepton::ExpressionProgram& energyExpression,
+        const Lepton::ExpressionProgram& forceExpression, const vector<string>& parameterNames, map<string, double> globalParameters) :
+        energyExpression(energyExpression), forceExpression(forceExpression), paramNames(parameterNames), globalParameters(globalParameters) {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomAngleIxn::ReferenceCustomAngleIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceCustomAngleIxn destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceCustomAngleIxn::~ReferenceCustomAngleIxn( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomAngleIxn::~ReferenceCustomAngleIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Calculate Custom Angle Ixn
+
+   @param atomIndices      atom indices of atom participating in bond
+   @param atomCoordinates  atom coordinates
+   @param parameters       parameters values
+   @param forces           force array (forces added to input values)
+   @param energiesByBond   energies by bond: energiesByBond[bondIndex]
+   @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
+
+   @return ReferenceForce::DefaultReturn;
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceCustomAngleIxn::calculateBondIxn( int* atomIndices,
+                                                RealOpenMM** atomCoordinates,
+                                                RealOpenMM* parameters,
+                                                RealOpenMM** forces,
+                                                RealOpenMM* energiesByBond,
+                                                RealOpenMM* energiesByAtom ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceCustomAngleIxn::calculateAngleIxn";
+
+   // ---------------------------------------------------------------------------------------
+
+   static const std::string methodName = "\nReferenceCustomAngleIxn::calculateAngleIxn";
+
+   static const RealOpenMM zero        = 0.0;
+   static const RealOpenMM one         = 1.0;
+
+   RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
+   map<string, double> variables = globalParameters;
+   for (int i = 0; i < (int) paramNames.size(); ++i)
+       variables[paramNames[i]] = parameters[i];
+
+   // ---------------------------------------------------------------------------------------
+
+   // Compute the angle between the three atoms
+
+   int atomAIndex = atomIndices[0];
+   int atomBIndex = atomIndices[1];
+   int atomCIndex = atomIndices[2];
+   ReferenceForce::getDeltaR(atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR[0]);
+   ReferenceForce::getDeltaR(atomCoordinates[atomCIndex], atomCoordinates[atomBIndex], deltaR[1]);
+   RealOpenMM pVector[3];
+   SimTKOpenMMUtilities::crossProductVector3(deltaR[0], deltaR[1], pVector);
+   RealOpenMM rp = SQRT(DOT3(pVector, pVector));
+   if (rp < 1.0e-06)
+      rp = (RealOpenMM) 1.0e-06;
+   RealOpenMM dot = DOT3(deltaR[0], deltaR[1]);
+   RealOpenMM cosine = dot/SQRT((deltaR[0][ReferenceForce::R2Index]*deltaR[1][ReferenceForce::R2Index]));
+   RealOpenMM angle;
+   if (cosine >= one)
+      angle = zero;
+   else if (cosine <= -one)
+      angle = PI_M;
+   else
+      angle = ACOS(cosine);
+   variables["theta"] = angle;
+
+   // Compute the force and energy, and apply them to the atoms.
+   
+   RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables);
+   RealOpenMM dEdR = (RealOpenMM) forceExpression.evaluate(variables);
+   RealOpenMM termA =  dEdR/(deltaR[0][ReferenceForce::R2Index]*rp);
+   RealOpenMM termC = -dEdR/(deltaR[1][ReferenceForce::R2Index]*rp);
+
+   RealOpenMM deltaCrossP[3][3];
+   SimTKOpenMMUtilities::crossProductVector3(deltaR[0], pVector, deltaCrossP[0]);
+   SimTKOpenMMUtilities::crossProductVector3(deltaR[1], pVector, deltaCrossP[2]);
+
+   for (int ii = 0; ii < 3; ii++) {
+      deltaCrossP[0][ii] *= termA;
+      deltaCrossP[2][ii] *= termC;
+      deltaCrossP[1][ii]  = -(deltaCrossP[0][ii]+deltaCrossP[2][ii]);
+   }
+
+   // accumulate forces
+
+   for (int jj = 0; jj < 3; jj++) {
+      for (int ii = 0; ii < 3; ii++) {
+         forces[atomIndices[jj]][ii] += deltaCrossP[jj][ii];
+      }
+   }
+
+   // accumulate energies
+
+   updateEnergy(energy, energiesByBond, 3, atomIndices, energiesByAtom);
+
+   return ReferenceForce::DefaultReturn;
+}
+

platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.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 __ReferenceCustomAngleIxn_H__
+#define __ReferenceCustomAngleIxn_H__
+
+#include "ReferenceBondIxn.h"
+#include "lepton/ExpressionProgram.h"
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceCustomAngleIxn : public ReferenceBondIxn {
+
+   private:
+      Lepton::ExpressionProgram energyExpression;
+      Lepton::ExpressionProgram forceExpression;
+      std::vector<std::string> paramNames;
+      std::map<std::string, double> globalParameters;
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceCustomAngleIxn(const Lepton::ExpressionProgram& energyExpression, const Lepton::ExpressionProgram& forceExpression,
+                              const std::vector<std::string>& parameterNames, std::map<std::string, double> globalParameters);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceCustomAngleIxn( );
+
+      /**---------------------------------------------------------------------------------------
+
+         Calculate Custom Angle Ixn
+
+         @param atomIndices      two bond indices
+         @param atomCoordinates  atom coordinates
+         @param parameters       parameter values
+         @param forces           force array (forces added)
+         @param energiesByBond   energies by bond: energiesByBond[bondIndex]
+         @param energiesByAtom   energies by atom: energiesByAtom[atomIndex]
+
+         --------------------------------------------------------------------------------------- */
+
+      int calculateBondIxn( int* atomIndices, RealOpenMM** atomCoordinates,
+                            RealOpenMM* parameters, RealOpenMM** forces,
+                            RealOpenMM* energiesByBond, RealOpenMM* energiesByAtom ) const;
+
+
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // _ReferenceCustomAngleIxn___

platforms/reference/tests/TestReferenceCustomAngleForce.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) 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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the reference implementation of CustomAngleForce.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/CustomAngleForce.h"
+#include "openmm/HarmonicAngleForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testAngles() {
+    ReferencePlatform platform;
+
+    // Create a system using a CustomAngleForce.
+
+    System customSystem;
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    customSystem.addParticle(1.0);
+    CustomAngleForce* custom = new CustomAngleForce("scale*k*(theta-theta0)^2");
+    custom->addPerAngleParameter("theta0");
+    custom->addPerAngleParameter("k");
+    custom->addGlobalParameter("scale", 0.5);
+    vector<double> parameters(2);
+    parameters[0] = 1.5;
+    parameters[1] = 0.8;
+    custom->addAngle(0, 1, 2, parameters);
+    customSystem.addForce(custom);
+
+    // Create an identical system using a HarmonicAngleForce.
+
+    System harmonicSystem;
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    harmonicSystem.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    HarmonicAngleForce* harmonic = new HarmonicAngleForce();
+    harmonic->addAngle(0, 1, 2, 1.5, 0.8);
+    harmonicSystem.addForce(harmonic);
+
+    // Set the atoms in various positions, and verify that both systems give identical forces and energy.
+
+    init_gen_rand(0);
+    vector<Vec3> positions(3);
+    VerletIntegrator integrator1(0.01);
+    VerletIntegrator integrator2(0.01);
+    for (int i = 0; i < 10; i++) {
+        Context c1(customSystem, integrator1, platform);
+        Context c2(harmonicSystem, integrator2, platform);
+        positions[0] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
+        positions[1] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
+        positions[2] = Vec3(5.0*genrand_real2(), 5.0*genrand_real2(), 5.0*genrand_real2());
+        c1.setPositions(positions);
+        c2.setPositions(positions);
+        State s1 = c1.getState(State::Forces | State::Energy);
+        State s2 = c2.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = s1.getForces();
+        ASSERT_EQUAL_VEC(s1.getForces()[0], s2.getForces()[0], TOL);
+        ASSERT_EQUAL_VEC(s1.getForces()[1], s2.getForces()[1], TOL);
+        ASSERT_EQUAL_VEC(s1.getForces()[2], s2.getForces()[2], TOL);
+        ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), TOL);
+    }
+}
+
+int main() {
+    try {
+        testAngles();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+
+