<jupyter_start><jupyter_text>GIAO 的 RHF 磁化率数值导数计算方式简述 > 创建时间：2020-08-30 在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算 GIAO 的 RHF 数值磁化率的程序。该文档大量参考 PySCF 的代码 [magnetizability/rhf.py](https://github.com/pyscf/pyscf/blob/master/pyscf/prop/magnetizability/rhf.py) 与 [nmr/rhf.py](https://github.com/pyscf/pyscf/blob/master/pyscf/prop/nmr/rhf.py)。一篇公式记号比较清晰的文章是 Laasner, Blum, et al. [^Laasner-Blum.arXiv.2018]。与上一篇文档一样，我们的讨论中所使用到的分子体系 `mol` 会是非对称的氨分子，并且取用最小基组。其 RHF 计算放在实例 `mf`，而磁化率计算实例会放在 `mf_mag`。在这篇文档中，我们仍然需要保留规范原点 (Gauge Origin) 的概念。规范原点会取在坐标原点上 `coord_orig`。<jupyter_code>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<jupyter_output><empty_output><jupyter_text>其自洽场能量为<jupyter_code>mf = scf.RHF(mol).run()
mf.e_tot<jupyter_output><empty_output><jupyter_text>我们定义磁化率计算类为 `mf_mag`。其磁化率张量 $\xi_{ts}$ 为：<jupyter_code>mf_mag = magnetizability.RHF(mf)
mf_mag.kernel()<jupyter_output><empty_output><jupyter_text>但留意到，与上一份文档不同地，如果我们对分子作平移操作 (譬如下述平移操作是将原点移动 $(x, y, z) = (10, -10, 5) \, \mathsf{a.u.}$)，其磁化率仍然保持不变：<jupyter_code>mol_trans = mol.copy()
mol_trans.set_geom_(mol.atom_coords() + np.array([[10, -10, 5]]), unit="AU")
mol_trans.build()
magnetizability.RHF(scf.RHF(mol_trans).run()).kernel()<jupyter_output><empty_output><jupyter_text>GIAO 规范不变原子轨道  基本概念 我们从上一篇文档中，得知对分子作平移操作后，磁化率可能会改变。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 [^Pople-Pople.MP.1958] eq 2.5 附近的讨论。事实上，规范 (原点) **不变** 原子轨道的称呼可能是不恰当的。Pople [^Pople-Pople.DFS.1962.34] 曾在文章里使用了规范 (原点) **依赖** 原子轨道 (Gauge Dependent Atomic Orbital)；这是因为 $\phi^\mathrm{GIAO}_\mu (\boldsymbol{r})$ 实际上是依规范原点位置不同而不同的 (注意 $\boldsymbol{R}_\mu$ 的定义)。当然，我们仍然沿用 GIAO 的约定俗成。  实例分析：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)$<jupyter_code>ni = dft.numint.NumInt()

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

ao = ni.eval_ao(mol, grids.coords)
ao.shape<jupyter_output><empty_output><jupyter_text>**数值格点积分** 我们预先已经知道，第 $\mu = 3$ 根基轨道是第 0 个原子 (N 原子)，第 $\nu = 6$ 跟轨道是第 2 个原子 (H 原子) (按照 0-index 计数)，那么向量 $\boldsymbol{R}_{\mu \nu} = \boldsymbol{R}_{36} = \boldsymbol{R}_{6} - \boldsymbol{R}_{3}$：<jupyter_code>R_36 = (mol.atom_coord(2) - coord_orig) - (mol.atom_coord(0) - coord_orig)
R_36<jupyter_output><empty_output><jupyter_text>上述代码中，看起来在 `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$$<jupyter_code>- 0.5j * np.einsum("g, g, g, gr -> r", grids.weights, ao[:, 3], ao[:, 6], np.cross(R_36, grids.coords))<jupyter_output><empty_output><jupyter_text>上述结果是一个纯虚数向量，其三个维度分别相当于 $(\mathscr{B}_x, \mathscr{B}_y, \mathscr{B}_z)$ 的维度。 最后，我们指出，在 PySCF/libcint 中，上述积分的虚部负值可以通过输入 `int1e_igovlp` 字符串实现。其意义可以参考 [auto_intor.cl](https://github.com/sunqm/libcint/blob/abf6948fa17e5b4ecbd26de05bf4b1d7b2b2fe3c/scripts/auto_intor.clL17) 与 [README](https://github.com/sunqm/libcint/blob/master/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}$$<jupyter_code>mol.intor("int1e_igovlp")[:, 3, 6]<jupyter_output><empty_output><jupyter_text>PySCF 导出是实 (反对称) 矩阵 $i \langle \boldsymbol{U}_\mathrm{g} \mu | \nu \rangle$；因此在实际使用该积分时，我们就需要乘以 $- i$，让该矩阵变为复厄米矩阵，成为真正的 GIAO 变换算符一阶的矩阵形式 $\langle \boldsymbol{U}_\mathrm{g} \mu | \nu \rangle$。<jupyter_code>- 1j * mol.intor("int1e_igovlp")[:, 3, 6]<jupyter_output><empty_output><jupyter_text>后文我们也会用 $\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{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}$$ 我们仍然拿 $(\mu, \nu) = (3, 6)$ 的情形讨论。我们应当注意到，$\langle \mu | \boldsymbol{r} | \nu \rangle$ 与偶极积分形式几乎一致 (依照不同的定义方式，偶极积分可能是该值或其相反数)，可以用字符串表示为 `int1e_r`。因此，$\nabla_{\boldsymbol{\mathscr{B}}} S_{36}^\mathrm{GIAO}$ 还可以程序化为<jupyter_code>- 0.5j * np.cross(R_36, mol.intor("int1e_r")[:, 3, 6])<jupyter_output><empty_output><jupyter_text>**微扰的 GIAO 重叠矩阵并非“规范不变”** 我们知道“规范不变”的意义是，作为最终结果的磁化率在任意的规范原点下值相同。但这并不意味着中间过程的矩阵或数值结果也相同。重叠矩阵就不满足这种“规范不变”性质：<jupyter_code>- 1j * mol.intor("int1e_igovlp")[:, 3, 6]
- 1j * mol_trans.intor("int1e_igovlp")[:, 3, 6]<jupyter_output><empty_output><jupyter_text>PySCF 所使用的默认规范原点就是坐标系原点，因此我们这里都没有严格地使用 `with mol.with_common_orig(coord_orig)` 语句。当然，由于 GIAO 下的磁化率应当不受规范原点的选取变化而受到影响，因此仅从最终结果上来说，也不需要刻意地规定规范原点。  Core Hamiltonian 的程序实现  一阶 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{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}$$ 其中，积分字符与公式表达之间的关系为- `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)}$ 的表达式可以用下述程序表示：<jupyter_code>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)
)<jupyter_output><empty_output><jupyter_text>PySCF 中有对应的函数，生成 GIAO 下的一阶 Core Hamiltonian 矩阵 (事实上上述程序块就是从下述函数中获得的)：<jupyter_code>np.allclose(
    hcore_1,
    1j * nmr.rhf.make_h10(mol, dm_guess)
)<jupyter_output><empty_output><jupyter_text>最后，我们指出，关于库伦积分，实际上应当表示为$$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)}$ 相乘并求和后必为零值。<jupyter_code>np.allclose(mol.intor("int2e_ig1"), - mol.intor("int2e_ig1").swapaxes(-3, -4))<jupyter_output><empty_output><jupyter_text>在程序初步编写时，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$ 的角标时，只能在竖线左或者竖线右相互交换，不能跨线)。  二阶 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{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}$$ 其中，$\mathbf{T}$ 是向量算符的转置，这种转置会生成矩阵形式的磁化率。 落实到程序中，则可以写为$$\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}$$ 其中，积分字符与公式表达之间的关系为- `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)$<jupyter_code>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<jupyter_output><empty_output><jupyter_text>因此，二阶 Core Hamiltonian `hcore_2` $h_{ts \mu \nu}^{(2)}$ 的表达式可以用下述程序表示：<jupyter_code>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)
)<jupyter_output><empty_output><jupyter_text>重叠矩阵的程序实现 在 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$<jupyter_code>ovlp_1 = - 1j * mol.intor("int1e_igovlp")
ovlp_2 = 0.5 * mol.intor('int1e_ggovlp').reshape(3, 3, mol.nao, mol.nao)<jupyter_output><empty_output><jupyter_text>数值导数求磁化率 与上一篇文档一样，我们已经获得了原子轨道形式的一阶、二阶 Core Hamiltonian、重叠积分。通过在 `scf.RHF` 的自洽场过程中重载 (override) 函数 `get_hcore` 与 `get_ovlp`，就可以得到受外磁场微扰的分子能量 $E_\mathrm{tot} (\mathscr{B}_x, \mathscr{B}_y, \mathscr{B}_z)$：<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>我们可以用与上一篇文档相同的代码实现二阶梯度：$$\xi_{ts} = - \frac{\partial^2 E_\mathrm{tot} (\boldsymbol{\mathscr{B}})}{\partial \mathscr{B}_t \partial \mathscr{B}_s}$$<jupyter_code>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<jupyter_output><empty_output><jupyter_text>我们最后与 PySCF 所给出的结果进行核验：<jupyter_code>mf_mag.kernel()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.21029968997588702, 'fraction_numerical': 0.023711103456194742, 'mean_word_length': 4.662873862158648, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 95, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '122272', 'license': 'MIT', 'n_tokens_mistral': 9261, 'n_tokens_neox': 9452, 'n_words': 1301, 'path': 'source/QC_Notes/Prop_Series/Mag_GIAO_NumDeriv.ipynb', 'repo_name': 'ajz34/ajz34.readthedocs.io', 'size': '17418'}