hv-fluctuation-os

Jarzynski's equality and Crooks' theorem applied to Langevin dynamics of canonical toy systems. NumPy-only engine. Any potential with a tunable control parameter λ can be plugged in.

Keywords: fluctuation theorems, Jarzynski, Crooks, Langevin, nonequilibrium thermodynamics, numpy-only.

Size: ~22 KB source, no weights. Runtime: ~40 s for the full benchmark on CPU. Dependencies: NumPy only.

The two theorems

Jarzynski's equality (1997). For a system driven from equilibrium at control-parameter value λ_i to λ_f over a finite time:

⟨exp(−β·W)⟩  =  exp(−β·ΔF)

where W is the work along one trajectory and ΔF is the equilibrium free-energy difference between λ_f and λ_i. The average is over many nonequilibrium trajectories.

Crooks' fluctuation theorem (1999). For forward and reverse driving:

P_F(W) / P_R(−W)  =  exp(β·(W − ΔF))

Taking logs:

log P_F(W) − log P_R(−W)  =  β·W − β·ΔF

A plot of the log-ratio vs W is a straight line with slope β and intercept −β·ΔF. The intercept gives ΔF directly.

The engine

component	role
System	potential U(x, λ), force, ∂U/∂λ, free energy
LinearSchedule	λ(t) linear in time over [0, T]
SinusoidalSchedule	λ(t) sinusoidal over the first half-period
sample_equilibrium	draw samples from the Boltzmann distribution
sample_forward	Euler–Maruyama Langevin; returns work W per trajectory
sample_forward_both	forward and reverse in one call
jarzynski_estimate	naive and bias-corrected ΔF estimates
crooks_test	fit of log-ratio vs W; returns slope and ΔF

Headline numbers

Harmonic trap (ΔF = 0 exactly, work dissipative):

n_traj	βΔF_naive	βΔF_bc
100	+0.114	+0.121
500	−0.153	−0.149
2000	−0.074	−0.073
10000	−0.021	−0.021

Converges to 0 as expected. The work itself has ⟨W⟩ = 0.64, σ_W = 1.14 — all dissipation.

Tunable quadratic (ΔF = ln(4)/2 ≈ 0.6931):

n_traj	βΔF_naive	⟨W⟩	exact
100	+0.676	0.733	0.693
1000	+0.717	0.781	0.693
10000	+0.696	0.759	0.693

Second law: ⟨W⟩ ≥ ΔF everywhere. Jarzynski recovers ΔF to 0.4% at n=10000.

Crooks test (n=5000 forward and reverse):

quantity	fit	exact
slope	+0.911	+1.000
ΔF from intercept	+0.658	+0.693

ΔF recovered to within 0.036 absolute error, R² = 0.82 for the linear fit.

Convergence scaling (std over 10 seeds):

n_traj	std	1/√N
100	0.038	0.100
500	0.013	0.045
2000	0.008	0.022
10000	0.002	0.010

Faster than 1/√N. The observed decay is closer to 1/N at large n, likely because the exponential average is dominated by a few trajectories whose contribution becomes smoother with more samples.

Double-well tilt (λ: 0.0 → 1.5, barrier crossing):

n_traj	βΔF_naive	⟨W⟩	σ_W	exact ΔF
200	−2.326	−0.920	2.29	−2.290
1000	−2.296	−0.897	2.27	−2.290
5000	−2.306	−0.881	2.30	−2.290

The work distribution is broad (σ_W ≈ 2.3) because barrier crossings create rare low-W trajectories. The naive estimator still recovers ΔF to 0.7% because the exponential average is dominated by those rare events.

17/17 consistency checks pass.

How to use

from hv_fluctuation_os import (
    MovingHarmonicTrap, TunableQuadratic, DoubleWellTilted,
    LinearSchedule, sample_equilibrium, sample_forward,
    sample_forward_both, jarzynski_estimate, crooks_test,
)

# Pick a system
system = TunableQuadratic()

# Define the driving schedule
schedule = LinearSchedule(lam_start=1.0, lam_end=4.0, T_total=5.0)
beta = 1.0

# Sample equilibrium initial conditions
rng = np.random.default_rng(0)
x0 = sample_equilibrium(system, schedule.lam_start, beta,
                         n_traj=10000, rng=rng)

# Drive forward; get work per trajectory
W = sample_forward(system, schedule, beta, x0,
                    dt=1e-3, n_steps=5000, rng=rng)

# Estimate ΔF from Jarzynski
j = jarzynski_estimate(W, beta)
print(f"ΔF = {j['delta_F_naive']:.4f}  (bc: {j['delta_F_bc']:.4f})")

# Or use Crooks with forward + reverse
W_f, W_r = sample_forward_both(system, schedule, beta,
                                 n_traj=5000, dt=1e-3, n_steps=5000,
                                 n_equil=0, rng=rng)
crooks = crooks_test(W_f, W_r, beta, n_bins=25)
print(f"slope = {crooks['slope']:.3f}  ΔF = {crooks['delta_F']:.4f}")