10. XYG3 型密度泛函

这一节我们讨论 XYG3 型密度泛函 (XYG3 type of Double Hybrid density functional, xDH)。后续文档的目标就是推导 XYG3 型泛函的一阶梯度与二阶梯度。

[1]:

from pyscf import scf, gto, dft, mp
import numpy as np
from functools import partial
import warnings

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

warnings.filterwarnings("ignore")
np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
np.set_printoptions(5, linewidth=150, suppress=True)

[2]:

mol = Mol_H2O2().mol
grids = Mol_H2O2().gen_grids()

10.1. 量化软件的 XYG3 计算

10.1.1. Gaussian 计算

在一个内部版本的 Gaussian 中，我们可以获得 H2O2 分子的 XYG3 下，6-31G 基组、(99, 590) 格点的非冻核近似的计算结果；这个结果的输入卡与 formchk 文件也已经在 pyxdh 的库中。调取方式如下：

[3]:

ref_fchk = FormchkInterface(resource_filename("pyxdh", "Validation/gaussian/H2O2-XYG3-force.fchk"))
ref_fchk.total_energy()

[3]:

-151.1962822786802

10.1.2. PySCF 计算

为了实现 XYG3 的计算，我们需要对 B2PLYP 与非自洽泛函的计算过程有大致印象；但我们将为了方便后文，重新定义记号。

记号说明

为了记号简便，在一篇文档讨论确定的泛函时 (这篇文档讨论 XYG3 泛函)，会将泛函名称忽略不写。

XYG3 包含 PT2 与普通非自洽泛函的贡献部分，分别记作 \(c_\mathrm{c} E_\mathrm{c, PT2}\) 与 \(E_\mathrm{xc, n}\)。其中，\(c_\mathrm{c}\) 是 PT2 贡献的系数。
普通非自洽泛函 \(E_\mathrm{xc, n}\) 可以分为 HF 型交换能 \(c_\mathrm{x}^\mathrm{n} E_\mathrm{x, exact}\) 与纯 GGA 交换相关能 \(E_\mathrm{GGA, n}\)。
相应地，自洽泛函 \(E_\mathrm{xc}\) 可以分为 HF 型交换能 \(c_\mathrm{x} E_\mathrm{x, exact}\) 与纯 GGA 交换相关能 \(E_\mathrm{GGA}\)。

对于 XYG3 而言，自洽泛函 scf_eng 是 B3LYP：

[4]:

scf_eng = dft.RKS(mol)
scf_eng.conv_tol = 1e-11
scf_eng.conv_tol_grad = 1e-9
scf_eng.xc = "B3LYPg"
scf_eng.grids = grids
scf_eng.kernel()

[4]:

-151.37754356054216

其非自洽泛函的形式根据文献 [ZXG09] 式 (12)，为

\[E_\mathrm{xc, xDH} = E_\mathrm{xc, LDA} + c_1 (E_\mathrm{x, exact} - E_\mathrm{x, LDA}) + c_2 \Delta E_\mathrm{x, GGA} + c_3 (E_\mathrm{c, PT2} - E_\mathrm{c, LDA}) + c_4 \Delta E_\mathrm{c, GGA}\]

上式是一个普遍的双杂化泛函的框架。对于 XYG3 泛函本身而言，

\[c_1 = c_\mathrm{x}^\mathrm{n} = 0.8033, \quad c_2 = 0.2107, \quad c_3 = c_\mathrm{c} = 0.3211\]

因此有

\[\begin{split}\begin{align} E_\mathrm{xc, XYG3} &= 0.3211 E_\mathrm{c, PT2} + E_\mathrm{xc, n} \\ E_\mathrm{xc, n} &= 0.8033 E_\mathrm{x, exact} + E_\mathrm{GGA, n} \\ E_\mathrm{GGA, n} &= 0.2107 E_\mathrm{x, B88} - 0.0140 E_\mathrm{x, LDA} + 0.6789 E_\mathrm{c, LYP} \end{align}\end{split}\]

通过上述关系，我们可以定义与 \(E_\mathrm{xc, n}\) 有关的 nc_eng：

[5]:

nc_eng = dft.RKS(mol)
nc_eng.xc = "0.8033*HF - 0.0140*LDA + 0.2107*B88, 0.6789*LYP"
nc_eng.grids = grids

以及与 \(E_\mathrm{c, PT2}\) 有关的 mp2_eng 和 PT2 相关系数 c_c \(\mathrm{c_c}\)：

[6]:

c_c = 0.3211
mp2_eng = mp.MP2(scf_eng)
mp2_eng.run()

[6]:

<pyscf.mp.mp2.MP2 at 0x7f7f3f05d828>

我们将非自洽普通交换相关能量与 PT2 部分能量求和，得到 XYG3 能量：

[7]:

e_XYG3 = nc_eng.energy_tot(dm=scf_eng.make_rdm1()) + c_c * mp2_eng.e_corr
e_XYG3

[7]:

-151.1962818850459

我们可以与 Gaussian 计算结果进行比较：

[8]:

np.allclose(e_XYG3, ref_fchk.total_energy())

[8]:

True

10.1.3. pyxdh 计算

pyxdh 计算 XYG3 的单点能与 HF-B3LYP、B2PLYP 仍然是很类似的。但我们需要注意到，我们同时需要像 HF-B3LYP 一样给出自洽泛函部分与非自洽泛函部分；同时需要指定 PT2 相关系数 \(c_\mathrm{c} = 0.3211\)。

[9]:

config = {
    "scf_eng": scf_eng,
    "nc_eng": nc_eng,
    "cc": 0.3211
}
xyg3h = GradXDH(config)

随后就可以立即获得 XYG3 的总能量了：

[10]:

xyg3h.eng

[10]:

-151.1962818850459

显然上述的 XYG3 确实是基本正确的：

[11]:

np.allclose(xyg3h.eng, ref_fchk.total_energy())

[11]:

True

10.2. 参考文献

[ZXG09]

Y. Zhang, X. Xu, and W. A. Goddard. Doubly hybrid density functional for accurate descriptions of nonbond interactions, thermochemistry, and thermochemical kinetics. Proc. Natl. Acad. Sci. U.S.A., 106(13):4963–4968, mar 2009. doi:10.1073/pnas.0901093106.