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 泛函本身是由下述项构成:
泛函分为两部分,上式的最后一项接近于 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