Add TIP4P-EW

devtools/conda-envs/test_env.yaml

@@ -7,5 +7,6 @@ dependencies:
   - openff-toolkit >=0.14.3
   - pytest
   - pytest-cov
+  - pytest-xdist
   - pandas
   - click

openforcefields/offxml/tip4p_ew-1.0.0.offxml

@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+<SMIRNOFF version="0.3" aromaticity_model="OEAroModel_MDL">
+    <vdW version="0.4" potential="Lennard-Jones-12-6" combining_rules="Lorentz-Berthelot" scale12="0.0" scale13="0.0" scale14="0.5" scale15="1.0" cutoff="9.0 * angstrom ** 1" switch_width="1.0 * angstrom ** 1" periodic_method="cutoff" nonperiodic_method="no-cutoff">
+        <Atom smirks="[#1]-[#8X2H2+0:1]-[#1]" epsilon="0.680946 * kilojoule_per_mole ** 1" id="n-tip4p-ew-O" sigma="0.316435 * nanometer ** 1"></Atom>
+        <Atom smirks="[#1:1]-[#8X2H2+0]-[#1]" epsilon="0.0 * kilocalorie_per_mole ** 1" id="n-tip4p-ew-H" sigma="1.0 * nanometer ** 1"></Atom>
+    </vdW>
+    <LibraryCharges version="0.3">
+        <LibraryCharge smirks="[#1]-[#8X2H2+0:1]-[#1]" charge1="0.0 * elementary_charge ** 1" id="q-tip4p-ew-O"></LibraryCharge>
+        <LibraryCharge smirks="[#1:1]-[#8X2H2+0]-[#1]" charge1="0.0 * elementary_charge ** 1" id="q-tip4p-ew-H"></LibraryCharge>
+    </LibraryCharges>
+    <Electrostatics version="0.4" scale12="0.0" scale13="0.0" scale14="0.8333333333" scale15="1.0" cutoff="9.0 * angstrom ** 1" switch_width="0.0 * angstrom ** 1" periodic_potential="Ewald3D-ConductingBoundary" nonperiodic_potential="Coulomb" exception_potential="Coulomb"></Electrostatics>
+    <Constraints version="0.3">
+        <Constraint smirks="[#1:1]-[#8X2H2+0:2]-[#1]" id="c-tip4p-ew-H-O" distance="0.09572 * nanometer ** 1"></Constraint>
+        <Constraint smirks="[#1:1]-[#8X2H2+0]-[#1:2]" id="c-tip4p-ew-H-O-H" distance="0.15139006545247014 * nanometer ** 1"></Constraint>
+    </Constraints>
+    <VirtualSites version="0.3" exclusion_policy="parents">
+        <VirtualSite smirks="[#1:2]-[#8X2H2+0:1]-[#1:3]" epsilon="0.0 * kilocalorie_per_mole ** 1" type="DivalentLonePair" match="once" distance="-0.0125 * nanometer ** 1" outOfPlaneAngle="0.0 * degree ** 1" inPlaneAngle="None" charge_increment1="0.0 * elementary_charge ** 1" charge_increment2="0.52422 * elementary_charge ** 1" charge_increment3="0.52422 * elementary_charge ** 1" sigma="1.0 * angstrom ** 1" name="EP"></VirtualSite>
+    </VirtualSites>
+</SMIRNOFF>

openforcefields/offxml/tip4p_ew.offxml

@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+<SMIRNOFF version="0.3" aromaticity_model="OEAroModel_MDL">
+    <vdW version="0.4" potential="Lennard-Jones-12-6" combining_rules="Lorentz-Berthelot" scale12="0.0" scale13="0.0" scale14="0.5" scale15="1.0" cutoff="9.0 * angstrom ** 1" switch_width="1.0 * angstrom ** 1" periodic_method="cutoff" nonperiodic_method="no-cutoff">
+        <Atom smirks="[#1]-[#8X2H2+0:1]-[#1]" epsilon="0.680946 * kilojoule_per_mole ** 1" id="n-tip4p-ew-O" sigma="0.316435 * nanometer ** 1"></Atom>
+        <Atom smirks="[#1:1]-[#8X2H2+0]-[#1]" epsilon="0.0 * kilocalorie_per_mole ** 1" id="n-tip4p-ew-H" sigma="1.0 * nanometer ** 1"></Atom>
+    </vdW>
+    <LibraryCharges version="0.3">
+        <LibraryCharge smirks="[#1]-[#8X2H2+0:1]-[#1]" charge1="0.0 * elementary_charge ** 1" id="q-tip4p-ew-O"></LibraryCharge>
+        <LibraryCharge smirks="[#1:1]-[#8X2H2+0]-[#1]" charge1="0.0 * elementary_charge ** 1" id="q-tip4p-ew-H"></LibraryCharge>
+    </LibraryCharges>
+    <Electrostatics version="0.4" scale12="0.0" scale13="0.0" scale14="0.8333333333" scale15="1.0" cutoff="9.0 * angstrom ** 1" switch_width="0.0 * angstrom ** 1" periodic_potential="Ewald3D-ConductingBoundary" nonperiodic_potential="Coulomb" exception_potential="Coulomb"></Electrostatics>
+    <Constraints version="0.3">
+        <Constraint smirks="[#1:1]-[#8X2H2+0:2]-[#1]" id="c-tip4p-ew-H-O" distance="0.09572 * nanometer ** 1"></Constraint>
+        <Constraint smirks="[#1:1]-[#8X2H2+0]-[#1:2]" id="c-tip4p-ew-H-O-H" distance="0.15139006545247014 * nanometer ** 1"></Constraint>
+    </Constraints>
+    <VirtualSites version="0.3" exclusion_policy="parents">
+        <VirtualSite smirks="[#1:2]-[#8X2H2+0:1]-[#1:3]" epsilon="0.0 * kilocalorie_per_mole ** 1" type="DivalentLonePair" match="once" distance="-0.0125 * nanometer ** 1" outOfPlaneAngle="0.0 * degree ** 1" inPlaneAngle="None" charge_increment1="0.0 * elementary_charge ** 1" charge_increment2="0.52422 * elementary_charge ** 1" charge_increment3="0.52422 * elementary_charge ** 1" sigma="1.0 * angstrom ** 1" name="EP"></VirtualSite>
+    </VirtualSites>
+</SMIRNOFF>

openforcefields/tests/test_water_models.py

@@ -171,28 +171,37 @@ def test_opc(water_molecule):
     )
 
 
-@pytest.mark.skip(reason="Skipping in first pass")
 def test_tip4p_ew(water_molecule):
-    interchange = ForceField("tip4p-ew_1.0.0.offxml").create_interchange(
-        water_molecule,
-    )
-
-    openmm_topology = to_openmm_topology(interchange)
-
-    reference = OpenMMForceField("tip4pew.xml").createSystem(
-        openmm_topology,
+    from openmm.app import Modeller
+    omm_water = water_molecule.to_openmm()
+    omm_ff = OpenMMForceField("tip4pew.xml")
+    mod = Modeller(omm_water, water_molecule.get_positions().to_openmm())
+    mod.addExtraParticles(omm_ff)
+    reference = omm_ff.createSystem(
+        mod.getTopology(),
         constraints=HAngles,
         rigidWater=True,
     )
 
+    interchange = ForceField("tip4p_ew-1.0.0.offxml").create_interchange(
+        water_molecule,
+    )
+    system = interchange.to_openmm()
+
     compare_water_systems(
-        reference, interchange.to_openmm(combine_nonbonded_forces=True)
+        reference,
+        system,
+        {
+            "charge": 1e-10 * openmm.unit.elementary_charge,
+            "sigma": 1e-10 * openmm.unit.nanometer,
+            "epsilon": 1e-10 * openmm.unit.kilojoule_per_mole,
+        },
     )
 
 
 @pytest.mark.skip(reason="Skipping in first pass")
 def test_tip5p(water_molecule):
-    interchange = ForceField("tip5p_1.0.0.offxml").create_openmm_system(
+    interchange = ForceField("tip5p-1.0.0.offxml").create_openmm_system(
         water_molecule,
     )
 
@@ -211,7 +220,14 @@ def test_tip5p(water_molecule):
 
 @pytest.mark.parametrize(
     "water_model,pattern",
-    [("tip3p", "^tip3p(?!.*fb)"), ("tip3p_fb", "^tip3p_fb"), ("tip4p_fb", "^tip4p_fb"), ("opc3", "^opc3"), ("opc", "^opc(?!3)")],
+    [
+        ("tip3p", "^tip3p(?!.*fb)"),
+        ("tip3p_fb", "^tip3p_fb"),
+        ("tip4p_fb", "^tip4p_fb"),
+        ("tip4p_ew", "^tip4p_ew"),
+        ("opc3", "^opc3"),
+        ("opc", "^opc(?!3)"),
+    ],
 )
 def test_most_recent_version_match(water_model, pattern):
     import re
@@ -320,12 +336,19 @@ def test_ion_parameter_assignment(water_molecule):
             parameter_was_used
         ), f"The ion LibraryCharge parameter with smirks {key} was not assigned"
 
-def test_water_model_is_compatible_with_mainline():
+@pytest.mark.parametrize(
+    "water_model",
+    [
+        "tip3p.offxml",
+        "tip3p_fb.offxml",
+        "tip4p_fb.offxml",
+        "tip4p_ew.offxml",
+        "opc3.offxml",
+        "opc.offxml",
+    ]
+)
+def test_water_model_is_compatible_with_mainline(water_model):
     """Ensure that the latest water model FF is compatible with the latest main-line FF"""
     # Since we don't have a way to get the most recent mainline FF, be sure
     # to occasionally update the first FF listed here
-    ForceField('openff-2.1.0.offxml', 'tip3p.offxml')
-    ForceField('openff-2.1.0.offxml', 'tip3p_fb.offxml')
-    ForceField('openff-2.1.0.offxml', 'tip4p_fb.offxml')
-    ForceField('openff-2.1.0.offxml', 'opc3.offxml')
-    ForceField('openff-2.1.0.offxml', 'opc.offxml')
+    ForceField('openff-2.1.0.offxml', water_model)

scripts/write_tip4p_ew.py

@@ -0,0 +1,105 @@
+"""
+Write TIP3P-EW parameters into a SMIRNOFF force field. Based on
+https://github.com/openmm/openmm/blob/116aed3927066b0a53eba929110d73f3daew64bd/wrappers/python/openmm/app/data/tip4pew.xml
+"""
+from pathlib import Path
+
+import pandas
+from openff.toolkit.typing.engines.smirnoff.forcefield import ForceField
+from openff.toolkit.typing.engines.smirnoff.parameters import (
+    ConstraintHandler,
+    ElectrostaticsHandler,
+    LibraryChargeHandler,
+    VirtualSiteHandler,
+    vdWHandler,
+)
+from openff.units import unit
+from packaging import version
+
+VERSION = version.Version("1.0.0")
+OFFXML_PATH = Path("openforcefields", "offxml")
+
+tip4p_ew = ForceField()
+
+tip4p_ew_electrostatics = ElectrostaticsHandler(version=0.4, scale14=0.8333333333)
+tip4p_ew_library = LibraryChargeHandler(version=0.3)
+tip4p_ew_vdw = vdWHandler(version=0.4)
+tip4p_ew_constraints = ConstraintHandler(version=0.3)
+tip4p_ew_virtual_sites = VirtualSiteHandler(version=0.3)
+
+tip4p_ew_vdw.add_parameter(
+    {
+        "smirks": "[#1]-[#8X2H2+0:1]-[#1]",
+        "epsilon": unit.Quantity(0.680946, unit.kilojoule_per_mole),
+        "sigma": unit.Quantity(0.316435, unit.nanometer),
+        "id": "n-tip4p-ew-O",
+    }
+)
+tip4p_ew_vdw.add_parameter(
+    {
+        "smirks": "[#1:1]-[#8X2H2+0]-[#1]",
+        "epsilon": unit.Quantity(0.0, unit.kilocalorie_per_mole),
+        "sigma": unit.Quantity(1.0, unit.nanometer),
+        "id": "n-tip4p-ew-H",
+    }
+)
+
+tip4p_ew_library.add_parameter(
+    {
+        "smirks": "[#1]-[#8X2H2+0:1]-[#1]",
+        "charge1": unit.Quantity(0.0, unit.elementary_charge),
+        "id": "q-tip4p-ew-O",
+    }
+)
+tip4p_ew_library.add_parameter(
+    {
+        "smirks": "[#1:1]-[#8X2H2+0]-[#1]",
+        "charge1": unit.Quantity(0.0, unit.elementary_charge),
+        "id": "q-tip4p-ew-H",
+    }
+)
+
+# Virtual site distance = 0.08984267127345 * 2 * 0.09572 * cos(1.82421813418 / 2)
+tip4p_ew_virtual_sites.add_parameter(
+    {
+        "type": "DivalentLonePair",
+        "smirks": "[#1:2]-[#8X2H2+0:1]-[#1:3]",
+        "match": "once",
+        "name": "EP",
+        "distance": unit.Quantity(-0.0125, unit.nanometer),
+        "sigma": unit.Quantity(1.0, unit.angstrom),
+        "epsilon": unit.Quantity(0.0, unit.kilocalorie_per_mole),
+        "outOfPlaneAngle": unit.Quantity(0.0, unit.degree),
+        "charge_increment1": unit.Quantity(0.0, unit.elementary_charge),
+        "charge_increment2": unit.Quantity(0.52422, unit.elementary_charge),
+        "charge_increment3": unit.Quantity(0.52422, unit.elementary_charge),
+    }
+)
+
+tip4p_ew_constraints.add_parameter(
+    {
+        "smirks": "[#1:1]-[#8X2H2+0:2]-[#1]",
+        "id": "c-tip4p-ew-H-O",
+        "distance": unit.Quantity(0.09572, unit.nanometer),
+    }
+)
+# H-H distance = 2 * 0.09572 * sin(1.82421813418 / 2)
+tip4p_ew_constraints.add_parameter(
+    {
+        "smirks": "[#1:1]-[#8X2H2+0]-[#1:2]",
+        "id": "c-tip4p-ew-H-O-H",
+        "distance": unit.Quantity(0.15139006545247014, unit.nanometer),
+    }
+)
+
+for handler in [
+    tip4p_ew_vdw,
+    tip4p_ew_library,
+    tip4p_ew_electrostatics,
+    tip4p_ew_constraints,
+    tip4p_ew_virtual_sites,
+]:
+    tip4p_ew.register_parameter_handler(handler)
+
+tip4p_ew.to_file(Path(OFFXML_PATH, "tip4p_ew.offxml"))
+tip4p_ew.to_file(Path(OFFXML_PATH, f"tip4p_ew-{VERSION}.offxml"))