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 and Wolff-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.

• 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.

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.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"}, 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"}, default="metropolis") – Update rule. Metropolis flips single spins; Wolff grows and flips whole clusters, which decorrelates much faster near \(T_c\).

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

Returns:

int or None – For "wolff", the size of the last cluster flipped; 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.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.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.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.