From 61901fc6430646eb452df1129aefa41462e69d04 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Wed, 29 Apr 2026 14:45:54 +0100
Subject: [PATCH 01/12] Add support for long-range LJ dispersion correction.

---
 src/somd2/config/_config.py | 18 ++++++++++++++++++
 src/somd2/runner/_base.py   | 10 ++++++++++
 2 files changed, 28 insertions(+)

diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py
index fce678e..7dae0aa 100644
--- a/src/somd2/config/_config.py
+++ b/src/somd2/config/_config.py
@@ -151,6 +151,7 @@ def __init__(
         gcmc_radius="4 A",
         gcmc_bulk_sampling_probability=0.1,
         gcmc_tolerance=0.0,
+        use_dispersion_correction=False,
         rest2_scale=1.0,
         rest2_selection=None,
         softcore_form="zacharias",
@@ -438,6 +439,12 @@ def __init__(
             of acceptance for a move. This can be used to exclude low probability candidates
             that can cause instabilities or crashes for the MD engine.
 
+        use_dispersion_correction: bool
+            Whether to use the long-range dispersion correction for LJ interactions.
+            When True, the correction is evaluated analytically via a CustomVolumeForce
+            and cached per lambda state, avoiding expensive recomputation on every
+            lambda change. Default False.
+
         rest2_scale: float, list(float)
             The scaling factor for Replica Exchange with Solute Tempering (REST) simulations.
             This is the factor by which the temperature of the solute is scaled with respect to
@@ -605,6 +612,7 @@ def __init__(
         self.gcmc_radius = gcmc_radius
         self.gcmc_bulk_sampling_probability = gcmc_bulk_sampling_probability
         self.gcmc_tolerance = gcmc_tolerance
+        self.use_dispersion_correction = use_dispersion_correction
         self.rest2_scale = rest2_scale
         self.rest2_selection = rest2_selection
         self.restart = restart
@@ -2286,6 +2294,16 @@ def gcmc_tolerance(self, gcmc_tolerance):
             raise ValueError("'gcmc_tolerance' must be greater than or equal to 0.0")
         self._gcmc_tolerance = gcmc_tolerance
 
+    @property
+    def use_dispersion_correction(self):
+        return self._use_dispersion_correction
+
+    @use_dispersion_correction.setter
+    def use_dispersion_correction(self, use_dispersion_correction):
+        if not isinstance(use_dispersion_correction, bool):
+            raise TypeError("'use_dispersion_correction' must be of type 'bool'")
+        self._use_dispersion_correction = use_dispersion_correction
+
     @property
     def rest2_scale(self):
         return self._rest2_scale
diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index a6c637c..c3748b4 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -192,6 +192,16 @@ def __init__(self, system, config):
             self._config.fix_perturbable_zero_sigmas
         )
 
+        # Long-range dispersion correction.
+        self._config._extra_args["use_dispersion_correction"] = (
+            self._config.use_dispersion_correction
+        )
+
+        # GCMC LRC map options.
+        if self._config.gcmc and self._config.use_dispersion_correction:
+            self._config._extra_args["use_gcmc_lrc"] = True
+            self._config._extra_args["num_gcmc_waters"] = self._config.gcmc_num_waters
+
         # If specified, use the Taylor soft-core form.
         if self._config.softcore_form == "taylor":
             self._config._extra_args["use_taylor_softening"] = True

From d3061d3e8b71334763fc51ef9f9ec578833a9809 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 30 Apr 2026 09:21:36 +0100
Subject: [PATCH 02/12] Warn when use_dispersion_correction=True and there is
 no space.

---
 src/somd2/runner/_base.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index c3748b4..571707e 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -301,6 +301,13 @@ def __init__(self, system, config):
         except:
             self._has_water = False
 
+        # Warn if dispersion correction is requested but can't be applied.
+        if self._config.use_dispersion_correction and not self._has_space:
+            _logger.warning(
+                "'use_dispersion_correction=True' has no effect: the system "
+                "has no periodic space. The option will be ignored."
+            )
+
         # Check the end state constraints.
         self._check_end_state_constraints()
 

From 3c75e5793341ef8305c00574d1b096333ffb04c8 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 30 Apr 2026 10:08:45 +0100
Subject: [PATCH 03/12] Add Beutler soft-core form for ABFE.

---
 src/somd2/_utils/_schedules.py | 261 ++++++++++++++++++++++++++
 src/somd2/config/_config.py    | 332 +++++----------------------------
 src/somd2/runner/_base.py      |  12 +-
 3 files changed, 315 insertions(+), 290 deletions(-)
 create mode 100644 src/somd2/_utils/_schedules.py

diff --git a/src/somd2/_utils/_schedules.py b/src/somd2/_utils/_schedules.py
new file mode 100644
index 0000000..9a22ce5
--- /dev/null
+++ b/src/somd2/_utils/_schedules.py
@@ -0,0 +1,261 @@
+######################################################################
+# SOMD2: GPU accelerated alchemical free-energy engine.
+#
+# Copyright: 2023-2026
+#
+# Authors: The OpenBioSim Team <team@openbiosim.org>
+#
+# SOMD2 is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# SOMD2 is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with SOMD2. If not, see <http://www.gnu.org/licenses/>.
+#####################################################################
+
+__all__ = [
+    "annihilate",
+    "decouple",
+    "ring_break_morph",
+    "reverse_ring_break_morph",
+]
+
+
+def annihilate():
+    """
+    Build the ABFE lambda schedule using decharge → annihilate with constant epsilon.
+
+    Annihilation removes ALL non-bonded interactions (including intramolecular LJ
+    between non-bonded pairs), matching the GROMACS ABFE protocol. Epsilon is held
+    constant at its real-atom value throughout the annihilate stage so that the
+    (1-alpha) prefactor for the Buetler soft-core provides the sole LJ decay pathway,
+    giving true single (1-lambda) scaling consistent with GROMACS Beutler.
+
+    Returns
+    -------
+
+    schedule : sire.legacy.CAS.LambdaSchedule
+        The lambda schedule.
+    """
+    from sire.cas import LambdaSchedule as _LambdaSchedule
+
+    # Start with the standard decouple schedule and modify the stages and
+    # equations as needed. This will be folded into Sire in future, but
+    # we will use this approach for prototyping.
+    s = _LambdaSchedule.standard_decouple()
+
+    s.remove_stage("decouple")
+
+    s.add_stage("decharge", equation=s.initial())
+    s.set_equation(
+        stage="decharge",
+        lever="charge",
+        equation=s.lam() * s.final() + s.initial() * (1 - s.lam()),
+    )
+    s.set_equation(stage="decharge", force="restraint", equation=s.lam() * s.final())
+
+    s.add_stage(
+        "annihilate",
+        equation=(-s.lam() + 1) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(stage="annihilate", lever="charge", equation=s.final())
+    s.set_equation(stage="annihilate", force="restraint", equation=s.final())
+    s.set_equation(stage="annihilate", lever="epsilon", equation=s.initial())
+
+    return s
+
+
+def decouple():
+    """
+    Build the ABFE lambda schedule using decharge → decouple with constant epsilon.
+
+    Decoupling removes only INTERMOLECULAR non-bonded interactions; intramolecular
+    terms are preserved via kappa=0 on ghost/ghost and ghost-14 forces. Epsilon is
+    held constant throughout the decouple stage (see annihilate for rationale).
+
+    Returns
+    -------
+
+    schedule : sire.legacy.CAS.LambdaSchedule
+        The lambda schedule.
+    """
+    from sire.cas import LambdaSchedule as _LambdaSchedule
+
+    # Start with the standard decouple schedule and modify the stages and
+    # equations as needed. This will be folded into Sire in future, but
+    # we will use this approach for prototyping.
+    s = _LambdaSchedule.standard_decouple()
+
+    s.set_equation(stage="decouple", lever="restraint", equation=s.final())
+    s.set_equation(stage="decouple", lever="kappa", force="ghost/ghost", equation=0)
+    s.set_equation(stage="decouple", lever="kappa", force="ghost-14", equation=0)
+    s.set_equation(stage="decouple", lever="charge", equation=s.final())
+    s.set_equation(stage="decouple", lever="epsilon", equation=s.initial())
+
+    s.prepend_stage("decharge", s.initial())
+    s.set_equation(
+        stage="decharge",
+        lever="charge",
+        equation=s.lam() * s.final() + s.initial() * (1 - s.lam()),
+    )
+    s.set_equation(stage="decharge", force="ghost/ghost", equation=s.initial())
+    s.set_equation(stage="decharge", force="ghost-14", equation=s.initial())
+    s.set_equation(
+        stage="decharge", lever="kappa", force="ghost/ghost", equation=-s.lam() + 1
+    )
+    s.set_equation(
+        stage="decharge", lever="kappa", force="ghost-14", equation=-s.lam() + 1
+    )
+    s.set_equation(stage="decharge", lever="restraint", equation=s.initial() * s.lam())
+
+    return s
+
+
+def ring_break_morph():
+    """
+    Build a lambda schedule for ring-breaking perturbations.
+
+    Three stages: potential_swap → restraints_off → morph.
+
+    Returns
+    -------
+
+    schedule : sire.legacy.CAS.LambdaSchedule
+        The lambda schedule.
+    """
+    from sire.cas import LambdaSchedule as _LambdaSchedule
+
+    s = _LambdaSchedule.standard_morph()
+
+    s.prepend_stage("restraints_off", s.initial())
+    s.set_equation(stage="restraints_off", lever="morse_soft", equation=1 - s.lam())
+    s.set_equation(stage="restraints_off", lever="morse_hard", equation=0)
+    s.set_equation(stage="restraints_off", lever="bond_k", equation=s.final())
+    s.set_equation(stage="restraints_off", lever="bond_length", equation=s.final())
+    s.set_equation(
+        stage="restraints_off",
+        lever="angle_k",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(
+        stage="restraints_off",
+        lever="angle_size",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(
+        stage="restraints_off",
+        lever="torsion_k",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(
+        stage="restraints_off",
+        lever="torsion_phase",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+
+    s.prepend_stage("potential_swap", s.initial())
+    s.set_equation(stage="potential_swap", lever="morse_hard", equation=1 - s.lam())
+    s.set_equation(stage="potential_swap", lever="morse_soft", equation=0 + s.lam())
+    s.set_equation(
+        stage="potential_swap",
+        lever="bond_k",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(
+        stage="potential_swap",
+        lever="bond_length",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(stage="potential_swap", lever="angle_k", equation=s.initial())
+    s.set_equation(stage="potential_swap", lever="angle_size", equation=s.initial())
+    s.set_equation(stage="potential_swap", lever="torsion_k", equation=s.initial())
+    s.set_equation(stage="potential_swap", lever="torsion_phase", equation=s.initial())
+
+    s.set_equation(stage="morph", lever="morse_hard", equation=0)
+    s.set_equation(stage="morph", lever="morse_soft", equation=0)
+    s.set_equation(stage="morph", lever="bond_k", equation=s.final())
+    s.set_equation(stage="morph", lever="bond_length", equation=s.final())
+    s.set_equation(stage="morph", lever="angle_k", equation=s.final())
+    s.set_equation(stage="morph", lever="angle_size", equation=s.final())
+    s.set_equation(stage="morph", lever="torsion_k", equation=s.final())
+    s.set_equation(stage="morph", lever="torsion_phase", equation=s.final())
+
+    return s
+
+
+def reverse_ring_break_morph():
+    """
+    Build a lambda schedule for reverse ring-breaking perturbations.
+
+    Three stages: morph → bonded_perturb → potential_swap.
+
+    Returns
+    -------
+
+    schedule : sire.legacy.CAS.LambdaSchedule
+        The lambda schedule.
+    """
+    from sire.cas import LambdaSchedule as _LambdaSchedule
+
+    s = _LambdaSchedule.standard_morph()
+
+    s.set_equation(stage="morph", lever="morse_hard", equation=0)
+    s.set_equation(stage="morph", lever="morse_soft", equation=0)
+    s.set_equation(stage="morph", lever="bond_k", equation=s.initial())
+    s.set_equation(stage="morph", lever="bond_length", equation=s.initial())
+    s.set_equation(stage="morph", lever="angle_k", equation=s.initial())
+    s.set_equation(stage="morph", lever="angle_size", equation=s.initial())
+    s.set_equation(stage="morph", lever="torsion_k", equation=s.initial())
+    s.set_equation(stage="morph", lever="torsion_phase", equation=s.initial())
+
+    s.append_stage("bonded_perturb", s.final())
+    s.set_equation(stage="bonded_perturb", lever="morse_soft", equation=0 + s.lam())
+    s.set_equation(stage="bonded_perturb", lever="morse_hard", equation=0)
+    s.set_equation(stage="bonded_perturb", lever="bond_k", equation=s.initial())
+    s.set_equation(stage="bonded_perturb", lever="bond_length", equation=s.initial())
+    s.set_equation(
+        stage="bonded_perturb",
+        lever="angle_k",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(
+        stage="bonded_perturb",
+        lever="angle_size",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(
+        stage="bonded_perturb",
+        lever="torsion_k",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(
+        stage="bonded_perturb",
+        lever="torsion_phase",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+
+    s.append_stage("potential_swap", s.final())
+    s.set_equation(stage="potential_swap", lever="morse_hard", equation=0 + s.lam())
+    s.set_equation(stage="potential_swap", lever="morse_soft", equation=1 - s.lam())
+    s.set_equation(
+        stage="potential_swap",
+        lever="bond_k",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(
+        stage="potential_swap",
+        lever="bond_length",
+        equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
+    )
+    s.set_equation(stage="potential_swap", lever="angle_k", equation=s.final())
+    s.set_equation(stage="potential_swap", lever="angle_size", equation=s.final())
+    s.set_equation(stage="potential_swap", lever="torsion_k", equation=s.final())
+    s.set_equation(stage="potential_swap", lever="torsion_phase", equation=s.final())
+
+    return s
diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py
index 7dae0aa..6d20615 100644
--- a/src/somd2/config/_config.py
+++ b/src/somd2/config/_config.py
@@ -70,9 +70,11 @@ class Config:
             "charge_scaled_morph",
             "ring_break_morph",
             "reverse_ring_break_morph",
+            "annihilate",
+            "decouple",
         ],
         "log_level": [level.lower() for level in _logger._core.levels],
-        "softcore_form": ["zacharias", "taylor"],
+        "softcore_form": ["zacharias", "taylor", "beutler"],
     }
 
     # A dictionary of nargs for the various options.
@@ -103,7 +105,6 @@ def __init__(
         lambda_schedule="standard_morph",
         charge_scale_factor=0.2,
         swap_end_states=False,
-        coulomb_power=0.0,
         shift_coulomb="1 A",
         shift_delta="1.5 A",
         restraints=None,
@@ -156,6 +157,7 @@ def __init__(
         rest2_selection=None,
         softcore_form="zacharias",
         taylor_power=1,
+        beutler_alpha=0.5,
         output_directory="output",
         restart=False,
         use_backup=False,
@@ -234,10 +236,6 @@ def __init__(
         swap_end_states: bool
             Whether to swap the end states of the alchemical system.
 
-        coulomb_power : float
-            Power to use for the soft-core Coulomb interaction. This is used
-            to soften the electrostatic interaction.
-
         shift_coulomb : str
             The soft-core shift-coulomb parameter. This is used to soften the
             Coulomb interaction.
@@ -465,14 +463,20 @@ def __init__(
             be applied to protein mutations.
 
         softcore_form: str
-            The soft-core potential form to use for alchemical interactions. This can be
-            either "zacharias" or "taylor". The default is "zacharias".
+            The soft-core potential form to use for alchemical interactions. Valid
+            options are "zacharias" (default), "taylor", and "beutler". The Beutler
+            form is recommended for ABFE calculations.
 
         taylor_power: int
             The power to use for the alpha term in the Taylor soft-core LJ expression,
             i.e. sig6 = sigma^6 / (alpha^m * sigma^6 + r^6). Must be between 0 and 4.
             The default is 1. Only used when softcore_form is "taylor".
 
+        beutler_alpha: float
+            The dimensionless scale factor for the r^6 shift in the Beutler soft-core
+            form. Must be >= 0. The default is 0.5. Only used when softcore_form is
+            "beutler".
+
         output_directory: str
             Path to a directory to store output files.
 
@@ -566,7 +570,6 @@ def __init__(
         self.lambda_schedule = lambda_schedule
         self.charge_scale_factor = charge_scale_factor
         self.swap_end_states = swap_end_states
-        self.coulomb_power = coulomb_power
         self.shift_coulomb = shift_coulomb
         self.shift_delta = shift_delta
         self.restraints = restraints
@@ -619,6 +622,7 @@ def __init__(
         self.use_backup = use_backup
         self.softcore_form = softcore_form
         self.taylor_power = taylor_power
+        self.beutler_alpha = beutler_alpha
         self.somd1_compatibility = somd1_compatibility
         self.pert_file = pert_file
         self.auto_fix_minimise = auto_fix_minimise
@@ -1075,279 +1079,29 @@ def lambda_schedule(self, lambda_schedule):
                     self._lambda_schedule = _LambdaSchedule.charge_scaled_morph(0.2)
                     self._lambda_schedule_name = "charge_scaled_morph"
                 elif lambda_schedule == "ring_break_morph":
-                    self._lambda_schedule = _LambdaSchedule.standard_morph()
-                    self._lambda_schedule.prepend_stage(
-                        "restraints_off", self._lambda_schedule.initial()
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="restraints_off",
-                        lever="morse_soft",
-                        equation=1 - self._lambda_schedule.lam(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="restraints_off", lever="morse_hard", equation=0
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="restraints_off",
-                        lever="bond_k",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="restraints_off",
-                        lever="bond_length",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="restraints_off",
-                        lever="angle_k",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="restraints_off",
-                        lever="angle_size",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="restraints_off",
-                        lever="torsion_k",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="restraints_off",
-                        lever="torsion_phase",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-
-                    self._lambda_schedule.prepend_stage(
-                        "potential_swap", self._lambda_schedule.initial()
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="morse_hard",
-                        equation=1 - self._lambda_schedule.lam(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="morse_soft",
-                        equation=0 + self._lambda_schedule.lam(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="bond_k",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="bond_length",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="angle_k",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="angle_size",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="torsion_k",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="torsion_phase",
-                        equation=self._lambda_schedule.initial(),
+                    from .._utils._schedules import (
+                        ring_break_morph as _ring_break_morph,
                     )
 
-                    self._lambda_schedule.set_equation(
-                        stage="morph", lever="morse_hard", equation=0
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph", lever="morse_soft", equation=0
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="bond_k",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="bond_length",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="angle_k",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="angle_size",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="torsion_k",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="torsion_phase",
-                        equation=self._lambda_schedule.final(),
-                    )
+                    self._lambda_schedule = _ring_break_morph()
                     self._lambda_schedule_name = "ring_break_morph"
                 elif lambda_schedule == "reverse_ring_break_morph":
-                    self._lambda_schedule = _LambdaSchedule.standard_morph()
-                    self._lambda_schedule.set_equation(
-                        stage="morph", lever="morse_hard", equation=0
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph", lever="morse_soft", equation=0
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="bond_k",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="bond_length",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="angle_k",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="angle_size",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="torsion_k",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="morph",
-                        lever="torsion_phase",
-                        equation=self._lambda_schedule.initial(),
+                    from .._utils._schedules import (
+                        reverse_ring_break_morph as _reverse_ring_break_morph,
                     )
 
-                    self._lambda_schedule.append_stage(
-                        "bonded_perturb", self._lambda_schedule.final()
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="bonded_perturb",
-                        lever="morse_soft",
-                        equation=0 + self._lambda_schedule.lam(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="bonded_perturb", lever="morse_hard", equation=0
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="bonded_perturb",
-                        lever="bond_k",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="bonded_perturb",
-                        lever="bond_length",
-                        equation=self._lambda_schedule.initial(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="bonded_perturb",
-                        lever="angle_k",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="bonded_perturb",
-                        lever="angle_size",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="bonded_perturb",
-                        lever="torsion_k",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="bonded_perturb",
-                        lever="torsion_phase",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-
-                    self._lambda_schedule.append_stage(
-                        "potential_swap", self._lambda_schedule.final()
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="morse_hard",
-                        equation=0 + self._lambda_schedule.lam(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="morse_soft",
-                        equation=1 - self._lambda_schedule.lam(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="bond_k",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="bond_length",
-                        equation=(1 - self._lambda_schedule.lam())
-                        * self._lambda_schedule.initial()
-                        + self._lambda_schedule.lam() * self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="angle_k",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="angle_size",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="torsion_k",
-                        equation=self._lambda_schedule.final(),
-                    )
-                    self._lambda_schedule.set_equation(
-                        stage="potential_swap",
-                        lever="torsion_phase",
-                        equation=self._lambda_schedule.final(),
-                    )
+                    self._lambda_schedule = _reverse_ring_break_morph()
                     self._lambda_schedule_name = "reverse_ring_break_morph"
+                elif lambda_schedule == "annihilate":
+                    from .._utils._schedules import annihilate as _annihilate
+
+                    self._lambda_schedule = _annihilate()
+                    self._lambda_schedule_name = "annihilate"
+                elif lambda_schedule == "decouple":
+                    from .._utils._schedules import decouple as _decouple
+
+                    self._lambda_schedule = _decouple()
+                    self._lambda_schedule_name = "decouple"
                 else:
                     try:
                         self._lambda_schedule = self._from_hex(lambda_schedule)
@@ -1393,19 +1147,6 @@ def swap_end_states(self, swap_end_states):
             raise ValueError("'swap_end_states' must be of type 'bool'")
         self._swap_end_states = swap_end_states
 
-    @property
-    def coulomb_power(self):
-        return self._coulomb_power
-
-    @coulomb_power.setter
-    def coulomb_power(self, coulomb_power):
-        if not isinstance(coulomb_power, float):
-            try:
-                coulomb_power = float(coulomb_power)
-            except Exception:
-                raise ValueError("'coulomb_power' must be a of type 'float'")
-        self._coulomb_power = coulomb_power
-
     @property
     def shift_coulomb(self):
         return self._shift_coulomb
@@ -2381,6 +2122,21 @@ def taylor_power(self, taylor_power):
             raise ValueError("'taylor_power' must be between 0 and 4")
         self._taylor_power = taylor_power
 
+    @property
+    def beutler_alpha(self):
+        return self._beutler_alpha
+
+    @beutler_alpha.setter
+    def beutler_alpha(self, beutler_alpha):
+        if not isinstance(beutler_alpha, float):
+            try:
+                beutler_alpha = float(beutler_alpha)
+            except Exception:
+                raise ValueError("'beutler_alpha' must be of type 'float'")
+        if beutler_alpha < 0.0:
+            raise ValueError("'beutler_alpha' must be >= 0")
+        self._beutler_alpha = beutler_alpha
+
     @property
     def use_backup(self):
         return self._use_backup
diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index 571707e..0743756 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -202,10 +202,19 @@ def __init__(self, system, config):
             self._config._extra_args["use_gcmc_lrc"] = True
             self._config._extra_args["num_gcmc_waters"] = self._config.gcmc_num_waters
 
-        # If specified, use the Taylor soft-core form.
+        # Set the soft-core form.
         if self._config.softcore_form == "taylor":
             self._config._extra_args["use_taylor_softening"] = True
             self._config._extra_args["taylor_power"] = self._config.taylor_power
+        elif self._config.softcore_form == "beutler":
+            schedule_name = self._config._lambda_schedule_name
+            if schedule_name not in (None, "annihilate", "decouple"):
+                raise ValueError(
+                    "The Beutler soft-core form is only supported with the 'annihilate' "
+                    "or 'decouple' lambda schedules, or a custom schedule."
+                )
+            self._config._extra_args["use_beutler_softening"] = True
+            self._config._extra_args["beutler_alpha"] = self._config.beutler_alpha
 
         # We're running in SOMD1 compatibility mode.
         if self._config.somd1_compatibility:
@@ -803,7 +812,6 @@ def __init__(self, system, config):
 
         # Common kwargs shared by both dynamics and GCMC sampling.
         self._common_kwargs = {
-            "coulomb_power": self._config.coulomb_power,
             "cutoff": self._config.cutoff,
             "cutoff_type": self._config.cutoff_type,
             "platform": self._config.platform,

From 4c51884b8612e4d472cbb22e0e424760586eeb02 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 30 Apr 2026 12:08:54 +0100
Subject: [PATCH 04/12] Add warning for non-ring-breaking connectivity changes.

---
 src/somd2/runner/_base.py | 34 ++++++++++++++++++++++++++++++++++
 1 file changed, 34 insertions(+)

diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index 0743756..14ff545 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -178,6 +178,40 @@ def __init__(self, system, config):
         # Link properties to the lambda = 0 end state.
         self._system = _sr.morph.link_to_reference(self._system)
 
+        # Whether this is a ring-breaking schedule.
+        if (
+            self._config._lambda_schedule_name is not None
+            and "ring_breaking" in self._config._lambda_schedule_name
+        ):
+            self._is_ring_breaking = True
+        else:
+            self._is_ring_breaking = False
+
+        # Check to see if the end-state connectivities are the same.
+        if not self._is_ring_breaking:
+            for mol in self._system["property is_perturbable"].molecules():
+                has_end_state_connectivity = False
+                try:
+                    # The molecule will have two connectivity properties if
+                    # the merge detected a change in connectivity.
+                    c0 = mol.property("connectivity0")
+                    c1 = mol.property("connectivity1")
+                    has_end_state_connectivity = True
+                except:
+                    # No connectivity change detected.
+                    has_end_state_connectivity = False
+                    pass
+
+                # Check the connectivities regardless.
+                if has_end_state_connectivity:
+                    if c0 != c1:
+                        msg = (
+                            "End-state connectivities are different. If this is a ring-breaking "
+                            "perturbation, please set 'lambda_schedule_name' to 'ring_breaking'."
+                        )
+                        _logger.warning(msg)
+                        break
+
         # Set the default configuration options.
 
         # Restrict the atomic properties used to define light atoms when

From 82f3b952948a1c47a41788ead6455b08b77677b9 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 30 Apr 2026 14:48:04 +0100
Subject: [PATCH 05/12] Remove redundant option.

---
 src/somd2/runner/_base.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index 14ff545..9bce2b8 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -325,9 +325,7 @@ def __init__(self, system, config):
             # Angle optimisation can sometimes fail.
             except Exception as e1:
                 try:
-                    self._system, self._modifications = modify(
-                        self._system, optimise_angles=False
-                    )
+                    self._system, self._modifications = modify(self._system)
                 except Exception as e2:
                     msg = f"Unable to apply modifications to ghost atom bonded terms: {e1}; {e2}"
                     _logger.error(msg)

From 9dad004c2e43f2ae0fce22f1ec57d30d50d51381 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Fri, 1 May 2026 11:34:11 +0100
Subject: [PATCH 06/12] Add support for LRC scaling.

---
 src/somd2/_utils/_schedules.py | 54 ++++++++++++++++++++++++++--------
 src/somd2/config/_config.py    |  8 ++---
 src/somd2/runner/_base.py      | 11 +++++++
 3 files changed, 55 insertions(+), 18 deletions(-)

diff --git a/src/somd2/_utils/_schedules.py b/src/somd2/_utils/_schedules.py
index 9a22ce5..0ba0801 100644
--- a/src/somd2/_utils/_schedules.py
+++ b/src/somd2/_utils/_schedules.py
@@ -27,15 +27,22 @@
 ]
 
 
-def annihilate():
+def annihilate(fix_epsilon=True):
     """
-    Build the ABFE lambda schedule using decharge → annihilate with constant epsilon.
+    Build the ABFE lambda schedule using decharge → annihilate.
 
     Annihilation removes ALL non-bonded interactions (including intramolecular LJ
-    between non-bonded pairs), matching the GROMACS ABFE protocol. Epsilon is held
-    constant at its real-atom value throughout the annihilate stage so that the
-    (1-alpha) prefactor for the Buetler soft-core provides the sole LJ decay pathway,
-    giving true single (1-lambda) scaling consistent with GROMACS Beutler.
+    between non-bonded pairs).
+
+    Parameters
+    ----------
+    fix_epsilon : bool, optional
+        If True (default), epsilon is held constant at its real-atom value
+        throughout the annihilate stage so that the (1-alpha) prefactor of the
+        Beutler soft-core provides the sole LJ decay pathway.  The ghost-LRC
+        force is then explicitly scaled to zero over the stage to compensate.
+        If False, epsilon is scaled normally from initial to final and the LRC
+        follows naturally.
 
     Returns
     -------
@@ -66,18 +73,33 @@ def annihilate():
     )
     s.set_equation(stage="annihilate", lever="charge", equation=s.final())
     s.set_equation(stage="annihilate", force="restraint", equation=s.final())
-    s.set_equation(stage="annihilate", lever="epsilon", equation=s.initial())
+
+    if fix_epsilon:
+        s.set_equation(stage="annihilate", lever="epsilon", equation=s.initial())
+        s.set_equation(
+            stage="annihilate",
+            force="ghost-lrc",
+            lever="lrc_scale",
+            equation=1 - s.lam(),
+        )
 
     return s
 
 
-def decouple():
+def decouple(fix_epsilon=True):
     """
-    Build the ABFE lambda schedule using decharge → decouple with constant epsilon.
+    Build the ABFE lambda schedule using decharge → decouple.
 
     Decoupling removes only INTERMOLECULAR non-bonded interactions; intramolecular
-    terms are preserved via kappa=0 on ghost/ghost and ghost-14 forces. Epsilon is
-    held constant throughout the decouple stage (see annihilate for rationale).
+    terms are preserved via kappa=0 on ghost/ghost and ghost-14 forces.
+
+    Parameters
+    ----------
+    fix_epsilon : bool, optional
+        If True (default), epsilon is held constant at its real-atom value
+        throughout the decouple stage (see annihilate for rationale).  The
+        ghost-LRC force is then explicitly scaled to zero over the stage.
+        If False, epsilon is scaled normally and the LRC follows naturally.
 
     Returns
     -------
@@ -96,7 +118,15 @@ def decouple():
     s.set_equation(stage="decouple", lever="kappa", force="ghost/ghost", equation=0)
     s.set_equation(stage="decouple", lever="kappa", force="ghost-14", equation=0)
     s.set_equation(stage="decouple", lever="charge", equation=s.final())
-    s.set_equation(stage="decouple", lever="epsilon", equation=s.initial())
+
+    if fix_epsilon:
+        s.set_equation(stage="decouple", lever="epsilon", equation=s.initial())
+        s.set_equation(
+            stage="decouple",
+            force="ghost-lrc",
+            lever="lrc_scale",
+            equation=1 - s.lam(),
+        )
 
     s.prepend_stage("decharge", s.initial())
     s.set_equation(
diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py
index 6d20615..9a605b7 100644
--- a/src/somd2/config/_config.py
+++ b/src/somd2/config/_config.py
@@ -1093,14 +1093,10 @@ def lambda_schedule(self, lambda_schedule):
                     self._lambda_schedule = _reverse_ring_break_morph()
                     self._lambda_schedule_name = "reverse_ring_break_morph"
                 elif lambda_schedule == "annihilate":
-                    from .._utils._schedules import annihilate as _annihilate
-
-                    self._lambda_schedule = _annihilate()
+                    self._lambda_schedule = None
                     self._lambda_schedule_name = "annihilate"
                 elif lambda_schedule == "decouple":
-                    from .._utils._schedules import decouple as _decouple
-
-                    self._lambda_schedule = _decouple()
+                    self._lambda_schedule = None
                     self._lambda_schedule_name = "decouple"
                 else:
                     try:
diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index 9bce2b8..8029f71 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -250,6 +250,17 @@ def __init__(self, system, config):
             self._config._extra_args["use_beutler_softening"] = True
             self._config._extra_args["beutler_alpha"] = self._config.beutler_alpha
 
+        # Build deferred schedules now that the softcore form is known.
+        fix_epsilon = self._config.softcore_form == "beutler"
+        if self._config._lambda_schedule_name == "annihilate":
+            from .._utils._schedules import annihilate as _annihilate
+
+            self._config._lambda_schedule = _annihilate(fix_epsilon=fix_epsilon)
+        elif self._config._lambda_schedule_name == "decouple":
+            from .._utils._schedules import decouple as _decouple
+
+            self._config._lambda_schedule = _decouple(fix_epsilon=fix_epsilon)
+
         # We're running in SOMD1 compatibility mode.
         if self._config.somd1_compatibility:
             from .._utils._somd1 import make_compatible

From e4de0a9699ac2458724a962fb5679094c0b728cb Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Tue, 5 May 2026 09:42:29 +0100
Subject: [PATCH 07/12] Disable LRC for vacuum simulations.

---
 src/somd2/runner/_base.py | 35 +++++++++++++++++------------------
 1 file changed, 17 insertions(+), 18 deletions(-)

diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index 8029f71..885cedd 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -212,6 +212,23 @@ def __init__(self, system, config):
                         _logger.warning(msg)
                         break
 
+        # Check for a periodic space.
+        self._has_space = self._check_space()
+
+        # Check for water.
+        try:
+            # The search will fail if there are no water molecules.
+            water = self._system["water"].molecules()
+            self._has_water = True
+        except:
+            self._has_water = False
+
+        # Warn if dispersion correction is requested but can't be applied.
+        if self._config.use_dispersion_correction and not self._has_water:
+            msg = "Cannot use dispersion correction for vacuum simulations. Disabling!"
+            _logger.warning(msg)
+            self._config.use_dispersion_correction = False
+
         # Set the default configuration options.
 
         # Restrict the atomic properties used to define light atoms when
@@ -342,24 +359,6 @@ def __init__(self, system, config):
                     _logger.error(msg)
                     raise RuntimeError(msg)
 
-        # Check for a periodic space.
-        self._has_space = self._check_space()
-
-        # Check for water.
-        try:
-            # The search will fail if there are no water molecules.
-            water = self._system["water"].molecules()
-            self._has_water = True
-        except:
-            self._has_water = False
-
-        # Warn if dispersion correction is requested but can't be applied.
-        if self._config.use_dispersion_correction and not self._has_space:
-            _logger.warning(
-                "'use_dispersion_correction=True' has no effect: the system "
-                "has no periodic space. The option will be ignored."
-            )
-
         # Check the end state constraints.
         self._check_end_state_constraints()
 

From d111e032b83f655618b366c746a41ade46b73430 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Tue, 5 May 2026 09:42:58 +0100
Subject: [PATCH 08/12] Add functionality for skipping/warning deprecated
 options.

---
 src/somd2/io/_io.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/src/somd2/io/_io.py b/src/somd2/io/_io.py
index e845e6a..3ccb9eb 100644
--- a/src/somd2/io/_io.py
+++ b/src/somd2/io/_io.py
@@ -34,8 +34,15 @@
 import pyarrow as _pa
 import pyarrow.parquet as _pq
 import pandas as _pd
+import warnings as _warnings
 import yaml as _yaml
 
+# Options that have been removed from the config. Any of these found in a YAML
+# config file will be silently dropped after emitting a deprecation warning.
+_REMOVED_OPTIONS = {
+    "coulomb_power": "'coulomb_power' has been removed and will be ignored.",
+}
+
 
 def dataframe_to_parquet(df, metadata, filepath=None, filename=None):
     """
@@ -142,6 +149,11 @@ def yaml_to_dict(path):
     except Exception as e:
         raise ValueError(f"Could not load YAML file: {e}")
 
+    for key, msg in _REMOVED_OPTIONS.items():
+        if key in d:
+            _warnings.warn(msg)
+            d.pop(key)
+
     return d
 
 

From 3fecef7458c811bee8d9582a257c743dc1aebc4a Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Tue, 5 May 2026 09:48:39 +0100
Subject: [PATCH 09/12] Remove coulomb_power option from unit test.

---
 tests/runner/test_restart.py | 7 -------
 1 file changed, 7 deletions(-)

diff --git a/tests/runner/test_restart.py b/tests/runner/test_restart.py
index bfb5fd0..6c1e2f6 100644
--- a/tests/runner/test_restart.py
+++ b/tests/runner/test_restart.py
@@ -123,13 +123,6 @@ def test_restart(mols, request):
         with pytest.raises(ValueError):
             runner_constraints = Runner(mols, Config(**config_diffconstraint))
 
-        config_diffcoulombpower = config_new.copy()
-        config_diffcoulombpower["runtime"] = "36fs"
-        config_diffcoulombpower["coulomb_power"] = 0.5
-
-        with pytest.raises(ValueError):
-            runner_coulombpower = Runner(mols, Config(**config_diffcoulombpower))
-
         config_diffcutofftype = config_new.copy()
         config_diffcutofftype["runtime"] = "36fs"
         config_diffcutofftype["cutoff_type"] = "rf"

From 55971646f547524b60a1421713dac2a4deb51107 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Tue, 5 May 2026 11:36:32 +0100
Subject: [PATCH 10/12] Fix schedule name check.

---
 src/somd2/runner/_base.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index 885cedd..c69c9ec 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -181,7 +181,7 @@ def __init__(self, system, config):
         # Whether this is a ring-breaking schedule.
         if (
             self._config._lambda_schedule_name is not None
-            and "ring_breaking" in self._config._lambda_schedule_name
+            and "ring_break" in self._config._lambda_schedule_name
         ):
             self._is_ring_breaking = True
         else:

From db7bdeefb13202bfee5bc1bc34df1537f17a76d4 Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 7 May 2026 09:13:26 +0100
Subject: [PATCH 11/12] Reverse LambdaSchedule when swap_end_states=True. [ref
 #78]

---
 src/somd2/runner/_base.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index c69c9ec..d78dbdf 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -906,6 +906,17 @@ def __init__(self, system, config):
         else:
             self._gcmc_kwargs = None
 
+        # Reverse the lambda schedule when swapping end states so that the
+        # schedule progresses from the perturbed end state to the reference.
+        if self._config.swap_end_states:
+            self._dynamics_kwargs["schedule"] = self._dynamics_kwargs[
+                "schedule"
+            ].reverse()
+            if self._gcmc_kwargs is not None:
+                self._gcmc_kwargs["lambda_schedule"] = self._gcmc_kwargs[
+                    "lambda_schedule"
+                ].reverse()
+
         # Limit the number of CPU threads available to Sire when running in parallel.
         if self._is_gpu:
             # First get the total number of threads that are available to Sire.

From 5d410f480ee918fb6c195650f3105d40d7aebded Mon Sep 17 00:00:00 2001
From: Lester Hedges <lester.hedges@gmail.com>
Date: Thu, 7 May 2026 15:10:06 +0100
Subject: [PATCH 12/12] Let loch handle it's own schedule reversal.

---
 src/somd2/runner/_base.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index d78dbdf..9e9f2c6 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -908,14 +908,11 @@ def __init__(self, system, config):
 
         # Reverse the lambda schedule when swapping end states so that the
         # schedule progresses from the perturbed end state to the reference.
+        # (The GCMC schedule is reversed inside loch itself.)
         if self._config.swap_end_states:
             self._dynamics_kwargs["schedule"] = self._dynamics_kwargs[
                 "schedule"
             ].reverse()
-            if self._gcmc_kwargs is not None:
-                self._gcmc_kwargs["lambda_schedule"] = self._gcmc_kwargs[
-                    "lambda_schedule"
-                ].reverse()
 
         # Limit the number of CPU threads available to Sire when running in parallel.
         if self._is_gpu: