6. CP-HF 方程相关计算

在这一节，我们会了解实数 Restricted Hartree-Fock 参考态下的 CP-HF 方程给的解法；这会作为下一节解决 xDH 型泛函，特别是 XYG3 的 CP-KS 方程的前置学习．通过一阶 CP-HF 方程，我们可以求解 HF 方法的极化率，以及 MP2 方法的偶极矩．

补充说明

通常来说，解旋转矩阵的方程称为 CP 方程，而在 MP2 中解单电子密度的方程称为 Z-Vector 方程．在之后，这两个问题将统称为 CP 方程．

同时，这一节的 CP 方程始终是一阶的．

6.1. 准备工作

6.1.1. 分子计算

# 环境搭建
import psi4
import numpy as np

# 引入 DIIS 模块
import sys
sys.path.append("include")
from DIIS_helper import DIIS_helper

# 简化矩阵输出
np.set_printoptions(8, linewidth=100, suppress=True)

# 输出文件
psi4.set_output_file("output.dat", False)

# 设置内存 0.5 GB
psi4.set_memory(int(5e8))

# 设置分子坐标
mol = psi4.geometry("""
    O  0.000000000000    -0.000000000000    -0.079135765807
    H  0.000000000000     0.707106781187     0.627971015380
    H  0.000000000000    -0.707106781187     0.627971015380
    symmetry c1
""")

# 设置计算选项
psi4.set_options({
    'basis':               '6-31g',
    'scf_type':            'pk',
    'mp2_type':            'conv',
    'e_convergence':        1e-8,
    'd_convergence':        1e-8,
})

# 波函数信息
scf_e, scf_wfn = psi4.energy("HF", molecule=mol, return_wfn=True)

# 验证结果
psi4.compare_values(scf_e, -75.9697009555, 6, 'HF Energy')

        HF Energy.........................................................PASSED

True

6.1.2. 中间变量

# 电子积分引擎
mints = psi4.core.MintsHelper(scf_wfn.basisset())
# 积分
T = np.asarray(mints.ao_kinetic())  # AO 基组动能积分
V = np.asarray(mints.ao_potential())  # AO 基组电子-核势能积分
eri = np.asarray(mints.ao_eri())  # AO 基组双电子排斥积分

# 数值量
nocc = scf_wfn.nalpha()
nbf = mints.nbf()
nmo = scf_wfn.nmo()
nvir = nmo - nocc

# HF 导出变量
F = np.asarray(scf_wfn.Fa())  # AO 基组 Fock 矩阵
D = np.asarray(scf_wfn.Da())  # AO 基组单电子密度
C = np.asarray(scf_wfn.Ca())  # AO <-> MO 轨道系数 C_{up}
Co = np.asarray(scf_wfn.Ca_subset('AO', 'OCC'))  # 占据轨道系数
Cv = np.asarray(scf_wfn.Ca_subset('AO', 'VIR'))  # 未占轨道系数
e = np.asarray(scf_wfn.epsilon_a())  # 轨道能级
eo = np.asarray(scf_wfn.epsilon_a_subset('AO', 'OCC'))  # 占据轨道能级
ev = np.asarray(scf_wfn.epsilon_a_subset('AO', 'VIR'))  # 未占轨道能级
dip_psi4 = mints.ao_dipole()
dip = np.array([np.asarray(mat) for mat in dip_psi4])  # AO 基组偶极矩积分

# MP2 导出变量
# 轨道能之差的张量 D_{ij}^{ab} -> d[i, a, j, b]
d_ovov = (eo.reshape(-1, 1, 1, 1) - ev.reshape(-1, 1, 1) + eo.reshape(-1, 1) - ev)
# 全 MO 基组 ERI，定义与上一节不同 <pq|rs> = (pr|qs) = g_{pq}^{rs} = g_{rs}^{pq} -> g[p, r, q, s]
g_pqrs = np.einsum("up, vr, uvkl, kq, ls -> prqs", C, C, eri, C, C, optimize=True)
g_ovov = g_pqrs[:nocc, nocc:, :nocc, nocc:]  # <ij|ab> = g_{ij}^{ab}
t_ovov = g_ovov / d_ovov  # 轨道对激发振幅

# 验证 MP2 能量
mp2_corr_my = (t_ovov ** 2 * d_ovov + 0.5 * (t_ovov - t_ovov.swapaxes(1,3)) ** 2 * d_ovov).sum()
psi4.compare_values(mp2_corr_my, -0.1343346885, 6, 'Constructed MP2 Correlation')

        Constructed MP2 Correlation.......................................PASSED

True

6.2. Hartree-Fock 极化率理论回顾

6.2.1. 记号重新定义

在这里，我们会使用 Yamaguchi 书 中的许多结论，但会对其中一些公式中的记号作重新定义．下面的表格是对必要的记号所作的声明．

意义	Yamaguchi	重新定义	代码调用
坐标轴	\(f, g\)
电场分量	\(F, G\)	\(E_f, E_g\)
偶极分量算符	\(r_f\)	\(-f\)
AO 偶极积分	\(d_{\mu \nu}^f\)	\(\mu_{\mu \nu}^f\)	dip[f][u, v]
任意物理量	\(a\)	\(\xi\)
MO ERI	\((pq \vert rs)\)	\(g_{pr}^{qs}\)	g[p, q, r, s]
	\(A_{pq,rs}\)	\(A_{pr}^{qs}\)	A[p, q, r, s]
	\(B_{0,pq}^a\)	\(B_{pq}^f\)	B[p, q]
电荷	\(e\)	\(1\)

同时，在轨道上，\(i, j, \cdots\) 代表占据轨道，\(a, b, \cdots\) 代表空间轨道；除此之外的大部分记号与 Yamaguchi 书中的保持基本一致．

6.2.2. 偶极矩基础回顾

在上一节中，我们已经计算了 Hartree-Fock 的偶极矩；在这里我们会从梯度的角度作回顾．我们在上一节中，考虑的 Hamiltonian 是电子与原子核之间的：

\begin{equation} \hat H^\textrm{tot} = \hat H \end{equation}

而如果现在存在关于外加物理量 \(\xi_1, \xi_2\) 的微扰，则 Hamiltonian 的形式应该更变为

\begin{equation} \hat H^\textrm{tot} = \hat H + \hat H_{\xi_1}' + \hat H_{\xi_2}' + \hat H_{\xi_1 \xi_1}'' + \hat H_{\xi_2 \xi_2}'' + \hat H_{\xi_1 \xi_2}'' + \cdots \end{equation}

现在的外加微扰是电场 \(\boldsymbol{E} = (E_x, E_y, E_z)\)，并且暂时忽略其他物理量微扰与高阶微扰，于是可以得到

\begin{equation} \hat H^\textrm{tot} = \hat H + \hat H_{\boldsymbol{E}}' = \hat F - \boldsymbol{E} \cdot \hat{\boldsymbol{\mu}} \end{equation}

其中，我们将 Fock 算符 \(\hat F\) 作为主要贡献，其余项作为微扰．其中，

\begin{equation} \hat{\boldsymbol{\mu}} = (- x, - y, - z) \end{equation}

之所以取负号，是因为电子本身是原子单位下负一电荷的．由于偶极矩算符只涉及一个电子坐标，因此它是单电子算符．

提示

只考虑外加电场与偶极的 Hamiltonian 算符实际上只有一项算符贡献，并不存在算符上的高阶微扰．我们将来会提到极化率，它是偶极矩对电场的导数，但极化率不反应在算符上，即 \(\nabla_{\boldsymbol{E}} \otimes \hat{\boldsymbol{\mu}} = \boldsymbol{0}\)；而其实反映在能量对电场的梯度上．这之后会详细提及．

提示

这里提及的偶极矩是电子云所产生的偶极；不考虑核电荷产生的偶极，也不考虑因坐标选取不同而产生的偶极大小的变化．但这不会影响极化率，因为极化率是偶极对电场的导数；所有额外产生的相对于波函数变化而呈常数的贡献项在电场的导数下便为零．这会反映在之后的极化率公式中．

6.2.3. Hartree-Fock 偶极矩

对于任意方法所产生的波函数 \(\Psi\) 而言，偶极矩的定义为

\begin{equation} \boldsymbol{\mu} = - \nabla_{\boldsymbol{E}} \langle \Psi | \hat H^\textrm{tot} | \Psi \rangle \end{equation}

对于 Hartree-Fock，其波函数所导出的偶极矩定义为

\begin{equation} \boldsymbol{\mu}^\textsf{HF} = - \nabla_{\boldsymbol{E}} \langle \Psi^\textsf{HF} | \hat H^\textrm{tot} | \Psi^\textsf{HF} \rangle \end{equation}

或者，我们说上述偶极 \(\boldsymbol{\mu}^\textsf{HF} = (\mu_x^\textsf{HF}, \mu_y^\textsf{HF}, \mu_z^\textsf{HF})\) 在 \(f\) 上的分量为 \(\mu_f^\textsf{HF}\) 为

\begin{equation} \mu_f^\textsf{HF} = - \partial_{E_f} \langle \Psi^\textsf{HF} | \hat H^\textrm{tot} | \Psi^\textsf{HF} \rangle \end{equation}

根据 Yamaguchi (17.37) 的直接推论，我们可以得到：

\begin{equation} \mu_f^\textsf{HF} = 2 \mu^f_{ii} = 2 \mu^f_{\mu \nu} P^\textsf{HF}_{\mu \nu} \end{equation}

下面的代码就是计算电子云本身所贡献的偶极矩．

# 计算偶极矩
dip_hf = 2 * np.einsum("fuv, uv -> f", dip, D)
dip_hf = dip_hf.round(decimals=6)
dip_hf

# 将结果截断到 6 位小数
dip_hf = dip_hf.round(decimals=6)
dip_hf

array([-0.      ,  0.      , -0.049403])

6.2.4. 极化率基础回顾

对于任意方法所产生的波函数 \(\Psi\) 而言，偶极矩的定义为

\begin{equation} \boldsymbol{\alpha} = \nabla_\boldsymbol{E} \otimes \boldsymbol{\mu} = \nabla_\boldsymbol{E} \otimes \nabla_\boldsymbol{E} \langle \Psi | \hat H^\textrm{tot} | \Psi \rangle \end{equation}

很自然地，Hartree-Fock 的偶极矩定义为

\begin{equation} \boldsymbol{\alpha}^\textsf{HF} = \nabla_\boldsymbol{E} \otimes \nabla_\boldsymbol{E} \langle \Psi^\textsf{HF} | \hat H^\textrm{tot} | \Psi^\textsf{HF} \rangle \end{equation}

或者，等价地，

\begin{equation} \alpha_{fg}^\textsf{HF} = \partial_{E_f} \partial_{E_g} \langle \Psi^\textsf{HF} | \hat H^\textrm{tot} | \Psi^\textsf{HF} \rangle = \partial_{E_g} \mu_f^\textsf{HF} \end{equation}

根据 Yamaguchi (17.54) 的结论，Hartree-Fock 的极化率计算公式非常简洁：

\begin{equation} \alpha_{fg}^\textsf{HF} = -4 U_{ai}^g h_{ai}^f \end{equation}

但其中 \(U_{ai}^g\) 不太容易计算，也是计算偶极矩中程序复杂、时间消耗量大的部分．它需要通过 CP-HF 方程给出．

6.2.5. CP-HF 方程

CP-HF 方程可以通过对 Fock 矩阵的导数 \(\partial_\xi F_{pq} = 0\) 确定，其非冗余的，即非占-占据部分的反对称旋转矩阵 \(U_{ai}^\xi\) 的 CP-HF 方程根据 Yamaguchi (10.21) 或 (X.1)，形式为

\begin{equation} (\varepsilon_i - \varepsilon_a) U_{ai}^\xi - A_{ij}^{ab} U_{bj}^\xi = F_{ai}^\xi - \frac{1}{2} S_{ai}^\xi \varepsilon_i - S_{kl}^\xi A_{il}^{ak} \end{equation}

其中，

\begin{equation} A_{pr}^{qs} = 4 (pq|rs) - (pr|qs) - (ps|qr) = 4 g_{pr}^{qs} - g_{pq}^{rs} - g_{pq}^{sr} \end{equation}

对于 CP-HF 方程的右方，常记为 \(B_{ai}^\xi\)．由于对于电场分量 \(\xi = E_g\) 时，由于 \(S_{pq}^g = 0\)，因此 CP-HF 方程的右侧形式大大简化：

\begin{equation} B_{ai}^g = F_{ai}^g = h_{ai}^g \end{equation}

故而电场分量 \(E_g\) 的 CP-HF 为

\begin{equation} (\varepsilon_i - \varepsilon_a) U_{ai}^g - A_{ij}^{ab} U_{bj}^g = h_{ai}^g \end{equation}

6.3. 偶极的 CP-HF 方程

6.3.1. 直接解法

我们现在关心 CP-HF 方程的左侧．等式左侧可以写为两个张量乘积的形式：

\begin{equation} - A'{}_{ij}^{ab} U_{bj}^g = h_{ai}^g \end{equation}

其中，

\begin{equation} A'{}_{ij}^{ab} = \delta_{ab} \delta_{ij} (\varepsilon_b - \varepsilon_j) + A_{ij}^{ab} \end{equation}

这就将问题化归为求解线性方程了．

从实现的角度来讲，我们可以更改角标，使得 CP-HF 方程更适合传统线性方程组问题．我们将上述方程写作

\begin{equation} - A'_{ai, bj} U_{bj}^g = h_{ai}^g \end{equation}

并且将 \(U_{bj}^g\) 与 \(h_{ai}^g\) 看作向量，把 \(A'_{ai, bj}\) 看作矩阵．这样我们就可以用比较方便的方式写出 \(U_{bj}^g\) 的解析结果：

\begin{equation} U_{bj}^g = - (A')^{-1}_{bj, ai} h_{ai}^g \end{equation}

代入偶极矩公式，可以立即得到

\begin{equation} \alpha_{fg}^\textsf{HF} = -4 U_{ai}^g h_{ai}^f = 4 h_{bj}^g (A')^{-1}_{bj, ai} h_{ai}^f = 4 \mu_{bj}^g (A')^{-1}_{bj, ai} \mu_{ai}^f \end{equation}

下面我们就计算极化率．首先我们给出 \(\mathbf{A}'\) 张量，并将其四脚标转换为两角标：

# 先给出 A 的形式，角标 pqrs
A = 4 * g_pqrs - np.einsum("prqs -> pqrs", g_pqrs) - np.einsum("psqr -> pqrs", g_pqrs)

# 随后向 A_prime 上加轨道能之差
#   对于 A'_pq^rs，以及 MO 基组下的 Fock 矩阵 Fmo (对角矩阵)
#   delta_rs delta_pq e_s = delta_pq Fmo_rs
#   注意到变量 F 在上面程序中，为 AO 基组，因此 Fmo = C.T @ F @ C
#   delta 函数用单元矩阵的生成函数 numpy.eye 就可以做
A_prime = A + np.einsum("pq, rs -> prqs", np.eye(F.shape[0]), C.T @ F @ C)
#   delta_rs delta_pq e_q = delta_rs Fmo_pq
A_prime -= np.einsum("pq, rs -> prqs", C.T @ F @ C, np.eye(F.shape[0]))

# 随后，因为我们关心的是 ijab 而非 pqrs，我们取出这个矩阵的 ijab 的部分 (A_prime[i, a, j, b])
A_prime = A_prime[:nocc, nocc:, :nocc, nocc:]

# 最后，我们将矩阵压平到二维矩阵，从而允许矩阵的求逆
# 在压平前，我们需要注意到我们需要把 A_prime[i, a, j, b] -> A_prime[ai, jb]
#   占据与非占轨道的顺序刚好反了一下，因此需要先转置矩阵
A_prime = A_prime.transpose(1,0,3,2)
A_prime = A_prime.reshape(nocc * nvir, nocc * nvir)

其次，我们给出 MO 基组下的单电子矩阵的导数 \(\boldsymbol{\mu}^f\)，并将其压到一维向量 (\(\boldsymbol{\mu}^g\) 是相同的)：

# dip 是原子轨道基组下的偶极积分，需要转为分子轨道
# 同时需要知道，我们只需要非占-占据部分即可
dip_mo = [None] * 3
for f in range(3):
    dip_mo[f] = (Cv.T @ dip[f] @ Co).ravel()
dip_mo = np.array(dip_mo)

最后简单代入 \(\alpha_{fg}^\textsf{HF} = 4 \mu_{bj}^g (A')^{-1}_{bj, ai} \mu_{ai}^f\) 即可．

alpha_hf = 4 * np.einsum("gp, pq, fq -> fg", dip_mo, np.linalg.inv(A_prime), dip_mo)
# 截断到第六位小数的偶极矩
alpha_hf.round(decimals=6)

array([[ 1.32196 ,  0.      ,  0.      ],
       [ 0.      ,  7.086627, -0.      ],
       [ 0.      , -0.      ,  6.05264 ]])

6.3.2. 自洽迭代与 DIIS 迭代解法

上面的 CP-HF 直接解法的好处是程序化比较方便，公式表达也比较直观；但缺点是计算消耗量大，涉及到矩阵求逆；同时，需要储存四脚标的 \(\mathbf{A}'\) 张量．尽管我们在直接解法里可以通过求解线性方程组问题来避免矩阵求逆，但四脚标张量的存储仍然是不可避免的；这仍然无法满足 CP-HF 实际应用的需求．因此，下面简单介绍 CP-HF 方程的迭代求解方法．

我们知道 CP-HF 方程写为

\begin{equation} (\varepsilon_i - \varepsilon_a) U_{ai}^g - A_{ij}^{ab} U_{bj}^g = h_{ai}^g \end{equation}

那么其形式可以变为

\begin{equation} U_{ai}^g = \frac{1}{\varepsilon_i - \varepsilon_a} \left( A_{ij}^{ab} U_{bj}^g + h_{ai}^g \right) = \frac{1}{\varepsilon_i - \varepsilon_a} \left( A_{ij}^{ab} U_{bj}^g - \mu_{ai}^g \right) \end{equation}

这就成为一个迭代方程了．乍一看我们仍然需要 \(A_{ij}^{ab}\)，但实际上并不需要显式地写出该变量．我们知道 \(A_{ij}^{ab}\) 是由 ERI 构成的，那么就应当可以通过 JK 引擎给出 \(A_{ij}^{ab} U_{bj}^g\) 的值．回顾 JK 积分的定义：

\begin{align} J_{\mu \nu} [D_{\mu \nu}] &= (\mu \nu | \kappa \lambda) C_{\kappa i}^\mathrm{left} C_{\lambda i}^\mathrm{right} \\ K_{\mu \nu} [D_{\mu \nu}] &= (\mu \kappa | \nu \lambda) C_{\kappa i}^\mathrm{left} C_{\lambda i}^\mathrm{right} \end{align}

同时注意到，实际上 \(D_{\mu \nu} = C_{\kappa i}^\mathrm{left} C_{\lambda i}^\mathrm{right}\) 就是求广义上的 Coulomb、Exchange 积分时所需要代入的密度；这个密度也是广义的．在当前的 CP-HF 方程中，我们可以将这个密度定为

\begin{equation} D_{\mu \nu} = C_{\mu b} U_{bj}^g C_{\nu j} \end{equation}

那么我们在调用这个引擎的时候，只需要定义左右矢为

\begin{align} C_{\kappa b}^\mathrm{left} &= C_{\kappa b} \\ C_{\lambda b}^\mathrm{right} &= U_{bj}^g C_{\lambda j} \end{align}

那么，最终 \(A_{ij}^{ab} U_{bj}^g\) 的结果则为

\begin{align} A_{ij}^{ab} U_{bj} &= (4 g_{ij}^{ab} - g_{ia}^{jb} - g_{ia}^{bj}) U_{bj}^g \\ &= C_{\mu i} C_{\nu a} \big( 4 (\mu \nu | \lambda \kappa) - (\mu \lambda | \nu \kappa) - (\mu \kappa | \nu \lambda) \big) C_{\kappa b} U_{bj}^g C_{\lambda j} \\ &= C_{\mu i} C_{\nu a} (4 J_{\mu \nu} [D_{\mu \nu}] - K_{\nu \mu} [D_{\mu \nu}] - K_{\mu \nu} [D_{\mu \nu}]) \end{align}

下面我们就进行迭代计算的程序实现．首先，我们需要初始化 JK 引擎，并给出 \(\varepsilon_i - \varepsilon_a\) 矩阵：

# JK 引擎初始化
jk = psi4.core.JK.build(scf_wfn.basisset())
jk.initialize()

# 定义左矢为 Cv，右矢为空矩阵
# 我们将来可以通过更改底层数据来改变右矢，这里只是声明而已
C_right_list = []
for g in range(3):
    jk.C_left_add(psi4.core.Matrix.from_array(Cv))
    mat = psi4.core.Matrix(nbf, nvir)
    C_right_list.append(np.asarray(mat))
    jk.C_right_add(mat)

# dvo_ai = e_i - e_a
d_vo = - ev.reshape(-1, 1) + eo

随后我们构建初猜，并建立迭代所需要的旧矩阵．初猜的构建方式即

\begin{equation} - (\varepsilon_i - \varepsilon_a)^{-1} \mu_{ai}^g = - (\varepsilon_i - \varepsilon_a)^{-1} \mu_{\mu \nu}^g C_{\mu a} C_{\nu i} \end{equation}

# 分子轨道下的偶极矩矩阵
dip_vo = np.einsum("guv, ua, vi -> gai", dip, Cv, Co, optimize=True)
# 初猜
U = - dip_vo / d_vo
# 旧矩阵
U_old = np.copy(U)

我们可以使用 DIIS 加速迭代．这使用了 Psi4NumPy 中的代码．在目前这个例子中，如果使用 DIIS 加速，则 11 次迭代后收敛；若使用普通自洽收敛，则 38 次迭代后收敛．这里不详细叙述 DIIS 迭代原理．

USE_DIIS = True
diis = []
if USE_DIIS:
    for _ in range(3):
        diis.append(DIIS_helper())

随后我们就可以进行迭代计算了．迭代的收敛方式是初猜与旧矩阵之间的差模小于一定阈值．

%%time
# 给定最大迭代次数与收敛限
MAX_ITER = 100
CONV = 1.e-8

for it in range(1, MAX_ITER + 1):
    # 首先更新 JK 引擎的右矢为 U_bj C_lj -> Cright_bl
    # 这里用到使用列表、以及使用 [:] 进行向量视图数据替换的技巧
    for g in range(3):
        C_right_list[g][:] = Co @ U[g].T

    # 随后计算 J[D]、K[D] 积分，并且将它们储存为临时变量
    # 这之后的计算都假定在同一个坐标取向 g 下进行
    jk.compute()
    for g in range(3):
        J = np.asarray(jk.J()[g])
        K = np.asarray(jk.K()[g])

        # 可以更新 U 矩阵了
        # Unew_ai^g = [ C_ui C_va (4 * JD_uv - KD_vu - KD_uv) - mu_ai^g ] / (e_i - e_a)
        U[g] = Cv.T @ (4 * J - K.T - K) @ Co - dip_vo[g]
        U[g] /= d_vo

    # 若使用了 DIIS 加速，则执行下述插值
    if USE_DIIS:
        for g in range(3):
            diis[g].add(U[g], U[g] - U_old[g])
            U[g] = diis[g].extrapolate()

    # 检查收敛情况
    rms = np.linalg.norm(U - U_old)
    # print('CPHF Iteration {:3d}: RMS = {:14.10f}'.format(it, rms))

    # 判断是否收敛
    if (rms < CONV):
        print("CPHF Converged in {:3d} iterations!".format(it))
        break
    else:
        U_old = np.copy(U)

CPHF Converged in  11 iterations!
CPU times: user 17.4 ms, sys: 275 µs, total: 17.7 ms
Wall time: 26.6 ms

最终的结果只要简单地代入 \(\alpha_{fg}^\textsf{HF} = -4 U_{ai}^g h_{ai}^f = 4 U_{ai}^g \mu_{ai}^f\) 即可．

alpha_hf = np.einsum("gai, fai -> fg", U, dip_vo) * 4
alpha_hf.round(decimals=6)

array([[ 1.32196 ,  0.      ,  0.      ],
       [ 0.      ,  7.086627, -0.      ],
       [ 0.      , -0.      ,  6.05264 ]])

6.4. MP2 偶极矩理论回顾

6.4.1. MP2 偶极矩公式

回顾上面对偶极矩的叙述，MP2 的偶极矩应该定义为

\begin{equation} \boldsymbol{\mu}^\textsf{MP2} = - \nabla_{\boldsymbol{E}} E^\textsf{MP2} = \boldsymbol{\mu}^\textsf{HF} + \boldsymbol{\mu}^{(2)} = \boldsymbol{\mu}^\textsf{HF} - \nabla_{\boldsymbol{E}} \langle \Psi^\textsf{HF} | \hat H^\textrm{tot} | \Psi^{(2)} \rangle \end{equation}

但其偶极矩不能简单地通过波函数作用在偶极算符上直接给出，而需要仔细地对 MP2 相关能进行导数计算．关于这部分的详细推导，见 Aikens, et al. 2003．在这里只列出重要的结论．

对于偶极矩，即电场导数而言，许多能量导数中的项为零；因此，MP2 相关能部分对偶极矩的贡献可以写为简单的弛豫密度 \(\mathbf{P}^\textsf{(2)}\) 与 AO 基组偶极矩矩阵的乘积：

\begin{equation} \mu_f^{(2)} = 2 P_{\mu \nu}^{(2)} \mu_{\mu \nu}^f \end{equation}

因此，弛豫密度 \(\mathbf{P}^{(2)}\) 的确定将是后面的重点．

6.4.2. 弛豫密度的拆分

现在我们所讨论的弛豫密度会在分子轨道基组下．其与原子轨道基组的转换是非常显然的：

\begin{equation} P_{\mu \nu}^{(2)} = C_{\mu p} P_{pq}^{(2)} C_{\nu q} \end{equation}

分子轨道基组分为占据与非占据的部分：

\begin{equation} P_{pq}^{(2)} = \begin{pmatrix} P_{ij}^{(2)} & P_{ia}^{(2)} \\ P_{ai}^{(2)} & P_{ab}^{(2)} \end{pmatrix} \end{equation}

其中，占据-占据部分与非占-非占部分的结果可以立即得到，但非占-占据部分的结果不是很显然．我们先把占据-占据与非占-非占部分的结果书写如下 (Aikens, eq.144-145，但这里定义的密度是 \(\alpha\) 或 \(\beta\) 自旋密度，并不是两者之和，因此会少两倍；后面的矩阵同样如此)：

\begin{align} P_{ij}^{(2)} &= - T_{ik}^{ab} t_{jk}^{ab} \\ P_{ab}^{(2)} &= T_{ij}^{ac} t_{ij}^{bc} \end{align}

其中引入了新的符号 \(T_{ij}^{ab}\)，它代表的是闭壳层下的简写张量：

\begin{equation} T_{ij}^{ab} = 2 t_{ij}^{ab} - t_{ij}^{ba} = t_{ij}^{ab} + \bar t_{ij}^{ab} \end{equation}

我们将这两个弛豫密度的部分写到代码中．

# 定义闭壳层的 T 张量
T_ovov = 2 * t_ovov - t_ovov.swapaxes(1,3)

# 占据-占据弛豫密度
P2_oo = - np.einsum("iakb, jakb -> ij", T_ovov, t_ovov)
# 非占-非占弛豫密度
P2_vv = np.einsum("iajc, ibjc -> ab", T_ovov, t_ovov)

6.4.3. 非占-占据部分弛豫密度

这部分的弛豫密度是通过 CP-HF 方程给出的 (Aikens, eq.163)：

\begin{equation} (\varepsilon_i - \varepsilon_a) P_{ai}^{(2)} - A_{ij}^{ab} P_{bj}^{(2)} = L_{ai} \end{equation}

该 CP-HF 方程与上面求 Hartree-Fock 极化率所使用的 CP-HF 方程的解法是相同的，只是其中的旋转矩阵变为了占据-非占部分弛豫密度，响应量变为了拉格朗日量．拉格朗日量的定义为 (Aikens, eq.159)

\begin{equation} L_{ai} = \frac{1}{2} P_{jk}^{(2)} A_{aj}^{ik} + \frac{1}{2} P_{bc}^{(2)} A_{ab}^{ic} - T_{jk}^{ab} g_{ib}^{jk} + T_{ij}^{bc} g_{aj}^{bc} \end{equation}

我们先写出 MP2 的拉格朗日矩阵量：

L_vo = 0.5 * np.einsum("jk, aijk -> ai", P2_oo, A[nocc:, :nocc, :nocc, :nocc])
L_vo += 0.5 * np.einsum("bc, aibc -> ai", P2_vv, A[nocc:, :nocc, nocc:, nocc:])
L_vo += - np.einsum("jakb, ijbk -> ai", T_ovov, g_pqrs[:nocc, :nocc, nocc:, :nocc])
L_vo += np.einsum("ibjc, abjc -> ai", T_ovov, g_pqrs[nocc:, nocc:, :nocc, nocc:])

6.4.4. CP-HF 方程求解与最终弛豫密度

仿照刚才的过程，我们很快就能得到 MP2 的弛豫密度的非占-占据部分，进而整个 MP2 弛豫密度：

# JK 引擎
jk = psi4.core.JK.build(scf_wfn.basisset())
jk.initialize()
jk.C_left_add(psi4.core.Matrix.from_array(Cv))
mat = psi4.core.Matrix(nbf, nvir)
C_right = np.asarray(mat)
jk.C_right_add(mat)

# 初猜与收敛设定
d_vo = - ev.reshape(-1, 1) + eo
P2_vo = L_vo / d_vo
P2_vo_old = np.copy(P2_vo)
USE_DIIS = True
diis = []
if USE_DIIS:
    diis = DIIS_helper()
MAX_ITER = 100
CONV = 1.e-8

# CP-HF 方程
for it in range(1, MAX_ITER + 1):
    C_right[:] = Co @ P2_vo.T
    jk.compute()
    J = np.asarray(jk.J()[0])
    K = np.asarray(jk.K()[0])
    P2_vo = Cv.T @ (4 * J - K.T - K) @ Co + L_vo
    P2_vo /= d_vo

    # 若使用了 DIIS 加速，则执行下述插值
    if USE_DIIS:
        diis.add(P2_vo, P2_vo - P2_vo_old)
        P2_vo = diis.extrapolate()

    # 检查收敛情况
    rms = np.linalg.norm(P2_vo - P2_vo_old)
    # print('CPHF Iteration {:3d}: RMS = {:14.10f}'.format(it, rms))

    # 判断是否收敛
    if (rms < CONV):
        print("CPHF Converged in {:3d} iterations!".format(it))
        break
    else:
        P2_vo_old = np.copy(P2_vo)

CPHF Converged in  10 iterations!

我们就最终获得了弛豫密度 \(P_{pq}^{(2)}\)．将它在代码中实现：

P2_pq = np.block([
    [P2_oo, P2_vo.T],
    [P2_vo, P2_vv]
])

其在 AO 基组下的结果就为

P2 = C @ P2_pq @ C.T

6.4.5. 自然轨道占据数

我们可以拿 MO 基组下的弛豫密度简单地作 MP2 的自然轨道．自然轨道相对于分子轨道的占据数与系数可以通过对角化同时求出：

# 定义自洽场密度
D_pq = np.zeros_like(D)
for i in range(5):
    D_pq[i, i] = 1
# 对自洽场密度与弛豫密度的和进行对角化
no_occ, no_coef = np.linalg.eigh(P2_pq + D_pq)
# 从而得到自然轨道的占据数
no_occ[::-1] * 2

array([1.999957  , 1.99015143, 1.98160379, 1.97443089, 1.97182765, 0.02646937, 0.02370363,
       0.01771884, 0.00974191, 0.00262307, 0.00142795, 0.00023002, 0.00011445])

6.4.6. MP2 相关能对偶极矩的贡献

得到弛豫密度后，简单地代入 MP2 偶极矩公式即可．

\begin{equation} \mu_f^{(2)} = 2 P_{\mu \nu}^{(2)} \mu_{\mu \nu}^f \end{equation}

# 计算偶极矩
dip_mp2 = 2 * np.einsum("fuv, uv -> f", dip, P2)
dip_mp2 = dip_mp2.round(decimals=6)
dip_mp2

# 将结果截断到 6 位小数
dip_mp2 = dip_mp2.round(decimals=6)
dip_mp2

array([-0.     , -0.     , -0.05608])

我们可以将 \(E_z\) 分量的 MP2 相关能偶极贡献与 Gaussian 的输出结果进行比对．Gaussian 输入卡 HF, MP2

psi4.compare_values(dip_mp2[2], 1.0715445 - 1.1276241, 6, 'MP2 corr Dipole')

        MP2 corr Dipole...................................................PASSED

True