diff --git a/src/somd2/_utils/_schedules.py b/src/somd2/_utils/_schedules.py
new file mode 100644
index 0000000..0ba0801
--- /dev/null
+++ b/src/somd2/_utils/_schedules.py
@@ -0,0 +1,291 @@
+######################################################################
+# 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(fix_epsilon=True):
+    """
+    Build the ABFE lambda schedule using decharge → annihilate.
+
+    Annihilation removes ALL non-bonded interactions (including intramolecular LJ
+    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
+    -------
+
+    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())
+
+    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(fix_epsilon=True):
+    """
+    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.
+
+    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
+    -------
+
+    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())
+
+    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(
+        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 fce678e..9a605b7 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,
@@ -151,10 +152,12 @@ 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",
         taylor_power=1,
+        beutler_alpha=0.5,
         output_directory="output",
         restart=False,
         use_backup=False,
@@ -233,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.
@@ -438,6 +437,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
@@ -458,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.
 
@@ -559,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
@@ -605,12 +615,14 @@ 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
         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
@@ -1067,279 +1079,25 @@ 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(),
+                    from .._utils._schedules import (
+                        ring_break_morph as _ring_break_morph,
                     )
 
-                    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(),
-                    )
-
-                    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(),
-                    )
-
-                    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(),
+                    from .._utils._schedules import (
+                        reverse_ring_break_morph as _reverse_ring_break_morph,
                     )
 
-                    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":
+                    self._lambda_schedule = None
+                    self._lambda_schedule_name = "annihilate"
+                elif lambda_schedule == "decouple":
+                    self._lambda_schedule = None
+                    self._lambda_schedule_name = "decouple"
                 else:
                     try:
                         self._lambda_schedule = self._from_hex(lambda_schedule)
@@ -1385,19 +1143,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
@@ -2286,6 +2031,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
@@ -2363,6 +2118,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/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
 
 
diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py
index a6c637c..9e9f2c6 100644
--- a/src/somd2/runner/_base.py
+++ b/src/somd2/runner/_base.py
@@ -178,6 +178,57 @@ 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_break" 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
+
+        # 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
@@ -192,10 +243,40 @@ def __init__(self, system, config):
             self._config.fix_perturbable_zero_sigmas
         )
 
-        # If specified, use the Taylor soft-core form.
+        # 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
+
+        # 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
+
+        # 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:
@@ -272,25 +353,12 @@ 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)
                     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
-
         # Check the end state constraints.
         self._check_end_state_constraints()
 
@@ -786,7 +854,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,
@@ -839,6 +906,14 @@ 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.
+        # (The GCMC schedule is reversed inside loch itself.)
+        if self._config.swap_end_states:
+            self._dynamics_kwargs["schedule"] = self._dynamics_kwargs[
+                "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.
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"