GIAO 的 RHF 核磁 (屏蔽) 共振常数 (NMR) 数值导数#
创建时间：2020-08-31
在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算 GIAO 的 RHF 数值核磁 (屏蔽) 共振常数 (Nuclear Magnetic Resonance constant, NMR) 的程序。该文档大量参考 PySCF 的代码 nmr/rhf.py。
与上一篇文档一样，我们的讨论中所使用到的分子体系 mol 会是非对称的氨分子，并且取用最小基组。其 RHF 计算放在实例 mf，而 NMR 计算实例会放在 mf_nmr。
from pyscf import gto, scf, dft
from pyscf.prop import nmr
from pyscf.data import nist
from scipy import constants
import numpy as np
np.set_printoptions(precision=5, linewidth=150, suppress=True)
mol = gto.Mole()
mol.atom = """
N 0. 0. 0.
H 0. 1. 0.2
H 0.1 0.3 1.5
H 0.9 0.4 -.2
"""
mol.basis = "STO-3G"
mol.verbose = 0
mol.build()
coord_orig = np.zeros(3)
nocc, nao, nmo, natm = mol.nelec[0], mol.nao, mol.nao, mol.natm
其自洽场能量为
mf = scf.RHF(mol).run()
mf.e_tot
-55.25354051468657
核磁屏蔽张量 (Shielding Constant) \(\sigma_{ts}^A\) 可以表示如下 (维度 \((A, t, s)\))：
\[
\sigma_{ts}^A = \frac{\partial^2 E_\mathrm{tot}}{\partial \mathscr{B}_t \partial \mu_{A_s}}
\]
mf_nmr = nmr.RHF(mf)
mf_nmr.kernel()
array([[[305.25641, 22.67382, 12.29714],
[-22.38193, 180.11491, -2.84064],
[ 13.65836, -28.07439, 383.02872]],
[[ 25.68677, -0.27932, -0.31477],
[ -0.84315, 35.70697, -4.26016],
[ -1.06334, 2.07184, 28.17535]],
[[ 24.06192, -0.12066, 0.81476],
[ 0.99821, 24.76595, -0.86457],
[ -0.74954, 0.2574 , 31.65924]],
[[ 39.48533, 4.28363, -7.95068],
[ 0.89957, 26.69818, -0.9659 ],
[ -4.40617, -0.48147, 28.90276]]])
上述结果是以 ppm 为单位的表达。
微扰原子矩阵的表达#
外磁场微扰的一阶 Core Hamiltonian#
关于该微扰矩阵，我们已经在前两篇文档中有比较多的说明了。
\[\begin{split}
\begin{split}\begin{align}
h_{\mu \nu}^{\mathscr{B}_t}
&=
\frac{1}{2} \langle \mu | \hat l_t | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu}
+ \langle U_\mathrm{g}^t \mu | \hat t | \nu \rangle
+ \langle U_\mathrm{g}^t \mu | \hat v_\mathrm{nuc} | \nu \rangle
\\ & \quad
+ \sum_{\kappa \lambda} ( U_\mathrm{g}^t \mu \nu | \kappa \lambda ) D_{\kappa \lambda}^{(0)}
- \frac{1}{2} \sum_{\kappa \lambda} ( U_\mathrm{g}^t \mu \lambda | \kappa \nu ) D_{\kappa \lambda}^{(0)}
- \frac{1}{2} \sum_{\kappa \lambda} ( U_\mathrm{g}^t \kappa \nu | \mu \lambda ) D_{\kappa \lambda}^{(0)}
\end{align}\end{split}
\end{split}\]
由于在 PySCF 中，函数 nmr.rhf.make_h10 就是专门用于生成该矩阵的，因此我们就用下述代码生成 hcore_1_B \(h_{\mu \nu}^{\mathscr{B}_t}\) 表示该矩阵。
dm_guess = mf.make_rdm1()
hcore_1_B = 1j * nmr.rhf.make_h10(mol, dm_guess)
外磁场微扰的一阶重叠矩阵#
该微扰矩阵也已经在前两篇文档中有比较多的说明了。
\[
S_{\mu \nu}^{\mathscr{B}_t} = \langle U_g^t \mu | \nu \rangle
\]
我们用下述代码生成 ovlp_1_B \(S_{\mu \nu}^{\mathscr{B}_t}\)：
ovlp_1_B = - 1j * mol.intor("int1e_igovlp")
核磁偶极的一阶 Core Hamiltonian#
核磁偶极 \(\mu_{A_t}\) 所产生的一阶算符贡献可以表达为
\[
\hat h {}^{(1)} (\boldsymbol{\mu}_A)
= - i \alpha^2 \boldsymbol{\mu}_A \cdot \left( \boldsymbol{\nabla} \frac{1}{\boldsymbol{r - \boldsymbol{R}_A}} \times \boldsymbol{\nabla} \right)
= - i \alpha^2 \boldsymbol{\mu}_A \cdot \left( \frac{\boldsymbol{r} - \boldsymbol{R}_A}{|\boldsymbol{r} - \boldsymbol{R}_A|^3} \times \boldsymbol{\nabla} \right)
\]
其中，\(\boldsymbol{R}_A\) 表示原子 \(A\) 的核坐标，\(\alpha\) 表示精细结构常数，\(1/\alpha \simeq 137\)。该常数可以从 PySCF 中获得，也可以从 SciPy 中获得。注意等式左边的 \(\boldsymbol{\mu_A}\) 看作是外加的核磁偶极大小，而等号右边的 \(\mu\) 表示原子轨道，两者意义不同；等号左边角标 \(\boldsymbol{A}\) 表示原子核坐标向量，而之前两篇文档中的 \(A\) 在很多文章或教材中表示 \(\frac{1}{2} \boldsymbol{B} \times \boldsymbol{r}\)。
1 / constants.alpha
137.0359990836958
但为了方便，在最后核算结果之前，我们会暂且将 \(\alpha\) 当作 1 来处理。
在 PySCF 中，实现上述过程的积分字符是 int1e_ia01p；但其使用需要告知 gto.mole.intor 函数以其 \(1 / \boldsymbol{r}\) 的规范原点位置具体处在哪个原子核中心。
\[
h_{\mu \nu}^{\mu_{A_t}} = - i \langle \mu | \boldsymbol{\nabla} \frac{1}{\boldsymbol{r}} \times \boldsymbol{\nabla} | \nu \rangle_{\text{Gauge of } \boldsymbol{r} \rightarrow \boldsymbol{R}_A}
\]
我们用下述代码生成 hcore_1_m \(h_{\mu \nu}^{\mu_{A_t}}\)：(维度为 \((A, t, \mu, \nu)\))
hcore_1_m = np.zeros((natm, 3, nao, nao), dtype=np.complex128)
for atom_idx in range(natm):
with mol.with_rinv_orig(mol.atom_coord(atom_idx)):
hcore_1_m[atom_idx] = - 1j * mol.intor("int1e_ia01p")
磁场与核磁偶极的二阶 Core Hamiltonian#
磁场与核磁偶极之间的算符乘积会产生二阶算符贡献项：
\[\begin{split}
\begin{align}
\hat h {}^{(2)} (\boldsymbol{\mathscr{B}}, \boldsymbol{\mu}_A) | \nu \rangle
&= \frac{\alpha^2}{2} \boldsymbol{\mathscr{B}}^\mathrm{T} \left( (\boldsymbol{r} - \boldsymbol{R}_\nu) \cdot \boldsymbol{\nabla} \frac{1}{\boldsymbol{r} - \boldsymbol{R}_A} - (\boldsymbol{r} - \boldsymbol{R}_\nu) \boldsymbol{\nabla} \frac{1}{(\boldsymbol{r} - \boldsymbol{R}_A)^\mathrm{T}} \right) \boldsymbol{\mu}_A | \nu \rangle \\
&\quad + \alpha^2 \boldsymbol{\mathscr{B}}^\mathrm{T} \boldsymbol{U}_\mathrm{g} \left( \boldsymbol{\nabla} \frac{1}{\boldsymbol{r} - \boldsymbol{R}_A} \times \boldsymbol{\hat p} \right)^\mathrm{T} \boldsymbol{\mu}_A | \nu \rangle
\end{align}
\end{split}\]
去除精细结构常数 \(\alpha\) 的贡献后，其矩阵的表达形式则是
\[\begin{split}
\begin{align}
h_{\mu \nu}^{A_s t} \mathscr{B}_t \mu_{A_s}
&= \langle \mu | - \frac{1}{2} \frac{(t - t_\nu) (s - s_A)}{|\boldsymbol{r} - \boldsymbol{R}_A|^3} | \nu \rangle
- \delta_{ts} \langle \mu | - \frac{1}{2} \sum_{w} \frac{(w - w_\nu) (w - w_A)}{|\boldsymbol{r} - \boldsymbol{R}_A|^3} | \nu \rangle \\
&\quad + \langle U_\mathrm{g}^t \mu | \left( \boldsymbol{\nabla} \frac{1}{\boldsymbol{r} - \boldsymbol{R}_A} \times \boldsymbol{\hat p} \right)_s | \nu \rangle
\end{align}
\end{split}\]
我们用下述代码生成 hcore_2 \(h_{\mu \nu}^{A_s t}\)：(维度为 \((A, t, s, \mu, \nu)\)，注意维度 \(t\) 对应外磁场，而 \(A, s\) 对应核磁偶极)
hcore_2 = np.zeros((natm, 3, 3, nao, nao))
for atom_idx in range(natm):
with mol.with_rinv_origin(mol.atom_coord(atom_idx)):
hcore_2[atom_idx] += mol.intor("int1e_giao_a11part").reshape((3, 3, nao, nao))
hcore_2[atom_idx] -= np.einsum("ts, uv -> tsuv", np.eye(3), mol.intor("int1e_giao_a11part").reshape((3, 3, nao, nao)).trace(axis1=0, axis2=1))
hcore_2[atom_idx] += mol.intor("int1e_a01gp").reshape((3, 3, nao, nao))
数值导数求 NMR 核磁屏蔽张量#
随后我们就可以通过数值梯度求核磁屏蔽张量 \(\sigma_{ts}^A\) (维度 \((A, t, s)\))：
\[
\sigma_{ts}^A = \frac{\partial^2 E_\mathrm{tot}}{\partial \mathscr{B}_t \partial \mu_{A_s}}
\]
在此之前，我们仍然需要构造一个通过更改 get_hcore Core Hamiltonian 与 get_ovlp 重叠矩阵的 PySCF 自洽场实例，以施加外场获得能量的函数 eng_nmr_field。其输入的参数 dev_xyz_B 是三维外加磁场大小 (对应维度 \(t\))，dev_xyz_m 是三维外加核磁偶极大小 (对应维度 \(s\))，atom_idx 是原子序号 (对应维度 \(A\))。
\[\begin{split}
\begin{align}
h_{\mu \nu} &= h_{\mu \nu}^{(0)} + \mathscr{B}_t h_{\mu \nu}^{\mathscr{B}_t} + \mu_{A_s} h_{\mu \nu}^{\mu_{A_s}} + \mathscr{B}_t \mu_{A_s} h_{\mu \nu}^{\mathscr{B}_t \mu_{A_s}} \\
S_{\mu \nu} &= S_{\mu \nu}^{(0)} + \mathscr{B}_t S_{\mu \nu}^{\mathscr{B}_t}
\end{align}
\end{split}\]
dm_guess = mf.make_rdm1()
def eng_nmr_field(dev_xyz_B, dev_xyz_m, atom_idx):
mf = scf.RHF(mol)
def get_hcore(mol_=mol):
hcore_total = np.asarray(scf.rhf.get_hcore(mol_), dtype=np.complex128)
hcore_total += np.einsum("tuv, t -> uv", hcore_1_B, dev_xyz_B)
hcore_total += np.einsum("tuv, t -> uv", hcore_1_m[atom_idx], dev_xyz_m)
hcore_total += np.einsum("tsuv, t, s -> uv", hcore_2[atom_idx], dev_xyz_B, dev_xyz_m)
return hcore_total
def get_ovlp(mol_):
ovlp_total = np.asarray(scf.rhf.get_ovlp(mol_), dtype=np.complex128)
ovlp_total += np.einsum("tuv, t -> uv", ovlp_1_B, dev_xyz_B)
return ovlp_total
mf.get_hcore = get_hcore
mf.get_ovlp = get_ovlp
return mf.kernel(dm=dm_guess)
随后使用与前文类似的数值梯度方式就能得到核磁屏蔽常数了；所使用的数值差分大小对于外磁场与核磁偶极均为 \(10^{-4}\) a.u.。
interval = 1e-4
num_nmr = np.zeros((natm, 3, 3))
for atom_idx in range(natm):
for t in range(3):
for s in range(3):
dev_xyzs_B, dev_xyzs_m = np.zeros((2, 3)), np.zeros((2, 3))
dev_xyzs_B[0, t] = dev_xyzs_m[0, s] = -interval
dev_xyzs_B[1, t] = dev_xyzs_m[1, s] = interval
num_nmr[atom_idx, t, s] = (
+ eng_nmr_field(dev_xyzs_B[0], dev_xyzs_m[0], atom_idx)
- eng_nmr_field(dev_xyzs_B[1], dev_xyzs_m[0], atom_idx)
- eng_nmr_field(dev_xyzs_B[0], dev_xyzs_m[1], atom_idx)
+ eng_nmr_field(dev_xyzs_B[1], dev_xyzs_m[1], atom_idx)
) / (4 * interval**2)
num_nmr
array([[[ 5.73238, 0.42579, 0.23092],
[-0.4203 , 3.38236, -0.05334],
[ 0.25651, -0.52722, 7.19284]],
[[ 0.48237, -0.00525, -0.00591],
[-0.01584, 0.67054, -0.08 ],
[-0.01997, 0.03891, 0.5291 ]],
[[ 0.45185, -0.00226, 0.01529],
[ 0.01875, 0.46508, -0.01624],
[-0.01408, 0.00484, 0.59453]],
[[ 0.74149, 0.08044, -0.14931],
[ 0.01689, 0.50136, -0.01814],
[-0.08274, -0.00904, 0.54276]]])
留意到我们之前一直都使用去除结构精细常数的结果，因此我们需要乘以 \(\alpha^2\)。同时由于单位是 ppm，因此最终我们需要乘以 \(10^6 \alpha^2\)：
num_nmr * constants.alpha**2 * 10**6
array([[[305.25688, 22.67379, 12.2969 ],
[-22.38171, 180.11542, -2.84064],
[ 13.65948, -28.0751 , 383.02872]],
[[ 25.68659, -0.27935, -0.3146 ],
[ -0.84328, 35.70701, -4.26026],
[ -1.06338, 2.07184, 28.17533]],
[[ 24.06171, -0.12057, 0.81443],
[ 0.99832, 24.76594, -0.86458],
[ -0.74969, 0.25765, 31.65949]],
[[ 39.48546, 4.28353, -7.95088],
[ 0.89964, 26.69806, -0.9661 ],
[ -4.40616, -0.48158, 28.9029 ]]])