3.2. RHF 自洽场计算

从这一节开始，我们将真正地开始量化方法的学习。最为基础、但也最为重要的量化方法是 RHF 自洽场计算。我们依靠 RHF 的实现过程来了解公式记号，以及电子积分、基组在实际问题中的应用。

同时，PySCF 作为易于扩充的量化程序包，实现了非常庞大的量化方法；RHF 是其中的一个。我们将简单地了解 PySCF 对 RHF 的实现、以及其中的一些函数接口。

import numpy as np
import scipy
from pyscf import scf, gto

from pkg_resources import resource_filename
from pyxdh.Utilities import FormchkInterface
from pyxdh.Utilities.test_molecules import Mol_H2O2

from functools import partial
np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])

np.set_printoptions(5, linewidth=150, suppress=True)

任务 (1)

1. 上面的代码块是引入外部程序以及作文档初始化的代码块；请解释上述代码的每一行的意义。

   不同的文档会有不同的初始化代码块，即使这些代码块可能看起来一样。请在阅读一份新的文档之前检查第一个代码块与其它文档是否有不同。

3.2.1. 量化软件的自洽场计算

3.2.1.1. PySCF 计算

我们先作一个 PySCF 的自洽场计算。我们以后会一直使用下面的 \(C_1\) 对称的 6-31G 基组的双氧水分子 mol 作为范例：

mol = gto.Mole()
mol.atom = """
O  0.0  0.0  0.0
O  0.0  0.0  1.5
H  1.0  0.0  0.0
H  0.0  0.7  1.0
"""
mol.basis = "6-31G"
mol.verbose = 0
mol.build()

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

而 PySCF 的 RHF 自洽场可以通过下述代码实现：

scf_eng = scf.RHF(mol)
scf_eng.conv_tol = 1e-12
scf_eng.conv_tol_grad = 1e-10
scf_eng.kernel()

-150.58503378083853

注意

1. PySCF 的自洽场、CP-HF 方程等迭代求解在失败的情况下，也不会抛出异常。如果将 mol.verbose 设到默认值，就可以看到一些警告信息。或者，我们通过下述语句来判断自洽场过程是否确实收敛。

   由于刚才给出的自洽场收敛条件过于苛刻，因此即使自洽场并没有收敛，其结果仍然可以用于定量分析。

scf_eng.converged

False

注意

在将来，为了简化文档的初始化代码，我们会使用 Mol_H2O2 类来初始化 mol 变量与 DFT 的格点变量 grids。其调用方式通过下述代码单元给出；其中，mol 变量的参数与方才定义的 mol 完全一致，是 6-31G 的 \(C_1\) 对称的分子；而 grids 参数则是

• 格点数量：径向 99，球面 590 的 Lebedev 格点；
• 截断模式：NWChem 模式 (PySCF 默认)
• 权重分配模式：Stratmann 模式 (Gaussian 默认)

mol = Mol_H2O2().mol
grids = Mol_H2O2().gen_grids()

3.2.1.2. Gaussian 计算

我们也可以使用 Gaussian 计算得到相同的结果。Gaussian 的计算结果储存在 pyxdh 库的资源文件夹下，其调用方式是：

ref_fchk = FormchkInterface(resource_filename("pyxdh", "Validation/gaussian/H2O2-HF-freq.fchk"))

上述语句会读取资源文件 H2O2-HF-freq.fchk。生成上述 formchk 文件的输入卡则在 H2O2-HF-freq.gjf；其内容是

%chk=H2O2-HF-freq
#p RHF/6-31G nosymm SCF(VeryTight) Freq

H2O2 HF Frequency

0 1
O  0.0  0.0  0.0
O  0.0  0.0  1.5
H  1.0  0.0  0.0
H  0.0  0.7  1.0

Gaussian 计算给出的 RHF 自洽场能量用下述方法调取：

ref_fchk.total_energy()

-150.5850337850876

可见，PySCF 所计算得到的 RHF 能量与 Gaussian 几乎一样。

3.2.1.3. pyxdh 计算

pyxdh 也可以进行 RHF 计算；但实际上，pyxdh 的 RHF 计算就是通过 PySCF 的接口实现的。这里仅仅介绍 pyxdh 程序在计算 RHF 时的接口，以为后续梯度的计算作准备。

我们已经在 B2PLYP 型泛函极化率的计算 一小节中介绍了 pyxdh 的计算方式，这里实际上是类似的，区别仅仅是没有引入 DFT：

from pyxdh.DerivOnce import GradSCF
config = {"scf_eng": scf_eng}
scfh = GradSCF(config)

scfh.eng

-150.58503378083853

其中，引入 GradSCF 与 DipoleSCF 在当前没有讨论梯度性质时是是没有区别的。可以发现其 RHF 给出的能量与 PySCF 和 Gaussian 的是近乎相同的。

事实上，pyxdh 仅仅是读入了 PySCF 已经给出的现成的 scf_eng 类，它不进行实际的自洽场计算。pyxdh 的便利之处是提供一些梯度计算中常用到的矩阵与张量，这些便利之处在以后才会体现。

3.2.2. 小型自洽场程序

3.2.2.1. 程序代码

自洽场程序会是通常量化课程的大作业。在这里，我们借用 PySCF 已经提供的函数接口来实现小的自洽场程序。这个程序是根据 Szabo and Ostlund [SO96] 书籍上的指示 (p146) 进行最简单的 SCF 程序编写。

提醒

1. 由于该分子已经难以通过零密度初猜来得到能量，因此这里暂且利用了 PySCF 所提供的默认初猜。

# Block 1

S = mol.intor("int1e_ovlp")
HC = mol.intor("int1e_kin") + mol.intor("int1e_nuc")
eri = mol.intor("int2e")
X = scipy.linalg.fractional_matrix_power(S, -0.5)
# X = np.linalg.inv(np.linalg.cholesky(S).T)

natm = mol.natm
nmo = nao = mol.nao
nocc = mol.nelec[0]
so = slice(0, nocc)

任务 (2)

1. (可选) 我们在上面的代码单元中已经生成了电子积分，随后我们的工作仅仅是给出一个 SCF 算法。这无外乎 Python 与 numpy 的代码书写。你完全可以尝试不看下面的代码，自己先试写一个 SCF 过程。这可以是物化研究生一年级的量化程序大作业。你可能必须要一个更好的初猜；零密度初猜似乎对双氧水分子不适用。密度矩阵的初猜可以通过以下代码得到：

   D = scf_eng.get_init_guess()
   

   上面的代码单元中，积分 int1e_ovlp、int1e_kin、int2e 的意义在 电子积分 一小节中提到；int1e_nuc 表示的是电子与原子核的外势积分。

   提示：你可以了解 np.linalg.eigh 函数的意义，并思考它可能对 SCF 过程的编写有何帮助。

# Block 2

A_t = mol.atom_coords()
r_ABt = A_t[:, None, :] - A_t[None, :, :]
r_AB = np.linalg.norm(r_ABt, axis=-1)
r_AB += np.diag(np.ones(mol.natm) * np.inf)
A_charge = mol.atom_charges()[:, None] * mol.atom_charges()
E_nuc = 0.5 * (A_charge / r_AB).sum()
print("Neucleus energy   ", E_nuc, " a.u.")

Neucleus energy    37.884674408641274  a.u.

# Block 3

D = scf_eng.get_init_guess()
D_old = np.zeros((nao, nao))
count = 0

while (not np.allclose(D, D_old)):
    if count > 500:
        raise ValueError("SCF not converged!")
    count += 1
    D_old = D
    F = HC + np.einsum("uvkl, kl -> uv", eri, D) - 0.5 * np.einsum("ukvl, kl -> uv", eri, D)
    Fp = X.T @ F @ X
    e, Cp = np.linalg.eigh(Fp)
    C = X @ Cp
    D = 2 * C[:, so] @ C[:, so].T

E_elec = (HC * D).sum() + 0.5 * np.einsum("uvkl, uv, kl ->", eri, D, D) - 0.25 * np.einsum("ukvl, uv, kl ->", eri, D, D)
E_tot = E_elec + E_nuc

print("SCF Converged in  ", count, " loops")
print("Electronic energy ", E_elec, " a.u.")
print("Total energy      ", E_tot, " a.u.")
print("----------------- ")
print("Energy allclose   ", np.allclose(E_tot, scf_eng.e_tot))
print("Density allclose  ", np.allclose(D, scf_eng.make_rdm1()))

SCF Converged in   124  loops
Electronic energy  -188.46970818948094  a.u.
Total energy       -150.58503378083967  a.u.
-----------------
Energy allclose    True
Density allclose   True

从 Total energy 的输出来看，我们能发现我们给出了正确的 RHF 能量。下面我们对代码进行说明。

3.2.2.2. 记号说明

在 Block 1 中，我们定义了三个电子积分、\(X_{\mu \nu}\) 矩阵以及与维度有关的量。除去分子轨道数，其余都是只与分子和基组有关的量。而一般来说，只要没有原子轨道线性依赖的情况，一般的程序都会定义分子轨道数与原子轨道基组数一致。

• S \(S_{\mu \nu}\)，或 int1e_ovlp 指交换积分，其在 Szabo (3.136) 定义

• HC \(h_{\mu \nu}\) 为动能积分 int1e_kin 与核排斥积分 int1e_nuc 之和，其在 Szabo (3.149) 定义

• eri \((\mu \nu | \kappa \lambda)\) 或 int2e 指双电子互斥积分，其在 Szabo (Table 2.2) 定义，采用 Chemistry Convention

• natm 表示原子数

• nmo 表示分子轨道数，以后默认与原子轨道数相等，但一般地，根据表达式总应当能区分我们应该采用原子轨道还是分子轨道

• nao 表示原子轨道数

• nocc 为占据轨道数；以后会出现 nvir，为非占轨道数

• so 为占据分子轨道分割；以后会出现 sv 为非占分子轨道分割，以及 sa 为全部分子轨道分割

• X 只在自洽场过程中出现，以后将不再使用；但会对该记号赋予新的意义 (密度的 U 偏导有关量)。在这一节中，其表达式为 \(X_{\mu \nu}\)，并满足关系式 (Szabo 3.165)

  \[X_{\kappa \mu} S_{\kappa \lambda} X_{\lambda \nu} = \delta_{\mu \nu}\]

记号说明

• \(\mu, \nu, \kappa, \lambda\) 代表原子轨道
• \(i, j, k, l\) 代表分子轨道，但出于程序编写需要 \(k, l\) 尽量不应与 \(\kappa, \lambda\) 同时出现
• \(a, b, c, d\) 代表非据轨道
• \(p, q, r, s, m\) 代表全分子轨道，但 \(r, s\) 的使用需要尽量避免，因与下述坐标分量记号冲突
• \(t, s, r, w, x\) 代表坐标分量；一般 \(t, s\) 特指原子坐标分量，\(r, w, x\) 特指电子坐标分量；坐标分量的三种可能取向是 \(x, y, z\) 方向
• \(A, B, M\) 代表原子；其中 \(M\) 一般是被求和的角标
• \(\boldsymbol{A}, \boldsymbol{B}, \boldsymbol{M}\) 代表原子坐标的三维向量，区别于普通斜体字母
• \(A_t, B_s\) 代表原子坐标的坐标分量，区别于 \(\boldsymbol{A}, \boldsymbol{B}\) 作为向量，也区别于 \(t, s\) 单纯地是坐标分量

任务 (3)

1. 尝试验证关系式 \(X_{\kappa \mu} S_{\kappa \lambda} X_{\lambda \nu} = \delta_{\mu \nu}\)。

在代码单元 Block 2 中，我们计算了核排斥能。其中，

• r_AB \(r_{AB}\) 表示原子间的欧氏距离，在 构型量输出 一小节中提及：

  \[\begin{split}\begin{equation} r_{AB} = \begin{cases} \Vert \boldsymbol{A} - \boldsymbol{B} \Vert_2 & A \neq B \\ + \infty & A = B \end{cases} \end{equation}\end{split}\]

• A_charge 表示量原子的电荷乘积：

  \[Z_{AB} = Z_A Z_B\]

• E_nuc 为原子核排斥能，以 a.u. 为单位：

  \[E_\mathrm{nuc} = \frac{1}{2} Z_{AB} r_{AB}^{-1}\]

3.2.2.3. 代码说明

随后我们考虑实际执行计算的 Block 3。

• Line 3

  D = scf_eng.get_init_guess()
  

  是除了电子积分外唯一使用 PySCF 的代码，它给一个合理的初猜 D \(D_{\mu \nu}\)。

• Line 12

  F = HC + np.einsum("uvkl, kl -> uv", eri, D) - 0.5 * np.einsum("ukvl, kl -> uv", eri, D)
  

  定义了 Fock 矩阵 F

  \[F_{\mu \nu} [D_{\kappa \lambda}] = h_{\mu \nu} + (\mu \nu | \kappa \lambda) D_{\kappa \lambda} - \frac{1}{2} (\mu \lambda| \kappa \nu) D_{\kappa \lambda}\]

• Line 16

  D = 2 * C[:, so] @ C[:, so].T
  

  通过使用占据轨道分割 so，更新了密度矩阵 D

  \[D_{\mu \nu} = 2 C_{\mu i} C_{\nu i}\]

• Line 18

  E_elec = (HC * D).sum() + 0.5 * np.einsum("uvkl, uv, kl ->", eri, D, D) - 0.25 * np.einsum("ukvl, uv, kl ->", eri, D, D)
  

  使用 SCF 收敛后的密度计算总能量 E_elec

  \[E_\mathrm{elec} [D_{\kappa \lambda}] = h_{\mu \nu} D_{\mu \nu} + \frac{1}{2} D_{\mu \nu} (\mu \nu | \kappa \lambda) D_{\kappa \lambda} - \frac{1}{4} D_{\mu \nu} (\mu \lambda| \kappa \nu) D_{\kappa \lambda}\]

至此，我们基本了解了 RHF 的实现过程。下面我们讨论一些更为细节的问题。

3.2.3. Hamiltonian Core 积分详述

我们刚才提到，在 \(h_{\mu \nu}\) 中，有动能的贡献量 \(t_{\mu \nu} = \langle \mu | \hat t | \nu \rangle\) 与核排斥能的贡献量 \(v_{\mu \nu} = \langle \mu | \hat v_\mathrm{nuc} | \nu \rangle\)。这两种积分可以通过更为底层的方式获得；特别是对于核排斥能的贡献量的理解，将会直接地影响到以后对 Hamiltonian Core 的原子核坐标梯度、二阶梯度的理解。

3.2.3.1. 动能积分

动能积分可以写为

\[t_{\mu \nu} = \langle \mu | \hat t | \nu \rangle = - \frac{1}{2} \phi_\mu \cdot (\partial_r^2 \phi_\nu) = - \frac{1}{2} \phi_\mu \phi_{r r \nu}\]

记号说明

• \(\phi\) 统一代表原子轨道函数，以电子坐标为自变量
• \(\phi_\mu\) 代表原子轨道 \(\mu\) 所对应的原子轨道函数
• \(\phi_{r \mu} = \partial_r \phi_\mu\) 代表原子轨道在电子坐标分量 \(r\) 下的偏导数
• \(\phi_{r w \mu} = \partial_r \partial_w \phi_\mu\) 代表原子轨道在电子坐标分量 \(r\) 与 \(w\) 下的二阶偏导数
• \(\boldsymbol{r}\) 作为加粗的 r 代表电子坐标；区别于电子坐标分量 \(r\) 是一维变量，\(\boldsymbol{r}\) 为三维向量

一般来说，如果一个表达式看起来是函数表达式，那么我们默认对其进行积分或者格点求和．譬如上式若不使用 Einstein Summation，则表达结果是是

\[t_{\mu \nu} = - \frac{1}{2} \int \phi_{\mu} (\boldsymbol{r}) \nabla_{\boldsymbol{r}}^2 \phi_{\nu} (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r}\]

在 PySCF 的积分引擎中，一个积分选项是生成关于 \((r, w, \mu, \nu)\) 的 AO 积分张量 int1e_ipipovlp \(\langle \partial_r \partial_w \mu | \nu \rangle = \phi_{r w \mu} \phi_\nu\)；我们可以对上述张量在 \(w = r\) 的情形求和，转置 \(\mu, \nu\)，并乘以系数 \(-0.5\)，就得到了动能积分 \(t_{\mu \nu}\) 了：

int1e_ipipovlp = mol.intor("int1e_ipipovlp").reshape((3, 3, nao, nao))
np.allclose(
    - 0.5 * (int1e_ipipovlp.diagonal(axis1=0, axis2=1).sum(axis=2)).T,
    mol.intor("int1e_kin")
)

True

任务 (4)

1. 上述代码单元中使用的是 sum(axis=2)，为什么？使用 sum(axis=0) 是否正确？使用 sum(axis=-1) 是否正确？

2. 我们还可以用另一种方法生成动能积分．现定义 int1e_ipovlpip 为 \(\langle \partial_r \mu | \partial_w \nu \rangle\)，请解释下述代码块为何输出 True？

   提示 1：算符 \(\partial_r\) 是，厄米算符还是反厄米算符？为什么？

   提示 2：动能算符为何是厄米算符？

   提示 3：上述代码中，如果不转置 \(\mu, \nu\)，即上述代码块的第三行末尾不加 .T，结果是否正确？

   对这些问题的了解将会允许我们更清晰地理解 AO 积分的对称性，辅助验证程序与公式的正确性，并辅助我们推导与核排斥势有关的导数．

int1e_ipovlpip = mol.intor("int1e_ipovlpip").reshape((3, 3, nao, nao))
np.allclose(
    0.5 * (int1e_ipovlpip.diagonal(axis1=0, axis2=1).sum(axis=-1)),
    mol.intor("int1e_kin"))

True

3.2.3.2. 核排斥势积分

势能积分可以写为

\[v_{\mu \nu} = \langle \mu | \frac{- Z_M}{|\boldsymbol{r} - \boldsymbol{M}|} | \nu \rangle = \langle \mu | \frac{- Z_M}{|\boldsymbol{r}|} | \nu \rangle_{\boldsymbol{r} \rightarrow M} = \left( \frac{- Z_M}{|\boldsymbol{r}|} \phi_\mu \phi_\nu \right)_{\boldsymbol{r} \rightarrow M}\]

记号说明

• 下标 \(\boldsymbol{r} \rightarrow M\) 代表电子积分的原点取在原子 \(M\) 的坐标上．

在 PySCF 的积分引擎中，\(\langle \mu | \frac{1}{|\boldsymbol{r}|} | \nu \rangle\) 的积分选项是 int1e_rinv；但其积分原点仍然是 \((0, 0, 0)\)。为了让特定原子坐标成为原点，PySCF 的一个便利函数是 gto.Mole.with_rinv_as_nucleus；它通过传入原子序号，将积分时的原点坐标更变为当前原子的坐标；但除了特定积分以外，分子的所有性质，包括坐标，都保持不变。

v = np.zeros((nao, nao))
for M in range(natm):
    with mol.with_rinv_as_nucleus(M):
        v += - mol.atom_charge(M) * mol.intor("int1e_rinv")
np.allclose(v, mol.intor("int1e_nuc"))

True

3.2.4. RHF 能量实现参考

在这里以及以后，实现参考一节将会展示不同的实现手段；这可能包括使用 PySCF 的高级函数，或者使用我们手写的 Python 脚本；并与当前的计算结果进行对照。

这份笔记的初衷有二：其一是记录非自洽 DFT 的计算方式；其二是尽可能只使用 PySCF 的积分、泛函与基组库，但不使用高级函数来构建我们的工作；即使使用高级函数，这些高级函数也已经通过其它实现参考用底层方法得以说明。

提示

一些不太容易编写，或者与效率有很强关联的程序，我们可能只叙述其原理，但最终还是会使用 PySCF 的库函数。SCF 过程和导出量、双电子积分函数、以及 pyscf.scf.cphf.solve 函数将会是其中几个例子。

3.2.4.1. 分子轨道系数 \(C_{\mu p}\)

PySCF 实现

PySCF 通过 mo_coeff 方法给出 \(C_{\mu p}\)：

np.allclose(np.abs(scf_eng.mo_coeff), np.abs(C), atol=1e-6, rtol=1e-4)

True

pyxdh 实现

pyxdh 可以通过 C 属性给出 \(C_{\mu p}\)：

np.allclose(np.abs(scfh.C), np.abs(C), atol=1e-6, rtol=1e-4)

True

任务 (5)

1. 如果对上述两个代码块去除 np.abs，即不对不同方法的 \(C_{\mu p}\) 取绝对值进行比较，是否还能给出 True 的结果？为什么？

2. 系数矩阵 \(C_{\mu p}\) 与重叠矩阵 \(S_{\mu \nu}\) 之间存在特殊的关系。请给出两者在某种矩阵乘法下给出单元矩阵的关系式。

   提示：可以利用 \(X_{\kappa \mu} S_{\kappa \lambda} X_{\lambda \nu} = \delta_{\mu \nu}\)。

pyxdh 还直接给出占据轨道、非占轨道下的分子轨道系数 Co、Cv；他们是根据占据轨道、非占轨道的分割 so、sv 给出的：

print(scfh.Co.shape)
print(scfh.Cv.shape)
print(np.allclose(scfh.Co, scfh.C[:, scfh.so]))
print(np.allclose(scfh.Cv, scfh.C[:, scfh.sv]))

(22, 9)
(22, 13)
True
True

3.2.4.2. 电子态密度 \(D_{\mu \nu}\)

\[D_{\mu \nu} = 2 C_{\mu i} C_{\nu i}\]

PySCF 实现

PySCF 通过 make_rdm1 成员函数给出：

np.allclose(scf_eng.make_rdm1(), D)

True

pyxdh 实现

pyxdh 可以通过 D 属性给出：

np.allclose(scfh.D, D)

True

3.2.4.3. Hamiltonian Core 积分 \(h_{\mu \nu}\)

\[h_{\mu \nu} = t_{\mu \nu} + v_{\mu \nu}\]

PySCF 实现

PySCF 可以通过 get_hcore 成员函数给出 \(h_{\mu \nu}\)：

scf_eng.get_hcore.__func__

<function pyscf.scf.hf.SCF.get_hcore(self, mol=None)>

np.allclose(scf_eng.get_hcore(), HC)

True

pyxdh 实现

pyxdh 可以通过 H_0_ao 属性给出 \(h_{\mu \nu}\)：

np.allclose(scfh.H_0_ao, HC)

True

同时，H_0_mo 方法给出 MO 基组下的 \(h_{pq}\)：

\[h_{pq} = C_{\mu p} h_{\mu \nu} C_{\nu q}\]

np.allclose(scfh.H_0_mo, scfh.C.T @ scfh.H_0_ao @ scfh.C)

True

3.2.4.4. 库伦积分 \(J_{\mu \nu}[X_{\kappa \lambda}]\)

\[J_{\mu \nu}[R_{\kappa \lambda}] = (\mu \nu | \kappa \lambda) R_{\kappa \lambda}\]

这里的 \(R_{\mu \nu}\) 代表的是任意矩阵；尽管在 SCF 过程中使用到库伦与交换积分处使用的是密度矩阵，但库伦积分的使用方式可以更为广泛，譬如使用广义密度而非电子态密度替代 \(R_{\kappa \lambda}\)。

PySCF 实现

PySCF 可以通过 get_j 成员函数实现库伦积分：

scf_eng.get_j.__func__

<function pyscf.scf.hf.SCF.get_j(self, mol=None, dm=None, hermi=1)>

R = np.random.random((nao, nao))
np.allclose(
    scf_eng.get_j(dm=R),
    np.einsum("uvkl, kl -> uv", mol.intor("int2e"), R)
)

True

3.2.4.5. 交换积分 \(K_{\mu \nu}[R_{\kappa \lambda}]\)

\[K_{\mu \nu}[R_{\kappa \lambda}] = (\mu \kappa | \nu \lambda) R_{\kappa \lambda}\]

PySCF 实现

PySCF 可以通过 get_k 成员函数实现交换积分：

注意

PySCF 中，交换积分对代入的 AO 基组广义密度矩阵有对称性要求。一般来说，我们以后工作中碰到的广义密度矩阵都是对称矩阵，因此 hermi 选项可以不设置。

scf_eng.get_k.__func__

<function pyscf.scf.hf.SCF.get_k(self, mol=None, dm=None, hermi=1)>

R = np.random.random((nao, nao))
[np.allclose(
    scf_eng.get_k(dm=R, hermi=hermi),
    np.einsum("ukvl, kl -> uv", mol.intor("int2e"), R)
) for hermi in [0, 1]]

[True, False]

R = np.random.random((nao, nao))
R += R.T
[np.allclose(
    scf_eng.get_k(dm=R, hermi=hermi),
    np.einsum("ukvl, kl -> uv", mol.intor("int2e"), R)
) for hermi in [0, 1]]

[True, True]

3.2.4.6. Fock 矩阵 \(F_{\mu \nu}[R_{\kappa \lambda}]\)

\[F_{\mu \nu}[R_{\kappa \lambda}] = h_{\mu \nu} + J_{\mu \nu}[R_{\kappa \lambda}] - \frac{1}{2} K_{\mu \nu}[R_{\kappa \lambda}]\]

PySCF 实现

PySCF 可以通过 get_fock 成员函数实现 Fock 矩阵；其可选参数之一是代入的矩阵 \(R_{\mu \nu}\)，默认下该矩阵为 \(D_{\mu \nu}\)：

注意

PySCF 中，Fock 矩阵同样对代入的广义密度矩阵有对称性要求。一般来说，我们也只处理对称矩阵。

scf_eng.get_fock.__func__

<function pyscf.scf.hf.get_fock(mf, h1e=None, s1e=None, vhf=None, dm=None, cycle=-1, diis=None, diis_start_cycle=None, level_shift_factor=None, damp_factor=None)>

R = np.random.random((nao, nao))
R += R.T
print(np.allclose(
    scf_eng.get_fock(dm=R),
    scf_eng.get_hcore() + scf_eng.get_j(dm=R) - 0.5 * scf_eng.get_k(dm=R)
))
print(np.allclose(scf_eng.get_fock(), F))

True
True

pyxdh 实现

pyxdh 可以通过 F_0_ao 属性给出 Fock 矩阵，但这个 Fock 矩阵是代入的是电子态密度的确定的矩阵 \(F_{\mu \nu} = F_{\mu \nu}[D_{\kappa \lambda}]\)：

np.allclose(
    scfh.F_0_ao,
    scf_eng.get_fock()
)

True

同时，F_0_mo 属性给出代入电子态密度的 MO 基组的 \(F_{pq}\)：

\[F_{pq} = C_{\mu p} F_{\mu \nu} C_{\nu q}\]

np.allclose(scfh.F_0_mo, scfh.C.T @ scfh.F_0_ao @ scfh.C)

True

3.2.4.7. 原子核排斥能 \(E_\mathrm{nuc}\)

\[E_\mathrm{nuc} = \frac{1}{2} Z_{AB} r_{AB}^{-1}\]

PySCF 实现

PySCF 可以通过 energy_nuc 成员函数实现排斥能：

np.allclose(scf_eng.energy_nuc(), E_nuc)

True

3.2.4.8. 电子态能量 \(E_\mathrm{elec}[R_{\mu \nu}]\)

\[E_\mathrm{elec}[R_{\mu \nu}] = (h_{\mu \nu} + \frac{1}{2} J_{\mu \nu} [R_{\kappa \lambda}] - \frac{1}{4} K_{\mu \nu} [R_{\kappa \lambda}]) R_{\mu \nu}\]

PySCF 实现

PySCF 可以通过 energy_elec 一般有两个返回值，前者是体系总能量；而后者是双电子积分能量：

scf_eng.energy_elec.__func__

<function pyscf.scf.hf.energy_elec(mf, dm=None, h1e=None, vhf=None)>

R = np.random.random((nao, nao))
R += R.T
print(np.allclose(
    scf_eng.energy_elec(dm=R)[0],
    ((scf_eng.get_hcore() + 0.5 * scf_eng.get_j(dm=R) - 0.25 * scf_eng.get_k(dm=R)) * R).sum()
))
print(np.allclose(
    scf_eng.energy_elec(dm=R)[1],
    ((0.5 * scf_eng.get_j(dm=R) - 0.25 * scf_eng.get_k(dm=R)) * R).sum()
))

True
True

如果代入的 \(R_{\mu \nu}\) 是电子态密度 \(D_{\mu \nu}\)，那么返回的能量将会是方才计算的双氧水电子态能量：

np.allclose(scf_eng.energy_elec()[0], E_elec)

True

3.2.4.9. 体系总能量 \(E_\mathrm{tot}[R_{\mu \nu}]\)

\[E_\mathrm{tot}[R_{\mu \nu}] = E_\mathrm{elec}[R_{\mu \nu}] + E_\mathrm{nuc}[R_{\mu \nu}]\]

PySCF 实现

PySCF 可以通过 energy_tot 成员函数实现体系总能量：

scf_eng.energy_tot.__func__

<function pyscf.scf.hf.energy_tot(mf, dm=None, h1e=None, vhf=None)>

R = np.random.random((nao, nao))
R += R.T
print(np.allclose(
    scf_eng.energy_tot(dm=R),
    scf_eng.energy_elec(dm=R)[0] + scf_eng.energy_nuc()
))
print(np.allclose(scf_eng.energy_tot(), E_tot))

True
True

pyxdh 实现

pyxdh 实例 scfh 的属性 eng 可以给出体系总能量，但不能代入任意 \(R_{\mu \nu}\) 所给出的能量：

np.allclose(scfh.eng, E_tot)

True

注意

PySCF 中，Fock 矩阵同样对代入的广义密度矩阵有对称性要求。一般来说，我们也只处理对称矩阵。

3.2.4.10. 轨道能量 \(\varepsilon_p\)

PySCF 实现

PySCF 可以通过 mo_energy 方法给出轨道能：

np.allclose(scf_eng.mo_energy, e)

True

pyxdh 实现

pyxdh 可以通过 e 属性给出：

np.allclose(scfh.e, e)

True

类似于系数矩阵 \(C_{\mu p}\)，pyxdh 也给出占据与非占的轨道能 eo、ev：

print(scfh.eo.shape)
print(scfh.ev.shape)
print(np.allclose(scfh.eo, scfh.e[scfh.so]))
print(np.allclose(scfh.ev, scfh.e[scfh.sv]))

(9,)
(13,)
True
True

事实上，轨道能就是 MO 基组下 Fock 矩阵 (作为对角矩阵) 的对角元：

np.allclose(np.diag(e), scfh.F_0_mo)

True

3.2.4.11. 轨道占据数

PySCF 实现

轨道占据数是指分子轨道的电子占据数量；对于 RHF 而言，占据轨道的每根轨道占据数为 2，而非占轨道的每根轨道占据数为 0。

scf_eng.mo_occ

array([2., 2., 2., 2., 2., 2., 2., 2., 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

3.2.5. 任务参考解答

3.2.5.1. 任务 (1)

3.2.5.1.1. 任务 (1.1)

• np 是 numpy，我们的整个张量计算都依靠 numpy 实现。
• scipy 用于计算稍复杂的矩阵与向量问题。尽管大多数时候我们不会使用 scipy；但这里我们需要使用矩阵的分数幂次以获得 \((\mathbf{S}^{-1/2})_{\mu \nu}\)。
• gto 是 PySCF 中用于给出基组、电子积分等信息的程序包，它们服务于各种量化计算。
• scf 是 PySCF 中计算自洽场的程序包；RHF 是其中一种自洽场，但 scf 还支持 UHF、ROHF、X2C 等自洽场计算。它是后续的各种性质和后自洽场量化方法的基础程序。
• resource_filename 是 Python 中用于给出程序库中的资源文件路径的函数；这是一种相对 Pythonic 的做法，可以避免使用有歧义的相对路径，以及避免跨系统平台可能会产生的问题。
• FormchkInterface 是 pyxdh 的辅助类，用于读取 Gaussian 程序所输出的 formchk 文件。在验证结果与单元测试中，该类与 resource_filename 会经常使用。
• Mol_H2O2 是 pyxdh 的辅助类，用于生成 H2O2 分子和一些预定义的自洽场方法。以后我们的学习都将基于这个分子。
• partial 可以用于重新定义函数的参数表。我们在 np.einsum 的使用：中转内存设置 一节中对此作过说明。
• np.set_printoptions 是控制 numpy 输出格式的语句；表明一般情况下，浮点数使用普通小数 (而非科学计数) 表示，且显示 5 位小数，每行显示 150 个字符。

3.2.5.2. 任务 (3)

3.2.5.2.1. 任务 (3.1)

np.allclose(X.T @ S @ X, np.eye(nao))

True

但请不要写为下述代码，尽管仍然会给出 True 的结果：

np.allclose(X @ S @ X, np.eye(nao))

True

之所以这么说，是因为 \(X_{\mu \nu}\) 的定义可以不唯一；使用 scipy.linalg.fractional_matrix_power 给出的 \(\mathbf{X} = \mathbf{S}^{-1/2}\) 是对称的。但如果使用 np.linalg.cholesky 与 np.linalg.inv 给出 \(\mathbf{X}\)，即是通过 Cholesky 分解 \(\mathbf{S} = \mathbf{L} \mathbf{L}^\dagger\) 给出的 \(\mathbf{X} = \mathbf{L^\dagger}^{-1}\)。

X_cd = np.linalg.inv(np.linalg.cholesky(S).T)
np.allclose(X_cd.T @ S @ X_cd, np.eye(nao))

True

但这时如果用错误的代码，就会给出 False 的结果：

np.allclose(X_cd @ S @ X_cd, np.eye(nao))

False

np.allclose(X_cd @ S @ X_cd.T, np.eye(nao))

False

读者可以尝试使用 Cholesky 分解给出的 X_cd \(X_{\mu \nu}\)，在这个双氧水分子的例子中，使用这种 Cholesky 分解比矩阵幂次少去将近 \(3/4\) 的迭代次数，看起来更为高效。

3.2.5.3. 任务 (4)

3.2.5.3.1. 任务 (4.1)

首先，我们回顾一下 int1e_ipipovlp 的维度：

int1e_ipipovlp.shape

(3, 3, 22, 22)

上述张量是 \(\phi_{r w \mu} \phi_\nu\)，其四个维度分别是 \(r, w, \mu, \nu\)。

随后，我们对 int1e_ipipovlp 取对角元，即令 \(w = r\)，其结果是 \(\phi_{r r \mu} \phi_\nu\)，其维度信息是

int1e_ipipovlp.diagonal(axis1=0, axis2=1).shape

(22, 22, 3)

注意这三个维度是 \(\mu, \nu, r\)；常识下我们可能会认为维度应该是 \(r, \mu, \nu\)，但之所以 \(r\) 会在最后一个维度，是因为 numpy 默认下会将被对角化的维度向最后放置。因此，动能积分的被求和的维度应当是第 3 维度，即 axis=2。若是其他维度，则被求和的角标将会是原子轨道角标。

在 numpy 中，axis=-1 代表最后一个维度，因此 axis=-1 也是正确的；axis=-1 有时是比 axis=2 更为推荐的用法，这可能在维度不确定的张量处理时给出普遍的结果。

3.2.5.3.2. 任务 (4.2)

提示 1 答 \(\partial_r\) 是反厄米算符。我们可以从 \(\langle \mu | \partial_r | \nu \rangle = \phi_{\mu} \phi_{r \nu}\) 矩阵的反对称性 \(\phi_{\mu} \phi_{r \nu} = - \phi_{\nu} \phi_{r \mu}\) 来验证：

int1e_ipovlp = mol.intor("int1e_ipovlp")
np.allclose(int1e_ipovlp, -int1e_ipovlp.swapaxes(-1, -2))

True

同时指出，\(\langle \mu | \partial_r | \nu \rangle\) 并非是零张量：

np.linalg.norm(int1e_ipovlp)

7.062131136920612

我们知道，动量算符是 \(- i \partial_r\) 是厄米算符，这与 \(\partial_r\) 为反厄米算符也是等价的。

之所以动量算符是厄米算符，或者等价地 \(\partial_r\) 为反厄米算符，可以参考 Levine [Lev13] (sec 7.2, p.158) 的说明。或者，可以参考下述简单的说明。根据分部积分原则，

\[\int \phi_\mu \partial_r \phi_\nu \, \mathrm{d} r = (\phi_\mu \phi_\nu) |_{r \rightarrow -\inf}^{+\inf} - \int \phi_\nu \partial_r \phi_\mu \, \mathrm{d} r\]

根据波函数的定义，\(\phi_\mu, \phi_\nu\) 应当满足 $r \rightarrow `:nbsphinx-math:infty $ 时，值趋近于零的性质；因此 :math:`(phi_mu phi_nu) |_{r rightarrow -inf}^{+inf} = 0；因此，

\[\int \phi_\mu \partial_r \phi_\nu \, \mathrm{d} r = - \int \phi_\nu \partial_r \phi_\mu \, \mathrm{d} r\]

即证明了 \(\partial_r\) 的反厄米性质。

提示 2 答 显然，我们知道动能算符 \(- \frac{1}{2} \partial_r^2\) 是厄米算符。

int1e_kin = mol.intor("int1e_kin")
np.allclose(int1e_kin, int1e_kin.swapaxes(-1, -2))

True

这是因为两个反厄米算符的乘积是厄米算符。普遍地，对于两个反厄米算符 \(\hat A, \hat B\)，

\[\begin{split}\begin{align} \langle \mu | \hat A \hat B | \nu \rangle &= \langle \mu | \hat A | \hat B \nu \rangle = - \langle \hat B \nu | \hat A | \mu \rangle^* \\ &= - \langle \hat B \nu | \hat A \mu \rangle^* = - \langle \hat A \mu | \hat B \nu \rangle \\ &= - \langle \hat A \mu | \hat B | \nu \rangle = \langle \nu | \hat B | \hat A \mu \rangle^* \\ &= \langle \nu | \hat B \hat A | \mu \rangle^* \end{align}\end{split}\]

对于 \(\hat B = \hat A\)，上式便能给出 \(\hat A^2\) 为厄米算符的结论；若问题只涉及实数，那么 \(\langle \mu | \hat A^2 | \nu \rangle = \langle \nu | \hat A^2 | \mu \rangle\)，因此动能矩阵是对称矩阵。

提示 3 答 结果正确；这可以根据动能算符是厄米算符的性质便立即给出。

int1e_ipovlpip = mol.intor("int1e_ipovlpip").reshape((3, 3, nao, nao))
np.allclose(
    0.5 * (int1e_ipovlpip.diagonal(axis1=0, axis2=1).sum(axis=-1)),
    mol.intor("int1e_kin"))

True

原问题 答

\[\begin{split}\begin{align} \langle \partial_r \mu | \partial_w \nu \rangle &= \langle \partial_r \mu | \partial_w | \nu \rangle = - \langle \nu | \partial_w | \partial_r \mu \rangle \\ &= - \langle \nu | \partial_w \partial_r | \mu \rangle \end{align}\end{split}\]

这里利用到 \(\partial_w\) 为反厄米算符以及问题不涉及复数的性质。这里有两种继续推导的方法：其一是令 \(r = w\)，因此

\[\langle \partial_r \mu | \partial_r \nu \rangle = - \langle \nu | \partial_r^2 | \mu \rangle = - \langle \mu | \partial_r^2 | \nu \rangle\]

另一种证明方法是利用对易性质 \([\partial_r, \partial_w]\)，因此

\[\langle \partial_r \mu | \partial_w \nu \rangle = - \langle \nu | \partial_w \partial_r | \mu \rangle = - \langle \nu | \partial_r \partial_w | \mu \rangle\]

利用 \(\langle \mu | \hat A \hat B | \nu \rangle = \langle \nu | \hat B \hat A | \mu \rangle\) 的结论，我们可以立即得到

\[\langle \partial_r \mu | \partial_w \nu \rangle = - \langle \mu | \partial_r \partial_w | \nu \rangle = - \langle \mu | \partial_w \partial_r | \nu \rangle\]

因此，\(\partial_r \partial_w\) 尽管两个算符并不一定相同，但确实是厄米算符。因此，\(\langle \partial_r \mu | \partial_r \nu \rangle\) 是关于 \(\mu, \nu\) 对称的。

最后，代入 \(\hat T = - \frac{1}{2} \partial_r^2\)，就可以得到以 \(\langle \partial_r \mu | \partial_w \nu \rangle\) 构造的动能矩阵 \(\langle \mu | \hat T | \nu \rangle\)，即原题的代码了。

3.2.5.4. 任务 (5)

3.2.5.4.1. 任务 (5.1)

分子轨道系数在列上 (分子的其中一根轨道系数) 可以取相反数，因此上述代码单元中使用了绝对值进行比较。取相反数一般不会影响到各种 AO 与 MO 基组矩阵的值。

3.2.5.4.2. 任务 (5.2)

关系式是

\[C_{\mu p} S_{\mu \nu} C_{\nu q} = \delta_{pq}\]

np.allclose(C.T @ S @ C, np.eye(nmo))

True

首先，我们指出，在 SCF 迭代过程中使用到的临时矩阵 Cp \(C'_{p \mu}\) (定义与性质参见 Szabo (3.175-178)) 因为是通过 \(\mathbf{F}'\) 对角化得到，因此具有以下正交矩阵的性质 (\({\mathbf{C}'}^\dagger \mathbf{C}' = \mathbf{1}\) 与 \({\mathbf{C}'}^{-1} = {\mathbf{C}'}^\dagger\))：

print(np.allclose(Cp.T @ Cp, np.eye(nmo)))
print(np.allclose(Cp.T, np.linalg.inv(Cp)))

True
True

因此，

\[\mathbf{C}^\dagger \mathbf{S} \mathbf{C} = (\mathbf{X} \mathbf{C}')^\dagger \mathbf{S} (\mathbf{X} \mathbf{C}') = {\mathbf{C}'}^\dagger (\mathbf{X}^\dagger \mathbf{S} \mathbf{X}) \mathbf{C}' = {\mathbf{C}'}^\dagger \mathbf{C}' = \mathbf{1}\]

将矩阵计算 \(\mathbf{C}^\dagger \mathbf{S} \mathbf{C} = \mathbf{1}\) 化为矩阵元素的运算后，便成为 \(C_{\mu p} S_{\mu \nu} C_{\nu q} = \delta_{pq}\)。

3.2.6. 参考文献

[Lev13]

Ira N. Levine. Quantum Chemistry (7th Edition). Pearson, 2013. ISBN 9780321803450. URL: https://www.amazon.com/Quantum-Chemistry-7th-Ira-Levine/dp/0321803450.

[SO96]

Attila Szabo and Neil S. Ostlund. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books on Chemistry). Dover Publications, 1996. ISBN 978-0486691862. URL: https://www.amazon.com/Modern-Quantum-Chemistry-Introduction-Electronic/dp/0486691861.