Added AMOEBA 2018 force field (#3299)

* Adding support for new AMOEBA features

* Support modern method of specifying in-plane angles

* Implemented stretch-torsions

* Implemented angle-torsions

* More AMOEBA fixes

* Bug fix

* Converted AMOEBA 2018 force field

* Added documentation for AMOEBA 2018

* Added a missing file for tests

devtools/forcefield-scripts/processTinkerForceField.py

@@ -35,16 +35,19 @@ import os.path
 # biotype    2006    CA      "Calcium Ion"                     256
 # biotype    2007    CL      "Chloride Ion"                    258
 
-ions    = { 'Li+' :  ['LI', 249],
-            'Na+' :  ['NA', 250],
-            'K+'  :  ['K',  251],
-            'Rb+' :  ['RB', 252],  
-            'Cs+' :  ['CS', 253], 
-            'Be+' :  ['BE', 254],
-            'Mg+' :  ['MG', 255],
-            'Ca+' :  ['CA', 256],
-            'Zn+' :  ['ZN', 257],
-            'Cl-' :  ['Cl', 258]
+ions    = { 'Li+' :  ['LI', 351],
+            'Na+' :  ['NA', 352],
+            'K+'  :  ['K',  353],
+            'Rb+' :  ['RB', 354],  
+            'Cs+' :  ['CS', 355], 
+            'Be2' :  ['BE', 356],
+            'Mg2' :  ['MG', 357],
+            'Ca2' :  ['CA', 358],
+            'Zn2' :  ['ZN', 359],
+            'F-'  :  ['F',  360],
+            'Cl-' :  ['Cl', 361],
+            'Br-' :  ['Br', 362],
+            'I-'  :  ['I',  363]
            }
 
 atomTypes                    = {}
@@ -149,9 +152,9 @@ def buildProteinResidue( residueDict, atoms, bondInfo, abbr, loc, tinkerLookupNa
                     residueDict[abbr]['atoms'][bondedAtom]['bonds'][atom] = 1
                     residueDict[abbr]['atoms'][atom]['bonds'][bondedAtom] = 1
                 else:
-                    print "Error: bonded atom=%s not in residue=%s" % ( atom, abbr )
+                    print("Error: bonded atom=%s not in residue=%s" % ( atom, abbr ))
         else:
-            print "Error: bonded atom=%s nt in residue=%s" % ( atom, abbr )
+            print("Error: bonded atom=%s nt in residue=%s" % ( atom, abbr ))
 
     return
 
@@ -204,7 +207,7 @@ def copyProteinResidue( residue ):
 def buildResidueDict( residueXmlFileName ):
 
     residueTree = etree.parse(residueXmlFileName)
-    print "Read %s" % (residueXmlFileName)
+    print("Read %s" % (residueXmlFileName))
     root        = residueTree.getroot()
     residueDict = dict()
 
@@ -262,7 +265,7 @@ def buildResidueDict( residueXmlFileName ):
                     residueDict[cResidueName]['tinkerLookupName']                 = 'C-Terminal ' + 'HIS ' + sub
                 else:
                     residueDict[cResidueName]['tinkerLookupName']                 = 'C-Terminal ' + abbr + ' ' + sub
-                print "tinkerLookupName %s %s" % ( abbr, residueDict[cResidueName]['tinkerLookupName'])
+                print("tinkerLookupName %s %s" % ( abbr, residueDict[cResidueName]['tinkerLookupName']))
             else:
                 residueDict[cResidueName]['tinkerLookupName']                 = 'C-Terminal ' + abbr
             residueDict[cResidueName]['atoms']['OXT']                         = copyAtom( residueDict[abbr]['atoms']['O'] )
@@ -296,10 +299,10 @@ def buildResidueDict( residueXmlFileName ):
 
             buildProteinResidue( residueDict, atoms, bondInfo, abbr, loc, tinkerName, 1, residueName, type )
 
-    print "Start Lookup XML FFFFinal\n\n"
+    print("Start Lookup XML FFFFinal\n\n")
     printXml = 1
     if( printXml ):
-        print "<Residues>"
+        print("<Residues>")
         for resName in sorted( residueDict.keys() ):
             if( 'include' in residueDict[resName] and residueDict[resName]['include'] ):
                 type         = residueDict[resName]['type']
@@ -307,13 +310,13 @@ def buildResidueDict( residueXmlFileName ):
                 tinkerLookupName   = residueDict[resName]['tinkerLookupName']
                 fullName     = residueDict[resName]['residueName']
                 outputString = """  <Residue abbreviation="%s" loc="%s" type="%s" tinkerLookupName="%s" fullName="%s">""" % (resName, loc, type, tinkerLookupName, fullName )
-                print "%s" % outputString
+                print("%s" % outputString)
 
                 atomsInfo    = residueDict[resName]['atoms']
                 for atomName in sorted( atomsInfo.keys() ):
                     tinkerLookupName = atomsInfo[atomName]['tinkerLookupName']
                     outputString = """  <Atom name="%s" tinkerLookupName="%s" />""" % (atomName, tinkerLookupName)
-                    print "%s" % outputString
+                    print("%s" % outputString)
 
                 includedBonds = dict()
                 for atomName in sorted( atomsInfo.keys() ):
@@ -327,9 +330,9 @@ def buildResidueDict( residueXmlFileName ):
                                 includedBonds[bondedAtom] = dict()
                             includedBonds[atomName][bondedAtom] = 1
                             includedBonds[bondedAtom][atomName] = 1
-                            print "%s" % outputString
-                print "</Residue>"
-        print "</Residues>"
+                            print("%s" % outputString)
+                print("</Residue>")
+        print("</Residues>")
 
     return residueDict
 
@@ -363,13 +366,13 @@ def setBioTypes( bioTypes, residueDict ):
             if lookupName in bioTypes:
                 res['atoms'][atom]['type'] = bioTypes[lookupName][3]
             else:
-                print "For %s lookupName=%s not in biotype" % (atom,lookupName)
+                print("For %s lookupName=%s not in biotype" % (atom,lookupName))
                 if( 'parent' in res ):
                     lookupName =  res['atoms'][atom]['tinkerLookupName']  + '_' +  res['parent']['tinkerLookupName']
                     if( lookupName in bioTypes ):
                         res['atoms'][atom]['type'] = bioTypes[lookupName][3]
                     else:
-                        print "Missing lookupName=%s from biotype" % (lookupName)
+                        print("Missing lookupName=%s from biotype" % (lookupName))
     return 0
 
 #=============================================================================================
@@ -460,12 +463,16 @@ forces                                   = {}
 recognizedForces                         = {}
 recognizedForces['bond']                 = 1
 recognizedForces['angle']                = 1
+recognizedForces['anglep']               = 1
 recognizedForces['strbnd']               = 1
 recognizedForces['ureybrad']             = 1
 recognizedForces['opbend']               = 1
 recognizedForces['torsion']              = 1
 recognizedForces['pitors']               = 1
+recognizedForces['strtors']              = 1
+recognizedForces['angtors']              = 1
 recognizedForces['vdw']                  = 1
+recognizedForces['vdwpr']                = 1
 recognizedForces['polarize']             = 1
 recognizedForces['tortors']              = addTorTor
 recognizedForces['multipole']            = addMultipole
@@ -573,7 +580,7 @@ while lineIndex < len( allLines ):
         scalars[fields[0]] = fields[1]
         lineIndex += 1
     else:
-        print "Field %s not recognized: line=<%s>" % ( fields[0], allLines[lineIndex] )
+        print("Field %s not recognized: line=<%s>" % ( fields[0], allLines[lineIndex] ))
         lineIndex += 1
 
 #=============================================================================================
@@ -589,12 +596,12 @@ setBioTypes( bioTypes, residueDict )
 tinkerXmlFileName           = scalars['forcefield']
 tinkerXmlFileName          += '.xml'
 tinkerXmlFile               = open( tinkerXmlFileName, 'w' )
-print "Opened %s." % (tinkerXmlFileName)
+print("Opened %s." % (tinkerXmlFileName))
 
 gkXmlFileName              = scalars['forcefield']
 gkXmlFileName             += '_gk.xml'
 gkXmlFile                  = open( gkXmlFileName, 'w' )
-print "Opened %s." % (gkXmlFileName)
+print("Opened %s." % (gkXmlFileName))
 
 today = datetime.date.today().isoformat()
 sourceFile = os.path.basename(sys.argv[1])
@@ -680,7 +687,7 @@ for resname, res in sorted(residueDict.items()):
             type  = res['atoms'][atom]['type']
             typeI = int( type )
             if( typeI < 0 ):
-                print "Error: type=%s for atom=%s of residue=%s" % (type, atom, resname)
+                print("Error: type=%s for atom=%s of residue=%s" % (type, atom, resname))
             tag  = "   <Atom name=\"%s\" type=\"%s\" />" % (atom, type)
             atomIndex[atom]  = atomCount
             atomCount       += 1
@@ -755,7 +762,7 @@ tinkerXmlFile.write('  </Residue>\n')
 
 # ions
  
-for ion,ionInfo in ions.iteritems():
+for ion,ionInfo in ions.items():
     outputString  = """  <Residue name="%s">\n""" % (ionInfo[0])
     outputString += """   <Atom name="%s" type="%s"/>\n""" % (ionInfo[0], str(ionInfo[1]))
     outputString += """  </Residue>\n""" 
@@ -792,20 +799,20 @@ if( isAmoeba ):
     sextic          = float(scalars['angle-sextic']) 
     outputString    = """ <AmoebaAngleForce angle-cubic="%s" angle-quartic="%s" angle-pentic="%s" angle-sextic="%s">""" %( str(cubic), str(quartic), str(pentic), str(sextic) )
     tinkerXmlFile.write( "%s\n" % (outputString ) )
-    angles          = forces['angle']
     radian          = 57.2957795130
     radian2         = 4.184/(radian*radian)
-    for angle in angles:
-       k            = float(angle[3])*radian2
-       outputString = """  <Angle class1="%s" class2="%s" class3="%s" k="%s" angle1="%s" """ % (angle[0], angle[1], angle[2], str(k), angle[4] ) 
-       if( len(angle) > 5 ):
-           outputString += """  angle2="%s" """ % (angle[5])
-
-       if( len(angle) > 6 ):
-           outputString += """  angle3="%s" """ % (angle[6])
-       outputString += " /> "
-
-       tinkerXmlFile.write( "%s\n" % (outputString ) )
+    for set in ['angle', 'anglep']:
+        for angle in forces[set]:
+           k            = float(angle[3])*radian2
+           outputString = '  <Angle class1="%s" class2="%s" class3="%s" k="%s" angle1="%s"' % (angle[0], angle[1], angle[2], str(k), angle[4] ) 
+           if( len(angle) > 5 ):
+               outputString += ' angle2="%s"' % (angle[5])
+    
+           if( len(angle) > 6 ):
+               outputString += ' angle3="%s"' % (angle[6])
+           outputString += ' inPlane="%s"/>' % (set == 'anglep')
+    
+           tinkerXmlFile.write( "%s\n" % (outputString ) )
     tinkerXmlFile.write( " </AmoebaAngleForce>\n" )
 
 #=============================================================================================
@@ -823,9 +830,12 @@ if( isAmoeba ):
     opbends         = forces['opbend']
     radian2         = 4.184/(radian*radian)
     for opbend in opbends:
-       k            = float(opbend[4])*radian2
-       outputString = """  <Angle class1="%s" class2="%s" class3="%s" class4="%s" k="%s"/>""" % (opbend[0], opbend[1], opbend[2],  opbend[3], str(k))
-       tinkerXmlFile.write( "%s\n" % (outputString ) )
+        k            = float(opbend[4])*radian2
+        for i in range(4):
+            if opbend[i] == '0':
+                opbend[i] = ''
+        outputString = """  <Angle class1="%s" class2="%s" class3="%s" class4="%s" k="%s"/>""" % (opbend[0], opbend[1], opbend[2],  opbend[3], str(k))
+        tinkerXmlFile.write( "%s\n" % (outputString ) )
     tinkerXmlFile.write( " </AmoebaOutOfPlaneBendForce>\n" )
 
 #=============================================================================================
@@ -869,6 +879,26 @@ if( isAmoeba ):
        tinkerXmlFile.write( "%s\n" % (outputString ) )
     tinkerXmlFile.write( " </AmoebaPiTorsionForce>\n" )
 
+#=============================================================================================
+
+    # Stretch torsion
+
+    tinkerXmlFile.write(' <AmoebaStretchTorsionForce>\n')
+    for torsion in forces['strtors']:
+        v = [float(x)*10*4.184 for x in torsion[4:]]
+        tinkerXmlFile.write(f'  <Torsion class1="{torsion[0]}" class2="{torsion[1]}" class3="{torsion[2]}" class4="{torsion[3]}" v11="{v[0]}" v12="{v[1]}" v13="{v[2]}" v21="{v[3]}" v22="{v[4]}" v23="{v[5]}" v31="{v[6]}" v32="{v[7]}" v33="{v[8]}"/>\n')
+    tinkerXmlFile.write(' </AmoebaStretchTorsionForce>\n')
+
+#=============================================================================================
+
+    # Angle torsion
+
+    tinkerXmlFile.write(' <AmoebaAngleTorsionForce>\n')
+    for torsion in forces['angtors']:
+        v = [float(x)*4.184 for x in torsion[4:]]
+        tinkerXmlFile.write(f'  <Torsion class1="{torsion[0]}" class2="{torsion[1]}" class3="{torsion[2]}" class4="{torsion[3]}" v11="{v[0]}" v12="{v[1]}" v13="{v[2]}" v21="{v[3]}" v22="{v[4]}" v23="{v[5]}"/>\n')
+    tinkerXmlFile.write(' </AmoebaAngleTorsionForce>\n')
+
 #=============================================================================================
 
     # AmoebaStretchBendForce
@@ -906,7 +936,7 @@ if( isAmoeba ):
        outputString  = """  <TorsionTorsionGrid grid="%s" nx="%s" ny="%s" >""" % (str(index), torInfo[5], torInfo[6] )
        tinkerXmlFile.write( "%s\n" % (outputString ) )
        for (gridIndex, gridEntry) in enumerate(grid):
-           print "Gxx %d  %s" % ( gridIndex, str(gridEntry) )
+           print("Gxx %d  %s" % ( gridIndex, str(gridEntry) ))
            if( len( gridEntry ) > 5 ):
                f   = float( gridEntry[2] )*4.184
                fx  = float( gridEntry[3] )*4.184
@@ -930,15 +960,19 @@ if( isAmoeba ):
     outputString         = """ <AmoebaVdwForce type="%s" radiusrule="%s" radiustype="%s" radiussize="%s" epsilonrule="%s" vdw-13-scale="%s" vdw-14-scale="%s" vdw-15-scale="%s" >""" % (
          scalars['vdwtype'], scalars['radiusrule'], scalars['radiustype'], scalars['radiussize'], scalars['epsilonrule'], scalars['vdw-13-scale'], scalars['vdw-14-scale'], scalars['vdw-15-scale'] )
     tinkerXmlFile.write( "%s\n" % (outputString ) )
-    vdws                 = forces['vdw']
-    for vdw in vdws:
+    for vdw in forces['vdw']:
        sigma             = float(vdw[1])*0.1
        epsilon           = float(vdw[2])*4.184
        if( len(vdw) > 3 ): 
            reduction = vdw[3]
        else:
            reduction = 1.0
-       outputString      = """  <Vdw class="%s" sigma="%s" epsilon="%s" reduction="%s" /> """ % (vdw[0], str(sigma), str(epsilon), str(reduction) )
+       outputString      = """  <Vdw class="%s" sigma="%s" epsilon="%s" reduction="%s"/>""" % (vdw[0], str(sigma), str(epsilon), str(reduction))
+       tinkerXmlFile.write( "%s\n" % (outputString ) )
+    for pair in forces['vdwpr']:
+       sigma             = float(pair[2])*0.1
+       epsilon           = float(pair[3])*4.184
+       outputString      = """  <Pair class1="%s" class2="%s" sigma="%s" epsilon="%s"/>""" % (pair[0], pair[1], str(sigma), str(epsilon))
        tinkerXmlFile.write( "%s\n" % (outputString ) )
     tinkerXmlFile.write( " </AmoebaVdwForce>\n" )
 
@@ -1013,11 +1047,11 @@ if( isAmoeba ):
           
        outputString     += "/>"
        tinkerXmlFile.write( "%s\n" % (outputString ) )
-       print m[polarize[0]]
+       print(m[polarize[0]])
     for t in sorted(m):
         for k in m[t]:
             if t not in m[k]:
-                print t, k
+                print(t, k)
 
     tinkerXmlFile.write( " </AmoebaMultipoleForce>\n" )
 
@@ -1037,10 +1071,10 @@ if( isAmoeba ):
     # radii are set in forcefield.py
 
     for type in sorted( atomTypes ):
-        print "atom type=%s  %s" % ( str(type), str(atomTypes[type]) )
+        print("atom type=%s  %s" % ( str(type), str(atomTypes[type]) ))
 
     for type in sorted( bioTypes ):
-        print "bio type=%s  %s" % ( str(type), str(bioTypes[type]) )
+        print("bio type=%s  %s" % ( str(type), str(bioTypes[type]) ))
 
     multipoleArray       = forces['multipole']
     for multipoleInfo in multipoleArray:
@@ -1067,9 +1101,9 @@ if( isAmoeba ):
            elif( element == 'Fe' ):
                shct = 0.88
            else:
-               print "Warning no overlap scale factor for type=%d element=%s" % (type, element)
+               print("Warning no overlap scale factor for type=%d element=%s" % (type, element))
        else:
-           print "Warning no overlap scale factor for type=%d " % (type)
+           print("Warning no overlap scale factor for type=%d " % (type))
 
        outputString      = """  <GeneralizedKirkwood type="%s" charge="%s" shct="%s"  /> """ % ( axisInfo[0], multipoles[0],  str(shct) )
        gkXmlFile.write( "%s\n" % (outputString ) )

docs-source/usersguide/application/02_running_sims.rst

@@ -649,25 +649,28 @@ dielectric constant, and :math:`T` is the temperature in Kelvin.
 AMOEBA
 ------
 
-The AMOEBA polarizable force field provides parameters for proteins, water, and ions.
+The AMOEBA polarizable force field provides parameters for proteins, nucleic acids, water, and ions.
 
 .. tabularcolumns:: |l|L|
 
 =============================  ================================================================================
 File                           Parameters
 =============================  ================================================================================
-:file:`amoeba2013.xml`         AMOEBA 2013\ :cite:`Shi2013`
-:file:`amoeba2013_gk.xml`      Generalized Kirkwood solvation model\ :cite:`Schnieders2007` for use with AMOEBA 2013 force field
+:file:`amoeba2018.xml`         AMOEBA 2018\ :cite:`Shi2013`
+:file:`amoeba2018_gk.xml`      Generalized Kirkwood solvation model\ :cite:`Schnieders2007` for use with AMOEBA 2018 force field
+:file:`amoeba2013.xml`         AMOEBA 2013.  This force field is deprecated.  It is
+                               recommended to use AMOEBA 2018 instead.
+:file:`amoeba2013_gk.xml`      Generalized Kirkwood solvation model for use with AMOEBA 2013 force field
 :file:`amoeba2009.xml`         AMOEBA 2009\ :cite:`Ren2002`.  This force field is deprecated.  It is
-                               recommended to use AMOEBA 2013 instead.
+                               recommended to use AMOEBA 2018 instead.
 :file:`amoeba2009_gk.xml`      Generalized Kirkwood solvation model for use with AMOEBA 2009 force field
 =============================  ================================================================================
 
-For explicit solvent simulations, just include the single file :file:`amoeba2013.xml`.
+For explicit solvent simulations, just include the single file :file:`amoeba2018.xml`.
 AMOEBA also supports implicit solvent using a Generalized Kirkwood model.  To enable
-it, also include :file:`amoeba2013_gk.xml`.
+it, also include :file:`amoeba2018_gk.xml`.
 
-The older AMOEBA 2009 force field is provided only for backward compatibility, and is not
+The older AMOEBA 2009 and 2013 force fields are provided only for backward compatibility, and are not
 recommended for most simulations.
 
 CHARMM Polarizable Force Field

wrappers/python/tests/TestForceField.py

@@ -1282,7 +1282,7 @@ class AmoebaTestForceField(unittest.TestCase):
         forcefield = ForceField('amoeba2013.xml', 'amoeba2013_gk.xml')
         system = forcefield.createSystem(pdb.topology, polarization='direct')
         integrator = VerletIntegrator(0.001)
-        context = Context(system, integrator)
+        context = Context(system, integrator, Platform.getPlatformByName('Reference'))
         context.setPositions(pdb.positions)
         state1 = context.getState(getForces=True)
         with open('systems/alanine-dipeptide-amoeba-forces.xml') as input:
@@ -1291,6 +1291,55 @@ class AmoebaTestForceField(unittest.TestCase):
             diff = norm(f1-f2)
             self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-3)
 
+    def computeAmoeba18Energies(self, filename):
+        pdb = PDBFile(filename)
+        forcefield = ForceField('amoeba2018.xml')
+        system = forcefield.createSystem(pdb.topology, polarization='mutual', mutualInducedTargetEpsilon=1e-5)
+        for i, f in enumerate(system.getForces()):
+            f.setForceGroup(i)
+        integrator = VerletIntegrator(0.001)
+        context = Context(system, integrator, Platform.getPlatformByName('Reference'))
+        context.setPositions(pdb.positions)
+        energies = {}
+        for i, f in enumerate(system.getForces()):
+            state = context.getState(getEnergy=True, groups={i})
+            energies[f.getName()] = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+        return energies
+
+    def test_Amoeba18BPTI(self):
+        """Test that AMOEBA18 computes energies correctly for BPTI."""
+        energies = self.computeAmoeba18Energies('systems/bpti.pdb')
+
+        # Compare to values computed with Tinker.
+
+        self.assertAlmostEqual(290.2445, energies['AmoebaBond'], 4)
+        self.assertAlmostEqual(496.4300, energies['AmoebaAngle']+energies['AmoebaInPlaneAngle'], 4)
+        self.assertAlmostEqual(51.2913, energies['AmoebaOutOfPlaneBend'], 4)
+        self.assertAlmostEqual(5.7695, energies['AmoebaStretchBend'], 4)
+        self.assertAlmostEqual(75.6890, energies['PeriodicTorsionForce'], 4)
+        self.assertAlmostEqual(19.3364, energies['AmoebaPiTorsion'], 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)
+
+    def test_Amoeba18DNA(self):
+        """Test that AMOEBA18 computes energies correctly for DNA."""
+        energies = self.computeAmoeba18Energies('systems/dna.pdb')
+
+        # Compare to values computed with Tinker.
+
+        self.assertAlmostEqual(688.7794, energies['AmoebaBond'], 4)
+        self.assertAlmostEqual(430.5887, energies['AmoebaAngle']+energies['AmoebaInPlaneAngle'], 4)
+        self.assertAlmostEqual(10.5715, energies['AmoebaOutOfPlaneBend'], 4)
+        self.assertAlmostEqual(0.7773, energies['AmoebaStretchBend'], 4)
+        self.assertAlmostEqual(77.0482, energies['PeriodicTorsionForce'], 4)
+        self.assertAlmostEqual(56.5654, energies['AmoebaPiTorsion'], 4)
+        self.assertAlmostEqual(-0.5457, energies['AmoebaStretchTorsion'], 4)
+        self.assertAlmostEqual(-1.6044, energies['AmoebaAngleTorsion'], 4)
+        self.assertAlmostEqual(106.4062, energies['AmoebaVdwForce'], 4)
+        self.assertAlmostEqual(90.5715-73.6051, energies['AmoebaMultipoleForce'], 4)
+        self.assertAlmostEqual(1385.5529, sum(list(energies.values())), 4)
 
 if __name__ == '__main__':
     unittest.main()