Add long-range dispersion correction, Beutler soft-core, and LambdaSchedule.reverse()

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<int> 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;
+}

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;
 

doc/source/changelog.rst

@@ -67,6 +67,23 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * 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.
+
+* 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 <https://github.com/openbiosim/sire/compare/2025.3.0...2025.4.0>`__ - February 2026
 ---------------------------------------------------------------------------------------------
 

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

doc/source/tutorial/part07/03_ghosts.rst

@@ -55,8 +55,22 @@ 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).
 
-There are two different soft-core potentials available. The default is
-the Zacharias potential, while the second is the Taylor potential.
+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 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
 -------------------
@@ -68,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]
 
@@ -103,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``
@@ -132,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]
 
@@ -169,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``
@@ -187,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
 -------------
 

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

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,