From d10c7775cbb3c0c2f4f1fb0fff30503b2dc332e3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 27 Apr 2026 09:22:45 +0100 Subject: [PATCH 01/21] Add support for caching long-range correction coefficients. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 9 ++++++--- .../SireOpenMM/sire_to_openmm_system.cpp | 19 +++++++++---------- 2 files changed, 15 insertions(+), 13 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index b23f8ab81..b4432b68f 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1925,14 +1925,17 @@ double LambdaLever::setLambda(OpenMM::Context &context, // update the parameters in the context for forces whose parameters changed if (has_changed_cljff) { + if (ghost_ghostff or ghost_nonghostff) + context.setParameter("lambda", std::round(lambda_value * 1e5) / 1e5); + if (cljff) - cljff->updateParametersInContext(context); + cljff->updateParametersInContext(context, true); if (ghost_ghostff) - ghost_ghostff->updateParametersInContext(context); + ghost_ghostff->updateParametersInContext(context, true); if (ghost_nonghostff) - ghost_nonghostff->updateParametersInContext(context); + ghost_nonghostff->updateParametersInContext(context, true); } if (ghost_14ff and has_changed_ghost14ff) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index dea3aa6f0..f4935b843 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -352,14 +352,13 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res } const auto energy_expression = QString( - "e_total;" - "e_total = rho * e_morse + e_repulsion;" - "e_morse = de * (1 - exp(-alpha * delta))^2;" - "e_repulsion = e_rep * (r_sigma / r)^r_pow;" - "alpha = sqrt(k / (2 * de));" - "delta = (r - r0)") - .toStdString(); - + "e_total;" + "e_total = rho * e_morse + e_repulsion;" + "e_morse = de * (1 - exp(-alpha * delta))^2;" + "e_repulsion = e_rep * (r_sigma / r)^r_pow;" + "alpha = sqrt(k / (2 * de));" + "delta = (r - r0)") + .toStdString(); auto *restraintff = new OpenMM::CustomBondForce(energy_expression); restraintff->setName("MorsePotentialRestraintForce"); @@ -1527,15 +1526,15 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, ghost_ghostff->addPerParticleParameter("two_sqrt_epsilon"); ghost_ghostff->addPerParticleParameter("alpha"); ghost_ghostff->addPerParticleParameter("kappa"); + ghost_ghostff->addGlobalParameter("lambda", 0.00000); ghost_nonghostff->addPerParticleParameter("q"); ghost_nonghostff->addPerParticleParameter("half_sigma"); ghost_nonghostff->addPerParticleParameter("two_sqrt_epsilon"); ghost_nonghostff->addPerParticleParameter("alpha"); ghost_nonghostff->addPerParticleParameter("kappa"); + ghost_nonghostff->addGlobalParameter("lambda", 0.00000); - // this will be slow if switched on, as it needs recalculating - // for every change in parameters ghost_ghostff->setUseLongRangeCorrection(use_dispersion_correction); ghost_nonghostff->setUseLongRangeCorrection(use_dispersion_correction); From 535f0d08321fe10a8473be84a8be98ac3b1b1065 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 28 Apr 2026 10:02:18 +0100 Subject: [PATCH 02/21] Split ghost Coulomb and LJ interactions into separate forces. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 57 ++-- wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 2 +- .../SireOpenMM/sire_to_openmm_system.cpp | 243 +++++++++++------- 3 files changed, 187 insertions(+), 115 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index b4432b68f..009a09d15 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1291,8 +1291,10 @@ double LambdaLever::setLambda(OpenMM::Context &context, // get copies of the forcefields in which the parameters will be changed auto qmff = this->getForce("qmff", system); auto cljff = this->getForce("clj", system); - auto ghost_ghostff = this->getForce("ghost/ghost", system); - auto ghost_nonghostff = this->getForce("ghost/non-ghost", system); + auto ghost_ghost_coulombff = this->getForce("ghost/ghost-coulomb", system); + auto ghost_ghost_ljff = this->getForce("ghost/ghost-lj", system); + auto ghost_nonghost_coulombff = this->getForce("ghost/non-ghost-coulomb", system); + auto ghost_nonghost_ljff = this->getForce("ghost/non-ghost-lj", system); auto ghost_14ff = this->getForce("ghost-14", system); auto bondff = this->getForce("bond", system); auto angff = this->getForce("angle", system); @@ -1300,7 +1302,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, auto cmapff = this->getForce("cmap", system); // we know if we have peturbable ghost atoms if we have the ghost forcefields - const bool have_ghost_atoms = (ghost_ghostff != 0 or ghost_nonghostff != 0); + const bool have_ghost_atoms = (ghost_ghost_ljff != 0 or ghost_nonghost_ljff != 0); // whether the constraints have changed bool have_constraints_changed = false; @@ -1320,6 +1322,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, // track whether parameters actually changed for each force, so we only // call updateParametersInContext when necessary bool has_changed_cljff = false; + bool has_changed_coulombff = false; + bool has_changed_ljff = false; bool has_changed_ghost14ff = false; bool has_changed_bondff = false; bool has_changed_angff = false; @@ -1491,7 +1495,11 @@ double LambdaLever::setLambda(OpenMM::Context &context, const int nparams = morphed_charges.count(); // Detect whether any CLJ or ghost-14 parameters changed - has_changed_cljff |= rest2_changed || cache.hasChanged("clj", "charge") || cache.hasChanged("clj", "sigma") || cache.hasChanged("clj", "epsilon") || cache.hasChanged("clj", "alpha") || cache.hasChanged("clj", "kappa") || cache.hasChanged("clj", "charge_scale") || cache.hasChanged("clj", "lj_scale") || cache.hasChanged("ghost/ghost", "charge") || cache.hasChanged("ghost/ghost", "sigma") || cache.hasChanged("ghost/ghost", "epsilon") || cache.hasChanged("ghost/ghost", "alpha") || cache.hasChanged("ghost/ghost", "kappa") || cache.hasChanged("ghost/non-ghost", "charge") || cache.hasChanged("ghost/non-ghost", "sigma") || cache.hasChanged("ghost/non-ghost", "epsilon") || cache.hasChanged("ghost/non-ghost", "alpha") || cache.hasChanged("ghost/non-ghost", "kappa"); + has_changed_cljff |= rest2_changed || cache.hasChanged("clj", "charge") || cache.hasChanged("clj", "sigma") || cache.hasChanged("clj", "epsilon") || cache.hasChanged("clj", "alpha") || cache.hasChanged("clj", "kappa") || cache.hasChanged("clj", "charge_scale") || cache.hasChanged("clj", "lj_scale"); + + has_changed_coulombff |= rest2_changed || cache.hasChanged("ghost/ghost", "charge") || cache.hasChanged("ghost/ghost", "alpha") || cache.hasChanged("ghost/ghost", "kappa") || cache.hasChanged("ghost/non-ghost", "charge") || cache.hasChanged("ghost/non-ghost", "alpha") || cache.hasChanged("ghost/non-ghost", "kappa"); + + has_changed_ljff |= rest2_changed || cache.hasChanged("ghost/ghost", "sigma") || cache.hasChanged("ghost/ghost", "epsilon") || cache.hasChanged("ghost/ghost", "alpha") || cache.hasChanged("ghost/non-ghost", "sigma") || cache.hasChanged("ghost/non-ghost", "epsilon") || cache.hasChanged("ghost/non-ghost", "alpha"); has_changed_ghost14ff |= rest2_changed || cache.hasChanged("ghost-14", "charge") || cache.hasChanged("ghost-14", "sigma") || cache.hasChanged("ghost-14", "epsilon") || cache.hasChanged("ghost-14", "alpha") || cache.hasChanged("ghost-14", "kappa") || cache.hasChanged("ghost-14", "charge_scale") || cache.hasChanged("ghost-14", "lj_scale"); @@ -1535,7 +1543,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, else if (custom_params[4] > 1) custom_params[4] = 1; - ghost_ghostff->setParticleParameters(start_index + j, custom_params); + ghost_ghost_coulombff->setParticleParameters(start_index + j, custom_params); + ghost_ghost_ljff->setParticleParameters(start_index + j, custom_params); // reduced charge custom_params[0] = sqrt_scale * morphed_nonghost_charges[j]; @@ -1560,7 +1569,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, else if (custom_params[4] > 1) custom_params[4] = 1; - ghost_nonghostff->setParticleParameters(start_index + j, custom_params); + ghost_nonghost_coulombff->setParticleParameters(start_index + j, custom_params); + ghost_nonghost_ljff->setParticleParameters(start_index + j, custom_params); if (is_from_ghost or is_to_ghost) { @@ -1923,19 +1933,30 @@ double LambdaLever::setLambda(OpenMM::Context &context, } // update the parameters in the context for forces whose parameters changed - if (has_changed_cljff) + if (has_changed_cljff and cljff) + cljff->updateParametersInContext(context, true); + + if (has_changed_coulombff or has_changed_ljff) { - if (ghost_ghostff or ghost_nonghostff) + // Update the lambda cache key before any ghost force updateParametersInContext. + if (ghost_ghost_ljff or ghost_nonghost_ljff) context.setParameter("lambda", std::round(lambda_value * 1e5) / 1e5); - if (cljff) - cljff->updateParametersInContext(context, true); - - if (ghost_ghostff) - ghost_ghostff->updateParametersInContext(context, true); + if (has_changed_coulombff) + { + if (ghost_ghost_coulombff) + ghost_ghost_coulombff->updateParametersInContext(context, true); + if (ghost_nonghost_coulombff) + ghost_nonghost_coulombff->updateParametersInContext(context, true); + } - if (ghost_nonghostff) - ghost_nonghostff->updateParametersInContext(context, true); + if (has_changed_ljff) + { + if (ghost_ghost_ljff) + ghost_ghost_ljff->updateParametersInContext(context, true); + if (ghost_nonghost_ljff) + ghost_nonghost_ljff->updateParametersInContext(context, true); + } } if (ghost_14ff and has_changed_ghost14ff) @@ -1991,8 +2012,10 @@ double LambdaLever::setLambda(OpenMM::Context &context, // record which named forces had parameters changed in this call last_changed_forces["clj"] = has_changed_cljff; - last_changed_forces["ghost/ghost"] = has_changed_cljff; - last_changed_forces["ghost/non-ghost"] = has_changed_cljff; + last_changed_forces["ghost/ghost-coulomb"] = has_changed_coulombff; + last_changed_forces["ghost/ghost-lj"] = has_changed_ljff; + last_changed_forces["ghost/non-ghost-coulomb"] = has_changed_coulombff; + last_changed_forces["ghost/non-ghost-lj"] = has_changed_ljff; last_changed_forces["ghost-14"] = has_changed_ghost14ff; last_changed_forces["bond"] = has_changed_bondff; last_changed_forces["angle"] = has_changed_angff; diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index e0c0e8133..0ab059933 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -1490,7 +1490,7 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) CODELOC); this->alphas = QVector(cljs.count(), 0.0); - this->kappas = QVector(cljs.count(), 0.0); + this->kappas = QVector(cljs.count(), 1.0); this->perturbed->alphas = this->alphas; this->perturbed->kappas = this->kappas; diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index f4935b843..06e0b7b84 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -475,8 +475,10 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, const double internal_to_k = (1 * SireUnits::kcal_per_mol / (SireUnits::angstrom2)).to(SireUnits::kJ_per_mol / (SireUnits::nanometer2)); auto cljff = lambda_lever.getForce("clj", system); - auto ghost_ghostff = lambda_lever.getForce("ghost/ghost", system); - auto ghost_nonghostff = lambda_lever.getForce("ghost/non-ghost", system); + auto ghost_ghost_coulombff = lambda_lever.getForce("ghost/ghost-coulomb", system); + auto ghost_ghost_ljff = lambda_lever.getForce("ghost/ghost-lj", system); + auto ghost_nonghost_coulombff = lambda_lever.getForce("ghost/non-ghost-coulomb", system); + auto ghost_nonghost_ljff = lambda_lever.getForce("ghost/non-ghost-lj", system); std::vector custom_params = {1.0, 0.0, 0.0}; // Define null parameters used to add these particles to the ghost forces (5 total) @@ -515,14 +517,16 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, cljff->addParticle(0, 0, 0); } - if (ghost_ghostff != 0) + if (ghost_ghost_coulombff != 0) { - ghost_ghostff->addParticle(custom_clj_params); + ghost_ghost_coulombff->addParticle(custom_clj_params); + ghost_ghost_ljff->addParticle(custom_clj_params); } - if (ghost_nonghostff != 0) + if (ghost_nonghost_coulombff != 0) { - ghost_nonghostff->addParticle(custom_clj_params); + ghost_nonghost_coulombff->addParticle(custom_clj_params); + ghost_nonghost_ljff->addParticle(custom_clj_params); } } @@ -538,20 +542,22 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, cljff->addException(anchor_index, atom_index, 0, 0, 0, true); } - if (ghost_ghostff != 0) + if (ghost_ghost_coulombff != 0) { // make sure to exclude interactions between // the atom being positionally restrained and // the anchor - ghost_ghostff->addExclusion(anchor_index, atom_index); + ghost_ghost_coulombff->addExclusion(anchor_index, atom_index); + ghost_ghost_ljff->addExclusion(anchor_index, atom_index); } - if (ghost_nonghostff != 0) + if (ghost_nonghost_coulombff != 0) { // make sure to exclude interactions between // the atom being positionally restrained and // the anchor - ghost_nonghostff->addExclusion(anchor_index, atom_index); + ghost_nonghost_coulombff->addExclusion(anchor_index, atom_index); + ghost_nonghost_ljff->addExclusion(anchor_index, atom_index); } restraintff->addBond(anchor_index, atom_index, custom_params); @@ -1296,8 +1302,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, /// OpenMM::CustomBondForce *ghost_14ff = 0; - OpenMM::CustomNonbondedForce *ghost_ghostff = 0; - OpenMM::CustomNonbondedForce *ghost_nonghostff = 0; + OpenMM::CustomNonbondedForce *ghost_ghost_coulombff = 0; + OpenMM::CustomNonbondedForce *ghost_ghost_ljff = 0; + OpenMM::CustomNonbondedForce *ghost_nonghost_coulombff = 0; + OpenMM::CustomNonbondedForce *ghost_nonghost_ljff = 0; if (any_perturbable) { @@ -1397,7 +1405,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, }; // see below for the description of this energy expression - std::string nb14_expression, clj_expression; + std::string nb14_expression, coul_expression, lj_expression; if (use_taylor_softening) { @@ -1463,18 +1471,21 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // 138.9354558466661 is the constant needed to get energies in // kJ mol-1 given the units of charge (|e|) and distance (nm) // - clj_expression = QString("coul_nrg+lj_nrg;" - "coul_nrg=138.9354558466661*q1*q2*(((%1)/sqrt((%2*max_alpha)+r_safe^2))-(max_kappa/r_safe));" - "lj_nrg=two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" - "sig6=(sigma^6)/(%3*sigma^6 + r_safe^6);" - "r_safe=max(r, 0.001);" - "max_kappa=max(kappa1, kappa2);" - "max_alpha=max(alpha1, alpha2);" - "sigma=half_sigma1+half_sigma2;") - .arg(coulomb_power_expression("max_alpha", coulomb_power)) - .arg(shift_coulomb) - .arg(taylor_power_expression("max_alpha", taylor_power)) - .toStdString(); + coul_expression = QString("138.9354558466661*q1*q2*(((%1)/sqrt((%2*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + "r_safe=max(r, 0.001);" + "max_kappa=max(kappa1, kappa2);" + "max_alpha=max(alpha1, alpha2);") + .arg(coulomb_power_expression("max_alpha", coulomb_power)) + .arg(shift_coulomb) + .toStdString(); + + lj_expression = QString("two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" + "sig6=(sigma^6)/(%1*sigma^6 + r_safe^6);" + "r_safe=max(r, 0.001);" + "max_alpha=max(alpha1, alpha2);" + "sigma=half_sigma1+half_sigma2;") + .arg(taylor_power_expression("max_alpha", taylor_power)) + .toStdString(); } else { @@ -1501,72 +1512,94 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // 138.9354558466661 is the constant needed to get energies in // kJ mol-1 given the units of charge (|e|) and distance (nm) // - clj_expression = QString("coul_nrg+lj_nrg;" - "coul_nrg=138.9354558466661*q1*q2*(((%1)/sqrt((%2*max_alpha)+r_safe^2))-(max_kappa/r_safe));" - "lj_nrg=two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" - "sig6=(sigma^6)/(((sigma*delta) + r_safe^2)^3);" - "delta=%3*max_alpha;" - "r_safe=max(r, 0.001);" - "max_kappa=max(kappa1, kappa2);" - "max_alpha=max(alpha1, alpha2);" - "sigma=half_sigma1+half_sigma2;") - .arg(coulomb_power_expression("max_alpha", coulomb_power)) - .arg(shift_coulomb) - .arg(shift_delta.to(SireUnits::nanometer)) - .toStdString(); - } - - ghost_ghostff = new OpenMM::CustomNonbondedForce(clj_expression); - ghost_ghostff->setName("GhostGhostNonbondedForce"); - ghost_nonghostff = new OpenMM::CustomNonbondedForce(clj_expression); - ghost_nonghostff->setName("GhostNonGhostNonbondedForce"); + coul_expression = QString("138.9354558466661*q1*q2*(((%1)/sqrt((%2*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + "r_safe=max(r, 0.001);" + "max_kappa=max(kappa1, kappa2);" + "max_alpha=max(alpha1, alpha2);") + .arg(coulomb_power_expression("max_alpha", coulomb_power)) + .arg(shift_coulomb) + .toStdString(); - ghost_ghostff->addPerParticleParameter("q"); - ghost_ghostff->addPerParticleParameter("half_sigma"); - ghost_ghostff->addPerParticleParameter("two_sqrt_epsilon"); - ghost_ghostff->addPerParticleParameter("alpha"); - ghost_ghostff->addPerParticleParameter("kappa"); - ghost_ghostff->addGlobalParameter("lambda", 0.00000); + lj_expression = QString("two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" + "sig6=(sigma^6)/(((sigma*delta) + r_safe^2)^3);" + "delta=%1*max_alpha;" + "r_safe=max(r, 0.001);" + "max_alpha=max(alpha1, alpha2);" + "sigma=half_sigma1+half_sigma2;") + .arg(shift_delta.to(SireUnits::nanometer)) + .toStdString(); + } - ghost_nonghostff->addPerParticleParameter("q"); - ghost_nonghostff->addPerParticleParameter("half_sigma"); - ghost_nonghostff->addPerParticleParameter("two_sqrt_epsilon"); - ghost_nonghostff->addPerParticleParameter("alpha"); - ghost_nonghostff->addPerParticleParameter("kappa"); - ghost_nonghostff->addGlobalParameter("lambda", 0.00000); + ghost_ghost_coulombff = new OpenMM::CustomNonbondedForce(coul_expression); + ghost_ghost_coulombff->setName("GhostGhostCoulombForce"); + ghost_ghost_ljff = new OpenMM::CustomNonbondedForce(lj_expression); + ghost_ghost_ljff->setName("GhostGhostLJForce"); + ghost_nonghost_coulombff = new OpenMM::CustomNonbondedForce(coul_expression); + ghost_nonghost_coulombff->setName("GhostNonGhostCoulombForce"); + ghost_nonghost_ljff = new OpenMM::CustomNonbondedForce(lj_expression); + ghost_nonghost_ljff->setName("GhostNonGhostLJForce"); + + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + { + ff->addPerParticleParameter("q"); + ff->addPerParticleParameter("half_sigma"); + ff->addPerParticleParameter("two_sqrt_epsilon"); + ff->addPerParticleParameter("alpha"); + ff->addPerParticleParameter("kappa"); + ff->addGlobalParameter("lambda", 0.00000); + } - ghost_ghostff->setUseLongRangeCorrection(use_dispersion_correction); - ghost_nonghostff->setUseLongRangeCorrection(use_dispersion_correction); + // Coulomb forces must not use LRC: the soft-core 1/sqrt(alpha+r^2)-1/r + // expression decays as 1/r^3 when alpha>0, making the LRC integral + // (proportional to int r^2 * 1/r^3 dr = int 1/r dr) log-divergent. + ghost_ghost_coulombff->setUseLongRangeCorrection(false); + ghost_nonghost_coulombff->setUseLongRangeCorrection(false); + // LJ soft-core decays as 1/r^6, so LRC is well-defined. + ghost_ghost_ljff->setUseLongRangeCorrection(use_dispersion_correction); + ghost_nonghost_ljff->setUseLongRangeCorrection(use_dispersion_correction); if (ffinfo.hasCutoff()) { if (ffinfo.space().isPeriodic()) { - ghost_ghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffPeriodic); - ghost_nonghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffPeriodic); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffPeriodic); } else { - ghost_ghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffNonPeriodic); - ghost_nonghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffNonPeriodic); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffNonPeriodic); } - ghost_ghostff->setCutoffDistance(ffinfo.cutoff().to(SireUnits::nanometers)); - ghost_nonghostff->setCutoffDistance(ffinfo.cutoff().to(SireUnits::nanometers)); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->setCutoffDistance(ffinfo.cutoff().to(SireUnits::nanometers)); } else { - ghost_ghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::NoCutoff); - ghost_nonghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::NoCutoff); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->setNonbondedMethod(OpenMM::CustomNonbondedForce::NoCutoff); } - ghost_ghostff->setForceGroup(force_group_counter); - lambda_lever.setForceIndex("ghost/ghost", system.addForce(ghost_ghostff)); - lambda_lever.setForceGroup("ghost/ghost", force_group_counter++); + ghost_ghost_coulombff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ghost/ghost-coulomb", system.addForce(ghost_ghost_coulombff)); + lambda_lever.setForceGroup("ghost/ghost-coulomb", force_group_counter++); + + ghost_ghost_ljff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ghost/ghost-lj", system.addForce(ghost_ghost_ljff)); + lambda_lever.setForceGroup("ghost/ghost-lj", force_group_counter++); + + ghost_nonghost_coulombff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ghost/non-ghost-coulomb", system.addForce(ghost_nonghost_coulombff)); + lambda_lever.setForceGroup("ghost/non-ghost-coulomb", force_group_counter++); - ghost_nonghostff->setForceGroup(force_group_counter); - lambda_lever.setForceIndex("ghost/non-ghost", system.addForce(ghost_nonghostff)); - lambda_lever.setForceGroup("ghost/non-ghost", force_group_counter++); + ghost_nonghost_ljff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ghost/non-ghost-lj", system.addForce(ghost_nonghost_ljff)); + lambda_lever.setForceGroup("ghost/non-ghost-lj", force_group_counter++); ghost_14ff->setForceGroup(force_group_counter); lambda_lever.setForceIndex("ghost-14", system.addForce(ghost_14ff)); @@ -1687,11 +1720,13 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // add a perturbable molecule, recording the start index // for each of the forcefields - start_indicies.reserve(7); + start_indicies.reserve(9); start_indicies.insert("clj", start_index); - start_indicies.insert("ghost/ghost", start_index); - start_indicies.insert("ghost/non-ghost", start_index); + start_indicies.insert("ghost/ghost-coulomb", start_index); + start_indicies.insert("ghost/ghost-lj", start_index); + start_indicies.insert("ghost/non-ghost-coulomb", start_index); + start_indicies.insert("ghost/non-ghost-lj", start_index); // the start index for this molecules first bond, angle or // torsion parameters will be however many of these @@ -1781,8 +1816,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, custom_params[4] = kappas_data[j]; // Add the particle to the ghost and nonghost forcefields - ghost_ghostff->addParticle(custom_params); - ghost_nonghostff->addParticle(custom_params); + ghost_ghost_coulombff->addParticle(custom_params); + ghost_ghost_ljff->addParticle(custom_params); + ghost_nonghost_coulombff->addParticle(custom_params); + ghost_nonghost_ljff->addParticle(custom_params); real_atoms.append(atom_index); @@ -1861,8 +1898,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, custom_params[3] = 0.0; // kappa - is 0 for non-ghost atoms custom_params[4] = 0.0; - ghost_ghostff->addParticle(custom_params); - ghost_nonghostff->addParticle(custom_params); + ghost_ghost_coulombff->addParticle(custom_params); + ghost_ghost_ljff->addParticle(custom_params); + ghost_nonghost_coulombff->addParticle(custom_params); + ghost_nonghost_ljff->addParticle(custom_params); non_ghost_atoms.append(atom_index); } } @@ -1933,8 +1972,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // kappa custom_params[4] = kappas_data[mol.molinfo.nAtoms() + k]; - ghost_ghostff->addParticle(custom_params); - ghost_nonghostff->addParticle(custom_params); + ghost_ghost_coulombff->addParticle(custom_params); + ghost_ghost_ljff->addParticle(custom_params); + ghost_nonghost_coulombff->addParticle(custom_params); + ghost_nonghost_ljff->addParticle(custom_params); const bool vs_to_ghost = mol.to_ghost_idxs.contains(parent_idx); const bool vs_from_ghost = mol.from_ghost_idxs.contains(parent_idx); @@ -1956,8 +1997,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, { // Add to ghost FFs if necessary custom_params = {vs_charge, 1.0, 0.0, 0.0, 0.0}; - ghost_ghostff->addParticle(custom_params); - ghost_nonghostff->addParticle(custom_params); + ghost_ghost_coulombff->addParticle(custom_params); + ghost_ghost_ljff->addParticle(custom_params); + ghost_nonghost_coulombff->addParticle(custom_params); + ghost_nonghost_ljff->addParticle(custom_params); non_ghost_atoms.append(atom_index); } } @@ -2133,14 +2176,16 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, /// Finally tell the ghost forcefields about the ghost and non-ghost /// interaction groups, so that they can correctly calculate the /// ghost/ghost and ghost/non-ghost energies - if (ghost_ghostff != 0 and ghost_nonghostff != 0) + if (ghost_ghost_ljff != 0 and ghost_nonghost_ljff != 0) { // set up the interaction groups - ghost / non-ghost // ghost / ghost std::set ghost_atoms_set(ghost_atoms.begin(), ghost_atoms.end()); std::set non_ghost_atoms_set(non_ghost_atoms.begin(), non_ghost_atoms.end()); - ghost_ghostff->addInteractionGroup(ghost_atoms_set, ghost_atoms_set); - ghost_nonghostff->addInteractionGroup(ghost_atoms_set, non_ghost_atoms_set); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff}) + ff->addInteractionGroup(ghost_atoms_set, ghost_atoms_set); + for (auto ff : {ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->addInteractionGroup(ghost_atoms_set, non_ghost_atoms_set); } // see if we want to remove COM motion @@ -2320,10 +2365,11 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // we need to make sure that the list of exclusions in // the NonbondedForce match those in the CustomNonbondedForces - if (ghost_ghostff != 0) + if (ghost_ghost_ljff != 0) { - ghost_ghostff->addExclusion(boost::get<0>(p), boost::get<1>(p)); - ghost_nonghostff->addExclusion(boost::get<0>(p), boost::get<1>(p)); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->addExclusion(boost::get<0>(p), boost::get<1>(p)); } } @@ -2351,10 +2397,11 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, cljff->addException(vs0_index, start_index + a, 0.0, 1, 0, false); - if (ghost_ghostff != 0) + if (ghost_ghost_ljff != 0) { - ghost_ghostff->addExclusion(vs0_index, start_index + a); - ghost_nonghostff->addExclusion(vs0_index, start_index + a); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->addExclusion(vs0_index, start_index + a); } for (int v1 = v0 + 1; v1 < atom_vs.size(); ++v1) @@ -2363,10 +2410,11 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, cljff->addException(vs0_index, vs1_index, 0.0, 1, 0, false); - if (ghost_ghostff != 0) + if (ghost_ghost_ljff != 0) { - ghost_ghostff->addExclusion(vs0_index, vs1_index); - ghost_nonghostff->addExclusion(vs0_index, vs1_index); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->addExclusion(vs0_index, vs1_index); } } } @@ -2396,8 +2444,9 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // and if the two atoms are in the same molecule if (mol_from == mol_to) { - ghost_ghostff->addExclusion(from_ghost_idx, to_ghost_idx); - ghost_nonghostff->addExclusion(from_ghost_idx, to_ghost_idx); + for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, + ghost_nonghost_coulombff, ghost_nonghost_ljff}) + ff->addExclusion(from_ghost_idx, to_ghost_idx); cljff->addException(from_ghost_idx, to_ghost_idx, 0.0, 1e-9, 1e-9, false); } From 02d8ed4f430f090f01ef1244273bb5deac268118 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 28 Apr 2026 14:14:23 +0100 Subject: [PATCH 03/21] Only need to use preserveLongRangeCorrection with CustomNonbondedForce. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 009a09d15..19c2e9096 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1934,7 +1934,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, // update the parameters in the context for forces whose parameters changed if (has_changed_cljff and cljff) - cljff->updateParametersInContext(context, true); + cljff->updateParametersInContext(context); if (has_changed_coulombff or has_changed_ljff) { From 20257f8bd76e81c935d7a3e10fbe42f704ec8db2 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 28 Apr 2026 14:52:49 +0100 Subject: [PATCH 04/21] Revert debugging update. --- wrapper/Convert/SireOpenMM/openmmmolecule.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 0ab059933..e0c0e8133 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -1490,7 +1490,7 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) CODELOC); this->alphas = QVector(cljs.count(), 0.0); - this->kappas = QVector(cljs.count(), 1.0); + this->kappas = QVector(cljs.count(), 0.0); this->perturbed->alphas = this->alphas; this->perturbed->kappas = this->kappas; From b293f3aef294e50e82f118ddebef6e0f50283671 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 28 Apr 2026 15:30:11 +0100 Subject: [PATCH 05/21] Add internal caching of CustomNonbondedForce LRC coefficients. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 81 +++++++++++++++++++ wrapper/Convert/SireOpenMM/lambdalever.h | 11 +++ .../SireOpenMM/sire_to_openmm_system.cpp | 20 ++++- 3 files changed, 109 insertions(+), 3 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 19c2e9096..2e4a41c64 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1959,6 +1959,86 @@ double LambdaLever::setLambda(OpenMM::Context &context, } } + // Update the ghost LJ dispersion correction via the CustomVolumeForce. + // At r > rc the soft-core shift is negligible, so the standard closed-form + // LJ tail integral applies. Results are cached per lambda state. + auto ghost_lrc_ff = this->getForce("ghost-lrc", system); + if (ghost_lrc_ff != nullptr && ghost_ghost_ljff != nullptr && ghost_nonghost_ljff != nullptr) + { + const qint64 lam_key = qRound64(lambda_value * 1e5); + double lrc_coeff = 0.0; + + if (this->lrc_coeff_cache.contains(lam_key)) + { + lrc_coeff = this->lrc_coeff_cache[lam_key]; + } + else + { + const double cutoff = ghost_ghost_ljff->getCutoffDistance(); + const double rc3 = cutoff * cutoff * cutoff; + const double rc9 = rc3 * rc3 * rc3; + const double four_pi = 4.0 * M_PI; + + // Interaction group sets: ghost atoms and non-ghost atoms. + std::set ghost_set, dummy_set, nonghost_set; + ghost_ghost_ljff->getInteractionGroupParameters(0, ghost_set, dummy_set); + ghost_nonghost_ljff->getInteractionGroupParameters(0, dummy_set, nonghost_set); + + // Cache half_sigma and two_sqrt_epsilon for each ghost atom. + QHash> ghost_params; + for (int i : ghost_set) + { + std::vector p; + ghost_ghost_ljff->getParticleParameters(i, p); + ghost_params[i] = {p[1], p[2]}; // half_sigma, two_sqrt_epsilon + } + + // Cache non-ghost params. + QHash> nonghost_params; + for (int j : nonghost_set) + { + std::vector p; + ghost_nonghost_ljff->getParticleParameters(j, p); + nonghost_params[j] = {p[1], p[2]}; + } + + // ghost-ghost unique pairs (i < j). + for (auto it_i = ghost_set.cbegin(); it_i != ghost_set.cend(); ++it_i) + { + const auto &pi = ghost_params[*it_i]; + auto it_j = it_i; + for (++it_j; it_j != ghost_set.cend(); ++it_j) + { + const auto &pj = ghost_params[*it_j]; + const double sig = pi.first + pj.first; + const double sig2 = sig * sig; + const double sig6 = sig2 * sig2 * sig2; + const double eps_pair = pi.second * pj.second; + lrc_coeff += four_pi * eps_pair * sig6 * (sig6 / (9.0 * rc9) - 1.0 / (3.0 * rc3)); + } + } + + // ghost-nonghost all pairs. + for (auto it_i = ghost_set.cbegin(); it_i != ghost_set.cend(); ++it_i) + { + const auto &pi = ghost_params[*it_i]; + for (int j : nonghost_set) + { + const auto &pj = nonghost_params[j]; + const double sig = pi.first + pj.first; + const double sig2 = sig * sig; + const double sig6 = sig2 * sig2 * sig2; + const double eps_pair = pi.second * pj.second; + lrc_coeff += four_pi * eps_pair * sig6 * (sig6 / (9.0 * rc9) - 1.0 / (3.0 * rc3)); + } + } + + this->lrc_coeff_cache[lam_key] = lrc_coeff; + } + + context.setParameter("lrc_coeff", lrc_coeff); + } + if (ghost_14ff and has_changed_ghost14ff) ghost_14ff->updateParametersInContext(context); @@ -2016,6 +2096,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, last_changed_forces["ghost/ghost-lj"] = has_changed_ljff; last_changed_forces["ghost/non-ghost-coulomb"] = has_changed_coulombff; last_changed_forces["ghost/non-ghost-lj"] = has_changed_ljff; + last_changed_forces["ghost-lrc"] = has_changed_ljff; last_changed_forces["ghost-14"] = has_changed_ghost14ff; last_changed_forces["bond"] = has_changed_bondff; last_changed_forces["angle"] = has_changed_angff; diff --git a/wrapper/Convert/SireOpenMM/lambdalever.h b/wrapper/Convert/SireOpenMM/lambdalever.h index 3167c01bb..f9a9de6d3 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.h +++ b/wrapper/Convert/SireOpenMM/lambdalever.h @@ -214,6 +214,11 @@ namespace SireOpenMM * so we can detect when it actually changes. Initialised to -1 as a * sentinel meaning "never been set". */ mutable double last_qmff_lam; + + /** Cache of pre-computed ghost LJ dispersion coefficients keyed by + * rounded lambda (qRound64(lambda * 1e5)). Populated on first visit + * to each lambda state and reused on warm passes. */ + mutable QHash lrc_coeff_cache; }; #ifndef SIRE_SKIP_INLINE_FUNCTION @@ -242,6 +247,12 @@ namespace SireOpenMM return "OpenMM::CustomCVForce"; } + template <> + inline QString _get_typename() + { + return "OpenMM::CustomVolumeForce"; + } + template <> inline QString _get_typename() { diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 06e0b7b84..eb9d2f7e9 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1555,9 +1555,12 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // (proportional to int r^2 * 1/r^3 dr = int 1/r dr) log-divergent. ghost_ghost_coulombff->setUseLongRangeCorrection(false); ghost_nonghost_coulombff->setUseLongRangeCorrection(false); - // LJ soft-core decays as 1/r^6, so LRC is well-defined. - ghost_ghost_ljff->setUseLongRangeCorrection(use_dispersion_correction); - ghost_nonghost_ljff->setUseLongRangeCorrection(use_dispersion_correction); + // LJ soft-core LRC is handled analytically via a CustomVolumeForce + // rather than OpenMM's numerical integrator, because the standard + // LJ tail (r > rc, soft-core shift negligible) has a closed-form + // expression and this allows the result to be cached per lambda state. + ghost_ghost_ljff->setUseLongRangeCorrection(false); + ghost_nonghost_ljff->setUseLongRangeCorrection(false); if (ffinfo.hasCutoff()) { @@ -1601,6 +1604,17 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, lambda_lever.setForceIndex("ghost/non-ghost-lj", system.addForce(ghost_nonghost_ljff)); lambda_lever.setForceGroup("ghost/non-ghost-lj", force_group_counter++); + // Analytic LJ LRC: E = lrc_coeff / V, updated each lambda step via + // a cached closed-form sum over interaction-group pairs. + if (use_dispersion_correction && ffinfo.hasCutoff() && ffinfo.space().isPeriodic()) + { + auto ghost_lrc_ff = new OpenMM::CustomVolumeForce("lrc_coeff/v"); + ghost_lrc_ff->addGlobalParameter("lrc_coeff", 0.0); + ghost_lrc_ff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ghost-lrc", system.addForce(ghost_lrc_ff)); + lambda_lever.setForceGroup("ghost-lrc", force_group_counter++); + } + ghost_14ff->setForceGroup(force_group_counter); lambda_lever.setForceIndex("ghost-14", system.addForce(ghost_14ff)); lambda_lever.setForceGroup("ghost-14", force_group_counter++); From 4dc6bc9dc7571d31cde178f7e85df4ffd9465142 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 28 Apr 2026 15:40:18 +0100 Subject: [PATCH 06/21] No longer need preserveLongRangeCorrection parameter. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 2e4a41c64..012df8a61 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1945,17 +1945,17 @@ double LambdaLever::setLambda(OpenMM::Context &context, if (has_changed_coulombff) { if (ghost_ghost_coulombff) - ghost_ghost_coulombff->updateParametersInContext(context, true); + ghost_ghost_coulombff->updateParametersInContext(context); if (ghost_nonghost_coulombff) - ghost_nonghost_coulombff->updateParametersInContext(context, true); + ghost_nonghost_coulombff->updateParametersInContext(context); } if (has_changed_ljff) { if (ghost_ghost_ljff) - ghost_ghost_ljff->updateParametersInContext(context, true); + ghost_ghost_ljff->updateParametersInContext(context); if (ghost_nonghost_ljff) - ghost_nonghost_ljff->updateParametersInContext(context, true); + ghost_nonghost_ljff->updateParametersInContext(context); } } From 356901a1ca05315e064f49bb248160d23970584b Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Apr 2026 14:17:14 +0100 Subject: [PATCH 07/21] Update XML parser to handle new custom forces. --- src/sire/morph/_xml.py | 57 ++++++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 30 deletions(-) diff --git a/src/sire/morph/_xml.py b/src/sire/morph/_xml.py index 19c3c864e..777d2482a 100644 --- a/src/sire/morph/_xml.py +++ b/src/sire/morph/_xml.py @@ -20,7 +20,8 @@ def evaluate_xml_force(mols, xml, force): force : str The name of the custom force to evaluate. Options are: - "ghost-ghost", "ghost-nonghost", "ghost-14". + "ghost-ghost-lj", "ghost-ghost-coulomb", + "ghost-nonghost-lj", "ghost-nonghost-coulomb", "ghost-14". Returns ------- @@ -35,8 +36,6 @@ def evaluate_xml_force(mols, xml, force): The Lennard-Jones energies for each atom pair. """ - from math import sqrt - import xml.etree.ElementTree as ET import sys @@ -78,21 +77,28 @@ def evaluate_xml_force(mols, xml, force): ) # Validate the force name. - if not force in ["ghostghost", "ghostnonghost", "ghost14"]: + valid = [ + "ghostghostlj", + "ghostghostcoulomb", + "ghostnonghostlj", + "ghostnonghostcoulomb", + "ghost14", + ] + if force not in valid: raise ValueError( - "'force' must be one of 'ghost-ghost', 'ghost-nonghost', or 'ghost-14'." + "'force' must be one of 'ghost-ghost-lj', 'ghost-ghost-coulomb', " + "'ghost-nonghost-lj', 'ghost-nonghost-coulomb', or 'ghost-14'." ) - # Create the name and index based on the force type. - if force == "ghostghost": - name = "GhostGhostNonbondedForce" - elif force == "ghostnonghost": - name = "GhostNonGhostNonbondedForce" - elif force == "ghost14": - name = "Ghost14BondForce" - - # Get the root of the XML tree. - root = tree.getroot() + # Map sanitised name to the OpenMM force name in the XML. + _force_name_map = { + "ghostghostlj": "GhostGhostLJForce", + "ghostghostcoulomb": "GhostGhostCoulombForce", + "ghostnonghostlj": "GhostNonGhostLJForce", + "ghostnonghostcoulomb": "GhostNonGhostCoulombForce", + "ghost14": "Ghost14BondForce", + } + name = _force_name_map[force] # Loop over the forces until we find the named CustomNonbondedForce. is_found = False @@ -163,11 +169,8 @@ def evaluate_xml_force(mols, xml, force): atom_i = atoms[i] # Set the parameters for this particle. - setattr(module, parameters[0] + "1", float(particle_i.get("param1"))) - setattr(module, parameters[1] + "1", float(particle_i.get("param2"))) - setattr(module, parameters[2] + "1", float(particle_i.get("param3"))) - setattr(module, parameters[3] + "1", float(particle_i.get("param4"))) - setattr(module, parameters[4] + "1", float(particle_i.get("param5"))) + for k, param in enumerate(parameters): + setattr(module, param + "1", float(particle_i.get(f"param{k + 1}"))) # Loop over particles in set2. for y in range(len(set2)): @@ -190,11 +193,8 @@ def evaluate_xml_force(mols, xml, force): atom_j = atoms[j] # Set the parameters for this particle. - setattr(module, parameters[0] + "2", float(particle_j.get("param1"))) - setattr(module, parameters[1] + "2", float(particle_j.get("param2"))) - setattr(module, parameters[2] + "2", float(particle_j.get("param3"))) - setattr(module, parameters[3] + "2", float(particle_j.get("param4"))) - setattr(module, parameters[4] + "2", float(particle_j.get("param5"))) + for k, param in enumerate(parameters): + setattr(module, param + "2", float(particle_j.get(f"param{k + 1}"))) # Get the distance between the particles. r = measure(atom_i, atom_j).to(nanometer) @@ -239,11 +239,8 @@ def evaluate_xml_force(mols, xml, force): atom_j = atoms[int(bond.get("p2"))] # Set the parameters for this bond. - setattr(module, parameters[0], float(bond.get("param1"))) - setattr(module, parameters[1], float(bond.get("param2"))) - setattr(module, parameters[2], float(bond.get("param3"))) - setattr(module, parameters[3], float(bond.get("param4"))) - setattr(module, parameters[4], float(bond.get("param5"))) + for k, param in enumerate(parameters): + setattr(module, param, float(bond.get(f"param{k + 1}"))) # Get the distance between the particles. r = measure(atom_i, atom_j).to(nanometer) From b27c346ff7358cd6693df4793d6fe03620e39fc6 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Apr 2026 14:23:56 +0100 Subject: [PATCH 08/21] Add full LRC handling via CustomVolume forces. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 73 ++++++++ wrapper/Convert/SireOpenMM/lambdalever.h | 12 ++ .../SireOpenMM/sire_to_openmm_system.cpp | 166 +++++++++++++++++- 3 files changed, 248 insertions(+), 3 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 012df8a61..a92ba62b8 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -438,6 +438,11 @@ bool LambdaLever::wasForceChanged(const QString &name) const return it.value(); } +void LambdaLever::setGCMCWaterAtoms(const QVector &atoms) +{ + gcmc_water_atoms = QSet(atoms.begin(), atoms.end()); +} + boost::tuple get_exception(int atom0, int atom1, int start_index, double coul_14_scl, double lj_14_scl, @@ -2039,6 +2044,72 @@ double LambdaLever::setLambda(OpenMM::Context &context, context.setParameter("lrc_coeff", lrc_coeff); } + // 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); + if (background_lrc_ff != nullptr && cljff != nullptr) + { + const qint64 lam_key = qRound64(lambda_value * 1e5); + double lrc_coeff = 0.0; + + if (this->background_lrc_coeff_cache.contains(lam_key)) + { + lrc_coeff = this->background_lrc_coeff_cache[lam_key]; + } + else + { + const double cutoff = cljff->getCutoffDistance(); + const double rc3 = cutoff * cutoff * cutoff; + const double rc9 = rc3 * rc3 * rc3; + const double four_pi = 4.0 * M_PI; + + // Classify particles by (sigma, epsilon); ghost atoms (epsilon=0), + // virtual sites, and GCMC water atoms are skipped. + std::map, int> class_counts; + for (int i = 0; i < cljff->getNumParticles(); ++i) + { + double charge, sigma, epsilon; + cljff->getParticleParameters(i, charge, sigma, epsilon); + if (epsilon == 0.0) + continue; + if (!gcmc_water_atoms.isEmpty() && gcmc_water_atoms.contains(i)) + continue; + class_counts[{sigma, epsilon}]++; + } + + // Diagonal pairs (same class, unique i(n) * (n - 1) * 0.5; + if (n_pairs == 0.0) + continue; + const double sig2 = key.first * key.first; + const double sig6 = sig2 * sig2 * sig2; + const double eps_pair = 4.0 * key.second; + lrc_coeff += n_pairs * four_pi * eps_pair * sig6 * (sig6 / (9.0 * rc9) - 1.0 / (3.0 * rc3)); + } + + // Off-diagonal pairs (class i != class j). + for (auto it1 = class_counts.cbegin(); it1 != class_counts.cend(); ++it1) + { + auto it2 = it1; + for (++it2; it2 != class_counts.cend(); ++it2) + { + const double sigma_ij = 0.5 * (it1->first.first + it2->first.first); + const double eps_pair = 4.0 * std::sqrt(it1->first.second * it2->first.second); + const double n_pairs = static_cast(it1->second) * it2->second; + const double sig2 = sigma_ij * sigma_ij; + const double sig6 = sig2 * sig2 * sig2; + lrc_coeff += n_pairs * four_pi * eps_pair * sig6 * (sig6 / (9.0 * rc9) - 1.0 / (3.0 * rc3)); + } + } + + this->background_lrc_coeff_cache[lam_key] = lrc_coeff; + } + + context.setParameter("lrc_background", lrc_coeff); + } + if (ghost_14ff and has_changed_ghost14ff) ghost_14ff->updateParametersInContext(context); @@ -2097,6 +2168,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, last_changed_forces["ghost/non-ghost-coulomb"] = has_changed_coulombff; last_changed_forces["ghost/non-ghost-lj"] = has_changed_ljff; last_changed_forces["ghost-lrc"] = has_changed_ljff; + 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["bond"] = has_changed_bondff; last_changed_forces["angle"] = has_changed_angff; diff --git a/wrapper/Convert/SireOpenMM/lambdalever.h b/wrapper/Convert/SireOpenMM/lambdalever.h index f9a9de6d3..0e2cbe494 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.h +++ b/wrapper/Convert/SireOpenMM/lambdalever.h @@ -165,6 +165,8 @@ namespace SireOpenMM QStringList getForceNames() const; bool wasForceChanged(const QString &name) const; + void setGCMCWaterAtoms(const QVector &atoms); + protected: void updateRestraintInContext(OpenMM::Force &ff, double rho, OpenMM::Context &context) const; @@ -219,6 +221,16 @@ namespace SireOpenMM * rounded lambda (qRound64(lambda * 1e5)). Populated on first visit * to each lambda state and reused on warm passes. */ mutable QHash lrc_coeff_cache; + + /** Cache of pre-computed background (NonbondedForce) LJ dispersion + * coefficients keyed by rounded lambda. Mirrors lrc_coeff_cache but + * covers all non-ghost atoms in the clj force. */ + mutable QHash background_lrc_coeff_cache; + + /** OpenMM atom indices that belong to GCMC water molecules. When + * non-empty these atoms are excluded from background-lrc (their LRC + * is handled by the gcmc-lrc CustomVolumeForce instead). */ + QSet gcmc_water_atoms; }; #ifndef SIRE_SKIP_INLINE_FUNCTION diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index eb9d2f7e9..01445a07e 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1184,9 +1184,21 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, use_dispersion_correction = map["use_dispersion_correction"].value().asABoolean(); } - // note that this will be very slow for perturbable systems, as - // it needs recalculating for every change of lambda - cljff->setUseDispersionCorrection(use_dispersion_correction); + bool is_gcmc = false; + int num_gcmc_waters = 0; + + if (map.specified("use_gcmc_lrc")) + { + is_gcmc = map["use_gcmc_lrc"].value().asABoolean(); + } + if (is_gcmc && map.specified("num_gcmc_waters")) + { + num_gcmc_waters = map["num_gcmc_waters"].value().asAnInteger(); + } + + // LRC for the NonbondedForce is handled analytically via a CustomVolumeForce + // (background-lrc) updated each lambda step, so we always disable it here. + cljff->setUseDispersionCorrection(false); // set the non-bonded cutoff type and length based on // the infomation in ffinfo @@ -1620,6 +1632,32 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, lambda_lever.setForceGroup("ghost-14", force_group_counter++); } + // Analytic LRC for the NonbondedForce (all non-ghost atoms): E = lrc_background / V, + // updated each lambda step via a cached closed-form class-pair sum. + if (use_dispersion_correction && ffinfo.hasCutoff() && ffinfo.space().isPeriodic()) + { + auto background_lrc_ff = new OpenMM::CustomVolumeForce("lrc_background/v"); + background_lrc_ff->addGlobalParameter("lrc_background", 0.0); + background_lrc_ff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("background-lrc", system.addForce(background_lrc_ff)); + lambda_lever.setForceGroup("background-lrc", force_group_counter++); + } + + // GCMC water LRC: E = (n_w * lrc_w_solute + n_w*(n_w-1) * lrc_ww_half) / V. + // lrc_w_solute and lrc_ww_half are pre-computed at setup; n_w is updated by + // the GCMC sampler at each insertion/deletion move. + if (is_gcmc && use_dispersion_correction && ffinfo.hasCutoff() && ffinfo.space().isPeriodic()) + { + auto gcmc_lrc_ff = new OpenMM::CustomVolumeForce( + "(n_w * lrc_w_solute + n_w * (n_w - 1) * lrc_ww_half) / v"); + gcmc_lrc_ff->addGlobalParameter("n_w", 0.0); + gcmc_lrc_ff->addGlobalParameter("lrc_w_solute", 0.0); + gcmc_lrc_ff->addGlobalParameter("lrc_ww_half", 0.0); + gcmc_lrc_ff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("gcmc-lrc", system.addForce(gcmc_lrc_ff)); + lambda_lever.setForceGroup("gcmc-lrc", force_group_counter++); + } + // Stage 4 is complete. We now have all(*) of the forces we need to run // a perturbable simulation. (*) well, we will define the restraint // forces in a much later stage after the particles have been added. @@ -2202,6 +2240,128 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, ff->addInteractionGroup(ghost_atoms_set, non_ghost_atoms_set); } + // Register GCMC water atom indices with the lambda lever and pre-compute the + // fixed LRC coefficients (lrc_w_solute and lrc_ww_half) for the gcmc-lrc force. + if (is_gcmc && use_dispersion_correction && ffinfo.hasCutoff() && ffinfo.space().isPeriodic()) + { + auto gcmc_lrc_ff = lambda_lever.getForce("gcmc-lrc", system); + if (gcmc_lrc_ff != nullptr) + { + const double cutoff = cljff->getCutoffDistance(); + const double rc3 = cutoff * cutoff * cutoff; + const double rc9 = rc3 * rc3 * rc3; + const double four_pi = 4.0 * M_PI; + + // All GCMC waters (real + virtual buffer) are excluded from the background + // LRC and tracked by n_w instead. Any water can be swapped in or out. + const auto water_result = mols.search("water"); + QSet water_mol_nums; + for (const auto &view : water_result.views()) + water_mol_nums.insert(view.data().number()); + + // Collect OpenMM atom indices for all water molecules. + QVector water_atom_indices; + for (int i = 0; i < nmols; ++i) + { + if (!water_mol_nums.contains(mols[i].number())) + continue; + const int mol_start = start_indexes[i]; + const int mol_natoms = openmm_mols_data[i].masses.count(); + for (int j = mol_start; j < mol_start + mol_natoms; ++j) + water_atom_indices.append(j); + } + + lambda_lever.setGCMCWaterAtoms(water_atom_indices); + + // Pre-compute lrc_w_solute (per active water molecule, interaction with + // all non-water atoms: protein, ligand, ions) and lrc_ww_half (per active + // water-molecule pair, halved). FEP ghost atoms have epsilon=0 and + // contribute zero automatically. + // n_w starts at n_all_waters - num_gcmc_waters (initially active count). + QSet water_set(water_atom_indices.begin(), water_atom_indices.end()); + + // Collect (sigma, epsilon) for water atom types and non-water solute atoms. + // Real waters are also in the GCMC pool so solute = protein + ligand + ions. + std::map, int> water_class_counts; + std::map, int> solute_class_counts; + for (int i = 0; i < cljff->getNumParticles(); ++i) + { + double charge, sigma, epsilon; + cljff->getParticleParameters(i, charge, sigma, epsilon); + if (epsilon == 0.0) + continue; + if (water_set.contains(i)) + water_class_counts[{sigma, epsilon}]++; + else + solute_class_counts[{sigma, epsilon}]++; + } + + // lrc_ww_half: half the LRC coefficient for one water-molecule pair. + // Sum over all water-atom-type pairs within a molecule pair. + // Each molecule pair contributes once (factor 1/2 already in the name). + const int n_water_mols = water_mol_nums.size(); + double lrc_ww = 0.0; + // diagonal water-water class pairs (same type within a water molecule pair) + for (const auto &[key, n] : water_class_counts) + { + // n atoms of this type spread across n_water_mols molecules: + // per-mol count = n / n_water_mols + const double per_mol = static_cast(n) / n_water_mols; + const double n_pairs = per_mol * per_mol; + const double sig2 = key.first * key.first; + const double sig6 = sig2 * sig2 * sig2; + const double eps_pair = 4.0 * key.second; + lrc_ww += n_pairs * four_pi * eps_pair * sig6 * (sig6 / (9.0 * rc9) - 1.0 / (3.0 * rc3)); + } + // off-diagonal water-water class pairs + for (auto it1 = water_class_counts.cbegin(); it1 != water_class_counts.cend(); ++it1) + { + auto it2 = it1; + for (++it2; it2 != water_class_counts.cend(); ++it2) + { + const double per_mol_1 = static_cast(it1->second) / n_water_mols; + const double per_mol_2 = static_cast(it2->second) / n_water_mols; + const double sigma_ij = 0.5 * (it1->first.first + it2->first.first); + const double eps_pair = 4.0 * std::sqrt(it1->first.second * it2->first.second); + const double sig2 = sigma_ij * sigma_ij; + const double sig6 = sig2 * sig2 * sig2; + lrc_ww += 2.0 * per_mol_1 * per_mol_2 * four_pi * eps_pair * sig6 * (sig6 / (9.0 * rc9) - 1.0 / (3.0 * rc3)); + } + } + const double lrc_ww_half = 0.5 * lrc_ww; + + // lrc_w_solute: LRC coefficient for one water molecule with all solute atoms. + double lrc_w_solute = 0.0; + for (const auto &[wkey, wn] : water_class_counts) + { + const double per_mol_w = static_cast(wn) / n_water_mols; + for (const auto &[skey, sn] : solute_class_counts) + { + const double sigma_ij = 0.5 * (wkey.first + skey.first); + const double eps_pair = 4.0 * std::sqrt(wkey.second * skey.second); + const double sig2 = sigma_ij * sigma_ij; + const double sig6 = sig2 * sig2 * sig2; + lrc_w_solute += per_mol_w * sn * four_pi * eps_pair * sig6 * (sig6 / (9.0 * rc9) - 1.0 / (3.0 * rc3)); + } + } + + // Update the CustomVolumeForce default parameter values so they are + // correct when the OpenMM Context is created. + // n_w starts at the number of initially active waters (total - buffer). + const double n_w_initial = static_cast(n_water_mols - num_gcmc_waters); + for (int p = 0; p < gcmc_lrc_ff->getNumGlobalParameters(); ++p) + { + const auto &name = gcmc_lrc_ff->getGlobalParameterName(p); + if (name == "n_w") + gcmc_lrc_ff->setGlobalParameterDefaultValue(p, n_w_initial); + else if (name == "lrc_w_solute") + gcmc_lrc_ff->setGlobalParameterDefaultValue(p, lrc_w_solute); + else if (name == "lrc_ww_half") + gcmc_lrc_ff->setGlobalParameterDefaultValue(p, lrc_ww_half); + } + } + } + // see if we want to remove COM motion const auto com_remove_prop = map["com_reset_frequency"]; From f5e4edf14e9d491e28f5f4827e8357810f67d197 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Apr 2026 14:24:51 +0100 Subject: [PATCH 09/21] Add unit test for LRC. --- tests/convert/test_openmm_lrc.py | 102 +++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 tests/convert/test_openmm_lrc.py diff --git a/tests/convert/test_openmm_lrc.py b/tests/convert/test_openmm_lrc.py new file mode 100644 index 000000000..ffb7f23b2 --- /dev/null +++ b/tests/convert/test_openmm_lrc.py @@ -0,0 +1,102 @@ +""" +Validate that the CustomVolumeForce-based LRC in the perturbable OpenMM system +gives the same total energy as a non-perturbable end-state system that uses the +standard NonbondedForce dispersion correction. + +Tests use merged_ethane_methanol (merged_molecule.s3) with the 5 nearest waters +so they run quickly while still exercising a periodic cutoff system. +""" + +import pytest +import sire as sr + + +def _get_end_state(mol, state, remove_state): + c = mol.cursor() + for key in c.keys(): + if key.endswith(state): + c[key.removesuffix(state)] = c[key] + del c[key] + elif key.endswith(remove_state): + del c[key] + c["is_perturbable"] = False + return c.commit() + + +def _build_systems(mols, platform): + space = mols.space() + + c = mols.cursor() + c["molidx 0"]["coordinates"] = c["molidx 0"]["coordinates0"] + c["molidx 0"]["coordinates1"] = c["molidx 0"]["coordinates0"] + mols = c.commit() + + merge = mols[0] + water = mols["closest 5 waters to molidx 0"] + + mols_pert = merge + water + mols_ref = _get_end_state(merge, "0", "1") + water + mols_pert_end = _get_end_state(merge, "1", "0") + water + + l = sr.cas.LambdaSchedule() + l.add_stage("morph", (1 - l.lam()) * l.initial() + l.lam() * l.final()) + + map = { + "platform": platform, + "schedule": l, + "constraint": "h-bonds-not-perturbed", + "include_constrained_energies": True, + "dynamic_constraints": False, + "use_dispersion_correction": True, + "space": space, + } + + omm_pert = sr.convert.to(mols_pert, "openmm", map=map) + omm_ref = sr.convert.to(mols_ref, "openmm", map=map) + omm_pert_end = sr.convert.to(mols_pert_end, "openmm", map=map) + + return omm_pert, omm_ref, omm_pert_end + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_lrc_lambda0_matches_reference_end_state( + merged_ethane_methanol, openmm_platform +): + """ + Perturbable system at lambda=0 must give the same energy as a freshly built + non-perturbable reference (lambda=0) end state with standard NonbondedForce LRC. + """ + omm_pert, omm_ref, _ = _build_systems( + merged_ethane_methanol.clone(), openmm_platform + ) + + omm_pert.set_lambda(0.0) + nrg_pert = omm_pert.get_energy().value() + nrg_ref = omm_ref.get_energy().value() + + assert nrg_pert == pytest.approx(nrg_ref, rel=1e-3) + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_lrc_lambda1_matches_perturbed_end_state( + merged_ethane_methanol, openmm_platform +): + """ + Perturbable system at lambda=1 must give the same energy as a freshly built + non-perturbable perturbed (lambda=1) end state with standard NonbondedForce LRC. + """ + omm_pert, _, omm_pert_end = _build_systems( + merged_ethane_methanol.clone(), openmm_platform + ) + + omm_pert.set_lambda(1.0) + nrg_pert = omm_pert.get_energy().value() + nrg_pert_end = omm_pert_end.get_energy().value() + + assert nrg_pert == pytest.approx(nrg_pert_end, rel=1e-3) From 758101f15ddf5c17cea51291e3faa3babef5b8b7 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Apr 2026 14:25:23 +0100 Subject: [PATCH 10/21] Add paragraph on LRC handling. --- doc/source/tutorial/part07/03_ghosts.rst | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/doc/source/tutorial/part07/03_ghosts.rst b/doc/source/tutorial/part07/03_ghosts.rst index 2fafe665e..70ba66b9c 100644 --- a/doc/source/tutorial/part07/03_ghosts.rst +++ b/doc/source/tutorial/part07/03_ghosts.rst @@ -55,6 +55,19 @@ Ghost atoms are then added to three custom OpenMM Forces: interaction and subtracts this from the total (to remove the real-space contribution that was calculated in the standard OpenMM NonBondedForce). +When ``use_dispersion_correction=True`` and the system uses a periodic +cutoff, two additional ``CustomVolumeForce`` objects are added: + +* A **background LRC** force (``lrc_background/v``) that evaluates the + LJ long-range correction for all non-ghost atoms analytically. This + replaces the built-in dispersion correction of the ``NonbondedForce``, + which would otherwise be recomputed on every λ change. The coefficient + is cached per λ state by the lambda lever. + +* A **ghost LRC** force (``lrc_coeff/v``) that evaluates the LJ + long-range correction for all ghost–ghost and ghost–non-ghost + interactions analytically. Its coefficient is also cached per λ state. + There are two different soft-core potentials available. The default is the Zacharias potential, while the second is the Taylor potential. From 4975dd28052825bca9e2d0a28f81a7a7c86d0340 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Apr 2026 15:55:55 +0100 Subject: [PATCH 11/21] Remove redundant parameters from CustomNonbondedForce objects. --- tests/convert/test_openmm_rest2.py | 54 ++++----- wrapper/Convert/SireOpenMM/lambdalever.cpp | 89 ++++++++------- .../SireOpenMM/sire_to_openmm_system.cpp | 106 ++++++++++-------- 3 files changed, 141 insertions(+), 108 deletions(-) diff --git a/tests/convert/test_openmm_rest2.py b/tests/convert/test_openmm_rest2.py index 7ac498a52..af4338044 100644 --- a/tests/convert/test_openmm_rest2.py +++ b/tests/convert/test_openmm_rest2.py @@ -91,32 +91,34 @@ def test_rest2(mols, rest2_selection, excluded_atoms, request): if is_perturbable: # Find the ghost/ghost nonbonded interaction. for force in omm_system.getForces(): - if force.getName() == "GhostGhostNonbondedForce": - break + if force.getName() == "GhostGhostCoulombForce": + ghost_ghost_coulomb_force = force + elif force.getName() == "GhostGhostLJForce": + ghost_ghost_lj_force = force # Store the initial parameters. ghost_ghost_params_initial = [] excluded_ghost_ghost_indices = [] - for i in range(force.getNumParticles()): - charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( - force.getParticleParameters(i) - ) + for i in range(ghost_ghost_coulomb_force.getNumParticles()): + charge, _, _ = ghost_ghost_coulomb_force.getParticleParameters(i) + _, two_sqrt_epsilon, _ = ghost_ghost_lj_force.getParticleParameters(i) ghost_ghost_params_initial.append((charge, two_sqrt_epsilon)) if i in excluded_atoms: excluded_ghost_ghost_indices.append(i) # Find the ghost/non-ghost nonbonded interaction. for force in omm_system.getForces(): - if force.getName() == "GhostNonGhostNonbondedForce": - break + if force.getName() == "GhostNonGhostCoulombForce": + ghost_non_ghost_coulomb_force = force + elif force.getName() == "GhostNonGhostLJForce": + ghost_non_ghost_lj_force = force # Store the initial parameters. ghost_non_ghost_params_initial = [] excluded_ghost_non_ghost_indices = [] - for i in range(force.getNumParticles()): - charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( - force.getParticleParameters(i) - ) + for i in range(ghost_non_ghost_coulomb_force.getNumParticles()): + charge, _, _ = ghost_non_ghost_coulomb_force.getParticleParameters(i) + _, two_sqrt_epsilon, _ = ghost_non_ghost_lj_force.getParticleParameters(i) ghost_non_ghost_params_initial.append((charge, two_sqrt_epsilon)) if i in excluded_atoms: excluded_ghost_non_ghost_indices.append(i) @@ -151,32 +153,34 @@ def test_rest2(mols, rest2_selection, excluded_atoms, request): for i in range(force.getNumExceptions()): exception_params_modified.append(force.getExceptionParameters(i)[-3::2]) - # Find the ghost/ghost nonbonded interaction. + # Find the ghost/ghost nonbonded forces. for force in omm_system.getForces(): - if force.getName() == "GhostGhostNonbondedForce": - break + if force.getName() == "GhostGhostCoulombForce": + ghost_ghost_coulomb_force = force + elif force.getName() == "GhostGhostLJForce": + ghost_ghost_lj_force = force # Handle custom forces for pertubable molecules. if is_perturbable: # Store the modified parameters. ghost_ghost_params_modified = [] - for i in range(force.getNumParticles()): - charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( - force.getParticleParameters(i) - ) + for i in range(ghost_ghost_coulomb_force.getNumParticles()): + charge, _, _ = ghost_ghost_coulomb_force.getParticleParameters(i) + _, two_sqrt_epsilon, _ = ghost_ghost_lj_force.getParticleParameters(i) ghost_ghost_params_modified.append((charge, two_sqrt_epsilon)) # Find the ghost/non-ghost nonbonded interaction. for force in omm_system.getForces(): - if force.getName() == "GhostNonGhostNonbondedForce": - break + if force.getName() == "GhostNonGhostCoulombForce": + ghost_non_ghost_coulomb_force = force + elif force.getName() == "GhostNonGhostLJForce": + ghost_non_ghost_lj_force = force # Store the modified parameters. ghost_non_ghost_params_modified = [] - for i in range(force.getNumParticles()): - charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( - force.getParticleParameters(i) - ) + for i in range(ghost_non_ghost_coulomb_force.getNumParticles()): + charge, _, _ = ghost_non_ghost_coulomb_force.getParticleParameters(i) + _, two_sqrt_epsilon, _ = ghost_non_ghost_lj_force.getParticleParameters(i) ghost_non_ghost_params_modified.append((charge, two_sqrt_epsilon)) # Store the scaling factor. diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index a92ba62b8..fd10875b2 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1312,9 +1312,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, // whether the constraints have changed bool have_constraints_changed = false; - // whether any CMAP map parameters were set (tracked to defer updateParametersInContext) - - std::vector custom_params = {0.0, 0.0, 0.0, 0.0, 0.0}; + std::vector custom_params_coul = {0.0, 0.0, 0.0}; + std::vector custom_params_lj = {0.0, 0.0, 0.0}; if (qmff != 0) { @@ -1526,56 +1525,68 @@ double LambdaLever::setLambda(OpenMM::Context &context, } // reduced charge - custom_params[0] = sqrt_scale * morphed_ghost_charges[j]; + custom_params_coul[0] = sqrt_scale * morphed_ghost_charges[j]; + // alpha + custom_params_coul[1] = morphed_ghost_alphas[j]; + // kappa + custom_params_coul[2] = morphed_ghost_kappas[j]; + // half_sigma - custom_params[1] = 0.5 * morphed_ghost_sigmas[j]; + custom_params_lj[0] = 0.5 * morphed_ghost_sigmas[j]; // two_sqrt_epsilon - custom_params[2] = 2.0 * sqrt_scale * std::sqrt(morphed_ghost_epsilons[j]); + custom_params_lj[1] = 2.0 * sqrt_scale * std::sqrt(morphed_ghost_epsilons[j]); // alpha - custom_params[3] = morphed_ghost_alphas[j]; - // kappa - custom_params[4] = morphed_ghost_kappas[j]; + custom_params_lj[2] = morphed_ghost_alphas[j]; // clamp alpha between 0 and 1 - if (custom_params[3] < 0) - custom_params[3] = 0; - else if (custom_params[3] > 1) - custom_params[3] = 1; + if (custom_params_coul[1] < 0) + { + custom_params_coul[1] = 0; + custom_params_lj[2] = 0; + } + else if (custom_params_coul[1] > 1) + { + custom_params_coul[1] = 1; + custom_params_lj[2] = 1; + } // clamp kappa between 0 and 1 - if (custom_params[4] < 0) - custom_params[4] = 0; - else if (custom_params[4] > 1) - custom_params[4] = 1; + if (custom_params_coul[2] < 0) + custom_params_coul[2] = 0; + else if (custom_params_coul[2] > 1) + custom_params_coul[2] = 1; - ghost_ghost_coulombff->setParticleParameters(start_index + j, custom_params); - ghost_ghost_ljff->setParticleParameters(start_index + j, custom_params); + ghost_ghost_coulombff->setParticleParameters(start_index + j, custom_params_coul); + ghost_ghost_ljff->setParticleParameters(start_index + j, custom_params_lj); // reduced charge - custom_params[0] = sqrt_scale * morphed_nonghost_charges[j]; + custom_params_coul[0] = sqrt_scale * morphed_nonghost_charges[j]; + // alpha + custom_params_coul[1] = morphed_nonghost_alphas[j]; + // kappa + custom_params_coul[2] = morphed_nonghost_kappas[j]; + // half_sigma - custom_params[1] = 0.5 * morphed_nonghost_sigmas[j]; + custom_params_lj[0] = 0.5 * morphed_nonghost_sigmas[j]; // two_sqrt_epsilon - custom_params[2] = 2.0 * sqrt_scale * std::sqrt(morphed_nonghost_epsilons[j]); + custom_params_lj[1] = 2.0 * sqrt_scale * std::sqrt(morphed_nonghost_epsilons[j]); // alpha - custom_params[3] = morphed_nonghost_alphas[j]; - // kappa - custom_params[4] = morphed_nonghost_kappas[j]; + custom_params_lj[2] = morphed_nonghost_alphas[j]; // clamp alpha between 0 and 1 - if (custom_params[3] < 0) - custom_params[3] = 0; - else if (custom_params[3] > 1) - custom_params[3] = 1; - - // clamp kappa between 0 and 1 - if (custom_params[4] < 0) - custom_params[4] = 0; - else if (custom_params[4] > 1) - custom_params[4] = 1; + if (custom_params_coul[1] < 0) + { + custom_params_coul[1] = 0; + custom_params_lj[2] = 0; + } + else if (custom_params_coul[1] > 1) + { + custom_params_coul[1] = 1; + custom_params_lj[2] = 1; + } - ghost_nonghost_coulombff->setParticleParameters(start_index + j, custom_params); - ghost_nonghost_ljff->setParticleParameters(start_index + j, custom_params); + ghost_nonghost_coulombff->setParticleParameters(start_index + j, custom_params_coul); + ghost_nonghost_ljff->setParticleParameters(start_index + j, custom_params_lj); if (is_from_ghost or is_to_ghost) { @@ -1995,7 +2006,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, { std::vector p; ghost_ghost_ljff->getParticleParameters(i, p); - ghost_params[i] = {p[1], p[2]}; // half_sigma, two_sqrt_epsilon + ghost_params[i] = {p[0], p[1]}; // half_sigma, two_sqrt_epsilon } // Cache non-ghost params. @@ -2004,7 +2015,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, { std::vector p; ghost_nonghost_ljff->getParticleParameters(j, p); - nonghost_params[j] = {p[1], p[2]}; + nonghost_params[j] = {p[0], p[1]}; } // ghost-ghost unique pairs (i < j). diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 01445a07e..6147f05da 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -482,7 +482,7 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, std::vector custom_params = {1.0, 0.0, 0.0}; // Define null parameters used to add these particles to the ghost forces (5 total) - std::vector custom_clj_params = {0.0, 0.0, 0.0, 0.0, 0.0}; + std::vector custom_clj_params = {0.0, 0.0, 0.0}; // we need to add all of the positions as anchor particles for (const auto &restraint : atom_restraints) @@ -1551,14 +1551,19 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, ghost_nonghost_ljff = new OpenMM::CustomNonbondedForce(lj_expression); ghost_nonghost_ljff->setName("GhostNonGhostLJForce"); - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) + for (auto ff : {ghost_ghost_coulombff, ghost_nonghost_coulombff}) { ff->addPerParticleParameter("q"); + ff->addPerParticleParameter("alpha"); + ff->addPerParticleParameter("kappa"); + ff->addGlobalParameter("lambda", 0.00000); + } + + for (auto ff : {ghost_ghost_ljff, ghost_nonghost_ljff}) + { ff->addPerParticleParameter("half_sigma"); ff->addPerParticleParameter("two_sqrt_epsilon"); ff->addPerParticleParameter("alpha"); - ff->addPerParticleParameter("kappa"); ff->addGlobalParameter("lambda", 0.00000); } @@ -1689,9 +1694,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, QHash idx_to_pert_idx; QHash idx_to_qm_idx; - // just a holder for all of the custom parameters for the - // ghost forces (prevents us having to continually re-allocate it) - std::vector custom_params = {0.0, 0.0, 0.0, 0.0, 0.0}; + // holders for all of custom parameters for the ghost forces (prevents us + // having to continually re-allocate it) + std::vector custom_params_coul = {0.0, 0.0, 0.0}; + std::vector custom_params_lj = {0.0, 0.0, 0.0}; // the sets of particle indexes for the ghost atoms and non-ghost atoms QVector ghost_atoms; @@ -1857,21 +1863,24 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, } // reduced_q - custom_params[0] = charge; + custom_params_coul[0] = charge; + // alpha + custom_params_coul[1] = alphas_data[j]; + // kappa + custom_params_coul[2] = kappas_data[j]; + // half_sigma - custom_params[1] = 0.5 * boost::get<1>(clj); + custom_params_lj[0] = 0.5 * boost::get<1>(clj); // two_sqrt_epsilon - custom_params[2] = 2.0 * std::sqrt(boost::get<2>(clj)); + custom_params_lj[1] = 2.0 * std::sqrt(boost::get<2>(clj)); // alpha - custom_params[3] = alphas_data[j]; - // kappa - custom_params[4] = kappas_data[j]; + custom_params_lj[2] = alphas_data[j]; // Add the particle to the ghost and nonghost forcefields - ghost_ghost_coulombff->addParticle(custom_params); - ghost_ghost_ljff->addParticle(custom_params); - ghost_nonghost_coulombff->addParticle(custom_params); - ghost_nonghost_ljff->addParticle(custom_params); + ghost_ghost_coulombff->addParticle(custom_params_coul); + ghost_ghost_ljff->addParticle(custom_params_lj); + ghost_nonghost_coulombff->addParticle(custom_params_coul); + ghost_nonghost_ljff->addParticle(custom_params_lj); real_atoms.append(atom_index); @@ -1940,20 +1949,24 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // forcefields if there are any perturbable molecules if (any_perturbable) { - // reduced charge - custom_params[0] = boost::get<0>(clj); + // reduced_q + custom_params_coul[0] = boost::get<0>(clj); + // alpha - is zero for non-ghost atoms + custom_params_coul[1] = 0.0; + // kappa - is 0 for non-ghost atoms + custom_params_coul[2] = 0.0; + // half_sigma - custom_params[1] = 0.5 * boost::get<1>(clj); + custom_params_lj[0] = 0.5 * boost::get<1>(clj); // two_sqrt_epsilon - custom_params[2] = 2.0 * std::sqrt(boost::get<2>(clj)); + custom_params_lj[1] = 2.0 * std::sqrt(boost::get<2>(clj)); // alpha - is zero for non-ghost atoms - custom_params[3] = 0.0; - // kappa - is 0 for non-ghost atoms - custom_params[4] = 0.0; - ghost_ghost_coulombff->addParticle(custom_params); - ghost_ghost_ljff->addParticle(custom_params); - ghost_nonghost_coulombff->addParticle(custom_params); - ghost_nonghost_ljff->addParticle(custom_params); + custom_params_lj[2] = 0.0; + + ghost_ghost_coulombff->addParticle(custom_params_coul); + ghost_ghost_ljff->addParticle(custom_params_lj); + ghost_nonghost_coulombff->addParticle(custom_params_coul); + ghost_nonghost_ljff->addParticle(custom_params_lj); non_ghost_atoms.append(atom_index); } } @@ -2007,27 +2020,31 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, system.setVirtualSite(atom_index, new_vs); // Add to forcefield, depending on whether the system is perturbable - // Note that VS with LJ parameters are currently not supported, so epsilon and sigma are hard-coded to 0 in all cases + // Note that VS with LJ parameters are currently not supported, + // so epsilon and sigma are hard-coded to 0 in all cases double vs_charge = mol.vs_charges.at(k).asADouble(); cljff->addParticle(vs_charge, 1.0, 0.0); if (any_perturbable and mol.isPerturbable()) { // reduced_q - custom_params[0] = vs_charge; + custom_params_coul[0] = vs_charge; + // alpha + custom_params_coul[1] = alphas_data[mol.molinfo.nAtoms() + k]; + // kappa + custom_params_coul[2] = kappas_data[mol.molinfo.nAtoms() + k]; + // half_sigma - custom_params[1] = 1.0; + custom_params_lj[0] = 1.0; // two_sqrt_epsilon - custom_params[2] = 0.0; + custom_params_lj[1] = 0.0; // alpha - custom_params[3] = alphas_data[mol.molinfo.nAtoms() + k]; - // kappa - custom_params[4] = kappas_data[mol.molinfo.nAtoms() + k]; + custom_params_lj[2] = alphas_data[mol.molinfo.nAtoms() + k]; - ghost_ghost_coulombff->addParticle(custom_params); - ghost_ghost_ljff->addParticle(custom_params); - ghost_nonghost_coulombff->addParticle(custom_params); - ghost_nonghost_ljff->addParticle(custom_params); + ghost_ghost_coulombff->addParticle(custom_params_coul); + ghost_ghost_ljff->addParticle(custom_params_lj); + ghost_nonghost_coulombff->addParticle(custom_params_coul); + ghost_nonghost_ljff->addParticle(custom_params_lj); const bool vs_to_ghost = mol.to_ghost_idxs.contains(parent_idx); const bool vs_from_ghost = mol.from_ghost_idxs.contains(parent_idx); @@ -2048,11 +2065,12 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, else if (any_perturbable) { // Add to ghost FFs if necessary - custom_params = {vs_charge, 1.0, 0.0, 0.0, 0.0}; - ghost_ghost_coulombff->addParticle(custom_params); - ghost_ghost_ljff->addParticle(custom_params); - ghost_nonghost_coulombff->addParticle(custom_params); - ghost_nonghost_ljff->addParticle(custom_params); + custom_params_coul = {vs_charge, 0.0, 0.0}; + custom_params_lj = {1.0, 0.0, 0.0}; + ghost_ghost_coulombff->addParticle(custom_params_coul); + ghost_ghost_ljff->addParticle(custom_params_lj); + ghost_nonghost_coulombff->addParticle(custom_params_coul); + ghost_nonghost_ljff->addParticle(custom_params_lj); non_ghost_atoms.append(atom_index); } } From 0772799fc98ae2c5f472a28cca460544dbe5e207 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 30 Apr 2026 09:59:41 +0100 Subject: [PATCH 12/21] Add Beutler soft-core form for ABFE. --- doc/source/changelog.rst | 10 ++ doc/source/cheatsheet/openmm.rst | 14 ++- doc/source/tutorial/part07/03_ghosts.rst | 85 +++++++++++-- src/sire/mol/__init__.py | 18 +-- src/sire/mol/_dynamics.py | 2 - src/sire/mol/_minimisation.py | 2 - src/sire/system/_system.py | 12 +- .../SireOpenMM/sire_to_openmm_system.cpp | 117 ++++++++++++------ 8 files changed, 174 insertions(+), 86 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 7b66d5970..f43950e00 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -67,6 +67,16 @@ organisation on `GitHub `__. * Added per-stage weights to :class:`~sire.cas.LambdaSchedule`, allowing stages to occupy unequal fractions of lambda space (e.g. ``add_morph_stage("morph", weight=2.0)``). +* Added support for the Beutler et al. softcore potential (``use_beutler_softening``, + ``beutler_alpha`` map options), currently intended for ABFE calculations. Removed the + ``coulomb_power`` parameter, which was invalid for any non-zero value. + +* Added analytic LJ long-range correction (LRC) for periodic simulations + (``use_dispersion_correction`` map option). A background LRC is computed for all + non-ghost atoms. When ghost atoms are present, a separate ghost LRC covers ghost–ghost + and ghost–non-ghost interactions. Both are updated efficiently via the lambda lever + without recomputing the full dispersion correction on every λ change. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/doc/source/cheatsheet/openmm.rst b/doc/source/cheatsheet/openmm.rst index 9f46bdd04..005a70773 100644 --- a/doc/source/cheatsheet/openmm.rst +++ b/doc/source/cheatsheet/openmm.rst @@ -122,9 +122,9 @@ Available keys and allowable values are listed below. | | ``bonds-h-angles-not-perturbed`` and | | | ``bonds-h-angles-not-heavy-perturbed | +------------------------------+----------------------------------------------------------+ -| coulomb_power | The coulomb power parameter used by the softening | -| | potential used to soften electrostatic interactions | -| | involving ghost atoms. This defaults to 0. | +| beutler_alpha | The alpha parameter for the Beutler softcore LJ | +| | potential. Controls the width of the softening in | +| | r^6-space. This defaults to 0.5. | +------------------------------+----------------------------------------------------------+ | cutoff | Size of the non-bonded cutoff, e.g. | | | ``7.5*sr.units.angstrom`` | @@ -226,14 +226,18 @@ Available keys and allowable values are listed below. | use_dispersion_correction | Whether or not to use the dispersion correction to | | | deal with cutoff issues. This is very expensive. | +------------------------------+----------------------------------------------------------+ +| use_beutler_softening | Whether or not to use the Beutler softcore potential to | +| | soften interactions involving ghost atoms. Currently | +| | designed for ABFE calculations only. This defaults to | +| | False. If True, overrides taylor and zacharias softening.| ++------------------------------+----------------------------------------------------------+ | use_taylor_softening | Whether or not to use the taylor algorithm to soften | | | interactions involving ghost atoms. This defaults to | | | False. | +------------------------------+----------------------------------------------------------+ | use_zacharias_softening | Whether or not to use the zacharias algorithm to soften | | | interactions involving ghost atoms. This defaults to | -| | True. Note that one of zacharias or taylor softening | -| | must be True, with zacharias taking precedence. | +| | True. | +------------------------------+----------------------------------------------------------+ Higher level API diff --git a/doc/source/tutorial/part07/03_ghosts.rst b/doc/source/tutorial/part07/03_ghosts.rst index 70ba66b9c..4ff7f4e4b 100644 --- a/doc/source/tutorial/part07/03_ghosts.rst +++ b/doc/source/tutorial/part07/03_ghosts.rst @@ -68,8 +68,9 @@ cutoff, two additional ``CustomVolumeForce`` objects are added: long-range correction for all ghost–ghost and ghost–non-ghost interactions analytically. Its coefficient is also cached per λ state. -There are two different soft-core potentials available. The default is -the Zacharias potential, while the second is the Taylor potential. +There are three different soft-core potentials available. The default is +the Zacharias potential, while the second is the Taylor potential, and +the third is the Beutler potential. Zacharias softening ------------------- @@ -81,7 +82,7 @@ It is based on the following electrostatic and Lennard-Jones potentials: .. math:: - V_{\text{elec}}(r) = q_i q_j \left[ \frac{(1 - \alpha)^n}{\sqrt{r^2 + \delta_\text{coulomb}}} - \frac{\kappa}{r} \right] + V_{\text{elec}}(r) = q_i q_j \left[ \frac{1}{\sqrt{r^2 + \delta_\text{coulomb}}} - \frac{\kappa}{r} \right] V_{\text{LJ}}(r) = 4\epsilon \left[ \frac{\sigma^{12}}{(\delta_\text{LJ} \sigma + r^2)^6} - \frac{\sigma^6}{(\delta_\text{LJ} \sigma + r^2)^3} \right] @@ -116,10 +117,6 @@ The soft-core parameters are: state, and 1 in the perturbed state. These values can be perturbed via the ``alpha`` lever in the λ-schedule. -* ``n`` is the "coulomb power", and is set to 0 by default. It can be - any integer between 0 and 4. It is set via ``coulomb_power`` map - parameter. - * ``shift_coulomb`` and ``shift_LJ`` are the so-called "shift delta" parameters, which are specified individually for the coulomb and LJ\ potentials. They are set via the ``shift_coulomb`` and ``shift_delta`` @@ -145,7 +142,7 @@ It is based on the following electrostatic and Lennard-Jones potentials: .. math:: - V_{\text{elec}}(r) = q_i q_j \left[ \frac{(1 - \alpha)^n}{\sqrt{r^2 + \delta^2}} - \frac{\kappa}{r} \right] + V_{\text{elec}}(r) = q_i q_j \left[ \frac{1}{\sqrt{r^2 + \delta^2}} - \frac{\kappa}{r} \right] V_{\text{LJ}}(r) = 4\epsilon \left[ \frac{\sigma^{12}}{(\alpha^m \sigma^6 + r^6)^2} - \frac{\sigma^6}{\alpha^m \sigma^6 + r^6} \right] @@ -182,10 +179,6 @@ The soft-core parameters are: any integer between 0 and 4. It is set via ``taylor_power`` map parameter. -* ``n`` is the "coulomb power", and is set to 0 by default. It can be - any integer between 0 and 4. It is set via ``coulomb_power`` map - parameter. - * ``shift_coulomb`` is the so-called "shift delta" parameters, which are specified only for the coulomb potential. This is set via the ``shift_coulomb`` @@ -200,6 +193,74 @@ The soft-core parameters are: intramolecular electrostatic interactions, when the "hard" interaction would not be calculated in the NonbondedForce. +Beutler softening +----------------- + +This is the third soft-core potential, based on Beutler et al., +*Chem. Phys. Lett.*, 1994. You can use it by setting the map option +``use_beutler_softening`` to True. + +.. note:: + + The Beutler potential is currently designed for absolute binding free + energy (ABFE) calculations only. It is not recommended for relative + free energy (RBFE) calculations. + +It is based on the following electrostatic and Lennard-Jones potentials: + +.. math:: + + V_{\text{elec}}(r) = q_i q_j \left[ \frac{1}{\sqrt{r^2 + \delta^2}} - \frac{\kappa}{r} \right] + + V_{\text{LJ}}(r) = (1 - \alpha) \cdot 4\epsilon \left[ \frac{\sigma^{12}}{(\beta \sigma^6 \alpha + r^6)^2} - \frac{\sigma^6}{\beta \sigma^6 \alpha + r^6} \right] + +where + +.. math:: + + \delta = \alpha \times \text{shift_coulomb} + +and + +.. math:: + + \alpha = \max(\alpha_i, \alpha_j) + + \kappa = \max(\kappa_i, \kappa_j) + +The parameters ``r``, ``q_i``, ``q_j``, ``\epsilon``, and ``\sigma`` +are the standard parameters for the electrostatic and Lennard-Jones +potentials. + +The soft-core parameters are: + +* ``α_i`` and ``α_j`` control the amount of "softening" of the + electrostatic and LJ interactions. A value of 0 means no softening + (fully hard), while a value of 1 means fully soft. Ghost atoms which + disappear as a function of λ have a value of α of 1 in the + reference state, and 0 in the perturbed state. Ghost atoms which appear + as a function of λ have a value of α of 0 in the reference + state, and 1 in the perturbed state. These values can be perturbed + via the ``alpha`` lever in the λ-schedule. + +* ``β`` is the Beutler alpha parameter that controls the softening of the + LJ interaction. It is set via the ``beutler_alpha`` map parameter and + defaults to 0.5. + +* ``shift_coulomb`` is the "shift delta" parameter for the electrostatic + potential. It is set via the ``shift_coulomb`` map parameter and + defaults to 1 Å. + +* ``κ_i`` and ``κ_j`` are the "hard" electrostatic parameters, + which control whether or not to calculate the "hard" electrostatic + interaction to subtract from the total energy and force (thus cancelling + out the double-counting of this interaction from the NonbondedForce). + By default, these are always equal to 1. You can perturb these via the + ``kappa`` lever in the λ-schedule, e.g. if you want to decouple the + intramolecular electrostatic interactions, when the "hard" interaction + would not be calculated in the NonbondedForce. + + Good practice ------------- diff --git a/src/sire/mol/__init__.py b/src/sire/mol/__init__.py index 20a0d5ca8..77409ac90 100644 --- a/src/sire/mol/__init__.py +++ b/src/sire/mol/__init__.py @@ -971,7 +971,7 @@ def __fixed__dihedral__(obj, idx=None, idx1=None, idx2=None, idx3=None, map=None raise KeyError("There is no matching dihedral in this view.") elif len(dihedrals) > 1: raise KeyError( - "More than one dihedral matches. Number of " f"matches is {len(dihedrals)}." + f"More than one dihedral matches. Number of matches is {len(dihedrals)}." ) return dihedrals[0] @@ -986,7 +986,7 @@ def __fixed__improper__(obj, idx=None, idx1=None, idx2=None, idx3=None, map=None raise KeyError("There is no matching improper in this view.") elif len(impropers) > 1: raise KeyError( - "More than one improper matches. Number of " f"matches is {len(impropers)}." + f"More than one improper matches. Number of matches is {len(impropers)}." ) return impropers[0] @@ -1573,7 +1573,6 @@ def _dynamics( vacuum=None, shift_delta=None, shift_coulomb=None, - coulomb_power=None, restraints=None, fixed=None, platform=None, @@ -1735,11 +1734,6 @@ def _dynamics( softening potential that smooths the creation and deletion of ghost atoms during a potential. This defaults to 1.0 A. - coulomb_power: int - The coulomb power parmeter that controls the electrostatic - softening potential that smooths the creation and deletion - of ghost atoms during a potential. This defaults to 0. - restraints: sire.mm.Restraints or list[sire.mm.Restraints] A single set of restraints, or a list of sets of restraints that will be applied to the atoms during @@ -1978,7 +1972,6 @@ def _dynamics( lambda_value=lambda_value, shift_delta=shift_delta, shift_coulomb=shift_coulomb, - coulomb_power=coulomb_power, swap_end_states=swap_end_states, ignore_perturbations=ignore_perturbations, restraints=restraints, @@ -2003,7 +1996,6 @@ def _minimisation( vacuum=None, shift_delta=None, shift_coulomb=None, - coulomb_power=None, platform=None, device=None, precision=None, @@ -2087,11 +2079,6 @@ def _minimisation( softening potential that smooths the creation and deletion of ghost atoms during a potential. This defaults to 1.0 A. - coulomb_power: int - The coulomb power parmeter that controls the electrostatic - softening potential that smooths the creation and deletion - of ghost atoms during a potential. This defaults to 0. - restraints: sire.mm.Restraints or list[sire.mm.Restraints] A single set of restraints, or a list of sets of restraints that will be applied to the atoms during @@ -2211,7 +2198,6 @@ def _minimisation( ignore_perturbations=ignore_perturbations, shift_delta=shift_delta, shift_coulomb=shift_coulomb, - coulomb_power=coulomb_power, restraints=restraints, fixed=fixed, map=map, diff --git a/src/sire/mol/_dynamics.py b/src/sire/mol/_dynamics.py index 975ce11ea..65feefb4a 100644 --- a/src/sire/mol/_dynamics.py +++ b/src/sire/mol/_dynamics.py @@ -1613,7 +1613,6 @@ def __init__( ignore_perturbations=None, shift_delta=None, shift_coulomb=None, - coulomb_power=None, restraints=None, fixed=None, qm_engine=None, @@ -1643,7 +1642,6 @@ def __init__( if shift_coulomb is not None: _add_extra(extras, "shift_coulomb", u(shift_coulomb)) - _add_extra(extras, "coulomb_power", coulomb_power) _add_extra(extras, "restraints", restraints) _add_extra(extras, "fixed", fixed) _add_extra(extras, "qm_engine", qm_engine) diff --git a/src/sire/mol/_minimisation.py b/src/sire/mol/_minimisation.py index b6c89b9fe..088256d36 100644 --- a/src/sire/mol/_minimisation.py +++ b/src/sire/mol/_minimisation.py @@ -22,7 +22,6 @@ def __init__( ignore_perturbations=None, shift_delta=None, shift_coulomb=None, - coulomb_power=None, restraints=None, fixed=None, ): @@ -46,7 +45,6 @@ def __init__( if shift_coulomb is not None: _add_extra(extras, "shift_coulomb", u(shift_coulomb)) - _add_extra(extras, "coulomb_power", coulomb_power) _add_extra(extras, "restraints", restraints) _add_extra(extras, "fixed", fixed) diff --git a/src/sire/system/_system.py b/src/sire/system/_system.py index 0846f535b..37a7175ee 100644 --- a/src/sire/system/_system.py +++ b/src/sire/system/_system.py @@ -440,11 +440,6 @@ def minimisation(self, *args, **kwargs): creation and deletion of ghost atoms during a potential. This defaults to 2.0 A. - coulomb_power: int - The coulomb power parmeter that controls the electrostatic - softening potential that smooths the creation and deletion - of ghost atoms during a potential. This defaults to 0. - vacuum: bool Whether or not to run the simulation in vacuum. If this is set to `True`, then the simulation space automatically be @@ -626,11 +621,6 @@ def dynamics(self, *args, **kwargs): creation and deletion of ghost atoms during a potential. This defaults to 2.0 A. - coulomb_power: int - The coulomb power parmeter that controls the electrostatic - softening potential that smooths the creation and deletion - of ghost atoms during a potential. This defaults to 0. - restraints: sire.mm.Restraints or list[sire.mm.Restraints] A single set of restraints, or a list of sets of restraints that will be applied to the atoms during @@ -948,7 +938,7 @@ def set_energy_trajectory(self, trajectory, map=None): if trajectory.what() != "SireMaths::EnergyTrajectory": raise TypeError( - f"You cannot set a {type(trajectory)} as an " "energy trajectory!" + f"You cannot set a {type(trajectory)} as an energy trajectory!" ) self._system.set_property(traj_propname.source(), trajectory) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 6147f05da..6704e051a 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1380,29 +1380,23 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, use_taylor_softening = not map["use_zacharias_softening"].value().asABoolean(); } - int coulomb_power = 0; + // use_beutler_softening overrides taylor/zacharias if set + bool use_beutler_softening = false; - if (map.specified("coulomb_power")) + if (map.specified("use_beutler_softening")) { - coulomb_power = map["coulomb_power"].value().asAnInteger(); + use_beutler_softening = map["use_beutler_softening"].value().asABoolean(); } - if (coulomb_power < 0) - coulomb_power = 0; - else if (coulomb_power > 4) - coulomb_power = 4; + double beutler_alpha = 0.5; - auto coulomb_power_expression = [](const QString &alpha, int power) + if (map.specified("beutler_alpha")) { - if (power == 0) - return QString("1"); - else if (power == 1) - return QString("(1-%1)").arg(alpha); - else if (power == 2) - return QString("(1-%1)*(1-%1)").arg(alpha); - else - return QString("(1-%1)^%2").arg(alpha).arg(power); - }; + beutler_alpha = map["beutler_alpha"].value().asADouble(); + } + + if (beutler_alpha < 0.0) + beutler_alpha = 0.0; auto taylor_power_expression = [](const QString &alpha, int power) { @@ -1419,15 +1413,32 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // see below for the description of this energy expression std::string nb14_expression, coul_expression, lj_expression; - if (use_taylor_softening) + if (use_beutler_softening) { + // Beutler et al., Chem. Phys. Lett., 1994 + // V_{LJ}(r) = (1-alpha) * 4 epsilon [ + // sigma^12 / (beutler_alpha*sigma^6*alpha + r^6)^2 + // - sigma^6 / (beutler_alpha*sigma^6*alpha + r^6) ] + // V_{coul}(r) = q_i q_j / 4 pi eps_0 sqrt(delta + r^2) + // delta = shift_coulomb^2 * alpha nb14_expression = QString( "coul_nrg+lj_nrg;" - "coul_nrg=138.9354558466661*q*(((%1)/sqrt((%2*alpha)+r_safe^2))-(kappa/r_safe));" + "coul_nrg=138.9354558466661*q*((1/sqrt((%1*alpha)+r_safe^2))-(kappa/r_safe));" + "lj_nrg=(1-alpha)*four_epsilon*sig6*(sig6-1);" + "sig6=(sigma^6)/(%2*sigma^6*alpha + r_safe^6);" + "r_safe=max(r, 0.001);") + .arg(shift_coulomb) + .arg(beutler_alpha) + .toStdString(); + } + else if (use_taylor_softening) + { + nb14_expression = QString( + "coul_nrg+lj_nrg;" + "coul_nrg=138.9354558466661*q*((1/sqrt((%1*alpha)+r_safe^2))-(kappa/r_safe));" "lj_nrg=four_epsilon*sig6*(sig6-1);" - "sig6=(sigma^6)/(%3*sigma^6 + r_safe^6);" + "sig6=(sigma^6)/(%2*sigma^6 + r_safe^6);" "r_safe=max(r, 0.001);") - .arg(coulomb_power_expression("alpha", coulomb_power)) .arg(shift_coulomb) .arg(taylor_power_expression("alpha", taylor_power)) .toStdString(); @@ -1436,12 +1447,11 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, { nb14_expression = QString( "coul_nrg+lj_nrg;" - "coul_nrg=138.9354558466661*q*(((%1)/sqrt((%2*alpha)+r_safe^2))-(kappa/r_safe));" + "coul_nrg=138.9354558466661*q*((1/sqrt((%1*alpha)+r_safe^2))-(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=%3*alpha;") - .arg(coulomb_power_expression("alpha", coulomb_power)) + "delta=%2*alpha;") .arg(shift_coulomb) .arg(shift_delta.to(SireUnits::nanometer)) .toStdString(); @@ -1460,17 +1470,51 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // periodic boundaries or cutoffs ghost_14ff->setUsesPeriodicBoundaryConditions(false); - if (use_taylor_softening) + if (use_beutler_softening) + { + // Beutler et al., Chem. Phys. Lett., 1994 + // + // V_{LJ}(r) = (1-alpha) * 4 epsilon [ + // sigma^12 / (beutler_alpha*sigma^6*alpha + r^6)^2 + // - sigma^6 / (beutler_alpha*sigma^6*alpha + r^6) ] + // + // V_{coul}(r) = q_i q_j / 4 pi eps_0 sqrt(delta + r^2) + // + // delta = shift_coulomb^2 * alpha + // + // Note that we supply half_sigma and two_sqrt_epsilon to save some + // cycles + // + // Note also that we subtract the normal coulomb energy as this + // is calculated during the standard NonbondedForce + // + // 138.9354558466661 is the constant needed to get energies in + // kJ mol-1 given the units of charge (|e|) and distance (nm) + // + coul_expression = QString("138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + "r_safe=max(r, 0.001);" + "max_kappa=max(kappa1, kappa2);" + "max_alpha=max(alpha1, alpha2);") + .arg(shift_coulomb) + .toStdString(); + + lj_expression = QString("(1-max_alpha)*two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" + "sig6=(sigma^6)/(%1*sigma^6*max_alpha + r_safe^6);" + "r_safe=max(r, 0.001);" + "max_alpha=max(alpha1, alpha2);" + "sigma=half_sigma1+half_sigma2;") + .arg(beutler_alpha) + .toStdString(); + } + else if (use_taylor_softening) { - // this uses the following potentials - // Zacharias and McCammon, J. Chem. Phys., 1994, and also, - // Michel et al., JCTC, 2007 - // LJ is Rich Taylor's softcore LJ + // Zacharias and McCammon, J. Chem. Phys., 1994, and also, + // Michel et al., JCTC, 2007; LJ is Rich Taylor's softcore LJ // // V_{LJ}(r) = 4 epsilon [ (sigma^12 / (alpha^m sigma^6 + r^6)^2) - // (sigma^6 / (alpha^m sigma^6 + r^6) ) ] // - // V_{coul}(r) = (1-alpha)^n q_i q_j / 4 pi eps_0 (delta+r^2)^(1/2) + // V_{coul}(r) = q_i q_j / 4 pi eps_0 (delta+r^2)^(1/2) // // delta = shift_coulomb^2 * alpha // @@ -1483,11 +1527,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // 138.9354558466661 is the constant needed to get energies in // kJ mol-1 given the units of charge (|e|) and distance (nm) // - coul_expression = QString("138.9354558466661*q1*q2*(((%1)/sqrt((%2*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + coul_expression = QString("138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" "r_safe=max(r, 0.001);" "max_kappa=max(kappa1, kappa2);" "max_alpha=max(alpha1, alpha2);") - .arg(coulomb_power_expression("max_alpha", coulomb_power)) .arg(shift_coulomb) .toStdString(); @@ -1501,16 +1544,15 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, } else { - // this uses the following potentials - // Zacharias and McCammon, J. Chem. Phys., 1994, and also, - // Michel et al., JCTC, 2007 + // Zacharias and McCammon, J. Chem. Phys., 1994, and also, + // Michel et al., JCTC, 2007 // // V_{LJ}(r) = 4 epsilon [ ( sigma^12 / (delta*sigma + r^2)^6 ) - // ( sigma^6 / (delta*sigma + r^2)^3 ) ] // // delta = shift_delta * alpha // - // V_{coul}(r) = (1-alpha)^n q_i q_j / 4 pi eps_0 (delta+r^2)^(1/2) + // V_{coul}(r) = q_i q_j / 4 pi eps_0 (delta+r^2)^(1/2) // // delta = shift_coulomb^2 * alpha // @@ -1524,11 +1566,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // 138.9354558466661 is the constant needed to get energies in // kJ mol-1 given the units of charge (|e|) and distance (nm) // - coul_expression = QString("138.9354558466661*q1*q2*(((%1)/sqrt((%2*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + coul_expression = QString("138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" "r_safe=max(r, 0.001);" "max_kappa=max(kappa1, kappa2);" "max_alpha=max(alpha1, alpha2);") - .arg(coulomb_power_expression("max_alpha", coulomb_power)) .arg(shift_coulomb) .toStdString(); From ca6ab9039740213d7bc0d82529a9cc842a94115e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 30 Apr 2026 16:42:58 +0100 Subject: [PATCH 13/21] Name LRC forces. --- wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 6704e051a..7efb17247 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1668,6 +1668,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, { auto ghost_lrc_ff = new OpenMM::CustomVolumeForce("lrc_coeff/v"); ghost_lrc_ff->addGlobalParameter("lrc_coeff", 0.0); + ghost_lrc_ff->setName("GhostLRCForce"); ghost_lrc_ff->setForceGroup(force_group_counter); lambda_lever.setForceIndex("ghost-lrc", system.addForce(ghost_lrc_ff)); lambda_lever.setForceGroup("ghost-lrc", force_group_counter++); @@ -1684,6 +1685,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, { auto background_lrc_ff = new OpenMM::CustomVolumeForce("lrc_background/v"); background_lrc_ff->addGlobalParameter("lrc_background", 0.0); + background_lrc_ff->setName("BackgroundLRCForce"); background_lrc_ff->setForceGroup(force_group_counter); lambda_lever.setForceIndex("background-lrc", system.addForce(background_lrc_ff)); lambda_lever.setForceGroup("background-lrc", force_group_counter++); @@ -1699,6 +1701,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, gcmc_lrc_ff->addGlobalParameter("n_w", 0.0); gcmc_lrc_ff->addGlobalParameter("lrc_w_solute", 0.0); gcmc_lrc_ff->addGlobalParameter("lrc_ww_half", 0.0); + gcmc_lrc_ff->setName("GCMCLRCForce"); gcmc_lrc_ff->setForceGroup(force_group_counter); lambda_lever.setForceIndex("gcmc-lrc", system.addForce(gcmc_lrc_ff)); lambda_lever.setForceGroup("gcmc-lrc", force_group_counter++); From 31fd5d859c7f5bde0386be936fd93fecbec1a6a2 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 1 May 2026 11:31:44 +0100 Subject: [PATCH 14/21] Add lrc_scale lever. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 8 ++++++++ wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp | 3 ++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index fd10875b2..312146d7d 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -2053,6 +2053,14 @@ double LambdaLever::setLambda(OpenMM::Context &context, } context.setParameter("lrc_coeff", lrc_coeff); + + // lrc_scale defaults to 1.0 (no effect). Schedules that fix epsilon + // (e.g. Beutler softcore) should set a force-specific equation for + // lever "lrc_scale" on force "ghost-lrc" to scale it to zero as the + // ghost is annihilated/decoupled. + double lrc_scale = this->lambda_schedule.morph( + "ghost-lrc", "lrc_scale", 1.0, 1.0, lambda_value); + context.setParameter("lrc_scale", lrc_scale); } // Update the NonbondedForce (background) LRC via its own CustomVolumeForce. diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 7efb17247..1081db3ad 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1666,8 +1666,9 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // a cached closed-form sum over interaction-group pairs. if (use_dispersion_correction && ffinfo.hasCutoff() && ffinfo.space().isPeriodic()) { - auto ghost_lrc_ff = new OpenMM::CustomVolumeForce("lrc_coeff/v"); + auto ghost_lrc_ff = new OpenMM::CustomVolumeForce("lrc_coeff*lrc_scale/v"); ghost_lrc_ff->addGlobalParameter("lrc_coeff", 0.0); + ghost_lrc_ff->addGlobalParameter("lrc_scale", 1.0); ghost_lrc_ff->setName("GhostLRCForce"); ghost_lrc_ff->setForceGroup(force_group_counter); lambda_lever.setForceIndex("ghost-lrc", system.addForce(ghost_lrc_ff)); From 2a5b8f8c6450481a47b8b48cd850b0bf6bb61d98 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 5 May 2026 09:28:50 +0100 Subject: [PATCH 15/21] Refactor soft-core expression definitions to avoid duplication. --- .../SireOpenMM/sire_to_openmm_system.cpp | 81 ++++++------------- 1 file changed, 23 insertions(+), 58 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 1081db3ad..e52cd19ce 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1470,6 +1470,25 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // periodic boundaries or cutoffs ghost_14ff->setUsesPeriodicBoundaryConditions(false); + // Coulomb soft-core expression is identical for all three soft-core forms + // (Beutler, Taylor, and Zacharias/Michel). The energy is: + // + // V_{coul}(r) = q_i q_j / 4 pi eps_0 * sqrt(delta + r^2) + // + // delta = shift_coulomb^2 * max(alpha_i, alpha_j) + // + // We use max(alpha) so that the interaction is soft whenever either + // particle is a ghost. We subtract the regular Coulomb (kappa/r) term + // because the standard NonbondedForce already contributes it. + // 138.9354558466661 converts from |e|^2/nm to kJ mol-1. + // + coul_expression = QString("138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + "r_safe=max(r, 0.001);" + "max_kappa=max(kappa1, kappa2);" + "max_alpha=max(alpha1, alpha2);") + .arg(shift_coulomb) + .toStdString(); + if (use_beutler_softening) { // Beutler et al., Chem. Phys. Lett., 1994 @@ -1478,26 +1497,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // sigma^12 / (beutler_alpha*sigma^6*alpha + r^6)^2 // - sigma^6 / (beutler_alpha*sigma^6*alpha + r^6) ] // - // V_{coul}(r) = q_i q_j / 4 pi eps_0 sqrt(delta + r^2) - // - // delta = shift_coulomb^2 * alpha - // - // Note that we supply half_sigma and two_sqrt_epsilon to save some - // cycles - // - // Note also that we subtract the normal coulomb energy as this - // is calculated during the standard NonbondedForce + // half_sigma and two_sqrt_epsilon are supplied to save cycles. // - // 138.9354558466661 is the constant needed to get energies in - // kJ mol-1 given the units of charge (|e|) and distance (nm) - // - coul_expression = QString("138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" - "r_safe=max(r, 0.001);" - "max_kappa=max(kappa1, kappa2);" - "max_alpha=max(alpha1, alpha2);") - .arg(shift_coulomb) - .toStdString(); - lj_expression = QString("(1-max_alpha)*two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" "sig6=(sigma^6)/(%1*sigma^6*max_alpha + r_safe^6);" "r_safe=max(r, 0.001);" @@ -1514,26 +1515,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // V_{LJ}(r) = 4 epsilon [ (sigma^12 / (alpha^m sigma^6 + r^6)^2) - // (sigma^6 / (alpha^m sigma^6 + r^6) ) ] // - // V_{coul}(r) = q_i q_j / 4 pi eps_0 (delta+r^2)^(1/2) - // - // delta = shift_coulomb^2 * alpha - // - // Note that we supply half_sigma and two_sqrt_epsilon to save some - // cycles + // half_sigma and two_sqrt_epsilon are supplied to save cycles. // - // Note also that we subtract the normal coulomb energy as this - // is calculated during the standard NonbondedForce - // - // 138.9354558466661 is the constant needed to get energies in - // kJ mol-1 given the units of charge (|e|) and distance (nm) - // - coul_expression = QString("138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" - "r_safe=max(r, 0.001);" - "max_kappa=max(kappa1, kappa2);" - "max_alpha=max(alpha1, alpha2);") - .arg(shift_coulomb) - .toStdString(); - lj_expression = QString("two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" "sig6=(sigma^6)/(%1*sigma^6 + r_safe^6);" "r_safe=max(r, 0.001);" @@ -1552,27 +1535,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // // delta = shift_delta * alpha // - // V_{coul}(r) = q_i q_j / 4 pi eps_0 (delta+r^2)^(1/2) - // - // delta = shift_coulomb^2 * alpha - // - // Note that we pre-calculate delta as a forcefield parameter, - // and also supply half_sigma and two_sqrt_epsilon to save some - // cycles - // - // Note also that we subtract the normal coulomb energy as this - // is calculated during the standard NonbondedForce + // half_sigma and two_sqrt_epsilon are supplied to save cycles. // - // 138.9354558466661 is the constant needed to get energies in - // kJ mol-1 given the units of charge (|e|) and distance (nm) - // - coul_expression = QString("138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" - "r_safe=max(r, 0.001);" - "max_kappa=max(kappa1, kappa2);" - "max_alpha=max(alpha1, alpha2);") - .arg(shift_coulomb) - .toStdString(); - lj_expression = QString("two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" "sig6=(sigma^6)/(((sigma*delta) + r_safe^2)^3);" "delta=%1*max_alpha;" @@ -1613,6 +1577,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // (proportional to int r^2 * 1/r^3 dr = int 1/r dr) log-divergent. ghost_ghost_coulombff->setUseLongRangeCorrection(false); ghost_nonghost_coulombff->setUseLongRangeCorrection(false); + // LJ soft-core LRC is handled analytically via a CustomVolumeForce // rather than OpenMM's numerical integrator, because the standard // LJ tail (r > rc, soft-core shift negligible) has a closed-form From 499dd22374836de29a1281c3ac2d58c7b206d8a5 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 6 May 2026 11:51:29 +0100 Subject: [PATCH 16/21] Update method name for clarity. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 4 ++-- wrapper/Convert/SireOpenMM/lambdalever.h | 2 +- wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 312146d7d..38c1ef3a8 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -438,7 +438,7 @@ bool LambdaLever::wasForceChanged(const QString &name) const return it.value(); } -void LambdaLever::setGCMCWaterAtoms(const QVector &atoms) +void LambdaLever::setWaterAtoms(const QVector &atoms) { gcmc_water_atoms = QSet(atoms.begin(), atoms.end()); } @@ -1977,7 +1977,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, // Update the ghost LJ dispersion correction via the CustomVolumeForce. // At r > rc the soft-core shift is negligible, so the standard closed-form - // LJ tail integral applies. Results are cached per lambda state. + // LJ tail integral applies. Results are cached per lambda state. auto ghost_lrc_ff = this->getForce("ghost-lrc", system); if (ghost_lrc_ff != nullptr && ghost_ghost_ljff != nullptr && ghost_nonghost_ljff != nullptr) { diff --git a/wrapper/Convert/SireOpenMM/lambdalever.h b/wrapper/Convert/SireOpenMM/lambdalever.h index 0e2cbe494..927a24b03 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.h +++ b/wrapper/Convert/SireOpenMM/lambdalever.h @@ -165,7 +165,7 @@ namespace SireOpenMM QStringList getForceNames() const; bool wasForceChanged(const QString &name) const; - void setGCMCWaterAtoms(const QVector &atoms); + void setWaterAtoms(const QVector &atoms); protected: void updateRestraintInContext(OpenMM::Force &ff, double rho, diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index e52cd19ce..bd6186d1f 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -2299,7 +2299,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, water_atom_indices.append(j); } - lambda_lever.setGCMCWaterAtoms(water_atom_indices); + lambda_lever.setWaterAtoms(water_atom_indices); // Pre-compute lrc_w_solute (per active water molecule, interaction with // all non-water atoms: protein, ligand, ions) and lrc_ww_half (per active From 457aff5e4402bccf7c7ff5fad9be47f093ef5688 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 6 May 2026 14:39:06 +0100 Subject: [PATCH 17/21] Recombine ghost nonbonded forces to improve performance. --- src/sire/morph/_xml.py | 18 +- tests/convert/test_openmm_rest2.py | 54 ++- wrapper/Convert/SireOpenMM/lambdalever.cpp | 147 +++---- wrapper/Convert/SireOpenMM/lambdalever.h | 8 +- .../SireOpenMM/sire_to_openmm_system.cpp | 379 ++++++++---------- 5 files changed, 257 insertions(+), 349 deletions(-) diff --git a/src/sire/morph/_xml.py b/src/sire/morph/_xml.py index 777d2482a..29bee2a1d 100644 --- a/src/sire/morph/_xml.py +++ b/src/sire/morph/_xml.py @@ -20,8 +20,7 @@ def evaluate_xml_force(mols, xml, force): force : str The name of the custom force to evaluate. Options are: - "ghost-ghost-lj", "ghost-ghost-coulomb", - "ghost-nonghost-lj", "ghost-nonghost-coulomb", "ghost-14". + "ghost-ghost", "ghost-nonghost", "ghost-14". Returns ------- @@ -78,24 +77,19 @@ def evaluate_xml_force(mols, xml, force): # Validate the force name. valid = [ - "ghostghostlj", - "ghostghostcoulomb", - "ghostnonghostlj", - "ghostnonghostcoulomb", + "ghostghost", + "ghostnonghost", "ghost14", ] if force not in valid: raise ValueError( - "'force' must be one of 'ghost-ghost-lj', 'ghost-ghost-coulomb', " - "'ghost-nonghost-lj', 'ghost-nonghost-coulomb', or 'ghost-14'." + "'force' must be one of 'ghost-ghost', 'ghost-nonghost', or 'ghost-14'." ) # Map sanitised name to the OpenMM force name in the XML. _force_name_map = { - "ghostghostlj": "GhostGhostLJForce", - "ghostghostcoulomb": "GhostGhostCoulombForce", - "ghostnonghostlj": "GhostNonGhostLJForce", - "ghostnonghostcoulomb": "GhostNonGhostCoulombForce", + "ghostghost": "GhostGhostNonbondedForce", + "ghostnonghost": "GhostNonGhostNonbondedForce", "ghost14": "Ghost14BondForce", } name = _force_name_map[force] diff --git a/tests/convert/test_openmm_rest2.py b/tests/convert/test_openmm_rest2.py index af4338044..7ac498a52 100644 --- a/tests/convert/test_openmm_rest2.py +++ b/tests/convert/test_openmm_rest2.py @@ -91,34 +91,32 @@ def test_rest2(mols, rest2_selection, excluded_atoms, request): if is_perturbable: # Find the ghost/ghost nonbonded interaction. for force in omm_system.getForces(): - if force.getName() == "GhostGhostCoulombForce": - ghost_ghost_coulomb_force = force - elif force.getName() == "GhostGhostLJForce": - ghost_ghost_lj_force = force + if force.getName() == "GhostGhostNonbondedForce": + break # Store the initial parameters. ghost_ghost_params_initial = [] excluded_ghost_ghost_indices = [] - for i in range(ghost_ghost_coulomb_force.getNumParticles()): - charge, _, _ = ghost_ghost_coulomb_force.getParticleParameters(i) - _, two_sqrt_epsilon, _ = ghost_ghost_lj_force.getParticleParameters(i) + for i in range(force.getNumParticles()): + charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( + force.getParticleParameters(i) + ) ghost_ghost_params_initial.append((charge, two_sqrt_epsilon)) if i in excluded_atoms: excluded_ghost_ghost_indices.append(i) # Find the ghost/non-ghost nonbonded interaction. for force in omm_system.getForces(): - if force.getName() == "GhostNonGhostCoulombForce": - ghost_non_ghost_coulomb_force = force - elif force.getName() == "GhostNonGhostLJForce": - ghost_non_ghost_lj_force = force + if force.getName() == "GhostNonGhostNonbondedForce": + break # Store the initial parameters. ghost_non_ghost_params_initial = [] excluded_ghost_non_ghost_indices = [] - for i in range(ghost_non_ghost_coulomb_force.getNumParticles()): - charge, _, _ = ghost_non_ghost_coulomb_force.getParticleParameters(i) - _, two_sqrt_epsilon, _ = ghost_non_ghost_lj_force.getParticleParameters(i) + for i in range(force.getNumParticles()): + charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( + force.getParticleParameters(i) + ) ghost_non_ghost_params_initial.append((charge, two_sqrt_epsilon)) if i in excluded_atoms: excluded_ghost_non_ghost_indices.append(i) @@ -153,34 +151,32 @@ def test_rest2(mols, rest2_selection, excluded_atoms, request): for i in range(force.getNumExceptions()): exception_params_modified.append(force.getExceptionParameters(i)[-3::2]) - # Find the ghost/ghost nonbonded forces. + # Find the ghost/ghost nonbonded interaction. for force in omm_system.getForces(): - if force.getName() == "GhostGhostCoulombForce": - ghost_ghost_coulomb_force = force - elif force.getName() == "GhostGhostLJForce": - ghost_ghost_lj_force = force + if force.getName() == "GhostGhostNonbondedForce": + break # Handle custom forces for pertubable molecules. if is_perturbable: # Store the modified parameters. ghost_ghost_params_modified = [] - for i in range(ghost_ghost_coulomb_force.getNumParticles()): - charge, _, _ = ghost_ghost_coulomb_force.getParticleParameters(i) - _, two_sqrt_epsilon, _ = ghost_ghost_lj_force.getParticleParameters(i) + for i in range(force.getNumParticles()): + charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( + force.getParticleParameters(i) + ) ghost_ghost_params_modified.append((charge, two_sqrt_epsilon)) # Find the ghost/non-ghost nonbonded interaction. for force in omm_system.getForces(): - if force.getName() == "GhostNonGhostCoulombForce": - ghost_non_ghost_coulomb_force = force - elif force.getName() == "GhostNonGhostLJForce": - ghost_non_ghost_lj_force = force + if force.getName() == "GhostNonGhostNonbondedForce": + break # Store the modified parameters. ghost_non_ghost_params_modified = [] - for i in range(ghost_non_ghost_coulomb_force.getNumParticles()): - charge, _, _ = ghost_non_ghost_coulomb_force.getParticleParameters(i) - _, two_sqrt_epsilon, _ = ghost_non_ghost_lj_force.getParticleParameters(i) + for i in range(force.getNumParticles()): + charge, half_sigma, two_sqrt_epsilon, alpha, kappa = ( + force.getParticleParameters(i) + ) ghost_non_ghost_params_modified.append((charge, two_sqrt_epsilon)) # Store the scaling factor. diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index 38c1ef3a8..f30799558 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -438,7 +438,7 @@ bool LambdaLever::wasForceChanged(const QString &name) const return it.value(); } -void LambdaLever::setWaterAtoms(const QVector &atoms) +void LambdaLever::setGCMCWaterAtoms(const QVector &atoms) { gcmc_water_atoms = QSet(atoms.begin(), atoms.end()); } @@ -1296,10 +1296,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, // get copies of the forcefields in which the parameters will be changed auto qmff = this->getForce("qmff", system); auto cljff = this->getForce("clj", system); - auto ghost_ghost_coulombff = this->getForce("ghost/ghost-coulomb", system); - auto ghost_ghost_ljff = this->getForce("ghost/ghost-lj", system); - auto ghost_nonghost_coulombff = this->getForce("ghost/non-ghost-coulomb", system); - auto ghost_nonghost_ljff = this->getForce("ghost/non-ghost-lj", system); + 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 bondff = this->getForce("bond", system); auto angff = this->getForce("angle", system); @@ -1307,13 +1305,12 @@ double LambdaLever::setLambda(OpenMM::Context &context, auto cmapff = this->getForce("cmap", system); // we know if we have peturbable ghost atoms if we have the ghost forcefields - const bool have_ghost_atoms = (ghost_ghost_ljff != 0 or ghost_nonghost_ljff != 0); + const bool have_ghost_atoms = (ghost_ghostff != 0 or ghost_nonghostff != 0); // whether the constraints have changed bool have_constraints_changed = false; - std::vector custom_params_coul = {0.0, 0.0, 0.0}; - std::vector custom_params_lj = {0.0, 0.0, 0.0}; + std::vector custom_params = {0.0, 0.0, 0.0, 0.0, 0.0}; if (qmff != 0) { @@ -1326,8 +1323,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, // track whether parameters actually changed for each force, so we only // call updateParametersInContext when necessary bool has_changed_cljff = false; - bool has_changed_coulombff = false; - bool has_changed_ljff = false; + bool has_changed_ghostff = false; bool has_changed_ghost14ff = false; bool has_changed_bondff = false; bool has_changed_angff = false; @@ -1501,9 +1497,7 @@ double LambdaLever::setLambda(OpenMM::Context &context, // Detect whether any CLJ or ghost-14 parameters changed has_changed_cljff |= rest2_changed || cache.hasChanged("clj", "charge") || cache.hasChanged("clj", "sigma") || cache.hasChanged("clj", "epsilon") || cache.hasChanged("clj", "alpha") || cache.hasChanged("clj", "kappa") || cache.hasChanged("clj", "charge_scale") || cache.hasChanged("clj", "lj_scale"); - has_changed_coulombff |= rest2_changed || cache.hasChanged("ghost/ghost", "charge") || cache.hasChanged("ghost/ghost", "alpha") || cache.hasChanged("ghost/ghost", "kappa") || cache.hasChanged("ghost/non-ghost", "charge") || cache.hasChanged("ghost/non-ghost", "alpha") || cache.hasChanged("ghost/non-ghost", "kappa"); - - has_changed_ljff |= rest2_changed || cache.hasChanged("ghost/ghost", "sigma") || cache.hasChanged("ghost/ghost", "epsilon") || cache.hasChanged("ghost/ghost", "alpha") || cache.hasChanged("ghost/non-ghost", "sigma") || cache.hasChanged("ghost/non-ghost", "epsilon") || cache.hasChanged("ghost/non-ghost", "alpha"); + has_changed_ghostff |= rest2_changed || cache.hasChanged("ghost/ghost", "charge") || cache.hasChanged("ghost/ghost", "sigma") || cache.hasChanged("ghost/ghost", "epsilon") || cache.hasChanged("ghost/ghost", "alpha") || cache.hasChanged("ghost/ghost", "kappa") || cache.hasChanged("ghost/non-ghost", "charge") || cache.hasChanged("ghost/non-ghost", "sigma") || cache.hasChanged("ghost/non-ghost", "epsilon") || cache.hasChanged("ghost/non-ghost", "alpha") || cache.hasChanged("ghost/non-ghost", "kappa"); has_changed_ghost14ff |= rest2_changed || cache.hasChanged("ghost-14", "charge") || cache.hasChanged("ghost-14", "sigma") || cache.hasChanged("ghost-14", "epsilon") || cache.hasChanged("ghost-14", "alpha") || cache.hasChanged("ghost-14", "kappa") || cache.hasChanged("ghost-14", "charge_scale") || cache.hasChanged("ghost-14", "lj_scale"); @@ -1525,68 +1519,54 @@ double LambdaLever::setLambda(OpenMM::Context &context, } // reduced charge - custom_params_coul[0] = sqrt_scale * morphed_ghost_charges[j]; - // alpha - custom_params_coul[1] = morphed_ghost_alphas[j]; - // kappa - custom_params_coul[2] = morphed_ghost_kappas[j]; - + custom_params[0] = sqrt_scale * morphed_ghost_charges[j]; // half_sigma - custom_params_lj[0] = 0.5 * morphed_ghost_sigmas[j]; + custom_params[1] = 0.5 * morphed_ghost_sigmas[j]; // two_sqrt_epsilon - custom_params_lj[1] = 2.0 * sqrt_scale * std::sqrt(morphed_ghost_epsilons[j]); + custom_params[2] = 2.0 * sqrt_scale * std::sqrt(morphed_ghost_epsilons[j]); // alpha - custom_params_lj[2] = morphed_ghost_alphas[j]; + custom_params[3] = morphed_ghost_alphas[j]; + // kappa + custom_params[4] = morphed_ghost_kappas[j]; // clamp alpha between 0 and 1 - if (custom_params_coul[1] < 0) - { - custom_params_coul[1] = 0; - custom_params_lj[2] = 0; - } - else if (custom_params_coul[1] > 1) - { - custom_params_coul[1] = 1; - custom_params_lj[2] = 1; - } + if (custom_params[3] < 0) + custom_params[3] = 0; + else if (custom_params[3] > 1) + custom_params[3] = 1; // clamp kappa between 0 and 1 - if (custom_params_coul[2] < 0) - custom_params_coul[2] = 0; - else if (custom_params_coul[2] > 1) - custom_params_coul[2] = 1; + if (custom_params[4] < 0) + custom_params[4] = 0; + else if (custom_params[4] > 1) + custom_params[4] = 1; - ghost_ghost_coulombff->setParticleParameters(start_index + j, custom_params_coul); - ghost_ghost_ljff->setParticleParameters(start_index + j, custom_params_lj); + ghost_ghostff->setParticleParameters(start_index + j, custom_params); // reduced charge - custom_params_coul[0] = sqrt_scale * morphed_nonghost_charges[j]; - // alpha - custom_params_coul[1] = morphed_nonghost_alphas[j]; - // kappa - custom_params_coul[2] = morphed_nonghost_kappas[j]; - + custom_params[0] = sqrt_scale * morphed_nonghost_charges[j]; // half_sigma - custom_params_lj[0] = 0.5 * morphed_nonghost_sigmas[j]; + custom_params[1] = 0.5 * morphed_nonghost_sigmas[j]; // two_sqrt_epsilon - custom_params_lj[1] = 2.0 * sqrt_scale * std::sqrt(morphed_nonghost_epsilons[j]); + custom_params[2] = 2.0 * sqrt_scale * std::sqrt(morphed_nonghost_epsilons[j]); // alpha - custom_params_lj[2] = morphed_nonghost_alphas[j]; + custom_params[3] = morphed_nonghost_alphas[j]; + // kappa + custom_params[4] = morphed_nonghost_kappas[j]; // clamp alpha between 0 and 1 - if (custom_params_coul[1] < 0) - { - custom_params_coul[1] = 0; - custom_params_lj[2] = 0; - } - else if (custom_params_coul[1] > 1) - { - custom_params_coul[1] = 1; - custom_params_lj[2] = 1; - } + if (custom_params[3] < 0) + custom_params[3] = 0; + else if (custom_params[3] > 1) + custom_params[3] = 1; + + // clamp kappa between 0 and 1 + if (custom_params[4] < 0) + custom_params[4] = 0; + else if (custom_params[4] > 1) + custom_params[4] = 1; - ghost_nonghost_coulombff->setParticleParameters(start_index + j, custom_params_coul); - ghost_nonghost_ljff->setParticleParameters(start_index + j, custom_params_lj); + ghost_nonghostff->setParticleParameters(start_index + j, custom_params); if (is_from_ghost or is_to_ghost) { @@ -1952,34 +1932,19 @@ double LambdaLever::setLambda(OpenMM::Context &context, if (has_changed_cljff and cljff) cljff->updateParametersInContext(context); - if (has_changed_coulombff or has_changed_ljff) + if (has_changed_ghostff) { - // Update the lambda cache key before any ghost force updateParametersInContext. - if (ghost_ghost_ljff or ghost_nonghost_ljff) - context.setParameter("lambda", std::round(lambda_value * 1e5) / 1e5); - - if (has_changed_coulombff) - { - if (ghost_ghost_coulombff) - ghost_ghost_coulombff->updateParametersInContext(context); - if (ghost_nonghost_coulombff) - ghost_nonghost_coulombff->updateParametersInContext(context); - } - - if (has_changed_ljff) - { - if (ghost_ghost_ljff) - ghost_ghost_ljff->updateParametersInContext(context); - if (ghost_nonghost_ljff) - ghost_nonghost_ljff->updateParametersInContext(context); - } + if (ghost_ghostff) + ghost_ghostff->updateParametersInContext(context); + if (ghost_nonghostff) + ghost_nonghostff->updateParametersInContext(context); } // Update the ghost LJ dispersion correction via the CustomVolumeForce. // At r > rc the soft-core shift is negligible, so the standard closed-form // LJ tail integral applies. Results are cached per lambda state. auto ghost_lrc_ff = this->getForce("ghost-lrc", system); - if (ghost_lrc_ff != nullptr && ghost_ghost_ljff != nullptr && ghost_nonghost_ljff != nullptr) + if (ghost_lrc_ff != nullptr && ghost_ghostff != nullptr && ghost_nonghostff != nullptr) { const qint64 lam_key = qRound64(lambda_value * 1e5); double lrc_coeff = 0.0; @@ -1990,23 +1955,23 @@ double LambdaLever::setLambda(OpenMM::Context &context, } else { - const double cutoff = ghost_ghost_ljff->getCutoffDistance(); + const double cutoff = ghost_ghostff->getCutoffDistance(); const double rc3 = cutoff * cutoff * cutoff; const double rc9 = rc3 * rc3 * rc3; const double four_pi = 4.0 * M_PI; // Interaction group sets: ghost atoms and non-ghost atoms. std::set ghost_set, dummy_set, nonghost_set; - ghost_ghost_ljff->getInteractionGroupParameters(0, ghost_set, dummy_set); - ghost_nonghost_ljff->getInteractionGroupParameters(0, dummy_set, nonghost_set); + ghost_ghostff->getInteractionGroupParameters(0, ghost_set, dummy_set); + ghost_nonghostff->getInteractionGroupParameters(0, dummy_set, nonghost_set); // Cache half_sigma and two_sqrt_epsilon for each ghost atom. QHash> ghost_params; for (int i : ghost_set) { std::vector p; - ghost_ghost_ljff->getParticleParameters(i, p); - ghost_params[i] = {p[0], p[1]}; // half_sigma, two_sqrt_epsilon + ghost_ghostff->getParticleParameters(i, p); + ghost_params[i] = {p[1], p[2]}; // half_sigma, two_sqrt_epsilon } // Cache non-ghost params. @@ -2014,8 +1979,8 @@ double LambdaLever::setLambda(OpenMM::Context &context, for (int j : nonghost_set) { std::vector p; - ghost_nonghost_ljff->getParticleParameters(j, p); - nonghost_params[j] = {p[0], p[1]}; + ghost_nonghostff->getParticleParameters(j, p); + nonghost_params[j] = {p[1], p[2]}; } // ghost-ghost unique pairs (i < j). @@ -2182,11 +2147,9 @@ double LambdaLever::setLambda(OpenMM::Context &context, // record which named forces had parameters changed in this call last_changed_forces["clj"] = has_changed_cljff; - last_changed_forces["ghost/ghost-coulomb"] = has_changed_coulombff; - last_changed_forces["ghost/ghost-lj"] = has_changed_ljff; - last_changed_forces["ghost/non-ghost-coulomb"] = has_changed_coulombff; - last_changed_forces["ghost/non-ghost-lj"] = has_changed_ljff; - last_changed_forces["ghost-lrc"] = has_changed_ljff; + last_changed_forces["ghost/ghost"] = has_changed_ghostff; + last_changed_forces["ghost/non-ghost"] = has_changed_ghostff; + last_changed_forces["ghost-lrc"] = has_changed_ghostff; last_changed_forces["background-lrc"] = has_changed_cljff; last_changed_forces["gcmc-lrc"] = false; last_changed_forces["ghost-14"] = has_changed_ghost14ff; diff --git a/wrapper/Convert/SireOpenMM/lambdalever.h b/wrapper/Convert/SireOpenMM/lambdalever.h index 927a24b03..4fd03aabf 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.h +++ b/wrapper/Convert/SireOpenMM/lambdalever.h @@ -165,7 +165,7 @@ namespace SireOpenMM QStringList getForceNames() const; bool wasForceChanged(const QString &name) const; - void setWaterAtoms(const QVector &atoms); + void setGCMCWaterAtoms(const QVector &atoms); protected: void updateRestraintInContext(OpenMM::Force &ff, double rho, @@ -218,16 +218,16 @@ namespace SireOpenMM mutable double last_qmff_lam; /** Cache of pre-computed ghost LJ dispersion coefficients keyed by - * rounded lambda (qRound64(lambda * 1e5)). Populated on first visit + * rounded lambda (qRound64(lambda * 1e5)). Populated on first visit * to each lambda state and reused on warm passes. */ mutable QHash lrc_coeff_cache; /** Cache of pre-computed background (NonbondedForce) LJ dispersion - * coefficients keyed by rounded lambda. Mirrors lrc_coeff_cache but + * coefficients keyed by rounded lambda. Mirrors lrc_coeff_cache but * covers all non-ghost atoms in the clj force. */ mutable QHash background_lrc_coeff_cache; - /** OpenMM atom indices that belong to GCMC water molecules. When + /** OpenMM atom indices that belong to GCMC water molecules. When * non-empty these atoms are excluded from background-lrc (their LRC * is handled by the gcmc-lrc CustomVolumeForce instead). */ QSet gcmc_water_atoms; diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index bd6186d1f..642197e37 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -475,14 +475,12 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, const double internal_to_k = (1 * SireUnits::kcal_per_mol / (SireUnits::angstrom2)).to(SireUnits::kJ_per_mol / (SireUnits::nanometer2)); auto cljff = lambda_lever.getForce("clj", system); - auto ghost_ghost_coulombff = lambda_lever.getForce("ghost/ghost-coulomb", system); - auto ghost_ghost_ljff = lambda_lever.getForce("ghost/ghost-lj", system); - auto ghost_nonghost_coulombff = lambda_lever.getForce("ghost/non-ghost-coulomb", system); - auto ghost_nonghost_ljff = lambda_lever.getForce("ghost/non-ghost-lj", system); + auto ghost_ghostff = lambda_lever.getForce("ghost/ghost", system); + auto ghost_nonghostff = lambda_lever.getForce("ghost/non-ghost", system); std::vector custom_params = {1.0, 0.0, 0.0}; - // Define null parameters used to add these particles to the ghost forces (5 total) - std::vector custom_clj_params = {0.0, 0.0, 0.0}; + // Null parameters for anchor particles added to the ghost forces {q, half_sigma, two_sqrt_epsilon, alpha, kappa} + std::vector custom_clj_params = {0.0, 0.0, 0.0, 0.0, 0.0}; // we need to add all of the positions as anchor particles for (const auto &restraint : atom_restraints) @@ -517,16 +515,14 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, cljff->addParticle(0, 0, 0); } - if (ghost_ghost_coulombff != 0) + if (ghost_ghostff != 0) { - ghost_ghost_coulombff->addParticle(custom_clj_params); - ghost_ghost_ljff->addParticle(custom_clj_params); + ghost_ghostff->addParticle(custom_clj_params); } - if (ghost_nonghost_coulombff != 0) + if (ghost_nonghostff != 0) { - ghost_nonghost_coulombff->addParticle(custom_clj_params); - ghost_nonghost_ljff->addParticle(custom_clj_params); + ghost_nonghostff->addParticle(custom_clj_params); } } @@ -542,22 +538,20 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, cljff->addException(anchor_index, atom_index, 0, 0, 0, true); } - if (ghost_ghost_coulombff != 0) + if (ghost_ghostff != 0) { // make sure to exclude interactions between // the atom being positionally restrained and // the anchor - ghost_ghost_coulombff->addExclusion(anchor_index, atom_index); - ghost_ghost_ljff->addExclusion(anchor_index, atom_index); + ghost_ghostff->addExclusion(anchor_index, atom_index); } - if (ghost_nonghost_coulombff != 0) + if (ghost_nonghostff != 0) { // make sure to exclude interactions between // the atom being positionally restrained and // the anchor - ghost_nonghost_coulombff->addExclusion(anchor_index, atom_index); - ghost_nonghost_ljff->addExclusion(anchor_index, atom_index); + ghost_nonghostff->addExclusion(anchor_index, atom_index); } restraintff->addBond(anchor_index, atom_index, custom_params); @@ -1314,10 +1308,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, /// OpenMM::CustomBondForce *ghost_14ff = 0; - OpenMM::CustomNonbondedForce *ghost_ghost_coulombff = 0; - OpenMM::CustomNonbondedForce *ghost_ghost_ljff = 0; - OpenMM::CustomNonbondedForce *ghost_nonghost_coulombff = 0; - OpenMM::CustomNonbondedForce *ghost_nonghost_ljff = 0; + OpenMM::CustomNonbondedForce *ghost_ghostff = 0; + OpenMM::CustomNonbondedForce *ghost_nonghostff = 0; if (any_perturbable) { @@ -1411,7 +1403,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, }; // see below for the description of this energy expression - std::string nb14_expression, coul_expression, lj_expression; + std::string nb14_expression, clj_expression; if (use_beutler_softening) { @@ -1470,25 +1462,6 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // periodic boundaries or cutoffs ghost_14ff->setUsesPeriodicBoundaryConditions(false); - // Coulomb soft-core expression is identical for all three soft-core forms - // (Beutler, Taylor, and Zacharias/Michel). The energy is: - // - // V_{coul}(r) = q_i q_j / 4 pi eps_0 * sqrt(delta + r^2) - // - // delta = shift_coulomb^2 * max(alpha_i, alpha_j) - // - // We use max(alpha) so that the interaction is soft whenever either - // particle is a ghost. We subtract the regular Coulomb (kappa/r) term - // because the standard NonbondedForce already contributes it. - // 138.9354558466661 converts from |e|^2/nm to kJ mol-1. - // - coul_expression = QString("138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" - "r_safe=max(r, 0.001);" - "max_kappa=max(kappa1, kappa2);" - "max_alpha=max(alpha1, alpha2);") - .arg(shift_coulomb) - .toStdString(); - if (use_beutler_softening) { // Beutler et al., Chem. Phys. Lett., 1994 @@ -1499,133 +1472,143 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // // half_sigma and two_sqrt_epsilon are supplied to save cycles. // - lj_expression = QString("(1-max_alpha)*two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" - "sig6=(sigma^6)/(%1*sigma^6*max_alpha + r_safe^6);" - "r_safe=max(r, 0.001);" - "max_alpha=max(alpha1, alpha2);" - "sigma=half_sigma1+half_sigma2;") - .arg(beutler_alpha) - .toStdString(); + clj_expression = QString("coul_nrg+lj_nrg;" + "coul_nrg=138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + "lj_nrg=(1-max_alpha)*two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" + "sig6=(sigma^6)/(%2*sigma^6*max_alpha + r_safe^6);" + "r_safe=max(r, 0.001);" + "max_kappa=max(kappa1, kappa2);" + "max_alpha=max(alpha1, alpha2);" + "sigma=half_sigma1+half_sigma2;") + .arg(shift_coulomb) + .arg(beutler_alpha) + .toStdString(); } else if (use_taylor_softening) { - // Zacharias and McCammon, J. Chem. Phys., 1994, and also, - // Michel et al., JCTC, 2007; LJ is Rich Taylor's softcore LJ + // this uses the following potentials + // Zacharias and McCammon, J. Chem. Phys., 1994, and also, + // Michel et al., JCTC, 2007 + // LJ is Rich Taylor's softcore LJ // // V_{LJ}(r) = 4 epsilon [ (sigma^12 / (alpha^m sigma^6 + r^6)^2) - // (sigma^6 / (alpha^m sigma^6 + r^6) ) ] // - // half_sigma and two_sqrt_epsilon are supplied to save cycles. + // V_{coul}(r) = q_i q_j / 4 pi eps_0 (delta+r^2)^(1/2) + // + // delta = shift_coulomb^2 * alpha + // + // Note that we supply half_sigma and two_sqrt_epsilon to save some + // cycles + // + // Note also that we subtract the normal coulomb energy as this + // is calculated during the standard NonbondedForce + // + // 138.9354558466661 is the constant needed to get energies in + // kJ mol-1 given the units of charge (|e|) and distance (nm) // - lj_expression = QString("two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" - "sig6=(sigma^6)/(%1*sigma^6 + r_safe^6);" - "r_safe=max(r, 0.001);" - "max_alpha=max(alpha1, alpha2);" - "sigma=half_sigma1+half_sigma2;") - .arg(taylor_power_expression("max_alpha", taylor_power)) - .toStdString(); + clj_expression = QString("coul_nrg+lj_nrg;" + "coul_nrg=138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + "lj_nrg=two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" + "sig6=(sigma^6)/(%2*sigma^6 + r_safe^6);" + "r_safe=max(r, 0.001);" + "max_kappa=max(kappa1, kappa2);" + "max_alpha=max(alpha1, alpha2);" + "sigma=half_sigma1+half_sigma2;") + .arg(shift_coulomb) + .arg(taylor_power_expression("max_alpha", taylor_power)) + .toStdString(); } else { - // Zacharias and McCammon, J. Chem. Phys., 1994, and also, - // Michel et al., JCTC, 2007 + // this uses the following potentials + // Zacharias and McCammon, J. Chem. Phys., 1994, and also, + // Michel et al., JCTC, 2007 // // V_{LJ}(r) = 4 epsilon [ ( sigma^12 / (delta*sigma + r^2)^6 ) - // ( sigma^6 / (delta*sigma + r^2)^3 ) ] // // delta = shift_delta * alpha // - // half_sigma and two_sqrt_epsilon are supplied to save cycles. + // V_{coul}(r) = q_i q_j / 4 pi eps_0 (delta+r^2)^(1/2) + // + // delta = shift_coulomb^2 * alpha + // + // Note that we pre-calculate delta as a forcefield parameter, + // and also supply half_sigma and two_sqrt_epsilon to save some + // cycles // - lj_expression = QString("two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" - "sig6=(sigma^6)/(((sigma*delta) + r_safe^2)^3);" - "delta=%1*max_alpha;" - "r_safe=max(r, 0.001);" - "max_alpha=max(alpha1, alpha2);" - "sigma=half_sigma1+half_sigma2;") - .arg(shift_delta.to(SireUnits::nanometer)) - .toStdString(); - } - - ghost_ghost_coulombff = new OpenMM::CustomNonbondedForce(coul_expression); - ghost_ghost_coulombff->setName("GhostGhostCoulombForce"); - ghost_ghost_ljff = new OpenMM::CustomNonbondedForce(lj_expression); - ghost_ghost_ljff->setName("GhostGhostLJForce"); - ghost_nonghost_coulombff = new OpenMM::CustomNonbondedForce(coul_expression); - ghost_nonghost_coulombff->setName("GhostNonGhostCoulombForce"); - ghost_nonghost_ljff = new OpenMM::CustomNonbondedForce(lj_expression); - ghost_nonghost_ljff->setName("GhostNonGhostLJForce"); - - for (auto ff : {ghost_ghost_coulombff, ghost_nonghost_coulombff}) - { - ff->addPerParticleParameter("q"); - ff->addPerParticleParameter("alpha"); - ff->addPerParticleParameter("kappa"); - ff->addGlobalParameter("lambda", 0.00000); - } - - for (auto ff : {ghost_ghost_ljff, ghost_nonghost_ljff}) - { - ff->addPerParticleParameter("half_sigma"); - ff->addPerParticleParameter("two_sqrt_epsilon"); - ff->addPerParticleParameter("alpha"); - ff->addGlobalParameter("lambda", 0.00000); - } - - // Coulomb forces must not use LRC: the soft-core 1/sqrt(alpha+r^2)-1/r - // expression decays as 1/r^3 when alpha>0, making the LRC integral - // (proportional to int r^2 * 1/r^3 dr = int 1/r dr) log-divergent. - ghost_ghost_coulombff->setUseLongRangeCorrection(false); - ghost_nonghost_coulombff->setUseLongRangeCorrection(false); - - // LJ soft-core LRC is handled analytically via a CustomVolumeForce - // rather than OpenMM's numerical integrator, because the standard - // LJ tail (r > rc, soft-core shift negligible) has a closed-form - // expression and this allows the result to be cached per lambda state. - ghost_ghost_ljff->setUseLongRangeCorrection(false); - ghost_nonghost_ljff->setUseLongRangeCorrection(false); + // Note also that we subtract the normal coulomb energy as this + // is calculated during the standard NonbondedForce + // + // 138.9354558466661 is the constant needed to get energies in + // kJ mol-1 given the units of charge (|e|) and distance (nm) + // + clj_expression = QString("coul_nrg+lj_nrg;" + "coul_nrg=138.9354558466661*q1*q2*((1/sqrt((%1*max_alpha)+r_safe^2))-(max_kappa/r_safe));" + "lj_nrg=two_sqrt_epsilon1*two_sqrt_epsilon2*sig6*(sig6-1);" + "sig6=(sigma^6)/(((sigma*delta) + r_safe^2)^3);" + "delta=%2*max_alpha;" + "r_safe=max(r, 0.001);" + "max_kappa=max(kappa1, kappa2);" + "max_alpha=max(alpha1, alpha2);" + "sigma=half_sigma1+half_sigma2;") + .arg(shift_coulomb) + .arg(shift_delta.to(SireUnits::nanometer)) + .toStdString(); + } + + ghost_ghostff = new OpenMM::CustomNonbondedForce(clj_expression); + ghost_ghostff->setName("GhostGhostNonbondedForce"); + ghost_nonghostff = new OpenMM::CustomNonbondedForce(clj_expression); + ghost_nonghostff->setName("GhostNonGhostNonbondedForce"); + + ghost_ghostff->addPerParticleParameter("q"); + ghost_ghostff->addPerParticleParameter("half_sigma"); + ghost_ghostff->addPerParticleParameter("two_sqrt_epsilon"); + ghost_ghostff->addPerParticleParameter("alpha"); + ghost_ghostff->addPerParticleParameter("kappa"); + + ghost_nonghostff->addPerParticleParameter("q"); + ghost_nonghostff->addPerParticleParameter("half_sigma"); + ghost_nonghostff->addPerParticleParameter("two_sqrt_epsilon"); + ghost_nonghostff->addPerParticleParameter("alpha"); + ghost_nonghostff->addPerParticleParameter("kappa"); + + // LRC for the ghost soft-core is handled analytically via a CustomVolumeForce + // (Coulomb has no well-defined LRC; LJ tail is handled by the ghost-lrc force). + ghost_ghostff->setUseLongRangeCorrection(false); + ghost_nonghostff->setUseLongRangeCorrection(false); if (ffinfo.hasCutoff()) { if (ffinfo.space().isPeriodic()) { - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffPeriodic); + ghost_ghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffPeriodic); + ghost_nonghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffPeriodic); } else { - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffNonPeriodic); + ghost_ghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffNonPeriodic); + ghost_nonghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::CutoffNonPeriodic); } - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->setCutoffDistance(ffinfo.cutoff().to(SireUnits::nanometers)); + ghost_ghostff->setCutoffDistance(ffinfo.cutoff().to(SireUnits::nanometers)); + ghost_nonghostff->setCutoffDistance(ffinfo.cutoff().to(SireUnits::nanometers)); } else { - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->setNonbondedMethod(OpenMM::CustomNonbondedForce::NoCutoff); + ghost_ghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::NoCutoff); + ghost_nonghostff->setNonbondedMethod(OpenMM::CustomNonbondedForce::NoCutoff); } - ghost_ghost_coulombff->setForceGroup(force_group_counter); - lambda_lever.setForceIndex("ghost/ghost-coulomb", system.addForce(ghost_ghost_coulombff)); - lambda_lever.setForceGroup("ghost/ghost-coulomb", force_group_counter++); - - ghost_ghost_ljff->setForceGroup(force_group_counter); - lambda_lever.setForceIndex("ghost/ghost-lj", system.addForce(ghost_ghost_ljff)); - lambda_lever.setForceGroup("ghost/ghost-lj", force_group_counter++); + ghost_ghostff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ghost/ghost", system.addForce(ghost_ghostff)); + lambda_lever.setForceGroup("ghost/ghost", force_group_counter++); - ghost_nonghost_coulombff->setForceGroup(force_group_counter); - lambda_lever.setForceIndex("ghost/non-ghost-coulomb", system.addForce(ghost_nonghost_coulombff)); - lambda_lever.setForceGroup("ghost/non-ghost-coulomb", force_group_counter++); - - ghost_nonghost_ljff->setForceGroup(force_group_counter); - lambda_lever.setForceIndex("ghost/non-ghost-lj", system.addForce(ghost_nonghost_ljff)); - lambda_lever.setForceGroup("ghost/non-ghost-lj", force_group_counter++); + ghost_nonghostff->setForceGroup(force_group_counter); + lambda_lever.setForceIndex("ghost/non-ghost", system.addForce(ghost_nonghostff)); + lambda_lever.setForceGroup("ghost/non-ghost", force_group_counter++); // Analytic LJ LRC: E = lrc_coeff / V, updated each lambda step via // a cached closed-form sum over interaction-group pairs. @@ -1704,10 +1687,9 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, QHash idx_to_pert_idx; QHash idx_to_qm_idx; - // holders for all of custom parameters for the ghost forces (prevents us - // having to continually re-allocate it) - std::vector custom_params_coul = {0.0, 0.0, 0.0}; - std::vector custom_params_lj = {0.0, 0.0, 0.0}; + // just a holder for all of the custom parameters for the + // ghost forces (prevents us having to continually re-allocate it) + std::vector custom_params = {0.0, 0.0, 0.0, 0.0, 0.0}; // the sets of particle indexes for the ghost atoms and non-ghost atoms QVector ghost_atoms; @@ -1788,13 +1770,11 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // add a perturbable molecule, recording the start index // for each of the forcefields - start_indicies.reserve(9); + start_indicies.reserve(7); start_indicies.insert("clj", start_index); - start_indicies.insert("ghost/ghost-coulomb", start_index); - start_indicies.insert("ghost/ghost-lj", start_index); - start_indicies.insert("ghost/non-ghost-coulomb", start_index); - start_indicies.insert("ghost/non-ghost-lj", start_index); + start_indicies.insert("ghost/ghost", start_index); + start_indicies.insert("ghost/non-ghost", start_index); // the start index for this molecules first bond, angle or // torsion parameters will be however many of these @@ -1873,24 +1853,19 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, } // reduced_q - custom_params_coul[0] = charge; - // alpha - custom_params_coul[1] = alphas_data[j]; - // kappa - custom_params_coul[2] = kappas_data[j]; - + custom_params[0] = charge; // half_sigma - custom_params_lj[0] = 0.5 * boost::get<1>(clj); + custom_params[1] = 0.5 * boost::get<1>(clj); // two_sqrt_epsilon - custom_params_lj[1] = 2.0 * std::sqrt(boost::get<2>(clj)); + custom_params[2] = 2.0 * std::sqrt(boost::get<2>(clj)); // alpha - custom_params_lj[2] = alphas_data[j]; + custom_params[3] = alphas_data[j]; + // kappa + custom_params[4] = kappas_data[j]; // Add the particle to the ghost and nonghost forcefields - ghost_ghost_coulombff->addParticle(custom_params_coul); - ghost_ghost_ljff->addParticle(custom_params_lj); - ghost_nonghost_coulombff->addParticle(custom_params_coul); - ghost_nonghost_ljff->addParticle(custom_params_lj); + ghost_ghostff->addParticle(custom_params); + ghost_nonghostff->addParticle(custom_params); real_atoms.append(atom_index); @@ -1959,24 +1934,18 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // forcefields if there are any perturbable molecules if (any_perturbable) { - // reduced_q - custom_params_coul[0] = boost::get<0>(clj); - // alpha - is zero for non-ghost atoms - custom_params_coul[1] = 0.0; - // kappa - is 0 for non-ghost atoms - custom_params_coul[2] = 0.0; - + // reduced charge + custom_params[0] = boost::get<0>(clj); // half_sigma - custom_params_lj[0] = 0.5 * boost::get<1>(clj); + custom_params[1] = 0.5 * boost::get<1>(clj); // two_sqrt_epsilon - custom_params_lj[1] = 2.0 * std::sqrt(boost::get<2>(clj)); + custom_params[2] = 2.0 * std::sqrt(boost::get<2>(clj)); // alpha - is zero for non-ghost atoms - custom_params_lj[2] = 0.0; - - ghost_ghost_coulombff->addParticle(custom_params_coul); - ghost_ghost_ljff->addParticle(custom_params_lj); - ghost_nonghost_coulombff->addParticle(custom_params_coul); - ghost_nonghost_ljff->addParticle(custom_params_lj); + custom_params[3] = 0.0; + // kappa - is 0 for non-ghost atoms + custom_params[4] = 0.0; + ghost_ghostff->addParticle(custom_params); + ghost_nonghostff->addParticle(custom_params); non_ghost_atoms.append(atom_index); } } @@ -2038,23 +2007,18 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, if (any_perturbable and mol.isPerturbable()) { // reduced_q - custom_params_coul[0] = vs_charge; - // alpha - custom_params_coul[1] = alphas_data[mol.molinfo.nAtoms() + k]; - // kappa - custom_params_coul[2] = kappas_data[mol.molinfo.nAtoms() + k]; - + custom_params[0] = vs_charge; // half_sigma - custom_params_lj[0] = 1.0; + custom_params[1] = 1.0; // two_sqrt_epsilon - custom_params_lj[1] = 0.0; + custom_params[2] = 0.0; // alpha - custom_params_lj[2] = alphas_data[mol.molinfo.nAtoms() + k]; + custom_params[3] = alphas_data[mol.molinfo.nAtoms() + k]; + // kappa + custom_params[4] = kappas_data[mol.molinfo.nAtoms() + k]; - ghost_ghost_coulombff->addParticle(custom_params_coul); - ghost_ghost_ljff->addParticle(custom_params_lj); - ghost_nonghost_coulombff->addParticle(custom_params_coul); - ghost_nonghost_ljff->addParticle(custom_params_lj); + ghost_ghostff->addParticle(custom_params); + ghost_nonghostff->addParticle(custom_params); const bool vs_to_ghost = mol.to_ghost_idxs.contains(parent_idx); const bool vs_from_ghost = mol.from_ghost_idxs.contains(parent_idx); @@ -2075,12 +2039,9 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, else if (any_perturbable) { // Add to ghost FFs if necessary - custom_params_coul = {vs_charge, 0.0, 0.0}; - custom_params_lj = {1.0, 0.0, 0.0}; - ghost_ghost_coulombff->addParticle(custom_params_coul); - ghost_ghost_ljff->addParticle(custom_params_lj); - ghost_nonghost_coulombff->addParticle(custom_params_coul); - ghost_nonghost_ljff->addParticle(custom_params_lj); + custom_params = {vs_charge, 1.0, 0.0, 0.0, 0.0}; + ghost_ghostff->addParticle(custom_params); + ghost_nonghostff->addParticle(custom_params); non_ghost_atoms.append(atom_index); } } @@ -2256,16 +2217,14 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, /// Finally tell the ghost forcefields about the ghost and non-ghost /// interaction groups, so that they can correctly calculate the /// ghost/ghost and ghost/non-ghost energies - if (ghost_ghost_ljff != 0 and ghost_nonghost_ljff != 0) + if (ghost_ghostff != 0 and ghost_nonghostff != 0) { // set up the interaction groups - ghost / non-ghost // ghost / ghost std::set ghost_atoms_set(ghost_atoms.begin(), ghost_atoms.end()); std::set non_ghost_atoms_set(non_ghost_atoms.begin(), non_ghost_atoms.end()); - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff}) - ff->addInteractionGroup(ghost_atoms_set, ghost_atoms_set); - for (auto ff : {ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->addInteractionGroup(ghost_atoms_set, non_ghost_atoms_set); + ghost_ghostff->addInteractionGroup(ghost_atoms_set, ghost_atoms_set); + ghost_nonghostff->addInteractionGroup(ghost_atoms_set, non_ghost_atoms_set); } // Register GCMC water atom indices with the lambda lever and pre-compute the @@ -2299,7 +2258,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, water_atom_indices.append(j); } - lambda_lever.setWaterAtoms(water_atom_indices); + lambda_lever.setGCMCWaterAtoms(water_atom_indices); // Pre-compute lrc_w_solute (per active water molecule, interaction with // all non-water atoms: protein, ligand, ions) and lrc_ww_half (per active @@ -2567,11 +2526,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // we need to make sure that the list of exclusions in // the NonbondedForce match those in the CustomNonbondedForces - if (ghost_ghost_ljff != 0) + if (ghost_ghostff != 0) { - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->addExclusion(boost::get<0>(p), boost::get<1>(p)); + ghost_ghostff->addExclusion(boost::get<0>(p), boost::get<1>(p)); + ghost_nonghostff->addExclusion(boost::get<0>(p), boost::get<1>(p)); } } @@ -2599,11 +2557,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, cljff->addException(vs0_index, start_index + a, 0.0, 1, 0, false); - if (ghost_ghost_ljff != 0) + if (ghost_ghostff != 0) { - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->addExclusion(vs0_index, start_index + a); + ghost_ghostff->addExclusion(vs0_index, start_index + a); + ghost_nonghostff->addExclusion(vs0_index, start_index + a); } for (int v1 = v0 + 1; v1 < atom_vs.size(); ++v1) @@ -2612,11 +2569,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, cljff->addException(vs0_index, vs1_index, 0.0, 1, 0, false); - if (ghost_ghost_ljff != 0) + if (ghost_ghostff != 0) { - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->addExclusion(vs0_index, vs1_index); + ghost_ghostff->addExclusion(vs0_index, vs1_index); + ghost_nonghostff->addExclusion(vs0_index, vs1_index); } } } @@ -2646,9 +2602,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // and if the two atoms are in the same molecule if (mol_from == mol_to) { - for (auto ff : {ghost_ghost_coulombff, ghost_ghost_ljff, - ghost_nonghost_coulombff, ghost_nonghost_ljff}) - ff->addExclusion(from_ghost_idx, to_ghost_idx); + ghost_ghostff->addExclusion(from_ghost_idx, to_ghost_idx); + ghost_nonghostff->addExclusion(from_ghost_idx, to_ghost_idx); cljff->addException(from_ghost_idx, to_ghost_idx, 0.0, 1e-9, 1e-9, false); } From 22c36d6d2855fea03e8f9ebf5fa8b391c30c5f4a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 6 May 2026 15:54:10 +0100 Subject: [PATCH 18/21] Add unit test for GCMC water LRC. --- tests/convert/test_openmm_lrc.py | 73 ++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/tests/convert/test_openmm_lrc.py b/tests/convert/test_openmm_lrc.py index ffb7f23b2..f575b353c 100644 --- a/tests/convert/test_openmm_lrc.py +++ b/tests/convert/test_openmm_lrc.py @@ -5,6 +5,10 @@ Tests use merged_ethane_methanol (merged_molecule.s3) with the 5 nearest waters so they run quickly while still exercising a periodic cutoff system. + +A separate test exercises the GCMC water LRC (GCMCLRCForce) using the ala_mols +system (alanine-dipeptide in water). GCMC is faked by decrementing n_w by hand +and verifying the energy change against the analytical formula. """ import pytest @@ -100,3 +104,72 @@ def test_lrc_lambda1_matches_perturbed_end_state( nrg_pert_end = omm_pert_end.get_energy().value() assert nrg_pert == pytest.approx(nrg_pert_end, rel=1e-3) + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_gcmc_lrc(ala_mols, openmm_platform): + """ + Test the GCMCLRCForce by faking a GCMC move: decrement n_w by 1 and verify + the energy change matches the analytical formula. + + E(n_w) = (n_w * lrc_w_solute + n_w * (n_w - 1) * lrc_ww_half) / V + + dE = E(n_w - 1) - E(n_w) = -(lrc_w_solute + 2 * (n_w - 1) * lrc_ww_half) / V + """ + # Build with GCMC LRC enabled, reserving 1 buffer water so n_w starts at + # n_waters - 1. + d = ala_mols.dynamics( + platform=openmm_platform, + map={ + "use_dispersion_correction": True, + "use_gcmc_lrc": True, + "num_gcmc_waters": 1, + }, + ) + d.set_lambda(0.0) + + context = d.context() + omm_system = context.getSystem() + + # Locate the GCMCLRCForce and its force group. + gcmc_force = None + for force in omm_system.getForces(): + if force.getName() == "GCMCLRCForce": + gcmc_force = force + break + assert gcmc_force is not None, "GCMCLRCForce not found in system" + fg = gcmc_force.getForceGroup() + + # Read the pre-computed LRC coefficients and initial n_w from the context. + state = context.getState(getParameters=True, getPositions=True) + params = state.getParameters() + n_w = params["n_w"] + lrc_w_solute = params["lrc_w_solute"] + lrc_ww_half = params["lrc_ww_half"] + + # The LJ tail is attractive so both coefficients must be negative. + assert lrc_w_solute < 0 + assert lrc_ww_half < 0 + assert n_w > 0 + + # Get the GCMC LRC energy at the initial n_w. + nrg1 = ( + context.getState(getEnergy=True, groups=(1 << fg)).getPotentialEnergy()._value + ) + + # "Turn off" one water by decrementing n_w. + context.setParameter("n_w", n_w - 1) + nrg2 = ( + context.getState(getEnergy=True, groups=(1 << fg)).getPotentialEnergy()._value + ) + + # Compute the expected energy change analytically. + # Box vectors are in nm; volume in nm^3 matches the lrc_coeff units (kJ/mol·nm^3). + box = state.getPeriodicBoxVectors() + V = box[0][0]._value * box[1][1]._value * box[2][2]._value + + expected_delta = -(lrc_w_solute + 2 * (n_w - 1) * lrc_ww_half) / V + assert nrg2 - nrg1 == pytest.approx(expected_delta, rel=1e-5) From 298f29bcb5c005e087d9bd21f105c4a02ea95653 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 6 May 2026 16:20:30 +0100 Subject: [PATCH 19/21] Improve formatting. --- wrapper/Convert/SireOpenMM/lambdalever.cpp | 36 +++++++++++++++++++--- 1 file changed, 31 insertions(+), 5 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp index f30799558..1f8023ee1 100644 --- a/wrapper/Convert/SireOpenMM/lambdalever.cpp +++ b/wrapper/Convert/SireOpenMM/lambdalever.cpp @@ -1495,11 +1495,37 @@ double LambdaLever::setLambda(OpenMM::Context &context, const int nparams = morphed_charges.count(); // Detect whether any CLJ or ghost-14 parameters changed - has_changed_cljff |= rest2_changed || cache.hasChanged("clj", "charge") || cache.hasChanged("clj", "sigma") || cache.hasChanged("clj", "epsilon") || cache.hasChanged("clj", "alpha") || cache.hasChanged("clj", "kappa") || cache.hasChanged("clj", "charge_scale") || cache.hasChanged("clj", "lj_scale"); - - has_changed_ghostff |= rest2_changed || cache.hasChanged("ghost/ghost", "charge") || cache.hasChanged("ghost/ghost", "sigma") || cache.hasChanged("ghost/ghost", "epsilon") || cache.hasChanged("ghost/ghost", "alpha") || cache.hasChanged("ghost/ghost", "kappa") || cache.hasChanged("ghost/non-ghost", "charge") || cache.hasChanged("ghost/non-ghost", "sigma") || cache.hasChanged("ghost/non-ghost", "epsilon") || cache.hasChanged("ghost/non-ghost", "alpha") || cache.hasChanged("ghost/non-ghost", "kappa"); - - has_changed_ghost14ff |= rest2_changed || cache.hasChanged("ghost-14", "charge") || cache.hasChanged("ghost-14", "sigma") || cache.hasChanged("ghost-14", "epsilon") || cache.hasChanged("ghost-14", "alpha") || cache.hasChanged("ghost-14", "kappa") || cache.hasChanged("ghost-14", "charge_scale") || cache.hasChanged("ghost-14", "lj_scale"); + // clang-format off + has_changed_cljff |= rest2_changed + || cache.hasChanged("clj", "charge") + || cache.hasChanged("clj", "sigma") + || cache.hasChanged("clj", "epsilon") + || cache.hasChanged("clj", "alpha") + || cache.hasChanged("clj", "kappa") + || cache.hasChanged("clj", "charge_scale") + || cache.hasChanged("clj", "lj_scale"); + + has_changed_ghostff |= rest2_changed + || cache.hasChanged("ghost/ghost", "charge") + || cache.hasChanged("ghost/ghost", "sigma") + || cache.hasChanged("ghost/ghost", "epsilon") + || cache.hasChanged("ghost/ghost", "alpha") + || cache.hasChanged("ghost/ghost", "kappa") + || cache.hasChanged("ghost/non-ghost", "charge") + || cache.hasChanged("ghost/non-ghost", "sigma") + || cache.hasChanged("ghost/non-ghost", "epsilon") + || cache.hasChanged("ghost/non-ghost", "alpha") + || cache.hasChanged("ghost/non-ghost", "kappa"); + + has_changed_ghost14ff |= rest2_changed + || cache.hasChanged("ghost-14", "charge") + || cache.hasChanged("ghost-14", "sigma") + || cache.hasChanged("ghost-14", "epsilon") + || cache.hasChanged("ghost-14", "alpha") + || cache.hasChanged("ghost-14", "kappa") + || cache.hasChanged("ghost-14", "charge_scale") + || cache.hasChanged("ghost-14", "lj_scale"); + // clang-format on if (have_ghost_atoms) { From 0a43fe704ab002912168964de04e8609348c5907 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 May 2026 09:07:49 +0100 Subject: [PATCH 20/21] Add method for reversing a LambdaSchedule. [ref OpenBioSim/somd2#78] --- corelib/src/libs/SireCAS/lambdaschedule.cpp | 53 ++++++++++++++ corelib/src/libs/SireCAS/lambdaschedule.h | 2 + doc/source/changelog.rst | 7 ++ tests/cas/test_lambdaschedule.py | 80 +++++++++++++++++++++ wrapper/CAS/LambdaSchedule.pypp.cpp | 8 +++ 5 files changed, 150 insertions(+) diff --git a/corelib/src/libs/SireCAS/lambdaschedule.cpp b/corelib/src/libs/SireCAS/lambdaschedule.cpp index 8464ff085..bcca60fa1 100644 --- a/corelib/src/libs/SireCAS/lambdaschedule.cpp +++ b/corelib/src/libs/SireCAS/lambdaschedule.cpp @@ -28,6 +28,8 @@ #include "lambdaschedule.h" +#include "SireCAS/identities.h" +#include "SireCAS/symbolexpression.h" #include "SireCAS/values.h" #include "SireBase/console.h" @@ -1670,3 +1672,54 @@ QVector LambdaSchedule::morph(const QString &force, return morphed; } + +/** Return a new LambdaSchedule that is the reverse of this schedule. + * The stages are reversed in order, and within each stage's equations + * the lambda symbol is replaced by (1 - lambda) and the initial and + * final symbols are swapped simultaneously. This flips the schedule + * about its midpoint so that the end state becomes the start state + * and vice versa. + * + * The invariant is: + * reversed.morph(force, lever, initial, final, λ) + * == original.morph(force, lever, final, initial, 1-λ) + */ +LambdaSchedule LambdaSchedule::reverse() const +{ + if (this->isNull()) + return *this; + + // Simultaneous substitution: λ → (1-λ), initial ↔ final + const Identities ids( + SymbolExpression(lambda_symbol, 1.0 - Expression(lambda_symbol)), + SymbolExpression(initial_symbol, Expression(final_symbol)), + SymbolExpression(final_symbol, Expression(initial_symbol))); + + auto transform = [&](const Expression &e) -> Expression + { + auto r = e.substitute(ids); + if (r == default_morph_equation) + r = default_morph_equation; + return r; + }; + + LambdaSchedule result(*this); + + std::reverse(result.stage_names.begin(), result.stage_names.end()); + std::reverse(result.default_equations.begin(), result.default_equations.end()); + std::reverse(result.stage_equations.begin(), result.stage_equations.end()); + std::reverse(result.stage_weights.begin(), result.stage_weights.end()); + + for (int i = 0; i < result.nStages(); ++i) + { + result.default_equations[i] = transform(result.default_equations[i]); + + for (auto &eq : result.stage_equations[i]) + eq = transform(eq); + } + + for (auto &mol_sched : result.mol_schedules) + mol_sched = mol_sched.reverse(); + + return result; +} diff --git a/corelib/src/libs/SireCAS/lambdaschedule.h b/corelib/src/libs/SireCAS/lambdaschedule.h index b2ff13222..5ee79ecef 100644 --- a/corelib/src/libs/SireCAS/lambdaschedule.h +++ b/corelib/src/libs/SireCAS/lambdaschedule.h @@ -233,6 +233,8 @@ namespace SireCAS double clamp(double lambda_value) const; + LambdaSchedule reverse() const; + protected: int find_stage(const QString &stage) const; diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index f43950e00..69834fabf 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -77,6 +77,13 @@ organisation on `GitHub `__. and ghost–non-ghost interactions. Both are updated efficiently via the lambda lever without recomputing the full dispersion correction on every λ change. +* Added a ``reverse()`` method to :class:`~sire.cas.LambdaSchedule` that returns a new + schedule flipped about its midpoint. Stages are reversed in order and each equation is + transformed by substituting ``λ → (1-λ)`` and swapping ``initial`` ↔ ``final`` + simultaneously. The invariant is that + ``reversed.morph(force, lever, initial, final, λ)`` equals + ``original.morph(force, lever, final, initial, 1-λ)``. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/tests/cas/test_lambdaschedule.py b/tests/cas/test_lambdaschedule.py index add7ee969..1f411d222 100644 --- a/tests/cas/test_lambdaschedule.py +++ b/tests/cas/test_lambdaschedule.py @@ -547,3 +547,83 @@ def test_stage_weight_get_stage_boundaries(): # Midpoint of stage "b": lambda=0.25 + 1.5/4 = 0.625 → local lam=0.5 assert l.get_lambda_in_stage(0.625) == pytest.approx(0.5) + + +def test_reverse_stage_order(): + """Reversing reverses the order of stage names.""" + l = sr.cas.LambdaSchedule() + l.add_stage("a", l.lam()) + l.add_stage("b", l.lam()) + l.add_stage("c", l.lam()) + + r = l.reverse() + assert r.get_stages() == ["c", "b", "a"] + + +def test_reverse_stage_weights(): + """Reversing preserves stage weights in reversed order.""" + l = sr.cas.LambdaSchedule() + l.add_stage("a", l.lam(), weight=1.0) + l.add_stage("b", l.lam(), weight=2.0) + l.add_stage("c", l.lam(), weight=3.0) + + r = l.reverse() + assert r.get_stage_weights() == pytest.approx([3.0, 2.0, 1.0]) + + +def test_reverse_invariant(): + """reversed.morph(init, fin, λ) == original.morph(fin, init, 1-λ).""" + import random + + random.seed(42) + + for schedule in [ + sr.cas.LambdaSchedule.standard_morph(), + sr.cas.LambdaSchedule.charge_scaled_morph(), + sr.cas.LambdaSchedule.standard_decouple(), + ]: + r = schedule.reverse() + for _ in range(50): + lam_val = random.uniform(0.0, 1.0) + assert r.morph( + initial=0.0, final=1.0, lambda_value=lam_val + ) == pytest.approx( + schedule.morph(initial=1.0, final=0.0, lambda_value=1.0 - lam_val), + abs=1e-10, + ) + + +def test_reverse_standard_morph_is_identity(): + """Reversing a standard morph schedule gives the same behaviour (it is symmetric).""" + import random + + random.seed(42) + l = sr.cas.LambdaSchedule.standard_morph() + r = l.reverse() + + for _ in range(50): + lam_val = random.uniform(0.0, 1.0) + assert l.morph(initial=0.0, final=1.0, lambda_value=lam_val) == pytest.approx( + r.morph(initial=0.0, final=1.0, lambda_value=lam_val), abs=1e-10 + ) + + +def test_reverse_twice_is_identity(): + """Reversing a schedule twice recovers the original behaviour.""" + import random + + random.seed(42) + + for schedule in [ + sr.cas.LambdaSchedule.standard_morph(), + sr.cas.LambdaSchedule.charge_scaled_morph(), + sr.cas.LambdaSchedule.standard_decouple(), + ]: + rr = schedule.reverse().reverse() + for _ in range(50): + lam_val = random.uniform(0.0, 1.0) + assert schedule.morph( + initial=0.0, final=1.0, lambda_value=lam_val + ) == pytest.approx( + rr.morph(initial=0.0, final=1.0, lambda_value=lam_val), abs=1e-10 + ) diff --git a/wrapper/CAS/LambdaSchedule.pypp.cpp b/wrapper/CAS/LambdaSchedule.pypp.cpp index 011114f5e..fc64963de 100644 --- a/wrapper/CAS/LambdaSchedule.pypp.cpp +++ b/wrapper/CAS/LambdaSchedule.pypp.cpp @@ -520,6 +520,14 @@ void register_LambdaSchedule_class() LambdaSchedule_exposer.def( "removeStage", removeStage_function_value, (bp::arg("stage")), bp::release_gil_policy(), "Remove the stage stage"); } + { //::SireCAS::LambdaSchedule::reverse + + typedef ::SireCAS::LambdaSchedule (::SireCAS::LambdaSchedule::*reverse_function_type)() const; + reverse_function_type reverse_function_value(&::SireCAS::LambdaSchedule::reverse); + + LambdaSchedule_exposer.def( + "reverse", reverse_function_value, bp::release_gil_policy(), "Return a new LambdaSchedule that is the reverse of this schedule.\n The stages are reversed in order and in each stages equations\n the lambda symbol is replaced by (1 - lambda) and the initial and\n final symbols are swapped simultaneously. The invariant is:\n reversed.morph(force, lever, initial, final, lam) ==\n original.morph(force, lever, final, initial, 1-lam)\n"); + } { //::SireCAS::LambdaSchedule::setConstant typedef ::SireCAS::Symbol (::SireCAS::LambdaSchedule::*setConstant_function_type)(::QString const &, double); From d022e0bbe064c0a8ec3303a9001ecc608e5a6b7a Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 May 2026 12:18:28 +0100 Subject: [PATCH 21/21] Autodetect Visual Studio generator version. --- setup.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 5e12e62f0..f10ae77da 100644 --- a/setup.py +++ b/setup.py @@ -271,6 +271,31 @@ def add_default_cmake_defs(cmake_defs, ncores): cmake_defs.append(["SIRE_DISABLE_AVX512F=ON"]) +def _detect_vs_generator(): + # Map the installed VS major version to the CMake generator name. + # VS 2019=16, VS 2022=17, VS 2026=18 (and beyond). + _VS_GENERATORS = { + 18: "Visual Studio 18 2026", + 17: "Visual Studio 17 2022", + 16: "Visual Studio 16 2019", + 15: "Visual Studio 15 2017", + } + try: + result = subprocess.run( + ["vswhere.exe", "-nologo", "-latest", "-property", "installationVersion"], + capture_output=True, + text=True, + ) + if result.returncode == 0 and result.stdout.strip(): + major = int(result.stdout.strip().split(".")[0]) + for version in sorted(_VS_GENERATORS, reverse=True): + if major >= version: + return _VS_GENERATORS[version] + except Exception: + pass + return "Visual Studio 17 2022" + + def make_cmd(ncores, install=False): if is_windows: action = "INSTALL" if install else "ALL_BUILD" @@ -661,7 +686,7 @@ def install(ncores: int = 1, npycores: int = 1): action = args.action[0] if is_windows and (args.generator is None or len(args.generator) == 0): - args.generator = [["Visual Studio 17 2022"]] + args.generator = [[_detect_vs_generator()]] args.architecture = [["x64"]] elif is_macos: # fix compile bug when INSTALL_NAME_TOOL is not set