diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index f2d0be394..86aade73b 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -76,6 +76,18 @@ organisation on `GitHub `__. structures. The indices are membership tested against every exception and torsion in the system, which was quadratic for large REST2 regions, e.g. proteins. +* Fixed the REST2 region mask not covering off-site charges (virtual sites), which + are appended to the OpenMM system as extra particles after the atoms of each + molecule. This caused an out-of-bounds read when applying the REST2 scaling to a + perturbable molecule with virtual sites. Each virtual site now inherits the REST2 + flag of its parent atom, so its charge is scaled along with it. + +* Fixed off-site charges (virtual sites) on non-ghost atoms of a perturbable molecule + being added to neither the ghost nor the non-ghost interaction group of the softcore + forces. Their interaction with any ghost atom was evaluated by the standard + ``NonbondedForce`` alone, i.e. without softening, and without the subtraction of the + coulomb energy that the softcore replaces. + `2026.1.0 `__ - June 2026 ----------------------------------------------------------------------------------------- diff --git a/tests/convert/test_openmm_vsites.py b/tests/convert/test_openmm_vsites.py index fd1208c58..9a331aa22 100644 --- a/tests/convert/test_openmm_vsites.py +++ b/tests/convert/test_openmm_vsites.py @@ -1,3 +1,5 @@ +from math import isclose + import sire as sr import pytest @@ -152,6 +154,199 @@ def test_vsite_pertubation(ethane_12dichloroethane, openmm_platform): assert expected_charge == nb_charge.value_in_unit(unit.elementary_charge) +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +@pytest.mark.parametrize( + ["rest2_selection", "rest2_vsites"], + [ + # No selection, so the whole perturbable molecule is the REST2 region + # and both virtual sites are scaled. + (None, [0, 1]), + # Only the parent atom of the first virtual site is in the REST2 + # region, so only that virtual site is scaled. + ("atomidx 0,1,2", [0]), + ], +) +def test_vsite_rest2( + ethane_12dichloroethane, openmm_platform, rest2_selection, rest2_vsites +): + # Do virtual sites inherit the REST2 flag of their parent atom? + mols = ethane_12dichloroethane + + # Just dichloroethane + mol = mols[0] + + # Set vsite properties. The parent atoms are 0 and 3, i.e. the first + # index of each vs_indices list. + vsite_dict = { + "0": { + "vs_indices": [0, 1, 2], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + }, + "1": { + "vs_indices": [3, 2, 1], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + }, + } + + parents_dict = {str(atom_i): [] for atom_i in range(mol.num_atoms())} + for v, vs in enumerate(vsite_dict): + parent = vsite_dict[vs]["vs_indices"][0] + parents_dict[str(parent)].append(v) + + n_virtual_sites = len(vsite_dict) + vs_charges0 = [0.1, 0.1] + vs_charges1 = [0.2, 0.2] + + cursor = mol.cursor() + cursor.set("n_virtual_sites", n_virtual_sites) + cursor.set("vs_charges0", vs_charges0) + cursor.set("vs_charges1", vs_charges1) + cursor.set("virtual_sites", vsite_dict) + cursor.set("parents", parents_dict) + mol = cursor.commit() + + mol = sr.morph.link_to_reference(mol) + + d = mol.dynamics( + lambda_value=0.0, + platform=openmm_platform, + rest2_selection=rest2_selection, + ) + + nb_force = next( + force + for force in d.context().getSystem().getForces() + if force.getName() == "NonbondedForce" + ) + + # Store the unscaled charges at the same lambda value, so that the + # comparison isolates the REST2 scaling from the lambda lever. + d.set_lambda(0.0, rest2_scale=1.0) + charges_initial = [ + nb_force.getParticleParameters(i)[0]._value + for i in range(nb_force.getNumParticles()) + ] + + # Update the REST2 scaling factor. + d.set_lambda(0.0, rest2_scale=2.0) + + # Store the scaling factor. + scale = 0.5 + + for vs_index in range(n_virtual_sites): + parent = vsite_dict[str(vs_index)]["vs_indices"][0] + + # The virtual sites are appended after the atoms of the molecule. + for i in [parent, mol.num_atoms() + vs_index]: + charge = nb_force.getParticleParameters(i)[0]._value + + if vs_index in rest2_vsites: + assert isclose(charge, charges_initial[i] * scale**0.5) + else: + assert isclose(charge, charges_initial[i]) + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_vsite_ghost_interaction_groups(solvated_neopentane_methane, openmm_platform): + # Are virtual sites placed in the same softcore interaction group as their + # parent atom? + mols = solvated_neopentane_methane + + mol0 = mols[0] + + # Atom 0 (C1) becomes a ghost at lambda = 1, while atom 1 (C2) is present + # in both end states, so the two virtual sites cover both cases. + vsite_dict = { + "0": { + "vs_indices": [0, 1, 2], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + }, + "1": { + "vs_indices": [1, 2, 3], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + }, + } + + parents_dict = {str(atom_i): [] for atom_i in range(mol0.num_atoms())} + for v, vs in enumerate(vsite_dict): + parent = vsite_dict[vs]["vs_indices"][0] + parents_dict[str(parent)].append(v) + + n_virtual_sites = len(vsite_dict) + vs_charges0 = [0.1, 0.1] + vs_charges1 = [0.2, 0.2] + + cursor = mol0.cursor() + cursor.set("n_virtual_sites", n_virtual_sites) + cursor.set("vs_charges0", vs_charges0) + cursor.set("vs_charges1", vs_charges1) + cursor.set("virtual_sites", vsite_dict) + cursor.set("parents", parents_dict) + mol0 = cursor.commit() + mols.update(mol0) + + mols = sr.morph.link_to_reference(mols) + + d = mols.dynamics(lambda_value=0.0, platform=openmm_platform) + system = d.context().getSystem() + + forces = {force.getName(): force for force in system.getForces()} + + # The ghost/ghost group is the ghost atoms paired with themselves, and the + # ghost/non-ghost group is the ghost atoms paired with everything else. + ghosts, ghosts_again = forces[ + "GhostGhostNonbondedForce" + ].getInteractionGroupParameters(0) + assert set(ghosts) == set(ghosts_again) + + ghosts_check, non_ghosts = forces[ + "GhostNonGhostNonbondedForce" + ].getInteractionGroupParameters(0) + assert set(ghosts_check) == set(ghosts) + + ghosts = set(ghosts) + non_ghosts = set(non_ghosts) + + # Every particle, including every virtual site, must be in exactly one of + # the two groups. A particle in neither would be left hard. + assert ghosts.isdisjoint(non_ghosts) + assert ghosts | non_ghosts == set(range(system.getNumParticles())) + + # The virtual sites are appended after the atoms of the molecule, which is + # the first molecule in the system. + for vs_index in range(n_virtual_sites): + parent = vsite_dict[str(vs_index)]["vs_indices"][0] + particle = mol0.num_atoms() + vs_index + + if parent in ghosts: + assert particle in ghosts + else: + assert particle in non_ghosts + + # Make sure that both cases were actually covered, i.e. that the choice of + # parent atoms above still holds for this perturbation. + assert vsite_dict["0"]["vs_indices"][0] in ghosts + assert vsite_dict["1"]["vs_indices"][0] in non_ghosts + + @pytest.mark.skipif( "openmm" not in sr.convert.supported_formats(), reason="openmm support is not available", diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 78b53999f..8d2d66f84 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -549,6 +549,24 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, is_rest2 = QVector(nats, true); } + // Virtual sites are appended as extra particles after the atoms, so the + // mask must cover them too. Each virtual site inherits the REST2 flag of + // its parent atom. + if (this->has_vs) + { + is_rest2.resize(nats + this->n_vs); + + for (int i = 0; i < nats; ++i) + { + auto atom_vs = this->vs_parents.property(std::to_string(i).c_str()).asAnArray(); + + for (int vs = 0; vs < atom_vs.size(); ++vs) + { + is_rest2[nats + atom_vs.at(vs).asAnInteger()] = is_rest2[i]; + } + } + } + if (nats <= 0) { return; diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 3116ecfab..ec91683f8 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -2274,6 +2274,13 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, ghost_atoms.append(atom_index); from_ghost_idxs.append(atom_index); } + else + { + // the parent isn't a ghost, so this virtual site must + // join the non-ghost group, else its interaction with + // any ghost atom would be left hard + non_ghost_atoms.append(atom_index); + } } else if (any_perturbable) {