Merge pull request #826 from peastman/angles

Fixed an incorrect unit conversion when constraining angles with AMOEBA

wrappers/python/simtk/openmm/app/forcefield.py

@@ -34,6 +34,7 @@ __version__ = "1.0"
 import os
 import itertools
 import xml.etree.ElementTree as etree
+import math
 from math import sqrt, cos
 import simtk.openmm as mm
 import simtk.unit as unit
@@ -2101,7 +2102,7 @@ class AmoebaAngleGenerator:
                     hit = 1
                     if isConstrained and self.k[i] != 0.0:
                         angleDict['idealAngle'] = self.angle[i][0]
-                        addAngleConstraint(angle, self.angle[i][0], data, sys)
+                        addAngleConstraint(angle, self.angle[i][0]*math.pi/180.0, data, sys)
                     elif self.k[i] != 0:
                         lenAngle = len(self.angle[i])
                         if (lenAngle > 1):
@@ -2179,7 +2180,7 @@ class AmoebaAngleGenerator:
                     hit = 1
                     angleDict['idealAngle'] = self.angle[i][0]
                     if (isConstrained and self.k[i] != 0.0):
-                        addAngleConstraint(angle, self.angle[i][0], data, sys)
+                        addAngleConstraint(angle, self.angle[i][0]*math.pi/180.0, data, sys)
                     else:
                         force.addAngle(angle[0], angle[1], angle[2], angle[3], self.angle[i][0], self.k[i])
                     break

wrappers/python/tests/TestForceField.py

@@ -5,6 +5,7 @@ from simtk.openmm import *
 from simtk.unit import *
 import simtk.openmm.app.element as elem
 import simtk.openmm.app.forcefield as forcefield
+import math
 
 class TestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method."""
@@ -208,7 +209,7 @@ class AmoebaTestForceField(unittest.TestCase):
         """
 
         self.pdb1 = PDBFile('systems/amoeba-ion-in-water.pdb')
-        self.forcefield1 = ForceField('amoeba2009.xml')
+        self.forcefield1 = ForceField('amoeba2013.xml')
         self.topology1 = self.pdb1.topology
 
 
@@ -253,6 +254,21 @@ class AmoebaTestForceField(unittest.TestCase):
                 if isinstance(force, AmoebaVdwForce):
                     self.assertEqual(useDispersionCorrection, force.getUseDispersionCorrection())
 
+    def test_RigidWater(self):
+        """Test that AMOEBA creates rigid water with the correct geometry."""
+
+        system = self.forcefield1.createSystem(self.pdb1.topology, rigidWater=True)
+        constraints = dict()
+        for i in range(system.getNumConstraints()):
+            p1,p2,dist = system.getConstraintParameters(i)
+            if p1 < 3:
+                constraints[(min(p1,p2), max(p1,p2))] = dist.value_in_unit(nanometers)
+        hoDist = 0.09572
+        hohAngle = 108.50*math.pi/180.0
+        hohDist = math.sqrt(2*hoDist**2 - 2*hoDist**2*math.cos(hohAngle))
+        self.assertAlmostEqual(constraints[(0,1)], hoDist)
+        self.assertAlmostEqual(constraints[(0,2)], hoDist)
+        self.assertAlmostEqual(constraints[(1,2)], hohDist)
 
 if __name__ == '__main__':
     unittest.main()