physicskit.statphys

physicskit.statphys hub · History · Examples · API reference

physicskit.statphys: interactive, computational statistical mechanics.

physicskit.statphys brings the core breakthroughs of statistical mechanics, lattice physics, critical phenomena, molecular dynamics, and disordered systems into a visual, computational Python framework:

• Lattice models (physicskit.statphys.chapters.ising_lattice): the 2D Ising model with Metropolis, Wolff and Swendsen-Wang cluster dynamics, the q-state Potts model, and the XY model’s topological Kosterlitz-Thouless transition.

• Disordered systems (physicskit.statphys.chapters.spin_glass): the Edwards-Anderson Ising spin glass, with quenched bond disorder, frustration, and the replica-overlap order parameter.

• Molecular dynamics (physicskit.statphys.chapters.molecular_dynamics): a Velocity-Verlet-integrated Lennard-Jones gas that relaxes toward the Maxwell-Boltzmann distribution, with live tracking of Boltzmann’s H-function.

• Diffusion (physicskit.statphys.chapters.random_walk): random walk ensembles illustrating the Einstein relation and the central limit theorem.

• Stochastic dynamics (physicskit.statphys.chapters.langevin, physicskit.statphys.utils.fokker_planck): overdamped Brownian motion and underdamped Langevin particles checking the Einstein relation \(D = k_BT/\gamma\), the Ornstein-Uhlenbeck process, and a 1D Fokker-Planck solver for the probability density itself.

• Critical phenomena (physicskit.statphys.chapters.percolation, physicskit.statphys.chapters.renormalization, physicskit.statphys.chapters.sandpile): site/bond percolation with Hoshen-Kopelman cluster labeling, Kadanoff block-spin renormalization group coarse-graining, and the self-organized-critical Bak-Tang-Wiesenfeld sandpile.

• Irreversibility (physicskit.statphys.chapters.ehrenfest_urn): the Ehrenfest urn model, the simplest system reconciling reversible microscopic dynamics with an observed thermodynamic arrow of time.

• Replica symmetry breaking (physicskit.statphys.chapters.spin_glass): the infinite-range Sherrington-Kirkpatrick spin glass and its numerically observable replica-overlap distribution.

• Interface growth (physicskit.statphys.chapters.kpz_growth): the Kardar-Parisi-Zhang universality class via Restricted Solid-On-Solid deposition.

• Nonequilibrium work (physicskit.statphys.chapters.nonequilibrium_work): Jarzynski’s equality, recovering exact equilibrium free energies from irreversible, finite-speed work measurements.

• Flat-histogram sampling (physicskit.statphys.chapters.wang_landau): Wang-Landau estimation of the Ising density of states, giving the thermodynamics at every temperature from one run.

• Quantum Monte Carlo (physicskit.statphys.chapters.path_integral): path-integral Monte Carlo of a particle in a harmonic or anharmonic well, as a classical ring polymer in imaginary time.

Examples

>>> from physicskit.statphys.chapters import Ising2D
>>> model = Ising2D(L=32, seed=0)
>>> round(float(model.T_C), 3)
2.269
>>> model.sweep(beta=1.0 / model.T_C, n_sweeps=100)

class physicskit.statphys.BTWSandpile(L=64, threshold=4, seed=None)

Bases: object

The Bak-Tang-Wiesenfeld sandpile on an L x L grid with open (absorbing) boundaries.

Grains are added one at a time to a chosen (or random) site; whenever a site’s height reaches the toppling threshold, it topples, distributing one grain to each orthogonal neighbor (grains that fall off the grid’s edge are permanently lost). A single added grain can trigger a chain reaction of topplings – an avalanche – whose size has no characteristic scale once the pile has self-organized to its critical state.

Parameters:

• L (int, default=64) – Linear grid size.

• threshold (int, default=4) – Toppling threshold (4 is the standard value on the square lattice).

• seed (int, optional) – Seed for the random drop-site sequence.

Variables:

heights (ndarray of shape (L, L)) – Current pile heights; always < threshold everywhere immediately after a call to add_grain() or run().

add_grain(site=None)

Add one grain and relax the pile to stability.

Parameters:

site (tuple of int, optional) – (i, j) grid coordinates to drop the grain at. A uniformly random site is chosen if not given.

Returns:

int – Avalanche size: the number of individual topplings triggered by this grain (0 if it did not trigger any toppling).

run(n_grains, warmup=0)

Drop many grains at random sites, recording each one’s avalanche size.

Parameters:

• n_grains (int) – Number of grains to add and record after the warmup period.

• warmup (int, default=0) – Grains added first, to let the pile reach its self-organized critical state, whose avalanche sizes are discarded.

Returns:

ndarray of shape (n_grains,), dtype int64 – Avalanche size triggered by each recorded grain.

Examples

>>> pile = BTWSandpile(L=20, seed=0)
>>> sizes = pile.run(n_grains=500, warmup=500)
>>> sizes.shape
(500,)

total_grains()

Total number of grains currently on the pile.

class physicskit.statphys.BlockSpinRG(L=64, T=None, J=1.0, n_equil_sweeps=500, seed=None)

Bases: object

Iterated 2x2 majority-rule block-spin coarse-graining of an Ising configuration.

Parameters:

• L (int, default=64) – Linear size of the initial lattice. Should be a power of 2 times the largest block factor you intend to reach, e.g. L=64 supports up to 6 coarse-graining steps.

• T (float, default=None) – Temperature at which to equilibrate the initial configuration, in units of J/kB. Defaults to the exact critical temperature (Ising2D.T_C) if not given, i.e. the critical fixed point.

• J (float, default=1.0) – Ising coupling constant used to equilibrate the initial configuration.

• n_equil_sweeps (int, default=500) – Wolff cluster updates used to equilibrate the initial fine-grained configuration at temperature T before coarse-graining begins.

• seed (int, optional) – Seed for the initial equilibration and for majority-rule tie breaking.

Variables:

initial_spins (ndarray of shape (L, L)) – The equilibrated fine-grained configuration coarse-graining starts from.

coarse_grain_step(spins)

Apply one 2x2 majority-rule block-spin transformation.

Each non-overlapping \(2 \times 2\) block of the input is replaced by a single spin equal to the sign of the block’s sum; an exact tie (sum of zero) is broken by an independent fair coin flip.

Parameters:

spins (ndarray of shape (L, L)) – Input configuration. L must be even.

Returns:

ndarray of shape (L // 2, L // 2) – Coarse-grained spin configuration.

iterate(n_steps=None)

Repeatedly coarse-grain the initial configuration.

Parameters:

n_steps (int, optional) – Number of coarse-graining steps. Defaults to floor(log2(L)) - 1, i.e. as many as keep the lattice at least \(2 \times 2\).

Returns:

list of ndarray – The sequence of configurations [initial_spins, step_1, step_2, ...], of shrinking size \(L, L/2, L/4, \ldots\)

Examples

>>> rg = BlockSpinRG(L=16, T=100.0, seed=0)
>>> grids = rg.iterate(n_steps=3)
>>> [g.shape[0] for g in grids]
[16, 8, 4, 2]

static order_parameter(spins)

Magnitude of the mean magnetization per site, \(|\langle s \rangle|\).

A convenient scalar summary of a grid in the RG flow: it drifts to 1 under iteration when the starting temperature is below \(T_c\) (flow to the ordered fixed point), to 0 when above (flow to the disordered fixed point), and hovers at an intermediate, scale-invariant value at \(T_c\).

Parameters:

spins (ndarray) – A spin configuration.

Returns:

float

class physicskit.statphys.BrownianMotion(n_particles=1000, dim=1, gamma=1.0, kT=1.0, force=0.0, stiffness=0.0, seed=None)

Bases: object

Overdamped Langevin (Brownian) particles, optionally driven and trapped.

Each coordinate obeys

\[dx = \frac{F - kx}{\gamma}\,dt + \sqrt{\frac{2k_BT}{\gamma}}\,dW,\]

with constant force \(F\), trap stiffness \(k\), and friction \(\gamma\) (the inverse mobility). Free particles (\(F = k = 0\)) diffuse with \(\langle \Delta x^2\rangle = 2Dt\) per dimension and \(D = k_BT/\gamma\) (Einstein 1905); a constant force adds a drift \(F/\gamma\) without changing the spread. Euler-Maruyama is exact in distribution for \(k = 0\), since the drift is then constant.

Parameters:

• n_particles (int, default=1000) – Number of independent particles.

• dim (int, default=1) – Spatial dimension.

• gamma (float, default=1.0) – Friction coefficient.

• kT (float, default=1.0) – Temperature (\(k_B = 1\)).

• force (float or array_like of shape (dim,), default=0.0) – Constant external force.

• stiffness (float, default=0.0) – Harmonic trap stiffness \(k\) (trap centered at the origin).

• seed (int, optional) – Seed for the noise generator.

Variables:

• times (ndarray or None) – Recorded times, set by run().

• positions (ndarray of shape (n_frames + 1, n_particles, dim) or None) – Recorded positions, set by run().

Examples

>>> bm = BrownianMotion(n_particles=4000, dim=2, gamma=2.0, kT=1.0, seed=0)
>>> t, x = bm.run(t_max=5.0, dt=0.01, n_frames=50)
>>> x.shape
(51, 4000, 2)
>>> round(bm.diffusion_coefficient, 3), round(bm.measured_diffusion_coefficient(), 1)
(0.5, 0.5)

property diffusion_coefficient: float

Einstein’s \(D = k_BT/\gamma\).

mean_squared_displacement()

Ensemble mean-squared displacement from the initial positions, summed over dimensions.

Returns:

• times (ndarray)

• msd (ndarray) – \(\langle |\mathbf x(t) - \mathbf x(0)|^2\rangle\); for free particles, \(2\,d\,D t\).

measured_diffusion_coefficient()

Diffusion coefficient from the growth of the displacement variance.

Fits \(\mathrm{Var}[\Delta x] = 2Dt\) (per dimension, through the origin) so that a constant-force drift does not bias it.

Return type:

float

Returns:

float

measured_mobility()

Mobility from the mean drift velocity along the applied force.

Return type:

float

Returns:

float – \(\langle \Delta\mathbf x\rangle\cdot\hat{\mathbf F} / (|\mathbf F|\, t)\) at the last frame; nan without a force.

property mobility: float

drift velocity per unit force.

Type:

Mobility \(\mu = 1/\gamma\)

run(t_max, dt, n_frames=200, x0=None)

Integrate every particle from x0 to t_max.

Parameters:

• t_max (float) – Total time (rounded to a whole number of steps per frame).

• dt (float) – Time step.

• n_frames (int, default=200) – Number of recorded frames after the initial one.

• x0 (array_like of shape (n_particles, dim), optional) – Initial positions; defaults to the origin.

Returns:

• times (ndarray of shape (n_frames + 1,))

• positions (ndarray of shape (n_frames + 1, n_particles, dim))

class physicskit.statphys.EdwardsAndersonSpinGlass2D(L=32, J=1.0, kB=1.0, seed=None)

Bases: object

The 2D +/-J Edwards-Anderson Ising spin glass on a periodic L x L lattice.

Parameters:

• L (int, default=32) – Linear lattice size.

• J (float, default=1.0) – Bond magnitude; each bond independently takes the value +J or -J with equal probability (the “bimodal” Edwards-Anderson glass).

• kB (float, default=1.0) – Boltzmann constant.

• seed (int, optional) – Seed for the quenched disorder realization, the initial spin configuration, and the Metropolis dynamics.

Variables:

• spins (ndarray of shape (L, L)) – Current spin configuration, values in {-1, +1}.

• J_down (J_right,) – Quenched bond couplings; see physicskit.statphys.core.monte_carlo.metropolis_sweep_spin_glass() for their indexing convention. Fixed for the model’s lifetime.

edwards_anderson_order_parameter(beta, n_equil=200, n_measure=200, measure_every=1)

Estimate \(\langle q^2 \rangle\), the mean-squared Edwards-Anderson replica overlap.

Two independently initialized and thermalized replicas, sharing the same quenched bonds, are compared via their spin overlap \(q = \frac{1}{N}\sum_i s_i^{(1)} s_i^{(2)}\). Even though each replica’s spins keep fluctuating, below the spin-glass transition they fluctuate around the same disorder-selected frozen pattern, so \(\langle q^2 \rangle\) stays nonzero; above it, the replicas decorrelate and \(\langle q^2 \rangle \to 0\).

Parameters:

• beta (float) – Inverse temperature.

• n_equil (int, default=200) – Equilibration sweeps for each replica.

• n_measure (int, default=200) – Number of overlap measurements.

• measure_every (int, default=1) – Sweeps between successive measurements.

Returns:

float – Estimated \(\langle q^2 \rangle \in [0, 1]\).

energy()

Total energy \(-\sum_{\langle i,j\rangle} J_{ij} s_i s_j\) of the current configuration.

frustration_density()

Fraction of elementary plaquettes that are frustrated.

A plaquette is frustrated when the product of its four bond couplings is negative – no spin configuration can satisfy all four bonds at once. This is a fixed geometric property of the quenched disorder, independent of the current spin configuration.

Returns:

float

sweep(beta, n_sweeps=1)

Advance the spins by n_sweeps Metropolis-Hastings sweeps (bonds are fixed).

Parameters:

• beta (float) – Inverse temperature.

• n_sweeps (int, default=1) – Number of sweeps to perform.

class physicskit.statphys.EhrenfestUrn(n_balls=100, n_left_init=None, seed=None)

Bases: object

The Ehrenfest urn model of N labeled balls exchanged between two boxes.

At each time step, one of the n_balls balls (chosen uniformly at random) is moved from whichever box it currently occupies to the other.

Parameters:

• n_balls (int, default=100) – Total number of balls, N.

• n_left_init (int, optional) – Initial number of balls in the left box. Defaults to n_balls (all balls start in one box, a maximally non-equilibrium state).

• seed (int, optional) – Seed for the ball-selection sequence.

Variables:

n_left (int) – Current number of balls in the left box.

entropy()

Boltzmann entropy \(S/k_B = \ln W\) of the current macrostate.

\(W = \binom{N}{n_{\text{left}}}\) counts the ball-labelings consistent with the current left-box occupancy; the binomial coefficient is evaluated via the log-gamma function (scipy.special.gammaln) to stay numerically stable for large \(N\).

Returns:

float

run(n_steps)

Run the urn for n_steps and record its history.

Parameters:

n_steps (int) – Number of steps to simulate.

Returns:

dict of str -> ndarray – Keys "t", "n_left", "entropy", each of length n_steps + 1 (including the initial state).

Examples

>>> urn = EhrenfestUrn(n_balls=50, seed=0)
>>> history = urn.run(200)
>>> history["n_left"].shape
(201,)
>>> int(history["n_left"][0])
50

step()

Move one uniformly random ball to the other box.

Returns:

int – The updated n_left.

class physicskit.statphys.Ising2D(L=32, J=1.0, kB=1.0, seed=None)

Bases: object

The 2D ferromagnetic Ising model on a periodic \(L \times L\) lattice.

Parameters:

• L (int, default=32) – Linear lattice size.

• J (float, default=1.0) – Ferromagnetic coupling constant.

• kB (float, default=1.0) – Boltzmann constant, in whatever unit system J and T are expressed in.

• seed (int, optional) – Seed for reproducible initial spin configurations and Metropolis dynamics.

Variables:

spins (ndarray of shape (L, L)) – Current spin configuration, values in {-1, +1}.

Examples

>>> model = Ising2D(L=16, seed=0)
>>> model.sweep(beta=1.0, n_sweeps=50)
>>> -1.0 <= model.magnetization() / model.n_sites <= 1.0
True

property T_C

Onsager’s exact critical temperature, \(2J / (k_B \ln(1 + \sqrt{2})) \approx 2.269 J / k_B\).

energy()

Total energy of the current configuration, \(-J \sum_{\langle i,j\rangle} s_i s_j\).

magnetization()

Total magnetization \(M = \sum_i s_i\) of the current configuration.

reset(ordered=False)

Reinitialize the spin configuration.

Parameters:

ordered (bool, default=False) – If True, start from a fully aligned (all +1) configuration (a low-temperature / \(T=0\) state). Otherwise draw an independent random \pm 1 configuration (a \(T=\infty\) state).

run_temperature_sweep(temperatures, n_equil=200, n_measure=500, algorithm='metropolis', measure_every=1)

Sweep over a range of temperatures, measuring equilibrium thermodynamics at each.

At every temperature the lattice is first equilibrated for n_equil sweeps, then sampled every measure_every sweeps for n_measure measurements, from which the specific heat and magnetic susceptibility are estimated via their fluctuation-dissipation relations. The final configuration at each temperature seeds the next (simulated annealing), which shortens equilibration as the sweep proceeds.

Parameters:

• temperatures (array_like) – Temperatures to sample, typically ordered from high to low to exploit annealing.

• n_equil (int, default=200) – Equilibration sweeps discarded before measuring at each temperature.

• n_measure (int, default=500) – Number of measurements collected at each temperature.

• algorithm ({"metropolis", "wolff", "swendsen-wang"}, default="metropolis") – Update rule, see sweep().

• measure_every (int, default=1) – Number of sweeps between successive measurements.

Returns:

dict of str -> ndarray – Keys "T", "E", "M", "C_v", "chi", each an array aligned with temperatures. "E" and "M" are per-site averages; "C_v" and "chi" are per-site.

Examples

>>> model = Ising2D(L=12, seed=1)
>>> result = model.run_temperature_sweep(
...     [3.5, 2.269, 1.0], n_equil=20, n_measure=20
... )
>>> sorted(result.keys())
['C_v', 'E', 'M', 'T', 'chi']

sweep(beta, algorithm='metropolis', n_sweeps=1)

Advance the lattice by one or more Monte Carlo sweeps.

Parameters:

• beta (float) – Inverse temperature \(1/(k_B T)\).

• algorithm ({"metropolis", "wolff", "swendsen-wang"}, default="metropolis") – Update rule. Metropolis flips single spins; Wolff grows and flips one cluster, and Swendsen-Wang partitions the whole lattice into clusters and flips each with probability 1/2. Both cluster rules decorrelate much faster near \(T_c\).

• n_sweeps (int, default=1) – Number of sweeps (Metropolis) or cluster updates (Wolff, Swendsen-Wang) to perform.

Returns:

int or None – For "wolff", the size of the last cluster flipped; for "swendsen-wang", the number of clusters in the last update; None for "metropolis".

class physicskit.statphys.JarzynskiHarmonicTrap(k=1.0, gamma=1.0, kB=1.0, T=1.0, seed=None)

Bases: object

A Brownian particle in a harmonic trap dragged at finite speed, verifying the Jarzynski equality.

A single overdamped particle sits in the potential \(U(x, \lambda) = \frac{1}{2} k (x - \lambda)^2\), whose center \(\lambda\) is the control parameter. run_protocol() drags \(\lambda\) linearly from one value to another over a finite time \(\tau\), for many independent realizations starting from thermal equilibrium, integrating the overdamped Langevin equation

\[\gamma\, dx = -k(x - \lambda)\, dt + \sqrt{2\gamma k_B T}\, dW_t,\]

and accumulating each trajectory’s thermodynamic work \(W = \int \frac{\partial U}{\partial \lambda}\, \dot\lambda\, dt = -\int k(x - \lambda)\, d\lambda\). Because the trap’s stiffness never changes – only its center moves – the true equilibrium free energy is identical before and after the protocol, \(\Delta F = 0\) exactly, by translational invariance. This makes the harmonic trap the cleanest possible test case: any protocol speed dissipates heat and drives \(\langle W \rangle\) strictly positive (the second law), while the exponential average of the same work samples should recover \(\Delta F = 0\) regardless of how fast or slow the drag is.

Parameters:

• k (float, default=1.0) – Trap stiffness.

• gamma (float, default=1.0) – Friction (drag) coefficient.

• kB (float, default=1.0) – Boltzmann constant.

• T (float, default=1.0) – Temperature.

• seed (int, optional) – Seed for the initial equilibrium sampling and the thermal noise.

jarzynski_free_energy_estimate(work)

Estimate \(\Delta F = -k_B T \ln \langle e^{-\beta W} \rangle\) from work samples.

Uses the log-sum-exp trick for numerical stability, since individual \(e^{-\beta W}\) terms can otherwise under- or overflow long before the average itself does.

Parameters:

work (array_like) – Work samples, e.g. from run_protocol().

Returns:

float – Estimated free energy difference. For JarzynskiHarmonicTrap, the true value is exactly 0 regardless of the protocol.

run_protocol(lambda_0=0.0, lambda_1=3.0, tau=5.0, n_steps=1000, n_trajectories=5000)

Drag the trap center from lambda_0 to lambda_1 over time tau, recording each trajectory’s work.

Every trajectory starts from an independent sample of the equilibrium Boltzmann distribution at the initial trap position, \(x_0 \sim \mathcal{N}(\lambda_0, k_B T / k)\), then is integrated forward by Euler-Maruyama while the trap center advances linearly in time.

Parameters:

• lambda_0 (float, default=0.0, 3.0) – Initial and final trap-center positions.

• lambda_1 (float, default=0.0, 3.0) – Initial and final trap-center positions.

• tau (float, default=5.0) – Protocol duration. Smaller tau drags the system further from equilibrium and dissipates more heat, but should not change the Jarzynski free-energy estimate.

• n_steps (int, default=1000) – Number of Euler-Maruyama integration steps.

• n_trajectories (int, default=5000) – Number of independent realizations. The exponential average underlying the Jarzynski estimator is dominated by rare, low-work trajectories, so it converges slowly (more trajectories are needed for a longer or more violent protocol).

Returns:

ndarray of shape (n_trajectories,) – Accumulated work \(W\) for each realization.

class physicskit.statphys.KPZInterface(L=256, seed=None)

Bases: object

A 1D Restricted Solid-On-Solid growing interface on a periodic lattice of L sites.

Parameters:

• L (int, default=256) – Number of sites (periodic boundary conditions).

• seed (int, optional) – Seed for the random deposition-site sequence.

Variables:

• heights (ndarray of shape (L,)) – Current interface height profile, initialized flat at zero.

• time (float) – Elapsed time, in units of Monte Carlo sweeps (one sweep = L attempted deposition events).

grow(n_sweeps=1)

Advance the interface by n_sweeps Monte Carlo sweeps of RSOS deposition.

Parameters:

n_sweeps (float, default=1) – Number of sweeps (L deposition attempts each) to perform.

mean_height()

Mean interface height \(\bar h\), which grows linearly in time on average.

run_growth_curve(t_max, n_points=40, log_spaced=True)

Grow the interface while recording the width at a sequence of checkpoint times.

The width is expected to grow as \(w(t) \sim t^{\beta}\) with \(\beta = 1/3\) at early times, crossing over to a length-limited saturation \(w_{\text{sat}} \sim L^{\alpha}\) with \(\alpha = 1/2\) once \(t\) exceeds the correlation time \(\sim L^{z}\) (\(z = 3/2\)).

Parameters:

• t_max (float) – Final time (in sweeps) to grow to.

• n_points (int, default=40) – Number of checkpoints to record.

• log_spaced (bool, default=True) – Use logarithmically spaced checkpoints (recommended for resolving the early-time power-law growth regime) rather than linearly spaced ones.

Returns:

• times (ndarray of shape (n_points,))

• widths (ndarray of shape (n_points,))

width()

Interface width \(w = \sqrt{\langle (h - \bar h)^2 \rangle}\), the key KPZ observable.

Returns:

float

class physicskit.statphys.LangevinDynamics(n_particles=1000, dim=1, mass=1.0, gamma=1.0, kT=1.0, force=0.0, stiffness=0.0, seed=None)

Bases: object

Underdamped Langevin particles: inertia, friction, and thermal noise.

\[dx = v\,dt, \qquad m\,dv = (F - kx - \gamma v)\,dt + \sqrt{2\gamma k_BT}\,dW\]

(Langevin 1908). The noise amplitude is tied to the friction by the fluctuation-dissipation relation, so every velocity component thermalizes to \(\langle v^2\rangle = k_BT/m\) on the time scale \(m/\gamma\), and at times long compared with it the particles diffuse with Einstein’s \(D = k_BT/\gamma\).

Integrated by default with the BAOAB splitting (baoab_integrate()), which stays at the bath temperature for weak friction. Plain Euler-Maruyama (method="euler_maruyama") heats a trapped particle of frequency \(\omega\) to roughly \(k_BT\,\gamma/(\gamma - m\omega^2\Delta t)\), a large error once \(\gamma/m\) is comparable to \(\omega^2\Delta t\).

Parameters:

• n_particles (int, default=1000) – Number of independent particles.

• dim (int, default=1) – Spatial dimension.

• mass (float, default=1.0) – Particle mass \(m\).

• gamma (float, default=1.0) – Friction coefficient.

• kT (float, default=1.0) – Bath temperature (\(k_B = 1\)).

• force (float or array_like of shape (dim,), default=0.0) – Constant external force.

• stiffness (float, default=0.0) – Harmonic trap stiffness.

• seed (int, optional) – Seed for the noise generator.

Examples

>>> ld = LangevinDynamics(n_particles=5000, mass=1.0, gamma=1.0, kT=0.5, seed=1)
>>> t, x, v = ld.run(t_max=10.0, dt=0.005, n_frames=20)
>>> round(ld.kinetic_temperature(), 1)
0.5

property diffusion_coefficient: float

Long-time Einstein diffusion coefficient \(k_BT/\gamma\) (free particles).

kinetic_temperature(frame=-1)

Kinetic temperature \(m\langle v^2\rangle\) per component at one frame.

Parameters:

frame (int) – Frame index.

Return type:

float

Returns:

float

mean_squared_displacement()

Ensemble mean-squared displacement, summed over dimensions.

Returns:

• times (ndarray)

• msd (ndarray)

run(t_max, dt, n_frames=200, x0=None, v0=None, method='baoab')

Integrate every particle to t_max.

Parameters:

• t_max (float) – Total time.

• dt (float) – Time step, which should be small compared with \(m/\gamma\).

• n_frames (int, default=200) – Number of recorded frames after the initial one.

• x0 (array_like of shape (n_particles, dim), optional) – Initial positions and velocities; default zero (at rest at the origin).

• v0 (array_like of shape (n_particles, dim), optional) – Initial positions and velocities; default zero (at rest at the origin).

• method ({"baoab", "euler_maruyama"}, default="baoab") – Integration scheme.

Returns:

• times (ndarray of shape (n_frames + 1,))

• positions, velocities (ndarray of shape (n_frames + 1, n_particles, dim))

property velocity_relaxation_time: float

Momentum relaxation time \(m/\gamma\).

class physicskit.statphys.LennardJonesGas(n_particles=100, box_size=20.0, temperature_init=1.0, mass=1.0, dt=0.005, epsilon=1.0, sigma=1.0, cutoff=2.5, initial_velocity_distribution='uniform', kB=1.0, seed=None)

Bases: object

A 2D N-particle Lennard-Jones gas in a periodic square box.

Particles are initialized on a lattice (to avoid the divergent Lennard-Jones repulsion of random overlapping placements) and given velocities drawn from an arbitrary, generally non-equilibrium initial distribution, then integrated with the symplectic Velocity Verlet scheme.

Parameters:

• n_particles (int, default=100) – Number of particles. Rounded up internally so that ceil(sqrt(n_particles))**2 particles fit on a square starting lattice; the requested count is used for all physical averages.

• box_size (float, default=20.0) – Side length of the square periodic simulation box.

• temperature_init (float, default=1.0) – Nominal temperature used to set the initial velocity scale (see initial_velocity_distribution).

• mass (float, default=1.0) – Common particle mass.

• dt (float, default=0.005) – Integration time step.

• epsilon (float, default=1.0) – Lennard-Jones well depth.

• sigma (float, default=1.0) – Lennard-Jones length scale.

• cutoff (float, default=2.5) – Interaction cutoff radius.

• initial_velocity_distribution ({"maxwell_boltzmann", "uniform", "delta"}, default="uniform") – Shape of the initial (generally non-equilibrium) velocity distribution. "uniform" and "delta" are useful starting points for watching relaxation toward Maxwell-Boltzmann via the H-theorem; "maxwell_boltzmann" starts already at equilibrium.

• kB (float, default=1.0) – Boltzmann constant.

• seed (int, optional) – Seed for the initial configuration.

Variables:

• positions (ndarray of shape (N, 2))

• velocities (ndarray of shape (N, 2))

• forces (ndarray of shape (N, 2)) – Forces at the current positions, kept up to date across calls to step() so each step costs one force evaluation.

• time (float) – Elapsed simulation time.

h_function(bins=40, v_max=None)

Estimate Boltzmann’s H-function from the instantaneous speed histogram.

\[H(t) = \int_0^\infty f(v, t) \ln \frac{f(v, t)}{2\pi v} \, dv\]

estimated by binning the current particle speeds into a normalized density histogram. For an isotropic 2D gas, \(f(v)/(2\pi v)\) is the velocity-space density, so this equals Boltzmann’s \(\int f \ln f \, d^2v\); dropping the \(2\pi v\) Jacobian would give a functional minimized by a half-Gaussian in speed rather than by the Maxwell-Boltzmann (Rayleigh) distribution. \(H\) decreases monotonically (on average) as an arbitrary initial velocity distribution relaxes toward the Maxwell-Boltzmann equilibrium, which minimizes \(H\) at fixed energy – the microscopic root of the second law of thermodynamics.

Parameters:

• bins (int, default=40) – Number of histogram bins.

• v_max (float, optional) – Upper edge of the speed histogram. Defaults to 1.5 * speeds().max().

Returns:

float – The estimated value of \(H(t)\).

kinetic_energy()

Total kinetic energy \(\sum_i \frac{1}{2} m v_i^2\).

potential_energy()

Total Lennard-Jones potential energy at the current configuration.

run(n_steps, steps_per_record=10, bins=40)

Integrate forward while recording a time series of thermodynamic observables.

Parameters:

• n_steps (int) – Total number of integration steps to perform.

• steps_per_record (int, default=10) – Number of integration steps between recorded samples.

• bins (int, default=40) – Histogram bins used for each h_function() estimate.

Returns:

dict of str -> ndarray – Keys "time", "H", "temperature", "kinetic_energy", "potential_energy", "total_energy", each a 1D array sampled every steps_per_record steps (including the initial state).

Examples

>>> gas = LennardJonesGas(n_particles=64, box_size=15.0, seed=0)
>>> history = gas.run(n_steps=200, steps_per_record=50)
>>> history["H"].shape[0]
5

speeds()

Instantaneous particle speeds \(|v_i|\).

step(n_steps=1)

Advance the simulation by n_steps Velocity Verlet integration steps.

Parameters:

n_steps (int, default=1) – Number of integration steps to perform.

temperature()

Instantaneous kinetic temperature from equipartition, \(T = \langle m v^2 \rangle / (d\, k_B)\).

With \(d=2\) spatial dimensions and 2 degrees of freedom removed for the fixed center-of-mass momentum.

total_energy()

Total (kinetic + potential) energy, conserved up to integration error.

class physicskit.statphys.OrnsteinUhlenbeck(theta=1.0, mu=0.0, sigma=1.0, n_paths=1000, seed=None)

Bases: object

An ensemble of Ornstein-Uhlenbeck processes, \(dX = \theta(\mu - X)\,dt + \sigma\,dW\).

The unique stationary, Gaussian, Markov process (Doob 1942). Started at \(x_0\) it has

\[\langle X_t\rangle = \mu + (x_0 - \mu)e^{-\theta t}, \qquad \mathrm{Var}\,X_t = \frac{\sigma^2}{2\theta}(1 - e^{-2\theta t}),\]

and in the stationary state the autocovariance is \(\frac{\sigma^2}{2\theta}e^{-\theta|\tau|}\) (G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev. 36, 823 (1930)). It is the velocity of a free Langevin particle (\(\theta = \gamma/m\), \(\sigma^2/2\theta = k_BT/m\)) and the position of an overdamped particle in a harmonic trap.

Parameters:

• theta (float, default=1.0) – Relaxation rate \(\theta > 0\).

• mu (float, default=0.0) – Long-time mean.

• sigma (float, default=1.0) – Noise amplitude.

• n_paths (int, default=1000) – Number of independent realizations.

• seed (int, optional) – Seed for the noise generator.

Examples

>>> ou = OrnsteinUhlenbeck(theta=2.0, sigma=1.0, n_paths=20000, seed=0)
>>> t, X = ou.run(t_max=5.0, dt=0.01, x0=3.0)
>>> round(ou.stationary_variance, 3), round(float(X[-1].var()), 2)
(0.25, 0.25)

autocovariance(tau)

Stationary autocovariance \(\frac{\sigma^2}{2\theta}e^{-\theta|\tau|}\).

Parameters:

tau (array_like)

Returns:

ndarray

mean(t, x0)

Exact ensemble mean at time t from a common start x0.

Parameters:

• t (array_like)

• x0 (float)

Returns:

ndarray

run(t_max, dt, n_frames=200, x0=None, method='exact')

Sample every path from x0 to t_max.

Parameters:

• t_max (float) – Total time.

• dt (float) – Time step.

• n_frames (int, default=200) – Number of recorded frames after the initial one.

• x0 (float or array_like of shape (n_paths,), optional) – Initial values; by default drawn from the stationary distribution.

• method ({"exact", "euler_maruyama"}, default="exact") – "exact" uses the exact Gaussian transition \(X_{t+\Delta t} = \mu + (X_t-\mu)e^{-\theta\Delta t} + \sigma\sqrt{(1-e^{-2\theta\Delta t})/2\theta}\;\mathcal N(0,1)\) (Gillespie, Phys. Rev. E 54, 2084 (1996)), free of time-step error; "euler_maruyama" integrates the SDE, with stationary variance \(\sigma^2/(\theta(2 - \theta\Delta t))\).

Returns:

• times (ndarray of shape (n_frames + 1,))

• paths (ndarray of shape (n_frames + 1, n_paths))

property stationary_variance: float

\(\sigma^2 / (2\theta)\).

variance(t)

Exact ensemble variance at time t from a deterministic start.

Parameters:

t (array_like)

Returns:

ndarray

class physicskit.statphys.PathIntegralParticle(beta, n_beads=32, omega=1.0, lam=0.0, seed=None)

Bases: object

Ring-polymer path-integral Monte Carlo for \(V(x) = \tfrac12\omega^2x^2 + \lambda x^4\).

Parameters:

• beta (float) – Inverse temperature (\(\hbar = m = k_B = 1\)).

• n_beads (int, default=32) – Number of imaginary-time slices \(P\).

• omega (float, default=1.0) – Harmonic frequency.

• lam (float, default=0.0) – Quartic anharmonicity \(\lambda \ge 0\).

• seed (int, optional) – Seed for reproducible sampling.

Variables:

path (ndarray of shape (n_beads,)) – Current bead positions.

Examples

>>> pimc = PathIntegralParticle(beta=2.0, n_beads=16, seed=0)
>>> out = pimc.run(n_equil=500, n_measure=4000)
>>> abs(out["x2"] - harmonic_x2_primitive(2.0, 16)) < 0.05
True

run(n_equil=1000, n_measure=10000)

Equilibrate, then sample the ring polymer.

Parameters:

• n_equil (int, default=1000) – Discarded sweeps.

• n_measure (int, default=10000) – Measured sweeps (one sample per sweep of \(P\) bead moves).

Returns:

dict – "x2" and "energy" (means of \(\langle x^2\rangle\) and of the virial energy estimator), their per-sweep series "x2_series" and "energy_series", and the bead-move "acceptance" rate.

class physicskit.statphys.Percolation2D(L=64, p=0.5, mode='site', seed=None)

Bases: object

Site or bond percolation on a periodic-free (open-boundary) square lattice.

Parameters:

• L (int, default=64) – Linear lattice size.

• p (float, default=0.5) – Occupation probability: fraction of sites (site percolation) or bonds (bond percolation) that are open.

• mode ({"site", "bond"}, default="site") – Percolation type. The square-lattice thresholds are \(p_c \approx 0.592746\) (site, known only numerically) and \(p_c = 0.5\) (bond, exact by self-duality).

• seed (int, optional) – Seed for reproducible lattice realizations.

Variables:

labels (ndarray of shape (L, L)) – Cluster labels of the current realization, from hoshen_kopelman().

P_C_BOND = 0.5

Exact threshold for square-lattice bond percolation (self-duality).

P_C_SITE = 0.592746

Known threshold for square-lattice site percolation (numerical estimate).

cluster_size_distribution()

Sizes of every cluster in the current realization.

At the percolation threshold, the number of clusters of size s follows a power law \(n_s \sim s^{-\tau}\) with the Fisher exponent \(\tau = 187/91 \approx 2.055\) for 2D percolation – the geometric analogue of the divergent susceptibility seen at a thermal critical point.

Returns:

ndarray – Cluster sizes, in sites, one entry per cluster (including the spanning cluster if one exists).

fractal_dimension(n_trials=20)

Estimate the fractal dimension of the largest cluster at \(p = p_c\).

A single lattice size gives no direct mass-versus-L scaling, so instead this estimates \(d_f\) from the mass-radius relation of the largest cluster in each realization, \(d_f = \log(\text{cluster size}) / \log(R_g)\), where \(R_g\) is the cluster’s radius of gyration, averaged over trials. The 2D percolation universality class predicts \(d_f = 91/48 \approx 1.896\).

Parameters:

n_trials (int, default=20) – Number of critical realizations to average over.

Returns:

float – Estimated fractal dimension.

generate(p=None)

Draw a fresh random realization and label its clusters.

Parameters:

p (float, optional) – Occupation probability to use; defaults to self.p. If given, also updates self.p.

Returns:

ndarray of shape (L, L) – The updated cluster labels array.

largest_cluster_size()

Number of (real) sites in the largest cluster of the current realization.

property p_c

Threshold occupation probability for the current mode.

spanning_labels()

Cluster labels of the current realization that span top to bottom.

For bond percolation the doubled lattice is used internally, but returned labels index the real (even-index) sites.

Returns:

ndarray – Labels present on both the first and last row of real sites.

spanning_probability(p_values, n_trials=100)

Monte Carlo estimate of the spanning probability \(P_{\text{span}}(p)\).

For each probability in p_values, generates n_trials independent realizations and reports the fraction that contain a top-to-bottom spanning cluster. As \(L \to \infty\) this curve sharpens into a step function at \(p_c\).

Parameters:

• p_values (array_like) – Occupation probabilities to sample.

• n_trials (int, default=100) – Independent realizations per probability.

Returns:

• p_values (ndarray)

• P_span (ndarray) – Estimated spanning probability at each p.

Examples

>>> perc = Percolation2D(L=20, seed=0)
>>> p, P = perc.spanning_probability([0.3, 0.593, 0.8], n_trials=20)
>>> bool(P[0] < P[2])
True

spans()

Whether the current realization has a top-to-bottom spanning cluster.

class physicskit.statphys.PottsModel2D(L=32, q=3, J=1.0, kB=1.0, seed=None)

Bases: object

The 2D ferromagnetic q-state Potts model on a periodic \(L \times L\) lattice.

Each site carries a discrete state \(s_i \in \{0, \ldots, q-1\}\), and the Hamiltonian rewards neighboring sites for matching: \(H = -J \sum_{\langle i,j \rangle} \delta(s_i, s_j)\). The Ising model is the special case \(q = 2\) (with a rescaled coupling). On the square lattice the transition is second order for \(q \le 4\) and first order for \(q > 4\).

Parameters:

• L (int, default=32) – Linear lattice size.

• q (int, default=3) – Number of Potts states.

• J (float, default=1.0) – Coupling constant.

• kB (float, default=1.0) – Boltzmann constant.

• seed (int, optional) – Seed for reproducible dynamics.

Variables:

spins (ndarray of shape (L, L)) – Current state configuration, values in {0, ..., q-1}.

property T_C

Exact critical temperature of the square-lattice Potts model, \(J / (k_B \ln(1 + \sqrt{q}))\).

energy()

Total energy \(-J \sum_{\langle i,j\rangle} \delta(s_i, s_j)\) of the current configuration.

order_parameter()

Potts order parameter, the rescaled fraction of sites in the majority state.

\[m = \frac{q \cdot n_{\max} / N - 1}{q - 1}\]

where \(n_{\max}\) is the largest population among the q states. m is 0 in the fully disordered phase and 1 in the fully ordered phase.

Returns:

float

run_temperature_sweep(temperatures, n_equil=200, n_measure=500, measure_every=1)

Sweep over temperatures, measuring energy, order parameter, and specific heat.

See Ising2D.run_temperature_sweep() for the equilibration and measurement protocol; this follows the same pattern.

Parameters:

• temperatures (array_like) – Temperatures to sample.

• n_equil (int, default=200) – Equilibration sweeps at each temperature.

• n_measure (int, default=500) – Measurements collected at each temperature.

• measure_every (int, default=1) – Sweeps between successive measurements.

Returns:

dict of str -> ndarray – Keys "T", "E", "m", "C_v".

sweep(beta, n_sweeps=1)

Advance the lattice by n_sweeps Metropolis-Hastings sweeps.

Parameters:

• beta (float) – Inverse temperature.

• n_sweeps (int, default=1) – Number of sweeps to perform.

class physicskit.statphys.RandomWalk(n_walkers=500, n_steps=1000, dim=2, kind='lattice', step_std=1.0, seed=None)

Bases: object

An ensemble of independent random walkers used to study diffusion.

Parameters:

• n_walkers (int, default=500) – Number of independent walkers in the ensemble.

• n_steps (int, default=1000) – Number of steps each walker takes.

• dim (int, default=2) – Spatial dimension.

• kind ({"lattice", "gaussian"}, default="lattice") – "lattice" walkers take unit steps along a randomly chosen axis (a simple-cubic lattice walk); "gaussian" walkers take continuous, independent Gaussian steps in every dimension at every time step (a discrete-time model of Brownian motion).

• step_std (float, default=1.0) – Per-axis standard deviation of a single step. Only used for kind="gaussian".

• seed (int, optional) – Seed for reproducible trajectories.

Variables:

trajectories (ndarray of shape (n_steps + 1, n_walkers, dim) or None) – Walker positions at every time, populated by run().

Examples

>>> walk = RandomWalk(n_walkers=200, n_steps=200, dim=2, seed=0)
>>> t, msd = walk.mean_squared_displacement()
>>> bool(msd[-1] > msd[10])
True

diffusion_coefficient()

Estimate the diffusion coefficient D via the Einstein relation.

Fits \(\text{MSD}(t) = 2 d D t\) (a line through the origin, least-squares) to the ensemble MSD, where \(d\) is dim.

Returns:

float – Estimated diffusion coefficient.

final_displacement_histogram(axis=0, bins=40)

Density histogram of the final-time displacement along one axis.

A demonstration of the central limit theorem: regardless of the (highly non-Gaussian) single-step distribution, the sum of many independent steps approaches a Gaussian.

Parameters:

• axis (int, default=0) – Which spatial axis to histogram.

• bins (int, default=40) – Number of histogram bins.

Returns:

• centers (ndarray of shape (bins,)) – Bin center positions.

• density (ndarray of shape (bins,)) – Probability density in each bin.

mean_squared_displacement()

Ensemble-averaged mean squared displacement \(\text{MSD}(t) = \langle |r(t)|^2 \rangle\).

Returns:

• t (ndarray of shape (n_steps + 1,)) – Time (in steps), starting at 0.

• msd (ndarray of shape (n_steps + 1,)) – Mean squared displacement from the (shared) origin at each time.

run()

Simulate every walker’s full trajectory.

Returns:

ndarray of shape (n_steps + 1, n_walkers, dim) – The same array stored in trajectories.

class physicskit.statphys.SherringtonKirkpatrick(N=128, J=1.0, kB=1.0, seed=None)

Bases: object

The infinite-range Sherrington-Kirkpatrick Ising spin glass on N fully-connected spins.

Every pair of spins is coupled by an independent Gaussian random bond,

\[H = -\sum_{i<j} J_{ij}\, s_i s_j, \qquad J_{ij} \sim \mathcal{N}\!\left(0, \frac{J^2}{N}\right),\]

where the \(1/N\) variance scaling is what keeps the energy extensive as \(N \to \infty\). Rather than attempt Parisi’s analytic replica-symmetry-breaking solution, overlap_distribution() measures its numerically observable fingerprint directly: the distribution \(P(q)\) of the replica overlap \(q = \frac{1}{N}\sum_i s_i^{(1)} s_i^{(2)}\), pooled over many independent disorder realizations. A single sharp peak at \(q=0\) signals the trivial (replica-symmetric) paramagnetic phase; a broad, non-Gaussian, non-self-averaging \(P(q)\) is the finite-size signature of the many-valley spin-glass phase that Parisi’s solution describes exactly in the \(N \to \infty\) limit.

Parameters:

• N (int, default=128) – Number of spins. Note one sweep costs \(O(N^2)\), since every spin couples to every other one.

• J (float, default=1.0) – Coupling scale; bond strengths are drawn from \(\mathcal{N}(0, J^2/N)\).

• kB (float, default=1.0) – Boltzmann constant.

• seed (int, optional) – Seed for the quenched disorder, initial spins, and dynamics.

Variables:

• spins (ndarray of shape (N,)) – Current spin configuration, values in {-1, +1}.

• J (ndarray of shape (N, N)) – Quenched, symmetric coupling matrix with zero diagonal. Fixed for the model’s lifetime.

energy()

Total energy \(-\sum_{i<j} J_{ij} s_i s_j\) of the current configuration.

overlap_distribution(beta, n_disorder=30, n_equil=300, n_measure=100, measure_every=1)

Sample the replica-overlap distribution P(q), pooled over independent disorder realizations.

For each of n_disorder independent bond realizations, two replicas are initialized randomly, thermalized under that shared disorder, and their overlap \(q = \frac{1}{N}\sum_i s_i^{(1)} s_i^{(2)}\) recorded repeatedly. Pooling samples across many disorder realizations (rather than just many measurements of one realization) is essential here: unlike the Edwards-Anderson model’s single well-defined \(\langle q^2 \rangle\), the SK spin glass’s overlap distribution is sample-dependent (it does not self-average) below the transition, and it is exactly this sample-to-sample variability that Parisi’s solution organizes into a hierarchy of pure states.

Parameters:

• beta (float) – Inverse temperature.

• n_disorder (int, default=30) – Number of independent quenched-disorder realizations to sample.

• n_equil (int, default=300) – Equilibration sweeps per replica, for each disorder realization.

• n_measure (int, default=100) – Number of overlap measurements per disorder realization.

• measure_every (int, default=1) – Sweeps between successive measurements.

Returns:

ndarray of shape (n_disorder * n_measure,) – Pooled overlap samples \(q \in [-1, 1]\).

sweep(beta, n_sweeps=1)

Advance the spins by n_sweeps Metropolis-Hastings sweeps (bonds are fixed).

Parameters:

• beta (float) – Inverse temperature.

• n_sweeps (int, default=1) – Number of sweeps to perform.

class physicskit.statphys.WangLandauIsing(L=8, J=1.0, seed=None)

Bases: object

Wang-Landau estimate of the density of states of the 2D Ising model.

The walk proposes single-spin flips accepted with probability \(\min(1, g(E_{\text{old}})/g(E_{\text{new}}))\) and raises \(\ln g\) of the current level by \(\ln f\) after every proposal. Whenever the visit histogram is flat (every visited level within flatness of the mean), it is reset and \(\ln f\) is halved, until \(\ln f\) falls below log_f_final.

Parameters:

• L (int, default=8) – Linear lattice size.

• J (float, default=1.0) – Ferromagnetic coupling.

• seed (int, optional) – Seed for the initial configuration and the walk.

Examples

>>> wl = WangLandauIsing(L=4, seed=0)
>>> E, log_g = wl.run(log_f_final=1e-3)
>>> round(float(np.exp(log_g[0])), 6)  # normalized to the two ground states
2.0

run(log_f_initial=1.0, log_f_final=1e-06, flatness=0.8, sweeps_per_check=1000)

Run the walk to convergence and return the normalized \(\ln g(E)\).

Parameters:

• log_f_initial (float, default=1.0) – Starting modification factor \(\ln f_0\).

• log_f_final (float, default=1e-6) – Stop once \(\ln f\) falls below this.

• flatness (float, default=0.8) – A histogram is flat when its minimum over the accessible levels is at least this fraction of its mean.

• sweeps_per_check (int, default=1000) – Sweeps between flatness checks.

Returns:

• energies (ndarray) – Energy levels \(E_k\).

• log_g (ndarray) – \(\ln g(E)\), normalized so that \(g(E_0) = 2\) (the two ground states), with -inf on the inaccessible levels.

class physicskit.statphys.XYModel2D(L=32, J=1.0, kB=1.0, seed=None)

Bases: object

The 2D classical XY model on a periodic \(L \times L\) lattice.

Each site carries a planar spin angle \(\theta_i \in [0, 2\pi)\), with Hamiltonian \(H = -J \sum_{\langle i,j \rangle} \cos(\theta_i - \theta_j)\). The Mermin-Wagner theorem forbids conventional long-range order in 2D, but the model still hosts the topological Kosterlitz-Thouless (KT) transition at \(T_{\text{KT}} \approx 0.893\, J / k_B\): bound vortex-antivortex pairs at low temperature unbind into a free vortex plasma above \(T_{\text{KT}}\).

Parameters:

• L (int, default=32) – Linear lattice size.

• J (float, default=1.0) – Coupling constant.

• kB (float, default=1.0) – Boltzmann constant.

• seed (int, optional) – Seed for reproducible dynamics.

Variables:

theta (ndarray of shape (L, L)) – Current spin angles, in radians.

property T_KT

Approximate (numerically established) Kosterlitz-Thouless transition temperature, \(\approx 0.893\, J / k_B\).

energy()

Total energy \(-J \sum_{\langle i,j\rangle} \cos(\theta_i - \theta_j)\) of the current configuration.

magnetization_vector()

Net planar magnetization vector \((\sum_i \cos\theta_i, \sum_i \sin\theta_i)\).

sweep(beta, n_sweeps=1, delta=1.0)

Advance the lattice by n_sweeps Metropolis-Hastings sweeps.

Parameters:

• beta (float) – Inverse temperature.

• n_sweeps (int, default=1) – Number of sweeps to perform.

• delta (float, default=1.0) – Maximum trial angular step (radians); tune towards an acceptance rate near 50% for efficient sampling.

vortex_count(threshold=0.5)

Count vortices and antivortices by thresholding the plaquette vorticity.

Parameters:

threshold (float, default=0.5) – Minimum |vorticity| for a plaquette to be counted as hosting a topological charge.

Returns:

• n_vortices (int) – Number of plaquettes with vorticity near +1.

• n_antivortices (int) – Number of plaquettes with vorticity near -1.

vorticity()

Topological charge of every plaquette; see physicskit.statphys.core.monte_carlo.xy_plaquette_vorticity().

Returns:

ndarray of shape (L, L) – Plaquette vorticity, close to an integer everywhere and nonzero (\(\pm 1\)) at vortex/antivortex cores.

physicskit.statphys.canonical_from_density_of_states(energies, log_g, temperatures, n_sites, kB=1.0)

Canonical energy and specific heat per site from a density of states.

\[\langle E^n \rangle_\beta = \frac{\sum_E E^n g(E) e^{-\beta E}}{\sum_E g(E) e^{-\beta E}}, \qquad C_v = \frac{\langle E^2 \rangle - \langle E \rangle^2}{k_B T^2 N}.\]

Parameters:

• energies (array_like) – Energy levels.

• log_g (array_like) – \(\ln g(E)\) on those levels; -inf marks an empty level.

• temperatures (array_like) – Temperatures at which to evaluate.

• n_sites (int) – Number of sites, to report per-site quantities.

• kB (float, default=1.0) – Boltzmann constant.

Returns:

dict of str -> ndarray – Keys "T", "E" (per site), "C_v" (per site) and "log_Z".

Examples

>>> E, g = ising_density_of_states_exact(2)
>>> out = canonical_from_density_of_states(E, np.log(g), [1e-3], n_sites=4)
>>> float(out["E"][0])
-2.0

physicskit.statphys.einstein_diffusion_coefficient(kT, gamma)

Einstein relation \(D = \mu k_B T = k_B T / \gamma\).

(A. Einstein, Ann. Phys. 322, 549 (1905), eq. for \(D\) in §4.)

Parameters:

• kT (float) – Temperature in energy units (\(k_B = 1\)).

• gamma (float) – Friction coefficient, the inverse of the mobility \(\mu\).

Return type:

float

Returns:

float – Diffusion coefficient.

Examples

>>> einstein_diffusion_coefficient(kT=2.0, gamma=4.0)
0.5

physicskit.statphys.fokker_planck_1d(p0, x, drift, diffusion, t_max, dt, n_frames=100, theta=0.5)

Evolve a 1D probability density under the Fokker-Planck equation.

Uses the finite-volume operator of fokker_planck_operator() with the \(\theta\)-scheme \((I - \theta\Delta t L)p^{n+1} = (I + (1-\theta)\Delta t L)p^n\): Crank-Nicolson (second order, \(\theta = 1/2\)) by default, or backward Euler (\(\theta = 1\), strongly damping) for rough initial data. Both are unconditionally stable; the matrix is factored once.

Parameters:

• p0 (array_like, shape (n,)) – Initial density on the grid (normalized to unit mass on return).

• x (array_like, shape (n,)) – Uniform grid of cell centers; the walls sit half a cell outside the end points.

• drift (float, array_like, or callable) – \(A(x)\) and \(D(x)\), as in fokker_planck_operator().

• diffusion (float, array_like, or callable) – \(A(x)\) and \(D(x)\), as in fokker_planck_operator().

• t_max (float) – Final time (rounded to a whole number of steps per frame).

• dt (float) – Time step.

• n_frames (int, default=100) – Number of recorded frames after the initial one.

• theta (float, default=0.5) – Implicitness, in [0.5, 1].

Returns:

• times (ndarray, shape (n_frames + 1,))

• p (ndarray, shape (n_frames + 1, n)) – Density at each frame, each with \(\sum p\,\Delta x = 1\).

Examples

Free diffusion spreads a Gaussian as \(\sigma^2(t) = \sigma_0^2 + 2Dt\):

>>> x = np.linspace(-10, 10, 401)
>>> p0 = np.exp(-x**2 / 2)
>>> t, p = fokker_planck_1d(p0, x, 0.0, 0.5, t_max=2.0, dt=0.01, n_frames=4)
>>> var = np.sum(x**2 * p[-1]) * (x[1] - x[0])
>>> round(float(var), 3)  # 1 + 2 * 0.5 * 2
3.0

physicskit.statphys.fokker_planck_operator(x, drift, diffusion)

Finite-volume Fokker-Planck operator \(L\) with \(dp/dt = Lp\) and zero-flux walls.

Each grid point is the center of a cell of width \(\Delta x\). The flux through the face between cells \(i\) and \(i+1\) is

\[J_{i+1/2} = A_{i+1/2}\,\frac{p_i + p_{i+1}}{2} - \frac{D_{i+1}p_{i+1} - D_i p_i}{\Delta x},\]

with \(J = 0\) through the two outer walls, so \(\sum_i p_i\,\Delta x\) is conserved exactly. Central differencing of the drift needs a cell Péclet number \(|A|\Delta x / D\) below about 2 to keep \(p \ge 0\).

Parameters:

• x (array_like, shape (n,)) – Uniform grid of cell centers.

• drift (float, array_like of shape (n,), or callable) – \(A(x)\). A callable is evaluated at the cell faces; an array is averaged onto them.

• diffusion (float, array_like of shape (n,), or callable) – \(D(x) > 0\), evaluated at the cell centers.

Return type:

csr_matrix

Returns:

scipy.sparse.csr_matrix, shape (n, n)

Examples

Columns sum to zero: probability is conserved.

>>> L = fokker_planck_operator(np.linspace(-1, 1, 11), lambda x: -x, 0.5)
>>> bool(np.allclose(np.asarray(L.sum(axis=0)).ravel(), 0.0))
True

physicskit.statphys.fokker_planck_stationary(x, drift, diffusion)

Zero-flux stationary solution of the 1D Fokker-Planck equation.

Setting \(J = 0\) gives

\[p_{\text{st}}(x) = \frac{N}{D(x)}\exp\!\int^x \frac{A(y)}{D(y)}\,dy\]

(Risken 1989, eq. 5.6). For \(A = -U'/\gamma\) and \(D = k_BT/\gamma\) this is the Boltzmann distribution \(e^{-U/k_BT}/Z\).

Parameters:

• x (array_like, shape (n,)) – Grid.

• drift (float, array_like of shape (n,), or callable) – \(A(x)\) and \(D(x)\), evaluated on the grid.

• diffusion (float, array_like of shape (n,), or callable) – \(A(x)\) and \(D(x)\), evaluated on the grid.

Return type:

ndarray

Returns:

ndarray, shape (n,) – Stationary density normalized by the trapezoid rule on x.

Examples

An Ornstein-Uhlenbeck process, \(A = -\theta x\), \(D = \sigma^2/2\), is Gaussian with variance \(\sigma^2/2\theta\):

>>> x = np.linspace(-6, 6, 2001)
>>> p = fokker_planck_stationary(x, lambda y: -2.0 * y, 0.5)
>>> round(float(np.sum(x**2 * p) * (x[1] - x[0])), 4)
0.25

physicskit.statphys.harmonic_x2_exact(beta, omega=1.0)

Exact thermal \(\langle x^2 \rangle\) of a quantum harmonic oscillator.

\[\langle x^2 \rangle = \frac{1}{2\omega} \coth\frac{\beta\omega}{2} \quad (\hbar = m = 1).\]

Parameters:

• beta (float) – Inverse temperature.

• omega (float, default=1.0) – Oscillator frequency.

Returns:

float

Examples

>>> round(float(harmonic_x2_exact(100.0)), 6)  # ground state: 1/(2 omega)
0.5

physicskit.statphys.harmonic_x2_primitive(beta, n_beads, omega=1.0)

Exact \(\langle x^2 \rangle\) of the \(P\)-bead primitive-action ring polymer.

For a harmonic well the discretized action is a Gaussian \(\tfrac12 x^T A x\) with \(A = \tfrac{P}{\beta}(2I - S - S^T) + \tfrac{\beta\omega^2}{P} I\) (\(S\) the cyclic shift), so \(\langle x^2 \rangle_P = \operatorname{tr}(A^{-1})/P\). This is the exact target of a PathIntegralParticle run with lam=0, and tends to harmonic_x2_exact() as \(P \to \infty\); the gap is the \(O(\tau^2)\) Trotter error.

Parameters:

• beta (float) – Inverse temperature.

• n_beads (int) – Number of imaginary-time slices \(P\).

• omega (float, default=1.0) – Oscillator frequency.

Returns:

float

Examples

>>> exact = harmonic_x2_exact(5.0)
>>> bool(abs(harmonic_x2_primitive(5.0, 200) - exact) < 1e-3 * exact)
True

physicskit.statphys.ising_density_of_states_exact(L, J=1.0)

Exact density of states of a periodic \(L \times L\) Ising lattice, by enumeration.

All \(2^{L^2}\) configurations are enumerated, so this is only practical up to \(L = 4\) (65536 states).

Parameters:

• L (int) – Linear lattice size, at most 4.

• J (float, default=1.0) – Ferromagnetic coupling.

Returns:

• energies (ndarray of shape (L*L + 1,)) – Energy levels \(E_k = -2JN + 4Jk\).

• g (ndarray of shape (L*L + 1,)) – Number of configurations at each level (zero for the inaccessible levels \(E = \pm(2JN - 4J)\)).

Examples

>>> E, g = ising_density_of_states_exact(2)
>>> int(g.sum()), int(g[0])
(16, 2)

physicskit.statphys.stokes_einstein_diffusion_coefficient(kT, viscosity, radius)

Stokes-Einstein diffusion coefficient of a sphere, \(D = k_B T / (6\pi\eta a)\).

Combines the Einstein relation with Stokes’ drag \(\gamma = 6\pi\eta a\) on a sphere of radius \(a\) in a fluid of viscosity \(\eta\).

Parameters:

• kT (float) – Thermal energy \(k_B T\) (in joules for SI inputs).

• viscosity (float) – Dynamic viscosity \(\eta\).

• radius (float) – Sphere radius \(a\).

Return type:

float

Returns:

float – Diffusion coefficient.

Examples

A 0.5 micron bead in room-temperature water diffuses about 0.43 square microns per second:

>>> D = stokes_einstein_diffusion_coefficient(1.380649e-23 * 293.0, 1.0e-3, 0.5e-6)
>>> round(D * 1e12, 2)
0.43