Add support for Foyer forcefields

devtools/conda-envs/test_env.yaml

@@ -15,6 +15,7 @@ dependencies:
   - requests-mock # For testing http requests.
   # smirnoff-plugins =0.0.3
   - coverage >=4.4
+  - foyer
 
     # Shims
   - pint =0.20.1

openff/evaluator/forcefield/__init__.py

@@ -1,5 +1,6 @@
 from .forcefield import (
     ForceFieldSource,
+    FoyerForceFieldSource,
     LigParGenForceFieldSource,
     SmirnoffForceFieldSource,
     TLeapForceFieldSource,
@@ -11,6 +12,7 @@
     SmirnoffForceFieldSource,
     LigParGenForceFieldSource,
     TLeapForceFieldSource,
+    FoyerForceFieldSource,
     ParameterGradient,
     ParameterGradientKey,
 ]

openff/evaluator/forcefield/forcefield.py

@@ -261,3 +261,43 @@ def __setstate__(self, state):
         self._cutoff = state["cutoff"]
         self._request_url = state["request_url"]
         self._download_url = state["download_url"]
+
+
+class FoyerForceFieldSource(ForceFieldSource):
+    """A wrapper around Foyer force fields"""
+
+    @property
+    def foyer_source(self):
+        """str: Foyer force field source."""
+        return self._foyer_source
+
+    @property
+    def cutoff(self):
+        """openff.evaluator.unit.Quantity: The non-bonded interaction cutoff."""
+        return self._cutoff
+
+    def __init__(self, foyer_source="", cutoff=0.9 * unit.nanometer):
+        """Constructs a new FoyerForceField Source
+
+        Parameters
+        ----------
+        foyer_source: str
+            'oplsaa' or a Foyer XML forcefield file
+        cutoff: openff.evaluator.unit.Quantity
+            The non-bonded interaction cutoff, default 0.9 nanometers.
+
+        Examples
+        --------
+        To create a source for the Foyer force field:
+
+        >>> foyer_source = FoyerForceFieldSource('oplsaa')
+        """
+        self._foyer_source = foyer_source
+        self._cutoff = cutoff
+
+    def __getstate__(self):
+        return {"foyer_source": self._foyer_source, "cutoff": self._cutoff}
+
+    def __setstate__(self, state):
+        self._foyer_source = state["foyer_source"]
+        self._cutoff = state["cutoff"]

openff/evaluator/protocols/forcefield.py

@@ -24,6 +24,7 @@
     SmirnoffForceFieldSource,
     TLeapForceFieldSource,
 )
+from openff.evaluator.forcefield.forcefield import FoyerForceFieldSource
 from openff.evaluator.forcefield.system import ParameterizedSystem
 from openff.evaluator.substances import Substance
 from openff.evaluator.utils.utils import (
@@ -87,6 +88,7 @@ def _append_system(existing_system, system_to_append, index_map=None):
             openmm.HarmonicAngleForce,
             openmm.PeriodicTorsionForce,
             openmm.NonbondedForce,
+            openmm.RBTorsionForce,
         ]
 
         number_of_appended_forces = 0
@@ -226,6 +228,15 @@ def _append_system(existing_system, system_to_append, index_map=None):
                         index_a + index_offset, index_b + index_offset, *parameters
                     )
 
+            elif isinstance(force_to_append, openmm.RBTorsionForce):
+                # Support for RBTorisionForce needed for OPLSAA, etc
+                for index in range(force_to_append.getNumTorsions()):
+                    torsion_params = force_to_append.getTorsionParameters(index)
+                    for i in range(4):
+                        torsion_params[i] = index_map[torsion_params[i]] + index_offset
+
+                    existing_force.addTorsion(*torsion_params)
+
             number_of_appended_forces += 1
 
         if number_of_appended_forces != system_to_append.getNumForces():
@@ -1052,3 +1063,57 @@ def _execute(self, directory, available_resources):
             )
 
         super(BuildTLeapSystem, self)._execute(directory, available_resources)
+
+
+@workflow_protocol()
+class BuildFoyerSystem(TemplateBuildSystem):
+    """Parameterize a set of molecules with a Foyer force field source"""
+
+    def _parameterize_molecule(self, molecule, force_field_source, cutoff):
+        """Parameterize the specified molecule.
+
+        Parameters
+        ----------
+        molecule: openff.toolkit.topology.Molecule
+            The molecule to parameterize.
+        force_field_source: FoyerForceFieldSource
+            The foyer source which describes which parameters to apply.
+
+        Returns
+        -------
+        openmm.System
+            The parameterized system.
+        """
+        from foyer import Forcefield as FoyerForceField
+        from openff.interchange import Interchange
+        from openff.toolkit import Topology
+
+        topology: Topology = molecule.to_topology()
+
+        force_field: FoyerForceField
+        if force_field_source.foyer_source.lower() == "oplsaa":
+            force_field = FoyerForceField(name="oplsaa")
+        else:
+            force_field = FoyerForceField(
+                forcefield_files=force_field_source.foyer_source
+            )
+
+        interchange = Interchange.from_foyer(topology=topology, force_field=force_field)
+
+        openmm_pdb_file = app.PDBFile(self.coordinate_file_path)
+        if openmm_pdb_file.topology.getPeriodicBoxVectors() is not None:
+            interchange.box = openmm_pdb_file.topology.getPeriodicBoxVectors()
+
+        openmm_system = interchange.to_openmm(combine_nonbonded_forces=False)
+
+        return openmm_system
+
+    def _execute(self, directory, available_resources):
+        force_field_source = ForceFieldSource.from_json(self.force_field_path)
+
+        if not isinstance(force_field_source, FoyerForceFieldSource):
+            raise ValueError(
+                "Only Foyer force field sources are supported by this protocol."
+            )
+
+        super(BuildFoyerSystem, self)._execute(directory, available_resources)

openff/evaluator/tests/test_protocols/test_forcefield.py

@@ -9,8 +9,10 @@
 from openff.toolkit.utils.rdkit_wrapper import RDKitToolkitWrapper
 
 from openff.evaluator.forcefield import LigParGenForceFieldSource, TLeapForceFieldSource
+from openff.evaluator.forcefield.forcefield import FoyerForceFieldSource
 from openff.evaluator.protocols.coordinates import BuildCoordinatesPackmol
 from openff.evaluator.protocols.forcefield import (
+    BuildFoyerSystem,
     BuildLigParGenSystem,
     BuildSmirnoffSystem,
     BuildTLeapSystem,
@@ -161,3 +163,69 @@ def download_callback(_, context):
         assign_parameters.substance = substance
         assign_parameters.execute(directory)
         assert path.isfile(assign_parameters.parameterized_system.system_path)
+
+
+def test_build_foyer_oplsaa_system():
+    force_field_source = FoyerForceFieldSource("oplsaa")
+    substance = Substance.from_components("C", "O", "CCC", "C1CCCC1")
+
+    with tempfile.TemporaryDirectory() as directory:
+        force_field_path = path.join(directory, "ff.json")
+
+        with open(force_field_path, "w") as file:
+            file.write(force_field_source.json())
+
+        build_coordinates = BuildCoordinatesPackmol("build_coordinates")
+        build_coordinates.max_molecules = 8
+        build_coordinates.substance = substance
+        build_coordinates.execute(directory)
+
+        assign_parameters = BuildFoyerSystem("assign_parameters")
+        assign_parameters.force_field_path = force_field_path
+        assign_parameters.coordinate_file_path = build_coordinates.coordinate_file_path
+        assign_parameters.substance = substance
+        assign_parameters.execute(directory)
+        assert path.isfile(assign_parameters.parameterized_system.system_path)
+
+
+def test_build_foyer_xml_system():
+    with tempfile.TemporaryDirectory() as directory:
+        force_field_source_path = path.join(directory, "ff.json")
+        force_field_xml_path = path.join(directory, "foyer_ff.xml")
+        force_field_source = FoyerForceFieldSource(force_field_xml_path)
+        substance = Substance.from_components("O")
+
+        with open(force_field_source_path, "w") as file:
+            file.write(force_field_source.json())
+
+        with open(force_field_xml_path, "w") as file:
+            file.write(
+                """<ForceField name="tip3p" version="0.0.3" combining_rule="geometric">
+ <AtomTypes>
+  <Type name="HW_tip3p" class="HW_tip3p" element="H" mass="1.008"/>
+  <Type name="OW_tip3p" class="OW_tip3p" element="O" mass="15.9994"/>
+ </AtomTypes>
+ <NonbondedForce coulomb14scale="0.5" lj14scale="0.5">
+  <Atom type="HW_tip3p" charge="0.417" sigma="1.0" epsilon="0.0"/>
+  <Atom type="OW_tip3p" charge="-0.834" sigma="0.315007" epsilon="0.63681228"/>
+ </NonbondedForce>
+ <HarmonicBondForce>
+  <Bond class1="OW_tip3p" class2="HW_tip3p" length="0.09572" k="1884.06"/>
+ </HarmonicBondForce>
+ <HarmonicAngleForce>
+  <Angle class1="HW_tip3p" class2="OW_tip3p" class3="HW_tip3p" angle="1.82421813" k="230.274"/>
+ </HarmonicAngleForce>
+</ForceField>"""
+            )
+
+        build_coordinates = BuildCoordinatesPackmol("build_coordinates")
+        build_coordinates.max_molecules = 8
+        build_coordinates.substance = substance
+        build_coordinates.execute(directory)
+
+        assign_parameters = BuildFoyerSystem("assign_parameters")
+        assign_parameters.force_field_path = force_field_source_path
+        assign_parameters.coordinate_file_path = build_coordinates.coordinate_file_path
+        assign_parameters.substance = substance
+        assign_parameters.execute(directory)
+        assert path.isfile(assign_parameters.parameterized_system.system_path)