コンテンツにスキップ
チュートリアル
Exact Diagonalization
ED-06 Full Diagonalization

ED-06 Full Diagonalization

このチュートリアルでは、fulldiag アプリケーションを使って小さな量子スピン系の完全なスペクトルを計算し、それを用いてすべての固有状態にわたる重み付き和として、量子モンテカルロの統計ノイズや DMRG の切断誤差を一切伴わない厳密な有限温度熱力学量——エネルギー、エントロピー、比熱、磁化率——を計算します。これは、ハミルトニアンの最低数個の固有状態だけを対象とする sparsediag を用いたこれまでのチュートリアルを補完するものです。ここではスペクトル全体が必要です。なぜなら、有限温度ではどれほどエネルギーの高い状態でもボルツマン重み e−En/Te^{-E_n/T} を担っているからです。

一次元スピン模型の熱力学

スピン1ハイゼンベルク鎖

最初のシミュレーションはスピン1ハイゼンベルク量子反強磁性鎖について行います。

H=J∑⟨i,j⟩Si⋅Sj,J=1, H = J \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j , \qquad J=1 ,

周期鎖上で定義され、その有限鎖の熱力学は最初に J.C. Bonner and M.E. Fisher, Phys. Rev. 135, A640 (1964) によって体系的に研究されました。

パラメータ

パラメータ意味値
LATTICE周期鎖chain lattice
MODEL量子スピン模型spin
local_S各サイトのスピン量子数1
J最近接結合1
L鎖長8
CONSERVED_QUANTUMNUMBERSHH をブロック対角化するために用いる対称性Sz

格子

このチュートリアルシリーズを通して使用してきたのと同じ周期的な chain lattice（ALPS格子ライブラリ参照）：

    J     J     J     J     J     J     J
o-------o-------o-------o-------o-------o-------o-------o
0       1       2       3       4       5       6       7
|_______________________________________________________|
                      J   (bond 7-0, periodic)

手法

熱力学量を求めるには、基底状態だけでなくスペクトル全体にわたってボルツマン重みを足し合わせる必要があるため、sparsediag の代わりに fulldiag を使用します。この S=1、L=8 鎖の完全なヒルベルト空間は次元 38=65613^8=6561 を持ちます。CONSERVED_QUANTUMNUMBERS=Sz を使うと、これがずっと小さな17個のブロック（最大のもの、すなわち Sz=0S_z=0 セクターでも高々 11071107 状態）に分解され、それぞれを直接対角化できます。これはより高速であるだけでなく、任意の温度で厳密な分配関数を再構成するのに十分です。

コマンドラインを用いる方法

パラメータファイル parm6a は、8サイトの一次元鎖上のスピン1ハイゼンベルク量子反強磁性体について全対角化を設定します。

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

CONSERVED_QUANTUMNUMBERS パラメータは、fulldiag がヒルベルト空間を不変部分空間に分解し、それぞれを個別に対角化するのに役立ちます。 標準的な一連のコマンドを使い、まず入力パラメータを XML に変換し、続いて fulldiag を使ってこの量子ハミルトニアンの完全なスペクトルを計算します。

parameter2xml parm6a
fulldiag --write-xml parm6a.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 ツールを使い、そのデータをお好みのプロットツールで表示できます。例えば、エネルギー密度対温度のデータを抽出するには次のようにします。

plot2text parm6a.task1.plot.energy.xml

同様の方法で、他の XML プロットファイルからもデータを抽出できます。 Grace がお好みのプロットツールであれば、plot2xmgr ツールを使って XML プロットファイルから直接 Grace のプロジェクトファイルを生成することもできます。例えば、エネルギー対温度の Grace プロジェクトファイルを生成するには次のようにします。

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

同様に、ツール plot2gp は Gnuplot スクリプトを生成し、plot2text はファイルをプレーンテキストに変換します。ただし、データ処理とプロットには Python を使う方法が推奨されます。

Python を用いる方法

Python でシミュレーションを設定・実行するには、スクリプト tutorial6a.py を使用します。このスクリプトの最初の部分では必要なモジュールをインポートし、続いて入力パラメータを Python の辞書のリストとして準備し、シミュレーションを実行します。

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

出力データ

同じハミルトニアンを独立に対角化すると、次のようなサイトあたりの熱力学量（エネルギー E/LE/L、比熱 C/LC/L、均一磁化率 χ/L\chi/L）が得られます。

T/JT/JE/LE/LC/LC/Lχ/L\chi/L
0.1-1.4170.0340.0066
0.2-1.4070.1380.0567
0.5-1.3330.4010.1276
1.0-1.0640.5590.1693
2.0-0.6530.2780.1639
5.0-0.2740.0550.1010

比熱は T≈JT\approx J 付近に幅広い極大を持ち、磁化率は T≈1T\approx1–1.5 J1.5\,J 付近でピークとなった後、T=0T=0 に向かって急激に低下します——これは、ED-02 で計算した Haldane ギャップ Δ≈0.41 J\Delta\approx0.41\,J から期待される指数関数的な活性化挙動 χ∼e−Δ/T\chi\sim e^{-\Delta/T} の有限サイズにおける名残です。

演習

  • 磁化率の計算を、今度は L=9 について繰り返してみてください（数分かかります。急ぐ場合は L=7 を使ってください）。そして L=8 の結果と比較してください！parameter2xml から始まるすべてのステップを繰り返す必要があることを思い出してください。

  • 無限鎖に対する近似として有限サイズの結果を使ってよい温度範囲を大まかに見積もってください。

スピン1/2ハイゼンベルクラダー

ここで模型を二本脚ラダーに変更します。

H=J0∑legsSi⋅Sj+J1∑rungsSi⋅Sj, H = J_0 \sum_{\text{legs}} \mathbf{S}_i \cdot \mathbf{S}_j + J_1 \sum_{\text{rungs}} \mathbf{S}_i \cdot \mathbf{S}_j ,

長さ L=6（つまり6ラング、12サイト）とし、パラメータ LATTICE を “ladder” に、結合 J0 と J1 を1に設定します——これは ED-03 で使用したのと同じ格子です。

      o---J0---o---J0---o---J0---o---J0---o---J0---o     leg 0
      |        |        |        |        |        |
      J1       J1       J1       J1       J1       J1
      |        |        |        |        |        |
      o---J0---o---J0---o---J0---o---J0---o---J0---o     leg 1

パラメータは parm6b、Python スクリプトは tutorial6b.py にあります。上の parm6a と全く同じように実行します。

parameter2xml parm6b
fulldiag --write-xml parm6b.in.xml
fulldiag_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T 0.1 parm6b.task1.out.xml

parm6b は、鎖のパラメータの代わりに LATTICE="ladder"、J0=1、J1=1、{L = 6} を使っている以外は parm6a と全く同じです。

出力データ

この L=6 ラダー（脚は開放境界条件）を独立に対角化すると、有限サイズの一重項-三重項ギャップは Δ≈0.68 J\Delta\approx0.68\,J となり、これは下のヒントで引用されている無限ラダーの受け入れられている値 ≈0.5 J\approx 0.5\,J よりもやや大きくなります——これは、6ラングの系がギャップを持つ状態については熱力学極限からまだかなり離れていることを思い出させてくれます。対応するサイトあたりの熱力学量は次の通りです。

T/JT/JE/LE/LC/LC/Lχ/L\chi/L
0.1-0.5500.0130.0019
0.2-0.5440.1260.0294
0.5-0.4510.4240.1065
0.7-0.3680.3820.1204
1.0-0.2740.2520.1189
2.0-0.1370.0720.0876

演習

  • 比熱の最大値の位置について議論してください（ヒント：無限ラダーのギャップは約 J/2 です）！
  • この計算を7ラング（急ぐ場合は5ラング）について繰り返してください。
  • 有限サイズの結果が無限系に対する良い近似となる温度範囲について、この比較から何が推測できますか。

補足： 量子モンテカルロ（QMC）シミュレーションに詳しい方であれば、より大きな系を扱うことができ、したがって熱力学極限に対するより良い近似が得られることをご存知でしょう。（試してみてください！上記の例については、QMC で全対角化より良い結果を出すことは実はそれほど簡単ではないことに気づくかもしれません。）

一定の条件下では、厳密対角化は依然として最良の手法であり続けます。第一に、系が本質的に有限（かつおそらく小さい）場合、全対角化は即座に厳密な結果を与えます。第二に、厳密対角化は符号問題に悩まされることがないため、追加の制約なしにフェルミオン模型やフラストレート模型にも使用できます。この両方の条件は、第2部で議論する磁性分子の模型において同時に満たされています。

磁性分子の熱力学

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

結合した2つの二量体

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

次のような、結合した2つの二量体からなる系を考えます：J0J_0 で結合された2つのスピン S0S_0 からなる「上側」二量体、同じく J0J_0 で結合された2つのスピン S1S_1 からなる「下側」二量体、そして上側二量体のすべてのスピンが下側二量体のすべてのスピンと J1J_1 で結合されているものです。

H=J0(S1⋅S2+S3⋅S4)+J1(S1⋅S3+S1⋅S4+S2⋅S3+S2⋅S4)−h∑i=14Siz. H = J_0\left(\mathbf{S}_1\cdot\mathbf{S}_2 + \mathbf{S}_3\cdot\mathbf{S}_4\right) + J_1\left(\mathbf{S}_1\cdot\mathbf{S}_3+\mathbf{S}_1\cdot\mathbf{S}_4+\mathbf{S}_2\cdot\mathbf{S}_3+\mathbf{S}_2\cdot\mathbf{S}_4\right) - h\sum_{i=1}^4 S^z_i . 
        J0
   1 ---------- 2        (upper dimer, spin S0)
   |  \        /  |
   |   \  J1  /   |
   J1    \  /    J1
   |     /  \     |
   |   /      \   |
   3 ---------- 4        (lower dimer, spin S1)
        J0

これはまさに以下でコード化されているグラフです：サイト1、2（タイプ0）が上側二量体を、サイト3、4（タイプ1）が下側二量体を構成し、タイプ0のエッジが二量体内部のボンド（J0J_0）、タイプ1のエッジが4本の二量体間ボンド（J1J_1）です。

まず、この問題を表すグラフが必要です。これはファイル dd-graph.xml の以下のエントリで定義されます。

<LATTICES>
<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>
</LATTICES>

注：このファイルは <?xml?> 宣言や <!DOCTYPE> ヘッダーを含まない素の XML でなければなりません——ALPS のパーサーはこれらをサポートしていません。

次に、タイプ0とタイプ1のエッジにそれぞれ異なるハイゼンベルク交換 J0 と J1 を割り当てたいと思います。これはファイル model-dspin.xml の以下のエントリによって実現されます。

<MODELS>
<SITEBASIS name="spin">
  <PARAMETER name="local_spin" default="local_S"/>
  <PARAMETER name="local_S" default="1/2"/>
  <QUANTUMNUMBER name="S" min="local_spin" max="local_spin"/>
  <QUANTUMNUMBER name="Sz" min="-S" max="S"/>
  <OPERATOR name="Splus" matrixelement="sqrt(S*(S+1)-Sz*(Sz+1))">
    <CHANGE quantumnumber="Sz" change="1"/>
  </OPERATOR>
  <OPERATOR name="Sminus" matrixelement="sqrt(S*(S+1)-Sz*(Sz-1))">
    <CHANGE quantumnumber="Sz" change="-1"/>
  </OPERATOR>
  <OPERATOR name="Sz" matrixelement="Sz"/>
</SITEBASIS>
<BASIS name="spin">
  <SITEBASIS ref="spin">
    <PARAMETER name="local_spin" value="local_S#"/>
    <PARAMETER name="local_S#" value="1/2"/>
  </SITEBASIS>
  <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
</BASIS>
<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?> 宣言や <!DOCTYPE> ヘッダーを含まない素の XML でなければならず、<SITEBASIS> と <BASIS> の定義を含まなければなりません——これらがないと fulldiag はヒルベルト空間を構築できません。

なお、デフォルトの models.xml ファイル中の “spin” ハミルトニアンにすでに適切な定義が含まれているため、実際にはこの定義は必要ありません。とはいえ、上の例は、ハッシュ記号（#）を使ってタイプ n のボンドに交換定数 Jn を自動的に割り当てる方法を示しています。

ついでに、サイト1、2 と 3、4 にも異なるタイプを割り当てました。したがって、local_S0 と local_S1 の値をそれぞれ指定することで、上側二量体と下側二量体に異なる局所スピンを割り当てることができます。

ここでは、局所スピン S0=1S_0=1（上側二量体）と S1=1/2S_1=1/2（下側二量体）、J0=1J_0=1、J1=0.4J_1=0.4、温度 T=0.02T=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_LIBRARY と MODEL_LIBRARY に注目してください。これらはカスタムの格子とモデルを含むファイルを指しています。 この計算は次のコマンド列で実行されます。

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

ここでは磁場の範囲を指定するために、コマンドライン引数付きの fulldiag_evaluate を使用していることに注意してください。特にコマンドライン引数 –versus h は、これまでの例のような温度ではなく、磁場を x 軸に置きます。

手法

わずか (2⋅1+1)2×(2⋅12+1)2=36(2\cdot1+1)^2\times(2\cdot\tfrac12+1)^2=36 個の基底状態しかないため、この系は fulldiag にとって取るに足らないものです——この例の要点は計算の難しさではなく、実際の磁性分子データに必要とされるような、完全にカスタムな格子グラフとハミルトニアンを、パラメータファイルだけで組み立てる方法を示すことにあります。

出力データ

ここで使用する低温 T=0.02 J0T=0.02\,J_0 では、hh を上げるにつれて基底状態の全スピン StotalS_{\text{total}} が単位段階で増加し、各段階は以下の解析スペクトルから2つの準位が交差する点で起こります。

h/J0h/J_0⟨Stotalz⟩\langle S^z_{\text{total}}\rangle基底状態のセクター
0.0 – 0.80.00Stotal=0S_{\text{total}}=0
1.00.67交差（Stotal=0→1S_{\text{total}}=0\to1）
1.21.00Stotal=1S_{\text{total}}=1
1.51.99交差（Stotal=1→2S_{\text{total}}=1\to2）
1.8 – 2.22.00Stotal=2S_{\text{total}}=2
2.52.99交差（Stotal=2→3S_{\text{total}}=2\to3）
3.0 – 4.03.00Stotal=3S_{\text{total}}=3

これら3つの交差磁場は、下のヒントで与えられている解析的エネルギーから直接導かれます：Stotal=nS_{\text{total}}=n（Sz=nS^z=n における）の最低エネルギーと Stotal=n+1S_{\text{total}}=n+1（Sz=n+1S^z=n+1 における）の最低エネルギーを等しいと置くと、J0J_0 を単位として h1=1.0h_1=1.0、h2=1.4h_2=1.4、h3=2.4h_3=2.4 が得られます——これは Stotalz=0,1,2,3S^z_{\text{total}}=0,1,2,3 にプラトーを持つ磁化階段であり、小さいながらもゼロでない温度によってわずかに丸められているだけです。

問題

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

ヒント： 結合した S0=1S_0=1 と S1=1/2S_1=1/2 の2つの二量体のスペクトルは解析的に求めることができます。エネルギーは（いくつかの縮退を伴って）次の通りです。

  • 全スピン Stotal=0S_{\text{total}}=0：−11J0/4-11J_0/4、−3J0/4−2J1-3J_0/4-2J_1；
  • Stotal=1S_{\text{total}}=1：−7J0/4-7J_0/4、−3J0/4−J1-3J_0/4-J_1、5J0/4−3J15J_0/4-3J_1；
  • Stotal=2S_{\text{total}}=2：−3J0/4+J1-3J_0/4+J_1、J0/4J_0/4、5J0/4−J15J_0/4-J_1；
  • Stotal=3S_{\text{total}}=3：5J0/4+2J15J_0/4+2J_1。

分子錯体 V15V_{15}

最後の例は、分子ナノ磁石 V15\mathrm{V}_{15} のモデルです。これは15個の V4+\mathrm{V}^{4+} イオン（それぞれ S=1/2S=1/2）が、中心の三角形を挟む2つの六角形上に配置されたものです。その低エネルギー磁性は、2つの六角形周りの強い反強磁性交換と、それらを中心の三角形に結ぶより弱い交換によって支配されており、そのため低温では15個のスピンが実効的にフラストレートしたスピン1/2三角形のように振る舞います：スペクトルの残りの部分からよく分離した、2つの低エネルギー二重項と1つの四重項です（G. Chaboussant et al., Europhys. Lett. 59, 291 (2002)）。

        1---2                 outer hexagon (6 spins)
       /     \
      6       3
       \     /
        5---4
       /     \
      o       o                inner hexagon (6 spins),
       \     /                 rotated relative to the outer one
        o---o
          |
        triangle (3 spins)    weakly coupled to both hexagons

シミュレーションでは、このグラフのすべての辺に沿って等しいハイゼンベルク交換 JJ を持つという簡単化された仮定を置きます。カスタム格子とモデルファイルの設定は、上の結合二量体の例と全く同じ手順に従います——15個の頂点を持つ <GRAPH> と、六角形および三角形のボンドに対応する <EDGE> リスト、そして同じ dimerized spin（あるいは内蔵の spin）ハミルトニアンです。関連するグラフは v15-graph.xml で定義されています。シミュレーションの実行は上と同様に進みます。パラメータは parm6d、Python スクリプトは tutorial6d.py にあります。

手法

15個のスピン1/2サイトの完全なヒルベルト空間は次元 215=327682^{15}=32768 を持ちます。SzS_z セクターに分解した後でも、fulldiag は低温熱力学を正しく求めるために数千状態のブロックを厳密に対角化しなければなりません。これが、下の注記にある通り、この実行がこれまでの例よりも明らかに長い時間を要する理由です。

問題

  • 磁気磁化率の低温での振る舞いをどのように説明しますか。（ヒント：上で述べた実効的な低エネルギースピン1/2三角形について考え、磁場中でのその準位交差が χ(T)\chi(T) にどのように現れるかを考えてみてください。）

この計算はすでにいくらか大変で、時間がかかることに注意してください。ですので、結果を見に戻ってくる前に、休憩を取るのも良い考えかもしれません。

まとめ

全対角化は、並進不変な鎖やラダーから、完全にカスタムでフラストレートした磁性分子に至るまで、あらゆる有限量子スピン系について数値的に厳密な熱力学を与えます——その代償はヒルベルト空間がシステムサイズとともに指数関数的に増大することであり、これが、対称性を利用してもこの手法が数十スピン程度に限られる理由です。対照的に、基底状態と低エネルギー性質だけであれば sparsediag によってはるかに大きな系に到達できます。

追加演習

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

  • 正方格子上のハバード模型のパラメータファイルを設定してください。
    • 使いたい格子を見つけてください。正方格子でも鎖でも構いません。境界条件に注意してください！
    • モデルライブラリの中からハバード模型を見つけてください。その各項を理解するようにしてください。
    • 対称性をオンにしてください：（x 方向と y 方向の）運動量と粒子数が保存量です。
    • 検証可能な試験的パラメータ（例えば t=0 や U=0）を選び、シミュレーションを実行してください。
    • t’ はどのように導入しますか。
ED-05 ED Phase Transition