Merge pull request #1080 from openforcefield/allow-hmr-with-virtual-s…

…ites

Allow HMR with virtual sites in OpenMM

docs/releasehistory.md

@@ -16,19 +16,20 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 ### New features
 
 * #1053 Logs, at the level of `logging.INFO`, how charges are assigned by SMIRNOFF force fields to each atom and virtual site.
+* #1080 HMR is supported with OpenMM when virtual sites are present.
 
 ### Bug fixes
 
-* The `charge_from_molecules` argument must include only molecules that contain partial charges and are non-isomorphic with each other.
-* The `charge_from_molecules` argument as used by the OpenFF Toolkit is handled internally as `molecules_with_preset_charges`.
+* #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`.
 
 ### Performance improvements
 
 * #1097 Migrates version handling to `versioningit`, which should result in shorter import times.
 
 ### Documentation improvements
 
-* Documents charge assignment hierarchy in the user guide.
+* #1070 Documents charge assignment hierarchy in the user guide.
 
 ## 0.4.0 - 2024-11-04
 

docs/using/status.md

@@ -2,11 +2,11 @@
 
 | Engine | Known bugs | HMR | Bond constraints | 4-site water models | 5-site water models | Virtual sites on ligands | Proteins | Nucleic acids | Lipids | Carbohydrates |
 |:--|:--|:--|:--|:--|:--|:--|:--|:--|:--|:--|
-| OpenMM | [Link](https://github.com/openforcefield/openff-interchange/issues?q=is%3Aissue+is%3Aopen+label%3Aopenmm+label%3Abug) | Mostly<sup>2</sup> supported | Supported | Supported | Supported | Supported | Supported |
+| OpenMM | [Link](https://github.com/openforcefield/openff-interchange/issues?q=is%3Aissue+is%3Aopen+label%3Aopenmm+label%3Abug) | Supported | Supported | Supported | Supported | Supported | Supported |
 | GROMACS | [Link](https://github.com/openforcefield/openff-interchange/issues?q=is%3Aissue+is%3Aopen+label%3Agromacs+label%3Abug) | Mostly<sup>2</sup> supported | Partially supported | Partially supported | Supported | Partially supported | Supported |
 | Amber | [Link](https://github.com/openforcefield/openff-interchange/issues?q=is%3Aissue+is%3Aopen+label%3Aamber+label%3Abug) | Unsupported | Partially supported | Unsupported<sup>1</sup> | Unsupported<sup>1</sup> | Unsupported<sup>1</sup> | Supported |
 | LAMMPS | [Link](https://github.com/openforcefield/openff-interchange/issues?q=is%3Aissue+is%3Aopen+label%3Alammps+label%3Abug) | Unsupported | Partially supported | Not supported | Not supported | Not supported | Not tested |
 | All | | | | | | | <td colspan=3>No available SMIRNOFF force fields |
 
 1. Unable to find upstream documentation
-2. See caveats in `Interchange.to_openmm` and `Interchange.to_gromacs` docstrings.
+2. See caveats in [`Interchange.to_gromacs` docstring](openff.interchange.Interchange.to_gromacs).

openff/interchange/_tests/unit_tests/interop/openmm/test_hmr.py

@@ -1,9 +1,10 @@
 import random
 
 import pytest
-from openff.toolkit import Molecule, unit
+from openff.toolkit import Molecule, Topology, unit
 
-from openff.interchange.exceptions import NegativeMassError, UnsupportedExportError
+from openff.interchange._tests import MoleculeWithConformer
+from openff.interchange.exceptions import NegativeMassError
 
 
 @pytest.mark.parametrize("reversed", [False, True])
@@ -80,6 +81,60 @@ def test_mass_is_positive(sage):
         sage.create_interchange(Molecule.from_smiles("C").to_topology()).to_openmm(hydrogen_mass=5.0)
 
 
-def test_virtual_sites_unsupported(tip4p, water):
-    with pytest.raises(UnsupportedExportError, match="with virtual sites"):
-        tip4p.create_interchange(water.to_topology()).to_openmm(hydrogen_mass=2.0)
+@pytest.mark.parametrize("hydrogen_mass", [1.1, 3.0])
+def test_hmr_not_applied_to_tip4p(
+    tip4p,
+    water,
+    hydrogen_mass,
+):
+    pytest.importorskip("openmm")
+
+    system = tip4p.create_interchange(
+        Topology.from_molecules(2 * [water]),
+    ).to_openmm_system(
+        hydrogen_mass=hydrogen_mass,
+    )
+
+    masses = [system.getParticleMass(particle_index)._value for particle_index in range(system.getNumParticles())]
+
+    assert masses[0] == masses[3] == 15.99943
+    assert masses[1] == masses[2] == masses[4] == masses[5] == 1.007947
+    assert masses[6] == masses[7] == 0.0
+
+
+@pytest.mark.parametrize("hydrogen_mass", [1.1, 3.0])
+def test_hmr_with_ligand_virtual_sites(sage_with_bond_charge, hydrogen_mass):
+    """Test that a single-molecule sigma hole example runs"""
+    pytest.importorskip("openmm")
+    import openmm
+    import openmm.unit
+
+    molecule = MoleculeWithConformer.from_mapped_smiles("[H:5][C:2]([H:6])([C:3]([H:7])([H:8])[Cl:4])[Cl:1]")
+
+    # should work fine with 4 fs, but be conservative and just use 3 fs
+    simulation = sage_with_bond_charge.create_interchange(molecule.to_topology()).to_openmm_simulation(
+        integrator=openmm.VerletIntegrator(3.0 * openmm.unit.femtosecond),
+        hydrogen_mass=hydrogen_mass,
+    )
+
+    system = simulation.system
+
+    assert system.getParticleMass(0)._value == 35.4532
+    assert system.getParticleMass(3)._value == 35.4532
+
+    assert system.getParticleMass(1)._value == 12.01078 - 2 * (hydrogen_mass - 1.007947)
+    assert system.getParticleMass(2)._value == 12.01078 - 2 * (hydrogen_mass - 1.007947)
+
+    assert system.getParticleMass(4)._value == hydrogen_mass
+    assert system.getParticleMass(5)._value == hydrogen_mass
+    assert system.getParticleMass(6)._value == hydrogen_mass
+    assert system.getParticleMass(7)._value == hydrogen_mass
+
+    assert system.getParticleMass(8)._value == 0.0
+    assert system.getParticleMass(9)._value == 0.0
+
+    # make sure it can minimize and run briefly without crashing
+    # this is a single ligand in vacuum with boring parameters
+    simulation.minimizeEnergy(maxIterations=100)
+
+    simulation.runForClockTime(10 * openmm.unit.second)

openff/interchange/components/interchange.py

@@ -335,7 +335,7 @@ def to_gromacs(
             "foo.top", "foo.gro", and "foo_pointenergy.mdp".
         decimal: int, default=3
             The number of decimal places to use when writing the GROMACS coordinate file.
-        hydrogen_mass : PostitiveFloat, default=1.007947
+        hydrogen_mass : PositiveFloat, default=1.007947
             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.
@@ -403,7 +403,7 @@ def to_top(
         ----------
         file_path
             The path to the GROMACS topology file to write.
-        hydrogen_mass : PostitiveFloat, default=1.007947
+        hydrogen_mass : PositiveFloat, default=1.007947
             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.
@@ -535,10 +535,10 @@ def to_openmm_system(
             on a bond or angle that is fully constrained.
         ewald_tolerance : float, default=1e-4
             The value passed to `NonbondedForce.setEwaldErrorTolerance`
-        hydrogen_mass : PostitiveFloat, default=1.007947
+        hydrogen_mass : PositiveFloat, default=1.007947
             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.
+            not applied to any waters.
 
         Returns
         -------
@@ -633,10 +633,10 @@ def to_openmm_simulation(
             on a bond or angle that is fully constrained.
         ewald_tolerance : float, default=1e-4
             The value passed to `NonbondedForce.setEwaldErrorTolerance`
-        hydrogen_mass : PostitiveFloat, default=1.007947
+        hydrogen_mass : PositiveFloat, default=1.007947
             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.
+            not applied to any waters.
         additional_forces : Iterable[openmm.Force], default=tuple()
             Additional forces to be added to the system, e.g. barostats, that are not
             added by the force field.

openff/interchange/interop/openmm/__init__.py

@@ -187,11 +187,6 @@ def _apply_hmr(
     if abs(hydrogen_mass - 1.008) < 1e-3:
         return
 
-    if system.getNumParticles() != interchange.topology.n_atoms:
-        raise UnsupportedExportError(
-            "Hydrogen mass repartitioning with virtual sites present, even on rigid water, is not yet supported.",
-        )
-
     water = Molecule.from_smiles("O")
 
     def _is_water(molecule: Molecule) -> bool: