User Guide

This guide covers the physics implemented by montecarlo and how to apply it to your own graphs.

The Ising model

For a graph \(G = (V, E)\) with edge weights \(J_{ij}\), montecarlo represents an Ising model with Hamiltonian

\[H(\alpha) = \sum_{(i,j) \in E} J_{ij}\, s_i s_j + \sum_{i \in V} \mu_i s_i,\]

where \(s_i = +1\) if site \(i\) is “up” (bit value 1) and \(s_i = -1\) if it is “down” (bit value 0). Each configuration \(\alpha\) is represented by a BitString.

Representing configurations

The BitString class wraps a NumPy array of 0/1 entries and exposes utilities for flipping bits, counting up/down spins, and converting to and from a decimal integer index:

import montecarlo

bs = montecarlo.BitString(N=4)
bs.set_integer_config(11)   # binary 1011
print(bs)                   # "1011"
print(bs.on(), bs.off())    # 3, 1
bs.flip_site(0)             # now "0011"

Building a Hamiltonian

Hamiltonians are constructed from a NetworkX graph with weight attributes on each edge. Optional local fields can be set with set_mu:

import networkx as nx
import montecarlo

N = 6
G = nx.Graph()
G.add_nodes_from(range(N))
G.add_edges_from([(i, (i + 1) % N) for i in range(N)])
for e in G.edges:
    G.edges[e]['weight'] = 2.0

ham = montecarlo.IsingHamiltonian(G)
ham.set_mu([0.1] * N)

bs = montecarlo.BitString(N)
bs.set_config([1, 0, 1, 0, 1, 0])
print(ham.energy(bs))

Exact thermodynamic averages

Given the Boltzmann distribution

\[P(\alpha) = \frac{e^{-\beta E(\alpha)}}{Z}, \qquad Z = \sum_{\alpha'} e^{-\beta E(\alpha')},\]

an expectation value of an observable \(A\) is

\[\langle A \rangle = \sum_\alpha A(\alpha) P(\alpha).\]

The compute_average_values() method enumerates all \(2^N\) configurations and returns

\[\langle E \rangle, \quad \langle M \rangle, \quad \frac{\langle E^2 \rangle - \langle E \rangle^2}{T^2}, \quad \frac{\langle M^2 \rangle - \langle M \rangle^2}{T},\]

i.e. the average energy, magnetization, heat capacity, and magnetic susceptibility:

E, M, HC, MS = ham.compute_average_values(T=1.0)

This is exact, but its cost grows as \(2^N\) and becomes infeasible beyond roughly 20 sites.

Metropolis sampling

For larger systems, MonteCarlo performs Metropolis sampling. From the current configuration \(\alpha\), a single-site flip proposes a new configuration \(\beta\), accepted with probability

\[\begin{split}W(\alpha \to \beta) = \begin{cases} 1 & \text{if } E(\beta) \le E(\alpha) \\ e^{-\beta\, [E(\beta) - E(\alpha)]} & \text{otherwise.} \end{cases}\end{split}\]

The resulting Markov chain converges to the Boltzmann distribution, so expectation values reduce to a simple arithmetic average over collected samples:

\[\langle A \rangle \approx \frac{1}{M} \sum_{\alpha \in \text{trajectory}} A(\alpha).\]

import numpy as np

mc = montecarlo.MonteCarlo(ham)
E_samples, M_samples = mc.run(T=2.0, n_samples=10_000, n_burn=200)

Eavg, Estd = np.mean(E_samples), np.std(E_samples)
Mavg, Mstd = np.mean(M_samples), np.std(M_samples)

HC = Estd**2 / 2.0**2
MS = Mstd**2 / 2.0

The n_burn argument discards an initial number of sweeps to allow the chain to equilibrate before recording samples. Each “sweep” attempts a flip on every site once.

Temperature scans

A common workflow is to compute observables across a range of temperatures and locate the critical temperature, where the heat capacity and magnetic susceptibility peak:

import numpy as np
import matplotlib.pyplot as plt

Ts = np.linspace(0.1, 10.0, 100)
Es, Ms, HCs, MSs = [], [], [], []
for T in Ts:
    E, M, HC, MS = ham.compute_average_values(T)
    Es.append(E); Ms.append(M); HCs.append(HC); MSs.append(MS)

plt.plot(Ts, Es, label="Energy")
plt.plot(Ts, Ms, label="Magnetization")
plt.plot(Ts, HCs, label="Heat capacity")
plt.plot(Ts, MSs, label="Susceptibility")
plt.xlabel("Temperature"); plt.legend(); plt.show()

Tc = Ts[int(np.argmax(MSs))]
print(f"Critical temperature ~ {Tc:.2f}")