Skip to content

Switch to reversible three-stage ring-breaking schedule #158

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lohedges merged 1 commit into from
Jun 16, 2026
Merged

Switch to reversible three-stage ring-breaking schedule #158

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
97 changes: 38 additions & 59 deletions src/somd2/_utils/_schedules.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,14 +151,16 @@ def ring_break_morph():
"""
Build a lambda schedule for ring-breaking perturbations.

Four stages: potential_swap → restraints_off → ring_open → morph.
Three stages: potential_swap → restraints_off → morph.

The ring-break softcore kappa ramps 0→1 through ring_open and is fixed at 1
in morph. The ring-make equations mirror ring-break so that
``ring_break_morph().reverse()`` is the correct schedule for the ring-making
direction (used by :func:`reverse_ring_break_morph`). Because ring_break_morph
is only used for ring-breaking perturbations (no ring-make force present), the
ring-make equations have no effect on forward simulations.
During restraints_off the Morse restraint ramps off (morse_soft: 1→0) while
the ring-break softcore simultaneously ramps on (alpha: 1→0, kappa: 0→1),
equations mirror ring-break so that ``ring_break_morph().reverse()`` is the
providing a smooth handover with no gap between the two forces. The ring-make
correct schedule for the ring-making direction (used by
:func:`reverse_ring_break_morph`). Because ring_break_morph is only used for
ring-breaking perturbations (no ring-make force present), the ring-make
equations have no effect on forward simulations.

Returns
-------
Expand All @@ -169,41 +171,13 @@ def ring_break_morph():
from sire.cas import LambdaSchedule as _LambdaSchedule

s = _LambdaSchedule.standard_morph()
s.set_stage_weight("morph", 2)

# ring_open: Morse is already off; ring-break nonbonded interaction ramps
# on (alpha: 1→0, kappa: 0→1) while non-bonded terms stay at initial and
# bonded terms remain at final. The softcore interaction gently pushes the
# atoms into the open-chain geometry before the full nonbonded morph begins,
# improving HREX overlap at the ring-break boundary.
#
# ring-make equations mirror ring-break so the reversed schedule is correct:
# after .reverse(), ring-make kappa ramps 1→0 through this stage, matching
# what ring-break does here in the forward direction.
s.prepend_stage("ring_open", s.initial(), weight=1)
s.set_equation(stage="ring_open", lever="morse_hard", equation=0)
s.set_equation(stage="ring_open", lever="morse_soft", equation=0)
s.set_equation(stage="ring_open", lever="bond_k", equation=s.final())
s.set_equation(stage="ring_open", lever="bond_length", equation=s.final())
s.set_equation(stage="ring_open", lever="angle_k", equation=s.final())
s.set_equation(stage="ring_open", lever="angle_size", equation=s.final())
s.set_equation(stage="ring_open", lever="torsion_k", equation=s.final())
s.set_equation(stage="ring_open", lever="torsion_phase", equation=s.final())
s.set_equation(
stage="ring_open", force="ring-break", lever="alpha", equation=1 - s.lam()
)
s.set_equation(
stage="ring_open", force="ring-break", lever="kappa", equation=s.lam()
)
# ring-make mirrors ring-break so reversed schedule ramps ring-make 1→0 here.
s.set_equation(
stage="ring_open", force="ring-make", lever="alpha", equation=1 - s.lam()
)
s.set_equation(
stage="ring_open", force="ring-make", lever="kappa", equation=s.lam()
)

s.prepend_stage("restraints_off", s.initial(), weight=1)
# restraints_off [1/3, 2/3): Morse ramps off while ring-break softcore ramps
# on simultaneously (alpha: 1→0, kappa: 0→1). Bonded terms (angles, torsions)
# interpolate initial→final over the same stage. ring-make mirrors ring-break
# so that after .reverse(), the ring-make softcore ramps off as morse_soft ramps
# on in the reversed restraints_off stage, correct for ring-making perturbations.
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())
Expand All @@ -228,8 +202,20 @@ def ring_break_morph():
lever="torsion_phase",
equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(),
)
s.set_equation(
stage="restraints_off", force="ring-break", lever="alpha", equation=1 - s.lam()
)
s.set_equation(
stage="restraints_off", force="ring-break", lever="kappa", equation=s.lam()
)
s.set_equation(
stage="restraints_off", force="ring-make", lever="alpha", equation=1 - s.lam()
)
s.set_equation(
stage="restraints_off", force="ring-make", lever="kappa", equation=s.lam()
)

s.prepend_stage("potential_swap", s.initial(), weight=2)
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(
Expand All @@ -247,10 +233,9 @@ def ring_break_morph():
s.set_equation(stage="potential_swap", lever="torsion_k", equation=s.initial())
s.set_equation(stage="potential_swap", lever="torsion_phase", equation=s.initial())

# morph: standard nonbonded morphing. Ring-break is fixed at fully open
# (kappa=1, alpha=0) since geometry has already relaxed in ring_open.
# ring-make mirrors ring-break: kappa=1, alpha=0 so that .reverse() gives
# kappa=1 at lam=0 of the reversed morph stage (ring-making direction start).
# morph [2/3, 1]: standard nonbonded morphing with ring-break/ring-make fixed
# at fully open (kappa=1, alpha=0). ring-make mirrors ring-break so .reverse()
# gives kappa=1 at lam=0 of the reversed morph stage (ring-making start).
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())
Expand All @@ -264,21 +249,16 @@ def ring_break_morph():
s.set_equation(stage="morph", force="ring-make", lever="alpha", equation=0)
s.set_equation(stage="morph", force="ring-make", lever="kappa", equation=1)

# coul_kappa decouples Coulomb onset from LJ onset: zero throughout the
# bonded stages so the CLJ exception carries no charge while atoms are at
# covalent distances, then ramps 0→1 during morph only (where the LJ
# softcore has already separated the atoms). ring-make mirrors ring-break
# so that .reverse() gives the correct reversed schedule (coul_kappa ramps
# 1→0 through the reversed morph stage for the ring-making direction).
# coul_kappa: zero through both bonded stages so the CLJ exception carries no
# charge while atoms are at covalent distances; ramps 0→1 in morph only once
# the softcore has already separated the atoms. ring-make mirrors ring-break
# so .reverse() gives coul_kappa ramps 1→0 through the reversed morph stage.
s.set_equation(
stage="potential_swap", force="ring-break", lever="coul_kappa", equation=0
)
s.set_equation(
stage="restraints_off", force="ring-break", lever="coul_kappa", equation=0
)
s.set_equation(
stage="ring_open", force="ring-break", lever="coul_kappa", equation=0
)
s.set_equation(
stage="morph", force="ring-break", lever="coul_kappa", equation=s.lam()
)
Expand All @@ -288,7 +268,6 @@ def ring_break_morph():
s.set_equation(
stage="restraints_off", force="ring-make", lever="coul_kappa", equation=0
)
s.set_equation(stage="ring_open", force="ring-make", lever="coul_kappa", equation=0)
s.set_equation(
stage="morph", force="ring-make", lever="coul_kappa", equation=s.lam()
)
Expand All @@ -300,9 +279,9 @@ def reverse_ring_break_morph():
"""
Build a lambda schedule for ring-making perturbations (reverse ring-break).

Returns ``ring_break_morph().reverse()``: four stages in reversed order
(morph → ring_open → restraints_off → potential_swap) with all equations
reflected about λ=½ and initial/final end-states swapped.
Returns ``ring_break_morph().reverse()``: three stages in reversed order
(morph → restraints_off → potential_swap) with all equations reflected about
λ=½ and initial/final end-states swapped.

This schedule is correct for two equivalent use-cases:

Expand Down
73 changes: 37 additions & 36 deletions tests/schedules/test_ring_break.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,47 +129,49 @@ def test_reverse_has_ring_make_not_ring_break(reverse_dynamics):
# They are completely independent of Sire's energy formula and will continue to
# work correctly regardless of changes to the softcore implementation.

# ring_break_morph kappa/alpha points:
# λ=0.00 potential_swap start kappa=0, alpha=1
# λ=0.15 potential_swap mid kappa=0, alpha=1
# λ=1/3 restraints_off start kappa=0, alpha=1
# λ=0.45 restraints_off mid kappa=0, alpha=1
# λ=0.50 ring_open start kappa=0, alpha=1 (within-stage lam=0)
# λ=0.55 ring_open mid kappa=0.3, alpha=0.7
# λ=0.60 ring_open mid kappa=0.6, alpha=0.4
# λ=2/3 morph start kappa=1, alpha=0
# λ=0.85 morph mid kappa=1, alpha=0
# λ=1.00 morph end kappa=1, alpha=0
# ring_break_morph kappa/alpha points (3 equal stages: [0,1/3), [1/3,2/3), [2/3,1]):
# λ=0.00 potential_swap start kappa=0, alpha=1
# λ=0.15 potential_swap mid kappa=0, alpha=1
# λ=1/3 restraints_off start kappa=0, alpha=1 (within-stage lam=0)
# λ=0.45 restraints_off mid kappa=0.35, alpha=0.65 (within-stage lam=0.35)
# λ=0.50 restraints_off mid kappa=0.5, alpha=0.5 (within-stage lam=0.5)
# λ=0.55 restraints_off mid kappa=0.65, alpha=0.35 (within-stage lam=0.65)
# λ=0.60 restraints_off near end kappa=0.8, alpha=0.2 (within-stage lam=0.8)
# λ=2/3 morph start kappa=1, alpha=0
# λ=0.85 morph mid kappa=1, alpha=0
# λ=1.00 morph end kappa=1, alpha=0
_FWD_KAPPA_ALPHA = [
(0.00, 0.0, 1.0),
(0.15, 0.0, 1.0),
(1 / 3, 0.0, 1.0),
(0.45, 0.0, 1.0),
(0.50, 0.0, 1.0),
(0.55, 0.3, 0.7),
(0.60, 0.6, 0.4),
(0.45, 0.35, 0.65),
(0.50, 0.5, 0.5),
(0.55, 0.65, 0.35),
(0.60, 0.8, 0.2),
(2 / 3, 1.0, 0.0),
(0.85, 1.0, 0.0),
(1.00, 1.0, 0.0),
]

# reverse_ring_break_morph ring-make kappa/alpha points (mirror of forward):
# λ=0.00 morph start kappa=1, alpha=0
# λ=0.15 morph mid kappa=1, alpha=0
# λ=1/3 ring_open start kappa=1, alpha=0 (within-stage lam=0)
# λ=0.45 ring_open mid kappa=0.3, alpha=0.7 (within-stage lam=0.7)
# λ=0.50 restraints_off start kappa=0, alpha=1
# λ=0.60 restraints_off mid kappa=0, alpha=1
# λ=2/3 potential_swap start kappa=0, alpha=1
# λ=0.85 potential_swap mid kappa=0, alpha=1
# λ=1.00 potential_swap end kappa=0, alpha=1
# λ=0.00 reversed morph start kappa=1, alpha=0
# λ=0.15 reversed morph mid kappa=1, alpha=0
# λ=1/3 reversed restraints_off start kappa=1, alpha=0 (within-stage lam=0)
# λ=0.45 reversed restraints_off mid kappa=0.65, alpha=0.35
# λ=0.50 reversed restraints_off mid kappa=0.5, alpha=0.5
# λ=0.55 reversed restraints_off mid kappa=0.35, alpha=0.65
# λ=0.60 reversed restraints_off near end kappa=0.2, alpha=0.8
# λ=2/3 reversed potential_swap start kappa=0, alpha=1
# λ=0.85 reversed potential_swap mid kappa=0, alpha=1
# λ=1.00 reversed potential_swap end kappa=0, alpha=1
_REV_KAPPA_ALPHA = [
(0.00, 1.0, 0.0),
(0.15, 1.0, 0.0),
(1 / 3, 1.0, 0.0),
(0.45, 0.3, 0.7),
(0.50, 0.0, 1.0),
(0.60, 0.0, 1.0),
(0.45, 0.65, 0.35),
(0.50, 0.5, 0.5),
(0.55, 0.35, 0.65),
(0.60, 0.2, 0.8),
(2 / 3, 0.0, 1.0),
(0.85, 0.0, 1.0),
(1.00, 0.0, 1.0),
Expand All @@ -194,10 +196,10 @@ def test_reverse_has_ring_make_not_ring_break(reverse_dynamics):
]

# reverse_ring_break_morph ring-make coul_kappa points (initial=1, final=0):
# λ=0.00 morph start (reversed) coul_kappa=1.0
# λ=0.15 morph mid coul_kappa=0.55 (1 - 0.15*3)
# λ=1/3 morph end / ring_open coul_kappa=0.0
# λ=0.45–1.0 ring_open/restraints_off/potential_swap coul_kappa=0
# λ=0.00 reversed morph start coul_kappa=1.0
# λ=0.15 reversed morph mid coul_kappa=0.55 (1 - 0.15*3)
# λ=1/3 reversed morph end coul_kappa=0.0
# λ=0.45–1.0 restraints_off/potential_swap coul_kappa=0
_REV_COUL_KAPPA = [
(0.00, 1.0),
(0.15, 0.55),
Expand Down Expand Up @@ -293,10 +295,9 @@ def test_reverse_ring_break_morph_coul_kappa(lam, expected_coul_kappa):


@pytest.mark.parametrize("lam", [2 / 3, 1.0])
def test_ring_break_active_after_ring_open(forward_dynamics, lam):
def test_ring_break_active_in_morph(forward_dynamics, lam):
"""
Ring-break energy is clearly non-zero (kappa=1) at the end of ring_open
and throughout morph.
Ring-break energy is clearly non-zero (kappa=1) throughout the morph stage.
"""
e = _force_energy_kcal(forward_dynamics, lam, "ring-break")
assert abs(e) > _ACTIVE_THRESHOLD, (
Expand Down Expand Up @@ -343,8 +344,8 @@ def test_ring_make_inactive_at_lambda_one(reverse_dynamics):
#
# Test points span zero and non-zero energy regions:
# λ=0.0 → forward kappa=0, reverse at 1-λ=1.0 kappa=0 (both ≈0)
# λ=0.55 → forward ring_open (kappa=0.3), reverse ring_open at 0.45 (kappa=0.3)
# λ=2/3 → forward morph start (kappa=1), reverse ring_open start at 1/3 (kappa=1)
# λ=0.55 → forward restraints_off (kappa=0.65), reverse restraints_off at 0.45 (kappa=0.65)
# λ=2/3 → forward morph start (kappa=1), reverse restraints_off start at 1/3 (kappa=1)
# λ=0.85 → forward morph (kappa=1), reverse reversed-morph at 0.15 (kappa=1)
# λ=1.0 → forward morph end (kappa=1), reverse at 0.0 reversed-morph (kappa=1)

Expand Down
Loading