Merge pull request #973 from openforcefield/gmx-merge-atomtypes

Upstream GROMACS atom type merging

openff/interchange/_tests/unit_tests/interop/gromacs/export/test_export.py

@@ -193,8 +193,6 @@ def test_atom_names_pdb(self):
 @needs_gmx
 class TestGROMACS:
     @pytest.mark.slow
-    @pytest.mark.skip("from_top is not yet refactored for new Topology API")
-    @pytest.mark.parametrize("reader", ["intermol", "internal"])
     @pytest.mark.parametrize(
         "smiles",
         [
@@ -207,21 +205,27 @@ class TestGROMACS:
             # "C1COC(=O)O1",  # This adds an improper, i2
         ],
     )
-    def test_simple_roundtrip(self, sage, smiles, reader):
+    def test_simple_roundtrip(self, sage_unconstrained, smiles, monkeypatch):
+        # Skip if using RDKit conformer, which does weird stuff around 180 deg
+        pytest.importorskip("openeye.oechem")
+
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
         molecule = MoleculeWithConformer.from_smiles(smiles)
         molecule.name = molecule.to_hill_formula()
         topology = molecule.to_topology()
 
-        out = Interchange.from_smirnoff(force_field=sage, topology=topology)
+        out = sage_unconstrained.create_interchange(topology)
+
         out.box = [4, 4, 4]
-        out.positions = molecule.conformers[0]
+        out.positions = numpy.round(molecule.conformers[0], 2)
 
         out.to_top("out.top")
-        out.to_gro("out.gro")
+        out.to_gro("out.gro", decimal=3)
 
-        converted = Interchange.from_gromacs("out.top", "out.gro", reader=reader)
+        converted = Interchange.from_gromacs("out.top", "out.gro")
 
-        assert numpy.allclose(out.positions, converted.positions)
+        assert numpy.allclose(out.positions, converted.positions, atol=1e-3)
         assert numpy.allclose(out.box, converted.box)
 
         get_gromacs_energies(out).compare(
@@ -496,6 +500,121 @@ def test_common_boxes(self, pdb_file):
         )
 
 
+class TestMergeAtomTypes:
+    @pytest.mark.slow
+    @pytest.mark.parametrize(
+        "smiles",
+        [
+            "CC",  # Two identical carbons
+            "C1CCCCC1",  # Identical carbons in the ring
+            "C1[C@@H](N)C[C@@H](N)CC1",  # Identical carbons and nitrogens in the ring
+        ],
+    )
+    def test_energies_with_merging_atom_types(self, sage, smiles):
+        """
+        Tests #962
+        """
+        molecule = MoleculeWithConformer.from_smiles(smiles)
+        molecule.name = molecule.to_hill_formula()
+        topology = molecule.to_topology()
+
+        out = Interchange.from_smirnoff(force_field=sage, topology=topology)
+        out.box = [4, 4, 4]
+        out.positions = molecule.conformers[0]
+
+        get_gromacs_energies(out).compare(
+            get_gromacs_energies(out, _merge_atom_types=True),
+            tolerances={
+                "Bond": 0.002 * molecule.n_bonds * unit.kilojoule / unit.mol,
+                "Electrostatics": 0.05 * unit.kilojoule / unit.mol,
+            },
+        )
+
+    @pytest.mark.slow
+    @pytest.mark.parametrize(
+        "smiles",
+        [
+            "CC",  # Two identical carbons
+            "C1CCCCC1",  # Identical carbons in the ring
+            "C1[C@@H](N)C[C@@H](N)CC1",  # Identical carbons and nitrogens in the ring
+        ],
+    )
+    def test_simple_roundtrip_with_merging_atom_types(
+        self,
+        sage_unconstrained,
+        smiles,
+        monkeypatch,
+    ):
+        """
+        Tests #962
+        """
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+        molecule = MoleculeWithConformer.from_smiles(smiles)
+        molecule.name = molecule.to_hill_formula()
+        topology = molecule.to_topology()
+
+        out = sage_unconstrained.create_interchange(topology)
+        out.box = [4, 4, 4]
+        out.positions = molecule.conformers[0]
+
+        out.to_top("out.top")
+        out.to_top("out_merged.top", _merge_atom_types=True)
+        out.to_gro("out.gro", decimal=3)
+
+        converted = Interchange.from_gromacs("out.top", "out.gro")
+        converted_merged = Interchange.from_gromacs(
+            "out_merged.top",
+            "out.gro",
+        )
+
+        assert numpy.allclose(converted.positions, converted_merged.positions)
+        assert numpy.allclose(converted.box, converted_merged.box)
+
+        get_gromacs_energies(converted_merged).compare(
+            get_gromacs_energies(converted),
+            tolerances={
+                "Bond": 0.002 * molecule.n_bonds * unit.kilojoule / unit.mol,
+                "Electrostatics": 0.05 * unit.kilojoule / unit.mol,
+            },
+        )
+
+    @pytest.mark.parametrize(
+        "molecule_list",
+        [
+            ["CC", "CCO"],
+            ["CC", "CCCC"],
+            ["c1ccccc1", "CC", "CCO", "O"],
+            ["c1ccncc1", "n1cnccc1", "[nH]1cccc1"],
+        ],
+    )
+    def test_merge_atom_types_of_similar_molecules(
+        self,
+        molecule_list,
+        sage,
+    ):
+        pytest.importorskip("parmed")
+
+        topology = Topology.from_molecules(
+            [Molecule.from_smiles(smi) for smi in molecule_list],
+        )
+
+        out = sage.create_interchange(topology)
+
+        for merge in [True, False]:
+            out.to_top(f"{merge}.top", _merge_atom_types=merge)
+
+        not_merged = parmed.load_file("False.top")
+        merged = parmed.load_file("True.top")
+
+        # These are (n_atoms x 1) lists of parameters, which should be
+        # read from [ atoms ] section and cross-referenced to [ atomtypes ]
+        for attr in ["sigma", "epsilon", "charge"]:
+            assert [getattr(atom, attr) for atom in not_merged.atoms] == [
+                getattr(atom, attr) for atom in merged.atoms
+            ]
+
+
 @needs_gmx
 class TestGROMACSVirtualSites:
     @pytest.fixture

openff/interchange/components/interchange.py

@@ -421,6 +421,7 @@ def to_gromacs(
         prefix: str,
         decimal: int = 3,
         hydrogen_mass: float = 1.007947,
+        _merge_atom_types: bool = False,
     ):
         """
         Export this Interchange object to GROMACS files.
@@ -436,6 +437,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.
+        _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.
 
         """
         from openff.interchange.interop.gromacs.export._export import GROMACSWriter
@@ -447,13 +451,14 @@ def to_gromacs(
             gro_file=prefix + ".gro",
         )
 
-        writer.to_top()
+        writer.to_top(_merge_atom_types=_merge_atom_types)
         writer.to_gro(decimal=decimal)
 
     def to_top(
         self,
         file_path: Path | str,
         hydrogen_mass: float = 1.007947,
+        _merge_atom_types: bool = False,
     ):
         """
         Export this Interchange to a GROMACS topology file.
@@ -466,6 +471,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.
+        _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.
 
         """
         from openff.interchange.interop.gromacs.export._export import GROMACSWriter
@@ -474,7 +482,7 @@ def to_top(
         GROMACSWriter(
             system=_convert(self, hydrogen_mass=hydrogen_mass),
             top_file=file_path,
-        ).to_top()
+        ).to_top(_merge_atom_types=_merge_atom_types)
 
     def to_gro(self, file_path: Path | str, decimal: int = 3):
         """

openff/interchange/drivers/gromacs.py

@@ -51,6 +51,7 @@ def get_gromacs_energies(
     mdp: str = "auto",
     round_positions: int = 8,
     detailed: bool = False,
+    _merge_atom_types: bool = False,
 ) -> EnergyReport:
     """
     Given an OpenFF Interchange object, return single-point energies as computed by GROMACS.
@@ -67,6 +68,8 @@ def get_gromacs_energies(
         A decimal precision for the positions in the `.gro` file.
     detailed : bool, default=False
         If True, return a detailed report containing the energies of each term.
+    _merge_atom_types: bool, default=False
+        If True, energy should be computed with merging atom types.
 
     Returns
     -------
@@ -79,6 +82,7 @@ def get_gromacs_energies(
             interchange=interchange,
             mdp=mdp,
             round_positions=round_positions,
+            merge_atom_types=_merge_atom_types,
         ),
         detailed=detailed,
     )
@@ -88,11 +92,16 @@ def _get_gromacs_energies(
     interchange: Interchange,
     mdp: str = "auto",
     round_positions: int = 8,
+    merge_atom_types: bool = False,
 ) -> dict[str, unit.Quantity]:
     with tempfile.TemporaryDirectory() as tmpdir:
         with temporary_cd(tmpdir):
             prefix = "_tmp"
-            interchange.to_gromacs(prefix=prefix, decimal=round_positions)
+            interchange.to_gromacs(
+                prefix=prefix,
+                decimal=round_positions,
+                _merge_atom_types=merge_atom_types,
+            )
 
             if mdp == "auto":
                 mdconfig = MDConfig.from_interchange(interchange)