LM-01 ハニカム反強磁性体

カスタム格子の定義：ハニカム反強磁性体

このチュートリアルでは、ALPSの標準ライブラリに含まれていない格子形状を定義し、量子スピンモデルと組み合わせ、モンテカルロシミュレーションを実行する方法を学びます。具体例としてハニカム（六角）格子を使用します。これはグラフェン、Kitaevハニカムモデル、および多くのフラストレート磁性体の基礎となる構造です。

このチュートリアルを終えると、以下のことができるようになります：

ALPSの格子ライブラリXMLファイルをゼロから記述する。
シミュレーションパラメータファイルでそのファイルを参照する。
spinmcを使って自作格子上で古典Heisenbergモデルを実行する。
ALPSの組み込みハニカム格子と比較して定義の正確さを検証する。

背景：ハニカム格子

ハニカム格子は、三角ブラベー格子と2サイト基底からなる2次元格子です。各単位格子にはAサブ格子とBサブ格子と呼ばれる2つの非等価サイトが含まれ、各サイトはちょうど3つの最近接格子点を持ちます。 AサイトはBサイトとのみ結合し、その逆も同様であるため、ハニカム格子は二部格子であり、反強磁性秩序の舞台として適しています。

原始格子ベクトルは

\mathbf{a}_1 = a(1, 0), \qquad \mathbf{a}_2 = a\left(\tfrac{1}{2}, \tfrac{\sqrt{3}}{2}\right)

であり、 $a$ は格子定数です（以下では $a=1$ とします）。 2つの基底サイトの位置は分数座標で

\text{A}: \tfrac{2}{3}\mathbf{a}_1 - \tfrac{1}{3}\mathbf{a}_2 = \left(\tfrac{2}{3}, -\tfrac{1}{3}\right), \qquad \text{B}: \tfrac{1}{3}\mathbf{a}_1 + \tfrac{1}{3}\mathbf{a}_2 = \left(\tfrac{1}{3}, \tfrac{1}{3}\right)

と表されます。これらの座標は六角形の中心を単位格子の原点に置きます。

ステップ1：格子ライブラリファイルの作成

ALPSは<LATTICES>ルート要素を持つXMLファイルから格子定義を読み込みます。 ALPSパーサーが直接読み込めるよう、ファイルは裸のXMLでなければなりません（<?xml ...?>宣言や<!DOCTYPE>前文は不要です）。

以下の内容でmy_honeycomb.xmlというファイルを作成してください：

<LATTICES>

  <!-- 1. ブラベー格子：三角格子、基底ベクトルa1とa2 -->
  <LATTICE name="triangular" dimension="2">
    <BASIS>
      <VECTOR>1 0</VECTOR>
      <VECTOR>0.5 0.8660254037844387</VECTOR>
    </BASIS>
  </LATTICE>

  <!-- 2. 単位格子：2サイト（AサブとBサブ格子）と3本の最近接結合 -->
  <UNITCELL name="honeycomb" dimension="2" vertices="2">
    <VERTEX id="1" type="0"><COORDINATE>0.6666666666666666 -0.3333333333333333</COORDINATE></VERTEX>
    <VERTEX id="2" type="1"><COORDINATE>0.3333333333333333 0.3333333333333333</COORDINATE></VERTEX>
    <EDGE type="0"><SOURCE vertex="1"/><TARGET vertex="2"/></EDGE>
    <EDGE type="1"><SOURCE vertex="2"/><TARGET vertex="1" offset="0 1"/></EDGE>
    <EDGE type="2"><SOURCE vertex="1"/><TARGET vertex="2" offset="1 -1"/></EDGE>
  </UNITCELL>

  <!-- 3. 周期境界条件付きL x Wハニカム格子 -->
  <LATTICEGRAPH name="honeycomb L x W">
    <FINITELATTICE>
      <LATTICE ref="triangular"/>
      <PARAMETER name="W" default="L"/>
      <EXTENT dimension="1" size="L"/>
      <EXTENT dimension="2" size="W"/>
      <BOUNDARY type="periodic"/>
    </FINITELATTICE>
    <UNITCELL ref="honeycomb"/>
  </LATTICEGRAPH>

</LATTICES>

各要素の役割

<LATTICE name="triangular"> は基底ベクトルを指定して無限ブラベー格子を定義します。第2基底ベクトル $(0.5,\, \sqrt{3}/2)$ は第1ベクトルに対して60°傾いており、三角格子配置を生み出します。

<UNITCELL name="honeycomb"> は2つの基底サイトを配置し、それらを最近接格子点に結ぶ3本の辺を宣言します：

辺タイプ	始点	終点	意味
0	A（格子(0,0)）	B（格子(0,0)）	単位格子内の結合
1	B（格子(0,0)）	A（格子(0,1)）	$\mathbf{a}_2$ 方向の結合
2	A（格子(0,0)）	B（格子(1,−1)）	$\mathbf{a}_1 - \mathbf{a}_2$ 方向の結合

これら3本の辺により、各AサイトはBの3近接格子点を持ち、各BサイトはAの3近接格子点を持ちます。

頂点座標は単位格子の分数座標（ブラベー格子ベクトルの係数）で与えられます。

<LATTICEGRAPH name="honeycomb L x W"> は有限シミュレーションセルを組み立てます。 WのデフォルトをLに設定することで、パラメータファイルにL=4のみを指定すれば $4\times 4$ の単位格子配置、すなわち $4 \times 4 \times 2 = 32$ サイトが得られます。 <BOUNDARY type="periodic"/> は両次元をトーラスに巻きつけます。

ステップ2：パラメータファイルの作成

my_honeycomb.xmlと同じディレクトリにparm_honeycombというパラメータファイルを作成します：

LATTICE_LIBRARY="my_honeycomb.xml"
LATTICE="honeycomb L x W"
MODEL="Heisenberg"
J=-1
L=4
W=4
UPDATE="cluster"
THERMALIZATION=5000
SWEEPS=50000
{T=0.1;}
{T=0.3;}
{T=0.5;}
{T=0.8;}
{T=1.0;}
{T=1.5;}
{T=2.0;}
{T=3.0;}

重要なポイント：

LATTICE_LIBRARYでALPSにカスタム格子定義ファイルの場所を伝えます。
LATTICEは<LATTICEGRAPH name="...">で指定した名前と完全に一致させる必要があります。
MODEL="Heisenberg"は組み込みの古典Heisenbergモデルを使用します。 spinmcの規約ではJ=-1が強磁性結合となりますが、磁化率ピークの振る舞いは対称なので $J = \pm 1$ どちらでも使えます。
L=4、W=4は $4 \times 4$ 単位格子の配置（32サイト）を選択します。
UPDATE="cluster"はWolffクラスター更新を使用します。秩序温度付近では局所更新よりもはるかに効率的です。

ステップ3：シミュレーションの実行

パラメータファイルをXMLに変換してspinmcを実行します：

parameter2xml parm_honeycomb
spinmc --Tmin 5 --write-xml parm_honeycomb.in.xml

シミュレーションは8つのタスク（温度ごとに1つ）を順に実行し、 parm_honeycomb.task1.out.xmlからparm_honeycomb.task8.out.xmlに結果を書き出します。

ステップ4：結果の抽出

以下のPythonコードで出力ファイルから磁化率とエネルギー密度を読み込みます：

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

data = pyalps.loadMeasurements(
    pyalps.getResultFiles(prefix='parm_honeycomb'),
    ['Susceptibility', 'Energy'])

chi = pyalps.collectXY(data, x='T', y='Susceptibility')
energy = pyalps.collectXY(data, x='T', y='Energy')

plt.figure()
pyalps.plot.plot(chi)
plt.xlabel('Temperature $T/|J|$')
plt.ylabel('Susceptibility $\\chi$')
plt.title('Classical Heisenberg on honeycomb (L=4)')
plt.show()

期待される出力（モンテカルロのノイズにより数値は若干変動します）：

  T |          chi |       E/site
------------------------------------
0.1 |     0.111780 |    -1.401380
0.3 |     0.122365 |    -1.191363
0.5 |     0.134909 |    -0.941824
0.8 |     0.145397 |    -0.595319
1.0 |     0.140822 |    -0.478350
1.5 |     0.120825 |    -0.326241
2.0 |     0.105280 |    -0.247001
3.0 |     0.080478 |    -0.165197

磁化率 $\chi(T)$ は $T \approx 0.8\,|J|$ 付近でピークを示し、これは古典Heisenbergハニカムモデルにおける短距離反強磁性相関の発現と一致します。高温では $\chi$ はキュリー則 $\chi \propto 1/T$ に近づき、エネルギー密度はゼロに収束します。

ステップ5：組み込み格子との比較検証

ALPSには標準ライブラリとして"honeycomb lattice"の定義が含まれています。カスタムXMLファイルの正確さを素早く確認するために使用できます：

LATTICE="honeycomb lattice"
MODEL="Heisenberg"
J=-1
L=4
W=4
UPDATE="cluster"
THERMALIZATION=5000
SWEEPS=50000
{T=0.5;}
{T=1.0;}
{T=2.0;}

LATTICE_LIBRARYがないことに注目してください——組み込み格子には追加ファイルは不要です。このパラメータファイルでspinmcを実行すると、カスタム格子の結果と統計誤差（数パーセント）の範囲内で一致する磁化率が得られるはずです。

発展問題

配位数：ハニカム格子の配位数は $z = 3$ です。磁化率のピーク温度は正方格子（ $z = 4$ ）と比べてどうでしょうか？ LATTICE="honeycomb L x W"をLATTICE="square lattice"（追加ファイル不要）に置き換えて計算を繰り返してみてください。
有限サイズ効果：Lを4から6、8に増やしてみてください。磁化率とエネルギー密度はどう変化しますか？どの温度範囲で $4 \times 4$ セルが熱力学的極限の良い近似になりますか？
結合の異方性：単位格子には3種類の辺タイプ（0、1、2）が定義されています。各結合に異なる結合定数 $J_0$ 、 $J_1$ 、 $J_2$ を割り当てるようハミルトニアンを変更してください。 ALPSのスピンモデルではJ0、J1、…パラメータで結合タイプごとの強度を指定できます——モデル定義ドキュメントを参照してください。 $J_2 = 0$ （二量体鎖に縮退）にすると磁化率はどうなりますか？
量子モデル：MODEL="Heisenberg"をMODEL="spin"に変え、local_S=0.5とALGORITHM="loop"を追加し、実行ファイルをloopに切り替えてみてください。カスタムハニカム格子上で量子 $S = 1/2$ Heisenberg反強磁性モデルを実行できます。古典モデルと量子モデルの磁化率曲線を比較してください。

LM-02 カスタムモデルハミルトニアン