Wang–Landau Algorithm¶
Wang–Landau sampling estimates the density of states directly by forcing a broad random walk in energy space. It greatly reduces canonical trapping and critical-slowing-down problems and gives access to the free energy, entropy, and other thermodynamic quantities over a wide temperature range. The energy-space walk still has a finite tunneling and convergence time that grows with system size.
Consider first a four-site Ising model.
Many configurations can have the same energy. Their number is the density of states $g(E)$.
For the square Ising model, let $L$ be the linear lattice size and $N_s=L^2$ the number of spins. For even $L$ with periodic boundary conditions,
$$ E_{\min}=-2JN_s,\qquad E_{\max}=2JN_s. $$
The energies differ by multiples of $4J$, but the levels $E_{\min}+4J$ and $E_{\max}-4J$ are absent. The number of allowed energies is of order $N_s$, whereas the number of configurations is $2^{N_s}$. Many energy levels are therefore highly degenerate.
We use units with $J=k_B=1$.
A canonical simulation at temperature $T$ visits energy $E$ with probability
$$ P_T(E)\propto g(E)e^{-E/T}. $$
The idea of Wang–Landau sampling is to estimate $g(E)$ directly, without choosing a simulation temperature. If the density of states is known,
$$ Z(T)=\sum_E g(E)e^{-E/T},\qquad F(T)=-T\ln Z(T). $$
Assign every configuration $X$ the statistical weight
$$ w(X)\propto\frac{1}{g(E(X))}. $$
Since there are $g(E)$ configurations at energy $E$, the probability of visiting that energy is
$$ P_{\mathrm{WL}}(E)\propto g(E)\frac{1}{g(E)}=\text{constant}. $$
The resulting energy histogram is flat when the estimated density of states is accurate. Since $g(E)$ is initially unknown, the estimate $\widehat g(E)$ is modified continuously during the walk.
For a proposed move from energy $E_1$ to $E_2$, use
$$ P(E_1\rightarrow E_2)=\min\left[1,\frac{\widehat g(E_1)}{\widehat g(E_2)}\right]. $$
Moves toward a smaller estimated density of states are always accepted; moves toward a larger density are accepted less often.
Algorithm¶
Initialize $\ln\widehat g(E)=0$ and $H(E)=0$ for every allowed energy. Start with $\ln f=1$, so $f=e$.
Propose a random spin flip and accept it with the probability above.
Whether the proposal is accepted or rejected, update the energy currently occupied by the walker:
$$ \ln\widehat g(E_{\mathrm{current}})\leftarrow\ln\widehat g(E_{\mathrm{current}})+\ln f, $$
$$ H(E_{\mathrm{current}})\leftarrow H(E_{\mathrm{current}})+1. $$
Periodically test whether the histogram is sufficiently flat. With flatness parameter $p$, use
$$ \min_EH(E)>p\,\frac{1}{N_E}\sum_EH(E), $$
where $N_E$ is the number of allowed energy levels. Typical choices are $p=0.8$ or $0.9$.
When the histogram is flat, reduce the modification factor, $\ln f\leftarrow\ln f/2$, reset $H(E)=0$, and continue from the current configuration.
Stop when $\ln f$ is smaller than a chosen target, for example $10^{-6}$ for a lecture demonstration or $10^{-8}$ for a more accurate calculation.
A flat histogram at one stage is a convergence diagnostic, not a proof that $\widehat g(E)$ is exact. The statistical error decreases as $\ln f$ is reduced, and independent runs provide a practical estimate of the remaining uncertainty.
Primary references: Wang and Landau, Phys. Rev. Lett. 86, 2050 (2001) and Wang and Landau, Phys. Rev. E 64, 056101 (2001).
Why store $\ln g(E)$?¶
The density of states grows exponentially with the number of spins and quickly exceeds floating-point range. We therefore store the microcanonical entropy
$$ S(E)=\ln g(E). $$
The multiplicative Wang–Landau update becomes an additive update,
$$ g(E)\leftarrow g(E)f \quad\Longleftrightarrow\quad S(E)\leftarrow S(E)+\ln f. $$
Only ratios of the density of states enter the acceptance probability, so $S(E)$ is initially determined only up to an additive constant. After the simulation, we fix that constant using the known ground-state degeneracy $g(E_{\min})=2$.
Implementation of the Wang–Landau algorithm¶
import matplotlib.pyplot as plt
import numpy as np
from numba import jit
@jit(nopython=True)
def CEnergy(lattice):
"""Return the total energy of an L x L Ising lattice with J=1."""
L = lattice.shape[0]
energy = 0
for i in range(L):
for j in range(L):
spin = lattice[i,j]
field = (lattice[(i+1)%L,j] + lattice[i,(j+1)%L]
+ lattice[(i-1)%L,j] + lattice[i,(j-1)%L])
energy -= spin*field
return energy/2.0 # every bond was counted twice
def RandomL(L, seed=509):
"""Return a random lattice, corresponding to infinite temperature."""
rng = np.random.default_rng(seed)
return rng.choice((-1,1), size=(L,L))
def PrepareEnergies(L):
"""Allowed energies and a fast energy-to-index lookup for even L."""
n_spins = L*L
energies = np.arange(-2*n_spins, 2*n_spins+1, 4, dtype=int)
# A single overturned spin breaks four bonds, so E_min+4 is absent.
# By ferro/antiferromagnetic symmetry, E_max-4 is absent as well.
energies = np.delete(energies, [1,-2])
energy_min, energy_max = energies[0],energies[-1]
index_of_energy = -np.ones(energy_max-energy_min+1, dtype=int)
for index,energy in enumerate(energies):
index_of_energy[energy-energy_min] = index
return energies,index_of_energy,energy_min
energies,index_of_energy,energy_min = PrepareEnergies(8)
print("number of energy levels =",len(energies))
print("minimum energy =",energy_min)
print("allowed energies:")
print(energies)
number of energy levels = 63 minimum energy = -128 allowed energies: [-128 -120 -116 -112 -108 -104 -100 -96 -92 -88 -84 -80 -76 -72 -68 -64 -60 -56 -52 -48 -44 -40 -36 -32 -28 -24 -20 -16 -12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 96 100 104 108 112 116 120 128]
def WangLandau(L,flatness=0.8,lnf_final=1e-6,
max_steps=50_000_000,check_every=10_000,seed=509):
"""Estimate log g(E) for an even L x L Ising lattice."""
energies,index_of_energy,_ = PrepareEnergies(L)
lattice = RandomL(L,seed)
log_dos,histogram,lnf,steps,stages = RunWangLandau(
max_steps,energies,lattice,index_of_energy,
flatness,lnf_final,check_every,seed+1,
)
return energies,log_dos,histogram,lnf,steps,stages
@jit(nopython=True)
def RunWangLandau(max_steps,energies,lattice,index_of_energy,
flatness,lnf_final,check_every,seed):
"""Perform Wang--Landau updates until lnf reaches its target."""
np.random.seed(seed)
L = lattice.shape[0]
n_spins = L*L
energy = int(CEnergy(lattice))
energy_min = energies[0]
log_dos = np.zeros(len(energies))
histogram = np.zeros(len(energies))
lnf = 1.0
stages = 0
steps = 0
for step in range(max_steps):
site = int(np.random.random()*n_spins)
i,j = site//L,site%L
spin = lattice[i,j]
field = (lattice[(i+1)%L,j] + lattice[i,(j+1)%L]
+ lattice[(i-1)%L,j] + lattice[i,(j-1)%L])
energy_new = energy + int(2*spin*field)
old_index = index_of_energy[energy-energy_min]
new_index = index_of_energy[energy_new-energy_min]
log_acceptance = log_dos[old_index]-log_dos[new_index]
if log_acceptance >= 0.0 or np.random.random() < np.exp(log_acceptance):
lattice[i,j] = -spin
energy = energy_new
# Update the energy occupied after every attempted move.
current_index = index_of_energy[energy-energy_min]
histogram[current_index] += 1.0
log_dos[current_index] += lnf
steps = step+1
if steps%check_every == 0:
mean_histogram = np.mean(histogram)
if np.min(histogram) > flatness*mean_histogram:
lnf *= 0.5
stages += 1
print(steps,"steps: histogram is flat; lnf =",lnf)
if lnf <= lnf_final:
break
histogram[:] = 0.0
return log_dos,histogram,lnf,steps,stages
L = 8
flatness = 0.8
lnf_final = 1e-6
energies,lngE,histogram,lnf,steps,stages = WangLandau(
L,flatness=flatness,lnf_final=lnf_final,
max_steps=50_000_000,check_every=10_000,seed=509,
)
print(f"finished after {steps:,} attempted moves and {stages} flat stages")
if lnf > lnf_final:
print("The maximum number of steps was reached before convergence.")
70000 steps: histogram is flat; lnf = 0.5 90000 steps: histogram is flat; lnf = 0.25 140000 steps: histogram is flat; lnf = 0.125 180000 steps: histogram is flat; lnf = 0.0625 210000 steps: histogram is flat; lnf = 0.03125 270000 steps: histogram is flat; lnf = 0.015625 350000 steps: histogram is flat; lnf = 0.0078125 440000 steps: histogram is flat; lnf = 0.00390625 520000 steps: histogram is flat; lnf = 0.001953125 700000 steps: histogram is flat; lnf = 0.0009765625 970000 steps: histogram is flat; lnf = 0.00048828125 1190000 steps: histogram is flat; lnf = 0.000244140625 1430000 steps: histogram is flat; lnf = 0.0001220703125 2440000 steps: histogram is flat; lnf = 6.103515625e-05 3220000 steps: histogram is flat; lnf = 3.0517578125e-05 3850000 steps: histogram is flat; lnf = 1.52587890625e-05 4370000 steps: histogram is flat; lnf = 7.62939453125e-06 4900000 steps: histogram is flat; lnf = 3.814697265625e-06 5530000 steps: histogram is flat; lnf = 1.9073486328125e-06 5850000 steps: histogram is flat; lnf = 9.5367431640625e-07 finished after 5,850,000 attempted moves and 20 flat stages
print("unnormalized endpoint values:",lngE[0],lngE[-1])
unnormalized endpoint values: 1632.8566226959229 1632.7482376098633
# Fix the otherwise arbitrary normalization with g(E_min)=2.
lngE += np.log(2.0)-lngE[0]
def LogSumExp(values):
"""Return log(sum(exp(values))) without overflow."""
maximum = np.max(values)
return maximum+np.log(np.sum(np.exp(values-maximum)))
print("g(E_min) =",np.exp(lngE[0]),"(exact: 2)")
print("g(E_max) =",np.exp(lngE[-1]),"(exact: 2)")
print("log sum_E g(E) =",LogSumExp(lngE),
"(exact:",L*L*np.log(2.0),")")
g(E_min) = 1.99999999999989 (exact: 2) g(E_max) = 1.7945639982058836 (exact: 2) log sum_E g(E) = 44.45775009232398 (exact: 44.3614195558365 )
fig,ax = plt.subplots(figsize=(7,4))
ax.plot(energies/(L*L),lngE/(L*L),"o-",markersize=3)
ax.set(xlabel=r"energy per spin $E/N_s$",
ylabel=r"entropy per spin $\ln g(E)/N_s$")
ax.grid(alpha=0.3)
plt.show()
Thermodynamics from the density of states¶
Once $g(E)$ is known, the same simulation gives thermodynamics at every temperature:
$$ Z(T)=\sum_E g(E)e^{-E/T}, \qquad \langle E^n\rangle=\frac{1}{Z}\sum_E g(E)E^n e^{-E/T}. $$
For $N_s=L^2$ spins,
$$ e=\frac{\langle E\rangle}{N_s},\qquad c_V=\frac{\langle E^2\rangle-\langle E\rangle^2}{N_sT^2}, $$
$$ f=-\frac{T\ln Z}{N_s},\qquad s=\frac{\ln Z+\langle E\rangle/T}{N_s}. $$
Directly exponentiating $\ln g(E)-E/T$ can overflow or underflow. We therefore compute $\ln Z$ with the log-sum-exp identity. If $a_E=\ln g(E)-E/T$ and $a_{\max}=\max_Ea_E$, then
$$ \ln Z=a_{\max}+\ln\left[\sum_E e^{a_E-a_{\max}}\right]. $$
The normalized Boltzmann probabilities are then $p_E=e^{a_E-\ln Z}$. For a finite lattice, the two degenerate ground states give $e\rightarrow-2$ and $s\rightarrow\ln2/N_s$ as $T\rightarrow0$; the entropy per spin vanishes only in the thermodynamic limit. As $T\rightarrow\infty$, $e\rightarrow0$ and $s\rightarrow\ln2$.
def Thermodynamics(T,log_dos,energies,L):
"""Return energy, heat capacity, free energy, and entropy per spin."""
n_spins = L*L
log_weights = log_dos-energies/T
log_partition = LogSumExp(log_weights)
probabilities = np.exp(log_weights-log_partition)
mean_energy = np.sum(probabilities*energies)
mean_energy2 = np.sum(probabilities*energies**2)
heat_capacity = (mean_energy2-mean_energy**2)/(n_spins*T**2)
free_energy = -T*log_partition/n_spins
entropy = (log_partition+mean_energy/T)/n_spins
return mean_energy/n_spins,heat_capacity,free_energy,entropy
temperatures = np.linspace(0.05,5.0,500)
thermodynamics = np.array([
Thermodynamics(T,lngE,energies,L) for T in temperatures
])
labels = (r"$\langle E\rangle/N_s$",r"$c_V$",r"$F/N_s$",r"$S/N_s$")
fig,axes = plt.subplots(2,2,figsize=(10,7),sharex=True)
for column,(ax,label) in enumerate(zip(axes.flat,labels)):
ax.plot(temperatures,thermodynamics[:,column])
ax.axvline(2/np.log(1+np.sqrt(2)),color="0.5",linestyle="--",
label=r"$T_c$" if column == 1 else None)
ax.set(ylabel=label)
ax.grid(alpha=0.3)
axes[0,1].legend()
axes[1,0].set_xlabel("temperature T")
axes[1,1].set_xlabel("temperature T")
fig.tight_layout()
plt.show()
low_temperature = temperatures[0]
high_temperature = 100.0
low_T = Thermodynamics(low_temperature,lngE,energies,L)
high_T = Thermodynamics(high_temperature,lngE,energies,L)
print(f"low-T energy per spin: {low_T[0]:.6f} (exact limit: -2)")
print(f"low-T entropy per spin: {low_T[3]:.6f} "
f"(finite-size limit: {np.log(2)/(L*L):.6f})")
print(f"at T={high_temperature:g}, entropy per spin: {high_T[3]:.6f} "
f"(infinite-temperature limit: {np.log(2):.6f})")
low-T energy per spin: -2.000000 (exact limit: -2) low-T entropy per spin: 0.010830 (finite-size limit: 0.010830) at T=100, entropy per spin: 0.694553 (infinite-temperature limit: 0.693147)
Homework¶
- Estimate the statistical uncertainty in $S(E)=\ln g(E)$ for the $32\times32$ Ising model. Perform at least five independent Wang–Landau runs, using a different random seed in every run and the same $\ln f$ stopping criterion. Normalize each result with $S_i(E_{\min})=\ln2$ before comparing or averaging the runs.
For $R$ independent runs, calculate
$$ \overline{S}(E)=\frac{1}{R}\sum_{i=1}^R S_i(E), $$
and the sample standard deviation between runs,
$$ s_S(E)=\sqrt{\frac{1}{R-1}\sum_{i=1}^R\left[S_i(E)-\overline{S}(E)\right]^2}. $$
The standard error of the mean is $s_S(E)/\sqrt{R}$. Plot $S_i(E)-\overline{S}(E)$ and $s_S(E)$. State clearly whether an error bar represents the scatter of individual runs or the uncertainty in their mean.
- Beale's Exact Distribution of Energies in the Two-Dimensional Ising Model, Phys. Rev. Lett. 76, 78 (1996), gives exact finite-lattice densities of states. For the $32\times32$ lattice, the first ten values are
g_exact = np.array([
2,2048,4096,1057792,4218880,371621888,2191790080,
100903637504,768629792768,22748079183872,
],dtype=float)
Plot $\overline{S}(E)$ and $\ln g_{\mathrm{exact}}(E)$ for these energies. In a second panel, plot their difference together with the standard error $s_S(E)/\sqrt{R}$. Is the exact result statistically consistent with the mean of your independent runs?
- Use your best estimate of $g(E)$ to calculate the energy, heat capacity, free energy, and entropy as functions of temperature. Check the low- and high-temperature limits stated above.