ALPS 2 Tutorials:DMRG-01 DMRG/ja

Models: ハイゼンベルグスピン鎖

DMRGに関して、次のハミルトニアンで与えれる長さLの反強磁性ハイゼンベルグ鎖（スピン1/2,1）の２つのmodelを考えます。

$H = \sum_{i=1}^{L-1} \frac{1}{2} (S^+_i S^-_{i+1} + S^-_i S^+_{i+1}) + S^z_i S^z_{i+1}$ .

既に多くのチュートリアルで取り上げたこの２つのmodelsを選択した理由は、これらのモデルは類似性があるにもかかわらず、物理的な振る舞いがまったく異なり、DMRGアルゴリズムを学ぶうえで適していると考えられます。以下にこれらの物理特性を述べていきます。

スピン-1/2鎖

スピン1/2鎖の基底状態はBethe仮説によって完全に記述されます。すなわち、基底状態エネルギーは正確に知ることができます。熱力学的極限$L\rightarrow\infty$のサイトあたりのエネルギーは次式で与えられます。

$E_0 = 1/4 - \ln 2 = -0,4431471805599...$

基底状態のエネルギーは他の状態と比較しないと、それ自体はあまり意味の無いものですが、DMRG法のベンチマークとしては非常に有益な情報が得られます。基底状態は、熱力学的極限でのエネルギー差によって生じる励起状態とどれくらい異なるか、例えば、gapは消失するのかしないのか、といった情報に注目します。スピン1/2鎖では、gapは0です。

同時に、異なるサイト上での相関はどのようになっているか？無限長スピン1/2鎖は漸近的(i.e. for $|i-j| \rightarrow \infty$)に次のようになります。

$\langle S^z_i S^z_j \rangle \sim (-1)^{|i-j|} \frac{\sqrt{\ln|i-j|}}{|i-j|}$ .

これらのことは、スピン1/2鎖はcriticalを意味します。スピン間の反強磁性相関はそのpower lawに従い、距離とともに減少します。この場合、power lawの指数は$-1$になります。長鎖上でのDMRG計算によって求められるこの　対数補正用に平方根が追加されています。これは長鎖上でのDMRG計算によって確かめられます。この変数を用いたアルゴリズムは非常にゆっくりですが、値が増加してしまいますが、最初の段階では無視できるレベルです。

スピン-1鎖

長年、異なるスピン長に起因する定量的な差異に関して、スピン-1鎖は同じように振る舞うであろうと考えられてきました。１９８２年Duncan Haldaneによって根本的な差異があると報告されました。それは、等方性反強磁性Heisenberg鎖はスピン長に依存、すなわち半整数($S=1/2,3/2,...$)と整数($S=1$)で異なる振る舞いがあることが明らかになりました。それ以降、スピン-1鎖への注目が高まり、さまざまな研究がおこなわれています。

スピン-1/2鎖とは異なり、スピン-1鎖は解析に完全に計算できる特性はありません。定量的な表現になるように、完全に数値的に対応しなくてはいけません。

定常状態のエネルギーは次式で与えられます。

$E_0 = -1.401484039 ...$.

また、ギャップの存在がスピン-1/2とは異なり非常に重要な問題になってきます。熱力学的極限において、スピン-1のギャップは有限であり、次式で与えられます。

$\Delta = 0.41052$

5桁の精度です。

スピン-スピン相関の振る舞いの問題も、スピン-1/2のケースとは大きく異なります。相関は漸近的に求められます。($|i-j| \rightarrow \infty$)

$\langle S^z_i S^z_j \rangle \sim (-1)^{|i-j|} \frac{\exp (-|i-j|/\xi)}{\sqrt{|i-j|}}$ .

支配的な寄与は、長さのスケール$\xi$でおこる指数関数的な減衰です。この特徴的なケースのcorrelation lengthは数値的に$\xi=6.02$と知られています。分母の距離の平方根によって解析的な補正があります。しかし、指数関数的減衰にくらべ影響が小さいため、この補正長の計算は無視されます。もちろん補正長が大きくなれば、問題になってきます。

スピン-1鎖は有限ギャップのnon-critical量子系、指数関数的に減衰する相関のサンプル計算としてとりあげます。

チュートリアルの進め方

チュートリアルを通して学びたいことは、ALPS DMRGを用いて上記のすべての数量を求める計算方法です。あわせて、間違いやすい初歩的な問題点についても学びます。

Vive la difference ...

他の数値法と比べて最も重要な差異はDMRGは他の数値解法や厳密解は閉境界であるのに対して開境界条件を選ぶということです。サイト１と$L$に２つの鎖端が存在します。

コードの実行

概論

まず最初に、 DMRGアルゴリズムについて簡単に説明します。スピン長$S$、$d=2S+1$であらわされる次元$d$の局所空間を持つ一次元量子系が与えられると、Hilbert空間の次元はシステムサイズ$L$の指数的$d^L$に増大します。厳密対角化により、小さい系では、指数関数的に増大するHilbert空間で正しい値がえられます。量子モンテカルロ計算は、確率的に大空間をサンプリングすることで近似値が得られます。DMRGは別の方法でアプローチし、指数関数的に増大するHilbert空間の非常に小さな部分空間を確認します。

最初の制御パラメータは$D$です。これは部分空間を番号づけしたものです。DMRGはこのパラメータに対して単調に比例します。部分空間が大きくなれば、近似値もより正確になります。制限もあります。もし$D\rightarrow d^L$であれば、状態は表現できず、正確に求めることができません。しかし、このような大きな状態が計算できるのであれば、厳密対角化は優れた代替手段となるでしょう。DMRGの計算では$D$は、matrix dimensionやnumber of block statesなど、様々な値を表現するので注意が必要です。

2番目の制御パラメータはシステムサイズ$L$です。3つ目はDMRGアルゴリズムを詳しく理解してからでないとよく分からないでしょう。最良の状態近似を求めるためにDMRGは次の２つのステップをおこないます。

1. 最初のステップ（いわゆる無限系DMRG）では、アルゴリズムは、目的のシステムサイズ$L$になるまで、2,4,6の鎖長の解析を反復することで、適した部分空間を求めようとします。この手順は、反復毎に鎖を分割し、中央に２つの新しいサイトを挿入します。この手順はもちろん無限に継続し、このことが名前の由来になっています。しかし、非常に有意義な結果を期待するのは危険です。2番目に注意する点は、この手順はDMRGを扱うための鎖長も支持するということです。
2. 2番目のステップ（いわゆる有限系DMRG）では、短鎖で部分空間は、すべての量子ゆらぎや相関を考慮できないことを取扱います。何をするかというと、部分空間の質を向上させるため、さらなる反復計算が必要になってきます。DMRGで鎖の全サイトで見られるこのような反復処理をsweepと呼びます。このsweepの回数が、3番目に重要な制御パラメータです。このパラメータは小さいと目的の$D$は達成できないし、長すぎると計算時間が膨大になってしまいます。

最後の注意として、DMRG実行によって得られた正確な指標truncation errorについて考察します。簡単に述べると、DMRGのアルゴリズムの各ポイントは状態空間の指数関数的な増大の方向に１ステップ進行し、もしそのステップが適切で無いなら、 その固有状態に対応する量（固有値）の分布に関する密度行列の分析を用いて、どれくらい正確かどうか求めます。DMRGの近似は、統計的な重みを除外されます。いわいる打ち切り誤差が発生します。多くのDMRGプログラムでは、様々な工夫によってこの誤差は$10^{-12}$程度になっており、微々たるレベルです。このチュートリアルでは、ローカルな計算量（エネルギー、磁化など）の誤差はおおまかに打ち切り誤差に比例（通常それよりも少し大きい）することを理解することが重要です。

ALPS DMRGコードと制御パラメータ

ハミルトニアンや格子形状のようなインプットに加え、DMRGシミュレーションは特定のパラメータセットが必要になります。詳細を述べていきます。

DMRGパラメータ

NUMBER_EIGENVALUES
計算するための固有状態とエネルギーの数。デフォルトでは1です。、ギャップを計算するために2に設定する必要があります。

SWEEPS
有限サイズアルゴリズムのためにDMRGスイープ数。各掃引では、左から右へ半掃引、右から左半分掃引が含まれます。

NUM_WARMUP_STATES
DMRGのブロックを成長させる初期の状態の数。指定されていない場合、アルゴリズムはデフォルトとして２０状態を使用します。

STATES
DMRG状態の数は、各半分の掃引を続けます。ユーザーは、STATES、または1つMAXSTATESまたはNUMSTATES値の2*sweep異なる値のいずれかを指定する必要があります。

MAXSTATES

DMRG状態の最大数は保持されます。ユーザーは各ハーフ掃引、またはMAXSTATESまたはプログラムの基礎を育てるために使用するNUMSTATESためSTATESの値のいずれかを指定することもできます。プログラムは自動的にどのように多くの状態がMAXSTATESに達するまでSTATES/(2*SWEEPS)のステップで基本サイズを増加、各掃引のために使用することを決定します。

NUMSTATES
全掃引を保持するDMRGの状態の定数値

TRUNCATION_ERROR ユーザーは状態の数ではなく、シミュレーションのための範囲を設定することを選択できます。プログラムは、この範囲が条件を満たすためにどれくらいの数の状態が必要か自動的に判断します。この場合、サイズがシミュレーションで制御できない程、大きくなる可能性があるので注意が必要です。それは、前に説明したように、MAXSTATESまたはNUMSTATESのいずれかを使用して、制約としての状態の最大数を指定することが賢明です。

LANCZOS_TOLERANCE 対角化（Davidson/Lanczos）コードのための許容範囲。デフォルト値は10^ -7です。

CONSERVED_QUANTUMNUMBERS

量子数は関連のあるモデルにて保存されます。ブロック形式で行列を減らすためにコードで使用されます。特定の量子数が指定されていない場合、プログラムはグランドカノニカルで動作します。スピン鎖でSz_total宣言しない場合は、プログラムでは、DIM=2^Nの状態とHilbert空間を用いて実行されます。"canonical"で実行する（Sz_total=0を設定することにより）ことで、大幅に削減次元で空間で作業することによって、パフォーマンスが向上します。これを行う方法の例については、DMRGコードに含まれてPARMSファイルを見てみましょう。

適切なパラメータの選択

デフォルトの入力値を使用するようにお勧めできません。DMRGの収束は前計算、掃引数や保持する状態の最大数に強く影響を受けます。良い計算例は基底状態エネルギーや、状態の数の関数としての打ち切り誤差の収束です。これは特定のしきい値以下の誤差を維持するために必要な状態の数を示します。掃引数が十分であるかどうか判断するために、相関関係、またはそのようなスピンのSZの投影や、粒子密度などのローカル量の空間分布を観測することがお勧めです。 例えば、反転対称モデルでは、それらの観測値も対称であると予想できます。他に対称となるべき数値はエンタングルメント・エントロピーがあります。もしこの結果が反映されていなかったら、掃引が十分でない可能性が考えられます。ハミルトニアンがSzやNといった量子数を保持するのであれば、次元を減少させる部分空間のシミュレーションで実行するためこれらの値を修正することが可能です。計算時間の短縮や使用メモリの削減につながります。

基底状態エネルギー

定常状態$| \psi_0 \rangle$とそのエネルギー$E_0$について考えます。ここでは、2つのケースを区別する必要があります。まず、与えられた長さ$L$の鎖上の与えられたハミルトニアンの基底状態のエネルギーに興味があるかもしれません。第二に、我々は熱力学的極限におけるサイトあたりのエネルギー（またはボンドあたり）に興味があるかもしれません。

固定長基底状態エネルギー

長さ$L =32,64,96,128</ math>の鎖を考慮してください。スピン-1/2とスピン-1両方で、 <math>D=50,100,150,200,300$の状態数で定常状態のエネルギー計算の準備をおこないます。各長さについて、$D$の関数として打ち切り誤差と定常状態のエネルギーを集計してください。結果が正しく収束するように、掃引数を注意深く調整してみてください。

1. 各システム長やスピン長で、トータルエネルギーvs.打ち切り後のプロットを通して、トータルエネルギーと打ち切り誤差間の関係性を考えてみてください。
2. $D$でスピン-1/2、スピン-1でどのように収束するのか観測してください。

これらの２つのケースについて収束の違いがわかりますか？ "Hint:"観点は何でしょう？　大規模システムで、スピン-1では収束性に弱い影響しかありませんが、スピン-1/2では非常に強い依存性があります。 このことは、スピン-1鎖は相関長のセグメントによって影響を受けているのに対し、スピン-1/2鎖は臨界のため有限長のスケールがありません。

1. $D\rightarrow\infty$の極限での各鎖長の定常状態エネルギーを推定してみてください。

一次元ハイゼンベルグ鎖（S=1/2）　

シリアル計算

最初の例では、100の状態を維持し、32サイト、開境界条件とスピン-1/2ハイゼンベルク鎖のシミュレーションを設定しています。

コマンドラインでの実行

パラメータファイルはspin_one_halfです。

LATTICE="open chain lattice"
MODEL="spin" 
CONSERVED_QUANTUMNUMBERS="N,Sz" 
Sz_total=0
J=1
SWEEPS=4
NUMBER_EIGENVALUES=1
L=32 
{MAXSTATES=100}

次の用にコマンドし、dmrgアプリケーションを使用します。

 parameter2xml spin_one_half
 dmrg --write-xml spin_one_half.in.xml

全計算結果は出力ファイルspin_one_half.task1.out.xmlに格納されています。ブラウザー等で確認してください。

DMRGは４種類の掃引（左から右の４つの半分の掃引と右から左の４つの半分の掃引）を実行します。設定されたMAXSTATES=100に到達するまで、MAXSTATES/(2*SWEEPS)づつ増加していきます。

これらは収束のデフォルトオプションです。下記のスピンS=1の例に示すように、状態数をカスタマイズすることができます。

Pythonを使った実行

Pythonを使った実行には、spin_one_half.pyのスクリプトを参照してください。

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

parms = [ { 
       'LATTICE'                   : "open chain lattice", 
       'MODEL'                     : "spin",
       'CONSERVED_QUANTUMNUMBERS'  : 'N,Sz',
       'Sz_total'                  : 0,
       'J'                         : 1,
       'SWEEPS'                    : 4,
       'NUMBER_EIGENVALUES'        : 1,
       'L'                         : 32,
       'MAXSTATES'                 : 100
      } ]

input_file = pyalps.writeInputFiles('parm_spin_one_half',parms)
res = pyalps.runApplication('dmrg',input_file,writexml=True)

次に、DMRGによって計算された定常状態の特性をロードします。

data = pyalps.loadEigenstateMeasurements(pyalps.getResultFiles(prefix='parm_spin_one_half'))

プリントします。

for s in data[0]:
   print s.props['observable'], ' : ', s.y[0]

さらに、我々は、各反復ステップの詳細なデータを読み込むことができます

iter = pyalps.loadMeasurements(pyalps.getResultFiles(prefix='parm_spin_one_half'),
                               what=['Iteration Energy','Iteration Truncation Error'])

最後の計算結果からどのようにDMRGアルゴリズムが収束したか見ることができます。

plt.figure()
pyalps.pyplot.plot(iter[0][0])
plt.title('Iteration history of ground state energy (S=1/2)')
plt.ylim(-15,0)
plt.ylabel('$E_0$')
plt.xlabel('iteration')

plt.figure()
pyalps.plot.plot(iter[0][1])
plt.title('Iteration history of truncation error (S=1/2)')
plt.yscale('log')
plt.ylabel('error')
plt.xlabel('iteration')

plt.show()

Vistrailsを使った実行

Vistrailsを用いた計算はdmrg-01-dmrg.vtを使用します。"spin 1/2 iterations"を実行してください。

並列計算

コマンドラインでの実行

コマンドラインでの実行にはspin_one_half_multipleファイルを参照してください。

LATTICE="open chain lattice"
CONSERVED_QUANTUMNUMBERS="N,Sz"
MODEL="spin"
Sz_total=0
J=1
NUMBER_EIGENVALUES=1
SWEEPS=4
L=32
{ MAXSTATES=20 }
{ MAXSTATES=40 }
{ MAXSTATES=60 }

見て明らかなように、前の例との主な違いは、括弧内にエンコードされたパラメータで構成されています。 前と同じように、実行します。

 parameter2xml spin_one_half_multiple
 dmrg --write-xml spin_one_half_multiple.in.xml

このケースでは、アウトプットとして、spin_one_half_multiple.task#.out.xmlが生成されます。

Pythonを使った実行

Pythonを使った実行には、spin_one_half_multiple.pyのスクリプトを参照してください。

parms= []
for m in [20,40,60]:
   parms.append({ 
       'LATTICE'                   : "open chain lattice", 
       'MODEL'                     : "spin",
       'CONSERVED_QUANTUMNUMBERS'  : 'N,Sz',
       'Sz_total'                  : 0,
       'J'                         : 1,
       'SWEEPS'                    : 4,
       'NUMBER_EIGENVALUES'        : 1,
       'L'                         : 32,
       'MAXSTATES'                 : m
      })

パラメータファイルを作成後、dmrgアプリケーションを実行し、シングル実行と同じように計算結果をロードしてください。全計算結果から物理量のプロットが可能です。

for run in data:
    for s in run:
        print s.props['observable'], ' : ', s.y[0]

Vistrailsを使った実行

Vistrailsを用いた計算はdmrg-01-dmrg.vtをしようします。"spin 1/2 multiple"を実行してください。

一次元ハイゼンベルグ鎖（S=1）　

S = 1ハイゼンベルグ鎖が開境界条件に起因するいくつかの特別な対処が必要になります。上記で説明したように、スピンS=1/2、それらのそれぞれの持つ鎖の両端に2つのサイトを含める必要があります。また、新しい格子ファイルを定義しなくてはなりません。結局のところ、これらを簡単におこなう方法はありませんので、手動でおこなう必要があります。手順を簡単にするため、格子を生成するシンプルなPythonスクリプトbuild_lattice.pyを準備します。唯一の入力値は、格子のサイト数です。例えば、次の用にタイプします。

$ alpspython build_lattice.py 6

出力は次のようになります。

<LATTICES>
 <GRAPH name = "open chain lattice with special edges" dimension="1" vertices="6" edges="5">
  <VERTEX id="1" type="0"><COORDINATE>0</COORDINATE></VERTEX>
  <VERTEX id="2" type="1"><COORDINATE>2</COORDINATE></VERTEX>
  <VERTEX id="3" type="1"><COORDINATE>3</COORDINATE></VERTEX>
  <VERTEX id="4" type="1"><COORDINATE>4</COORDINATE></VERTEX>
  <VERTEX id="5" type="1"><COORDINATE>5</COORDINATE></VERTEX>
  <VERTEX id="6" type="0"><COORDINATE>6</COORDINATE></VERTEX>
  <EDGE source="1" target="2" id="1" type="0" vector="1"/>
  <EDGE source="2" target="3" id="2" type="0" vector="1"/>
  <EDGE source="3" target="4" id="3" type="0" vector="1"/>
  <EDGE source="4" target="5" id="4" type="0" vector="1"/>
  <EDGE source="5" target="6" id="5" type="0" vector="1"/>
 </GRAPH>
</LATTICES>

明らかなように、格子は、６つのバーテックスとその近傍をつなぐ辺を含む1次元グラフで定義されます。最初と最後のバーテックスはタイプ"0"でその他は"1"です。この格子の上にモデルを実装するために、この定義を使用しなければいけません。これらのバーテックスに様々な情報が設定されている必要があります。これらの設定をおこなうパラメータの指定方法は、次の通りです。

local_S0=0.5
local_S1=1

32サイト格子での実行です。

$ alpspython build_lattice.py 32 > my_lattice.xml

コマンドラインでの実行

最後です。こちらのパラメータファイルspin_oneを参照してください。

LATTICE_LIBRARY="my_lattice.xml"
LATTICE="open chain lattice with special edges"
MODEL="spin"
local_S0=0.5
local_S1=1
CONSERVED_QUANTUMNUMBERS="N,Sz"
Sz_total=0
J=1
SWEEPS=4
NUMBER_EIGENVALUES=1
MAXSTATES=100

明らかに、それは各システムのサイズに対してこのプロセスを繰り返すことが面倒なことになります。簡単な解決方法は自動的に実行するスクリプトを生成することです。単純な方法は、必要な格子ライブラリを定義することです。 my_lattices.xmlファイルを参照してください。$L=32,64,96,128,192$の格子情報があります。

LATTICE_LIBRARY="my_lattices.xml"
LATTICE="open chain lattice with special edges 32"

名前に格子サイズを入れています。

Pythonでの実行

Pythonを使った実行には、spin_one.pyのスクリプトを参照してください。

parms = [ { 
       'LATTICE_LIBRARY'           : 'my_lattice.xml',
       'LATTICE'                   : 'open chain lattice with special edges',
       'MODEL'                     : 'spin',
       'local_S0'                  : '0.5',
       'local_S1'                  : '1',
       'CONSERVED_QUANTUMNUMBERS'  : 'N,Sz',
       'Sz_total'                  : 0,
       'J'                         : 1,
       'SWEEPS'                    : 4,
       'NUMBER_EIGENVALUES'        : 1,
       'MAXSTATES'                 : 100
      } ]

spin_one_half.pyも参照してください。スクリプトは上述のものと同様です。

Vistrailsを使った実行

Vistrailsを用いた計算はdmrg-01-dmrg.vtを使用します。"spin 1/2 multiple"を実行してください。

多重実行

コマンドラインでの実行

次のようにスピンS = 1/2の場合と同じで、単一のパラメータ·ファイルに複数の実行を設定することができます。

LATTICE_LIBRARY="my_lattices.xml"
LATTICE="open chain lattice with special edges 32"
MODEL="spin"
local_S0=0.5
local_S1=1
CONSERVED_QUANTUMNUMBERS="N,Sz"
Sz_total=0
J=1 
NUMBER_EIGENVALUES=1 
SWEEPS=4
{ MAXSTATES=20 } 
{ MAXSTATES=40 }
{ MAXSTATES=60 }

Pythonを使った実行

Pythonを使った実行には、spin_one_multiple.pyのスクリプトを参照してください。

Vistrailsを使った実行

Vistrailsを用いた計算はdmrg-01-dmrg.vtを使用します。"spin 1 multiple"を実行してください。

サイト（ボンド）での基底状態エネルギー

ハミルトニアンをよく見れば、鎖長$L$のエネルギーは $L$サイト上ではなく、$L-1$であることが分かると思います。まず最初に最後におこなったシミュレーションと計算$e_0 = \frac{E_0(L)}{L-1}$を見てみます。概要に記載されている値に近い値が得られますか？間違っている基本的な仮説は何でしょうか？正しい方法は、鎖の中央の結合エネルギーを考慮し、開境界条件の影響を排除することです。それをおこなう方法は２つあります。

1. $L$と$L+2$の鎖長の定常状態のエネルギーを計算し、既出の鎖長で、結合あたりのエネルギー$e_0 = (E_0(L+2) - E_0 (L))/2$を計算してみてください。結果はどのようになりますか？
2. 低コストで、通常の方法としては、隣接するサイト間の相関を用います。

$e_0 = \frac{1}{2} (\langle S^+_i S^-_{i+1}\rangle + \langle S^-_i S^+_{i+1}\rangle ) + \langle S^z_i S^z_{i+1} \rangle$
$i,i+1$は鎖中心のサイトです。