diff --git a/src/somd2/_utils/_schedules.py b/src/somd2/_utils/_schedules.py index 1c4e48d..06985eb 100644 --- a/src/somd2/_utils/_schedules.py +++ b/src/somd2/_utils/_schedules.py @@ -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 ------- @@ -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()) @@ -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( @@ -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()) @@ -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() ) @@ -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() ) @@ -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: diff --git a/tests/schedules/test_ring_break.py b/tests/schedules/test_ring_break.py index bf6cac5..bfc0146 100644 --- a/tests/schedules/test_ring_break.py +++ b/tests/schedules/test_ring_break.py @@ -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), @@ -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), @@ -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, ( @@ -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)