From 44bd395ea50dcb9d81c3c5c561d3ce2cab6debf1 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 8 May 2026 10:27:34 +0100 Subject: [PATCH 1/4] Add ring-break and ring-make lever equations to ring-break schedules. --- src/somd2/_utils/_schedules.py | 41 ++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/somd2/_utils/_schedules.py b/src/somd2/_utils/_schedules.py index 0ba0801..04d42f2 100644 --- a/src/somd2/_utils/_schedules.py +++ b/src/somd2/_utils/_schedules.py @@ -216,6 +216,22 @@ def ring_break_morph(): s.set_equation(stage="morph", lever="torsion_k", equation=s.final()) s.set_equation(stage="morph", lever="torsion_phase", equation=s.final()) + # Ring-breaking bonds: softcore interaction grows from zero as the ring + # opens during the morph stage (alpha: 1→0, kappa: 0→1). + # Ring-making bonds: softcore interaction shrinks to zero as the ring + # closes during the morph stage (alpha: 0→1, kappa: 1→0). + # potential_swap and restraints_off stages use default (s.initial()): + # ring-break alpha=1.0/kappa=0.0 (no interaction, ring still bonded) + # ring-make alpha=0.0/kappa=1.0 (full interaction, ring not yet formed) + s.set_equation( + stage="morph", force="ring-break", lever="alpha", equation=1 - s.lam() + ) + s.set_equation(stage="morph", force="ring-break", lever="kappa", equation=s.lam()) + s.set_equation(stage="morph", force="ring-make", lever="alpha", equation=s.lam()) + s.set_equation( + stage="morph", force="ring-make", lever="kappa", equation=1 - s.lam() + ) + return s @@ -288,4 +304,29 @@ def reverse_ring_break_morph(): s.set_equation(stage="potential_swap", lever="torsion_k", equation=s.final()) s.set_equation(stage="potential_swap", lever="torsion_phase", equation=s.final()) + # morph stage (first): nonbonded-only changes; ring bonds still intact/absent. + # ring-break: alpha=1.0/kappa=0.0 throughout (no interaction, ring bonded at λ=0) + # ring-make: alpha=0.0/kappa=1.0 throughout (full interaction, ring absent at λ=0) + # bonded_perturb stage (second): ring bonds established/dissolved via Morse. + # ring-make softcore turns off as ring forms (alpha: 0→1, kappa: 1→0) + # ring-break softcore turns on as ring opens (alpha: 1→0, kappa: 0→1) + # potential_swap stage (last): Morse→harmonic swap; ring fully transitioned. + # defaults (s.final()) give ring-break alpha=0/kappa=1, ring-make alpha=1/kappa=0. + 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) + s.set_equation( + stage="bonded_perturb", force="ring-break", lever="alpha", equation=1 - s.lam() + ) + s.set_equation( + stage="bonded_perturb", force="ring-break", lever="kappa", equation=s.lam() + ) + s.set_equation( + stage="bonded_perturb", force="ring-make", lever="alpha", equation=s.lam() + ) + s.set_equation( + stage="bonded_perturb", force="ring-make", lever="kappa", equation=1 - s.lam() + ) + return s From 963eacb53bccddbfd6361c8a03c83da3dc0e01c5 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 8 May 2026 12:17:34 +0100 Subject: [PATCH 2/4] Add ring_open/ring_close stages to ring-break schedules. --- src/somd2/_utils/_schedules.py | 116 +++++++++++++++++++++------------ 1 file changed, 73 insertions(+), 43 deletions(-) diff --git a/src/somd2/_utils/_schedules.py b/src/somd2/_utils/_schedules.py index 04d42f2..bae6a68 100644 --- a/src/somd2/_utils/_schedules.py +++ b/src/somd2/_utils/_schedules.py @@ -151,7 +151,7 @@ def ring_break_morph(): """ Build a lambda schedule for ring-breaking perturbations. - Three stages: potential_swap → restraints_off → morph. + Four stages: potential_swap → restraints_off → ring_open → morph. Returns ------- @@ -163,6 +163,33 @@ def ring_break_morph(): s = _LambdaSchedule.standard_morph() + # 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. + s.prepend_stage("ring_open", s.initial()) + 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() + ) + s.set_equation( + stage="ring_open", force="ring-make", lever="alpha", equation=s.lam() + ) + s.set_equation( + stage="ring_open", force="ring-make", lever="kappa", equation=1 - s.lam() + ) + 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) @@ -207,6 +234,8 @@ 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. 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()) @@ -215,22 +244,10 @@ def ring_break_morph(): 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()) - - # Ring-breaking bonds: softcore interaction grows from zero as the ring - # opens during the morph stage (alpha: 1→0, kappa: 0→1). - # Ring-making bonds: softcore interaction shrinks to zero as the ring - # closes during the morph stage (alpha: 0→1, kappa: 1→0). - # potential_swap and restraints_off stages use default (s.initial()): - # ring-break alpha=1.0/kappa=0.0 (no interaction, ring still bonded) - # ring-make alpha=0.0/kappa=1.0 (full interaction, ring not yet formed) - s.set_equation( - stage="morph", force="ring-break", lever="alpha", equation=1 - s.lam() - ) - s.set_equation(stage="morph", force="ring-break", lever="kappa", equation=s.lam()) - s.set_equation(stage="morph", force="ring-make", lever="alpha", equation=s.lam()) - s.set_equation( - stage="morph", force="ring-make", lever="kappa", equation=1 - s.lam() - ) + 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 @@ -239,7 +256,7 @@ def reverse_ring_break_morph(): """ Build a lambda schedule for reverse ring-breaking perturbations. - Three stages: morph → bonded_perturb → potential_swap. + Four stages: morph → ring_close → bonded_perturb → potential_swap. Returns ------- @@ -251,6 +268,8 @@ def reverse_ring_break_morph(): s = _LambdaSchedule.standard_morph() + # 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()) @@ -259,7 +278,35 @@ def reverse_ring_break_morph(): 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()) + 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) + s.set_equation( + stage="ring_close", force="ring-make", lever="alpha", equation=s.lam() + ) + s.set_equation( + stage="ring_close", force="ring-make", lever="kappa", equation=1 - s.lam() + ) + # 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()) s.set_equation(stage="bonded_perturb", lever="morse_soft", equation=0 + s.lam()) s.set_equation(stage="bonded_perturb", lever="morse_hard", equation=0) @@ -285,6 +332,14 @@ def reverse_ring_break_morph(): 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() + ) + s.set_equation( + stage="bonded_perturb", force="ring-break", lever="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()) s.set_equation(stage="potential_swap", lever="morse_hard", equation=0 + s.lam()) @@ -304,29 +359,4 @@ def reverse_ring_break_morph(): s.set_equation(stage="potential_swap", lever="torsion_k", equation=s.final()) s.set_equation(stage="potential_swap", lever="torsion_phase", equation=s.final()) - # morph stage (first): nonbonded-only changes; ring bonds still intact/absent. - # ring-break: alpha=1.0/kappa=0.0 throughout (no interaction, ring bonded at λ=0) - # ring-make: alpha=0.0/kappa=1.0 throughout (full interaction, ring absent at λ=0) - # bonded_perturb stage (second): ring bonds established/dissolved via Morse. - # ring-make softcore turns off as ring forms (alpha: 0→1, kappa: 1→0) - # ring-break softcore turns on as ring opens (alpha: 1→0, kappa: 0→1) - # potential_swap stage (last): Morse→harmonic swap; ring fully transitioned. - # defaults (s.final()) give ring-break alpha=0/kappa=1, ring-make alpha=1/kappa=0. - 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) - s.set_equation( - stage="bonded_perturb", force="ring-break", lever="alpha", equation=1 - s.lam() - ) - s.set_equation( - stage="bonded_perturb", force="ring-break", lever="kappa", equation=s.lam() - ) - s.set_equation( - stage="bonded_perturb", force="ring-make", lever="alpha", equation=s.lam() - ) - s.set_equation( - stage="bonded_perturb", force="ring-make", lever="kappa", equation=1 - s.lam() - ) - return s From 18b66e189a67ab229b457828f7b7351069b6dea2 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 8 May 2026 13:36:06 +0100 Subject: [PATCH 3/4] Try weighting ring-breaking schedule stages. --- src/somd2/_utils/_schedules.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/somd2/_utils/_schedules.py b/src/somd2/_utils/_schedules.py index bae6a68..5207a84 100644 --- a/src/somd2/_utils/_schedules.py +++ b/src/somd2/_utils/_schedules.py @@ -162,13 +162,14 @@ 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. - s.prepend_stage("ring_open", s.initial()) + 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()) @@ -190,7 +191,7 @@ def ring_break_morph(): stage="ring_open", force="ring-make", lever="kappa", equation=1 - s.lam() ) - s.prepend_stage("restraints_off", s.initial()) + s.prepend_stage("restraints_off", s.initial(), weight=1) 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()) @@ -216,7 +217,7 @@ def ring_break_morph(): equation=(1 - s.lam()) * s.initial() + s.lam() * s.final(), ) - s.prepend_stage("potential_swap", s.initial()) + s.prepend_stage("potential_swap", s.initial(), weight=2) 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( @@ -267,6 +268,7 @@ def reverse_ring_break_morph(): 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). @@ -286,7 +288,7 @@ def reverse_ring_break_morph(): # 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()) + 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()) @@ -307,7 +309,7 @@ def reverse_ring_break_morph(): # 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()) + 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()) @@ -341,7 +343,7 @@ def reverse_ring_break_morph(): 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()) + 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( From 78a311de2bf4df50d8a5e4d995dce12d4e1b549e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 12 May 2026 16:09:22 +0100 Subject: [PATCH 4/4] Add unit test for ring-break/make forces. --- tests/conftest.py | 12 ++ tests/schedules/test_ring_break.py | 232 +++++++++++++++++++++++++++++ 2 files changed, 244 insertions(+) create mode 100644 tests/schedules/test_ring_break.py diff --git a/tests/conftest.py b/tests/conftest.py index db0a348..a4ee48c 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -86,3 +86,15 @@ def pert_rev_mols(): mols = sr.load_test_files("somd1_backward.prm7", "somd1_backward.rst7") pert_file = str(Path(__file__).parent / "inputs" / "backward.pert") return apply_pert(mols, pert_file) + + +@pytest.fixture(scope="session") +def syk_ring_break_mols(): + """ + Load the SYK 5035→5033 ring-breaking perturbation system. + + Reference state (λ=0): SYK-5035 with an intact ring containing a + breaking bond. Perturbed state (λ=1): SYK-5033, the open-chain analogue. + """ + mols = sr.load_test_files("syk_5035_5033.s3") + return sr.morph.link_to_reference(mols) diff --git a/tests/schedules/test_ring_break.py b/tests/schedules/test_ring_break.py new file mode 100644 index 0000000..6f589f2 --- /dev/null +++ b/tests/schedules/test_ring_break.py @@ -0,0 +1,232 @@ +""" +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 + +pytestmark = pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + 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 ────────────────────────────────────────────────────────────────── + + +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 + + mols = mols.clone() + + hard_restraints, mols = sr.restraints.morse_potential( + mols, + de="150 kcal/mol", + auto_parametrise=True, + direct_morse_replacement=True, + name="morse_hard", + ) + soft_restraints, _ = sr.restraints.morse_potential( + mols, + atoms0=hard_restraints[0].atom0(), + atoms1=hard_restraints[0].atom1(), + r0=hard_restraints[0].r0(), + k="125 kcal mol-1 A-2", + auto_parametrise=False, + de="50 kcal mol-1", + name="morse_soft", + ) + mols = make_compatible(mols) + + return mols.dynamics( + constraint="h_bonds", + perturbable_constraint="h_bonds_not_heavy_perturbed", + cutoff="10A", + cutoff_type="rf", + dynamic_constraints=True, + include_constrained_energies=False, + swap_end_states=swap_end_states, + map={ + "ghosts_are_light": True, + "check_for_h_by_max_mass": True, + "check_for_h_by_mass": False, + "check_for_h_by_element": False, + "check_for_h_by_ambertype": False, + "fix_perturbable_zero_sigmas": True, + "lambda_schedule": schedule, + "restraints": [hard_restraints, soft_restraints], + }, + ) + + +def _force_energy_kcal(d, lam, force_name): + """Return the energy (kcal/mol) for *force_name* at *lam*.""" + import openmm + + context = d.context() + d.set_lambda(lam, update_constraints=True) + grp = context._force_group_map[force_name] + state = context.getState(getEnergy=True, groups=(1 << grp)) + 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.""" + from somd2._utils._schedules import ring_break_morph + + return _build_dynamics( + syk_ring_break_mols, ring_break_morph(), swap_end_states=False + ) + + +@pytest.fixture(scope="module") +def reverse_dynamics(syk_ring_break_mols): + """Reverse ring-making dynamics: swap_end_states=True, reverse_ring_break_morph.""" + from somd2._utils._schedules import reverse_ring_break_morph + + return _build_dynamics( + syk_ring_break_mols, reverse_ring_break_morph(), swap_end_states=True + ) + + +# ── force-presence tests ────────────────────────────────────────────────────── + + +def test_forward_has_ring_break_not_ring_make(forward_dynamics): + """ + Forward direction: ring-break CustomBondForce is registered; ring-make is absent. + """ + fmap = forward_dynamics.context()._force_group_map + assert "ring-break" in fmap, "ring-break force group missing for forward direction" + assert "ring-make" not in fmap, ( + "ring-make force group should not exist for forward direction" + ) + + +def test_reverse_has_ring_make_not_ring_break(reverse_dynamics): + """ + Reverse direction (swap_end_states=True): ring-make CustomBondForce is + registered; ring-break is absent. + """ + fmap = reverse_dynamics.context()._force_group_map + assert "ring-make" in fmap, "ring-make force group missing for reverse direction" + assert "ring-break" not in fmap, ( + "ring-break force group should not exist for reverse direction" + ) + + +# ── forward schedule energy tests ───────────────────────────────────────────── + +# 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). + + +@pytest.mark.parametrize("lam", [0.0, 1 / 3, 0.5]) +def test_ring_break_inactive_before_ring_open(forward_dynamics, lam): + """ + Ring-break energy is near-zero (kappa=0) in potential_swap and + restraints_off stages and at the start of ring_open. + """ + 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)" + ) + + +@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_break_energy_increases_with_kappa(forward_dynamics): + """ + Energy at kappa=1 (λ=1, morph) substantially exceeds energy at kappa=0 (λ=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)" + ) + + +# ── reverse schedule energy 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). + + +def test_ring_make_active_at_lambda_zero(reverse_dynamics): + """ + Ring-make energy is non-zero at λ=0: the morph stage fixes kappa=1 + so the ring-make interaction is fully on from the start. + """ + e = _force_energy_kcal(reverse_dynamics, 0.0, "ring-make") + assert abs(e) > _ACTIVE_THRESHOLD, ( + f"ring-make energy {e:.4f} kcal/mol at λ=0 is below active threshold " + f"{_ACTIVE_THRESHOLD} kcal/mol (kappa should be 1 in morph stage)" + ) + + +@pytest.mark.parametrize("lam", [2 / 3, 1.0]) +def test_ring_make_inactive_after_ring_close(reverse_dynamics, lam): + """ + Ring-make energy is near-zero (kappa=0) in bonded_perturb and + potential_swap stages of reverse_ring_break_morph. + """ + 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)" + )