Merge pull request #883 from openforcefield/interop

Interoperability trunk

.github/workflows/beta.yaml

@@ -71,5 +71,6 @@ jobs:
     - name: Codecov
       uses: codecov/codecov-action@v4
       with:
+        token: ${{ secrets.CODECOV_TOKEN }}
         file: ./coverage.xml
         fail_ci_if_error: false

.github/workflows/ci.yaml

@@ -144,5 +144,6 @@ jobs:
     - name: Codecov
       uses: codecov/codecov-action@v4
       with:
+        token: ${{ secrets.CODECOV_TOKEN }}
         file: ./coverage.xml
         fail_ci_if_error: false

devtools/conda-envs/beta_env.yaml

@@ -12,7 +12,7 @@ dependencies:
   - pydantic =1
   - openmm >=7.6
   # OpenFF stack
-  - openff-toolkit >=0.14.5
+  - openff-toolkit >=0.15.2
   - openff-models
   - openff-nagl
   - openff-nagl-models

devtools/conda-envs/dev_env.yaml

@@ -11,7 +11,7 @@ dependencies:
   - pydantic
   - openmm
   # OpenFF stack
-  - openff-toolkit ==0.14.5
+  - openff-toolkit =0.15.2
   - openff-interchange-base
   - openff-models
   - smirnoff-plugins =2023

devtools/conda-envs/docs_env.yaml

@@ -9,7 +9,7 @@ dependencies:
   - pip
   - numpy =1
   - pydantic =1
-  - openff-toolkit-base =0.14.5
+  - openff-toolkit-base =0.15.2
   - openff-models
   - openmm >=7.6
   - mbuild

devtools/conda-envs/examples_env.yaml

@@ -10,7 +10,8 @@ dependencies:
   - pydantic =1
   - openmm >=7.6
   # OpenFF stack
-  - openff-toolkit =0.14.5
+  - openff-toolkit =0.15.2
+  - openff-models
   - openff-nagl ==0.3.1
   - openff-nagl-models ==0.1
   # Optional features

devtools/conda-envs/test_env.yaml

@@ -8,7 +8,7 @@ dependencies:
   - numpy
   - pydantic
   # OpenFF stack
-  - openff-toolkit-base =0.14.5
+  - openff-toolkit-base =0.15.2
   - openff-units
   - openff-models
   # Needs to be explicitly listed to not be dropped when AmberTools is removed

docs/releasehistory.md

@@ -13,15 +13,26 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 
 ## Current development
 
-### New Features
+### 0.3.19 - 2023-02-05
 
 * #867 Tags `PotentialKey.virtual_site_type` with the associated type provided by SMIRNOFF parameters.
+* #857 Tags `PotentialKey.associated_handler` when importing data from OpenMM.
 * #848 Raises a more useful error when `Interchange.minimize` is called while positions are not present.
 * #852 Support LJPME in OpenMM export.
+* #871 Re-introduces Foyer compatibility with version 0.12.1.
+* #883 Improve topology interoperability after importing data from OpenMM, using OpenFF Toolkit 0.15.2.
+* #883 Falls back to `Topology.visualize` in most cases.
 
 ### Bugfixes
 
 * #848 Fixes a bug in which `Interchange.minimize` erroneously appended virtual site positions to the `positions` attribute.
+* #883 Using `openff-models` 0.1.2, fixes parsing box information from OpenMM data.
+* #883 Skips writing unnecessary PDB file during visualization.
+* #883 Preserves atom metadata when roundtripping topologies with OpenMM.
+
+### Documentation improvements
+
+* #864 Updates installation instructions.
 
 ## 0.3.18 - 2023-11-16
 

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

@@ -305,13 +305,14 @@ def test_to_openmm_simulation(self, sage):
 
     @skip_if_missing("nglview")
     @skip_if_missing("openmm")
-    def test_visualize(self, sage):
+    @pytest.mark.parametrize("include_virtual_sites", [True, False])
+    def test_visualize(self, include_virtual_sites, tip4p, sage):
         import nglview
 
-        molecule = Molecule.from_smiles("CCO")
+        molecule = Molecule.from_smiles("O")
 
         out = Interchange.from_smirnoff(
-            force_field=sage,
+            force_field=tip4p if include_virtual_sites else sage,
             topology=molecule.to_topology(),
         )
 
@@ -324,7 +325,10 @@ def test_visualize(self, sage):
         molecule.generate_conformers(n_conformers=1)
         out.positions = molecule.conformers[0]
 
-        assert isinstance(out.visualize(), nglview.NGLWidget)
+        assert isinstance(
+            out.visualize(include_virtual_sites=include_virtual_sites),
+            nglview.NGLWidget,
+        )
 
 
 @skip_if_missing("openmm")

openff/interchange/_tests/unit_tests/interop/openmm/_import/test_import.py

@@ -1,8 +1,10 @@
+import copy
 import random
 
+import numpy
 import pytest
-from openff.units import unit
-from openff.utilities import has_package, skip_if_missing
+from openff.toolkit import Molecule, Topology, unit
+from openff.utilities import get_data_file_path, has_package, skip_if_missing
 
 from openff.interchange import Interchange
 from openff.interchange._tests import _BaseTest
@@ -49,6 +51,64 @@ def test_simple_roundtrip(self, monkeypatch, sage_unconstrained, ethanol):
             },
         )
 
+        assert isinstance(converted.box.m, numpy.ndarray)
+
+        # OpenMM seems to avoid using the built-in type
+        assert converted.box.m.dtype in (float, numpy.float32, numpy.float64)
+
+    def test_openmm_roundtrip_metadata(self, monkeypatch):
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+        # Make an example OpenMM Topology with metadata.
+        # Here we use OFFTK to make the OpenMM Topology, but this could just as easily come from another source
+        ethanol = Molecule.from_smiles("CCO")
+        benzene = Molecule.from_smiles("c1ccccc1")
+        for atom in ethanol.atoms:
+            atom.metadata["chain_id"] = "1"
+            atom.metadata["residue_number"] = "1"
+            atom.metadata["insertion_code"] = ""
+            atom.metadata["residue_name"] = "ETH"
+        for atom in benzene.atoms:
+            atom.metadata["chain_id"] = "1"
+            atom.metadata["residue_number"] = "2"
+            atom.metadata["insertion_code"] = "A"
+            atom.metadata["residue_name"] = "BNZ"
+        top = Topology.from_molecules([ethanol, benzene])
+
+        # Roundtrip the topology with metadata through openmm
+        interchange = Interchange.from_openmm(topology=top.to_openmm())
+
+        # Ensure that the metadata is the same
+        for atom in interchange.topology.molecule(0).atoms:
+            assert atom.metadata["chain_id"] == "1"
+            assert atom.metadata["residue_number"] == "1"
+            assert atom.metadata["insertion_code"] == ""
+            assert atom.metadata["residue_name"] == "ETH"
+        for atom in interchange.topology.molecule(1).atoms:
+            assert atom.metadata["chain_id"] == "1"
+            assert atom.metadata["residue_number"] == "2"
+            assert atom.metadata["insertion_code"] == "A"
+            assert atom.metadata["residue_name"] == "BNZ"
+
+    def test_openmm_native_roundtrip_metadata(self, monkeypatch):
+        """
+        Test that metadata is the same whether we load a PDB through OpenMM+Interchange vs. Topology.from_pdb.
+        """
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+        pdb = openmm.app.PDBFile(
+            get_data_file_path("ALA_GLY/ALA_GLY.pdb", "openff.interchange._tests.data"),
+        )
+        interchange = Interchange.from_openmm(topology=pdb.topology)
+        off_top = Topology.from_pdb(
+            get_data_file_path("ALA_GLY/ALA_GLY.pdb", "openff.interchange._tests.data"),
+        )
+        for roundtrip_atom, off_atom in zip(interchange.topology.atoms, off_top.atoms):
+            # off_atom's metadata also includes a little info about how the chemistry was
+            # assigned, so we remove this from the comparison
+            off_atom_metadata = copy.deepcopy(off_atom.metadata)
+            del off_atom_metadata["match_info"]
+            assert roundtrip_atom.metadata == off_atom_metadata
+
 
 @skip_if_missing("openmm")
 class TestConvertNonbondedForce:

openff/interchange/components/_viz.py

@@ -0,0 +1,27 @@
+import uuid
+from io import StringIO
+
+from openff.toolkit.utils._viz import TopologyNGLViewStructure
+
+from openff.interchange import Interchange
+
+
+class InterchangeNGLViewStructure(TopologyNGLViewStructure):
+    """Subclass of the toolkit's NGLView interface."""
+
+    def __init__(
+        self,
+        interchange: Interchange,
+        ext: str = "PDB",
+    ):
+        self.interchange = interchange
+        self.ext = ext.lower()
+        self.params: dict = dict()
+        self.id = str(uuid.uuid4())
+
+    def get_structure_string(self) -> str:
+        """Get structure as a string."""
+        with StringIO() as f:
+            self.interchange.to_pdb(f, include_virtual_sites=True)
+            structure_string = f.getvalue()
+        return structure_string

openff/interchange/components/interchange.py

@@ -282,7 +282,7 @@ def visualize(
         self,
         backend: str = "nglview",
         include_virtual_sites: bool = False,
-    ):
+    ) -> "nglview.NGLWidget":
         """
         Visualize this Interchange.
 
@@ -301,9 +301,42 @@ def visualize(
             The NGLWidget containing the visualization.
 
         """
+        from openff.toolkit.utils.exceptions import (
+            IncompatibleUnitError,
+            MissingConformersError,
+        )
+
         if backend == "nglview":
-            return self._visualize_nglview(include_virtual_sites=include_virtual_sites)
+            if include_virtual_sites:
+
+                return self._visualize_nglview(include_virtual_sites=True)
+
+            else:
+
+                # Interchange.topology might have its own positions;
+                # just use Interchange.positions
+                original_positions = self.topology.get_positions()
+
+                try:
+                    self.topology.set_positions(self.positions)
+                    widget = self.topology.visualize()
+                except (MissingConformersError, IncompatibleUnitError) as error:
+                    raise MissingPositionsError(
+                        "Cannot visualize system without positions.",
+                    ) from error
+
+                # but don't modify them long-term
+                # work around https://github.com/openforcefield/openff-toolkit/issues/1820
+                if original_positions is not None:
+                    self.topology.set_positions(original_positions)
+                else:
+                    for molecule in self.topology.molecules:
+                        molecule._conformers = None
+
+                return widget
+
         else:
+
             raise UnsupportedExportError
 
     @requires_package("nglview")
@@ -319,16 +352,35 @@ def _visualize_nglview(
         """
         import nglview
 
+        from openff.interchange.components._viz import InterchangeNGLViewStructure
+
         try:
-            self.to_pdb(
-                "_tmp_pdb_file.pdb",
-                include_virtual_sites=include_virtual_sites,
+            widget = nglview.NGLWidget(
+                InterchangeNGLViewStructure(
+                    interchange=self,
+                    ext="pdb",
+                ),
+                representations=[
+                    dict(type="unitcell", params=dict()),
+                ],
             )
+
         except MissingPositionsError as error:
             raise MissingPositionsError(
                 "Cannot visualize system without positions.",
             ) from error
-        return nglview.show_structure_file("_tmp_pdb_file.pdb")
+
+        widget.add_representation("line", sele="water")
+        widget.add_representation("spacefill", sele="ion")
+        widget.add_representation("cartoon", sele="protein")
+        widget.add_representation(
+            "licorice",
+            sele="not water and not ion and not protein",
+            radius=0.25,
+            multipleBond=False,
+        )
+
+        return widget
 
     def minimize(
         self,

openff/interchange/components/toolkit.py

@@ -2,13 +2,11 @@
 
 from typing import TYPE_CHECKING, Union
 
-import networkx as nx
-import numpy as np
-from openff.toolkit import Molecule, Topology
+import networkx
+import numpy
+from openff.toolkit import ForceField, Molecule, Quantity, Topology
 from openff.toolkit.topology._mm_molecule import _SimpleMolecule
-from openff.toolkit.typing.engines.smirnoff import ForceField
 from openff.toolkit.utils.collections import ValidatedList
-from openff.units import Quantity
 from openff.utilities.utilities import has_package
 
 if has_package("openmm") or TYPE_CHECKING:
@@ -66,7 +64,7 @@ def _validated_list_to_array(validated_list: "ValidatedList") -> "Quantity":
     from openff.units import unit
 
     unit_ = validated_list[0].units
-    return unit.Quantity(np.asarray([val.m for val in validated_list]), unit_)
+    return unit.Quantity(numpy.asarray([val.m for val in validated_list]), unit_)
 
 
 def _combine_topologies(topology1: Topology, topology2: Topology) -> Topology:
@@ -110,14 +108,24 @@ def _check_electrostatics_handlers(force_field: "ForceField") -> bool:
 
 def _simple_topology_from_openmm(openmm_topology: "openmm.app.Topology") -> Topology:
     """Convert an OpenMM Topology into an OpenFF Topology consisting **only** of so-called `_SimpleMolecule`s."""
-    # TODO: Residue metadata
     # TODO: Splice in fully-defined OpenFF `Molecule`s?
-    graph = nx.Graph()
 
+    graph = networkx.Graph()
+
+    # TODO: This is nearly identical to Topology._openmm_topology_to_networkx.
+    #  Should this method be replaced with a direct call to that?
     for atom in openmm_topology.atoms():
         graph.add_node(
             atom.index,
             atomic_number=atom.element.atomic_number,
+            name=atom.name,
+            residue_name=atom.residue.name,
+            # Note that residue number is mapped to residue.id here. The use of id vs. number varies in other packages
+            # and the convention for the OpenFF-OpenMM interconversion is recorded at
+            # https://docs.openforcefield.org/projects/toolkit/en/0.15.1/users/molecule_conversion.html
+            residue_number=atom.residue.id,
+            insertion_code=atom.residue.insertionCode,
+            chain_id=atom.residue.chain.id,
         )
 
     for bond in openmm_topology.bonds():
@@ -129,10 +137,11 @@ def _simple_topology_from_openmm(openmm_topology: "openmm.app.Topology") -> Topo
     return _simple_topology_from_graph(graph)
 
 
-def _simple_topology_from_graph(graph: nx.Graph) -> Topology:
+def _simple_topology_from_graph(graph: networkx.Graph) -> Topology:
+    """Convert a networkx Graph into an OpenFF Topology consisting only of `_SimpleMolecule`s."""
     topology = Topology()
 
-    for component in nx.connected_components(graph):
+    for component in networkx.connected_components(graph):
         subgraph = graph.subgraph(component)
         topology.add_molecule(_SimpleMolecule._from_subgraph(subgraph))