Merge remote-tracking branch 'upstream/main' into release-0.4.3

Author: mattwthompson

docs/releasehistory.md

@@ -17,8 +17,9 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 
 * #1174 Support Python 3.13.
 * #1216 Type labels can (optionally) be included in LAMMPS files.
+* #1219 Adds `SMIRNOFFElectrostaticsCollection.get_charge_array`, which currently only works with systems lacking virtual sites.
+* #1220 Adds `Interchange.set_positions_from_gro`.
 * #1250 Allow `Interchange.combine` to proceed when cutoffs differ by up to 1e-6 nanometers.
-* #1219 Adds `SMIRNOFFElectrostaticsCollection.get_charge_array`.
 
 ### Behavior changes
 

openff/interchange/_tests/data/mdp/em.mdp

@@ -0,0 +1,20 @@
+integrator              = steep
+nsteps                  = 10000
+nstxout                 = 0
+nstlog                  = 0
+nstenergy               = 0
+
+cutoff-scheme           = Verlet
+
+coulombtype             = PME
+rcoulomb                = 0.9
+verlet-buffer-tolerance = 0.005
+rlist                   = 0.9
+
+vdwtype                 = Cut-off
+vdw-modifier            = force-switch
+rvdw                    = 0.9
+rvdw-switch             = 0.8
+DispCorr                = EnerPres
+
+constraints             = h-bonds

openff/interchange/_tests/unit_tests/components/test_interchange.py

@@ -7,6 +7,7 @@
     ElectrostaticsHandler,
     ParameterHandler,
 )
+from openff.utilities import temporary_cd
 from openff.utilities.testing import skip_if_missing
 from pydantic import ValidationError
 
@@ -17,9 +18,10 @@
     needs_gmx,
     needs_lmp,
 )
-from openff.interchange.drivers import get_openmm_energies
+from openff.interchange.drivers import get_gromacs_energies, get_openmm_energies
 from openff.interchange.exceptions import (
     ExperimentalFeatureException,
+    InvalidPositionsError,
     InvalidTopologyError,
     MissingParameterHandlerError,
     MissingParametersError,
@@ -124,9 +126,7 @@ def test_atom_ordering(self):
 
         from openff.interchange import Interchange
         from openff.interchange.drivers import (
-            get_gromacs_energies,
             get_lammps_energies,
-            get_openmm_energies,
         )
 
         oplsaa = foyer.forcefields.load_OPLSAA()
@@ -442,6 +442,75 @@ def test_from_gromacs_called(self, simple_interchange, monkeypatch):
             gro_file="tmp_.gro",
         )
 
+    @needs_gmx
+    def test_set_positions_from_gro_same_coordinates(self, simple_interchange):
+        """Write and read back the same coordinates, make sure energies match."""
+        with temporary_cd():
+            simple_interchange.positions = simple_interchange.positions.round(3)
+            simple_interchange.to_gromacs(prefix="test12")
+
+            original_energies = get_gromacs_energies(simple_interchange)
+
+            simple_interchange.set_positions_from_gro(gro_file="test12.gro")
+
+            new_energies = get_gromacs_energies(simple_interchange)
+
+            # if energies match, positions must also match
+            assert original_energies.energies.keys() == new_energies.energies.keys()
+
+            for key in original_energies.energies:
+                assert numpy.allclose(
+                    original_energies[key].m,
+                    new_energies[key].m,
+                )
+
+    @needs_gmx
+    def test_set_positions_from_gro_minimize(self, simple_interchange):
+        """Write, minimize, and read new coordinates, make sure energy descreased."""
+        with temporary_cd():
+            simple_interchange.positions = simple_interchange.positions.round(3)
+
+            # with 0.001 nm precision, not guaranteed a difference between conformer-generation
+            # coordinates and MM-minimized coordinates, so just perturb a little bit
+            simple_interchange.positions[-1] += Quantity([0.03, 0.03, 0.03], "nanometer")
+            simple_interchange.to_gromacs(prefix="test13")
+
+            original_energies = get_gromacs_energies(simple_interchange)
+
+            mdp = get_test_file_path("mdp/em.mdp")
+
+            for cmd in [
+                f"gmx grompp -f {mdp} -c test13.gro -p test13.top -o test13.tpr",
+                "gmx mdrun -s test13.tpr -c final_coords.gro",
+                "echo 0 | gmx trjconv -f final_coords.gro -s test13.tpr -o final_coords_wrapped.gro -pbc whole",
+            ]:
+                subprocess.check_output(cmd, shell=True)
+
+            simple_interchange.set_positions_from_gro(gro_file="final_coords_wrapped.gro")
+
+            new_energies = get_gromacs_energies(simple_interchange)
+
+            assert new_energies.total_energy < original_energies.total_energy
+
+    @needs_gmx
+    def test_set_positions_from_gro_wrong_coordinates(self, sage):
+        with temporary_cd():
+            eth = sage.create_interchange(MoleculeWithConformer.from_smiles("CCO").to_topology())
+            eth.box = Quantity([4, 4, 4], unit.nanometer)
+            eth.to_gromacs(prefix="ethanol")
+            benzene = sage.create_interchange(MoleculeWithConformer.from_smiles("c1ccccc1").to_topology())
+            benzene.box = Quantity([4, 4, 4], unit.nanometer)
+            benzene.to_gromacs(prefix="benzene")
+
+            with pytest.raises(InvalidPositionsError, match="12.*9"):
+                eth.set_positions_from_gro(gro_file="benzene.gro")
+
+            with pytest.raises(
+                InvalidPositionsError,
+                match="9.*12",
+            ):
+                benzene.set_positions_from_gro(gro_file="ethanol.gro")
+
     @skip_if_missing("openmm")
     def test_minimize(self, simple_interchange):
         original_energy = get_openmm_energies(simple_interchange).total_energy

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

@@ -287,6 +287,20 @@ def test_get_charge_array(self, sage):
         )
 
 
+def test_get_charge_array_fails_if_virtual_sites_present(water, tip4p):
+    with pytest.raises(
+        NotImplementedError,
+        match="Not yet implemented with virtual sites",
+    ):
+        tip4p.create_interchange(water.to_topology())["Electrostatics"].get_charge_array(include_virtual_sites=True)
+
+    with pytest.raises(
+        NotImplementedError,
+        match="Not yet implemented when virtual sites are present",
+    ):
+        tip4p.create_interchange(water.to_topology())["Electrostatics"].get_charge_array(include_virtual_sites=False)
+
+
 def test_nonintegral_molecule_charge_error(sage, water):
     funky_charges = Quantity([0, 0, -5.5], "elementary_charge")
 

openff/interchange/components/interchange.py

@@ -27,6 +27,7 @@
 from openff.interchange.components.mdconfig import MDConfig
 from openff.interchange.components.potentials import Collection, _AnnotatedCollections
 from openff.interchange.exceptions import (
+    InvalidPositionsError,
     MissingParameterHandlerError,
     MissingPositionsError,
     UnsupportedExportError,
@@ -903,6 +904,33 @@ def from_gromacs(
             ),
         )
 
+    def set_positions_from_gro(
+        self,
+        gro_file: Path | str,
+    ):
+        """
+        Set the positions of this `Interchange` from a GROMACS coordinate `.gro` file.
+
+        Only the coordinates from the `.gro` file are used. No effort is made to ensure the topologies are compatible
+        with each other. This includes, for example, a lack of guarantee that the atom ordering in the `.gro` file
+        matches the atom ordering in the `Interchange` object.
+
+        `InvalidPositionsError` is raised if the number of rows in the coordinate array does not match the number of
+        atoms in the topology of this `Interchange`.
+        """
+        from openff.interchange.interop.gromacs._import._import import _read_coordinates
+
+        # should already be in nm, might not be necessary
+        coordinates = _read_coordinates(gro_file).to(unit.nanometer)
+
+        if coordinates.shape != (self.topology.n_atoms, 3):
+            raise InvalidPositionsError(
+                f"Coordinates in {gro_file} do not match the number of atoms in the topology. ",
+                f"Parsed coordinates have shape {coordinates.shape} but topology has {self.topology.n_atoms} atoms.",
+            )
+
+        self.positions = coordinates
+
     @classmethod
     def from_openmm(
         cls,

openff/interchange/exceptions.py

@@ -38,6 +38,12 @@ class InvalidBoxError(InterchangeException, ValueError):
     """
 
 
+class InvalidPositionsError(InterchangeException, ValueError):
+    """
+    Generic exception for errors with positions.
+    """
+
+
 class InvalidTopologyError(InterchangeException, ValueError):
     """
     Generic exception for errors reading chemical topology data.

openff/interchange/smirnoff/_nonbonded.py

@@ -330,10 +330,19 @@ def charges(
         return self._charges
 
     def get_charge_array(self, include_virtual_sites: bool = False) -> Quantity:
-        """Return a one-dimensional array-like of atomic charges, ordered topologically."""
+        """
+        Return a one-dimensional array-like of atomic charges, ordered topologically.
+
+        If virtual sites are present in the system, `NotImplementedError` is raised.
+        """
         if include_virtual_sites:
             raise NotImplementedError("Not yet implemented with virtual sites")
 
+        if VirtualSiteKey in {type(key) for key in self.key_map}:
+            raise NotImplementedError(
+                "Not yet implemented when virtual sites are present, even with `include_virtual_sites=False`.",
+            )
+
         return Quantity.from_list([q for _, q in sorted(self.charges.items(), key=lambda x: x[0].atom_indices)])
 
     def _get_charges(