ALPS 2 Tutorials:ED-06 FullDiagonalization/ja

From ALPS
Jump to: navigation, search

量子ハミルトニアンの完全対角化

このチュートリアルでは、fulldiag を使って量子系の全対角化を行い、 得られたスペクトルから熱力学的量を計算します。

1次元スピンモデルの熱力学

spin-1 ハイゼンベルグ鎖

最初に、S=1 反強磁性ハイゼンベルグ鎖の計算をおこないます。

コマンドラインでの実行

parm6aは、8サイトでのS=1 反強磁性ハイゼンベルグ鎖の完全対角化の入力設定ファイルです。

LATTICE="chain lattice"
MODEL="spin"
local_S = 1
J       = 1
CONSERVED_QUANTUMNUMBERS="Sz"
{L = 8}

CONSERVED_QUANTUMNUMBERS パラメータでハミルトニアンと交換する物理量を指定することで、 ヒルベルト空間を不変部分空間に分割して、それぞれ独立して対角化できます。

以下の標準的な計算手順に従い、fulldiag を用いて量子ハミルトニアンの計算をおこないます。

parameter2xml parm1
fulldiag --write-xml parm1.in.xml

出力ファイルには全固有値・固有ベクトルが記述されています。ここから、fulldiag_evaluate を用いて熱力学的量・磁気的量のグラフを作成するXML ファイルを作る事ができます。

fulldiag_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T  0.1  parm6a.task1.out.xml

次のXML形式のファイルが生成されます。

parm6a.task1.plot.energy.xml
parm6a.task1.plot.free_energy.xml
parm6a.task1.plot.entropy.xml
parm6a.task1.plot.specific_heat.xml
parm6a.task1.plot.uniform_susceptibility.xml
parm6a.task1.plot.magnetization.xml

fulldiag_evaluate によって生成されたXMLファイルから該当の計算結果を抽出するために、plot2text ツールを使用します。得られるデータはプレーンテキスト形式なので、お好きなツールで扱う事ができます。例えば、エネルギー密度と温度の2次元プロットは、次の用に記述します。

plot2text parm1.task1.plot.energy.xml

他の物理量も同様に変換できます。

プロットツールにGraceを使用する場合は、plot2xmgrツールを使用してGrace用のファイルを直接生成できます。エネルギーと温度のプロットでは、次のように入力してください。

plot2xmgr parm1.task1.plot.energy.xml > energy.agr

同様に、gnuplotではplot2gp、テキスト形式ではplot2textを使用してください。推奨はPythonを使用する方法です。

Pythonを用いた実行

Pythonで計算をする場合、Pythonスクリプトtutorial6a.pyを使います。 これまで通り、必要なモジュールをインポートし、 入力パラメータを設定し、計算をします。

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

parms = [{ 
         'LATTICE'                   : "chain lattice", 
         'MODEL'                     : "spin",
         'CONSERVED_QUANTUMNUMBERS'  : 'Sz',
         'local_S'                   : 1,
         'J'                         : 1,
         'L'                         : 8
       }]

input_file = pyalps.writeInputFiles('parm6a',parms)
res = pyalps.runApplication('fulldiag',input_file)

次に、全出力ファイルに対して解析プログラムを実行します。

data = pyalps.evaluateFulldiagVersusT(pyalps.getResultFiles(prefix='parm6a'),DELTA_T=0.1,\
                                      T_MIN=0.1, T_MAX=10.0)

最後に、すべての物理量を図示します。

for s in pyalps.flatten(data):
  plt.figure()
  plt.title("Antiferromagnetic Heisenberg chain")
  pyalps.pyplot.plot(s)
plt.show()

Vistrailsを用いた実行

Vistrailsを用いた計算の実行のチュートリアルed-06-fulldiag.vtを開き、"antiferromagnetic chain"のワークフローを見てください。"Execute"をクリックすると、入力の設定、計算、結果の出力を行います。

演習

  • L=9(計算には数分かかりますのでL=7でも構いません)として磁化率の計算を行い、L=8のときの結果と比較してください。コマンドラインツールで解析する場合、parameter2xml を行うのを忘れないようにしてください。
  • 有限サイズの計算を無限サイズの計算と直接みなせるような温度範囲を見積もってみてください。

spin-1/2ハイゼンベルグ梯子

次に、長さ6のS=1/2ハイゼンベルグ梯子モデルの計算をおこないます。LATTICEパラメータを"ladder"に、J0,J1を1に変更します。パラメターファイルはparm6bです。Pythonスクリプトはtutorial6b.py、Vistrailsのワークフローはed-06-fulldiag.vtをそれぞれ参照してください。

演習

  • 比熱のピークの位置を議論してください(無限梯子はおおよそJ/2のギャップがあります。)
  • L=7(時間がかかるようなら、L=5でも可)の計算をおこなってみてください。
  • 有限サイズ計算がそのまま無限系の結果の良い近似となるような温度範囲を見積もってください。この範囲をS=1 ハイゼンベルグ鎖のそれと比べた時、何が言えるでしょうか?

注意: 量子モンテカルロ計算(QMC)に精通している方なら、QMC ではより大きな系での計算が行えて、熱力学極限を得るためにより良い近似が得られることがわかるでしょう。 (しかし実際のところ今回の系では、厳密対角化でも十分な結果を与えていて、これを超える計算をすることはそんなに簡単ではありません。実際に計算して確かめてみてください!)

特定の条件下では厳密対角化が最もよい計算手法となります。まず、系が本質的に有限で(可能な限り小さい)あれば、完全対角化で厳密解を得られるでしょう。また、厳密対角化には負符号問題が現れないので、フェルミオン系やフラストレートスピン系でも、制限なく計算できます。どちらの条件も同時に満たすのが次章で述べる磁性分子モデルです。

磁性分子の熱力学計算

次に、fulldiag を用いて小さなスピンクラスターの完全なスペクトルを計算し、厳密な熱力学シミュレーションをおこないます。

二個のカップリングした二量体

格子(グラフ)とモデルのカスタム設定

ここでは2つのカップリングした二量体を考えます。

最初に、この問題を表現するグラフが必要です。次のファイルdd-graph.xmlを参照してください

<GRAPH name="double dimer" vertices="4">
 <VERTEX id="1" type="0"></VERTEX>
 <VERTEX id="2" type="0"></VERTEX>
 <VERTEX id="3" type="1"></VERTEX>
 <VERTEX id="4" type="1"></VERTEX>
 <EDGE type="0" source="1" target="2"/>
 <EDGE type="0" source="3" target="4"/>
 <EDGE type="1" source="1" target="3"/>
 <EDGE type="1" source="1" target="4"/>
 <EDGE type="1" source="2" target="3"/>
 <EDGE type="1" source="2" target="4"/>
</GRAPH>

また、2種類のボンド0 と1 とに、異なるハイゼンベルグ交換相互作用 J0 とJ1 とを指定する必要があります。次のファイルmodel-dspin.xmlを参照してください。

<HAMILTONIAN name="dimerized spin">
 <PARAMETER name="J" default="1"/>
 <PARAMETER name="h" default="0"/>
 <BASIS ref="spin"/>
 <SITETERM site="i">
 <PARAMETER name="h#" default="h"/>
    -h#*Sz(i)
 </SITETERM>
 <BONDTERM source="i" target="j">
 <PARAMETER name="J#" default="J"/>
    J#*Sz(i)*Sz(j)+J#/2*(Splus(i)*Sminus(j)+Sminus(i)*Splus(j))
 </BONDTERM>
</HAMILTONIAN>

デフォルトで使われるmodels.xml ファイルにある'Spin' という名前のモデルには既に適切な​​定義が含まれていますので、実際に自分で定義する必要としないことに注意してください。 しかしながら上記の例は、ハッシュ記号(#)を使うことで、 n 番目のボンドにの自動的に相互作用 Jn を割り当てられるということを示す教育的な例となっています。

また、バーテックスの1番,2番と3番,4番に異なるタイプを割り当てました。これによって、それぞれlocal_S0とlocal_S1に対応する値を指定することにより、上側と下側の二量体に異なる大きさのローカルスピンを割り当てることができます。

このモデル磁化過程を求める計算をおこないます。計算条件は、ローカルスピンS0=1 (upper dimer)、S1=1/2 (lower dimer)、J0=1, J1=0.4、温度はT=0.02です。

コマンドラインを用いた実行

計算パラメータはparm6cを参照してください。

LATTICE="double dimer"
MODEL="dimerized spin"
LATTICE_LIBRARY="dd-graph.xml"
MODEL_LIBRARY="model-dspin.xml"
local_S0=1
local_S1=1/2
J0      = 1
J1      = 0.4
h       = 0 
CONSERVED_QUANTUMNUMBERS="Sz"
{T = 0.02}

新しいパラメーターLATTICE_LIBRARYMODEL_LIBRARYで 格子ファイルとモデルファイルの設定をおこないます。

次の手順にしたがって計算をおこないます。

parameter2xml parm3
fulldiag parm3.in.xml
fulldiag_evaluate --H_MIN 0 --H_MAX 4 --DELTA_H 0.025 --versus h parm3.task1.out.xml

ここではfulldiag_evaluate のコマンドライン引数として磁場の範囲を指定していることに注意してください。特に、コマンドライン引数 --versus h は、x軸として温度ではなく磁場を設定します。

Python,Vistrailsを用いた実行

Python、Vistrailsでの実行でも、格子ファイルとモデルファイルを指定する以外は以前と同じように行えます。tutorialcb.pyed-06-fulldiag.vtを参照してください。

Vistrailsでは2つの手法があります。一つは、Vistrailsでlattice,modelファイルを生成する方法、もう一つは、webからファイルをダウンロードする方法です。

問題

  • 結果をプロットし、解釈してください!

ヒント:カップルしたS0=1 ダイマーとS1=1/2ダイマーのスペクトルは解析的にわかります。エネルギーは(縮退を無視して)次のようになります。

  • -11J0/4, -3J0/4-2J1 for total spin Stotal=0;
  • -7J0/4, -3J0/4-J1, 5J0/4-3J1 for Stotal=1;
  • -3J0/4+J1, J0/4, 5J0/4-J1 for Stotal=2;
  • 5J0/4+2J1 for Stotal=3.

分子錯体 V15

最後の例は、V15分子錯体です。15個の1/2スピンで構成されています。

V15.png

グラフは v15-graph.xmlで定義されます。

今回のシミュレーションでは、相互作用 J はすべてのボンドで等しいと仮定します。パラメータファイルはparm6d,Pythonスクリプトはtutorial6d.py、Vistrailsのワークフローはed-06-fulldiag.vtを参照してください。


問題

  • 磁化率の低温度での振る舞いが説明できますか?

この計算は少しが難易度が高めで、計算にかかる時間も(休憩が取れるぐらいに)大きくなります。

追加演習

正方格子上のハバード模型

  • 正方格子上のハバード模型を計算するためのパラメータファイルを設定しましょう。
    • あなたが使いたい格子を探しましょう。正方格子ではなくて一次元鎖でも構いません。境界条件に注意してください。
    • model.xml からハバード模型を探しましょう。それぞれの項の意味をよく確かめてください。
    • この計算でも対称性を利用できます。このモデルで保存されるのは運動量(x、y方向それぞれ)と粒子数です。
    • t=0やU=0といったパラメータを用いて試行計算してみてください。
    • 次近接相互作用 t' を導入するにはどのようにすればよいでしょうか?

© 2004-2010 by Andreas Honecker, Emanuel Gull, and Matthias Troyer