Revision as of 03:25, 19 March 2012

Tutorial 02: Hybridization Expansion CT-HYB

まず最初に、連続時間量子モンテカルロコードについて説明します。ハイブリッド拡張アルゴリズムCT-HYBです。DFMTのFig.11(Georges it et al.)を例にします。この図の６つの曲線は、系が冷却時にどのように反強磁性相になるのか示しています。系はhalf filling近傍の相互作用でのBethe格子のHubbard模型です。

CT-HYBシミュレーションは反復処理ごとに約1分間実行されます。このシミュレーションを実行するためのパラメータ·ファイルはディレクトリtutorials/dmft-02-hybridizationにあります。

全DMFTチュートリアルでは、Vistrailsを使って計算を始めるか、直接コマンドラインから実行させるか、どちらも可能です。vistrailsスクリプトは、以下に説明するように、自動的にパラメータを生成し、実行し、結果をプロットします。DMFTパラメータを手動で実行すると、例えば、tutorials/dmft-02-hybridizationの'beta_14_U3_tsqrt2'ディレクトリで、dmftを'/opt/alps/bin/dmft hybrid.param'のように実行すると、同じ結果が得られます。パラメータは次の通りです。

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

バンド構造や、格子タイプのパラメータが無いことに注意してください。デフォルトではBethe格子が採用されていますが、変更が可能です。（後述）

シミュレーションを実行します。

dmft hybrid.param

10回の自己無どう着逐次近似法を実行します。プログラムを実行しているディレクトリで、グリーン関数ファイル G_tau_?と自己のエネルギー(selfenergy_?)ならびに、グリーン関数の周波数G_omega_?が出力されます。サンプルのG_tauは２つの値を持ちます。スピンアップとスピンダウンの列です。\betaの値は、負の密度です。系が反強磁性相で、エラーバーの外側という点で異なります。異なる\betaでのグリーン関数のプロットやGeorges it et al.のFig.11の結果と比較するために次のPythonシェルスクリプトを使用することができます。チュートリアル07では、離散時間量子モンテカルロプログラムHirsch Fyeコードで同じ結果を再現します。パラメータは同じです。

サンプルtutorials/dmft-02-hybridizationにパラメータファイル(*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)

シミュレーション後、チュートリアル07のHirsch FyeやDFMTreviewの結果と比較してみてください。また、次のチュートリアルの相互作用展開の結果とも比較してみてください。また、再計算するために、G0OMEGA_INPUTを定義することによって再開する箇所を決めることができます。例えば、G0omga_outputをG0_omega_inputにコピーし、 パラメータファイルでG0OMEGA_INPUT = G0_omega_inputと規程し、再計算させます。http://alps.comp-phys.org/mediawiki/index.php/ALPS_2_Tutorials:DMFT-01_Hirsch-Fye Tutorial DMFT-01 Hirsch-Fye]で、このスクリプトを用いてグリーン関数をプロットすることによって反強磁性相の相転移を見ることができます。

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)を長くしたり、多くのCPUを使用することによって統計結果を改善することができます。MPIでの実行は、mpirun -np procs dmft parameter_fileとしてください。詳細はMPIのman等を参考にしてください。

DFMTの自己無どう着逐次近似法の収束状況を確認したいなら、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()

収束が難しいself energyが収束していることが確認されることが望ましいです。長時間シミュレーションが必要なことに注意してください。