3.3. MP2 计算

XYG3 型泛函从实现上，可以看作是将普通杂化泛函上引入一部分 MP2 的能量；因此这里对 MP2 的实现过程作初步的了解。

import numpy as np
from pyscf import scf, gto, mp

from pkg_resources import resource_filename
from pyxdh.Utilities import FormchkInterface
from pyxdh.Utilities.test_molecules import Mol_H2O2
from pyxdh.DerivOnce import GradMP2

from functools import partial
np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])

import warnings
warnings.filterwarnings("ignore")

np.set_printoptions(5, linewidth=150, suppress=True)

mol = Mol_H2O2().mol

3.3.1. 量化软件的 MP2 计算

3.3.1.1. PySCF 计算

PySCF 的 MP2 计算可以在给出 SCF 类 scf_eng 的前提下，通过下述代码给出 mp2_eng 以实现：

scf_eng = scf.RHF(mol)
scf_eng.kernel()

-150.5850337808368

mp2_eng = mp.MP2(scf_eng)
mp2_eng.kernel()[0]

-0.26901177599951515

上述的输出是体系的相关矫正能 \(E_\mathrm{MP2, c}\)；两者相加则可以得到总能量 \(E_\mathrm{MP2} = E_\mathrm{SCF} + E_\mathrm{MP2, c}\)：

scf_eng.e_tot + mp2_eng.e_corr

-150.8540455568363

3.3.1.2. Gaussian 计算

我们可以将上述计算结果与 Gaussian 结果 (输入卡，formchk 结果) 进行比对：

ref_fchk = FormchkInterface(resource_filename("pyxdh", "Validation/gaussian/H2O2-MP2-freq.fchk"))

通过 FormchkInterface 给出的是总能量：

ref_fchk.total_energy()

-150.8540455443488

若要获得其中的相关能 \(E_\mathrm{MP2, c}\)，需要直接读取其中的文件词条：

ref_fchk.key_to_value("MP2 Energy") - ref_fchk.key_to_value("SCF Energy")

-0.2690117592611898

3.3.1.3. pyxdh 计算

pyxdh 也可以进行 MP2 计算，但需要使用 DerivOnceMP2 的子类：

config = {"scf_eng": scf_eng}
mp2h = GradMP2(config)

mp2h.eng

-150.8540455568363

GradMP2 继承于 GradSCF，因此 GradSCF 类的所有属性都被 GradMP2 继承。

nmo = nao = mp2h.nmo
natm = mp2h.natm
nocc, nvir = mp2h.nocc, mp2h.nvir
so, sv, sa = mp2h.so, mp2h.sv, mp2h.sa

C, Co, Cv = mp2h.C, mp2h.Co, mp2h.Cv

3.3.2. MP2 相关能计算

3.3.2.1. 变量定义

在今后的文档中，我们会经常地使用其中的一些矩阵。在这里我们列举出以后程序中会常用到的变量名称和意义。

• nmo 分子电子数
• nao 原子轨道数
• natm 原子数
• nocc 占据轨道数
• nvir 未占轨道数
• so 占据轨道分割
• sv 未占轨道分割
• sa 全轨道分割

nmo = nao = mp2h.nmo
natm = mp2h.natm
nocc, nvir = mp2h.nocc, mp2h.nvir
so, sv, sa = mp2h.so, mp2h.sv, mp2h.sa

• C 系数矩阵 \(C_{\mu p}\)
• e 轨道能量 \(\varepsilon_p\)
• Co 占据轨道系数矩阵 \(C_{\mu i}\)
• Cv 未占轨道系数矩阵 \(C_{\mu a}\)
• eo 占据轨道能量 \(\varepsilon_i\)
• ev 未占轨道能量 \(\varepsilon_a\)
• D 密度矩阵 \(D_{\mu \nu}\)
• F_0_ao AO 基组 Fock 矩阵 \(F_{\mu \nu}\)
• F_0_mo MO 基组 Fock 矩阵 \(F_{pq}\) (为对角阵)
• H_0_ao AO 基组 Hamiltonian Core 矩阵 \(h_{\mu \nu}\)
• H_0_mo MO 基组 Hamiltonian Core 矩阵 \(h_{pq}\)
• eri0_ao AO 基组双电子互斥积分 \((\mu \nu | \kappa \lambda)\)
• eri0_mo MO 基组双电子互斥积分 \((pq | rs)\)
• mo_occ 轨道占据数

C, Co, Cv = mp2h.C, mp2h.Co, mp2h.Cv
e, eo, ev = mp2h.e, mp2h.eo, mp2h.ev
D = mp2h.D
F_0_ao, F_0_mo = mp2h.F_0_ao, mp2h.F_0_mo
H_0_ao, H_0_mo = mp2h.H_0_ao, mp2h.H_0_mo
eri0_ao, eri0_mo = mp2h.eri0_ao, mp2h.eri0_mo
mo_occ = mp2h.mo_occ

3.3.2.2. 实际 MP2 计算

事实上刚刚我们已经把 MP2 中计算量最大的部分，即 MO 基组的原子轨道 eri0_mo 已经生成出来了。我们回顾一下该矩阵是如何生成的。AO 原子轨道 eri0_ao 为四维张量 \((\mu \nu | \kappa \lambda)\)；那么 MO 原子轨道 eri0_mo 的表达式是

\[(pq|rs) = C_{\mu p} C_{\nu q} (\mu \nu | \kappa \lambda) C_{\kappa r} C_{\lambda s}\]

np.allclose(
    np.einsum("up, vq, uvkl, kr, ls -> pqrs", C, C, eri0_ao, C, C),
    eri0_mo
)

True

在 RHF 下，MP2 计算表现为 (Szabo, (6.74))

\[E_\mathrm{MP2, c} = \frac{(ia|jb) \big( 2 (ia|jb) - (ib|ja) \big)}{\varepsilon_i - \varepsilon_a + \varepsilon_j - \varepsilon_b} = T_{ij}^{ab} t_{ij}^{ab} D_{ij}^{ab}\]

其中，

\[\begin{split}\begin{align} D_{ij}^{ab} &= \varepsilon_i - \varepsilon_a + \varepsilon_j - \varepsilon_b \\ t_{ij}^{ab} &= \frac{(ia|jb)}{D_{ij}^{ab}} \\ T_{ij}^{ab} &= 2 t_{ij}^{ab} - t_{ij}^{ba} \end{align}\end{split}\]

我们发现，实际上我们不需要全部轨道的 MO 基组张量，只需要其中的 (占据, 非占, 占据, 非占) 部分；因此，我们定义下述的张量：

• D_iajb \(D_{ij}^{ab}\)
• t_iajb \(t_{ij}^{ab}\)
• T_iajb \(T_{ij}^{ab}\)

这些张量不能定义在全分子轨道下，因为如果推广 \(D_{ij}^{ab}\) 到 \(D_{pq}^{rs}\)，那么遇到类似于 \(D_{ij}^{ij} = 0\) 的情形，\(t_{ij}^{ij}\) 则表示为非零的 \((ii|jj)\) 与零值的 \(D_{ij}^{ij}\) 相除；因此会引起许多程序上的问题。

D_iajb = mp2h.D_iajb
t_iajb = mp2h.t_iajb
T_iajb = mp2h.T_iajb

我们可以验证从 pyxdh 给出的这些变量与上述公式给出的结果是相同的：

print(np.allclose(eo[:, None, None, None] - ev[None, :, None, None] + eo[None, None, :, None] - ev[None, None, None, :], D_iajb))
print(np.allclose(np.einsum("ui, va, uvkl, kj, lb -> iajb", Co, Cv, eri0_ao, Co, Cv) / D_iajb, t_iajb))
print(np.allclose(2 * t_iajb - t_iajb.swapaxes(1, 3), T_iajb))

True
True
True

在以后，我们还可能为了公式表达便利，会使用 \(g_{pq}^{rs}\) 表示 \((pq|rs)\)；以及 \(G_{pq}^{rs}\) 表示 \(2 g_{pq}^{rs} - g_{pq}^{sr}\)。

通过上述对变量的定义，相关能 \(E_\mathrm{MP2, c}\) 可以通过简单的张量相乘求和给出：

np.allclose((T_iajb * t_iajb * D_iajb).sum(), mp2_eng.e_corr)

True