Ising on a graph

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

• Name: graph_ising
• Field: physics/statistical_mechanics
• Formalism: Monte Carlo sampling
• References: Argun et al. 2021; Dorogovtsev, Goltsev, and Mendes 2008

Factory

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

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