Merge pull request #1119 from openforcefield/write-itp

Optionally write GROMACS topologies as multiple `.itp` files

docs/using/output.md

@@ -34,6 +34,8 @@ interchange.to_gromacs("mysim")  # Produces the same three files
 
 Note that the MDP file generated is configured for a single-point energy calculation and must be modified to run other simulations.
 
+By default, the topology is written to as a monolithic file, which can be large. To split this into separate files with the `#include "molecule.itp"` convention, use `monolithic=False`. This produces a functionally equivalent topology file which splits non-bonded interactions into a file `mysim_nonbonded.itp` and each molecule's parameters in a separate file named according to the `Molecule.name` attribute.
+
 ## LAMMPS
 
 An [`Interchange`] object can be written to LAMMPS data and run input files with [`Interchange.to_lammps()`]

examples/protein_ligand/protein_ligand.ipynb

@@ -566,7 +566,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "system_intrcg.to_gromacs(prefix=\"out\")"
+    "system_intrcg.to_gromacs(prefix=\"out\", monolithic=False)"
    ]
   },
   {
@@ -703,7 +703,7 @@
  "metadata": {
   "category": "tutorial",
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -717,12 +717,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
-  },
-  "vscode": {
-   "interpreter": {
-    "hash": "86c9b142c8dc60dd36d17e2a57efabbd2ed015b9d3db80dd77f3e0894d5aea85"
-   }
+   "version": "3.11.11"
   }
  },
  "nbformat": 4,

openff/interchange/_tests/_gromacs.py

@@ -0,0 +1,16 @@
+"""Helpers for testing GROMACS interoperability."""
+
+from openff.utilities import temporary_cd
+
+from openff.interchange import Interchange
+from openff.interchange.drivers import get_gromacs_energies
+
+
+def gmx_monolithic_vs_itp(state: Interchange):
+    with temporary_cd():
+        get_gromacs_energies(state, _monolithic=True).compare(get_gromacs_energies(state, _monolithic=False))
+
+        # TODO: More helpful handling of failures, i.e.
+        #         * Detect differences in positions
+        #         * Detect differences in box vectors
+        #         * Detect differences in non-bonded settings

openff/interchange/_tests/conftest.py

@@ -287,6 +287,21 @@ def gbsa_force_field() -> ForceField:
     )
 
 
+@pytest.fixture
+def ff14sb() -> ForceField:
+    return ForceField("ff14sb_off_impropers_0.0.3.offxml")
+
+
+@pytest.fixture
+def ligand():
+    return MoleculeWithConformer.from_smiles("CC[C@@](/C=C\\[H])(C=C)O", allow_undefined_stereo=True)
+
+
+@pytest.fixture
+def caffeine():
+    return MoleculeWithConformer.from_smiles("Cn1cnc2c1c(=O)n(C)c(=O)n2C")
+
+
 @pytest.fixture
 def basic_top() -> Topology:
     topology = MoleculeWithConformer.from_smiles("C").to_topology()
@@ -469,27 +484,33 @@ def ethanol_top(ethanol):
 
 
 @pytest.fixture
-def mainchain_ala():
-    molecule = Molecule.from_file(
-        get_data_file_path("proteins/MainChain_ALA.sdf", "openff.toolkit"),
+def alanine_dipeptide() -> Topology:
+    return Topology.from_pdb(
+        get_data_file_path(
+            "proteins/MainChain_ALA_ALA.pdb",
+            "openff.toolkit",
+        ),
     )
-    molecule._add_default_hierarchy_schemes()
-    molecule.perceive_residues()
-    molecule.perceive_hierarchy()
-
-    return molecule
 
 
 @pytest.fixture
-def mainchain_arg():
-    molecule = Molecule.from_file(
-        get_data_file_path("proteins/MainChain_ARG.sdf", "openff.toolkit"),
-    )
-    molecule._add_default_hierarchy_schemes()
-    molecule.perceive_residues()
-    molecule.perceive_hierarchy()
+def mainchain_ala() -> Molecule:
+    return Topology.from_pdb(
+        get_data_file_path(
+            "proteins/MainChain_ALA.pdb",
+            "openff.toolkit",
+        ),
+    ).molecule(0)
 
-    return molecule
+
+@pytest.fixture
+def mainchain_arg() -> Molecule:
+    return Topology.from_pdb(
+        get_data_file_path(
+            "proteins/MainChain_ARG.pdb",
+            "openff.toolkit",
+        ),
+    ).molecule(0)
 
 
 @pytest.fixture

openff/interchange/_tests/interoperability_tests/gromacs/test_systems.py

@@ -0,0 +1,24 @@
+from openff.toolkit import Quantity
+
+from openff.interchange._tests._gromacs import gmx_monolithic_vs_itp
+
+
+def test_ligand_vacuum(caffeine, sage_unconstrained, monkeypatch):
+    monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+    topology = caffeine.to_topology()
+    topology.box_vectors = Quantity([4, 4, 4], "nanometer")
+
+    gmx_monolithic_vs_itp(sage_unconstrained.create_interchange(topology))
+
+
+def test_water_dimer(water_dimer, tip3p, monkeypatch):
+    monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+    gmx_monolithic_vs_itp(tip3p.create_interchange(water_dimer))
+
+
+def test_alanine_dipeptide(alanine_dipeptide, ff14sb, monkeypatch):
+    monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+    gmx_monolithic_vs_itp(ff14sb.create_interchange(alanine_dipeptide))

openff/interchange/components/interchange.py

@@ -1,5 +1,6 @@
 """An object for storing, manipulating, and converting molecular mechanics data."""
 
+import tempfile
 import warnings
 from collections.abc import Iterable
 from pathlib import Path
@@ -323,6 +324,7 @@ def to_gromacs(
         prefix: str,
         decimal: int = 3,
         hydrogen_mass: PositiveFloat = 1.007947,
+        monolithic: bool = True,
         _merge_atom_types: bool = False,
     ):
         """
@@ -339,6 +341,9 @@ def to_gromacs(
             The mass to use for hydrogen atoms if not present in the topology. If non-trivially different
             than the default value, mass will be transferred from neighboring heavy atoms. Note that this is currently
             not applied to any waters and is unsupported when virtual sites are present.
+        monolithic: bool, default=False
+            Whether the topology file should be monolithic (True) or reference individual .itp files (False). Note that
+            these individual .itp files rely on ad hoc atom types and cannot be transferred between systems.
         _merge_atom_types: bool, default = False
             The flag to define behaviour of GROMACSWriter. If True, then similar atom types will be merged.
             If False, each atom will have its own atom type.
@@ -360,7 +365,7 @@ def to_gromacs(
             gro_file=prefix + ".gro",
         )
 
-        writer.to_top(_merge_atom_types=_merge_atom_types)
+        writer.to_top(monolithic=monolithic, _merge_atom_types=_merge_atom_types)
         writer.to_gro(decimal=decimal)
 
         self.to_mdp(prefix + "_pointenergy.mdp")
@@ -394,6 +399,7 @@ def to_top(
         self,
         file_path: Path | str,
         hydrogen_mass: PositiveFloat = 1.007947,
+        monolithic: bool = True,
         _merge_atom_types: bool = False,
     ):
         """
@@ -407,6 +413,9 @@ def to_top(
             The mass to use for hydrogen atoms if not present in the topology. If non-trivially different
             than the default value, mass will be transferred from neighboring heavy atoms. Note that this is currently
             not applied to any waters and is unsupported when virtual sites are present.
+        monolithic: bool, default=False
+            Whether the topology file should be monolithic (True) or reference individual .itp files (False). Note that
+            these individual .itp files rely on ad hoc atom types and cannot be transferred between systems.
         _merge_atom_types: book, default=False
             The flag to define behaviour of GROMACSWriter. If True, then similar atom types will be merged.
             If False, each atom will have its own atom type.
@@ -423,7 +432,11 @@ def to_top(
         GROMACSWriter(
             system=_convert(self, hydrogen_mass=hydrogen_mass),
             top_file=file_path,
-        ).to_top(_merge_atom_types=_merge_atom_types)
+            gro_file=tempfile.NamedTemporaryFile(suffix=".gro").file.name,
+        ).to_top(
+            monolithic=monolithic,
+            _merge_atom_types=_merge_atom_types,
+        )
 
     def to_gro(self, file_path: Path | str, decimal: int = 3):
         """
@@ -456,6 +469,7 @@ def to_gro(self, file_path: Path | str, decimal: int = 3):
         # TODO: Write the coordinates without the full conversion
         GROMACSWriter(
             system=_convert(self),
+            top_file=tempfile.NamedTemporaryFile(suffix=".top").file.name,
             gro_file=file_path,
         ).to_gro(decimal=decimal)
 

openff/interchange/drivers/gromacs.py

@@ -52,6 +52,7 @@ def get_gromacs_energies(
     round_positions: int = 8,
     detailed: bool = False,
     _merge_atom_types: bool = False,
+    _monolithic: bool = True,
 ) -> EnergyReport:
     """
     Given an OpenFF Interchange object, return single-point energies as computed by GROMACS.
@@ -83,6 +84,7 @@ def get_gromacs_energies(
             mdp=mdp,
             round_positions=round_positions,
             merge_atom_types=_merge_atom_types,
+            monolithic=_monolithic,
         ),
         detailed=detailed,
     )
@@ -93,13 +95,15 @@ def _get_gromacs_energies(
     mdp: str = "auto",
     round_positions: int = 8,
     merge_atom_types: bool = False,
+    monolithic: bool = True,
 ) -> dict[str, Quantity]:
     with tempfile.TemporaryDirectory() as tmpdir:
         with temporary_cd(tmpdir):
             prefix = "_tmp"
             interchange.to_gromacs(
                 prefix=prefix,
                 decimal=round_positions,
+                monolithic=monolithic,
                 _merge_atom_types=merge_atom_types,
             )
 

openff/interchange/interop/gromacs/export/_export.py

@@ -1,5 +1,6 @@
 import pathlib
 import warnings
+from typing import IO
 
 import numpy
 from openff.toolkit import unit
@@ -23,30 +24,42 @@ class GROMACSWriter(_BaseModel):
     """Thin wrapper for writing GROMACS systems."""
 
     system: GROMACSSystem
-    top_file: pathlib.Path | str | None = None
-    gro_file: pathlib.Path | str | None = None
+    top_file: pathlib.Path | str
+    gro_file: pathlib.Path | str
 
-    def to_top(self, _merge_atom_types: bool = False):
+    def to_top(self, monolithic: bool = True, _merge_atom_types: bool = False):
         """Write a GROMACS topology file."""
-        if self.top_file is None:
-            raise ValueError("No TOP file specified.")
-
         with open(self.top_file, "w") as top:
             self._write_defaults(top)
+
+            if monolithic:
+                atomtypes_file_object: IO[str] = top
+            else:
+                prefix = str(self.top_file).split(".top")[0]
+                atomtypes_file_object = open(
+                    f"{prefix}_atomtypes.itp",
+                    "w",
+                )
+                top.write(f'#include "{prefix}_atomtypes.itp"\n')
+
             mapping_to_reduced_atom_types = self._write_atomtypes(
-                top,
-                _merge_atom_types,
+                top=atomtypes_file_object,
+                merge_atom_types=_merge_atom_types,
             )
 
             self._write_moleculetypes(
-                top,
-                mapping_to_reduced_atom_types,
-                _merge_atom_types,
+                top=top,
+                monolithic=monolithic,
+                mapping_to_reduced_atom_types=mapping_to_reduced_atom_types,
+                merge_atom_types=_merge_atom_types,
             )
 
             self._write_system(top)
             self._write_molecules(top)
 
+            if not monolithic:
+                atomtypes_file_object.close()
+
     def to_gro(self, decimal: int = 3):
         """Write a GROMACS coordinate file."""
         if self.gro_file is None:
@@ -55,7 +68,7 @@ def to_gro(self, decimal: int = 3):
         with open(self.gro_file, "w") as gro:
             self._write_gro(gro, decimal)
 
-    def _write_defaults(self, top):
+    def _write_defaults(self, top: IO[str]):
         top.write("[ defaults ]\n")
         top.write("; nbfunc\tcomb-rule\tgen-pairs\tfudgeLJ\tfudgeQQ\n")
 
@@ -68,7 +81,7 @@ def _write_defaults(self, top):
             f"{self.system.coul_14:8.6f}\n\n",
         )
 
-    def _write_atomtypes(self, top, merge_atom_types: bool) -> dict[str, bool | str]:
+    def _write_atomtypes(self, top: IO[str], merge_atom_types: bool) -> dict[str, bool | str]:
         top.write("[ atomtypes ]\n")
         top.write(
             ";type, bondingtype, atomic_number, mass, charge, ptype, sigma, epsilon\n",
@@ -156,30 +169,46 @@ def _get_new_entry_name(atom_type_list) -> str:
 
     def _write_moleculetypes(
         self,
-        top,
+        top: IO[str],
+        monolithic: bool,
         mapping_to_reduced_atom_types,
         merge_atom_types: bool,
     ):
         for molecule_name, molecule_type in self.system.molecule_types.items():
-            top.write("[ moleculetype ]\n")
+            # this string needs to be something that plays nicely in file paths
+            # and also works as GROMACS's label for the moleculetype "name"
+            canonicalized_name = molecule_name.translate({ord(c): "_" for c in r' \/:<>"|?*'})
 
-            top.write(
-                f"{molecule_name.replace(' ', '_')}\t{molecule_type.nrexcl:10d}\n\n",
-            )
+            if monolithic:
+                molecule_file = top
+            else:
+                prefix = str(self.top_file).split(".top")[0]
+                molecule_file = open(
+                    file=f"{prefix}_{canonicalized_name}.itp",
+                    mode="w",
+                )
+                top.write(f'#include "{prefix}_{canonicalized_name}.itp"\n')
+
+            molecule_file.write("[ moleculetype ]\n")
+
+            molecule_file.write(f"{canonicalized_name}\t{molecule_type.nrexcl:10d}\n\n")
 
             self._write_atoms(
-                top,
+                molecule_file,
                 molecule_type,
                 mapping_to_reduced_atom_types,
                 merge_atom_types,
             )
-            self._write_pairs(top, molecule_type)
-            self._write_bonds(top, molecule_type)
-            self._write_angles(top, molecule_type)
-            self._write_dihedrals(top, molecule_type)
-            self._write_settles(top, molecule_type)
-            self._write_virtual_sites(top, molecule_type)
-            self._write_exclusions(top, molecule_type)
+            self._write_pairs(molecule_file, molecule_type)
+            self._write_bonds(molecule_file, molecule_type)
+            self._write_angles(molecule_file, molecule_type)
+            self._write_dihedrals(molecule_file, molecule_type)
+            self._write_settles(molecule_file, molecule_type)
+            self._write_virtual_sites(molecule_file, molecule_type)
+            self._write_exclusions(molecule_file, molecule_type)
+
+            if not monolithic:
+                molecule_file.close()
 
         top.write("\n")