From d00f7655e68b3039bbbd0bbc8def2cab643f871d Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 8 May 2026 10:27:09 +0100 Subject: [PATCH 1/4] Add softcore CustomBondForce for ring-breaking and ring-making pairs. --- doc/source/changelog.rst | 2 + wrapper/Convert/SireOpenMM/lambdalever.cpp | 90 +++++++ wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 30 +++ wrapper/Convert/SireOpenMM/openmmmolecule.h | 24 +- .../SireOpenMM/sire_to_openmm_system.cpp | 223 +++++++++++++++++- 5 files changed, 362 insertions(+), 7 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index c71a60ee9..053f4b6a8 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -86,6 +86,8 @@ organisation on `GitHub `__. * Add support for using a switching function for QM/MM simulations. +* Add softcore ``CustomBondForce`` for ring-breaking and ring-making pairs. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 1f8023ee1..ab780c13c 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1299,6 +1299,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, auto ghost_ghostff = this->getForce("ghost/ghost", system); auto ghost_nonghostff = this->getForce("ghost/non-ghost", system); auto ghost_14ff = this->getForce("ghost-14", system); + auto ring_breaking_ff = this->getForce("ring-break", system); + auto ring_making_ff = this->getForce("ring-make", system); auto bondff = this->getForce("bond", system); auto angff = this->getForce("angle", system); auto dihff = this->getForce("torsion", system); @@ -1330,6 +1332,25 @@ double LambdaLever::setLambda(OpenMM::Context &context, bool has_changed_dihff = false; bool has_changed_cmap = false; + // Pre-compute ring-break/make kappa values so the per-mol exception update + // loop and the later global-parameter block both use the same values. + const double rb_kappa = (ring_breaking_ff != nullptr) + ? std::max(0.0, std::min(1.0, this->lambda_schedule.morph( + "ring-break", "kappa", 0.0, 1.0, lambda_value))) + : 0.0; + const double rb_alpha = (ring_breaking_ff != nullptr) + ? std::max(0.0, std::min(1.0, this->lambda_schedule.morph( + "ring-break", "alpha", 1.0, 0.0, lambda_value))) + : 1.0; + const double rm_kappa = (ring_making_ff != nullptr) + ? std::max(0.0, std::min(1.0, this->lambda_schedule.morph( + "ring-make", "kappa", 1.0, 0.0, lambda_value))) + : 1.0; + const double rm_alpha = (ring_making_ff != nullptr) + ? std::max(0.0, std::min(1.0, this->lambda_schedule.morph( + "ring-make", "alpha", 0.0, 1.0, lambda_value))) + : 0.0; + // change the parameters for all of the perturbable molecules for (int i = 0; i < this->perturbable_mols.count(); ++i) { @@ -1658,6 +1679,12 @@ double LambdaLever::setLambda(OpenMM::Context &context, scale = rest2_scale; } + // ring-breaking/making pairs have idx=-1 (their CLJ + // exception is fixed at 1e-9 and must not be updated; + // the ring force handles the interaction via global params) + if (boost::get<0>(idxs[j]) == -1) + continue; + // don't set LJ terms for ghost atoms if (atom0_is_ghost or atom1_is_ghost) { @@ -1725,6 +1752,45 @@ double LambdaLever::setLambda(OpenMM::Context &context, } } + // Update CLJ exceptions for ring-breaking pairs: scale charge product by + // rb_kappa (0 at bonded end, 1 at nonbonded end). LJ stays at 1e-9 + // so the ring-break CustomBondForce provides the full LJ correction. + if (cljff != nullptr and ring_breaking_ff != nullptr) + { + const auto rb_exc = perturbable_mol.getExceptionIndicies("ring-break"); + for (int j = 0; j < rb_exc.count(); ++j) + { + const int clj_idx = boost::get<0>(rb_exc[j]); + if (clj_idx < 0) + continue; + const int bond_idx = boost::get<1>(rb_exc[j]); + int a0, a1; + std::vector bp; + ring_breaking_ff->getBondParameters(bond_idx, a0, a1, bp); + cljff->setExceptionParameters(clj_idx, a0, a1, + rb_kappa * bp[0], 1e-9, 1e-9); + has_changed_cljff = true; + } + } + + if (cljff != nullptr and ring_making_ff != nullptr) + { + const auto rm_exc = perturbable_mol.getExceptionIndicies("ring-make"); + for (int j = 0; j < rm_exc.count(); ++j) + { + const int clj_idx = boost::get<0>(rm_exc[j]); + if (clj_idx < 0) + continue; + const int bond_idx = boost::get<1>(rm_exc[j]); + int a0, a1; + std::vector bp; + ring_making_ff->getBondParameters(bond_idx, a0, a1, bp); + cljff->setExceptionParameters(clj_idx, a0, a1, + rm_kappa * bp[0], 1e-9, 1e-9); + has_changed_cljff = true; + } + } + // update all of the perturbable constraints if (update_constraints) { @@ -2054,6 +2120,28 @@ double LambdaLever::setLambda(OpenMM::Context &context, context.setParameter("lrc_scale", lrc_scale); } + // Update ring-breaking/making softcore force global parameters using the + // values pre-computed before the per-mol loop (rb_alpha/kappa, rm_alpha/kappa). + bool has_changed_ring_breaking_ff = false; + if (ring_breaking_ff != nullptr) + { + has_changed_ring_breaking_ff = + (rb_alpha != context.getParameter("ring_break_alpha") || + rb_kappa != context.getParameter("ring_break_kappa")); + context.setParameter("ring_break_alpha", rb_alpha); + context.setParameter("ring_break_kappa", rb_kappa); + } + + bool has_changed_ring_making_ff = false; + if (ring_making_ff != nullptr) + { + has_changed_ring_making_ff = + (rm_alpha != context.getParameter("ring_make_alpha") || + rm_kappa != context.getParameter("ring_make_kappa")); + context.setParameter("ring_make_alpha", rm_alpha); + context.setParameter("ring_make_kappa", rm_kappa); + } + // Update the NonbondedForce (background) LRC via its own CustomVolumeForce. // Ghost atoms have epsilon=0 in cljff so they contribute nothing naturally. auto background_lrc_ff = this->getForce("background-lrc", system); @@ -2179,6 +2267,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, last_changed_forces["background-lrc"] = has_changed_cljff; last_changed_forces["gcmc-lrc"] = false; last_changed_forces["ghost-14"] = has_changed_ghost14ff; + last_changed_forces["ring-break"] = has_changed_ring_breaking_ff; + last_changed_forces["ring-make"] = has_changed_ring_making_ff; last_changed_forces["bond"] = has_changed_bondff; last_changed_forces["angle"] = has_changed_angff; last_changed_forces["torsion"] = has_changed_dihff; diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index e0c0e8133..aa19a231f 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -294,6 +294,36 @@ OpenMMMolecule::OpenMMMolecule(const Molecule &mol, std::swap(map0, map1); } + // Read ring-breaking/making bond pairs from molecule properties, + // swapping them if end states are inverted so that the members + // always reflect the λ=0/λ=1 convention of the (possibly swapped) + // end states. + auto read_ring_pairs = [&](const QString &prop_name) + { + QVector> pairs; + if (mol.hasProperty(prop_name)) + { + const auto &flat = mol.property(prop_name) + .asA() + .toVector(); + pairs.reserve(flat.count() / 2); + for (int k = 0; k + 1 < flat.count(); k += 2) + pairs.append(QPair(flat[k], flat[k + 1])); + } + return pairs; + }; + + if (swap_end_states) + { + this->ring_breaking_pairs = read_ring_pairs("ring_making_bonds"); + this->ring_making_pairs = read_ring_pairs("ring_breaking_bonds"); + } + else + { + this->ring_breaking_pairs = read_ring_pairs("ring_breaking_bonds"); + this->ring_making_pairs = read_ring_pairs("ring_making_bonds"); + } + // save this perturbable map - this will help us set // new properties from the results of dynamics, e.g. // updating coordinates after minimisation diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.h b/wrapper/Convert/SireOpenMM/openmmmolecule.h index 79af7a268..39a1e2621 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.h +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.h @@ -3,17 +3,17 @@ #include -#include "SireMol/moleculeinfo.h" -#include "SireMol/core.h" #include "SireMol/atom.h" +#include "SireMol/core.h" +#include "SireMol/moleculeinfo.h" #include "SireMol/selector.hpp" -#include "SireMM/mmdetail.h" -#include "SireMM/excludedpairs.h" #include "SireMM/amberparams.h" -#include "SireMM/bond.h" #include "SireMM/angle.h" +#include "SireMM/bond.h" #include "SireMM/dihedral.h" +#include "SireMM/excludedpairs.h" +#include "SireMM/mmdetail.h" #include @@ -34,7 +34,7 @@ namespace SireOpenMM }; Type type; - int vsite_idx; // molecule-local index of the virtual site atom + int vsite_idx; // molecule-local index of the virtual site atom int p1_idx, p2_idx, p3_idx; // molecule-local parent indices (O, H1, H2) // Weights for ThreeParticleAverageSite: pos = w1*p1 + w2*p2 + w3*p3 @@ -208,6 +208,18 @@ namespace SireOpenMM */ QSet from_ghost_idxs; + /** Pairs of atom indices (molecule-local) whose bond is present + * at λ=0 but absent at λ=1 (ring-breaking). Already swapped if + * swap_end_states was active at construction time. + */ + QVector> ring_breaking_pairs; + + /** Pairs of atom indices (molecule-local) whose bond is absent + * at λ=0 but present at λ=1 (ring-making). Already swapped if + * swap_end_states was active at construction time. + */ + QVector> ring_making_pairs; + /** What type of constraint to use */ qint32 constraint_type; diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 642197e37..3562fb95c 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1125,6 +1125,9 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, } // check to see if there are any perturbable molecules + bool any_ring_breaking = false; + bool any_ring_making = false; + if (not ignore_perturbations) { for (int i = 0; i < nmols; ++i) @@ -1132,8 +1135,13 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, if (openmm_mols_data[i].isPerturbable()) { any_perturbable = true; - break; } + + if (not openmm_mols_data[i].ring_breaking_pairs.isEmpty()) + any_ring_breaking = true; + + if (not openmm_mols_data[i].ring_making_pairs.isEmpty()) + any_ring_making = true; } } @@ -1308,6 +1316,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, /// OpenMM::CustomBondForce *ghost_14ff = 0; + OpenMM::CustomBondForce *ring_breaking_ff = 0; + OpenMM::CustomBondForce *ring_making_ff = 0; OpenMM::CustomNonbondedForce *ghost_ghostff = 0; OpenMM::CustomNonbondedForce *ghost_nonghostff = 0; @@ -1449,6 +1459,114 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, .toStdString(); } + // Ring-breaking/making softcore expressions: same functional form as + // ghost-14 but using global parameters so a single alpha/kappa value + // is shared across all bonds in each force and can be driven by the + // schedule without per-bond tracking infrastructure. + std::string rb_expression, rm_expression; + const bool need_rb = any_ring_breaking or any_ring_making; + + // ring_break_kappa prefactor on Coulomb ensures zero interaction at the + // bonded end (kappa=0) and a clean handoff to the CLJ exception at the + // nonbonded end (kappa=1, softcore correction = kappa*(1/r - kappa/r) = 0). + // ring_make_kappa plays the same role for the opposite direction. + if (need_rb and use_beutler_softening) + { + rb_expression = QString( + "coul_nrg+lj_nrg;" + "coul_nrg=ring_break_kappa*138.9354558466661*q*((1/sqrt((%1*ring_break_alpha)+r_safe^2))-(ring_break_kappa/r_safe));" + "lj_nrg=(1-ring_break_alpha)*four_epsilon*sig6*(sig6-1);" + "sig6=(sigma^6)/(%2*sigma^6*ring_break_alpha + r_safe^6);" + "r_safe=max(r, 0.001);") + .arg(shift_coulomb) + .arg(beutler_alpha) + .toStdString(); + rm_expression = QString( + "coul_nrg+lj_nrg;" + "coul_nrg=ring_make_kappa*138.9354558466661*q*((1/sqrt((%1*ring_make_alpha)+r_safe^2))-(ring_make_kappa/r_safe));" + "lj_nrg=(1-ring_make_alpha)*four_epsilon*sig6*(sig6-1);" + "sig6=(sigma^6)/(%2*sigma^6*ring_make_alpha + r_safe^6);" + "r_safe=max(r, 0.001);") + .arg(shift_coulomb) + .arg(beutler_alpha) + .toStdString(); + } + else if (use_taylor_softening) + { + rb_expression = QString( + "coul_nrg+lj_nrg;" + "coul_nrg=ring_break_kappa*138.9354558466661*q*((1/sqrt((%1*ring_break_alpha)+r_safe^2))-(ring_break_kappa/r_safe));" + "lj_nrg=four_epsilon*sig6*(sig6-1);" + "sig6=(sigma^6)/(%2*sigma^6 + r_safe^6);" + "r_safe=max(r, 0.001);") + .arg(shift_coulomb) + .arg(taylor_power_expression("ring_break_alpha", taylor_power)) + .toStdString(); + rm_expression = QString( + "coul_nrg+lj_nrg;" + "coul_nrg=ring_make_kappa*138.9354558466661*q*((1/sqrt((%1*ring_make_alpha)+r_safe^2))-(ring_make_kappa/r_safe));" + "lj_nrg=four_epsilon*sig6*(sig6-1);" + "sig6=(sigma^6)/(%2*sigma^6 + r_safe^6);" + "r_safe=max(r, 0.001);") + .arg(shift_coulomb) + .arg(taylor_power_expression("ring_make_alpha", taylor_power)) + .toStdString(); + } + else + { + rb_expression = QString( + "coul_nrg+lj_nrg;" + "coul_nrg=ring_break_kappa*138.9354558466661*q*((1/sqrt((%1*ring_break_alpha)+r_safe^2))-(ring_break_kappa/r_safe));" + "lj_nrg=four_epsilon*sig6*(sig6-1);" + "sig6=(sigma^6)/(((sigma*delta)+r_safe^2)^3);" + "r_safe=max(r, 0.001);" + "delta=%2*ring_break_alpha;") + .arg(shift_coulomb) + .arg(shift_delta.to(SireUnits::nanometer)) + .toStdString(); + rm_expression = QString( + "coul_nrg+lj_nrg;" + "coul_nrg=ring_make_kappa*138.9354558466661*q*((1/sqrt((%1*ring_make_alpha)+r_safe^2))-(ring_make_kappa/r_safe));" + "lj_nrg=four_epsilon*sig6*(sig6-1);" + "sig6=(sigma^6)/(((sigma*delta)+r_safe^2)^3);" + "r_safe=max(r, 0.001);" + "delta=%2*ring_make_alpha;") + .arg(shift_coulomb) + .arg(shift_delta.to(SireUnits::nanometer)) + .toStdString(); + } + + // ring_break_alpha=1 initially: fully soft at the bonded (ring-closed) + // end state so the pair interaction grows from zero as lambda moves into + // the morph stage. ring_break_kappa=0 initially: no coulomb contribution + // at the bonded end (pair is excluded there). + if (any_ring_breaking) + { + ring_breaking_ff = new OpenMM::CustomBondForce(rb_expression); + ring_breaking_ff->setName("RingBreakingBondForce"); + ring_breaking_ff->addGlobalParameter("ring_break_alpha", 1.0); + ring_breaking_ff->addGlobalParameter("ring_break_kappa", 0.0); + ring_breaking_ff->addPerBondParameter("q"); + ring_breaking_ff->addPerBondParameter("sigma"); + ring_breaking_ff->addPerBondParameter("four_epsilon"); + ring_breaking_ff->setUsesPeriodicBoundaryConditions(false); + } + + // ring_make_alpha=0 initially: hard at the nonbonded (ring-open) end + // so the pair interacts normally there. ring_make_kappa=1 initially: + // full coulomb at the nonbonded end. + if (any_ring_making) + { + ring_making_ff = new OpenMM::CustomBondForce(rm_expression); + ring_making_ff->setName("RingMakingBondForce"); + ring_making_ff->addGlobalParameter("ring_make_alpha", 0.0); + ring_making_ff->addGlobalParameter("ring_make_kappa", 1.0); + ring_making_ff->addPerBondParameter("q"); + ring_making_ff->addPerBondParameter("sigma"); + ring_making_ff->addPerBondParameter("four_epsilon"); + ring_making_ff->setUsesPeriodicBoundaryConditions(false); + } + ghost_14ff = new OpenMM::CustomBondForce(nb14_expression); ghost_14ff->setName("Ghost14BondForce"); @@ -1626,6 +1744,20 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, ghost_14ff->setForceGroup(force_group_counter); lambda_lever.setForceIndex("ghost-14", system.addForce(ghost_14ff)); lambda_lever.setForceGroup("ghost-14", force_group_counter++); + + if (ring_breaking_ff != 0) + { + ring_breaking_ff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ring-break", system.addForce(ring_breaking_ff)); + lambda_lever.setForceGroup("ring-break", force_group_counter++); + } + + if (ring_making_ff != 0) + { + ring_making_ff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ring-make", system.addForce(ring_making_ff)); + lambda_lever.setForceGroup("ring-make", force_group_counter++); + } } // Analytic LRC for the NonbondedForce (all non-ghost atoms): E = lrc_background / V, @@ -2410,6 +2542,27 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, const bool is_perturbable = any_perturbable and mol.isPerturbable(); + // Build local sets of ring-breaking/making pairs (molecule-local + // indices) for fast lookup during exception processing. + QSet rb_pairs_local, rm_pairs_local; + // Per-exception index arrays for ring-break/make CLJ exceptions. + // Filled alongside exception_idxs; (-1,-1) for non-ring pairs. + QVector> rb_exception_idxs; + QVector> rm_exception_idxs; + int rb_bond_count = 0; + int rm_bond_count = 0; + if (is_perturbable) + { + for (const auto &p : mol.ring_breaking_pairs) + rb_pairs_local.insert(IndexPair(p.first, p.second)); + for (const auto &p : mol.ring_making_pairs) + rm_pairs_local.insert(IndexPair(p.first, p.second)); + rb_exception_idxs = QVector>( + mol.exception_params.count(), boost::make_tuple(-1, -1)); + rm_exception_idxs = QVector>( + mol.exception_params.count(), boost::make_tuple(-1, -1)); + } + if (is_perturbable) { exception_idxs = QVector>(mol.exception_params.count(), @@ -2461,6 +2614,9 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, int idx = -1; int nbidx = -1; + const bool is_ring_breaking = rb_pairs_local.contains(IndexPair(atom0, atom1)); + const bool is_ring_making = rm_pairs_local.contains(IndexPair(atom0, atom1)); + if (atom0_is_ghost or atom1_is_ghost) { // don't include the LJ term, as this is calculated @@ -2506,6 +2662,65 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, excluded_ghost_pairs.insert(IndexPair(boost::get<0>(p), boost::get<1>(p))); } } + else if (is_ring_breaking or is_ring_making) + { + // Add CLJ exception with near-zero charges initially; the + // lambda lever will scale these by rb_kappa / rm_kappa each + // step, growing the hard Coulomb as the ring opens/closes. + // LJ stays at 1e-9: the ring-break/make CustomBondForce + // provides the full softcore LJ. + idx = cljff->addException(boost::get<0>(p), boost::get<1>(p), + 1e-9, 1e-9, 1e-9, true); + // nbidx stays -1 (no ghost-14 bond for ring-break pairs) + + if (is_ring_breaking and ring_breaking_ff != 0) + { + // CLJ parameters from the nonbonded end state (λ=1, + // perturbed), where the bond is absent and the pair + // interacts normally (full scale factors). + auto pp = mol.perturbed->getException( + atom0, atom1, start_index, 1.0, 1.0); + std::vector params_rb = { + boost::get<2>(pp), + boost::get<3>(pp), + 4.0 * boost::get<4>(pp)}; + if (params_rb[0] == 0) + params_rb[0] = 1e-9; + if (params_rb[1] == 0) + params_rb[1] = 1e-9; + ring_breaking_ff->addBond(boost::get<0>(p), + boost::get<1>(p), + params_rb); + rb_exception_idxs[j] = boost::make_tuple(idx, rb_bond_count); + ++rb_bond_count; + } + else if (is_ring_making and ring_making_ff != 0) + { + // CLJ parameters from the nonbonded end state (λ=0, + // reference), where the bond is absent. + auto pp = mol.getException( + atom0, atom1, start_index, 1.0, 1.0); + std::vector params_rm = { + boost::get<2>(pp), + boost::get<3>(pp), + 4.0 * boost::get<4>(pp)}; + if (params_rm[0] == 0) + params_rm[0] = 1e-9; + if (params_rm[1] == 0) + params_rm[1] = 1e-9; + ring_making_ff->addBond(boost::get<0>(p), + boost::get<1>(p), + params_rm); + rm_exception_idxs[j] = boost::make_tuple(idx, rm_bond_count); + ++rm_bond_count; + } + + // Store idx in exception_idxs so the main loop guard + // (`if (boost::get<0>(idxs[j]) == -1) continue`) still + // skips these pairs — we handle them separately below. + // Reset idx to -1 so the guard works correctly. + idx = -1; + } else { idx = cljff->addException(boost::get<0>(p), boost::get<1>(p), @@ -2538,6 +2753,12 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, auto pert_idx = idx_to_pert_idx.value(i, openmm_mols.count() + 1); lambda_lever.setExceptionIndicies(pert_idx, "clj", exception_idxs); + if (rb_bond_count > 0) + lambda_lever.setExceptionIndicies(pert_idx, + "ring-break", rb_exception_idxs); + if (rm_bond_count > 0) + lambda_lever.setExceptionIndicies(pert_idx, + "ring-make", rm_exception_idxs); lambda_lever.setConstraintIndicies(pert_idx, constraint_idxs); } From 39289607d19420356c2dcbe0294011d02a9eaba2 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 8 May 2026 14:26:46 +0100 Subject: [PATCH 2/4] Add missing REST2 scaling. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index ab780c13c..fb0152acf 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1753,8 +1753,11 @@ double LambdaLever::setLambda(OpenMM::Context &context, } // Update CLJ exceptions for ring-breaking pairs: scale charge product by - // rb_kappa (0 at bonded end, 1 at nonbonded end). LJ stays at 1e-9 - // so the ring-break CustomBondForce provides the full LJ correction. + // rb_kappa (0 at bonded end, 1 at nonbonded end) and by rest2_scale when + // both atoms are REST2 solute atoms. LJ stays at 1e-9 so the ring-break + // CustomBondForce provides the full LJ correction. + // Note: getBondParameters always returns the original unscaled q (bp[0]) + // because setBondParameters is never called on ring_breaking_ff here. if (cljff != nullptr and ring_breaking_ff != nullptr) { const auto rb_exc = perturbable_mol.getExceptionIndicies("ring-break"); @@ -1767,8 +1770,12 @@ double LambdaLever::setLambda(OpenMM::Context &context, int a0, a1; std::vector bp; ring_breaking_ff->getBondParameters(bond_idx, a0, a1, bp); + double q_scale = rb_kappa; + if (perturbable_mol.isRest2(a0 - start_index) and + perturbable_mol.isRest2(a1 - start_index)) + q_scale *= rest2_scale; cljff->setExceptionParameters(clj_idx, a0, a1, - rb_kappa * bp[0], 1e-9, 1e-9); + q_scale * bp[0], 1e-9, 1e-9); has_changed_cljff = true; } } @@ -1785,8 +1792,12 @@ double LambdaLever::setLambda(OpenMM::Context &context, int a0, a1; std::vector bp; ring_making_ff->getBondParameters(bond_idx, a0, a1, bp); + double q_scale = rm_kappa; + if (perturbable_mol.isRest2(a0 - start_index) and + perturbable_mol.isRest2(a1 - start_index)) + q_scale *= rest2_scale; cljff->setExceptionParameters(clj_idx, a0, a1, - rm_kappa * bp[0], 1e-9, 1e-9); + q_scale * bp[0], 1e-9, 1e-9); has_changed_cljff = true; } } From e0d5942323dc384089584be18d3cba6f19d4d9f3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 12 May 2026 16:34:20 +0100 Subject: [PATCH 3/4] Fix ring-break/make CLJ exception to use morphed charges from CLJ. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 31 +++++++++++----------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index fb0152acf..2eb49512f 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1752,12 +1752,11 @@ double LambdaLever::setLambda(OpenMM::Context &context, } } - // Update CLJ exceptions for ring-breaking pairs: scale charge product by - // rb_kappa (0 at bonded end, 1 at nonbonded end) and by rest2_scale when - // both atoms are REST2 solute atoms. LJ stays at 1e-9 so the ring-break - // CustomBondForce provides the full LJ correction. - // Note: getBondParameters always returns the original unscaled q (bp[0]) - // because setBondParameters is never called on ring_breaking_ff here. + // Update CLJ exceptions for ring-breaking pairs: the exception charge is + // rb_kappa * q_a0(λ) * q_a1(λ), where q_a0/q_a1 are read from the CLJ + // particle parameters that were already morphed (and REST2-scaled via + // sqrt_scale) earlier in this function. LJ stays at 1e-9; the ring-break + // CustomBondForce provides the full LJ correction. if (cljff != nullptr and ring_breaking_ff != nullptr) { const auto rb_exc = perturbable_mol.getExceptionIndicies("ring-break"); @@ -1770,12 +1769,12 @@ double LambdaLever::setLambda(OpenMM::Context &context, int a0, a1; std::vector bp; ring_breaking_ff->getBondParameters(bond_idx, a0, a1, bp); - double q_scale = rb_kappa; - if (perturbable_mol.isRest2(a0 - start_index) and - perturbable_mol.isRest2(a1 - start_index)) - q_scale *= rest2_scale; + double q_a0, sig_a0, eps_a0; + double q_a1, sig_a1, eps_a1; + cljff->getParticleParameters(a0, q_a0, sig_a0, eps_a0); + cljff->getParticleParameters(a1, q_a1, sig_a1, eps_a1); cljff->setExceptionParameters(clj_idx, a0, a1, - q_scale * bp[0], 1e-9, 1e-9); + rb_kappa * q_a0 * q_a1, 1e-9, 1e-9); has_changed_cljff = true; } } @@ -1792,12 +1791,12 @@ double LambdaLever::setLambda(OpenMM::Context &context, int a0, a1; std::vector bp; ring_making_ff->getBondParameters(bond_idx, a0, a1, bp); - double q_scale = rm_kappa; - if (perturbable_mol.isRest2(a0 - start_index) and - perturbable_mol.isRest2(a1 - start_index)) - q_scale *= rest2_scale; + double q_a0, sig_a0, eps_a0; + double q_a1, sig_a1, eps_a1; + cljff->getParticleParameters(a0, q_a0, sig_a0, eps_a0); + cljff->getParticleParameters(a1, q_a1, sig_a1, eps_a1); cljff->setExceptionParameters(clj_idx, a0, a1, - q_scale * bp[0], 1e-9, 1e-9); + rm_kappa * q_a0 * q_a1, 1e-9, 1e-9); has_changed_cljff = true; } } From b3d954afd88cc173bc25ecb66d93558bc42e026b Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 13 May 2026 10:57:37 +0100 Subject: [PATCH 4/4] Exclude ghost atoms from ring-breaking forces to avoid double counting. --- wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 3562fb95c..0bde46e0a 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -2614,11 +2614,17 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, int idx = -1; int nbidx = -1; - const bool is_ring_breaking = rb_pairs_local.contains(IndexPair(atom0, atom1)); - const bool is_ring_making = rm_pairs_local.contains(IndexPair(atom0, atom1)); + bool is_ring_breaking = rb_pairs_local.contains(IndexPair(atom0, atom1)); + bool is_ring_making = rm_pairs_local.contains(IndexPair(atom0, atom1)); if (atom0_is_ghost or atom1_is_ghost) { + // don't add ring-breaking/making forces for pairs involving ghost atoms, + // since the GhostNonbondedForce already provides a softcore interaction + // for these pairs. + is_ring_breaking = false; + is_ring_making = false; + // don't include the LJ term, as this is calculated // elsewhere - note that we need to use 1e-9 to // make sure that OpenMM doesn't eagerly remove