5. “安全的”MP2 Hessian 与“轨道旋转”效应的消除

我们在上一节中，使用了 \(U_{pq}^\mathbb{A}\) 来表达 MP2 的 Hessian。但我们都知道，\(U_{ij}^\mathbb{A}\) 与 \(U_{ab}^\mathbb{A}\) 由于存在奇点，因此应当用“轨道旋转”过后的 \(\mathscr{U}_{ij}^\mathbb{A}\) 与 \(\mathscr{U}_{ab}^\mathbb{A}\) 替代；这个过程中，很容易出现问题，也难以用直接对一阶导数作数值导数的方法，来验证二阶导数采用“轨道旋转”过后的 U 矩阵的正确性。这一节，我们就省视这种更变的合理性，并深入地理解 MP2 的 Hessian 推导。

这一节并不作公式的推演。尽管这一节的内容有它的意义，但读者不一定有必要完整地了解整个过程。读者只需要知道，我们在这篇文档中，验证了“轨道旋转”在二阶梯度求取中的可行性，以及在实际应用过程中，注意到普通 \(U_{pq}^\mathbb{A}\) 与 \(\mathscr{U}_{pq}^\mathbb{A}\) 的区别就可以了。

未完成文档

这份文档由于需要太过于长的公式推导过程，因此关于弛豫密度的“轨道旋转”问题就没有继续进行下去。

关于“轨道旋转”问题，可能的另一种讨论思路是使用非正则 RHF 参考态的 Λ-CCSD 方程退化到 MP2 的情形。

5.1. 准备工作

%matplotlib notebook

from pyscf import gto, scf, dft, lib, grad, hessian
from pyscf.scf import cphf
import numpy as np
from functools import partial
import warnings
from matplotlib import pyplot as plt
from pyxdh.Utilities import NucCoordDerivGenerator, DipoleDerivGenerator, NumericDiff
from pyxdh.DerivOnce import GradMP2
from pyxdh.DerivTwice import HessMP2, HessSCF

np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
np.allclose = partial(np.allclose, atol=1e-6, rtol=1e-4)
np.set_printoptions(5, linewidth=150, suppress=True)
warnings.filterwarnings("ignore")

mol = gto.Mole()
mol.atom = """
O  0.0  0.0  0.0
O  0.0  0.0  1.5
H  1.0  0.0  0.0
H  0.0  0.7  1.0
"""
mol.basis = "6-31G"
mol.verbose = 0
mol.build()

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

我们将轨道“未经旋转”的计算实例用变量 gradh_nr 与 hessh_nr 表示；而“经过旋转”的计算实例用 gradh 与 hessh 表示。

gradh_nr = GradMP2({"scf_eng": scf.RHF(mol), "cphf_tol": 1e-12, "rotation": False})
hessh_nr = HessMP2({"deriv_A": gradh_nr, "deriv_B": gradh_nr, "rotation": False})

gradh = GradMP2({"scf_eng": scf.RHF(mol), "cphf_tol": 1e-12})
hessh = HessMP2({"deriv_A": gradh, "deriv_B": gradh})

一些与“轨道旋转”无关的变量定义如下：

nmo, nao, natm, nocc, nvir = gradh.nao, gradh.nao, gradh.natm, gradh.nocc, gradh.nvir
so, sv, sa = gradh.so, gradh.sv, gradh.sa
mol_slice = gradh.mol_slice
C, Co, Cv, e, eo, ev, D = gradh.C, gradh.Co, gradh.Cv, gradh.e, gradh.eo, gradh.ev, gradh.D
H_0_ao, S_0_ao, eri0_ao, F_0_ao = gradh.H_0_ao, gradh.S_0_ao, gradh.eri0_ao, gradh.F_0_ao
H_0_mo, S_0_mo, eri0_mo, F_0_mo = gradh.H_0_mo, gradh.S_0_mo, gradh.eri0_mo, gradh.F_0_mo
H_1_ao, S_1_ao, eri1_ao, F_1_ao = gradh.H_1_ao, gradh.S_1_ao, gradh.eri1_ao, gradh.F_1_ao
H_1_mo, S_1_mo, eri1_mo, F_1_mo = gradh.H_1_mo, gradh.S_1_mo, gradh.eri1_mo, gradh.F_1_mo
Ax0_Core, B_1, U_1_vo = gradh.Ax0_Core, gradh.B_1, gradh.U_1_vo
H_2_ao, S_2_ao, eri2_ao, F_2_ao = hessh.H_2_ao, hessh.S_2_ao, hessh.eri2_ao, hessh.F_2_ao
H_2_mo, S_2_mo, eri2_mo, F_2_mo = hessh.H_2_mo, hessh.S_2_mo, hessh.eri2_mo, hessh.F_2_mo
Ax1_Core = gradh.Ax1_Core

D_iajb, t_iajb, T_iajb, D_r, W_I, L =  gradh.D_iajb, gradh.t_iajb, gradh.T_iajb, gradh.D_r, gradh.W_I, gradh.L

而一阶 U 矩阵则如以前的文档一样，分为“未经旋转”的 U_1 \(U_{pq}^\mathbb{A}\) 与“经过旋转”的 U_1_nr \(\mathscr{U}_{pq}^\mathbb{A}\)。

U_1, U_1_nr = gradh.U_1, gradh_nr.U_1

最后，我们能发现，使用“经过旋转”与“未经旋转”的 U 矩阵下，MP2 Hessian 矩阵的结果是完全一致的：

np.allclose(hessh.E_2, hessh_nr.E_2)

True

5.2. MP2 Hessian 相关能贡献回顾

上一篇文档，我们提到，

\[\begin{split}\begin{align} \partial_\mathbb{A} \partial_\mathbb{B} E_\mathrm{MP2, c} &= \partial_\mathbb{B} D_{pq}^\mathrm{MP2, oo-vv} B_{pq}^\mathbb{A} + \mathtt{RHS}_{ai}^\mathbb{B} U_{ai}^\mathbb{A} + D_{pq}^\mathrm{MP2} \cdot \partial_\mathbb{B} B_{pq}^\mathbb{A} \\ &\quad + \partial_\mathbb{B} W_{pq}^\mathrm{MP2} [\mathrm{I}] \cdot S_{pq}^\mathbb{A} + W_{pq}^\mathrm{MP2} [\mathrm{I}] \cdot \partial_\mathbb{B} S_{pq}^\mathbb{A} \\ &\quad + 2 \partial_\mathbb{B} T_{ij}^{ab} \cdot (ia | jb)^\mathbb{A} + 2 T_{ij}^{ab} \cdot \partial_\mathbb{B} (ia | jb)^\mathbb{A} \end{align}\end{split}\]

E_2_MP2_Contrib = (
    # D_r * B
    + np.einsum("Bpq, Apq -> AB", hessh.pdB_D_r_oovv, B_1)
    + np.einsum("Bai, Aai -> AB", hessh.RHS_B, U_1[:, sv, so])
    + np.einsum("pq, ABpq -> AB", D_r, hessh.pdB_B_A)
    # W_I * S
    + np.einsum("Bpq, Apq -> AB", hessh.pdB_W_I, S_1_mo)
    + np.einsum("pq, ABpq -> AB", W_I, hessh.pdB_S_A_mo)
    # T * g
    + 2 * np.einsum("Biajb, Aiajb -> AB", gradh.pdA_T_iajb, eri1_mo[:, so, sv, so, sv])
    + 2 * np.einsum("iajb, ABiajb -> AB", T_iajb, hessh.pdB_pdpA_eri0_iajb)
)
np.allclose(E_2_MP2_Contrib, hessh._get_E_2_MP2_Contrib())

True

我们会发现，上述的计算使用到的是“旋转”后的 \(\mathscr{U}_{pq}^\mathbb{A}\) 矩阵，但仍然能达到与“未旋转”的 \(U_{pq}^\mathbb{A}\) 一样的效果。我们说，下述程序可以验证之：

np.allclose(hessh._get_E_2_MP2_Contrib(), hessh_nr._get_E_2_MP2_Contrib())

True

除此之外，每一个大分项贡献都是等价的。譬如说弛豫密度部分的贡献而言，

\[\partial_\mathbb{A} \partial_\mathbb{B} E_\mathrm{MP2, c} \leftarrow \partial_\mathbb{B} D_{pq}^\mathrm{MP2, oo-vv} B_{pq}^\mathbb{A} + \mathtt{RHS}_{ai}^\mathbb{B} U_{ai}^\mathbb{A} + D_{pq}^\mathrm{MP2} \cdot \partial_\mathbb{B} B_{pq}^\mathbb{A}\]

np.allclose(
    # D_r * B, rotation
    + np.einsum("Bpq, Apq -> AB", hessh.pdB_D_r_oovv, B_1)
    + np.einsum("Bai, Aai -> AB", hessh.RHS_B, U_1[:, sv, so])
    + np.einsum("pq, ABpq -> AB", D_r, hessh.pdB_B_A),
    # D_r * B, no rotation
    + np.einsum("Bpq, Apq -> AB", hessh_nr.pdB_D_r_oovv, B_1)
    + np.einsum("Bai, Aai -> AB", hessh_nr.RHS_B, U_1_nr[:, sv, so])
    + np.einsum("pq, ABpq -> AB", D_r, hessh_nr.pdB_B_A)
)

True

上述的结论对于其它的两项 (\(W_{pq}^\mathrm{MP2} [\mathrm{I}]\) 贡献项与双电子密度贡献项)。但具体落到小分项 \(\partial_\mathbb{B} D_{pq}^\mathrm{MP2, oo-vv} B_{pq}^\mathbb{A}\) 而言，则就有差异了：

np.allclose(
    # D_r * B, rotation
    + np.einsum("Bpq, Apq -> AB", hessh.pdB_D_r_oovv, B_1),
    # D_r * B, no rotation
    + np.einsum("Bpq, Apq -> AB", hessh_nr.pdB_D_r_oovv, B_1)
)

False

因此，若要验证“旋转”后的 U 矩阵确实能给出正确的二阶梯度贡献大小，仍然需要作不少细致的分析。我们后文就分别对三项贡献作推敲。

5.3. 双电子密度项“旋转不变”性质

我们先从最容易讨论的一项开始。

\[\begin{split}\begin{align} \partial_\mathbb{A} \partial_\mathbb{B} E_\mathrm{MP2, c} &\leftarrow 2 \partial_\mathbb{B} T_{ij}^{ab} \cdot (ia | jb)^\mathbb{A} + 2 T_{ij}^{ab} \cdot \partial_\mathbb{B} (ia | jb)^\mathbb{A} \\ &= \partial_\mathbb{B} t_{ij}^{ab} \cdot \big( 2 (ia | jb)^\mathbb{A} - (ib | ja)^\mathbb{A} \big) + 2 T_{ij}^{ab} \cdot \partial_\mathbb{B} (ia | jb)^\mathbb{A} \end{align}\end{split}\]

这里开始的文档可能会一改以前文档的风格。以前的文档都是从公式出发编写程序，并回头验证公式的正确性；但在这里，我们会通过正确的程序，反推出公式与我们所需要的性质。

我们首先拆分 \(\partial_\mathbb{B} t_{ij}^{ab}\) 与 \(\partial_\mathbb{B} (ia | jb)^\mathbb{A}\)，其中 \(\partial_\mathbb{B} t_{ij}^{ab}\) 仅仅展开到 \(\partial_\mathbb{B} F_{pq}\) 的程度，因为再进行下一级展开，会使得代码太过冗杂。

np.allclose(
    # + 2 * np.einsum("Biajb, Aiajb -> AB", gradh.pdA_t_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Biajb, iajb, Aiajb -> AB", eri1_mo[:, so, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bmi, majb, iajb, Aiajb -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bmj, iamb, iajb, Aiajb -> AB", U_1[:, sa, so], eri0_mo[so, sv, sa, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bma, imjb, iajb, Aiajb -> AB", U_1[:, sa, sv], eri0_mo[so, sa, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bmb, iajm, iajb, Aiajb -> AB", U_1[:, sa, sv], eri0_mo[so, sv, so, sa], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bkj, iakb, iajb, Aiajb -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bca, icjb, iajb, Aiajb -> AB", gradh.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bcb, iajc, iajb, Aiajb -> AB", gradh.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    # + 2 * np.einsum("iajb, ABiajb -> AB", hessh.T_iajb, hessh.pdB_pdpA_eri0_iajb)
    + 2 * np.einsum("iajb, ABiajb -> AB", T_iajb, eri2_mo[:, :, so, sv, so, sv])
    + 2 * np.einsum("iajb, Bmi, Amajb -> AB", T_iajb, U_1[:, sa, so], eri1_mo[:, sa, sv, so, sv])
    + 2 * np.einsum("iajb, Bmj, Aiamb -> AB", T_iajb, U_1[:, sa, so], eri1_mo[:, so, sv, sa, sv])
    + 2 * np.einsum("iajb, Bma, Aimjb -> AB", T_iajb, U_1[:, sa, sv], eri1_mo[:, so, sa, so, sv])
    + 2 * np.einsum("iajb, Bmb, Aiajm -> AB", T_iajb, U_1[:, sa, sv], eri1_mo[:, so, sv, so, sa]),
    # True value
    + 2 * np.einsum("Biajb, Aiajb -> AB", gradh_nr.pdA_T_iajb, eri1_mo[:, so, sv, so, sv])
    + 2 * np.einsum("iajb, ABiajb -> AB", T_iajb, hessh_nr.pdB_pdpA_eri0_iajb)
)

True

我们将上式中包含 \(\mathscr{U}_{mi}^\mathbb{B}\) 的项，与根据“旋转”后的 U 矩阵所作的导数 \(\partial_\mathbb{B} F_{ki}\) 的项提出，我们验证一下，这些项在使用“旋转”与“未旋转”的 U 矩阵下，结果是否相同。

np.allclose(
    # No rotation code counterpart
    + 2 * np.einsum("Bmi, majb, iajb, Aiajb -> AB", U_1_nr[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", gradh_nr.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("iajb, Bmi, Amajb -> AB", T_iajb, U_1_nr[:, sa, so], eri1_mo[:, sa, sv, so, sv]),
    # True value: last code block line 4, 8, 14
    + 2 * np.einsum("Bmi, majb, iajb, Aiajb -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("iajb, Bmi, Amajb -> AB", T_iajb, U_1[:, sa, so], eri1_mo[:, sa, sv, so, sv])
)

True

这就意味着我们可以将所有包含 \(i\) 角标的 U 矩阵导数是否正确的问题单独提出。关于其它的角标 (即 \(a, j, b\))，也可以作同样的处理；这里就不再展开。剩下的项都只包含 Skeleton 导数，因此不可能因为 U 矩阵旋转与否而产生变化。因此，我们只要讨论清楚上面代码为何给出 True 的结果，就等于完成了关于 U 矩阵旋转性质的推敲了。

我们再对 \(\partial_\mathbb{B} F_{ki}\) 作展开。

np.allclose(
    + 2 * np.einsum("Bmi, majb, iajb, Aiajb -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", F_1_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bmk, mi, kajb, iajb, Aiajb -> AB", U_1[:, sa, so], F_0_mo[sa, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bmi, mk, kajb, iajb, Aiajb -> AB", U_1[:, sa, so], F_0_mo[sa, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", Ax0_Core(so, so, sa, so)(U_1[:, sa, so]), t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("iajb, Bmi, Amajb -> AB", T_iajb, U_1[:, sa, so], eri1_mo[:, sa, sv, so, sv]),
    # True value: last code block line 4, 8, 14
    + 2 * np.einsum("Bmi, majb, iajb, Aiajb -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("iajb, Bmi, Amajb -> AB", T_iajb, U_1[:, sa, so], eri1_mo[:, sa, sv, so, sv])
)

True

我们这里指出，\(A_{ki, ml} U_{ml}^\mathbb{B} = A_{ki, ml} \mathscr{U}_{ml}^\mathbb{B}\)，若等式两边对 \(m, l\) 角标求和。

np.allclose(
    Ax0_Core(so, so, sa, so)(U_1[:, sa, so]),
    Ax0_Core(so, so, sa, so)(U_1_nr[:, sa, so])
)

True

于是，我们可以将代码拆分为使用了 \(\mathscr{U}_{ki}^\mathbb{B}\)，与 \(\mathscr{U}_{ci}^\mathbb{B}\) 或 \(U_{ml}^\mathbb{B}\) 的两类：

np.allclose(
    # unsafe if use not-rotated U
    + 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eri0_mo[so, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bik, i, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, k, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bki, iajb, Akajb -> AB", U_1[:, so, so], T_iajb, eri1_mo[:, so, sv, so, sv])
    # safe if use not-rotated U
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", F_1_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bci, cajb, iajb, Aiajb -> AB", U_1[:, sv, so], eri0_mo[sv, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", Ax0_Core(so, so, sa, so)(U_1[:, sa, so]), t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("iajb, Bci, Acajb -> AB", T_iajb, U_1[:, sv, so], eri1_mo[:, sv, sv, so, sv])
    ,
    # True value: last code block line 4, 8, 14
    + 2 * np.einsum("Bmi, majb, iajb, Aiajb -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("iajb, Bmi, Amajb -> AB", T_iajb, U_1[:, sa, so], eri1_mo[:, sa, sv, so, sv])
)

True

我们注意到上述代码中利用到了 \(\mathscr{U}_{ki}^\mathbb{B}\) 与 \(\mathscr{U}_{ik}^\mathbb{B}\)。我们不妨将 \(\mathscr{U}_{ik}^\mathbb{B}\) 用 \(- S_{ki}^\mathbb{A} - \mathscr{U}_{ki}^\mathbb{A}\) 替代。同时我们注意到方才的代码里第 6 行使用的是 \((ka|jb)^\mathbb{A}\)，我们将该行的 \(k, i\) 角标互换，更换为 \((ia|jb)^\mathbb{A}\)。

np.allclose(
    # unsafe if use not-rotated U
    + 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eri0_mo[so, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bki, i, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, k, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, Aiajb -> AB", U_1[:, so, so], T_iajb, eri1_mo[:, so, sv, so, sv])
    # safe if use not-rotated U
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", F_1_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bci, cajb, iajb, Aiajb -> AB", U_1[:, sv, so], eri0_mo[sv, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", Ax0_Core(so, so, sa, so)(U_1[:, sa, so]), t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("iajb, Bci, Acajb -> AB", T_iajb, U_1[:, sv, so], eri1_mo[:, sv, sv, so, sv])
    + 2 * np.einsum("Bki, i, kajb, iajb, Aiajb -> AB", S_1_mo[:, so, so], eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bik, kajb, Aiajb -> AB", S_1_mo[:, so, so], T_iajb, eri1_mo[:, so, sv, so, sv])
    ,
    # True value: last code block line 4, 8, 14
    + 2 * np.einsum("Bmi, majb, iajb, Aiajb -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("iajb, Bmi, Amajb -> AB", T_iajb, U_1[:, sa, so], eri1_mo[:, sa, sv, so, sv])
)

True

随后，我们发现，与 \(\mathscr{U}_{ki}^\mathbb{B}\) 有关的项的总和恰好为零：

np.abs(
    + 2 * np.einsum("Bki, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eri0_mo[so, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("Bki, i, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, k, kajb, iajb, Aiajb -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("Bki, kajb, Aiajb -> AB", U_1[:, so, so], T_iajb, eri1_mo[:, so, sv, so, sv])
).sum()

6.570593390529218e-17

我们把所有与 \(\mathscr{U}_{ki}^\mathbb{B}\) 排除，看看剩下的项，也恰好为零：

np.abs(
    + 2 * np.einsum("kajb, iajb, Aiajb -> Aki", eri0_mo[so, sv, so, sv], 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    + 2 * np.einsum("i, kajb, iajb, Aiajb -> Aki", eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("k, kajb, iajb, Aiajb -> Aki", eo, t_iajb, 1 / D_iajb, 2 * gradh.eri1_mo[:, so, sv, so, sv] - gradh.eri1_mo[:, so, sv, so, sv].swapaxes(-1, -3))
    - 2 * np.einsum("kajb, Aiajb -> Aki", T_iajb, eri1_mo[:, so, sv, so, sv])
).sum()

1.7243299317277436e-15

我们对上式要作一定程度的更变：

np.abs(
    + 2 * np.einsum("kajb, kajb, iajb, Aiajb -> Aki", T_iajb, D_iajb, 1 / D_iajb, gradh.eri1_mo[:, so, sv, so, sv])
    + 2 * np.einsum("i, kajb, iajb, Aiajb -> Aki", eo, T_iajb, 1 / D_iajb, gradh.eri1_mo[:, so, sv, so, sv])
    - 2 * np.einsum("k, kajb, iajb, Aiajb -> Aki", eo, T_iajb, 1 / D_iajb, gradh.eri1_mo[:, so, sv, so, sv])
    - 2 * np.einsum("kajb, Aiajb -> Aki", T_iajb, eri1_mo[:, so, sv, so, sv])
).sum()

1.4084491772561986e-15

np.abs(
    + 2 * np.einsum("kajb, kajb, iajb, Aiajb -> Aki", T_iajb, D_iajb, 1 / D_iajb, gradh.eri1_mo[:, so, sv, so, sv])
    + 2 * np.einsum("i, kajb, iajb, Aiajb -> Aki", eo, T_iajb, 1 / D_iajb, gradh.eri1_mo[:, so, sv, so, sv])
    - 2 * np.einsum("k, kajb, iajb, Aiajb -> Aki", eo, T_iajb, 1 / D_iajb, gradh.eri1_mo[:, so, sv, so, sv])
    - 2 * np.einsum("kajb, Aiajb -> Aki", T_iajb, eri1_mo[:, so, sv, so, sv])
).sum()

1.4084491772561986e-15

我们发现上式可以写为 (对 \(j, a, b\) 角标求和)

\[\left[ \frac{T_{kj}^{ab} D_{kj}^{ab} + (\varepsilon_i - \varepsilon_k) T_{kj}^{ab}}{D_{ij}^{ab}} - T_{kj}^{ab} \right] (ia|jb)^\mathbb{A} = \frac{T_{kj}^{ab}}{D_{ij}^{ab}} (D_{kj}^{ab} - D_{ij}^{ab} + \varepsilon_i - \varepsilon_k) (ia|jb)^\mathbb{A}\]

我们注意到 \(D_{kj}^{ab} - D_{ij}^{ab} + \varepsilon_i - \varepsilon_k = 0\)。因此，上式确实为零。

至此，我们就完成了对双电子密度项 U 矩阵“旋转不变”的推敲。

5.4. \(W_{pq}^\mathrm{MP2} [\mathrm{I}]\) 贡献项“旋转不变”性质

我们随后以下面的一个占据-占据块的例子，讨论下述贡献项的“旋转不变”性质：

\[\partial_\mathbb{A} \partial_\mathbb{B} E_\mathrm{MP2, c} \leftarrow \partial_\mathbb{B} W_{ij}^\mathrm{MP2} [\mathrm{I}] \cdot S_{ij}^\mathbb{A} + W_{ij}^\mathrm{MP2} [\mathrm{I}] \cdot \partial_\mathbb{B} S_{ij}^\mathbb{A}\]

我们首先对上式的导数作初步展开。第一项展开到 \(\partial_\mathbb{B} T_{ij}^{ab}\) 级别，第二项完全展开：

np.allclose(
    - 2 * np.einsum("Biakb, jakb, Aij -> AB", gradh.pdA_T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("iakb, Bjakb, Aij -> AB", T_iajb, gradh.pdA_eri0_mo[:, so, sv, so, sv], S_1_mo[:, so, so])
    + np.einsum("ij, ABij -> AB", W_I[so, so], S_2_mo[:, :, so, so])
    + np.einsum("Bmi, ij, Amj -> AB", U_1[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    + np.einsum("Bmj, ij, Ami -> AB", U_1[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    ,
    # True value
    + np.einsum("Bij, Aij -> AB", hessh.pdB_W_I[:, so, so], S_1_mo[:, so, so])
    + np.einsum("ij, ABij -> AB", W_I[so, so], hessh.pdB_S_A_mo[:, :, so, so])
)

True

我们进一步对第一项展开，直到 \(\partial_\mathbb{B} F_{pq}\) 级别：

np.allclose(
    - 2 * np.einsum("Biakb, iakb, jakb, Aij -> AB", eri1_mo[:, so, sv, so, sv], 1 / D_iajb, 2 * eri0_mo[so, sv, so, sv] - eri0_mo[so, sv, so, sv].swapaxes(-1, -3), S_1_mo[:, so, so])
    - 2 * np.einsum("Bmi, makb, iakb, jakb, Aij -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iamb, iakb, jakb, Aij -> AB", U_1[:, sa, so], eri0_mo[so, sv, sa, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bma, imkb, iakb, jakb, Aij -> AB", U_1[:, sa, sv], eri0_mo[so, sa, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmb, iakm, iakb, jakb, Aij -> AB", U_1[:, sa, sv], eri0_mo[so, sv, so, sa], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", gradh.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", gradh.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("iakb, Bjakb, Aij -> AB", T_iajb, eri1_mo[:, so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmj, iakb, makb, Aij -> AB", U_1[:, sa, so], T_iajb, eri0_mo[sa, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iakb, jamb, Aij -> AB", U_1[:, sa, so], T_iajb, eri0_mo[so, sv, sa, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bma, iakb, jmkb, Aij -> AB", U_1[:, sa, sv], T_iajb, eri0_mo[so, sa, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmb, iakb, jakm, Aij -> AB", U_1[:, sa, sv], T_iajb, eri0_mo[so, sv, so, sa], S_1_mo[:, so, so])
    + np.einsum("ij, ABij -> AB", W_I[so, so], S_2_mo[:, :, so, so])
    + np.einsum("Bmi, ij, Amj -> AB", U_1[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    + np.einsum("Bmj, ij, Ami -> AB", U_1[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    ,
    # True value
    + np.einsum("Bij, Aij -> AB", hessh_nr.pdB_W_I[:, so, so], S_1_mo[:, so, so])
    + np.einsum("ij, ABij -> AB", W_I[so, so], hessh_nr.pdB_S_A_mo[:, :, so, so])
)

True

5.4.1. 涉及角标 \(a, b\) 的项

我们单列出涉及角标 \(a, b\) 的项，并且将这些项在使用“旋转”的 U 矩阵与未旋转的 U 矩阵作对比：

np.allclose(
    # not rotated
    - 2 * np.einsum("Bma, imkb, iakb, jakb, Aij -> AB", U_1[:, sa, sv], eri0_mo[so, sa, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmb, iakm, iakb, jakb, Aij -> AB", U_1[:, sa, sv], eri0_mo[so, sv, so, sa], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", gradh.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", gradh.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bma, iakb, jmkb, Aij -> AB", U_1[:, sa, sv], T_iajb, eri0_mo[so, sa, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmb, iakb, jakm, Aij -> AB", U_1[:, sa, sv], T_iajb, eri0_mo[so, sv, so, sa], S_1_mo[:, so, so])
    ,
    # rotated
    - 2 * np.einsum("Bma, imkb, iakb, jakb, Aij -> AB", U_1_nr[:, sa, sv], eri0_mo[so, sa, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmb, iakm, iakb, jakb, Aij -> AB", U_1_nr[:, sa, sv], eri0_mo[so, sv, so, sa], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bma, iakb, jmkb, Aij -> AB", U_1_nr[:, sa, sv], T_iajb, eri0_mo[so, sa, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmb, iakb, jakm, Aij -> AB", U_1_nr[:, sa, sv], T_iajb, eri0_mo[so, sv, so, sa], S_1_mo[:, so, so])
)

True

随后我们将上式的 \(\partial_\mathbb{B} F_{pq}\) 项作展开，并且排除所有非占-占据与占据-非占部分的 U 矩阵贡献项 (这些 U 矩阵贡献不可能因为“旋转”而改变)。我们还同时更改了上式中第 7, 8 行的下角标，并将类如 \(t_{ik}^{cb} T_{jk}^{ab}\) 统一换成 \(T_{ik}^{cb} t_{jk}^{ab}\)：

np.allclose(
    # not rotated
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb * D_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb * D_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", F_1_mo[:, sv, sv], T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, c, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bac, a, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", Ax0_Core(sv, sv, sa, so)(U_1[:, sa, so]), T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", F_1_mo[:, sv, sv], T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, c, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bbc, b, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", Ax0_Core(sv, sv, sa, so)(U_1[:, sa, so]), T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bac, ickb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bbc, iakc, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    ,
    # rotated
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", U_1_nr[:, sv, sv], eri0_mo[so, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", U_1_nr[:, sv, sv], eri0_mo[so, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, iakb, jckb, Aij -> AB", U_1_nr[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakb, jakc, Aij -> AB", U_1_nr[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
)

True

我们将上式中与轨道“旋转”无关的项与有关的项分列，并且使 U 矩阵的 \(c\) 角标始终放在前面：

np.allclose(
    # not rotated
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb * D_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb * D_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, c, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bca, a, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, c, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bcb, b, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bca, ickb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bcb, iakc, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    # rotation unrelated
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", F_1_mo[:, sv, sv], T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", F_1_mo[:, sv, sv], T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", Ax0_Core(sv, sv, sa, so)(U_1[:, sa, so]), T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", Ax0_Core(sv, sv, sa, so)(U_1[:, sa, so]), T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bca, a, ickb, iakb, jakb, Aij -> AB", S_1_mo[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bcb, b, iakc, iakb, jakb, Aij -> AB", S_1_mo[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bca, ickb, jakb, Aij -> AB", S_1_mo[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bcb, iakc, jakb, Aij -> AB", S_1_mo[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    ,
    # rotated
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", U_1_nr[:, sv, sv], eri0_mo[so, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", U_1_nr[:, sv, sv], eri0_mo[so, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, sv, sv], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, iakb, jckb, Aij -> AB", U_1_nr[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakb, jakc, Aij -> AB", U_1_nr[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
)

True

我们发现，与轨道旋转有关的所有项的总和，为零：

np.abs(
    - 2 * np.einsum("Bca, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb * D_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb * D_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bca, c, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bca, a, ickb, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bcb, c, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bcb, b, iakc, iakb, jakb, Aij -> AB", U_1[:, sv, sv], ev, T_iajb, 1 / D_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bca, ickb, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bcb, iakc, jakb, Aij -> AB", U_1[:, sv, sv], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
).sum()

7.09647114885409e-17

关于上一个代码块为何为零，可以参考上一段对双电子密度“旋转不变”性质的讨论。

5.4.2. 涉及角标 \(i, j, k\) 的项

首先，我们将所有涉及角标 \(i, j, k\) 且与轨道“旋转”有关的项单列出来：

np.allclose(
    # not rotated
    - 2 * np.einsum("Bmi, makb, iakb, jakb, Aij -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iamb, iakb, jakb, Aij -> AB", U_1[:, sa, so], eri0_mo[so, sv, sa, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", gradh.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmj, iakb, makb, Aij -> AB", U_1[:, sa, so], T_iajb, eri0_mo[sa, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iakb, jamb, Aij -> AB", U_1[:, sa, so], T_iajb, eri0_mo[so, sv, sa, sv], S_1_mo[:, so, so])
    + np.einsum("Bmi, ij, Amj -> AB", U_1[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    + np.einsum("Bmj, ij, Ami -> AB", U_1[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    ,
    # rotated
    - 2 * np.einsum("Bmi, makb, iakb, jakb, Aij -> AB", U_1_nr[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iamb, iakb, jakb, Aij -> AB", U_1_nr[:, sa, so], eri0_mo[so, sv, sa, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmj, iakb, makb, Aij -> AB", U_1_nr[:, sa, so], T_iajb, eri0_mo[sa, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iakb, jamb, Aij -> AB", U_1_nr[:, sa, so], T_iajb, eri0_mo[so, sv, sa, sv], S_1_mo[:, so, so])
    + np.einsum("Bmi, ij, Amj -> AB", U_1_nr[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    + np.einsum("Bmj, ij, Ami -> AB", U_1_nr[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
)

True

我们将上述的 \(\partial_\mathbb{B} F_{pq}\) 作展开：

np.allclose(
    # not rotated
    - 2 * np.einsum("Bmi, makb, iakb, jakb, Aij -> AB", U_1[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iamb, iakb, jakb, Aij -> AB", U_1[:, sa, so], eri0_mo[so, sv, sa, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", F_1_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bil, i, lakb, iakb, jakb, Aij -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, l, lakb, iakb, jakb, Aij -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", Ax0_Core(so, so, sa, so)(U_1[:, sa, so]), t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", F_1_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bkl, k, ialb, iakb, jakb, Aij -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, l, ialb, iakb, jakb, Aij -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", Ax0_Core(so, so, sa, so)(U_1[:, sa, so]), t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmj, iakb, makb, Aij -> AB", U_1[:, sa, so], T_iajb, eri0_mo[sa, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iakb, jamb, Aij -> AB", U_1[:, sa, so], T_iajb, eri0_mo[so, sv, sa, sv], S_1_mo[:, so, so])
    + np.einsum("Bmi, ij, Amj -> AB", U_1[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    + np.einsum("Bmj, ij, Ami -> AB", U_1[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    ,
    # rotated
    - 2 * np.einsum("Bmi, makb, iakb, jakb, Aij -> AB", U_1_nr[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iamb, iakb, jakb, Aij -> AB", U_1_nr[:, sa, so], eri0_mo[so, sv, sa, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmj, iakb, makb, Aij -> AB", U_1_nr[:, sa, so], T_iajb, eri0_mo[sa, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iakb, jamb, Aij -> AB", U_1_nr[:, sa, so], T_iajb, eri0_mo[so, sv, sa, sv], S_1_mo[:, so, so])
    + np.einsum("Bmi, ij, Amj -> AB", U_1_nr[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    + np.einsum("Bmj, ij, Ami -> AB", U_1_nr[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
)

True

我们将所有与轨道“旋转”无关的项单列出来，并且作一些角标更换，以及利用 \(\mathscr{U}_{pq}^\mathbb{B} + \mathscr{U}_{qp}^\mathbb{B} + S_{pq}^\mathbb{B} = 0\) 的性质，在 U 矩阵中尽量让角标 \(l\) 提前，就得到：

np.allclose(
    # not rotated
    - 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", U_1[:, so, so], eri0_mo[so, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", U_1[:, so, so], eri0_mo[so, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bli, i, lakb, iakb, jakb, Aij -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, l, lakb, iakb, jakb, Aij -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Blk, k, ialb, iakb, jakb, Aij -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, l, ialb, iakb, jakb, Aij -> AB", U_1[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Blk, iakb, jalb, Aij -> AB", U_1[:, so, so], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, jakb, Aij -> AB", U_1[:, so, so], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    # rotation-unrelated
    - 2 * np.einsum("Bci, cakb, iakb, jakb, Aij -> AB", U_1[:, sv, so], eri0_mo[sv, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bck, iacb, iakb, jakb, Aij -> AB", U_1[:, sv, so], eri0_mo[so, sv, sv, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", F_1_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", F_1_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", Ax0_Core(so, so, sa, so)(U_1[:, sa, so]), t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", Ax0_Core(so, so, sa, so)(U_1[:, sa, so]), t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bcj, iakb, cakb, Aij -> AB", U_1[:, sv, so], T_iajb, eri0_mo[sv, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bck, iakb, jacb, Aij -> AB", U_1[:, sv, so], T_iajb, eri0_mo[so, sv, sv, sv], S_1_mo[:, so, so])
    + np.einsum("Bci, ij, Acj -> AB", U_1[:, sv, so], W_I[so, so], S_1_mo[:, sv, so])
    + np.einsum("Bcj, ij, Aci -> AB", U_1[:, sv, so], W_I[so, so], S_1_mo[:, sv, so])
    + 2 * np.einsum("Bjl, iakb, lakb, Aij -> AB", S_1_mo[:, so, so], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bli, i, lakb, iakb, jakb, Aij -> AB", S_1_mo[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bkl, k, ialb, iakb, jakb, Aij -> AB", S_1_mo[:, so, so], eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, jakb, Aij -> AB", S_1_mo[:, so, so], T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
    ,
    # rotated
    - 2 * np.einsum("Bmi, makb, iakb, jakb, Aij -> AB", U_1_nr[:, sa, so], eri0_mo[sa, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iamb, iakb, jakb, Aij -> AB", U_1_nr[:, sa, so], eri0_mo[so, sv, sa, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Bli, lakb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("Blk, ialb, iakb, jakb, Aij -> AB", gradh_nr.pdA_F_0_mo[:, so, so], t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("Bmj, iakb, makb, Aij -> AB", U_1_nr[:, sa, so], T_iajb, eri0_mo[sa, sv, so, sv], S_1_mo[:, so, so])
    - 2 * np.einsum("Bmk, iakb, jamb, Aij -> AB", U_1_nr[:, sa, so], T_iajb, eri0_mo[so, sv, sa, sv], S_1_mo[:, so, so])
    + np.einsum("Bmi, ij, Amj -> AB", U_1_nr[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
    + np.einsum("Bmj, ij, Ami -> AB", U_1_nr[:, sa, so], W_I[so, so], S_1_mo[:, sa, so])
)

True

随后我们就可以分离角标 \(i, k\) 分别进行讨论。其中，对于关于角标 \(i\) 的四行，我们去除其中关于 \(\mathscr{U}_{li}^\mathbb{B}\) 的部分 (即不对该矩阵进行张量缩并)，得到的仍然是零张量：

np.abs(
    - 2 * np.einsum("lakb, iakb, jakb, Aij -> Ali", eri0_mo[so, sv, so, sv], 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    - 2 * np.einsum("i, lakb, iakb, jakb, Aij -> Ali", eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("l, lakb, iakb, jakb, Aij -> Ali", eo, t_iajb, 1 / D_iajb, T_iajb * D_iajb, S_1_mo[:, so, so])
    + 2 * np.einsum("lakb, jakb, Aij -> Ali", T_iajb, eri0_mo[so, sv, so, sv], S_1_mo[:, so, so])
).sum()

6.082105753059246e-16

其为零的缘由在前文已经有所描述了。

我们就不再进行更多推导了。我们相信，对于角标 \(k\) 而言，也有相同的情况。我们方才仅仅对

\[\partial_\mathbb{A} \partial_\mathbb{B} E_\mathrm{MP2, c} \leftarrow \partial_\mathbb{B} W_{ij}^\mathrm{MP2} [\mathrm{I}] \cdot S_{ij}^\mathbb{A} + W_{ij}^\mathrm{MP2} [\mathrm{I}] \cdot \partial_\mathbb{B} S_{ij}^\mathbb{A}\]

的情况作了推导；而对于 W 矩阵在非占-占据、非占-非占的情况下，相信有着相同的结果。

5.5. 弛豫密度贡献项“旋转不变”性质

\[\partial_\mathbb{A} \partial_\mathbb{B} E_\mathrm{MP2, c} \leftarrow \partial_\mathbb{B} D_{pq}^\mathrm{MP2, oo-vv} B_{pq}^\mathbb{A} + \mathtt{RHS}_{ai}^\mathbb{B} U_{ai}^\mathbb{A} + D_{pq}^\mathrm{MP2} \cdot \partial_\mathbb{B} B_{pq}^\mathbb{A}\]

我们可以将弛豫密度的贡献项分为占据-占据、非占-非占、非占-占据贡献项三种。下面我们只对占据-占据提供讨论的思路；详细的讨论确实过于繁杂，这里就不再详述了。

\[\partial_\mathbb{A} \partial_\mathbb{B} E_\mathrm{MP2, c} \leftarrow \partial_\mathbb{B} D_{ij}^\mathrm{MP2} B_{ij}^\mathbb{A} + D_{ij}^\mathrm{MP2} \cdot \partial_\mathbb{B} B_{ij}^\mathbb{A}\]

我们首先，判断轨道旋转与未旋转的情况下，对相关能二阶梯度的贡献是否相等：

np.allclose(
    # D_r * B, rotation
    + np.einsum("Bij, Aij -> AB", hessh.pdB_D_r_oovv[:, so, so], B_1[:, so, so])
    + np.einsum("ij, ABij -> AB", D_r[so, so], hessh.pdB_B_A[:, :, so, so])
    ,
    # D_r * B, no rotation
    + np.einsum("Bij, Aij -> AB", hessh_nr.pdB_D_r_oovv[:, so, so], B_1[:, so, so])
    + np.einsum("ij, ABij -> AB", D_r[so, so], hessh_nr.pdB_B_A[:, :, so, so])
)

True

我们首先对 \(\partial_\mathbb{B} D_{ij}^\mathrm{MP2} B_{ij}^\mathbb{A}\) 与 \(\partial_\mathbb{B} B_{ij}^\mathbb{A}\) 作初步展开，得到：

np.allclose(
    # pd D_r * B, rotation
    - 2 * np.einsum("Biakb, jakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, F_1_mo[:, so, so])
    + 2 * np.einsum("Biakb, jakb, Aij, j -> AB", gradh.pdA_t_iajb, T_iajb, S_1_mo[:, so, so], eo)
    + np.einsum("Biakb, jakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, Ax0_Core(so, so, so, so)(S_1_mo[:, so, so]))
    - 2 * np.einsum("Bjakb, iakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, F_1_mo[:, so, so])
    + 2 * np.einsum("Bjakb, iakb, Aij, j -> AB", gradh.pdA_t_iajb, T_iajb, S_1_mo[:, so, so], eo)
    + np.einsum("Bjakb, iakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, Ax0_Core(so, so, so, so)(S_1_mo[:, so, so]))
    # D_r * pd B, rotation
    - 2 *  np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, hessh.pdB_F_A_mo[:, :, so, so])
    + 2 * np.einsum("iakb, jakb, ABij, j -> AB", t_iajb, T_iajb, hessh.pdB_S_A_mo[:, :, so, so], eo)
    + 2 * np.einsum("iakb, jakb, Ami, Bmj -> AB", t_iajb, T_iajb, S_1_mo[:, :, so], gradh.pdA_F_0_mo[:, :, so])
    + np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, Ax0_Core(so, so, so, so)(hessh.pdB_S_A_mo[:, :, so, so]))
    + np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, Ax1_Core(so, so, so, so)(S_1_mo[:, so, so]).swapaxes(0, 1))
    + 2 * np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, Ax0_Core(so, so, sa, so)(np.einsum("Bml, Akl -> ABmk", U_1[:, :, so], S_1_mo[:, so, so])))
    + np.einsum("iakb, jakb, Aim, Bmj -> AB", t_iajb, T_iajb, Ax0_Core(so, sa, so, so)(S_1_mo[:, so, so]), U_1[:, :, so])
    + np.einsum("iakb, jakb, Amj, Bmi -> AB", t_iajb, T_iajb, Ax0_Core(sa, so, so, so)(S_1_mo[:, so, so]), U_1[:, :, so])
    ,
    # D_r * B, no rotation
    + np.einsum("Bij, Aij -> AB", hessh_nr.pdB_D_r_oovv[:, so, so], B_1[:, so, so])
    + np.einsum("ij, ABij -> AB", D_r[so, so], hessh_nr.pdB_B_A[:, :, so, so]),
)

True

我们将其中的三种贡献项分开；我们知道

\[B_{pq}^\mathbb{A} = F_{pq}^\mathbb{A} - S_{pm}^\mathbb{A} F_{qm} - \frac{1}{2} A_{pq, kl} S_{kl}^\mathbb{A}\]

这三种贡献项我们就拆分如下：

np.allclose(
    # rotation: Fock part
    - 2 * np.einsum("Biakb, jakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, F_1_mo[:, so, so])
    - 2 * np.einsum("Bjakb, iakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, F_1_mo[:, so, so])
    - 2 * np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, hessh.pdB_F_A_mo[:, :, so, so])
    # rotation: S part
    + 2 * np.einsum("Biakb, jakb, Aij, j -> AB", gradh.pdA_t_iajb, T_iajb, S_1_mo[:, so, so], eo)
    + 2 * np.einsum("Bjakb, iakb, Aij, j -> AB", gradh.pdA_t_iajb, T_iajb, S_1_mo[:, so, so], eo)
    + 2 * np.einsum("iakb, jakb, ABij, j -> AB", t_iajb, T_iajb, hessh.pdB_S_A_mo[:, :, so, so], eo)
    + 2 * np.einsum("iakb, jakb, Ami, Bmj -> AB", t_iajb, T_iajb, S_1_mo[:, :, so], gradh.pdA_F_0_mo[:, :, so])
    # rotation: Ax0.S part, subscript ij
    + np.einsum("Biakb, jakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, Ax0_Core(so, so, so, so)(S_1_mo[:, so, so]))
    + np.einsum("Bjakb, iakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, Ax0_Core(so, so, so, so)(S_1_mo[:, so, so]))
    + np.einsum("iakb, jakb, Aim, Bmj -> AB", t_iajb, T_iajb, Ax0_Core(so, sa, so, so)(S_1_mo[:, so, so]), U_1[:, :, so])
    + np.einsum("iakb, jakb, Amj, Bmi -> AB", t_iajb, T_iajb, Ax0_Core(sa, so, so, so)(S_1_mo[:, so, so]), U_1[:, :, so])
    # rotation: Ax0.S part, subscript kl
    + np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, Ax0_Core(so, so, so, so)(hessh.pdB_S_A_mo[:, :, so, so]))
    + np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, Ax1_Core(so, so, so, so)(S_1_mo[:, so, so]).swapaxes(0, 1))
    + 2 * np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, Ax0_Core(so, so, sa, so)(np.einsum("Bml, Akl -> ABmk", U_1[:, :, so], S_1_mo[:, so, so])))
    ,
    # D_r * B, no rotation
    + np.einsum("Bij, Aij -> AB", hessh_nr.pdB_D_r_oovv[:, so, so], B_1[:, so, so])
    + np.einsum("ij, ABij -> AB", D_r[so, so], hessh_nr.pdB_B_A[:, :, so, so]),
)

True

其中，每一部分的贡献项都是“旋转不变”的。拿 Fock 矩阵贡献项举例：

np.allclose(
    - 2 * np.einsum("Biakb, jakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, F_1_mo[:, so, so])
    - 2 * np.einsum("Bjakb, iakb, Aij -> AB", gradh.pdA_t_iajb, T_iajb, F_1_mo[:, so, so])
    - 2 *  np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, hessh.pdB_F_A_mo[:, :, so, so])
    ,
    - 2 * np.einsum("Biakb, jakb, Aij -> AB", gradh_nr.pdA_t_iajb, T_iajb, F_1_mo[:, so, so])
    - 2 * np.einsum("Bjakb, iakb, Aij -> AB", gradh_nr.pdA_t_iajb, T_iajb, F_1_mo[:, so, so])
    - 2 *  np.einsum("iakb, jakb, ABij -> AB", t_iajb, T_iajb, hessh_nr.pdB_F_A_mo[:, :, so, so])
)

True