TST: Test interaction of preset charges and virtual sites

openff/interchange/_tests/conftest.py

@@ -84,7 +84,7 @@ def sage_with_sigma_hole(sage):
             smirks="[#6:1]-[#17:2]",
             distance=1.4 * unit.angstrom,
             type="BondCharge",
-            match="once",
+            match="all_permutations",
             charge_increment1=0.1 * unit.elementary_charge,
             charge_increment2=0.2 * unit.elementary_charge,
         ),

openff/interchange/_tests/unit_tests/smirnoff/test_create.py

@@ -22,6 +22,7 @@
     _upconvert_vdw_handler,
     library_charge_from_molecule,
 )
+from openff.interchange.warnings import PresetChargesAndVirtualSitesWarning
 
 
 def _get_interpolated_bond_k(bond_handler) -> float:
@@ -264,6 +265,76 @@ def test_error_multiple_isomorphic_molecules(self, sage, ethanol, reversed_ethan
                 charge_from_molecules=[ethanol, reversed_ethanol],
             )
 
+    def test_virtual_site_charge_increments_applied_after_preset_charges_water(
+        self,
+        tip4p,
+        water,
+    ):
+        """
+        Test that charge increments to/from virtual sites are still applied after preset charges are set.
+
+        This is funky user input, since preset charges override library charges, which are important for water.
+        """
+        water.partial_charges = Quantity(
+            [-2.0, 1.0, 1.0],
+            "elementary_charge",
+        )
+
+        with pytest.warns(PresetChargesAndVirtualSitesWarning):
+            out = tip4p.create_interchange(
+                water.to_topology(),
+                charge_from_molecules=[water],
+            )
+
+        # This dict is probably ordered O H H VS
+        charges = out["Electrostatics"].charges
+
+        # tip4p_fb.offxml is meant to result in charges of -1.0517362213526 (VS) 0 (O) and 0.5258681106763 (H)
+        # the force field has 0.0 for all library charges (skipping AM1-BCC), so preset charges screw this up.
+        # the resulting charges, if using preset charges of [-2, 1, 1], should be the result of adding together
+        # [-2, 1, 1, 0] (preset charges) and
+        # [0, 1.5258681106763001, 1.5258681106763001, -1.0517362213526] (charge increments)
+
+        numpy.testing.assert_allclose(
+            [charge.m for charge in charges.values()],
+            [-2.0, 1.5258681106763001, 1.5258681106763001, -1.0517362213526],
+        )
+
+        # The cache seems to leak charges between tests
+        out["Electrostatics"]._charges.clear()
+
+    def test_virtual_site_charge_increments_applied_after_preset_charges_ligand(
+        self,
+        sage_with_sigma_hole,
+    ):
+        ligand = Molecule.from_mapped_smiles("[C:1]([Cl:2])([Cl:3])([Cl:4])[Cl:5]")
+
+        ligand.partial_charges = Quantity(
+            [1.2, -0.3, -0.3, -0.3, -0.3],
+            "elementary_charge",
+        )
+
+        out = sage_with_sigma_hole.create_interchange(
+            ligand.to_topology(),
+            charge_from_molecules=[ligand],
+        )
+
+        # This dict is probably ordered C Cl Cl Cl Cl VS VS VS VS
+        charges = out["Electrostatics"].charges
+
+        # this fake sigma hole parameter shifts 0.1 e from the carbon and 0.2 e from the chlorine, so the
+        # resulting charges should be like
+        # [1.2, -0.3, -0.3, -0.3, -0.3] +
+        # [0.4, 0.2, 0.2, 0.2, 0.2, -0.3, -0.3, -0.3, -0.3]
+
+        numpy.testing.assert_allclose(
+            [charge.m for charge in charges.values()],
+            [1.6, -0.1, -0.1, -0.1, -0.1, -0.3, -0.3, -0.3, -0.3],
+        )
+
+        # The cache seems to leak charges between tests
+        out["Electrostatics"]._charges.clear()
+
 
 @skip_if_missing("openmm")
 class TestPartialBondOrdersFromMolecules:

openff/interchange/smirnoff/_create.py

@@ -1,3 +1,5 @@
+import warnings
+
 from openff.toolkit import ForceField, Molecule, Quantity, Topology
 from openff.toolkit.typing.engines.smirnoff import ParameterHandler
 from openff.toolkit.typing.engines.smirnoff.plugins import load_handler_plugins
@@ -26,6 +28,7 @@
     SMIRNOFFProperTorsionCollection,
 )
 from openff.interchange.smirnoff._virtual_sites import SMIRNOFFVirtualSiteCollection
+from openff.interchange.warnings import PresetChargesAndVirtualSitesWarning
 
 _SUPPORTED_PARAMETER_HANDLERS: set[str] = {
     "Constraints",
@@ -141,6 +144,14 @@ def _create_interchange(
 
     _check_supported_handlers(force_field)
 
+    if molecules_with_preset_charges is not None and "VirtualSites" in force_field.registered_parameter_handlers:
+        warnings.warn(
+            "Preset charges were provided (via `charge_from_molecules`) alongside a force field that includes "
+            "virtual site parameters. Note that virtual sites will be applied charges from the force field and "
+            "cannot be given preset charges.",
+            PresetChargesAndVirtualSitesWarning,
+        )
+
     # interchange = Interchange(topology=topology)
     # or maybe
     interchange = Interchange(topology=validate_topology(topology))

openff/interchange/warnings.py

@@ -21,5 +21,9 @@ class ForceFieldModificationWarning(UserWarning):
     """Warning for when a ForceField is modified."""
 
 
+class PresetChargesAndVirtualSitesWarning(UserWarning):
+    """Warning when possibly using preset charges and virtual sites together."""
+
+
 class InterchangeCombinationWarning(UserWarning):
     """Warning for when combining Interchange objects."""