diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 58905392e..789d30ee8 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -87,6 +87,9 @@ organisation on `GitHub `__. * Fixed ``sire.restraints.boresch_search()`` failing for charged ligands, since alchemical ions used to keep the charge constant were counted as additional perturbable molecules. +* Fixed MM non-bonded interactions not being removed between QM atoms in different + molecules, e.g. when a QM water reacts with a QM ligand. + `2026.2.0 `__ - September 2026 ---------------------------------------------------------------------------------------------- diff --git a/tests/qm/test_qm.py b/tests/qm/test_qm.py index 22f61747c..b43fddbbb 100644 --- a/tests/qm/test_qm.py +++ b/tests/qm/test_qm.py @@ -660,6 +660,76 @@ def callback(numbers_qm, charges_mm, xyz_qm, xyz_mm, cell=None, idx_mm=None): assert lever.was_force_changed("qmff"), "qmff should be changed when lambda changes" +def test_qm_cross_molecule_exceptions(ala_mols): + """ + Make sure that the MM non-bonded interactions between QM atoms in + different molecules are scaled from fully on at lambda = 0 to fully + off at lambda = 1. + """ + + import math + import openmm + + def callback(numbers_qm, charges_mm, xyz_qm, xyz_mm, cell=None, idx_mm=None): + return (0.0, xyz_qm, xyz_mm) + + mols = ala_mols.clone() + qm_mols, engine = sr.qm.create_engine( + mols, mols["molidx 0 or molidx 1"].atoms(), callback, callback=None + ) + + d = qm_mols.dynamics( + timestep="1fs", constraint="none", qm_engine=engine, platform="cpu" + ) + + ctx = d._d._omm_mols + system = ctx.getSystem() + + cljff = [f for f in system.getForces() if isinstance(f, openmm.NonbondedForce)][0] + + num_atoms0 = qm_mols[0].num_atoms() + num_atoms1 = qm_mols[1].num_atoms() + + cross_pairs = { + (i, j) + for i in range(num_atoms0) + for j in range(num_atoms0, num_atoms0 + num_atoms1) + } + + def get_exceptions(): + exceptions = {} + for k in range(cljff.getNumExceptions()): + i, j, q, sigma, epsilon = cljff.getExceptionParameters(k) + pair = (min(i, j), max(i, j)) + if pair in cross_pairs: + exceptions[pair] = ( + q.value_in_unit(openmm.unit.elementary_charge**2), + sigma.value_in_unit(openmm.unit.nanometer), + epsilon.value_in_unit(openmm.unit.kilojoule_per_mole), + ) + return exceptions + + # The context is initialised at lambda = 1, so the pairs are switched off. + exceptions = get_exceptions() + assert set(exceptions.keys()) == cross_pairs + for q, sigma, epsilon in exceptions.values(): + assert q == pytest.approx(0.0) + assert epsilon == pytest.approx(0.0, abs=1e-8) + + # At lambda = 0 the pairs should have the full MM interaction. + ctx.set_lambda(0.0) + exceptions = get_exceptions() + for (i, j), (q, sigma, epsilon) in exceptions.items(): + q0, sig0, eps0 = cljff.getParticleParameters(i) + q1, sig1, eps1 = cljff.getParticleParameters(j) + q0 = q0.value_in_unit(openmm.unit.elementary_charge) + q1 = q1.value_in_unit(openmm.unit.elementary_charge) + eps0 = eps0.value_in_unit(openmm.unit.kilojoule_per_mole) + eps1 = eps1.value_in_unit(openmm.unit.kilojoule_per_mole) + assert q == pytest.approx(q0 * q1) + assert epsilon == pytest.approx(math.sqrt(eps0 * eps1), abs=1e-8) + + def test_property_bases(ala_mols): """ Make sure that the SireOpenMM property classes are exposed as diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index d271e4a0b..3aa75bddf 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -279,7 +279,8 @@ LambdaLever::LambdaLever(const LambdaLever &other) perturbable_maps(other.perturbable_maps), lambda_cache(other.lambda_cache), last_rest2_scale(-1.0), - last_qmff_lam(-1.0) + last_qmff_lam(-1.0), + qm_exception_idxs(other.qm_exception_idxs) { } @@ -299,6 +300,7 @@ LambdaLever &LambdaLever::operator=(const LambdaLever &other) start_indices = other.start_indices; perturbable_maps = other.perturbable_maps; lambda_cache = other.lambda_cache; + qm_exception_idxs = other.qm_exception_idxs; Property::operator=(other); } @@ -312,7 +314,8 @@ bool LambdaLever::operator==(const LambdaLever &other) const lambda_schedule == other.lambda_schedule and perturbable_mols == other.perturbable_mols and start_indices == other.start_indices and - perturbable_maps == other.perturbable_maps; + perturbable_maps == other.perturbable_maps and + qm_exception_idxs == other.qm_exception_idxs; } bool LambdaLever::operator!=(const LambdaLever &other) const @@ -2047,6 +2050,47 @@ double LambdaLever::setLambda(OpenMM::Context &context, } } + // scale the MM interactions between QM atoms in different molecules + // from fully on (MM, lambda = 0) to fully off (QM, lambda = 1) + if (cljff != nullptr and not qm_exception_idxs.isEmpty()) + { + const double coul_scale = this->lambda_schedule.morph( + "clj", "charge_scale", 1.0, 0.0, lambda_value); + const double lj_scale = this->lambda_schedule.morph( + "clj", "lj_scale", 1.0, 0.0, lambda_value); + + for (const auto &idx : qm_exception_idxs) + { + int atom0, atom1; + double old_charge, old_sigma, old_epsilon; + cljff->getExceptionParameters(idx, atom0, atom1, + old_charge, old_sigma, old_epsilon); + + double q0, sig0, eps0; + double q1, sig1, eps1; + cljff->getParticleParameters(atom0, q0, sig0, eps0); + cljff->getParticleParameters(atom1, q1, sig1, eps1); + + const double charge = coul_scale * q0 * q1; + double sigma = 0.5 * (sig0 + sig1); + double epsilon = lj_scale * std::sqrt(eps0 * eps1); + + // stop OpenMM from turning this exception into an exclusion + if (charge == 0 and epsilon == 0) + { + sigma = 1e-9; + epsilon = 1e-9; + } + + if (charge != old_charge or sigma != old_sigma or epsilon != old_epsilon) + { + cljff->setExceptionParameters(idx, atom0, atom1, + charge, sigma, epsilon); + has_changed_cljff = true; + } + } + } + // update the parameters in the context for forces whose parameters changed if (has_changed_cljff and cljff) cljff->updateParametersInContext(context); @@ -2683,6 +2727,14 @@ void LambdaLever::setConstraintIndicies(int mol_idx, const QVector &cons this->perturbable_mols[mol_idx].setConstraintIndicies(constraint_idxs); } +/** Set the indices of the NonbondedForce exceptions between QM atoms + * in different molecules + */ +void LambdaLever::setQMExceptionIndicies(const QVector &exception_idxs) +{ + this->qm_exception_idxs = exception_idxs; +} + /** Return all of the property maps used to find the perturbable properties * of the perturbable molecules. This is indexed by molecule number */ diff --git a/wrapper/Convert/SireOpenMM/lambdalever.h b/wrapper/Convert/SireOpenMM/lambdalever.h index 886e436db..bea447220 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.h +++ b/wrapper/Convert/SireOpenMM/lambdalever.h @@ -139,6 +139,8 @@ namespace SireOpenMM void setConstraintIndicies(int idx, const QVector &constraint_idxs); + void setQMExceptionIndicies(const QVector &exception_idxs); + void setSchedule(const SireCAS::LambdaSchedule &schedule); void addLever(const QString &lever_name); @@ -231,6 +233,10 @@ namespace SireOpenMM * non-empty these atoms are excluded from background-lrc (their LRC * is handled by the gcmc-lrc CustomVolumeForce instead). */ QSet gcmc_water_atoms; + + /** Indices of the NonbondedForce exceptions between QM atoms + * in different molecules */ + QVector qm_exception_idxs; }; #ifndef SIRE_SKIP_INLINE_FUNCTION diff --git a/wrapper/Convert/SireOpenMM/qmmm.h b/wrapper/Convert/SireOpenMM/qmmm.h index ba662949b..5e2437bb1 100644 --- a/wrapper/Convert/SireOpenMM/qmmm.h +++ b/wrapper/Convert/SireOpenMM/qmmm.h @@ -32,6 +32,8 @@ #include "OpenMM.h" #include "openmm/Force.h" +#include + #include "sireglobal.h" #include "SireBase/property.h" @@ -54,6 +56,9 @@ namespace SireOpenMM //! Set the lambda weighting factor. virtual void setLambda(double lambda) = 0; + //! Get the indices of the QM atoms. + virtual QVector getAtoms() const = 0; + protected: virtual OpenMM::ForceImpl *createImpl() const = 0; }; @@ -74,7 +79,7 @@ namespace SireOpenMM } //! Create a QM force object. - virtual QMForce* createForce() const=0; + virtual QMForce *createForce() const = 0; //! Get a null QM engine. static const NullQMEngine &null(); @@ -96,7 +101,7 @@ namespace SireOpenMM const char *what() const; //! Create a QM force object. - QMForce* createForce() const; + QMForce *createForce() const; }; } diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 37d03840d..86b5eeac7 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -3034,6 +3034,67 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, } } + // add exceptions for QM atoms in different molecules, since these can't + // be removed via each molecule's intrascale. These are scaled by the lever. + if (qmff != 0 and cljff != 0) + { + const auto qm_atoms = qmff->getAtoms(); + + QVector qm_mol_idxs(qm_atoms.count()); + + for (int i = 0; i < qm_atoms.count(); ++i) + { + int mol_idx = 0; + + while (mol_idx + 1 < nmols and start_indexes[mol_idx + 1] <= qm_atoms[i]) + mol_idx++; + + qm_mol_idxs[i] = mol_idx; + } + + QVector qm_exception_idxs; + + for (int i = 0; i < qm_atoms.count(); ++i) + { + for (int j = i + 1; j < qm_atoms.count(); ++j) + { + if (qm_mol_idxs[i] == qm_mol_idxs[j]) + continue; + + const int atom0 = qm_atoms[i]; + const int atom1 = qm_atoms[j]; + + double q0, sig0, eps0; + double q1, sig1, eps1; + cljff->getParticleParameters(atom0, q0, sig0, eps0); + cljff->getParticleParameters(atom1, q1, sig1, eps1); + + const double charge = q0 * q1; + double sigma = 0.5 * (sig0 + sig1); + double epsilon = std::sqrt(eps0 * eps1); + + // stop OpenMM from turning this exception into an exclusion + if (charge == 0 and epsilon == 0) + { + sigma = 1e-9; + epsilon = 1e-9; + } + + qm_exception_idxs.append(cljff->addException(atom0, atom1, + charge, sigma, + epsilon, false)); + + if (ghost_ghostff != 0) + { + ghost_ghostff->addExclusion(atom0, atom1); + ghost_nonghostff->addExclusion(atom0, atom1); + } + } + } + + lambda_lever.setQMExceptionIndicies(qm_exception_idxs); + } + // Stage 6 is complete. We have set up all of the exceptions. The // total energy / force calculated for the system should now be // correct.