6. RHF 极化率

上一节，我们已经讨论过了 RHF 的核坐标二阶导数 (即 Hessian) 的计算。这一节，我们简单描述电场强度的二阶导数，即极化率。我们不打算对该物理量作其物理意义的描述。

事实上，Hessian 在数学中就相当于对向量作二阶导数；但在计算化学中，它通常特指对原子核坐标的二阶导数，因此我们不称极化率矩阵为 Hessian 矩阵。

from pyscf import gto, scf, lib, prop
import numpy as np
from pyxdh.Utilities import FormchkInterface, DipoleDerivGenerator, NumericDiff
from pyxdh.Utilities.test_molecules import Mol_H2O2
from pyxdh.DerivOnce import DipoleSCF
from pyxdh.DerivTwice import PolarSCF
import warnings

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

6.1. 量化软件计算 RHF 极化率

6.1.1. Gaussian 极化率

在上一篇文档，我们已经给出从 Gaussian 频率分析的程序输出得到极化率的过程。尽管极化率严格来说不是频率分析的结果，但一般来说，从频率分析过程中的矩阵导出极化率是不太消耗计算量的；因此，Gaussian 会在频率分析过程中给出极化率。

fchk = FormchkInterface("assets/H2O2-freq.fch")
fchk.polarizability()

array([[ 6.58142, -0.0841 , -1.45378],
       [-0.0841 ,  4.26836,  0.39969],
       [-1.45378,  0.39969, 17.89033]])

可以见到，极化率是与分子大小无关的 \((3, 3)\) 维度的矩阵。极化率一般写为 \(\alpha_{ts}\)，其中 \(\alpha\) 是极化率量的表示符号，下标 \(t, s\) 都表示三维电子坐标的分量。

6.1.2. PySCF 极化率

PySCF 可以通过其 prop 库进行计算；它需要一个 scf.RHF 的计算实例 scf_eng，用于生成计算极化率的计算实例 scf_polar：

molh = Mol_H2O2()
mol = molh.mol
scf_eng = molh.hf_eng.run()
scf_polar = prop.polarizability.rhf.Polarizability(scf_eng).run()

极化率的结果可以用 polarizability 成员函数得到。

scf_polar.polarizability()

array([[ 6.58142, -0.0841 , -1.45378],
       [-0.0841 ,  4.26835,  0.39969],
       [-1.45378,  0.39969, 17.89033]])

可见 Gaussian 与 PySCF 的结果非常相近：

np.allclose(scf_polar.polarizability(), fchk.polarizability())

True

6.1.3. pyxdh 极化率

pyxdh 也提供了计算 RHF 的函数。其调用方式比较类似于上一篇文档提到的 HessSCF 的使用方法。首先，我们需要先定义 DipoleSCF 的实例 diph 计算偶极矩：

diph = DipoleSCF({"scf_eng": scf_eng})
diph.E_1

array([ 0.88992,  0.66299, -0.29469])

随后，我们将 diph 代入到 PolarSCF 的实例化过程中，得到极化率实例 polh：

polh = PolarSCF({"deriv_A": diph, "deriv_B": diph})
- polh.E_2

array([[ 6.58142, -0.0841 , -1.45379],
       [-0.0841 ,  4.26835,  0.39969],
       [-1.45379,  0.39969, 17.89032]])

需要注意，上述的 E_2 property 调用后还需要乘以 -1，才能得到极化率结果。我们拿该结果与 Gaussian 核对：

np.allclose(- polh.E_2, fchk.polarizability())

True

6.2. 数值导数得到极化率

6.2.1. 三点差分得到极化率

现在我们假设需要得到极化率的第 1 行 \(\alpha_{ty}\)；那么，我们对 Hamiltonian Core 作如下的变化：

\[\hat h (F) = \hat t + \hat v_\mathrm{nuc} + F y\]

其中，\(F\) 为微扰的外场强度。这个微扰外场就相当于三点差分过程中 \(x - h\) 与 \(x + h\) 的逼近参数 \(h\)。上式的 \(y\) 表示的是我们外加电场的分量取向；由于我们是求其中一个取向为 \(y\) 方向的极化率，那么我们对解析偶极矩的 \(y\) 方向作数值求导即可。作为逼近参数的微扰外加电场大小是 \(10^{-4}\)，单位为原子单位。

def mf_func(t, f):
    mf = scf.RHF(mol)
    mf.conv_tol = 1e-10
    mf.get_hcore = lambda mol_: scf.rhf.get_hcore(mol_) - f * mol_.intor("int1e_r")[t]
    return mf.run()

scf_eng_p1 = mf_func(1,  1e-4)
scf_eng_m1 = mf_func(1, -1e-4)
- (scf_eng_p1.dip_moment(unit="A.U.") - scf_eng_m1.dip_moment(unit="A.U.")) / 2e-4

Dipole moment(X, Y, Z, A.U.):  0.88992,  0.66256, -0.29473
Dipole moment(X, Y, Z, A.U.):  0.88991,  0.66342, -0.29465

array([-0.0841 ,  4.26836,  0.39967])

上面的程序使用了一次负号；这个负号是因为两次负电荷的的负号累加导致，因此需要消去一个负号。

我们已经通过三点差分得到了极化率中关于 \(y\) 取向的一行了；那么剩下两个取向也会是非常容易获得的了。

6.2.2. pyxdh 数值求导

pyxdh 的数值求导机制就是上述过程，即需要一个生成更变了 Hamiltonian Core 的计算实例 mf_func 以实例化 DipoleDerivGenerator，随后在 NumDiff 中对偶极矩作三点差分：

generator = DipoleDerivGenerator(mf_func, interval=1e-6)
diff = NumericDiff(generator, lambda mf: mf.dip_moment(unit="A.U.", verbose=0))

那么，极化率的值就可以导出如下：

- diff.derivative

array([[ 6.58142, -0.0841 , -1.45377],
       [-0.0841 ,  4.26836,  0.39967],
       [-1.45379,  0.39968, 17.89025]])

我们可以验证上述极化率是否与 Gaussian 相等：

np.allclose(- diff.derivative, fchk.polarizability(), atol=1e-6, rtol=1e-4)

True