Changed forcefield to not use dispersion corrections for CHARMM

devtools/forcefield-scripts/processCharmmForceField.py

@@ -377,7 +377,7 @@ for values in sorted(uniqueCmaps, key=lambda x: uniqueCmaps[x]):
 for map in cmaps:
     print('   <Torsion map="%d" class1="%s" class2="%s" class3="%s" class4="%s" class5="%s"/>' % (uniqueCmaps[map.values], map.classes[0], map.classes[1], map.classes[2], map.classes[3], map.classes[4]))
 print(' </CMAPTorsionForce>')
-print(' <NonbondedForce coulomb14scale="1.0" lj14scale="1.0">')
+print(' <NonbondedForce coulomb14scale="1.0" lj14scale="1.0" useDispersionCorrection="False">')
 for type in atomTypes:
     print('  <Atom type="%d" charge="%g" sigma="1.0" epsilon="1.0"/>' % (type.type, type.charge))
 print(' </NonbondedForce>')
@@ -464,6 +464,7 @@ print("""]
 
 customNonbondedForce = mm.CustomNonbondedForce('4*eps*((sig/r)^12-(sig/r)^6); eps=epsilon(type1, type2); sig=sigma(type1, type2)')
 customNonbondedForce.setNonbondedMethod(min(nonbondedForce.getNonbondedMethod(), 2))
+customNonbondedForce.setUseLongRangeCorrection(False)
 customNonbondedForce.addTabulatedFunction('epsilon', mm.Discrete2DFunction(numAtomClasses, numAtomClasses, epsilon))
 customNonbondedForce.addTabulatedFunction('sigma', mm.Discrete2DFunction(numAtomClasses, numAtomClasses, sigma))
 customNonbondedForce.addPerParticleParameter('type')

wrappers/python/simtk/openmm/app/data/charmm36.xml

@@ -102920,7 +102920,7 @@
     <Torsion map="6" type1="C" type2="NH1" type3="CT" type4="C" type5="NH1"/>
     <Torsion map="7" type1="C" type2="N" type3="CT2" type4="C" type5="N"/>
   </CMAPTorsionForce>
-  <NonbondedForce coulomb14scale="1.0" lj14scale="1.0">
+  <NonbondedForce coulomb14scale="1.0" lj14scale="1.0" useDispersionCorrection="False">
     <UseAttributeFromResidue name="charge"/>
     <Atom epsilon="0.0" sigma="1.0" type="H"/>
     <Atom epsilon="0.0" sigma="1.0" type="HC"/>
@@ -103335,7 +103335,7 @@
     <Atom epsilon="0.0" sigma="1.0" type="NE"/>
     <Atom epsilon="0.0" sigma="1.0" type="DUM"/>
   </NonbondedForce>
-  <LennardJonesForce lj14scale="1.0">
+  <LennardJonesForce lj14scale="1.0" useDispersionCorrection="False">
     <Atom epsilon="0.192464" sigma="0.04000135244450124" type="H"/>
     <Atom epsilon="0.192464" sigma="0.04000135244450124" type="HC"/>
     <Atom epsilon="0.092048" sigma="0.2351972615890496" type="HA"/>

wrappers/python/simtk/openmm/app/data/charmm_polar_2013.xml

@@ -7200,7 +7200,7 @@
    <Torsion map="2" class1="CD2O1A" class2="ND3A3" class3="CD31C" class4="CD2O1A" class5="ND2A2"/>
    <Torsion map="2" class1="CD2O1A" class2="ND3A3" class3="CD31C" class4="CD2O1A" class5="ND3A3"/>
  </CMAPTorsionForce>
- <NonbondedForce coulomb14scale="1.0" lj14scale="1.0">
+ <NonbondedForce coulomb14scale="1.0" lj14scale="1.0" useDispersionCorrection="False">
   <Atom type="0" charge="1.71639" sigma="1.0" epsilon="1.0"/>
   <Atom type="1" charge="-1.11466" sigma="1.0" epsilon="1.0"/>
   <Atom type="2" charge="0.55733" sigma="1.0" epsilon="1.0"/>
@@ -11370,6 +11370,7 @@ sigma14 = [
 customNonbondedForce = mm.CustomNonbondedForce('4*eps*((sig/r)^12-(sig/r)^6); eps=epsilon(type1, type2); sig=sigma(type1, type2)')
 customNonbondedForce.setNonbondedMethod(min(nonbondedForce.getNonbondedMethod(), 2))
+customNonbondedForce.setUseLongRangeCorrection(False)
 customNonbondedForce.setCutoffDistance(nonbondedForce.getCutoffDistance())
 customNonbondedForce.addTabulatedFunction('epsilon', mm.Discrete2DFunction(numAtomClasses, numAtomClasses, epsilon))
 customNonbondedForce.addTabulatedFunction('sigma', mm.Discrete2DFunction(numAtomClasses, numAtomClasses, sigma))

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

@@ -37,6 +37,7 @@ import os
 import itertools
 import xml.etree.ElementTree as etree
 import math
+import warnings
 from math import sqrt, cos
 from copy import deepcopy
 from collections import defaultdict
@@ -2273,10 +2274,11 @@ class NonbondedGenerator(object):
 
     SCALETOL = 1e-5
 
-    def __init__(self, forcefield, coulomb14scale, lj14scale):
+    def __init__(self, forcefield, coulomb14scale, lj14scale, useDispersionCorrection):
         self.ff = forcefield
         self.coulomb14scale = coulomb14scale
         self.lj14scale = lj14scale
+        self.useDispersionCorrection = useDispersionCorrection
         self.params = ForceField._AtomTypeParameters(forcefield, 'NonbondedForce', 'Atom', ('charge', 'sigma', 'epsilon'))
 
     def registerAtom(self, parameters):
@@ -2285,15 +2287,31 @@ class NonbondedGenerator(object):
     @staticmethod
     def parseElement(element, ff):
         existing = [f for f in ff._forces if isinstance(f, NonbondedGenerator)]
+        if 'useDispersionCorrection' in element.attrib:
+            useDispersionCorrection = bool(eval(element.attrib.get('useDispersionCorrection')))
+        else:
+            useDispersionCorrection = None
         if len(existing) == 0:
-            generator = NonbondedGenerator(ff, float(element.attrib['coulomb14scale']), float(element.attrib['lj14scale']))
+            generator = NonbondedGenerator(
+                ff,
+                float(element.attrib['coulomb14scale']),
+                float(element.attrib['lj14scale']),
+                useDispersionCorrection
+            )
             ff.registerGenerator(generator)
         else:
             # Multiple <NonbondedForce> tags were found, probably in different files.  Simply add more types to the existing one.
             generator = existing[0]
-            if abs(generator.coulomb14scale - float(element.attrib['coulomb14scale'])) > NonbondedGenerator.SCALETOL or \
-                    abs(generator.lj14scale - float(element.attrib['lj14scale'])) > NonbondedGenerator.SCALETOL:
+            if (abs(generator.coulomb14scale - float(element.attrib['coulomb14scale'])) > NonbondedGenerator.SCALETOL
+                or abs(generator.lj14scale - float(element.attrib['lj14scale'])) > NonbondedGenerator.SCALETOL
+            ):
                 raise ValueError('Found multiple NonbondedForce tags with different 1-4 scales')
+            if (
+                    generator.useDispersionCorrection is not None
+                    and useDispersionCorrection is not None
+                    and generator.useDispersionCorrection != useDispersionCorrection
+            ):
+                raise ValueError('Found multiple NonbondedForce tags with different useDispersionCorrection settings.')
         generator.params.parseDefinitions(element)
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
@@ -2317,7 +2335,18 @@ class NonbondedGenerator(object):
         if 'ewaldErrorTolerance' in args:
             force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
         if 'useDispersionCorrection' in args:
-            force.setUseDispersionCorrection(bool(args['useDispersionCorrection']))
+            lrc_keyword = bool(args['useDispersionCorrection'])
+            if self.useDispersionCorrection is not None and lrc_keyword != self.useDispersionCorrection:
+                warnings.warn(
+                    "Conflicting settings for useDispersionCorrection in createSystem() and forcefield file. "
+                    "Using the one specified in createSystem()."
+                )
+            force.setUseDispersionCorrection(lrc_keyword)
+        elif self.useDispersionCorrection is not None:
+            force.setUseDispersionCorrection(self.useDispersionCorrection)
+        else:
+            # by default
+            force.setUseDispersionCorrection(True)
         sys.addForce(force)
 
     def postprocessSystem(self, sys, data, args):
@@ -2334,10 +2363,11 @@ parsers["NonbondedForce"] = NonbondedGenerator.parseElement
 class LennardJonesGenerator(object):
     """A NBFix generator to construct the L-J force with NBFIX implemented as a lookup table"""
 
-    def __init__(self, forcefield, lj14scale):
+    def __init__(self, forcefield, lj14scale, useDispersionCorrection):
         self.ff = forcefield
         self.nbfixTypes = {}
         self.lj14scale = lj14scale
+        self.useDispersionCorrection = useDispersionCorrection
         self.ljTypes = ForceField._AtomTypeParameters(forcefield, 'LennardJonesForce', 'Atom', ('sigma', 'epsilon'))
 
     def registerNBFIX(self, parameters):
@@ -2356,14 +2386,28 @@ class LennardJonesGenerator(object):
     @staticmethod
     def parseElement(element, ff):
         existing = [f for f in ff._forces if isinstance(f, LennardJonesGenerator)]
+        if 'useDispersionCorrection' in element.attrib:
+            useDispersionCorrection = bool(eval(element.attrib.get('useDispersionCorrection')))
+        else:
+            useDispersionCorrection = None
         if len(existing) == 0:
-            generator = LennardJonesGenerator(ff, float(element.attrib['lj14scale']))
+            generator = LennardJonesGenerator(
+                ff,
+                float(element.attrib['lj14scale']),
+                useDispersionCorrection=useDispersionCorrection
+            )
             ff.registerGenerator(generator)
         else:
             # Multiple <LennardJonesForce> tags were found, probably in different files
             generator = existing[0]
             if abs(generator.lj14scale - float(element.attrib['lj14scale'])) > NonbondedGenerator.SCALETOL:
                 raise ValueError('Found multiple LennardJonesForce tags with different 1-4 scales')
+            if (
+                    generator.useDispersionCorrection is not None
+                    and useDispersionCorrection is not None
+                    and generator.useDispersionCorrection != useDispersionCorrection
+            ):
+                raise ValueError('Found multiple LennardJonesForce tags with different useDispersionCorrection settings.')
         for LJ in element.findall('Atom'):
             generator.registerLennardJones(LJ.attrib)
         for Nbfix in element.findall('NBFixPair'):
@@ -2434,12 +2478,24 @@ class LennardJonesGenerator(object):
         if args['switchDistance'] is not None:
             self.force.setUseSwitchingFunction(True)
             self.force.setSwitchingDistance(args['switchDistance'])
+        if 'useDispersionCorrection' in args:
+            lrc_keyword = bool(args['useDispersionCorrection'])
+            if self.useDispersionCorrection is not None and lrc_keyword != self.useDispersionCorrection:
+                warnings.warn(
+                    "Conflicting settings for useDispersionCorrection in createSystem() and forcefield file. "
+                    "Using the one specified in createSystem()."
+                )
+            self.force.setUseLongRangeCorrection(lrc_keyword)
+        elif self.useDispersionCorrection is not None:
+            self.force.setUseLongRangeCorrection(self.useDispersionCorrection)
+        else:
+            # by default
+            self.force.setUseLongRangeCorrection(True)
 
         # Add the particles
 
         for atom in data.atoms:
             self.force.addParticle((typeToMergedType[data.atomType[atom]],))
-        self.force.setUseLongRangeCorrection(bool(args.get('useLongRangeCorrection', False)))
         self.force.setCutoffDistance(nonbondedCutoff)
         sys.addForce(self.force)
 

wrappers/python/tests/TestForceField.py

@@ -6,6 +6,7 @@ from simtk.unit import *
 import simtk.openmm.app.element as elem
 import simtk.openmm.app.forcefield as forcefield
 import math
+import textwrap
 try:
     from cStringIO import StringIO
 except ImportError:
@@ -50,16 +51,77 @@ class TestForceField(unittest.TestCase):
                                 for f in forces))
 
     def test_DispersionCorrection(self):
-        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
-
-        for useDispersionCorrection in [True, False]:
-            system = self.forcefield1.createSystem(self.pdb1.topology,
-                                                   nonbondedCutoff=2*nanometer,
-                                                   useDispersionCorrection=useDispersionCorrection)
+        """Test to make sure that the dispersion/long-range correction is set properly."""
+        top = Topology()
+        chain = top.addChain()
 
+        for lrc in (True, False):
+            xml = textwrap.dedent(
+                """
+                <ForceField>
+                 <LennardJonesForce lj14scale="0.3" useDispersionCorrection="{lrc}">
+                  <Atom type="A" sigma="1" epsilon="0.1"/>
+                  <Atom type="B" sigma="2" epsilon="0.2"/>
+                  <NBFixPair type1="A" type2="B" sigma="2.5" epsilon="1.1"/>
+                 </LennardJonesForce>
+                 <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5" useDispersionCorrection="{lrc2}">
+                  <Atom type="A" sigma="0.315" epsilon="0.635"/>
+                 </NonbondedForce>
+                </ForceField>
+                """
+            )
+            ff = ForceField(StringIO(xml.format(lrc=lrc, lrc2=lrc)))
+            system = ff.createSystem(top)
+            checked_nonbonded = False
+            checked_custom = False
             for force in system.getForces():
                 if isinstance(force, NonbondedForce):
-                    self.assertEqual(useDispersionCorrection, force.getUseDispersionCorrection())
+                    self.assertEqual(force.getUseDispersionCorrection(), lrc)
+                    checked_nonbonded = True
+                elif isinstance(force, CustomNonbondedForce):
+                    self.assertEqual(force.getUseLongRangeCorrection(), lrc)
+                    checked_custom = True
+            self.assertTrue(checked_nonbonded and checked_custom)
+
+            # check that the keyword argument overwrites xml input
+            lrc_kwarg = not lrc
+            with warnings.catch_warnings(record=True) as w:
+                warnings.simplefilter("always")
+                system2 = ff.createSystem(top, useDispersionCorrection=lrc_kwarg)
+                self.assertTrue(len(w) == 2)
+                assert "conflict" in str(w[-1].message).lower()
+            checked_nonbonded = False
+            checked_custom = False
+            for force in system2.getForces():
+                if isinstance(force, NonbondedForce):
+                    self.assertEqual(force.getUseDispersionCorrection(), lrc_kwarg)
+                    checked_nonbonded = True
+                elif isinstance(force, CustomNonbondedForce):
+                    self.assertEqual(force.getUseLongRangeCorrection(), lrc_kwarg)
+                    checked_custom = True
+            self.assertTrue(checked_nonbonded and checked_custom)
+
+            # check that no warning is generated when useDispersionCorrection is not in the xml file
+            xml = textwrap.dedent(
+                """
+                <ForceField>
+                 <LennardJonesForce lj14scale="0.3">
+                  <Atom type="A" sigma="1" epsilon="0.1"/>
+                  <Atom type="B" sigma="2" epsilon="0.2"/>
+                  <NBFixPair type1="A" type2="B" sigma="2.5" epsilon="1.1"/>
+                 </LennardJonesForce>
+                 <NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
+                  <Atom type="A" sigma="0.315" epsilon="0.635"/>
+                 </NonbondedForce>
+                </ForceField>
+                """
+            )
+            ff = ForceField(StringIO(xml))
+            system = ff.createSystem(top)
+            for lrc_kwarg in [True, False]:
+                with warnings.catch_warnings():
+                    warnings.simplefilter("error")
+                    system2 = ff.createSystem(top, useDispersionCorrection=lrc_kwarg)
 
     def test_Cutoff(self):
         """Test to make sure the nonbondedCutoff parameter is passed correctly."""
@@ -227,7 +289,7 @@ class TestForceField(unittest.TestCase):
         angles = forcefield.HarmonicAngleGenerator(ff)
         angles.registerAngle({'class1':'HW', 'class2':'OW', 'class3':'HW', 'angle':1.82421813418*radians, 'k':836.8*kilojoules_per_mole/radian})
         ff.registerGenerator(angles)
-        nonbonded = forcefield.NonbondedGenerator(ff, 0.833333, 0.5)
+        nonbonded = forcefield.NonbondedGenerator(ff, 0.833333, 0.5, True)
         nonbonded.registerAtom({'type':'tip3p-O', 'charge':-0.834, 'sigma':0.31507524065751241*nanometers, 'epsilon':0.635968*kilojoules_per_mole})
         nonbonded.registerAtom({'type':'tip3p-H', 'charge':0.417, 'sigma':1*nanometers, 'epsilon':0*kilojoules_per_mole})
         ff.registerGenerator(nonbonded)
@@ -712,7 +774,7 @@ class TestForceField(unittest.TestCase):
    <Atom name="SOD" type="SOD"/>
   </Residue>
  </Residues>
- <LennardJonesForce lj14scale="1.0">
+ <LennardJonesForce lj14scale="1.0" useDispersionCorrection="False">
   <Atom type="CLA" sigma="0.404468018036" epsilon="0.6276"/>
   <Atom type="SOD" sigma="0.251367073323" epsilon="0.1962296"/>
   <NBFixPair type1="CLA" type2="SOD" sigma="0.33239431" epsilon="0.350933"/>
@@ -940,5 +1002,6 @@ class AmoebaTestForceField(unittest.TestCase):
             diff = norm(f1-f2)
             self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-3)
 
+
 if __name__ == '__main__':
     unittest.main()