Revision as of 06:34, 19 March 2012

Single site DMFTのNeel転移

このチュートリアルでは、DFMTのFig.11Georges it et al.を再現します。この6つの曲線は、half filling近傍で相互作用を持つBethe格子上でのHubbard模型を表しています。温度が下がると反強磁性相になります。

このサンプル計算は、コマンドライン、Pythonスクリプト、Vistrailsでの実行に対応しています。Vistrailは、自動的にパラメータファイルを生成し、計算を実行し、結果をプロットします。DMFTパラメーターを手動で設定する場合は、tutorials/dmft-0j-xxx内の'beta_14_U3_tsqrt2'ディレクトリを参照してください。'/opt/alps/bin/dmft xxx.param'で実行してください。

Hybridization 拡張CT-HYB

Hybridization拡張アルゴリズムCT-HYB、まず連続時間量子モンテカルロコードを実行させます。CT-HYBシミュレーションは、反復処理毎に約１分実行します。実行用のパラメータファイルはディレクトリtutorials/dmft-02-hybridizationにあります。

主なパラメータは次の通りです。

SEED = 0; //Monte Carlo Random Number Seed 
THERMALIZATION = 1000;  Thermalization Sweeps 
SWEEPS = 1000000; Total Sweeps to be computed 
MAX_TIME = 60;  Maximum time to run the simulation 
BETA = 12.;  Inverse temperature 
SITES = 1;  This is a single site DMFT simulation, so Sites is 1 
N = 16;  Number of time slices (you will see that this parameter is rather small) 
NMATSUBARA = 500;  The number of Matsubara frequencies 
U = 3;  Interaction energy 
t = 1;  hopping parameter. For the Bethe lattice considered here $W=2D=4t$
MU = 0;  Chemical potential 
H = 0;  Magnetic field 
SYMMETRIZATION = 0;  We are not enforcing a paramagnetic self consistency condition 
SOLVER = Hybridization;  The Hybridization solver

次のようにタイプし、シミュレーションを実行させます。

dmft hybrid.param

このサンプルプログラムでは10回の繰り返し計算をおこないます。計算が実行されると、出力ディレクトリに、グリーン関数ファイルG_tau_?、自己エネルギー(selfenergy_?)、周波数空間でのグリーン関数G_omega_?が生成されます。G_tauは２つの値をもちます。spin-upとspin-downです。\beta値は負の密度です。次に示すPythonスクリプトを用いて計算し、他の\betaでのグリーン関数をプロットしたり、Georges it et al.のFig.11と比較してみてください。

Hirsch-Fye のセクションで、離散時間量子モンテカルロコード：Hirsch Fyeコードを用いて再現計算をおこないます。パラメータは同じです。

計算実行後、tutorials/dmft-02-hybridizationディレクトリにパラメータファイル(called *tsqrt2)が生成されます。DFMTシミュレーションをVistrailsで行う場合は、dmft-02-hybridization.vt、Python用のスクリプトは tutorial2a.pyになります。

import pyalps
import numpy as np
import matplotlib.pyplot as plt
import pyalps.pyplot

#prepare the input parameters
parms=[]
for b in [6.,8.,10.,12.,14.,16.]:
    parms.append(
            {
              'ANTIFERROMAGNET'     : 1,
              'CONVERGED'           : 0.03,
              'F'                   : 10,
              'FLAVORS'             : 2,
              'H'                   : 0,
              'H_INIT'              : 0.2,
              'MAX_IT'              : 10,
              'MAX_TIME'            : 60,
              'MU'                  : 0,
              'N'                   : 1000,
              'NMATSUBARA'          : 1000,
              'N_FLIP'              : 0,
              'N_MEAS'              : 10000,
              'N_MOVE'              : 0,
              'N_ORDER'             : 50,
              'N_SHIFT'             : 0,
              'OMEGA_LOOP'          : 1,
              'OVERLAP'             : 0,
              'SEED'                : 0,
              'SITES'               : 1,
              'SOLVER'              : 'Hybridization',
              'SYMMETRIZATION'      : 0,
              'TOLERANCE'           : 0.01,
              'U'                   : 3,
              't'                   : 0.707106781186547,
              'SWEEPS'              : 100000000,
              'THERMALIZATION'      : 1000,
              'BETA'                : b,
              'CHECKPOINT'          : 'solverdump_beta_'+str(b)+'.task1.out.h5',
              'G0TAU_INPUT'         : 'G0_tau_input_beta_'+str(b)
            }
        )

#write the input file and run the simulation
for p in parms:
    input_file = pyalps.writeParameterFile('parm_beta_'+str(p['BETA']),p)
    res = pyalps.runDMFT(input_file)

これらのシミュレーションが実行後、 Hirsch-Fyeセクションの Hirsch FyeやDMFTの結果と比較してみてい下さい。また、 Interaction Expansion CT-INTのセクションの 相互作用拡張の結果とも比較してみてください。G0OMEGA_INPUTを設定することで計算のスタートポイントを定義することができます。例えば、G0omga_output を G0_omega_inputにコピーし specify G0OMEGA_INPUT = G0_omega_inputと定義してみてください。

flavors=parms[0]['FLAVORS']
listobs=[]   
for f in range(0,flavors):
    listobs.append('Green_'+str(f))

ll=pyalps.load.Hdf5Loader()
data = ll.ReadMeasurementFromFile(pyalps.getResultFiles(pattern='parm_beta_*h5'), respath='/simulation/results/G_tau', measurements=listobs, verbose=True)
for d in pyalps.flatten(data):
    d.x = d.x*d.props["BETA"]/float(d.props["N"])
    d.props['label'] = r'$\beta=$'+str(d.props['BETA'])
plt.figure()
plt.xlabel(r'$\tau$')
plt.ylabel(r'$G(\tau)$')
plt.title('Hubbard model on the Bethe lattice')
pyalps.pyplot.plot(data)
plt.legend()
plt.show()

計算結果は比較的ノイズが多いことに気付くかと思います。理由は、測定結果を描写するのに、この様な高温度では展開数が非常に小さいからです。計算時間(MAX_TIME)を長くしたり、MPIを使って並列計算をおこなったりして計算の精度を向上させることができます。試してみてください。

DMFTの自己無憧着の収束を確認したい場合は、tutorial2b.pyを用いて各繰り返し計算でのグリーン関数をプロットしてみてください。

ll=pyalps.load.Hdf5Loader()
for p in parms:
    data = ll.ReadDMFTIterations(pyalps.getResultFiles(pattern='parm_beta_'+str(p['BETA'])+'.h5'), measurements=listobs, verbose=True)
    grouped = pyalps.groupSets(pyalps.flatten(data), ['iteration'])
    nd=[]
    for group in grouped:
        r = pyalps.DataSet()
        r.y = np.array(group[0].y)
        r.x = np.array([e*group[0].props['BETA']/float(group[0].props['N']) for e in group[0].x])
        r.props = group[0].props
        r.props['label'] = r.props['iteration']
        nd.append( r )
    plt.figure()
    plt.xlabel(r'$\tau$')
    plt.ylabel(r'$G(\tau)$')
    plt.title(r'\beta = %.4s' %nd[0].props['BETA'])
    pyalps.pyplot.plot(nd)
    plt.legend()
    plt.show()

自己エネルギーの収束性を確認する最も適した方法です。グリーン関数、自己エネルギーの収束には長いシミュレーションが必要なことに注意してください。

相互作用展開CT-INT

CT-INTコードを元にした Hybridization Expansion CT-HYB での計算例について述べます。このコードは相互作用エネルギーの展開をおこないます。パラメータファイルはtutorials/dmft-03-interactionディレクトリ内にあります。Pythonスクリプト、Vistrails用のファイルは、それぞれtutorial3a.py、tutorial3b.py、 Vistrails fileです。

Hirsch Fye

離散時間モンテカルロコードHirsch Fyeの計算をおこない、連続時間との比較をおこないます。Hirsch Fyeアルゴリズムの詳細はこちらを参照してください。このコードはオープンソースとして公開されています。さらなる情報はBlümer's PhDにあります。多くのアルゴリズムが開発され（Alvarez08、 Nukala09など）、連続時間アルゴリズムが主流となっています。

Hirsch Fyeシミュレーションのチュートリアルサンプルはtutorials/dmft-01-hirschfyeディレクトリ内を参照してください。

パラメータファイルな主な部分は、次の通りです。

SEED = 0; //Monte Carlo Random Number Seed 
THERMALIZATION = 10000;  Thermalization Sweeps 
SWEEPS = 1000000; Total Sweeps to be computed 
MAX_TIME = 60;  Maximum time to run the simulation 
BETA = 12.;  Inverse temperature 
SITES = 1;  This is a single site DMFT simulation, so Sites is 1 
N = 16;  Number of time slices (you will see that this parameter is rather small) 
NMATSUBARA = 500;  The number of Matsubara frequencies 
U = 3;  Interaction energy 
t = 1;  hopping parameter. For the Bethe lattice considered here $W=2D=4t$
MU = 0;  Chemical potential 
H = 0;  Magnetic field 
SYMMETRIZATION = 0;  We are not enforcing a paramagnetic self consistency condition 
SOLVER = /opt/alps/bin/hirschfye;  The path to the external Hirsch Fye solver

計算を実行させます。

dmft hirschfye.param

もしくは、Pythonスクリプトtutorial1a.pyを使用してください。

import pyalps
import numpy as np
import matplotlib.pyplot as plt
import pyalps.pyplot

#prepare the input parameters 
parms=[]
for b in [6.,8.,10.,12.,14.,16.]: 
    parms.append(
            { 
              'ANTIFERROMAGNET'     : 1,
              'CONVERGED'           : 0.03,
              'FLAVORS'             : 2,
              'H'                   : 0,
              'MAX_IT'              : 10,
              'MAX_TIME'            : 60,
              'MU'                  : 0,
              'N'                   : 16,
              'NMATSUBARA'          : 500, 
              'OMEGA_LOOP'          : 1,
              'SEED'                : 0, 
              'SITES'               : 1,
              'SOLVER'              : '/opt/alps/bin/hirschfye',
              'SYMMETRIZATION'      : 0,
              'TOLERANCE'           : 0.01,
              'U'                   : 3,
              't'                   : 0.707106781186547,
              'SWEEPS'              : 1000000,
              'THERMALIZATION'      : 10000,
              'BETA'                : b,
              'G0OMEGA_INPUT'       : 'G0_omega_input_beta_'+str(b),
              'BASENAME'            : 'hirschfye.param'
            }
        )

#write the input file and run the simulation
for p in parms:
    input_file = pyalps.writeParameterFile('parm_beta_'+str(p['BETA']),p)
    res = pyalps.runDMFT(input_file)

コードは10回の反復計算をおこないます。出力ディレクトリに、グリーン関数ファイルG_tau_?、自己エネルギー(selfenergy_?)、周波数空間でのグリーン関数G_omega_?が生成されます。G_tauは２つの値をもちます。spin-upとspin-downです。\beta値は負の密度です。次に示すPythonスクリプトを用いて計算し、他の\betaでのグリーン関数をプロットしたり、Georges it et al.のFig.11と比較してみてください。

flavors=parms[0]['FLAVORS']
listobs=[]   
for f in range(0,flavors):
    listobs.append('Green_'+str(f))

ll=pyalps.load.Hdf5Loader()
data = ll.ReadMeasurementFromFile(pyalps.getResultFiles(pattern='parm_beta_*h5'), respath='/simulation/results/G_tau', measurements=listobs, verbose=True)
for d in pyalps.flatten(data):
    d.x = d.x*d.props["BETA"]/float(d.props["N"])
    d.props['label'] = r'$\beta=$'+str(d.props['BETA'])
plt.figure()
plt.xlabel(r'$\tau$')
plt.ylabel(r'$G(\tau)$')
plt.title('Hubbard model on the Bethe lattice')
pyalps.pyplot.plot(data)
plt.legend()
plt.show()

離散時間の方法では、Hirsch Fye法は誤差が顕著にあらわれます。大きな<tt>Nを用いて計算してください。また、DMFTシミュレーションではほとんど収束された入力ファイルを用いてシミュレーションをおこなっています。G0_omega_inputファイルを削除して、フリーな状態から計算して、収束状況を確認してみてください。

DMFTシミュレーションはVistrailsの利用が可能です。dmft-01-hirschfye.vtを参照してください。