2.5. pyxdh 使用介绍

pyxdh 库是一个迷你量化库,可以用不多的几行命令直接获得 xDH 型双杂化泛函 (目前只能闭壳层、非冻核,下略) 的一阶梯度,以及 B2-PLYP 型双杂化泛函的二阶梯度。同时,它还能提供一些便利的中间变量的输出,可以帮助我们更好地理解梯度的程序实现过程;我们以后会非常经常地使用这个工具。

这一节我们会讨论 pyxdh 库的最基本调用方法。更为细节的调用将在后文呈现。

警告

pyxdh 没有经过严格的测评,目前也没有任何同行评议。在这份警告撤销之前,请不要在正式发表的论文中使用此处的做法作为 XYG3 及其导数性质的计算方法。对于其它方法,譬如 MP2、双杂化泛函等性质,也请在正式发表论文中使用成熟的量化软件。

[64]:
import warnings
warnings.filterwarnings("ignore")

2.5.1. 示例:XYG3 核坐标梯度

[94]:
import numpy as np
from pyscf import gto, dft, lib, grad
from pyxdh.DerivOnce import GradXDH

2.5.1.1. 分子与格点构建

下面两个代码块实际上是 mol PySCF 的分子构建:

[65]:
mol = gto.Mole()
mol.atom = """
O  0.0  0.0  0.0
O  0.0  0.0  1.5
H  1.0  0.0  0.0
H  0.0  0.7  1.0
"""
mol.basis = "6-31G"
mol.verbose = 0
mol.build()
[65]:
<pyscf.gto.mole.Mole at 0x7fc926372278>

我们打算进行格点积分;如果要使用 grids 自定义的格点,则在 PySCF 中可作如下的定义:

[66]:
grids = dft.Grids(mol)
grids.atom_grid = (99, 590)
grids.build()
[66]:
<pyscf.dft.gen_grid.Grids at 0x7fc92631ec50>

2.5.1.2. B3LYP 自洽泛函与 XYG3 非自洽泛函的定义

XYG3 的自洽场泛函,或者说为 XYG3 泛函提供密度与轨道的参考态是 B3LYP 泛函。B3LYP 泛函的能量求取类 scf_eng 在 PySCF 可以写为

[67]:
scf_eng = dft.RKS(mol)
scf_eng.xc = "B3LYPg"
scf_eng.grids = grids

XYG3 泛函本身是由下述项构成:

\[E_\mathrm{xc}^\mathrm{XYG3} = 0.8033 E_\mathrm{x}^\mathrm{exact} - 0.0140 E_\mathrm{x}^\mathrm{LDA} + 0.2107 E_\mathrm{x}^\mathrm{B88} + 0.6789 E_\mathrm{c}^\mathrm{LYP} + 0.3211 E_\mathrm{c}^\mathrm{PT2}\]

泛函分为两部分,上式的最后一项接近于 MP2;而前四项是一般的杂化泛函贡献项,因此前四项可以归于普通的泛函类 nc_eng

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

2.5.1.3. XYG3 能量与泛函梯度计算

所有 pyxdh 包的类都由配置字典 config 实例化。XYG3 型泛函可以通过三个变量给出定义:scf_eng 自洽场泛函、nc_eng 不包含 PT2 的非自洽泛函、以及 PT2 系数 (对于 XYG3 而言为 0.3211)。代入 pyxdh.DerivOnce.GradXDH,便得到可以计算能量与梯度的 XYG3 的类 grad_xDH

[69]:
config = {
    "scf_eng": scf_eng,
    "nc_eng": nc_eng,
    "cc": 0.3211
}
grad_xDH = GradXDH(config)

XYG3 的单点能 (原子单位) 可以通过如下方式生成:

[70]:
grad_xDH.eng
[70]:
-151.1962818434804

而梯度 (原子单位) 则可以通过下述方式生成:

[71]:
grad_xDH.E_1
[71]:
array([[-0.03967538,  0.06717703,  0.14149365],
       [ 0.00876854,  0.15758362, -0.17123915],
       [ 0.01226317,  0.01305055,  0.03179645],
       [ 0.01864365, -0.23781121, -0.00205102]])

2.5.2. 示例:B2PLYP 型泛函极化率

2.5.2.1. B2PLYP 型泛函极化率的计算

下面一个例子是计算 B2PLYP 的极化率。为了计算极化率,我们先引入以下两个类:

[31]:
from pyxdh.DerivOnce import DipoleMP2
from pyxdh.DerivTwice import PolarMP2

分子 mol 与格点 grids 选取都与上面 XYG3 的梯度计算相同:

[32]:
b2plyp_eng = dft.RKS(mol)
b2plyp_eng.xc = "0.53*HF + 0.47*B88, 0.73*LYP"
b2plyp_eng.grids = grids

我们首先给出 B2PLYP 的偶极矩助手 b2plyp_dip_helper 类;其中,cc = 0.27 是 B2PLYP 中 PT2 部分的系数。偶极矩的计算类似于梯度,可以通过 E_1 给出:

[33]:
config = {
    "scf_eng": b2plyp_eng,
    "cc": 0.27
}
b2plyp_dip_helper = DipoleMP2(config)
b2plyp_dip_helper.E_1
[33]:
array([ 0.83235837,  0.60533609, -0.34817672])

通过上述定义的偶极矩助手 b2plyp_dip_helper,我们可以给出极化率助手 b2plyp_polar_helper

[34]:
config = {
    "deriv_A": b2plyp_dip_helper,
    "deriv_B": b2plyp_dip_helper,
}
b2plyp_polar_helper = PolarMP2(config)

极化率助手所给出的极化率可以通过 E_2 给出;需要注意对于极化率,需要对 E_2 给出的值取负值才是通常的极化率结果:

[35]:
- b2plyp_polar_helper.E_2
[35]:
array([[ 6.89985042, -0.11067359, -1.07620158],
       [-0.11067359,  4.74839861,  0.25707374],
       [-1.07620158,  0.25707374, 14.38297718]])

2.5.2.2. pyxdh 所提供的标准结果

pyxdh 中,B2PLYP 的极化率是 pyxdh.DerivTwice.polar_mp2 文件的其中一个 pytest 测试例。这里我们简单了解一下文件的使用方式。我们首先引入以下两个对象:

[17]:
from pkg_resources import resource_filename
from pyxdh.Utilities import FormchkInterface

我们已经预先通过 Gaussian 计算了双氧水分子 B2PLYP 的频率结果;这个结果通过 b2plyp_formchk 储存下来:

[40]:
b2plyp_formchk = FormchkInterface(resource_filename("pyxdh", "Validation/gaussian/H2O2-B2PLYP-freq.fchk"))

通过调用其中的方法,就可以获得偶极、极化率等必要的信息:

[41]:
b2plyp_formchk.dipole()
[41]:
array([ 0.83235859,  0.60533611, -0.34817773])
[42]:
b2plyp_formchk.polarizability()
[42]:
array([[ 6.89984471, -0.11067149, -1.07619714],
       [-0.11067149,  4.74839444,  0.25707124],
       [-1.07619714,  0.25707124, 14.3829714 ]])

最后,我们可以通过 np.allclose 判断我们计算的偶极、极化率是否与 Gaussian 给出的结果大致一致:

[43]:
np.allclose(
    b2plyp_dip_helper.E_1, b2plyp_formchk.dipole(),
    atol=1e-6, rtol=1e-4
)
[43]:
True
[44]:
np.allclose(
    - b2plyp_polar_helper.E_2, b2plyp_formchk.polarizability(),
    atol=1e-6, rtol=1e-4
)
[44]:
True