LM-02 自定义模型哈密顿量

定义自定义模型哈密顿量：自旋、费米子与玻色子

在教程 LM-01 中，你学习了如何定义自定义晶格。本教程将更进一步，介绍如何编写模型库 XML 文件，从基本构建块——局域 Hilbert 空间和相互作用项——出发，定义一个量子哈密顿量。这使你能够模拟任意格点模型，而不仅限于 ALPS 默认提供的少数几种。

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

描述四种 XML 构建块：SITEBASIS、QUANTUMNUMBER、OPERATOR 和 HAMILTONIAN。
从零开始构建带单离子各向异性的自旋-1 Heisenberg 链，并用 sparsediag 提取自旋能隙。
为无自旋费米子和软核玻色子编写类似的模型文件。
理解不同粒子类型之间矩阵元公式的差异。

ALPS 模型库 XML 格式

模型库文件是一个不含 <?xml ...?> 前导声明的裸 XML 文档，根元素为 <MODELS>。一个完整的模型由三个步骤构建：

第一步 — SITEBASIS：局域 Hilbert 空间

<SITEBASIS name="my basis">
  <PARAMETER name="P" default="1"/>              <!-- 可调参数，例如自旋 S -->
  <QUANTUMNUMBER name="Q" min="-P" max="P"
                 type="half_integer"/>           <!-- 守恒量子数           -->
  <OPERATOR name="Raise"
            matrixelement="sqrt(P*(P+1)-Q*(Q+1))">  <!-- 非对角升算符 -->
    <CHANGE quantumnumber="Q" change="1"/>
  </OPERATOR>
  <OPERATOR name="Lower"
            matrixelement="sqrt(P*(P+1)-Q*(Q-1))">
    <CHANGE quantumnumber="Q" change="-1"/>
  </OPERATOR>
  <OPERATOR name="Diag" matrixelement="Q"/>      <!-- 对角算符            -->
</SITEBASIS>

关键规则：

matrixelement 是一个代数表达式，在 <CHANGE> 作用之前的量子数值处求值（即在源态上求值）。
<CHANGE quantumnumber="Q" change="±1"/> 将 Q 改变指定量；ALPS 自动强制 min/max 边界，并在范围外将矩阵元置零。
type="half_integer" 允许半整数量子数值（自旋 $S=\tfrac{1}{2}$ 、 $\tfrac{3}{2}$ 、… 时需要）。

第二步 — BASIS：量子数约束

BASIS 是 SITEBASIS 与 HAMILTONIAN 之间的连接器，声明哪些量子数守恒以及每个扇区对应哪个模拟参数。

<BASIS name="my basis">
  <SITEBASIS ref="my basis"/>
  <CONSTRAINT quantumnumber="Q" value="Q_total"/>  <!-- 将总 Q 固定为参数 Q_total -->
</BASIS>

第三步 — HAMILTONIAN：格点项与键项

<HAMILTONIAN name="my model">
  <PARAMETER name="J" default="1"/>
  <BASIS ref="my basis"/>             <!-- 引用 BASIS，而非 SITEBASIS -->
  <SITETERM site="i">
    D*Diag(i)*Diag(i)                <!-- 格点表达式                 -->
  </SITETERM>
  <BONDTERM source="i" target="j">
    J*Diag(i)*Diag(j)                <!-- 键表达式                   -->
  </BONDTERM>
</HAMILTONIAN>

常见错误： 在 <HAMILTONIAN> 内写 <SITEBASIS ref="..."/> 会导致 ALPS 报错 “unexpected element”。请务必在 <HAMILTONIAN> 内使用 <BASIS ref="..."/>，并单独定义包裹 <SITEBASIS> 的 <BASIS> 元素。

完整文件的结构为：

<MODELS>
  <!-- SITEBASIS 元素（局域 Hilbert 空间） -->
  <!-- BASIS 元素（约束声明）             -->
  <!-- HAMILTONIAN 元素（相互作用项）     -->
</MODELS>

在参数文件中通过以下方式引用：

MODEL_LIBRARY = "my_models.xml"
MODEL         = "my model"

示例一 — 自旋：带单离子各向异性的 S=1 链

物理背景

哈密顿量为

H = J \sum_{\langle i,j \rangle} \mathbf{S}_i \cdot \mathbf{S}_j + D \sum_i \bigl(S^z_i\bigr)^2

其中 $J > 0$ 为反铁磁交换耦合， $D$ 为单离子各向异性。用升降算符改写 $\mathbf{S}_i \cdot \mathbf{S}_j$ ：

\mathbf{S}_i \cdot \mathbf{S}_j = S^z_i S^z_j + \tfrac{1}{2}\bigl(S^+_i S^-_j + S^-_i S^+_j\bigr).

在 $D = 0$ 时，该链处于Haldane 相：一个具有体自旋能隙 $\Delta \approx 0.41\,J$ 的有能隙对称保护拓扑相（Haldane 1983；White & Huse 1993）。随着 $D$ 增大，能隙在量子临界点 $D_c \approx 0.97\,J$ 处关闭，系统进入大- $D$ 相（其基态接近乘积态 $|0, 0, \ldots, 0\rangle$ ，即所有格点均处于 $S^z = 0$ 态）。

哈密顿量骨架（开放链，4 个格点示意）：

  D       D       D       D
  o --J-- o --J-- o --J-- o

D  = 单离子各向异性  D*(Sz)^2      [SITETERM]
J  = Heisenberg 交换              [BONDTERM]
h  = 均匀外场       -h*Sz         [SITETERM，默认为 0]

晶格选择

本示例使用 ALPS 内置的 chain lattice（周期边界条件），链长 $L=8$ 。周期边界条件下，sparsediag 可同时利用平移对称性按总动量 $k$ 进行块对角化，从而更清晰地呈现体相能隙。关于所有可用晶格几何，参见 ALPS 晶格库。

方法选择

sparsediag（稀疏 Lanczos 精确对角化）给出精确本征值，适用于 Hilbert 空间维数 $\lesssim 10^6$ 的体系。对于 $L=8$ 、 $S=1$ 、 $S^z_{\text{tot}}=0$ 扇区，基矢数仅约 1.3 万——Lanczos 方法可在几秒内完成。本例采用精确对角化的原因在于：我们需要对四位有效数字精度的能量差（能隙）随 $D$ 的变化规律进行分析。

第一步：编写模型库文件

新建 my_models.xml：

<MODELS>

  <!-- ===========================================================
       格点基：通用自旋-S
         * local_S 默认为 1（S=1 链）；设置 local_S=0.5 可得 S=1/2 链
         * 量子数：S（固定）和 Sz（从 -S 到 S）
         * 算符：Splus、Sminus（升降）、Sz（对角）
  ============================================================ -->
  <SITEBASIS name="my spin">
    <PARAMETER name="local_S" default="1"/>

    <QUANTUMNUMBER name="S"  min="local_S" max="local_S" type="half_integer"/>
    <QUANTUMNUMBER name="Sz" min="-S" max="S" type="half_integer"/>

    <!-- S^+ 将 Sz 升高 1；矩阵元由 Wigner-Eckart 定理给出 -->
    <OPERATOR name="Splus" matrixelement="sqrt(S*(S+1)-Sz*(Sz+1))">
      <CHANGE quantumnumber="Sz" change="1"/>
    </OPERATOR>

    <!-- S^- 将 Sz 降低 1 -->
    <OPERATOR name="Sminus" matrixelement="sqrt(S*(S+1)-Sz*(Sz-1))">
      <CHANGE quantumnumber="Sz" change="-1"/>
    </OPERATOR>

    <!-- Sz 为对角算符：矩阵元等于 Sz 的当前值 -->
    <OPERATOR name="Sz" matrixelement="Sz"/>
  </SITEBASIS>

  <!-- BASIS 将格点基与总-Sz 约束关联 -->
  <BASIS name="anisotropic spin chain">
    <SITEBASIS ref="my spin"/>
    <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
  </BASIS>

  <!-- ===========================================================
       哈密顿量：Heisenberg 交换 + 单离子各向异性
         H = J * sum_{<ij>} (Sz(i)*Sz(j) + 0.5*(S+(i)*S-(j) + S-(i)*S+(j)))
           + D * sum_i  Sz(i)^2
           - h * sum_i  Sz(i)      [可选均匀外场]
  ============================================================ -->
  <HAMILTONIAN name="anisotropic spin chain">
    <PARAMETER name="J" default="1"/>
    <PARAMETER name="D" default="0"/>
    <PARAMETER name="h" default="0"/>

    <BASIS ref="anisotropic spin chain"/>

    <!-- 格点项：单离子各向异性 + 外场 -->
    <SITETERM site="i">
      D*Sz(i)*Sz(i) - h*Sz(i)
    </SITETERM>

    <!-- 键项：以 S+/S-/Sz 写出的 Heisenberg 交换 -->
    <BONDTERM source="i" target="j">
      J*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j) + Sminus(i)*Splus(j)))
    </BONDTERM>
  </HAMILTONIAN>

</MODELS>

各元素的作用

<SITEBASIS name="my spin"> 定义局域 Hilbert 空间。自旋- $S$ 格点有 $2S+1$ 个态，以 $S^z \in \{-S,\, -S+1,\, \ldots,\, S\}$ 标记。默认 local_S=1，每个格点有三个态： $|-1\rangle$ 、 $|0\rangle$ 、 $|{+}1\rangle$ 。

<PARAMETER name="local_S"> 设置自旋大小。 ALPS 要求自旋量子数（S）与可调参数名称不同，以避免运行时名称冲突；遵循 ALPS 惯例，可调参数命名为 local_S，量子数命名为 S。

<QUANTUMNUMBER name="S" min="local_S" max="local_S" ...> 将 $S$ 固定为 local_S 参数的值，使每个格点携带相同的总自旋。 <QUANTUMNUMBER name="Sz" ...> 将 $S^z$ 注册为流动量子数，ALPS 可用其对哈密顿量进行块对角化。

<BASIS name="..."> 包裹格点基并添加 <CONSTRAINT>，将总 $S^z$ 固定为模拟参数 Sz_total。在 <HAMILTONIAN> 内写 <BASIS ref="..."/>（而非 <SITEBASIS ref="..."/>）；若使用错误标签，ALPS 将报错。

<OPERATOR name="Splus" matrixelement="..."> 定义升算符 $S^+$ 。 matrixelement 给出 $\langle S,\, S^z{+}1 \mid S^+ \mid S,\, S^z \rangle = \sqrt{S(S+1)-S^z(S^z+1)}$ ，在 $S^z$ 的源值处求值。子元素 <CHANGE quantumnumber="Sz" change="1"/> 告知 ALPS $S^+$ 将 $S^z$ 升高 1；当 $S^z = S$ 时，ALPS 自动将矩阵元置零。

<SITETERM site="i"> 指定格点哈密顿量。表达式 D*Sz(i)*Sz(i) 在格点 $i$ 上两次计算对角算符 $S^z_i$ ，得到单离子各向异性 $D(S^z_i)^2$ 。

<BONDTERM source="i" target="j"> 指定双格点键相互作用。乘积 Splus(i)*Sminus(j) 表示 $S^+_i \otimes S^-_j$ ；ALPS 自动构建张量积矩阵。系数 $\tfrac{1}{2}$ 来自 $S^+S^- = 2(S^xS^x + S^yS^y)$ 。

第二步：编写参数文件

要计算自旋能隙 $\Delta = E_1(S^z_{\text{tot}}=1) - E_0(S^z_{\text{tot}}=0)$ ，需要在每个 $S^z_{\text{tot}}$ 扇区分别运行 sparsediag。新建 parm_spin1：

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "chain lattice"
MODEL                    = "anisotropic spin chain"
L                        = 8
J                        = 1
NUMBER_EIGENVALUES       = 1
CONSERVED_QUANTUMNUMBERS = "Sz"

{D=0.0; Sz_total=0;}
{D=0.0; Sz_total=1;}
{D=1.0; Sz_total=0;}
{D=1.0; Sz_total=1;}
{D=2.0; Sz_total=0;}
{D=2.0; Sz_total=1;}

要点说明：

LATTICE="chain lattice" — ALPS 标准周期链；无需 LATTICE_LIBRARY。周期边界条件下，sparsediag 还利用平移对称性（逐个全动量扇区对角化），能更清晰地呈现体相 Haldane 能隙。
CONSERVED_QUANTUMNUMBERS="Sz" 激活块对角化；ALPS 只对角化 Sz_total 等于所请求值的块。
共创建六个任务：在三个 $D$ 值下各两个扇区（ $S^z_{\text{tot}} = 0,\, 1$ ）。
NUMBER_EIGENVALUES=1 仅请求每个动量块中的最低本征值。

第三步：运行模拟

parameter2xml parm_spin1
sparsediag --write-xml parm_spin1.in.xml

sparsediag 依次运行六个任务，将结果写入 parm_spin1.task1.out.xml 到 parm_spin1.task6.out.xml。

第四步：提取自旋能隙

import xml.etree.ElementTree as ET, glob

# 周期边界条件下，sparsediag 对每个全动量扇区分别对角化，
# 每个扇区写一个 <EIGENVALUES> 块。
# 取各扇区最低本征值的最小值即为各 Sz 扇区的基态能量。
results = {}
for fname in sorted(glob.glob('parm_spin1.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    D  = float(params['D'])
    Sz = int(params['Sz_total'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results.setdefault(D, {})[Sz] = min(means)

print(f"{'D/J':>5}  {'E(Sz=0)':>12}  {'E(Sz=1)':>12}  {'能隙 Δ/J':>9}  相")
print('─' * 62)
for D in sorted(results):
    E0 = results[D][0]
    E1 = results[D][1]
    gap = E1 - E0
    phase = 'Haldane 相' if D < 0.97 else '大-D 相'
    print(f"{D:5.1f}  {E0:12.5f}  {E1:12.5f}  {gap:9.5f}  {phase}")

预期输出（L=8，周期链， $J=1$ ）：

  D/J      E(Sz=0)      E(Sz=1)    能隙 Δ/J  相
──────────────────────────────────────────────────────────────
  0.0   -11.33696   -10.74340    0.59356  Haldane 相
  1.0    -7.00947    -6.62505    0.38442  大-D 相
  2.0    -4.33528    -3.66645    0.66883  大-D 相

在 $D = 0$ 时， $L=8$ 下能隙 $\Delta \approx 0.59\,J$ ；体相值为 $0.41\,J$ ，有限尺寸修正使其偏大（见延伸问题 1）。在 $D \approx J$ 附近，能隙缩小至 $0.38\,J$ ——量子相变（ $D_c \approx 0.97\,J$ ）的有限尺寸前兆。在 $D = 2J$ 时，系统已进入大- $D$ 相，能隙增大至 $0.67\,J$ 。

示例二 — 费米子：带近邻排斥的无自旋费米子链

无自旋费米子是最简单的费米子模型：每个格点要么空（ $|0\rangle$ ）要么被占据（ $|1\rangle$ ），没有自旋自由度。哈密顿量为

H = -t \sum_{\langle i,j \rangle} \bigl(c^\dagger_i c_j + c^\dagger_j c_i\bigr) + V \sum_{\langle i,j \rangle} n_i n_j - \mu \sum_i n_i,

描述自由跳跃能带（由 $t$ 控制）、近邻密度-密度排斥 $V$ 以及化学势 $\mu$ 。当 $V > 2t$ 时，基态为电荷密度波；当 $V < 2t$ 时，为 Luttinger 液体（金属态）。

哈密顿量骨架（开放链，4 个格点示意）：

  -μ        -μ        -μ        -μ
   o ------- o ------- o ------- o
      -t, V      -t, V      -t, V

-μ = 化学势    -μ*n           [SITETERM]
-t = 跳跃振幅  -t*(c†c+h.c.)  [BONDTERM]
 V = 密度排斥   V*n*n          [BONDTERM]

将以下内容添加到 my_models.xml（在 <MODELS> 元素内）：

  <!-- ===========================================================
       格点基：无自旋费米子
         * 占据数 N ∈ {0, 1}（由 max="1" 实现泡利不相容）
         * type="fermion" 告知 ALPS 对一维系统自动应用 Jordan-Wigner 串
  ============================================================ -->
  <SITEBASIS name="my spinless fermion" type="fermion">
    <QUANTUMNUMBER name="N" min="0" max="1"/>

    <!-- c† 产生一个费米子：⟨1|c†|0⟩ = 1 -->
    <OPERATOR name="cdag" matrixelement="1">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>

    <!-- c 湮灭一个费米子：⟨0|c|1⟩ = 1 -->
    <OPERATOR name="c" matrixelement="1">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>

    <!-- 粒子数算符：对角，矩阵元等于当前 N -->
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="spinless fermion chain">
    <SITEBASIS ref="my spinless fermion"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <!-- ===========================================================
       哈密顿量：无自旋费米子链
         H = -t*(cdag(i)*c(j)+h.c.) + V*n(i)*n(j) - mu*n(i)
  ============================================================ -->
  <HAMILTONIAN name="spinless fermion chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="V"  default="0"/>
    <PARAMETER name="mu" default="0"/>

    <BASIS ref="spinless fermion chain"/>

    <SITETERM site="i">
      -mu*n(i)
    </SITETERM>

    <BONDTERM source="i" target="j">
      -t*(cdag(i)*c(j) + cdag(j)*c(i)) + V*n(i)*n(j)
    </BONDTERM>
  </HAMILTONIAN>

为何 `type="fermion"` 至关重要

费米子算符满足反对易关系： $\{c_i, c^\dagger_j\} = \delta_{ij}$ 。在一维中，这由 Jordan-Wigner 变换精确处理，它将费米子映射为带有非局域符号串的硬核玻色子。在 SITEBASIS 上设置 type="fermion" 指示 ALPS 在费米子算符跨越中间格点跳跃时自动插入该符号串。若不设置，符号将会出错，谱也将不正确。

晶格选择

本示例使用 open chain lattice（开放边界条件），链长 $L=10$ 。开放边界条件避免了相互作用费米子在动量扇区对角化时出现的符号复杂性，是一维链 ED 的标准选择。关于其他链变体，参见 ALPS 晶格库。

方法选择

sparsediag 在此精确高效：对于 $L=10$ 、 $N=5$ ，Hilbert 空间维数仅为 $\binom{10}{5}=252$ ——几乎即时完成。设置 NUMBER_EIGENVALUES=4 可同时获得基态和最低中性激发态。

参数文件

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "open chain lattice"
MODEL                    = "spinless fermion chain"
L                        = 10
t                        = 1
mu                       = 0
NUMBER_EIGENVALUES       = 4
CONSERVED_QUANTUMNUMBERS = "N"

{V=0.0; N_total=5;}
{V=1.0; N_total=5;}
{V=2.0; N_total=5;}
{V=3.0; N_total=5;}

参数	含义	使用值
`L`	链长（格点数）	10
`t`	跳跃振幅	1
`mu`	化学势	0（半填充）
`V`	近邻排斥	0–3（扫描）
`N_total`	总粒子数	5（半填充）

N_total=5 为半填充（每两个格点一个费米子），此时电荷密度波不稳定性最强。

运行模拟

parameter2xml parm_fermion
sparsediag --write-xml parm_fermion.in.xml

提取并显示结果

import xml.etree.ElementTree as ET, glob

results = {}
for fname in sorted(glob.glob('parm_fermion.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    V  = float(params['V'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results[V] = sorted(means)

print(f"{'V/t':>5}  {'E₀':>10}  {'E₁':>10}  {'E₁−E₀':>8}")
print('─' * 42)
for V in sorted(results):
    Es = results[V]
    print(f"{V:5.1f}  {Es[0]:10.5f}  {Es[1]:10.5f}  {Es[1]-Es[0]:8.5f}")

预期输出（L=10，开放链， $t=1$ ， $N=5$ ）：

  V/t        E₀          E₁       E₁−E₀
──────────────────────────────────────────
  0.0   -6.02667    -5.45742    0.56926
  1.0   -5.00000    -4.30391    0.69609
  2.0   -4.26419    -3.50162    0.76257
  3.0   -3.74967    -2.97936    0.77031

小结

半填充扇区内的激发能隙随 $V$ 单调增大：排斥势使电荷扇区变硬，提高了中性（粒子-空穴）激发的能量代价。要测量真正的电荷能隙 $\Delta_c = E_0(N+1) + E_0(N-1) - 2E_0(N)$ （在 $V < 2t$ 时为零，在 $V > 2t$ 时打开），需在参数文件中增加 N_total=4 和 N_total=6 的任务（参见延伸问题 3）。

示例三 — 玻色子：软核 Bose-Hubbard 链

Bose-Hubbard 模型允许每个格点存在多个玻色子（最多 Nmax 个），描述光学晶格实验中的超流–Mott 绝缘体转变：

H = -t \sum_{\langle i,j \rangle} \bigl(b^\dagger_i b_j + b^\dagger_j b_i\bigr) + \frac{U}{2} \sum_i n_i(n_i - 1) - \mu \sum_i n_i.

格点排斥 $U$ 惩罚双占据（及更高占据），在 $U \gg t$ 且公度填充时，基态为具有粒子-空穴能隙的 Mott 绝缘体。

哈密顿量骨架（开放链，4 个格点示意）：

   U,-μ       U,-μ       U,-μ       U,-μ
    o -------- o -------- o -------- o
          -t          -t          -t

 U  = 格点排斥    (U/2)*n*(n-1)  [SITETERM]
-μ  = 化学势      -μ*n           [SITETERM]
-t  = 跳跃振幅    -t*(b†b+h.c.)  [BONDTERM]

将以下内容添加到 my_models.xml：

  <!-- ===========================================================
       格点基：软核玻色子
         * 占据数 N ∈ {0, 1, …, Nmax}（默认 Nmax=3）
         * b† 和 b 带有来自玻色对易关系 [b, b†] = 1 的 sqrt 因子
  ============================================================ -->
  <SITEBASIS name="my boson">
    <PARAMETER name="Nmax" default="3"/>

    <QUANTUMNUMBER name="N" min="0" max="Nmax"/>

    <!-- b† 将 N 升高 1；矩阵元 ⟨N+1|b†|N⟩ = sqrt(N+1) -->
    <OPERATOR name="bdag" matrixelement="sqrt(N+1)">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>

    <!-- b 将 N 降低 1；矩阵元 ⟨N-1|b|N⟩ = sqrt(N) -->
    <OPERATOR name="b" matrixelement="sqrt(N)">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>

    <!-- 粒子数算符 -->
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="Bose-Hubbard chain">
    <SITEBASIS ref="my boson"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <!-- ===========================================================
       哈密顿量：软核 Bose-Hubbard 链
         H = -t*(bdag(i)*b(j)+h.c.) + U/2*n(i)*(n(i)-1) - mu*n(i)
  ============================================================ -->
  <HAMILTONIAN name="Bose-Hubbard chain">
    <PARAMETER name="t"   default="1"/>
    <PARAMETER name="U"   default="0"/>
    <PARAMETER name="mu"  default="0"/>

    <BASIS ref="Bose-Hubbard chain"/>

    <SITETERM site="i">
      U/2*n(i)*(n(i)-1) - mu*n(i)
    </SITETERM>

    <BONDTERM source="i" target="j">
      -t*(bdag(i)*b(j) + bdag(j)*b(i))
    </BONDTERM>
  </HAMILTONIAN>

矩阵元对比：玻色子、费米子与自旋

粒子类型	产生算符	湮灭算符
自旋（ $S^+$ ）	$\sqrt{S(S+1)-S^z(S^z+1)}$	$\sqrt{S(S+1)-S^z(S^z-1)}$
软核玻色子（ $b^\dagger$ ）	$\sqrt{N+1}$	$\sqrt{N}$
硬核玻色子 / 无自旋费米子（ $c^\dagger$ ）	$1$	$1$

软核玻色子中的 $\sqrt{N+1}$ 因子反映了向已占据模式的受激发射增强——玻色-爱因斯坦统计的精髓。硬核玻色子（占据数上限为 1）和无自旋费米子共享矩阵元 $= 1$ ；其区别在于费米子的反对易关系，通过 SITEBASIS 上的 type="fermion" 编码。

晶格选择

本示例使用 open chain lattice，链长 $L=6$ 。对于 $N_{\max}=3$ 、 $L=6$ 的体系， $N=6$ 扇区的 Hilbert 空间维数可控，适合 Lanczos ED。关于用于研究二维超流-Mott 相图的正方晶格等其他几何，参见 ALPS 晶格库。

方法选择

sparsediag 在任意 $U/t$ 下均给出精确结果，可直接获得固定 $N$ 扇区中的最低本征值。量子蒙特卡洛（worm）更适用于较大体系，但 ED 是建立基准值并在小 $L$ 下进行定量比较的理想选择。

参数文件

MODEL_LIBRARY            = "my_models.xml"
LATTICE                  = "open chain lattice"
MODEL                    = "Bose-Hubbard chain"
Nmax                     = 3
L                        = 6
t                        = 1
NUMBER_EIGENVALUES       = 2
CONSERVED_QUANTUMNUMBERS = "N"

{U=0.0;  mu=0; N_total=6;}
{U=2.0;  mu=0; N_total=6;}
{U=5.0;  mu=0; N_total=6;}
{U=10.0; mu=0; N_total=6;}

参数	含义	使用值
`L`	链长（格点数）	6
`Nmax`	每格点最大玻色子数	3
`t`	跳跃振幅	1
`U`	格点排斥	0–10（扫描）
`mu`	化学势	0
`N_total`	总玻色子数	6（单位填充）

N_total=6 给出单位填充（每个格点一个玻色子），此时 Mott 能隙最大。

运行模拟

parameter2xml parm_boson
sparsediag --write-xml parm_boson.in.xml

提取并显示结果

import xml.etree.ElementTree as ET, glob

results = {}
for fname in sorted(glob.glob('parm_boson.task*.out.xml')):
    tree = ET.parse(fname)
    root = tree.getroot()
    params = {p.get('name'): p.text.strip() for p in root.findall('.//PARAMETER')}
    U  = float(params['U'])
    means = [float(m.text)
             for m in root.findall('.//SCALAR_AVERAGE[@name="Energy"]/MEAN')]
    if means:
        results[U] = sorted(means)

print(f"{'U/t':>5}  {'E₀(N=6)':>12}  {'E₁(N=6)':>12}  {'E₁−E₀':>8}")
print('─' * 48)
for U in sorted(results):
    Es = results[U]
    print(f"{U:5.1f}  {Es[0]:12.5f}  {Es[1]:12.5f}  {Es[1]-Es[0]:8.5f}")

预期输出（L=6，开放链， $t=1$ ， $N_{\max}=3$ ， $N=6$ ）：

  U/t      E₀(N=6)      E₁(N=6)     E₁−E₀
────────────────────────────────────────────────
  0.0   -10.26758     -9.55861     0.70897
  2.0    -6.64995     -5.60418     1.04577
  5.0    -3.71673     -1.39313     2.32360
 10.0    -1.96696      4.35707     6.32403

小结

单位填充扇区内的激发能隙随 $U$ 迅速增大：在 $U=0$ （超流态）时，能隙仅由动能决定（ $\approx 0.71\,t$ ）；当 $U=10\,t$ 时，格点排斥使任何粒子-空穴激发极为代价高昂（ $\approx 6.3\,t$ ），标志着深度 Mott 绝缘体特征。要测量真正的单粒子能隙 $\Delta = E_0(N+1) + E_0(N-1) - 2E_0(N)$ （在超流相中为零，在 Mott 相中打开），可在参数文件中增加 N_total=5 和 N_total=7 的任务（参见延伸问题 4）。

完整的 `my_models.xml`

将三种格点基和哈密顿量合并到一个文件：

<MODELS>

  <!-- ===== 自旋 ===== -->
  <SITEBASIS name="my spin">
    <PARAMETER name="local_S" default="1"/>
    <QUANTUMNUMBER name="S"  min="local_S" max="local_S" type="half_integer"/>
    <QUANTUMNUMBER name="Sz" min="-S" max="S" type="half_integer"/>
    <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="anisotropic spin chain">
    <SITEBASIS ref="my spin"/>
    <CONSTRAINT quantumnumber="Sz" value="Sz_total"/>
  </BASIS>

  <HAMILTONIAN name="anisotropic spin chain">
    <PARAMETER name="J" default="1"/>
    <PARAMETER name="D" default="0"/>
    <PARAMETER name="h" default="0"/>
    <BASIS ref="anisotropic spin chain"/>
    <SITETERM site="i">D*Sz(i)*Sz(i) - h*Sz(i)</SITETERM>
    <BONDTERM source="i" target="j">
      J*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j) + Sminus(i)*Splus(j)))
    </BONDTERM>
  </HAMILTONIAN>


  <!-- ===== 无自旋费米子 ===== -->
  <SITEBASIS name="my spinless fermion" type="fermion">
    <QUANTUMNUMBER name="N" min="0" max="1"/>
    <OPERATOR name="cdag" matrixelement="1">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>
    <OPERATOR name="c" matrixelement="1">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="spinless fermion chain">
    <SITEBASIS ref="my spinless fermion"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <HAMILTONIAN name="spinless fermion chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="V"  default="0"/>
    <PARAMETER name="mu" default="0"/>
    <BASIS ref="spinless fermion chain"/>
    <SITETERM site="i">-mu*n(i)</SITETERM>
    <BONDTERM source="i" target="j">
      -t*(cdag(i)*c(j) + cdag(j)*c(i)) + V*n(i)*n(j)
    </BONDTERM>
  </HAMILTONIAN>


  <!-- ===== 软核玻色子 ===== -->
  <SITEBASIS name="my boson">
    <PARAMETER name="Nmax" default="3"/>
    <QUANTUMNUMBER name="N" min="0" max="Nmax"/>
    <OPERATOR name="bdag" matrixelement="sqrt(N+1)">
      <CHANGE quantumnumber="N" change="1"/>
    </OPERATOR>
    <OPERATOR name="b" matrixelement="sqrt(N)">
      <CHANGE quantumnumber="N" change="-1"/>
    </OPERATOR>
    <OPERATOR name="n" matrixelement="N"/>
  </SITEBASIS>

  <BASIS name="Bose-Hubbard chain">
    <SITEBASIS ref="my boson"/>
    <CONSTRAINT quantumnumber="N" value="N_total"/>
  </BASIS>

  <HAMILTONIAN name="Bose-Hubbard chain">
    <PARAMETER name="t"  default="1"/>
    <PARAMETER name="U"  default="0"/>
    <PARAMETER name="mu" default="0"/>
    <BASIS ref="Bose-Hubbard chain"/>
    <SITETERM site="i">U/2*n(i)*(n(i)-1) - mu*n(i)</SITETERM>
    <BONDTERM source="i" target="j">
      -t*(bdag(i)*b(j) + bdag(j)*b(i))
    </BONDTERM>
  </HAMILTONIAN>

</MODELS>

思考与延伸

能隙与系统尺寸的关系：在 $D=0$ 时，对 $L = 6, 8, 10, 12$ 运行自旋-1 示例。画出 $\Delta(L)$ 关于 $1/L$ 的图。能隙是否外推到体相 Haldane 值 $\Delta_\infty \approx 0.41\,J$ ？开放与周期边界条件的有限尺寸修正有何不同？
自旋- $\tfrac{1}{2}$ 与自旋-1 的比较：在 parm_spin1 中添加 local_S=0.5（SITEBASIS 使用 local_S 作为可调参数）。 $S = \tfrac{1}{2}$ 的 Heisenberg 模型是无能隙的（Luttinger 液体）。验证随 $L \to \infty$ 能隙外推到零，这与 $S = 1$ 情况形成对比。
无自旋费米子中的电荷密度波能隙：用 $V = 0, 1, 2, 3$ 运行 parm_fermion。通过添加 N_total=4 和 N_total=6 的任务，提取电荷能隙 $\Delta_c = E_0(N_{\text{tot}}+1) + E_0(N_{\text{tot}}-1) - 2E_0(N_{\text{tot}})$ 。在 $V/t$ 的哪个值处电荷密度波能隙打开？
玻色子中的 Mott 能隙：计算单位填充下 Bose-Hubbard 链的单粒子能隙随 $U/t$ 的变化。与无自旋费米子的电荷密度波能隙比较——两者都是电荷能隙，但物理起源不同（分别由排斥驱动和动能驱动）。
自定义晶格 + 自定义模型：将 my_models.xml 与 LM-01 中的蜂窝晶格结合。对小集团使用 sparsediag 在蜂窝晶格上运行自旋-1 Heisenberg 模型。单离子各向异性 $D$ 是否在阻挫蜂窝晶格上稳定出与链上不同的相？
添加次近邻键：向 anisotropic spin chain 哈密顿量中添加第二个 type="1" 的 <BONDTERM>：
```
<BONDTERM type="1" source="i" target="j">
  J2*(Sz(i)*Sz(j) + 0.5*(Splus(i)*Sminus(j)+Sminus(i)*Splus(j)))
</BONDTERM>
```
使用 nnn chain lattice，研究次近邻耦合 $J_2$ 对 Haldane 能隙的影响。参见晶格操作指南以了解该晶格中棱类型 1 的定义方式。

LM-01 蜂窝反铁磁体

LM-02 自定义模型哈密顿量

定义自定义模型哈密顿量：自旋、费米子与玻色子

ALPS 模型库 XML 格式

第一步 — SITEBASIS：局域 Hilbert 空间

第二步 — BASIS：量子数约束

第三步 — HAMILTONIAN：格点项与键项

示例一 — 自旋：带单离子各向异性的 S=1 链

物理背景

晶格选择

方法选择

第一步：编写模型库文件

各元素的作用

第二步：编写参数文件

第三步：运行模拟

第四步：提取自旋能隙

示例二 — 费米子：带近邻排斥的无自旋费米子链

为何 type="fermion" 至关重要

晶格选择

方法选择

参数文件

运行模拟

提取并显示结果

小结

示例三 — 玻色子：软核 Bose-Hubbard 链

矩阵元对比：玻色子、费米子与自旋

晶格选择

方法选择

参数文件

运行模拟

提取并显示结果

小结

完整的 my_models.xml

思考与延伸

为何 `type="fermion"` 至关重要

完整的 `my_models.xml`