Flowers for Tuesday

Computational Physics

HF Method ＆ Post HF

下载 Markdown

参考教材

[1] Jos Thijssen , Computational Physics

[2] Attila Szabo, Neil S. Ostlund , Modern Quantum Chemistry

那么，开始吧

1. Hartree Fock方法

HF方法算是很基础的近似方法，用正常方式讲就有点太无聊了，所以不妨直接从二次量子化来写。

1.1 二次量子化

我们知道单体波函数可以写成能量本征态基底{φi}\{ \varphi_i \}的展开

φ=∑cnφn\varphi = \sum c_n \varphi_n

这意味着单体波函数的算符O^\hat{O}（例如动能算符 ∑p^a22m\sum \frac{\hat{p}_a^2}{2m}、外势能算符 ∑V(x^a)\sum V(\hat{x}_a)）其实可以根据本征态进行矩阵元展开

oij=⟨φi∣o^∣φj⟩=∫φi∗(x)o^φj(x) dxo_{ij} = \langle \varphi_i \vert{} \hat{o} \vert{} \varphi_j \rangle = \int \varphi_i^*(x) \hat{o} \varphi_j(x) \, dx

从这个角度看，算符 o^\hat{o} 作用在 ∣φj⟩\vert{}\varphi_j\rangle 上可依据完备性条件展开

o^∣φj⟩=∑i∣φi⟩⟨φi∣o^∣φj⟩=∑ioij∣φi⟩\hat{o} \vert{}\varphi_j\rangle = \sum_i \vert{}\varphi_i\rangle \langle \varphi_i \vert{} \hat{o} \vert{} \varphi_j \rangle = \sum_i o_{ij} \vert{}\varphi_i\rangle

这件事有一个等价的理解角度，也就是相当于将状态 jj 湮灭，并以振幅 oijo_{ij} 产生状态 ii。因此用二次量子化的写法就是

O^=∑i,joijai†aj\hat{O} = \sum_{i, j} o_{ij} a_i^\dagger a_j

所以说，二次量子化就是把线性代数这件事做得更加彻底一点，把波函数彻底变成了基底下的一串矢量。

另一件有趣的事情是，多体波函数的希尔伯特空间是可以被多组完备的单体波函数基底所组合展开的。

例如一个二元函数f(x,y)f(x,y)，我们有两组完备基底{ui(x)}\{u_i(x)\}和{vi(y)}\{v_i(y)\}，那么有

f(x,y)=∑cijui(x)vj(y)f(x,y)=\sum c_{ij}u_i(x)v_j(y)

证明就是先把某个变量视作常数，再完备展开即可。这件事就可以自然拓展到多体波函数，例如两个粒子的波函数就可以用各自的轨道进行展开

Ψ(r1,r2)=∑n=1∞∑m=1∞cnmϕn(r1)ϕm(r2)\Psi(\mathbf{r}_1, \mathbf{r}_2) = \sum_{n=1}^\infty \sum_{m=1}^\infty c_{nm} \phi_n(\mathbf{r}_1) \phi_m(\mathbf{r}_2)

用更线性代数的语言，就是多体波函数存在于单体波函数张成的直积空间里，对于 NN 粒子系统，直积空间为

HN=H1⊗N=H1⊗H1⊗⋯⊗H1\mathcal{H}_N = \mathcal{H}_1^{\otimes N} = \mathcal{H}_1 \otimes \mathcal{H}_1 \otimes \dots \otimes \mathcal{H}_1

不难发现这个空间对真实波函数来说是扩大的，原因是真实体系存在玻色和费米子的正反对称约束，定义对称化算符 S^+\hat{S}_+ 与反对称化算符 S^−\hat{S}_-

S^±=1N!∑σ∈SN(±1)sgn(σ)P^σ\hat{S}_\pm = \frac{1}{N!} \sum_{\sigma \in S_N} (\pm 1)^{\text{sgn}(\sigma)} \hat{P}_\sigma

其中 sgn(σ)\text{sgn}(\sigma) 为置换的奇偶性（奇置换为 11，偶置换为 00）。玻色子空间 HN+\mathcal{H}_N^+ 和费米子空间 HN−\mathcal{H}_N^- 分别为

HN±=S^±(H1⊗N)\mathcal{H}_N^\pm = \hat{S}_\pm (\mathcal{H}_1^{\otimes N})

例如单纯从数学上看，无相互作用体系下，单体波函数之间的简单乘积就是多体波函数的本征态，但是满足反对称的只有Slater行列式的形式

ψ(x1,…,xN)=1N!det⁡(φk1(x1)…φkN(x1)⋮⋱⋮φk1(xN)…φkN(xN))\psi(x_1, \dots, x_N) = \frac{1}{\sqrt{N!}} \det \begin{pmatrix} \varphi_{k_1}(x_1) & \dots & \varphi_{k_N}(x_1) \\ \vdots & \ddots & \vdots \\ \varphi_{k_1}(x_N) & \dots & \varphi_{k_N}(x_N) \end{pmatrix}

二次量子化定义 Fock 空间 F±\mathcal{F}_\pm 为所有 NN 粒子子空间的直和

F±=⨁N=0∞HN±=H0⊕H1⊕H2±⊕H3±⊕…\mathcal{F}_\pm = \bigoplus_{N=0}^{\infty} \mathcal{H}_N^\pm = \mathcal{H}_0 \oplus \mathcal{H}_1 \oplus \mathcal{H}_2^\pm \oplus \mathcal{H}_3^\pm \oplus \dots

于是可以说，满足反对称性下的NN粒子的直积空间的基底就是任选NN个轨道组成的Slater行列式

∣Φ0⟩=a1†a2†…aN†∣0⟩\vert{}\Phi_0\rangle = a_1^\dagger a_2^\dagger \dots a_N^\dagger \vert{}0\rangle

如果对单粒子轨道做一个截断，例如认为上限是MM个轨道，那么对应的NN粒子波函数的基底数量自然就是(MN)=M!N!(M−N)!\binom{M}{N} = \frac{M!}{N!(M-N)!}个。

不难发现如果存在相互作用，Slater行列式就不再是多体波函数的本征态，但可证明它们仍然是多体波函数的完备基底，这是之后一切讨论的基础。

1.2 二次量子化下的HF方法

HF方法的本质，是选取一组最佳的单粒子轨道{φi(r)}\{\varphi_i(\mathbf{r})\}，使得真实多体波函数在其基态Slater行列式对应的基底上的投影最大。

而根据量子力学的变分原理，这件事的等价表述就是使得基态Slater行列式对应的基底能量值最小。于是我们来试着求一下这个能量关于轨道的泛函表达式。

假设分子系统中有 NN 个电子，我们在单粒子空间中选择一组正交归一的已知基函数 {χp(r)}p=1M\{\chi_p(\mathbf{r})\}_{p=1}^M。利用电子的产生算符 ap†a_p^\dagger 和湮灭算符 aqa_q，系统的哈密顿量可写为

H^=∑p,qhpqap†aq+12∑p,q,r,sgpqrsap†aq†asar\hat{H} = \sum_{p,q} h_{pq} a_p^\dagger a_q + \frac{1}{2} \sum_{p,q,r,s} g_{pqrs} a_p^\dagger a_q^\dagger a_s a_r

其中

单电子积分 hpqh_{pq}（动能 + 原子核吸引势）

hpq=∫χp∗(r)(−12∇2−∑AZA∣r−RA∣)χq(r) d3rh_{pq} = \int \chi_p^*(\mathbf{r}) \left( -\frac{1}{2}\nabla^2 - \sum_A \frac{Z_A}{\vert{}\mathbf{r} - \mathbf{R}_A\vert{}} \right) \chi_q(\mathbf{r}) \, d^3r

双电子库仑积分 gpqrsg_{pqrs}

gpqrs=∬χp∗(r1)χq∗(r2)χr(r1)χs(r2)∣r1−r2∣ d3r1d3r2g_{pqrs} = \iint \frac{\chi_p^*(\mathbf{r}_1) \chi_q^*(\mathbf{r}_2) \chi_r(\mathbf{r}_1) \chi_s(\mathbf{r}_2)}{\vert{}\mathbf{r}_1 - \mathbf{r}_2\vert{}} \, d^3r_1 d^3r_2

然后来计算基态基底对应的能量，E0=⟨Φ0∣H^∣Φ0⟩E_0 = \langle \Phi_0 \vert{} \hat{H} \vert{} \Phi_0 \rangle，这里不加推导地给出 Hartree-Fock 能量表达式

EHF=∑i=1Nhii+12∑i,j=1N(gijij−gijji)E_{\text{HF}} = \sum_{i=1}^N h_{ii} + \frac{1}{2} \sum_{i,j=1}^N \left( g_{ijij} - g_{ijji} \right)

为了找到使 EHFE_{\text{HF}} 最小的最佳轨道，我们需要对正交轨道基底做无穷小旋转变换，这等价于

ap†→ap†+∑qϵpqaq†a_p^\dagger \to a_p^\dagger + \sum_q \epsilon_{pq} a_q^\dagger

对算符变化的变分为 0 (δEHF=0\delta E_{\text{HF}} = 0)，这一操作直接导出Fock 算符

f^=∑p,qfpqap†aq\hat{f} = \sum_{p,q} f_{pq} a_p^\dagger a_q

其中 Fock 矩阵元 fpqf_{pq} 表达为

fpq=hpq+∑i=1N(gpiqi−gpiiq)f_{pq} = h_{pq} + \sum_{i=1}^N \left( g_{piqi} - g_{piiq} \right)

对 Fock 矩阵对角化，就得到矩阵形式的 Hartree-Fock 方程

FC=SCE\mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C} \mathbf{\mathbf{\mathcal{E}}}

通过自洽场方法迭代求解即可。

不过由于含相互作用的多体波函数无法被Slater行列式完备表达，所以这个方法注定只能是一个最基础的近似。后续的一系列Post HF方法会以Φ0\Phi_0作为0级近似进一步地优化。

2. Post HF方法

2.1 构型相互作用方法（Configuration Interaction, CI）

多体希尔伯特空间中的一个任意向量形式可以用不同NN粒子轨道的斯莱特行列式展开（二次量子化下对应的就是不同的激发态）

∣Ψ⟩=C0∣Φ0⟩+∑aiCia∣Φia⟩+∑abijCijab∣Φijab⟩+…\vert{}\Psi\rangle = C_0 \vert{}\Phi_0\rangle + \sum_{ai} C_i^a \vert{}\Phi_i^a\rangle + \sum_{abij} C_{ij}^{ab} \vert{}\Phi_{ij}^{ab}\rangle + \dots

∣Φ0⟩\vert{}\Phi_0\rangle为参考态（0阶激发），∣Φia⟩=aa†ai∣Φ0⟩\vert{}\Phi_i^a\rangle = a_a^\dagger a_i \vert{}\Phi_0\rangle为单激发行列式（Singles, S），∣Φijab⟩=aa†ab†ajai∣Φ0⟩\vert{}\Phi_{ij}^{ab}\rangle = a_a^\dagger a_b^\dagger a_j a_i \vert{}\Phi_0\rangle为双激发行列式（Doubles, D），以此类推。

所以一个优化HF方法的最基本想法就是放弃单行列式的限制，直接将波函数展开为全希尔伯特空间中所有正交斯莱特行列式的线性叠加。

当然由于正交轨道数理论上是无限的，所以还是得先对轨道数做一个截断，假设上限为MM。

根据线性变分原理，我们希望寻找一组展开系数 C={C0,Cia,Cijab,… }T\mathbf{C} = \{C_0, C_i^a, C_{ij}^{ab}, \dots\}^T，使得体系的能量期望值极小化

ECI=⟨ΨCI∣H^∣ΨCI⟩⟨ΨCI∣ΨCI⟩E_{\text{CI}} = \frac{\langle \Psi_{\text{CI}} \vert{} \hat{H} \vert{} \Psi_{\text{CI}} \rangle}{\langle \Psi_{\text{CI}} \vert{} \Psi_{\text{CI}} \rangle}

利用拉格朗日乘子法 δ(⟨ΨCI∣H^∣ΨCI⟩−E(⟨ΨCI∣ΨCI⟩−1))=0\delta \left( \langle \Psi_{\text{CI}} \vert{} \hat{H} \vert{} \Psi_{\text{CI}} \rangle - E (\langle \Psi_{\text{CI}} \vert{} \Psi_{\text{CI}} \rangle - 1) \right) = 0，直接导出 CI 矩阵特征值方程

HC=EC\mathbf{H} \mathbf{C} = E \mathbf{C}

剩下的数值计算步骤其实就没什么好说的，毕竟轨道都由HF算出来了。（emm理论上跳过HF，直接用已知正交轨道组求CI也是可行的，因为CI其实是一个独立的方法）不过考虑到实际计算量，肯定不会算所有的激发态，一般取到二阶激发即可。

2.2 耦合簇理论 （Coupled-Cluster, CC）

与CI这种很容易get到思路的方法不同，CC让人不太能搞懂它是怎么被想出来的

其核心想法是相似变换，把多体波函数 ∣ΨCC⟩\vert{}\Psi_{\text{CC}}\rangle 表达为参考态 ∣Φ0⟩\vert{}\Phi_0\rangle 上的指数激发算符

∣ΨCC⟩=eT^∣Φ0⟩\vert{}\Psi_{\text{CC}}\rangle = e^{\hat{T}} \vert{}\Phi_0\rangle

于是

H^eT^∣Φ0⟩=ECCeT^∣Φ0⟩\hat{H} e^{\hat{T}} \vert{}\Phi_0\rangle = E_{\text{CC}} e^{\hat{T}} \vert{}\Phi_0\rangle
Hˉ∣Φ0⟩=ECC∣Φ0⟩,其中 Hˉ≡e−T^H^eT^\bar{H} \vert{}\Phi_0\rangle = E_{\text{CC}} \vert{}\Phi_0\rangle, \quad \text{其中 } \bar{H} \equiv e^{-\hat{T}} \hat{H} e^{\hat{T}}

其中簇算符T^\hat{T} 可以展开为激发算符的线性组合：T^=T^1+T^2+T^3+⋯+T^N\hat{T} = \hat{T}_1 + \hat{T}_2 + \hat{T}_3 + \dots + \hat{T}_N

例如单激发算符 T^1\hat{T}_1

T^1=∑i,atia aa†ai\hat{T}_1 = \sum_{i,a} t_i^a \, a_a^\dagger a_i

双激发算符 T^2\hat{T}_2 （考虑到重复计算）

T^2=14∑i,j,a,btijab aa†ab†ajai\hat{T}_2 = \frac{1}{4} \sum_{i,j,a,b} t_{ij}^{ab} \, a_a^\dagger a_b^\dagger a_j a_i

其中的簇振幅 tia,tijabt_i^a, t_{ij}^{ab} 就是需要求解的参数。

为什么这么展开是可行的呢，似乎激发算符族TiT_i并不是一个完备的算符展开基底。这里更严格的说法是，所有从Φ0→Φ\Phi_0 \rightarrow \Phi的变换算符，都可以用激发算符族TiT_i展开，这件事是显然的。（CI本质上就在告诉我们这件事）

如果T^\hat{T}不在{Ti}\{T_i\}张成的空间里，eT^e^{\hat{T}}显然也不会，所以可证。

那么，利用泰勒级数将指数算符展开

eT^=I^+T^+12!T^2+13!T^3+…e^{\hat{T}} = \hat{I} + \hat{T} + \frac{1}{2!}\hat{T}^2 + \frac{1}{3!}\hat{T}^3 + \dots

如果考虑最常用的 CCSD 方法（截断到 T^=T^1+T^2\hat{T} = \hat{T}_1 + \hat{T}_2），将展开式写出来

eT^1+T^2=I^+T^1+(T^2+12T^12)+(T^1T^2+16T^13)+(12T^22+12T^12T^2+124T^14)+…e^{\hat{T}_1 + \hat{T}_2} = \hat{I} + \hat{T}_1 + \left( \hat{T}_2 + \frac{1}{2}\hat{T}_1^2 \right) + \left( \hat{T}_1 \hat{T}_2 + \frac{1}{6}\hat{T}_1^3 \right) + \left( \frac{1}{2}\hat{T}_2^2 + \frac{1}{2}\hat{T}_1^2 \hat{T}_2 + \frac{1}{24}\hat{T}_1^4 \right) + \dots

会发现这依然是一种扩展基底的方式，但如果直接将 eT^∣Φ0⟩e^{\hat{T}}\vert{}\Phi_0\rangle 代入变分原理会导致多项式级数不终止，还不如直接用CI。

CC采用的是非变分的投影法，根据 BCH 公式，相似变换后的哈密顿量 Hˉ\bar{H} 可以展开为嵌套对易子

Hˉ=H^+[H^,T^]+12![[H^,T^],T^]+13![[[H^,T^],T^],T^]+14![[[[H^,T^],T^],T^],T^]\bar{H} = \hat{H} + [\hat{H}, \hat{T}] + \frac{1}{2!} [[\hat{H}, \hat{T}], \hat{T}] + \frac{1}{3!} [[[\hat{H}, \hat{T}], \hat{T}], \hat{T}] + \frac{1}{4!} [[[[\hat{H}, \hat{T}], \hat{T}], \hat{T}], \hat{T}]

因为哈密顿量 H^\hat{H} 最多只包含双体算符（4 个费米子生成湮灭算符），嵌套对易子级数会在第 4 阶严格截断。所以某种意义上可以认为CC是针对特定哈密顿量形式的特定方法。

于是现在我们有了代数闭合的方程

Hˉ∣Φ0⟩=E∣Φ0⟩\bar{H} \vert{}\Phi_0\rangle = E \vert{}\Phi_0\rangle

为了解出未知量（T^\hat{T} 中的振幅系数 tia,tijab,…t_i^a, t_{ij}^{ab}, \dots）以及能量 EE，可以直接对基底 {⟨Φ0∣,⟨Φia∣,⟨Φijab∣,… }\{\langle \Phi_0\vert{}, \langle \Phi_i^a\vert{}, \langle \Phi_{ij}^{ab}\vert{}, \dots\} 进行左投影。

其中对参考态的投影得到能量方程（由于 H^\hat{H} 最多是两体算符，超过双激发的对易项在与 ⟨Φ0∣\langle \Phi_0\vert{} 叠合时矩阵元为 0）

⟨Φ0∣Hˉ∣Φ0⟩=E⟨Φ0∣Φ0⟩\langle \Phi_0 \vert{} \bar{H} \vert{}\Phi_0\rangle = E \langle \Phi_0 \vert{} \Phi_0 \rangle
ECC=⟨Φ0∣Hˉ∣Φ0⟩=⟨Φ0∣H^+[H^,T^]+12![[H^,T^],T^]∣Φ0⟩E_{\text{CC}} = \langle \Phi_0 \vert{} \bar{H} \vert{}\Phi_0\rangle = \langle \Phi_0 \vert{} \hat{H} + [\hat{H}, \hat{T}] + \frac{1}{2!} [[\hat{H}, \hat{T}], \hat{T}] \vert{}\Phi_0\rangle

对其余激发态投影可得

⟨ΦK∣Hˉ∣Φ0⟩=E⟨ΦK∣Φ0⟩\langle \Phi_K \vert{} \bar{H} \vert{}\Phi_0\rangle = E \langle \Phi_K \vert{} \Phi_0 \rangle

于是有严格的振幅方程组

{单激发方程 (S):⟨Φia∣H^+[H^,T^]+12![[H^,T^],T^]+13![[[H^,T^],T^],T^]∣Φ0⟩=0双激发方程 (D):⟨Φijab∣H^+[H^,T^]+12![[H^,T^],T^]+13![[[H^,T^],T^],T^]+14![[[[H^,T^],T^],T^],T^]∣Φ0⟩=0高阶激发方程:⟨Φijk…abc…∣Hˉ∣Φ0⟩=0\begin{cases} \mathbf{\text{单激发方程 (S):}} & \langle \Phi_i^a \vert{} \hat{H} + [\hat{H}, \hat{T}] + \frac{1}{2!} [[\hat{H}, \hat{T}], \hat{T}] + \frac{1}{3!} [[[\hat{H}, \hat{T}], \hat{T}], \hat{T}] \vert{}\Phi_0\rangle = 0 \\ \mathbf{\text{双激发方程 (D):}} & \langle \Phi_{ij}^{ab} \vert{} \hat{H} + [\hat{H}, \hat{T}] + \frac{1}{2!} [[\hat{H}, \hat{T}], \hat{T}] + \frac{1}{3!} [[[\hat{H}, \hat{T}], \hat{T}], \hat{T}] + \frac{1}{4!} [[[[\hat{H}, \hat{T}], \hat{T}], \hat{T}], \hat{T}] \vert{}\Phi_0\rangle = 0 \\ \mathbf{\text{高阶激发方程:}} & \langle \Phi_{ijk\dots}^{abc\dots} \vert{} \bar{H} \vert{}\Phi_0\rangle = 0 \end{cases}

CCSD截断就是强行让高阶项的振幅为0，然后根据残余的方程组解出低阶项的振幅系数。最终得到簇算符形式。

一个有趣的事情是，我们的方程根本没有涉及到变分，这意味着我们只要求能量为定值，也就是对应的可以是任何一个本征态而不一定是基态。事实也的确如此，CC的方程的解是不唯一的，一般都是迭代求解出一组特解。但很容易猜到，弱相关体系下，我们的基态特解一定很接近解空间的零点，也就是靠近参考态Φ0\Phi_0

3. 与DFT的一点对比

下面是个人的一些思考，我们知道HF和DFT都采用了单Slater行列式的形式来表达多体波函数。但二者的思想其实有本质性的差异。

HF方法正如前面讨论的，是真实多体波函数在限定形式下的最佳近似，所以其实质就是投影，并没有改变真实的物理环境。

而DFT是假设一个理论上的势场加一个无相关的体系，其基态波函数的密度能恰好与真实体系的基态波函数密度一致。无相关体系下单Slater行列式无疑是完备的，所以DFT的一切问题都转移到了这个理论势场的求解，进而转移到交换相关泛函形式的求解。

当然一个不可回避的问题是，我们怎么知道这个势场是存在的呢？

HK定理虽然确保了映射的唯一性，但似乎没有确保其存在，也就是vv的可构性未知。别急，这件事其实是有人证明的，详见Lieb的Density functionals for coulomb systems

HF Method ＆ Post HF · 图 1

啊没错就是严格求解了一维 Hubbard 模型的那位Lieb。

继续阅读 · 计算物理

密度矩阵重整化群 / Density Matrix Renormalization Group
← 返回笔记目录回到顶部 ↑