6. 磁性质数值导数 (2)：RHF 的 GIAO 磁化率#

创建时间：2020-08-30

危险

该文档的数值导数策略有误。数值导数应当要对双电子排斥积分 (ERIs) 作改动，而非将 Fock 矩阵的贡献纳入 Core Hamiltonian 中。

近日该文档将会作修订。

在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算 GIAO 的 RHF 数值磁化率的程序。该文档大量参考 PySCF 的代码 magnetizability/rhf.py 与 nmr/rhf.py。一篇公式记号比较清晰的文章是 Laasner, Blum, et al. [1]。

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

在这篇文档中，我们仍然需要保留规范原点 (Gauge Origin) 的概念。规范原点会取在坐标原点上 coord_orig。

from pyscf import gto, scf, dft
from pyscf.prop import nmr, magnetizability
import numpy as np

np.set_printoptions(precision=5, linewidth=150, suppress=True)
np.random.seed(0)

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 = mol.nelec[0], mol.nao, mol.nao

其自洽场能量为

mf = scf.RHF(mol).run()
mf.e_tot

-55.253540514686556

我们定义磁化率计算类为 mf_mag。其磁化率张量 \(\xi_{ts}\) 为：

mf_mag = magnetizability.RHF(mf)
mf_mag.kernel()

array([[-3.72497, -0.02158, -0.07514],
       [-0.02158, -2.88877,  0.20345],
       [-0.07514,  0.20345, -3.57485]])

但留意到，与上一份文档不同地，如果我们对分子作平移操作 (譬如下述平移操作是将原点移动 \((x, y, z) = (10, -10, 5) \, \mathsf{a.u.}\))，其磁化率仍然保持不变：

mol_trans = mol.copy()
mol_trans.set_geom_(mol.atom_coords() + np.array([[10, -10, 5]]), unit="AU")
mol_trans.build()

<pyscf.gto.mole.Mole at 0x7fa032502d00>

magnetizability.RHF(scf.RHF(mol_trans).run()).kernel()

array([[-3.72498, -0.0216 , -0.07515],
       [-0.0216 , -2.8888 ,  0.20345],
       [-0.07515,  0.20345, -3.57482]])

6.1. GIAO 规范不变原子轨道#

6.1.1. 基本概念#

我们从上一篇文档中，得知对分子作平移操作后，磁化率可能会改变。

GIAO 规范不变原子轨道 (Gauge Invariant Atomic Orbital) 是一种解决方案，它可以让分子在平移操作后，保证磁化率不变。事实上，分子旋转操作下磁化率也不变，但这会额外引入三维旋转矩阵，为了方便我们就不讨论分子旋转的情形。

GIAO 的基本思路是对于外磁场引入的微扰，其波函数更换为

\[ \phi^\mathrm{GIAO}_\mu (\boldsymbol{r}) = e^{- \frac{i}{2} (\boldsymbol{R}_\mu \times \boldsymbol{r}) \cdot \boldsymbol{\mathscr{B}}} \phi_\mu (\boldsymbol{r}) \]

其中，\(\phi_\mu (\boldsymbol{r})\) 是普通的原子轨道基组，\(\boldsymbol{R}_\mu\) 是原子轨道 \(\mu\) 作为 Gaussian 基组的中心坐标相对于规范原点 (Gauge Origin) 向量，\(\boldsymbol{r}\) 是电子坐标，\(\boldsymbol{\mathscr{B}}\) 是外加微扰磁场。根据其定义，我们知道，GIAO 方法必须要使用原子轨道基组下使用。需要注意，这是一个复函数轨道。

使用该转换后的原子轨道就可以达到平移操作后磁化率不变的结果。关于这一点的证明，参考 Pople [2] eq 2.5 附近的讨论。

事实上，规范 (原点) 不变原子轨道的称呼可能是不恰当的。Pople [3] 曾在文章里使用了规范 (原点) 依赖原子轨道 (Gauge Dependent Atomic Orbital)；这是因为 \(\phi^\mathrm{GIAO}_\mu (\boldsymbol{r})\) 实际上是依规范原点位置不同而不同的 (注意 \(\boldsymbol{R}_\mu\) 的定义)。当然，我们仍然沿用 GIAO 的约定俗成。

6.1.2. 实例分析：GIAO 重叠矩阵在外磁场下的微扰#

概念与程序准备

我们拿一个具体的例子来说明 GIAO 的计算；最简单的例子是重叠矩阵：

\[ S_{\mu \nu}^\mathrm{GIAO} = \langle \mu | e^{- \frac{i}{2} (\boldsymbol{R}_{\mu \nu} \times \boldsymbol{r}) \cdot \boldsymbol{\mathscr{B}}} | \nu \rangle \]

其中，\(\boldsymbol{R}_{\mu \nu} = \boldsymbol{R}_\nu - \boldsymbol{R}_\mu\)，它表示原子轨道 \(\mu\) 与 \(\nu\) 作为 Gaussian 基组的中心坐标的向量之差。之所以会出现 \(- \boldsymbol{R}_\mu\) 的负号，是因为 \(\phi^\mathrm{GIAO}_\mu (\boldsymbol{r})\) 出现在左矢时，其虚数使得在作复共轭时，指数项应当乘以负号。

我们在实际计算中，只关心其一阶算符的矩阵形式：

\[ \nabla_{\boldsymbol{\mathscr{B}}} S_{\mu \nu}^\mathrm{GIAO} = \langle \mu | - \frac{i}{2} \boldsymbol{R}_{\mu \nu} \times \boldsymbol{r} | \nu \rangle \]

其分量形式不太容易写出，但我们下面会用程序来具体地求取其分量。

我们先使用格点积分来计算矩阵。借助 PySCF 的 DFT 格点，我们写出下述代码：

ni 格点积分引擎
grids (50, 194) 大小的格点
ao \(\phi_\mu (\boldsymbol{r}_g)\) 即处于格点 \(\boldsymbol{r}_g\) 上 \(\phi_\mu\) 基轨道的函数值，维度表示为 \((g, \mu)\)

ni = dft.numint.NumInt()

grids = dft.Grids(mol)
grids.atom_grid = (50, 194)
grids.build()

ao = ni.eval_ao(mol, grids.coords)
ao.shape

(26896, 8)

数值格点积分

我们预先已经知道，第 \(\mu = 3\) 根基轨道是第 0 个原子 (N 原子)，第 \(\nu = 6\) 跟轨道是第 2 个原子 (H 原子) (按照 0-index 计数)，那么向量 \(\boldsymbol{R}_{\mu \nu} = \boldsymbol{R}_{36} = \boldsymbol{R}_{6} - \boldsymbol{R}_{3}\)：

R_36 = (mol.atom_coord(2) - coord_orig) - (mol.atom_coord(0) - coord_orig)
R_36

array([0.18897, 0.56692, 2.83459])

上述代码中，看起来在 coord_orig 上有多余的代码；这是因为 \(\boldsymbol{R}_\mu\) 本身应当是原子核坐标相对于规范原点的距离；这里我们的规范原点恰好是原点 coord_orig。

下面我们就求取积分值 \(\nabla_{\boldsymbol{\mathscr{B}}} S_{36}^\mathrm{GIAO}\)：

\[ \nabla_{\boldsymbol{\mathscr{B}}} S_{36}^\mathrm{GIAO} = \int \phi_3 (\boldsymbol{r}) \left( - \frac{i}{2} \boldsymbol{R}_{36} \times \boldsymbol{r} \right) \phi_6 (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} \simeq - \frac{i}{2} \sum_g w_g \phi_3 (\boldsymbol{r}_g) \phi_6 (\boldsymbol{r}_g) \cdot \boldsymbol{R}_{36} \times \boldsymbol{r}_g \]

- 0.5j * np.einsum("g, g, g, gr -> r", grids.weights, ao[:, 3], ao[:, 6], np.cross(R_36, grids.coords))

array([0.+0.22657j, 0.-0.j     , 0.-0.0151j ])

上述结果是一个纯虚数向量，其三个维度分别相当于 \((\mathscr{B}_x, \mathscr{B}_y, \mathscr{B}_z)\) 的维度。

最后，我们指出，在 PySCF/libcint 中，上述积分的虚部负值可以通过输入 int1e_igovlp 字符串实现。其意义可以参考 auto_intor.cl 与 README。留意在 libcint 中，这类积分被表示为 \(i \langle \boldsymbol{U}_\mathrm{g} \mu | \nu \rangle\)，其 \(\boldsymbol{U}_\mathrm{g} = (U_\mathrm{g}^x, U_\mathrm{g}^y, U_\mathrm{g}^z)\) 的意义相当于三维向量算符

\[ \boldsymbol{U}_\mathrm{g} = - \frac{i}{2} \boldsymbol{R}_{\mu \nu} \times \boldsymbol{r} \]

mol.intor("int1e_igovlp")[:, 3, 6]

array([-0.22658, -0.     ,  0.01511])

PySCF 导出是实 (反对称) 矩阵 \(i \langle \boldsymbol{U}_\mathrm{g} \mu | \nu \rangle\)；因此在实际使用该积分时，我们就需要乘以 \(- i\)，让该矩阵变为复厄米矩阵，成为真正的 GIAO 变换算符一阶的矩阵形式 \(\langle \boldsymbol{U}_\mathrm{g} \mu | \nu \rangle\)。

- 1j * mol.intor("int1e_igovlp")[:, 3, 6]

array([0.+0.22658j, 0.+0.j     , 0.-0.01511j])

后文我们也会用 \(\boldsymbol{U}_\mathrm{g}\) 来简化一阶导数下的 GIAO 变换算符。

当然，我们需要知道这并不是一个很严格的写法，因为 \(\boldsymbol{U}_\mathrm{g}\) 的等式右边与 \(\mu, \nu\) 有关；因此，使用这种简写时一定需要留意哪些轨道与该 GIAO 变换算符一阶导数有关。拿双电子积分 \(( U_\mathrm{g}^t \mu \lambda | \kappa \nu )\) 来说，由于 \(U_\mathrm{g}^t\) 作用在 \(\mu\) 上，并且 \(\lambda\) 与 \(\mu\) 在积分过程中使用相同的电子坐标，因此 \(U_\mathrm{g}^t\) 的作用对象就是 \(\mu\) 与 \(\lambda\)。

从偶极积分出发给出微扰的 GIAO 重叠矩阵

实际上，

\[\begin{split} \begin{align} \nabla_{\boldsymbol{\mathscr{B}}} S_{\mu \nu}^\mathrm{GIAO} &= \langle \mu | - \frac{i}{2} \boldsymbol{R}_{\mu \nu} \times \boldsymbol{r} | \nu \rangle \\ &= - \frac{i}{2} \boldsymbol{R}_{\mu \nu} \times \langle \mu | \boldsymbol{r} | \nu \rangle \end{align} \end{split}\]

我们仍然拿 \((\mu, \nu) = (3, 6)\) 的情形讨论。我们应当注意到，\(\langle \mu | \boldsymbol{r} | \nu \rangle\) 与偶极积分形式几乎一致 (依照不同的定义方式，偶极积分可能是该值或其相反数)，可以用字符串表示为 int1e_r。因此，\(\nabla_{\boldsymbol{\mathscr{B}}} S_{36}^\mathrm{GIAO}\) 还可以程序化为

- 0.5j * np.cross(R_36, mol.intor("int1e_r")[:, 3, 6])

array([0.+0.22658j, 0.-0.j     , 0.-0.01511j])

微扰的 GIAO 重叠矩阵并非“规范不变”

我们知道“规范不变”的意义是，作为最终结果的磁化率在任意的规范原点下值相同。但这并不意味着中间过程的矩阵或数值结果也相同。重叠矩阵就不满足这种“规范不变”性质：

- 1j * mol.intor("int1e_igovlp")[:, 3, 6]

array([0.+0.22658j, 0.+0.j     , 0.-0.01511j])

- 1j * mol_trans.intor("int1e_igovlp")[:, 3, 6]

array([0.-0.52235j, 0.-0.65814j, 0.+0.16645j])

PySCF 所使用的默认规范原点就是坐标系原点，因此我们这里都没有严格地使用 with mol.with_common_orig(coord_orig) 语句。当然，由于 GIAO 下的磁化率应当不受规范原点的选取变化而受到影响，因此仅从最终结果上来说，也不需要刻意地规定规范原点。

6.2. Core Hamiltonian 的程序实现#

6.2.1. 一阶 Hamiltonian Core#

前一篇文档中，所有的一阶微扰是

\[ \hat h {}^{(1)} (\boldsymbol{\mathscr{B}}) = \frac{1}{2} \boldsymbol{\mathscr{B}} \cdot \boldsymbol{r} \times \boldsymbol{\hat{p}} \]

但当考虑到 GIAO，一阶微扰算符则应当写作

\[ \hat h {}^{(1)} (\boldsymbol{\mathscr{B}}, \boldsymbol{U}_\mathrm{g}) | \nu \rangle = \frac{1}{2} \boldsymbol{\mathscr{B}} \cdot (\boldsymbol{r} - \boldsymbol{R}_\nu) \times \boldsymbol{\hat{p}} + \boldsymbol{\mathscr{B}} \cdot \boldsymbol{U}_\mathrm{g} \hat f {}^{(0)} | \nu \rangle \]

需要注意到，这里使用到了零阶 Fock 算符 \(\hat f {}^{(0)}\)，该算符需要代入未受外场微扰的密度矩阵才能获得，而不是单纯地由分子构型与基组就能决定的。该项贡献的来源是 \(\langle \Psi^{(0), \mathrm{GIAO}} | \hat H^{(0)} | \Psi^{(0), \mathrm{GIAO}} \rangle\)。

因此，顺磁项对应的积分可以公式表达为

\[\begin{split} \begin{align} h_{t \mu \nu}^{(1)} &= \frac{1}{2} \langle \mu | \hat l_t | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu} \\ & \quad + \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}\]

其中，积分字符与公式表达之间的关系为

int1e_cg_irxp \(i \langle \mu | \hat l_t | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu}\)
int1e_igkin \(i \langle \boldsymbol{U}_\mathrm{g} \mu | \hat t | \nu \rangle\)
int1e_ignuc \(i \langle \boldsymbol{U}_\mathrm{g} \mu | \hat v_\mathrm{nuc} | \nu \rangle\)
int2e_ig1 \(i ( \boldsymbol{U}_\mathrm{g} \mu \nu | \kappa \lambda )\)

需要留意 \(i\)，因为这会使得上式中很多项的正负号在编写程序时是相反的。

因此，一阶 Core Hamiltonian hcore_1 \(h_{t \mu \nu}^{(1)}\) 的表达式可以用下述程序表示：

dm_guess = mf.make_rdm1()
hcore_1 = 1j * (
    - 0.5 * mol.intor("int1e_giao_irjxp")
    - mol.intor("int1e_igkin")
    - mol.intor("int1e_ignuc")
    - np.einsum("tuvkl, kl -> tuv", mol.intor("int2e_ig1"), dm_guess)
    + 0.5 * np.einsum("tulkv, kl -> tuv", mol.intor("int2e_ig1"), dm_guess)
    + 0.5 * np.einsum("tkvul, kl -> tuv", mol.intor("int2e_ig1"), dm_guess)
)

PySCF 中有对应的函数，生成 GIAO 下的一阶 Core Hamiltonian 矩阵 (事实上上述程序块就是从下述函数中获得的)：

np.allclose(
    hcore_1,
    1j * nmr.rhf.make_h10(mol, dm_guess)
)

True

最后，我们指出，关于库伦积分，实际上应当表示为

\[ h_{t \mu \nu}^{(1)} \leftarrow - i \sum_{\kappa \lambda} ( U_\mathrm{g}^t \mu \nu | \kappa \lambda ) D_{\kappa \lambda}^{(0)} - i \sum_{\kappa \lambda} ( \mu \nu | U_\mathrm{g}^t \kappa \lambda ) D_{\kappa \lambda}^{(0)} \]

但由于反对称性质 \(( \mu \nu | U_\mathrm{g}^t \kappa \lambda ) = - ( \mu \nu | U_\mathrm{g}^t \lambda \kappa )\) (反映到复数的情形其实是复共轭)，导致它与对称的零阶密度矩阵 \(D_{\kappa \lambda}^{(0)}\) 相乘并求和后必为零值。

np.allclose(mol.intor("int2e_ig1"), - mol.intor("int2e_ig1").swapaxes(-3, -4))

True

在程序初步编写时，ERI 积分上很容易遇到角标顺序应该怎么写的问题。这是因为在实对称矩阵情形下，\(\sum_{\kappa \lambda} (\mu \kappa | \nu \lambda) D_{\kappa \lambda}\) 与 \(\sum_{\kappa \lambda} (\mu \lambda | \kappa \nu) D_{\kappa \lambda}\) 是完全相同的，但当 \(\boldsymbol{U}_g\) 作用于 ERI 后，这种性质就不满足了。一个比较容易达成正确程序编写与公式推导的技巧是，保证 \(\mu, \kappa\) 处在 ERI 积分的复共轭位上，而 \(\nu, \lambda\) 处在普通位上 (或者用物理的记号来讲，当需要交换 \(\langle \mu \kappa | \nu \lambda \rangle\) 的角标时，只能在竖线左或者竖线右相互交换，不能跨线)。

6.2.2. 二阶 Hamiltonian Core#

在前一篇文档中，所有的二阶微扰是

\[ \hat h {}^{(2)} (\boldsymbol{\mathscr{B}}) = \frac{1}{8} \big( \boldsymbol{\mathscr{B}}^2 \boldsymbol{r}^2 - (\boldsymbol{\mathscr{B}} \cdot \boldsymbol{r}) (\boldsymbol{\mathscr{B}} \cdot \boldsymbol{r})^\mathrm{T} \big) \]

但当考虑到 GIAO，二阶微扰算符则应当写作

\[\begin{split} \begin{align} \hat h {}^{(2)} (\boldsymbol{\mathscr{B}}, \boldsymbol{U}_\mathrm{g}) | \nu \rangle &= \frac{1}{8} \big( \boldsymbol{\mathscr{B}}^2 (\boldsymbol{r} - \boldsymbol{R}_\nu)^2 - \boldsymbol{\mathscr{B}}^\mathrm{T} \big( (\boldsymbol{r} - \boldsymbol{R}_\nu) (\boldsymbol{r} - \boldsymbol{R}_\nu)^\mathrm{T} \big) \boldsymbol{\mathscr{B}} | \nu \rangle \\ & \quad + \frac{i}{4} \boldsymbol{\mathscr{B}}^\mathrm{T} \big( (\boldsymbol{r} - \boldsymbol{R}_\nu) \times \boldsymbol{\hat{p}} \big) \boldsymbol{U}_\mathrm{g}^\mathrm{T} \boldsymbol{\mathscr{B}} | \nu \rangle + \frac{i}{4} \boldsymbol{\mathscr{B}}^\mathrm{T} \boldsymbol{U}_\mathrm{g} \big( (\boldsymbol{r} - \boldsymbol{R}_\nu) \times \boldsymbol{\hat{p}} \big)^\mathrm{T} \boldsymbol{\mathscr{B}} | \nu \rangle \\ & \quad + \frac{1}{2} \boldsymbol{\mathscr{B}}^\mathrm{T} \big( \boldsymbol{U}_\mathrm{g} \boldsymbol{U}_\mathrm{g}^\mathrm{T} \big) \boldsymbol{\mathscr{B}} \hat f^{(0)} | \nu \rangle \end{align} \end{split}\]

其中，\(\mathbf{T}\) 是向量算符的转置，这种转置会生成矩阵形式的磁化率。

落实到程序中，则可以写为

\[\begin{split} \begin{align} h_{ts \mu \nu}^{(2)} \cdot \mathscr{B}_t \mathscr{B}_s &= \frac{1}{8} \big( \delta_{ts} \langle \mu | x^2 + y^2 + z^2 | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu} - \langle \mu | ts | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu} \big) \\ &\quad + \frac{1}{4} \langle U_g^t \mu | \hat l_s | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu} + \frac{1}{4} \langle U_g^s \mu | \hat l_t | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu} \\ &\quad + \frac{1}{2} \langle U_g^t U_g^s \mu | \hat t | \nu \rangle + \frac{1}{2} \langle U_g^t U_g^s \mu | \hat v_\mathrm{nuc} | \nu \rangle \\ &\quad + \frac{1}{4} \sum_{\kappa \lambda} \big( (U_g^t U_g^s \mu \nu | \kappa \lambda) + (U_g^t \mu \nu | U_g^s \kappa \lambda) \big) D_{\kappa \lambda}^{(0)} + \frac{1}{4} \sum_{\kappa \lambda} \big( (U_g^t U_g^s \kappa \lambda | \mu \nu) + (U_g^t \kappa \lambda | U_g^s \mu \nu) \big) D_{\kappa \lambda}^{(0)} \\ &\quad - \frac{1}{8} \sum_{\kappa \lambda} \big( (U_g^t U_g^s \mu \lambda | \kappa \nu) + (U_g^t \mu \lambda | U_g^s \kappa \nu) \big) D_{\kappa \lambda}^{(0)} - \frac{1}{8} \sum_{\kappa \lambda} \big( (U_g^t U_g^s \kappa \nu | \mu \lambda) + (U_g^t \kappa \nu | U_g^s \mu \lambda) \big) D_{\kappa \lambda}^{(0)} \\ \end{align} \end{split}\]

其中，积分字符与公式表达之间的关系为

int1e_rr_origj \(\langle \mu | ts | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu}\)
int1e_grjxp \(\langle U_g^t \mu | \hat l_s | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu}\)
int1e_ggkin \(\langle U_g^t \mu | \hat t | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu}\)
int1e_ggnuc \(\langle U_g^s \mu | \hat v_\mathrm{nuc} | \nu \rangle_{\mathrm{Gauge} \rightarrow \boldsymbol{R}_\nu}\)
int2e_gg1 \((U_g^t U_g^s \mu \nu | \kappa \lambda)\)
int2e_g1g2 \((U_g^t \mu \nu | U_g^s \kappa \lambda)\)

int1e_rr_origj = mol.intor("int1e_rr_origj").reshape(3, 3, nao, nao)
int1e_grjxp = mol.intor('int1e_grjxp').reshape(3, 3, nao, nao)
int1e_ggkin = mol.intor('int1e_ggkin').reshape(3, 3, nao, nao)
int1e_ggnuc = mol.intor('int1e_ggnuc').reshape(3, 3, nao, nao)
int2e_gg1   = mol.intor("int2e_gg1")  .reshape(3, 3, nao, nao, nao, nao)
int2e_g1g2  = mol.intor("int2e_g1g2") .reshape(3, 3, nao, nao, nao, nao)
int2e_gg    = int2e_gg1 + int2e_g1g2

因此，二阶 Core Hamiltonian hcore_2 \(h_{ts \mu \nu}^{(2)}\) 的表达式可以用下述程序表示：

hcore_2 = (
    + 1/8 * (np.einsum("ts, uv -> tsuv", np.eye(3), int1e_rr_origj.diagonal(0, 0, 1).sum(-1)) - int1e_rr_origj)
    + 1/4 * mol.intor('int1e_grjxp').reshape(3, 3, nao, nao)
    + 1/4 * mol.intor('int1e_grjxp').reshape(3, 3, nao, nao).swapaxes(0, 1)
    + 1/2 * mol.intor('int1e_ggkin').reshape(3, 3, nao, nao)
    + 1/2 * mol.intor('int1e_ggnuc').reshape(3, 3, nao, nao)
    + 1/4 * np.einsum("tsuvkl, kl -> tsuv", int2e_gg, dm_guess)
    + 1/4 * np.einsum("tskluv, kl -> tsuv", int2e_gg, dm_guess)
    - 1/8 * np.einsum("tsulkv, kl -> tsuv", int2e_gg, dm_guess)
    - 1/8 * np.einsum("tskvul, kl -> tsuv", int2e_gg, dm_guess)
)

6.3. 重叠矩阵的程序实现#

在 GIAO 下，除了 Hamiltonian Core 会发生改变，重叠矩阵也同样会产生变化。

一阶重叠矩阵程序中表示为

\[ S_{t \mu \nu}^{(1)} = \langle U_g^t \mu | \nu \rangle \]

二阶重叠矩阵程序中表示为

\[ S_{t \mu \nu}^{(2)} = \frac{1}{2} \langle U_g^t U_g^s \mu | \nu \rangle \]

其中，积分字符与公式表达之间的关系为

int1e_igovlp \(i \langle U_g^t \mu | \nu \rangle\)
int1e_ggovlp \(\langle U_g^t U_g^s \mu | \nu \rangle\)

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

ovlp_2 = 0.5 * mol.intor('int1e_ggovlp').reshape(3, 3, mol.nao, mol.nao)

6.4. 数值导数求磁化率#

与上一篇文档一样，我们已经获得了原子轨道形式的一阶、二阶 Core Hamiltonian、重叠积分。通过在 scf.RHF 的自洽场过程中重载 (override) 函数 get_hcore 与 get_ovlp，就可以得到受外磁场微扰的分子能量 \(E_\mathrm{tot} (\mathscr{B}_x, \mathscr{B}_y, \mathscr{B}_z)\)：

dm_guess = mf.make_rdm1()

def hcore_mag_field(dev_xyz):
    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, dev_xyz)
        hcore_total += np.einsum("tsuv, t, s -> uv", hcore_2, dev_xyz, dev_xyz)
        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, dev_xyz)
        ovlp_total += np.einsum("tsuv, t, s -> uv", ovlp_2, dev_xyz, dev_xyz)
        return ovlp_total
    mf.get_hcore = get_hcore
    mf.get_ovlp  = get_ovlp
    return mf.kernel(dm=dm_guess)

我们可以用与上一篇文档相同的代码实现二阶梯度：

\[ \xi_{ts} = - \frac{\partial^2 E_\mathrm{tot} (\boldsymbol{\mathscr{B}})}{\partial \mathscr{B}_t \partial \mathscr{B}_s} \]

eng_origin = hcore_mag_field((0, 0, 0))
interval = 1e-4
num_polar = np.zeros((3, 3))
for t in range(3):
    for s in range(3):
        if t != s:
            dev_xyzs = np.zeros((4, 3))
            dev_xyzs[0, t] = dev_xyzs[0, s] = dev_xyzs[1, t] = dev_xyzs[2, s] =  interval
            dev_xyzs[3, t] = dev_xyzs[3, s] = dev_xyzs[2, t] = dev_xyzs[1, s] = -interval
            num_polar[t, s] = (
                + hcore_mag_field(dev_xyzs[0])
                - hcore_mag_field(dev_xyzs[1])
                - hcore_mag_field(dev_xyzs[2])
                + hcore_mag_field(dev_xyzs[3])
            ) / (4 * interval**2)
        else:
            dev_xyzs = np.zeros((2, 3))
            dev_xyzs[0, t], dev_xyzs[1, t] = interval, -interval
            num_polar[t, t] = (
                + hcore_mag_field(dev_xyzs[0])
                + hcore_mag_field(dev_xyzs[1])
                - eng_origin * 2
            ) / (interval ** 2)
num_polar *= -1

num_polar

array([[-3.72496, -0.02158, -0.07515],
       [-0.02158, -2.88877,  0.20344],
       [-0.07514,  0.20345, -3.57484]])

我们最后与 PySCF 所给出的结果进行核验：

mf_mag.kernel()

array([[-3.72497, -0.02158, -0.07514],
       [-0.02158, -2.88877,  0.20345],
       [-0.07514,  0.20345, -3.57485]])

磁性质数值导数 (2)：RHF 的 GIAO 磁化率

文档目录

6. 磁性质数值导数 (2)：RHF 的 GIAO 磁化率#

6.1. GIAO 规范不变原子轨道#

6.1.1. 基本概念#

6.1.2. 实例分析：GIAO 重叠矩阵在外磁场下的微扰#

6.2. Core Hamiltonian 的程序实现#

6.2.1. 一阶 Hamiltonian Core#

6.2.2. 二阶 Hamiltonian Core#

6.3. 重叠矩阵的程序实现#

6.4. 数值导数求磁化率#