Update AMOEBA XML, .prm files, and parsing of Tinker files (#5086)

* Correct AmoebaAngleTorsion in test_Amoeba18Nucleic of TestForceField

* Update processTinkerForceField to handle latest .prm Tinker files

* Update amoeba2018 XML files

* Update amoeba2013 XML files

* Update amoeba2009 XML files

* Adapt addTorTor to new format in the .prm files

* Fix TorsionTorsion

* Also update the total energy in test_Amoeba18Nucleic

* Update amoebabio18.prm

* Fix nucleic acid test energies

* Correct AmoebaAngleTorsionForce params

* Add new addTorTor to TinkerFiles

* Revert unit fix

* Change to .pdb file which Tinker likes

* Update test_Amoeba18BPTI

* Remove trailing zeros from XML files

* Leave trailing zeros only on 2018 ff

* New element names in 2018

* More digits for surfaceAreaFactor

* More digits for surfaceAreaFactor

* More digits

* Remove debugging print

* Add support to 2009 and 2013 AMOEBA ffs to processTinkerForceField.py

* Add FF specific residues XML files

* Delete old residuesFinal.xml

* Update AMOEBA XML ffs

* Update FFs

* Fix some formatting issues

* Fix "." in scientific notation

* Remove old assertions

wrappers/python/openmm/app/tinkerfiles.py

@@ -268,13 +268,25 @@ class TinkerFiles:
                 forces["tortors"] = []
             fields = allLines[lineIndex]
             tortorInfo = fields[1:]
-            lastGridLine = lineIndex + int(fields[6]) * int(fields[7])
+            nx = int(fields[6])
+            ny = int(fields[7])
+            totalGridPoints = nx * ny
             grid = []
-            while lineIndex < lastGridLine:
-                lineIndex += 1
-                grid.append(allLines[lineIndex])
-            forces["tortors"].append([tortorInfo, grid])
-            return lineIndex
+            currentLineIndex = lineIndex + 1
+            gridPointsProcessed = 0
+            while gridPointsProcessed < totalGridPoints and currentLineIndex < len(allLines):
+                lineFields = allLines[currentLineIndex]
+                i = 0
+                while i + 3 <= len(lineFields) and gridPointsProcessed < totalGridPoints:
+                    angle1 = lineFields[i]
+                    angle2 = lineFields[i + 1] 
+                    f = lineFields[i + 2]
+                    grid.append([angle1, angle2, f])
+                    gridPointsProcessed += 1
+                    i += 3
+                currentLineIndex += 1
+            forces['tortors'].append( [ tortorInfo, grid ] )
+            return (currentLineIndex - 1)
 
         RECOGNIZED_FORCES = {
             "bond": 1,

wrappers/python/tests/TestForceField.py

@@ -1953,7 +1953,24 @@ class AmoebaTestForceField(unittest.TestCase):
         return energies
 
     def test_Amoeba18BPTI(self):
-        """Test that AMOEBA18 computes energies correctly for BPTI."""
+        """
+        Test that AMOEBA18 computes energies correctly for BPTI.
+        
+        Total Potential Energy :               -259.8636 Kcal/mole
+
+        Energy Component Breakdown :           Kcal/mole        Interactions
+
+        Bond Stretching                         290.2445              906
+        Angle Bending                           496.4300             1626
+        Stretch-Bend                              5.7695             1455
+        Out-of-Plane Bend                        51.2913              597
+        Torsional Angle                          75.6890             2391
+        Pi-Orbital Torsion                       19.3364              109
+        Torsion-Torsion                         -32.6689                6
+        Van der Waals                           383.8705           394854
+        Atomic Multipoles                     -1325.1825           394854
+        Polarization                           -224.6434           394854
+        """
         energies = self.computeAmoeba18Energies('systems/bpti.pdb')
 
         # Compare to values computed with Tinker.
@@ -1966,8 +1983,8 @@ class AmoebaTestForceField(unittest.TestCase):
         self.assertAlmostEqual(19.3364, energies['AmoebaPiTorsionForce'], 4)
         self.assertAlmostEqual(-32.6689, energies['AmoebaTorsionTorsionForce'], 4)
         self.assertAlmostEqual(383.8705, energies['AmoebaVdwForce'], 4)
-        self.assertAlmostEqual(-1323.5640-225.3660, energies['AmoebaMultipoleForce'], 2)
-        self.assertAlmostEqual(-258.9676, sum(list(energies.values())), 2)
+        self.assertAlmostEqual(-1325.1825-224.6434, energies['AmoebaMultipoleForce'], 2)
+        self.assertAlmostEqual(-259.8636, sum(list(energies.values())), 2)
 
     def test_Amoeba18Nucleic(self):
         """Test that AMOEBA18 computes energies correctly for DNA and RNA."""
@@ -1985,7 +2002,7 @@ class AmoebaTestForceField(unittest.TestCase):
         self.assertAlmostEqual(-5.0402, energies['AmoebaAngleTorsionForce'], 4)
         self.assertAlmostEqual(187.1103, energies['AmoebaVdwForce'], 4)
         self.assertAlmostEqual(1635.1289-236.1484, energies['AmoebaMultipoleForce'], 3)
-        self.assertAlmostEqual(3146.3046, sum(list(energies.values())), 3)
+        self.assertAlmostEqual(3146.3045, sum(list(energies.values())), 3)
 
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/TestTinkerFiles.py

@@ -175,7 +175,7 @@ class TestTinkerFiles(unittest.TestCase):
 
         Intermolecular Energy :                 896.3435 Kcal/mole
 
-        Total Potential Energy :               3151.2568 Kcal/mole
+        Total Potential Energy :               3146.3045 Kcal/mole
 
         Energy Component Breakdown :           Kcal/mole        Interactions
 
@@ -186,7 +186,7 @@ class TestTinkerFiles(unittest.TestCase):
         Torsional Angle                         166.7233             2197
         Pi-Orbital Torsion                       57.2066              142
         Stretch-Torsion                          -4.2538               68
-        Angle-Torsion                            -0.0880              112
+        Angle-Torsion                            -5.0402              112
         Van der Waals                           187.1103           291451
         Atomic Multipoles                      1635.1289           291451
         Polarization                           -236.1484           291451
@@ -229,10 +229,10 @@ class TestTinkerFiles(unittest.TestCase):
         self.assertEnergyEqual(166.7233, energies["PeriodicTorsionForce"])
         self.assertEnergyEqual(57.2066, energies["AmoebaPiTorsionForce"])
         self.assertEnergyEqual(-4.2538, energies["AmoebaStretchTorsionForce"])
-        self.assertEnergyEqual(-0.0880, energies["AmoebaAngleTorsionForce"], 1e-4)
+        self.assertEnergyEqual(-5.0402, energies["AmoebaAngleTorsionForce"], 1e-3)
         self.assertEnergyEqual(187.1103, energies["AmoebaVdwForce"])
         self.assertEnergyEqual(1635.1289 - 236.1484, energies["AmoebaMultipoleForce"])
-        self.assertEnergyEqual(3151.2568, sum(list(energies.values())))
+        self.assertEnergyEqual(3146.3045, sum(list(energies.values())))
 
     def test_Amoeba13ForcesImplicit(self):
         """Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""