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
的表达式是
[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))
其中,
我们发现,实际上我们不需要全部轨道的 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