Merge pull request #1134 from openforcefield/from-openmm-water-box

Create blank collections when importing only water

docs/releasehistory.md

@@ -22,6 +22,7 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 
 * #1070 The `charge_from_molecules` argument must include only molecules that contain partial charges and are non-isomorphic with each other.
 * #1070 The `charge_from_molecules` argument as used by the OpenFF Toolkit is handled internally as `molecules_with_preset_charges`.
+* #1134 Fixes a bug in which importing OpenMM systems containing only rigid water would crash.
 
 ### Performance improvements
 

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

@@ -2,8 +2,9 @@
 import random
 
 import numpy
+import openmm
 import pytest
-from openff.toolkit import Molecule, Quantity, Topology, unit
+from openff.toolkit import ForceField, Molecule, Quantity, Topology, unit
 from openff.utilities import get_data_file_path, has_package, skip_if_missing
 
 from openff.interchange import Interchange
@@ -249,6 +250,33 @@ def test_fill_in_rigid_water_parameters(self, water_dimer, monkeypatch):
             },
         )
 
+    def test_only_constrained_water(self, water_dimer):
+        water_dimer.box_vectors = Quantity([4, 4, 4], "nanometer")
+
+        interchange = ForceField("openff-2.2.1.offxml").create_interchange(water_dimer)
+
+        simulation = interchange.to_openmm_simulation(integrator=openmm.LangevinIntegrator(300, 1, 0.001))
+        system = simulation.system
+
+        for index in range(system.getNumForces()):
+            if isinstance(system.getForce(index), openmm.HarmonicBondForce):
+                break
+
+        system.removeForce(index)
+
+        for index in range(system.getNumForces()):
+            if isinstance(system.getForce(index), openmm.HarmonicAngleForce):
+                break
+
+        system.removeForce(index)
+
+        interchange2 = Interchange.from_openmm(
+            system=simulation.system,
+            topology=simulation.topology,
+        )
+
+        assert interchange2.topology.n_bonds == interchange.topology.n_bonds
+
 
 @skip_if_missing("openmm")
 class TestProcessTopology:

openff/interchange/interop/openmm/_import/_import.py

@@ -122,7 +122,12 @@ def from_openmm(
     # TODO: Does this run through the Interchange.box validator?
     interchange.box = _box_vectors
 
-    if interchange.topology.n_bonds > len(interchange.collections["Bonds"].key_map):
+    try:
+        num_physics_bonds = len(interchange["Bonds"].key_map)
+    except LookupError:
+        num_physics_bonds = 0
+
+    if interchange.topology.n_bonds > num_physics_bonds:
         # There are probably missing (physics) bonds from rigid waters. The topological
         # bonds are probably processed correctly.
         _fill_in_rigid_water_bonds(interchange)
@@ -348,6 +353,12 @@ def _fill_in_rigid_water_bonds(interchange: "Interchange"):
     from openff.interchange.components.potentials import Potential
     from openff.interchange.models import AngleKey, BondKey, PotentialKey
 
+    if "Bonds" not in interchange.collections:
+        interchange.collections.update({"Bonds": BondCollection()})
+
+    if "Angles" not in interchange.collections:
+        interchange.collections.update({"Angles": AngleCollection()})
+
     simple_water = _SimpleMolecule.from_molecule(Molecule.from_smiles("O"))
 
     rigid_water_bond_key = PotentialKey(id="rigid_water", associated_handler="Bonds")