From b46b16bb99c7a6e689af4709c2c866f22c5f84ce Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 4 Sep 2026 12:26:30 +0100 Subject: [PATCH] Account for v-site charges when checking for charge change. --- CHANGELOG.md | 1 + src/somd2/runner/_base.py | 11 +++++++++++ tests/runner/test_alchemical_ions.py | 22 ++++++++++++++++++++++ 3 files changed, 34 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 112246d..648817d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ Changelog * Add a `precision` option for GPU platforms, defaulting to `mixed` [#191](https://github.com/OpenBioSim/somd2/pull/191). * Add support for generating Morse restraints for ring-breaking perturbations. * Remove the unused `kappa` lever equations from the ring-breaking/making lambda schedules. +* Account for off-site charges (virtual sites) when computing the charge difference between the end states. They are held as a molecule property rather than on the atoms, so a charge-preserving perturbation could appear to change charge and be given spurious alchemical ions. [2026.1.0](https://github.com/openbiosim/somd2/compare/2025.1.0...2026.1.0) - Jun 2026 -------------------------------------------------------------------------------------- diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py index 2deb608..809e189 100644 --- a/src/somd2/runner/_base.py +++ b/src/somd2/runner/_base.py @@ -1344,6 +1344,17 @@ def _get_charge_difference(system): reference = _sr.morph.link_to_reference(system).charge().value() perturbed = _sr.morph.link_to_perturbed(system).charge().value() + # Off-site charges are held as a molecule property, not on the atoms, + # so they are invisible to charge() above. + try: + vsite_mols = system.molecules("property vs_charges0") + except KeyError: + vsite_mols = [] + + for mol in vsite_mols: + reference += sum(float(x) for x in mol.property("vs_charges0")) + perturbed += sum(float(x) for x in mol.property("vs_charges1")) + return perturbed - reference def _save_alchemical_ion_indices(self, mol_indices): diff --git a/tests/runner/test_alchemical_ions.py b/tests/runner/test_alchemical_ions.py index 2c5d7fe..a255372 100644 --- a/tests/runner/test_alchemical_ions.py +++ b/tests/runner/test_alchemical_ions.py @@ -175,3 +175,25 @@ def test_alchemical_ion_abfe_schedule(schedule_name, ethane_methanol_ions): ion_schedule = schedule.get_molecule_schedule(ion_idx) assert ion_schedule.get_stages() == ["morph"] assert ion_schedule.to_string() == LambdaSchedule.standard_morph().to_string() + + +def test_charge_difference_with_virtual_sites(ethane_methanol): + """ + Off-site charges are held as a molecule property rather than on the atoms, + so they are invisible to charge() and must be added separately. + """ + + mols = ethane_methanol.clone() + + # No virtual sites, so the end states have the same charge. + assert math.isclose(Runner._get_charge_difference(mols), 0.0, abs_tol=1e-6) + + # Give the perturbable molecule off-site charges that differ by one unit + # between the end states. + mol = mols.molecules("property is_perturbable")[0] + cursor = mol.cursor() + cursor.set("vs_charges0", [0.5, 0.5]) + cursor.set("vs_charges1", [1.0, 1.0]) + mols.update(cursor.commit()) + + assert math.isclose(Runner._get_charge_difference(mols), 1.0, abs_tol=1e-6)