PeaPods

_{A 2D XY model quenched to T = 0.1: vortex–antivortex pairs, where every color meets, annihilate until one domain remains (script).}

A Python library for simulating Ising and XY spin systems with modern Monte Carlo methods. The core simulation loop is written in Rust (via PyO3) for performance, with a thin Python wrapper for ease of use. Currently focused on spin glass simulation, where the combination of parallel tempering, replica cluster moves, and overlap statistics is of particular interest.

Reproduces published results

• 3D ±J spin glass. The spin-glass correlation lengths ξ/L of L = 4 to 12 (3840 disorder samples each) cross at T = 1.204(17), 1.178(19) and 1.146(18) for (L, 2L) = (4, 8), (5, 10), (6, 12). With the Janus correction exponent the crossings extrapolate to T_c = 1.085(42), against Janus's 1.1019(29) (Baity-Jesi et al. 2013); ξ/L at the crossings rises from 0.56 toward their 0.652 as L grows (validation/spin_glass_3d/, about 5 hours on 96 CPU cores).
• 2D XY at the BKT point. At β = 1.1199 the helicity modulus, ξ/L and χ for L = 16 and 32 agree with Hasenbusch (2005), Table 1, within 1.9 combined standard errors (validation/xy_finite_size.py).
• Ising ferromagnets. Binder cumulants of square and triangular lattices (exact T_c) and simple-cubic, BCC and FCC lattices cross at the known T_c on every CI run (tests/binder_crossings.py).

Scripts and references for each result are listed in validation/.

Features

• Ising ferromagnets and spin glasses on periodic Bravais lattices (hypercubic, triangular, or any custom neighbor offsets)
• XY vectors with signed couplings, site/bond dilution, embedded SW/Wolff and optional overrelaxation
• Physical XY thermodynamics, helicity, structure factors, uniform correlation lengths, optional correlations and angle winding
• Arbitrary, bimodal (±J), or Gaussian coupling distributions
• Multiple replicas with overlap statistics for spin glass order parameters

The following algorithms are currently supported:

• Single-spin flips (Metropolis and Gibbs sampling)
• Swendsen-Wang cluster updates
• Wolff cluster updates
• Parallel tempering
• Houdayer isoenergetic cluster move (replica cluster move for spin glasses)
• Jörg move (stochastic overlap cluster move)
• CMR move (Chayes-Machta-Redner blue-bond overlap cluster move)
• Replica Monte Carlo (Swendsen-Wang τ-cluster move between adjacent temperatures)
• Passive full-graph FK and overlap-cluster observations with per-disorder statistics
• Full-ladder parallel tempering diagnostics and reproducible optional seeds

Quickstart

Run simulations directly from the terminal:

# 2D ferromagnet with cluster updates and parallel tempering
peapods simulate --shape 32 32 --temp-min 1.5 --temp-max 3.0 \
    --n-sweeps 5000 --cluster-interval 1 --pt-interval 1 --n-replicas 2

# 3D spin glass with Houdayer ICM
peapods simulate --shape 8 8 8 --couplings bimodal \
    --temp-min 0.8 --temp-max 1.4 --n-temps 24 --n-sweeps 10000 \
    --pt-interval 1 --overlap-cluster-update-interval 1 --n-replicas 4

# Save full results to .npz
peapods simulate --shape 16 16 --temp-min 1.5 --temp-max 3.0 \
    --n-sweeps 5000 --pt-interval 1 --n-replicas 2 -o results.npz

For full control over geometry and parameters, use the Python API directly:

import numpy as np
from peapods import Ising

# Triangular lattice ferromagnet (T_c = 4/ln(3) ≈ 3.641)
tri = Ising((32, 32), temperatures=np.linspace(3.0, 4.2, 32),
            n_replicas=2, neighbor_offsets=[[1, 0], [0, 1], [1, -1]])
tri.sample(n_sweeps=5000, sweep_mode="metropolis",
           cluster_update_interval=1, pt_interval=1)
print(tri.binder_cumulant)

Pass an explicit seed for deterministic couplings and dynamics. Full SW-style graphs can be measured without acting on spins:

model = Ising((32, 32), couplings="bimodal", temperatures=np.linspace(1, 3, 16),
              n_replicas=2, seed=42)
result = model.sample(n_sweeps=5000, cluster_update_interval=10,
                      cluster_mode="sw", cluster_action="observe",
                      pt_interval=1, pt_schedule="full_ladder")
fk = result["per_disorder"]["cluster_observations"]["fk"]

Autocorrelation diagnostics use the exact bounded-memory ring backend by default. Set autocorrelation_backend="fft" together with autocorrelation_max_lag for faster evaluation when retaining the full measurement history and its higher memory use are acceptable.

More complete examples:

• validation/spin_glass_3d/: spin-glass correlation length and T_c crossings, end to end
• validation/plot/overlap_csd.py: CMR blue-cluster size distribution of the 2D ±J spin glass
• skills/simulate/sweep_config.toml: a peapods sweep --config template
• The agent skills index every model, lattice, move and output

XY models

The Rust and Python APIs support periodic hypercubic XY systems with finite signed couplings. See the XY guide for conventions, dilution, block errors, examples, and bounded comparison with published finite-size measurements. CLI and sweep-framework XY integration are deferred.

from peapods import XY

xy = XY((8, 8), temperatures=[0.6, 1.0, 2.0], seed=42)
result = xy.sample(1024, collect_blocks=True)
print(result["helicity_modulus"])

Installation

We recommend uv for package management. Bare pip works too (if you don't mind using a worse management tool).

uv venv
source .venv/bin/activate
uv pip install peapods

Pre-built wheels are available for Linux (x86_64, aarch64), macOS (Intel, Apple Silicon), and Windows (x86_64).

Agent skills

skills/ holds Agent Skills that index peapods for coding agents: simulate (models, lattices, moves, outputs), analyze (equilibration, errors, critical points), benchmark and develop. In Claude Code, install them as a plugin:

/plugin marketplace add PeaBrane/peapods
/plugin install peapods@peapods

Other agents can load the SKILL.md files directly.

Building from source

Requires a Rust toolchain and maturin:

uv pip install maturin numpy
maturin develop --release