Add a test for AmoebaAngleForce. There's a problem here, though, since the input

angles are in degrees (unlike *every* other angle force out there), and
stripUnits automatically reduces all angles to radians if they come in with
units.

The approach here is to *slightly* change the API, so that
AmoebaAngleForce.addAngle interprets input angles <2*pi as radians, and >2*pi as
degrees. This is heuristic, but should work in every case out in the wild so
far.

I've also updated the documentation to reflect this behavior, and fixed the
units attached to the return value of AmoebaAngleForce.getAngleParameters() to
return degrees instead of radians.

plugins/amoeba/openmmapi/include/openmm/AmoebaAngleForce.h

@@ -127,7 +127,7 @@ public:
      * @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 length        the angle measured in degrees
+     * @param length        the equilibrium angle, measured in degrees if >2*pi, otherwise interpreted as radians
      * @param quadratic k   the quadratic force constant for the angle, measured in kJ/mol/radian^2
      * @return the index of the angle that was added
      */
@@ -140,7 +140,7 @@ public:
      * @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 length        the equilibrium angle, measured in degress
+     * @param length        the equilibrium angle, measured in degrees if >2*pi, otherwise interpreted as radians
      * @param quadratic k   the quadratic force constant for the angle, measured in kJ/mol/radian^2
      */
     void getAngleParameters(int index, int& particle1, int& particle2, int& particle3, double& length, double& quadraticK) const;
@@ -152,7 +152,7 @@ public:
      * @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 length        the equilibrium angle, measured in degrees
+     * @param length        the equilibrium angle, measured in degrees if >2*pi, otherwise interpreted as radians
      * @param quadratic k   the quadratic force constant for the angle, measured in kJ/mol/radian^2
      */
     void setAngleParameters(int index, int particle1, int particle2, int particle3, double length, double quadraticK);

plugins/amoeba/openmmapi/src/AmoebaAngleForce.cpp

@@ -29,6 +29,7 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include <cmath>
 #include "openmm/Force.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/AmoebaAngleForce.h"
@@ -41,6 +42,7 @@ AmoebaAngleForce::AmoebaAngleForce() {
 }
 
 int AmoebaAngleForce::addAngle(int particle1, int particle2, int particle3,  double length, double quadraticK) {
+    length = length < 2*M_PI ? length * 180.0/M_PI : length;
     angles.push_back(AngleInfo(particle1, particle2, particle3, length, quadraticK));
     return angles.size()-1;
 }
@@ -59,7 +61,7 @@ void AmoebaAngleForce::setAngleParameters(int index, int particle1, int particle
     angles[index].particle1  = particle1;
     angles[index].particle2  = particle2;
     angles[index].particle3  = particle3;
-    angles[index].length     = length;
+    angles[index].length     = length < 2*M_PI ? length * 180.0/M_PI : length;
     angles[index].quadraticK = quadraticK;
 }
 

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -218,11 +218,11 @@ UNITS = {
 ("AmoebaGeneralizedKirkwoodForce",       "getProbeRadius")                                :  ( 'unit.nanometer', ()),
 ("AmoebaGeneralizedKirkwoodForce",       "getSurfaceAreaFactor")                          :  ( '(unit.nanometer*unit.nanometer)/unit.kilojoule_per_mole',()),
 
-("AmoebaAngleForce",             "getAmoebaGlobalAngleCubic")             :  ( None,()),
-("AmoebaAngleForce",             "getAmoebaGlobalAngleQuartic")           :  ( None,()),
-("AmoebaAngleForce",             "getAmoebaGlobalAnglePentic")            :  ( None,()),
-("AmoebaAngleForce",             "getAmoebaGlobalAngleSextic")            :  ( None,()),
-("AmoebaAngleForce",             "getAngleParameters")                            :  ( None, (None, None, None, 'unit.radian', 'unit.kilojoule_per_mole/(unit.radian*unit.radian)')),
+("AmoebaAngleForce",             "getAmoebaGlobalAngleCubic")             :  ( '1/unit.radian',()),
+("AmoebaAngleForce",             "getAmoebaGlobalAngleQuartic")           :  ( '1/unit.radian**2',()),
+("AmoebaAngleForce",             "getAmoebaGlobalAnglePentic")            :  ( '1/unit.radian**3',()),
+("AmoebaAngleForce",             "getAmoebaGlobalAngleSextic")            :  ( '1/unit.radian**4',()),
+("AmoebaAngleForce",             "getAngleParameters")                            :  ( None, (None, None, None, 'unit.degree', 'unit.kilojoule_per_mole/(unit.radian*unit.radian)')),
 
 ("AmoebaBondForce",              "getAmoebaGlobalBondCubic")              :  ( '1/unit.nanometer',()),
 ("AmoebaBondForce",              "getAmoebaGlobalBondQuartic")            :  ( '1/unit.nanometer**2',()),

wrappers/python/tests/TestAPIUnits.py

@@ -547,6 +547,49 @@ class TestAPIUnits(unittest.TestCase):
 
     def testAmoebaAngleForce(self):
         """ Tests the AmoebaAngleForce API features """
+        force = AmoebaAngleForce()
+        force.setAmoebaGlobalAngleCubic(1.0)
+        force.setAmoebaGlobalAngleQuartic(2.0/radians**2)
+        force.setAmoebaGlobalAnglePentic(3.0/degrees**3)
+        force.setAmoebaGlobalAngleSextic(4.0)
+
+        self.assertEqual(force.getAmoebaGlobalAngleCubic(), 1.0/radians)
+        self.assertEqual(force.getAmoebaGlobalAngleQuartic(), 2.0/radians**2)
+        self.assertAlmostEqualUnit(force.getAmoebaGlobalAnglePentic(), 3.0/degrees**3)
+        self.assertEqual(force.getAmoebaGlobalAngleSextic(), 4.0/radians**4)
+
+        force.addAngle(0, 1, 2, math.pi*radians, 1.5*kilocalories_per_mole/radians**2)
+        force.addAngle(1, 2, 3, 180*degrees, 1.5*kilocalories_per_mole/radians**2)
+        force.addAngle(2, 3, 4, 109.4, 1.5)
+
+        self.assertEqual(force.getNumAngles(), 3)
+
+        i, j, k, t, tk = force.getAngleParameters(0)
+        self.assertEqual(i, 0)
+        self.assertEqual(j, 1)
+        self.assertEqual(k, 2)
+        self.assertAlmostEqualUnit(t, math.pi*radians)
+        self.assertIs(t.unit, degree)
+        self.assertAlmostEqualUnit(tk, 1.5*kilocalories_per_mole/radians**2)
+        self.assertIs(tk.unit, kilojoules_per_mole/radians**2)
+
+        i, j, k, t, tk = force.getAngleParameters(1)
+        self.assertEqual(i, 1)
+        self.assertEqual(j, 2)
+        self.assertEqual(k, 3)
+        self.assertAlmostEqualUnit(t, 180*degrees)
+        self.assertIs(t.unit, degree)
+        self.assertAlmostEqualUnit(tk, 1.5*kilocalories_per_mole/radians**2)
+        self.assertIs(tk.unit, kilojoules_per_mole/radians**2)
+
+        i, j, k, t, tk = force.getAngleParameters(2)
+        self.assertEqual(i, 2)
+        self.assertEqual(j, 3)
+        self.assertEqual(k, 4)
+        self.assertAlmostEqualUnit(t, 109.4*degrees)
+        self.assertIs(t.unit, degree)
+        self.assertAlmostEqualUnit(tk, 1.5*kilojoules_per_mole/radians**2)
+        self.assertIs(tk.unit, kilojoules_per_mole/radians**2)
 
 if __name__ == '__main__':
     unittest.main()