Sherrington-Kirkpatrick spin glass

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

• Name: sherrington_kirkpatrick
• Field: physics/statistical_mechanics
• Formalism: jump process
• References: Sherrington and Kirkpatrick 1975; Bertin 2021; Mezard, Parisi, and Virasoro 1987

Factory

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

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