5. RHF 原子核坐标二阶梯度

前面我们提及,对 RHF 的能量作一阶梯度,可以求出一些分子的性质。譬如,原子核坐标的一阶梯度 \(E^{A_t}\),可以得到分子自身结构所产生的张力 (分子力);对参考电荷所产生的电场的导数,就能得到分子偶极矩。从这一节开始,我们简单地讨论能量量的二阶梯度计算。

我们首先先是比较易于求导数的计算原子核坐标二阶导数 \(E^{A_t B_s}\)。核坐标二阶导数的最重要的意义在于,对于稳定构象分子而言,可以求取其分子频率。

提醒

我们下面使用以后经常使用的非对称双氧水分子。但该分子并非处于稳定构象,因此我们后文计算所得的分子频率并非是有物理意义的。后文所给的计算过程仅仅是演示而已。

[1]:
from pyscf import gto, scf, lib, hessian
import numpy as np
from pyxdh.Utilities import FormchkInterface, NucCoordDerivGenerator, NumericDiff
from pyxdh.Utilities.test_molecules import Mol_H2O2
from pyxdh.DerivOnce import GradSCF
from pyxdh.DerivTwice import HessSCF
import warnings

warnings.filterwarnings("ignore")
np.set_printoptions(5, linewidth=150, suppress=True)

5.1. 量化程序的频率计算

5.1.1. Gaussian 的频率分析

我们可以写如下的 输入卡,并得到 输出文件fch 文件

[2]:
with open("assets/H2O2-freq.gjf", "r") as f:
    print(f.read())
%chk=H2O2-freq
#p RHF/6-31G Freq NoSymm

H2O2 Frequency Analysis

0 1
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

通过输出文件,我们可以得到如下与频率分析 (或高阶导数) 有关的量:

[3]:
fchk = FormchkInterface("assets/H2O2-freq.fch")
  • Hessian 矩阵,维度为 \((3 n_\mathrm{atom}, 3 n_\mathrm{atom})\),可以用于计算分子频率,单位为 a.u.:
[4]:
fchk.hessian()
[4]:
array([[ 0.36765, -0.01096, -0.02986, -0.02036,  0.0064 ,  0.03848, -0.4029 ,  0.00214, -0.02579,  0.0556 ,  0.00242,  0.01717],
       [-0.01096,  0.02901,  0.11159,  0.0047 ,  0.08453, -0.11579,  0.00851, -0.03718,  0.0048 , -0.00226, -0.07637, -0.0006 ],
       [-0.02986,  0.11159,  0.47024, -0.00243,  0.00961, -0.33099,  0.02687,  0.0038 , -0.03383,  0.00542, -0.125  , -0.10542],
       [-0.02036,  0.0047 , -0.00243, -0.07793, -0.00283, -0.04145, -0.00102, -0.0056 ,  0.04462,  0.09931,  0.00372, -0.00074],
       [ 0.0064 ,  0.08453,  0.00961, -0.00283,  0.66306, -0.43734,  0.00034,  0.00409, -0.00816, -0.00392, -0.75168,  0.43588],
       [ 0.03848, -0.11579, -0.33099, -0.04145, -0.43734,  0.426  , -0.00318,  0.00431, -0.04919,  0.00616,  0.54882, -0.04582],
       [-0.4029 ,  0.00851,  0.02687, -0.00102,  0.00034, -0.00318,  0.41067, -0.01219, -0.02918, -0.00675,  0.00333,  0.00549],
       [ 0.00214, -0.03718,  0.0038 , -0.0056 ,  0.00409,  0.00431, -0.01219,  0.02907,  0.00724,  0.01565,  0.00402, -0.01535],
       [-0.02579,  0.0048 , -0.03383,  0.04462, -0.00816, -0.04919, -0.02918,  0.00724,  0.0954 ,  0.01035, -0.00389, -0.01238],
       [ 0.0556 , -0.00226,  0.00542,  0.09931, -0.00392,  0.00616, -0.00675,  0.01565,  0.01035, -0.14815, -0.00947, -0.02193],
       [ 0.00242, -0.07637, -0.125  ,  0.00372, -0.75168,  0.54882,  0.00333,  0.00402, -0.00389, -0.00947,  0.82402, -0.41993],
       [ 0.01717, -0.0006 , -0.10542, -0.00074,  0.43588, -0.04582,  0.00549, -0.01535, -0.01238, -0.02193, -0.41993,  0.16362]])

Hessian 矩阵 (或者张量) 其实就是分子能量对所有原子核坐标的三个分量的二次导数 \(E^{A_t B_s} = \frac{\partial^2 E}{\partial A_t \partial B_s}\)。上述矩阵的第一维度表示 \(A_t\),第二维度表示 \(B_s\)

  • 分子频率,对于非线性分子而言是 \(3 n_\mathrm{atom} - 6\) 维度;但该值只能从 out 文件得到而不能从 fch 文件给出,单位为 1/cm:
[5]:
with open("assets/H2O2-freq.out", "r") as f:
    for idx, line in enumerate(f.readlines()):
        if "Frequencies" in line:
            print("line {:4d}:".format(idx + 1), line[:-1])
line  574:  Frequencies -- -1580.6089             -1218.3809              1370.6206
line  588:  Frequencies --  1647.5426              3389.8666              5347.9015
  • 偶极矩的核坐标导数,维度为 \((3 n_\mathrm{atom}, 3)\),可以用于计算红外光谱峰强度 (单位为 a.u.):
[6]:
fchk.dipolederiv()
[6]:
array([[-0.2343 ,  0.01785,  0.16617],
       [ 0.07423, -0.4948 ,  0.00587],
       [ 0.04888, -0.08356, -0.50397],
       [-0.41785,  0.01883, -0.16946],
       [ 0.00548, -0.32846, -0.15515],
       [ 0.02287,  0.03304, -0.09475],
       [ 0.21503,  0.00182, -0.03122],
       [-0.04486,  0.44397,  0.02416],
       [-0.03371, -0.01025,  0.25124],
       [ 0.43712, -0.0385 ,  0.03451],
       [-0.03485,  0.37929,  0.12513],
       [-0.03804,  0.06077,  0.34748]])
  • 红外光谱强度,对于非线性分子而言是 \(3 n_\mathrm{atom} - 6\) 维度,与分子频率一一对应,单位为 km/mol (千米每摩尔):
[7]:
with open("assets/H2O2-freq.out", "r") as f:
    for idx, line in enumerate(f.readlines()):
        if "IR Inten" in line:
            print("line {:4d}:".format(idx + 1), line[:-1])
line  577:  IR Inten    --   195.2349               105.4141                99.7388
line  591:  IR Inten    --    17.5360                47.6775               105.1116
  • 极化率,维度为 \((3, 3)\),单位为 a.u.:
[8]:
fchk.polarizability()
[8]:
array([[ 6.58142, -0.0841 , -1.45378],
       [-0.0841 ,  4.26836,  0.39969],
       [-1.45378,  0.39969, 17.89033]])

上面这五个导出量中,

  • Hessian 矩阵、分子频率是能量的二阶核坐标导数量的导出结果;
  • 偶极矩的核坐标导数、红外光谱强度是能量的一阶核坐标与一阶电场到数量的导出结果;
  • 极化率是能量的二阶核坐标导数量的导出结果。

我们将会分为三篇文档来介绍这三种类型的导出量。这篇文档,我们会具体地给出 Hessian 矩阵的计算,并且借助外部程序计算分子频率。

5.1.2. PySCF 计算 RHF Hessian 矩阵

我们首先定义自洽场计算实例 scf_eng。由于非对称双氧水分子在测评和文档中都经常使用,我们可以很方便地活用 Mol_H2O2 的代码生成 RHF 类 scf.RHF 实例:

[9]:
molh = Mol_H2O2()
mol = molh.mol
scf_eng = molh.hf_eng.run()
scf_eng.e_tot
[9]:
-150.58503378083688

为了后文的便利,我们补充定义原子数 natm \(n_\mathrm{atom}\) 与 Hessian 作为矩阵时的大小 dhess \(3 n_\mathrm{atom}\)

[10]:
natm = mol.natm
dhess = natm * 3

PySCF 中,Hessian 的计算可以通过如下代码实现:

[11]:
scf_hess = hessian.RHF(scf_eng).run()

Hessian 储存在 de 变量中,其维度并非是我们常用的 \((A, t, B, s)\)\((4, 3, 4, 3)\) 的大小,而是 \((A, B, t, s)\)\((4, 4, 3, 3)\) 的大小:

[12]:
scf_hess.de.shape
[12]:
(4, 4, 3, 3)

如果我们要将 PySCF 的 Hessian 能与 Gaussian 的核对是否一致,我们需要将 Hessian 张量的中间两个维度转置:

[13]:
scf_hess.de.swapaxes(1, 2).reshape(dhess, dhess)
[13]:
array([[ 0.36765, -0.01096, -0.02986, -0.02036,  0.0064 ,  0.03848, -0.4029 ,  0.00214, -0.02579,  0.0556 ,  0.00242,  0.01717],
       [-0.01096,  0.02901,  0.11159,  0.0047 ,  0.08453, -0.11579,  0.00851, -0.03718,  0.0048 , -0.00226, -0.07637, -0.0006 ],
       [-0.02986,  0.11159,  0.47024, -0.00243,  0.00961, -0.33099,  0.02687,  0.0038 , -0.03383,  0.00542, -0.125  , -0.10542],
       [-0.02036,  0.0047 , -0.00243, -0.07793, -0.00283, -0.04145, -0.00102, -0.0056 ,  0.04462,  0.09931,  0.00372, -0.00074],
       [ 0.0064 ,  0.08453,  0.00961, -0.00283,  0.66306, -0.43734,  0.00034,  0.00409, -0.00816, -0.00392, -0.75168,  0.43588],
       [ 0.03848, -0.11579, -0.33099, -0.04145, -0.43734,  0.426  , -0.00318,  0.00431, -0.04919,  0.00616,  0.54882, -0.04582],
       [-0.4029 ,  0.00851,  0.02687, -0.00102,  0.00034, -0.00318,  0.41067, -0.01219, -0.02918, -0.00675,  0.00333,  0.00549],
       [ 0.00214, -0.03718,  0.0038 , -0.0056 ,  0.00409,  0.00431, -0.01219,  0.02907,  0.00724,  0.01565,  0.00402, -0.01535],
       [-0.02579,  0.0048 , -0.03383,  0.04462, -0.00816, -0.04919, -0.02918,  0.00724,  0.0954 ,  0.01035, -0.00389, -0.01238],
       [ 0.0556 , -0.00226,  0.00542,  0.09931, -0.00392,  0.00616, -0.00675,  0.01565,  0.01035, -0.14815, -0.00947, -0.02193],
       [ 0.00242, -0.07637, -0.125  ,  0.00372, -0.75168,  0.54882,  0.00333,  0.00402, -0.00389, -0.00947,  0.82402, -0.41993],
       [ 0.01717, -0.0006 , -0.10542, -0.00074,  0.43588, -0.04582,  0.00549, -0.01535, -0.01238, -0.02193, -0.41993,  0.16362]])

上述矩阵是对称矩阵了,我们可以看看它是否与 Gaussian 的结果吻合:

[14]:
np.allclose(scf_hess.de.swapaxes(1, 2).reshape(dhess, dhess), fchk.hessian())
[14]:
False

看似是不吻合的。但如果我们稍稍放低一些判断标准,将绝对值误差 atol 容忍到 \(10^{-6}\),或相对值误差容忍到 \(10^{-4}\),就能认为 PySCF 的计算结果与 Gaussian 接近了。

[15]:
np.allclose(scf_hess.de.swapaxes(1, 2).reshape(dhess, dhess), fchk.hessian(), atol=1e-6, rtol=1e-4)
[15]:
True

我们以后一般也沿用上述的评判标准,判断两矩阵或张量是否相等。

5.1.3. pyxdh 计算 RHF Hessian 矩阵

pyxdh 也提供 RHF 的 Hessian 计算。我们要首先给出其梯度辅助类 GradSCF 的实例 grdh

[16]:
config = {"scf_eng": scf_eng}
grdh = GradSCF(config)
grdh.E_1
[16]:
array([[-0.06727,  0.06951,  0.0961 ],
       [ 0.01291,  0.14195, -0.11756],
       [ 0.03423,  0.01409,  0.03949],
       [ 0.02013, -0.22555, -0.01803]])

随后通过上述的实例 grdh 构建 Hessian 辅助类 HessSCF 的实例 hessh

[17]:
config = {"deriv_A": grdh, "deriv_B": grdh}
hessh = HessSCF(config)
hessh.E_2
[17]:
array([[ 0.36765, -0.01096, -0.02986, -0.02036,  0.0064 ,  0.03848, -0.4029 ,  0.00214, -0.02579,  0.0556 ,  0.00242,  0.01717],
       [-0.01096,  0.02901,  0.11159,  0.0047 ,  0.08453, -0.11579,  0.00851, -0.03718,  0.0048 , -0.00226, -0.07637, -0.0006 ],
       [-0.02986,  0.11159,  0.47024, -0.00243,  0.00961, -0.33099,  0.02687,  0.0038 , -0.03383,  0.00542, -0.125  , -0.10542],
       [-0.02036,  0.0047 , -0.00243, -0.07793, -0.00283, -0.04145, -0.00102, -0.0056 ,  0.04462,  0.09931,  0.00372, -0.00074],
       [ 0.0064 ,  0.08453,  0.00961, -0.00283,  0.66306, -0.43734,  0.00034,  0.00409, -0.00816, -0.00392, -0.75168,  0.43588],
       [ 0.03848, -0.11579, -0.33099, -0.04145, -0.43734,  0.426  , -0.00318,  0.00431, -0.04919,  0.00616,  0.54882, -0.04582],
       [-0.4029 ,  0.00851,  0.02687, -0.00102,  0.00034, -0.00318,  0.41067, -0.01219, -0.02918, -0.00675,  0.00333,  0.00549],
       [ 0.00214, -0.03718,  0.0038 , -0.0056 ,  0.00409,  0.00431, -0.01219,  0.02907,  0.00724,  0.01565,  0.00402, -0.01535],
       [-0.02579,  0.0048 , -0.03383,  0.04462, -0.00816, -0.04919, -0.02918,  0.00724,  0.0954 ,  0.01035, -0.00389, -0.01238],
       [ 0.0556 , -0.00226,  0.00542,  0.09931, -0.00392,  0.00616, -0.00675,  0.01565,  0.01035, -0.14815, -0.00947, -0.02193],
       [ 0.00242, -0.07637, -0.125  ,  0.00372, -0.75168,  0.54882,  0.00333,  0.00402, -0.00389, -0.00947,  0.82402, -0.41993],
       [ 0.01717, -0.0006 , -0.10542, -0.00074,  0.43588, -0.04582,  0.00549, -0.01535, -0.01238, -0.02193, -0.41993,  0.16362]])

我们可以验证上述 Hessian 矩阵是否与 Gaussian 相等:

[18]:
np.allclose(hessh.E_2, fchk.hessian(), atol=1e-6, rtol=1e-4)
[18]:
True

5.2. 数值导数求取 Hessian

5.2.1. Hessian 矩阵中单个值的计算

事实上,Hessian 就是能量值的二阶导数构成的矩阵。我们拿第 1 个氧原子的 \(z\) 轴分量、与第 1 个氢原子的 \(x\) 轴分量的 Hessian 矩阵值来举例:

[19]:
fchk.hessian()[2, 6]
[19]:
0.0268686422

之所以索引是 \((2, 6)\),是因为第一个氧原子占用索引 0, 1, 2,其 \(z\) 轴分量则是索引 2;而第 1 个氢原子占用索引 6, 7, 8,其 \(x\) 轴分量则是索引 6。

我们指出,Hessian 矩阵具有对称性,即 \(E^{A_t B_s} = E^{B_s A_t}\),或者我们也能发现下述矩阵值与上面的值是一样的:

[20]:
fchk.hessian()[6, 2]
[20]:
0.0268686422

我们之前已经会三点差分的一阶导数了,事实上求取二阶导数也是相同的。我们首先定义三点差分计算中需要使用到的 \(x - h\) 的点与 \(x + h\) 的点 (分子) mol_m1, mol_p1。这里的 \(x\) 相当于分子的原始坐标,\(h\) 相当于第 1 个氢原子 \(x\) 分量求导所用的偏移量。这里采用的偏移量 (逼近参数) 是 \(10^{-4}\),单位 Bohr。

[21]:
def gen_H2O2(coord):
    """
    Generate H2O2 molecule (with basis 6-31G)
    """
    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()
    mol.set_geom_(coord * lib.param.BOHR)
    return mol.build()
[22]:
coord_orig = mol.atom_coords()
coord_m1 = coord_orig.copy()
coord_m1[2, 0] -= 1e-4
coord_p1 = coord_orig.copy()
coord_p1[2, 0] += 1e-4
[23]:
mol_m1 = gen_H2O2(coord_m1)
mol_p1 = gen_H2O2(coord_p1)

随后,我们可以对上述用于三点差分的分子计算其分子力 \(E^{A_t}\)

[24]:
grad_m1 = scf.RHF(mol_m1).run().nuc_grad_method().run().de
grad_p1 = scf.RHF(mol_p1).run().nuc_grad_method().run().de

我们对上述分子力的第 1 个氧原子 (索引 0) \(z\) 坐标分量 (索引 2) 的值作三点差分导数计算 (相当于 \(E^{A_t B_s} = \frac{\partial E_{A_t}}{\partial B_s}\)):

[25]:
(grad_p1[0, 2] - grad_m1[0, 2]) / (2e-4)
[25]:
0.02686858719513907

我们就会发现,上述的值与 Hessian 矩阵中对应的值是相等的:

[26]:
fchk.hessian()[2, 6]
[26]:
0.0268686422

事实上,我们也可以对分子力的所有值作三点差分:

[27]:
(grad_p1 - grad_m1).flatten() / (2e-4)
[27]:
array([-0.4029 ,  0.00851,  0.02687, -0.00102,  0.00034, -0.00318,  0.41067, -0.01218, -0.02918, -0.00675,  0.00333,  0.00549])

这其实与 Hessian 关于第 1 个氢原子的 \(x\) 轴导数部分完全一致:

[28]:
fchk.hessian()[:, 6]
[28]:
array([-0.4029 ,  0.00851,  0.02687, -0.00102,  0.00034, -0.00318,  0.41067, -0.01219, -0.02918, -0.00675,  0.00333,  0.00549])

我们对其中一个坐标分量作数值导数,就可以得到 Hessian 矩阵的一行。很容易想到,如果我们对所有分子坐标分量作导数,那么完整的 Hessian 矩阵就能获得了。至此,我们就描述好了数值导数计算 Hessian 的原理。

5.2.2. pyxdh 数值梯度助手计算 Hessian

我们以前介绍过,使用 pyxdh 的 NucCoordDerivGenerator 类进行能量的核坐标导数 \(\frac{\partial E}{\partial A_t}\);事实上,这个类原则上可以帮助实现任意维度张量的导数,譬如我们现在需要计算 \(\frac{\partial E^{A_t}}{\partial B_s}\)

相对于之前的文档,这里在生成 NucCoordDerivGenerator 实例时,lambda 函数输入仍然是分子实例,但将 lambda 的输出更改为 pyxdh.grad.RHF 类型作为计算实例。相应的, NumericDiff 实例的 lambda 函数也要更改成输入 pyxdh.grad.RHF 类型,输出分子的梯度矢量。

[29]:
generator = NucCoordDerivGenerator(mol, lambda mol_: scf.RHF(mol_).run().nuc_grad_method().run())
diff = NumericDiff(generator, lambda mf: mf.de.flatten())

最后,我们求取梯度,就得到 Hessian 矩阵:

[30]:
diff.derivative
[30]:
array([[ 0.36765, -0.01096, -0.02986, -0.02036,  0.0064 ,  0.03848, -0.4029 ,  0.00214, -0.02579,  0.0556 ,  0.00242,  0.01717],
       [-0.01096,  0.02901,  0.11159,  0.0047 ,  0.08453, -0.11579,  0.00851, -0.03718,  0.0048 , -0.00226, -0.07637, -0.0006 ],
       [-0.02985,  0.11158,  0.47024, -0.00244,  0.00962, -0.33099,  0.02687,  0.0038 , -0.03383,  0.00542, -0.125  , -0.10542],
       [-0.02036,  0.0047 , -0.00243, -0.07793, -0.00283, -0.04145, -0.00102, -0.0056 ,  0.04462,  0.09931,  0.00372, -0.00074],
       [ 0.0064 ,  0.08453,  0.00961, -0.00282,  0.66305, -0.43734,  0.00034,  0.00409, -0.00816, -0.00391, -0.75168,  0.43588],
       [ 0.03848, -0.11578, -0.33099, -0.04145, -0.43735,  0.42601, -0.00318,  0.00431, -0.04919,  0.00616,  0.54882, -0.04582],
       [-0.4029 ,  0.00851,  0.02687, -0.00102,  0.00034, -0.00318,  0.41067, -0.01219, -0.02918, -0.00675,  0.00333,  0.00549],
       [ 0.00214, -0.03718,  0.0038 , -0.0056 ,  0.00409,  0.00431, -0.01218,  0.02907,  0.00724,  0.01565,  0.00402, -0.01535],
       [-0.02579,  0.0048 , -0.03383,  0.04462, -0.00816, -0.04919, -0.02918,  0.00724,  0.0954 ,  0.01035, -0.00389, -0.01238],
       [ 0.0556 , -0.00226,  0.00542,  0.09931, -0.00392,  0.00616, -0.00675,  0.01565,  0.01035, -0.14815, -0.00947, -0.02193],
       [ 0.00242, -0.07637, -0.125  ,  0.00372, -0.75168,  0.54882,  0.00333,  0.00402, -0.00389, -0.00948,  0.82402, -0.41993],
       [ 0.01717, -0.0006 , -0.10542, -0.00073,  0.43588, -0.04582,  0.00549, -0.01535, -0.01238, -0.02193, -0.41993,  0.16362]])

上面的数值梯度与 Gaussian 结果出入稍大 (这可能与收敛判标有关);我们再降低一些判定条件,将绝对值条件降为 \(10^{-5}\),则可以判定数值梯度与 Gaussian 梯度接近等同:

[31]:
np.allclose(diff.derivative, fchk.hessian(), atol=1e-5, rtol=1e-4)
[31]:
True

5.2.3. 通过数值导数计算分子振动频率

我们之前提及,Hessian 的计算的一个很重要的意义是计算分子的振动频率。在这里,我们就不详细讨论如何计算频率。

我们引入一个自编的 Python 脚本 freqanal.py,该脚本能帮助我们进行频率计算:

[32]:
from freqanal import FreqAnal

其输入参量是原子质量列表 (单位 AMU)、分子坐标 (单位 Bohr)、以及 Hessian (单位 a.u.)。为了与 Gaussian 的分子频率作核对,我们需要额外定义较为精确的原子质量。

[33]:
mol_weight = np.array([15.99491, 15.99491, 1.00783, 1.00783])
freqanal = FreqAnal(mol_weight=mol_weight, mol_coord=mol.atom_coords(), hessian=diff.derivative)
freqanal.freq
[33]:
array([-1580.60283, -1218.3735 ,  1370.61697,  1647.5344 ,  3389.85998,  5347.88146])

我们可以将其与 Gaussian 输出的频率值作核对:

[34]:
with open("assets/H2O2-freq.out", "r") as f:
    for idx, line in enumerate(f.readlines()):
        if "Frequencies" in line:
            print("line {:4d}:".format(idx + 1), line[:-1])
line  574:  Frequencies -- -1580.6089             -1218.3809              1370.6206
line  588:  Frequencies --  1647.5426              3389.8666              5347.9015