第12章 · 连续空间中的现实模拟#
12.1 引言与动机#
在前几章中,我们介绍了高效的变分蒙特卡罗和投影蒙特卡罗技术,这些技术能够非常精确地描述格点模型的基态性质。将电子坐标限制在给定格点上已被证明非常有用且成功,可以描述许多现象,包括安德森局域化、近藤屏蔽以及莫特绝缘体的存在。在数值计算中,格点具有巨大优势,因为 (i) 构型易于枚举,最重要的是 (ii) 存在天然的短程截断,使我们无需进一步近似即可处理模型。在计算方法中,连续极限代表一种数学抽象,由于所有数值计算的基础是有限精度算术,计算机无法真正捕捉连续极限。因此,我们可以想象连续空间是通过考虑间距为 $a$ 的格点并取极限 $a 0$ 得到的。在这种情况下,电子通过长度为 $a$ 的离散跳跃进行运动,并通过原始的库仑势相互作用。这种方法被称为"格点正则化",其优点是,对于任何给定的格点间距 $a$,系统可以直接在计算机上模拟。然后,通过取足够小的格点间距,我们就能获得连续哈密顿量的正确物理性质。
然而,在变分方法中,每当波函数用少数几个参数表示时(McMillan, 1965; Ceperley et al., 1977),格点正则化会引入不必要的近似,而不会提高算法的效率。尽管如此,我们可能会尝试使用格点正则化来考虑在网格上对波函数进行系统参数化。不幸的是,当需要对电子轨道进行完全优化时,这种方法不适用于现实材料的第一性原理计算;事实上,格点正则化会引入太多参数,目前无法用随机技术处理。例如,考虑一个包含200个原子(对应相同数量的电子和100个分子轨道)的液态氢系统,位于一个 $100 \times 100 \times 100$ 的网格上;这种情况将意味着 $10^{8}$ 个变分参数,任何随机优化方法都无法处理:尽管近年来取得了显著进展,但目前这个数字不能超过 $\approx 10^{4}$ (Neuscamman et al., 2012)。因此,与格点模型一样,变分态可以写成一个非关联部分(如斯莱特行列式)和一个关联因子(如Jastrow因子)的乘积,前者体现了正确的反对称费米子性质,后者包含了电子间的关联。波函数的斯莱特部分通常使用确定性平均场方法获得,其中通过使用非常廉价且高效的自洽方法可以处理具有大量变分参数的问题。然后,斯莱特部分在存在Jastrow因子的情况下不改变,而Jastrow因子可以通过合适的参数化进行优化(Filippi and Umrigar, 1996; Wagner et al., 2009; Needs et al., 2010; Kim et al., 2012)。
在这里,我们希望遵循这样一种思路:由平均场方法给出的斯莱特行列式并非最优解。事实上,只要平均场不能正确描述该系统,这一过程就无法捕捉到正确的低能行为。例如,在这种方法中,我们永远无法描述仅由排斥相互作用驱动的超导现象。相反,即使斯莱特行列式由一个小基组构建,但若在其与Jastrow因子同时存在的情况下进行优化,引入关联(如Jastrow)因子也能实现对基态的定性正确描述。在这种情况下,优化过程可以映射为一个由变分参数张成的缩减子空间中的迭代线性问题。我们注意到,一个参数数量有限的优化波函数对应于某个局域哈密顿量的基态,对于一个好的尝试波函数,该哈密顿量“足够接近”原始哈密顿量。相比之下,变分方法可能本质上存在偏差,因为它无法描述波函数中未被明确考虑的性质。尽管如此,我们认为,即使存在这一局限性,变分方法仍然极其重要,因为目前还不存在一个能解决电子问题的“黑箱”工具。对于给定的系统,我们应该能够构建合理的波函数并获得定量结果。然后,通过与实验进行直接且明确的比较(无需可调参数),可以验证这一方法,因为变分波函数是通过最小化能量期望值(即通过从头算过程)唯一确定的。此外,即使变分方法与实验结果对比不佳,它仍然是有用的,因为这种情况将证伪波函数的形式,从而激发在其中引入新的物理要素。本着这一精神,我们举几个成功的例子,如Anderson(1987)为描述高温超导引入的共振价键态,以及Laughlin(1983)描述分数量子霍尔效应的波函数。
12.2 局域轨道变分波函数#
在这里,我们讨论用于连续情况的关联波函数的一般形式。我们将$N_{e}$个电子的坐标表示为$\bf{X} = \{\mathbf{r}_{1} , \ldots . . . \mathbf{r}_{N_{e}} \}$。在下文中,我们将考虑自旋向上和自旋向下的电子数目相等(即$N_{e} / 2$对自旋相反的电子)。为实用起见,我们可以固定每个电子的自旋,使得前(后)$N_{e} / 2$个粒子具有自旋向上(自旋向下)。在此选择下,自旋向上和自旋向下电子的坐标分别记为$\{\mathbf{r}_{i}^{\uparrow} \}$和$\{\mathbf{r}_{i}^{\downarrow} \}$(其中$i = 1 , \ldots , N_{e} / 2$)。然后,根据泡利原理,变分波函数$\Psi ( \mathbf{x} )$的轨道部分(与此自旋配置相关联)对于所有同自旋粒子的交换必须反对称(而对于交换反自旋粒子则无对称性要求)。多体态的变分尝试波函数写作两个项的乘积:
$$ \boldsymbol{\Psi} ( \mathbf{x} ) = \mathcal{I} ( \mathbf{x} ) F ( \mathbf{x} ) ,\tag{12.1} $$其中,(正的)Jastrow因子$\mathcal{I} ( \mathbf{x} )$在粒子坐标的所有置换下是对称的,并以两体赝势$u ( {\mathbf{r}} , {\mathbf{r}}^{\prime} )$的形式给出:
$$ \mathcal{T} ( \mathbf{x} ) = \exp \left[ \sum_{i < j} u ( \mathbf{r}_{i} , \mathbf{r}_{j} ) \right] ,\tag{12.2} $$而行列式部分$F ( \mathbf{x} )$负责强制执行同自旋粒子间的反对称性,由配对函数$f ( \mathbf{r} , \mathbf{r}^{\prime} )$(定义在配置空间中)给出:
$$ F ( \mathbf{x} ) = \operatorname*{det} f \left( \mathbf{r}_{i}^{\uparrow} , \mathbf{r}_{j}^{\downarrow} \right) .\tag{12.3} $$接下来,我们将限制在单重态波函数的情形下,这意味着配对函数是对称的,即$f ( \mathbf{r} , \mathbf{r}^{\prime} ) = f ( \mathbf{r}^{\prime} , \mathbf{r} )$(在大多数情况下,也取为实数)。
这个行列式由$N_{e} / 2 \times N_{e} / 2$的配对函数矩阵计算得到,称为反对称化双电子幂(AGP)。我们注意到,对于特定的配对函数选择(见下文),非关联部分$F ( \mathbf{x} )$可以简化为一个斯莱特行列式,但它允许我们获得更一般的量子态,例如BCS超导体。此外,重要的是要强调,类似于在格点模型中的情况,Jastrow因子是描述电子间库仑排斥的最简单项,它引入了当两个电子非常接近时对非关联$F ( \mathbf{x} )$的“重整化”。这种效应无法由任何单电子变分Ansatz考虑,因为在这种情况下,每个电子表现为一个独立粒子,受到由所有其他电子平均效应决定的平均场势的影响。通过使用上述表达式,可以实现标准的变分蒙特卡罗方法,从而评估能量和相关函数的期望值,计算时间与$O ( N_{e}^{3} )$成比例(Ceperley等人,1977)。
这里,我们详细描述如何通过使用一组原子轨道来方便地展开两个函数$u ( \mathbf{r} , \mathbf{r} )$和$f ( \mathbf{r} , \mathbf{r}^{\prime} )$(分别出现在Jastrow因子和行列式部分的定义中)。为此,我们考虑一个原子基组$\{\phi_{p} ( \mathbf{r} - \mathbf{R}_{I} ) \}$,其中每个元素都有一个轨道指标$p$,并依赖于一个给定指标I,表示轨道所居中的离子。这里,$p = 1 , \ldots , N_{\mathrm{orb}} ( I )$且$I = 1 , \ldots N_{\mathrm{ion}} , N_{\mathrm{orb}} ( I )$和$N_{\mathrm{ion}}$分别表示与给定离子I相关的轨道总数和离子数。以下,为简化符号,我们使用单一指标$\mu$来表示轨道$p$和原子中心I,从而得到$\{\phi_{\mu} ( \mathbf{r} ) \}$。在这种紧凑符号下,离子位置和轨道指标都是$\mu$的整数函数,即$I = I ( \mu )$和$p = p ( \mu )$。
原子基组不一定是正交归一的。实际上,出于实用目的,选择简单的(例如斯莱特或高斯型,见下文)局域轨道是方便的。在这种方案下,我们通过下式定义重叠矩阵S:
$$ S_{\mu , \nu} = \int \mathbf{dr} \phi_{\mu} ( \mathbf{r} ) \phi_{\nu} ( \mathbf{r} ) .\tag{12.4} $$在实际中,对于开放系统中的高斯型轨道,矩阵S可以精确计算;然而,我们经验证,通过使用有限格点网格并随后数值计算相应的积分,也可以得到非常好的近似。后一种方法非常通用,允许我们考虑斯莱特轨道(甚至更一般的轨道)。我们还强调,在现代超级计算机上,S的数值计算非常高效,通常只需忽略不计的计算量。
然后,配对函数可以一般性地写为:
$$ f ( {\bf r} , {\bf r}^{\prime} ) = \sum_{\mu , \nu} f_{\mu , \nu} \phi_{\mu} ( {\bf r} ) \phi_{\nu} ( {\bf r}^{\prime} ) ,\tag{12.5} $$其中 $f_{\mu , \nu}$ 定义一个实对称矩阵。需要指出,只要原子基组足够大,原则上可以精确表示任意(可归一化的)函数 $f ( \mathbf{r} , \mathbf{r}^{\prime} )$。类似地,关联项 $u ( \mathbf{r}_{i} , \mathbf{r}_{j} )$ 也可以在一组可能不同的原子轨道 $\{\chi_{\mu} ( \mathbf{r} ) \}$ 中展开,并沿用上述紧凑记号。然而,为了加快向完备基组极限的收敛(换言之,在小基组内令人满意地参数化Jastrow项),必须满足所谓尖点条件,这些条件由库仑势在短距离处的发散性所隐含。这样,当 $\mathrm{~ \bf ~ r ~} \to \mathrm{~ \bf ~ r ~}^{\prime}$(电子-电子)或 $\mathrm{~ \bf ~ r ~} \to \mathrm{~ \bf ~ R ~}_{I}$(电子-离子)时,主导的奇异行为得以精确满足,即 $u ( {\bf r , r^{\prime}} ) \approx | {\bf r - r^{\prime}} | / 2$ 和 $u ( {\bf r} , {\bf r^{\prime}} ) \approx - Z_{I} | {\bf r} - {\bf R}_{I} | / ( N_{e} - 1 )$(其中 $Z_{I}$ 表示原子I的电荷数)。因此,关联因子 $u ( {\mathbf{r}} , {\mathbf{r}}^{\prime} )$ 的一般形式写为(这里我们假设 $N_{e} > 1$):
$$ u ( \mathbf{r} , \mathbf{r}^{\prime} ) = \left[ \frac{u_{\mathrm{ei}} ( \mathbf{r} ) + u_{\mathrm{ei}} ( \mathbf{r}^{\prime} )} {N_{e} - 1} \right] + u_{\mathrm{ee}} ( \mathbf{r} , \mathbf{r}^{\prime} ) ,\tag{12.6} $$其中括号内的前两项描述电子-离子相互作用,第三项与真正的电子-电子关联有关:
$$ u_{\mathrm{ei}} ( {\bf r} ) = u_{\mathrm{ei}}^{\mathrm{cusp}} ( {\bf r} ) + \sum_{\nu} w_{\nu} \chi_{\nu} ( {\bf r} ) ,\tag{12.7} $$$$ u_{\mathrm{ee}} ( {\bf r} , {\bf r^{\prime}} ) = u_{\mathrm{ee}}^{\mathrm{cusp}} ( {\bf r} , {\bf r^{\prime}} ) + \sum_{\mu , \nu} u_{\mu , \nu} \chi_{\mu} ( {\bf r} ) \chi_{\nu} ( {\bf r^{\prime}} ) ;\tag{12.8} $$这里,$u_{\mathrm{ei}}^{\mathrm{cusp}} ( \mathbf{r} )$ 和 $u_{\mathrm{ee}}^{\mathrm{cusp}} ( \mathbf{r} , \mathbf{r}^{\prime} )$ 分别被假设满足电子-电子和电子-离子尖点条件。最常用的始终有界的简单表达式为(Foulkes等人,2001):
$$ u_{\mathrm{ei}}^{\mathrm{cusp}} ( \mathbf{r} ) = - \sum_{I} \frac{Z_{I} | \mathbf{r} - \mathbf{R}_{I} |} {1 + \sqrt{2 Z_{I}} b_{\mathrm{ei}} | \mathbf{r} - \mathbf{R}_{I} |} ,\tag{12.9} $$$$ u_{\mathrm{ee}}^{\mathrm{cusp}} ( {\mathbf{r}} , {\mathbf{r}^{\prime}} ) = \frac{| \mathbf{r} - \mathbf{r}^{\prime} |} {2 ( 1 + b_{\mathrm{ee}} | \mathbf{r} - \mathbf{r}^{\prime} | )} ,\tag{12.10} $$其中 $b_{\mathrm{ei}}$ 和 $b_{\mathrm{ee}}$ 是两个独立的(正)变分参数。此外,$u_{\mu , \nu}$ 是一个对称矩阵,$w_{\nu}$ 是一个向量,它们描述了Jastrow因子在给定(有限)原子基组中的变分自由度。将单体贡献乘以 $1 / ( N_{e} - 1 )$ 的原因源于以下恒等式:
$$ \exp \left\{\sum_{i < j} \left[ \frac{u_{\mathrm{ei}} ( \mathbf{r}_{i} ) + u_{\mathrm{ei}} ( \mathbf{r}_{j}^{\prime} )} {N_{e} - 1} \right] \right\} = \exp \left[ \sum_{i} u_{\mathrm{ei}} ( \mathbf{r}_{i} ) \right] ,\tag{12.11} $$因为上式的右边是式(12.2)中单体Jastrow因子常用的表达式。注意,通过选取足够大的基组并采用合适的 $u_{\mu , \nu}$ 和 $w_{\nu}$ 值,可以实现对单体项和双体项的一般性完整描述。
实际上,我们可以对Jastrow因子和配对函数采用相同的基组。然而,通常Jastrow项(包括电子-离子和电子-电子项)的基组维数 $( D_{u}$)可以远小于配对函数所用的基组维数 $( D_{f} )$,因为用较小的 $D_{u}$ 值即可获得各种化学性质的收敛结果。综上所述,我们需要优化的总参数数量为:
$$ N_{\mathrm{par}} = \frac{D_{f} ( D_{f} + 1 ) + D_{u} ( D_{u} + 1 )} {2} + D_{u} .\tag{12.12} $$幸运的是,与量子化学计算中的情况类似,即使存在电子关联,一个非常小的基组也足以相当精确地描述电子波函数。实际上,问题的主要能量尺度是强烈的电子-离子吸引,这对于 $Z_{I}$ 较大的情况构成了主要贡献,并且可以通过仅包含少数元素的适当原子基组非常精确地表示。
12.2.1 原子轨道#
在此,我们讨论原子轨道 $\{\phi_{\mu} ( \mathbf{r} ) \}$ 的形式,这些轨道以离子 $\mathbf{R}_{I}$ 的位 $I$ 置为中心,并通过径向向量 $\mathbf{r} - \mathbf{R}_{I}$ 表示,该向量连接电子位置 r 与离子 $\mathbf{R}_{I}$ 的位置。请注意,在最一般情况下,特定原子 I 上定义有多个轨道(这构成了离子 I 上的原子基组)。基本对象由径向部分和角向部分决定:径向部分由高斯型或斯莱特型函数给出,角向部分则由角动量 l 及其沿 z 轴的投影 m 表征。一方面,高斯型形式为:
$$ \phi_{l , \pm | m | , I}^{G} ( \mathbf{r} ; \zeta ) \propto | \mathbf{r} - \mathbf{R}_{I} | ^{l} e^{- \zeta | \mathbf{r} - \mathbf{R}_{I} | ^{2}} \left[ Y_{l , m , I} ( \pmb{\Omega} ) \pm Y_{l , - m , I} ( \pmb{\Omega} ) \right] ;\tag{12.13} $$另一方面,斯莱特型形式为:
$$ \phi_{l , \pm | m | , I}^{S} ( \mathbf{r} ; \boldsymbol{\zeta} ) \propto | \mathbf{r} - \mathbf{R}_{I} | ^{l} e^{- \zeta | \mathbf{r} - \mathbf{R}_{I} |} \left[ Y_{l , m , I} ( \pmb{\Omega} ) \pm Y_{l , - m , I} ( \pmb{\Omega} ) \right] ;\tag{12.14} $$在这两种情况下,$Y_{l , m , I} ( \Omega )$ 表示围绕 $\mathbf{R}_{I}$ 的球谐函数(取决于单位向量);对于角动量 l(及其沿 z 轴的投影 m)的每个值,斯莱特型和高斯型轨道均取决于(变分)参数 $\zeta$ 。通过将 $Y_{l , m , I} ( \Omega )$ 与其复共轭 $Y_{l , - m , I} ( \pmb{\Omega} )$ (对于 $m \ne 0$)进行线性组合,我们得到一个实的原子基组;$m = 0$ 的情况无需任何对称化,因为 $Y_{l , 0 , I}$ 已经是实的。这些轨道可以用简单的笛卡尔形式表示,这种形式可以高效地实现。
最一般的收缩原子轨道 $\phi_{\mu}^{\mathrm{contr}} ( \mathbf{r} )$ 可以写为高斯型或斯莱特型轨道(记为 $X = G \text{或} S$)的线性组合,每个轨道都与角动量 l 相关联,并由一个额外的索引 k 指定。那么,变分参数 $\zeta$ 可能依赖于 l 和 k;最终收缩轨道由下式给出:
$$ \phi_{\mu}^{\mathrm{contr}} ( {\bf r} ) = \sum_{k , l , m} c_{l , m , I}^{k} \phi_{l , m , I}^{X_{k , l}} ( {\bf r} , \zeta_{k , l} ) .\tag{12.15} $$在量子化学中,这种线性组合被称为收缩,用于降低原子基的维度以描述强局域化的原子轨道(例如1s),尽管混合不同角动量的态并不常见。然而,在量子化学中,缩减原子基组并不关键,主要的计算优势是通过使用少量分子轨道实现的。相比之下,在量子蒙特卡罗方法中,最小化变分参数的数量极为重要,因此必须尽可能降低原子基维度 $D_{f}$。为此,建议在由初等函数构成的大基组中,通过优化形如式(12.15)的若干独立原子轨道来降低 $D_{f}$。这样,只需少量收缩轨道便可达到收敛甚至精确的结果。因此,这类轨道是广义杂化轨道,对于以紧凑高效的方式描述波函数具有根本重要性(Sorella 等,2015)。
与标准电子结构方法中的常见做法一致(Dovesi 等,2014),高斯轨道和斯莱特轨道均可通过考虑 $L_{x} \times L_{y} \times L_{z}$ 超胞推广到周期性体系。这可以方便地通过向量 ${\bf L}_{n_{x} , n_{y} , n_{z}} = ( n_{x} L_{x} , n_{y} L_{y} , n_{z} L_{z} )$(标记所有可能的周期性镜像)来实现,并对一般扭转边界条件(由角度 $| \theta_{x} | , | \theta_{y} |$ 和 $| \theta_{z} | \le \pi$ 确定)求和:
$$ \phi_{l , m , l}^{G_{\mathrm{PBC}} , S_{\mathrm{PBC}}} ( \mathbf{r} ; \zeta ) = \sum_{n_{x} , n_{y} , n_{z}} \phi_{l , m , l}^{G , S} \left( \mathbf{r} + \mathbf{L}_{n_{x} , n_{y} , n_{z}} ; \zeta \right) e^{i ( n_{x} \theta_{x} + n_{y} \theta_{y} + n_{z} \theta_{z} )} .\tag{12.16} $$高斯因子或指数因子的存在保证了上述级数快速收敛。这一新基满足扭转边界条件:
$$ \phi_{l , m , l}^{G_{\mathrm{PBC}} , S_{\mathrm{PBC}}} ( \mathbf{r} + \mathbf{L}_{n_{x} , n_{y} , n_{z}} ; \zeta ) = e^{- i ( n_{x} \theta_{x} + n_{y} \theta_{y} + n_{z} \theta_{z} )} \phi_{l , m , l}^{G_{\mathrm{PBC}} , S_{\mathrm{PBC}}} ( \mathbf{r} ; \zeta ) ,\tag{12.17} $$这对研究材料的体性质非常重要。例如,通过使用适当选择的扭转角(Baldereschi, 1973)或对所有可能角度平均物理可观测量(Lin 等,2001),可以平滑且快速地收敛到热力学极限。
同样的过程可应用于Jastrow因子的轨道,其中 $\phi_{l , \pm | m | , I}^{S , G} ( \mathbf{r} ; \zeta ) \stackrel{-} {\to} \chi_{l , \pm | m | , I}^{S , G} ( \mathbf{r} ; \zeta )$。然而,由于Jastrow因子仅耦合局域密度,我们仅限于 $\theta_{x} = \theta_{y} = \theta_{z} = 0$ 的周期情况,这意味着 $\chi_{l , m , I}^{G_{\mathrm{PBC}} , S_{\mathrm{PBC}}} ( \mathbf{r} + \mathbf{L}_{n_{x} , n_{y} , n_{z}} ; \zeta ) = \chi_{l , m , I}^{G_{\mathrm{PBC}} , S_{\mathrm{PBC}}} ( \mathbf{r} ; \zeta )$。最后,式(12.9)和(12.10)中定义的单体和双体函数 $u_{\mathrm{ei}}^{\mathrm{cusp}} ( \mathbf{r} )$ 和 $u_{\mathrm{ee}}^{\mathrm{cusp}} ( \mathbf{r} , \mathbf{r}^{\prime} )$ 是缓慢衰减的函数,因此用欧几里得距离的周期性推广来表达更为方便:
$$ | \mathbf{r} | \to{\sqrt{\left[ {\frac{L_{x}} {\pi}} \sin \left( {\frac{\pi x} {L_{x}}} \right) \right]^{2} + \left[ {\frac{L_{y}} {\pi}} \sin \left( {\frac{\pi y} {L_{y}}} \right) \right]^{2} + \left[ {\frac{L_{z}} {\pi}} \sin \left( {\frac{\pi z} {L_{z}}} \right) \right]^{2}}} .\tag{12.18} $$12.2.2 分子轨道#
现在,我们考虑另一种不同的泛函形式来表示由配对函数 $f ( \mathbf{r} , \mathbf{r}^{\prime} )$ 所描述的非关联态。通过这种方式,AGP波函数与Slater行列式之间的关系将变得非常透明。(实且对称的)配对函数可被视为连续空间中的一个线性算符;因此,可以通过寻找满足下式的本征函数 $\Phi_{\alpha} ( \mathbf{r} )$ 来对其进行对角化:
$$ \int \mathbf{dr}^{\prime} f ( \mathbf{r} , \mathbf{r}^{\prime} ) \Phi_{\alpha} ( \mathbf{r}^{\prime} ) = \lambda_{\alpha} \Phi_{\alpha} ( \mathbf{r} ) ,\tag{12.19} $$其中 $\lambda_{\alpha}$ 是对应的本征值。在维数为 $D_{f}$ 的有限基组 $\{\phi_{\mu} ( \mathbf{r} ) \}$ 中,本征函数可写为:
$$ \Phi_{\alpha} ( \mathbf{r} ) = \sum_{\mu} P_{\mu , \alpha} \phi_{\mu} ( \mathbf{r} ) .\tag{12.20} $$其中矩阵 P 的元素由广义本征值方程确定,该方程是通过将此表达式以及配对函数展开式 (12.5) 代入方程 (12.19) 得到的。由此可得:
$$ \mathbf{fSP} = \mathbf{P} \pmb{\Lambda} ,\tag{12.21} $$其中 $\pmb{\Lambda}$ 是包含本征值 $\lambda_{\alpha}$ 的对角矩阵,S 是方程 (12.4) 的重叠矩阵,f 是方程 (12.5) 中元素为 $f_{\mu , \nu}$ 的矩阵。注意,作为广义本征值问题的结果,矩阵 P 定义了一组正交归一态:
$$ \int \mathbf{dr} \Phi_{\alpha} ( \mathbf{r} ) \Phi_{\beta} ( \mathbf{r} ) = \delta_{\alpha , \beta} ,\tag{12.22} $$这导致 $\mathbf{P}^{T} \mathbf{S} \mathbf{P} = \mathbf{1}$。单粒子波函数 $\{\Phi_{\alpha} ( \mathbf{r} ) \}$ 通常被称为分子轨道。注意,根据矩阵 $\mathbf{P}_{\mathrm{:}}$ 的性质,分子轨道可以是局域的,也可以是扩展的。
最后,由方程 (12.19),配对函数可以用分子轨道表示:
$$ f ( \mathbf{r} , \mathbf{r}^{\prime} ) = \sum_{\alpha = 1}^{D_{f}} \lambda_{\alpha} \Phi_{\alpha} ( \mathbf{r} ) \Phi_{\alpha} ( \mathbf{r}^{\prime} ) .\tag{12.23} $$值得注意的是,当非零本征值 $\lambda_{\alpha}$ 的个数等于电子对数 $N_{e} / 2$ 时,方程 (12.3) 中的行列式 $F ( \mathbf{x} )$ 可以写成三个维数为 $N_{e} / 2 \times N_{e} / 2$ 的方阵的乘积:
$$ \operatorname*{det} f ( \mathbf{r}_{i}^{\uparrow} , \mathbf{r}_{j}^{\downarrow} ) = \operatorname*{det} ( \mathbf{P}_{\uparrow} \tilde{\Lambda} \mathbf{P}_{\downarrow}^{T} ) = \left( \prod_{\alpha = 1}^{N_{e} / 2} \lambda_{\alpha} \right) \times \operatorname*{det} \mathbf{P}_{\uparrow} \times \operatorname*{det} \mathbf{P}_{\downarrow}\tag{12.24} $$其中 $\tilde{\Lambda}$ 是仅包含非零本征值的(对角)方阵;矩阵 $\mathbf{P}_{\uparrow}$ 和 $\mathbf{P}_{\downarrow}$ 由下式给出:
$$ ( P_{\uparrow} )_{i , \alpha} = \Phi_{\alpha} ( \mathbf{r}_{i}^{\uparrow} ) ,\tag{12.25} $$$$ ( P_{\downarrow} )_{i , \alpha} = \Phi_{\alpha} ( \mathbf{r}_{i}^{\downarrow} ) .\tag{12.26} $$因此,在这一特殊情况下,我们得到了与斯莱特行列式标准表达式一致的形式(无关常数,即非零特征值的乘积除外),即 AGP 波函数简化为一个具有 $N_{e} / 2$ 个分子轨道的斯莱特行列式。换言之,在 AGP 形式体系中定义的分子轨道,与广泛用于定义斯莱特行列式的分子轨道相一致。
在数值计算角度,将简单行列式表示为具有受限分子轨道数的 AGP 态有重要优势。最简单的一点是,只需考虑一个 $N_{e} / 2 \times N_{e} / 2$ 矩阵的行列式,即可通过式 (12.3) 计算波函数;相比之下,在斯莱特行列式的标准表示中,我们需要计算两个相同维度的行列式。另一个优势在于,可以简单地利用对称性,例如晶体中的平移不变性。此时,只需要求对应于相同类型轨道的配对矩阵 $\mathbf{f}$ 仅依赖于离子位置之差,即 $f_{\mu , \nu} = f_{\mu^{\prime} , \nu^{\prime}} \mathrm{if} | \mathbf{R}_{I ( \mu^{\prime} )} - \mathbf{R}_{I ( \nu^{\prime} )} | = | \mathbf{R}_{I ( \mu )} - \mathbf{R}_{I ( \nu )} |$,这表示一个非常简单的线性约束。
最后,我们还注意到,即使 AGP 态与斯莱特行列式完全等价,对 Jastrow 因子和分子轨道的联合优化也可能产生性质上不同的波函数,即其化学和物理性质不同于无关联的 Hartree-Fock 或密度泛函理论所得到的结果。具有 $N_{e} / 2$ 个分子轨道的完全优化 Jastrow-斯莱特波函数,能够精确描述原子(Foulkes 等人,2001),可恢复约 90% 的相关能(定义为通过实验或高质量计算得到的精确估计值与最佳 Hartree-Fock 值之间的能量差)。通过释放对分子轨道的约束,我们获得了更大的变分自由度;显然,无约束的 AGP 波函数可以改进 Hartree-Fock 方法,尤其是在同时包含 Jastrow 因子的情况下,能够提高化学键描述的精度,甚至产生新的定性性质。特别地,AGP 态允许电子间显式配对,因此可以描述超导体;而在存在 Jastrow 因子的情况下,还可以描述 Mott(共振价键)绝缘体,正如安德森(1987)在高温超导体发现后不久所提出的。
12.3 变分波函数的尺寸一致性#
量子化学中的一个基本概念是所谓的尺寸一致性(size consistency),即当我们计算一个可明确区分出两个空间上分离的区域A和B的分子系统的能量时,组合系统 $E_{A + B}$ 的总能量应等于孤立系统 $E_{A}$ 和 $E_{B}$ 的能量之和,即 $E_{A + B} = E_{A} + E_{B}$。尺寸一致性对化学键有至关重要的影响,因此被认为是完全变分优化波函数最重要的性质之一。众所周知,一般来说,AGP 波函数不具有尺寸一致性,即便简单的 Hartree-Fock 态也不受此问题影响。因此,AGP 态在化学中已被弃用。在 AGP 方法中,尺寸不一致性源自对给定粒子数的投影。例如,考虑 $\mathrm{Be}_{2}$ 分子,其中区域 A 和 B 各包含一个 Be 原子;那么,单个原子的 AGP 波函数正确地只包含两对电子。然而,在分子中,电子数的约束只作用于总数为四对的情况。因此,发生成键时,$A + B$ 系统的 AGP 波函数也允许这样的情况:靠近一个原子的区域有三对甚至四对电子,而靠近另一个原子的区域只有一对或零对电子。当分子被拉伸,两个原子处于极大距离时,这种效应不会产生两个独立的原子 AGP 波函数,因为投影总电子数并不等同于投影每个原子内的电子数。因此,相应的能量将远高于两个独立 Be 原子的能量,尺寸一致性因此得不到满足。
在此,我们证明,当 AGP 波函数与 Jastrow 因子结合时,尺寸一致性得以恢复,因为后者具有变分自由度,允许对不同空间区域进行部分数投影(Sorella 等人,2007;Neuscamman,2012)。为此,我们可以考虑一个 Jastrow 因子,它将所有在区域 A 和 B 中分别不恰好有 $N_{A}$ 和 $N_{B}$ 个粒子的电子构型投影掉:
$$ \mathcal{I} ( \mathbf{x} ) \propto \exp \{- \frac{C} {2} [ N_{A} ( \mathbf{x} ) - N_{A} ]^{2} + A B \} ,\tag{12.27} $$其中 $\begin{array} {r} {N_{A , B} ( {\bf x} ) = \sum_{i} g_{A , B} ( {\bf r}_{i} )} \end{array}$,且 $g_{A} ( \mathbf{r} ) = 1[ g_{B} ( \mathbf{r} ) = 1]$ 若 $\mathbf{r} \in A( \mathbf{r} \in B)$,否则为零,C 为足够大的常数。这种形式的 Jastrow 因子可以通过取如下特定形式的赝势得到:
$$ u ( {\bf r , r^{\prime}} ) = - C \{g_{A} ( {\bf r} ) g_{A} ( {\bf r^{\prime}} ) - \frac{N_{A} - 1 / 2} {N_{e} - 1} [ g_{A} ( {\bf r} ) + g_{A} ( {\bf r^{\prime}} ) ] + A B \} .\tag{12.28} $$因此,只要Jastrow因子的基组足够大,能够表示式(12.28)的形式,对全局波函数的完整优化就能保证总能量低于两个碎片(当区域A和B充分分离时)的能量,从而保证了大小一致性。上述规则的一个例外情况是两个碎片具有非零自旋:此时,由单个Slater行列式构建的波函数往往不具有大小一致性(例如两个氧原子的情况,每个独立原子自旋为$S = 1$,而$\mathrm{O}_{2}$分子自旋为$S = 0$)。因此,只要整个系统不具有与两个碎片相容的最大自旋,限制性Hartree-Fock行列式就不具有大小一致性。
12.4 变分波函数的优化#
优化配对函数的一个直接方法是将式(12.5)中的所有矩阵元$f_{\mu , \nu}$视为独立的变分参数。这将是一种无约束优化。然而,这种方法并不总是高效的,因为变分能量对长程配对项的依赖性可能非常弱,尤其是在绝缘相中。将配对函数显式地定义在局域基组上的主要优势在于可以利用关联的局域性。当geminal函数$f ( \mathbf{r} , \mathbf{r}^{\prime} )$由特征长度$\xi$描述时,减少变分参数数量的一种有效方法是,在优化中只考虑连接距离小于某个截断半径$R_{\mathrm{cut}}$的轨道的矩阵元$f_{\mu , \nu}$。对于给定的$R_{\mathrm{cut}}$,矩阵$f_{\mu , \nu}$的每一列(或每一行)中只有固定数量的非零元素,因此参数总数与总原子基组维度成线性关系,导致变分参数数量$N_{\mathrm{par}} \propto D_{f} \propto N_{e}$。与对所有矩阵元$f_{\mu , \nu}$进行独立优化的无约束最小化相比,这种方法高效得多。当物理问题中存在有限关联长度时,通过增大$R_{\mathrm{cut}}$,精度应呈指数级提高。例如,如果geminal矩阵的秩等于电子数的一半$N_{e} / 2$,则geminal函数对应于密度矩阵$\rho ( {\bf r} , {\bf r^{\prime}} )$。那么,在绝缘相中,我们有:
$$ \operatorname*{lim}_{| \mathbf{r} - \mathbf{r}^{\prime} | \to \infty} \rho ( \mathbf{r} , \mathbf{r}^{\prime} ) \propto \exp \left( - \frac{| \mathbf{r} - \mathbf{r}^{\prime} |} {\xi} \right) ,\tag{12.29} $$这证明了引入$R_{\mathrm{cut}}$以高效减少参数数量的合理性。不幸的是,当考虑金属相时,这种方法变得非常低效,因为在三维空间中:
$$ \operatorname*{lim}_{| \mathbf{r} - \mathbf{r}^{\prime} | \to \infty} | \rho ( \mathbf{r} , \mathbf{r}^{\prime} ) | \propto \frac{1} {| \mathbf{r} - \mathbf{r}^{\prime} | ^{2}} ,\tag{12.30} $$这意味着随着$R_{\mathrm{cut}}$的增大,收敛速度非常缓慢。密度矩阵缓慢衰减的原因是费米面的存在,它限制了部分填充能带内的占据数。
在此,我们定义一种约束优化程序,它能够同时处理绝缘体和金属行为,即使在存在有限截断$R_{\mathrm{cut}}$的情况下也是如此。这一目标可以通过将配对函数的秩固定为给定数n来实现,即只在geminal函数的定义中取n个分子轨道:
$$ f^{c} ( \mathbf{r} , \mathbf{r}^{\prime} ) = \sum_{\alpha = 1}^{n} \lambda_{\alpha} \Phi_{\alpha} ( \mathbf{r} ) \Phi_{\alpha} ( \mathbf{r}^{\prime} ) ,\tag{12.31} $$其中分子轨道 $\{\Phi_{\alpha} ( \mathbf{r} ) \}$ ($\alpha = 1 , \ldots , n$)定义了“占据态”。在这一过程中,通过施加上述 n 个占据态的约束,可以“重构”费米面,从而允许我们恢复缓慢衰减的密度矩阵。通过将配对函数的优化约束为有限数量 n 个分子轨道的形式,我们获得了三个主要优势。第一个是配对函数的双重表示,这使我们能够将矩阵元 $f_{\mu , \nu}$ 视为变分参数,同时满足固定秩 n(即固定数量的 n 个优化分子轨道)的约束;第二个优势是,在优化过程中,我们可以通过取来自同一局域基组中定义的平均场方法(如 Hartree-Fock 或密度泛函理论)的分子轨道 $\Phi_{\alpha} ( \mathbf{r} )$ 来初始化配对函数(即矩阵 $f_{\mu , \nu}$);第三个优势是,我们可以在受约束的方法中实现局域性,从而极大减少变分参数的数量,而不会对计算精度造成太大影响,因为配对函数的长程部分通过限制分子轨道数量(即强制存在费米面)被自动考虑了。
在讨论技术细节之前,我们想提一下,在晶格上我们使用了不同的方法来最小化变分参数(见第5章和第6章)。在那里,我们考虑了一种由具有有限作用范围的等效(BCS)哈密顿量给出的参数化(即,仅对应于距离小于 $R_{\mathrm{cut}}$ 的参数);该等效哈密顿量的基态隐含地定义了配对函数,在金属或无能隙超导体的情况下,该配对函数是长程的。在连续系统中,实现基于辅助哈密顿量的方法比较困难,而直接处理配对函数要容易得多。
对于固定秩 n 的受约束配对函数,通过进行无穷小变化并保持相同秩,我们得到:
$$ \delta f^{c} ( {\bf r} , {\bf r}^{\prime} ) = \sum_{\alpha = 1}^{n} \delta \lambda_{\alpha} \Phi_{\alpha} ( {\bf r} ) \Phi_{\alpha} ( {\bf r}^{\prime} ) + \lambda_{\alpha} \left[ \delta \Phi_{\alpha} ( {\bf r} ) \Phi_{\alpha} ( {\bf r}^{\prime} ) + \Phi_{\alpha} ( {\bf r} ) \delta \Phi_{\alpha} ( {\bf r}^{\prime} ) \right] .\tag{12.32} $$然后,我们可以很容易地验证这种修改满足以下条件(通过使用原始分子轨道的正交归一化得到):
$$ \int \mathrm{d} \mathbf{r}_{1} \int \mathrm{d} \mathbf{r}_{2} \left[ \delta ( \mathbf{r} - \mathbf{r}_{1} ) - \Pi ( \mathbf{r} , \mathbf{r}_{1} ) \right] \delta f^{c} ( \mathbf{r}_{1} , \mathbf{r}_{2} ) \left[ \delta ( \mathbf{r}^{\prime} - \mathbf{r}_{2} ) - \Pi ( \mathbf{r}_{2} , \mathbf{r}^{\prime} ) \right] = 0 ,\tag{12.33} $$其中投影到占据态上的投影算子 $\Pi ( \mathbf{r} , \mathbf{r}^{\prime} )$ 由下式给出:
$$ \Pi ( {\bf r} , {\bf r}^{\prime} ) = \sum_{\alpha = 1}^{n} \Phi_{\alpha} ( {\bf r} ) \Phi_{\alpha} ( {\bf r}^{\prime} ) .\tag{12.34} $$在矩阵表示中,对应于式 (12.5) 中的有限基组 Ansatz,式 (12.33) 变为:
$$ ( \mathbf{1} - \mathbf{L} ) \delta \mathbf{f}^{c} ( \mathbf{1} - \mathbf{R} ) = 0\tag{12.35} $$其中 L 和 R 是两个 $D_{f} \times D_{f}$ 矩阵,由于原始基组 $\{\phi_{\mu} ( \mathbf{r} ) \}$ 是非正交的,它们由下式给出:
$$ \mathbf{R} = ( \mathbf{S} \mathbf{P} ) \mathbf{P}^{T} ,\tag{12.36} $$$$ \mathbf{L} = ( \mathbf{R} )^{T} = \mathbf{P} ( \mathbf{S} \mathbf{P} )^{T} ;\tag{12.37} $$这里,$D_{f} \times n$ 矩阵 P 由式(12.20)定义,而 $D_{f} \times D_{f}$ 矩阵 S 是定义基组的重叠矩阵,参见式(12.4)。在实际操作中,由于维度 n 通常远小于 $D_{f}$,仅使用两个矩阵 P 和 SP 会更方便;实际上,所有涉及 L 和 R 的线性运算在利用矩形矩阵 P 和 SP 进行时都会快得多。
至此,我们能够以一种简单的方式实现约束。假设我们对配对函数进行无约束变化 δf,那么,通过考虑下式,我们可以从该变化中移除所有(在线性阶上)影响矩阵秩的贡献:
$$ \delta \mathbf{f}^{c} = \delta \mathbf{f} - ( \mathbf{1} - \mathbf{L} ) \delta \mathbf{f} ( \mathbf{1} - \mathbf{R} ) .\tag{12.38} $$该方程表示无约束变化与约束变化之间的线性映射。需要强调的是,这个映射可以看作是对应能函数变化的一种简单参数化,它被约束为固定数量的 n 个分子轨道。此外,即使通过引入 $R_{\mathrm{{cut}}}$ 将 f 限制在有限范围内,应用式(12.38)也可能在对应能函数中诱导出长程尾部。这种方法有助于表示无隙相,如金属和节超导体,并凸显了即使在金属行为的情况下,$R_{\mathrm{{cut}}}$ 收敛速度也可能非常快的原因。
在应用第6章描述的优化技术时,我们必须计算波函数对对应能矩阵 f 的任意变化的导数。我们用如下公式表示无约束导数:
$$ D_{\mu , \nu} = \frac{\partial \ln | F ( \mathbf{x} ) |} {\partial f_{\mu , \nu}} .\tag{12.39} $$相对于约束参数的总变化由 $\mathrm{Tr} [ \mathbf{D} ( \delta \mathbf{f}^{c} )^{T} ]$ 给出,其中 $\delta \mathbf{f}^{c}$ 由式(12.38)给出。因此,我们可以定义约束导数,即相对于受约束参数的变化:
$$ D_{\mu , \nu}^{c} = \left[ \mathbf{D} - \left( \mathbf{1} - \mathbf{L}^{T} \right) \mathbf{D} \left( \mathbf{1} - \mathbf{R}^{T} \right) \right]_{\mu , \nu} ,\tag{12.40} $$这可以通过标准且高效的线性代数运算获得。
总结来说,配对函数的约束优化方案如下:
计算式(12.39)的导数矩阵 D,不考虑任何约束。
利用当前分子轨道应用式(12.40)的投影。通过利用所有涉及的矩阵都可以用小型矩形矩阵表示这一事实,可以实现复杂度为 $O ( n D_{f}^{2} )$ 的非常方便的计算。
在积累足够多的统计信息后,应用最速下降法(或更复杂的方案,如第6章所述)得到 $\delta f_{\mu , \nu}$ ∝ $D_{\mu , \nu}^{c}$。
使用式(12.38)获取约束后的 δfc。
根据下式更新参数:
只考虑截断半径 $R_{\mathrm{cut}}$ 允许的那些轨道。
最后,对角化矩阵 f,仅保留绝对值最大的特征值对应的 n 个特征向量。在此对角化之后,定义新的分子轨道,并重新计算相应的投影矩阵 L 和 R(如果在优化收缩轨道参数时也重新计算了重叠矩阵 S,则一并处理)。
重复上述所有步骤,直至收敛。
上述算法非常稳定且具有通用性,因为它即使在准粒子能量在费米能量处简并(这是金属行为的典型特征)时也能工作。实际上,在该方法中,特征值与准粒子能量无关,仅代表分子轨道的平均占据数。在表 12.1 中,我们报告了在 $r_{s} = 1.32$ 的体心立方晶体中 128 个氢原子(即 128 个电子)的情况。
为了实现对最小值的稳定收敛,在每一步定义新矩阵 f 时,所考虑的非零特征值必须显著非零,因为如果它们过于接近零,则在数值精度内确定矩阵的秩可能会出现问题。特别地,当 $n = N_{e} / 2$ 时,我们可以在优化的每次迭代中将最大的 n 个特征值重置为 $\lambda_{\alpha} = 1$(因为如前所述,对于 Slater 行列式,特征值 $\lambda_{\alpha}$ 仅定义一个无关紧要的归一化常数),从而继续进行非常稳定的优化。
12.4.1 铍分子和氟分子#
这里,我们展示了前面章节描述的优化方案在 $\mathrm{F}_{2}$ 和 $\mathrm{Be}_{2}$ 分子上的两个示例(Marchi 等人,2009)。让我们从前者开始,它相对容易处理。在此,行列式由每个离子上的 $5s5p2d$ 高斯轨道基组构建而成(导致 $D_{f} = 60$)。Hartree-Fock 方法(使用 $n = 9$ 个分子轨道),辅以 Jastrow 因子,在结合能上存在相当大的误差($E_{\text{binding}} \approx 39 \, \text{mH}$,而估计的精确值约为 62 mH),但平衡距离近似得很好($R_{\text{eq}} \approx 2.67 \, \text{a.u.}$,这与估计的精确值一致),见图 12.1。整体能量随两个氟原子间相对距离 R 的变化在 $R \approx 4.5 \, \text{a.u.}$ 处出现一个非物理的最大值。此外,在大距离下,这种描述总单态波函数无法恢复两个 $S = 1/2$ 氟原子的尺寸一致结果。
表 12.1. 在 $r_{s} = 1.32$ 的体心立方晶体中,128 个氢原子在不同 $R_{\mathrm{cut}}$ 下的优化能量(a.u.)。使用了扭曲边界条件,其中 $\theta_{x} = \theta_{y} = 2\pi/6$ 且 $\theta_{z} = \pi$,参见式 (12.17)。原子基组通过收缩式 (12.15) 中的两个杂化轨道获得,这些轨道在 3s1p 基组中展开,包含三个 s 波(具有不同指数 $\zeta_{k, l}$)和一个 p 波高斯轨道。为了体现小 $R_{\mathrm{cut}}$ 时获得的显著精度,我们还报告了相对于 $R_{\mathrm{cut}} = 0$ 情况的关联能差异(使用重加权技术计算)。同时报告了使用密度泛函理论(DFT)计算得到的 Slater 行列式,并配以优化后的 Jastrow 因子(每个原子有两个 s 波高斯轨道)的情况。注意到,即使 $R_{\mathrm{cut}} = 0$ 的情况,相比后者也已有约 90% 的改进。
| $R_{\mathrm{cut}}$ | AGP 参数数量 | 总能 | 关联能差异 |
| DFT | 0 | -64.945(2) | 0.064(2) |
| 0 | 11 | -65.009(1) | 0.0 |
| 2.5 | 43 | –65.010(2) | -0.0018(1) |
| 3 | 67 | –65.017(2) | -0.0064(1) |
| 4 | 115 | –65.014(1) | −0.0072(2) |
| 5 | 243 | –65.015(1) | −0.0072(2) |
| ∞ | 519 | –65.015(2) | -0.0071(2) |
通过增加分子轨道数 $n$ 可以显著改进结果。事实上,当 $n = 10$ 时,加入 Jastrow 因子后的 AGP 波函数在结合能和平衡距离上都给出了非常精确的结果(即 $E_{\mathrm{结合能}} \approx 61\ \mathrm{mH}$),参见图12.1。此外,在这种情况下还恢复了大小一致性。进一步增加分子轨道数,我们并没有显著改善 $n = 10$ 的结果。我们强调这样一个事实:通过使用单个行列式,可以达到与使用数千个行列式的量子化学方法相当(甚至更好)的精度。
双原子分子中有一个有趣的例子——铍二聚体,它对大多数数值技术来说都是一个真正的挑战。


图12.1 $\mathrm{F}_{2}$(上)和 Be2(下)分子结合能的结果。添加了 Jastrow 因子的 AGP(Hartree-Fock)行列式用实心(空心)圆表示。估计的精确结果由实线表示。$\mathrm{F}_{2} \ ( \mathrm{Be}_{2} )$ 分子的 Hartree-Fock 状态具有 $n = 9 ( n = 4 )$ 个分子轨道,这些轨道由每个离子上的高斯轨道基组 $5s5p2d$(高斯型和斯莱特型混合的 $5s2p1d$ 基组)构成。$\Gamma_{2} ~ ( \mathrm{Be}_{2} )$ 分子的 AGP 状态具有 $n = 10 ( n = 32 )$,并使用与 Hartree-Fock 波函数相同的原子轨道。参数 $\zeta_{k , l}$ 经过优化以提高有限基组的精度。
这里的困难来自较小的结合能,其中弱范德瓦尔斯力也起着重要作用。事实上,简单的 Hartree-Fock 方法无法束缚该分子,而其他基于密度泛函理论的自洽方法则严重高估了结合能。我们考虑一个行列式,每个离子上使用了高斯型和斯莱特型混合的 $5s2p1d$ 基组(导致 $D_{f} = 32$)。值得注意的是,通过考虑 $n = 4$ 个分子轨道(即 AGP 函数对应于 Hartree-Fock 状态的情况),我们在结合能(例如,$E_{\mathrm{结合能}} \approx 4.1\ \mathrm{mH}$,可与估计的精确值约 4.2 mH 相比)和平衡位置($R_{\mathrm{eq}} \approx 4.75\ \mathrm{a.u.}$)上都得到了出色的结果,参见图12.1。然而,当包含更多分子轨道时,总能量的精度有所提高,但整体的能量分散性变差,导致键长更大而结合能更小。这种结果是因为在这种情况下,关联的 AGP 状态对于单个原子(它达到了相关能的 99.9%)比对化学键的描述要有效得多。
12.5 晶格正则化扩散蒙特卡洛#
在本节中,我们将考虑如下形式的电子哈密顿量(以原子单位表示):
$$ \mathcal{H} = - \frac{1} {2} \nabla^{2} + V ( \mathbf{x} ) + V_{\mathrm{ii}} ,\tag{12.42} $$其中 $\begin{array} {r} {\nabla^{2} = \sum_{i = 1}^{N_{e}} \nabla_{i}^{2}} \end{array}$ 是定义总动能 $\nabla_{i}^{2}$ 的拉普拉斯算子,$\nabla_{i}^{2}$ 是作用于单个电子(坐标为 $\mathbf{r}_{i}$)的相应拉普拉斯算子,而 $V ( \mathbf{x} ) = V_{\mathrm{ee}} ( \mathbf{x} ) + V_{\mathrm{ei}} ( \mathbf{x} )$ 是标准的多体势,包括电子-电子相互作用:
$$ V_{\mathrm{ee}} ( \mathbf{x} ) = \sum_{i < j} {\frac{1} {| \mathbf{r}_{i} - \mathbf{r}_{j} |}} ,\tag{12.43} $$以及电子-离子相互作用:
$$ V_{\mathrm{ei}} ( {\bf x} ) = - \sum_{i} \nu_{\mathrm{ei}} ( {\bf r}_{i} ) ,\tag{12.44} $$其中
$$ \nu_{\mathrm{ei}} ( \mathbf{r}_{i} ) = \sum_{I} {\frac{Z_{I}} {| \mathbf{r}_{i} - \mathbf{R}_{I} |}}\tag{12.45} $$表示其对位于 $\mathbf{r}_{i}$ 的电子的贡献;最后,与电子坐标无关的离子-离子项 $V_{\mathrm{ii}}$ 由下式给出:
$$ V_{\mathrm{ii}} = \sum_{I < J} {\frac{Z_{I} Z_{J}} {| \mathbf{R}_{I} - \mathbf{R}_{J} |}} ;\tag{12.46} $$这里,离子坐标 $\mathbf{R}_{I}$($I = 1 , \ldots , N_{\mathrm{ion}}$)在玻恩-奥本海默近似下只是给定的参数,因此 $V_{\mathrm{ii}}$ 表示电子能量的一个常数偏移。
在连续模型中,能量谱具有任意大的正本征值(例如,具有大动量的平面波);因此,为了从给定试探函数 $| \Psi_{0} \rangle$ 中过滤出基态 $| \Upsilon_{0} \rangle$,无法应用第1.7节的幂法,而引入指数算符更为方便:
$$ \operatorname*{lim}_{\tau \to \infty} e^{- \mathcal{H} \tau} | \Psi_{0} \rangle \propto | \Upsilon_{0} \rangle .\tag{12.47} $$然而,对于大的 $\tau$ 值,很难精确计算这个投影算符。传统方法是处理短时间 $\Delta \tau$ 的极限,其中Trotter分解产生一个随 $\Delta \tau \to 0$ 消失的小误差(Trotter, 1959; Suzuki, 1976a,b),并多次应用近似形式。在这种方法中,存在一些通常没有严格证明的近似,特别是对于费米子。实际上,扩散蒙特卡洛方法及其相关的固定节点近似(Kalos et al., 1974; Ceperley and Alder, 1980; Foulkes et al., 2001)的少数方面仍有争议(例如,误差如何随 $\Delta \tau \to 0$ 消失,或者在使用所谓的赝势去除核心电子时,能否严格定义变分上界性质)。
在此,我们采用一种依赖晶格正则化的替代观点。我们从拉普拉斯算符 $\nabla^{2}$(定义每个电子的动能)的近似开始,使用格点间距为 a 的有限网格:
$$ \begin{array} {c} {{\nabla_{a}^{2} g ( x , y , z ) = \displaystyle \frac{1} {a^{2}} \left[ g ( x + a , y , z ) + g ( x - a , y , z ) - 2 g ( x , y , z ) \right]}} \\ {{{}}} \\ {{+ \displaystyle \frac{1} {a^{2}} \left[ g ( x , y + a , z ) + g ( x , y - a , z ) - 2 g ( x , y , z ) \right]}} \\ {{{}}} \\ {{{} + \displaystyle \frac{1} {a^{2}} \left[ g ( x , y , z + a ) + g ( x , y , z - a ) - 2 g ( x , y , z ) \right] ,}} \end{array}\tag{12.48} $$其中 $g ( x , y , z )$ 是定义在三维空间 $( x , y , z )$ 中的一般函数。对于 $a \to 0$,线性算符 $\nabla_{a}^{2}$ 收敛到精确拉普拉斯算符,对于足够规则函数 g,误差按 $a^{2}$ 量级缩放。这一事实可以通过对函数 $g ( x \pm a , y , z )$(以及 y 和 z 分量类似)进行泰勒展开,并注意到所有奇数阶导数消失,留下一个 $O ( a^{2} )$ 的残余项轻松得到。
在下文中,根据式(12.48),我们可以将 $\nabla_{i , a}^{2}$ 定义为作用在单个电子坐标 $\mathbf{r}_{i}$ 上的拉普拉斯算符。哈密顿量 $\mathcal{H}$ 的晶格正则化对应于定义一个厄米算符 ${\mathcal{H}}_{a}$,使得当 $a \to 0$ 时 $\mathcal{H}_{a} \to \mathcal{H}$;然后,可以通过稳定模拟数值计算 ${\mathcal{H}}_{a}$ 的性质,并恢复 $a \to 0$ 时连续极限的性质。上述要求绝非显而易见;事实上,$V_{\mathrm{ei}} ( \mathbf{x} )$ 是一个无限吸引势,可能导致出现无下界的负能量。这一事实引导我们,当电子坐标过于接近离子位置时,定义一种电子-离子相互作用的正则化:
$$ \nu_{\mathrm{ei}}^{a} ( \mathbf{r}_{i} ) = \sum_{I} \frac{Z_{I}} {\mathrm{Max} ( | \mathbf{r}_{i} - \mathbf{R}_{I} | , a )} ,\tag{12.49} $$这显然在短距离处截断了库仑奇点。因此,电子-离子势被修改为:
$$ V_{\mathrm{ei}}^{a} ( {\bf x} ) = - \sum_{i} \nu_{\mathrm{ei}}^{a} ( {\bf r}_{i} ) .\tag{12.50} $$容易验证,式(12.49)的正则化引入的能量误差与拉普拉斯算子离散化引入的误差同阶,因为近似成立$( \mathrm{i} . \mathrm{e} . , \ | \mathbf{r} - \mathbf{R}_{I} | \ < \ a )$区域的体积正比于$a^{3}$,且在此区域中,原始电子-离子势与正则化势的平均差异为${\cal O} ( 1 / a )$。我们还注意到,电子-电子项$V_{\mathrm{ee}} ( \mathbf{x} )$中的正发散不需要特别处理,因为无穷大的排斥势不会影响有限基态能量的存在。
晶格正则化的一个自然定义对应简单立方晶格${\bf r}_{l , m , n} = ( l a , m a , n a )$,其中l、m和n为任意整数。然而,这一方案会随着a的变化产生不可控的误差,因为对于不太简单的系统,离子位置与晶格点之间的最小距离表现得极不规则。为了改进当$a \ \ 0$时的收敛性并定义有效的计算方案,我们注意到算子$\nabla_{a}^{2}$不一定定义在规则晶格上。例如,考虑周期性超胞$L \times L \times L$中的电子。那么,每当晶格间距a与盒子边长L不可通约时(即$L / a$为无理数,或在任何计算机实现中非常接近无理数),$\nabla_{a}^{2}$将产生电子“跳跃”,从而遍及地填充整个周期性超胞的连续空间。一般而言,在开放系统中同样可以产生相同的效应。通过在拉普拉斯算子的定义中使用两个不可通约的晶格间距,例如$\nabla_{a}^{2} p \nabla_{a}^{2} + ( 1 - p ) \nabla_{a^{\prime}}^{2}$,其中$0 < p < 1$且$a^{\prime} / a$固定为无理数,可以定义作用在原始哈密顿量同一连续空间中的${\mathcal{H}}_{a}$,其显著优势是改善了当$a 0$时外推的平滑性(Casula et al., 2005)。在下文中,即使哈密顿量${\mathcal{H}}_{a}$作用于连续空间,我们仍继续将这种方法称为“晶格正则化”。这是因为我们吸收了晶格最重要的优势:能够列举${\mathcal{H}}_{a}$的非对角矩阵元(见下文)以及引入由晶格间距a给出的自然截断。
总之,原始哈密顿量的晶格正则化给出:
$$ \mathcal{H}_{a} = - \frac{1} {2} \nabla_{a}^{2} + V^{a} ( \mathbf{x} ) + V_{\mathrm{ii}} ,\tag{12.51} $$其中势项由下式给出:
$$ V^{a} ( {\mathbf{x}} ) = V_{\mathrm{ee}} ( {\mathbf{x}} ) + V_{\mathrm{ei}}^{a} ( {\mathbf{x}} ) .\tag{12.52} $$我们强调$\mathcal{H}_{a} = \mathcal{H} + O ( a^{2} )$,以及正如我们之前所解释的,$\mathcal{H}$和${\mathcal{H}}_{a}$均作用于同一希尔伯特空间。这种方法最重要的优势是,我们可以将已经讨论过的晶格模型技术应用于正则化哈密顿量${\mathcal{H}}_{a}$。确实,所有这些技术都基于一个事实:将哈密顿量${\mathcal{H}}_{a}$作用于给定构型$| x \rangle$上,会得到有限个$( 6 N_{e} + 1 )$新构型:
$$ \mathcal{H}_{a} | x \rangle = \left[ V^{a} ( {\mathbf{x}} ) + \frac{3 N_{e}} {a^{2}} \right] | x \rangle - \frac{1} {2 a^{2}} \sum_{j = 1}^{6 N_{e}} | x_{j} \rangle ,\tag{12.53} $$其中右端的第二项来源于,根据式(12.48),对于每个笛卡尔分量,我们必须将每个电子位置改变 $\pm a$;在这种情况下,每个单独的改变将产生一个新构型 $\left| x_{j} \right\rangle$。
然后,类似于格点情形(见第10章),固定节点近似可以应用于 ${\mathcal{H}}_{a}$。这里,固定节点哈密顿量 $\mathcal{H}_{a}^{\mathrm{FN} , \gamma}$ 可以通过引导函数 $\Psi_{G} ( \mathbf{x} )$ 来定义,该函数用于固定基态的节点面,见式(10.34)。在该近似下,${\mathcal{H}}_{a}$ 中导致符号问题的非对角矩阵元乘以 $-\gamma$,其中 $\gamma$ 为任意正数,而对角元则按如下方式变化:
$$ V^{a} ( {\bf x} ) \rightarrow V^{a} ( {\bf x} ) + ( 1 + \gamma ) V_{\mathrm{sf}}^{a} ( {\bf x} ) ,\tag{12.54} $$其中
$$ V_{\mathrm{sf}}^{a} ( {\bf x} ) = - \frac{1} {2 a^{2}} \sum_{j : s ( {\bf x} , {\bf x}_{j} ) < 0} \frac{\Psi_{G} ( {\bf x}_{j} )} {\Psi_{G} ( {\bf x} )} ,\tag{12.55} $$其中 $s ( \mathbf{x} , \mathbf{x}_{j} ) = \Psi_{G} ( \mathbf{x}_{j} ) \Psi_{G} ( \mathbf{x} )$。讨论该方法在 $a \to 0$ 极限下的推论是很有趣的。由于条件 $\Psi_{G} ( \mathbf{x}_{j} ) \Psi_{G} ( \mathbf{x} ) < 0$ 仅在非常靠近节面(即距离 $\approx a$ 处,其中 $\Psi_{G} ( \mathbf{x} ) = 0$)时才成立,因此 $\mathcal{H}_{a}^{\mathrm{FN} , \gamma}$ 的矩阵元几乎总是与 ${\mathcal{H}}_{a}$ 的矩阵元一致。出于同样的原因,势能 $V_{\mathrm{sf}}^{a} ( {\bf x} )$ 仅在节面附近非零。因此,我们得到,对于 $a \to 0$,固定节点哈密顿量将通过节面处的边界条件与正则化哈密顿量区分开来,在该边界条件下,$\mathcal{\dot{H}}_{a}^{\mathrm{FN} , \gamma}$ 的固定节点基态波函数为零。事实上,势能 $V_{\mathrm{sf}}^{a} ( {\bf x} )$ 在 $a \to 0$ 时发散,而节面上不为零的波函数在 $a \to 0$ 时,其 $\mathcal{H}_{a}^{\mathrm{FN} , \gamma}$ 的期望值将为无穷大(正数)。因此,在这个极限下,我们清楚地看到,具有晶格正则化的固定节点近似等价于一个边界条件,并提供了与导引函数节点相容的最佳变分能量。请注意,这一事实是在没有使用所谓的铺砌定理(Ceperley,1991)的情况下得到的,因为对于 $\gamma > 0$ 和 $a > 0$,晶格上的固定节点蒙特卡罗采样总是遍历的,并且 $\mathcal{\hat{H}}_{a}^{\mathrm{FN} , \bar{\gamma}}$ 的基态在所有节点袋中都有定义。特别地,$\gamma = 1$ 的模拟是最有效的,而 $\gamma = 0$ 的情况则不是遍历的,因为节点袋无法被穿越。
然后,固定节点方法可以在格林函数蒙特卡罗算法中实现;然而,由于哈密顿量的谱无上界,有必要考虑已在第8.4节中描述的连续时间极限。通过这种方式,我们得到了近似的基态,它满足由导引函数 $\Psi_{G} ( \mathbf{x} )$ 给出的固定节点约束,类似于标准扩散蒙特卡罗框架(Casula 等人,2005)。我们指出,传统的扩散蒙特卡罗技术存在关于 $\Delta \tau$ 的时间步长误差,而当前的公式则引入了晶格空间误差;两种方法都具有相同的上界性质,并分别在 $\Delta \tau \to 0$ 和 $a \to 0$ 的极限下,收敛到相同的投影固定节点能量。
12.6 晶格正则化的改进方案#
不幸的是,尽管可以随着 $a \to 0$ 获得非常平滑的收敛,但这种简单方法在实践中并不非常有用,因为误差通常很大,这意味着需要极其小的 a 值才能获得精确的结果。为了改进 $a \to 0$ 时的收敛性,我们应该考虑一种替代方法,其中拉普拉斯算符仍然使用式(12.48)进行离散化,但采用势能算符的不同版本,即 $V ( \mathbf{x} ) \to V_{\mathrm{zv}}^{a} ( \mathbf{x} )$。实际上,受晶格上固定节点近似的启发,我们可以通过要求原始哈密顿量的局部能量与正则化哈密顿量的局部能量一致来修改势能算符:
$$ \frac{\langle x | \mathcal{H}_{a} | \Psi_{G} \rangle} {\langle x | \Psi_{G} \rangle} = \frac{\langle x | \mathcal{H} | \Psi_{G} \rangle} {\langle x | \Psi_{G} \rangle} .\tag{12.56} $$假设该等式对于最佳引导函数 $| \Psi_{G} \rangle$(即对应最低变分能量的引导函数)成立。这是一个有用的条件,因为它满足了零方差性质,即使在存在晶格正则化的情况下也是如此:如果引导函数与精确本征态重合,那么 ${\mathcal{H}}_{a}$ 的局域能量就等于精确本征值,相应的固定节点能量也将是精确的,无论 $^{a}$ 的大小如何。于是,式 (12.56) 隐含了 $V_{\mathrm{zv}}^{a} ( \mathbf{x} )$ 的唯一确定。实际上,如果我们假设拉普拉斯算符通过替换 $\nabla^{2} \to \nabla_{a}^{2}$ 进行离散化,那么 $V_{\mathrm{ZV}}^{a} ( \mathbf{x} )$ 由下式给出:
$$ V_{\mathrm{zv}}^{a} ( \mathbf{x} ) = V ( \mathbf{x} ) + \frac{\left( \nabla_{a}^{2} - \nabla^{2} \right) \Psi_{G} ( \mathbf{x} )} {2 \Psi_{G} ( \mathbf{x} )} ,\tag{12.57} $$该式表明,当 $a \to 0$ 时 $V_{\mathrm{zv}}^{a} ( {\bf x} ) \to V ( {\bf x} )$。这一选择最初由 Casula 等人(2005)提出。在他们的工作中,模拟非常稳定,并且极限 $a \to 0$ 可以很容易地提取出来;不幸的是,这些计算对于费米子通常定义不佳,因为哈密顿量 ${\mathcal{H}}_{a}$ 对于 $a > 0$ 并不总是具有有限的基态能量。实际上,正则化势 $V_{\mathrm{ZV}}^{a} ( \mathbf{x} )$ 在引导函数 $\Psi_{G} ( \mathbf{x} )$ 的节点面附近可能取无界的负值,而晶格正则化的动能(其矩阵元总是有限的)无法补偿这些负值。然而,由于采用统计方法很难选择导致不稳定性的节点面构型,使用有限计算资源的模拟仍可能产生合理的结果。尽管如此,不能先验地排除出现不稳定性的可能性,因此找到解决这种潜在灾难的方法至关重要。类似的情况也常见于标准的扩散蒙特卡罗技术(Foulkes 等人,2001;Umrigar 等人,1993),尤其是在使用赝势时。在后一种情况中,对于给定的 $\Delta \tau$,尽管普遍认为只要 $\Delta \tau$ 足够小并采用适当的截断,模拟就应该是稳定的,但目前尚无通用的补救措施。在晶格正则化方法中,情况更加透明,并且可以找到简单的补救措施,我们将在下文描述。
为了定义一种对于任意有限 $^{a}$ 值都严格稳定且稳健的方法,我们需要修改电子-离子相互作用。具体来说,我们注意到式 (12.57) 的正则化实际上对应于通过替换式 (12.45) 中的电子-离子势来实现修改:
$$ \nu_{\mathrm{ei}} ( {\bf r}_{i} ) \to \nu_{\mathrm{zv} , i}^{a} ( \mathbf{x} ) ,\tag{12.58} $$其中
$$ \nu_{\mathrm{zv} , i}^{a} ( \mathbf{x} ) = \nu_{\mathrm{ei}} ( \mathbf{r}_{i} ) + \frac{\left( \nabla_{i , a}^{2} - \nabla_{i}^{2} \right) \Psi_{G} ( \mathbf{x} )} {2 \Psi_{G} ( \mathbf{x} )} .\tag{12.59} $$然后,对于满足电子-离子尖点条件(如第 12.2 节所述)的良好引导函数,上述替换抵消了吸引势的大部分奇异性;然而,在引导函数 $\Psi_{G} ( \mathbf{x} )$ 的节点面上,我们仍然可能得到无界的负值。一种稳妥地保护模拟免受引导函数节点面所有可能不稳定性的方法是,对每个电子 $i$ 的电子-离子势进行如下修改:
$$ \begin{array} {r} {\nu_{\mathrm{max} , i}^{a} ( \mathbf{x} ) = \mathrm{Max} \left[ \nu_{\mathrm{zv} , i}^{a} ( \mathbf{x} ) , \nu_{\mathrm{ei}}^{a} ( \mathbf{r}_{i} ) \right] ,} \end{array}\tag{12.60} $$其中 $\nu_{\mathrm{ei}}^{a} ( \mathbf{r}_{i} )$ 由式(12.49)定义。在此框架下,正则化势能可避免节点面附近可能出现的失稳问题,并保证对于任意 $a > 0$ 都存在有限的基态能量。由于条件 $\nu_{\mathrm{ei}}^{a} ( \mathbf{r}_{i} ) > \nu_{\mathrm{zv} , i}^{a} ( \mathbf{x} )$ 经常成立,因此我们可以改进这一选择,仅当电子非常接近节点面时应用条件(12.60)。为此,我们仅在电子在离散拉普拉斯算子作用下穿越节点面时应用(单电子)条件(12.60):
$$ \begin{array} {r} {\nu_{\mathrm{opt} , i}^{a} ( \mathbf{x} ) = \left\{\begin{array} {l l} {\nu_{\mathrm{max} , i}^{a} ( \mathbf{x} )} & {\mathrm{if} \Psi_{G} ( \mathbf{x} ) \Psi_{G} ( \mathbf{x}_{j} ) < 0 \mathrm{for} \mathrm{at} \mathrm{least} \mathrm{one} \mathbf{x}_{j} ,} \\ {\nu_{\mathrm{zv} , i}^{a} ( \mathbf{x} )} & {\mathrm{otherwise} ;} \end{array} \right.} \end{array}\tag{12.61} $$这里,$\mathbf{X}_{j}$ 表示由电子 $i$ 的离散拉普拉斯算子 $\nabla_{a , i}^{2}$ 生成的可能构型。注意,对于有限的晶格间距 $a$,上述条件可能会遗漏发生双重符号变化的情况。然而,当 $a$ 值较小时,这种情况的影响可以忽略不计。图12.2报告了前述方案在一维链上放置的18个等间距氢原子(最近邻原子距离 $R = 2 a.u.$)情形下的计算结果。这里,引导函数由优化后的Jastrow因子(采用3s2p1d基组)应用于从密度泛函计算(采用4s2p1d基组)得到的Slater行列式构成。最简方法,即式(12.57)不经额外正则化的方案,已给出相当准确的结果。通过施加式(12.60)的修正,我们得到了高得多的能量,表明这种方案过于严格。最终,采用式(12.61)的方案得到了非常准确且平滑的能量,表明该方法为晶格正则化扩散蒙特卡洛方法提供了最佳近似。
只要引导函数具有尺寸一致性,这种晶格正则化方案就完全满足尺寸一致性。事实上,考虑一个大距离下可分解的波函数 $\Psi_{G} ( \mathbf{x} )$:
$$ \Psi_{G} ( \mathbf{x} ) \approx \Psi_{G}^{A} ( \mathbf{x}_{A} ) \Psi_{G}^{B} ( \mathbf{x}_{B} ) ,\tag{12.62} $$其中我们假设相关电子构型 $\mathbf{x}$ 由区域A(B)内对应于构型 $\mathbf{X}_{A} ~ \left( \mathbf{X}_{B} \right)$ 的 $N_{A}$($N_{B}$)个电子得到。

图12.2 不同正则化方案下,一维链上18个氢原子(即18个电子)在距离 $R = 2 a.u.$ 处的总能量结果。此处使用的引导函数由优化后的Jastrow因子(采用 $3 s 2 p 1 d$ 基组)应用于从密度泛函计算(采用 $4 s 2 p 1 d$ 基组)得到的Slater行列式构成。报告了式(12.57)的简单近似(空心三角形)、式(12.60)给出的结果(空心正方形)以及式(12.61)给出的结果(实心圆)。$a 0$ 外推能量为 $E = -10.1925(2)$ 哈特里。
那么,式(12.59)的正则化势能将独立作用在两个区域中:
$$ \begin{array} {r} {\nu_{z \mathrm{v} , i}^{a} ( \mathbf{x} ) = \left\{\begin{array} {l l} {\displaystyle \nu_{ei} ( \mathbf{r}_{i} ) + \frac{\left( \nabla_{i , a}^{2} - \nabla_{i}^{2} \right) \Psi_{G}^{A} ( \mathbf{x}_{A} )} {2 \Psi_{G}^{A} ( \mathbf{x}_{A} )}} & {\mathrm{~ i f ~} \mathbf{r}_{i} \in A ,} \\ {\displaystyle \nu_{ei} ( \mathbf{r}_{i} ) + \frac{\left( \nabla_{i , a}^{2} - \nabla_{i}^{2} \right) \Psi_{G}^{B} ( \mathbf{x}_{B} )} {2 \Psi_{G}^{B} ( \mathbf{x}_{B} )}} & {\mathrm{~ i f ~} \mathbf{r}_{i} \in B ,} \end{array} \right.} \end{array}\tag{12.63} $$从而得到:
$$ \nu_{\mathrm{max} , i}^{a} ( \mathbf{x} ) = \nu_{\mathrm{max} , A , i}^{a} ( \mathbf{x}_{A} ) + \nu_{\mathrm{max} , B , i}^{a} ( \mathbf{x}_{B} ) .\tag{12.64} $$这个事实解释了为什么我们采用了一种正则化方法,该方法独立考虑每个单电子贡献,参见方程(12.60)。相比之下,一种采用 $V_{\mathrm{max}}^{a} ( \mathbf{x} ) = \mathrm{Max} \left[ V_{\mathrm{zv}}^{a} ( \mathbf{x} ) , V^{a} ( \mathbf{x} ) \right] )$ 的正则化方法(其中 $V^{a} ( {\bf x} )$ 由方程(12.50)给出)是不满足大小一致性的,因为 Max 操作所隐含的全局关系必然会耦合 A 和 B 两个区域,即使它们相距很远。
最后,晶格正则化方法可以自然地扩展到赝势的使用,赝势通常被认为可以消除来自芯电子的巨大能量尺度。赝势具有一些项,它们通过非局域贡献简单地叠加到 $\nu_{ei} ( {\bf r}_{i} )$ 上,而利用晶格算符可以对这种非局域贡献进行自然的描述。事实上,它们可以被安全地视为需要包含在正则化拉普拉斯算子 $\nabla_{a}^{2}$ 中的额外跳跃项。因此,应用前面描述的固定节点方案(The fixed-node scheme described before)是直接的,并且能够保持其能量的上界性质。此外,使用软赝势(它能够消除电子-离子奇异性)会得到一个光滑的局域势,该势不需要进行正则化,因此方程(12.49)的正则化步骤也无需应用。
附录 计算机生成的伪随机数#
在大多数情况下,计算机生成的随机数并非真正的随机数,而是伪随机数,这意味着它们是通过使用精确的确定性算法以可预测的方式生成的。事实上,尽管存在纯随机的算法(基于计算机中的某些大气噪声),但在许多情况下,使用能够确定可重复计算的代码要好得多。这是因为能够运行相同的随机数序列有助于调试。此外,基于计算机中某些热噪声的随机数生成器通常比确定性算法慢得多,这也提供了另一个避免使用它们的理由。所有产生伪随机数的例程都是复杂的算法,它们将大量数论中的巧妙技巧浓缩在几百行代码中。这些算法的输出看起来是随机的,但实际上并非如此。事实上,当我们(在相同的初始条件下)运行它们两次时,它们总会给出相同的输出。在实践中,这些算法生成一个数列 $r_{0} , r_{1} , r_{2} , . . .$ ,在每次迭代中,$r_{k}$ 是由前一个数的某个复杂函数给出的。最简单的例子是线性同余生成器,其形式为:
$$ r_{k} = f ( r_{k - 1} ) .\tag{A.1} $$一个简单的(不应再使用的)选择由下式给出:
$$ i_{k} = ( a i_{k - 1} + b ) \mod \mathrm{M} ,\tag{A.2} $$其中 $a$ 和 $b$ 是合适的整数(生成器的质量关键取决于它们),而 $M = 2^{31} - 1$(用于32位生成器)。然后通过取 $r_{k} = i_{k} / M$ 来获得 $[0, 1)$ 区间内的实数。从一个给定的种子(这里即 $i_{0}$ 的初始值)开始,整个伪随机数链被迭代地生成。这里的主要问题是,该序列在(至多)$2^{31} - 1$ 次迭代后会完全相同地重复。从这个意义上说,滞后斐波那契生成器给出了更好的结果:
$$ i_{k} = i_{k - p} \bigotimes i_{k - q} \mod{\mathrm{M}} ,\tag{A.3} $$其中 $\otimes$ 表示任意二元运算,$p$ 和 $q$ 是固定整数。根据所选的二元运算及整数 $p$ 和 $q$,我们可以实现周期极长的生成器,例如周期可达 $10^{400}$ 甚至更长。
在处理伪随机数时,我们必须牢记:由于这些数字来自确定性过程,相邻数字之间总是存在相关性。然而,关键在于,在不知道算法的情况下,它们至少在"粗粒化"层面上必须看起来几乎是随机的。更准确地说,条件概率 $\omega ( i_{k} | i_{k - 1} )$ 应尽可能平坦,这意味着获得某个给定 $i_{k}$ 的概率应与上一步获得 $i_{k - 1}$ 这一事实无关。对于像式 (A.2) 这样的算法,这显然不成立,因为在第 $k - 1$ 步的结果完全决定了下一步的结果。然而,一旦我们对区间 [0, 1) 进行粗粒化处理,就可能出现一个常数化的 $\omega ( i_{k} | i_{k - 1} )$。此外,对于滞后斐波那契生成器,只要其周期远大于 $2^{31} - 1 \approx 2 \times 10^{9}$,那么在序列中,相同的 $i_{k - 1}$ 将对应多个不同的 $i_{k}$ 值,从而导致概率实际上均匀分布。为了示例化这一概念,我们将区间 [0, 1) 划分为 $L$ 个宽度为 $1 / L$ 的小子区间,生成 $M$ 个伪随机数,并计算当前一个数落在 $[ ( m - 1 ) / L , m / L ]$ 时,后一个数落在 $[ ( n - 1 ) / L , n / L ]$ 的次数 $N ( n | m )$。于是,估计的条件概率由 $\begin{array} {r} {\overline{{\omega}} ( n | m ) = N ( n | m ) / \sum_{n} N ( n | m )} \end{array}$ 给出,这本身也是一个伪随机数。对于好的生成器,根据中心极限定理,${\overline{{\omega}}} ( n | m )$ 的概率分布是均值为 $1 / L$、方差为 $1 / M$ 的高斯分布(${\overline{{\omega}}} ( n | m )$ 的均值是 $1 / L$,方差是 $1 / L + O ( 1 / L^{2} )$;每个 m 的采样次数平均为 $O ( M / L )$)。在图 A.1 中,我们报告了 $L = 1000$ 和 $M = 10^{9}$ 的情况,其中 ${\overline{{\omega}}} ( n | m )$ 的分布函数完美地拟合为高斯分布。

图 A.1 区间 [0, 1) 被划分为宽度为 $1 / L$ 的子区间,其中 $L = 1000$。评估了两个后续数字分别落在子区间 $[ ( n - 1 ) / L , n / L ]$ 和 $[ ( m {-} 1 ) / L , m / L ]$ 的条件概率 ${\overline{{\omega}}} ( n | m )$,并对 $M = 10^{9}$ 个样本进行了平均。图中显示了 $\overline{{\omega}} ( n | m )$ 的分布函数。
一般来说,随机数生成器必须满足两个主要要求(Knuth, 1997):
具有长周期。事实上,每个确定性生成器最终都会循环,目标是使周期尽可能长。
具有良好的统计质量。在这方面,生成器的输出应与真正的随机数生成器的输出几乎无法区分,并且在数字序列中不应表现出任何相关性。
生成器质量差可能会破坏蒙特卡洛应用的结果(Ferrenberg et al., 1992);因此,生成器能够通过一组经验测试非常重要,例如由 George Marsaglia 开发的 Diehard 测试集,或由 Pierre L'Ecuyer 和 Richard Simard 引入的 TestU01 库。
索引#
组态相互作用, 20 约束路径量子蒙特卡洛, 219 连续时间, 179, 196 关联采样, 59 关联能, 243 关联时间, 73 协方差, 45 累积概率, 43 尖点条件, 239 延迟更新, 121, 227 直接采样, 60 分布函数 伯努利分布, 46 二项分布, 46 高斯分布, 46 泊松分布, 46 均匀分布, 46 能量导数, 134 误差棒 分块分析, 83 自助法, 78 误差传播, 77 刀切法, 81 费米子波函数 行列式, 112 行列式的快速更新, 116 普法夫量的快速更新, 126 格林函数, 119, 128 普法夫量, 123 固定节点近似, 204, 209, 255 固定节点能量 期望值, 212 基态值, 211 上界性质, 211 固定节点误差, 212 固定节点哈密顿量, 210, 255 固定相位近似, 204 福克-普朗克方程, 90 趋近平衡, 92 薛定谔方程, 91 格拉姆-施密特正交化, 222 格林函数蒙特卡洛 累积权重, 173, 176, 184 分支过程, 181 连续时间, 179 前向行走, 176, 186 导引函数, 173 重要性采样, 173 主方程, 169 混合平均, 177 符号问题, 168 单游走子技术, 170 古茨威勒因子, 17 古茨威勒投影算符, 18, 107 古茨威勒波函数, 17 哈丹-沙斯特里波函数, 23 哈特里-福克波函数, 4, 15 海森堡模型, 10 海特勒-伦敦波函数, 3 赫尔曼-费曼定理, 212 海森矩阵, 148 哈伯德模型(玻色子), 12 哈伯德模型(费米子), 7 哈伯德-斯特拉托诺维奇变换, 218 虚时传播, 35, 202, 205, 214 重要性采样, 61 刀切法, 81 贾斯特罗因子 单体(连续空间中的电子-离子), 239 两体(连续空间中的电子-电子), 239 两体(晶格上的电子-电子), 19, 107 贾斯特罗波函数, 18, 107 贾斯特罗-斯莱特波函数, 20, 107, 237 朗之万动力学 加速动力学, 96 细致平衡, 92 离散化, 87 离散化误差, 95 一阶方程, 85, 140 主方程, 90 二阶方程, 86 大偏差, 55 晶格正则化, 235 电子-离子势, 253 拉普拉斯算子, 253 尺寸一致性, 258 零方差性质, 256 劳克林波函数, 25, 237 光锥, 163 线性方法, 132, 147, 248 “强”零方差性质, 149 广义本征值方程, 149 局域能量, 105, 172, 175, 191 局域估计量, 104 马尔可夫链 辅助场量子蒙特卡洛, 223 定义, 66 细致平衡, 70 遍历性, 70 格林函数蒙特卡洛, 169, 200 主方程, 67 周期性, 70 可约性, 70 爬行动物蒙特卡洛, 191 平稳解, 70 均值, 44 梅特罗波利斯-黑斯廷斯算法, 74 接受概率, 74 试验概率, 74 分子轨道, 243 莫特绝缘体, 9 莫特相变, 10 多行列式波函数, 20 牛顿-拉弗森优化, 132 节面, 202 优化 线性方法, 147 重加权技术, 132 最速下降法, 137 随机重构, 139 配对(孪生)函数, 243 粒子-空穴变换, 115, 229 路径积分蒙特卡洛, 190, 220 泡利面, 202 佩龙-弗罗贝尼乌斯定理, 72, 199, 210 普法夫量, 123 幂法, 34 概率 条件概率, 42 定义(频率学派), 41 联合概率, 42 边缘概率, 42 实验的可重复性, 40, 52 概率密度, 43 投影技术, 34, 167, 189, 214, 253 伪随机数 定义, 261 斐波那契(滞后)生成器, 261 周期性, 261 量子淬火, 161 随机过程, 66 随机变量, 42 实时传播, 156 爬行动物蒙特卡洛 弹跳算法, 195 连续时间, 196 重要性采样, 190 主方程, 195 伪配分函数, 189 符号问题, 190 简单更新, 191 共振价键波函数, 5, 29, 237, 244 重加权技术, 59 采样, 56 连续分布, 64 离散分布, 62 符号问题, 199 尺寸一致性, 244 尺寸广延性, 31 斯莱特行列式, 4 稀疏哈密顿量矩阵, 35 标准差, 45 最速下降法, 88, 132, 137, 248 随机重构, 132, 139, 248 协方差性质, 142 预处理, 143 投影, 146 正则化, 143, 147 尺度不变正则化, 144 信噪比, 144
索引#
热化时间, 73, 172, 176, 193 含时哈特里-福克, 156 含时变分蒙特卡洛 连续时间极限, 159 能量守恒, 160 方程, 158 模守恒, 160 平稳作用量, 160 Trotter 近似, 216 Trotter 误差, 180, 196, 216 范德蒙行列式, 22 方差, 45 变分哈密顿量, 209 变分蒙特卡洛 期望值, 104 局部基组, 108 局域能量, 105 局部估计量, 104 零方差性质, 105 变分原理 关联性精度, 14 能量精度, 13 一般性质, 13
波函数度量, 139 大数弱定律, 52 白噪声, 85 Wick 定理, 119, 128, 224
参考文献#
Allen, M. P., and Tildesley, D. J. 1987. Computer Simulation of Liquids. Oxford University Press.
Anderson, H. 1986. Metropolis, Monte Carlo and the MANIAC. Los Alamos Science, 14, 96–108.
Anderson, J. B. 1975. Random walk simulation of the Schrodinger equation:¨ $\mathrm{He}_{3}^{+}$ . J. Chem. Phys., 63, 1499–1503.
Anderson, J. B. 1976. Quantum chemistry by random walk. J. Chem. Phys., 65, 4121–4127.
Anderson, P. W. 1987. The resonating valence bond state in ${\mathrm{La}}_{2} {\mathrm{CuO}}_{4}$ and superconductivity. Science, 235, 1196–1198.
Anderson, P. W., Baskaran, G., Zou, Z., and Hsu, T. 1987. Resonating-valence-bond theory of phase transitions and superconductivity in ${\mathrm{La}}_{2} {\mathrm{CuO}}_{4}$ -based compounds. Phys. Rev. Lett., 58, 2790–2793.
Aoki, H., Tsuji, N., Eckstein, M., Kollar, M., Oka, T., and Werner, P. 2014. Non-equilibrium dynamical mean-field theory and its applications. Rev. Mod. Phys., 86, 779–837.
Arovas, D., Schrieffer, J. R., and Wilczek, F. 1984. Fractional statistics and the quantum Hall effect. Phys. Rev. Lett., 53, 722–723.
Bajdich, M., Mitas, L., Drobny, G., Wagner, L. K., and Schmidt, K. E. 2006. Pfaffian pairing wave functions in electronic-structure quantum Monte Carlo simulations. Phys. Rev. Lett., 96, 130201.
Bajdich, M., Mitas, L., Wagner, L. K., and Schmidt, K. E. 2008. Pfaffian pairing and backflow wave functions for electronic-structure quantum Monte Carlo methods. Phys. Rev. B, 77, 115112.
Baldereschi, A. 1973. Mean-value point in the Brillouin zone. Phys. Rev. B, 7, 5212–5215.
Bardeen, J., Cooper, L. N., and Schrieffer, J. R. 1957. Theory of superconductivity. Phys. Rev., 108, 1175–1204.
Baroni, S., and Moroni, S. 1998. Reptation quantum Monte Carlo. arXiv:cond-mat/ 9808213.
Bartlett, R. J. 1981. Many-body perturbation theory and coupled cluster theory for electron correlation in molecules. Ann. Rev. Phys. Chem., 32, 359–401.
Baskaran, G., and Anderson, P. W. 1988. Gauge theory of high-temperature superconductors and strongly correlated Fermi systems. Phys. Rev. B, 37, 580–583.
Bethe, H. 1931. Zur Theorie der metalle. I. eigenwerte und eigenfunktionen der linearen atomkette. Z. Phys., 71, 205–226.
Bouchaud, J. P., Georges, A., and Lhuillier, C. 1988. Pair wave functions for strongly correlated fermions and their determinantal representation. J. Phys. (Paris), 49, 553–559.
Calandra Buonaura, M., and Sorella, S. 1998. Numerical study of the two-dimensional Heisenberg model using a Green function Monte Carlo technique with a fixed number of walkers. Phys. Rev. B, 57, 11446–11456.
Capello, M., Becca, F., Fabrizio, M., Sorella, S., and Tosatti, E. 2005. Variational description of Mott insulators. Phys. Rev. Lett., 94, 026406.
Capello, M., Becca, F., Yunoki, S., and Sorella, S. 2006. Unconventional metal-insulator transition in two dimensions. Phys. Rev. B, 73, 245116.
Capello, M., Becca, F., Fabrizio, M., and Sorella, S. 2007. Superfluid to Mott-insulator transition in Bose-Hubbard models. Phys. Rev. Lett., 99, 056402.
Capello, M., Becca, F., Fabrizio, M., and Sorella, S. 2008. Mott transition in bosonic systems: insights from the variational approach. Phys. Rev. B, 77, 144517.
Capriotti, L., Becca, F., Parola, A., and Sorella, S. 2001. Resonating valence bond wave functions for strongly frustrated spin systems. Phys. Rev. Lett., 87, 097201.
Carleo, G., Becca, F., Schiro, M., and Fabrizio, M. 2012. Localization and glassy dynamics of many-body quantum system. Sci. Rep., 2, 243.
Carleo, G., Becca, F., Sanchez-Palencia, L., Sorella, S., and Fabrizio, M. 2014. Light-cone effect and supersonic correlations in one- and two-dimensional bosonic superfluids. Phys. Rev. A, 89, 031602.
Casula, M., Filippi, C., and Sorella, S. 2005. Diffusion Monte Carlo method with lattice regularization. Phys. Rev. Lett., 95, 100201.
Ceperley, D., Chester, G. V., and Kalos, M.H. 1977. Monte Carlo simulation of a manyfermion study. Phys. Rev. B, 16, 3081–3099.
Ceperley, D. M. 1991. Fermion nodes. J. Stat. Phys., 63, 1237–1267.
Ceperley, D. M., and Alder, B. J. 1980. Ground state of the electron gas by a stochastic method. Phys. Rev. Lett., 45, 566–569.
Corboz, P., Rice, T. M., and Troyer, M. 2014. Competing states in the t−J model: uniform d-wave state versus stripe state. Phys. Rev. Lett., 113, 046402.
Daley, A. J., Kollath, C., Schollwock, U., and Vidal, G. 2004. Time-dependent density- ¨ matrix renormalization-group using adaptive effective Hilbert spaces. J. Stat. Mech.: Theor. Exp., P04005.
Dirac, P. A. M. 1930. Note on exchange phenomena in the Thomas atom. Math. Proc. Cambridge Philos. Soc., 26, 376–385.
Dongarra, J., and Sullivan, F. 2000. Guest editors’ introduction: the top 10 algorithms. Comput. Sci Eng., 2, 22–23.
Dovesi, R., Orlando, R., Erba, A., Zicovich-Wilson, C. M., Civalleri, B., Casassa, S., Maschio, L., Ferrabone, M., De La Pierre, M., D’Arco, P., Noel, Y., Caus ¨ a, M., R ` erat, ´ M., and Kirtman, B. 2014. CRYSTAL14: a program for the ab initio investigation of crystalline solids. Int. J. Quantum Chem., 114, 1287–1317.
Eichenberger, D., and Baeriswyl, D. 2007. Superconductivity and antiferromagnetism in the two-dimensional Hubbard model: a variational study. Phys. Rev. B, 76, 180504.
Eisert, J., Cramer, M., and Plenio, M. B. 2010. Area laws for the entanglement entropy. Rev. Mod. Phys., 82, 277–306.
Faddeev, L. D., and Takhtajan, L. A. 1981. What is the spin of a spin wave? Phys. Lett. A, 85, 375–377.
Fahy, S., and Hamann, D. R. 1991. Diffusive behavior of states in the Hubbard-Stratonovich transformation. Phys. Rev. B, 43, 765–779.
Fano, G., Ortolani, F., and Colombo, E. 1986. Configuration-interaction calculations on the fractional quantum Hall effect. Phys. Rev. B, 34, 2670–2680.
Fazekas, P. 1999. Lecture Notes on Electron Correlation and Magnetism. World Scientific.
Fazekas, P., and Anderson, P. W. 1974. On the ground state properties of the anisotropic triangular antiferromagnet. Phil. Mag., 30, 423–440.
Ferrenberg, A. M., Landau, D. P., and Wong, Y. J. 1992. Monte Carlo simulations: Hidden errors from “good” random number generators. Phys. Rev. Lett., 69, 3382–3384.
Fetter, A. L., and Walecka, J. D. 2003. Quantum Theory of Many-Particle Systems. Dover Publications Inc.
Feynman, R. P. 1954. Atomic theory of the two-fluid model of liquid Helium. Phys. Rev., 94, 262–277.
Feynman, R. P., and Cohen, M. 1956. Energy spectrum of the excitations in liquid Helium. Phys. Rev., 102, 1189–1204.
Filippi, C., and Umrigar, C. 1996. Multiconfiguration wave functions for quantum Monte Carlo calculations of first-row diatomic molecules. J. Chem. Phys., 105, 213–226.
Foulkes, W. M. C., Mitas, L., Needs, R. J., and Rajagopal, G. 2001. Quantum Monte Carlo simulations of solids. Rev. Mod. Phys., 73, 33–83.
Gnedenko, B. V. 2014. The Theory of Probability. Martino Fine Books.
Gnedenko, B. V., and Kolmogorov, A. N. 1954. Limit Distributions for Sums of Independent Random Variables. Addison-Wesley.
Gros, C. 1988. Superconductivity in correlated wave functions. Phys. Rev. B, 38, 931–934.
Gros, C. 1989. Physics of projected wave functions. Ann. Phys., 189, 53–88.
Gros, C., Joynt, R., and Rice, T. M. 1987. Antiferromagnetic correlations in almostlocalized Fermi liquids. Phys. Rev. B, 36, 381–393.
Gubernatis, J., Kawashima, N., and Werner, P. 2016. Quantum Monte Carlo methods: algorithms for lattice models. Cambridge University Press.
Gutzwiller, M. C. 1963. Effect of correlation on the ferromagnetism of transition metals. Phys. Rev. Lett., 10, 159–162.
Haldane, F. D. M. 1983. Fractional quantization of the Hall effect: a hierarchy of incompressible quantum fluid states. Phys. Rev. Lett., 51, 605–608.
Haldane, F. D. M. 1988. Exact Jastrow-Gutzwiller resonating-valence-bond ground state of the spin-1/2 antiferromagnetic Heisenberg chain with $1 / r^{2}$ exchange. Phys. Rev. Lett., 60, 635–638.
Haldane, F. D. M. 1991. “Spinon gas” description of the $S = 1 / 2$ Heisenberg chain with inverse-square exchange: exact spectrum and thermodynamics. Phys. Rev. Lett., 66, 1529–1532.
Haldane, F. D. M., and Rezayi, E. H. 1985. Periodic Laughlin-Jastrow wave functions for the fractional quantized Hall effect. Phys. Rev. B, 60, 2529–2531.
Harju, A., Barbiellini, B., Siljamaki, S., Nieminen, R. M., and Ortiz, G. 1997. Stochastic ¨ gradient approximation: an efficient method to optimize many-body wave functions. Phys. Rev. Lett., 79, 1173–1177.
Hastings, W. K. 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57, 97–109.
Heitler, W., and London, F. 1927. Wechselwirkung neutraler atome und homoopolare¨ bindung nach der quantenmechanik. Z. Phys., 44, 455–472.
Hirsch, J. E. 1985. Two-dimensional Hubbard model: numerical simulation study. Phys. Rev. B, 31, 4403–4419.
Hu, W.-J., Becca, F., and Sorella, S. 2012. Absence of static stripes in the two-dimensional t−J model determined using an accurate and systematic quantum Monte Carlo approach. Phys. Rev. B, 85, 081110.
Hubbard, J. 1959. Calculation of partition functions. Phys. Rev. Lett., 3, 77–78.
Hubbard, J. 1963. Electron correlations in narrow energy bands. Proc. Royal Soc. of London, 276, 238–257.
Ido, K., Ohgoe, T., and Imada, M. 2015. Time-dependent many-variable variational Monte Carlo method for non-equilibrium strongly correlated electron systems. Phys. Rev. B, 92, 245106.
Imada, M., Fujimori, A., and Tokura, Y. 1998. Metal-insulator transitions. Rev. Mod. Phys., 70, 1039–1263.
Iqbal, Y., Becca, F., and Poilblanc, D. 2011. Projected wave function study of $Z_{2}$ spin liquids on the kagome lattice for the spin-1/2 quantum Heisenberg antiferromagnet. Phys. Rev. B, 84, 020407.
Jain, J. K. 2012. Composite Fermions. Cambridge University Press.
Jastrow, R. 1955. Many-body problem with strong forces. Phys. Rev., 98, 1479–1484.
Kalos, M. H., Levesque, D., and Verlet, L. 1974. Helium at zero temperature with hardsphere and other forces. Phys. Rev. A, 9, 2178–2195.
Kanamori, J. 1963. Electron correlation and ferromagnetism of transition metals. Prog. Theor. Phys., 30, 275–289.
Kaneko, R., Tocchio, L. F., Valent´ı, R., Becca, F., and Gros, C. 2016. Spontaneous symmetry breaking in correlated wave functions. Phys. Rev. B, 93, 125127.
Kim, J., Esler, K. P., McMinis, J., Morales, M. A., Clark, B. K., Shulenburger, L., and D. M., Ceperley. 2012. Hybrid algorithms in quantum Monte Carlo. J. Phys.: Conf. Ser., 402, 012008.
Kivelson, S. A., Rokhsar, D. S., and Sethna, J. P. 1987. Topology of the resonating valence-bond state: solitons and high-Tc superconductivity. Phys. Rev. B, 35, 8865– 8868.
Knuth, D. 1997. The Art of Computer Programming. Addison-Wesley.
Krauth, W. 2006. Statistical Mechanics: Algorithms and Computations. Oxford University Press.
Krauth, W., Caffarel, M., and Bouchaud, J. 1992. Gutzwiller wave function for a model of strongly interacting bosons. Phys. Rev. B, 45, 3137–3140.
Kwon, Y., Ceperley, D. M., and Martin, R. M. 1993. Effects of three-body and backflow correlations in the two-dimensional electron gas. Phys. Rev. B, 48, 12037–12046.
Kwon, Y., Ceperley, D. M., and Martin, R. M. 1998. Effects of backflow correlation in the three-dimensional electron gas: quantum Monte Carlo study. Phys. Rev. B, 58, 6800–6806.
Laughlin, R. B. 1983. Anomalous quantum Hall effect: an incompressible quantum fluid with fractionally charged excitations. Phys. Rev. Lett., 50, 1395–1398.
Lee, P. A., Nagaosa, N., and Wen, X.-G. 2006. Doping a Mott insulator: physics of hightemperature superconductivity. Rev. Mod. Phys., 78, 17–85.
Liang, S., Doucot, B., and Anderson, P. W. 1988. Some new variational resonating-valencebond-type wave functions for the spin-1/2 antiferromagnetic Heisenberg model on a square lattice. Phys. Rev. Lett., 61, 365–368.
Lieb, E. H., and Wu, F. Y. 1968. Absence of Mott transition in an exact solution of the short-range, one-band model in one dimension. Phys. Rev. Lett., 20, 1445–1448.
Lin, C., Zong, F. H., and Ceperley, D. M. 2001. Twist-averaged boundary conditions in continuum quantum Monte Carlo algorithms. Phys. Rev. E, 64, 016702.
Loh, E. Y., Gubernatis, J. E., Scalettar, R. T., White, S. R., Scalapino, D. J., and Sugar, R. L. 1990. Sign problem in the numerical simulation of many-electron systems. Phys. Rev. B, 41, 9301–9307.
Marchi, M., Azadi, S., Casula, M., and Sorella, S. 2009. Resonating valence bond wave function with molecular orbitals: Application to first-row molecules. J. Chem. Phys., 131, 154116.
Marshall, W. 1955. Antiferromagnetism. Proc. R. Soc. London Ser., A 232, 48–68.
Martin, R. M. 2004. Electronic Structure: Basic Theory and Practical Methods. Cambridge University Press.
Martin, R.M., Reining, L., and Ceperley, D.M. 2016. Interacting Electrons: Theory and Computational Approaches. Cambridge University Press.
Mazzola, G., and Sorella, S. 2017. Accelerating ab initio molecular dynamics and probing the weak dispersive forces in dense liquid Hydrogen. Phys. Rev. Lett., 118, 015703.
McMillan, W. L. 1965. Ground State of Liquid He4. Phys. Rev., 138, A442–A451.
Metropolis, N. 1987. The beginning of the Monte Carlo method. Los Alamos Science, 15, 125–130.
Metropolis, N., and Ulam, S. 1949. The Monte Carlo method. J. Am. Stat. Ass., 44, 335–341.
Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., and Teller, E. 1957. Equations of state calculations by fast computing machines. J. Chem. Phys., 21, 1087–1092.
Meyer, C. D. 2000. Matrix Analysis and Applied Linear Algebra. SIAM.
Moore, G., and Read, N. 1991. Nonabelions in the fractional quantum Hall effect. Nucl. Phys. B, 360, 362–396.
Moroni, S., and Baroni, S. 1999. Reptation quantum Monte Carlo: A method for unbiased ground-state averages and imaginary-time correlations. Phys. Rev. Lett., 82, 4745–4748.
Moskowitz, J. W., Schmidt, K. E., Lee, M. A., and Kalos, M. H. 1982. A new look at correlation energy in atomic and molecular systems. II: the application of the Green’s function Monte Carlo method to LiH. J. Chem. Phys., 77, 349–355.
Mott, N. F. 1949. The basis of the electron theory of metals, with special reference to the transition metals. Proc. Phys. Soc. (London), 62, 416–422.
Mott, N. F. 1990. Metal-Insulator Transitions. Taylor and Francis.
Nagaoka, Y. 1966. Ferromagnetism in a narrow, almost half-filled s band. Phys. Rev., 147, 392–405.
Needs, R. J., Towler, M. D., Drummond, N. D., and Lopez R ´ ´ıos, P. 2010. Continuum variational and diffusion quantum Monte Carlo calculations. J. Phys.: Condens. Matter, 22, 023201.
Nelson, E. 1966. Derivation of the Schrodinger equation from Newtonian mechanics.¨ Phys. Rev., 150, 1079–1085.
Neuscamman, E. 2012. Size consistency error in the antisymmetric geminal power wave function can be completely removed. Phys. Rev. Lett., 109, 203001.
Neuscamman, E., Umrigar, C. J., and Chan, G.K.-L. 2012. Optimizing large parameter sets in variational quantum Monte Carlo. Phys. Rev. B, 85, 045103.
Nightingale, M. P., and Melik-Alaverdian, V. 2001. Optimization of ground- and excitedstate wave functions and van der Waals clusters. Phys. Rev. Lett., 87, 043401.
Norris, J. R. 1997. Markov Chains. Cambridge Series in Statistical and Probabilistic Mathematics.
Nozieres, P. 1964. Theory of Interacting Fermi Systems. New York: W.A. Benjamin.
Ortiz, G., Ceperley, D. M., and Martin, R. M. 1993. New stochastic method for systems with broken time-reversal symmetry: 2D fermions in a magnetic field. Phys. Rev. Lett., 71, 2777–2780.
Oshikawa, M., and Senthil, T. 2006. Fractionalization, topological order, and quasiparticle statistics. Phys. Rev. Lett., 96, 060601.
Pandharipande, V. R., and Itoh, N. 1973. Effective mass of 3He in liquid 4He. Phys. Rev. A, 8, 2564–2566.
Parisi, G. 1984. Prolegomena to any future computer evaluation of the QCD mass spectrum. Pages 531–541 of: Progress in Gauge Field Theory. NATO ASI Series, vol. 115. Springer.
Parisi, G., and Wu, Y.-S. 1981. Perturbation theory without gauge fixing. Sci. Sinica, 24, 483–496.
Pauling, L. 1960. The Nature of the Chemical Bond. Cornell University Press.
Pierleoni, C., and Ceperley, D. M. 2005. Computational methods in coupled electronion Monte Carlo simulations. ChemPhysChem, 6, 1872–1878.
Pitaevskii, L., and Stringari, S. 1991. Uncertainly principle, quantum fluctuations, and broken symmetries. J. Low Temp. Phys., 85, 377–388.
Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. 2007. Numerical Recipes 3rd edition: The Art of Scientific Computing. Cambridge University Press.
Read, N., and Chakraborty, B. 1989. Statistics of the excitations of the resonating-valencebond state. Phys. Rev. B, 40, 7133–7140.
Read, N., and Rezayi, E. 1999. Beyond paired quantum Hall states: Parafermions and incompressible states in the first excited Landau level. Phys. Rev. B, 59, 8084–8092.
Reger, J. D., and Young, A. P. 1988. Monte Carlo simulations of the spin-1/2 Heisenberg antiferromagnet on a square lattice. Phys. Rev. B, 37, 5978–5981.
Reynolds, P. J., Ceperley, D. M., Alder, B. J., and Lester, Jr. 1982. Fixednode quantum Monte Carlo for molecules. J. Chem. Phys., 77, 5593–5603.
Ring, P., and Schuck, P. 2004. The Nuclear Many-body Problem. Springer-Verlag.
Rokhsar, D. S., and Kotliar, B. G. 1991. Gutzwiller projection for bosons. Phys. Rev. B, 44, 10328–10332.
Schmidt, K. E., and Pandharipande, V. R. 1979. New variational wave function for liquid 3He. Phys. Rev. B, 19, 2504–2519.
Schollwock, U. 2005. The density-matrix renormalization group.¨ Rev. Mod. Phys., 77, 259–315.
Schollwock, U. 2011. The density-matrix renormalization group in the age of matrix¨ product states. Ann. Phys., 326, 96–192.
Schrieffer, J. R. 1964. Theory of Superconductivity. New York: W.A. Benjamin.
Senthil, T., and Fisher, M. P. A. 2000. Z2 gauge theory of electron fractionalization in strongly correlated systems. Phys. Rev. B, 62, 7850–7881.
Shastry, B. S. 1988. Exact solution of an S = 1/2 Heisenberg antiferromagnetic chain with long-ranged interactions. Phys. Rev. Lett., 60, 639–642.
Slater, J. C. 1930. The electronic structure of metals. Rev. Mod. Phys., 6, 209–280.
Sorella, S. 1998. Green function Monte Carlo with stochastic reconfiguration. Phys. Rev. Lett., 80, 4558–4561.
Sorella, S. 2001. Generalized Lanczos algorithm for variational quantum Monte Carlo. Phys. Rev. B, 64, 024512.
Sorella, S. 2002. Effective Hamiltonian approach for strongly correlated lattice models. arXiv:cond-mat/0201388.
Sorella, S. 2005. Wave function optimization in the variational Monte Carlo method. Phys. Rev. B, 71, 241103.
Sorella, S., Baroni, S., Car, R., and Parrinello, M. 1989. A novel technique for the simulation of interacting fermion systems. Europhys. Lett., 8, 663–668.
Sorella, S., Martins, G. B., Becca, F., Gazza, C., Capriotti, L., Parola, A., and Dagotto, E. 2002. Superconductivity in the two-dimensional t−J model. Phys. Rev. Lett., 88, 117002.
Sorella, S., Casula, M., and Rocca, D. 2007. Weak binding between two aromatic rings: Feeling the van der Waals attraction by quantum Monte Carlo methods. J. Chem. Phys., 127, 014105.
Sorella, S., Devaux, N., Dagrada, M., Mazzola, G., and Casula, M. 2015. Geminal embedding scheme for optimal atomic basis set construction in correlated calculations. J. Chem. Phys., 143, 244112.
Stewart, G. R. 1984. Heavy-fermion systems. Rev. Mod. Phys., 56, 755–787.
Stratonovich, R. L. 1957. A method for the computation of quantum distribution functions. Doklady Akad. Nauk S.S.S.R., 115, 1097–1100.
Sutherland, B. 1971. Exact results for a quantum many-body problem in one dimension. Phys. Rev. A, 4, 2019–2021.
Sutherland, B. 1975. Model for a multicomponent quantum system. Phys. Rev. B, 12, 3795–3805.
Suzuki, M. 1976a. Generalized Trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems. Commun. Math. Phys., 51, 183–190.
Suzuki, M. 1976b. Relationship between d-dimensional quantal spin systems and (d + 1)- dimensional Ising systems. Prog. Theor. Phys., 56, 1454–1469.
Szabo, A., and Ostlund, N. S. 1996. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory. Dover Publications Inc.
Tasaki, H. 1998. From Nagaoka’s ferromagnetism to flat-band ferromagnetism and beyond. Prog. Theor. Phys., 99, 489–548.
ten Haaf, D. F. B., van Bemmel, H. J. M., van Leeuwen, J. M. J., van Saarloos, W., and Ceperley, D. M. 1995. Proof of upper bound in fixed-node Monte Carlo for lattice fermions. Phys. Rev. B, 51, 13039–13045.
Tocchio, L. F., Becca, F., Parola, A., and Sorella, S. 2008. Role of backflow correlations for the nonmagnetic phase of the t−t Hubbard model. Phys. Rev. B, 78, 041101.
Tocchio, L. F., Becca, F., and Gros, C. 2011. Backflow correlations in the Hubbard model: An efficient tool for the study of the metal-insulator transition and the large-U limit. Phys. Rev. B, 83, 195138.
Toulouse, J., and Umrigar, C. J. 2007. Optimization of quantum Monte Carlo wave functions by energy minimization. J. Chem. Phys., 126, 084102.
Trivedi, N., and Ceperley, D. M. 1989. Green-function Monte Carlo study of quantum antiferromagnets. Phys. Rev. B, 40, 2737–2740.
Trivedi, N., and Ceperley, D. M. 1990. Ground-state correlations of quantum antiferromagnets: a Green-function Monte Carlo study. Phys. Rev. B, 41, 4552–4569.
Trotter, H. F. 1959. On the product of semi-groups of operators. Proc. Am. Math. Soc., 10, 545–551.
Tuckerman, M. E. 2010. Statistical Mechanics: Theory and Molecular Simulation. Oxford University Press.
Umrigar, C. J., and Filippi, C. 2005. Energy and variance optimization of many-body wave functions. Phys. Rev. Lett., 94, 150201.
Umrigar, C. J., Wilson, K. G., and Wilkins, J. W. 1988. Optimized trial wave functions for quantum Monte Carlo calculations. Phys. Rev. Lett., 60, 1719–1722.
Umrigar, C. J., Nightingale, M. P., and Runge, K. J. 1993. A diffusion Monte Carlo algorithm with very small timestep errors. J. Chem. Phys., 99, 2865–2890.
Umrigar, C. J., Toulouse, J., Filippi, C., Sorella, S., and Hennig, R. G. 2007. Alleviation of the fermion-sign problem by optimization of many-body wave functions. Phys. Rev. Lett., 98, 110201.
Wagner, L. K., Bajdich, M., and Mitas, L. 2009. QWalk: A quantum Monte Carlo program for electronic structure. J. Comp. Phys., 228, 3390–3404.
Wen, X.-G. 1991. Topological orders and Chern-Simons theory in strongly correlated quantum liquid. Int. J. Mod. Phys. B, 5, 1641–1648.
Wen, X.-G., and Niu, Q. 1990. Ground-state degeneracy of the fractional quantum Hall states in the presence of a random potential and on high-genus Riemann surfaces. Phys. Rev. B, 41, 9377–9396.
White, S. R. 1992. Density matrix formulation for quantum renormalization groups. Phys. Rev. Lett., 69, 2863–2866.
White, S. R, and Feiguin, A. 2004. Real-time evolution using the density matrix renormalization group. Phys. Rev. Lett., 93, 076401.
White, S. R., and Scalapino, D.J. 1998. Density matrix renormalization group study of the striped phase in the 2D t−J model. Phys. Rev. Lett., 80, 1272–1275.
Wigner, E., and Seitz, F. 1934. On the constitution of metallic Sodium. II. Phys. Rev., 46, 509–524.
Yokoyama, H., and Shiba, H. 1987a. Variational Monte Carlo studies of Hubbard model. I. J. Phys. Soc. Jpn., 56, 1490–1506.
Yokoyama, H., and Shiba, H. 1987b. Variational Monte Carlo studies of Hubbard model. II. J. Phys. Soc. Jpn., 56, 3582–3592.
Yokoyama, H., and Shiba, H. 1990. Variational Monte Carlo studies of Hubbard model. III. Intersite correlation effects. J. Phys. Soc. Jpn., 59, 3669–3686.
Young, P. 2012. Everything you wanted to know about data analysis and fitting but were afraid to ask. arXiv:1210.3781.
Zhang, F. C., and Rice, T. M. 1988. Effective Hamiltonian for the superconducting Cu oxides. Phys. Rev. B, 37, 3759–3761.
Zhang, S., Carlson, J., and Gubernatis, J. E. 1995. Constrained path quantum Monte Carlo method for fermion ground states. Phys. Rev. Lett., 74, 3652–3655.
Acceptance-rejection method, 62 Anti-symmetric geminal power Constrained optimization, 246 Definition, 237 Area law, 31 Atomic orbitals, 238 Contracted, 241 Gaussian type, 240 Periodic systems, 241 Slater type, 240 Twisted boundary conditions, 241 Auxiliary-field quantum Monte Carlo Attractive Hubbard model, 228 Backward propagation, 221 Forward propagation, 221 Pseudo-partition function, 220 Repulsive Hubbard model, 229 Sequential updates, 223 Sign problem, 221 Backflow correlations, 21, 122 Bayes formula, 43 BCS Hamiltonian, 27, 115, 125 BCS wave function, 26, 125 Block analysis, 83 Boltzmann distribution, 86, 90 Bootstrap method, 78 Bosonic wave function Condensed state, 110 Fast update of the Jastrow factor, 111 Permanent, 111 Box-Muller trick, 65 Brownian motion, 86 Central limit theorem, 53 Changing random variables, 47 Characteristic function, 45 Chebyshev’s inequality, 48 Cholesky decomposition, 222 Classical mapping, 33