MC-07 相变

教程资源

蒙特卡洛模拟

MC-07 相变

本教程的目标是学习如何从有限尺寸模拟中探测二级相变，所用的经典例子是二维伊辛模型。有限系统上不可能出现真正的相变，但有限尺寸模拟会显示出清晰的前兆特征——结合有限尺寸标度，就能精确确定相变的临界温度和普适类。

由于二维伊辛模型可以精确求解，关于它相变的几乎一切都已为人所知。这里我们假装并不知道这些结果，重新确定临界点的位置和临界指数，以此说明相关方法。获得精确的估计需要很长的模拟时间，因此我们一开始就在后台启动细网格的运行，先分析粗网格的结果。

后台模拟

现在就启动细网格模拟，让它在你阅读本教程其余部分时持续运行。

命令行

下载 parm7b 并运行：

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

--Tmin 10 选项将检查点间隔设为 10 秒，这样模拟可以安全地中断和续算。

Python

tutorial7b.py 的第一部分会设置并启动同样的模拟：

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)

粗扫描：定位相变

我们首先在小系统上做一次粗略的温度扫描，以确定临界区域的大致位置。

命令行

下载 parm7a 并运行：

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

Python

使用 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：两次检查点之间的最短实际耗时（秒）

磁化强度、磁化率与比热

伊辛相变的序参量是每格点的磁化强度 $m$ 。在有限系统上，由于对称性有 $\langle m \rangle = 0$ ，因此我们改为绘制平均绝对值 $\langle |m| \rangle$ 。载入结果并作图：

pyalps.evaluateSpinMC(pyalps.getResultFiles(prefix='parm7a'))  # 由原始蒙特卡洛输出计算导出量（磁化率、比热、Binder 累积量）
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')

磁化强度从高温下的零值上升到低 $T$ 时的饱和值。随着系统尺寸增大，上升变得更陡峭，但直接读出精确的相变温度仍然很困难。

要更精确地定位相变，需要考察涨落。连通磁化率 $\chi = \beta N (\langle m^2 \rangle - \langle |m| \rangle^2)$ 刻画磁化强度的涨落，并在临界点处（于热力学极限下）发散：

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

在 $T = 2.2$ – $2.4$ 附近出现一个明显的峰，且峰值随系统尺寸增大而增长。比热 $c_v = \beta^2 N (\langle e^2 \rangle - \langle e \rangle^2)$ （其中 $e$ 是每格点的内能）在同一区间内显示出类似但不那么显著的峰：

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

$T \approx 2.3$ 附近的比热峰随系统尺寸只有很弱的增长——远慢于磁化率——这表明临界指数 $\alpha$ 接近零。

Binder 累积量

Binder 累积量为定位临界温度提供了更锐利的工具。定义比值

U_4 = \frac{\langle m^4 \rangle}{\langle m^2 \rangle^2},

它在完全有序的低温相中等于 1，而对于高斯型（高温）序参量分布则等于 3。在临界点处它的取值是普适的：对所有系统尺寸都相同，因此不同 $L$ 的曲线都会在 $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()

由曲线的交点可以读出 $T_c \in [2.2, 2.3]$ 。第二累积量 $U_2 = \langle m^2 \rangle / \langle |m| \rangle^2$ 具有相同的交叉性质；用已经载入的 binder_u2 数据把它画出来，留作练习。

精确确定 $T_c$ 与临界指数

利用 parm7b 中更细的温度网格和更大的系统尺寸，我们可以更准确地提取 $T_c$ 和临界指数。

载入结果：

pyalps.evaluateSpinMC(pyalps.getResultFiles(prefix='parm7b'))  # 由原始蒙特卡洛输出计算导出量
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 累积量的交叉

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

$L = 32$ 、 $48$ 和 $64$ 三条曲线的交点给出 $T_c$ 的精细估计。

数据坍缩与关联长度指数 $\nu$

有限尺寸标度预言 Binder 累积量服从如下标度形式：

U_4 = f\!\left(L^{1/\nu}\,\frac{T - T_c}{T_c}\right),

其中 $f$ 是普适函数。把横轴平移为 $(T - T_c)/T_c$ 后，各条曲线应当在零点附近相交。再乘以 $L^a$ ，当 $a = 1/\nu$ 时所有曲线应坍缩到同一条主曲线上。尝试取接近 1 的 $a$ 值：

Tc = ...   # 你由 Binder 交点得到的估计值
a  = ...   # 尝试接近 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()

当坍缩效果良好时，关联长度指数为 $\nu = 1/a$ 。这种确定 $\nu$ 的方法与其他所有临界指数无关。

磁化率与比热的峰值标度

有限尺寸标度还预言了在 $T_c$ 处磁化率和比热的峰值如何随系统尺寸增长：

\chi(T_c) \sim L^{\gamma/\nu}, \qquad C_v(T_c) \sim L^{\alpha/\nu}.

注意峰的位置会随 $L$ 略有漂移，遵循 $T_c(L) = T_c + A\,L^{-1/\nu}$ ，因此你既可以在由 Binder 交点得到的真实 $T_c$ 处取峰值，也可以在每个 $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()

磁化率峰值拟合

提取峰值并对 $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()

拟合得到的指数即 $\gamma/\nu$ 。利用超标度关系 $\gamma/\nu = 2 - \eta$ ，即可读出反常量纲 $\eta$ 。

比热峰值拟合

对比热重复同样的分析：

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

注意 $\alpha$ 可以为正也可以为负：在连续相变处，比热未必以幂律形式发散。

磁化率与磁化强度的数据坍缩

完整的标度形式为：

\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),

其中 $g$ 和 $h$ 是普适函数。利用你从 Binder 坍缩得到的 $\nu$ 和 $T_c$ 估计值，调节 $2-\eta$ 使磁化率曲线坍缩：

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 = ...   # 你的估计值
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()

再调节 $\beta/\nu$ 使磁化强度曲线坍缩：

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 = ...   # 你的估计值
for d in magnetization_abs:
    l   = d.props['L']
    d.y = d.y / pow(float(l), -beta_over_nu)  # 除以 L^(-β/ν) 等于乘以 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()

精确值与比较

二维伊辛模型的精确临界指数为：

\nu = 1, \quad \eta = \tfrac{1}{4}, \quad \beta = \tfrac{1}{8}, \quad \alpha = 0.

关于 $\alpha = 0$ ：比热是对数发散而非幂律发散，因此幂律拟合会给出一个接近零但不严格为零的指数。

若要得到更精确的估计，可以把临界温度固定为精确值 $T_c = 2/\ln(1+\sqrt{2}) \approx 2.2692$ ，并模拟更大的系统尺寸。在 $T_c$ 固定的情况下，提取 $\nu$ 的一个有效办法是利用 Binder 累积量的导数。在 $T_c$ 处，有限尺寸标度给出 $dU_4/dT \sim L^{1/\nu}$ ，因此对该导数关于 $L$ 作幂律拟合即可直接得到 $\nu$ 。这个导数既可以通过对蒙特卡洛数据作数值微分获得，也可以在模拟过程中作为热力学平均量直接测量（这需要良好的统计精度）。

问题

Binder 累积量曲线在什么温度相交？你的估计值与精确值 $T_c \approx 2.2692$ 相比如何？
你对 $\nu$ 、 $\eta$ 、 $\beta$ 和 $\alpha$ 的数值估计与精确值相比如何？
从 Binder 累积量的数据坍缩中你得到的 $\nu$ 是多少？与精确值 $\nu = 1$ 比较。
磁化率的峰随系统尺寸强烈增长，而比热的峰增长得慢得多。这对 $\alpha$ 意味着什么？
利用 $\gamma/\nu = 2 - \eta$ ，你得到的 $\eta$ 是多少？与精确值 $\eta = 1/4$ 比较。
尝试增大最大系统尺寸的 SWEEPS。这如何改善数据坍缩的质量？
（进阶）数值计算 $dU_4/dT$ 并验证 $L^{1/\nu}$ 的标度关系。

MC-06 量子王-朗道 MC-08 量子相变

MC-07 相变

后台模拟

命令行

Python

粗扫描：定位相变

命令行

Python

磁化强度、磁化率与比热

Binder 累积量

精确确定 TcT_cTc​ 与临界指数

Binder 累积量的交叉

数据坍缩与关联长度指数 ν\nuν

磁化率与比热的峰值标度

磁化率峰值拟合

比热峰值拟合

磁化率与磁化强度的数据坍缩

精确值与比较

问题

精确确定 $T_c$ 与临界指数

数据坍缩与关联长度指数 $\nu$