6. `pyxdh` 使用介绍

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

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

警告

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

[1]:

import warnings
warnings.filterwarnings("ignore")

6.1. 示例：XYG3 核坐标梯度

[2]:

import numpy as np
from pyscf import gto, dft, lib, grad
from pyxdh.DerivOnce import GradXDH

6.1.1. 分子与格点构建

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

[3]:

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()

[3]:

<pyscf.gto.mole.Mole at 0x7f89241f2be0>

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

[4]:

grids = dft.Grids(mol)
grids.atom_grid = (99, 590)
grids.build()

[4]:

<pyscf.dft.gen_grid.Grids at 0x7f89586269b0>

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

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

[5]:

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：

[6]:

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

6.1.3. XYG3 能量与泛函梯度计算

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

[7]:

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

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

[8]:

grad_xDH.eng

[8]:

-151.1962818434803

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

[9]:

grad_xDH.E_1

[9]:

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]])

6.2. 示例：B2PLYP 型泛函极化率

6.2.1. B2PLYP 型泛函极化率的计算

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

[10]:

from pyxdh.DerivOnce import DipoleMP2
from pyxdh.DerivTwice import PolarMP2

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

[11]:

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 给出：

[12]:

config = {
    "scf_eng": b2plyp_eng,
    "cc": 0.27
}
b2plyp_dip_helper = DipoleMP2(config)
b2plyp_dip_helper.E_1

[12]:

array([ 0.83235837,  0.60533609, -0.34817672])

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

[13]:

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

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

[14]:

- b2plyp_polar_helper.E_2

[14]:

array([[ 6.89985042, -0.11067359, -1.07620158],
       [-0.11067359,  4.74839861,  0.25707374],
       [-1.07620158,  0.25707374, 14.38297718]])

6.2.2. `pyxdh` 所提供的标准结果

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

[15]:

from pkg_resources import resource_filename
from pyxdh.Utilities import FormchkInterface

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

[16]:

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

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

[17]:

b2plyp_formchk.dipole()

[17]:

array([ 0.83235859,  0.60533611, -0.34817773])

[18]:

b2plyp_formchk.polarizability()

[18]:

array([[ 6.89984471, -0.11067149, -1.07619714],
       [-0.11067149,  4.74839444,  0.25707124],
       [-1.07619714,  0.25707124, 14.3829714 ]])

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

[19]:

np.allclose(
    b2plyp_dip_helper.E_1, b2plyp_formchk.dipole(),
    atol=1e-6, rtol=1e-4
)

[19]:

True

[20]:

np.allclose(
    - b2plyp_polar_helper.E_2, b2plyp_formchk.polarizability(),
    atol=1e-6, rtol=1e-4
)

[20]:

True

6.3. 示例：XYG3 型泛函极化率

为了计算 XYG3 型泛函极化率，我们先引入以下两个类：

[21]:

from pyxdh.DerivOnce import DipoleXDH
from pyxdh.DerivTwice import PolarXDH

首先，我们给出 XYG3 型泛函自洽场与非自洽场泛函部分的 PySCF 能量计算类 scf_eng, nc_eng：

[22]:

scf_eng = dft.RKS(mol)
scf_eng.xc = "B3LYPg"
scf_eng.grids = grids

[23]:

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

与梯度类的构造相似地，我们给出 XYG3 的偶极矩类 dip_xDH：

[24]:

config = {
    "scf_eng": scf_eng,
    "nc_eng": nc_eng,
    "cc": 0.3211
}
dip_xDH = DipoleXDH(config)

与 B2PLYP 的极化率构造相似地，我们给出 XYG3 的极化率类 polar_xDH：

[25]:

config = {
    "deriv_A": dip_xDH,
    "deriv_B": dip_xDH,
}
polar_xDH = PolarXDH(config)

XYG3 的极化率也可以通过 E_2 属性的负值给出：

[26]:

- polar_xDH.E_2

[26]:

array([[ 6.87997982, -0.1021484 , -1.09976624],
       [-0.1021484 ,  4.7171979 ,  0.29678172],
       [-1.09976624,  0.29678172, 14.75690205]])

这就获得了 XYG3 的极化率了。最后我们可以通过预置在 pyxdh 库中，通过数值梯度给出的 XYG3 极化率参考值 ref_polar，来验证我们所计算的 - polar_xDH.E_2 确实是正确的极化率值：

[27]:

import pickle
with open(resource_filename("pyxdh", "Validation/numerical_deriv/xdh_polarizability_xyg3.dat"), "rb") as f:
    ref_polar = pickle.load(f)["polarizability"]

[28]:

np.allclose(- polar_xDH.E_2, ref_polar, atol=1e-7, rtol=1e-5)

[28]:

True

6. pyxdh 使用介绍