MC-07 相転移

このチュートリアルの目的は、2次元イジング模型という古典的な例を用いて、有限サイズのシミュレーションから二次相転移を検出する方法を学ぶことです。有限系では真の相転移は起こりませんが、有限サイズのシミュレーションには明確な前兆が現れます。これを有限サイズスケーリングと組み合わせることで、臨界温度と転移の普遍性クラスを精密に決定できます。

2次元イジング模型の相転移は厳密に解けるため、ほぼすべてが分かっています。ここでは手法を説明するために、臨界点の位置と臨界指数があたかも未知であるかのように求め直します。精密な推定には長時間のシミュレーションが必要なので、細かい温度グリッドの計算を最初にバックグラウンドで開始し、先に粗いグリッドの計算を解析します。

バックグラウンドでのシミュレーション

チュートリアルの残りの部分を進めている間に計算が走るよう、細かいグリッドのシミュレーションを今すぐ開始します。

コマンドライン

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'))  # 生のモンテカルロ出力から派生量（帯磁率、比熱、ビンダーキュムラント）を計算する
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$ がゼロに近いことを示しています。

ビンダーキュムラント

ビンダーキュムラントは、臨界温度を特定するためのより鋭い道具を与えてくれます。次の比を定義します。

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]$ と読み取れます。第2キュムラント $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'])

ビンダーキュムラントの交差

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$

有限サイズスケーリングは、ビンダーキュムラントが次のスケーリング形に従うことを予言します。

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$ の値を 1 の近くで試してください。

Tc = ...   # ビンダー交差から得たあなたの推定値
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}.

ピーク位置は $T_c(L) = T_c + A\,L^{-1/\nu}$ に従って $L$ とともにわずかにずれることに注意してください。したがってピークの評価は、ビンダー交差から得た真の $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$ は普遍関数です。ビンダーコラプスから得た $\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()

厳密値との比較

2次元イジング模型の厳密な臨界指数は次のとおりです。

\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$ を取り出す有効な方法はビンダーキュムラントの微分を使うことです。 $T_c$ において、有限サイズスケーリングは $dU_4/dT \sim L^{1/\nu}$ を意味するので、この微分を $L$ の関数としてべき則フィットすれば $\nu$ が直接得られます。この微分は、モンテカルロデータの数値微分によって求めることもできますし、シミュレーション中に熱力学的平均として測定することもできます（後者には良い統計が必要です）。

問題

ビンダーキュムラントの曲線はどの温度で交差しますか。あなたの推定値は厳密値 $T_c \approx 2.2692$ とどの程度一致しますか。
$\nu$ 、 $\eta$ 、 $\beta$ 、 $\alpha$ について、あなたの数値的な推定値は厳密値とどう比較できますか。
ビンダーキュムラントのデータコラプスから、 $\nu$ としてどんな値が得られますか。厳密値 $\nu = 1$ と比較してください。
帯磁率のピークは系のサイズとともに大きく増大する一方、比熱のピークははるかにゆっくりとしか増大しません。このことは $\alpha$ について何を意味しますか。
$\gamma/\nu = 2 - \eta$ を使うと、 $\eta$ としてどんな値が得られますか。厳密値 $\eta = 1/4$ と比較してください。
最も大きな系のサイズについて SWEEPS を増やしてみてください。データコラプスの質はどう改善されますか。
（発展） $dU_4/dT$ を数値的に計算し、 $L^{1/\nu}$ のスケーリングを確認してください。

MC-06 量子 Wang-Landau MC-08 量子相転移

MC-07 相転移

バックグラウンドでのシミュレーション

コマンドライン

Python

粗いスキャン: 相転移の位置を大まかに探す

コマンドライン

Python

磁化、帯磁率、比熱

ビンダーキュムラント

TcT_cTc​ と臨界指数の精密決定

ビンダーキュムラントの交差

データコラプスと相関長指数 ν\nuν

帯磁率と比熱のピークのスケーリング

帯磁率ピークのフィット

比熱ピークのフィット

帯磁率と磁化のデータコラプス

厳密値との比較

問題

$T_c$ と臨界指数の精密決定

データコラプスと相関長指数 $\nu$