7. 磁性质数值导数 (3)：RMP2 的 GIAO 核磁 (屏蔽) 共振常数 (NMR)#

创建时间：2020-08-31；修订前的文档对原子轨道积分的磁化率梯度叙述上有错误。

修订时间：2022-02-09

在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算 GIAO 的 RHF 数值核磁 (屏蔽) 共振常数 (Nuclear Magnetic Resonance constant, NMR) 的程序。该文档大量参考 PySCF 的代码 nmr/rhf.py。

该文档经过修订。之所以发现先前文档的错误，是因为读到最近的一篇数值 NMR 工作[1]。目前版本的文档应当能处理类似于 MP2、dRPA@HF、CCSD 等 Post-HF 方法。DFT 方法的 NMR 暂不在本文档的讨论范畴。

警告

该文档已经修订完毕，但使用到了先前的一份文档对 GIAO 的公式、记号与讨论。先前的文档还未修订完毕。

之前犯的主要问题是，在解析导数求解时，将 Fock 贡献归入 Core Hamiltonian 是程序实现上非常方便的；但原理上，Fock 贡献要拆分为单电子与双电子部分。因此，在求数值导数时，自洽场部分的双电子积分 (ERIs) 也应是受了磁场影响的。

与之前的文档一样，我们的讨论中所使用到的分子体系 mol 会是非对称的氨分子，并且取用最小基组。其 RHF 计算放在实例 mf，而 NMR 计算实例会放在 mf_nmr。

7.1. 准备工作#

import warnings
warnings.filterwarnings("ignore")

from pyscf import gto, scf, mp
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.253540514686556

RHF 的核磁屏蔽张量 (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 为单位的表达。不过我们在这篇文档打算考虑 RMP2 的核磁计算；这里的 RHF NMR 张量只是用于演示而已。

7.2. 微扰原子矩阵的表达#

7.2.1. 外磁场微扰的一阶 Core Hamiltonian#

关于该微扰矩阵，即是统合磁效应在 Core Hamiltonian 的作用、以及 GIAO 轨道对动能与核静电势能的一阶贡献考虑进来：

\[ \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 \end{align*} \]

hcore_1_B = - 1j * (
    + 0.5 * mol.intor('int1e_giao_irjxp', 3)
    + mol.intor('int1e_ignuc', 3)
    + mol.intor('int1e_igkin', 3))

7.2.2. 外磁场微扰的一阶重叠矩阵#

\[ S_{\mu \nu}^{\mathscr{B}_t} = \langle U_\mathrm{g}^t \mu | \nu \rangle \]

我们用下述代码生成 ovlp_1_B \(S_{\mu \nu}^{\mathscr{B}_t}\)：

ovlp_1_B = - 1j * mol.intor("int1e_igovlp")

7.2.3. 外磁场微扰的一阶 ERI 张量#

\[ (\mu \nu | \kappa \lambda)^{\mathscr{B}_t} = (U_\mathrm{g}^t \mu \nu | \kappa \lambda) + (\mu \nu | U_\mathrm{g}^t \kappa \lambda) \]

我们用下述代码生成 eri_1_B \((\mu \nu | \kappa \lambda)^{\mathscr{B}_t}\)：

eri_1_B = -1j * (
    + np.einsum("tuvkl -> tuvkl", mol.intor('int2e_ig1'))
    + np.einsum("tkluv -> tuvkl", mol.intor('int2e_ig1')))

7.2.4. 核磁偶极的一阶 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")

7.2.5. 磁场与核磁偶极的二阶 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))

7.3. 数值导数求 RMP2 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 重叠矩阵、_eri 双电子积分的 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}} + o(\mathscr{B}_t \mu_{A_s}) \\ S_{\mu \nu} &= S_{\mu \nu}^{(0)} + \mathscr{B}_t S_{\mu \nu}^{\mathscr{B}_t} + o(\mathscr{B}_t) \\ (\mu \nu | \kappa \lambda) &= (\mu \nu | \kappa \lambda)^{(0)} + (\mu \nu | \kappa \lambda)^{\mathscr{B}_t} + o(\mathscr{B}_t) \end{align} \end{split}\]

在跑完自洽场后，立即简单地执行 MP2，就可以得到受外磁场与核磁偶极扰动的 MP2 能量了。

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
    mf._eri = mol.intor("int2e") + np.einsum("tuvkl, t -> uvkl", eri_1_B, dev_xyz_B)
    mf.run()
    mf_mp = mp.MP2(mf).run()
    return mf_mp.e_tot

随后使用与前文类似的数值梯度方式就能得到核磁屏蔽常数了；所使用的数值差分大小对于外磁场与核磁偶极均为 \(3 \times 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.63279,  0.41743,  0.26518],
        [-0.11122,  3.81459, -0.13199],
        [ 0.38918, -0.53118,  7.02831]],

       [[ 0.47397, -0.00892,  0.00238],
        [-0.02979,  0.68584, -0.06142],
        [-0.0151 ,  0.0332 ,  0.52897]],

       [[ 0.44192, -0.0033 ,  0.01888],
        [ 0.0113 ,  0.45509, -0.01335],
        [-0.01723,  0.00653,  0.59529]],

       [[ 0.73128,  0.08291, -0.13904],
        [ 0.03982,  0.50304, -0.0239 ],
        [-0.06955, -0.01146,  0.54013]]])

留意到我们之前一直都使用去除结构精细常数的结果，因此我们需要乘以 \(\alpha^2\)。同时由于单位是 ppm，因此最终我们需要乘以 \(10^6 \alpha^2\)：

num_nmr_ppm = num_nmr * constants.alpha**2 * 10**6
num_nmr_ppm

array([[[299.95362,  22.2286 ,  14.12107],
        [ -5.92259, 203.13221,  -7.02871],
        [ 20.72447, -28.28604, 374.26729]],

       [[ 25.2397 ,  -0.47482,   0.12685],
        [ -1.58641,  36.52216,  -3.27057],
        [ -0.80421,   1.76797,  28.16842]],

       [[ 23.53299,  -0.17577,   1.00526],
        [  0.60148,  24.23406,  -0.71087],
        [ -0.91735,   0.3478 ,  31.69988]],

       [[ 38.9417 ,   4.41492,  -7.40411],
        [  2.12046,  26.78774,  -1.27294],
        [ -3.70343,  -0.61013,  28.76271]]])

7.4. RMP2 NMR 核磁屏蔽常数#

最终可以用于汇报的核磁屏蔽常数需要经过 \(3 \times 3\) 矩阵对角化给出。对于其中一个原子而言，其核磁屏蔽张量可以表示为

\[\begin{split} \boldsymbol{\sigma} = \begin{pmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \end{pmatrix} = X \begin{pmatrix} \sigma_{xx}' & 0 & 0 \\ 0 & \sigma_{yy}' & 0 \\ 0 & 0 & \sigma_{zz}' \end{pmatrix} X^\dagger \end{split}\]

其中，\((\sigma_{xx}', \sigma_{yy}', \sigma_{zz}')\) 是屏蔽张量的本征值，不随规范原点变化或坐标旋转而变化；而屏蔽张量 \(\boldsymbol{\sigma}\) 本身是随规范原点或坐标旋转变化而可以改变的。因此，汇报 NMR 数据时，应当要汇报与本征值有关的结果。

一般的 NMR 谱会打出同性核磁屏蔽常数 \(\sigma_\text{iso} = \frac{1}{3} (\sigma_{xx}' + \sigma_{yy}' + \sigma_{zz}')\)。有时异性屏蔽常数 \(\sigma_\text{aniso} = \sigma_{zz}' - \frac{1}{2} (\sigma_{xx}' + \sigma_{yy}')\) 也会使用到[2]。我们可以用下面的程序，对 MP2/STO-3G 的 NMR 屏蔽常数与Gaussian 结果 作对照。

# Gaussian results
! grep "Isotropic" NH3-nmr.out | tail -n 4

N    Isotropic =   292.5252   Anisotropy =   130.4969
H    Isotropic =    29.9780   Anisotropy =    10.0555
H    Isotropic =    26.4878   Anisotropy =     7.8257
H    Isotropic =    31.4998   Anisotropy =    15.9451

def proc_nmr_tensor(tsr):
    tsr = (tsr + tsr.T) / 2
    eigs = np.linalg.eigvalsh(tsr)
    eigs.sort()
    nmr_iso = eigs.sum() / 3
    nmr_ani = eigs[2] - (eigs[0] + eigs[1]) / 2
    return nmr_iso, nmr_ani

# Numerical derivative results
for tsr in num_nmr_ppm:
    print("Isotropic = {:8.4f}  Anisotropy = {:8.4f}".format(*proc_nmr_tensor(tsr)))

Isotropic = 292.4510  Anisotropy = 130.5998
Isotropic =  29.9768  Anisotropy =  10.0491
Isotropic =  26.4890  Anisotropy =   7.8233
Isotropic =  31.4974  Anisotropy =  15.9435