MC-07 Phase Transition
The goal of this tutorial is to learn how to detect a second-order phase transition from finite-size simulations, using the classic example of the two-dimensional Ising model. No true phase transition can occur on a finite system, but finite-size simulations show clear precursor signs that — combined with finite-size scaling — allow a precise determination of the critical temperature and the universality class of the transition.
Almost everything is known about the phase transition in the 2D Ising model since it is exactly solvable. Here we recover the location of the critical point and the critical exponents as if they were unknown, to illustrate the methods. Extracting precise estimates requires long simulations, so we start the fine-grid run in the background at the outset and analyse the coarse-grid run first.
Background simulations
Start the fine-grid simulations now so they run while you work through the rest of the tutorial.
Command line
Download parm7b and run:
parameter2xml parm7b
spinmc --Tmin 10 parm7b.in.xml &The --Tmin 10 flag sets a checkpoint interval of 10 seconds, allowing the simulation to be safely interrupted and resumed.
Python
The first part of tutorial7b.py sets up and launches the same simulation:
import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
import pyalps.fit_wrapper as fw
parms = []
for l in [32, 48, 64]:
for t in [2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.30, 2.31, 2.32, 2.33, 2.34, 2.35]:
parms.append(
{
'LATTICE' : "square lattice",
'T' : t,
'J' : 1,
'THERMALIZATION' : 5000,
'SWEEPS' : 150000,
'UPDATE' : "cluster",
'MODEL' : "Ising",
'L' : l
}
)
input_file = pyalps.writeInputFiles('parm7b', parms)
pyalps.runApplication('spinmc', input_file, Tmin=5)Rough scan: locating the phase transition
We first make a coarse temperature scan on small systems to locate the critical region.
Command line
Download parm7a and run:
parameter2xml parm7a
spinmc --Tmin 5 parm7a.in.xmlPython
Using tutorial7a.py:
import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
parms = []
for l in [4, 8, 16]:
for t in [5.0, 4.5, 4.0, 3.5, 3.0, 2.9, 2.8, 2.7]:
parms.append(
{
'LATTICE' : "square lattice",
'T' : t,
'J' : 1,
'THERMALIZATION' : 1000,
'SWEEPS' : 400000,
'UPDATE' : "cluster",
'MODEL' : "Ising",
'L' : l
}
)
for t in [2.6, 2.5, 2.4, 2.3, 2.2, 2.1, 2.0, 1.9, 1.8, 1.7, 1.6, 1.5, 1.2]:
parms.append(
{
'LATTICE' : "square lattice",
'T' : t,
'J' : 1,
'THERMALIZATION' : 1000,
'SWEEPS' : 40000,
'UPDATE' : "cluster",
'MODEL' : "Ising",
'L' : l
}
)
input_file = pyalps.writeInputFiles('parm7a', parms)
pyalps.runApplication('spinmc', input_file, Tmin=5) # Tmin: minimum wall time in seconds between checkpointsMagnetization, susceptibility, and specific heat
The order parameter of the Ising transition is the magnetization per site $m$. On a finite system $\langle m \rangle = 0$ by symmetry, so we plot the average absolute value $\langle |m| \rangle$ instead. Load the results and plot them:
pyalps.evaluateSpinMC(pyalps.getResultFiles(prefix='parm7a')) # computes derived quantities (susceptibility, specific heat, Binder cumulant) from raw MC output
data = pyalps.loadMeasurements(
pyalps.getResultFiles(prefix='parm7a'),
['|Magnetization|', 'Connected Susceptibility', 'Specific Heat', 'Binder Cumulant', 'Binder Cumulant U2']
)
magnetization_abs = pyalps.collectXY(data, x='T', y='|Magnetization|', foreach=['L'])
connected_susc = pyalps.collectXY(data, x='T', y='Connected Susceptibility', foreach=['L'])
spec_heat = pyalps.collectXY(data, x='T', y='Specific Heat', foreach=['L'])
binder_u4 = pyalps.collectXY(data, x='T', y='Binder Cumulant', foreach=['L'])
binder_u2 = pyalps.collectXY(data, x='T', y='Binder Cumulant U2', foreach=['L'])
plt.figure()
pyalps.plot.plot(magnetization_abs)
plt.xlabel('Temperature $T$')
plt.ylabel('Magnetization $|m|$')
plt.title('2D Ising model')The magnetization rises from zero at high temperature to its saturation value at low $T$. The upturn sharpens with increasing system size, but the precise transition temperature is hard to read off directly.
To locate the transition more precisely, examine fluctuations. The connected susceptibility $\chi = \beta N (\langle m^2 \rangle - \langle |m| \rangle^2)$ measures magnetization fluctuations and diverges at the critical point:
plt.figure()
pyalps.plot.plot(connected_susc)
plt.xlabel('Temperature $T$')
plt.ylabel('Connected Susceptibility $\chi_c$')
plt.title('2D Ising model')A marked peak appears around $T = 2.2$–$2.4$, and the peak grows with system size. The specific heat $c_v = \beta^2 N (\langle e^2 \rangle - \langle e \rangle^2)$, where $e$ is the internal energy per site, shows a similar but less pronounced peak in the same range:
plt.figure()
pyalps.plot.plot(spec_heat)
plt.xlabel('Temperature $T$')
plt.ylabel('Specific Heat $c_v$')
plt.title('2D Ising model')
plt.show()The specific heat peak near $T \approx 2.3$ grows only weakly with system size — much more slowly than the susceptibility — indicating that the exponent $\alpha$ is close to zero.
Binder cumulant
The Binder cumulant provides a sharper tool for locating the critical temperature. Define the ratio
$$U_4 = \frac{\langle m^4 \rangle}{\langle m^2 \rangle^2},$$
which equals 1 in the fully ordered low-temperature phase and 3 for a Gaussian (high-temperature) order-parameter distribution. At the critical point its value is universal: it is the same for all system sizes, so curves for different $L$ all cross at $T_c$.
plt.figure()
pyalps.plot.plot(binder_u4)
plt.xlabel('Temperature $T$')
plt.ylabel('Binder Cumulant $U_4$')
plt.title('2D Ising model')
plt.show()From the crossing of the curves we can read off $T_c \in [2.2, 2.3]$.
The second cumulant $U_2 = \langle m^2 \rangle / \langle |m| \rangle^2$ has the same crossing property; plotting it from the already-loaded binder_u2 data is left as an exercise.
Precise determination of $T_c$ and critical exponents
Using the finer temperature grid and larger system sizes from parm7b, we can extract $T_c$ and the critical exponents more accurately.
Load the results:
pyalps.evaluateSpinMC(pyalps.getResultFiles(prefix='parm7b')) # computes derived quantities from raw MC output
data = pyalps.loadMeasurements(
pyalps.getResultFiles(prefix='parm7b'),
['|Magnetization|', 'Connected Susceptibility', 'Specific Heat', 'Binder Cumulant', 'Binder Cumulant U2']
)
magnetization_abs = pyalps.collectXY(data, x='T', y='|Magnetization|', foreach=['L'])
connected_susc = pyalps.collectXY(data, x='T', y='Connected Susceptibility', foreach=['L'])
spec_heat = pyalps.collectXY(data, x='T', y='Specific Heat', foreach=['L'])
binder_u4 = pyalps.collectXY(data, x='T', y='Binder Cumulant', foreach=['L'])
binder_u2 = pyalps.collectXY(data, x='T', y='Binder Cumulant U2', foreach=['L'])Binder cumulant crossing
plt.figure()
pyalps.plot.plot(binder_u4)
plt.xlabel('Temperature $T$')
plt.ylabel('Binder Cumulant $U_4$')
plt.title('2D Ising model')
plt.show()The crossing of the $L = 32$, $48$, and $64$ curves gives a refined estimate of $T_c$.
Data collapse and the correlation-length exponent $\nu$
Finite-size scaling predicts that the Binder cumulant obeys the scaling form
$$U_4 = f!\left(L^{1/\nu},\frac{T - T_c}{T_c}\right),$$
where $f$ is a universal function. Shifting the horizontal axis to $(T - T_c)/T_c$ should make the curves cross near zero. Multiplying by $L^a$ should then collapse all curves onto a single master curve when $a = 1/\nu$. Try values of $a$ close to 1:
Tc = ... # your estimate from the Binder crossing
a = ... # try values near 1.0
for d in binder_u4:
d.x -= Tc
d.x = d.x / Tc
l = d.props['L']
d.x = d.x * pow(float(l), a)
plt.figure()
pyalps.plot.plot(binder_u4)
plt.xlabel('$L^{1/\\nu}(T-T_c)/T_c$')
plt.ylabel('Binder Cumulant $U_4$')
plt.title('2D Ising model')
plt.show()When the collapse is good, the correlation-length exponent is $\nu = 1/a$. This determination of $\nu$ is independent of all other critical exponents.
Peak scaling of susceptibility and specific heat
Finite-size scaling also predicts how the peak values of the susceptibility and specific heat grow with system size at $T_c$:
$$\chi(T_c) \sim L^{\gamma/\nu}, \qquad C_v(T_c) \sim L^{\alpha/\nu}.$$
Note that the peak positions drift slightly with $L$ as $T_c(L) = T_c + A,L^{-1/\nu}$, so you can evaluate the peaks either at the true $T_c$ from the Binder crossing or at the numerically determined peak temperature for each $L$.
plt.figure()
pyalps.plot.plot(connected_susc)
plt.xlabel('Temperature $T$')
plt.ylabel('Connected Susceptibility $\chi_c$')
plt.title('2D Ising model')
plt.figure()
pyalps.plot.plot(spec_heat)
plt.xlabel('Temperature $T$')
plt.ylabel('Specific Heat $c_v$')
plt.title('2D Ising model')
plt.show()Susceptibility peak fit
Extract the peak values and fit a power law in $L$:
cs_mean = []
for q in connected_susc:
cs_mean.append(np.array([d.mean for d in q.y]))
peak_cs = pyalps.DataSet()
peak_cs.props = pyalps.dict_intersect([q.props for q in connected_susc])
peak_cs.y = np.array([np.max(q) for q in cs_mean])
peak_cs.x = np.array([q.props['L'] for q in connected_susc])
sel = np.argsort(peak_cs.x)
peak_cs.y = peak_cs.y[sel]
peak_cs.x = peak_cs.x[sel]
pars = [fw.Parameter(1), fw.Parameter(1)]
f = lambda self, x, pars: pars[0]() * np.power(x, pars[1]())
fw.fit(None, f, pars, peak_cs.y, peak_cs.x)
gamma_nu = pars[1].get()
plt.figure()
plt.plot(peak_cs.x, f(None, peak_cs.x, pars))
pyalps.plot.plot(peak_cs)
plt.xlabel('System Size $L$')
plt.ylabel('Connected Susceptibility $\\chi_c(T_c)$')
plt.title('2D Ising model, $\\gamma/\\nu$ = %.4s' % gamma_nu)
plt.show()The fitted exponent gives $\gamma/\nu$. Using the hyperscaling relation $\gamma/\nu = 2 - \eta$, you can read off the anomalous dimension $\eta$.
Specific heat peak fit
Repeat the same analysis for the specific heat:
sh_mean = []
for q in spec_heat:
sh_mean.append(np.array([d.mean for d in q.y]))
peak_sh = pyalps.DataSet()
peak_sh.props = pyalps.dict_intersect([q.props for q in spec_heat])
peak_sh.y = np.array([np.max(q) for q in sh_mean])
peak_sh.x = np.array([q.props['L'] for q in spec_heat])
sel = np.argsort(peak_sh.x)
peak_sh.y = peak_sh.y[sel]
peak_sh.x = peak_sh.x[sel]
pars = [fw.Parameter(1), fw.Parameter(1)]
f = lambda self, x, pars: pars[0]() * np.power(x, pars[1]())
fw.fit(None, f, pars, peak_sh.y, peak_sh.x)
alpha_nu = pars[1].get()
plt.figure()
plt.plot(peak_sh.x, f(None, peak_sh.x, pars))
pyalps.plot.plot(peak_sh)
plt.xlabel('System Size $L$')
plt.ylabel('Specific Heat $c_v(T_c)$')
plt.title('2D Ising model, $\\alpha/\\nu$ = %.4s' % alpha_nu)
plt.show()Note that $\alpha$ can be positive or negative: the specific heat need not diverge as a power law at a continuous transition.
Data collapse for susceptibility and magnetization
The full scaling forms are
$$\chi = L^{2-\eta}, g!\left(L^{1/\nu},\frac{T-T_c}{T_c}\right), \qquad |m| = L^{-\beta/\nu}, h!\left(L^{1/\nu},\frac{T-T_c}{T_c}\right),$$
where $g$ and $h$ are universal functions. Using your estimates of $\nu$ and $T_c$ from the Binder collapse, tune $2-\eta$ to collapse the susceptibility curves:
for d in connected_susc:
d.x -= Tc
d.x = d.x / Tc
l = d.props['L']
d.x = d.x * pow(float(l), a)
two_minus_eta = ... # your estimate
for d in connected_susc:
l = d.props['L']
d.y = d.y / pow(float(l), two_minus_eta)
plt.figure()
pyalps.plot.plot(connected_susc)
plt.xlabel('$L^{1/\\nu}(T-T_c)/T_c$')
plt.ylabel('$L^{-(2-\\eta)}\\chi_c$')
plt.title('2D Ising model')
plt.show()And tune $\beta/\nu$ to collapse the magnetization curves:
for d in magnetization_abs:
d.x -= Tc
d.x = d.x / Tc
l = d.props['L']
d.x = d.x * pow(float(l), a)
beta_over_nu = ... # your estimate
for d in magnetization_abs:
l = d.props['L']
d.y = d.y / pow(float(l), -beta_over_nu) # dividing by L^(-β/ν) = multiplying by L^(β/ν)
plt.figure()
pyalps.plot.plot(magnetization_abs)
plt.xlabel('$L^{1/\\nu}(T-T_c)/T_c$')
plt.ylabel('$L^{\\beta/\\nu}|m|$')
plt.title('2D Ising model')
plt.show()Exact values and comparison
The exact critical exponents of the 2D Ising model are
$$\nu = 1, \quad \eta = \tfrac{1}{4}, \quad \beta = \tfrac{1}{8}, \quad \alpha = 0.$$
For $\alpha = 0$: the specific heat diverges logarithmically rather than as a power law, so a power-law fit returns an exponent close to zero but not exactly zero.
For more precise estimates, fix the critical temperature to the exact value $T_c = 2/\ln(1+\sqrt{2}) \approx 2.2692$ and simulate larger system sizes. With $T_c$ fixed, a useful way to extract $\nu$ is via the derivative of the Binder cumulant. At $T_c$, finite-size scaling implies $dU_4/dT \sim L^{1/\nu}$, so a power-law fit of this derivative as a function of $L$ gives $\nu$ directly. The derivative can be obtained by numerical differentiation of the Monte Carlo data, or measured as a thermodynamic average during the simulation (which requires good statistics).
Questions
- At what temperature do the Binder cumulant curves cross? How does your estimate compare with the exact value $T_c \approx 2.2692$?
- How do your numerical estimates of $\nu$, $\eta$, $\beta$, and $\alpha$ compare with the exact values?
- What value of $\nu$ do you extract from the Binder cumulant data collapse? Compare with the exact value $\nu = 1$.
- The susceptibility peak grows strongly with system size while the specific heat peak grows much more slowly. What does this imply about $\alpha$?
- Using $\gamma/\nu = 2 - \eta$, what value of $\eta$ do you obtain? Compare with the exact value $\eta = 1/4$.
- Try increasing
SWEEPSfor the largest system sizes. How does this improve the quality of the data collapse? - (Advanced) Compute $dU_4/dT$ numerically and verify the $L^{1/\nu}$ scaling.