LM-01 蜂窝反铁磁体

定义自定义晶格：蜂窝反铁磁体

本教程将介绍如何定义一个不在 ALPS 标准库中的晶格几何，将其与量子自旋模型关联，并在此基础上运行蒙特卡洛模拟。我们以蜂窝（六角）晶格为具体示例——它是石墨烯、Kitaev 蜂窝模型以及众多阻挫磁体的基础结构。

完成本教程后，你将能够：

从零开始编写 ALPS 晶格库 XML 文件。
在模拟参数文件中引用该文件。
使用 spinmc 在自定义晶格上运行经典 Heisenberg 模型。
与 ALPS 内置蜂窝晶格对比，验证自定义定义的正确性。

背景：蜂窝晶格

蜂窝晶格是一种二维晶格，具有三角 Bravais 格子和两格点基元。每个原胞包含两个不等价格点——通常称为 A 和 B 子格——每个格点恰好有三个最近邻。由于 A 格点只与 B 格点相连，反之亦然，蜂窝晶格是二分的，天然适合承载反铁磁有序。

原始格矢为

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

其中 $a$ 为晶格常数（下文取为 1）。两个基元格点的位置为

\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),

以原胞分数坐标表示，即以 Bravais 格矢为基的坐标系。这些坐标将六边形中心置于原胞原点。

第一步：编写晶格库文件

ALPS 从具有 <LATTICES> 根元素的 XML 文件中读取晶格定义。该文件必须为裸 XML——不含 <?xml ...?> 声明和 <!DOCTYPE> 前导声明——以便 ALPS 解析器直接读取。

新建文件 my_honeycomb.xml，内容如下：

<LATTICES>

  <!-- 1. Bravais 格子：三角格，基矢为 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"> 通过给出基矢定义无限 Bravais 格子。第二个基矢 $(0.5,\, \sqrt{3}/2)$ 相对于第一个基矢倾斜 60°，产生三角铺排。

<UNITCELL name="honeycomb"> 放置两个基元格点，并声明将它们连接到最近邻的三条棱：

棱类型	源	目标	含义
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$ 方向的键

这三条棱使每个 A 格点有三个 B 邻居，每个 B 格点有三个 A 邻居。

格点坐标以原胞分数坐标（即 Bravais 格矢的系数）给出。

<LATTICEGRAPH name="honeycomb L x W"> 组合出有限模拟元胞。将 W 默认设为 L 意味着在参数文件中只指定 L=4，就会产生 $4\times 4$ 的原胞排列，共 $4 \times 4 \times 2 = 32$ 个格点。 <BOUNDARY type="periodic"/> 将两个维度都卷成环面。

第二步：编写参数文件

在与 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 模型。 J=-1 在 spinmc 约定下为铁磁耦合；磁化率峰值行为是对称的，可使用 $J = \pm 1$ 。
L=4、W=4 选择 $4 \times 4$ 原胞排列（32 个格点）。
UPDATE="cluster" 使用 Wolff 团簇更新，在有序温度附近比局部更新高效得多。

第三步：运行模拟

将参数文件转换为 XML 格式并运行 spinmc：

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

模拟依次运行八个任务（每个温度一个），结果写入 parm_honeycomb.task1.out.xml 到 parm_honeycomb.task8.out.xml。

第四步：提取结果

下面的 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$ ，而能量密度趋向零。

第五步：与内置晶格对比验证

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$ 的元胞是热力学极限的良好近似？
键各向异性：原胞定义了三种棱类型（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 自定义模型哈密顿量