最近在写CC。因为CC公式大多比较长，而且推起来很麻烦，所以我就打算整理一下公式再写几篇博客，以防之后忘记。
Coupled Cluster（耦合簇，CC）理论是一种量子化学计算方法，它能够恰当地考虑电子的动态相关，从而能够高精度地计算一些体系的相关能。CC 最早被核物理学家所提出，主要用来计算核反应体系的一些性质1
作为本系列博客的第一篇文章，本文将会介绍CC形式理论，并会给出CC能量与振幅所满足的一般方程。在阅读本系列文章之前，建议读者至少要有一定的量子力学基础。
1. 组态波函数的二次量子化表示 对于一特定的 n 电子波函数，可用符号 ∣ p , q , ⋯ r , s ⟩ \ket{p,q,\cdots r,s} ∣ p , q , ⋯ r , s ⟩ p p p q q q
∣ Ψ ⟩ = ∑ C p q ⋯ r s ∣ p , q , ⋯ r , s ⟩ \begin{equation} \ket{\Psi} = \sum C_{pq\cdots rs}\ket{p,q,\cdots r,s} \end{equation} ∣ Ψ ⟩ = ∑ C pq ⋯ rs ∣ p , q , ⋯ r , s ⟩ 如果我们选定一个态作为参考态（如 HF 基态），把其他态看作参考态的激发，上式就可以写成（在本文中，i j k l ijkl ijk l a b c d abcd ab c d p q r s pqrs pq rs
∣ Ψ ⟩ = ∣ Ψ 0 ⟩ + ∑ i a c i a ∣ Ψ i a ⟩ + ∑ i a j b c i j a b ∣ Ψ i j a b ⟩ + ⋯ \begin{equation} \ket{\Psi} = \ket{\Psi_0}+\sum_{ia}c_i^a\ket{\Psi_i^a}+\sum_{iajb}c_{ij}^{ab}\ket{\Psi_{ij}^{ab}}+\cdots \end{equation} ∣ Ψ ⟩ = ∣ Ψ 0 ⟩ + ia ∑ c i a ∣ Ψ i a ⟩ + iajb ∑ c ij ab ∣ Ψ ij ab ⟩ + ⋯ 其中 ∣ Ψ 0 ⟩ \ket{\Psi_0} ∣ Ψ 0 ⟩ ∣ Ψ i a ⟩ \ket{\Psi_i^a} ∣ Ψ i a ⟩ i i i a a a CI 展开 。
1.2 组态波函数的二次量子化表示 写出组态波函数二次量子化表示
∣ p , q , ⋯ r , s ⟩ = a ^ s † a ^ r † ⋯ a ^ q † a ^ p † ∣ ⟩ \begin{equation} \ket{p,q,\cdots r,s} = \hat{a}^\dagger_s\hat{a}^\dagger_r\cdots\hat{a}^\dagger_q\hat{a}^\dagger_p\ket{} \end{equation} ∣ p , q , ⋯ r , s ⟩ = a ^ s † a ^ r † ⋯ a ^ q † a ^ p † ∣ ⟩ 式中 a ^ † \hat{a}^\dagger a ^ † a ^ \hat{a} a ^
[ a ^ p † , a ^ q † ] + = a ^ p † a ^ q † + a ^ q † a ^ p † = 0 [ a ^ p , a ^ q ] + = a ^ p a ^ + a ^ q a ^ p = 0 [ a ^ p † , a ^ q ] + = a ^ p † a ^ q + a ^ q a ^ p † = δ p q \begin{equation} \begin{aligned} [\hat{a}^\dagger_p,\hat{a}^\dagger_q]_+ &= \hat{a}^\dagger_p\hat{a}^\dagger_q+\hat{a}^\dagger_q\hat{a}^\dagger_p=0\\ [\hat{a}_p,\hat{a}_q]_+ &= \hat{a}_p\hat{a}+\hat{a}_q\hat{a}_p=0\\ [\hat{a}^\dagger_p,\hat{a}_q]_+ &= \hat{a}^\dagger_p\hat{a}_q+\hat{a}_q\hat{a}^\dagger_p=\delta_{pq} \end{aligned} \end{equation} [ a ^ p † , a ^ q † ] + [ a ^ p , a ^ q ] + [ a ^ p † , a ^ q ] + = a ^ p † a ^ q † + a ^ q † a ^ p † = 0 = a ^ p a ^ + a ^ q a ^ p = 0 = a ^ p † a ^ q + a ^ q a ^ p † = δ pq 通过上述对易关系，我们可以验证一下费米子体系的泡利不相容原理
∣ p , q , ⋯ r , s ⟩ = a ^ s † a ^ r † ⋯ a ^ q † a ^ p † ∣ ⟩ = − a ^ s † a ^ r † ⋯ a ^ p † a ^ q † ∣ ⟩ = ∣ q , p , ⋯ r , s ⟩ \begin{equation} \ket{p,q,\cdots r,s} = \hat{a}^\dagger_s\hat{a}^\dagger_r\cdots\hat{a}^\dagger_q\hat{a}^\dagger_p\ket{} = -\hat{a}^\dagger_s\hat{a}^\dagger_r\cdots\hat{a}^\dagger_p\hat{a}^\dagger_q\ket{} = \ket{q,p,\cdots r,s} \end{equation} ∣ p , q , ⋯ r , s ⟩ = a ^ s † a ^ r † ⋯ a ^ q † a ^ p † ∣ ⟩ = − a ^ s † a ^ r † ⋯ a ^ p † a ^ q † ∣ ⟩ = ∣ q , p , ⋯ r , s ⟩ 即交换两电子位置，波函数变号。也可以通过产生湮灭算符的对易关系和真空态的一些性质，来计算一些积分，比如用于计算两个双电子态的重叠积分
⟨ p , q ∣ r , s ⟩ = ⟨ ∣ a ^ p a ^ q a ^ s † a ^ r † ∣ ⟩ = ⟨ ∣ a ^ p ( δ q s − a ^ s † a ^ q ) a ^ r † ∣ ⟩ = δ q s ⟨ ∣ a ^ p a ^ r † ∣ ⟩ − ⟨ ∣ a ^ p a ^ s † a ^ q a ^ r † ∣ ⟩ = δ p r δ q s − δ q r ⟨ ∣ a ^ p a ^ s † ∣ ⟩ = δ p r δ q s − δ p s δ q r \begin{equation} \begin{aligned} \braket{p,q|r,s} &= \bra{}\hat{a}_p\hat{a}_q\hat{a}^\dagger_s\hat{a}^\dagger_r\ket{} = \bra{}\hat{a}_p(\delta_{qs}-\hat{a}^\dagger_s\hat{a}_q)\hat{a}^\dagger_r\ket{}\\ &= \delta_{qs}\bra{}\hat{a}_p\hat{a}^\dagger_r\ket{}-\bra{}\hat{a}_p\hat{a}^\dagger_s\hat{a}_q\hat{a}^\dagger_r\ket{}\\ &= \delta_{pr}\delta_{qs}-\delta_{qr}\bra{}\hat{a}_p\hat{a}^\dagger_s\ket{}\\ &= \delta_{pr}\delta_{qs}-\delta_{ps}\delta_{qr} \end{aligned} \end{equation} ⟨ p , q ∣ r , s ⟩ = ⟨ ∣ a ^ p a ^ q a ^ s † a ^ r † ∣ ⟩ = ⟨ ∣ a ^ p ( δ q s − a ^ s † a ^ q ) a ^ r † ∣ ⟩ = δ q s ⟨ ∣ a ^ p a ^ r † ∣ ⟩ − ⟨ ∣ a ^ p a ^ s † a ^ q a ^ r † ∣ ⟩ = δ p r δ q s − δ q r ⟨ ∣ a ^ p a ^ s † ∣ ⟩ = δ p r δ q s − δ p s δ q r 这就是双电子情况下不同组态之间的正交归一关系。通过二次量子化语言，我们可以在不写出组态波函数的显式表达式的情况下，仅用代数方法便可求得各种积分。但对于实际体系，由于电子数非常多，直接使用（6）式中所示方法就需要处理大量产生湮灭算符的乘积，这将会变得非常困难，有没有方法能够避免这个问题呢？
2. 费米真空与正则序 2.1 费米真空与 q-产生湮灭算符 答案当然是有的。回顾波函数的 CI 展开式（2），我们可以把任意组态看作参考态的激发，即将参考态看作真空态，将其他组态写成产生湮灭算符作用在参考态上的形式，即
∣ Ψ 0 ⟩ = ∣ Ψ 0 ⟩ ∣ Ψ i a ⟩ = a ^ a † a ^ i ∣ Ψ 0 ⟩ ∣ Ψ i j a b ⟩ = a ^ a † a ^ b † a ^ i a ^ j ∣ Ψ 0 ⟩ … \begin{equation} \begin{aligned} \ket{\Psi_0} &= \ket{\Psi_0} \\ \ket{\Psi_i^a} &= \hat{a}^\dagger_a\hat{a}_i\ket{\Psi_0}\\ \ket{\Psi_{ij}^{ab}} &= \hat{a}^\dagger_a\hat{a}^\dagger_b\hat{a}_i\hat{a}_j\ket{\Psi_0}\\ & \ldots \end{aligned} \end{equation} ∣ Ψ 0 ⟩ ∣ Ψ i a ⟩ ∣ Ψ ij ab ⟩ = ∣ Ψ 0 ⟩ = a ^ a † a ^ i ∣ Ψ 0 ⟩ = a ^ a † a ^ b † a ^ i a ^ j ∣ Ψ 0 ⟩ … 通过把真空态与一部分产生算符合并成一个参考态，我们可以极大地减少要处理的产生湮灭算符的个数。如果选取的参考态为 HF 基态，我们就称该参考态为费米真空 。
在将真空态换为费米真空后，产生湮灭算符的一些性质也会发生些许改变。原本湮灭算符作用在真空态上会等于零，即 a ^ p ∣ ⟩ = 0 \hat{a}_p\ket{}=0 a ^ p ∣ ⟩ = 0
a ^ a ∣ Ψ 0 ⟩ = 0 a ^ i † ∣ Ψ 0 ⟩ = 0 \begin{equation} \begin{aligned} \hat{a}_a\ket{\Psi_0}=0\\ \hat{a}^\dagger_i\ket{\Psi_0}=0 \end{aligned} \end{equation} a ^ a ∣ Ψ 0 ⟩ = 0 a ^ i † ∣ Ψ 0 ⟩ = 0 可以看到，在选择费米真空作为真空态后，占据轨道的产生算符会表现出‘湮灭算符‘的性质，占据轨道的湮灭算符会表现出‘产生算符’的性质。因此，在选择费米真空作为真空态后，我们可以将占据轨道的产生算符与非占据轨道的湮灭算符归为一类‘湮灭算符’、将占据轨道的湮灭算符与非占据轨道的产生算符归为一类‘产生算符’。
为了与通常意义下的产生湮灭算符做区别，我们通常会把占据轨道的产生算符与非占据轨道的湮灭算符称为 q q q q q q q q q
[ a ^ p † , a ^ q † ] + = a ^ p † a ^ q † + a ^ q † a ^ p † = 0 [ a ^ p , a ^ q ] + = a ^ p a ^ + a ^ q a ^ p = 0 [ a ^ i † , a ^ j ] + = a ^ i † a ^ j + a ^ j a ^ i † = δ i j [ a ^ a † , a ^ b ] + = a ^ a † a ^ b + a ^ b a ^ a † = δ a b [ a ^ i † , a ^ a ] + = a ^ i † a ^ a + a ^ a a ^ i † = 0 [ a ^ a † , a ^ i ] + = a ^ a † a ^ i + a ^ i a ^ a † = 0 \begin{equation} \begin{aligned} [\hat{a}^\dagger_p,\hat{a}^\dagger_q]_+ &= \hat{a}^\dagger_p\hat{a}^\dagger_q+\hat{a}^\dagger_q\hat{a}^\dagger_p=0\\ [\hat{a}_p,\hat{a}_q]_+ &= \hat{a}_p\hat{a}+\hat{a}_q\hat{a}_p=0\\ [\hat{a}^\dagger_i,\hat{a}_j]_+ &= \hat{a}^\dagger_i\hat{a}_j+\hat{a}_j\hat{a}^\dagger_i=\delta_{ij}\\ [\hat{a}^\dagger_a,\hat{a}_b]_+ &= \hat{a}^\dagger_a\hat{a}_b+\hat{a}_b\hat{a}^\dagger_a=\delta_{ab}\\ [\hat{a}^\dagger_i,\hat{a}_a]_+ &= \hat{a}^\dagger_i\hat{a}_a+\hat{a}_a\hat{a}^\dagger_i=0\\ [\hat{a}^\dagger_a,\hat{a}_i]_+ &= \hat{a}^\dagger_a\hat{a}_i+\hat{a}_i\hat{a}^\dagger_a=0 \end{aligned} \end{equation} [ a ^ p † , a ^ q † ] + [ a ^ p , a ^ q ] + [ a ^ i † , a ^ j ] + [ a ^ a † , a ^ b ] + [ a ^ i † , a ^ a ] + [ a ^ a † , a ^ i ] + = a ^ p † a ^ q † + a ^ q † a ^ p † = 0 = a ^ p a ^ + a ^ q a ^ p = 0 = a ^ i † a ^ j + a ^ j a ^ i † = δ ij = a ^ a † a ^ b + a ^ b a ^ a † = δ ab = a ^ i † a ^ a + a ^ a a ^ i † = 0 = a ^ a † a ^ i + a ^ i a ^ a † = 0 下面我们来讨论 q q q q q q q q q q q q Partical-Hole 表述 。
以一张表来归纳真空与费米真空之间的异同
真空 费米真空 真空态 ∣ ⟩ \ket{} ∣ ⟩ ∣ Ψ 0 ⟩ \ket{\Psi_0} ∣ Ψ 0 ⟩ 产生算符 a ^ p † \hat{a}^\dagger_p a ^ p † a ^ a † , a ^ i \hat{a}^\dagger_a,\hat{a}_i a ^ a † , a ^ i 湮灭算符 a ^ p \hat{a}_p a ^ p a ^ i † , a ^ a \hat{a}^\dagger_i,\hat{a}_a a ^ i † , a ^ a 产湮算符的物理意义 产生（湮灭）电子 占据轨道的算符为产生（湮灭）空穴，非占据轨道为产生（湮灭）电子 非零对易关系 [ a ^ p † , a ^ q ] + = δ p q [\hat{a}^\dagger_p,\hat{a}_q]_+ =\delta_{pq} [ a ^ p † , a ^ q ] + = δ pq [ a ^ i † , a ^ j ] + = δ i j , [ a ^ a † , a ^ b ] + = δ a b \begin{aligned}[\hat{a}^\dagger_i,\hat{a}_j]_+ = \delta_{ij},[\hat{a}^\dagger_a,\hat{a}_b]_+ = \delta_{ab}\end{aligned} [ a ^ i † , a ^ j ] + = δ ij , [ a ^ a † , a ^ b ] + = δ ab
用一个例题来加深对这些概念的理解
例1. 令 H ^ = ∑ h p q a ^ p † a ^ q \hat{H} = \sum h_{pq}\hat{a}^\dagger_p\hat{a}_q H ^ = ∑ h pq a ^ p † a ^ q H ^ \hat{H} H ^
答:
⟨ Ψ 0 ∣ H ^ ∣ Ψ 0 ⟩ = ⟨ Ψ 0 ∣ ∑ h p q a ^ p † a ^ q ∣ Ψ 0 ⟩ = ∑ h p q ⟨ Ψ 0 ∣ a ^ p † a ^ q ∣ Ψ 0 ⟩ = ∑ h i j ⟨ Ψ 0 ∣ a ^ i † a ^ j ∣ Ψ 0 ⟩ + ∑ h i a ⟨ Ψ 0 ∣ a ^ i † a ^ a ∣ Ψ 0 ⟩ + ∑ h a i ⟨ Ψ 0 ∣ a ^ a † a ^ i ∣ Ψ 0 ⟩ + ∑ h a b ⟨ Ψ 0 ∣ a ^ a † a ^ b ∣ Ψ 0 ⟩ \begin{equation} \begin{aligned} \bra{\Psi_0}\hat{H}\ket{\Psi_0} &= \bra{\Psi_0}\sum h_{pq}\hat{a}^\dagger_p\hat{a}_q\ket{\Psi_0}\\ &=\sum h_{pq} \bra{\Psi_0}\hat{a}^\dagger_p\hat{a}_q\ket{\Psi_0}\\ &=\sum h_{ij} \bra{\Psi_0}\hat{a}^\dagger_i\hat{a}_j\ket{\Psi_0}+\sum h_{ia} \bra{\Psi_0}\hat{a}^\dagger_i\hat{a}_a\ket{\Psi_0}\\ &+\sum h_{ai} \bra{\Psi_0}\hat{a}^\dagger_a\hat{a}_i\ket{\Psi_0}+\sum h_{ab} \bra{\Psi_0}\hat{a}^\dagger_a\hat{a}_b\ket{\Psi_0} \end{aligned} \end{equation} ⟨ Ψ 0 ∣ H ^ ∣ Ψ 0 ⟩ = ⟨ Ψ 0 ∣ ∑ h pq a ^ p † a ^ q ∣ Ψ 0 ⟩ = ∑ h pq ⟨ Ψ 0 ∣ a ^ p † a ^ q ∣ Ψ 0 ⟩ = ∑ h ij ⟨ Ψ 0 ∣ a ^ i † a ^ j ∣ Ψ 0 ⟩ + ∑ h ia ⟨ Ψ 0 ∣ a ^ i † a ^ a ∣ Ψ 0 ⟩ + ∑ h ai ⟨ Ψ 0 ∣ a ^ a † a ^ i ∣ Ψ 0 ⟩ + ∑ h ab ⟨ Ψ 0 ∣ a ^ a † a ^ b ∣ Ψ 0 ⟩ 其中
⟨ Ψ 0 ∣ a ^ i † a ^ j ∣ Ψ 0 ⟩ = ⟨ Ψ 0 ∣ δ i j − a ^ j a ^ i † ∣ Ψ 0 ⟩ = δ i j ⟨ Ψ 0 ∣ a ^ i † a ^ a ∣ Ψ 0 ⟩ = 0 ⟨ Ψ 0 ∣ a ^ a † a ^ i ∣ Ψ 0 ⟩ = − ⟨ Ψ 0 ∣ a ^ i a ^ a † ∣ Ψ 0 ⟩ = 0 ⟨ Ψ 0 ∣ a ^ a † a ^ b ∣ Ψ 0 ⟩ = 0 \begin{equation} \begin{aligned} &\bra{\Psi_0}\hat{a}^\dagger_i\hat{a}_j\ket{\Psi_0}=\bra{\Psi_0}\delta_{ij}-\hat{a}_j\hat{a}^\dagger_i\ket{\Psi_0}=\delta_{ij}\\ &\bra{\Psi_0}\hat{a}^\dagger_i\hat{a}_a\ket{\Psi_0}=0\\ &\bra{\Psi_0}\hat{a}^\dagger_a\hat{a}_i\ket{\Psi_0}=-\bra{\Psi_0}\hat{a}_i\hat{a}^\dagger_a\ket{\Psi_0}=0\\ &\bra{\Psi_0}\hat{a}^\dagger_a\hat{a}_b\ket{\Psi_0}=0 \end{aligned} \end{equation} ⟨ Ψ 0 ∣ a ^ i † a ^ j ∣ Ψ 0 ⟩ = ⟨ Ψ 0 ∣ δ ij − a ^ j a ^ i † ∣ Ψ 0 ⟩ = δ ij ⟨ Ψ 0 ∣ a ^ i † a ^ a ∣ Ψ 0 ⟩ = 0 ⟨ Ψ 0 ∣ a ^ a † a ^ i ∣ Ψ 0 ⟩ = − ⟨ Ψ 0 ∣ a ^ i a ^ a † ∣ Ψ 0 ⟩ = 0 ⟨ Ψ 0 ∣ a ^ a † a ^ b ∣ Ψ 0 ⟩ = 0 将（11）代入（10）中，可得
⟨ Ψ 0 ∣ H ^ ∣ Ψ 0 ⟩ = ∑ h i j δ i j = ∑ h i i \begin{equation} \bra{\Psi_0}\hat{H}\ket{\Psi_0} = \sum h_{ij}\delta_{ij}=\sum h_{ii} \end{equation} ⟨ Ψ 0 ∣ H ^ ∣ Ψ 0 ⟩ = ∑ h ij δ ij = ∑ h ii 2.2 算符的正则序与wick定理 从（6）式和（11）式的计算过程来看，无论是真空还是费米真空，计算这类含产生湮灭算符的积分的主要步骤实际上是通过对易关系，将原本的产生湮灭算符的乘积给化成左边是产生算符右边是湮灭算符的形式，容易证明，这种左产生右湮灭形式的算符作用到真空态上会等于零。我们将这种左产生又湮灭形式的称为算符的 正则序（normal-ordering） 形式，定义正则序算符，其作用是将作用到的算符变换为左产生右湮灭的形式，用符号 { ⋯ } v \{\cdots\}_v { ⋯ } v
{ X ^ Y ^ Z ^ ⋯ U ^ V ^ W ^ } v = ( − 1 ) P ( U ^ X ^ ⋯ Y ^ ⏟ 产生算符 V ^ ⋯ W ^ Z ^ ⏟ 湮灭算符 ) \begin{equation} \{\hat{X}\hat{Y}\hat{Z}\cdots\hat{U}\hat{V}\hat{W}\}_v = (-1)^P(\underbrace{\hat{U}\hat{X}\cdots\hat{Y}}_{\text{产生算符}}\underbrace{\hat{V}\cdots\hat{W}\hat{Z}}_{\text{湮灭算符}}) \end{equation} { X ^ Y ^ Z ^ ⋯ U ^ V ^ W ^ } v = ( − 1 ) P ( 产生算符 U ^ X ^ ⋯ Y ^ 湮灭算符 V ^ ⋯ W ^ Z ^ ) 上式中，P P P
在介绍wick定理之前，我们先介绍一个概念——算符的收缩 。两个算符 A ^ , B ^ \hat{A},\hat{B} A ^ , B ^ A ^ ∙ B ^ ∙ \hat{A}^\bullet\hat{B}^\bullet A ^ ∙ B ^ ∙ 2
A ^ ∙ B ^ ∙ = A ^ B ^ − { A ^ B ^ } v \begin{equation} \hat{A}^\bullet\hat{B}^\bullet=\hat{A}\hat{B}-\{\hat{A}\hat{B}\}_v \end{equation} A ^ ∙ B ^ ∙ = A ^ B ^ − { A ^ B ^ } v 即一个算符的收缩等于算符本身减去算符的正则序。根据正则序的定义式（13）和收缩的定义式（14），我们可以计算后面需要用到的 q q q q q q
a ^ i † ∙ a ^ j ∙ = δ i j a ^ a ∙ a ^ b † ∙ = δ a b \begin{equation} \begin{aligned} \hat{a}^{\dagger\bullet}_i\hat{a}_j^\bullet = \delta_{ij}\\ \hat{a}^{\bullet}_a\hat{a}_b^{\dagger\bullet} = \delta_{ab} \end{aligned} \end{equation} a ^ i †∙ a ^ j ∙ = δ ij a ^ a ∙ a ^ b †∙ = δ ab 在定义了算符的收缩后，我们就能够介绍了。wick定理可以表述为：一个算符可以等于该算符的正则序与对该算符所有可能的收缩的正则序之和，即
A ^ B ^ ⋯ C ^ D ^ = { A ^ B ^ ⋯ C ^ D ^ } v + { A ^ ∙ B ^ ⋯ C ^ ∙ D ^ } v + 所有可能的单收缩的正则序 + { A ^ ∙ B ^ ∙ ∙ ⋯ C ^ ∙ D ^ ∙ ∙ } v + 所有可能的双收缩的正则序 + ⋯ \begin{equation} \begin{aligned} \hat{A}\hat{B}\cdots\hat{C}\hat{D} &= \{\hat{A}\hat{B}\cdots\hat{C}\hat{D}\}_v\\ &+\{\hat{A}^\bullet\hat{B}\cdots\hat{C}^\bullet\hat{D}\}_v+\text{所有可能的单收缩的正则序}\\ &+\{\hat{A}^\bullet\hat{B}^{\bullet\bullet}\cdots\hat{C}^\bullet\hat{D}^{\bullet\bullet}\}_v+\text{所有可能的双收缩的正则序}\\ &+\cdots \end{aligned} \end{equation} A ^ B ^ ⋯ C ^ D ^ = { A ^ B ^ ⋯ C ^ D ^ } v + { A ^ ∙ B ^ ⋯ C ^ ∙ D ^ } v + 所有可能的单收缩的正则序 + { A ^ ∙ B ^ ∙∙ ⋯ C ^ ∙ D ^ ∙∙ } v + 所有可能的双收缩的正则序 + ⋯ wick定理的一个简单证明过程可参见3
3. CC形式理论 3.1 精确波函数的簇展开 下面我们将从波函数的 CI 展开出发，逐步引入耦合簇理论的一些基本概念。
使用上一节引入的费米真空与产生湮灭算符，我们可以把精确波函数的 CI 展开式（2）写成一个激发算符作用在 HF 基态上的形式，即
∣ Ψ ⟩ = ∣ Ψ 0 ⟩ + ∑ i a c i a a ^ a † a ^ i ∣ Ψ 0 ⟩ + ∑ i a j b c i j a b a ^ a † a ^ b † a ^ i a ^ j ∣ Ψ 0 ⟩ + ⋯ = ( 1 + ∑ i a c i a a ^ a † a ^ i + ∑ i a j b c i j a b a ^ a † a ^ b † a ^ i a ^ j + ⋯ ) ∣ Ψ 0 ⟩ = C ^ ∣ Ψ 0 ⟩ \begin{equation} \begin{aligned} \ket{\Psi} &= \ket{\Psi_0}+\sum_{ia} c_i^a\hat{a}^\dagger_a\hat{a}_i\ket{\Psi_0}+\sum_{iajb}c_{ij}^{ab}\hat{a}^\dagger_a\hat{a}^\dagger_b\hat{a}_i\hat{a}_j\ket{\Psi_0}+\cdots\\ &= (1+\sum_{ia} c_i^a\hat{a}^\dagger_a\hat{a}_i+\sum_{iajb}c_{ij}^{ab}\hat{a}^\dagger_a\hat{a}^\dagger_b\hat{a}_i\hat{a}_j+\cdots)\ket{\Psi_0}\\ &= \hat{C}\ket{\Psi_0} \end{aligned} \end{equation} ∣ Ψ ⟩ = ∣ Ψ 0 ⟩ + ia ∑ c i a a ^ a † a ^ i ∣ Ψ 0 ⟩ + iajb ∑ c ij ab a ^ a † a ^ b † a ^ i a ^ j ∣ Ψ 0 ⟩ + ⋯ = ( 1 + ia ∑ c i a a ^ a † a ^ i + iajb ∑ c ij ab a ^ a † a ^ b † a ^ i a ^ j + ⋯ ) ∣ Ψ 0 ⟩ = C ^ ∣ Ψ 0 ⟩ 其中
C ^ = 1 + ∑ i a c i a a ^ a † a ^ i + ∑ i a j b c i j a b a ^ a † a ^ b † a ^ i a ^ j + ⋯ \begin{equation} \begin{aligned} \hat{C}=1+\sum_{ia} c_i^a\hat{a}^\dagger_a\hat{a}_i+\sum_{iajb}c_{ij}^{ab}\hat{a}^\dagger_a\hat{a}^\dagger_b\hat{a}_i\hat{a}_j+\cdots \end{aligned} \end{equation} C ^ = 1 + ia ∑ c i a a ^ a † a ^ i + iajb ∑ c ij ab a ^ a † a ^ b † a ^ i a ^ j + ⋯ 这就是 CI 展开的二次量子化表示。引入簇算符 T ^ \hat{T} T ^
T ^ = ∑ n = 1 ∞ T ^ n \begin{equation} \hat{T} = \sum_{n=1}^\infty\hat{T}_n \end{equation} T ^ = n = 1 ∑ ∞ T ^ n 其中
T ^ 1 = ∑ t i a a ^ a † a ^ i T ^ 2 = ∑ t i j a b a ^ a † a ^ b † a ^ i a ^ j ⋯ \begin{equation} \begin{aligned} \hat{T}_1 &= \sum t_i^a\hat{a}^\dagger_a\hat{a}_i\\ \hat{T}_2 &= \sum t_{ij}^{ab}\hat{a}^\dagger_a\hat{a}^\dagger_b\hat{a}_i\hat{a}_j\\ &\cdots \end{aligned} \end{equation} T ^ 1 T ^ 2 = ∑ t i a a ^ a † a ^ i = ∑ t ij ab a ^ a † a ^ b † a ^ i a ^ j ⋯ 其中各项系数一般被称为 CC 振幅 。令
C ^ = 1 + T ^ + T ^ 2 2 ! + ⋯ = e T ^ \begin{equation} \hat{C} = 1+\hat{T}+\frac{\hat{T}^2}{2!}+\cdots=e^{\hat{T}} \end{equation} C ^ = 1 + T ^ + 2 ! T ^ 2 + ⋯ = e T ^ 可以由上式计算出算符 C ^ \hat{C} C ^ T ^ \hat{T} T ^
c i a = t i a c i j a b = t i j a b + t i a t j b − t i b t j a c i j k a b c = t i j k a b c + t i a t j k b c − t j a t i k b c + t k a t i j b c − t i b t j k a c + t j b t i k a c − t k b t i j a c + t i c t j k a b − t j c t i k a b + t k c t i j a b + t i a t j b t k c − t i a t j c t k b − t i b t j a t k c + t i b t j c t k a + t i c t j a t k b − t i c t j b t k a … \begin{equation} \begin{aligned} c_i^a = &t_i^a\\ c_{ij}^{ab} = &t_{ij}^{ab}+t_i^at_j^b-t_i^bt_j^a\\ c_{ijk}^{abc} = &t_{ijk}^{abc}+t_i^at_{jk}^{bc}-t_j^at_{ik}^{bc}+t_k^at_{ij}^{bc}-t_i^bt_{jk}^{ac}+t_j^bt_{ik}^{ac}-t_k^bt_{ij}^{ac}+t_i^ct_{jk}^{ab}-t_j^ct_{ik}^{ab}\\ &+t_k^ct_{ij}^{ab}+t_i^at_j^bt_k^c-t_i^at_j^ct_k^b-t_i^bt_j^at_k^c+t_i^bt_j^ct_k^a+t_i^ct_j^at_k^b-t_i^ct_j^bt_k^a\\ \ldots \end{aligned} \end{equation} c i a = c ij ab = c ijk ab c = … t i a t ij ab + t i a t j b − t i b t j a t ijk ab c + t i a t jk b c − t j a t ik b c + t k a t ij b c − t i b t jk a c + t j b t ik a c − t k b t ij a c + t i c t jk ab − t j c t ik ab + t k c t ij ab + t i a t j b t k c − t i a t j c t k b − t i b t j a t k c + t i b t j c t k a + t i c t j a t k b − t i c t j b t k a 由此可知，算符 C ^ \hat{C} C ^ T ^ \hat{T} T ^
∣ Ψ ⟩ = e T ^ ∣ Ψ 0 ⟩ \begin{equation} \ket{\Psi} = e^{\hat{T}}\ket{\Psi_0} \end{equation} ∣ Ψ ⟩ = e T ^ ∣ Ψ 0 ⟩ 这就是精确波函数的耦合簇展开。从上述推导过程中可以看到，精确波函数的簇展开与 CI 展开实际上是等价的。将波函数展开成（17）的计算方法被称为组态相互作用（Configure Interaction，CI）理论 ，而将波函数展开成（23）的计算方法被称为耦合簇（Coupled Cluster，CC）理论 。
3.2 CC能量方程与振幅方程 在上一小节，我们废了老大劲才从CI展开式（17）过渡到CC展开式（23），一些初学者可能会感到困惑：为什么要把原本的激发算符给转化成指数算符（23）呢，这种指数形式的算符处理起来与原来的激发算符相比更为简单吗？乍一看这一通操作好像把原来展开式变得更抽象更复杂了，但实际上，将算符转化成指数形式将会大大简化后续的推导。在这一小节中我们将会推导CC能量方程与振幅方程，在推导过程中我们就能体会到指数算符的便利之处。
在将精确波函数写成式（23）的形式后，因为参考态一般是给定的，我们只需做个 HF 单点计算就能得到 HF 基态能量以及基态轨道系数和密度矩阵，所以实际上我们只需要知道振幅的值就能够确定精确波函数的表达式，进而计算能量以及一些其他的我们感兴趣的物理量，从而确定体系的所有性质。
下面我们就来推导CC能量和振幅所满足的方程。将（23）代入定态薛定谔方程，可得
H ^ e T ^ ∣ Ψ 0 ⟩ = E e T ^ Ψ 0 \begin{equation} \hat{H}e^{\hat{T}}\ket{\Psi_0} = Ee^{\hat{T}}\Psi_0 \end{equation} H ^ e T ^ ∣ Ψ 0 ⟩ = E e T ^ Ψ 0 两边同时左乘 e − T ^ e^{-\hat{T}} e − T ^
e − T ^ H ^ e T ^ ∣ Ψ 0 ⟩ = H ^ c c ∣ Ψ 0 ⟩ = E ∣ Ψ 0 ⟩ \begin{equation} e^{-\hat{T}}\hat{H}e^{\hat{T}}\ket{\Psi_0}=\hat{H}_{cc}\ket{\Psi_0}=E\ket{\Psi_0} \end{equation} e − T ^ H ^ e T ^ ∣ Ψ 0 ⟩ = H ^ cc ∣ Ψ 0 ⟩ = E ∣ Ψ 0 ⟩ 其中
H ^ c c = e − T ^ H ^ e T ^ \begin{equation} \hat{H}_{cc}=e^{-\hat{T}}\hat{H}e^{\hat{T}} \end{equation} H ^ cc = e − T ^ H ^ e T ^ 是哈密顿算符 H ^ \hat{H} H ^ ⟨ Ψ 0 ∣ \bra{\Psi_0} ⟨ Ψ 0 ∣
E = ⟨ Ψ 0 ∣ H ^ c c ∣ Ψ 0 ⟩ \begin{equation} E = \bra{\Psi_0}\hat{H}_{cc}\ket{\Psi_0} \end{equation} E = ⟨ Ψ 0 ∣ H ^ cc ∣ Ψ 0 ⟩ 同理，两边同时左乘 ⟨ Ψ i a ∣ \bra{\Psi_i^a} ⟨ Ψ i a ∣ t 1 t_1 t 1
0 = ⟨ Ψ i a ∣ H ^ c c ∣ Ψ 0 ⟩ \begin{equation} 0 = \bra{\Psi_i^a}\hat{H}_{cc}\ket{\Psi_0} \end{equation} 0 = ⟨ Ψ i a ∣ H ^ cc ∣ Ψ 0 ⟩ 同理，两边同时左乘 ⟨ Ψ i j a b ∣ \bra{\Psi_{ij}^{ab}} ⟨ Ψ ij ab ∣ t 2 t_2 t 2
0 = ⟨ Ψ i j a b ∣ H ^ c c ∣ Ψ 0 ⟩ \begin{equation} 0 = \bra{\Psi_{ij}^{ab}}\hat{H}_{cc}\ket{\Psi_0} \end{equation} 0 = ⟨ Ψ ij ab ∣ H ^ cc ∣ Ψ 0 ⟩ 以此类推，最终我们可以得到任意阶振幅所满足的方程，即
{ ⟨ Ψ i a ∣ H ^ c c ∣ Ψ 0 ⟩ = 0 ⟨ Ψ i j a b ∣ H ^ c c ∣ Ψ 0 ⟩ = 0 ⟨ Ψ i j k a b c ∣ H ^ c c ∣ Ψ 0 ⟩ = 0 ⋯ \begin{equation} \begin{cases} \bra{\Psi_i^a}\hat{H}_{cc}\ket{\Psi_0} & =0\\ \bra{\Psi_{ij}^{ab}}\hat{H}_{cc}\ket{\Psi_0} &=0\\ \bra{\Psi_{ijk}^{abc}}\hat{H}_{cc}\ket{\Psi_0} &=0\\ \qquad\cdots \end{cases} \end{equation} ⎩ ⎨ ⎧ ⟨ Ψ i a ∣ H ^ cc ∣ Ψ 0 ⟩ ⟨ Ψ ij ab ∣ H ^ cc ∣ Ψ 0 ⟩ ⟨ Ψ ijk ab c ∣ H ^ cc ∣ Ψ 0 ⟩ ⋯ = 0 = 0 = 0 现在能量方程（27）和振幅方程（30）中只剩下 H ^ c c \hat{H}_{cc} H ^ cc H ^ c c \hat{H}_{cc} H ^ cc
e − A B e A = B + [ B , A ] + 1 2 ! [ [ B , A ] , A ] + 1 3 ! [ [ [ B , A ] , A ] , A ] + ⋯ \begin{equation} e^{-A}Be^{A} = B+[B,A]+\frac{1}{2!}[[B,A],A]+\frac{1}{3!}[[[B,A],A],A]+\cdots \end{equation} e − A B e A = B + [ B , A ] + 2 ! 1 [[ B , A ] , A ] + 3 ! 1 [[[ B , A ] , A ] , A ] + ⋯ 将上式代入方程（26）中，可得
H ^ c c = e − T ^ H ^ e T ^ = H ^ + [ H ^ , T ^ ] + 1 2 ! [ [ H ^ , T ^ ] , T ^ ] + 1 3 ! [ [ [ H ^ , T ^ ] , T ^ ] , T ^ ] + ⋯ \begin{equation} \hat{H}_{cc}=e^{-\hat{T}}\hat{H}e^{\hat{T}}= \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}]+\cdots \end{equation} H ^ cc = e − T ^ H ^ e T ^ = H ^ + [ H ^ , T ^ ] + 2 ! 1 [[ H ^ , T ^ ] , T ^ ] + 3 ! 1 [[[ H ^ , T ^ ] , T ^ ] , T ^ ] + ⋯ 这就是 H ^ c c \hat{H}_{cc} H ^ cc
3.3 讨论 通过上面的推导过程，我们就能很清楚地看到指数算符的优越性。如果我们直接将 CI 的激发算符（18）代入定态薛定谔方程，我们就会发现我们很难直接求得算符 C ^ \hat{C} C ^ 4 H ^ c c \hat{H}_{cc} H ^ cc 5
如果我们将 CI 展开式（2）和簇算符（19）展开至无穷阶，即使用Full CI或Full CC，那么 CI 和 CC 等价，这里上文已经证明过，但是在实际计算中，在大多数情况下，随着展开式阶数的增加，所需的计算时间会呈几何式增加，所以我们一般会对展开式进行截断，从而节省计算时间。对于 CI ，我们一般截断展开式（2），如 CISD 代表在展开式（2）中只保留基态、单激发组态和双激发组态，双激发以上的组态会被舍弃，即令
∣ Ψ ⟩ = ∣ Ψ 0 ⟩ + ∑ i a c i a ∣ Ψ i a ⟩ + ∑ i a j b c i j a b ∣ Ψ i j a b ⟩ \begin{equation} \ket{\Psi} = \ket{\Psi_0}+\sum_{ia}c_i^a\ket{\Psi_i^a}+\sum_{iajb}c_{ij}^{ab}\ket{\Psi_{ij}^{ab}} \end{equation} ∣ Ψ ⟩ = ∣ Ψ 0 ⟩ + ia ∑ c i a ∣ Ψ i a ⟩ + iajb ∑ c ij ab ∣ Ψ ij ab ⟩ 而在 CC 中，我们一般会截断簇算符的表达式（19），如 CCSD 代表在簇算符中只保留单激发和双激发部分，其他激发算符会被舍弃，即令
T ^ = T ^ 1 + T ^ 2 ∣ Ψ ⟩ = e T ^ ∣ Ψ 0 ⟩ = e T ^ 1 + T ^ 2 ∣ Ψ 0 ⟩ \begin{equation} \begin{aligned} \hat{T} &= \hat{T}_1+\hat{T}_2\\ \ket{\Psi} &=e^{\hat{T}}\ket{\Psi_0}= e^{\hat{T}_1+\hat{T}_2}\ket{\Psi_0} \end{aligned} \end{equation} T ^ ∣ Ψ ⟩ = T ^ 1 + T ^ 2 = e T ^ ∣ Ψ 0 ⟩ = e T ^ 1 + T ^ 2 ∣ Ψ 0 ⟩ 在这里，一些读者可能就要说了：在对簇算符进行截断后，式（27）（30）不还是包含 H ^ c c \hat{H}_{cc} H ^ cc
在截断的 CI 中，我们仅保留低阶的激发组态，把高阶激发组态全部舍弃，即假设高阶激发组态的组态系数远小于低阶激发组态的组态系数，这种截断方案其实并不妥当，因为这个假设并不一定成立，高阶组态系数未必一定会比低阶组态系数小，所以使用截断 CI计算得到的相关能精度十分有限。而截断的 CC 并不存在这个问题。由于在截断簇算符后，H ^ c c \hat{H}_{cc} H ^ cc
附录A. Baker–Campbell–Hausdorff 公式的简单证明 这里的证明步骤参考了6 U ( x ) U(x) U ( x )
U ( x ) = e − A x B e A x \begin{equation} U(x) = e^{-Ax}Be^{Ax} \end{equation} U ( x ) = e − A x B e A x 其中 x x x x x x U ( x ) U(x) U ( x )
d d x U = − A e − A x B e A x + e − A x B A e A x \begin{equation} \frac{d}{dx}U=-Ae^{-Ax}Be^{Ax}+e^{-Ax}BAe^{Ax} \end{equation} d x d U = − A e − A x B e A x + e − A x B A e A x 由于 算符 A A A 7
d d x U = − A U + U A = [ U , A ] \begin{equation} \frac{d}{dx}U=-AU+UA=[U,A] \end{equation} d x d U = − A U + U A = [ U , A ] 这就是算符函数 U ( x ) U(x) U ( x ) U ( 0 ) = B U(0)=B U ( 0 ) = B
U ( x ) = B + ∫ 0 x [ U ( t ) , A ] d t \begin{equation} U(x) = B+\int_0^x[U(t),A]dt \end{equation} U ( x ) = B + ∫ 0 x [ U ( t ) , A ] d t 将 U ( x ) U(x) U ( x ) x x x
左边 = U 0 + U 1 x + 1 2 ! U 2 x 2 + ⋯ 右边 = B + ∫ 0 x U ( t ) A d t − ∫ 0 x A U ( t ) d t = B + [ U 0 , A ] x + 1 2 ! [ U 1 , A ] x 2 + ⋯ \begin{equation} \begin{aligned} \text{左边}&=U_0+U_1x+\frac{1}{2!}U_2x^2+\cdots\\ \text{右边}&=B+\int_0^xU(t)Adt-\int_0^xAU(t)dt\\ &=B+[U_0,A]x+\frac{1}{2!}[U_1,A]x^2+\cdots \end{aligned} \end{equation} 左边 右边 = U 0 + U 1 x + 2 ! 1 U 2 x 2 + ⋯ = B + ∫ 0 x U ( t ) A d t − ∫ 0 x A U ( t ) d t = B + [ U 0 , A ] x + 2 ! 1 [ U 1 , A ] x 2 + ⋯ 对比两边系数，可得
U 0 = B U 1 = [ U 0 , A ] = [ B , A ] U 2 = [ U 1 , A ] = [ [ B , A ] , A ] ⋯ \begin{equation} \begin{aligned} U_0 &= B\\ U_1 &=[U_0,A]=[B,A]\\ U_2 &=[U_1,A]=[[B,A],A]\\ &\cdots \end{aligned} \end{equation} U 0 U 1 U 2 = B = [ U 0 , A ] = [ B , A ] = [ U 1 , A ] = [[ B , A ] , A ] ⋯ 将系数代入（39）中，即可得算符函数 U ( x ) U(x) U ( x )
U ( x ) = U 0 + U 1 x + 1 2 ! U 2 x 2 + ⋯ = e − A x B e A x = B + [ B , A ] x + 1 2 ! [ [ B , A ] , A ] x 2 + 1 3 ! [ [ [ B , A ] , A ] , A ] x 3 + ⋯ \begin{equation} \begin{aligned} U(x)&=U_0+U_1x+\frac{1}{2!}U_2x^2+\cdots=e^{-Ax}Be^{Ax} \\ &= B+[B,A]x+\frac{1}{2!}[[B,A],A]x^2+\frac{1}{3!}[[[B,A],A],A]x^3+\cdots \end{aligned} \end{equation} U ( x ) = U 0 + U 1 x + 2 ! 1 U 2 x 2 + ⋯ = e − A x B e A x = B + [ B , A ] x + 2 ! 1 [[ B , A ] , A ] x 2 + 3 ! 1 [[[ B , A ] , A ] , A ] x 3 + ⋯ 令 x = 1 x=1 x = 1
U ( 1 ) = U 0 + U 1 + 1 2 ! U 2 + ⋯ = e − A B e A = B + [ B , A ] + 1 2 ! [ [ B , A ] , A ] + 1 3 ! [ [ [ B , A ] , A ] , A ] + ⋯ \begin{equation} \begin{aligned} U(1)&=U_0+U_1+\frac{1}{2!}U_2+\cdots=e^{-A}Be^{A} \\ &= B+[B,A]+\frac{1}{2!}[[B,A],A]+\frac{1}{3!}[[[B,A],A],A]+\cdots \end{aligned} \end{equation} U ( 1 ) = U 0 + U 1 + 2 ! 1 U 2 + ⋯ = e − A B e A = B + [ B , A ] + 2 ! 1 [[ B , A ] , A ] + 3 ! 1 [[[ B , A ] , A ] , A ] + ⋯ 