Merge pull request #906 from openforcefield/virtual-site-virtual-site-exclusions

Add (virtual site)-to-(virtual site) exceptions

Author: mattwthompson

docs/releasehistory.md

@@ -11,6 +11,12 @@ 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.21 - 2023-02-20
+
+* #906 Fixes a bug in which intramolecular interactions between virtual sites were not properly excluded with OpenMM.
+* #901 `Interchange.from_openmm` now requires the `system` argument.
+* #903 The Python API of LAMMPS is now internally used for LAMMPS energy calculations.
+
 ## 0.3.20 - 2023-02-12
 
 * #891 Adds support for hydrogen mass repartitioning (HMR) in GROMACS export. Note that this implementaiton never modifies masses in waters and requires the system contains no virtual sites.

openff/interchange/_tests/conftest.py

@@ -101,6 +101,33 @@ def sage_with_trivalent_nitrogen():
     return sage_210
 
 
+@pytest.fixture()
+def sage_with_off_center_hydrogen(sage):
+    virtual_sites = sage.get_parameter_handler("VirtualSites")
+
+    # Add a virtual site for an off-center hydrogen, see issue #905
+    # this differs by JH's example by adding an arbitrary charge increment in
+    # order to test the electrostatics of virtual site pairs that interact
+    # by 1-4 interactions (as determined by their parent relationship)
+    virtual_sites.add_parameter(
+        parameter_kwargs={
+            "smirks": "[#1:1]~[*:2]",
+            "epsilon": Quantity(0.0157, "kilocalorie / mole"),
+            "type": "BondCharge",
+            "match": "all_permutations",
+            "distance": Quantity(-0.1, "angstrom"),
+            "outOfPlaneAngle": None,
+            "inPlaneAngle": None,
+            "charge_increment1": Quantity(0.12345, "elementary_charge"),
+            "charge_increment2": Quantity(0, "elementary_charge"),
+            "sigma": Quantity(1.069078461768407, "angstrom"),
+            "name": "OFF_SITE_HYDROGEN",
+        },
+    )
+
+    return sage
+
+
 @pytest.fixture()
 def _simple_force_field():
     # TODO: Create a minimal force field for faster tests

openff/interchange/_tests/interoperability_tests/test_openmm.py

@@ -569,78 +569,6 @@ def test_valence_term_paticle_index_offsets(self, water, tip5p):
             assert isinstance(out.getVirtualSite(index), openmm.VirtualSite)
 
 
-@skip_if_missing("openmm")
-class TestOpenMMVirtualSiteExclusions:
-    def test_tip5p_num_exceptions(self, water):
-        tip5p = ForceField(get_test_file_path("tip5p.offxml"))
-
-        out = Interchange.from_smirnoff(tip5p, [water]).to_openmm(
-            combine_nonbonded_forces=True,
-        )
-
-        # In a TIP5P water    expected exceptions include (total 10)
-        #
-        # V(3)  V(4)          Oxygen to hydrogens and particles (4)
-        #    \ /                - (0, 1), (0, 2), (0, 3), (0, 4)
-        #     O(0)            Hyrogens to virtual particles (4)
-        #    / \                - (1, 3), (1, 4), (2, 3), (2, 4)
-        # H(1)  H(2)          Hydrogens and virtual particles to each other (2)
-        #                       - (1, 2), (3, 4)
-
-        for force in out.getForces():
-            if isinstance(force, openmm.NonbondedForce):
-                assert force.getNumExceptions() == 10
-
-    def test_dichloroethane_exceptions(self, sage):
-        """Test a case in which a parent's 1-4 exceptions must be 'imported'."""
-        from openff.toolkit._tests.mocking import VirtualSiteMocking
-
-        # This molecule has heavy atoms with indices (1-indexed) CL1, C2, C3, Cl4,
-        # resulting in 1-4 interactions between the Cl-Cl pair and some Cl-H pairs
-        dichloroethane = Molecule.from_mapped_smiles(
-            "[Cl:1][C:2]([H:5])([H:6])[C:3]([H:7])([H:8])[Cl:4]",
-        )
-
-        # This parameter pulls 0.1 and 0.2e from Cl (parent) and C, respectively, and has
-        # LJ parameters of 4 A, 3 kJ/mol
-        parameter = VirtualSiteMocking.bond_charge_parameter("[Cl:1]-[C:2]")
-
-        handler = VirtualSiteHandler(version="0.3")
-        handler.add_parameter(parameter=parameter)
-
-        sage.register_parameter_handler(handler)
-
-        system = Interchange.from_smirnoff(sage, [dichloroethane]).to_openmm(
-            combine_nonbonded_forces=True,
-        )
-
-        assert system.isVirtualSite(8)
-        assert system.isVirtualSite(9)
-
-        non_bonded_force = [
-            f for f in system.getForces() if isinstance(f, openmm.NonbondedForce)
-        ][0]
-
-        for exception_index in range(non_bonded_force.getNumExceptions()):
-            p1, p2, q, sigma, epsilon = non_bonded_force.getExceptionParameters(
-                exception_index,
-            )
-            if p2 == 8:
-                # Parent Cl, adjacent C and its bonded H, and the 1-3 C
-                if p1 in (0, 1, 2, 4, 5):
-                    assert q._value == epsilon._value == 0.0
-                # 1-4 Cl or 1-4 Hs
-                if p1 in (3, 6, 7):
-                    for value in (q, sigma, epsilon):
-                        assert value._value != 0, (q, sigma, epsilon)
-            if p2 == 9:
-                if p1 in (3, 1, 2, 6, 7):
-                    assert q._value == epsilon._value == 0.0
-                if p1 in (0, 4, 5):
-                    for value in (q, sigma, epsilon):
-                        assert value._value != 0, (q, sigma, epsilon)
-
-
 @skip_if_missing("openmm")
 class TestToOpenMMTopology:
     def test_num_virtual_sites(self, water, tip4p):