Merge pull request #1039 from openforcefield/openmm-reaction-field

Support reaction field electrostatics in OpenMM

docs/releasehistory.md

@@ -11,6 +11,10 @@ Dates are given in YYYY-MM-DD format.
 
 Please note that all releases prior to a version 1.0.0 are considered pre-releases and many API changes will come before a stable release.
 
+## 0.3.30 - 2024-08
+
+* #1039 Updates support of "cutoff" electrostatics in `.to_openmm` to better reflect what OpenMM supports. Set `"reaction-field"` to force the use of `CutoffPeriodic`, provided the vdW and electrostatic cutoff distances match. The potential/method `"cutoff"` is no longer supported but may be re-added in the future.
+
 ## 0.3.29 - 2024-08-01
 
 * #1023 Fixes a bug in which non-bonded parameter lookup sometimes crashed when virtual sites were present.

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

@@ -2,7 +2,10 @@
 from openff.toolkit import Molecule, unit
 from openff.utilities.testing import skip_if_missing
 
-from openff.interchange.exceptions import UnsupportedCutoffMethodError
+from openff.interchange.exceptions import (
+    UnsupportedCutoffMethodError,
+    UnsupportedExportError,
+)
 
 
 @skip_if_missing("openmm")
@@ -19,20 +22,98 @@ def test_ljpme_nonperiodic(self, sage):
             interchange.to_openmm(combine_nonbonded_forces=False)
 
     @pytest.mark.parametrize("periodic", [True, False])
-    def test_reaction_field(self, sage, periodic):
+    @pytest.mark.parametrize("combine", [True, False])
+    def test_hard_cutoff(self, sage, periodic, combine):
         interchange = sage.create_interchange(Molecule.from_smiles("CC").to_topology())
 
         if periodic:
             interchange.box = [4, 4, 4] * unit.nanometer
-            interchange["Electrostatics"].periodic_potential = "reaction-field"
+            interchange["Electrostatics"].periodic_potential = "cutoff"
         else:
-            interchange["Electrostatics"].nonperiodic_potential = "reaction-field"
+            interchange["Electrostatics"].nonperiodic_potential = "cutoff"
 
         with pytest.raises(
             UnsupportedCutoffMethodError,
-            match="Reaction field electrostatics not supported. ",
+            match="does not support.*Consider using",
         ):
-            interchange.to_openmm(combine_nonbonded_forces=False)
+            interchange.to_openmm(combine_nonbonded_forces=combine)
+
+
+@skip_if_missing("openmm")
+class TestCutoffElectrostatics:
+    @pytest.mark.parametrize("periodic", [True, False])
+    @pytest.mark.parametrize("combine", [True, False])
+    def test_export_reaction_field_electrostatics(
+        self,
+        sage,
+        basic_top,
+        periodic,
+        combine,
+    ):
+
+        import openmm
+        import openmm.unit
+
+        out = sage.create_interchange(basic_top)
+
+        if periodic:
+            out["Electrostatics"].periodic_potential = "reaction-field"
+        else:
+            out["Electrostatics"].nonperiodic_potential = "reaction-field"
+            out.box = None
+
+        assert (out.box is not None) == periodic
+
+        out["Electrostatics"].cutoff = out["vdW"].cutoff
+
+        if combine and not periodic:
+            with pytest.raises(
+                UnsupportedCutoffMethodError,
+                match="Combination of",
+            ):
+                system = out.to_openmm(combine_nonbonded_forces=combine)
+            return
+
+        system = out.to_openmm(combine_nonbonded_forces=combine)
+
+        for force in system.getForces():
+            if isinstance(force, openmm.NonbondedForce):
+
+                assert 0.9 == force.getCutoffDistance() / openmm.unit.nanometer
+
+                if periodic:
+                    expected_force = openmm.NonbondedForce.CutoffPeriodic
+                else:
+                    expected_force = openmm.NonbondedForce.CutoffNonPeriodic
+
+                assert expected_force == force.getNonbondedMethod()
+
+                del force
+
+                break
+        else:
+            pytest.fail("Found no `NonbondedForce`")
+
+    def test_vdw_electrostatics_cutoff_mismatch(
+        self,
+        sage,
+        basic_top,
+    ):
+        import random
+
+        out = sage.create_interchange(basic_top)
+
+        out["Electrostatics"].periodic_potential = "reaction-field"
+        out["Electrostatics"].cutoff = out["vdW"].cutoff
+
+        for key in ("vdW", "Electrostatics"):
+            out[key].cutoff *= random.random()
+
+        with pytest.raises(
+            UnsupportedExportError,
+            match="cutoffs must match",
+        ):
+            out.to_openmm(combine_nonbonded_forces=True)
 
 
 @skip_if_missing("openmm")

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

@@ -330,7 +330,6 @@ def test_tip5p_num_exceptions(self, water, tip5p, combine, n_molecules):
                     # Safeguard against some of the behavior seen in #919
                     for index in range(num_exceptions):
                         p1, p2, *_ = force.getExceptionParameters(index)
-                        print(p1, p2)
 
                         if sorted([p1, p2]) == [0, 3]:
                             raise Exception(

openff/interchange/interop/openmm/_nonbonded.py

@@ -3,6 +3,7 @@
 """
 
 import itertools
+import warnings
 from collections import defaultdict
 from typing import DefaultDict, NamedTuple, Optional
 
@@ -25,6 +26,7 @@
     TopologyKey,
     VirtualSiteKey,
 )
+from openff.interchange.warnings import NonbondedSettingsWarning
 
 if has_package("openmm"):
     import openmm
@@ -272,7 +274,7 @@ def _prepare_input_data(interchange: "Interchange") -> _NonbondedData:
         mixing_rule=mixing_rule,
         electrostatics_collection=electrostatics,
         electrostatics_method=electrostatics_method,
-        periodic=interchange.box is None,
+        periodic=interchange.box is not None,
     )
 
 
@@ -303,6 +305,14 @@ def _create_single_nonbonded_force(
     non_bonded_force.setName("Nonbonded force")
     system.addForce(non_bonded_force)
 
+    # This limitation is independent of periodicity and vdW method, so check it first
+    if data.electrostatics_method == "cutoff":
+        raise UnsupportedCutoffMethodError(
+            "OpenMM does not support electrostatics with a hard cutoff (and no reaction field "
+            'modification). Consider using `"reaction-field"` to get force(s) with CutoffPeriodic '
+            "or CutoffNonperiodic.",
+        )
+
     if interchange.box is None:
         if (data.vdw_method in ("no-cutoff", None)) and (
             data.electrostatics_method in ("Coulomb", None)
@@ -342,21 +352,23 @@ def _create_single_nonbonded_force(
             non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.LJPME)
             non_bonded_force.setEwaldErrorTolerance(ewald_tolerance)
 
-        elif data["vdw_method"] == data["electrostatics_method"] == "cutoff":
-            if data["vdw_cutoff"] != data["electrostatics_collection"].cutoff:
+        elif (data.vdw_method == "cutoff") and (
+            data.electrostatics_method == "reaction-field"
+        ):
+            if data.vdw_cutoff != data.electrostatics_collection.cutoff:
                 raise UnsupportedExportError(
-                    "If using cutoff vdW and electrostatics, cutoffs must match.",
+                    "If using cutoff vdW and reaction-field electrostatics, cutoffs must match.",
                 )
 
             non_bonded_force.setNonbondedMethod(openmm.NonbondedForce.CutoffPeriodic)
             non_bonded_force.setCutoffDistance(
-                to_openmm_quantity(data["vdw_cutoff"]),
+                to_openmm_quantity(data.vdw_cutoff),
             )
 
         else:
             raise UnsupportedCutoffMethodError(
                 f"Combination of non-bonded cutoff methods {data.vdw_method} (vdW) and "
-                "{data.electrostatics_method} (Electrostatics) not currently supported or invalid with "
+                f"{data.electrostatics_method} (Electrostatics) not currently supported or invalid with "
                 f"`combine_nonbonded_forces=True` and `.box={interchange.box}`.",
             )
 
@@ -735,9 +747,9 @@ def _create_multiple_nonbonded_forces(
             system.addForce(force)
 
     if vdw_force is not None and electrostatics_force is not None:
-        if (vdw_force.getNonbondedMethod() > 0) ^ (
-            electrostatics_force.getNonbondedMethod() > 0
-        ):
+        vdw_uses = vdw_force.usesPeriodicBoundaryConditions()
+        elec_uses = electrostatics_force.usesPeriodicBoundaryConditions()
+        if vdw_uses != elec_uses:
             raise UnsupportedCutoffMethodError(
                 "When using `openmm.CustomNonbondedForce`, vdW and electrostatics cutoff methods "
                 "must agree on whether or not periodic boundary conditions should be used. "
@@ -874,9 +886,24 @@ def _create_electrostatics_force(
                 ].append(force_index)
 
     if data.electrostatics_method == "reaction-field":
-        raise UnsupportedCutoffMethodError(
-            "Reaction field electrostatics not supported. If this use case is important to you, "
-            "please raise an issue describing the scope of functionality you would like to use.",
+        warnings.warn(
+            "OpenMM implements reaction field electrostatics with CutoffPeriodic and CutoffNonperiodic. "
+            "See http://docs.openmm.org/latest/userguide/theory/02_standard_forces.html?highlight="
+            "reaction%20field#coulomb-interaction-with-cutoff",
+            NonbondedSettingsWarning,
+        )
+
+        if data.periodic:
+            electrostatics_force.setNonbondedMethod(
+                openmm.NonbondedForce.CutoffPeriodic,
+            )
+        else:
+            electrostatics_force.setNonbondedMethod(
+                openmm.NonbondedForce.CutoffNonPeriodic,
+            )
+
+        electrostatics_force.setCutoffDistance(
+            to_openmm_quantity(data.electrostatics_collection.cutoff),
         )
 
     elif data.electrostatics_method == _PME:
@@ -898,8 +925,11 @@ def _create_electrostatics_force(
 
     elif data.electrostatics_method == "cutoff":
         raise UnsupportedCutoffMethodError(
-            "OpenMM does not support electrostatics with a hard cutoff.",
+            "OpenMM does not support electrostatics with a hard cutoff (and no reaction field "
+            'modification). Consider using `"reaction-field"` to get force(s) with CutoffPeriodic '
+            "or CutoffNonperiodic.",
         )
+
     else:
         raise UnsupportedCutoffMethodError(
             f"Electrostatics method {data.electrostatics_method} not supported",
@@ -912,6 +942,7 @@ def _create_electrostatics_force(
         openff_openmm_particle_map,
         parent_virtual_particle_mapping,
     )
+
     return electrostatics_force
 
 

openff/interchange/warnings.py

@@ -11,3 +11,7 @@ class SwitchingFunctionNotImplementedWarning(UserWarning):
 
 class MissingPositionsWarning(UserWarning):
     """Warning for when positions are likely needed but missing."""
+
+
+class NonbondedSettingsWarning(UserWarning):
+    """Warning for when details of nonbonded implementations get messy."""