Merge pull request #665 from openforcefield/write-settles

Write `[ settles ]`

docs/releasehistory.md

@@ -31,6 +31,7 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 * #649 Removes the use of `pkg_resources`, which is deprecated.
 * #660 Moves the contents of `openff.interchange.components.foyer` to `openff.interchange.foyer` while maintaining existing import paths.
 * #663 Improves the performance of `Interchange.to_prmtop`.
+* #665 Properly write `[ settes ]` directive in GROMACS files.
 
 ## 0.3.0 - 2023-04-10
 

openff/interchange/__init__.py

@@ -2,7 +2,7 @@
 import importlib
 from types import ModuleType
 
-from openff.interchange._version import get_versions  # type: ignore
+from openff.interchange._version import get_versions
 from openff.interchange.components.interchange import Interchange
 
 # Handle versioneer

openff/interchange/_tests/__init__.py

@@ -285,7 +285,7 @@ def acetaldehyde(self):
 
     @pytest.fixture()
     def water(self):
-        return Molecule.from_mapped_smiles("[H:1][O:2][H:3]")
+        return Molecule.from_mapped_smiles("[H:2][O:1][H:3]")
 
 
 HAS_GROMACS = _find_gromacs_executable() is not None

openff/interchange/_tests/unit_tests/smirnoff/test_create.py

@@ -55,7 +55,7 @@ def test_sage_tip3p_charges(self, water, sage):
         out = Interchange.from_smirnoff(force_field=sage, topology=[water])
         found_charges = [v.m for v in out["Electrostatics"].charges.values()]
 
-        assert numpy.allclose(found_charges, [0.417, -0.834, 0.417])
+        assert numpy.allclose(found_charges, [-0.834, 0.417, 0.417])
 
     def test_infer_positions(self, sage):
         from openff.toolkit.tests.create_molecules import create_ethanol

openff/interchange/_tests/unit_tests/smirnoff/test_gromacs.py

@@ -1,9 +1,17 @@
 """Test SMIRNOFF-GROMACS conversion."""
+import pytest
 from openff.toolkit import Molecule
+from openff.units import unit
 
 from openff.interchange import Interchange
 from openff.interchange._tests import _BaseTest
-from openff.interchange.smirnoff._gromacs import _convert
+from openff.interchange.interop.gromacs.models.models import GROMACSMolecule
+from openff.interchange.smirnoff._gromacs import (
+    _convert,
+    _convert_angles,
+    _convert_bonds,
+    _convert_settles,
+)
 
 
 class TestConvert(_BaseTest):
@@ -25,3 +33,86 @@ def test_residue_names(self, sage):
         for molecule_type in system.molecule_types.values():
             for atom in molecule_type.atoms:
                 assert atom.residue_name == "LIG"
+
+
+class TestSettles(_BaseTest):
+    @pytest.fixture()
+    def tip3p_interchange(self, tip3p, water):
+        return tip3p.create_interchange(water.to_topology())
+
+    def test_catch_other_water_ordering(self, tip3p):
+        molecule = Molecule.from_mapped_smiles("[H:1][O:2][H:3]")
+        interchange = tip3p.create_interchange(molecule.to_topology())
+
+        with pytest.raises(Exception, match="OHH"):
+            _convert_settles(
+                GROMACSMolecule(name="foo"),
+                interchange.topology.molecule(0),
+                interchange,
+            )
+
+    def test_convert_settles(self, tip3p_interchange):
+        molecule = GROMACSMolecule(name="foo")
+
+        _convert_settles(
+            molecule,
+            tip3p_interchange.topology.molecule(0),
+            tip3p_interchange,
+        )
+
+        assert len(molecule.settles) == 1
+
+        settle = molecule.settles[0]
+
+        assert settle.first_atom == 1
+        assert settle.hydrogen_hydrogen_distance.m_as(unit.angstrom) == pytest.approx(
+            1.5139006545247014,
+        )
+        assert settle.oxygen_hydrogen_distance.m_as(unit.angstrom) == pytest.approx(
+            0.9572,
+        )
+
+        assert molecule.exclusions[0].first_atom == 1
+        assert molecule.exclusions[0].other_atoms == [2, 3]
+        assert molecule.exclusions[1].first_atom == 2
+        assert molecule.exclusions[1].other_atoms == [3]
+
+    def test_convert_no_settles_unconstrained_water(self, tip3p_interchange):
+        tip3p_interchange.collections["Constraints"].key_map = dict()
+
+        molecule = GROMACSMolecule(name="foo")
+
+        _convert_settles(
+            molecule,
+            tip3p_interchange.topology.molecule(0),
+            tip3p_interchange,
+        )
+
+        assert len(molecule.settles) == 0
+
+    def test_convert_no_settles_no_constraints(self, tip3p_interchange):
+        tip3p_interchange.collections.pop("Constraints")
+
+        molecule = GROMACSMolecule(name="foo")
+
+        _convert_settles(
+            molecule,
+            tip3p_interchange.topology.molecule(0),
+            tip3p_interchange,
+        )
+
+        assert len(molecule.settles) == 0
+
+    def test_no_bonds_or_angles_if_settle(self, tip3p_interchange):
+        molecule = GROMACSMolecule(name="foo")
+
+        for function in [_convert_settles, _convert_bonds, _convert_angles]:
+            function(
+                molecule,
+                tip3p_interchange.topology.molecule(0),
+                tip3p_interchange,
+            )
+
+        assert len(molecule.settles) == 1
+        assert len(molecule.angles) == 0
+        assert len(molecule.bonds) == 0

openff/interchange/_version.py

@@ -1,5 +1,3 @@
-# type: ignore # noqa
-
 # This file helps to compute a version number in source trees obtained from
 # git-archive tarball (such as those provided by githubs download-from-tag
 # feature). Distribution tarballs (built by setup.py sdist) and build
@@ -16,6 +14,7 @@
 import re
 import subprocess
 import sys
+from typing import Any
 
 
 def get_keywords():
@@ -53,8 +52,8 @@ class NotThisMethod(Exception):
     """Exception raised if a method is not valid for the current scenario."""
 
 
-LONG_VERSION_PY = {}
-HANDLERS = {}
+LONG_VERSION_PY: dict = dict()
+HANDLERS: dict = {}
 
 
 def register_vcs_handler(vcs, method):  # decorator
@@ -506,7 +505,7 @@ def render(pieces, style):
     }
 
 
-def get_versions():
+def get_versions() -> dict[str, Any]:
     """Get version information or return default if unable to do so."""
     # I am in _version.py, which lives at ROOT/VERSIONFILE_SOURCE. If we have
     # __file__, we can work backwards from there to the root. Some

openff/interchange/smirnoff/_gromacs.py

@@ -1,3 +1,4 @@
+import itertools
 from typing import Optional
 
 from openff.toolkit.topology.molecule import Atom, Molecule
@@ -16,15 +17,17 @@
     GROMACSAngle,
     GROMACSAtom,
     GROMACSBond,
+    GROMACSExclusion,
     GROMACSMolecule,
     GROMACSPair,
+    GROMACSSettles,
     GROMACSSystem,
     LennardJonesAtomType,
     PeriodicImproperDihedral,
     PeriodicProperDihedral,
     RyckaertBellemansDihedral,
 )
-from openff.interchange.models import TopologyKey
+from openff.interchange.models import BondKey, TopologyKey
 
 
 def _convert(interchange: Interchange) -> GROMACSSystem:
@@ -171,12 +174,12 @@ def _convert(interchange: Interchange) -> GROMACSSystem:
             else:
                 raise NotImplementedError()
 
+        _convert_settles(molecule, unique_molecule, interchange)
         _convert_bonds(molecule, unique_molecule, interchange)
         _convert_angles(molecule, unique_molecule, interchange)
         # pairs
         _convert_dihedrals(molecule, unique_molecule, interchange)
-        # settles?
-        # constraints?
+        # other constraints?
 
         system.molecule_types[unique_molecule.name] = molecule
 
@@ -199,6 +202,9 @@ def _convert_bonds(
     unique_molecule: "Molecule",
     interchange: Interchange,
 ):
+    if len(molecule.settles) > 0:
+        return
+
     collection = interchange["Bonds"]
 
     for bond in unique_molecule.bonds:
@@ -246,6 +252,9 @@ def _convert_angles(
     unique_molecule: "Molecule",
     interchange: Interchange,
 ):
+    if len(molecule.settles) > 0:
+        return
+
     collection = interchange["Angles"]
 
     for angle in unique_molecule.angles:
@@ -407,3 +416,75 @@ def _convert_dihedrals(
                             multiplicity=int(params["periodicity"]),
                         ),
                     )
+
+
+def _convert_settles(
+    molecule: GROMACSMolecule,
+    unique_molecule: "Molecule",
+    interchange: Interchange,
+):
+    if "Constraints" not in interchange.collections:
+        return
+
+    if not unique_molecule.is_isomorphic_with(Molecule.from_smiles("O")):
+        return
+
+    if unique_molecule.atom(0).atomic_number != 8:
+        raise Exception(
+            "Writing `[ settles ]` assumes water is ordered as OHH. Please raise an issue "
+            "if you would benefit from this assumption changing.",
+        )
+
+    topology_atom_indices = [
+        interchange.topology.atom_index(atom) for atom in unique_molecule.atoms
+    ]
+
+    constraint_lengths = set()
+
+    for atom_pair in itertools.combinations(topology_atom_indices, 2):
+        key = BondKey(atom_indices=atom_pair)
+
+        if key not in interchange["Constraints"].key_map:
+            return
+
+        try:
+            constraint_lengths.add(
+                interchange["Bonds"]
+                .potentials[interchange["Bonds"].key_map[key]]
+                .parameters["length"],
+            )
+        # KeyError (subclass of LookupErrorR) when this BondKey is not found in the bond collection
+        # LookupError for the Interchange not having a bond collection
+        # in either case, look to the constraint collection for this distance
+        except LookupError:
+            constraint_lengths.add(
+                interchange["Constraints"]
+                .constraints[interchange["Constraints"].key_map[key]]
+                .parameters["distance"],
+            )
+
+    if len(constraint_lengths) != 2:
+        raise RuntimeError(
+            "Found three unique constraint lengths in constrained water.",
+        )
+
+    molecule.settles.append(
+        GROMACSSettles(
+            first_atom=1,  # TODO: documentation unclear on if this is first or oxygen
+            oxygen_hydrogen_distance=min(constraint_lengths),
+            hydrogen_hydrogen_distance=max(constraint_lengths),
+        ),
+    )
+
+    molecule.exclusions.append(
+        GROMACSExclusion(
+            first_atom=1,
+            other_atoms=[2, 3],
+        ),
+    )
+    molecule.exclusions.append(
+        GROMACSExclusion(
+            first_atom=2,
+            other_atoms=[3],
+        ),
+    )

openff/interchange/smirnoff/_valence.py

@@ -287,7 +287,7 @@ class SMIRNOFFConstraintCollection(SMIRNOFFCollection):
     expression: Literal[""] = ""
     constraints: dict[
         PotentialKey,
-        bool,
+        Potential,
     ] = dict()  # should this be named potentials for consistency?
 
     @classmethod
@@ -390,7 +390,7 @@ def store_constraints(
                     "distance": distance,
                 },
             )
-            self.constraints[potential_key] = potential  # type: ignore[assignment]
+            self.constraints[potential_key] = potential
 
 
 class SMIRNOFFAngleCollection(SMIRNOFFCollection, AngleCollection):