9. B2PLYP 能量计算

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

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

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

9.1.1. Gaussian 计算

输入卡formchk 结果

[3]:
ref_fchk = FormchkInterface(resource_filename("pyxdh", "Validation/gaussian/H2O2-B2PLYP-freq.fchk"))
ref_fchk.total_energy()
[3]:
-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 [Grimme-Grimme.JCP.2006.124] 的泛函定义。根据原文的公式 (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 定义为

[4]:
scf_eng = dft.RKS(mol)
scf_eng.grids = grids
scf_eng.xc = "0.53*HF + 0.47*B88, 0.73*LYP"
scf_eng.kernel()
[4]:
-151.11160929386716

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

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

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

[6]:
b2ph.eng
[6]:
-151.20399686033448
[7]:
np.allclose(b2ph.eng, ref_fchk.total_energy())
[7]:
True

9.1.3. PySCF 计算

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

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

[8]:
mp2_eng = mp.MP2(scf_eng)
mp2_eng.run()
[8]:
<pyscf.mp.mp2.MP2 at 0x7fad1b55feb8>

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

[9]:
mp2_eng.e_tot
[9]:
-151.45378546596837

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

[10]:
scf_eng.e_tot + 0.27 * mp2_eng.e_corr
[10]:
-151.20399686033448

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

[11]:
np.allclose(scf_eng.e_tot + 0.27 * mp2_eng.e_corr, ref_fchk.total_energy())
[11]:
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.