5. 磁性质数值导数 (1)：RHF 的非 GIAO 磁化率#
创建时间：2020-08-27
在这篇文档中，我们会讨论使用 PySCF 以及其作为 libcint 的接口，计算非 GIAO 的 RHF 数值磁化率的程序。该文档大量参考 PySCF 的代码 magnetizability/rhf.py 与 nmr/rhf.py。一些公式符号参考 Atkins, Friedman [1]。
我们的讨论中所使用到的分子体系 mol 会是非对称的氨分子，并且取用最小基组。规范原点 (Gauge Origin) 会取在坐标原点上 coord_orig。其 RHF 计算放在实例 mf，而磁化率计算实例会放在 mf_mag。
from pyscf import gto, scf
from pyscf.prop import nmr, magnetizability
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)
其自洽场能量为
mf = scf.RHF(mol).run()
mf.e_tot
-55.253540514686556
其磁化率张量 \(\xi_{ts}\) 为 (其中，\(t, s \in \{ x, y, z \}\) 表示三个坐标方向，需要注意这里选择了规范原点为坐标原点，若选取其它坐标则会得到非常不同的结果)
mf_mag = magnetizability.RHF(mf)
mf_mag.gauge_orig = coord_orig
mf_mag.kernel()
array([[-4.94475, 0.21773, -0.08268],
[ 0.21773, -4.27801, 0.49885],
[-0.08268, 0.49885, -4.15348]])
5.1. 基础概念#
5.1.1. 分子能量作为外加微扰量的函数#
我们指出，磁化率可以看作是分子处在某一恒定外加磁场 \(\boldsymbol{\mathscr{B}}\) 下 (作为三维矢量)，所产生的能量变化的表征：
以一般物理的约定俗成而言，对于外加的磁场微扰 (Atkins and Friedman, eq 13.34)
其中，\(\boldsymbol{\xi}\) 是二维对称矩阵 (或称张量，如之前代码所展示)。我们使用了 \(\boldsymbol{\xi}\) (Atkins and Friedman, eq 13.34, termed as magnetizability) 而非 \(\boldsymbol{\chi}\) (Atkins and Friedman, eq 13.3c, termed as magnetic susceptibility) 来表示磁化率。因此，磁化率本身可以表示为 (矩阵元的形式与向量形式)
5.1.2. 哈密顿算符作为外加微扰量的算符#
能量可以通过波函数在哈密顿算符的变分极小值处的期望获得：
其中，
上述算符是体系的多电子总哈密顿算符；而 \(\hat h (\boldsymbol{\mathscr{B}})\) 则是单电子的 Core Hamiltonian 算符；\(\hat V_\mathrm{ee}\) 为电子互斥算符，\(\hat V_\mathrm{NN}\) 为原子核互斥算符。需要注意，由于我们不使用 GIAO，因此 \(\hat V_\mathrm{ee}\) 就是普通的电子互斥算符，不受外场 \(\boldsymbol{\mathscr{B}}\) 干扰；但使用 GIAO 的情况下，可能需要额外考虑这部分贡献。
\(\hat h {}^{(0)}\) 是没有外加场的算符 (这与自洽场计算过程所用到的算符相同)。其余的算符则为 (Atkins and Friedman, eq 13.26, eq 13.29)
其中，
其中，
前一项被称为顺磁项 (Paramagnetic)，后一项称为抗磁项 (Diamagnetic)。\(\Psi^{(0)}\) 是指未微扰的体系哈密顿算符 \(\hat H {}^{(0)}\) 的本征态，\(\Psi^{(1)} (\boldsymbol{\mathscr{B}})\) 则是一阶微扰的波函数；其解析的求取方法是在程序中表示为 U 矩阵，通过 CP-HF 方程求取；我们在这里不会对解析方法作说明，但了解这两项的区分是有帮助的。
5.2. Core Hamiltonian 的程序实现#
我们知道，PySCF 中，在自洽场实例中更改 Core Hamiltonian 的类方法函数 (method function) 就可以实现外场微扰下的能量计算。这在 pyxdh 偶极矩的计算 文档 中有所说明。在这里我们也要做类似的工作。
5.2.1. 顺磁项#
在 PySCF 中，顺磁项 \(\hat h {}^{(1)} (\boldsymbol{\mathscr{B}}) = \frac{1}{2} \boldsymbol{\mathscr{B}} \cdot \boldsymbol{r} \times \boldsymbol{\hat p}\) 有其对应的积分 hcore_1 (\(h_{t \mu \nu}^{(1)}\)，需要注意它不包含作为标量的 \(\mathscr{B}_t\))
hcore_1 = - 0.5 * mol.intor("int1e_cg_irxp") * 1j
hcore_1.shape, hcore_1.dtype
((3, 8, 8), dtype('complex128'))
上述的程序看起来会有些奇怪，因为这里出现了复数。我们需要分段对其作解释。
积分字符
我们使用到了积分字符 int1e_cg_irxp。关于这段字符，其意义需要通过 auto_intor.cl 程序了解：
'("int1e_cg_irxp" (#C(0 1) \| rc cross p))
其右侧是积分的具体形式，说明在 README 文件中，意义为
其维度是 \((t, \mu, \nu)\)，但其第一个维度是通过向量叉乘给出，因此它与 \(\boldsymbol{r}\) 或 \(\boldsymbol{p}\) 的维度不是直接相关的。如果我们令动量算符 \(\boldsymbol{\hat l} = \boldsymbol{r} \times \boldsymbol{\hat p}\)，那么可以将上述积分写为
反对称性厄米性
我们应当留意到 \(\mathtt{int1e\_cg\_irxp}_{t \mu \nu}\) 是一个反对称矩阵 (即随 \(\mu, \nu\) 交换成相反值)
np.allclose(mol.intor("int1e_cg_irxp"), - mol.intor("int1e_cg_irxp").swapaxes(-1, -2))
True
这是由于 \(\nabla\) 算符本身是一个反对称算符。但是需要留意到，动量算符在此基础上乘以了虚数单位 \(- i\)，因此，该矩阵是厄米的，即其转置后的共轭是其本身。我们所定义的 \(h_{t \mu \nu}^{(1)}\) 就具有这样的性质：
np.allclose(hcore_1, hcore_1.swapaxes(-1, -2).conj())
True
因此，我们会说 hcore_1 \(h_{t \mu \nu}^{(1)}\) 是厄米的。
5.2.2. 抗磁项#
抗磁项 \(\hat h {}^{(2)} (\boldsymbol{\mathscr{B}}) = \frac{1}{8} \big( \boldsymbol{\mathscr{B}}^2 \boldsymbol{r}^2 - (\boldsymbol{\mathscr{B}} \cdot \boldsymbol{r})^2 \big)\) 需要一些技巧生成。PySCF 中可以生成张量 int1e_rr：
我们定义 hcore_2 \(h_{ts \mu \nu}^{(2)}\) 为
with mol.with_common_orig(coord_orig):
int1e_rr = mol.intor("int1e_rr").reshape(3, 3, mol.nao, mol.nao)
hcore_2 = 1/8 * (np.einsum("ts, uv -> tsuv", np.eye(3), int1e_rr.diagonal(0, 0, 1).sum(-1)) - int1e_rr)
并且，上述张量具有下述性质：
5.2.3. Core Hamiltonian 程序实现#
我们最后可以编写外加磁场微扰下的 Core Hamiltonian，以及在此微扰下的分子体系能量。为了加快计算速度，我们会使用为微扰的自洽场密度作为初猜 dm_guess。Core Hamiltonian 表达式为
dm_guess = mf.make_rdm1()
def hcore_mag_field(dev_xyz):
mf = scf.RHF(mol)
def hcore(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
mf.get_hcore = hcore
return mf.kernel(dm=dm_guess)
上述函数的参数 t, s 表示坐标方向分量，dev_t, dev_s 表示外加微扰大小，单位为 a.u.。
譬如，若在 \(x\) 方向的磁场上施加 \(\mathscr{B}_x = 1 \, \mathsf{a.u.}\)，而 \(y\) 方向上施加 \(\mathscr{B}_y = 2 \, \mathsf{a.u.}\) (即 \(\boldsymbol{\mathscr{B}} = (\mathscr{B}_x, \mathscr{B}_y, \mathscr{B}_z) = (1, 2, 0) \, \mathsf{a.u.}\))，那么下述程序会给出该自洽场能量 \(E_\mathrm{tot} (\mathscr{B}_x, \mathscr{B}_y, \mathscr{B}_z)\)：
hcore_mag_field(0, 1, 1, 2)
-57.997213868466154
5.3. 数值导数求取磁化率#
我们已经有了求取 \(E_\mathrm{tot} (\mathscr{B}_x, \mathscr{B}_y, \mathscr{B}_z)\) 的程序了，接下来就可以进行数值导数计算。数值导数的计算公式可以简单地使用三点差分法；对于被求导量 \(x :\neq y\)，有 (当 \(h\) 足够小时)
而对被求导量相同的情形，有
下面的程序就依照上述两个公式进行二阶数值导数求取。求导的原点取在 \((\mathscr{B}_x, \mathscr{B}_y, \mathscr{B}_z) = (0, 0, 0)\) 即不受外磁场影响的情形的自洽场能量 eng_origin，差分大小为 interval \(h = 10^{-3} \, \mathsf{a.u.}\)。需要注意，根据约定俗成，
因此求取得到的磁化率 num_polar \(\xi_{ts}\) 需要乘以 -1。
eng_origin = hcore_mag_field((0, 0, 0))
interval = 1e-3
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[2])
- eng_origin * 2
) / (interval ** 2)
num_polar *= -1
num_polar
array([[-4.94475, 0.21773, -0.08268],
[ 0.21773, -4.27801, 0.49885],
[-0.08268, 0.49885, -4.15348]])
我们再与 PySCF 的解析结果作对照：
mf_mag.kernel()
array([[-4.94475, 0.21773, -0.08268],
[ 0.21773, -4.27801, 0.49885],
[-0.08268, 0.49885, -4.15348]])