ED-06 Full Diagonalization

在本教程中，我们将使用 fulldiag 应用程序计算小型量子自旋系统的完整能谱，并用它计算精确的有限温度热力学量——能量、熵、比热和磁化率——作为对所有本征态的加权求和，完全没有量子蒙特卡罗方法的统计噪声，也没有 DMRG 的截断误差。这补充了前面的教程：那些教程使用 sparsediag 只针对哈密顿量的最低若干本征态；而这里我们需要整个能谱，因为在任意有限温度下，无论能量多高的态都携带玻尔兹曼权重 $e^{-E_n/T}$ 。

一维自旋模型的热力学

S=1 海森堡链

我们的第一个模拟针对的是 S=1 海森堡量子反铁磁链，

H = J \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j , \qquad J=1 ,

定义在一条周期链上，其有限链热力学最早由 J.C. Bonner and M.E. Fisher, Phys. Rev. 135, A640 (1964) 系统地研究。

参数

参数	含义	取值
`LATTICE`	周期链	`chain lattice`
`MODEL`	量子自旋模型	`spin`
`local_S`	每个格点的自旋量子数	1
`J`	最近邻耦合	1
`L`	链长	8
`CONSERVED_QUANTUMNUMBERS`	用于将 $H$ 分块对角化的对称性	`Sz`

格子

与本系列教程中一直使用的相同的周期性 chain lattice（参见 ALPS 格子库）：

    J     J     J     J     J     J     J
o-------o-------o-------o-------o-------o-------o-------o
0       1       2       3       4       5       6       7
|_______________________________________________________|
                      J   (bond 7-0, periodic)

方法

热力学量需要对整个能谱（而不仅仅是基态）求玻尔兹曼权重之和，因此这里使用 fulldiag 而非 sparsediag。这个 S=1、L=8 链的完整希尔伯特空间维数为 $3^8=6561$ ；使用 CONSERVED_QUANTUMNUMBERS=Sz 将其分解为 17 个小得多的块（最大的块，即 $S_z=0$ 子空间，最多有 $1107$ 个态），每个块都被直接对角化，这既更快，又足以在任意温度下重构出精确的配分函数。

使用命令行

参数文件 parm6a 为一维 8 格点上的 S=1 海森堡量子反铁磁体设置了全对角化：

LATTICE="chain lattice"
MODEL="spin"
local_S = 1
J       = 1
CONSERVED_QUANTUMNUMBERS="Sz"
{L = 8}

CONSERVED_QUANTUMNUMBERS 参数帮助 fulldiag 将希尔伯特空间分解为不变子空间，并分别对其进行对角化。按照标准的命令序列，你可以先将输入参数转换为 XML，然后使用 fulldiag 计算该量子哈密顿量的完整能谱：

parameter2xml parm6a
fulldiag --write-xml parm6a.in.xml

现在输出文件中包含了所有本征矢量的结果，你可以使用 fulldiag_evaluate 为热力学和磁学可观测量生成 XML 绘图文件：

fulldiag_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T  0.1  parm6a.task1.out.xml

这将生成以下 XML 绘图文件：

parm6a.task1.plot.energy.xml
parm6a.task1.plot.free_energy.xml
parm6a.task1.plot.entropy.xml
parm6a.task1.plot.specific_heat.xml
parm6a.task1.plot.uniform_susceptibility.xml
parm6a.task1.plot.magnetization.xml

要从 fulldiag_evaluate 生成的 XML 绘图文件中提取计算结果，可以使用 plot2text 工具，然后用你喜欢的绘图工具查看这些数据。例如，要提取能量密度随温度变化的数据，可以使用

plot2text parm6a.task1.plot.energy.xml

以类似的方式，你可以从其他 XML 绘图文件中提取数据。如果 Grace 是你喜欢的绘图工具，你也可以使用 plot2xmgr 工具直接从 XML 绘图文件生成 Grace 项目文件。例如，要生成能量随温度变化的 Grace 项目文件，可以使用

plot2xmgr parm6a.task1.plot.energy.xml > energy.agr

类似地，工具 plot2gp 生成 Gnuplot 脚本，plot2text 将文件转换为纯文本。不过，进行数据处理和绘图更推荐使用 Python。

使用 Python

要在 Python 中设置并运行该模拟，我们使用脚本 tutorial6a.py。该脚本的前几部分导入所需模块，然后将输入参数准备为一个 Python 字典列表，并运行模拟：

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
parms = [{ 
        'LATTICE'                   : "chain lattice", 
        'MODEL'                     : "spin",
        'CONSERVED_QUANTUMNUMBERS'  : 'Sz',
        'local_S'                   : 1,
        'J'                         : 1,
        'L'                         : 8
    }]

input_file = pyalps.writeInputFiles('parm6a',parms)
res = pyalps.runApplication('fulldiag',input_file)

接下来我们对所有输出文件运行评估程序：

data = pyalps.evaluateFulldiagVersusT(pyalps.getResultFiles(prefix='parm6a'),DELTA_T=0.1, T_MIN=0.1, T_MAX=10.0)

最后显示所有图像：

for s in pyalps.flatten(data):
    plt.figure()
    plt.title("Antiferromagnetic Heisenberg chain")
    pyalps.plot.plot(s)
    plt.show()

输出数据

独立对角化同一个哈密顿量，得到以下每格点的热力学量（能量 $E/L$ 、比热 $C/L$ 、均匀磁化率 $\chi/L$ ）：

$T/J$	$E/L$	$C/L$	$\chi/L$
0.1	-1.417	0.034	0.0066
0.2	-1.407	0.138	0.0567
0.5	-1.333	0.401	0.1276
1.0	-1.064	0.559	0.1693
2.0	-0.653	0.278	0.1639
5.0	-0.274	0.055	0.1010

比热在 $T\approx J$ 附近有一个宽阔的最大值，而磁化率在 $T\approx1$ – $1.5\,J$ 附近出现峰值，随后向 $T=0$ 急剧下降——这是在 ED-02 中计算出的 Haldane 能隙 $\Delta\approx0.41\,J$ 所预期的指数激活行为 $\chi\sim e^{-\Delta/T}$ 的有限尺寸残留。

练习

将磁化率的计算重复一遍，但改为 L=9（这会花几分钟时间；如果你不耐烦，可以用 L=7），并与 L=8 的结果进行比较！注意所有从 parameter2xml 开始的步骤都必须重复。
粗略估计一下，在什么温度范围内可以将有限尺寸结果用作无限链的近似。

S=1/2 海森堡梯子

现在我们把模型改为一个双腿梯子，

H = J_0 \sum_{\text{legs}} \mathbf{S}_i \cdot \mathbf{S}_j + J_1 \sum_{\text{rungs}} \mathbf{S}_i \cdot \mathbf{S}_j ,

长度 L=6（即 6 个梯级，12 个格点），将参数 LATTICE 设为 “ladder”，耦合 J0 和 J1 都设为 1——这与 ED-03 中使用的是同一个格子：

      o---J0---o---J0---o---J0---o---J0---o---J0---o     leg 0
      |        |        |        |        |        |
      J1       J1       J1       J1       J1       J1
      |        |        |        |        |        |
      o---J0---o---J0---o---J0---o---J0---o---J0---o     leg 1

参数在 parm6b 中，Python 脚本在 tutorial6b.py 中。运行方式与上面的 parm6a 完全相同：

parameter2xml parm6b
fulldiag --write-xml parm6b.in.xml
fulldiag_evaluate --T_MIN 0.1 --T_MAX 10 --DELTA_T 0.1 parm6b.task1.out.xml

parm6b 除了用 LATTICE="ladder"、J0=1、J1=1 和 {L = 6} 取代链的参数之外，其余都与 parm6a 完全相同。

输出数据

独立对角化这个 L=6（开放腿边界）的梯子，得到有限尺寸单态-三重态能隙 $\Delta\approx0.68\,J$ ，略高于下方提示中给出的公认无限梯子值 $\approx 0.5\,J$ ——这提醒我们，对于带能隙的态而言，6 梯级系统离热力学极限仍相当远。相应的每格点热力学量为：

$T/J$	$E/L$	$C/L$	$\chi/L$
0.1	-0.550	0.013	0.0019
0.2	-0.544	0.126	0.0294
0.5	-0.451	0.424	0.1065
0.7	-0.368	0.382	0.1204
1.0	-0.274	0.252	0.1189
2.0	-0.137	0.072	0.0876

练习

讨论比热最大值所在的位置（提示：无限梯子的能隙约为 J/2）！
对 7 个梯级重复该计算（如果你不耐烦，可以用 5 个梯级）。
通过比较，你能对有限尺寸结果能够很好地近似无限系统的温度范围得出什么结论？

备注： 如果你熟悉量子蒙特卡罗（QMC）模拟，你会知道可以处理更大的系统，从而得到对热力学极限更好的近似。（试试看！你可能会发现，对于上面的例子，用 QMC 做得比全对角化更好其实并不那么容易。）

在某些条件下，精确对角化无疑仍然是首选方法。首先，如果系统本身就是有限的（并且可能很小），全对角化可以立即给出精确结果。其次，精确对角化不会遇到符号问题，因此也可以毫无附加限制地用于费米子模型或阻挫模型。这两个条件在第二部分将要讨论的磁性分子模型中同时满足。

磁性分子的热力学

现在我们将使用 fulldiag，通过计算完整能谱来模拟小型自旋团簇的精确热力学。

两个耦合的二聚体

建立自定义格子/图和模型

考虑如下这样一个由两个耦合二聚体组成的系统：“上"二聚体由两个自旋 $S_0$ 通过 $J_0$ 耦合而成，“下"二聚体由两个自旋 $S_1$ 同样通过 $J_0$ 耦合而成，而上二聚体的每个自旋都通过 $J_1$ 与下二聚体的每个自旋相耦合，

H = J_0\left(\mathbf{S}_1\cdot\mathbf{S}_2 + \mathbf{S}_3\cdot\mathbf{S}_4\right) + J_1\left(\mathbf{S}_1\cdot\mathbf{S}_3+\mathbf{S}_1\cdot\mathbf{S}_4+\mathbf{S}_2\cdot\mathbf{S}_3+\mathbf{S}_2\cdot\mathbf{S}_4\right) - h\sum_{i=1}^4 S^z_i .

        J0
   1 ---------- 2        (upper dimer, spin S0)
   |  \        /  |
   |   \  J1  /   |
   J1    \  /    J1
   |     /  \     |
   |   /      \   |
   3 ---------- 4        (lower dimer, spin S1)
        J0

这正是下面所编码的图：格点 1、2（类型 0）构成上二聚体，格点 3、4（类型 1）构成下二聚体，类型 0 的边是二聚体内部的键（ $J_0$ ），类型 1 的边是四条二聚体间的键（ $J_1$ ）。

首先，我们需要一个图来表示这个问题。这由文件 dd-graph.xml 中的以下条目定义：

<LATTICES>
<GRAPH name="double dimer" vertices="4">
<VERTEX id="1" type="0"></VERTEX>
<VERTEX id="2" type="0"></VERTEX>
<VERTEX id="3" type="1"></VERTEX>
<VERTEX id="4" type="1"></VERTEX>
<EDGE type="0" source="1" target="2"/>
<EDGE type="0" source="3" target="4"/>
<EDGE type="1" source="1" target="3"/>
<EDGE type="1" source="1" target="4"/>
<EDGE type="1" source="2" target="3"/>
<EDGE type="1" source="2" target="4"/>
</GRAPH>
</LATTICES>

注意：该文件必须使用不带 <?xml?> 声明或 <!DOCTYPE> 头部的纯 XML——ALPS 的解析器不支持这些。

接下来我们还想为类型 0 和类型 1 的边分别赋予不同的海森堡交换 J0 和 J1。这通过文件 model-dspin.xml 中的以下条目实现：

<MODELS>
<SITEBASIS name="spin">
  <PARAMETER name="local_spin" default="local_S"/>
  <PARAMETER name="local_S" default="1/2"/>
  <QUANTUMNUMBER name="S" min="local_spin" max="local_spin"/>
  <QUANTUMNUMBER name="Sz" min="-S" max="S"/>
  <OPERATOR name="Splus" matrixelement="sqrt(S*(S+1)-Sz*(Sz+1))">
    <CHANGE quantumnumber="Sz" change="1"/>
  </OPERATOR>
  <OPERATOR name="Sminus" matrixelement="sqrt(S*(S+1)-Sz*(Sz-1))">
    <CHANGE quantumnumber="Sz" change="-1"/>
  </OPERATOR>
  <OPERATOR name="Sz" matrixelement="Sz"/>
</SITEBASIS>
<BASIS name="spin">
  <SITEBASIS ref="spin">
    <PARAMETER name="local_spin" value="local_S#"/>
    <PARAMETER name="local_S#" value="1/2"/>
  </SITEBASIS>
  <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
</BASIS>
<HAMILTONIAN name="dimerized spin">
<PARAMETER name="J" default="1"/>
<PARAMETER name="h" default="0"/>
<BASIS ref="spin"/>
<SITETERM site="i">
<PARAMETER name="h#" default="h"/>
    -h#*Sz(i)
</SITETERM>
<BONDTERM source="i" target="j">
<PARAMETER name="J#" default="J"/>
    J#*Sz(i)*Sz(j)+J#/2*(Splus(i)*Sminus(j)+Sminus(i)*Splus(j))
</BONDTERM>
</HAMILTONIAN>
</MODELS>

注意：该文件必须使用不带 <?xml?> 声明或 <!DOCTYPE> 头部的纯 XML，并且必须包含 <SITEBASIS> 和 <BASIS> 定义——否则 fulldiag 无法构造希尔伯特空间。

请注意，我们实际上并不需要这个定义，因为默认的 models.xml 文件中 “spin” 哈密顿量已经包含了合适的定义。尽管如此，上面的例子说明了如何使用井号（#）自动为类型 n 的键赋予交换常数 Jn。

顺便说一下，我们还为格点 1、2 和 3、4 赋予了不同的类型。因此，我们能够通过分别指定 local_S0 和 local_S1 的值，为上二聚体和下二聚体赋予不同的局域自旋。

现在我们将计算局域自旋 $S_0=1$ （上二聚体）和 $S_1=1/2$ （下二聚体）、 $J_0=1$ 、 $J_1=0.4$ ，温度 $T=0.02$ 下的磁化曲线。

使用命令行

参数文件 parm6c 设置了该模拟：

LATTICE="double dimer"
MODEL="dimerized spin"
LATTICE_LIBRARY="dd-graph.xml"
MODEL_LIBRARY="model-dspin.xml"
local_S0=1
local_S1=1/2
J0      = 1
J1      = 0.4
h       = 0 
CONSERVED_QUANTUMNUMBERS="Sz"
{T = 0.02}

注意新出现的参数 LATTICE_LIBRARY 和 MODEL_LIBRARY，它们指向包含我们自定义格子和模型的文件。该计算通过以下命令序列执行：

parameter2xml parm6c
fulldiag --write-xml parm6c.in.xml
fulldiag_evaluate --H_MIN 0 --H_MAX 4 --DELTA_H 0.025 --versus h parm6c.task1.out.xml

请注意，这里我们使用带命令行参数的 fulldiag_evaluate 来指定磁场范围。特别是命令行参数 –versus h 会把磁场（而不是像前面例子中的温度）放在 x 轴上。

方法

由于只有 $(2\cdot1+1)^2\times(2\cdot\tfrac12+1)^2=36$ 个基矢态，这个系统对 fulldiag 而言微不足道——这个例子的重点不在于计算难度，而在于展示如何完全通过参数文件，搭建出真实磁性分子数据所需的那种完全自定义的格子图和哈密顿量。

输出数据

在这里所用的低温 $T=0.02\,J_0$ 下，随着 $h$ 的增大，基态的总自旋 $S_{\text{total}}$ 以整数步长增加，每一步都发生在下方解析能谱中两条能级交叉之处：

$h/J_0$	$\langle S^z_{\text{total}}\rangle$	基态所在子空间
0.0 – 0.8	0.00	$S_{\text{total}}=0$
1.0	0.67	交叉（ $S_{\text{total}}=0\to1$ ）
1.2	1.00	$S_{\text{total}}=1$
1.5	1.99	交叉（ $S_{\text{total}}=1\to2$ ）
1.8 – 2.2	2.00	$S_{\text{total}}=2$
2.5	2.99	交叉（ $S_{\text{total}}=2\to3$ ）
3.0 – 4.0	3.00	$S_{\text{total}}=3$

这三个交叉场直接由下面提示中给出的解析能量得出：令 $S_{\text{total}}=n$ （在 $S^z=n$ 处）的最低能量与 $S_{\text{total}}=n+1$ （在 $S^z=n+1$ 处）的最低能量相等，以 $J_0$ 为单位分别得到 $h_1=1.0$ 、 $h_2=1.4$ 和 $h_3=2.4$ ——这是一个在 $S^z_{\text{total}}=0,1,2,3$ 处出现平台的磁化台阶，只是被这个虽小但非零的温度稍稍磨圆了。

问题

绘制并解释这个结果！

提示： 两个耦合的 $S_0=1$ 和 $S_1=1/2$ 二聚体的能谱可以解析求出。这些能量为（存在一些简并）：

总自旋 $S_{\text{total}}=0$ ： $-11J_0/4$ 、 $-3J_0/4-2J_1$ ；
$S_{\text{total}}=1$ ： $-7J_0/4$ 、 $-3J_0/4-J_1$ 、 $5J_0/4-3J_1$ ；
$S_{\text{total}}=2$ ： $-3J_0/4+J_1$ 、 $J_0/4$ 、 $5J_0/4-J_1$ ；
$S_{\text{total}}=3$ ： $5J_0/4+2J_1$ 。

分子团簇 $V_{15}$

最后一个例子是分子纳米磁体 $\mathrm{V}_{15}$ 的模型，其中 15 个 $\mathrm{V}^{4+}$ 离子（各自 $S=1/2$ ）分布在夹着一个中心三角形的两个六边形上。它的低能磁性主要由两个六边形周围强烈的反铁磁交换，以及将它们与中心三角形相连的较弱交换所主导，因此在低温下这 15 个自旋的行为等效于一个阻挫的自旋 1/2 三角形：两个低能二重态加上一个四重态，与谱的其余部分明显分离（G. Chaboussant et al., Europhys. Lett. 59, 291 (2002)）。

        1---2                 outer hexagon (6 spins)
       /     \
      6       3
       \     /
        5---4
       /     \
      o       o                inner hexagon (6 spins),
       \     /                 rotated relative to the outer one
        o---o
          |
        triangle (3 spins)    weakly coupled to both hexagons

在模拟中，我们做一个简化假设：沿该图所有边的海森堡交换 $J$ 都相等。搭建自定义格子和模型文件的方式与上面耦合二聚体的例子完全相同——一个具有 15 个顶点的 <GRAPH>，以及相应描述六边形和三角形键的 <EDGE> 列表，再加上同样的 dimerized spin（或内置的 spin）哈密顿量。相关的图定义在 v15-graph.xml 中。运行模拟的方式与前面相同。参数在 parm6d 中，Python 脚本在 tutorial6d.py 中。

方法

15 个自旋 1/2 格点的完整希尔伯特空间维数为 $2^{15}=32768$ ；即使按 $S_z$ 子空间分解，fulldiag 也必须精确对角化一个包含数千个态的块，才能得到正确的低温热力学——这正是为什么，正如下面的说明所指出的，这次运行明显比前面的例子花费更长时间。

问题

你如何解释磁化率在低温下的行为？（提示：想想上面提到的有效低能自旋 1/2 三角形，以及它自身在磁场中的能级交叉会如何体现在 $\chi(T)$ 中。）

请注意，这个计算已经有些吃力，需要花上一段时间。因此，不妨先休息一下，再回来查看结果。

总结

全对角化给出了从平移不变的链和梯子，到完全自定义、带有阻挫的磁性分子等任意有限量子自旋系统的数值精确热力学——代价是希尔伯特空间随系统尺寸呈指数增长，这正是为什么即使利用了对称性，该方法也只能局限于几十个自旋的规模，而 sparsediag 仅针对基态和低能性质就能达到远大得多的系统尺寸。

附加练习

正方格子上的哈伯德模型

为正方格子上的哈伯德模型建立一个参数文件。
- 找到你想使用的格子。你可以使用正方格子，也可以使用链。注意边界条件！
- 在模型库中找到哈伯德模型。确保你理解它的各项。
- 打开对称性：守恒的是（x 和 y 方向的）动量以及粒子数。
- 选取可以检验的试探参数（例如 t=0 或 U=0），并运行模拟。
- 你会如何引入 t’？

ED-05 ED Phase Transition

ED-06 Full Diagonalization

一维自旋模型的热力学

S=1 海森堡链

参数

格子

方法

使用命令行

使用 Python

输出数据

练习

S=1/2 海森堡梯子

输出数据

练习

磁性分子的热力学

两个耦合的二聚体

建立自定义格子/图和模型

使用命令行

方法

输出数据

问题

分子团簇 V15V_{15}V15​

方法

问题

总结

附加练习

正方格子上的哈伯德模型

分子团簇 $V_{15}$