9. B2PLYP 能量计算

在这一节中，我们将讨论 B2PLYP 泛函的能量计算。这将作为框架更大的 XYG3 型泛函文档的铺垫。

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

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
grids = Mol_H2O2().gen_grids()

9.1. 量化软件的 B2PLYP 能量计算

9.1.1. Gaussian 计算

输入卡，formchk 结果

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

-151.2039968828007

9.1.2. pyxdh 计算

在 pyxdh 中，B2PLYP 泛函的处理方式与 MP2 完全相同，即将 PySCF 的自洽场计算实例代入 GradMP2 中即可，但有两个区别。其一是代入的自洽场实例由 scf.RHF 改成 dft.RKS；其二是需要引入 PT2 型 (从程序角度上看可以与 MP2 等同，就如同 \(E_\mathrm{x, exact}\) 作为精确交换能可以等同于 HF 交换能一样) 相关能 \(E_\mathrm{c, PT2}\) 的系数 \(c_\mathrm{c}\)。

我们回顾 B2PLYP [Gri06] 的泛函定义。根据原文的公式 (1)，

\[E_\mathrm{xc, B2PLYP} = (1 - a_\mathrm{x}) E_\mathrm{x, B88} + a_x E_\mathrm{x, exact} + b E_\mathrm{c, LYP} + c E_\mathrm{c, PT2}\]

其中，

\[a_\mathrm{x} = 0.53, \quad c = 0.27, \quad b = 1 - c = 0.73\]

但上述泛函并非用于自洽场计算的泛函。其自洽场计算所使用的泛函 (Self-consistent, \(E_\mathrm{xc}\)) 是

\[E_\mathrm{xc, B2PLYP, SCF} = 0.53 E_\mathrm{x, exact} + 0.47 E_\mathrm{x, B88} + 0.73 E_\mathrm{c, LYP}\]

而在能量计算过程中，再补上其 PT2 贡献部分 \(0.27 E_\mathrm{c, PT2}\)。

pyxdh 计算首先需要代入自洽场的 PySCF 类实例；这个类实例 scf_eng 定义为

scf_eng = dft.RKS(mol)
scf_eng.grids = grids
scf_eng.xc = "0.53*HF + 0.47*B88, 0.73*LYP"
scf_eng.kernel()

-151.11160929386716

显然，上述能量与 Gaussian 给出的能量不太一样，因为这是没有补上 PT2 贡献的能量。现在我们代入 GradMP2 中：

config = {
    "scf_eng": scf_eng,
    "cc": 0.27
}
b2ph = GradMP2(config)

其中，配置字典 config 不仅包含了 PySCF 的自洽场计算实例，还包含了 \(c_\mathrm{c} = 0.27\)。与 MP2 计算一样地，通过 eng 属性可以给出 B2PLYP 体系总能量：

b2ph.eng

-151.20399686033448

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

True

9.1.3. PySCF 计算

尽管说 PySCF 没有默认的双杂化泛函的计算方式，API 文档或说明文档也没有提及 PySCF 是否可以计算双杂化泛函；但若了解双杂化泛函的形式，我们应当发现 PySCF 计算双杂化泛函的方式与 MP2 几乎完全相同；有所区别之处也仅仅在于 PT2 部分的系数，以及用的是 DFT 自洽场而不是 HF 自洽场。因此，PySCF 可以计算双杂化泛函的能量 而不需要借助其他工具。不过需要指出，PySCF 目前不支持对双杂化泛函的梯度性质计算。

自洽场我们已经通过 scf_eng 给出了；剩下的是计算 PT2 部分。仿照 MP2 的计算代码，

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

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

不过给出最后的结果 不能 使用下述代码：

mp2_eng.e_tot

-151.45378546596837

上述代码给出的是 \(E_\mathrm{xc, B2PLYP, SCF} + E_\mathrm{c, PT2}\)，而不是 \(E_\mathrm{xc, B2PLYP, SCF} + 0.27 E_\mathrm{c, PT2}\)。因此，B2PLYP 总能量应当通过下述代码给出：

scf_eng.e_tot + 0.27 * mp2_eng.e_corr

-151.20399686033448

可以验证上述结果与 Gaussian 的结果相等：

np.allclose(scf_eng.e_tot + 0.27 * mp2_eng.e_corr, ref_fchk.total_energy())

True

9.2. 参考文献

[Gri06]

Stefan Grimme. Semiempirical hybrid density functional with perturbative second-order correlation. J. Chem. Phys., 124(3):034108, jan 2006. doi:10.1063/1.2148954.