Fully connected Ising model with Glauber dynamics

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

• Name: curie_weiss
• Field: physics/statistical_mechanics
• Formalism: jump process
• References: Bertin 2021; Glauber 1963

Factory

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

Components

class FullyConnectedSpinElement(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`.

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

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

        spin:
            The initial spin, 1 or -1.

        J:
            The scale of the couplings.

        h:
            The external field.

        temperature:
            The temperature.

        couplings:
            The coupling matrix `C`, one row and one column per spin, or `None` for the uniform couplings `1 / N`.

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

    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

    def _field(self, s: SimpleNamespace) -> Expression:
        """Return the local field of a spin, without its own contribution."""
        if self.couplings is None:
            return s.J * (mean(s.spin) - s.spin / count()) + s.h
        return s.J * total(s.spin, weights=np.asarray(self.couplings, dtype=float)) + s.h

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the flips of up spins and of down spins at the Glauber rate.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        rate = 1 / (1 + exp(2 * s.spin * self._field(s) / s.temperature))
        return [JumpLaw(where(s.spin > 0, rate, 0), {s.spin: -2}), JumpLaw(where(s.spin < 0, rate, 0), {s.spin: 2})]