ALPS 2 Tutorials:MC-07 Phase Transition/ja

From ALPS
Revision as of 09:35, 12 March 2012 by Kota (talk | contribs) (Collapsing data)

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から低Tでの飽和値まで上昇することが明らかにわかります。システムサイズを大きくなるにつれて、上昇が鋭くなるにつれて、これらが起きる温度はそれほど明白ではありません。 次に、磁化の変動を見てみます。

ここで、磁化率\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を取り上げます。低温度で、この比率は、低温度では一定ですが、高温度ではパラメータのガウシアンゆらぎU_4は3となります。Binderキュムラントの特徴は、その値が、システムサイズから独立していることにあります。すなわち、温度に対するBinderキュムラントのプロットは、Tcの決定に非常に有益です。


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 cumulant交差の計算をおこないます。

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

このプロットから

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


現在、何をやっているかというと、臨界指数を得るため有名なテクニックで、スケーリング用のデータを調整します。この場合、相関距離の臨界指数を読み取ることが可能です。\nu = 1/aです。\nuを決定するためにBinder cumulantを考慮することは有益です。

臨界指数

比熱と関連感受性に関して、チュートリアルをおこないます。

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

バインダーキュムラント横断で得られたT_cとは少し異なりますTの値で接続された感受性と比熱ピークの両方。効果的な臨界温度T_cの(L)</ math>は一つのシステムサイズでドリフトすると予想されるピークの位置で定義することができます:<math>T_c(L) = T_c + A L^{-1/\nu}のような小さな違いは、小さな有限のシステムで期待されています。定数A T_cの(L)</ math>を定義するさまざまな方法とは異なる場合があります。臨界指数についてはどうですか?彼らは<math>T_c(L)で、またはバインダーキュムラントから得られたT_cのいずれかで、接続された感受性や比熱の値から読み取ることができます。一つには、次のスケーリングを期待しています。
\chi (T_c) \sim L^{\gamma/\nu} and C_v (T_c) \sim L^{\alpha / \nu} 例えば、別の曲線に接続された感受性の最大値を見つける L </ math>はの関数としてプロットし、べき乗則近似を試してみてください。注意<math>\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()

臨界指数の正確な推定値

2Dイジングの厳密な臨界指数は、次のとおりです。\nu=1、\eta=1/4、\beta=1/8はと\alpha=0は。この最後の指数上の特定の熱が発散していますが、対数、べき乗則、したがって、特定の熱指数が0次元ではありませんように。どのように正確な値を使用して見積もりを比較するのですか? このチュートリアルの時間スケールの間にはアクセスできませんでした大規模なシステムのサイズを、使用するときにいずれかの臨界指数と温度を取得する推定の質が大幅に増加します。また、彼らはまた、臨界温度の推定値の品質に依存しています。臨界指数で、より高精度、1つは臨界温度のT_c = 2 / \ln(1+\sqrt(2)) = 2.269\ldotsの正確な値を適用し、我々は演習としてあなたに残して大規模なシステムのサイズを、シミュレートすることができます。パラメータは、データの崩壊を通じて\nuをを決定する以前の方法として、温度の値が一つだけあるのでそれは役立ちません。バインダーのキュムラントは、上記のスケーリング形は\ NUの別のより直接的な測定が可能になります。派生的Tに対するバインダーのキュムラントは、1つは容易に認識しているT_cでのL^{1/\nu}としてデリバティブdU_4/dTに比例します。あなたはどちらのデータ、またはそれ以上の数値微分からこの誘導体を得ることができ、また、モンテカルロシミュレーションの間に熱力学的平均値として取得することができます。これは良いの統計情報を必要としないことに注意してください。いったんデータを取得した、指数\nuを、システムのサイズLをの関数としてのT_cで、このデリバティブdU_4/dT のべき乗則近似によって決定することができる。