4. MP2 计算

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

[1]:
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)
[2]:
mol = Mol_H2O2().mol

4.1. 量化软件的 MP2 计算

4.1.1. PySCF 计算

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

[3]:
scf_eng = scf.RHF(mol)
scf_eng.kernel()
[3]:
-150.5850337808368
[4]:
mp2_eng = mp.MP2(scf_eng)
mp2_eng.kernel()[0]
[4]:
-0.26901177599951515

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

[5]:
scf_eng.e_tot + mp2_eng.e_corr
[5]:
-150.8540455568363

4.1.2. Gaussian 计算

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

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

通过 FormchkInterface 给出的是总能量:

[7]:
ref_fchk.total_energy()
[7]:
-150.8540455443488

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

[8]:
ref_fchk.key_to_value("MP2 Energy") - ref_fchk.key_to_value("SCF Energy")
[8]:
-0.2690117592611898

4.1.3. pyxdh 计算

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

[9]:
config = {"scf_eng": scf_eng}
mp2h = GradMP2(config)
[10]:
mp2h.eng
[10]:
-150.8540455568363

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

[11]:
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

4.2. MP2 相关能计算

4.2.1. 变量定义

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

  • nmo 分子电子数
  • nao 原子轨道数
  • natm 原子数
  • nocc 占据轨道数
  • nvir 未占轨道数
  • so 占据轨道分割
  • sv 未占轨道分割
  • sa 全轨道分割
[12]:
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 轨道占据数
[13]:
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

4.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}\]
[14]:
np.allclose(
    np.einsum("up, vq, uvkl, kr, ls -> pqrs", C, C, eri0_ao, C, C),
    eri0_mo
)
[14]:
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}\) 相除;因此会引起许多程序上的问题。

[15]:
D_iajb = mp2h.D_iajb
t_iajb = mp2h.t_iajb
T_iajb = mp2h.T_iajb

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

[16]:
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}\) 可以通过简单的张量相乘求和给出:

[17]:
np.allclose((T_iajb * t_iajb * D_iajb).sum(), mp2_eng.e_corr)
[17]:
True