ALPS 2 Tutorials:DMRG-01 DMRG/ja

From ALPS
Jump to: navigation, search


ハイゼンベルグスピン鎖

密度行列くりこみ群法(Density Matrix Renormalization Group method, DMRG) の適用例として、次のハミルトニアンで与えられる長さLのスピン1/2反強磁性ハイゼンベルグ鎖とスピン1反強磁性ハイゼンベルグ模型を考えます。

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} .

既に多くのチュートリアルで取り上げたこの2つの模型を選択した理由は、これらは見た目がほとんど同じですが、その物理的な振る舞いがまったく異なり、DMRG アルゴリズムにとっては全く別の問題になっているからです。これらの模型を簡単にレビューします。

スピン1/2鎖

スピン1/2鎖の基底状態はBethe仮説を使うことで厳密解が得られます。すなわち、基底状態エネルギーは厳密に知ることができます。熱力学的極限L\rightarrow\inftyのサイトあたりのエネルギーは次式で与えられます。

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

基底状態のエネルギーは他の状態と比較しないと、それ自体はあまり意味の無いものですが、DMRG法のベンチマークとしては非常に有益です。 物理的に興味のあることとしては、基底状態と励起状態の間には有限のエネルギーギャップが熱力学的極限においても存在するかしないのか、と言ったことがらがあります。スピン1/2鎖ではギャップは存在しません。

それと同時に、異なるサイト上のスピンの間の相関がどうなっているのか、ということにも興味が持たれます。 無限長スピン1/2鎖では、スピンの相関関数は漸近的に(つまり長距離極限 |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) であることを意味します。つまり、スピン間の反強磁性相関は距離とともにベキ的に減少します。この場合、ベキの指数は-1になります。 この相関関数には距離の対数の平方根という形で対数補正がかかっています。 この補正項は非常に長い鎖に対する素晴らしいDMRG 計算によって、数値的にも確かめられていますが、 非常にゆっくりとしか変化しないため、最初のうちは無視してしまうことにします。

スピン-1鎖

長年、スピン1鎖はスピン1/2鎖と定性的には同じ振る舞いをするのだと考えられてきました。 そのため、1982年のDuncan Haldaneによる、反強磁性ハイゼンベルグ鎖はスピン長が半整数(S=1/2,3/2,...) か整数か(S=1,2,...) によってその振る舞いに根本的な差異があるという主張は、強い驚きをもって迎えられました。 それ以降、スピン1鎖への注目が高まり、さまざまな研究がおこなわれており、 DMRG はその流れの中で非常に大事な役割を果たすことになりました。

スピン1/2鎖とは異なり、スピン1鎖には解析解がありません。そのため、定量的な研究には数値計算が必須となります。

サイトあたりの基底エネルギーは

 E_0 = -1.401484039 ...

です。

また、エネルギーギャップの存在が非常に重要です。これがまさにスピン1 とスピン1/2との定性的な違いとなっています。 熱力学的極限において、スピン1鎖のギャップは有限であり、その値は5桁の精度では

 \Delta = 0.41052

となっています。

スピン-スピン相関の振る舞いも、スピン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でおこる指数関数的な減衰が支配的です。この相関長と呼ばれる長さスケールは、スピン1鎖では数値的に\xi=6.02だと知られています。 距離の平方根というベキ的な補正が分母にあります。 しかし、これは指数関数的減衰にくらべるとはるかにゆっくりとした変化なので、 相関長を求める際にはたいてい無視されます。もちろん相関長が十分に長い場合には、考慮にいれる必要が出てきます。

スピン-1鎖は有限のエネルギーギャップおよび指数関数的に減衰する相関関数を持った非臨界量子系の代表例です。 これから明らかになっていきますが、DMRG は臨界系よりも非臨界系の方により適した手法です。

チュートリアルの進め方

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

Vive la difference ...

他の数値法と比べて最も重要な差異に境界条件があります。 厳密対角化などの他の数値手法や解析計算では周期境界条件を用いるのに対して、 DMRG では開境界条件が好まれるということです。つまり、サイト1とLの2つの鎖端が存在します。

コードの実行

概論

まず最初に、 DMRGアルゴリズムについて簡単に説明します。状態空間の次元がd=2S+1である局所系をつないで一次元量子系を作ると、全体のHilbert空間の次元はシステムサイズLに対して指数的d^Lに増大します。厳密対角化では、この指数関数的に大きなHilbert空間を対角化することで数値的に厳密な結果を得られますが、扱える系の大きさに強い制限がかかるという問題があります。 これに量子モンテカルロ計算は、この大きな状態空間から確率的にサンプリングして統計的な近似値を得ることで、より大きな系の計算を可能にします。 DMRG は別の方法で大サイズ計算を目指します。指数関数的に増大するHilbert空間の中から、非常に小さな、しかし目的の状態(たいてい基底状態)を十分近似した状態を含むだろうと期待される部分空間を選びとっていきます。

最初の制御パラメータはDです。これは部分空間の大きさを制御するパラメータです。DMRGの精度はこのパラメータに対して単調に増加していきます。D が大きくなれば部分空間が大きくなり、近似値もより正確になります。 D には上限があります。もしD\rightarrow d^Lとすると、すべての状態を捨てずに取り扱うこととなり、厳密な結果を得ることができます。しかし、このような大きな状態空間を扱えるのであれば、厳密対角化がよりよい手法となります。 なお、歴史的な事情で、Dは、matrix dimensionnumber of block statesなど、様々な名前で呼ばれています。

2番目の制御パラメータはシステムサイズLです。

3つ目はDMRGアルゴリズムを詳しく見ていくことで現れます。DMRG では最良の状態近似を求めるために次の2つのステップをおこないます。

  1. 最初のステップ(いわゆる無限系DMRG)では、2スピンから始めて部分空間を選びとり、その結果を用いてスピン数を4,6,... と逐次的に増やしていき、目的のシステムサイズLを目指します。増やす際には鎖を分割し、中央に2つの新しいサイトを挿入します。この手順はもちろんL を増やすことで限度なく継続できて、このことが名前の由来になっています。とはいえ、これで熱力学的極限の意味のある結果を得られるとは期待しないでください。2番目に注意する点は、この手順から明らかなように、DMRG 計算では偶数システムサイズのほうがよく扱えるということです。
  2. 2番目のステップ(いわゆる有限系DMRG)では、ステップ1で作った部分空間は、長さL の系におけるすべての量子ゆらぎや相関をまだ考慮できていないことを改善します。このステップでは部分空間の質を向上させるため、さらなる反復計算を行います。反復処理1回で全サイトを見ていきますが、これをDMRG ではsweepと呼びます。このsweepの回数が、3番目に重要な制御パラメータです。このパラメータが小さいと与えたDにおいて期待される精度が達成できないし、かと言って長すぎると計算時間が膨大になってしまいます。


最後に、DMRG 計算の正確さを見積もる良い指標である打ち切り誤差(truncation error)について考察します。 DMRG では、各操作で状態空間の大きさは指数関数的に増えようとするのですが、仮に状態空間の大きさを増やさなかった時にどれだけの精度が保たれるのかを考えます。 この見積りではエネルギー固有状態および固有値から作った密度行列の解析を行います。 DMRG の近似の度合いは、捨てられた状態の持っていた状態重みを反映するので、これを打ち切り誤差と呼びます。 多くの計算では打ち切り誤差は10^{-12} 程度に小さな量となっています。 これDMRG が実際にはほとんど近似をしていないことを示し、DMRG が多くの成功を収めている証左となっています。 本チュートリアルの重要なゴールの一つとして、局所的な物理量(磁化やボンドエネルギーなど)の誤差が打ち切り誤差に対して大雑把に比例することを確認する、ということを設定しておきます。


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

DMRG への入力として、ハミルトニアンや格子形状の他に、制御パラメータが必要になります。

DMRGパラメータ

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

SWEEPS
有限サイズアルゴリズムで行うDMRG掃引の回数。各掃引は、左から右への半掃引と、右から左への半掃引からなります。

NUM_WARMUP_STATES
無限サイズアルゴリズムで用いる初期状態の数。指定されていない場合、アルゴリズムはデフォルトとして20状態を使用します。

STATES
各半掃引で残しておく状態数。 全半掃引分の値、つまり2*SWEEPS 個のSTATES か、1つのMAXSTATES か、1つのNUMSTATES、 これら3通りのどれか1つを指定する必要があります。

MAXSTATES
掃引時に残しておく状態数の最大値。各半掃引で残しておく状態数は、MAXSTATES/(2*SWEEPS) ずつ増えていくように自動的に計算されます。

NUMSTATES
すべての半掃引で共通の、残しておく状態数。

TRUNCATION_ERROR
目標打ち切り誤差。これを達成するだけに必要な状態数が自動的に計算されます。 誤差の見積もりによっては、必要な状態数が大きくなりすぎてシミュレーションが異常終了することに注意してください。MAXSTATES を与えて状態数に制限を課すことが賢明です。

LANCZOS_TOLERANCE
対角化(Davidson/Lanczos)の精度。デフォルト値は10^-7です。

CONSERVED_QUANTUMNUMBERS
保存される量子数を指定します。これに基づいて系のブロック対角化を行います。 量子数を指定しつつ、具体的な値を指定しない場合は、グランドカノニカル計算を行います。 例えばスピン1/2鎖で、Sz_total を使わないと、全Hilbert空間 (dim=2^N) で計算を行います。 例えばSz_total に具体的な値を設定する、例えばSz_total=0とすると、カノニカル計算が行われます。 この時、Sz_total=0 の部分空間のみで計算が進むため、計算性能が向上します。

パラメータをどうとるか

デフォルトのパラメータをそのまま使うことはお勧めできません。 DMRG の収束は状態数や掃引数などに強く影響されます。 パラメータを決める良い方法として、基底エネルギーの収束の仕方や打ち切り誤差の値が状態数に対してどのように変化するのかを見ることです。 ある精度を達成するためには状態数をどのぐらいにすればよいのかがわかります。

掃引数が十分であるかどうかを判断するためには、 相関関数や局所物理量、例えばスピンのz成分や粒子数密度などの 空間分布を見るのがよいでしょう。 たとえば、空間反転対称な系では、先に述べたような局所的な物理量もまた空間反転対称性を持つはずです。 他にもエンタングルメントエントロピーも同じ対称性を持ちます。 もしもこれらの計算結果が期待される対称性を持っていなかったならば、 十分な掃引が行われていなかった可能性があります(他には、相分離が起きている可能性もあります)。

ハミルトニアンが保存量をもつ場合、例えばハイゼンベルグ模型におけるスピンのz成分Sz やハバード模型における粒子数N などがある場合、 これらの値を固定することで部分空間のみで計算が完結するため、計算時間や使用メモリを抑えることができます。

基底状態エネルギー

最初に基底状態とそのエネルギーについて考えます。 (基底状態の)エネルギーについては、有限系全体のエネルギーと、(熱力学的極限における)サイトあたり、もしくはボンドあたりのエネルギーの2種類のエネルギーのどちらを考えるのかをはっきりさせておく必要があります。

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

長さ L = 32, 64, 96, 128 の一次元系で、 スピンの大きさは1/2, 1 のそれぞれについて、 状態数D=50,100,150,200,300 として 基底エネルギーを計算します。 それぞれの長さに対し、得られた打ち切り誤差と基底エネルギーをD の関数としてまとめます。 それぞれの計算で結果がちゃんと収束するように掃引数を注意深く調整しましょう。

  1. 各鎖長、各スピン長で、打ち切り誤差に対する全エネルギーをプロットしてください。
  2. D に対して物理量は収束しますが、その収束の仕方がシステムサイズの変化によってどう変わっていくかをスピン1/2, スピン1 のそれぞれで観察してみましょう。スピンの大きさによって、どのような違いが見受けられましたか? ヒント: スピン1/2 のほうがスピン1 よりも強い変化を示すはずです。スピン1 の系は非臨界で相関長が有限なので、系の性質は相関長程度の大きさのセグメントで決まるため、系を大きくしてもあまり変化はありません。一方スピン1/2 の系は臨界系で相関長は発散しており、特徴的な長さスケールが存在しません。この違いが収束の違いにつながります。
  3. D\rightarrow\inftyの極限での各鎖長の基底エネルギーを外挿してみてください。

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

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

コマンドラインでの実行

パラメータファイルは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}

入力XMLファイルを生成し、dmrg を実行します。

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

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

4回の掃引(左から右へ向かうの半掃引と右から左への半掃引をそれぞれ4回ずつ)を実行します。それぞれの半掃引で残す基底の数は、設定された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

入力パラメータをPython の辞書形式(のリスト)で作成し、

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
      } ]

XML ファイルを生成、dmrg に渡します。

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

次に、計算結果を読み出して、

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'])

これを使えば、計算がどのように収束しているのかを見ることができます。

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ファイルを使います。 ここでは、L=32 の系の計算を、保持状態数を変えていくつか実行します(チュートリアルなので、状態数は少なめにとってあります)。

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が計算タスクごとに生成されます(#=1,2,3)。

Pythonを使った実行

Pythonを使った実行には、spin_one_half_multiple.pyのスクリプトを使います。 入力パラメータとしてMAXSTATES の値のみが違う3つの辞書からなるリストを用意します。

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ハイゼンベルグ鎖では開境界条件に対処するためのいくつかの特別な対処が必要になります。 あとで説明するように、鎖の両端は2つの大きさS=1/2 のスピンにする必要があります。これには新たな格子を定義する必要がありますが、異なる長さの系について同時に表す方法が無く、それぞれの長さに対して別々の格子を手で定義する必要があります。簡単のために長さを入力すると対応する格子ファイルを生成するPython スクリプトを用意しました。 例えば長さ6の鎖を生成するには

$ 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>

です。 格子は、6つのサイト(バーテックス)と、隣接する2点をつなぐボンド(エッジ)からなる1次元グラフとして定義されます。 最初と最後のサイトはタイプ"0"でその他は"1"です。 この格子定義を用いて今回の1次元開放端S=1 ハイゼンベルグ模型を定義するために、 それぞれのサイトの上にあるスピン自由度をパラメータとして与える必要があります。

local_S0=0.5
local_S1=1

32 サイト系を定義して my_lattice.xml に保存します。

$ 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

実際の計算では、格子ファイルを生成し、このパラメータファイルで計算を行うという手順を各システムサイズで繰り返すという、明らかに面倒なことをする必要があります。 これに対する第一の単純な解法は、この処理を自動で行うスクリプトを組むことです。 もっと単純な第二の解法は、複数の大きさの格子定義に違う名前をつけて、1つの格子ファイルにまとめて保存することです。 L=32,64,96,128,192 の格子定義をまとめたファイル my_lattices.xml が同梱されているので、これを使ってください。 これを使えば、格子の名前としてシステムサイズを与えることができます。

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
      } ]

パラメータや格子ファイル名以外はスピン1/2 のスクリプト 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のスクリプトを使います。 スピン1/2 のものとパラメータをのぞいて同じです。

Vistrailsを使った実行

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

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

ハミルトニアンをよく見ると、エネルギーはL 個のサイトではなくてL-1 本のボンドに付随していることがわかります。 これを踏まえ、エネルギー密度の最初の(単純な)評価として

e_0 = \frac{E_0(L)}{L-1}

を採用してみます。 実際にこれを計算して、上に示した値と比較してみてください。 結果はどうでしたか?ダメだとしたら、一体何が間違っていたのでしょうか?

先ほどのエネルギー評価子は境界の影響を受けてしまいます。 この効果を取り除くためには、系の中心付近のボンドのエネルギーを直接見る必要があり、 それには2通りの方法があります。

  1. LL+2 の2つの系のエネルギーを計算して、

e_0 = (E_0(L+2) - E_0 (L))/2 をボンドあたりのエネルギーとします。

  1. もっと計算コストが軽くて、よく使われるのは、系の中心でのスピンの短距離相関を測ることでボンドエネルギーを直接評価する方法です。この手法は別のチュートリアルで詳しく取り扱います。

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