1. Gaussian 中的 PUHF/PMP2 结果的重新实现

创建时间：2019-08-31，最后修改：2019-09-01

在这一份笔记中，我们将使用 PySCF 的功能与 NumPy 重复 Gaussian 中计算的 PUHF 与 PMP2 能量结果；并对 PUHF 与 PMP2 的推导作简单的说明。

from pyscf import gto, scf, mp

1.1. 参考结果与体系定义

1.1.1. Gaussian 结果

在 Gaussian 中，我们使用以下输入卡可以得到 PUHF/PMP2 能量：

#p UMP2(Full)/6-31G nosymm

H2O

3 4
O  0.  0.  0.
H  1.  0.  0.
H  0.  1.  0.

对于上述分子，其中一些重要的输出结果是

1. \(E_\mathrm{UHF}\)：-73.0451423839

2. \(E_\mathrm{UMP2, corr}\): -0.02646719276

3. \(E_\mathrm{UMP2}\): -73.071609576661

4. \(\langle S_z \rangle\): 1.5

5. \(\langle S^{2(0)} \rangle\): 3.7531

6. \(\langle S^{2(0)} \rangle + \langle S^{2(1)} \rangle\): 3.7504

7. \(E_\mathrm{PUHF}\)：-73.046146318

8. \(E_\mathrm{PMP2}\): -73.072180589

输出文件参见 assets/PUHF_and_PMP2.out；其中，有效的数据可以通过下述的代码获得：

with open("assets/PUHF_and_PMP2.out", "r") as output:
    output_lines = output.read().split("\n")

for line_num, line_text in enumerate(output_lines):
    if any([keyword in line_text for keyword in
            ["SCF Done", "EUMP2", "<S**2>", "(S**2,1)", "E(PMP2)"]]) \
        and "Initial guess" not in line_text:
        print("line {:03d}: {}".format(line_num, line_text))

line 336:  SCF Done:  E(UHF) =  -73.0451423839     A.U. after   11 cycles
line 338:  <Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 1.5000 <S**2>= 3.7531 S= 1.5008
line 370:  (S**2,0)=  0.37531D+01           (S**2,1)=  0.37504D+01
line 371:  E(PUHF)=      -0.73046146318D+02        E(PMP2)=      -0.73072180589D+02
line 373:  E2 =    -0.2646719276D-01 EUMP2 =    -0.73071609576661D+02

我们的目标就是近乎准确无误地重复上述八个结果。

1.1.2. PySCF 体系定义

为了获得与 Gaussian 相同的结果，我们需要定义相同的分子与电荷、多重度环境：

mol = gto.Mole()
mol.atom = """
O 0. 0. 0.
H 1. 0. 0.
H 0. 1. 0.
"""
mol.charge = 3
mol.spin = 3
mol.basis = "6-31G"
mol.build()

<pyscf.gto.mole.Mole at 0x7fa4e5342160>

通过 PySCF 计算 UHF 能量：

scf_eng = scf.UHF(mol)
scf_eng.conv_tol = 1e-10
scf_eng.run();

converged SCF energy = -73.0451423536459  <S^2> = 3.7530824  2S+1 = 4.0015409

上述结果应当能与 \(E_\mathrm{UHF}\) 和 \(\langle S^{2(0)} \rangle\) 对应。\(\langle S_z \rangle = 1.5\) 几乎是显然的。不过，我们仍然不了解 \(\langle S^{2(0)} \rangle\) 是如何生成的。

通过 PySCF 计算 UMP2 能量：

mp2_eng = mp.UMP2(scf_eng)
mp2_eng.run();

E(UMP2) = -73.0716095455079  E_corr = -0.0264671918619837

上述结果应当能与 \(E_\mathrm{UMP2, corr}\) 和 \(E_\mathrm{UMP2}\) 对应。

因此，当前的问题将是回答：如何重复

1. \(\langle S^{2(0)} \rangle\): 3.7531

2. \(\langle S^{2(0)} \rangle + \langle S^{2(1)} \rangle\): 3.7504

3. \(E_\mathrm{PUHF}\)：-73.046146318

4. \(E_\mathrm{PMP2}\): -73.072180589

1.1.3. 部分变量定义

首先，我们遵从大多数量化文章中的记号

• \(i, j\) 代表占据分子轨道

• \(a, b\) 代表非占分子轨道

• \(p, q\) 代表任意分子轨道

• \(\alpha, \beta\) 代表任意原子轨道

Table 1. 分子相关变量

变量名	元素记号	意义与注解	标量或区间
nocc_a	\(n_\mathrm{occ}^\alpha\)	\(\alpha\) 自旋电子数	\(5\)
nocc_b	\(n_\mathrm{occ}^\beta\)	\(\beta\) 自旋电子数	\(2\)
N	\(N\)	总电子数	\(7\)
nmo	\(n_\mathrm{MO}\)	分子轨道数	\(13\)
nao	\(n_\mathrm{AO}\)	原子轨道数	\(13\)
S	\(S_{\mu \nu}\)	原子轨道重叠积分
so_a		\(\alpha\) 占据轨道分割	\([0, 5)\)
so_b		\(\beta\) 占据轨道分割	\([0, 2)\)
sv_a		\(\alpha\) 非占轨道分割	\([5, 13)\)
sv_b		\(\beta\) 非占轨道分割	\([2, 13)\)
Sx	\(S_x\)	\(x\) 分量自旋	\(0\)
Sy	\(S_y\)	\(y\) 分量自旋	\(0\)
Sz	\(S_z\)	\(z\) 分量自旋	\(3/2\)

Table 2. UHF 计算相关变量

变量名	元素记号	意义与注解
C_a	\(C_{\mu p}\)	\(\alpha\) 系数矩阵
C_b	\(C_{\mu \bar p}\)	\(\beta\) 系数矩阵
e_a	\(e_{p}\)	\(\alpha\) 轨道能
e_b	\(e_{\bar p}\)	\(\alpha\) 轨道能
eo_a	\(e_{i}\)	\(\beta\) 占据轨道能
eo_b	\(e_{\bar i}\)	\(\alpha\) 占据轨道能
ev_a	\(e_{a}\)	\(\alpha\) 非占轨道能
ev_b	\(e_{\bar a}\)	\(\beta\) 非占轨道能
D2_aa	\(D_{ij}^{ab}\)	\(\alpha, \alpha\) 轨道能差
D2_ab	\(D_{i \bar j}^{a \bar b}\)	\(\alpha, \beta\) 轨道能差
D2_bb	\(D_{\bar i \bar j}^{\bar a \bar b}\)	\(\beta, \beta\) 轨道能差

Table 3. UMP2 计算相关变量

变量名	元素记号	意义与注解
t2_aa	\(t_{ij}^{ab}\)	\(\alpha, \alpha\) MP2 激发系数
t2_ab	\(t_{i \bar j}^{a \bar b}\)	\(\alpha, \beta\) MP2 激发系数
t2_bb	\(t_{\bar i \bar j}^{\bar a \bar b}\)	\(\beta, \beta\) MP2 激发系数
D2_aa	\(D_{ij}^{ab}\)	\(\alpha, \alpha\) MP2 激发系数分母
D2_ab	\(D_{i \bar j}^{a \bar b}\)	\(\alpha, \beta\) MP2 激发系数分母
D2_bb	\(D_{\bar i \bar j}^{\bar a \bar b}\)	\(\beta, \beta\) MP2 激发系数分母

上述需要补充说明的公式有：

\[ S_z = \frac{1}{2} (n_\mathrm{occ}^\alpha - n_\mathrm{occ}^\beta) \]

\[ D_{i \bar j}^{a \bar b} = e_i + e_{\bar j} - e_a - e_{\bar b} \]

对于 MP2 激发系数分母，另外两种自旋情况的 \(D_{ij}^{ab}\) 与 \(D_{\bar i \bar j}^{\bar a \bar b}\) 也可以类似地生成。

# === Molecular
# --- Definition
nocc_a, nocc_b = mol.nelec
N = nocc_a + nocc_b
nmo = nao = mol.nao
S = mol.intor("int1e_ovlp")
# --- Derivative
so_a, so_b = slice(0, nocc_a), slice(0, nocc_b)
sv_a, sv_b = slice(nocc_a, nmo), slice(nocc_b, nmo)
Sx, Sy, Sz = 0, 0, 0.5 * (nocc_a - nocc_b)

# === UHF Calculation
# --- Definition
C_a, C_b = scf_eng.mo_coeff
e_a, e_b = scf_eng.mo_energy
# --- Derivative
eo_a, eo_b = e_a[so_a], e_b[so_b]
ev_a, ev_b = e_a[sv_a], e_b[sv_b]
D2_aa = eo_a[:, None, None, None] + eo_a[None, :, None, None] - ev_a[None, None, :, None] - ev_a[None, None, None, :]
D2_ab = eo_a[:, None, None, None] + eo_b[None, :, None, None] - ev_a[None, None, :, None] - ev_b[None, None, None, :]
D2_bb = eo_b[:, None, None, None] + eo_b[None, :, None, None] - ev_b[None, None, :, None] - ev_b[None, None, None, :]

# === MP2 Calculation
t2_aa, t2_ab, t2_bb = mp2_eng.t2

作为对四脚标张量性质的验证，我们计算 MP2 相关能 \(E_\mathrm{MP2, corr}\) 如下：

\[ E_\mathrm{MP2, corr} = \frac{1}{4} \sum_{ijab} (t_{ij}^{ab})^2 D_{ij}^{ab} + \frac{1}{4} \sum_{\bar i \bar j \bar a \bar b} (t_{\bar i \bar j}^{\bar a \bar b})^2 D_{i \bar j}^{a \bar b} + \sum_{i \bar j a \bar b} (t_{i \bar j}^{a\bar b})^2 D_{i \bar j}^{a \bar b} \]

(+ 0.25 * (t2_aa**2 * D2_aa).sum()
 + 0.25 * (t2_bb**2 * D2_bb).sum()
 + (t2_ab**2 * D2_ab).sum())

-0.02646719186198365

PySCF 所给出的 \(E_\mathrm{MP2, corr}\) 可以给出相同的结果：

mp2_eng.e_corr

-0.02646719186198366

1.2. \(\langle S^2 \rangle\) 相关计算

1.2.1. 分子轨道基组重叠矩阵 S_pq \(S_{p \bar q}\)

\[ S_{p \bar q} = \sum_{\mu \nu} C_{\mu p} S_{\mu \nu} C_{\nu \bar q} \]

若用量子力学记号，上述矩阵元素可能表示为

\[ S_{p \bar q} = \int \phi_p (\boldsymbol{r}) \phi_{\bar q} (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} \]

注意上述的积分是空间坐标的积分，不包含自旋部分的积分。

S_pq = C_a.T @ S @ C_b
S_pq.shape

(13, 13)

我们以后还会使用上述矩阵的占据-占据部分 S_ij \(S_{i \bar j}\)、占据-非占部分 S_ia \(S_{i \bar a}\) 与非占-占据部分 S_ai \(S_{a \bar i} = S_{\bar i a}\)：

S_ij, S_ia, S_ai = S_pq[so_a, so_b], S_pq[so_a, sv_b], S_pq[sv_a, so_b]
[S_ij.shape, S_ia.shape, S_ai.shape]

[(5, 2), (5, 11), (8, 2)]

1.2.2. S2_0 \(\langle S^{2(0)} \rangle\)

\(\langle S^{2(0)} \rangle\) 在程序中通常写为 <S^2> 或 <S**2>。在 Gaussian 计算 PUHF 处，还写为 (S**2,0)。这意味着是 UHF 波函数的 \(\langle S^2 \rangle_\mathrm{UHF}\)。相对地，UMP2 波函数给出的对 \(\langle S^2 \rangle\) 的矫正将记作 \(\langle S^{2(1)} \rangle\)。

参考 Chen and Schlegel [1] Table 1, \(0 \rightarrow 0\) 或等价地，Szabo and Ostlund [2] eq (2.271)

\[ \langle S^{2(0)} \rangle = \langle \Psi_0 | \hat S^2 | \Psi_0 \rangle = S_z (S_z + 1) + n_\mathrm{occ}^\beta - \sum_{i \bar j} (S_{i \bar j})^2 \]

S2_0 = Sz * (Sz + 1) + nocc_b - (S_ij**2).sum()
S2_0

3.7530823820890378

Gaussian 的参考值是 3.7531。

为了以后的记号便利，我们在这里定义 L

\[ L = \sum_{i \bar j} (S_{i \bar j})^2 \]

L = (S_ij**2).sum()

1.2.3. S2_1 \(\langle S^{2(1)} \rangle\)

\[\begin{split} \begin{align} \langle S^{2(1)} \rangle &= 2 \langle \Psi_0 | \hat S^2 | \Psi^{(1)} \rangle = 2 \sum_{i \bar j a \bar b} t_{i \bar j}^{a \bar b} \langle \Psi_0 | \hat S^2 | \Psi_{i \bar j}^{a \bar b} \rangle \\ &= - 2 \sum_{i \bar j a \bar b} t_{i \bar j}^{a \bar b} \langle i | \bar b \rangle \langle a | \bar j \rangle = - 2 \sum_{i \bar j a \bar b} t_{i \bar j}^{a \bar b} S_{i \bar b} S_{a \bar j} \end{align} \end{split}\]

上式的第一个等号是 Chen and Schlegel [1] eq (5) 所给出的；而第三个等号是 Table 1 \(0 \rightarrow \alpha \beta (i, a: \alpha; j, b: \beta)\) 给出的。

上式的推导中有一处关于 \(| \Psi^{(1)} \rangle\) 的展开的推导省略。我们知道

\[ | \Psi^{(1)} \rangle = \hat T_2 | \Psi_0 \rangle = \frac{1}{4} \sum_{ijab} t_{ij}^{ab} | \Psi_{ij}^{ab} \rangle + \frac{1}{4} \sum_{\bar i \bar j \bar a \bar b} t_{\bar i \bar j}^{\bar a \bar b} | \Psi_{\bar i \bar j}^{\bar a \bar b} \rangle + \sum_{i \bar j a \bar b} t_{i \bar j}^{a \bar b} | \Psi_{i \bar j}^{a \bar b} \rangle \]

但由于利用到 \(\langle 0 | \hat S^2 | \Psi_{ij}^{ab} \rangle = \langle 0 | \hat S^2 | \Psi_{\bar i \bar j}^{\bar a \bar b} \rangle = 0\)，因此在第二个等号时只从三个 \(| \Psi^{(1)} \rangle\) 中留下了一项。关于 \(\hat S^2\) 作用在 UHF 波函数与轨道下的性质，可以参考 Schlegel [3] eq (5) 的说明。

S2_1 = - 2 * (t2_ab * S_ia[:, None, None, :] * S_ai.T[None, :, :, None]).sum()
S2_1

-0.002658400959213991

因此，UMP2 矫正过的 \(\langle S^2 \rangle_\mathrm{UMP2} = \langle S^{2(0)} \rangle + \langle S^{2(1)} \rangle\) 的结果是

S2_0 + S2_1

3.7504239811298237

Gaussian 的参考值是 3.7504。

1.2.4. S4SD \(\texttt{S4SD}\)

S4SD 的表达式较为复杂，我们也直接使用 \(\texttt{S4SD}\) 而不用其他记号表示该项：

\[ \begin{align} \texttt{S4SD} = (n_\mathrm{occ}^\alpha - L) (n_\mathrm{occ}^\beta - L) + 2 L - 2 \sum_{i \bar j k \bar l} S_{i \bar j} S_{\bar j k} S_{k \bar l} S_{\bar l i} + \langle S^{2(0)} \rangle^2 \end{align} \]

S4SD = (nocc_a - L) * (nocc_b - L) + 2 * L - 2 * (S_ij @ S_ij.T @ S_ij @ S_ij.T).trace() + S2_0**2
S4SD

14.101029785519533

该表达式的来源可能是 Amos and Hall [4]。该文的 eq (7·02) 下方公式中，有通过稍高阶的投影而获得的 \(\langle S^2 \rangle\) 的计算方式

\[ \langle S^2 \rangle \simeq \langle S^{2(0)} \rangle + \frac{\texttt{S4SD} - \langle S^{2(0)} \rangle^2}{\langle S^{2(0)} \rangle - (S_z + 1) (S_z + 2)} \]

通过这种方式获得的 \(\langle S^2 \rangle\) 近似值可以相当精确，比 \(\langle S^{2(0)} \rangle + \langle S^{2(1)} \rangle\) 还要接近精确值 \(3.75\)：

S2_0 + (S4SD - S2_0**2) / (S2_0 - (Sz + 1) * (Sz + 2))

3.749999998117527

相信 \(\texttt{S4SD}\) 的存在意义是用于计算 Schlegel [3] 式 eq (24) 中的 \(\langle \tilde \Phi_1 | \tilde \Phi_1 \rangle = \langle \Phi_0 | A_{s + 1}^\dagger A_{s + 1} | \Phi_0 \rangle\)；但关于这一关系我还不确定是否正确。后面计算 PMP2 能量时会使用上 \(\texttt{S4SD}\)。

1.3. 自旋污染矫正能量计算

1.3.1. EPUHF \(E_\mathrm{PUHF}\)

根据 Schlegel [3] eq (22)，PUHF 能量可以表达为

\[ E_\mathrm{PUHF} = E_\mathrm{UHF} + \frac{1}{\langle \Psi_0 | \hat P_s | \Psi_0 \rangle} \sum_{i \bar j a \bar b} \langle \Psi_0 | \hat H | \Psi_{i \bar j}^{a \bar b} \rangle \langle \Psi_{i \bar j}^{a \bar b} | \hat P_s | \Psi_0 \rangle \]

其中，\(\hat P_s\) 算符称为 Löwdin 算符 [5] eq (7)，

\[ \hat P_s = \prod_{k \neq s}^{N / 2} \frac{\hat S^2 - k (k + 1)}{s (s + 1) - k (k + 1)} \]

相当于将自旋不纯的波函数纯化为自旋量子数为 \(s\) 的态。在实际使用中，通常使用 \(\hat A_{s + 1} \simeq \hat P_s\) 替代；关于这段讨论可以参考 Rossky and Karplus [6] section V.A 的讨论，而下面公式的形式参考 Schlegel [3] eq (14)；其中，\(s\) 一般取 \(S_z\)：

\[ \hat A_{s + 1} = \frac{\hat S^2 - (s + 1)(s + 2)}{\langle S^{2(0)} \rangle - (s + 1)(s + 2)} \]

关于 \(\hat A_{s + 1}\)，一个显然的性质是 \(\langle \Psi_0 | \hat A_{s + 1} | \Psi_0 \rangle = 1\)。

为了程序方便，定义下述临时变量 Y

\[ Y = \langle S^{2(0)} \rangle - (S_z + 1) (S_z + 2) \]

那么 D_EPUHF

\[\begin{split} \begin{align} \Delta E_\mathrm{PUHF} &= \sum_{i \bar j a \bar b} t_{i \bar j}^{a \bar b} D_{i \bar j}^{a \bar b} \cdot \langle \Psi_{i \bar j}^{a \bar b} | \frac{\hat S^2}{Y} | \Psi_0 \rangle \\ &= - \frac{1}{Y} \sum_{i \bar j a \bar b} t_{i \bar j}^{a \bar b} D_{i \bar j}^{a \bar b} S_{i \bar b} S_{\bar j a} \end{align} \end{split}\]

Y = S2_0 - (Sz + 1) * (Sz + 2)
D_EPUHF = - 1 / Y * (t2_ab * D2_ab * S_ia[:, None, None, :] * S_ai.T[None, :, :, None]).sum()
D_EPUHF

-0.0010039336320381543

因而 \(E_\mathrm{PUHF} = E_\mathrm{UHF} + \Delta E_\mathrm{PUHF}\)：

scf_eng.e_tot + D_EPUHF

-73.04614628727792

Gaussian 的参考值是 -73.046146318。

1.3.2. EPMP2 \(E_\mathrm{PMP2}\)

根据 Schlegel [3] eq (24)，PMP2 能量可以表达为

\[ \begin{align} \Delta E_\mathrm{PMP2} = \Delta E_\mathrm{PUHF} \left( 1 - \frac{\langle \Phi^{(1)} | \hat A_{s + 1} | \Psi_0 \rangle}{\langle \Phi_0 | \hat A_{s + 1}^\dagger \hat A_{s + 1} | \Psi_0 \rangle} \right) \end{align} \]

关于上式的分数项，分子部分可以写为

\[ \begin{align} \langle \Phi^{(1)} | \hat A_{s + 1} | \Psi_0 \rangle = \langle \Phi^{(1)} | \frac{\hat S^2}{Y} - \frac{(s + 1)(s + 2)}{Y} | \Psi_0 \rangle = \frac{1}{2} \frac{\langle S^{2(1)} \rangle}{Y} \end{align} \]

而关于分子项，参考在 \(\texttt{S4SD}\) 的讨论，

\[ \langle \Phi_0 | \hat A_{s + 1}^\dagger \hat A_{s + 1} | \Psi_0 \rangle \simeq \langle S^2 \rangle - \langle S^{2(0)} \rangle = \frac{\texttt{S4SD} - \langle S^{2(0)} \rangle^2}{Y^2} \]

但作者不能断定上述论断的正确性。

将分子、分母的结果代入 \(\Delta E_\mathrm{PMP2}\) 的算式中，可以得到 D_EPMP2

\[ \Delta E_\mathrm{PMP2} = \Delta E_\mathrm{PUHF} \left( 1 - \frac{1}{2} \frac{\langle S^{2(1)} \rangle \cdot Y}{\texttt{S4SD} - \langle S^{2(0)} \rangle^2} \right) \]

D_EPMP2 = D_EPUHF * (1 - 0.5 * S2_1 * Y / (S4SD - S2_0**2))
D_EPMP2

-0.0005710125302116336

因而 \(E_\mathrm{PMP2} = E_\mathrm{UMP2} + \Delta E_\mathrm{PMP2}\)：

mp2_eng.e_tot + D_EPMP2

-73.07218055803808

Gaussian 的参考值是 -73.072180589。

至此，我们已经完成了使用 PySCF 的功能与 NumPy 重复 Gaussian 的 PUHF、PMP2 的能量结果了。

1.4. 修订时间轴

• 2019/08/30 写完文档；文档基于 2019/08/13 的一份笔记。

• 2019/09/01 补充一部分推导。

---

[1] (1,2)

Chen, W.; Schlegel, H. B. Evaluation of S2 for Correlated Wave Functions and Spin Projection of Unrestricted Moller–Plesset Perturbation Theory. J. Chem. Phys. 1994, 101 (7), 5957–5968. doi: 10.1063/1.467312.

[2]

Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books on Chemistry); Dover Publications, 1996.

[3] (1,2,3,4,5)

Schlegel, H. B. Potential Energy Curves Using Unrestricted Mo/ller–Plesset Perturbation Theory with Spin Annihilation. J. Chem. Phys. 1986, 84 (8), 4530–4534. doi: 10.1063/1.450026.

[4]

Amos, A. T.; Hall, G. G. Single Determinant Wave Functions. Proc. R. Soc. Lond. A 1961, 263 (1315), 483–493. doi: 10.1098/rspa.1961.0175.

[5]

Lowdin, P.-O. Quantum Theory of Many-Particle Systems. III. Extension of the Hartree-Fock Scheme to Include Degenerate Systems and Correlation Effects. Phys. Rev. 1955, 97 (6), 1509–1520. 10.1103/physrev.97.1509.

[6]

Rossky, P. J.; Karplus, M. Spin Dependent Properties of Perturbed Wave Functions: An Analytic Comparison of the Exact, UHF, and Spin-Projected UHF States. J. Chem. Phys. 1980, 73 (12), 6196–6214. 10.1063/1.440115.