Skip to content

Barostat the reversible-scaling sweep at λp instead of the full pressure - #4

Closed
pmrv wants to merge 8 commits into
mainfrom
claude/confident-noether-83ei6e
Closed

pmrv wants to merge 8 commits into
mainfrom
claude/confident-noether-83ei6e

Conversation

@pmrv

@pmrv pmrv commented Sep 16, 2026 •

Copy link
Copy Markdown
Owner

Closes #3.

What was wrong

In mode: ts the sweep scales the potential to λU but pinned the barostat at the full target pressure p for every λ. Sampling λU at (T₀, λp) is the real system at (T₀/λ, p), so the scaled system has to be held at λp for the swept volumes to lie on the isobar that integrate_rs assumes when it adds its constant p·V term. pf = lf * pi was computed and only ever logged. At p = 0 the two targets coincide, which is why every existing golden passed.

Changes

mode: ts (_reversible_scaling_forward / _backward)

  • fix npt only ramps linearly over a run, so a new helper _rs_sweep_barostat re-issues the COM-constrained f1 fix immediately before the sweep run with iso p → lf·p (forward) and iso lf·p → p (backward), and re-attaches fix_modify f1 temp tcm, which the redefinition drops.
  • The backward leg's middle equilibration runs the lf-scaled potential, so it now sits at lf·p instead of p.
  • lambda_schedule: uniform_temperature is refused at finite pressure at input validation (p = 0 and npt: false stay allowed). λ is a hyperbola in the step under that schedule, and no LAMMPS barostat accepts a variable target, so a linear ramp cannot track λp. The schedule is deprecated upstream, so this refusal is the final behaviour rather than a stopgap.

mode: tscale and the pre-flight scan (the mirror-image slips flagged in the issue)

  • Both ramp the real thermostat with the unscaled potential, so their barostat now stays at p instead of ramping p → lf·p.
  • The tscale backward sweep held the thermostat at t0; it now ramps tf → t0 as its log message already claimed.

Tests / goldens

  • Goldens (all at p = 0) gain the re-issued fix lines and the corrected tscale backward thermostat; nothing else in the streams moved.
  • Tests at 50 kbar pin Pstart/Pstop of every NPT block in ts forward/backward (and that no run sits between the re-issued fix and the sweep), the uniform_temperature refusal, tscale and the prescan. They fail on main.

Housekeeping

  • pyscal3<4 pinned: 4.0.0 removed the System class calphy imports and broke test collection on main too.

Verification

MD, against independent references (bcc Fe SNAP, 5488 atoms, mode: ts 100 → 800 K, details in this comment):

reference unpatched patched
F(500 K, +5 GPa) −1.435703 eV/atom (mode: fe) −9.835 meV/atom off −0.151
V(500 K, +5 GPa) 9.3737 ų/atom (true NPT) −6.27 % −0.17 %
V over 100 → 800 K must expand contracted 8.1 % expands, monotonic
F(500 K, −2.5 GPa) −1.882265 eV/atom −0.591 meV/atom off −0.110

Residuals sit inside the runs' own RS dissipation (0.22 and 0.26 meV/atom, single iteration each). The LAMMPS logs confirm the sweep barostat ramps 50000 → 6250 and -25000 → -3125 bar while the equilibrations stay pinned. The fix also recovers the negative-pressure wide-window runs that previously died with Domain too large: that was the full pressure applied to a potential scaled to 0.125·U, not the λ-span itself.

An independent review with short 256-atom Cu EAM runs (+50 kbar heating, −20 kbar cooling) confirmed press tracks λp to within 1 % along the sweep, volumes match true NPT values, and the fix re-issue causes no transient beyond the one already present on main.

Not MD-tested: tscale, covered at the command-stream level.

Unit tests: tests/test_golden_streams.py and the full suite pass locally. One pre-existing failure is unrelated: test_damaged_sweep_file_does_not_abort_the_gather fails on main too under pandas 3 (NaN vs None).

Out of scope, filed separately: #5 (integrate_rs adds p·V under npt: False).

🤖 Generated with Claude Code

https://claude.ai/code/session_01JgUx1cRgji82mZBoX7UDCy

pmrv and others added 2 commits September 16, 2026 19:20
In mode ts the sweep scales the potential to λU but held the barostat at
the full target pressure p for every λ. Under the RS/NPT mapping,
sampling λU at (T0, λp) is the real system at (T0/λ, p), so the scaled
system has to be barostatted at λp for the recorded volumes to lie on the
isobar that integrate_rs assumes when it adds the constant p·V term. The
routine already computed pf = lf*p and never used it; at p = 0 the two
targets coincide, which is why the error went unnoticed. With |p| > 0 the
swept volumes were too small and the free energy off by up to ~10 meV/atom
in the reported campaign.

fix npt only ramps its target linearly over a run, so the sweep re-issues
the COM-constrained npt fix immediately before the sweep run with the
ramp p -> lf*p (forward) and lf*p -> p (backward), re-attaching the COM
temperature compute the redefinition drops. The middle equilibration of
the backward leg, which runs the lf-scaled potential, now also sits at
lf*p. The ramp matches λp exactly for the linear lambda schedule; for
uniform_temperature it only interpolates λp between the end points, so a
warning is logged at finite pressure and the docs say to use linear.

The mirror-image slips in the real-thermostat ramps are fixed too: mode
tscale and the pre-flight temperature scan ramped the barostat p -> lf*p
while ramping the real thermostat with the unscaled potential, where the
barostat has to stay at p. The tscale backward sweep also held the
thermostat at t0 instead of ramping tf -> t0 as its log claimed.

Goldens (all at p = 0) gain the re-issued fix lines and the corrected
tscale backward thermostat. New tests at 50 kbar pin the barostat targets
of every block in ts, tscale and the prescan.

Closes #3

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JgUx1cRgji82mZBoX7UDCy
pyscal3 4.0.0 (on PyPI since the last green run on main) removed the
top-level System class that calphy.input imports, so every test module
fails at collection with an unpinned install. Pin the 3.x series in
pyproject.toml and the conda environment files.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JgUx1cRgji82mZBoX7UDCy

pmrv commented Sep 16, 2026

Copy link
Copy Markdown
Owner Author

CI status on this PR, for the record:

All four test jobs were red on ea74642 because every test module failed at collection: pyscal3 4.0.0 shipped after the last green run on main and removed the top-level System that calphy.input imports. main would fail the same way today. Pushed 452f2ec pinning pyscal3<4 in pyproject.toml and the conda env files; the golden suite passes locally against the pinned install.

One more pre-existing failure will remain once collection works: tests/test_gather_results_sweep_data.py::test_damaged_sweep_file_does_not_abort_the_gather. It fails on main too under pandas 3.0.5 (what CI installs): pandas 3 infers a str dtype for the error_code column and converts the None of intact rows to NaN, so df.loc[good]["error_code"] is None is false. The code this PR touches is unrelated, so I have not folded a fix in. Proposed one-liner in gather_results (calphy/postprocessing.py:377), which preserves the documented "stays None" contract instead of loosening the test:

df = pd.DataFrame(data=datadict)
df["error_code"] = df["error_code"].astype(object).where(df["error_code"].notna(), None)

Happy to push it here or open a separate PR, whichever you prefer.


Generated by Claude Code

@pmrv
pmrv changed the base branch from main to master September 17, 2026 08:27
Comment thread calphy/phase.py Outdated
Review suggestion on the backward-sweep comment; applied to the matching
forward-sweep comment as well.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JgUx1cRgji82mZBoX7UDCy
@pmrv
pmrv changed the base branch from master to main September 17, 2026 09:13
@pmrv

pmrv commented Sep 18, 2026

Copy link
Copy Markdown
Owner Author

Verified by MD against independent references. The fix holds on both signs of pressure, and it also recovers runs that used to crash.

Tested at 452f2ec; branch head is now 1011763 and the diff between them is comment-only in phase.py (two ASCII-arrow lines, checked), so this applies to the current head.

System: bcc Fe, SNAP, 5488 atoms, 32 ranks, execution_mode: library, LAMMPS 11 Feb 2026. Two mode: ts runs sweeping 100 → 800 K, one at +5 GPa and one at −2.5 GPa, i.e. reruns of inputs the unpatched code got wrong.

The patched path is what ran

From the LAMMPS logs — the sweep barostat ramps, the equilibrations stay pinned:

run forward sweep backward sweep
+5 GPa iso 50000 → 6250 iso 6250 → 50000
−2.5 GPa iso -25000 → -3125 iso -3125 → -25000

6250 = 0.125 × 50000 = λ_min·p. The jobs also asserted calphy.__file__ pointed at the checkout before starting.

Results

References are independent of the sweep: a single-temperature mode: fe Frenkel–Ladd calculation at the same (T, p), and a mode: pscale integral from a separately validated p = 0 anchor. The two agree with each other to 0.011 meV/atom.

reference unpatched patched
F(500 K, +5 GPa) −1.435703 eV/atom −9.835 meV/atom off −0.151
V(500 K, +5 GPa) 9.3737 ų/atom (true NPT) −6.27 % −0.17 %
V over 100 → 800 K must expand contracted 8.1 % expands, monotonic
F(500 K, −2.5 GPa) −1.882265 eV/atom −0.591 meV/atom off −0.110

The volume row is the cleanest of these: it is a state function, so it needs no anchor, no alignment and no stitching convention, and mode: fe reports the true NPT equilibrium value directly.

It also retires a separate explanation

Chunk-A runs at negative pressure had been dying with Domain too large, and we had attributed that to λ-span — wide chunks scale the potential further, and the crashes correlated with λ_min.

It was the barostat. The −2.5 GPa run above uses the 100 → 800 K window that had produced no complete RS iteration unpatched; patched, it completed in 18 h. The mechanism is the same as the free-energy error: at λ_min = 0.125 the unpatched code applies the full −25000 bar to a potential scaled to 0.125·U, where this patch applies −3125. λ-span was a correlate of how much unphysical tension the cold end saw, not the cause.

Caveats on this evidence

  • Single iteration each (n_iterations: 1), chosen because the effect being measured is ~10 meV/atom and ~6 % in volume. So the residuals above, −0.151 and −0.110 meV/atom — both slightly low, same sign — sit inside the runs' own RS dissipation of 0.220 and 0.256 meV/atom. That is neither evidence of a remaining systematic nor evidence against one; it would take 3+ iterations to separate.
  • tscale is not tested here. The removal of the λp ramp and the backward-sweep thermostat fix both look right by inspection, but no MD was run on them.
  • lambda_schedule: uniform_temperature is not tested at finite pressure, so the interpolation caveat in the docstring and the runtime warning remain unexercised — our sweeps use linear.

One note on test design, if it is useful

The decisive gate was not +5 GPa, where the error is 9.8 meV/atom — a fix is hard to get wrong there. It was −2.5 GPa, where the error is 0.6 meV/atom: a patch that overshoots, or that gets the sense of λ backwards, passes the large-error test and fails the small one. Both predictions were fixed in writing before either job was submitted.

pmrv and others added 3 commits September 20, 2026 19:49
…t barostat chords

The uniform_temperature schedule makes T0/lambda linear in the MD step,
so lambda is a hyperbola in the step. No LAMMPS barostat accepts a
variable target (fix npt, nph, press/berendsen and press/langevin all
reject v_name) and fix npt only ramps linearly over one run, so the
previous commit could only interpolate lambda*p between the sweep end
points and warn: the barostat overshot lambda*p mid-sweep by 8 % for a
400 -> 700 K sweep, 33 % for 1000 -> 3000 K and 153 % for 100 -> 800 K.

Instead of approximating lambda*p by a chord, make the schedule itself
piecewise linear. _sweep_segments splits the sweep into up to
SWEEP_MAX_SEGMENTS sub-runs whose end points lie on the exact hyperbola;
_run_sweep redefines the lambda variable as ramp() between them before
each sub-run and re-issues the integrator fix with the chord
p*lambda_a -> p*lambda_b. Within a run lambda is then linear in the
step, so the barostat equals lambda*p at every step and the free-energy
integral, taken over the recorded lambda grid, is exact. The only
approximation left is in the sampling density, which deviates from
"uniform in T" by O(1/K^2), 0.2 % at K = 20 for 1000 -> 3000 K. Each
segment boundary zeroes the barostat velocity when the fix is
re-issued, so segments are kept at least ten damping times long; short
sweeps fall back to fewer segments, down to a single linear ramp.

The runtime warning is removed. The ts_uniform_temperature golden now
shows the five 1000-step segments of the 5000-step B6 sweep. New tests
check both sweep directions: contiguous lambda chain, every chord equal
to p*lambda at both ends, boundaries on the hyperbola, step counts
summing to the sweep length, and the segment-count bounds. The linear-
schedule tests additionally assert that no run sits between the
re-issued fix and the sweep run.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JgUx1cRgji82mZBoX7UDCy
…ith exact barostat chords"

This reverts commit 0793c34.
A barostatted reversible-scaling sweep has to hold the scaled system at
lambda*p, and fix npt can only ramp its target linearly over a run
(no LAMMPS barostat accepts a variable target). That equals lambda*p
for the linear schedule only; uniform_temperature makes lambda a
hyperbola in the step, so the barostat would overshoot lambda*p
mid-sweep by 8 % for 400 -> 700 K and 153 % for 100 -> 800 K.

Rather than run that knowingly wrong behind a warning, reject the
combination at input validation. p = 0 (where the ramp is 0 -> 0) and
npt: false stay allowed. A piecewise-linear segmented schedule that
makes the combination exact follows in a separate change.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JgUx1cRgji82mZBoX7UDCy

@pmrv pmrv left a comment

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

revise

Comment thread calphy/phase.py
Comment thread calphy/phase.py Outdated
Comment thread calphy/phase.py Outdated
…dule

Review: comments that point at each other across the file are unreadable.
Every comment this PR added now states its own reason in place, without
"see X" pointers. The tscale comment takes the reviewer's wording. The
reversible_scaling docstring gains a flowchart of the thermostat and
barostat targets through one cycle, and the helper docstring explains
why fixedpoint stays at the centre of mass.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JgUx1cRgji82mZBoX7UDCy

@pmrv pmrv left a comment

Copy link
Copy Markdown
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

revise

Comment thread calphy/phase.py Outdated
Comment thread calphy/phase.py Outdated
Review: neither belongs in the code. The fixedpoint keyword stays,
matching the other npt fixes in the sweep routines; the pressure
schedule chart lives in the PR discussion.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JgUx1cRgji82mZBoX7UDCy

pmrv commented Oct 7, 2026

Copy link
Copy Markdown
Owner Author

Closing: superseded by upstream. ICAMS/calphy PR ICAMS#283 (merged 2026-09-29, in main from 2.2.1) contains the same fix: the sweep barostat ramps λ·P in both directions via _install_sweep_barostat, the middle equilibration sits at lf·P, tscale and the pre-flight scan hold the real pressure with the backward thermostat ramping Tf → T0, and tests/test_sweep_barostat_targets.py pins the targets at 50 kbar. lambda_schedule was removed upstream in the same series. The MD verification above applies to that implementation as well, since the emitted command streams are identical apart from an explicit unfix f1 before the ramping fix.

Not covered upstream and still open: #5.


Generated by Claude Code

@pmrv pmrv closed this Oct 7, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

mode: ts — RS sweep pins the barostat at the full pressure instead of λp (the computed pf is dead)

2 participants