第4章 · 朗之万分子动力学#

4.1 引言#

本章介绍一种模拟有限体积 $\Omega$ 中给定数目 $N_{p}$ 个经典粒子在温度 T 下的最简单但稳健的方法，即用于计算所谓 NVT 系综（Allen and Tildesley, 1987; Tuckerman, 2010）中的物理量。该方法通过一阶朗之万动力学（first-order Langevin dynamics）实现，该动力学由一组随机微分方程的解定义：

其中 $\mathbf{R} ( t )$ 是一个（时间依赖的）D 维向量，其分量为 $R_{\alpha} \ ( \alpha = 1 , \ldots , D$ ，且 $D = N_{p} \times d$ ，d 为空间维度）。换句话说，这个 D 维向量是所有粒子位置的简洁表示，这些粒子通过经典势 $V ( \mathbf{R} )$ 相互作用，该势定义了确定性“力”的分量 $f_{\alpha} ( \mathbf{R} )$

最后，$\pmb{\eta} ( t )$ 是一个随机向量，其分量为 $\eta_{\alpha} ( t ) ~ ( \alpha = 1 , \ldots , D )$ ，表示具有零均值的随机噪声，且不同分量 $\alpha \neq \beta$ 以及不同时刻 $t \neq t^{\prime}$ 之间无相关性，即所谓的白噪声（white noise）：

其中 $\langle \dots \rangle$ 表示期望值。我们需要强调，一阶朗之万方程（4.1）与牛顿运动方程有本质区别，后者是将实际力与粒子加速度联系起来的二阶（确定性）微分方程。由牛顿方程生成的动力学可用于模拟固定能量 $E_{\ast}$ 下有限体积内的一团经典粒子，即所谓的 NVE 系综（Allen and Tildesley, 1987; Tuckerman, 2010）。一方面，一阶朗之万方程可视为一种特设方法，用于生成根据玻尔兹曼分布（Boltzmann distribution）达到平衡的构型（参见下文）。另一方面，它们可以从布朗运动（Brownian motion）的粗粒化描述中导出：在此描述中，流体中质量为 m 的小颗粒除了受到外力 $\mathbf{f} ( \mathbf{R} )$ 外，还受到摩擦系数为 $\gamma$ 的摩擦力以及由流体中随机密度涨落引起的随机力。这种方法导出了二阶朗之万动力学：

然后，在 $m \to 0$ 的极限下得到一阶朗之万方程（在此极限下，γ 仅设定时间尺度）。通常，此极限被称为二阶朗之万动力学的过阻尼 regime（over-damped regime）。

与常微分方程不同，从给定初始条件出发，随机微分方程并不提供唯一解，而是使可能的解具有统计意义：从 $t = t_{0}$ 时刻的 ${\bf R}_{0} \equiv{\bf R} ( t_{0} )$ 出发，根据噪声的具体实现，可能产生多条随机轨迹。因此，一个有意义的量是在时间 t 找到给定构型 R 的概率 $\mathcal{P} ( {\bf R} , t )$ ，初始条件为：

下面我们将展示，经过一段平衡时间后，$\mathcal{P} ( {\bf R} , t )$ 收敛到一个与初始条件 $\mathbf{R}_{0}$ 无关的平衡分布。特别地，平衡概率由玻尔兹曼分布给出：

其中 $\mathcal{Z}$ 是配分函数（partition function），即归一化常数：

在4.4节中，我们将展示如何从朗之万方程（4.1）出发，证明 $\mathcal{P} ( {\bf R} , t )$ 满足一个偏微分方程——福克-普朗克方程（Fokker-Planck equation），其稳态解确实是玻尔兹曼分布（4.8）。然而，在此之前，我们首先讨论如何求解给定的随机微分方程（4.1）。由于这些方程的非线性特性——势 $V ( \mathbf{R} )$ 通常包含粒子之间的相互作用势——解析解是不可能的。我们需要使用数值方法求解随机微分方程。

其中 $\mathcal{Z}$ 是配分函数，用于概率的归一化条件：

因此，微分方程（4.1）的解可用于在离散时间点 $t_{n}$ 上采样大量（N个）构型 $\mathbf{R}_{n}$，以表示正则分布 $\mathcal{P}_{\mathrm{eq}} ( \mathbf{R} )$，这使我们能够计算任意关联函数 ${\mathcal{O}} ( \mathbf{R} )$：

这与3.1节讨论的普通蒙特卡洛方法类似。

对于给定系统，没有通用规则说明使用分子动力学还是蒙特卡洛方法更方便。然而，在某些情况下，分子动力学是最佳选择。例如，在一些第一性原理（ab initio）方法中，经典势 $V(\mathbf{R})$ 并非简单形式，而是复杂算法（例如，表示玻恩-奥本海默近似下密度泛函理论得出的离子-离子相互作用）的结果。此时，计算所有力分量与执行单一能量计算（用于 Metropolis 接受步骤）所需的计算量相当。在这种情况下，分子动力学通常比蒙特卡洛更便捷，因为用相同的计算量就能一次性改变所有原子的位置。此外，无需采用限于单个或少数原子的局域移动来维持良好的接受率。相比之下，蒙特卡洛方法的优势在于除统计偏差外没有其他偏差，而分子动力学则始终需要引入时间离散化，这会产生系统性但可控的误差（我们将在后续讨论）。

4.2 离散时间朗之万动力学（Discrete-Time Langevin Dynamics）#

在计算机模拟中，需要时间离散化才能获得微分方程的解。因此，为了定义模拟有限温度下经典粒子的合适算法，我们在有限区间 $( t_{n} , t_{n + 1} )$ 上对方程（4.1）两边积分，其中 $t_{n} = t_{0} + \Delta n$ 是离散时间点（$t_{0}$ 是初始时间，$\Delta$ 表示离散时间步长）。这样，我们得到（在 $\Delta$ 的最低阶近似下）：

其中 ${\bf R}_{n} = {\bf R} ( t_{n} )$，${\bf f}_{n} = {\bf f} [ {\bf R} ( t_{n} ) ]$；这里，我们使用最低阶近似 $\Delta \mathbf{f}_{n}$ 来近似该区间内力的积分，因为在微小时间区间内力近似为常数，参见方程（4.15）和（4.16）。注意，这是上述积分中唯一的近似。方程（4.11）右边的积分可以通过注意到多个随机变量之和（因此也是积分）是高斯随机变量（参见2.7节）来估计。特别地，通过引入一个正态（高斯）分布的随机向量 ${\bf z}_{n}$（其均值为零，方差为单位矩阵），我们有：

其中右边的系数给出了积分的正确方差：

这里我们使用了方程（4.4）。此外，对于 $n \neq m$，由于被积函数始终为零，因此 $\langle z_{\alpha , n} z_{\beta , m} \rangle = 0$。

综合上述结果，我们可以写出离散时间朗之万方程的最终表达式：

这个迭代方案定义了新的坐标 $ \mathbf{R}_{n+1} $（在时刻 $ t_{n+1} $）如何由旧坐标 $ \mathbf{R}_n $（在时刻 $ t_n $）、离散化的力 $ \mathbf{f}_n $ 以及一组高斯变量 $ \mathbf{z}_n $ 给出。这个迭代过程恰好是一个马尔可夫过程，可以通过简单的迭代算法实现，因为只需要有一个算法能够计算给定 $ N $ 个经典粒子任意位置 $ \mathbf{R}_n $ 时的力即可。需要强调的是，当 $ \Delta 0 $ 时，噪声项（与 $ \sqrt{\Delta} $ 成正比）主导了确定性力项（与 $ \Delta $ 成正比）。注意，噪声项 $ \pmb{\eta}(t) $ 的存在使得朗之万方程的解变得不连续且不可微，原因如下：

这证伪了方程 (4.11) 中确定性力积分的近似：

由于方程 (4.15) 的结果，当 $ \Delta 0 $ 时，实际轨迹实际上是不确定的。尽管如此，我们将会证明方程 (4.14) 的时间离散化允许我们确定概率分布 $ \mathcal{P} ( \mathbf{R} , t ) $ 在 $ \Delta 0 $ 极限下的演化，且不存在不确定性。最后我们指出，在零温 $ T = 0 $ 的极限下，噪声项从运动方程中消失，因此算法变为确定性的。特别地，它退化为最速下降法（steepest descent method, Press 等人, 2007），在 $ n \infty $ 且 $ \Delta $ 足够小时，以确定性的方式给出势能 $ V ( \mathbf{R} ) $ 的最小值（或更一般地，局部最小值）。

4.3 从朗之万方程到福克-普朗克方程#

在这里，我们推导出方程 (4.14) 中离散化朗之万动力学的主方程，该方程定义了经典变量 $ \mathbf{R} $ 在有限时间步长 $ \Delta $ 下概率 $ \mathcal{P}_n ( \mathbf{R} ) $ 的演化。最重要的是，我们可以取 $ \Delta 0 $ 的极限，从而得到所谓的福克-普朗克方程（Fokker-Planck equation），它是朗之万方程 (4.1) 的主方程。

由于离散化的朗之万动力学定义了一个马尔可夫过程，概率 $ \mathcal{P}_n ( \mathbf{R} ) $ 完全由与该马尔可夫步相关的条件概率 $ K ( \mathbf{R}' | \mathbf{R} ) $ 决定：

条件概率可以通过注意到在方程 (4.14) 中，只有变量 $ \mathbf{z}_n $ 是随机的，而力是完全确定性的这一事实来确定。因此，给定 $ \mathbf{R} $，新变量 $ \mathbf{R}' $ 由确定性部分 $ \mathbf{R} + \Delta \mathbf{f}( \mathbf{R} ) $ 加上一个均值为零、方差为 $ 2 T \Delta $ 的正态分布随机噪声组成。于是有：

显然这个概率是归一化的，满足 $ \begin{array} {r} {\int \mathbf{dR}^{\prime} K ( \mathbf{R}^{\prime} | \mathbf{R} ) = 1} \end{array} $。将这种形式的条件概率代入主方程，我们得到：

现在，我们进行对 $ R_{\alpha} $ 的积分。因此，对于固定的 $ \mathbf{R}' $ 和 $ \mathbf{z} $，我们需要求解方程 (4.19) 中 δ 函数自变量的零点：

它表示一组关于 D 个未知数（$R_{\alpha}$）的 D 个非线性方程。这个问题显然难以求解；不过，我们可以通过对小量 $\Delta$ 进行系统展开来找到解：

这样，我们就能在主方程 (4.19) 中对 $R_{\alpha}$ 完成积分，并得到：

其中 $f_{\alpha}^{\prime} ( {\bf R}^{\prime} ) = \partial f_{\alpha} / \partial R_{\alpha}^{\prime}$，且在 $\mathbf{R}^{\prime} {\mathrm{:}}$ 处计算；此关系式精确到 ${\cal O} ( \Delta^{3 / 2} )$ 阶。进一步对 $\Delta$ 进行到首阶展开，我们得到：

最后，我们可以展开概率分布：

原则上，泰勒展开对于大的 $z_{\alpha}$ 值并不成立；然而，这不成问题，因为所有关于 $\{z_{\alpha} \}$ 的积分都主要贡献在 $| z_{\alpha} | ~ < ~ 1$ 的区域。将上述展开代入式 (4.23) 并完成对 $\{z_{\alpha} \}$ 的高斯积分，我们得到（将 $\mathbf{R}^{\prime}$ 改为 $\mathbf{R}$）：

现在，可以取极限 $\Delta \to 0$。确实，对于小的 $\Delta$：

这使我们导出了概率密度 $\mathcal{P} ( {\bf R} , t )$ 的福克-普朗克方程（Fokker-Planck equation）：

不难证明，式 (4.8) 中的玻尔兹曼分布（Boltzmann distribution）是福克-普朗克方程的一个定态解；事实上，我们可以将式 (4.27) 重写为：

其右端在玻尔兹曼分布下为零。福克-普朗克方程的另一个重要性质是右端为全散度。这正是概率归一化 $\int dR \mathcal{P} ( {\bf R} , t ) = N ( t ) = 1$ 的推论，该归一化代表方程的一个运动常数。事实上，对方程两边在给定体积上积分，并对右端应用高斯定理，可得：

其中

而 ${\bf n} ( {\bf R} )$ 是垂直于曲面 S 的单位向量。于是，${\bf A} ( {\bf R} , t )$ 在无穷远处为零（因为概率在无穷远处为零），这意味着 $N ( t )$ 与 t 无关。

4.4 福克-普朗克方程与量子力学#

在此，我们展示给定初始条件 $\mathcal{P} ( \mathbf{R} , t_{0} )$ 后，福克-普朗克方程如何计算任意时刻的概率 $\mathcal{P} ( {\bf R} , t )$。此外，我们还将证明，从任意初始概率出发，$\mathcal{P} ( {\bf R} , t )$ 在长时间 t 下会以指数方式趋近于平衡分布。

福克-普朗克方程与虚时薛定谔方程之间存在深刻联系 (Nelson, 1966; Parisi and Wu, 1981)。这一联系可通过将方程 (4.27) 的解写成如下形式得到：

其中 $\Upsilon_{0} ( \mathbf{R} ) = \sqrt{\mathscr{P}_{\mathrm{eq}} ( \mathbf{R} )}$ 表示一个归一化的量子态：

将上述 $\mathcal{P} ( {\bf R} , t )$ 的定义代入福克-普朗克方程 (4.27)，可得 $\Phi ( \mathbf{R} , t )$ 满足虚时薛定谔方程：

其中 $\mathcal{H}_{\mathrm{eff}}$ 是一个有效哈密顿量，由下式给出：

这里，逆（经典）温度扮演了粒子质量的角色。此外，$V_{\mathrm{eff}} ( \mathbf{R} )$ 是一个依赖于经典势 $V ( \mathbf{R} )$ 的有效势：

注意，在零温极限下，原始势 $V ( \mathbf{R} )$ 的极小值也是 $V_{\mathrm{eff}} ( \mathbf{R} )$ 的极小值；事实上，上式的第一项（在 $T \to 0$ 时占主导）在原始势 $V ( \mathbf{R} )$ 的极小值处为零。有效势可以重写为：

现在，$\Upsilon_{0} ( {\bf R} )$ 是 $\mathcal{H}_{\mathrm{eff}}$ 的一个本征态，其能量 $E_{0} = 0$，这可以通过直接检查方程 (4.34) 和 (4.36) 得到验证。由于它无节点（例如，$\mathcal{P}_{\mathrm{eq}} ( \mathbf{R} )$ 对所有构型 R 均为正），它也是有效哈密顿量的实际基态。于是，薛定谔方程以及相应的福克-普朗克方程的解，可以通过将初始条件按 $\mathcal{H}_{\mathrm{eff}}$ 的本征态 $\Upsilon_{n} ( \mathbf{R} )$ 展开，形式上给出闭式解：

其中

这意味着根据 $\mathcal{P} ( \mathbf{R} , t_{0} )$ 的归一化条件，$a_{0} = 1$。因此我们得到概率 $\mathcal{P} ( {\bf R} , t )$ 的完整演化：

因此，对于大时间 t，$\mathcal{P} ( {\bf R} , t )$ 指数收敛到平稳平衡分布 $\Upsilon_{0}^{2} ( {\bf R} ) = \mathcal{P}_{\mathrm{eq}} ( {\bf R} )$。平衡的特征时间 τ 由第一激发态的能隙倒数给出，即 $\tau = 1 / E_{1}$。

让我们通过证明由福克-普朗克方程生成的演化满足细致平衡条件来结束本节。薛定谔方程的形式解由下式给出：

因此，利用方程 (4.31)，概率密度的演化可写为：

这意味着条件概率由下式给出：

由于 $\mathcal{H}_{\mathrm{eff}}$ 是对称的，即 $\langle \mathbf{R}^{\prime} | e^{- \mathcal{H}_{\mathrm{eff}} t} | \mathbf{R} \rangle = \langle \mathbf{R} | e^{- \mathcal{H}_{\mathrm{eff}} t} | \mathbf{R}^{\prime} \rangle$，我们有：

这表明细致平衡条件确实被满足。我们要强调的是，时间离散化引入的误差会破坏细致平衡条件。在这个意义上，对于任何有限值的时间步长 $\Delta$，平衡分布并非由玻尔兹曼分布 (4.8) 给出，但会在 $\Delta 0$ 时退化为该分布。

4.4.1 谐振子情形的精确解#

为了理解时间离散化如何影响动力学，我们考虑一个一般的二次型势：

其中 $\mathbf{R}_{\mathrm{eq}}$ 是平衡位置。力常数矩阵 K 是对称的，因此可以通过酉矩阵 U 对角化：

其特征值记为 $\bar{K}_{\alpha}$，且 $\bar{K}_{1} \le \bar{K}_{2} \le \dots \le \bar{K}_{D}$，即对角矩阵 $\bar{\mathbf{K}}$ 的对角元素。方程 (4.35) 中考虑的有效势在这种情况下可以显式给出，有效哈密顿量 $\mathcal{H}_{\mathrm{eff}}$ 仍为谐振子形式，且有：

其中有效谐振子耦合为：

由于有效哈密顿量是二次型的，我们得到一个谐振子问题，其力常数矩阵如上所述，且所有自由度都具有质量 $m = 1 / 2 T$。因此，根据标准计算，动力学矩阵 K 的特征值提供了谐振子问题的本征模式（忽略一个能量偏移）：

其中 $n_{\alpha}$ 是非负整数（基态对应于所有 $\alpha$ 值下 $n_{\alpha} = 0$）。这里，我们利用了式 (4.46) 中的常数项是矩阵 K 的迹，即其所有特征值之和。该情况下的关联时间可以显式给出为：

其中 $\bar{K}_{1}$ 是谐振子势的最小特征值。值得注意的是，关联时间 $\tau$ 不依赖于温度 $T$。

在谐振子情形下，朗之万方程 (4.14) 的离散化版本可以显式求解，如下所述。实际上，将式 (4.14) 两边乘以 U 并引入简正模式坐标 ${\bf Q} = {\bf U} ( {\bf R} - {\bf R}_{\mathrm{eq}} )$，离散化的朗之万动力学变为：

其中 ${\bf y}_{n} = {\bf U} {\bf z}_{n}$ 是满足以下条件的正态分布随机变量：

由于各个简正模式之间不耦合，我们可以单独处理每个模式。对应于式 (4.50) 的主方程依赖于条件概率 $K ( Q_{\alpha}^{\prime} | Q_{\alpha} )$，该情况下的条件概率可以显式推导得出：

然后，我们得到概率分布：

其中

满足细致平衡条件：

事实上，利用式 (4.52)，我们得到：

这与分布 (4.53) 是一致的，其中 $K_{\alpha}^{\prime}$ 由式 (4.54) 给出。因此，将所有模式结合起来，平衡分布为：

关于这个结果的一个重要说明是，仅当对于所有 $\alpha$ 满足以下条件时，$\mathcal{P}_{\mathrm{eq}} ( Q_{\alpha} )$ 才有定义：

否则，朗之万迭代将产生无界的坐标值 $Q_{\alpha}$。因此，我们得出时间步长 $\Delta$ 必须满足以下条件：

其中 ${\bar{K}}_{D}$ 是力常数矩阵 K 的最大特征值。事实上，这是朗之万分子动力学的一个非常通用的结果，该方法仅当时间步长足够小时才稳定。

总之，需要强调以下几点：

_• 当 $\Delta \to 0$ 时，朗之万方程离散化的误差正确归零，因为在此极限下，精确分布由 $K_{\alpha}^{\prime} \to \bar{K}_{\alpha}$ 给出。注意，平衡分布确定中的相对误差，即弹簧常数 $\bar{K}_{\alpha}$ 的误差，与温度 $T$ 无关；因此，对于给定的目标精度，时间步长 $\Delta$ 可以保持与温度无关。

_• 在 $\Delta$ 为有限值时，离散化误差决定了（略为）不同的平衡分布，但该分布仍然保持相同的高斯形式。仅对于单一模式（即 $D = 1$），这种变化可以解释为有效温度的重新归一化，即 $\begin{array} {r} {T \to T / ( 1 - \frac{\bar{K}_{1} \Delta} {2} )} \end{array}$。

_• 从给定构型 R 生成一个独立构型 R' 所需的迭代次数 $n_{\mathrm{corr}}$ 由下式给出：

其中关联时间 τ 由式 (4.49) 确定，$\bar{K}_{\mathrm{cond}}$ 是矩阵 K 的条件数，即其最大与最小（非零）特征值之比。注意，对于谐势，$n_{\mathrm{corr}}$ 与温度无关。

4.5 加速朗之万动力学#

当条件数 $\bar{K}_{\mathrm{cond}}$ 很大时，矩阵是病态的，这意味着生成一个独立新构型所需的迭代次数极大。因此，该方法效率低下，需要采取一些技巧来加速算法 (Parisi, 1984)。一个简单的想法是在朗之万方程 (4.1) 中引入加速矩阵 S，使得：

其中噪声满足以下条件：

当 S 选取得尽可能接近 Hessian 矩阵 H 时，改进后的算法相对于原算法 (4.1)、(4.3) 和 (4.4) 有显著改进：

在下文中，我们考虑矩阵 S 不依赖于坐标 {R} 的情况，尽管在此意义上进行推广是可能的 (Mazzola and Sorella, 2017)。考虑方案 (4.61) 的优势可以通过利用低温极限来理解。事实上，当 $T \to 0$ 时，标准朗之万方程 (4.1) 退化为最速下降法，该方法从给定初始点 R 开始，提供了一种寻找势 V(R) 最近极小值的实用方法 (Press et al., 2007)。收敛所需的步数与条件数相关。相反，通过使用牛顿-拉夫逊方法，在谐情况下，无论条件数多大，我们都可以一步达到目标。因此，通过考虑 S 能紧密近似 Hessian 矩阵，由式 (4.61) 生成的动力学代表了牛顿-拉夫逊方法在有限温度下的推广，该方法最小化了热化时间/关联时间。

按照与 4.2 节相同的步骤，我们可以得到朗之万方程的离散化版本：

其中

之前的推导表明，在 $\Delta \to 0$ 的极限下，以下修正的福克-普朗克方程（modified Fokker-Planck equation）成立：

其平衡分布与之前相同，即式(4.8)中给出的玻尔兹曼分布（Boltzmann distribution）。

通过考虑谐振情形（harmonic case）可以体会到这种加速方案的效率。在足够低的温度下，热涨落被限制在势能最小值附近的构型，此时谐振是一个很好的近似。在这种情况下：

因此，通过选择 $\mathbf{S} = \mathbf{K}$，我们得到：

然后，取 $\Delta = 1$（如同牛顿-拉夫逊方法），式(4.65)大幅简化为：

这是一种生成新构型的高效方法，关联时间非常短，即 $n_{\mathrm{corr}} = 1$；实际上，上述方程根据下式设定一个新的独立构型：

注意到该平衡分布对应于谐振势能（式(4.44)）的正确（高斯）结果，只是温度项多了一个因子2，这是由式(4.70)中的 $\sqrt{2 T}$ 项引起的。这是时间离散化引入的近似的一个缺点。然而，在一般情况下，当势能非谐振时，该方法能够显著减少关联时间，是一种非常简单高效的方法，在进行一阶朗之万动力学分子动力学模拟之前值得考虑。图4.1展示了一个经典系统（64个粒子通过兰纳-琼斯势相互作用）的加速朗之万动力学示例，温度足够低，处于固相内部。

图4.1 对于由64个粒子通过兰纳-琼斯势相互作用的经典系统，展示了总能量随用于积分朗之万方程的离散化步长 $\Delta$ 的收敛情况。立方盒子边长为 $L = 7.8$，温度为 $T = 0.4$。报告了标准方案（式(4.14)，下三角）、加速动力学（式(4.65)，上三角）以及式(4.76)中 $A = 1$（实心圆）和 $A = 4$（实心方块）的结果。

对加速朗之万方程（式(4.61)）的离散化可以进一步改进，假设所选的加速矩阵S在平衡态下能很好地近似海森矩阵（Hessian） $\mathbf{K} \approx A \mathbf{S}$，只是存在一个整体常数A，可以通过经验调整以实现时间离散化中尽可能小的误差。因此，我们可以认为，即使势能不是谐振的，也有：

其中 $\mathbf{C} [ \mathbf{R} ( t ) ]$ 对 $\mathbf{R} ( t )$ 的依赖很弱——一旦我们忽略了平衡态下的非谐振贡献；那么，取 $\mathbf{C} [ \mathbf{R} ( t ) ] = A \mathbf{R}_{\mathrm{eq}}$ 是一个很好的近似。因此，我们转而考虑：

由于此微分方程是线性的，我们可以通过引入新变量 $\mathbf{R}^{\prime} ( t )$（满足 ${\bf R} ( t ) = {\bf X} ( t ) {\bf R}^{\prime} ( t )$，其中 $\mathbf{X} ( t )$ 是一个合适的含时矩阵）来精确积分（即在 $\Delta$ 下无误差）：

于是，在区间 $( t_{n} , t_{n + 1} )$ 上的结果为：

其中

该方案是对加速动力学的更精确积分，对于大的 $\Delta$ 值，它给出了式 (4.44) 中谐振势的精确动力学。事实上，当 $\Delta \infty$ 时，我们有：

这导致了温度 $T$ 下的正确玻尔兹曼分布（因为 $\langle z_{\alpha , n} z_{\beta , n} \rangle =$ $A K_{\alpha , \beta}^{- 1}$）：

在一般的非谐振情形中，我们可以将式 (4.76) 右端的第二项替换为 $\Delta_{A} \mathbf{S}^{- 1} \mathbf{f}_{n}$，从而得到迭代积分方案：

若非简谐项可忽略且 $\mathbf{K}$ 恰好等于 $\mathbf{A}\mathbf{S}$，则该式是精确的。注意，式 (4.81) 等价于在式 (4.65) 中将 $\Delta \to \Delta_{A}$、$T \to T \Delta_{2A} / \Delta_{A}$。实际中，若已知海森矩阵 $\mathbf{S}$ 的近似形式，可凭经验选择 $A$ 的值以最小化时间步长误差，并使用较大的 $\Delta$ 值。尽管对于一般势能 $V(\mathbf{R})$，误差在 $\Delta$ 上仍是线性的，但由于其仅依赖于非简谐项（前提是海森矩阵近似足够精确），误差通常很小。

图4：1 展示了用此方法处理 64 个经 Lennard-Jones 势相互作用的经典粒子所得结果。我们发现，$A=1$（对应精确海森矩阵）并非最优选择。这是因为在该情形下非简谐项很强，而经过重新标度的海森矩阵更能准确描述势能的整体形状。

第三部分

变分蒙特卡洛

变分蒙特卡洛 (Variational Monte Carlo)#

5.1 量子平均与统计采样#

本章将讨论变分蒙特卡洛方法的一般框架，以及其在格点上有相互作用的玻色子和费米子系统中的几个重要实现。采用此方法的主要优势在于第 1.4 节所展示的变分原理：任意量子态的能量总是受限于精确基态能量，这为我们提供了获得问题最佳可能解的途径。在多数情况下，通过对变分波函数进行合适的参数化，我们可以考虑多种不同的量子相（例如，金属、超导体和绝缘体）。通过优化参数，我们能达到能量最低的态，期望能抓住正确的基态行为。因此，变分波函数是一种灵活且有价值的方法，可用于深入了解那些无法用精确方法求解的模型的低能性质。相反，此方法的主要局限性在于它基于给定的形式 (Ansatz)，该形式可能引入所选参数化无法消除的相关偏差。

我们还注意到，变分波函数可以很容易地在一大类模型上定义和处理，无论相互作用的范围以及局域（即单点）希尔伯特空间的维数如何；该维数甚至可以是无限的，无需使用不可控的截断即可在有限维空间中工作。在这方面，变分蒙特卡洛优于其他方法，如密度矩阵重正化群 (density-matrix renormalization group) 或张量网络 (tensor-network) 方法，后者在相互作用范围和局域希尔伯特空间维数增加时，其复杂度会急剧增加 (White, 1992; Schollwöck, 2005, 2011)。

让我们从描述变分蒙特卡洛方法定义的一般框架开始。首先，我们在希尔伯特空间中固定一组完备基 $\{|x\rangle\}$，为简单起见，假设这些态是正交归一化的，满足：

那么，任意量子态 $|\Psi\rangle$ 可写为：

相应地，算子 $\mathcal{O}$ 在给定变分波函数 $|\Psi\rangle$ 上的期望值形式如下：

评估期望值的主要问题在于，求和中的组态数目随粒子数量呈指数增长。尽管通过采用少数守恒律（例如，粒子总数和/或沿量子化轴的总自旋分量的守恒）可以略微减小希尔伯特空间的维度，但其仍随粒子数量呈指数增长。因此，对于大系统，无法通过精确枚举所有组态来精确计算 $\langle\mathcal{O}\rangle$。然而，方程(5.3)可以改写成一种易于用标准蒙特卡洛方法处理的形式。实际上，我们有：

其中我们定义了算符 $\mathcal{O}$ 的局部估计量（local estimator）：

关键点在于

可以被解释为一个概率，因为它对于所有组态 $|x\rangle$ 都是非负的量，并且是归一化的，即 $\sum_{x} \mathcal{P} ( x ) = 1$。因此，计算算符 $\mathcal{O}$ 的量子平均值的问题可以转化为计算随机变量 $\mathcal{O}_{L} ( x )$ 在方程(5.6)给出的概率分布 $\mathcal{P} ( x )$ 下的平均值。特别地，如果我们考虑哈密顿量的期望值，局部估计量对应于所谓的局部能量（local energy），其定义为：

总之，可以定义一个随机算法（例如，马尔可夫过程），在其中生成一个组态序列 $\{\vert x_{n} \rangle \}$（例如，使用3.9节描述的梅特罗波利斯算法）。然后，在达到平衡后，这些组态根据所需的概率 $\mathcal{P} ( x )$ 分布。接着，量子期望值 $\langle\mathcal{O} \rangle$ 通过访问组态上随机变量 $\mathcal{O}_{L} ( x )$ 的平均值来评估：

最后，如3.11节所述，可以计算误差条。

5.2 零方差性质#

变分蒙特卡洛方法的一个重要特征是零方差性质（zero-variance property）。假设变分态 $| \Psi \rangle$ 恰好是 $\mathcal{H}$ 的一个精确本征态（不一定是基态），即 $\mathcal{H} | \Psi \rangle = E | \Psi \rangle$。那么，局部能量 $e_{L} ( x )$ 是常数：

因此，随机变量 $e_{L} ( x )$ 不依赖于 $| x \rangle$，这立即意味着其方差为零，而其平均值 $E$ 恰好等于精确本征值（换句话说，$e_{L} ( x )$ 不是一个随机变量）。显然，这是一个极端情况，对于一般关联模型来说非常罕见。然而，通常情况下，只要变分态 $| \Psi \rangle$ 趋近于一个精确本征态，$e_{L} ( x )$ 的方差就会减小。这一事实对于减少统计涨落和提高数值效率非常重要。零方差性质（zero-variance property）是量子期望值所独有的特征，而在经典计算中并不存在，因为经典计算中的可观测量具有热涨落。

最后，我们注意到局部能量平方的平均值对应于哈密顿量平方的量子平均值：

因此，随机变量 $e_{L} ( x )$ 的方差恰好等于哈密顿量在变分态 $| \Psi \rangle$ 上的量子方差：

5.3 Jastrow 和 Jastrow-Slater 波函数#

在量子平均值的蒙特卡洛评估中，参见方程 (5.4)，我们必须计算：

_• 不同构型概率的比值，这涉及到给定变分态与基组中两个构型之间交叠的比值：

这部分由 Metropolis 算法所需。

_• 局部估计量 ${\mathcal{O}}_{L} ( x )$，这进而涉及到基组态之间的交叠比值以及可观测量矩阵元的比值。例如，在考虑能量时，我们有：

直观看来，计算局域估计量似乎是一项极其困难的任务，因为它需要对多体希尔伯特空间中的所有状态进行求和；然而，得益于哈密顿量（或者类似的任何其他局域算符或关联函数）的局域性，实际上只有少数几项对求和有贡献。确实，给定组态$| x \rangle$，矩阵元$\langle x | \mathcal{H} | x^{\prime} \rangle$仅在$O ( L )$个组态$| x^{\prime} \rangle$上非零。以费米子哈伯德模型为例：使用局域基时，$| x \rangle$仅与少数其他组态相连，这些组态的差异在于一个电子从某个格点跳跃到其邻近格点；那么，此类过程的最大数量等于格点数$L$乘以键数再乘以2（考虑自旋）。因此，局域估计量的计算仅需少量操作，通常与格点/粒子数成正比。

如我们所见，变分蒙特卡洛方法的基本构建块是计算$\langle x | \Psi \rangle$，即变分态在基组中一般元素上的振幅。更精确地说，在整个马尔可夫过程中，只需计算这些重叠的比值。该计算必须在马尔可夫过程访问的每个组态上进行，因此必须尽可能快地完成。这给变分波函数的形式施加了某些约束。通常，费米子态需要计算行列式或Pfaffian，而玻色子波函数则需要计算永久（permanent）（矩阵$M$的永久与行列式的定义区别在于，不考虑排列的符号）。幸运的是，有快速（即多项式时间）算法来计算行列式和Pfaffian，因此我们可以将这些态视为电子系统的变分拟设（variational Ansätze）；相比之下，永久的计算仍然是一个指数困难的问题，这严重限制了在合理计算时间内能处理的玻色子数量。然而，所有玻色子都凝聚在单个态中的波函数不需要此类计算，并且可以在变分蒙特卡洛方法中高效处理。

在前述变分技术框架内，可以处理大量费米子和玻色子态，这些态可以在弱关联和强关联区域之间插值。接下来，我们将详细描述如何构造变分态以及如何计算马尔可夫过程中所需的基本构建块。我们将考虑第1章中介绍并讨论过的Jastrow（玻色子）或Jastrow-Slater（费米子）态：

其中$\left| {\Phi_{0}} \right.$是通用的无关联（玻色子或费米子）态，而$\mathcal{I}$是Jastrow因子，它考虑了电子关联。一个特别简单但重要的情况是$\mathcal{I}$包含密度-密度关联：

其中 $\nu_{i , j}$ 可取为变分参数，总数为 $L^{2}$。这里的赝势 $\nu_{i , j}$ 耦合了不同格点 i 和 j 上的密度，而非密度涨落（与式 (1.65) 不同）；对于粒子数守恒的系统，两种 Jastrow 因子的定义除一个无关的乘性因子外是一致的。对于平移不变模型，$\nu_{i , j}$ 仅取决于 i 和 j 之间的距离，因此参数数量可减少至 $O ( L )$。我们想指出的是，在蒙特卡洛方法中，可以精确处理（尽管仍存在统计误差）具有奇异 Jastrow 因子的极限情况，即 $\nu_{i , i} = \infty$，例如，消除所有多占据格点的 Gutzwiller 投影子（此情形与无限 Hubbard-U 相互作用相关）。实际上，这种情况可以通过构建一条仅访问满足该约束的构型 $| x \rangle$ 的马尔可夫链轻松实现。

在变分蒙特卡洛技术中考虑式 (5.14) 的波函数，其优势在于计算可以极其高效和快速。事实上，只要 Jastrow 因子在所选基下是对角的，我们就有：

其中 $\mathcal{I} ( x )$ 是针对构型 $| x \rangle$ 计算的 Jastrow 算子的值，即 $\mathcal{I} | x \rangle = \mathcal{I} ( x ) | x \rangle$。因此，给定电子构型，$\mathcal{I} ( x )$ 是一个数值，对于式 (5.15) 的 Jastrow 项，可以在 $O ( L^{2} )$ 次运算内求值。我们强调，为了实现多项式算法，Jastrow 因子必须只包含在基 $| x \rangle$ 下对角的算子，否则 ${\mathcal{I}} | x \rangle$ 会产生指数级数量的态，从而排斥对大型系统的任何计算。此外，$\langle x | \Phi_{0} \rangle$ 也可以在 $O ( L^{3} )$ 次运算内（对于费米子）或 $O ( 1 )$ 次运算内（对于凝聚玻色子）轻松求值，见第 5.5、5.6 和 5.7 节。

5.4 基组的选择#

在大多数情况下（但并非强制），基组选择为实空间中的乘积态。在具有 L 个格点的晶格上，我们首先必须对它们进行排序，$i = 1 , \ldots , L$，然后在每个格点上选取一个基组 $\{| \xi \rangle_{i} \}$，其维度取决于所考虑的模型。例如，在单带 Hubbard 模型中，$| \xi \rangle_{i}$ 可以是以下四种态中的任意一种：

即，格点 i 上的空态、两种单占据态（自旋沿量子化轴 $z$ 向上或向下）、以及双占据态。

对于玻色 Hubbard 模型，局部希尔伯特空间是无穷维的，同一个格点上可以有任意数量的玻色子：

最后，对于自旋为 S 的海森堡模型，局部基可取为 $S_{i}^{z} = - S , \ldots , S$ 这 $2 S + 1$ 个态：

一旦确定了局部希尔伯特空间，整个晶格基组 $\{\vert x \rangle \}$ 的一个一般元素由乘积态给出：

为后续使用，我们还定义了整个晶格的真空态为：

在格点上的每个位置取局部基矢的优势在于，哈密顿量通常是空间局域的，其矩阵元易于计算：相互作用项往往是对角的，例如Hubbard-U相互作用，或者是局域的，例如多带模型中的Hund耦合；而跳跃项通常是短程的，这导致对每个多体构型 $|x\rangle$ 而言，非零的矩阵元数量很少（即与粒子数成正比）。

在蒙特卡洛采样（即马尔可夫过程）过程中，方便的做法是从通过应用哈密顿量得到的那些构型中生成新的（提议）构型。这样，第 n 步的通用构型 $\left| x_{n} \right.$ 可能仅因一个整体符号（对于费米子系统）而与基组中的构型不同。在计算局域估计量 ${\mathcal{O}}_{L} ( x )$（例如局域能量 (5.7)）时，也会出现类似的问题。只要我们记录下可能出现的符号变化，这并不会构成问题。考虑在 L 个格点上有 $N_{e}$ 个无自旋电子的情况，构型

可以用一个长度为 L 的向量 κ 来描述，其第 i 个元素（若非零）表示在定义该构型的上述字符串中，对应产生算符的位置（从左到右）。例如，取 $N_{e} = 3$ 且 $L = 6$，构型 $c_{3}^{\dagger} c_{1}^{\dagger} c_{6}^{\dagger} | 0 \rangle$ 对应向量 (2, 0, 1, 0, 0, 3)。那么，应用跳跃项 $c_{R_{l}^{\prime}}^{\dagger} c_{R_{l}}$ 得到：

这表明字符串中唯一的变化出现在位置 $R_{l}$ 和 $R_{l}^{\prime}$。在上面的例子中，通过应用 $c_{2}^{\dagger} c_{1}$，我们得到 $c_{3}^{\dagger} c_{2}^{\dagger} c_{6}^{\dagger} | 0 \rangle$，向量 κ 变为 (0, 2, 1, 0, 0, 3)。在马尔可夫链中，存储并记录 κ 对于正确定义多体态是很有用的，参见第 5.6 和 5.7 节。

5.5 玻色子系统#

此处，我们详细描述如何高效计算不同构型上的波函数比值，以应用 Metropolis 算法，见方程 (5.12)。首先，我们考虑非相互作用部分 $\langle x | \Phi_{0} \rangle$ 的计算，然后讨论来自 $\mathcal{I} ( x )$ 的贡献。

5.5.1 非相互作用态的定义#

对于玻色子的非相互作用波函数，自然的选择是将所有粒子凝聚到一个单一态上。最简单的例子是考虑由下式定义的单粒子“轨道”态：

其中 $\{V_{i , \alpha} \}$ 是玻色子位于格点 $i$ 的泛振幅；$\alpha$ 是指定轨道的指标：例如，玻色子处于零动量态的情况由 $V_{i , \alpha} = 1 / \sqrt{L}$ 给出。那么，有 $N_{b}$ 个玻色子占据同一轨道的多体态为：

展开求和，我们可以将多体态重写为基于方程 (5.18) 的基组形式：

这使我们能够非常简单地表达重叠积分 $\langle x | \Phi_{0} \rangle$，其中构型 $| x \rangle \equiv | n_{1} , ..., n_{L} \rangle$ 在格点 $i$ 上包含 $n_{i}$ 个玻色子：

进而得到：

5.5.2 无相互作用态（Non-Interacting State）的快速计算#

对于无相互作用态，其形式为 $\langle n_{1} , \ldots , n_{L} | \Phi_{0} \rangle$ 由式(5.28)给出，两个不同构型的重叠比可以通过简单计算得到。例如，考虑仅因玻色子在格点 l 和 k 之间跳跃而不同的两个构型，即 $| x \rangle \equiv | n_{1} , \ldots , n_{k} , \ldots , n_{l} , \ldots , n_{L} \rangle$ 和 $| x^{\prime} \rangle \equiv | n_{1} , \dots , n_{k} + 1 , \dots , n_{l} - 1 , \dots , n_{L} \rangle$，我们有：

因此，在马尔可夫链中考虑单玻色子跳跃过程时，波函数无相互作用部分对 Metropolis 比率 (5.12) 的贡献非常简单，只需 $O ( 1 )$ 次运算即可完成。移动多个玻色子的情况可以通过直接推广上述分析得到。

我们想强调的是，如果玻色子不是凝聚的，而是占据不同的轨道，那么计算 $\langle x | \Phi_{0} \rangle$ 会变得极其复杂：

事实上，在这种情况下，重叠量由一个 $N_{b} \times N_{b}$ 矩阵的积和式（permanent）给出：

其中 $R_{j} {\mathrm{' s}} ( {\mathrm{with}} j = 1 , \dots , N_{b} )$ 取占据格点的值。注意，根据每个格点的占据数，$R_{j}$ 的值可能重复出现，例如 $\begin{array} {r} {\mathrm{e} . \mathrm{g} . , R_{j} = k} \end{array}$ 会出现 $n_{k}$ 次。计算积和式（或两个积和式的比）在计算上是指数级困难的，因此阻碍了我们对大 $N_{b}$ 值进行计算。

5.5.3 Jastrow 因子（Jastrow Factor）的快速计算#

现在考虑 Jastrow 因子的贡献。直接计算 $\mathcal{I} ( x )$ 需要 $O ( L^{2} )$ 次运算来得到式(5.15)的 Jastrow 因子，因此在计算 Metropolis 比率 (5.12) 中的比值 $\mathcal{I} ( \boldsymbol{x}^{\prime} ) / \mathcal{I} ( \boldsymbol{x} )$ 时具有相同的复杂度。然而，当两个构型仅由少数玻色子跳跃不同时，可以对比值进行快速计算，其运算量为 $O ( L )$。为此，考虑：

这里我们利用了赝势（pseudo-potential）$\nu_{i , j}$ 的对称性，即 $\nu_{i , j} = \nu_{j , i}$，和平移不变性，即 $\nu_{k , k} = \nu_{l , l}$。式(5.32)右边第二项与玻色子构型无关，可以在模拟开始时一次性计算。而第一项依赖于玻色子构型，该构型是通过马尔可夫链采样的。一旦我们计算并存储一个依赖于 $| x \rangle = | n_{1} , \dots , n_{k} , \dots , n_{l} , \dots , n_{L} \rangle$ 的 L 维向量，该比值的计算可以在 $O ( 1 )$ 次运算内完成：

然后，一旦新的候选态 $| x^{\prime} \rangle = | n_{1} , \dots , n_{k} + 1 , \dots , n_{l} - 1 , \dots , n_{L} \rangle$ 被接受，就必须更新向量 $T_{\mathrm{Jastrow}} ( j )$，这需要 $O ( L )$ 次运算。实际上，由于 $n_{i}^{\prime} = n_{i} + \delta_{i , k} - \delta_{i , l}$，我们有：

因此，我们可以在马尔可夫链开始时对所有格点 $j = 1 , \dots , L$ 从头计算 $T_{\mathrm{{Jastrow}}} ( j )$，然后在马尔可夫过程中每当接受一个新构型时，使用式 (5.34) 对其进行更新。作为安全习惯，建议定期从头重新计算 $T_{\mathrm{Jastrow}} ( j )$，因为快速更新的舍入误差会累积并引起数值误差。

总之，通过使用这种针对 Jastrow 因子的方案，单次移动算法的复杂度为 $O ( L )$。然而，通过进行这些局域更新，我们需执行 $O ( L )$ 次移动才能得到一个几乎独立的构型，这对于具有有限关联时间的系统而言。因此，玻色子代码的计算复杂度与格点/玻色子数目成二次关系。

5.6 行列式费米子系统#

现在转向由 Jastrow-Slater 波函数描述的费米子系统。由于 Jastrow 项的计算部分与之前的玻色子情况类似，我们在此不再重复。相同的快速更新方法也可用于费米子。相反，涉及无相互作用态 $\langle x | \Phi_{0} \rangle$ 计算的部分与玻色子情况完全不同。下面，我们描述如何设计一个高效算法来处理这一对象。

5.6.1 无相互作用态的定义#

任何无相互作用波函数都可以作为某个合适的二次型哈密顿量 $\mathcal{H}_{0}$ 的基态获得。首先，我们将自旋指标 $\sigma$ 和晶格格点 i 合并成一个从 1 到 2L 的单一指标 I：

然后，我们从最简单的无相互作用哈密顿量形式开始：

该哈密顿量仅包含跳跃项，也包括自旋沿 z 方向不守恒的过程，即 $I \ \leq \ L$ 且 $J \ > \ L$ 以及反之的项。简洁起见，式 (5.37) 的无相互作用哈密顿量可以写为：

其中

且

由于 $\mathcal{H}_{0}$ 与总电子数 $\begin{array} {r} {N_{e} = \sum_{I} {d_{I}^{\dagger} d_{I}}} \end{array}$ 对易，其本征态为单粒子轨道。实际上，$2 L \times 2 L$ 的矩阵 T 可以使用标准库（例如 LAPACK 例程）轻松对角化：

其中 U 是幺正矩阵（保持费米子算符反对易关系不变），$\mathbf{E} = \mathrm{diag} ( \varepsilon_{1} , \dots , \varepsilon_{2 L} )$ 是包含 $\mathbf{T}$ 的 2L 个本征值 $\varepsilon_{\alpha}$ 的对角矩阵，且 $\Phi^{\dagger} = ( \phi_{1}^{\dagger} , \dots , \phi_{2 L}^{\dagger} )$ 由 T 的本征向量定义：

现在，多体态 $\left| {\Phi_{0}} \right.$ 可以通过占据 $N_{e}$ 个最低能级轨道来构建：

根据第5.4节的讨论，在马尔可夫过程中访问到的、包含 $N_{e}$ 个电子的通用位型（generic configuration）为：

其中 $j = 1 , \ldots , N_{e}$ 同时包含自旋向上和向下的电子，$R_{j}$ 的取值范围为 1 到 2L：自旋向上电子的位置与格点编号一致，而自旋向下电子的位置则需要偏移 L。重叠积分 $\langle x | \Phi_{0} \rangle$ 于是给出为：

中方括号内的求和针对 $| x \rangle$ 中 $\{R_{j} \}$ 的所有可能排列；符号的出现源于费米子算符的反对易关系。于是，我们得到：

在构建多体态（5.43）时，我们必须注意构建一个唯一的多体态，即占据正确的能量最低轨道。当最高占据轨道与最低未占据轨道具有不同能量时（闭壳层位型，closed shell configuration），选择是唯一的。相反，也可能出现简并情况，导致无法唯一选择（开壳层位型，open shell configuration）。每当对无相互作用哈密顿量 $\mathcal{H}_{0}$ 进行数值对角化时，标准库提供的本征态不带确定的量子数（如动量），而是由简并轨道的通用线性组合给出（这种组合会不规则地依赖于计算所用的数值精度）。只要所有简并本征态都包含在多体态（5.43）中，这就不是问题，因为对矩阵 U 中列的任何线性组合都不会改变式（5.46）中行列式的值。相反，如果不包含所有简并本征态，就会产生问题，因为行列式将取决于所考虑的具体态的组合。因此，每当进行数值对角化时，我们必须验证是否出现了闭壳层位型。否则，可以通过构造具有确定量子数的合适轨道来克服能隙消失的问题，从而获得多体态的可重复模拟。

在实践中，无相互作用哈密顿量 $\mathcal{H}_{0}$ 的对角化必须在蒙特卡洛计算开始时进行；然后，我们需要存储 U 的缩减部分，即只保留对应于占据轨道的 $N_{e}$ 列。这是一个 $2 L \times N_{e}$ 矩阵：

然后，与一般电子位型的交叠由仅考虑对应电子位置 $\{R_{j} \}$ 的行（5.47）所得到的矩阵的行列式给出，从而得到一个 $N_{e} \times N_{e}$ 矩阵。

值得注意的是，如果无相互作用哈密顿量包含一个耦合自旋向上和向下的电子配对，则同样形式的框架也可使用。事实上，让我们考虑由下式描述的BCS哈密顿量：

其中我们引入了一个化学势 $\mu_{0}$ ，它在平均上固定了电子数。在这种情况下，总粒子数不守恒，单个粒子轨道的概念不再适用。实际上，BCS哈密顿量的基态自然是用配对函数（pairing function）来写的（见第1章）。尽管如此，我们可以对自旋向下的电子进行粒子-空穴变换：

然后，除常数项外，变换后的BCS哈密顿量具有方程（5.38）的形式，其中矩阵T由下式给出（i 和 $j \leq L )$：

由于在粒子-空穴变换后，粒子数（而非自旋的z分量）是守恒的，BCS哈密顿量的本征态可以表达为“轨道”形式，类似于方程（5.42）中的轨道，但沿z轴没有确定的自旋分量。

注意，粒子-空穴变换（5.49）和（5.50）将自旋向下态变为真空，反之亦然，因此：

其中 $| \tilde{0} \rangle_{i}$ 是f电子的真空，即 $f_{i , \uparrow} | \tilde{0} \rangle = f_{i , \downarrow} | \tilde{0} \rangle = 0$ 。通过使用费米子算符的反对易关系，新费米子算符的占据数 $u_{i , \sigma}$ 与原算符的占据数关系如下：

这意味着在位点i上的密度和磁化强度之间的关系由下式给出：

其中我们引入了 $u_{i}$ 和 $w_{i}^{z}$，分别作为变换后算符的密度和磁化强度。选择字母 u 和 w 是因为它们看起来分别像是倒过来的 n 和 m。因此，除了常数项之外，通过 Eq. (5.49) 和 Eq. (5.50) 的粒子-空穴变换（particle-hole transformation），密度变成了自旋的 z 分量，反之亦然。作为一个推论，我们得到，每当原始粒子的总磁化强度为零时，新粒子的总数等于格点数目 L。

5.6.2 快速计算行列式#

现在让我们展示，当两个组态 $| x \rangle$ 和 $| x^{\prime} \rangle$ 仅相差一个或少数几个电子跃迁时，如何高效地计算行列式之比。根据 Eq. (5.46)，重叠积分 $\langle x | \Phi_{0} \rangle$ 由 $\tilde{U}_{j , \alpha} \equiv U_{R_{j} , \alpha}$ 的行列式给出

该矩阵是通过仅取矩阵 U (Eq. (5.47)) 中与被占据格点对应的行得到的。注意，$\tilde{\mathbf{U}}$ 中行的实际顺序由 |x-( 中产生算符的顺序决定，见 Eq. (5.44)。在这方面，对于有自旋情况，向量 κ（见第5.4节）的推广可用于将 $\tilde{\mathbf{U}}$ 的行按正确顺序排列。

让我们开始考虑电子组态，其中 $| x^{\prime} \rangle$ 是通过将第 l 个电子从 $R_{l}$ 跃迁到 $R_{l}^{\prime}$ 从 |x- 得到的，即 $| x^{\prime} \rangle = d_{R_{l}^{\prime}}^{\dagger} d_{R_{l}} | x \rangle$

新的矩阵 $\tilde{\mathbf{U}}^{\prime}$ 将等于 ${\tilde{\mathbf{U}}}$，但第 l 行的元素从 $U_{R_{l} , \alpha}$ 变为 $U_{R_{l}^{\prime} , \alpha}$：

那么，仅相差一个单费米子跃迁的两个组态之比为：

用 K 表示第 l 个电子的新格点（即 $K \equiv R_{l}^{\prime}$），更新后的矩阵元可以写成紧凑形式：

其中我们定义了 $\nu_{\alpha}^{K , l} ~ \equiv ~ U_{K , \alpha} ~ - ~ \tilde{U}_{l , \alpha}$；这里索引 K 和 l 是固定的，因为它们分别指定了电子跃迁到的格点和电子索引。该方程可以重写为如下形式：

其中

因此，$\tilde{\mathbf{U}}^{\prime} = \tilde{\mathbf{U}} \mathbf{Q}$，这意味着计算 $\tilde{\mathbf{U}}^{\prime}$ 和 $\tilde{\mathbf{U}}$ 的行列式之比等价于计算 Q 的行列式：

极大的简化来源于矩阵 Q 的行列式可以很容易地计算。实际上，Q 具有特别简洁的形式，可以写成：

其中 $\boldsymbol{B_{\beta}} = \boldsymbol{\tilde{U}}_{\beta , l}^{- 1}$，$\mathcal{A}_{\alpha} = \nu_{\alpha}^{K , l}$。虽然该矩阵不是厄米矩阵，但矩阵 Q 的特征值可以从久期方程（secular equation）得到：

利用式 (5.71) 的显式形式，我们得到：

这意味着所有与 $\mathcal{A}_{\alpha}$ 正交的向量 $\upsilon_{\alpha}$ 都是特征值为 $\lambda = 1$ 的特征向量（共有 $N_{e} - 1$ 个这样的向量）；此外，$\upsilon_{\alpha} = B_{\alpha}$ 也是一个特征向量，其特征值为 $\begin{array} {r} {\lambda = 1 + \sum_{\alpha} \mathcal{A}_{\alpha} \mathcal{B}_{\alpha} = \sum_{\alpha} U_{K , \alpha} \tilde{U}_{\alpha , l}^{- 1}} \end{array}$。因此，我们有：

在模拟开始时存储了针对组态 |x- 的矩阵 $\tilde{\mathbf{U}}^{- 1}$ 后，此计算需要 $O ( N_{e} )$ 次运算，而不是计算行列式所需的 $O ( N_{e}^{3} )$ 次。然后，一旦新组态 $\vert x^{\prime} \rangle$ 在马尔可夫过程中被接受，就必须更新矩阵 $\tilde{\mathbf{U}}^{- 1}$。这可以在 $O ( N_{e}^{2} )$ 次运算中完成。事实上，我们有 $( \tilde{\mathbf{U}}^{\prime} )^{- 1} = \mathbf{Q}^{- 1} \tilde{\mathbf{U}}^{- 1}$，而矩阵 Q 的逆矩阵由下式给出（易于验证）：

那么，$( \tilde{\mathbf{U}}^{\prime} )^{- 1}$ 的更新后的矩阵元由下式给出：

这是一个用于更新矩阵 $\tilde{\mathbf{U}}^{- 1}$ 的封闭方程。

我们要强调，前述关于行列式比值计算和更新的结果可以进一步简化。实际上，在计算开始时，我们可以计算并存储一个 $2 L \times N_{e}$ 的矩阵 W，其元素由下式给出：

那么行列式比值 (5.74) 只需要 $O ( 1 )$ 次运算，因为它等价于取新格点（行）与执行跳跃过程的电子（列）所对应的元素：

W 的求值需要知道完整的矩阵 U（该矩阵已在模拟开始时计算并存储一次，不依赖于电子组态）和 $\tilde{\mathbf{U}}^{- 1}$（它依赖于组态 $| x \rangle$）。然后，可以给出一个简单的 W 更新方案。事实上，将式 (5.76) 两边乘以 $U_{I , \alpha}$ 并对 α 求和，我们得到：

其中，我们用到 det ${\bf Q} = {\cal W}_{K , l}$，这根据式 (5.78) 和式 (5.77) 的定义得出。由于每个矩阵元都需用 O(1) 次运算更新，总成本为 $O ( 2 L N_{e} )$。

通过使用由填充单粒子态构建的波函数（也包括进行粒子-空穴变换后具有BCS配对的情况），该算法可以用矩阵元 W 来表示，W 正是等时格林函数（equal-time Green's function），见式 (5.66)：

我们想要指出，对于更新后的（等时）格林函数，式(5.79)可以利用威克定理（Wick theorem）（具有不同的左矢和右矢）直接解释：

在这种情况下，我们必须考虑那些不包括反常（超导）收缩的收缩，这些收缩会消失，因为$| x \rangle$和$\left| {\Phi_{0}} \right.$都有固定的粒子数。

让我们通过将先前形式推广到多于一个电子跳跃的情况来结束这一部分，即$| x^{\prime} \rangle = d_{R_{l_{1}}^{\prime}}^{\dagger} d_{R_{l_{1}}} \cdot \cdot \cdot d_{R_{l_{m}}^{\prime}}^{\dagger} d_{R_{l_{m}}} | x \rangle$，从而导致对$\tilde{\mathbf{U}}$矩阵的m行进行修改；例如，对跳跃或自旋翻转过程就是这种情况。于是，式(5.67)被推广为：

和之前一样，指标$K_{r}$和$l_{r}$（对于$r = 1 , \ldots , m$）是固定的，因为它们分别指定了电子跳跃的位点和电子的指标。通过执行与之前相同的代数运算，我们得到：

其中矩阵Q现在具有以下形式：

其中$\begin{array} {r} {\mathcal{B}_{\beta}^{r} = \tilde{U}_{\beta , l_{r}}^{- 1}} \end{array}$和$\mathcal{A}_{\alpha}^{r} = \nu_{\alpha}^{K_{r} , l_{r}}$。和之前一样，$\mathbf{Q}$的行列式可以通过求解相应的本征值问题轻松计算：

这意味着所有与由$\mathcal{A}_{\alpha}^{r} \mathrm{{^{* {\Omega}_{\alpha}}}}$定义的子空间正交的向量都是$\lambda = 1$的本征向量；此外，$\begin{array} {r} {\upsilon_{\alpha} = \sum_{r = 1}^{m} x_{r} \mathcal{B}_{\alpha}^{r}} \end{array}$是一个本征向量，前提是系数$x_{r}$满足：

因此，$\mathbf{Q}$的m个非平凡本征值由$m \times m$矩阵的本征值给出：

两个行列式之比的最终表达式由下式给出：

同样在这种情况下，一旦移动被接受，我们必须更新矩阵$( \mathbf{\tilde{U}}^{\prime} )^{- 1} =$ $\mathbf{Q}^{- 1} \tilde{\mathbf{U}}^{- 1}$。和之前一样，可以求出Q矩阵的逆矩阵：

因此，我们得到：

然后，将前一方程的两边乘以 $U_{i , \alpha}$ 并对 α 求和，得到更新后的 $\mathbf{W}^{\prime}$：

其中

与Jastrow因子的快速更新类似，强烈建议定期从头重新计算矩阵W。其计算开销为 $O ( 2 L N_{e}^{2} )$，因此如果每隔 $O ( N_{e} )$ 次电子移动执行一次，该计算不会影响算法的复杂度。

5.6.3 延迟更新#

考虑只有一个电子被移动的简单更新。更新中的基本操作是所谓的秩1更新（rank-1 update），作用于一个 $2 L \times N_{e}$ 的矩阵，见式(5.79)，可写为：

其中 a 和 b 分别是长度为 2L 和 $N_{e}$ 的两个向量：

当矩阵规模较大以至于W无法放入处理器缓存时，这种更新操作在计算上可能效率低下。克服这一缺点的一种方法是在不丢失信息的情况下延迟矩阵W的更新。这可以通过存储一组左向量 $\mathbf{a}^{( p )}$ 和右向量 $\mathbf{b}^{( p )}$（$p = 1 , \ldots , m$），以及开始延迟更新的“初始”矩阵（记为 $\mathbf{W}^{0}$）来实现。于是，经过 m 次更新后的矩阵W为：

每当我们接受一个新组态，可以用 $\mathbf{W}^{0}$ 和之前的 $p = 1 , \ldots , m$ 的向量，通过少量运算计算出新的一对向量 $\mathbf{a}^{( m + 1 )}$ 和 $\mathbf{b}^{( m + 1 )}$，只需将式(5.96)代入式(5.94)和(5.95)即可。一旦矩阵W写成式(5.96)的形式，计算求和中因子所需的运算次数为 $O [ m ( 2 L + N_{e} ) ]$，对于 $m \ll L$ 的情况，这与完全更新相比可以忽略不计。

这样，我们可以找到一个最优的 $m_{\mathrm{max}}$，通过标准矩阵乘法计算完整矩阵W：

其中 A 和 B 分别是 $2 L \times m_{\mathrm{max}}$ 和 $N_{e} \times m_{\mathrm{max}}$ 的矩阵，分别由 $p = 1 , \ldots , m_{\mathrm{max}}$ 个向量 $\mathbf{a}^{( p )}$ 和 $\mathbf{b}^{( p )}$ 构成。之后，我们可以将整数 m 重新初始化为零，以新的 $\mathbf{W}^{0} = \mathbf{W}$ 继续一个新的延迟更新。这种更新过程的好处在于，经过 $m_{\mathrm{max}}$ 步的一个周期后，大部分计算量来自式 (5.97) 中的矩阵-矩阵乘积，这比式 (5.79) 中给出的对 W 进行 $m_{\mathrm{max}}$ 次秩1更新（rank-1 update）要高效得多（并且不受缓存限制）。对于大量电子，采用最优 $m_{\mathrm{max}}$ 值的延迟更新过程，可以将变分蒙特卡洛（variational Monte Carlo）代码的速度提升大约一个数量级。

5.6.4 反流关联（Backflow Correlations）#

在这里，我们将讨论在存在已引入格点问题的反流关联（backflow correlations）时，如何实现行列式部分的更新。我们注意到，在格点情形下，无法实现标准的反流关联（1.70）所涉及的连续坐标变化。因此，这种 ansatz 需要根据允许的格点位置，找到与轨道泰勒展开（Tocchio 等人，2008 和 Tocchio 等人，2011）相一致的通用离格反流坐标中轨道的表达式。实际上，我们考虑一个二次型哈密顿量（如第 5.6 节中讨论的那些），来构造无相互作用轨道 $\{U_{I , \alpha} \}$，见式 (5.42)。在最简单的处理方法中，反流关联作用于空穴-双占（holon-doublon）（最近邻）对上（Tocchio 等人，2008），通配构型 |x- 与反流波函数 $| \Phi_{0}^{b} \rangle$ 之间的重叠由带有反流修正的“关联”轨道构造：

其中，对于自旋向上（向下）的电子，$I = i ( I = i + L )$，J 和 $j$ 类似；$\langle j \rangle_{i}$ 表示与 i 最近邻的格点 j；$D_{i} \left( H_{i} \right)$ 是算子，若格点 i 是双占（空位）则给出 1，否则为 0；最后 $\eta_{0}$ 和 $\eta_{1}$ 是变分参数。然后，给定构型 $| x \rangle$ 的波函数通过取以下矩阵的行列式得到：

其中行对应电子占据的格点。

对于两个行列式之比，其构型分别为 $| x \rangle$ 和 $| x^{\prime} \rangle$，且相差一个或少数几个电子跃迁时，其快速更新可以像之前一样计算。

现在，唯一的区别在于，即使只有一个电子跃迁，新矩阵 $\tilde{\mathbf{U}}^{b \prime}$ 不仅会改变一行，而是会改变少数几行。实际上，移动一个电子可能会产生一个空位或双占位，从而导致式 (5.98) 中“关联”轨道的修改。因此，一般来说，我们必须使用通用更新，即改变多行；然而，在这种情况下，直接处理逆矩阵 $\tilde{\mathbf{U}}^{\prime}$ 更为方便，可以根据式 (5.90) 对其进行更新。

在实践中，我们可以考虑更一般的“关联”轨道形式，加入更多的项 (Tocchio et al., 2011)，但依然遵循相同的思路，即考虑依赖于多体构型 $| x \rangle$ 的无相互作用轨道的线性组合（在前一种情况下，对于近邻格点上有空穴和双占 (holon-doublon) 的构型，需要取这样的线性组合）。

5.7 基于 Pfaffian 的费米子系统#

在上一节中，我们解释了如何从由式 (5.37) 哈密顿量描述的无相互作用电子出发构建关联波函数，这自然引出了单粒子轨道的概念。这里，多体态的无关联部分由行列式给出，见式 (5.46)。通过对自旋向下的电子进行粒子-空穴变换，同样的方法也可以推广到包含不同自旋电子配对的 BCS 哈密顿量，见式 (5.48)。相反，当配对项耦合同向自旋的电子时，即使经过粒子-空穴变换，也无法定义轨道，此时无关联态用配对函数表示。进而，行列式并不代表最一般的无关联反对称态，而后者可以通过反对称矩阵的 Pfaffian 来定义。下面我们将证明，在最一般情况下，波函数的无关联部分由 Pfaffian 给出，它作为特例可以描述单重态配对态或 Slater 行列式。这种平均场波函数最早由 Bouchaud 等人 (1988) 引入，随后用于研究小原子和分子 (Bajdich et al., 2006, 2008)。我们将描述这种态以及进行局部更新的高效算法。

5.7.1 Pfaffian 的定义#

在讨论无相互作用态之前，我们先总结一下 Pfaffian 的定义和几个重要性质。考虑一个 $2 n \times 2 n$

的斜对称矩阵 A，即满足 $a_{i , j} = - a_{j , i}$ 的矩阵。A 的 Pfaffian 定义为反对称乘积：

其中求和遍历所有 $( 2 n - 1 ) ! !$ 种由 $\alpha \quad =$ $\{( i_{1} , j_{1} ) , \dotsc , ( i_{n} , j_{n} ) \}$ 定义的对分割，且满足 $i_{k} ~ < ~ j_{k}$ 和 $i_{1} < i_{2} < \cdots < i_{n}$；sign(α) 表示分割 α 对应的置换的奇偶性。例如，对于 $n = 2$

注意 det $\mathbf{A} = ( \mathrm{Pf} \mathbf{A} )^{2}$；因此，除了一个整体符号外，矩阵 Pfaffian 的计算可以归结为行列式的计算。

Pfaffian 满足以下关系式：

其中 A 和 $\mathbf{A}^{\prime}$ 分别为 $2 n \times 2 n$ 和 $2 m \times 2 n$ 的斜对称矩阵。此外：

其中 B 是任意 $2 n \times 2 n$ 矩阵。

最后，我们报告两个重要恒等式，将用于演示 Pfaffian 的快速更新。考虑两个可逆斜对称矩阵 A $( 2 n \times 2 n )$ 和 $\mathbf{C_{\lambda}} ( 2 m \times 2 m )$ 以及一个实矩阵 B $( 2 n \times 2 m )$。则我们有：

因此，利用式 (5.102) 和 (5.103)，我们得到：

5.7.2 非相互作用态的定义#

一种无法回溯到基于单粒子轨道之前形式化方法的情况是，当原始非相互作用哈密顿量包含同自旋配对，或同时存在自旋翻转跳跃和配对项时。确实，这里无法通过进行粒子-空穴变换来消除配对项。费米子系统最一般的二次型哈密顿量由广义 BCS 哈密顿量给出：

其中 $t_{i , j}^{\sigma , \tau} = \left( t_{j , i}^{\tau , \sigma} \right)^{*}$ 和 $\Delta_{j , i}^{\tau , \sigma} = - \Delta_{i , j}^{\sigma , \tau}$ 分别为跳跃项和配对项。

哈密顿量 $\mathcal{H}_{0}$ 可通过广义博戈留波夫变换（Bogoliubov transformation）对角化：

其中 $2 L \times 2 L$ 矩阵 u 和 v 通过施加以下条件确定：

以及正交性条件（该条件保持费米子算符的反对易关系）：

注意，式 (5.108) 允许成对的解，因为如果 $\Phi_{\alpha}$ 具有本征值 $\varepsilon_{\alpha}$，则 $\Phi_{\alpha}^{\dagger}$ 也将是带有本征值 $- \varepsilon_{\alpha}$ 的解。因此，我们可以限制只考虑具有非负本征值的本征向量。

在大多数情况下，非相互作用哈密顿量 (5.106) 的基态可以写成广义 BCS 波函数：

其中 $F_{I , J}$ 是反对称配对函数（pairing function），通过要求 $\left| {\Phi_{0}} \right.$ 被所有 $\varepsilon_{\alpha} \geq 0$ 的算符 $\Phi_{\alpha}$ 湮灭来定义，这导致：

与行列式情况类似，这里我们也必须注意定义唯一的多体波函数。因此，只应考虑没有零本征值 $\varepsilon_{\alpha}$ 的情况。此外，为了得到非奇异的配对函数 F，我们还必须要求矩阵 u 可逆。注意，奇异的 u 对应于存在占据的单粒子轨道。在这种情况下，多体态的构建必须包含这些未配对的轨道，这增加了算法的计算复杂度 (Bajdich et al., 2006, 2008)。在下面的快速更新中，我们将局限于非奇异情况。

$| \Phi_{0} \rangle$ 与式 (5.44) 的通用位型之间的重叠由下式给出：

然后，展开乘积的和，我们得到对 $| x \rangle$ 中费米子的 $( N_{e} - 1 ) ! !$ 种对配对的求和，其中每个贡献因子有 $( N_{e} / 2 ) !$ ！个恒等项。费米子算符的排序产生了定义 Pfaffian 所需的符号，见式 (5.100)，因此得到：

5.7.3 快速计算 Pfaffian#

现在我们展示如何有效计算两个 Pfaffian 的比值，其中两个位型 $| x \rangle$ 和 $| x^{\prime} \rangle$ 相差一个单电子跳跃。特别地，假设 $R_{l}$ 处的第 l 个电子改变其位置到 $R_{l}^{\prime}$，即 $| x^{\prime} \rangle = d_{R_{l}^{\prime}}^{\dagger} d_{R_{l}} | x \rangle$ 。用于计算与 $| x^{\prime} \rangle$ 重叠的新矩阵 $\mathbf{F^{\prime}}$ 仅在第 l 行和第 l 列与原矩阵 F 不同。根据式 (5.114)，重叠 $\langle x | \Phi_{0} \rangle$ 由以下矩阵的 Pfaffian 给出：

该矩阵可由原始矩阵 F 构建。与行列式情况一样，行（和列）的实际顺序由位型 $| x \rangle$ 中产生算符的顺序决定，见式 (5.44)。然后，可以将自旋情况下的向量 κ 的推广（见 5.4 节）用于选择正确的行列顺序。那么，重叠 $\langle x^{\prime} | \Phi_{0} \rangle$ 由以下矩阵的 Pfaffian 给出：

在紧凑形式下，我们有：

这可以化为形式 $\tilde{\mathbf{F}}^{\prime} = \tilde{\mathbf{F}} + \mathbf{BCB}^{T}$ ，其中 B 是一个 $N_{e} \times 2$ 矩阵：

以及 $\mathbf{C}$ 是一个 $2 \times 2$ 矩阵：

因此，利用等式(5.105)，我们得到Pfaffian比值的简洁表达式：

该表达式易于计算，因为它只需要计算一个 $2 \times 2$ 矩阵的Pfaffian，其中 $\mathbf{P} \mathbf{f} \mathbf{C}^{- 1} = - 1$。利用矩阵 $\mathbf{B}$ 的等式(5.118)，可以推导出其显式形式：

由于在模拟开始时已经计算了矩阵 $\tilde{\mathbf{F}}^{-1}$，此计算 $\tilde{\mathbf{F}}^{-1}$ 需要 $O(N_e)$ 次操作。与之前一样，一旦新的位型（configuration）$| x^{\prime} \rangle$ 在马尔可夫链中被接受，就必须更新矩阵 $\tilde{\mathbf{F}}^{-1}$。根据等式(5.117)，我们有：

其中

然后，可以遵循与第5.6节相同的策略方便地求出逆矩阵 $\tilde{\mathbf{F}}^{-1}$：

与第5.6节情况的重要区别在于，不能仅用矩阵 $\mathbf{G}$ 来写出更新式，还必须保留 $\tilde{\mathbf{F}}^{-1}$。事实上，将等式(5.125)两边乘以 $F_{R , R_{i}}$ 并对 $i$ 求和，我们得到：

其中我们使用了如下事实：

然后，利用等式(5.125)将 $\tilde{F}_{l , j}^{- 1 \prime}$ 表达为 $\tilde{F}_{l , j}^{- 1}$ 的函数，我们最终得到：

因此，通过考虑矩阵 $\mathbf{G}$ 和 $\tilde{\mathbf{F}}^{-1}$，并在每次提议的新位型被接受时使用等式(5.125)和(5.128)，可以执行快速更新（fast updating）。

在Pfaffian波函数中，除了标准的（等时）格林函数（Green's function）：

还存在反常的（超导）格林函数：

如前所述，式 (5.125) 和 (5.128) 的解释由 Wick 定理（当左矢和右矢不同时）给出，用于更新的格林函数 $G_{R , j}^{\prime}$ 和 $\tilde{F}_{i , j}^{- 1 \prime}$；在这种情况下，反常收缩也会出现，因为 $\left| {\Phi_{0}} \right.$ 没有固定的粒子数。

利用式 (5.105) 推广到多个电子跳跃的情形虽然直接，但非常繁琐。不过，我们可以通过连续应用单电子更新来绕过 m 电子更新。

5.8 能量与关联函数#

变分能量或任何其他可观测量（包括关联函数）都可以利用我们在前面章节中讨论的快速更新（fast updates）技巧轻松计算。实际上，根据式 (5.13)，局部能量或其他局部可观测量可以表示为波函数比值对乘以算子矩阵元 $\langle x | \mathcal{O} | x^{\prime} \rangle$ 的和。让我们考虑一个一般的 n 体算子：

其中 $d_{I}^{\dagger} \left( d_{I} \right)$ 在位点 $I = 1 , \ldots , 2 L$ 上产生（湮灭）一个电子（对于玻色子系统也有类似定义）。于是，式 (5.5) 的局部估计量由下式给出：

该方程右端的比值可以通过利用快速更新算法章节中讨论的结果轻松计算。

5.9 实际实现#

在此，我们概述变分蒙特卡洛算法实际实现中的重要步骤。

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

_• 生成基组 $| x \rangle$ 的一个随机态，该态存储于向量 iconf(L) 中，其元素给出了格点 $i = 1 , \ldots , L$ 上的局域组态。例如，对于费米子（单带）哈伯德模型，iconf(i) 可以取值为 0（空位点）、±1（一个自旋向上或向下的电子）或 2（双占据位点）。为了更改矩阵 $\tilde{\mathbf{U}}$ 的正确行（对于行列式）或矩阵 $\tilde{\mathbf{F}}$ 的正确行和列（对于 Pfaffian），还需存储一个向量 kel(2L)，其非零元素对于每个格点，给出了定义采样组态 $| x \rangle$（式 (5.44)）的字符串中产生算符的位置。对于玻色子哈伯德模型，iconf(i) 可以取所有非负整数值，且无需存储 kel(2L)。

_• 验证初始组态是非奇异的，即 $\langle x \vert \Psi_{\mathrm{J}} \rangle ~ \ne ~ 0$ 。在这方面，Jastrow 因子通常不会引起任何问题，我们只需检查 $\langle x | \Phi_{0} \rangle$ 是否为零即可。实际上，在处理费米子态时，对于某些 |x>，可能会出现 $\langle x | \Phi_{0} \rangle = 0$ 的情况。虽然这在模拟过程中不是问题（因为如果这些奇异组态被提议，Metropolis 算法不会接受它们），但用奇异组态初始化马尔可夫链将是个问题（因为 Metropolis 算法中的首次接受概率分母将为零）。注意，重叠 $\langle x | \Phi_{0} \rangle$ 的典型值随系统规模的增大呈指数级减小；然而，在快速更新算法中，我们永远不需要计算行列式或 Pfaffian 的实际值，而只需计算它们的比值。因此，在模拟开始时，我们只需确保矩阵 $\tilde{\mathbf{U}}$ 在无数值舍入的情况下是可逆的，以便构造等时格林函数。

_• 计算式 (5.33) 的表格，以执行 Jastrow 因子的快速更新，以及所有进行快速更新所必需的格林函数。对于行列式情况，只需式 (5.80) 的格林函数；而对于 Pfaffian 情况，则需要标准格林函数 (5.129) 以及反常收缩格林函数 (5.130) 和 (5.131)。

2. 使用 Metropolis 算法的马尔可夫过程。

_• 通过移动一个或少数粒子来提议一个新的（随机）组态 $| x^{\prime} \rangle$ 。

_• 计算新波函数与旧波函数之间的比值，参见式 (5.78)、(5.88) 或 (5.121)。

_• 根据 Metropolis 算法（式 (5.12)）接受或拒绝提议的组态。

_• 如果新组态被接受，则更新 Jastrow 因子的表格 (5.34) 以及所有格林函数，即对于行列式为式 (5.79)，对于 Pfaffian 为式 (5.125) 和 (5.128)。

3. 可观测量（Observables）的计算#

为了获得不相关的组态，可以每 $O(L)$ 步（一步对应提议移动一个或少数粒子的操作）计算一次可观测量（例如变分能量）。当然，计算观测量的频率取决于接受率，例如，如果接受率较小，则需要在重新计算观测量之前执行更多马尔可夫步，以使后续测量之间去相关。进一步的打包（binning）技术（见第 3.11 节）可用于减少相关性。

通常建议将观测量写入硬盘，然后进行统计分析，因为通常不知道相关时间（correlation time）。实际上，在蒙特卡洛模拟过程中计算观测量需要在模拟开始前即知道这个量。若不巧未能正确估计它，就必须重新开始模拟，而这可能带来巨大的计算开销。相比之下，后处理分析通常代价小得多，因为它只需要简单的操作。