Statistical mechanics

Statistical mechanics: gases, fluids, spin models, lattice processes, urns, master equations, and their measures.

BoltzmannAverageElement(label, log_partition=0.0, energy=0.0, entropy=0.0)

Bases: Element

Log partition function, mean energy, and entropy per spin of a group of levels.

Parameters:

Name	Type	Description	Default
label	str	The label of the element.	required
log_partition	Annotated[float, StateSpec()]	The initial guess of the log partition function, which the equations replace.	0.0
energy	Annotated[float, StateSpec()]	The initial guess of the mean energy, which the equations replace.	0.0
entropy	Annotated[float, StateSpec()]	The initial guess of the entropy per spin, which the equations replace.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import BoltzmannAverageElement
>>> BoltzmannAverageElement('ensemble').entropy
0.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_boltzmann_average.py

def __init__(
    self,
    label: str,
    log_partition: Annotated[float, StateSpec()] = 0.0,
    energy: Annotated[float, StateSpec()] = 0.0,
    entropy: Annotated[float, StateSpec()] = 0.0,
) -> None:
    super().__init__(label)
    self.log_partition = log_partition
    self.energy = energy
    self.entropy = entropy

BoltzmannSumInteraction(label, *, target, levels, temperature=1.0, n_spins=1.0)

Bases: Interaction

Boltzmann sums over a group of levels, which give the log partition function, mean energy, and entropy per spin.

The log partition function is ln Z = log sum_j exp(-E_j / T), the mean energy is U = sum_j E_j exp(-E_j / T - ln Z), and the entropy per spin is (ln Z + U / T) / N. The sums reduce in log space, so they do not overflow at any temperature.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the element that holds the averages.	required
levels	Annotated[str, RoleSpec()]	The path of the group of levels.	required
temperature	Annotated[float, ParameterSpec(domain=open(0, oo))]	The temperature.	1.0
n_spins	Annotated[float, ParameterSpec(domain=open(0, oo))]	The number of spins the entropy is divided by.	1.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import BoltzmannSumInteraction
>>> BoltzmannSumInteraction('sums', target='ensemble', levels='levels', temperature=0.5).temperature
0.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_boltzmann_sum.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    levels: Annotated[str, RoleSpec()],
    temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    n_spins: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label)
    self.target = target
    self.levels = levels
    self.temperature = temperature
    self.n_spins = n_spins

BondSamplingInteraction(label, *, target)

Bases: SpinSamplingInteraction

Metropolis sweeps of Ising spins with quenched couplings to their neighbours, of energy -s (sum J_k s_k + h).

The couplings J_up, J_left, J_right, and J_down are constant states of every spin, in the order of the offsets of the von Neumann neighbourhood, and the field h and the temperature are its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction, unique among the parts of its system.	required
target	Annotated[str, RoleSpec()]	The path of the lattice group of spins.	required

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import BondSamplingInteraction
>>> BondSamplingInteraction('sampling', target='spins').neighborhood
VonNeumannNeighborhood()

Source code in src/skcomplex/models/physics/statistical_mechanics/_bond_sampling.py

def __init__(self, label: str, *, target: Annotated[str, RoleSpec()]) -> None:
    super().__init__(label, target=target, state='spin', energy=self._energy, neighborhood=NEAREST)

BondSpinElement(label, visual=None, *, spin=1.0, J_up=1.0, J_left=1.0, J_right=1.0, J_down=1.0, h=0.0, temperature=1.0)

Bases: Element

Ising spin on a square lattice with quenched couplings to its four neighbours, of local energy -s sum J_k s_k.

The couplings are constant states of every site, so a bond between two sites is stored at both of them with the same value. A BondSamplingInteraction flips the spins by Metropolis Monte Carlo at the temperature T, with the energy -sum_<ij> J_ij s_i s_j - h sum_i s_i.

Parameters:

Name	Type	Description	Default
label	str	The label of the spin or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
spin	Annotated[float, StateSpec(domain=ZZ)]	The initial spin, 1 or -1.	1.0
J_up	Annotated[float, StateSpec()]	The coupling to the upper neighbour.	1.0
J_left	Annotated[float, StateSpec()]	The coupling to the left neighbour.	1.0
J_right	Annotated[float, StateSpec()]	The coupling to the right neighbour.	1.0
J_down	Annotated[float, StateSpec()]	The coupling to the lower neighbour.	1.0
h	Annotated[float, ParameterSpec()]	The external field.	0.0
temperature	Annotated[float, ParameterSpec(domain=open(0, oo))]	The temperature.	1.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import BondSpinElement
>>> BondSpinElement('spin', J_left=-1.0).J_left
-1.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_bond_spin.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    spin: Annotated[float, StateSpec(domain=ZZ)] = 1.0,
    J_up: Annotated[float, StateSpec()] = 1.0,
    J_left: Annotated[float, StateSpec()] = 1.0,
    J_right: Annotated[float, StateSpec()] = 1.0,
    J_down: Annotated[float, StateSpec()] = 1.0,
    h: Annotated[float, ParameterSpec()] = 0.0,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.spin = spin
    self.J_up = J_up
    self.J_left = J_left
    self.J_right = J_right
    self.J_down = J_down
    self.h = h
    self.temperature = temperature

ContactSiteElement(label, visual=None, *, infected=0.0, rate=3.0, tau=0.1)

Bases: Element

Site of the contact process, which infects a neighbour at the rate rate / 4 and recovers at rate 1.

One step is a time slice tau, in which an infected site recovers with probability 1 - exp(-tau) and a healthy site is infected with probability 1 - exp(-rate tau f), where f is the fraction of its infected neighbours.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
infected	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial state, 1 for infected and 0 for healthy.	0.0
rate	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The infection rate.	3.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time slice of one step.	0.1

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ContactSiteElement
>>> ContactSiteElement('site').rate
3.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_contact_site.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    infected: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
    rate: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
) -> None:
    super().__init__(label, visual=visual)
    self.infected = infected
    self.rate = rate
    self.tau = tau

DephasingInteraction(label, *, target, gamma=0.1)

Bases: Interaction

Dephasing of every spin of a chain, a dissipator with the collapse operator S^z_i at the rate gamma.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the spin chain.	required
gamma	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The dephasing rate.	0.1

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import DephasingInteraction
>>> DephasingInteraction('dephasing', target='chain', gamma=0.1).gamma
0.1

Source code in src/skcomplex/models/physics/statistical_mechanics/_dephasing.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
) -> None:
    super().__init__(label)
    self.target = target
    self.gamma = gamma

DiskElement(label, x=0.0)

Bases: Element

Disk of the plane at a position that no law moves.

Parameters:

Name	Type	Description	Default
label	str	The label of the disk or the group.	required
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The position of its centre.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import DiskElement
>>> DiskElement('disk', x=(1.0, 2.0)).x
(1.0, 2.0)

Source code in src/skcomplex/models/physics/statistical_mechanics/_disk.py

def __init__(
    self,
    label: str,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
) -> None:
    super().__init__(label)
    self.x = x

EhrenfestUrnElement(label, visual=None, *, n=0, N=40)

Bases: Element

First box of the Ehrenfest urn, which holds n of the N balls, one of which moves to the other box every step.

The count falls by one with the probability n / N that the drawn ball is in the box, and rises by one otherwise.

Parameters:

Name	Type	Description	Default
label	str	The label of the box.	required
visual	VisualSpec | None	How the box is drawn, or None for default_visual_.	None
n	Annotated[int, StateSpec(domain=NN)]	The initial number of balls in the box.	0
N	Annotated[int, ParameterSpec(domain=open(0, oo))]	The number of balls in both boxes.	40

Attributes:

Name	Type	Description
default_visual_		The visual of the element when it sets none.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import EhrenfestUrnElement
>>> EhrenfestUrnElement('urn', n=10).N
40

Source code in src/skcomplex/models/physics/statistical_mechanics/_urns.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    n: Annotated[int, StateSpec(domain=NN)] = 0,
    N: Annotated[int, ParameterSpec(domain=RealSet.open(0, oo))] = 40,
) -> None:
    super().__init__(label, visual=visual)
    self.n = n
    self.N = N

EnergyLevelElement(label, energy=0.0)

Bases: Element

Energy level of a configuration, held fixed.

Parameters:

Name	Type	Description	Default
label	str	The label of the level or the group.	required
energy	Annotated[float, StateSpec()]	The energy of the level.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import EnergyLevelElement
>>> EnergyLevelElement('level', energy=-1.5).energy
-1.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_energy_level.py

def __init__(self, label: str, energy: Annotated[float, StateSpec()] = 0.0) -> None:
    super().__init__(label)
    self.energy = energy

ExclusionHopsInteraction(label, *, sites, p=1.0, q=0.0)

Bases: TransferInteraction

Hops of a particle from every site of a line to an empty neighbour, to the right at rate p and left at q.

No particle hops past the ends of the line. It is a transfer of one particle along every edge of the line, to the right and back, from an occupied site to an empty one.

Parameters:

Name	Type	Description	Default
label	str	The label of the hops.	required
sites	Annotated[str, RoleSpec()]	The group of sites, in line order.	required
p	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of hops to the right.	1.0
q	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate of hops to the left.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ExclusionHopsInteraction
>>> ExclusionHopsInteraction('hops', sites='lane').p
1.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_exclusion_hops.py

def __init__(
    self,
    label: str,
    *,
    sites: Annotated[str, RoleSpec()],
    p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    q: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(
        label,
        source=sites,
        target=sites,
        state='n',
        amount=1.0,
        rate=self._get_right_rate,
        reverse_rate=self._get_left_rate,
        offset=1,
        boundary='fixed',
    )
    self.sites = sites
    self.p = p
    self.q = q

ExclusionSiteElement(label, visual=None, *, occupied=0.0, p=0.5, alpha=0.3, beta=0.7)

Bases: Element

Site of the totally asymmetric simple exclusion process, empty at 0 or holding one particle at 1.

It holds the probabilities of a hop, of an entry at the left end of an open line, and of an exit at its right end, which the transfers of the process read.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
occupied	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial occupation.	0.0
p	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability of a hop.	0.5
alpha	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability of an entry at the left end of an open line.	0.3
beta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability of an exit at the right end of an open line.	0.7

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ExclusionSiteElement
>>> ExclusionSiteElement('site').p
0.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_exclusion_site.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    occupied: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
    p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
    beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.7,
) -> None:
    super().__init__(label, visual=visual)
    self.occupied = occupied
    self.p = p
    self.alpha = alpha
    self.beta = beta

FullyConnectedSpinElement(label, visual=None, *, spin=1.0, J=1.0, h=0.0, temperature=0.5, couplings=None)

Bases: Element

Ising spin coupled to every other spin, flipped at the Glauber rate 1 / (1 + exp(2 s h / T)).

The local field is h = J (sum_{j != i} s_j) / N + h_ext of the Curie-Weiss energy -(J / N) sum_{i<j} s_i s_j, or h = J sum_j C_ij s_j + h_ext with a symmetric coupling matrix C of zero diagonal. The rate satisfies detailed balance with the Boltzmann distribution at the temperature T.

Parameters:

Name	Type	Description	Default
label	str	The label of the spin or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
spin	Annotated[float, StateSpec(domain=ZZ)]	The initial spin, 1 or -1.	1.0
J	Annotated[float, ParameterSpec()]	The scale of the couplings.	1.0
h	Annotated[float, ParameterSpec()]	The external field.	0.0
temperature	Annotated[float, ParameterSpec(domain=open(0, oo))]	The temperature.	0.5
couplings	ArrayLike | None	The coupling matrix C, one row and one column per spin, or None for the uniform couplings 1 / N.	None

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import FullyConnectedSpinElement
>>> FullyConnectedSpinElement('spin').temperature
0.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_fully_connected_spins.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    spin: Annotated[float, StateSpec(domain=ZZ)] = 1.0,
    J: Annotated[float, ParameterSpec()] = 1.0,
    h: Annotated[float, ParameterSpec()] = 0.0,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    couplings: npt.ArrayLike | None = None,
) -> None:
    super().__init__(label, visual=visual)
    self.spin = spin
    self.J = J
    self.h = h
    self.temperature = temperature
    self.couplings = couplings

IsingSpinElement(label, visual=None, *, spin=1.0, J=1.0, h=0.0, temperature=2.0)

Bases: Element

Ising spin, up at 1 or down at -1, with its coupling J, its field h, and its temperature.

An IsingSamplingInteraction flips the spins by Monte Carlo sweeps at the temperature of the spins.

Parameters:

Name	Type	Description	Default
label	str	The label of the spin or the lattice group.	required
visual	VisualSpec | None	How the spin is drawn, or None for default_visual_.	None
spin	Annotated[float, StateSpec(domain=ZZ)]	The initial spin, 1 or -1.	1.0
J	Annotated[float, ParameterSpec()]	The coupling between neighbouring spins.	1.0
h	Annotated[float, ParameterSpec()]	The external field.	0.0
temperature	Annotated[float, ParameterSpec(domain=open(0, oo))]	The temperature in units of the coupling.	2.0

Attributes:

Name	Type	Description
default_visual_		The visual of the spin when it sets none.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import IsingSpinElement
>>> IsingSpinElement('spin').temperature
2.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_ising_spin.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    spin: Annotated[float, StateSpec(domain=ZZ)] = 1.0,
    J: Annotated[float, ParameterSpec()] = 1.0,
    h: Annotated[float, ParameterSpec()] = 0.0,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
) -> None:
    super().__init__(label, visual=visual)
    self.spin = spin
    self.J = J
    self.h = h
    self.temperature = temperature

KineticTemperatureAnalysis(group=None, *, mass=1.0)

Bases: DistributionAnalysis

Kinetic temperature m <|v|²> / d of the agents of every frame in d dimensions, with k_B = 1.

It is the second moment of the speeds of the agents, the distribution analysis of kind 'moments', scaled by the mass over the dimension, as equipartition gives k_B T / 2 per degree of freedom. It reads the velocities of the agents or particles of a trajectory, or an array of velocities of one snapshot, one row per agent, or of one snapshot per frame along its first axis.

Parameters:

Name	Type	Description	Default
group	str | None	The path of the agent or particle group of a trajectory, or None for the first one.	None
mass	float	The mass of every agent.	1.0

Attributes:

Name	Type	Description
values_		The velocities, frames by agents by dimensions.
temperature_		The kinetic temperature of every frame.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import KineticTemperatureAnalysis
>>> KineticTemperatureAnalysis().check().prepare([[1.0, 0.0], [0.0, 1.0]]).run().temperature_
array([0.5])

Source code in src/skcomplex/models/physics/statistical_mechanics/_kinetic_temperature.py

def __init__(self, group: str | None = None, *, mass: float = 1.0) -> None:
    super().__init__("moments", group=group)
    self.mass = mass

MattisSpinElement(label, visual=None, *, spin=1.0, J_up=1.0, J_left=1.0, J_right=1.0, J_down=1.0, h=0.0, temperature=1.0, pattern=1.0)

Bases: BondSpinElement

Spin of the Mattis model, a bond spin that also holds its quenched pattern xi, with couplings J xi_i xi_j.

Parameters:

Name	Type	Description	Default
label	str	The label of the spin or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
spin	float	The initial spin, 1 or -1.	1.0
J_up	float	The coupling to the upper neighbour.	1.0
J_left	float	The coupling to the left neighbour.	1.0
J_right	float	The coupling to the right neighbour.	1.0
J_down	float	The coupling to the lower neighbour.	1.0
h	float	The external field.	0.0
temperature	float	The temperature.	1.0
pattern	Annotated[float, StateSpec(domain=ZZ)]	The quenched pattern xi, 1 or -1.	1.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import MattisSpinElement
>>> MattisSpinElement('spin', pattern=-1.0).pattern
-1.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_mattis_spin.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    spin: float = 1.0,
    J_up: float = 1.0,
    J_left: float = 1.0,
    J_right: float = 1.0,
    J_down: float = 1.0,
    h: float = 0.0,
    temperature: float = 1.0,
    pattern: Annotated[float, StateSpec(domain=ZZ)] = 1.0,
) -> None:
    super().__init__(
        label,
        visual=visual,
        spin=spin,
        J_up=J_up,
        J_left=J_left,
        J_right=J_right,
        J_down=J_down,
        h=h,
        temperature=temperature,
    )
    self.pattern = pattern

MaxwellBoltzmannAnalysis(group=None, *, mass=1.0, temperature=None)

Bases: FitAnalysis

Kolmogorov-Smirnov test of the speeds of every frame against the Maxwell-Boltzmann law of their dimension.

It is the fit of the named law 'maxwell' to the speeds, generalized to d dimensions: the speeds of an ideal gas at temperature T follow a chi law of d degrees of freedom and scale sqrt(T / m), the Maxwell law in three dimensions and the Rayleigh law in two. The temperature is the kinetic temperature of every frame unless it is given. It reads the velocities of the agents or particles of a trajectory, or an array of velocities of one snapshot, one row per agent, or of one snapshot per frame along its first axis.

Parameters:

Name	Type	Description	Default
group	str | None	The path of the agent or particle group of a trajectory, or None for the first one.	None
mass	float	The mass of every agent.	1.0
temperature	float | None	The temperature of the law, or None for the kinetic temperature of every frame.	None

Attributes:

Name	Type	Description
values_		The velocities, frames by agents by dimensions.
temperature_		The temperature of the tested law at every frame.
statistic_		The Kolmogorov-Smirnov statistic of every frame.
p_value_		The p-value of every frame.

Examples:

>>> import numpy as np
>>> from skcomplex.models.physics.statistical_mechanics import MaxwellBoltzmannAnalysis
>>> velocities = np.random.default_rng(0).standard_normal((2000, 3))
>>> bool(MaxwellBoltzmannAnalysis().check().prepare(velocities).run().p_value_[0] > 0.05)
True

Source code in src/skcomplex/models/physics/statistical_mechanics/_maxwell_boltzmann.py

def __init__(self, group: str | None = None, *, mass: float = 1.0, temperature: float | None = None) -> None:
    super().__init__("maxwell", group=group)
    self.mass = mass
    self.temperature = temperature

MicrostateElement(label, visual=None, *, P=0.0, U=0.0, beta=1.0, w=0.5)

Bases: Element

State of a discrete master equation on a ring, with a probability P and a potential U.

Every step the probability moves to each neighbouring state B from a state A with the Glauber probability w / (1 + exp(beta (U_B - U_A))), so P_A(t + 1) = P_A + sum_B T_AB P_B - sum_B T_BA P_A, and the stationary distribution is exp(-beta U) / Z.

Parameters:

Name	Type	Description	Default
label	str	The label of the state or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
P	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial probability.	0.0
U	Annotated[float, StateSpec()]	The potential.	0.0
beta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The inverse temperature.	1.0
w	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The attempt probability of a move to each neighbour, at most 1/2.	0.5

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import MicrostateElement
>>> MicrostateElement('state').w
0.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_microstate.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    P: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    U: Annotated[float, StateSpec()] = 0.0,
    beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    w: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.P = P
    self.U = U
    self.beta = beta
    self.w = w

OrderParameterElement(label, visual=None, *, phi=0.1, a=1.0, b=1.0, T=0.5, Tc=1.0, h=0.0, Gamma=1.0)

Bases: Element

Order parameter of Landau theory relaxing down its free energy, dphi/dt = -Gamma dF/dphi.

The free energy is F = a (T - Tc) phi^2 / 2 + b phi^4 / 4 - h phi.

Parameters:

Name	Type	Description	Default
label	str	The label of the order parameter.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
phi	Annotated[float, StateSpec()]	The initial order parameter.	0.1
a	Annotated[float, ParameterSpec(domain=open(0, oo))]	The coefficient of the quadratic term per unit temperature.	1.0
b	Annotated[float, ParameterSpec(domain=open(0, oo))]	The coefficient of the quartic term.	1.0
T	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The temperature.	0.5
Tc	Annotated[float, ParameterSpec(domain=open(0, oo))]	The critical temperature.	1.0
h	Annotated[float, ParameterSpec()]	The external field.	0.0
Gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The kinetic coefficient.	1.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import OrderParameterElement
>>> OrderParameterElement('magnet').Tc
1.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_order_parameter.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    phi: Annotated[float, StateSpec()] = 0.1,
    a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    b: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    T: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
    Tc: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    h: Annotated[float, ParameterSpec()] = 0.0,
    Gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.phi = phi
    self.a = a
    self.b = b
    self.T = T
    self.Tc = Tc
    self.h = h
    self.Gamma = Gamma

ParticleSiteElement(label, visual=None, *, n=0.0)

Bases: Element

Site of a lattice of particles, holding the number n of particles or quanta on it.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
n	Annotated[float, StateSpec(domain=NN)]	The initial number of particles.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ParticleSiteElement
>>> ParticleSiteElement('site', n=3.0).n
3.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_particle_site.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    n: Annotated[float, StateSpec(domain=NN)] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.n = n

PercolationSiteElement(label, visual=None, *, cover=0.0)

Bases: Element

Site of a percolation lattice: empty at 0, occupied at 1, burning at 2, or burnt at 3.

Fire spreads from burning sites to their occupied neighbours and burns out in one step, so a fire lit along one edge burns exactly the clusters that touch that edge.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
cover	Annotated[float, StateSpec(domain=NN)]	The initial cover, 0, 1, 2, or 3.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import PercolationSiteElement
>>> PercolationSiteElement('site').cover
0.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_percolation_site.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    cover: Annotated[float, StateSpec(domain=NN)] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.cover = cover

PottsSpinElement(label, visual=None, *, color=0.0, J=1.0, temperature=1.0)

Bases: Element

Spin of the q-state Potts model, in a state from 0 to q - 1, with its coupling J and its temperature.

A PottsSamplingInteraction redraws the states by Monte Carlo sweeps at the temperature of the spins.

Parameters:

Name	Type	Description	Default
label	str	The label of the spin or the lattice group.	required
visual	VisualSpec | None	How the spin is drawn, or None for default_visual_.	None
color	Annotated[float, StateSpec(domain=NN)]	The initial state.	0.0
J	Annotated[float, ParameterSpec()]	The coupling between neighbouring spins.	1.0
temperature	Annotated[float, ParameterSpec(domain=open(0, oo))]	The temperature in units of the coupling.	1.0

Attributes:

Name	Type	Description
default_visual_		The visual of the spin when it sets none.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import PottsSpinElement
>>> PottsSpinElement('spin').temperature
1.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_potts_spin.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    color: Annotated[float, StateSpec(domain=NN)] = 0.0,
    J: Annotated[float, ParameterSpec()] = 1.0,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.color = color
    self.J = J
    self.temperature = temperature

Q2RSpinElement(label, visual=None, *, spin=1.0, phase=0.0)

Bases: Element

Site of the Q2R automaton, a spin and a phase that says whether the site updates in the current step.

The local rule of Q2R flips the spins of the updating sites whose neighbours sum to zero, which leaves the Ising energy unchanged, and swaps the phase of every site, so the two sublattices update in turn and the dynamics is reversible.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the lattice group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
spin	Annotated[float, StateSpec(domain=ZZ)]	The initial spin, 1 or -1.	1.0
phase	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	The initial phase, 1 for updating and 0 for quiescent.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import Q2RSpinElement
>>> Q2RSpinElement('site').phase
0.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_q2r_spin.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    spin: Annotated[float, StateSpec(domain=ZZ)] = 1.0,
    phase: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.spin = spin
    self.phase = phase

ReactingParticlesElement(label, visual=None, *, n=100.0, kappa=2.0, nu=1.0, lam=1.0, volume=100.0, epsilon=0.0)

Bases: Element

Number n of particles in a volume V, with the reactions A -> 2A, A -> 0, and 2A -> A.

A particle branches at rate kappa, decays at rate nu, and every pair annihilates at rate lambda / V, so W(n + 1 | n) = kappa n + epsilon delta_n0 and W(n - 1 | n) = nu n + lambda n (n - 1) / V. The empty state absorbs when epsilon = 0.

Parameters:

Name	Type	Description	Default
label	str	The label of the particles.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
n	Annotated[float, StateSpec(domain=NN)]	The initial number of particles.	100.0
kappa	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The branching rate.	2.0
nu	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The decay rate.	1.0
lam	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The annihilation rate.	1.0
volume	Annotated[float, ParameterSpec(domain=open(0, oo))]	The volume V.	100.0
epsilon	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate at which a particle appears in the empty state.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ReactingParticlesElement
>>> ReactingParticlesElement('particles').kappa
2.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_absorbing_processes.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    n: Annotated[float, StateSpec(domain=NN)] = 100.0,
    kappa: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
    nu: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    lam: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    volume: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.n = n
    self.kappa = kappa
    self.nu = nu
    self.lam = lam
    self.volume = volume
    self.epsilon = epsilon

ReservoirInteraction(label, *, sites, inflow=0.0, outflow=0.0, exclusion=False)

Bases: TransferInteraction

Reservoir that adds one particle to each of its sites at rate inflow and removes one at rate outflow.

A particle leaves only an occupied site. With exclusion a particle enters only an empty site. It is a transfer of one particle from the environment into every site, with jumps back to the environment at the outflow rate.

Parameters:

Name	Type	Description	Default
label	str	The label of the reservoir.	required
sites	Annotated[list[str], RoleSpec(variadic=True)]	The sites in contact with the reservoir.	required
inflow	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate at which a particle enters a site.	0.0
outflow	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate at which a particle leaves an occupied site.	0.0
exclusion	bool	Whether a site holds at most one particle.	False

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ReservoirInteraction
>>> ReservoirInteraction('bath', sites=['lane[0]'], inflow=0.5).outflow
0.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_reservoir.py

def __init__(
    self,
    label: str,
    *,
    sites: Annotated[list[str], RoleSpec(variadic=True)],
    inflow: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    outflow: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    exclusion: bool = False,
) -> None:
    super().__init__(
        label,
        source='',
        target='',
        state=Source(),
        into='n',
        amount=1.0,
        rate=self._get_entry_rate,
        reverse_rate=self._get_exit_rate,
    )
    self.sites = sites
    self.inflow = inflow
    self.outflow = outflow
    self.exclusion = exclusion

SusceptibilityAnalysis(variables=None, *, group=None, size=1.0, temperatures=None, parameters=None)

Bases: DistributionAnalysis

Susceptibility N (<m²> - <|m|>²) of every series of samples of an order parameter, one per control value.

It is the distribution analysis of kind 'moments' of the absolute samples, whose variance it scales by the system size N and divides by the temperature when temperatures are given. Over a sweep of a control parameter the susceptibility peaks near the transition, and the peak is refined by a parabola through its neighbours.

Parameters:

Name	Type	Description	Default
variables	Sequence[str] | None	The variables of a trajectory to read, or None for every scalar series.	None
group	str | None	The path of the agent or particle group whose speeds are read, or None for the series of a trajectory.	None
size	float	The system size N that scales the variance.	1.0
temperatures	Sequence[float] | None	The temperature of every series that divides the susceptibility, or None to leave it undivided.	None
parameters	Sequence[float] | None	The control parameter value of every series, or None without a peak.	None

Attributes:

Name	Type	Description
values_		The samples, one per column.
labels_		The label of every series.
mean_		The mean absolute order parameter of every series.
variance_		The variance of the absolute order parameter of every series.
skewness_		The skewness of the absolute order parameter of every series.
kurtosis_		The excess kurtosis of the absolute order parameter of every series.
susceptibility_		The susceptibility of every series.
peak_		The control parameter value of the largest susceptibility, refined by a parabola through its neighbours, or None without control parameter values.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import SusceptibilityAnalysis
>>> susceptibility = SusceptibilityAnalysis(size=4).check().prepare([0.5, -1.0, 0.5, -1.0]).run()
>>> susceptibility.susceptibility_.tolist()
[0.25]

Source code in src/skcomplex/models/physics/statistical_mechanics/_susceptibility.py

def __init__(
    self,
    variables: Sequence[str] | None = None,
    *,
    group: str | None = None,
    size: float = 1.0,
    temperatures: Sequence[float] | None = None,
    parameters: Sequence[float] | None = None,
) -> None:
    super().__init__("moments", variables, group=group, absolute=True)
    self.size = size
    self.temperatures = temperatures
    self.parameters = parameters

TrapWalkerElement(label, visual=VisualSpec('points', {'position': 'x', 'color_by': 'E', 'color_range': (0.0, 5.0)}), *, x=0.5, E=1.0, T=0.5, tau0=1.0)

Bases: Element

Walker in a trap of depth E, which escapes after an exponential time of mean tau0 exp(E / T).

On escape the walker falls into a new trap of the fully connected landscape, at a uniformly random place in the unit square, whose depth is drawn from the density exp(-E).

Parameters:

Name	Type	Description	Default
label	str	The label of the walker or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('points', {'position': 'x', 'color_by': 'E', 'color_range': (0.0, 5.0)})
x	Annotated[ArrayLike, StateSpec(vector=True, role=position)]	The place of the trap.	0.5
E	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The depth of the trap.	1.0
T	Annotated[float, ParameterSpec(domain=open(0, oo))]	The temperature, in units of the mean depth.	0.5
tau0	Annotated[float, ParameterSpec(domain=open(0, oo))]	The escape time from a trap of zero depth.	1.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import TrapWalkerElement
>>> TrapWalkerElement('walker', T=0.5).T
0.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_traps.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('points', {'position': 'x', 'color_by': 'E', 'color_range': (0.0, 5.0)}),
    *,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.5,
    E: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    T: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    tau0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.E = E
    self.T = T
    self.tau0 = tau0

TreeNodeElement(label, visual=VisualSpec('points', {'position': 'x', 'radius': 0.04, 'color_by': 'n', 'color_range': (0.0, 2.0), 'states': (), 'trail': 0}), *, n=0.0, x=0.0)

Bases: ParticleSiteElement

Node of a tree that holds a number n of energy quanta, drawn at a fixed position.

Parameters:

Name	Type	Description	Default
label	str	The label of the node or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('points', {'position': 'x', 'radius': 0.04, 'color_by': 'n', 'color_range': (0.0, 2.0), 'states': (), 'trail': 0})
n	float	The initial number of quanta.	0.0
x	Annotated[ArrayLike, StateSpec(vector=True)]	The position of the node in the drawing of the tree.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import TreeNodeElement
>>> TreeNodeElement('root').n
0.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_tree_node.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec(
        'points',
        {'position': 'x', 'radius': 0.04, 'color_by': 'n', 'color_range': (0.0, 2.0), 'states': (), 'trail': 0},
    ),
    *,
    n: float = 0.0,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True)] = 0.0,
) -> None:
    super().__init__(label, visual=visual, n=n)
    self.x = x

WalkerElement(label, visual=VisualSpec('agents', {'color_by': 'stuck', 'color_range': (0.0, 1.0)}), *, stuck=0.0)

Bases: Element

Random walker of diffusion-limited aggregation, which sticks when it reaches a stuck walker.

Parameters:

Name	Type	Description	Default
label	str	The label of the agent or group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	VisualSpec('agents', {'color_by': 'stuck', 'color_range': (0.0, 1.0)})
stuck	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))]	Whether the walker belongs to the aggregate.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import WalkerElement
>>> WalkerElement('walker').stuck
0.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_walker.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = VisualSpec('agents', {'color_by': 'stuck', 'color_range': (0.0, 1.0)}),
    *,
    stuck: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.stuck = stuck

walk(agent, model)

Stick next to the aggregate, or step to a random neighbouring cell.

Parameters:

Name	Type	Description	Default
agent	WalkerView	The walker.	required
model	ModelView	The model.	required

Source code in src/skcomplex/models/physics/statistical_mechanics/_walker.py

def walk(self, agent: WalkerView, model: ModelView) -> None:
    """Stick next to the aggregate, or step to a random neighbouring cell.

    Args:
        agent:
            The walker.

        model:
            The model.
    """
    if agent.stuck:
        return
    if any(isinstance(other, WalkerView) and other.stuck for other in agent.neighbors()):
        agent.stuck = 1.0
        return
    cells = agent.neighbor_cells()
    agent.move_to(cells[int(model.rng.integers(len(cells)))])

WalkerView(agent, stepper)

Bases: AgentView

Walker of diffusion-limited aggregation as rules see it, with whether it belongs to the aggregate.

Parameters:

Name	Type	Description	Default
agent	Agent	The Mesa agent.	required
stepper	MesaSolver	The stepper that runs the model.	required

Source code in src/skcomplex/engine/solvers/_mesa.py

def __init__(self, agent: mesa.Agent, stepper: 'MesaSolver') -> None:
    self._agent = agent
    self._stepper = stepper

XXZExchangeInteraction(label, *, target, J=1.0, delta=1.0)

Bases: NeighborInteraction

Exchange J (S^x_i S^x_j + S^y_i S^y_j + Δ S^z_i S^z_j) between neighbouring spins.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the spin chain.	required
J	Annotated[float, ParameterSpec()]	The exchange coupling, antiferromagnetic when positive.	1.0
delta	Annotated[float, ParameterSpec()]	The anisotropy Δ, 1 for the isotropic chain and 0 for the XX chain.	1.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import XXZExchangeInteraction
>>> XXZExchangeInteraction('exchange', target='chain', delta=0.5).delta
0.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_xxz_exchange.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    J: Annotated[float, ParameterSpec()] = 1.0,
    delta: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label)
    self.target = target
    self.J = J
    self.delta = delta

XYSpinElement(label, visual=None, *, theta=0.0, J=1.0, temperature=0.5)

Bases: Element

Planar spin of the XY model, an angle theta, with its coupling J and its temperature.

An XYSamplingInteraction rotates the angles by Monte Carlo sweeps at the temperature of the spins.

Parameters:

Name	Type	Description	Default
label	str	The label of the spin or the lattice group.	required
visual	VisualSpec | None	How the spin is drawn, or None for default_visual_.	None
theta	Annotated[float, StateSpec(period=2 * pi)]	The initial angle.	0.0
J	Annotated[float, ParameterSpec()]	The coupling between neighbouring spins.	1.0
temperature	Annotated[float, ParameterSpec(domain=open(0, oo))]	The temperature in units of the coupling.	0.5

Attributes:

Name	Type	Description
default_visual_		The visual of the spin when it sets none.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import XYSpinElement
>>> XYSpinElement('spin').temperature
0.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_xy_spin.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    theta: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
    J: Annotated[float, ParameterSpec()] = 1.0,
    temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
) -> None:
    super().__init__(label, visual=visual)
    self.theta = theta
    self.J = J
    self.temperature = temperature

ZeemanFieldInteraction(label, *, target, h=0.0)

Bases: Interaction

Uniform field -h Σ S^z_i on the spins of a chain.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the spin chain.	required
h	Annotated[float, ParameterSpec()]	The field.	0.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ZeemanFieldInteraction
>>> ZeemanFieldInteraction('field', target='chain', h=0.2).h
0.2

Source code in src/skcomplex/models/physics/statistical_mechanics/_zeeman_field.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    h: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label)
    self.target = target
    self.h = h

ZeroRangeHopInteraction(label, *, link, p=1.0, q=0.0, b=4.0)

Bases: TransferInteraction

Hop of one particle along a link, forward at rate p u(n) and backward at rate q u(n').

The rate depends only on the occupation of the departure site, through u(n) = 1 + b / n. It is a transfer of one particle from the first site to the second, with jumps back at the reverse rate.

Parameters:

Name	Type	Description	Default
label	str	The label of the link.	required
link	Annotated[list[str], RoleSpec(variadic=True)]	The departure and the arrival site of a forward hop.	required
p	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The forward rate scale.	1.0
q	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The backward rate scale.	0.0
b	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the decrease of u(n) with the occupation, 0 for a constant rate.	4.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ZeroRangeHopInteraction
>>> ZeroRangeHopInteraction('hop', link=['sites[0]', 'sites[1]']).b
4.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_zero_range_hop.py

def __init__(
    self,
    label: str,
    *,
    link: Annotated[list[str], RoleSpec(variadic=True)],
    p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    q: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    b: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
) -> None:
    super().__init__(
        label,
        source='',
        target='',
        state='n',
        amount=1.0,
        rate=self._get_forward_rate,
        reverse_rate=self._get_backward_rate,
    )
    self.link = link
    self.p = p
    self.q = q
    self.b = b

ZeroRangeHopsInteraction(label, *, sites, p=1.0, q=0.0, b=4.0)

Bases: TransferInteraction

Hops of one particle from every site of a ring, to the right at rate p u(n) and to the left at rate q u(n).

The rate depends only on the occupation n of the departure site, through u(n) = 1 + b / n. It is a transfer of one particle along every edge of the ring, to the right and back.

Parameters:

Name	Type	Description	Default
label	str	The label of the hops.	required
sites	Annotated[str, RoleSpec()]	The group of sites, in ring order.	required
p	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate scale of hops to the right.	1.0
q	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate scale of hops to the left.	0.0
b	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the decrease of u(n) with the occupation, 0 for a constant rate.	4.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ZeroRangeHopsInteraction
>>> ZeroRangeHopsInteraction('hops', sites='sites').b
4.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_zero_range_hops.py

def __init__(
    self,
    label: str,
    *,
    sites: Annotated[str, RoleSpec()],
    p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    q: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    b: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
) -> None:
    super().__init__(
        label,
        source=sites,
        target=sites,
        state='n',
        amount=1.0,
        rate=self._get_right_rate,
        reverse_rate=self._get_left_rate,
        offset=1,
    )
    self.sites = sites
    self.p = p
    self.q = q
    self.b = b

ZoneElement(label, visual=None, *, N=0.0, Na=100.0, p=0.1, k=1.0)

Bases: Element

Small zone of a two-phase apparatus that exchanges particles with the rest of it.

Every one of the Na particles enters the zone at the rate k p and leaves it at the rate k (1 - p), where p = dV / V is the volume fraction of the zone, so the number in the zone is binomial with mean Na p.

Parameters:

Name	Type	Description	Default
label	str	The label of the zone or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
N	Annotated[float, StateSpec(domain=NN)]	The initial number of particles in the zone.	0.0
Na	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The number of particles in the apparatus and the zone.	100.0
p	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The volume fraction of the zone.	0.1
k	Annotated[float, ParameterSpec(domain=open(0, oo))]	The rate of exchange.	1.0

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ZoneElement
>>> ZoneElement('zone').Na
100.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_zone.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    N: Annotated[float, StateSpec(domain=NN)] = 0.0,
    Na: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 100.0,
    p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.N = N
    self.Na = Na
    self.p = p
    self.k = k

absorbing_birth_death(*, kappa=2.0, nu=1.0, lam=1.0, volume=100.0, density=1.0, epsilon=0.0)

Return the fully connected branching and annihilation process, active above kappa = nu until it dies out.

The deterministic density obeys d rho / dt = (kappa - nu) rho - lambda rho^2, with the active fixed point rho_0 = (kappa - nu) / lambda. With fluctuations the empty state absorbs every run, but above the transition the active state is metastable, with relative fluctuations sigma / rho_0 = (lambda / (kappa - nu)) sqrt(kappa / (lambda V)) and a lifetime that grows like exp(V phi(0)) / nu. Below the transition the particles die out in a time of order 1 / (nu - kappa). A small epsilon makes the empty state leave, as in exercise 3.3 of Bertin.

Parameters:

Name	Type	Description	Default
kappa	float	The branching rate.	2.0
nu	float	The decay rate.	1.0
lam	float	The annihilation rate.	1.0
volume	float	The volume V.	100.0
density	float	The initial density, rounded to a whole number of particles.	1.0
epsilon	float	The rate at which a particle appears in the empty state.	0.0

Returns:

Name	Type	Description
system	System	The model, with the number of particles particles.n.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import absorbing_birth_death
>>> absorbing_birth_death(volume=50.0).elements[0].n
50.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_absorbing_processes.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('absorbing state', 'phase transition', 'birth-death process', 'metastability', 'jump process'),
    references=('Bertin 2021', 'Hinrichsen 2000'),
    defaults={'until': 100.0, 'dt': 0.1},
)
def absorbing_birth_death(
    *,
    kappa: float = 2.0,
    nu: float = 1.0,
    lam: float = 1.0,
    volume: float = 100.0,
    density: float = 1.0,
    epsilon: float = 0.0,
) -> System:
    """Return the fully connected branching and annihilation process, active above `kappa = nu` until it dies out.

    The deterministic density obeys `d rho / dt = (kappa - nu) rho - lambda rho^2`, with the active fixed point
    `rho_0 = (kappa - nu) / lambda`. With fluctuations the empty state absorbs every run, but above the transition the
    active state is metastable, with relative fluctuations `sigma / rho_0 = (lambda / (kappa - nu)) sqrt(kappa /
    (lambda V))` and a lifetime that grows like `exp(V phi(0)) / nu`. Below the transition the particles die out in a
    time of order `1 / (nu - kappa)`. A small `epsilon` makes the empty state leave, as in exercise 3.3 of Bertin.

    Args:
        kappa:
            The branching rate.

        nu:
            The decay rate.

        lam:
            The annihilation rate.

        volume:
            The volume `V`.

        density:
            The initial density, rounded to a whole number of particles.

        epsilon:
            The rate at which a particle appears in the empty state.

    Returns:
        system:
            The model, with the number of particles `particles.n`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import absorbing_birth_death
        >>> absorbing_birth_death(volume=50.0).elements[0].n
        50.0
    """
    particles = ReactingParticlesElement(
        'particles',
        n=float(round(density * volume)),
        kappa=kappa,
        nu=nu,
        lam=lam,
        volume=volume,
        epsilon=epsilon,
    )
    return System(elements=[particles])

asep(n_sites=100, *, alpha=0.3, gamma=0.7, p=1.0, q=0.0, density=0.0, random_state=0)

Return the asymmetric simple exclusion process on an open line, in continuous time.

A particle hops to the empty right neighbour at rate p and to the empty left neighbour at rate q, enters the empty first site at rate alpha, and leaves the last site at rate gamma, the exit rate that Bertin writes gamma and others write beta. For p = 1 and q = 0 the current is alpha (1 - alpha) in the low-density phase alpha < 1/2, gamma > alpha, gamma (1 - gamma) in the high-density phase gamma < 1/2, gamma < alpha, and 1/4 in the maximal-current phase alpha, gamma > 1/2.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of sites.	100
alpha	float	The entry rate.	0.3
gamma	float	The exit rate.	0.7
p	float	The rate of hops to the right.	1.0
q	float	The rate of hops to the left.	0.0
density	float	The initial fraction of occupied sites.	0.0
random_state	int | Generator | None	The seed or generator of the initial occupations.	0

Returns:

Name	Type	Description
system	System	The model, with the occupations lane.n.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import asep
>>> len(asep(10).interactions)
3

Source code in src/skcomplex/models/physics/statistical_mechanics/_driven_particles.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('exclusion', 'nonequilibrium', 'phase transition', 'traffic', 'jump process'),
    references=('Derrida, Evans, Hakim, and Pasquier 1993', 'Bertin 2021'),
    defaults={'until': 500.0, 'dt': 2.0, 'backend': NumbaBackend()},
)
def asep(
    n_sites: int = 100,
    *,
    alpha: float = 0.3,
    gamma: float = 0.7,
    p: float = 1.0,
    q: float = 0.0,
    density: float = 0.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the asymmetric simple exclusion process on an open line, in continuous time.

    A particle hops to the empty right neighbour at rate `p` and to the empty left neighbour at rate `q`, enters the
    empty first site at rate `alpha`, and leaves the last site at rate `gamma`, the exit rate that Bertin writes
    `gamma` and others write `beta`. For `p = 1` and `q = 0` the current is `alpha (1 - alpha)` in the low-density
    phase `alpha < 1/2, gamma > alpha`, `gamma (1 - gamma)` in the high-density phase `gamma < 1/2, gamma < alpha`, and
    `1/4` in the maximal-current phase `alpha, gamma > 1/2`.

    Args:
        n_sites:
            The number of sites.

        alpha:
            The entry rate.

        gamma:
            The exit rate.

        p:
            The rate of hops to the right.

        q:
            The rate of hops to the left.

        density:
            The initial fraction of occupied sites.

        random_state:
            The seed or generator of the initial occupations.

    Returns:
        system:
            The model, with the occupations `lane.n`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import asep
        >>> len(asep(10).interactions)
        3
    """
    occupied = (np.random.default_rng(random_state).random(n_sites) < density).astype(float)
    lane = LatticeElement(
        'lane',
        element=ParticleSiteElement,
        space=LatticeSpace(n_sites, boundary='fixed'),
        values={'n': occupied},
    )
    interactions = [
        ExclusionHopsInteraction('hops', sites='lane', p=p, q=q),
        ReservoirInteraction('entry', sites=['lane[0]'], inflow=alpha, exclusion=True),
        ReservoirInteraction('exit', sites=[f'lane[{n_sites - 1}]'], outflow=gamma, exclusion=True),
    ]
    return System(elements=[lane], interactions=interactions)

contact_process(shape=(64, 64), *, rate=3.0, tau=0.1, infected=0.5, random_state=0)

Return the contact process on a periodic lattice, which survives above a critical infection rate near 1.65.

Below the critical rate the infection dies out and the lattice falls into the absorbing healthy state.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns.	(64, 64)
rate	float	The infection rate.	3.0
tau	float	The time slice of one step.	0.1
infected	float	The initial fraction of infected sites.	0.5
random_state	int | Generator | None	The seed or generator of the initial infections.	0

Returns:

Name	Type	Description
system	System	The model, with the infections lattice.infected.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import contact_process
>>> contact_process((10, 10), rate=1.0).elements[0].values['rate']
1.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_processes.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('lattice', 'stochastic', 'absorbing state', 'phase transition', 'epidemics'),
    references=('Harris 1974', 'Marro and Dickman 1999', 'Pastor-Satorras et al. 2015'),
    defaults={'until': 200, 'dt': 1},
)
def contact_process(
    shape: Sequence[int] = (64, 64),
    *,
    rate: float = 3.0,
    tau: float = 0.1,
    infected: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the contact process on a periodic lattice, which survives above a critical infection rate near 1.65.

    Below the critical rate the infection dies out and the lattice falls into the absorbing healthy state.

    Args:
        shape:
            The number of rows and columns.

        rate:
            The infection rate.

        tau:
            The time slice of one step.

        infected:
            The initial fraction of infected sites.

        random_state:
            The seed or generator of the initial infections.

    Returns:
        system:
            The model, with the infections `lattice.infected`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import contact_process
        >>> contact_process((10, 10), rate=1.0).elements[0].values['rate']
        1.0
    """
    shape = tuple(shape)
    states = (np.random.default_rng(random_state).random(shape) < infected).astype(float)
    space = LatticeSpace(shape)
    return System(
        elements=[
            LatticeElement(
                'lattice',
                element=ContactSiteElement,
                space=space,
                values={'infected': states, 'rate': rate, 'tau': tau},
            ),
        ],
    )

continuum_percolation(filling=1.1, *, side=60.0, radius=1.0, random_state=0)

Return overlapping disks at uniform positions in a square, whose clusters percolate near a filling of 1.13.

The number of disks is filling side^2 / (pi radius^2). The clusters of touching disks are read by ClusterAnalysis(radius=radius), and an Experiment sweeps filling and random_state for the spanning probability (Satz 2022 ch. 3).

Parameters:

Name	Type	Description	Default
filling	float	The filling factor, the total area of the disks over the area of the square.	1.1
side	float	The side of the square.	60.0
radius	float	The radius of a disk.	1.0
random_state	int | Generator | None	The seed or generator of the positions.	0

Returns:

Name	Type	Description
system	System	The model, with the positions disks.x.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import continuum_percolation
>>> continuum_percolation(1.0, side=10.0).elements[0].size
32

Source code in src/skcomplex/models/physics/statistical_mechanics/_percolation.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('percolation', 'disorder', 'clusters', 'static'),
    references=('Satz 2022', 'Mertens and Moore 2012'),
)
def continuum_percolation(
    filling: float = 1.1,
    *,
    side: float = 60.0,
    radius: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return overlapping disks at uniform positions in a square, whose clusters percolate near a filling of 1.13.

    The number of disks is `filling side^2 / (pi radius^2)`. The clusters of touching disks are read by
    `ClusterAnalysis(radius=radius)`, and an `Experiment` sweeps `filling` and `random_state` for the spanning
    probability (Satz 2022 ch. 3).

    Args:
        filling:
            The filling factor, the total area of the disks over the area of the square.

        side:
            The side of the square.

        radius:
            The radius of a disk.

        random_state:
            The seed or generator of the positions.

    Returns:
        system:
            The model, with the positions `disks.x`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import continuum_percolation
        >>> continuum_percolation(1.0, side=10.0).elements[0].size
        32
    """
    size = round(filling * side**2 / (np.pi * radius**2))
    positions = np.random.default_rng(random_state).uniform(0.0, side, (size, 2))
    visual = VisualSpec('points', {'position': 'x', 'radius': radius, 'states': ()})
    disks = GroupElement('disks', element=DiskElement, size=size, values={'x': positions}, visual=visual)
    return System(elements=[disks], space=ContinuousSpace(extent=[side, side]), time=StaticTime())

curie_weiss(n_spins=500, *, temperature=0.5, J=1.0, h=0.0, initial='random', random_state=0)

Return the fully connected Ising model with Glauber dynamics, ordered below the critical temperature T_c = J.

The energy is -(J / N) sum_{i<j} s_i s_j - h sum_i s_i. In the large system limit the magnetization m solves m = tanh((J m + h) / T), which has a nonzero root below T_c when h = 0. One unit of time is one attempted flip per spin. The spins are drawn in a row, which orders the drawing only.

Parameters:

Name	Type	Description	Default
n_spins	int	The number of spins.	500
temperature	float	The temperature, which can change during a run.	0.5
J	float	The coupling.	1.0
h	float	The external field.	0.0
initial	str	'random' for random spins or 'up' for all spins up.	'random'
random_state	int | Generator | None	The seed or generator of the random initial spins.	0

Returns:

Name	Type	Description
system	System	The model, with the spins spins.spin and the observable magnetization.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import curie_weiss
>>> curie_weiss(10).elements[0].check().prepare().size_
10

Source code in src/skcomplex/models/physics/statistical_mechanics/_fully_connected_spins.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('mean field', 'phase transition', 'spins', 'jump process'),
    references=('Bertin 2021', 'Glauber 1963'),
    defaults={'until': 30.0, 'dt': 0.1},
    observables=lambda system: [_magnetization()],
)
def curie_weiss(
    n_spins: int = 500,
    *,
    temperature: float = 0.5,
    J: float = 1.0,
    h: float = 0.0,
    initial: str = 'random',
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the fully connected Ising model with Glauber dynamics, ordered below the critical temperature `T_c = J`.

    The energy is `-(J / N) sum_{i<j} s_i s_j - h sum_i s_i`. In the large system limit the magnetization `m` solves
    `m = tanh((J m + h) / T)`, which has a nonzero root below `T_c` when `h = 0`. One unit of time is one attempted
    flip per spin. The spins are drawn in a row, which orders the drawing only.

    Args:
        n_spins:
            The number of spins.

        temperature:
            The temperature, which can change during a run.

        J:
            The coupling.

        h:
            The external field.

        initial:
            `'random'` for random spins or `'up'` for all spins up.

        random_state:
            The seed or generator of the random initial spins.

    Returns:
        system:
            The model, with the spins `spins.spin` and the observable `magnetization`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import curie_weiss
        >>> curie_weiss(10).elements[0].check().prepare().size_
        10
    """
    spins = _initial_spins(n_spins, initial, random_state)
    group = LatticeElement(
        'spins',
        element=FullyConnectedSpinElement,
        space=LatticeSpace(n_spins),
        values={'spin': spins, 'J': J, 'h': h, 'temperature': temperature},
        visual=SPIN_IMAGE,
    )
    return System(elements=[group])

diffusion_limited_aggregation(n_walkers=400, *, size=41, random_state=0)

Return diffusion-limited aggregation on a periodic grid, where random walkers stick to a seed at the center.

Walkers that touch the aggregate join it, so the tips of the aggregate catch most walkers and it grows into a branched fractal of dimension about 1.7.

Parameters:

Name	Type	Description	Default
n_walkers	int	The number of walkers, the seed included.	400
size	int	The side of the grid.	41
random_state	int | Generator | None	The seed or generator of the initial cells of the walkers.	0

Returns:

Name	Type	Description
system	System	The model, with the agent group walkers.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import diffusion_limited_aggregation
>>> diffusion_limited_aggregation(10, size=11).elements[0].check().prepare().size_
10

Source code in src/skcomplex/models/physics/statistical_mechanics/_processes.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('agents', 'aggregation', 'fractal', 'random walk'),
    references=('Witten and Sander 1981', 'Charbonneau 2017'),
    defaults={'until': 200, 'dt': 1},
)
def diffusion_limited_aggregation(
    n_walkers: int = 400,
    *,
    size: int = 41,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return diffusion-limited aggregation on a periodic grid, where random walkers stick to a seed at the center.

    Walkers that touch the aggregate join it, so the tips of the aggregate catch most walkers and it grows into a
    branched fractal of dimension about 1.7.

    Args:
        n_walkers:
            The number of walkers, the seed included.

        size:
            The side of the grid.

        random_state:
            The seed or generator of the initial cells of the walkers.

    Returns:
        system:
            The model, with the agent group `walkers`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import diffusion_limited_aggregation
        >>> diffusion_limited_aggregation(10, size=11).elements[0].check().prepare().size_
        10
    """
    rng = np.random.default_rng(random_state)
    middle = size // 2
    cells = rng.integers(0, size, (n_walkers, 2))
    cells[0] = middle
    stuck = np.zeros(n_walkers)
    stuck[0] = 1.0
    walkers = GroupElement(
        'walkers',
        element=WalkerElement,
        size=n_walkers,
        values={'stuck': stuck},
        location=[tuple(map(int, cell)) for cell in cells],
    )
    return System(
        elements=[walkers],
        space=LatticeSpace((size, size), boundary='periodic', neighborhood='von_neumann'),
    )

discrete_master_equation(n_states=50, *, beta=1.0, depth=2.0, w=0.5)

Return the discrete-time master equation of states on a ring, which relaxes to the Boltzmann distribution.

The potential is U = depth cos(2 pi i / n_states), the probability starts on the state of highest potential, and moves to neighbouring states with Glauber probabilities, which satisfy detailed balance.

Parameters:

Name	Type	Description	Default
n_states	int	The number of states.	50
beta	float	The inverse temperature.	1.0
depth	float	The amplitude of the potential.	2.0
w	float	The attempt probability of a move to each neighbour, at most 1/2.	0.5

Returns:

Name	Type	Description
system	System	The model, with the probabilities states.P and potentials states.U.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import discrete_master_equation
>>> discrete_master_equation(10).elements[0].check().prepare().size_
10

Source code in src/skcomplex/models/physics/statistical_mechanics/_master_equations.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('master equation', 'probability', 'detailed balance', 'lattice'),
    references=('Jensen 2023', 'van Kampen 2007'),
    defaults={'until': 500, 'dt': 10},
)
def discrete_master_equation(
    n_states: int = 50,
    *,
    beta: float = 1.0,
    depth: float = 2.0,
    w: float = 0.5,
) -> System:
    """Return the discrete-time master equation of states on a ring, which relaxes to the Boltzmann distribution.

    The potential is `U = depth cos(2 pi i / n_states)`, the probability starts on the state of highest potential,
    and moves to neighbouring states with Glauber probabilities, which satisfy detailed balance.

    Args:
        n_states:
            The number of states.

        beta:
            The inverse temperature.

        depth:
            The amplitude of the potential.

        w:
            The attempt probability of a move to each neighbour, at most 1/2.

    Returns:
        system:
            The model, with the probabilities `states.P` and potentials `states.U`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import discrete_master_equation
        >>> discrete_master_equation(10).elements[0].check().prepare().size_
        10
    """
    potential = depth * np.cos(2 * np.pi * np.arange(n_states) / n_states)
    probability = np.zeros(n_states)
    probability[0] = 1.0
    space = LatticeSpace(n_states)
    return System(
        elements=[
            LatticeElement(
                'states',
                element=MicrostateElement,
                space=space,
                values={'P': probability, 'U': potential, 'beta': beta, 'w': w},
            ),
        ],
    )

dissipative_zero_range(n_levels=5, *, branches=2, nu=1.0, alpha=0.5, beta_ext=1.0, epsilon0=1.0, dissipation=1.0)

Return the dissipative zero-range process on a tree, with energy injected at the root and lost at the leaves.

Level j of the tree has m^(j - 1) nodes, each with m branches to the next level. A quantum epsilon0 of energy crosses a branch between levels j and j + 1 at rate nu_j = nu k_j^alpha, k_j = m^(j - 1), in either direction from an occupied node. A reservoir at the inverse temperature beta_ext exchanges quanta with the root, adding one at rate nu exp(-beta_ext epsilon0) and removing one at rate nu, and every leaf loses a quantum at the rate dissipation. The stationary distribution is a product of geometric laws z_j^n, with z_j the solution of the linear balance of the currents nu_j (z_j - z_(j + 1)) at every level, and the injected flux is nu (exp(-beta_ext epsilon0) - z_1).

Parameters:

Name	Type	Description	Default
n_levels	int	The number of levels M.	5
branches	int	The number m of branches of every node above the leaves.	2
nu	float	The rate scale of transfers and of the exchange with the reservoir.	1.0
alpha	float	The exponent of the growth of the transfer rate with the pseudo wave number k_j.	0.5
beta_ext	float	The inverse temperature of the reservoir.	1.0
epsilon0	float	The quantum of energy.	1.0
dissipation	float	The rate at which every leaf loses a quantum.	1.0

Returns:

Name	Type	Description
system	System	The model, with the quanta nodes.n of the nodes in breadth-first order, from the root to the leaves.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import dissipative_zero_range
>>> system = dissipative_zero_range(3, branches=3)
>>> system.elements[0].size, len(system.interactions)
(13, 14)

Source code in src/skcomplex/models/physics/statistical_mechanics/_driven_particles.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('zero-range process', 'energy cascade', 'nonequilibrium', 'tree', 'jump process'),
    references=('Bertin 2021', 'Bertin 2005'),
    defaults={'until': 200.0, 'dt': 1.0, 'backend': NumbaBackend()},
)
def dissipative_zero_range(
    n_levels: int = 5,
    *,
    branches: int = 2,
    nu: float = 1.0,
    alpha: float = 0.5,
    beta_ext: float = 1.0,
    epsilon0: float = 1.0,
    dissipation: float = 1.0,
) -> System:
    """Return the dissipative zero-range process on a tree, with energy injected at the root and lost at the leaves.

    Level `j` of the tree has `m^(j - 1)` nodes, each with `m` branches to the next level. A quantum `epsilon0` of
    energy crosses a branch between levels `j` and `j + 1` at rate `nu_j = nu k_j^alpha`, `k_j = m^(j - 1)`, in either
    direction from an occupied node. A reservoir at the inverse temperature `beta_ext` exchanges quanta with the root,
    adding one at rate `nu exp(-beta_ext epsilon0)` and removing one at rate `nu`, and every leaf loses a quantum at the
    rate `dissipation`. The stationary distribution is a product of geometric laws `z_j^n`, with `z_j` the solution of
    the linear balance of the currents `nu_j (z_j - z_(j + 1))` at every level, and the injected flux is
    `nu (exp(-beta_ext epsilon0) - z_1)`.

    Args:
        n_levels:
            The number of levels `M`.

        branches:
            The number `m` of branches of every node above the leaves.

        nu:
            The rate scale of transfers and of the exchange with the reservoir.

        alpha:
            The exponent of the growth of the transfer rate with the pseudo wave number `k_j`.

        beta_ext:
            The inverse temperature of the reservoir.

        epsilon0:
            The quantum of energy.

        dissipation:
            The rate at which every leaf loses a quantum.

    Returns:
        system:
            The model, with the quanta `nodes.n` of the nodes in breadth-first order, from the root to the leaves.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import dissipative_zero_range
        >>> system = dissipative_zero_range(3, branches=3)
        >>> system.elements[0].size, len(system.interactions)
        (13, 14)
    """
    n_nodes = _level_start(n_levels + 1, branches)
    nodes = GroupElement('nodes', element=TreeNodeElement, size=n_nodes, values={'x': _tree_layout(n_levels, branches)})
    interactions: list[BaseInteraction] = []
    for level in range(1, n_levels):
        rate = nu * float(branches ** (level - 1)) ** alpha
        start, below = _level_start(level, branches), _level_start(level + 1, branches)
        for index in range(branches ** (level - 1)):
            for branch in range(branches):
                child = below + branches * index + branch
                link = [f'nodes[{start + index}]', f'nodes[{child}]']
                interactions.append(
                    ZeroRangeHopInteraction(f'transfer_{start + index}_{child}', link=link, p=rate, q=rate, b=0.0),
                )
    leaves = [f'nodes[{index}]' for index in range(_level_start(n_levels, branches), n_nodes)]
    interactions += [
        ReservoirInteraction('injection', sites=['nodes[0]'], inflow=nu * np.exp(-beta_ext * epsilon0), outflow=nu),
        ReservoirInteraction('dissipation', sites=leaves, outflow=dissipation),
    ]
    return System(elements=[nodes], interactions=interactions, space=ContinuousSpace(n_dim=2))

edwards_anderson(size=48, *, temperature=0.8, J=1.0, h=0.0, disorder='bimodal', random_state=0)

Return the Edwards-Anderson spin glass on a periodic square lattice, with quenched random bonds.

Every bond takes the coupling +J or -J with equal probability, or a Gaussian coupling of mean 0 and standard deviation J, drawn once. The magnetization stays near zero at every temperature, while at low temperature the spins freeze in random directions and the energy per spin approaches the ground-state energy, about -1.40 J for +J and -J bonds. One unit of time is one sweep of Metropolis updates.

Parameters:

Name	Type	Description	Default
size	int | Sequence[int]	The side of the square lattice, or its two sides.	48
temperature	float	The temperature, which can change during a run.	0.8
J	float	The scale of the couplings.	1.0
h	float	The external field.	0.0
disorder	str	'bimodal' for couplings +J and -J, or 'gaussian' for Gaussian couplings.	'bimodal'
random_state	int | Generator | None	The seed or generator of the couplings and of the random initial spins.	0

Returns:

Name	Type	Description
system	System	The model, with the spins spins.spin, the couplings spins.J_right and spins.J_down of the bonds to the right and down, and the observable magnetization.

Raises:

Type	Description
ValidationError	The disorder is neither 'bimodal' nor 'gaussian'.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import edwards_anderson
>>> sorted(set(edwards_anderson(8).elements[0].values['J_right'].ravel().tolist()))
[-1.0, 1.0]

Source code in src/skcomplex/models/physics/statistical_mechanics/_spin_glasses.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('spin glass', 'disorder', 'lattice', 'Monte Carlo', 'spins'),
    references=('Edwards and Anderson 1975', 'Bertin 2021', 'Binder and Young 1986'),
    defaults={'until': 500, 'dt': 5},
    observables=_magnetization_observables,
)
def edwards_anderson(
    size: int | Sequence[int] = 48,
    *,
    temperature: float = 0.8,
    J: float = 1.0,
    h: float = 0.0,
    disorder: str = 'bimodal',
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Edwards-Anderson spin glass on a periodic square lattice, with quenched random bonds.

    Every bond takes the coupling `+J` or `-J` with equal probability, or a Gaussian coupling of mean 0 and standard
    deviation `J`, drawn once. The magnetization stays near zero at every temperature, while at low temperature the
    spins freeze in random directions and the energy per spin approaches the ground-state energy, about `-1.40 J` for
    `+J` and `-J` bonds. One unit of time is one sweep of Metropolis updates.

    Args:
        size:
            The side of the square lattice, or its two sides.

        temperature:
            The temperature, which can change during a run.

        J:
            The scale of the couplings.

        h:
            The external field.

        disorder:
            `'bimodal'` for couplings `+J` and `-J`, or `'gaussian'` for Gaussian couplings.

        random_state:
            The seed or generator of the couplings and of the random initial spins.

    Returns:
        system:
            The model, with the spins `spins.spin`, the couplings `spins.J_right` and `spins.J_down` of the bonds to the
            right and down, and the observable `magnetization`.

    Raises:
        ValidationError:
            The disorder is neither `'bimodal'` nor `'gaussian'`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import edwards_anderson
        >>> sorted(set(edwards_anderson(8).elements[0].values['J_right'].ravel().tolist()))
        [-1.0, 1.0]
    """
    shape = _square(size)
    rng = np.random.default_rng(random_state)
    if disorder == 'bimodal':
        right, down = J * rng.choice([-1.0, 1.0], size=(2, *shape))
    elif disorder == 'gaussian':
        right, down = J * rng.standard_normal((2, *shape))
    else:
        error_msg = f"The disorder is 'bimodal' or 'gaussian', not {disorder!r}."
        raise ValidationError(error_msg)
    spins = rng.choice([-1.0, 1.0], size=shape)
    group = LatticeElement(
        'spins',
        element=BondSpinElement,
        space=LatticeSpace(shape),
        values={'spin': spins, 'h': h, 'temperature': temperature, **_bond_states(right, down)},
        visual=SPIN_IMAGE,
    )
    return System(elements=[group], interactions=[BondSamplingInteraction('sampling', target='spins')])

ehrenfest_urn(n_balls=40, *, in_first=None)

Return the Ehrenfest urn model, balls moved one at a time between two boxes, which relaxes to equal counts.

At every step one of the N numbered balls is drawn at random and moved to the other box, so the count n of the first box falls by one with probability n / N and rises by one otherwise. Started with all balls in the first box, the count relaxes to N / 2 with the exponential mean N / 2 + (n_0 - N / 2) (1 - 2 / N)^t, and then fluctuates around it: the chain is reversible, yet the relaxation looks irreversible. Its stationary distribution is the binomial distribution of N balls with probability one half.

Parameters:

Name	Type	Description	Default
n_balls	int	The number of balls N.	40
in_first	int | None	The initial number of balls in the first box, or None for all of them.	None

Returns:

Name	Type	Description
system	System	The model, with the count urn.n of balls in the first box.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ehrenfest_urn
>>> ehrenfest_urn(10).elements[0].n
10

Source code in src/skcomplex/models/physics/statistical_mechanics/_urns.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('markov chain', 'irreversibility', 'entropy', 'discrete map'),
    references=('Ehrenfest and Ehrenfest 1907', 'Barsky 2020', 'Kac 1947'),
    defaults={'until': 200, 'dt': 1},
)
def ehrenfest_urn(n_balls: int = 40, *, in_first: int | None = None) -> System:
    """Return the Ehrenfest urn model, balls moved one at a time between two boxes, which relaxes to equal counts.

    At every step one of the `N` numbered balls is drawn at random and moved to the other box, so the count `n` of the
    first box falls by one with probability `n / N` and rises by one otherwise. Started with all balls in the first
    box, the count relaxes to `N / 2` with the exponential mean `N / 2 + (n_0 - N / 2) (1 - 2 / N)^t`, and then
    fluctuates around it: the chain is reversible, yet the relaxation looks irreversible. Its stationary distribution
    is the binomial distribution of `N` balls with probability one half.

    Args:
        n_balls:
            The number of balls `N`.

        in_first:
            The initial number of balls in the first box, or `None` for all of them.

    Returns:
        system:
            The model, with the count `urn.n` of balls in the first box.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import ehrenfest_urn
        >>> ehrenfest_urn(10).elements[0].n
        10
    """
    n = n_balls if in_first is None else in_first
    return System(elements=[EhrenfestUrnElement('urn', n=n, N=n_balls)])

graph_ising(graph=None, *, temperature=2.0, J=1.0, h=0.0, initial='random', sampler='metropolis', random_state=0)

Return the Ising model on a graph, ordered below a critical temperature set by the degrees of the nodes.

On a random 4-regular graph the Bethe approximation gives the critical temperature 2 / ln 2. One unit of time is one sweep of colored updates.

Parameters:

Name	Type	Description	Default
graph	Graph | None	The graph whose nodes hold the spins, or None for a random 4-regular graph of 400 nodes.	None
temperature	float	The temperature in units of the coupling, which can change during a run.	2.0
J	float	The coupling between neighbouring spins.	1.0
h	float	The external field.	0.0
initial	str	'random' for random spins or 'up' for all spins up.	'random'
sampler	str	The acceptance rule, 'metropolis', 'glauber', or 'heat-bath'.	'metropolis'
random_state	int | Generator | None	The seed or generator of the default graph and the random initial spins.	0

Returns:

Name	Type	Description
system	System	The model, with the spins spins.spin in node order.

Examples:

>>> import networkx as nx
>>> from skcomplex.models.physics.statistical_mechanics import graph_ising
>>> graph_ising(nx.cycle_graph(10)).elements[0].check().prepare().size_
10

Source code in src/skcomplex/models/physics/statistical_mechanics/_graph_spins.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('network', 'Monte Carlo', 'phase transition', 'spins'),
    references=('Argun et al. 2021', 'Dorogovtsev, Goltsev, and Mendes 2008'),
    defaults={'until': 200, 'dt': 1},
)
def graph_ising(
    graph: nx.Graph | None = None,
    *,
    temperature: float = 2.0,
    J: float = 1.0,
    h: float = 0.0,
    initial: str = 'random',
    sampler: str = 'metropolis',
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Ising model on a graph, ordered below a critical temperature set by the degrees of the nodes.

    On a random 4-regular graph the Bethe approximation gives the critical temperature `2 / ln 2`. One unit of time is
    one sweep of colored updates.

    Args:
        graph:
            The graph whose nodes hold the spins, or `None` for a random 4-regular graph of 400 nodes.

        temperature:
            The temperature in units of the coupling, which can change during a run.

        J:
            The coupling between neighbouring spins.

        h:
            The external field.

        initial:
            `'random'` for random spins or `'up'` for all spins up.

        sampler:
            The acceptance rule, `'metropolis'`, `'glauber'`, or `'heat-bath'`.

        random_state:
            The seed or generator of the default graph and the random initial spins.

    Returns:
        system:
            The model, with the spins `spins.spin` in node order.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.models.physics.statistical_mechanics import graph_ising
        >>> graph_ising(nx.cycle_graph(10)).elements[0].check().prepare().size_
        10
    """
    rng = np.random.default_rng(random_state)
    if graph is None:
        graph = nx.random_regular_graph(DEGREE, NODES, seed=int(rng.integers(2**31)))
    size = graph.number_of_nodes()
    spins = np.ones(size) if initial == 'up' else rng.choice([-1.0, 1.0], size=size)
    space = NetworkSpace(graph)
    group = LatticeElement(
        'spins',
        element=IsingSpinElement,
        space=space,
        values={'spin': spins, 'J': J, 'h': h, 'temperature': temperature},
    )
    sampling = IsingSamplingInteraction(
        'sampling',
        target='spins',
        rule=sampler,
        schedule='colored',
        neighborhood=AdjacencyNeighborhood(),
    )
    return System(elements=[group], interactions=[sampling])

ideal_gas(n_molecules=500, *, kT=1.0, size=1.0, n_dim=2, mass=1.0, random_state=0, label='gas')

Return an ideal gas of non-interacting molecules in a reflecting box.

Positions are uniform in the box and velocity components are normal with variance kT / mass, with the mean velocity removed. The molecules are drawn as points colored by speed.

Parameters:

Name	Type	Description	Default
n_molecules	int	The number of molecules.	500
kT	float	The temperature in energy units.	1.0
size	float	The side of the box.	1.0
n_dim	int	The number of dimensions.	2
mass	float	The mass of every molecule.	1.0
random_state	int | Generator | None	The seed or generator of the initial positions and velocities.	0
label	str	The label of the molecule group. The box is labelled box.	'gas'

Returns:

Name	Type	Description
system	System	The gas.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ideal_gas
>>> group = ideal_gas(1000, random_state=0).elements[0]
>>> group.label, group.element.__name__, group.size
('gas', 'MassPointElement', 1000)

Source code in src/skcomplex/models/physics/statistical_mechanics/_gases.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('gas', 'ode', 'boundary'),
    references=('Maxwell 1860',),
    defaults={'until': 2.0, 'dt': 0.01},
)
def ideal_gas(
    n_molecules: int = 500,
    *,
    kT: float = 1.0,
    size: float = 1.0,
    n_dim: int = 2,
    mass: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
    label: str = 'gas',
) -> System:
    """Return an ideal gas of non-interacting molecules in a reflecting box.

    Positions are uniform in the box and velocity components are normal with variance `kT / mass`, with the mean
    velocity removed. The molecules are drawn as points colored by speed.

    Args:
        n_molecules:
            The number of molecules.

        kT:
            The temperature in energy units.

        size:
            The side of the box.

        n_dim:
            The number of dimensions.

        mass:
            The mass of every molecule.

        random_state:
            The seed or generator of the initial positions and velocities.

        label:
            The label of the molecule group. The box is labelled `box`.

    Returns:
        system:
            The gas.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import ideal_gas
        >>> group = ideal_gas(1000, random_state=0).elements[0]
        >>> group.label, group.element.__name__, group.size
        ('gas', 'MassPointElement', 1000)
    """
    rng = np.random.default_rng(random_state)
    velocities = rng.normal(0.0, np.sqrt(kT / mass), size=(n_molecules, n_dim))
    velocities -= velocities.mean(axis=0)
    visual = VisualSpec('points', {'position': 'x', 'radius': size / 200, 'color_by': 'x_d1', 'states': ()})
    molecules = GroupElement(
        label,
        element=MassPointElement,
        size=n_molecules,
        values={'x': rng.uniform(0.0, size, size=(n_molecules, n_dim)), 'x_d1': velocities, 'm': mass},
        visual=visual,
    )
    return System(
        elements=[molecules],
        space=ContinuousSpace(extent=[size] * n_dim, boundary='reflecting'),
    )

ising(size=64, *, temperature=2.0, J=1.0, h=0.0, initial='random', sampler='metropolis', random_state=0)

Return the two-dimensional Ising model on a torus, ordered below the critical temperature 2 / ln(1 + sqrt 2).

One unit of time is one sweep of checkerboard updates.

Parameters:

Name	Type	Description	Default
size	int | Sequence[int]	The side of the square lattice, or its shape.	64
temperature	float	The temperature in units of the coupling, which can change during a run.	2.0
J	float	The coupling between neighbouring spins.	1.0
h	float	The external field.	0.0
initial	str	'random' for random spins or 'up' for all spins up.	'random'
sampler	str	The acceptance rule, 'metropolis', 'glauber', or 'heat-bath'.	'metropolis'
random_state	int | Generator | None	The seed or generator of the random initial spins.	0

Returns:

Name	Type	Description
system	System	The model, with the spins spins.spin.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import ising
>>> ising(8, temperature=1.5).elements[0].values['temperature']
1.5

Source code in src/skcomplex/models/physics/statistical_mechanics/_spins.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('lattice', 'Monte Carlo', 'phase transition', 'spins'),
    references=('Ising 1925', 'Onsager 1944', 'Jensen 2023'),
    defaults={'until': 200, 'dt': 1},
)
def ising(
    size: int | Sequence[int] = 64,
    *,
    temperature: float = 2.0,
    J: float = 1.0,
    h: float = 0.0,
    initial: str = 'random',
    sampler: str = 'metropolis',
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the two-dimensional Ising model on a torus, ordered below the critical temperature `2 / ln(1 + sqrt 2)`.

    One unit of time is one sweep of checkerboard updates.

    Args:
        size:
            The side of the square lattice, or its shape.

        temperature:
            The temperature in units of the coupling, which can change during a run.

        J:
            The coupling between neighbouring spins.

        h:
            The external field.

        initial:
            `'random'` for random spins or `'up'` for all spins up.

        sampler:
            The acceptance rule, `'metropolis'`, `'glauber'`, or `'heat-bath'`.

        random_state:
            The seed or generator of the random initial spins.

    Returns:
        system:
            The model, with the spins `spins.spin`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import ising
        >>> ising(8, temperature=1.5).elements[0].values['temperature']
        1.5
    """
    shape = _square(size)
    spins = np.ones(shape) if initial == 'up' else np.random.default_rng(random_state).choice([-1.0, 1.0], size=shape)
    space = LatticeSpace(shape)
    group = LatticeElement(
        'spins',
        element=IsingSpinElement,
        space=space,
        values={'spin': spins, 'J': J, 'h': h, 'temperature': temperature},
    )
    return System(elements=[group], interactions=[IsingSamplingInteraction('sampling', target='spins', rule=sampler)])

landau_relaxation(phi=0.1, *, T=0.5, Tc=1.0, a=1.0, b=1.0, h=0.0, Gamma=1.0)

Return the relaxation of a Landau order parameter, to sqrt(a (Tc - T) / b) below Tc and to 0 above it.

The relaxation time 1 / (Gamma a |T - Tc|) diverges at the critical temperature, which is critical slowing down.

Parameters:

Name	Type	Description	Default
phi	float	The initial order parameter.	0.1
T	float	The temperature.	0.5
Tc	float	The critical temperature.	1.0
a	float	The coefficient of the quadratic term per unit temperature.	1.0
b	float	The coefficient of the quartic term.	1.0
h	float	The external field.	0.0
Gamma	float	The kinetic coefficient.	1.0

Returns:

Name	Type	Description
system	System	The model, with the order parameter magnet.phi.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import landau_relaxation
>>> landau_relaxation(T=2.0).elements[0].T
2.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_master_equations.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('ode', 'mean field', 'phase transition', 'critical slowing down'),
    references=('Landau 1937', 'Gros 2015'),
    defaults={'until': 20.0, 'dt': 0.1},
)
def landau_relaxation(
    phi: float = 0.1,
    *,
    T: float = 0.5,
    Tc: float = 1.0,
    a: float = 1.0,
    b: float = 1.0,
    h: float = 0.0,
    Gamma: float = 1.0,
) -> System:
    """Return the relaxation of a Landau order parameter, to `sqrt(a (Tc - T) / b)` below `Tc` and to 0 above it.

    The relaxation time `1 / (Gamma a |T - Tc|)` diverges at the critical temperature, which is critical slowing
    down.

    Args:
        phi:
            The initial order parameter.

        T:
            The temperature.

        Tc:
            The critical temperature.

        a:
            The coefficient of the quadratic term per unit temperature.

        b:
            The coefficient of the quartic term.

        h:
            The external field.

        Gamma:
            The kinetic coefficient.

    Returns:
        system:
            The model, with the order parameter `magnet.phi`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import landau_relaxation
        >>> landau_relaxation(T=2.0).elements[0].T
        2.0
    """
    return System(elements=[OrderParameterElement('magnet', phi=phi, a=a, b=b, T=T, Tc=Tc, h=h, Gamma=Gamma)])

lennard_jones_fluid(n_particles=400, *, density=0.5, temperature=None, cutoff=2.5, random_state=0)

Return Lennard-Jones particles in a periodic square box, optionally held at a temperature by a thermostat.

Parameters:

Name	Type	Description	Default
n_particles	int	The number of particles.	400
density	float	The number of particles per unit area.	0.5
temperature	float | None	The temperature of the Langevin thermostat, or None for constant energy. A thermostat makes the model stochastic, so it runs with the Euler-Maruyama solver instead of the default velocity Verlet.	None
cutoff	float	The distance beyond which pairs do not interact.	2.5
random_state	int | Generator | None	The seed or generator of the initial jitter and velocities.	0

Returns:

Name	Type	Description
system	System	The model, with the particle group fluid in two dimensions.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import lennard_jones_fluid
>>> lennard_jones_fluid(100).elements[0].check().prepare().size_
100

Source code in src/skcomplex/models/physics/statistical_mechanics/_gases.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('molecular dynamics', 'pair potential', 'ode'),
    references=('Lennard-Jones 1924', 'Verlet 1967'),
    defaults={'until': 0.5, 'dt': 0.005, 'solver': VelocityVerletSolver()},
)
def lennard_jones_fluid(
    n_particles: int = 400,
    *,
    density: float = 0.5,
    temperature: float | None = None,
    cutoff: float = 2.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return Lennard-Jones particles in a periodic square box, optionally held at a temperature by a thermostat.

    Args:
        n_particles:
            The number of particles.

        density:
            The number of particles per unit area.

        temperature:
            The temperature of the Langevin thermostat, or `None` for constant energy. A thermostat makes the model
            stochastic, so it runs with the Euler-Maruyama solver instead of the default velocity Verlet.

        cutoff:
            The distance beyond which pairs do not interact.

        random_state:
            The seed or generator of the initial jitter and velocities.

    Returns:
        system:
            The model, with the particle group `fluid` in two dimensions.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import lennard_jones_fluid
        >>> lennard_jones_fluid(100).elements[0].check().prepare().size_
        100
    """
    rng = np.random.default_rng(random_state)
    size = float(np.sqrt(n_particles / density))
    side = int(np.ceil(np.sqrt(n_particles)))
    spacing = size / side
    grid = np.stack(np.meshgrid(np.arange(side), np.arange(side)), axis=-1).reshape(-1, 2)[:n_particles]
    positions = np.mod((grid + 0.5) * spacing + rng.normal(0, 0.02 * spacing, (n_particles, 2)), size)
    velocities = rng.normal(0, np.sqrt(temperature or 1.0), (n_particles, 2))
    velocities -= velocities.mean(axis=0)
    interactions: list[BaseInteraction] = [
        LennardJonesInteraction('lj', first='fluid', second='fluid', cutoff=cutoff, divisor='m', period=size),
    ]
    if temperature is not None:
        interactions.append(
            LangevinThermostatInteraction('thermostat', target='fluid', temperature=temperature),
        )
    return System(
        elements=[
            GroupElement(
                'fluid',
                element=MassPointElement,
                size=n_particles,
                values={'x': positions, 'x_d1': velocities},
            ),
        ],
        interactions=interactions,
        space=ContinuousSpace(extent=(size, size), boundary='periodic'),
    )

local_magnetization(member)

Return the z component of one spin.

Parameters:

Name	Type	Description	Default
member	SimpleNamespace	The namespace of one spin.	required

Returns:

Name	Type	Description
operator	Expression	The local magnetisation.

Source code in src/skcomplex/models/physics/statistical_mechanics/_quantum_spins.py

def local_magnetization(member: SimpleNamespace) -> Expression:
    """Return the `z` component of one spin.

    Args:
        member:
            The namespace of one spin.

    Returns:
        operator:
            The local magnetisation.
    """
    return local_operator(member.psi, 'z')

mattis_spin_glass(size=48, *, temperature=1.5, J=1.0, h=0.0, random_state=0)

Return the Mattis spin glass on a periodic square lattice, with couplings J_ij = J xi_i xi_j of a pattern xi.

The gauge change sigma_i = xi_i s_i maps the model at h = 0 to the Ising model, so the overlap m_xi = sum_i xi_i s_i / N orders below the Ising critical temperature 2 J / ln(1 + sqrt 2) while the magnetization averages to zero over the disorder. One unit of time is one sweep of Metropolis updates.

Parameters:

Name	Type	Description	Default
size	int | Sequence[int]	The side of the square lattice, or its two sides.	48
temperature	float	The temperature, which can change during a run.	1.5
J	float	The coupling.	1.0
h	float	The external field.	0.0
random_state	int | Generator | None	The seed or generator of the pattern and of the random initial spins.	0

Returns:

Name	Type	Description
system	System	The model, with the spins spins.spin, the pattern spins.pattern, and the observables magnetization and overlap.

Examples:

>>> import numpy as np
>>> from skcomplex.models.physics.statistical_mechanics import mattis_spin_glass
>>> values = mattis_spin_glass(8).elements[0].values
>>> pattern = values['pattern'].reshape(8, 8)
>>> bool(np.array_equal(values['J_right'].reshape(8, 8), pattern * np.roll(pattern, -1, axis=1)))
True

Source code in src/skcomplex/models/physics/statistical_mechanics/_spin_glasses.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('spin glass', 'disorder', 'lattice', 'Monte Carlo', 'spins'),
    references=('Mattis 1976', 'Bertin 2021'),
    defaults={'until': 300, 'dt': 2},
    observables=_pattern_observables,
)
def mattis_spin_glass(
    size: int | Sequence[int] = 48,
    *,
    temperature: float = 1.5,
    J: float = 1.0,
    h: float = 0.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Mattis spin glass on a periodic square lattice, with couplings `J_ij = J xi_i xi_j` of a pattern `xi`.

    The gauge change `sigma_i = xi_i s_i` maps the model at `h = 0` to the Ising model, so the overlap
    `m_xi = sum_i xi_i s_i / N` orders below the Ising critical temperature `2 J / ln(1 + sqrt 2)` while the
    magnetization averages to zero over the disorder. One unit of time is one sweep of Metropolis updates.

    Args:
        size:
            The side of the square lattice, or its two sides.

        temperature:
            The temperature, which can change during a run.

        J:
            The coupling.

        h:
            The external field.

        random_state:
            The seed or generator of the pattern and of the random initial spins.

    Returns:
        system:
            The model, with the spins `spins.spin`, the pattern `spins.pattern`, and the observables `magnetization`
            and `overlap`.

    Examples:
        >>> import numpy as np
        >>> from skcomplex.models.physics.statistical_mechanics import mattis_spin_glass
        >>> values = mattis_spin_glass(8).elements[0].values
        >>> pattern = values['pattern'].reshape(8, 8)
        >>> bool(np.array_equal(values['J_right'].reshape(8, 8), pattern * np.roll(pattern, -1, axis=1)))
        True
    """
    shape = _square(size)
    rng = np.random.default_rng(random_state)
    pattern = rng.choice([-1.0, 1.0], size=shape)
    right = J * pattern * np.roll(pattern, -1, axis=1)
    down = J * pattern * np.roll(pattern, -1, axis=0)
    spins = rng.choice([-1.0, 1.0], size=shape)
    group = LatticeElement(
        'spins',
        element=MattisSpinElement,
        space=LatticeSpace(shape),
        values={'spin': spins, 'pattern': pattern, 'h': h, 'temperature': temperature, **_bond_states(right, down)},
        visual=VisualSpec('image', {'states': ()}),
    )
    return System(elements=[group], interactions=[BondSamplingInteraction('sampling', target='spins')])

potts(size=64, *, q=3, temperature=1.0, J=1.0, sampler='heat-bath', random_state=0)

Return the q-state Potts model on a torus, ordered below the critical temperature 1 / ln(1 + sqrt q).

Parameters:

Name	Type	Description	Default
size	int | Sequence[int]	The side of the square lattice, or its shape.	64
q	int	The number of states.	3
temperature	float	The temperature in units of the coupling.	1.0
J	float	The coupling between neighbouring spins.	1.0
sampler	str	The acceptance rule, 'metropolis', 'glauber', or 'heat-bath'.	'heat-bath'
random_state	int | Generator | None	The seed or generator of the random initial states.	0

Returns:

Name	Type	Description
system	System	The model, with the states spins.color.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import potts
>>> potts(8, q=4).interactions[0].q
4

Source code in src/skcomplex/models/physics/statistical_mechanics/_spins.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('lattice', 'Monte Carlo', 'phase transition', 'spins'),
    references=('Potts 1952', 'Wu 1982'),
    defaults={'until': 200, 'dt': 1},
)
def potts(
    size: int | Sequence[int] = 64,
    *,
    q: int = 3,
    temperature: float = 1.0,
    J: float = 1.0,
    sampler: str = 'heat-bath',
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the q-state Potts model on a torus, ordered below the critical temperature `1 / ln(1 + sqrt q)`.

    Args:
        size:
            The side of the square lattice, or its shape.

        q:
            The number of states.

        temperature:
            The temperature in units of the coupling.

        J:
            The coupling between neighbouring spins.

        sampler:
            The acceptance rule, `'metropolis'`, `'glauber'`, or `'heat-bath'`.

        random_state:
            The seed or generator of the random initial states.

    Returns:
        system:
            The model, with the states `spins.color`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import potts
        >>> potts(8, q=4).interactions[0].q
        4
    """
    shape = _square(size)
    colors = np.random.default_rng(random_state).integers(0, q, size=shape).astype(float)
    space = LatticeSpace(shape)
    group = LatticeElement(
        'spins',
        element=PottsSpinElement,
        space=space,
        values={'color': colors, 'J': J, 'temperature': temperature},
    )
    return System(
        elements=[group],
        interactions=[PottsSamplingInteraction('sampling', target='spins', q=q, rule=sampler)],
    )

q2r(shape=(64, 64), *, up=0.3, random_state=0)

Return the Q2R automaton, a reversible, energy-conserving Ising dynamics on a periodic lattice of even sides.

The spins start uncorrelated with a fraction up of spins up, and the sites of even coordinate sum update first.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of sites along each axis, all even.	(64, 64)
up	float	The probability that an initial spin is up.	0.3
random_state	int | Generator | None	The seed or generator of the initial spins.	0

Returns:

Name	Type	Description
system	System	The model, with the spins spins.spin and phases spins.phase.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import q2r
>>> q2r((8,)).elements[0].check().prepare().size_
8

Source code in src/skcomplex/models/physics/statistical_mechanics/_processes.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('lattice', 'cellular automaton', 'reversible', 'spins'),
    references=('Vichniac 1984', 'Pomeau 1984', 'Lindgren 2024'),
    defaults={'until': 100, 'dt': 1},
)
def q2r(
    shape: Sequence[int] = (64, 64),
    *,
    up: float = 0.3,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Q2R automaton, a reversible, energy-conserving Ising dynamics on a periodic lattice of even sides.

    The spins start uncorrelated with a fraction `up` of spins up, and the sites of even coordinate sum update first.

    Args:
        shape:
            The number of sites along each axis, all even.

        up:
            The probability that an initial spin is up.

        random_state:
            The seed or generator of the initial spins.

    Returns:
        system:
            The model, with the spins `spins.spin` and phases `spins.phase`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import q2r
        >>> q2r((8,)).elements[0].check().prepare().size_
        8
    """
    shape = tuple(shape)
    spins = np.where(np.random.default_rng(random_state).random(shape) < up, 1.0, -1.0)
    space = LatticeSpace(shape)
    group = LatticeElement(
        'spins',
        element=Q2RSpinElement,
        space=space,
        values={'spin': spins, 'phase': _checkerboard(shape)},
    )
    flip = LocalRuleInteraction(
        'flip',
        target='spins',
        state=('spin', 'phase'),
        rule=_flip_without_cost,
        neighborhood=VonNeumannNeighborhood(n_dim=len(shape), center=True),
    )
    return System(elements=[group], interactions=[flip])

quantum_heisenberg(n_sites=12, *, J=1.0, delta=1.0, h=0.0, gamma=0.0, boundary='periodic', initial='flip')

Return a chain of spins one half with the XXZ Hamiltonian, a field, and optional dephasing.

The Hamiltonian is J Σ (S^x_i S^x_j + S^y_i S^y_j + Δ S^z_i S^z_j) - h Σ S^z_i over neighbouring pairs, with ħ = 1. A dephasing rate gamma > 0 adds the dissipator S^z_i on every site, so the state is a density matrix that loses its coherence. One flipped spin in a ferromagnetic background spreads as a magnon, whose probability at distance d is J_d(J t)² in the XX chain, a light cone in the space-time image of S^z. The total magnetisation is conserved.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of spins.	12
J	float	The exchange coupling, antiferromagnetic when positive.	1.0
delta	float	The anisotropy Δ.	1.0
h	float	The field.	0.0
gamma	float	The dephasing rate.	0.0
boundary	str	'periodic' for a ring, or 'open' for a chain with two ends.	'periodic'
initial	str	'flip' for one flipped spin at site n_sites // 2, 'neel' for alternating spins, 'domain_wall' for the left half up and the right half down, or 'ground' for the ground state, which run gives the simulation as initial_state={'quantum_heisenberg.state': 'ground'}.	'flip'

Returns:

Name	Type	Description
system	System	The model, with the observables sz of every site, magnetization, and the energy.

Raises:

Type	Description
ValidationError	The boundary or the initial state is unknown.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import quantum_heisenberg
>>> quantum_heisenberg(6).elements[0].check().prepare().size_
6

Source code in src/skcomplex/models/physics/statistical_mechanics/_quantum_spins.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('quantum', 'spins', 'lattice', 'magnetism', 'open system'),
    references=('Heisenberg 1928', 'Mihailović et al. 2023'),
    defaults={'until': 8.0, 'dt': 0.05, 'compiler': OperatorCompiler()},
    initial_state=lambda initial, **_: {'quantum_heisenberg.state': 'ground'} if initial == 'ground' else None,
    observables=lambda system: [
        ObservableSpec('sz', local_magnetization, group='chain', reduce=_sites),
        ObservableSpec('magnetization', total_magnetization),
    ],
)
def quantum_heisenberg(
    n_sites: int = 12,
    *,
    J: float = 1.0,
    delta: float = 1.0,
    h: float = 0.0,
    gamma: float = 0.0,
    boundary: str = 'periodic',
    initial: str = 'flip',
) -> System:
    """Return a chain of spins one half with the XXZ Hamiltonian, a field, and optional dephasing.

    The Hamiltonian is `J Σ (S^x_i S^x_j + S^y_i S^y_j + Δ S^z_i S^z_j) - h Σ S^z_i` over neighbouring pairs, with
    `ħ = 1`. A dephasing rate `gamma > 0` adds the dissipator `S^z_i` on every site, so the state is a density matrix
    that loses its coherence. One flipped spin in a ferromagnetic background spreads as a magnon, whose probability at
    distance `d` is `J_d(J t)²` in the XX chain, a light cone in the space-time image of `S^z`. The total magnetisation
    is conserved.

    Args:
        n_sites:
            The number of spins.

        J:
            The exchange coupling, antiferromagnetic when positive.

        delta:
            The anisotropy `Δ`.

        h:
            The field.

        gamma:
            The dephasing rate.

        boundary:
            `'periodic'` for a ring, or `'open'` for a chain with two ends.

        initial:
            `'flip'` for one flipped spin at site `n_sites // 2`, `'neel'` for alternating spins, `'domain_wall'` for
            the left half up and the right half down, or `'ground'` for the ground state, which `run` gives the
            simulation as `initial_state={'quantum_heisenberg.state': 'ground'}`.

    Returns:
        system:
            The model, with the observables `sz` of every site, `magnetization`, and the energy.

    Raises:
        ValidationError:
            The boundary or the initial state is unknown.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import quantum_heisenberg
        >>> quantum_heisenberg(6).elements[0].check().prepare().size_
        6
    """
    if boundary not in BOUNDARIES or initial not in INITIAL_STATES:
        error_msg = (
            f'Unknown boundary {boundary!r} or initial state {initial!r}. Use a boundary of {", ".join(BOUNDARIES)} '
            f'and an initial state of {", ".join(INITIAL_STATES)}.'
        )
        raise ValidationError(error_msg)
    space = LatticeSpace(n_sites, boundary=BOUNDARIES[boundary])
    labels = _initial_labels(n_sites, initial)
    chain = LatticeElement('chain', element=SpinElement, space=space, values={'psi': labels})
    interactions: list[BaseInteraction] = [
        XXZExchangeInteraction('exchange', target='chain', J=J, delta=delta),
        ZeemanFieldInteraction('field', target='chain', h=h),
    ]
    if gamma > 0:
        interactions.append(DephasingInteraction('dephasing', target='chain', gamma=gamma))
    return System('quantum_heisenberg', elements=[chain], interactions=interactions)

random_energy_model(n_spins=10, *, temperature=1.0, random_state=0)

Return the random energy model, 2^N quenched Gaussian energies of variance N / 2 at a temperature.

The entropy per spin follows ln 2 - 1 / (4 T^2) above the glass temperature T_g = 1 / (2 sqrt(ln 2)) and freezes below it, where the Boltzmann weight falls on the few lowest levels (Bertin 2021 §1.5.3). The energies are drawn from random_state, so a disorder average sweeps it through GridDesign.

Parameters:

Name	Type	Description	Default
n_spins	int	The number of spins N.	10
temperature	float	The temperature.	1.0
random_state	int | Generator | None	The seed or generator of the energies.	0

Returns:

Name	Type	Description
system	System	The model, with the states ensemble.log_partition, ensemble.energy, and ensemble.entropy.

Examples:

>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.physics.statistical_mechanics import random_energy_model
>>> system = random_energy_model(12, temperature=2.0)
>>> trajectory = Simulation(system).check().prepare().run(1.0).trajectory_
>>> bool(abs(trajectory.values('ensemble.entropy')[0] - (np.log(2) - 1 / 16)) < 0.05)
True

Source code in src/skcomplex/models/physics/statistical_mechanics/_random_energy.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('disorder', 'glass transition', 'static'),
    references=('Derrida 1981', 'Bertin 2021'),
)
def random_energy_model(
    n_spins: int = 10,
    *,
    temperature: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the random energy model, `2^N` quenched Gaussian energies of variance `N / 2` at a temperature.

    The entropy per spin follows `ln 2 - 1 / (4 T^2)` above the glass temperature `T_g = 1 / (2 sqrt(ln 2))` and
    freezes below it, where the Boltzmann weight falls on the few lowest levels (Bertin 2021 §1.5.3). The energies are
    drawn from `random_state`, so a disorder average sweeps it through `GridDesign`.

    Args:
        n_spins:
            The number of spins `N`.

        temperature:
            The temperature.

        random_state:
            The seed or generator of the energies.

    Returns:
        system:
            The model, with the states `ensemble.log_partition`, `ensemble.energy`, and `ensemble.entropy`.

    Examples:
        >>> from skcomplex.engine.runtime import Simulation
        >>> from skcomplex.models.physics.statistical_mechanics import random_energy_model
        >>> system = random_energy_model(12, temperature=2.0)
        >>> trajectory = Simulation(system).check().prepare().run(1.0).trajectory_
        >>> bool(abs(trajectory.values('ensemble.entropy')[0] - (np.log(2) - 1 / 16)) < 0.05)
        True
    """
    size = 2**n_spins
    energies = np.random.default_rng(random_state).normal(0.0, np.sqrt(n_spins / 2), size)
    visual = VisualSpec('histograms', {'color_by': 'energy'})
    levels = GroupElement('levels', element=EnergyLevelElement, size=size, values={'energy': energies}, visual=visual)
    sums = BoltzmannSumInteraction(
        'sums',
        target='ensemble',
        levels='levels',
        temperature=temperature,
        n_spins=float(n_spins),
    )
    return System(elements=[levels, BoltzmannAverageElement('ensemble')], interactions=[sums], time=StaticTime())

sherrington_kirkpatrick(n_spins=200, *, temperature=0.5, J=1.0, h=0.0, random_state=0)

Return the Sherrington-Kirkpatrick spin glass with Glauber dynamics, frozen below the temperature T_c = J.

The energy is -sum_{i<j} J_ij s_i s_j - h sum_i s_i with quenched Gaussian couplings J_ij of mean 0 and variance J^2 / N, drawn once. Above T_c the energy per spin is -J^2 / (2 T), and below it the spins freeze in random directions, with no magnetization. The spins are drawn in a row, which orders the drawing only.

Parameters:

Name	Type	Description	Default
n_spins	int	The number of spins.	200
temperature	float	The temperature, which can change during a run.	0.5
J	float	The standard deviation of the couplings times the square root of the number of spins.	1.0
h	float	The external field.	0.0
random_state	int | Generator | None	The seed or generator of the couplings and of the random initial spins.	0

Returns:

Name	Type	Description
system	System	The model, with the spins spins.spin, the coupling matrix of unit scale in the values of the group, and the observable magnetization.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import sherrington_kirkpatrick
>>> couplings = sherrington_kirkpatrick(10).elements[0].values['couplings']
>>> bool((couplings == couplings.T).all()), float(couplings.trace())
(True, 0.0)

Source code in src/skcomplex/models/physics/statistical_mechanics/_fully_connected_spins.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('spin glass', 'disorder', 'mean field', 'spins', 'jump process'),
    references=('Sherrington and Kirkpatrick 1975', 'Bertin 2021', 'Mezard, Parisi, and Virasoro 1987'),
    defaults={'until': 50.0, 'dt': 0.5},
    observables=lambda system: [_magnetization()],
)
def sherrington_kirkpatrick(
    n_spins: int = 200,
    *,
    temperature: float = 0.5,
    J: float = 1.0,
    h: float = 0.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Sherrington-Kirkpatrick spin glass with Glauber dynamics, frozen below the temperature `T_c = J`.

    The energy is `-sum_{i<j} J_ij s_i s_j - h sum_i s_i` with quenched Gaussian couplings `J_ij` of mean 0 and
    variance `J^2 / N`, drawn once. Above `T_c` the energy per spin is `-J^2 / (2 T)`, and below it the spins freeze
    in random directions, with no magnetization. The spins are drawn in a row, which orders the drawing only.

    Args:
        n_spins:
            The number of spins.

        temperature:
            The temperature, which can change during a run.

        J:
            The standard deviation of the couplings times the square root of the number of spins.

        h:
            The external field.

        random_state:
            The seed or generator of the couplings and of the random initial spins.

    Returns:
        system:
            The model, with the spins `spins.spin`, the coupling matrix of unit scale in the values of the group, and
            the observable `magnetization`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import sherrington_kirkpatrick
        >>> couplings = sherrington_kirkpatrick(10).elements[0].values['couplings']
        >>> bool((couplings == couplings.T).all()), float(couplings.trace())
        (True, 0.0)
    """
    rng = np.random.default_rng(random_state)
    upper = np.triu(rng.normal(0.0, 1.0 / np.sqrt(n_spins), (n_spins, n_spins)), 1)
    spins = _initial_spins(n_spins, 'random', rng)
    group = LatticeElement(
        'spins',
        element=FullyConnectedSpinElement,
        space=LatticeSpace(n_spins),
        values={'spin': spins, 'J': J, 'h': h, 'temperature': temperature, 'couplings': upper + upper.T},
        visual=SPIN_IMAGE,
    )
    return System(elements=[group])

site_percolation(shape=(64, 64), *, p=0.6, random_state=0)

Return site percolation probed by a fire lit along the top row, which spans the lattice above p_c = 0.5927.

Sites are occupied with probability p, and the fire burns the clusters connected to the top row, one layer of sites per step.

Parameters:

Name	Type	Description	Default
shape	Sequence[int]	The number of rows and columns.	(64, 64)
p	float	The occupation probability.	0.6
random_state	int | Generator | None	The seed or generator of the occupations.	0

Returns:

Name	Type	Description
system	System	The model, with the covers lattice.cover.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import site_percolation
>>> site_percolation((10, 10)).elements[0].check().prepare().size_
100

Source code in src/skcomplex/models/physics/statistical_mechanics/_processes.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('lattice', 'percolation', 'phase transition', 'fire'),
    references=('Broadbent and Hammersley 1957', 'Stauffer and Aharony 1994', 'Charbonneau 2017'),
    defaults={'until': 200, 'dt': 1},
)
def site_percolation(
    shape: Sequence[int] = (64, 64),
    *,
    p: float = 0.6,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return site percolation probed by a fire lit along the top row, which spans the lattice above `p_c = 0.5927`.

    Sites are occupied with probability `p`, and the fire burns the clusters connected to the top row, one layer of
    sites per step.

    Args:
        shape:
            The number of rows and columns.

        p:
            The occupation probability.

        random_state:
            The seed or generator of the occupations.

    Returns:
        system:
            The model, with the covers `lattice.cover`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import site_percolation
        >>> site_percolation((10, 10)).elements[0].check().prepare().size_
        100
    """
    shape = tuple(shape)
    cover = (np.random.default_rng(random_state).random(shape) < p).astype(float)
    cover[0] *= BURNING
    space = LatticeSpace(shape, boundary='fixed')
    return System(
        elements=[LatticeElement('lattice', element=PercolationSiteElement, space=space, values={'cover': cover})],
    )

tasep(n_sites=200, *, p=0.5, alpha=0.3, beta=0.7, ring=False, density=0.5, random_state=0)

Return the totally asymmetric simple exclusion process with parallel update, on an open line or a ring.

Every step all particles whose right neighbour is empty hop to it with probability p at once, a transfer along every edge of the lane with one draw per edge. On an open line particles enter the first site with probability alpha and leave the last site with probability beta. On a ring at density rho the current is (1 - sqrt(1 - 4 p rho (1 - rho))) / 2. On an open line the entry and exit probabilities select the low-density, high-density, or maximal-current phase.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of sites.	200
p	float	The probability of a hop.	0.5
alpha	float	The probability of an entry at the left end of an open line.	0.3
beta	float	The probability of an exit at the right end of an open line.	0.7
ring	bool	Whether the line is a ring instead of an open line.	False
density	float	The initial fraction of occupied sites.	0.5
random_state	int | Generator | None	The seed or generator of the initial occupations.	0

Returns:

Name	Type	Description
system	System	The model, with the occupations lane.occupied and the hops lane.hops, and on an open line the entries entry.flow_0 and the exits exit.flow_0.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import tasep
>>> [interaction.label for interaction in tasep(20).interactions]
['hops', 'entry', 'exit']

Source code in src/skcomplex/models/physics/statistical_mechanics/_processes.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('lattice', 'stochastic', 'exclusion', 'traffic', 'nonequilibrium'),
    references=('Derrida, Domany, and Mukamel 1992', 'Schütz 2001', 'Helbing 2001'),
    defaults={'until': 500, 'dt': 1},
)
def tasep(
    n_sites: int = 200,
    *,
    p: float = 0.5,
    alpha: float = 0.3,
    beta: float = 0.7,
    ring: bool = False,
    density: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the totally asymmetric simple exclusion process with parallel update, on an open line or a ring.

    Every step all particles whose right neighbour is empty hop to it with probability `p` at once, a transfer along
    every edge of the lane with one draw per edge. On an open line particles enter the first site with probability
    `alpha` and leave the last site with probability `beta`. On a ring at density `rho` the current is
    `(1 - sqrt(1 - 4 p rho (1 - rho))) / 2`. On an open line the entry and exit probabilities select the low-density,
    high-density, or maximal-current phase.

    Args:
        n_sites:
            The number of sites.

        p:
            The probability of a hop.

        alpha:
            The probability of an entry at the left end of an open line.

        beta:
            The probability of an exit at the right end of an open line.

        ring:
            Whether the line is a ring instead of an open line.

        density:
            The initial fraction of occupied sites.

        random_state:
            The seed or generator of the initial occupations.

    Returns:
        system:
            The model, with the occupations `lane.occupied` and the hops `lane.hops`, and on an open line the
            entries `entry.flow_0` and the exits `exit.flow_0`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import tasep
        >>> [interaction.label for interaction in tasep(20).interactions]
        ['hops', 'entry', 'exit']
    """
    rng = np.random.default_rng(random_state)
    occupied = np.zeros(n_sites)
    occupied[rng.choice(n_sites, round(density * n_sites), replace=False)] = 1.0
    boundary = 'periodic' if ring else 'fixed'
    lane = LatticeElement(
        'lane',
        element=ExclusionSiteElement,
        space=LatticeSpace(n_sites, boundary=boundary),
        values={'occupied': occupied, 'p': p, 'alpha': alpha, 'beta': beta},
    )
    hops = TransferInteraction(
        'hops',
        source='lane',
        target='lane',
        state='occupied',
        amount=_hop,
        flow='hops',
        offset=1,
        boundary=boundary,
    )
    if ring:
        return System(elements=[lane], interactions=[hops], time=DiscreteTime())
    first, last = 'lane[0]', f'lane[{n_sites - 1}]'
    entry = TransferInteraction(
        'entry',
        source=first,
        target=first,
        state=Source(),
        into='occupied',
        amount=_enter,
    )
    leaving = TransferInteraction('exit', source=last, target=last, state='occupied', into=Sink(), amount=_leave)
    return System(elements=[lane], interactions=[hops, entry, leaving], time=DiscreteTime())

total_magnetization(s)

Return the sum of the z components of the spins of the chain.

Parameters:

Name	Type	Description	Default
s	SimpleNamespace	The namespace of the system.	required

Returns:

Name	Type	Description
operator	Expression	The total magnetisation.

Source code in src/skcomplex/models/physics/statistical_mechanics/_quantum_spins.py

def total_magnetization(s: SimpleNamespace) -> Expression:
    """Return the sum of the `z` components of the spins of the chain.

    Args:
        s:
            The namespace of the system.

    Returns:
        operator:
            The total magnetisation.
    """
    return sum(local_operator(member.psi, 'z') for member in s.chain)

trap_model(n_walkers=500, *, T=0.5, tau0=1.0, random_state=0)

Return the Bouchaud trap model, walkers that escape from traps of exponentially distributed depths.

The escape times tau0 exp(E / T) of depths with the density exp(-E) have the tail tau^(-1 - T). Below the glass temperature T = 1 their mean diverges, so the walkers fall into ever deeper traps, the number of escapes grows as t^T instead of linearly, and the system ages without reaching equilibrium.

Parameters:

Name	Type	Description	Default
n_walkers	int	The number of independent walkers.	500
T	float	The temperature, in units of the mean depth.	0.5
tau0	float	The escape time from a trap of zero depth.	1.0
random_state	int | Generator | None	The seed or generator of the first traps.	0

Returns:

Name	Type	Description
system	System	The model, with the places walkers.x, the depths walkers.E, and the escape count escapes.renewals.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import trap_model
>>> trap_model(10, T=0.8).elements[0].values['T']
0.8

Source code in src/skcomplex/models/physics/statistical_mechanics/_traps.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('trap model', 'renewal process', 'aging', 'glass'),
    references=('Bouchaud 1992', 'Monthus and Bouchaud 1996', 'Bertin 2021'),
    defaults={'until': 10000.0, 'dt': 50.0},
)
def trap_model(
    n_walkers: int = 500,
    *,
    T: float = 0.5,
    tau0: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Bouchaud trap model, walkers that escape from traps of exponentially distributed depths.

    The escape times `tau0 exp(E / T)` of depths with the density `exp(-E)` have the tail `tau^(-1 - T)`. Below the
    glass temperature `T = 1` their mean diverges, so the walkers fall into ever deeper traps, the number of escapes
    grows as `t^T` instead of linearly, and the system ages without reaching equilibrium.

    Args:
        n_walkers:
            The number of independent walkers.

        T:
            The temperature, in units of the mean depth.

        tau0:
            The escape time from a trap of zero depth.

        random_state:
            The seed or generator of the first traps.

    Returns:
        system:
            The model, with the places `walkers.x`, the depths `walkers.E`, and the escape count `escapes.renewals`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import trap_model
        >>> trap_model(10, T=0.8).elements[0].values['T']
        0.8
    """
    rng = np.random.default_rng(random_state)
    walkers = GroupElement(
        'walkers',
        element=TrapWalkerElement,
        size=n_walkers,
        values={'x': rng.random((n_walkers, 2)), 'E': rng.exponential(1.0, n_walkers), 'T': T, 'tau0': tau0},
    )
    return System(elements=[walkers], space=ContinuousSpace(n_dim=2))

two_phase_zone_exchange(n_zones=100, *, Na=100.0, p=0.1, k=1.0)

Return independent zones of a two-phase apparatus, each exchanging particles with the rest of the apparatus.

The zones start empty and their numbers of particles settle on the binomial distribution of mean Na p and variance Na p (1 - p).

Parameters:

Name	Type	Description	Default
n_zones	int	The number of independent zones.	100
Na	float	The number of particles in the apparatus and a zone.	100.0
p	float	The volume fraction of a zone.	0.1
k	float	The rate of exchange.	1.0

Returns:

Name	Type	Description
system	System	The model, with the numbers zones.N.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import two_phase_zone_exchange
>>> two_phase_zone_exchange(5).elements[0].check().prepare().size_
5

Source code in src/skcomplex/models/physics/statistical_mechanics/_master_equations.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('jump process', 'binomial', 'fluctuations', 'two-phase flow'),
    references=('Barsky 2020',),
    defaults={'until': 20.0, 'dt': 0.5},
)
def two_phase_zone_exchange(
    n_zones: int = 100,
    *,
    Na: float = 100.0,
    p: float = 0.1,
    k: float = 1.0,
) -> System:
    """Return independent zones of a two-phase apparatus, each exchanging particles with the rest of the apparatus.

    The zones start empty and their numbers of particles settle on the binomial distribution of mean `Na p` and
    variance `Na p (1 - p)`.

    Args:
        n_zones:
            The number of independent zones.

        Na:
            The number of particles in the apparatus and a zone.

        p:
            The volume fraction of a zone.

        k:
            The rate of exchange.

    Returns:
        system:
            The model, with the numbers `zones.N`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import two_phase_zone_exchange
        >>> two_phase_zone_exchange(5).elements[0].check().prepare().size_
        5
    """
    return System(
        elements=[GroupElement('zones', element=ZoneElement, size=n_zones, values={'Na': Na, 'p': p, 'k': k})],
    )

xy(size=64, *, temperature=0.5, J=1.0, width=1.0, random_state=0)

Return the XY model of planar spins on a torus, with vortices that bind below the Kosterlitz-Thouless point.

Parameters:

Name	Type	Description	Default
size	int | Sequence[int]	The side of the square lattice, or its shape.	64
temperature	float	The temperature in units of the coupling.	0.5
J	float	The coupling between neighbouring spins.	1.0
width	float	The largest change of a proposed rotation.	1.0
random_state	int | Generator | None	The seed or generator of the random initial angles.	0

Returns:

Name	Type	Description
system	System	The model, with the angles spins.theta.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import xy
>>> xy(8).interactions[0].width
1.0

Source code in src/skcomplex/models/physics/statistical_mechanics/_spins.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('lattice', 'Monte Carlo', 'vortices', 'spins'),
    references=('Kosterlitz and Thouless 1973', 'Jensen 2023'),
    defaults={'until': 200, 'dt': 1},
)
def xy(
    size: int | Sequence[int] = 64,
    *,
    temperature: float = 0.5,
    J: float = 1.0,
    width: float = 1.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the XY model of planar spins on a torus, with vortices that bind below the Kosterlitz-Thouless point.

    Args:
        size:
            The side of the square lattice, or its shape.

        temperature:
            The temperature in units of the coupling.

        J:
            The coupling between neighbouring spins.

        width:
            The largest change of a proposed rotation.

        random_state:
            The seed or generator of the random initial angles.

    Returns:
        system:
            The model, with the angles `spins.theta`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import xy
        >>> xy(8).interactions[0].width
        1.0
    """
    shape = _square(size)
    angles = np.random.default_rng(random_state).uniform(0, 2 * np.pi, size=shape)
    space = LatticeSpace(shape)
    group = LatticeElement(
        'spins',
        element=XYSpinElement,
        space=space,
        values={'theta': angles, 'J': J, 'temperature': temperature},
    )
    return System(elements=[group], interactions=[XYSamplingInteraction('sampling', target='spins', width=width)])

zero_range_process(n_sites=50, *, density=3.0, b=4.0, p=1.0, q=0.0, random_state=0)

Return the zero-range process on a ring, which condenses above the critical density 1 / (b - 2) for b > 2.

A particle leaves a site with n particles at rate u(n) = 1 + b / n, to the right with weight p and to the left with weight q. The stationary distribution factorizes with the weights f(n) = prod_k 1 / u(k), which decay as n^-b. Above the critical density a fluid at density 1 / (b - 2) coexists with one site that holds the excess particles.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of sites.	50
density	float	The number of particles per site.	3.0
b	float	The strength of the decrease of the departure rate with the occupation.	4.0
p	float	The rate scale of hops to the right.	1.0
q	float	The rate scale of hops to the left.	0.0
random_state	int | Generator | None	The seed or generator of the initial sites of the particles.	0

Returns:

Name	Type	Description
system	System	The model, with the occupations sites.n.

Examples:

>>> from skcomplex.models.physics.statistical_mechanics import zero_range_process
>>> system = zero_range_process(10, density=2.0)
>>> float(system.elements[0].values['n'].sum()), len(system.interactions)
(20.0, 1)

Source code in src/skcomplex/models/physics/statistical_mechanics/_driven_particles.py

@register_model(
    field='physics/statistical_mechanics',
    tags=('zero-range process', 'condensation', 'nonequilibrium', 'jump process'),
    references=('Evans and Hanney 2005', 'Bertin 2021'),
    defaults={'until': 2000.0, 'dt': 10.0, 'backend': NumbaBackend()},
)
def zero_range_process(
    n_sites: int = 50,
    *,
    density: float = 3.0,
    b: float = 4.0,
    p: float = 1.0,
    q: float = 0.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the zero-range process on a ring, which condenses above the critical density `1 / (b - 2)` for `b > 2`.

    A particle leaves a site with `n` particles at rate `u(n) = 1 + b / n`, to the right with weight `p` and to the
    left with weight `q`. The stationary distribution factorizes with the weights `f(n) = prod_k 1 / u(k)`, which
    decay as `n^-b`. Above the critical density a fluid at density `1 / (b - 2)` coexists with one site that holds the
    excess particles.

    Args:
        n_sites:
            The number of sites.

        density:
            The number of particles per site.

        b:
            The strength of the decrease of the departure rate with the occupation.

        p:
            The rate scale of hops to the right.

        q:
            The rate scale of hops to the left.

        random_state:
            The seed or generator of the initial sites of the particles.

    Returns:
        system:
            The model, with the occupations `sites.n`.

    Examples:
        >>> from skcomplex.models.physics.statistical_mechanics import zero_range_process
        >>> system = zero_range_process(10, density=2.0)
        >>> float(system.elements[0].values['n'].sum()), len(system.interactions)
        (20.0, 1)
    """
    occupations = _scatter(round(density * n_sites), n_sites, random_state)
    sites = LatticeElement('sites', element=ParticleSiteElement, space=LatticeSpace(n_sites), values={'n': occupations})
    return System(elements=[sites], interactions=[ZeroRangeHopsInteraction('hops', sites='sites', p=p, q=q, b=b)])