ALPS 2 Tutorials:MC-07 Phase Transition/ja

From ALPS
Jump to: navigation, search


最初に

このチュートリアルの目的は、有名な2次元イジングモデルを例にとり、有限サイズのシミュレーションから二次相転移を再現する方法を学ぶことです。有限系では真の相転移や臨界現象は起きませんが、前駆現象は明確にあらわれるので、有限サイズスケーリング法を用いることで、ユニバーサリティクラスを決定することができます。 相転移現象は2次元イジングモデルでは正確に解けることが知られています。このチュートリアルでは、臨界点並びに臨界指数の計算手順を説明します。このチュートリアルの計算は、ある程度の時間を要しますので、効率的に進めてください。

最初の計算をする前に、2番目の計算をあらかじめバックグラウンドで始めておくとよいでしょう。そのパラメータファイルはparm7bです。次のコマンドで実行開始してください。

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

Pythonでは、tutorial7b.py

import pyalps
import matplotlib.pyplot as plt
import pyalps.pyplot
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)

Vistrailを用いる場合は、mc-07-tutorial.vtを参照してください。

相転移のおおよその特定 

最初に、臨界点のおおよその範囲を特定するために、小さい系で温度による探索をおこないます。使用するパラメータファイルはparm7aです。

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

Pythonを用いる場合、スクリプトはtutorial7a.pyです。

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

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

磁化、帯磁率、比熱

低温におけるイジングモデルを記述する秩序変数は、サイトあたりの磁化m です。もちろん有限系では、平均値< m>は対称性により0となるので、磁化の平均絶対値に注目すべきです。 次の記述はPythonシェルによる、実行、データ解析、計算結果読み込みの一連の手順です。

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

Tに対する|Magnetization|のプロットをおこないます。

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

高温から低温へ行くと磁化が0からある飽和値まで上昇することが明らかにわかります。システムサイズを大きくするにつれて変化が鋭くなりますが、この変化が起きる温度はそれほど明白ではありません。より明らかにするために、磁化のゆらぎを見てみましょう。 そのために、帯磁率\chi = \beta .N .( < m^2>- <|m|>^2)を見てみます。Nは全スピン数です。

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. ( < e^2 >- < e >^2 ) (eはサイトあたりの内部エネルギー)をみることで、確かめることができます。

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

帯磁率のピークよりは弱いですが、同じ場所に比熱のピークを観測できます。比熱のピークが帯磁率のピークより弱い理由も、後で見ていきます。

Binderキュムラント

シミュレーションの結果、常に曲線の最大値(帯磁率や比熱のピーク)を容易に求められるとは限りません。また、前途のプロットでも見れたように、ピーク温度はシステムサイズによって変動することもあります。これらは後述する有限サイズスケーリング法によって対処が可能です。 相転移点を特定する他の効率的な方法にキュムラントを用いる方法があります。ここでBinderキュムラントと呼ばれる比率U_4=< m^4> /<m^2>^2を取り上げます。この量は、低温ではスピンがすべて揃うので1をとり、高温では秩序変数のガウシアンゆらぎから3となります。Binderキュムラントの特徴は、臨界点直上において、システムサイズによらずに、ユニバーサルなある値を取ることです。よって、サイズごとに温度に対するBinderキュムラントをプロットすると、その交点として転移点を決定できます。

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

このプロットから、相転移は T_c \in [2.2,2.3]の狭い範囲であると結論づけられます。また、別の課題として、キュムラントU_2 = < m^2> / <|m|>^2の場合も考えて見てください。

相転移場所のより精確な特定、臨界指数の評価、collapse plots

密な温度グリッド点を用い大規模系を計算することで、より正確な相転移の性質をシミュレーションすることができます。チュートリアルで用いるファイルはparm7bです。結果ファイルから観測量を取り出し解析します。

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

Bider キュムラントの交点を計算します。

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

プロットから転移点は得られましたか?どこにありましたか?さて、このプロットからは転移点以外にも様々な情報を得ることができます。有限サイズスケーリング理論からは、Binder キュムラントはU_4 = f (L^{1/\nu} (T-T_c)/T_c)のようにスケールすることが示されます。ここで、f はユニバーサルなスケーリング関数です。parm7b に用意された各システムサイズに対して、Binder キュムラントを(T-T_c)/T_c の関数としてプロットしてください。T_c には転移点の推定値を使います。この時、すべての曲線が横軸=0の近くで交わるはずです。次に、ある定数a を導入して、横軸にL^aをかけることで、すべての曲線が1本の曲線に乗るようにしてください。

Hint : a を1 の近くから探してみてください ...

以下の様なpython による解析コードを使うとよいでしょう。

#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']g
    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()

ここでやった、臨界領域で曲線を重ねる手法(collapse plots)は、臨界指数を得るための有名なテクニックとなっています。この場合、相関長の臨界指数を\nu = 1/aとして読み取ることができます。このようにBinder キュムラントを使うことで、\nuを他の臨界指数とは独立に得ることができます。他の物理量のスケーリングをする際には大抵、x 成分だけでなくy 成分に関してもL^bというスケーリングを行う必要が生まれます。

臨界指数

今度は、比熱と帯磁率に関する解析を行います。

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

帯磁率と比熱のピークは、Binder キュムラントから得られる転移点からすこしずれた位置に現れます。この小さな差は計算した系が小さいことから現れます。つまり、ピークの位置によって実効的な転移温度T_c(L)<math>を定義すると、その位置は真の転移点<math>T_cから有限サイズ効果によって、T_c(L) = T_c + A L^{-1/\nu} という形でずれます。AT_c(L) を決める手法に依存する定数です。 臨界指数は、帯磁率や比熱の、T_c(L) もしくはBinder 解析から得たT_c での値から読み取ることができます。これらは\chi (T_c) \sim L^{\gamma/\nu}C_v (T_c) \sim L^{\alpha / \nu}という形でスケーリングすることが期待されます。例として、帯磁率の各サイズに対する最大値を得てプロットし、べき乗則にしたがっていることを確かめてみましょう。 各臨界指数の間に成り立つ関係式、\gamma / \nu = 2 - \eta を用いて\etaの値を計算してください。

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

比熱に関しても同様の手順で解析をします。臨界指数\alpha の値はどうなりましたか? \alpha は一般に正にも負にもなりえることに注意してください。このことは、比熱は帯磁率と違い、連続相転移においても発散するとは限らないということを意味します。

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

練習として、L^{-\beta/\nu}としてのスケールする絶対磁化に対しても同様の解析をおこなうことができます。また、帯磁率や磁化についてもユニバーサル関数によるデータ解析をおこなってみてください。 hint: 対応するスケーリング関数は、\chi = L^{2-\eta} g ( L^{1/\nu} (T-T_c)/T_c))|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)/Tc:
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)/Tc
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()

臨界指数の正確な推定値

二次元イジングモデルの臨界指数の厳密解は\nu=1, \eta=1/4, \beta=1/8 and \alpha=0です。 最後の\alpha=0 ですが、これは比熱が発散するけれどもそれはベキ的ではなく対数的であるという事を意味します。この厳密解と、今までに計算で得られた値とを比べてみましょう。 今回は計算時間が長すぎるためにやらなかったような、より大規模な系で計算を行うことで、転移点や臨界指数の見積りの精度は大幅に上がります。さらに、臨界指数の見積り精度は、転移点の見積りの精度自体にも依存します。臨界指数をより精度よく得るためには、臨界点の厳密解T_c = 2 / \ln(1+\sqrt(2)) = 2.269\ldotsを用いて、より大きな系で計算する必要があります。これは練習問題とします。 前に行った、ユニバーサルスケーリング関数による\nu の計算は、あまり使いやすくありません。 別の方法として、先にやった帯磁率や比熱のように、Binder キュムラントの転移温度上でのスケーリング関数を用いることで、\nu を直接求めることができます。 Binder キュムラントの温度に関する導関数は、実際に微分することで容易にわかるように、転移点直上でL^{1/\nu} at T_cという形でスケールします。この量は、数値微分で求めることもできますし、よりよい方法として、モンテカルロ法によって直接求めることもできます。なお、この計算には十分なサンプル数が必要になります。 こうして転移点での微分量が得られたなら、そのシステムサイズ依存性をべき乗則でフィッティングすることで、臨界指数\nu が求まります。