HΦの演習問題 (解説)三澤 貴宏
東京大学物性研究所 特任研究員（PCoMS PI）

例題1: spin 1/2 dimer (full diag)
例題2: spin 1/2 chain (Lanczos)
例題3: J1-J2 Heisenberg model(Lanczos,TPQ)
例題4: Kitaev model (Lanczos,TPQ)
例題5: Hubbard chain (Lanczos,TPQ)
好きなものからやって下さい
ほとんどlaptop PCでできるはずです(TPQはちょっと重いですが…)
もちろん、自分のやりたい別の課題もやっても OKですBasic exercise !

いくつかのサンプルは以下のサイトにあります
https://github.com/yomichi/HPhi-samples

例題 1: 反強磁性 S = 1/2 Heisenberg dimer

最も簡単な量子系のひとつである 反強磁性 2 S = 1/2Heisenberg dimer model
H = S⃗1 · S⃗2

ただひとつの基底状態とみっつの縮退した励起状態とがある
基底状態は全スピン Stot = 0 で、エネルギーは E = −3/4
励起状態は全スピン Stot = 1 で、エネルギーは E = 1/4

最初の一歩としてこれを HPhi で解く

2以下特筆なければすべて反強磁性模型です©Y.Motoyama
2 / 17

例題 1: S = 1/2 Heisenberg dimer

入力ファイル SpinHalf.def
1
2
3
4
5
6model = " SpinGC "lattice = " chain ␣ lattice "method = " FullDiag "2S = 1
L = 2
J = 0.5

これを HPhi の Standard モードで処理するだけ
SpinGC の GC は Grand CanonicalJ = 0.5 なのは周期的境界条件のためボンドが 2 本あるから
1> HPhi -s SpinHalf . def余裕があれば, S=1, S=3/2…などをやってみましょう(2S=2などとするだけ)
Emin=-S(S+1), Emax=S^2になるはずです©Y.Motoyama
3 / 17

例題 2: S = 1/2 Heisenberg chain

入力ファイル StdFace.def
1234567

model = " Spin "lattice = " chain ␣ lattice "method = " Lanczos "2S = 12 Sz = 0
J = 1.0
L = 12Spin は Canonical
2Sz は 計算したい状態の全スピンの z 成分の 2 倍
今回は低エネルギー状態のみでいいので Lanczos 法を使う

1> HPhi -s StdFace . def©Y.Motoyama
6 / 17

例題 2-A: S = 1/2 Heisenberg chain

S = 1/2 Heisenberg chain
H=!L⃗

⃗

i=1 Si · Si+1
周期境界 SL+1 = S1 を課す（以下同様）

基底状態は朝永・Luttinger 流体で、ギャップレス
励起は線形分散を示す
有限系では状態の取れる波数が離散化されるため、エネルギーギャップが開く
分散が線形なので、∆ ∝ 1/L となる 3有限系のエネルギーギャップを HPhi で計算する

3

実際には対数補正がかかる©Y.Motoyama
5 / 17

例題 2: S = 1/2 Heisenberg chain

入力ファイル StdFace.def
12
3
4
5
67

model = " Spin "lattice = " chain ␣ lattice "method = " Lanczos "2S = 12 Sz = 0
J = 1.0
L = 12Spin は Canonical
2Sz は 計算したい状態の全スピンの z 成分の 2 倍
今回は低エネルギー状態のみでいいので Lanczos 法を使う

1> HPhi -s StdFace . def©Y.Motoyama
6 / 17例題 2: S = 1/2 Heisenberg chain
今回も大量の出力とともに一瞬で計算が終わる（はず）
必要な情報は output/zvo_Lanczos_Step.dat に出力されている 41
23
4> tail -n 3 output / zvo_Lanczos_Step . dat
stp =42 -5.3873909174 -5.0315434037-4.7773893336 -4.5693744101
stp =44 -5.3873909174 -5.0315434037-4.7773893337 -4.5693744108
stp =46 -5.3873909174 -5.0315434037-4.7773893337 -4.5693744108

これは各 “Lanczos step” で計算した「ちいさい固有値」が幾つか並んでいる

デフォルト設定では第一励起エネルギーの収束をチェックしている

これの左 2 つの差がエネルギーギャップ
4

一応ログ出力にもほぼ同じのが出ている©Y.Motoyama
7 / 17

例題 2: S = 1/2 Heisenberg chain
データ抽出の自動化

一回だけなら電卓に手入力でもいいけれど、超絶面倒くさいので自動化すべき
tail と awk でワンライナーを書くのが楽 5
12> tail -n 1 output / zvo_Lanczos_Step . dat| awk ’{ print ␣ $3 - $2 } ’
0.355848

長さ L を変えながらエネルギーギャップ ∆ を何度も計算して、ギャップの長さ依存性を見てみよう
ノート PC、特に仮想マシンでは L = 18, 20 ぐらいで止めておくのが賢明

それ以上のサイズは自習時間にがんばってください

L を変えて何度も計算するという操作も自動化を試みてみよう→ 次のページ

5

自分の手に馴染んでいるなら何使ってもいいです©Y.Motoyama
8 / 17

例題 2: S = 1/2 Heisenberg chain
入力の自動化

今回はひとつのパラメータを変えるだけなのでシェルで事足りる共通となるパラメータをくくりだして保存しておく(StdFace.common)
12
3
4
56

model = " Spin "lattice = " chain ␣ lattice "method = " Lanczos "2S = 12 Sz = 0
J = 1.0

別々の L 毎にこの共通ファイルをコピーして、足りないパラメータを追記すればよい
1234
5678910rm -f res . dat
for L in 10 12 14 16; docp StdFace . common StdFace . def
echo " L = $L " >> StdFace . def
HPhi -s StdFace . def
gap = $ (( tail - n1 \output / zvo_Lanczos_Step . dat |\awk ’{ print ␣ $3 - $2 } ’)echo $L $gap >> res . datdone

ヒアドキュメントを使うのもアリ。自分の好きなようにやりましょう©Y.Motoyama

9 / 17

例題 2: S = 1/2 Heisenberg chain
結果

0.5

a + b*x
Lanczos

0.4はギャップの閉じるぎりぎりのところにあり、log L の対数補正が現れるので結構数値解析が難しい\DeltaS = 1/2 Heisenberg chain

0.3
0.2
0.1
000.020.04

0.06

0.080.11/L

S = 1/2 XY chain (Jz = 0, Jx = 1.0) ではもっと綺麗にかける→ 練習問題・S=1 Heisenberg chainにするとギャップが見えるはず(Haldane gap)
→練習問題 [input fileで 2S=2 とするだけ]©Y.Motoyama

10 / 17

Basic exercise 2 !

例題3: J1-J2ハイゼンベルク模型
H = J1X

Si Sj + J 2hi,ji

G, HONG YAO, AND LEON BALENTSX

S i Sjhhi,jii

PHYSICAL REVIEW B 86, 024424 (201

最近接 J1,次近接J2

lattice.gpで描画可能

0ms(k0)0.3∆Tms(kx)

0.3∆S120.24012ms

80.2

0.1

0.0
0.0

Neel AFMStripe AFM

0.1

0.8

0.0
1.0

QSL

0.2

0.4

J2/J10.6∆T, ∆S

ture on the 3intermediatephaseof2 the J13-J201
on the possibility of symmetry breaking
f these priorstudiedhavesuggestedthat
1513141512
ate has VBS order. We note, however,results for the J1 -J2 model are based
891011 orechniques 11(such asseries
expansion
thods, or fixed node or related versions
apted to avoid
thewhich
4 sign 5 problem,
677
ased Monte Carlo in this system), or on
n of very smallsystems.
Some theoretical
01233
possibility of VBS order comes from theed quantum criticality,28 which predicts
1512131415
quantum phasetransition—adeconfinedoint (DQCP)—should occur between an
and a plaquette or columnar VBS state,owever, the existence of such a transition
imply that it occurs for the J1 -J2 model inis particular model even harbors a VBS
etical motivation for VBS order comes024424(2012)
FIG. 1. PRB(Color 86,online)
The ground-state phase diagram for t
1spin- 2 AFM Heisenberg J1 -J2 model on the square lattice,determined by accurate DMRG calculations on long cylinders wi
Ly up to 14. Changing the coupling parameter J2 /J1 , three differe
phases are found: Néel antiferromagnet (AFM), topological quantu
spin liquid (QSL), and stripe AFM phase. ms (k0 = (π,π )) [ms (kx(π,0))] denotes the staggered magnetization in the Néel AFM pha

J2/J1~0.5で非磁性の基底状態(スピン液体?)例題3: J1-J2ハイゼンベルク模型
H = J1

Xhi,jiSi Sj + J 2X

S i Sjhhi,jii

1. Lanczosでエネルギーを計算 (サイズ4 4位)
2. TPQで比熱を計算(J2/J1~0.5でどうなるか？)
3. (発展)余裕があればスピン相関も計算してみましょう
L =4
W = 4
model = "Spin"method = "Lanczos"
lattice = "square lattice"J = 2.0
J’ = 1.0
2Sz = 0can be a finite-size efI'ect. Then, we conjecture a phase diagram for this model as shown in Fig. 4(b). Since themagnetizations behave very smoothly we expect second-

例題3: 基底状態の答え

J1-J2 Heisenberg model, Ns=4×4, J1=2.0
E. Dagotto and A. Moreo, PRB (R) 39 , 4744 (1989)

TABLE I. Ground-state energy (Ep) and first excited-stateenergy (E~) per site (both singlets with zero momentum) of the2D Heisenberg model with frustration as a function of J2 on a4X41attice. The error is in the last digit.
Ep0.950-statez near-state). For(Es)nberg

1. 1001. 150
1.200
1.325
1.400
1.500
1.600
1.750

—1.065 978—1.047 189—1.047 183—1.051 792—1.089 305—1.127 716—1.188 546—1.254670—1.358437—1.0160—1.0254—1.0307—1.0380—1.0804—1.1169—1.1691—1.2233—1.3072

今ならPCで数秒で計算できる。

例題4: Kitaev model
H=Jx

XSix Sjx

XJyx bondSiz Sjzz bond

2010X

相図

Jx,Jy,Jzで相互作用
4Jzy bond3方向のそれぞれが5Siy Sjy

32A. Kitaev / Annals of Physics 321 (2006) 2–111

5Jz=1,

4

Jx= Jy=0

gapless

17131517161214161211791171068106

51351

4
12

130

2

401315

1713

14

16

12

lattice.gpで描画可能

Az

gappedB
Jx=1,
Jy= Jz=0Ax

Ay

Jy=1,Jx= Jz=0Fig. 5. Phase diagram of the model. The triangle is the section of the positive octant (Jx, Jy, Jz P 0) by
Jx + Jy + Jz = 1. The diagrams for the other octants are similar.

Annals of Physics 321, 2-111 (2016)is impossible, even if we introduce new terms in the Hamiltonian. On the other height copies of each phase (corresponding to diﬀerent sign combinations of Jx, Jy,
the same translational properties. It is unknown whether the eight copies of thephase are algebraically diﬀerent.
We now consider the zeros of the spectrum that exist in the gapless phase. Thetum q is defined modulo the reciprocal lattice, i.e., it belongs to a torus. We repre
momentum space by the parallelogram spanned by (q1, q2)—the basis dual to (nthe symmetric case (Jx = Jy = Jz) the zeros of the spectrum are given by

可解模型→スピン液体

例題4: Kitaev model
H=Jx

Xx bondSix Sjx

JyX

Siy Sjyy bondJz

XSiz Sjzz bond

1. Lanczosでエネルギーを計算 (サイズ18サイト位)
2. TPQで比熱を計算: マヨラナ粒子の兆候がみえるか？
3. (発展)次近接のスピン相関が厳密に0を確認
4. (発展)ハイゼンベルク項をたすとどうなるか？
5. (発展)磁場をかけて磁化の温度依存性から帯磁率が計算可能
W = 3
L = 3model = "SpinGC"
method = "Lanczos"
lattice = "Honeycomb"
J0x = -1.0
J0y = 0.0J0z = 0.0
J1x = 0.0
J1y = -1.0
J1z = 0.0
J2x = 0.0
J2y = 0.0J2z = -1.0

例題5: Hubbard chain
H=tX

hi,ji

†
(ci cj+ h.c.) + U

Xni" ni#i

1. Lanczosでエネルギー・二重占有度を計算 (サイズ8サイト位)
2. TPQで比熱・二重占有度を計算:
3. (発展) 全対角化でアンサンブル平均を計算してTPQと比較

L = 8
model = "FermionHubbard"
method = "Lanczos"
lattice = "chain"t = 1.0
U = 8.0
nelec = 82Sz = 0例題5: Hubbard chain (有限温度計算)
FullDiag, TPQ, Lanczosの比較
Hubbard model, L=8, U/t=8, half filling, Sz=02.5
20.25

TPQ
FullDiag
Lanczos

0.2

D/Ns

E/Ns

1.5
1
0.5

0.15
0.1
0.05

0-0.5
10-2

TPQ
FullDiag
Lanczos

10-1

10 0T/t

101

10

2

0
10-2

10-1

10 0T/t101

10 2

3つの手法は互いに一致その他の例
HPhi/samples/Standard/にHubbard模型、Heisenberg模型、Kitaev模型、近藤格子模型のStandard modeのStdFace.defがあります。

StdFace.defを適宜変えて遊んでみてください。
注意：-サイズは大きくし過ぎないこと(PCだとspin 1/2で24site, Hubbardで12サイトくらいが限界)-Lanczosはサイズが小さいと不安定です(数千次元以上を推奨)-TPQはサイズが小さいとうまくいきません(数千次元以上を推奨)

物性研スパコンを使ってみよう!
物性研システムB (sekirei)✓fat node: 1node (40 cores) 1TBのメモリ、2nodesまで使用可 → ~2TB
✓cpu node: 1node (24cores) 120GBのメモリ,144nodesまで使用可 →~17TB
潤沢なメモリ: spin 1/2 36 sites, Hubbard 18 sitesまでなら比較的簡単に計算可能 (5-10年程度前の最先端の計算)
TPQ法を使えば有限温度計算も可能！

お試しクラス(A) : 随時受付(100ポイント以下)
通常クラス(B,C,E): 年２回(6月,12月)受付
緊急クラス(D): 随時受け付け,ただ緊急性が必要物性研スパコンを使ってみよう!
物性研システムB (sekirei)にはインストール済み
以下のようなスクリプトをサブミットするだけで、大規模並列計算が可能#!/bin/sh#QSUB -queue F144cpu#QSUB -node 128#QSUB -mpi
128
#QSUB -omp
24
#QSUB -place pack
#QSUB -over false#PBS -l walltime=24:00:00#PBS -N HPhi
cd ${PBS_O_WORKDIR}
source /home/issp/materiapps/HPhi/HPhivars.sh
mpijob HPhi -s StdFace.defExpert mode !

How to use HΦ: What is Expert mode ?
HPhi -s StdFace.def
Standard mode:
計算に必要なファイルを自動生成
模型を指定するパラメーターファイル (3個)zInterAll.def,zTrans.def, zlocspn.def
計算条件を指定するパラメーターファイル (2個)
modpara.def,calcmod.def
相関関数を指定するパラメーターファイル (2個)
greenone.def, greentwo.def+ファイル名を列挙したファイル: namelist.def
Expert mode:計算に必要なファイルを自分で用意→Standard modeで生成したものを書き換えるのが楽How to use HΦ: What is Expert mode ?
Expert mode:計算に必要な以下のファイルを自分で用意→Standard modeで生成したものを書き換えるのが楽
模型を指定するパラメーターファイル (3個)zInterAll.def,zTrans.def, zlocspn.def
計算条件を指定するパラメーターファイル (2個)
modpara.def,calcmod.def
相関関数を指定するパラメーターファイル (2個)
greenone.def, greentwo.def
用意して以下のコマンドを実行
HPhi -e namelist.def

How to use HΦ: zInterall.def
模型を指定するパラメーターファイルの例
XH+ =i,j,k,lX

Iijkl 1 2 3 4 c†i 1 cj 2 c†k 3 cl 4

1, 2, 3, 4相互作用の総数
実部

iσ1

jσ2

kσ3l

虚部

σ4

「任意」の2体の相互作用を指定することができる(自動生成するようにしておくと便利)→任意の格子を指定することが出来る

How Exchange
to use HΦ:Expert mode
オンサイトクーロン相互作用をハミルトニアンに付け加えます(S = 1/2 の系
2.11
指定ファイル

使用可能、ver. 0.2 では MPI 非対応)。付け加える項は以下で与えられます
模型を指定する簡易版のパラメータファイル
Exchange カップリングをハミルトニアンに付け加えます(S = 1/2!
= )。電子系の場合にはUi ni↑ ni↓(4- CoulombIntra
用可能、ver. 0.2 では MPIH+
非対応
i!
=================================
NCoulombintra 2H+ =JijEx (c†i↑ cj↑ c†j↓ ci↓ + c†i↓ cj↓ c†j↑ ci↑ )にファイル例を記載します。
✓=================================
i,j===========Exchange=====================================================================
NCoulombIntra
604.0
付け加えられ、スピン系の場合には
14.0======================✏!
========i_0LocSpn_1IteElc ======+ −− +Ex
H+=J
(SS+S-Exchangeijiji Sj )======================
i,j0 =================================
4.000000
21 NExchange
4.000000+ −− +=================================付け加えられます。スピン系の(SS+Siji Sj ) を電子系の演算子で書
2 ===========Exchange==============
4.000000†† =================================††ccc+ccc34.000000
i↑ j↑ j↓ 0 i↓ 1 0.5i↓ j↓ j↑ ci↑ ) となることに注意して下さい。以下にファイ
2 0.5ます。
4 14.000000
指定が楽&高速化することも✓
5 4.000000✒======================

詳細はマニュアル参照✑

HΦを使いこなすには..
- 簡単なファイルの書き出し、読み込み&解析ができると計算後の解析の幅が格段に拡がります-もちろん、通常のプログラミング言語(fortran, C++/Cなど)でも十分できますが、スクリプト言語(python, perl, rubyなど)を覚えておくと便利です。sed、awkといったコマンドを使うのも楽ですが、拡張性が低いのが難点です-自分で気に入ったものを使うのが一番ですが、周りに詳しい人がいる言語を選択するのが無難でしょうか..
-計算物質科学の分野だと主流はpythonでしょうか..

あとはやりたいことをやってください- 近藤格子模型- カゴメ格子,三角格子…
解析したい模型などの相談・質問も大歓迎です!