diff --git a/src/somd2/_utils/_schedules.py b/src/somd2/_utils/_schedules.py index 5207a84..1c4e48d 100644 --- a/src/somd2/_utils/_schedules.py +++ b/src/somd2/_utils/_schedules.py @@ -153,6 +153,13 @@ def ring_break_morph(): Four stages: potential_swap → restraints_off → ring_open → 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. + Returns ------- @@ -169,6 +176,10 @@ def ring_break_morph(): # 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) @@ -184,11 +195,12 @@ def ring_break_morph(): 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=s.lam() + stage="ring_open", force="ring-make", lever="alpha", equation=1 - s.lam() ) s.set_equation( - stage="ring_open", force="ring-make", lever="kappa", equation=1 - s.lam() + stage="ring_open", force="ring-make", lever="kappa", equation=s.lam() ) s.prepend_stage("restraints_off", s.initial(), weight=1) @@ -237,6 +249,8 @@ def ring_break_morph(): # 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). 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()) @@ -247,118 +261,67 @@ def ring_break_morph(): s.set_equation(stage="morph", lever="torsion_phase", equation=s.final()) s.set_equation(stage="morph", force="ring-break", lever="alpha", equation=0) s.set_equation(stage="morph", force="ring-break", lever="kappa", equation=1) - s.set_equation(stage="morph", force="ring-make", lever="alpha", equation=1) - s.set_equation(stage="morph", force="ring-make", lever="kappa", equation=0) - - return s - - -def reverse_ring_break_morph(): - """ - Build a lambda schedule for reverse ring-breaking perturbations. - - Four stages: morph → ring_close → 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_stage_weight("morph", 2) - - # morph: standard nonbonded morphing. Ring-break fixed at fully open - # (kappa=1, alpha=0); ring-make fixed at full interaction (kappa=1, alpha=0). - 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.set_equation(stage="morph", force="ring-break", lever="alpha", equation=1) - s.set_equation(stage="morph", force="ring-break", lever="kappa", equation=0) s.set_equation(stage="morph", force="ring-make", lever="alpha", equation=0) s.set_equation(stage="morph", force="ring-make", lever="kappa", equation=1) - # ring_close: non-bonded terms fixed at final; ring-make interaction ramps - # off (alpha: 0→1, kappa: 1→0) to allow atoms to relax into ring geometry - # before Morse is applied. Symmetric counterpart to ring_open. - s.append_stage("ring_close", s.final(), weight=1) - s.set_equation(stage="ring_close", lever="morse_hard", equation=0) - s.set_equation(stage="ring_close", lever="morse_soft", equation=0) - s.set_equation(stage="ring_close", lever="bond_k", equation=s.initial()) - s.set_equation(stage="ring_close", lever="bond_length", equation=s.initial()) - s.set_equation(stage="ring_close", lever="angle_k", equation=s.initial()) - s.set_equation(stage="ring_close", lever="angle_size", equation=s.initial()) - s.set_equation(stage="ring_close", lever="torsion_k", equation=s.initial()) - s.set_equation(stage="ring_close", lever="torsion_phase", equation=s.initial()) - s.set_equation(stage="ring_close", force="ring-break", lever="alpha", equation=1) - s.set_equation(stage="ring_close", force="ring-break", lever="kappa", equation=0) + # 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). s.set_equation( - stage="ring_close", force="ring-make", lever="alpha", equation=s.lam() + stage="potential_swap", force="ring-break", lever="coul_kappa", equation=0 ) s.set_equation( - stage="ring_close", force="ring-make", lever="kappa", equation=1 - s.lam() + stage="restraints_off", force="ring-break", lever="coul_kappa", equation=0 ) - - # bonded_perturb: Morse soft ramps on; ring-make already off from ring_close. - # Ring-break softcore turns on as any ring-break bond opens (alpha: 1→0, - # kappa: 0→1); ring-make stays off (kappa=0, alpha=1). - s.append_stage("bonded_perturb", s.final(), weight=1) - 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(), + stage="ring_open", force="ring-break", lever="coul_kappa", equation=0 ) s.set_equation( - stage="bonded_perturb", - lever="angle_size", - equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(), + stage="morph", force="ring-break", lever="coul_kappa", equation=s.lam() ) s.set_equation( - stage="bonded_perturb", - lever="torsion_k", - equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(), + stage="potential_swap", force="ring-make", lever="coul_kappa", equation=0 ) s.set_equation( - stage="bonded_perturb", - lever="torsion_phase", - equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(), - ) - s.set_equation( - stage="bonded_perturb", force="ring-break", lever="alpha", equation=1 - s.lam() + 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="bonded_perturb", force="ring-break", lever="kappa", equation=s.lam() + stage="morph", force="ring-make", lever="coul_kappa", equation=s.lam() ) - s.set_equation(stage="bonded_perturb", force="ring-make", lever="alpha", equation=1) - s.set_equation(stage="bonded_perturb", force="ring-make", lever="kappa", equation=0) - - s.append_stage("potential_swap", s.final(), weight=2) - 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 + + +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. + + This schedule is correct for two equivalent use-cases: + + 1. A ring-making perturbation run with ``swap_end_states=False``: the + ring-make softcore force (kappa=1 at λ=0, ramping to 0) is controlled + directly by the ring-make lever equations. + 2. A ring-breaking perturbation run with ``swap_end_states=True`` (the + runner reverses the schedule automatically, yielding the same effective + schedule): the ring-make softcore — which now controls the original + ring-breaking bond after the end-state swap — is handled identically. + + The energy symmetry invariant holds for both cases: + ``E_ring_make_reverse(λ) == E_ring_break_forward(1-λ)`` at any fixed + geometry. + + Returns + ------- + + schedule : sire.legacy.CAS.LambdaSchedule + The lambda schedule. + """ + return ring_break_morph().reverse() diff --git a/tests/schedules/test_ring_break.py b/tests/schedules/test_ring_break.py index 6f589f2..044136b 100644 --- a/tests/schedules/test_ring_break.py +++ b/tests/schedules/test_ring_break.py @@ -1,29 +1,3 @@ -""" -Tests for ring-break / ring-make CustomBondForce setup and schedule behaviour. - -Checks: - - ``ring-break`` force is registered (not ``ring-make``) for the forward - direction (swap_end_states=False, ring_break_morph schedule). - - ``ring-make`` force is registered (not ``ring-break``) for the reverse - direction (swap_end_states=True, reverse_ring_break_morph schedule). - - Ring-break energy is near-zero when kappa=0 (potential_swap and - restraints_off stages) and clearly non-zero when kappa=1 (morph stage). - - Ring-make energy follows the symmetric pattern under reverse_ring_break_morph: - non-zero at λ=0 (morph, kappa=1) and near-zero at λ=1 (potential_swap, kappa=0). - -Stage boundaries for ring_break_morph (weights 2:1:1:2, total 6): - potential_swap λ ∈ [0, 1/3) ring-break kappa = 0 - restraints_off λ ∈ [1/3, 1/2) ring-break kappa = 0 - ring_open λ ∈ [1/2, 2/3) ring-break kappa ramps 0 → 1 - morph λ ∈ [2/3, 1] ring-break kappa = 1 (fixed) - -Stage boundaries for reverse_ring_break_morph (weights 2:1:1:2, total 6): - morph λ ∈ [0, 1/3) ring-make kappa = 1 (fixed) - ring_close λ ∈ [1/3, 1/2) ring-make kappa ramps 1 → 0 - bonded_perturb λ ∈ [1/2, 2/3) ring-make kappa = 0 - potential_swap λ ∈ [2/3, 1] ring-make kappa = 0 -""" - import pytest import sire as sr @@ -32,24 +6,14 @@ reason="openmm support is not available", ) -# Energy threshold (kcal/mol) separating "inactive" from "active" force states. -# Chosen conservatively: sp_scan shows kappa=0 energies up to ~0.07 kcal/mol -# and kappa=1 energies of at least ~0.19 kcal/mol at the ring_open/morph boundary. -_INACTIVE_THRESHOLD = 0.1 # abs(E) must be below this when kappa=0 -_ACTIVE_THRESHOLD = 0.1 # abs(E) must be above this when kappa=1 - - -# ── helpers ────────────────────────────────────────────────────────────────── +# Energy threshold (kcal/mol) for the "active" state: kappa=1 should give +# clearly non-zero CustomBondForce energy. +_ACTIVE_THRESHOLD = 0.1 def _build_dynamics(mols, schedule, swap_end_states): """ Construct a dynamics context for the ring-break system with Morse restraints. - - Mirrors the production setup in somd2_api_runner_dmr.py / sp_scan.py so - that the force groups match what a real simulation would create. - Uses reaction-field (rf) rather than PME since the test system is a - vacuum/non-periodic input with no periodic box vectors. """ from somd2._utils._somd1 import make_compatible @@ -106,9 +70,6 @@ def _force_energy_kcal(d, lam, force_name): return state.getPotentialEnergy().value_in_unit(openmm.unit.kilocalories_per_mole) -# ── fixtures ───────────────────────────────────────────────────────────────── - - @pytest.fixture(scope="module") def forward_dynamics(syk_ring_break_mols): """Forward ring-breaking dynamics: swap_end_states=False, ring_break_morph.""" @@ -121,7 +82,13 @@ def forward_dynamics(syk_ring_break_mols): @pytest.fixture(scope="module") def reverse_dynamics(syk_ring_break_mols): - """Reverse ring-making dynamics: swap_end_states=True, reverse_ring_break_morph.""" + """ + Reverse ring-making dynamics: swap_end_states=True, reverse_ring_break_morph. + + reverse_ring_break_morph() == ring_break_morph().reverse(), so this fixture + also implicitly tests the reversed schedule path used by the runner when a + ring-breaking perturbation is run with swap_end_states=True. + """ from somd2._utils._schedules import reverse_ring_break_morph return _build_dynamics( @@ -155,56 +122,187 @@ def test_reverse_has_ring_make_not_ring_break(reverse_dynamics): ) -# ── forward schedule energy tests ───────────────────────────────────────────── +# ── schedule kappa/alpha tests ──────────────────────────────────────────────── +# +# These tests verify kappa and alpha values by calling schedule.morph() directly, +# using the same initial/final values that lambdalever passes in production. +# 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 +_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), + (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 +_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), + (2 / 3, 0.0, 1.0), + (0.85, 0.0, 1.0), + (1.00, 0.0, 1.0), +] + +# ring_break_morph coul_kappa points (initial=0, final=1): +# λ=0.00–2/3 all pre-morph stages coul_kappa=0 +# λ=2/3 morph start coul_kappa=0 (within-stage lam=0) +# λ=0.85 morph mid coul_kappa=0.55 ((0.85-2/3)*3) +# λ=1.00 morph end coul_kappa=1.0 +_FWD_COUL_KAPPA = [ + (0.00, 0.0), + (0.15, 0.0), + (1 / 3, 0.0), + (0.45, 0.0), + (0.50, 0.0), + (0.55, 0.0), + (0.60, 0.0), + (2 / 3, 0.0), + (0.85, 0.55), + (1.00, 1.0), +] + +# 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 +_REV_COUL_KAPPA = [ + (0.00, 1.0), + (0.15, 0.55), + (1 / 3, 0.0), + (0.45, 0.0), + (0.50, 0.0), + (0.60, 0.0), + (2 / 3, 0.0), + (0.85, 0.0), + (1.00, 0.0), +] + + +@pytest.mark.parametrize("lam,expected_kappa,expected_alpha", _FWD_KAPPA_ALPHA) +def test_ring_break_morph_schedule(lam, expected_kappa, expected_alpha): + """ + ring_break_morph() produces the correct ring-break kappa and alpha at each λ. -# ring_break_morph: kappa=0 throughout potential_swap (λ<1/3) and -# restraints_off (λ<1/2); kappa ramps 0→1 through ring_open (1/2→2/3); -# kappa=1 fixed throughout morph (λ>2/3). + Uses lambdalever's initial/final values (kappa: 0→1, alpha: 1→0) to ensure + the test matches production behaviour exactly. + """ + from somd2._utils._schedules import ring_break_morph + + s = ring_break_morph() + kappa = s.morph("ring-break", "kappa", 0.0, 1.0, lam) + alpha = s.morph("ring-break", "alpha", 1.0, 0.0, lam) + assert abs(kappa - expected_kappa) < 1e-10, ( + f"ring-break kappa={kappa:.8f} at λ={lam:.4f}, expected {expected_kappa}" + ) + assert abs(alpha - expected_alpha) < 1e-10, ( + f"ring-break alpha={alpha:.8f} at λ={lam:.4f}, expected {expected_alpha}" + ) -@pytest.mark.parametrize("lam", [0.0, 1 / 3, 0.5]) -def test_ring_break_inactive_before_ring_open(forward_dynamics, lam): +@pytest.mark.parametrize("lam,expected_kappa,expected_alpha", _REV_KAPPA_ALPHA) +def test_reverse_ring_break_morph_schedule(lam, expected_kappa, expected_alpha): """ - Ring-break energy is near-zero (kappa=0) in potential_swap and - restraints_off stages and at the start of ring_open. + reverse_ring_break_morph() produces the correct ring-make kappa and alpha at each λ. + + Uses lambdalever's initial/final values (kappa: 1→0, alpha: 0→1) to ensure + the test matches production behaviour exactly. """ - e = _force_energy_kcal(forward_dynamics, lam, "ring-break") - assert abs(e) < _INACTIVE_THRESHOLD, ( - f"ring-break energy {e:.4f} kcal/mol at λ={lam:.4f} exceeds inactive " - f"threshold {_INACTIVE_THRESHOLD} kcal/mol (kappa should be 0)" + from somd2._utils._schedules import reverse_ring_break_morph + + s = reverse_ring_break_morph() + kappa = s.morph("ring-make", "kappa", 1.0, 0.0, lam) + alpha = s.morph("ring-make", "alpha", 0.0, 1.0, lam) + assert abs(kappa - expected_kappa) < 1e-10, ( + f"ring-make kappa={kappa:.8f} at λ={lam:.4f}, expected {expected_kappa}" + ) + assert abs(alpha - expected_alpha) < 1e-10, ( + f"ring-make alpha={alpha:.8f} at λ={lam:.4f}, expected {expected_alpha}" ) -@pytest.mark.parametrize("lam", [2 / 3, 1.0]) -def test_ring_break_active_after_ring_open(forward_dynamics, lam): +@pytest.mark.parametrize("lam,expected_coul_kappa", _FWD_COUL_KAPPA) +def test_ring_break_morph_coul_kappa(lam, expected_coul_kappa): """ - Ring-break energy is clearly non-zero (kappa=1) at the end of ring_open - and throughout morph. + ring_break_morph() pins coul_kappa=0 through all pre-morph stages and ramps + it 0→1 during morph only, so Coulomb only activates once atoms are separated. + + Uses lambdalever's initial/final values (coul_kappa: 0→1). """ - e = _force_energy_kcal(forward_dynamics, lam, "ring-break") - assert abs(e) > _ACTIVE_THRESHOLD, ( - f"ring-break energy {e:.4f} kcal/mol at λ={lam:.4f} is below active " - f"threshold {_ACTIVE_THRESHOLD} kcal/mol (kappa should be 1)" + from somd2._utils._schedules import ring_break_morph + + s = ring_break_morph() + coul_kappa = s.morph("ring-break", "coul_kappa", 0.0, 1.0, lam) + assert abs(coul_kappa - expected_coul_kappa) < 1e-10, ( + f"ring-break coul_kappa={coul_kappa:.8f} at λ={lam:.4f}, " + f"expected {expected_coul_kappa}" ) -def test_ring_break_energy_increases_with_kappa(forward_dynamics): +@pytest.mark.parametrize("lam,expected_coul_kappa", _REV_COUL_KAPPA) +def test_reverse_ring_break_morph_coul_kappa(lam, expected_coul_kappa): """ - Energy at kappa=1 (λ=1, morph) substantially exceeds energy at kappa=0 (λ=0). + reverse_ring_break_morph() ramps ring-make coul_kappa 1→0 through the + reversed morph stage and pins it to 0 in all subsequent stages. + + Uses lambdalever's initial/final values (coul_kappa: 1→0). """ - e0 = _force_energy_kcal(forward_dynamics, 0.0, "ring-break") - e1 = _force_energy_kcal(forward_dynamics, 1.0, "ring-break") - assert abs(e1) > abs(e0) + 0.5, ( - f"ring-break energy at λ=1 ({e1:.4f} kcal/mol) should be much larger " - f"than at λ=0 ({e0:.4f} kcal/mol)" + from somd2._utils._schedules import reverse_ring_break_morph + + s = reverse_ring_break_morph() + coul_kappa = s.morph("ring-make", "coul_kappa", 1.0, 0.0, lam) + assert abs(coul_kappa - expected_coul_kappa) < 1e-10, ( + f"ring-make coul_kappa={coul_kappa:.8f} at λ={lam:.4f}, " + f"expected {expected_coul_kappa}" ) -# ── reverse schedule energy tests ───────────────────────────────────────────── +# ── energy magnitude tests ──────────────────────────────────────────────────── -# reverse_ring_break_morph: ring-make kappa=1 fixed in morph (λ<1/3); -# kappa ramps 1→0 in ring_close (1/3→1/2); kappa=0 in bonded_perturb and -# potential_swap (λ>1/2). + +@pytest.mark.parametrize("lam", [2 / 3, 1.0]) +def test_ring_break_active_after_ring_open(forward_dynamics, lam): + """ + Ring-break energy is clearly non-zero (kappa=1) at the end of ring_open + and throughout morph. + """ + e = _force_energy_kcal(forward_dynamics, lam, "ring-break") + assert abs(e) > _ACTIVE_THRESHOLD, ( + f"ring-break energy {e:.4f} kcal/mol at λ={lam:.4f} is below active " + f"threshold {_ACTIVE_THRESHOLD} kcal/mol (kappa should be 1)" + ) def test_ring_make_active_at_lambda_zero(reverse_dynamics): @@ -219,14 +317,83 @@ def test_ring_make_active_at_lambda_zero(reverse_dynamics): ) -@pytest.mark.parametrize("lam", [2 / 3, 1.0]) -def test_ring_make_inactive_after_ring_close(reverse_dynamics, lam): +def test_ring_make_inactive_at_lambda_one(reverse_dynamics): + """ + Ring-make energy is near-zero at λ=1 (potential_swap end, kappa=0). + + At λ=1 the system is at the ring-open end state; the hard-hard correction + term in the CustomBondForce is small because the pair is at nonbonded + separation, so the absolute energy remains below the active threshold. + """ + e = _force_energy_kcal(reverse_dynamics, 1.0, "ring-make") + assert abs(e) < _ACTIVE_THRESHOLD, ( + f"ring-make energy {e:.4f} kcal/mol at λ=1 exceeds threshold " + f"{_ACTIVE_THRESHOLD} kcal/mol (kappa should be 0)" + ) + + +# ── energy symmetry tests ───────────────────────────────────────────────────── +# +# The invariant ring_break_morph().reverse() == reverse_ring_break_morph() means +# that the softcore kappa/alpha values at (forward, λ) and (reverse, 1-λ) are +# equal. Both forces act on the same bond (the original ring_breaking_bond, +# which swap_end_states=True maps to ring_making_pairs), so the energies must +# also match. The hard-hard correction appears identically on both sides and +# cancels in the comparison, making this test robust to formula changes. +# +# 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.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) + + +@pytest.mark.parametrize("lam", [0.0, 0.55, 2 / 3, 0.85, 1.0]) +def test_energy_symmetry_forward_reverse(forward_dynamics, reverse_dynamics, lam): """ - Ring-make energy is near-zero (kappa=0) in bonded_perturb and - potential_swap stages of reverse_ring_break_morph. + Single-point energy symmetry: E_ring_break_forward(λ) == E_ring_make_reverse(1-λ). + + Verifies that reverse_ring_break_morph() == ring_break_morph().reverse() and + that the mirrored kappa/alpha produce identical corrections on the same bond. """ - e = _force_energy_kcal(reverse_dynamics, lam, "ring-make") - assert abs(e) < _INACTIVE_THRESHOLD, ( - f"ring-make energy {e:.4f} kcal/mol at λ={lam:.4f} exceeds inactive " - f"threshold {_INACTIVE_THRESHOLD} kcal/mol (kappa should be 0)" + e_fwd = _force_energy_kcal(forward_dynamics, lam, "ring-break") + e_rev = _force_energy_kcal(reverse_dynamics, 1.0 - lam, "ring-make") + assert abs(e_fwd - e_rev) < 1e-4, ( + f"Energy symmetry broken at λ={lam:.4f}: " + f"ring-break forward = {e_fwd:.6f} kcal/mol, " + f"ring-make reverse(1-λ={1 - lam:.4f}) = {e_rev:.6f} kcal/mol, " + f"difference = {abs(e_fwd - e_rev):.2e} kcal/mol" ) + + +def test_schedule_symmetry(): + """ + reverse_ring_break_morph() must equal ring_break_morph().reverse(). + + Checks that the simplified implementation produces identical schedules by + comparing kappa values at a dense grid of lambda points using the default + initial/final values that lambdalever passes for ring-break kappa. + """ + from somd2._utils._schedules import ring_break_morph, reverse_ring_break_morph + + fwd = ring_break_morph() + rev = reverse_ring_break_morph() + rev_via_reverse = fwd.reverse() + + test_lambdas = [i / 20 for i in range(21)] + for lam in test_lambdas: + for force, lever, init, fin in [ + ("ring-break", "kappa", 0.0, 1.0), + ("ring-break", "alpha", 1.0, 0.0), + ("ring-break", "coul_kappa", 0.0, 1.0), + ("ring-make", "kappa", 1.0, 0.0), + ("ring-make", "alpha", 0.0, 1.0), + ("ring-make", "coul_kappa", 1.0, 0.0), + ]: + v_rev = rev.morph(force, lever, init, fin, lam) + v_rev2 = rev_via_reverse.morph(force, lever, init, fin, lam) + assert abs(v_rev - v_rev2) < 1e-12, ( + f"Mismatch for {force}/{lever} at λ={lam:.2f}: " + f"reverse_ring_break_morph={v_rev}, ring_break_morph().reverse()={v_rev2}" + )