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