Ising

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

• Name: ising
• Field: physics/statistical_mechanics
• Formalism: Monte Carlo sampling
• References: Ising 1925; Onsager 1944; Jensen 2023

Factory

@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)])

Components

class IsingSpinElement(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.

    Args:
        label:
            The label of the spin or the lattice group.

        visual:
            How the spin is drawn, or `None` for `default_visual_`.

        spin:
            The initial spin, 1 or -1.

        J:
            The coupling between neighbouring spins.

        h:
            The external field.

        temperature:
            The temperature in units of the coupling.

    Attributes:
        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
    """

    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

class IsingSamplingInteraction(SpinSamplingInteraction):
    """Monte Carlo sweeps of Ising spins `s = ±1` of local energy `-J s sum_j s_j - h s`, by flips.

    The coupling `J`, the field `h`, and the temperature are parameters of the spins.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        target:
            The path of the group of spins.

        state:
            The name of the state that holds the spin.

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

        schedule:
            `'checkerboard'` on a square lattice, `'colored'` on a network, or `'random-sequential'`.

        neighborhood:
            The neighbours of a spin, or `None` for the von Neumann neighbourhood. A network takes
            `AdjacencyNeighborhood()`.

    Examples:
        >>> from types import SimpleNamespace
        >>> from sage.all__sagemath_symbolics import var
        >>> from skcomplex.models.physics.common.spins import IsingSamplingInteraction
        >>> ising = IsingSamplingInteraction('ising', target='spins')
        >>> other = SimpleNamespace(spin=var('t'))
        >>> ising.energy(SimpleNamespace(target=SimpleNamespace(spin=var('s'), J=1, h=var('h'), other=other)))
        -h*s - s*neighbors(t, reduce='sum')
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        state: str = 'spin',
        rule: str = 'metropolis',
        schedule: str = 'checkerboard',
        neighborhood: BaseNeighborhood | None = None,
    ) -> None:
        super().__init__(
            label,
            target=target,
            state=state,
            energy=self._energy,
            rule=rule,
            schedule=schedule,
            neighborhood=neighborhood,
        )

    def _energy(self, namespace: SimpleNamespace) -> Expression:
        """Return the local energy `-J s sum_j s_j - h s` of a spin."""
        spin, other = getattr(namespace.target, self.state), getattr(namespace.target.other, self.state)
        return -namespace.target.J * spin * neighbors(other, over=self.neighborhood) - namespace.target.h * spin