第11章 · 辅助场量子蒙特卡洛#

11.1 引言#

哈特里-福克方法（或更一般地，任何平均场理论）最重要的优势之一是，可以通过简单的线性代数运算（可能包括小矩阵的对角化）得到完整且高效的数值解。然而，如第1章所讨论，这些基于独立电子波函数的技术通常不足以准确描述关联系统，例如哈伯德模型。在这方面，人们引入了投影方法，通过迭代应用合适的算子来滤除基态；第8章和第9章讨论了两个例子。在此，我们想介绍另一种投影技术，它特别适用于费米子哈伯德模型（Sorella et al., 1989），即所谓的辅助场量子蒙特卡洛（auxiliary-field quantum Monte Carlo, AFQMC）。该方法的主要任务是将多体问题简化为多个类似平均场的计算，这些计算可以通过在斯莱特行列式空间中使用随机方法来完成。这样，只要在所述变换过程中不出现符号问题（sign problem），统计误差就可以得到控制；值得注意的是，在许多有趣的情况下这确实是可能的。此外，这里的符号问题与第10.1节讨论的玻色-费米子不稳定性无关，因为在AFQMC中明确采样了反对称的斯莱特行列式，因此它通常比常规情况（例如在构型空间中定义的格林函数蒙特卡洛（Fahy and Hamann, 1991））更温和。

总体思路是通过应用虚时τ传播（见第1.7节）来滤除基态波函数：

与之前一样，必要条件是选择一个初始态 $| \Psi_{0} \rangle$，使得 $\langle \Upsilon_{0} | \Psi_{0} \rangle \neq 0$。首先，如第5.6节所做，我们将自旋指标σ和格点指标i压缩为单个指标I，其取值范围从1到2L：

然后，我们将专门考虑哈密顿量包含动能项和势能项的情况：

其中

这里，我们用以下符号表示格点 $i = 1 , \ldots , L$ 上自旋向上和向下电子的密度：

动能项用一个 $2 L \times 2 L$ 的矩阵 $\mathbf{K}$ 表示；在最简单的情况下，如果 i 和 $j$ 是最近邻，则 $K_{i , j} = K_{i + L , j + L} = - t$，否则为0。在粒子数守恒的一般情况下，可以将其视为一般的厄米矩阵。此外，在该框架下，也可以扩展到费米子场中任意双线性形式的情况，即不假设动能与粒子数对易。

在下文中，我们将初始态视为具有 $N_{e}$ 个电子的非相互作用波函数：

这一般情形由 $2 L \times N_{e}$ 矩阵 $\psi$ 描述，其中由 $\{\psi_{I , \alpha} \}$ 定义的轨道不具有确定的 自旋值（这包括类 BCS 态，前提是已对自旋向下电子进行过粒子-空穴变换）。相比之下，只要自旋向上分量和自旋向下分量之间没有混合，矩阵 $\psi$ 就可以写成两个块的形式，维度分别为 $L \times N_{e}^{\uparrow}$ 和 $L \times N_{e}^{\downarrow}$。

在没有相互作用项的情况下，方程 (11.1) 的投影技术可以用简单而高效的方式实现。实际上，单粒子传播子 exp $\textstyle ( - \tau{\boldsymbol{\mathcal{K}}} )$ 将一个无关联波函数变换为另一个仍然无关联的态：

其中演化后的轨道由下式给出：

这里，$\{\epsilon_{\beta} \}$ 是动能项 $\kappa$ 的本征值，矩阵 U 包含其本征矢量，即 $\begin{array} {r} {\phi_{\beta}^{\dagger} = \sum_{I} U_{I , \beta} d_{I}^{\dagger}} \end{array}$ 。因此，一个斯莱特行列式的虚时演化可以方便地写成原始 $2 L \times N_{e}$ 矩阵 $\psi$ 与 $2 L \times 2 L$ 传播矩阵 $\mathbf P ( \tau )$ 的矩阵-矩阵乘法：

其中

总之，我们证明了仅使用动能算子的虚时传播将一个斯莱特行列式变换为另一个斯莱特行列式；此外，这种传播是通过一个简单（且高效）的矩阵-矩阵乘法实现的。当相互作用项 $\nu$ 非零时，时间演化就不再如此简单了，因为将多体传播子 exp $( - \tau \mathcal{H} )$ 作用于一个斯莱特行列式并不会得到一个单一的斯莱特行列式。AFQMC 基于 Trotter 近似（Trotter, 1959; Suzuki, 1976a,b），这使我们能够将完整的演化 exp $( - \tau \mathcal{H} )$ 表达为一系列分别涉及动能项或势能项的不同的传播。然后，势能项通过使用所谓的 Hubbard-Stratonovich 变换（Hubbard, 1959; Stratonovich, 1957）来处理，该变换用一个依赖于经典辅助变量（这些变量需要通过随机技术积分掉）的单体传播叠加来取代多体项。

11.2 Trotter 近似#

让我们考虑方程 (11.1) 虚时演化中的一个微小步长 $\Delta \tau$；然后，我们可以定义一个接近原始哈密顿量的有效哈密顿量 $\overline{{\mathcal{H}}}$：

实际上，对于小的 $\Delta \tau$ 值，对上述方程两边进行直接展开，我们可以将 $\overline{{\mathcal{H}}}$ 评估为：

其中厄米算子 $\mathcal{O}_{\mathrm{Tr}}$ 由下式给出：

我们需要强调的是，算符 $\mathcal{O}_{\mathrm{Tr}}$ 是广延的（即它随格点数L标度变化）。事实上，给定算符 $\kappa$ 和 $\mathcal{V}$ 的实际形式（它们都是局域算符的和），对易子 $[ \mathcal{K} , \mathcal{V} ]$ 是广延的；更一般地，两个广延算符的对易子也是广延的。因此，我们在固定 $\Delta \tau$ 下用有效哈密顿量 $\overline{{\mathcal{H}}}$ 计算的所有性质，都与原始哈密顿量的性质相差一个 $\Delta \tau^{2}$ 量级的小误差。所以，与式(11.14)的 Trotter 分解相关的误差是完全可控的，因为它可以通过减小 $\Delta \tau$ 而系统地减小。此外，对于热力学量（例如每个格点的能量）的给定精度，我们在增加系统尺寸时不需要减小 Trotter 时间步长 $\Delta \tau$，即 Trotter 近似相关的误差是尺寸一致的（size consistent）。

11.3 Hubbard-Stratonovich 变换#

Hubbard-Stratonovich 变换允许我们将指数中包含密度（或自旋）中两体相互作用项的传播子，替换为更简单的单体传播子。这是以在每个格点 $i$ 上引入辅助经典场（变量）为代价的，这些变量下文记为 $\left\{\sigma_{i} \right\}$。由于 $\exp ( - \Delta \tau \mathcal{V} )$ 可以分解为 L 项，每一项作用在单个格点 i 上，我们可以独立地关注每个格点。让我们讨论通常用于排斥性哈伯德模型的退耦方法。最简单的例子是通过使用连续辅助场得到的：

其中 $\lambda^{2} = g$。事实上，在 $( n_{i , \uparrow} - n_{i , \downarrow} )$ 对角化的基底下，上式右边的项变成一个标准的高斯积分，其中 $( n_{i , \uparrow} - n_{i , \downarrow} )$ 可以仅被视为一个数。这个变换可以简化并推广到离散的，即伊辛（Ising）辅助场 $\sigma_{i} = \pm 1$，通过以下关系（Hirsch, 1985）：

其中现在双曲余弦 $\lambda$ 满足 $\cosh \lambda = \exp ( g / 2 )$。这个表达式可以很容易地通过展开 $\exp ( \lambda n_{i , \sigma} )$ 并注意到所有非零的 k 次幂贡献项为 $\lambda^{k} / k !$ 来证明（因为由于泡利原理，当 $k \geq 1$ 时 $n_{i , \sigma}^{k} = n_{i , \sigma}$）。这个事实导致：

因此，我们有：

(11.21)

这意味着式(11.18)确实在 $\cosh \lambda = \exp ( g / 2 )$ 时成立。

Hubbard-Stratonovich 变换非常有用，因为它将一个复杂的、指数中包含 $( n_{i , \uparrow} - n_{i , \downarrow} )$ 平方的两体算符，转化为依赖于辅助场的简单单体项的叠加。特别地，通过使用 $n_{i , \sigma}^{2} = n_{i , \sigma}$，多体传播子可以改写为适合 Hubbard-Stratonovich 变换的形式：

通过这种方式，对每个格点引入一个独立的伊辛场 $\sigma_{i} = \pm 1$，我们可以写出：

其中我们使用了 $g = U \Delta \tau$ 时的式(11.18)，这意味着：

在式(11.23)中，我们省略了一个无关紧要的常数 $\exp ( - U N_{e} \Delta \tau / 2 ) / 2^{L}$ ，因为我们假设总粒子数 $N_{e}$ 是固定的。因此，现在我们可以对式(11.9)中的无相互作用态应用多体算符 $\exp ( - \tau \mathcal{V} )$。事实上，通过定义：

使得：

我们得到：

其中演化后的轨道由下式给出：

这样，不仅动能项，就连势能的传播也可以通过一个 $2 L \times 2 L$ 的对角矩阵 $\mathbf{V} ( \lambda )$ 写成矩阵-矩阵相乘的形式：

其中

最后，我们想指出，还有几种其他方式进行哈伯德-斯特拉托诺维奇退耦。例如，两体算符可以写作：

那么，传播子 exp(−τ ) 包含具有负前因子的平方单体算符；因此，我们可以使用式(11.17)的变换，其中包含一个虚常数 $\lambda i \lambda$ 。这在量子蒙特卡洛计算中引入了符号问题（phase problem），可以通过使用某种近似来解决，例如所谓的约束路径量子蒙特卡洛（constrained-path quantum Monte Carlo）(Zhang et al., 1995)。值得注意的是，在半填充 $N_{e} = L$ 时，符号问题消失。据我们所知，尚未有系统性研究来确定对于给定模型最有效的哈伯德-斯特拉托诺维奇退耦方式，即在给定计算时间内最小化统计误差的方式。

11.4 路径积分表示#

现在，我们可以计算有效哈密顿量 ${\overline{{\mathcal{H}}}}$ 基态的所有性质，该哈密顿量可以任意接近精确哈密顿量，即误差随 $\Delta \tau^{2}$ 消失。为此，我们考虑伪配分函数：

其中 $2 T$ 是特罗特切片的数量，$\tau = 2 T \Delta \tau$ 。当 $T \to \infty$ 时，多体态

收敛到 ${\overline{{\mathcal{H}}}}$ 的精确基态 $| \overline{{\Upsilon}}_{0} \rangle$ 。因此，任何算符或关联函数 $\mathcal{O}$ 可以计算为：

接着，通过将式（11.34）的分子和分母中的每个特罗特切片都用式（11.23）的辅助场变换替换，可以得到路径积分表示：

其中，我们引入了单粒子传播子 $U_{\sigma} ( 2 T , 0 )$ 的紧凑定义，该传播子由辅助场变换（11.23）定义，这意味着对每个特罗特切片 $1 \leq q \leq 2 T$ 有一个相应的单体算符 $\begin{array} {r} {h_{\sigma}^{q} ( \lambda ) = \lambda \sum_{i} \sigma_{i}^{q} ( n_{i , \uparrow} - n_{i , \downarrow} )} \end{array}$；这里，场 $\sigma_{i}^{q} = \pm 1$ 除了与格点 i 关联的指标外，还获得了一个离散时间指标 $q$。于是，我们有：

由于 $U_{\sigma} ( 2 T , 0 )$ 是单粒子传播子，当它作用于一个斯莱特行列式时，会得到另一个斯莱特行列式。因此，对于伊辛场 $\{\sigma_{i}^{q} \}$ 的任何固定构型，量 $\langle \Psi_{0} | U_{\sigma} ( 2 T , 0 ) | \Psi_{0} \rangle$ 都可以进行数值计算。在式（11.34）的分子中也使用同样的场，我们得到：

然后，可以设计一个蒙特卡洛方案，通过计算两个随机变量平均值的比值来评估该方程的右端。特别地，我们考虑正权重：

这样，任何算符或关联函数的期望值由下式给出：

其中 $\langle \langle \dots \rangle \rangle$ 表示由权重（11.38）定义的分布上的统计平均，S(σ ) 是 $\langle \Psi_{0} | U_{\sigma} ( 2 T , 0 ) | \Psi_{0} \rangle$ 的符号，并且

注意，在 AFQMC 技术中，只要 $\langle \Psi_{0} | U_{\sigma} ( 2 T , 0 ) | \Psi_{0} \rangle$ 对于所有场构型 $\{\sigma_{i}^{q} \}$ 没有确定的符号，就会出现符号（或相位）问题。

一旦对于每个格点和时间切片 $q$ 给定场 $\sigma_{i}^{q}$，所有量（即随机变量和权重）都可以在多项式时间内计算。前向传播的状态 $| \overline{{{R}}}_{T} \rangle = U_{\sigma} ( T , 0 ) | \Psi_{0} \rangle$ 可以通过一个 $2 L \times N_{e}$ 矩阵写成与式（11.9）相同的形式：

其中 $\mathbf P ( \tau )$ 和 $\mathbf{V}_{q}$ 是式(11.13)和(11.30)中定义的 $2 L \times 2 L$ 矩阵，对应于每个时间切片 q，即具有场 $\{\sigma_{i}^{q} \}$。注意，矩阵 ${\bf R}_{T}$ 可以在时间切片上按 q 递增的顺序迭代计算，通过矩阵-矩阵乘法实现，总运算量为 $O ( 8 T L^{2} N_{e} )$：

初始条件为：

向后传播的态 $\langle \overline{{L}}_{T} | ~ = ~ \langle \Psi_{0} | U_{\sigma} ( 2 T , T )$ 可以用类似方式计算：

这也可以迭代计算：

初始条件为：

那么，两个波函数之间的重叠就是一个 $N_{e} \times N_{e}$ 矩阵的常规行列式：

这可以通过标准数值库以 $O ( N_{e}^{3} )$ 运算量计算。最后，式(11.40)的量子期望值 $\mathcal{O} ( \sigma )$ 也可以利用 Wick 定理在多项式时间内计算，如第11.5节所示。

11.4.1 稳定的虚时传播#

式(11.42)和(11.46)的虚时演化是在无限精度算术下定义的。然而，实际计算会受到截断误差的影响，这导致在将单粒子算符应用于斯莱特行列式并考虑大的虚时 $\tau$ 时出现不稳定问题。事实上，让我们考虑传播子的动能部分，其中演化轨道由式(11.11)给出。那么，因子 $\exp ( - \tau \epsilon_{\beta} )$ 会导致算法的指数不稳定性，因为所有轨道都将由对应于最低本征值的轨道主导；更精确地说，在双精度浮点数中 $1 + 10^{- 15} = 1$，因此，只要 $\exp [ - \tau ( \epsilon_{N_{e}} - \epsilon_{1} ) ] < 10^{- 15}$，那么在数值精度内将只有 $N_{e} - 1$ 个线性独立的轨道；这一事实导致斯莱特行列式消失。为了克服这个问题，我们必须取 $\Delta \tau$ 使得：

请注意，这个条件是尺寸一致的，因为在格点模型中，单粒子带宽在热力学极限下是有限的（例如，在最近邻跳跃 t 的 d 维情况下为 $4 d t$）。然后，在每次短时演化之后，传播轨道可以在不丢失斯莱特行列式信息的情况下进行正交化。实际上，通过使用矩阵 A 生成的任何线性变换：

斯莱特行列式除了一个整体常数 det A 外不会改变，该常数可以在整个传播过程中保存和更新。在这种背景下，格拉姆-施密特正交化（Gram-Schmidt orthogonalization）是用于此目的的最简单算法；然而，通过应用乔列斯基分解（Cholesky decomposition）可以获得更高效的方法。这里，考虑 $N_{e} \times N_{e}$ 重叠矩阵并进行分解：

其中 $\psi$ 现在是一个通用矩阵，L 是一个下三角矩阵。那么，由以下给出的新轨道：

将相互正交。实际上，我们有：

这些轨道将描述相同的行列式，但有一个整体常数因子 $\det \mathbf{L}$，由于该因子等于 $\mathbf{L}$ 的对角元素乘积，因此很容易计算。在本文描述的基态技术中，斯莱特行列式的归一化是一个无关紧要的常数，它不会出现在任何物理量中，因此可以忽略。正交化后，可以进行另一次稳定的传播，而不会遇到任何不稳定性问题。

11.5 顺序更新#

在本节中，我们描述一个简单的马尔可夫过程，可用于对公式 (11.38) 中的权重进行采样。特别地，对于每个时间切片 $q$，只改变位点 $k$ 上的单个伊辛变量。为了加速更新，重要的是在移动到下一个时间切片之前，先提出给定时间切片内所有位点的辅助场新值。如果在位点 $k$ 翻转场 $\sigma_{k}^{q} - \sigma_{k}^{q}$（意味着在位点上进行顺序循环），只有右波函数 $| R_{q} \rangle$ 会以简单的方式改变：

其中 $\lambda_{\uparrow} = \exp ( - 2 \lambda \sigma_{k}^{q} ) - 1$，$\lambda_{\downarrow} = \exp ( 2 \lambda \sigma_{k}^{q} ) - 1$。对于所提出的伊辛变量翻转，权重 $\mathcal{W} ( \sigma ) = | \langle L_{q} | R_{q} \rangle |$ 的变化以及相应的符号（或相位）$S ( \sigma )$ 可以通过以下比率更新：

该比率表示在两个不同的斯莱特行列式 $\langle L_{q} |$ 和 $| R_{q} \rangle$ 之间计算的多体关联函数。最后，接受概率由第 3.9 节描述的标准 Metropolis 算法获得；一种略有改进的方案，具有更小的关联时间和更快的平衡速度，由所谓的热浴算法（heat-bath algorithm）（Krauth, 2006）给出，其中提出的构型 $\sigma_{k}^{\prime}$ 以概率被接受：

比率 $\mathcal{R}_{k}^{q}$ 可以根据等时格林函数（equal-time Green's function）来计算，该函数由一个 $2 L \times 2 L$ 矩阵定义：

在最一般的情况下，即使对于非对角自旋分量，即 $I ~ \leq ~ L$ 且 $J \ > \ L$ 或 $I \ > \ L$ 且 $J \ \leq \ L$ 的情况，$( G_{q} )_{J , I}$ 也可能非零（一旦对自旋向下电子执行粒子-空穴变换，这些元素与反常平均值（anomalous averages）相关）。相反，当斯莱特行列式被分解为自旋向上和自旋向下分量时，$( G_{q} )_{J , I}$ 具有块对角形式。定义了上述格林函数后，所有可能的两个不同行列式之间的多体关联函数都可以使用维克定理（Wick theorem）（针对左矢（bra）和右矢（ket）状态不同的情况）进行计算。

现在让我们计算公式 (11.58) 中的格林函数。首先，两个斯莱特行列式之间的重叠由相应重叠矩阵的行列式给出，参见公式 (11.49)；例如，公式 (11.58) 右边分母由下式给出：

分子可以采用类似方式进行计算，注意到它可以写成两个具有$N_{e} + 1$个电子的斯莱特行列式之间的交叠：

其中

对应于交叠$\langle \Psi_{J} | \Psi_{I} \rangle$的新矩阵M 可以通过认识到额外$N_{e} + 1$个轨道以特别简单的方式进入来计算：

其中我们突出了额外行/列的贡献，而所有其他矩阵元保持与原始$N_{e} \times N_{e}$矩阵M完全相等。为了计算上述矩阵的行列式，我们利用一个性质：如果我们将最左边一列加上其他列的适当线性组合，行列式保持不变。这可以通过消除其所有值（除了对角元）的方式实现：线性组合的系数$x_{\beta}$由下式给出：

从而得到：

因此，我们得到：

将此结果与式(11.60)结合，我们得到：

这可以表示为紧凑形式：

该量可以在计算开始时使用标准线性代数运算从头计算，然后在每次接受蒙特卡洛移动时进行更新（见下文）。出于稳定性考虑，需要在选定时间切片（例如，对于$q = m \times p$，其中$m$是整数，$p$是固定的小整数）从头重新计算格林函数。

通过类似的推导，我们可以得到包含四个费米子算符的表达式：

这使得我们可以计算式(11.56)中的比值：

这里最后一项仅当行列式没有分解为自旋向上和自旋向下块时才存在。

一旦试探移动$| R_{q} \rangle \to | R_{q}^{\prime} \rangle$被接受，必须更新格林函数：

下面将证明，这无需从头重新计算，而只需简单应用维克定理即可完成。注意，如果斯莱特行列式对自旋向上和自旋向下电子是因子化的，那么方程(11.70)中的最后一项消失，格林函数的更新也会简化。若不假设斯莱特行列式因子化，我们得到的最终表达式为：

其中

其中 $\overline{{( G_{q} )}}_{J , I} = ( G_{q} )_{J , I} - \delta_{J , I}$ 。我们想指出，每当格林函数因子化为自旋向上和自旋向下分量时，可以采用一个更简单的更新方式，对各自旋部分独立进行：

其中现在的系数为：

总之，对于给定时间切片 $q$ 的更新方案可以在 $O ( L^{3} )$ 次运算内完成。为了从一个时间切片过渡到相邻切片，我们需要使用方程(11.47)的矩阵-矩阵乘法传播左边的斯莱特行列式 $\langle L_{q} |$ ，这同样最多需要 $O ( L^{3} )$ 次运算。而右边的斯莱特行列式 $| R_{q} \rangle$ 可以一开始就存储起来并保存在内存中（否则，也可以使用方程(11.43)反向传播，即 ${\bf R}_{q - 1} =$ $\mathbf{P} ( - \Delta \tau ) [ \mathbf{V}_{q} ( \lambda ) ]^{- 1} \mathbf{R}_{q} $ ）。格林函数同样可以通过类似的矩阵-矩阵运算，从时间切片 $q$ 反向传播到时间切片 $q - 1$ 。实际上，由方程(11.68)以及左右矩阵的传播，我们有：

因此，遍历所有时间切片一次最多需要 $O ( T L^{3} )$ 次运算。所以，该算法对于基态性质的计算特别高效，因为 $T \to \infty$ 的极限可以在仅与时间切片数成线性关系的代价下达到。

11.5.1 延迟更新#

更新过程方程(11.72)中的基本操作是对一个通用的 $2 L \times 2 L$ 矩阵进行秩-1更新：

其中 $a_{J}^{\sigma}$ 定义于方程(11.73)和(11.74)，而 $b_{I}^{\uparrow} ~ = ~ ( G_{q} )_{k , I}$ 和 $b_{I}^{\downarrow} ~ =$ $( G_{q} )_{k + L , I} $ 。正如第5.6节已经强调过的，当系统规模较大时，如果矩阵 $\mathbf{G}_{q}$ 不能完全装入处理器缓存，这类操作可能会在计算上效率低下。克服这一缺点的一种方法是延迟矩阵的更新，同时不丢失其信息。这可以通过存储一组左、右向量 $\mathbf{a}^{( p )}$ 和 $\mathbf{b}^{( p )}$（其中 $p = 1 , \ldots , 2 m$ ，包括“自旋”指标 $\sigma $ ）以及“初始”矩阵（记为 ${\bf G}_{{q}}^{0} .$ ）来实现，延迟更新即从该初始矩阵开始。那么，经过 m 次更新后的矩阵 $\mathbf{G}_{q}$ 由下式给出：

每次接受新构型时，新的向量 $\mathbf{a}^{( m + 1 )} , \ \mathbf{a}^{( m + 2 )} , \ \mathbf{b}^{( m + 1 )}$ 和 $\mathbf{b}^{( m + 2 )}$ 可通过少量操作，基于 ${\bf G}_{q}^{0}$ 和之前的向量（$p = 1 , \ldots , 2 m$）计算得出。注意，一旦矩阵 ${\bf G}_{q}$ 被写成式(11.81)的形式，计算求和项中各因子所需的操作量为 O(8mL)，当 $m \ll L$ 时，这相对于完全更新来说是微不足道的。

通过执行这种延迟更新（delayed update），我们可以找到最优的 $m_{\mathrm{max}}$，并通过一次标准的矩阵乘法来计算完整的矩阵 $\mathbf{G}_{q}$：

其中 A 和 B 是 $2 L \times 2 m_{\mathrm{max}}$ 的矩阵，分别由 $p = 1 , \ldots , 2 m_{\mathrm{max}}$ 的向量 $\mathbf{a}^{( p )}$ 和 $\mathbf{b}^{( p )}$ 构成。之后，我们可以用新的 ${\bf G}_{q}^{0} = {\bf G}_{q}$ 继续下一次延迟更新，并将整数 m 重置为零。这种更新过程的优势在于，经过一个 $2 m_{\mathrm{max}}$ 步的周期后，计算的主体部分由式(11.82)中的矩阵-矩阵乘法完成，这比 $2 m_{\mathrm{max}}$ 次秩-1更新（rank-1 update）高效得多（并且不受缓存内存大小的限制）。对于大量电子，延迟更新过程可以使蒙特卡洛代码的速度提高约一个数量级。

11.6 基态能量与关联函数#

基态能量或任何其他可观测量，包括关联函数，都可以通过使用第11.5节讨论的维克定理（Wick theorem）轻松计算。事实上，考虑一个一般的 n 体算符：

其中 $d_{I}^{\dagger} \left( d_{I} \right)$ 在位点 $I = 1 , \ldots , 2 L$ 上产生（湮灭）一个电子。那么，式(11.34)的关联函数通过计算式(11.40)给出的 $\mathcal{O} ( \sigma )$ 来评估，这可以利用如下事实：$| \overline{{R}}_{T} \rangle = U_{\sigma} ( T , 0 ) | \Psi_{0} \rangle$ 和 $\langle \overline{{L}}_{T} | = \langle \Psi_{0} | U_{\sigma} ( 2 T , T )$ 都是斯莱特行列式（Slater determinants），因此可以应用维克定理。

最后值得一提的是，基态能量（在 $T \to \infty$ 时得到）的误差以 ${\cal O} ( \Delta \tau^{4} )$ 的量级消失。实际上，鉴于式(11.14)的对称特罗特近似（symmetric Trotter approximation），有效哈密顿量 $\overline{{\mathcal{H}}}$ 与原始哈密顿量 $\mathcal{H}$ 的差别在于一个消失量为 $O ( \Delta \tau^{2} )$ 的微扰项，参见式(11.15)。$\overline{{\mathcal{H}}}$ 的基态可以通过式(11.33)的投影技术得到，其与 $\mathcal{H}$ 的精确基态之间可以利用标准微扰论建立联系，即 $| \overline{{\Upsilon}}_{0} \rangle = | \Upsilon_{0} \rangle + \Delta \tau^{2} | \Upsilon^{\prime} \rangle + O ( \Delta \tau^{4} )$（其中 $| \Upsilon^{\prime} \rangle$ 与 $| \Upsilon_{0} \rangle$ 正交）。因此，当计算 $\mathcal{H}$ 在 $| \overline{{\Upsilon}}_{0} \rangle$ 上的期望值时，主导的 $O ( \Delta \tau^{2} )$ 修正项消失，我们得到主导误差为 ${\cal O} ( \Delta \tau^{4} )$。

11.7 无符号问题的简单情形#

让我们讨论几个AFQMC技术不受符号问题（sign problem）影响的例子。首先，我们考虑吸引型哈伯德模型（attractive Hubbard model），即 $U < 0$ 的情况（以及任意实的动能项）。在这种情况下，可以使用一个与总密度 $\sigma_{i} ( n_{i , \uparrow} + n_{i , \downarrow} - 1 )$ 耦合的实场 $\sigma_{i}$ 进行辅助场变换。事实上，式(11.18)的离散变换可以修改为：

其中双曲余弦 $\lambda ~ = ~ \exp ( g / 2 )$，$g = | U | \Delta \tau$。然后，我们采用一个试探波函数 $\left| {\Psi_{0}} \right.$，该函数对两个自旋分量使用相同的实轨道，即 $| \Psi_{0} \rangle =$ $| \Psi_{0 , \uparrow} \rangle \otimes | \Psi_{0 , \downarrow} \rangle$，其中：

在这种情况下，单粒子传播子 $U_{\sigma} ( 2 T , 0 )$ 可以分解为上自旋和下自旋部分，即 $U_{\sigma} ( 2 T , 0 ) = U_{\sigma}^{\uparrow} ( 2 T , 0 ) \times U_{\sigma}^{\downarrow} ( 2 T , 0 )$，因为总传播对波函数的上自旋和下自旋分量的作用方式相同。其结果是，式(11.35)的赝配分函数 $\mathcal{Z}$ 变为：

因此，蒙特卡洛权重严格为正，因为它是实数的平方。总之，对于吸引型（自旋无关）相互作用，在任何晶格结构下，只要上自旋粒子数与下自旋粒子数严格相等，AFQMC 就不存在符号问题。

另一个没有符号问题的例子是半满填充下（即 $N_{e}^{\uparrow} + N_{e}^{\downarrow} = L$）二分格点上的排斥型哈伯德模型（例如，动能项仅耦合相反子格点 A 和 B）。实际上，我们总可以通过对下自旋粒子进行粒子-空穴变换（particle-hole transformation），将排斥型模型映射为吸引型模型：

其中，若格点 i 属于 A (B) 子格点，则 $s_{i} = 1 \ ( s_{i} = -1 )$；引入这一符号是为了不改变动能项。然后，下自旋电子的数量在粒子-空穴变换后发生改变（见第5.6节）：

这意味着"新"自旋向上电子的总数为 $N_{e}^{\uparrow}$，而"新"自旋向下电子的总数为 $L - N_{e}^{\downarrow}$。因此，基于吸引模型的结果，我们得出结论：当粒子-空穴变换后自旋向下电子数恰好等于自旋向上电子数，即 $N_{e}^{\downarrow} + N_{e}^{\uparrow} = L$ 时，不存在符号问题，而这正是半填充条件。

对于 $U / t = 4$ 的排斥哈伯德模型，图11.1给出了 $N_{e} = 42$ 个电子在50个格点上的结果。初始波函数 $\left| {\Psi_{0}} \right.$ 由无相互作用波函数给出，对应闭壳层构型（即最高占据能级与最低未占据能级之间存在有限尺寸能隙）。我们给出了基态能量、双占据数 $D = 1 / L \sum_{i} n_{i , \uparrow} n_{i , \downarrow}$ 以及模拟过程中的平均符号 $\langle \langle S ( \sigma ) \rangle \rangle$ 的结果。需要指出，在这种情况下，随着投影时间 τ（与特罗特分解相关联的离散化时间为 $\Delta \tau = 0.1 $）的变化，能快速收敛到精确值。我们想强调的是，由于闭壳层条件，平均符号在 $\tau \approx 2$ 之前保持接近1，然后非常缓慢地降至0，从而允许对物理量进行精确评估。

图11.1 哈伯德模型在50个格点、$N_{e} = 42$ 个电子、$U / t = 4$ 下的结果；特罗特近似中的时间离散化为 $\Delta \tau = 0.1$。所有量均显示为投影时间 τ 的函数：基态能量（上栏）、双占据数 $D = 1 / L \sum_{i} n_{i , \uparrow} n_{i , \downarrow}$（中栏）和平均符号（下栏）。

最后，我们想指出，当存在长程相互作用时（在二分格点上），无符号分解也是可能的。事实上，令：

那么，我们可以考虑多体传播子的如下哈伯德-斯特拉托诺维奇变换：

其中唯一的约束是 U 必须为正定矩阵。因此，每当动能项分别作用于自旋向上和自旋向下电子时，传播子便得以因子化。此外，通过执行方程(11.87)和(11.88)的粒子-空穴变换，我们立即得到自旋向上和自旋向下传播子互为复共轭，从而产生正权重：

11.8 实际实现#

在此，我们概述辅助场量子蒙特卡洛算法在实际实现中的关键步骤。

1. 计算开始时的初始化。

   _• 开始之前，我们必须确定特罗特切片的总数 $2 T$ 以及两次连续正交化（使用乔列斯基算法）之间的时间切片数 $p$。特别地，$p$ 应选择得足够小，以保证计算有足够的精度，这可以通过减小 $p$ 来系统地进行测试（$p = 1$ 是最精确的情况）。

   从每个位点和时间切片随机生成的伊辛场 $\sigma_{i}^{q} = $ ±1 开始。$2 L \times 2 L$ 矩阵 $\mathbf{P} ( \pm \Delta \tau )$ 和 $\mathbf{P} ( \pm \Delta \tau / 2 )$ 必须在开始时通过方程 (11.13) 计算并存储在内存中。而对角矩阵 ${\bf V}_{q} ( \lambda )$ 则无需存储，因为利用这些矩阵是对角矩阵的性质（一个对角矩阵乘以一个完整矩阵仅需 $4 L^{2}$ 次操作），所有计算都可以高效完成。

2. 初始传播。

   _• 使用方程 (11.43) 将 $\mathbf{R}_{q}$ 向前传播 $q = 1 , \ldots , 2 T$ 步，并出于稳定性考虑，每 $p$ 步正交化一次传播；在计算机内存中保存对应于 $q = k \times p$ 时波函数 $| R_{q} \rangle$ 的 $2 L \times N_{e}$ 矩阵（其中 $k$ 为整数）。到达最后一个区间 $q = 2 T$ 时，初始化方程 (11.48) 中反向传播的波函数 $\mathbf{L}_{2 T}$，并从头计算方程 (11.68) 中的格林函数 $\mathbf{G}_{2 T}$。

3. 使用梅特罗波利斯算法的马尔可夫过程。

   _• 在时间切片 q，开始对位点 $i = 1 , \ldots 2 L$ 进行遍历。每一步，我们可以使用方程 (11.70) 计算新权重与旧权重的比率。然后，可以使用梅特罗波利斯算法或热浴算法来生成马尔可夫过程。每当提议的构型被接受时，根据方程 (11.72) 更新格林函数。对于大规模计算，也可以使用延迟更新以获得更好的性能。

   _• 在对格点位点遍历结束时，我们使用方程 (11.47) 将反向传播状态从 $q$ 更改为 $q - 1$。注意，在传播 p 步之后，出于稳定性考虑，必须使用乔列斯基正交化。在此阶段，还需要根据方程 (11.79) 传播格林函数。只有在 $\mathbf{L}_{q}$ 被正交化时，我们才从头计算格林函数（使用初始传播中保存的 $\mathbf{R}_{q}$）。

   _• 回到初始传播并继续蒙特卡洛过程以累积统计量。注意，算法的一个小改进是反转扫描方向，即从 $q = 1$ 扫描到 $q = 2 T$，从而可以避免初始传播。这种方法通常可以节省大约10%的总计算量。

4. 可观测量的计算。

   可观测量（即总能量）可以在沿时间切片的扫描过程中计算。特别地，对于 $q \approx T$（即中间区间附近），我们得到方程 (11.33) 中状态 $| \Psi_{T} \rangle$ 上的期望值。在这种情况下，我们应该提醒，用于计算关联函数的格林函数 $\mathbf{G}_{q}$ 需要用动能传播半个时间步：

$$ \mathbf{G}{q} \to \mathbf{P} \left( \frac{\Delta \tau} {2} \right) \mathbf{G}{q} \mathbf{P} \left( - \frac{\Delta \tau} {2} \right) ,\tag{11.94}

$$

这对于以足够小的特罗特误差计算物理量尤其重要。

第五部分 高级专题