MC-07 Phase Transition

Introduction

The goal of this tutorial is to learn how to detect a second-order phase transition from finite-size simulations, using the venerable example of the 2d Ising model. No true phase transition can occur on a finite system. Nevertheless finite-size simulations show clear precursor signs of the phase transition which, combined with finite-size scaling, allow a precise determination of the universality class of a continuous phase transition.

Almost everything is known about the phase transition in the 2d Ising model since it is exactly solvable. In this tutotial, we will try to recover the location of the critical point, as well as critical exponents as if we would not know them in order to illustrate the methods. To make precise estimations requires quite some time. For this let us start the final simulation parameters in the background while doing the first part of the tutorial.

You can start the second simulation in the background with the parameter file parm7b and type:

parameter2xml parm7b
spinmc --Tmin 10  parm7b.in.xml &

or run the first part of the tutorial7b.py:

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
import pyalps.fit_wrapper as fw

prepare the input parameters

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

write the input file and run the simulation

input_file = pyalps.writeInputFiles('parm7b',parms) 
pyalps.runApplication('spinmc',input_file,Tmin=5)

Locate roughly the phase transition

First, let us make a rough temperature scan on small systems, in order to locate roughly the critical range. Use the parameter file parm7a and the command.

parameter2xml parm7a
spinmc --Tmin 5 parm7a.in.xml

Alternatively, you can run the simulations in python with the file tutorial7a.py:

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot

prepare the input parameters

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

write the input file and run the simulation

input_file = pyalps.writeInputFiles('parm7a',parms)
pyalps.runApplication('spinmc',input_file,Tmin=5)

Magnetization, susceptibility and specific heat

The order parameter describing the low-temperature phase in the Ising model is the magnetization per site m. Of course, on a finite-system, $\langle m\rangle$ is zero by symmetry since there is no true symmetry breaking, and one should rather look at the average absolute value of magnetization per site. In order to evaluate the data and load the observables we want to inspect run the following lines in the python shell:

pyalps.evaluateSpinMC(pyalps.getResultFiles(prefix='parm7a'))

load the observables and collect them as function of temperature T

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

and make a a plot of the |Magnetization| versus T:

plt.figure()
pyalps.plot.plot(magnetization_abs)
plt.xlabel('Temperature $T$')
plt.ylabel('Magnetization $|m|$')
plt.title('2D Ising model')

One can clearly see that the magnetization rises from its zero high-temperature value to its saturation value at low T. The temperature at which it happens is not so clear, as the upturn gets sharper as system size is increased. To have a clearer idea, let us look at fluctuations of the magnetization.

For this, let us consider the connected susceptibility $\chi = \beta N ( \langle m^2\rangle- \langle m\rangle^2)$, where N is the total number of spins and where we subtracted the square of the average absolute value of magnetization, hence the name connected.

plt.figure()
pyalps.plot.plot(connected_susc)
plt.xlabel('Temperature $T$')
plt.ylabel('Connected Susceptibility $\chi_c$')
plt.title('2D Ising model')

One observes a marked peak around T=2.2-2.4, where fluctuations are strongest. We note that the peaks tend to diverge with system size; as we will see later, this divergence can be characterized by a critical exponent.

This first rough estimation of the critical temperature is confirmed by looking at the behaviour of the specific heat, which characterizes fluctuations of energy : $C_v = \beta^2 N ( \langle e^2 \rangle - \langle e \rangle^2 )$, where $e$ is the internal energy per site.

plt.figure()
pyalps.plot.plot(spec_heat)
plt.xlabel('Temperature $T$')
plt.ylabel('Specific Heat $c_v$')
plt.title('2D Ising model')

We observe here a peak, albeit less marked, around the same values of temperature. We will account for the less pronounced divergence of the specific heat peak later too.

Binder Cumulant

Sometimes locating the maximum of a curve (such as the peaks in susceptibility or specific heat) is not always easy given the temperature grid used in the simulations. Moreover, and this can be seen in the previous plots, the peak temperature can drift with system size. This is also accounted for in the finite-size scaling theory outlined below.

An alternative efficient way of locating phase transitions is through the use of cumulants. Let us consider the ratio $U_4=\langle m^4\rangle /\angle m^2\rangle^2$, first introduced by Binder often called the Binder cumulant. At low temperature, this ratio tends to unity, whereas Gaussian fluctuations of the order parameter give a value 3 for $U_4$ in the high temperature phase. The non-trivial feature of the Binder cumulant resides in its value at the critical point, which can be shown to be universal and independent of the system size. Therefore, a plot of the Binder cumulant versus temperature allows a determination of Tc at the crossing points of the curves with different system sizes:

plt.figure()
pyalps.plot.plot(binder_u4)
plt.xlabel('Temperature $T$')
plt.ylabel('Binder Cumulant U4 $g$')
plt.title('2D Ising model')

From this last plot, we thus conclude that the phase transition is located in the finer range $T_c \in [2.2,2.3]$. One can also consider another cumulant $U_2 = \langle m^2\rangle / \langle |m|\rangle^2$, which is is left as an exercice.

Locate precisely the phase transition, first estimates of critical exponents, collapse plots

We can now try to perform a more precise determination of the nature of the phase transition using a finer grid of temperature points in the critical range T and using larger system sizes. For this, we will use the results of parm7b. As before we evaluate and load the specific observables from the result files./

pyalps.evaluateSpinMC(pyalps.getResultFiles(prefix='parm7b'))

load the observables and collect them as function of temperature T

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

Collapsing data

We start with performing the Binder cumulant crossing:

plt.figure()
pyalps.plot.plot(binder_u4)
plt.xlabel('Temperature $T$')
plt.ylabel('Binder Cumulant U4 $g$')
plt.title('2D Ising model')
plt.show()

What is the estimate of $T_c$ you obtain from this plot? There is more information in this plot. Indeed, the theory of finite size scaling indicates the following scaling form for the Binder cumulant $U_4 = f (L^{1/\nu} (T-T_c)/T_c)$, where $f$ is a universal function. For the different system sizes used in parm7b, plot the Binder cumulant as a function of $(T-T_c)/T_c$ where $T_c$ is the estimation you determined previously: all the curves should cross close to 0. Now try to find a good constant a such that when multiplied by $L^a$, all curves corresponding to system size $L$, collapse into a single master curve, i.e. the curve for $L=32$ should be multiplied by $32^a$, $L=48$ by $48^a$ etc.

Hint : Try a close to unity …

perform a data collapse of the Binder cumulant:

Tc=... #your estimate
a=...  #your estimate

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^a(T-T_c)/T_c$')
plt.ylabel('Binder Cumulant U4 $g$')
plt.title('2D Ising model')
plt.show()

What you are currently doing is a data collapse in the scaling regime, a famous technique to obtain critical exponents. In this case, you can read off the correlation length critical exponent : $\nu = 1/a$. The Binder cumulant collapse is pretty useful as it allows to determine $\nu$ independently of other critical exponents. A data collapse for other quantities often need to also scale the quantity on the y-axis by some other exponent $L^b$.

Critical exponents

Now consider the specific heat and connected susceptibility:

make a plot of the specific heat and connected susceptibility:

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()

Both connected susceptibility and specific heat peak at values of $T$ which are slightly different from the $T_c$ obtained with the Binder cumulant crossing. Such small differences are expected on small finite-systems : the effective critical temperature $T_c(L)$ that one can define at the position of the peaks is expected to drift with system size : $T_c(L) = T_c + A L^{-1/\nu}$. The constant $A$ can be different with different ways of defining $T_c(L)$. What about critical exponents? They can be read off from the value of the connected susceptibility or specific heat, either at $T_c(L)$ or at the $T_c$ obtained from the Binder cumulant. One expects the following scaling : $\chi (T_c) \sim L^{\gamma/\nu} and C_v (T_c) \sim L^{\alpha / \nu}$ For instance, locate the maximum of the connected susceptibility for the different curves, plot it as a function of $L$ and try a power-law fit. Note that by using relations between critical exponents $\gamma / \nu = 2 - \eta$. What value of $\eta$ do you obtain ?

make a fit of the connected susceptibility as a function of $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)
prefactor = pars[0].get()
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$ is %.4s' % gamma_nu)
plt.show()

Repeat the same procedure for the specific heat. What value of $\alpha$ do you obtain ? Note that $\alpha$ can in general be positive or negative, meaning that the specific heat need not to diverge at a continuous transition in contrast with the connected susceptibility.

make a fit of the specific heat as a function of $L$:

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)
prefactor = pars[0].get()
alpha_nu = pars[1].get()
plt.figure()
plt.plot(peak_sh.x, f(None, peak_sh.x, pars))
pyalps.plot.plot(peak_cs)
plt.xlabel('System Size $L$')
plt.ylabel('Specific Heat $c_v(T_c)$')
plt.title(r'2D Ising model, $\alpha$ is %.4s' % alpha_nu)
plt.show()

As an exercice, you can repeat the same analysis for the absolute magnetization, which scales as $L^{-\beta/\nu}$.

Also, try to perform a data collapse for the connected susceptibility and magnetization.

Hint : The corresponding scaling forms are $\chi = L^{2-\eta} g ( L^{1/\nu} (T-T_c)/T_c))$ and $|m| = L^{-\beta/\nu} h ( L^{1/\nu} (T-T_c)/T_c))$.

make a data collapse of the connected susceptibility as a function of $(T-Tc)/T_c$:

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_eta)

plt.figure()
pyalps.plot.plot(connected_susc)
plt.xlabel(' $L^a(T-T_c)/T_c$')
plt.ylabel(r'$L^{\gamma/\nu}\chi_c$')
plt.title('2D Ising model')
plt.show()

make a data collapse of the |magnetization| as a function of $(T-Tc)/T_c$

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)

plt.figure()
pyalps.plot.plot(magnetization_abs)
plt.xlabel(' $L^a(T-T_c)/T_c$')
plt.ylabel(r'$L^{-\beta/\nu} |m|$')
plt.title('2D Ising model')
plt.show()                                                                                                                                                                                  

Precise estimates of critical exponents

The exact critical exponents of the 2d Ising are : $\nu=1$, $\eta=1/4$, $\beta=1/8$ and $\alpha=0$. On this last exponent : the specific heat does diverge, but logarithmically, not as power-law and hence the specific heat exponent is 0 dimensional. How do you compare your estimates with the exact values ?

The quality of the estimates one obtains for the critical exponents and temperature increases considerably when using larger system sizes, which were not accessible during the time scale of this tutorial. Moreover, they also rely on the quality of the estimate of the critical temperature. For more precision on critical exponents, one can enforce the exact value of the critical temperature $T_c = 2 / \ln(1+\sqrt(2)) = 2.269\ldots$ and simulate larger system sizes, which we leave to you as an exercise. Since there is only one value of temperature as a parameter, the previous way of determining $\nu$ through a data collapse is not that helpful. The above scaling form for the Binder cumulant allows for another more direct determination of $\nu$. Derivating the Binder cumulant with respect to $T$, one easily sees that the derivative $dU_4/dT$ scales as $L^{1/\nu}$ at $T_c$. You can either obtain this derivative from a numerical differentiation of your data, or better, it can also be obtained as a thermodynamical average during the Monte Carlo simulations. Note that this requires good statistics. Once you have obtained the data, the exponent $\nu$ can be determined by a power-law fit of this derivative $dU_4/dT$ at T_c as a function of system size $L$ .