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.