2. 分子结构、基组、电子积分

对于真空下的分子，通过定态的、电子态的 Schrodinger 方程 \(\hat H \Psi = E \Psi\) 的波函数 \(\Psi\)，原则上我们可以获得所有我们所期望的分子性质。\(E\) 是能量或其它物理可观测量表征，因此 \(E\) 也是被求解量。

\(\Psi\) 与 \(E\) 的求解的第一步即是确定 \(\hat H\) 的形式。为了确定 \(\hat H\)，我们至少要做四件事：

理论上，\(\hat H\) 是以分子构型为变量的算符量。因此，分子结构是必要的成分。
从实践上，由于我们不可能写出 \(\Psi\) 的解析形式，因此需要使用近似的、截断的函数展开。这些函数称为基组 (basis)。
所有与算符 \(\hat H\) 相关的运算结果则是电子积分。
事实上，即使使用了截断的函数展开，使用 \(\hat H\) 的计算 \(\Psi\) 的代价仍然非常巨大。为此，我们需要使用近似的量化方法来处理问题。

其中，分子结构、电子积分都是能使用程序确切定义的量；基组一般来说使用约定俗成的定义即可 (譬如 6-31G、cc-pVTZ 等)；但量化方法现在仍然在积极地发展。我们假设前三个问题不再是困扰我们的问题；而量化方法则是我们最为关心的问题。

尽管分子结构、基组与电子积分并不是这一系列文档的核心问题，但为了解决量化方法问题，我们有必要对这些问题有所熟悉。在这一小节中，这里我们讨论 PySCF 下的分子结构、基组与电子积分的调用过程。这也可以看作是对 pyscf.gto.Mole 类的不完全的解释文档。这将是我们第一份讨论量化程序问题的文档。

[52]:

import numpy as np
from pyscf import gto, lib
from functools import partial

np.set_printoptions(5, suppress=True, linewidth=120)
np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])

2.1. 分子结构

2.1.1. 分子构建

在这一节、以及今后的文档中，我们会始终以下述结构的双氧水分子进行说明。

[3]:

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()

[3]:

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

上述指令中，

mol = gto.Mole() 是初始化一个 gto.Mole 类到变量 mol 中。该类保存了与分子和基组有关的绝大部分信息、以及电子积分的调用方式。
mol.atom = ... 通过较为简单的方式作分子的结构定义。默认情况下，长度的单位是 Angstrom；这也类同于 Gaussian 的定义方式。在 PySCF 中，还可以通过内坐标等方式定义分子。
mol.basis = "6-31G" 是定义分子统一使用 6-31G 基组。PySCF 支持对分子中不同原子给出不同的基组的解决方案；但在这系列文档中就从简考虑。
mol.verbose = 0 指定 PySCF 对这个分子计算过程中的输出控制。默认值为 3，但为了减少文档长度，我们尽量静默这部分输出。
mol.build() 对分子进行构建。尽管一般来说，PySCF 会在进行量化计算时也会对分子进行构建，但建议平时写脚本时手动补上这一句。

也有一些指令可能在实际实践中经常会使用到，但在整个文档中不作考虑：

mol.cart = False 指定是否使用三维笛卡尔坐标形式的基组。默认是 False，即使用球谐形式的基组。默认与 Gaussian 中定义 5D 7F 是相同的。需要留意，在 Gaussian 中，6-31G 系列基组在 Gaussian 中的默认是 6D 7F，意味着在 Gaussian 中 (以及历史上的文献)，6-31G 的 \(d\) 轨道与 \(f, g, h, i, ...\) 轨道是分开考虑的，且前者使用球谐基组，后者使用笛卡尔基组。
mol.spin 指定分子的自旋数；在整份文档中，由于我们只考虑闭壳层的情形，因此使用默认值 0。
mol.charge 指定分子的带电数；由于我们使用的是双氧水分子而非离子，因此使用默认值 0。

2.1.2. 构型量输出

我们知道，一般来说，分子体系的能量是由电子态能量 \(E_\mathrm{elec}\) 与原子核互斥能量 \(E_\mathrm{nuc}\) 的加和构成的。如果看 Gaussian 程序的输出，SCF Done 之后的能量便是如此计算出来的。电子态的能量需要通过量化方法获得，但原子核互斥能只需要通过简单的初中物理知识就能解决。下面，我们通过计算原子核互斥能来了解 gto.Mole 类的分子构型量的输出。

首先，原子核互斥能可以表示为

\[E_\mathrm{nuc} = \frac{1}{2} \frac{Z_A Z_B}{r_{AB}}\]

回顾 Einstein Summation，上面的等式右边应当要对 \(A, B\) 有所求和。其中，对角线特殊处理的原子间距离矩阵

\[\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, B\)：表示原子；对于双氧水，可以是两个氢、两个氧原子中的任意一个。
三维向量 \(\boldsymbol{A}, \boldsymbol{B}\)：表示原子三维笛卡尔坐标。以上述双氧水为例，第四个原子 (氢原子) 的坐标表示为 \((0, 0.7, 1.0)\)，向量单位为 Angstrom。
下角标 \(t, s, r, w\)：表示三维坐标分量，取值范围 \(\{ x, y, z \}\)。以如果第四个原子被记号 \(A\) 表示，那么 \((A_x, A_y, A_z) = (0, 0.7, 1.0)\)，向量单位为 Angstrom。上述公式并没有出现，但以后会非常经常地使用到。
电荷标量 \(Z_A\)：\(A\) 原子的核电荷数；正值，单位为 a.u.。若第二个原子 (氧原子) 被记号 \(A\) 表示，那么 \(Z_A = 8\)。

记号冲突

在以后，上下角标 \(A, B\) 还可能代表被求导的分子性质量。在不加说明的情况下，这种记号定义会产生混淆。以后的文档中，默认情况下会让 \(A, B\) 代表原子；代表分子性质量时会在文档开始处作说明。
\(r\) 在斜体的情况下，如果单独出现则表示坐标分量；但若是 \(r_t\) 的带有下标的形式，则表示电子坐标 \(\boldsymbol{r}\) 的三维分量 \((r_x, r_y, r_z)\) 之一。
\(w\) 在斜体的情况下，如果单独出现则表示坐标分量；但若是 \(w_g\) 的带有下标的形式，则表示 DFT 格点在 \(g\) 点处的权重。

为此，我们需要得到分子的核电荷量 \(Z_A\) 与原子核坐标量 \(\boldsymbol{A}\)。核电荷量 Z_A \(Z_A\) 可以通过下述方法获得：

[4]:

Z_A = mol.atom_charges()
Z_A

[4]:

array([8, 8, 1, 1], dtype=int32)

可以看到，这是一个关于 \(A\) 的四维整数向量，分别代表两个氧、两个氢原子的核电荷数。

原子核坐标量 A_t \(\boldsymbol{A}\) 事实上应当是关于原子 \(A\) 与其坐标分量 \(t\) 的二维张量。一般来说，以后会写作 \(A_t\)：

[5]:

A_t = mol.atom_coords()
A_t

[5]:

array([[0.        , 0.        , 0.        ],
       [0.        , 0.        , 2.83458919],
       [1.88972612, 0.        , 0.        ],
       [0.        , 1.32280829, 1.88972612]])

需要注意到，上述二维张量的单位是 Bohr 半径，或者等价地，a.u.。如果我们将上述矩阵乘以 Bohr 半径到 Angstrom 的转换系数 lib.param.BOHR，那么就与我们构造分子时所使用的坐标值相同了：

[6]:

A_t * lib.param.BOHR

[6]:

array([[0. , 0. , 0. ],
       [0. , 0. , 1.5],
       [1. , 0. , 0. ],
       [0. , 0.7, 1. ]])

原子间距离矩阵 r_AB \(r_{AB}\) 可以通过下述循环获得，其中 mol.natm 表示分子的原子数目，对于双氧水而言则是 4：

[7]:

r_AB = np.empty((mol.natm, mol.natm))
for A in range(mol.natm):
    for B in range(mol.natm):
        if A != B:
            r_AB[A, B] = np.linalg.norm(A_t[A] - A_t[B])
        else:
            r_AB[A, B] = np.infty
r_AB

[7]:

array([[       inf, 2.83458919, 1.88972612, 2.3067047 ],
       [2.83458919,        inf, 3.40675222, 1.62560388],
       [1.88972612, 3.40675222,        inf, 2.98193753],
       [2.3067047 , 1.62560388, 2.98193753,        inf]])

我们已经凑齐了所有计算 E_nuc \(E_\mathrm{nuc}\) 所需要的所有要素了，随后便是简单的整合：

[8]:

E_nuc = 0.5 * np.einsum("A, B, AB ->", Z_A, Z_A, 1 / r_AB)
E_nuc

[8]:

37.88467440864127

在 PySCF 中，存在函数 mol.energy_nuc 用于计算原子核排斥能 \(E_\mathrm{nuc}\)。我们以该函数来验证我们的结果是否正确：

[9]:

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

[9]:

True

任务 (1)

(可选) 我们刚才使用了双重循环获得 r_AB \(r_{AB}\)。但 Python 的循环效率低下，一般来说不建议使用显式的 for 循环。请尝试不使用 for，构造 r_AB；并分别对使用与不使用 for 循环的效率作评价。

需要指出，计算原子核排斥能不是非常耗时的计算，因此使用显示 for 循环也未尝不可；但这权当是练习。

2.2. 基组与电子积分

2.2.1. 基组的表示

尽管基组已经通过 mol.basis = "6-31G" 有所定义，但显然我们还需要进一步将其转化为数字，才能让程序进行计算。获得基组的数值信息可以通过弱保护变量 mol._basis 获得：

[10]:

mol._basis

[10]:

{'O': [[0,
   [5484.6717, 0.0018311],
   [825.23495, 0.0139501],
   [188.04696, 0.0684451],
   [52.9645, 0.2327143],
   [16.89757, 0.470193],
   [5.7996353, 0.3585209]],
  [0, [15.539616, -0.1107775], [3.5999336, -0.1480263], [1.0137618, 1.130767]],
  [0, [0.2700058, 1.0]],
  [1, [15.539616, 0.0708743], [3.5999336, 0.3397528], [1.0137618, 0.7271586]],
  [1, [0.2700058, 1.0]]],
 'H': [[0,
   [18.731137, 0.0334946],
   [2.8253937, 0.23472695],
   [0.6401217, 0.81375733]],
  [0, [0.1612778, 1.0]]]}

我们简单回顾线性组合的 Gaussian 基组 (Contracted Gaussian-Type Orbital, CGTO) 的相关知识。由于我们以后不回顾这部分内容，因此记号都是临时的，并且与 Szabo 书中的记号较为相似。

对于高斯基组 CGTO，根据 (Szabo, 3.212)，可以表示为 (下述等式 RHS (Reft-Hand Side of equation) 的 Einstein Summation 是关于 \(p\) 的求和)

\[\phi_\mu^\mathrm{CGTO} (\boldsymbol{r} - \boldsymbol{A}) = d_{p \mu} \phi_p^\mathrm{GTO} (\alpha_{p \mu}, \boldsymbol{r} - \boldsymbol{A})\]

上述记号中，\(\mu\) 代表用于量化计算的原子轨道，\(p\) 代表原始的 Gaussian 函数 (primitive)，\(d_{p \mu}\) 代表 primitive 对 CGTO 的贡献系数；\(\alpha_{p \mu}\) 表征的是 Gaussian 函数的形状，值愈大则 Gaussian 函数愈尖锐。根据 (Szabo, 3.203)，如果原子轨道对应的是角量子数为 1 的原子轨道 (\(s\) 轨道) (下述公式不能用于高角量子数的函数)，那么 primitive 函数

\[\phi_{s}^\mathrm{GTO} (\alpha_{p \mu}, \boldsymbol{r} - \boldsymbol{A}) = \left( \frac{2 \alpha_{p \mu}}{\pi} \right)^{3/4} \tilde g_{s} (\alpha_{p \mu}, \boldsymbol{r} - \boldsymbol{A}) = \left( \frac{2 \alpha_{p \mu}}{\pi} \right)^{3/4} \exp(- \alpha_{p \mu} \Vert \boldsymbol{r} - \boldsymbol{A} \Vert_2^2)\]

记号说明

三维向量 \(\boldsymbol{r}\) 表示电子坐标。
下标 \(\mu, \nu, \kappa, \lambda\) 表示原子轨道角标。
函数或格点 \(\phi\) 一般都指原子轨道。以后原子轨道的右上角会去除 \(\mathrm{CGTO}\) 与 \(\mathrm{GTO}\) 的记号，因为一般来说不会引起歧义；但这一节会保留该记号。

记号冲突

\(r\) 在斜体的情况下，如果单独出现则表示坐标分量；但若是 \(r_t\) 的带有下标的形式，则表示电子坐标 \(\boldsymbol{r}\) 的三维分量 \((r_x, r_y, r_z)\) 之一。
\(p, q\) 在这一节中表示 primitive 角标；但以后会用来表示任意分子轨道角标，因为以后不再涉及 primitive 函数问题。

以氢原子为例，每个氢原子相当于贡献了两个原子轨道基组；这两个原子轨道可以表示为 (为了讨论方便，暂时假定这个氢原子的坐标是 \(\boldsymbol{A} = (0, 0, 0)\))

\begin{align} \phi_\mu^\mathrm{CGTO} (\boldsymbol{r}) = & + 0.0334946 \times \phi_s^\mathrm{GTO} (18.731137, \boldsymbol{r}) \\ & + 0.23472695 \times \phi_s^\mathrm{GTO} (2.8253937, \boldsymbol{r}) \\ & + 0.81375733 \times \phi_s^\mathrm{GTO} (0.6401217, \boldsymbol{r}) \\ \phi_\nu^\mathrm{CGTO} (\boldsymbol{r}) = & + 1.0 \times \phi_s^\mathrm{GTO} (0.1612778, \boldsymbol{r}) \end{align}

对应到方才输出的 mol._basis 的结果，我们指出，\(\phi_\mu^\mathrm{CGTO} (\boldsymbol{r})\) 的线性系数 \(d_{p \mu}\) 与指数系数 \(\alpha_{p \mu}\) 就是下述输出：

[11]:

mol._basis["H"][0]

[11]:

[0, [18.731137, 0.0334946], [2.8253937, 0.23472695], [0.6401217, 0.81375733]]

其中，列表的第一个值为 0，表示该轨道是 \(s\) 轨道；我们会在氧原子的基组下看到一些 primitive 列表的第一个值是 1，表示那是 \(p\) 轨道。而列表后三个值是三个子列表，分别表示三个 primitive 的指数系数 \(\alpha_{p \mu}\) 与线性系数 \(d_{p \mu}\)。

对于 \(\phi_\nu^\mathrm{CGTO}\)，我们也能从 mol._basis 中获取相应的信息：

[12]:

mol._basis["H"][1]

[12]:

[0, [0.1612778, 1.0]]

任务 (2)

我们假定读者已经阅读过 Szabo 书本的 1, 2, 3, 6 章，并对书中的 Appendix A 有所了解；这里的一些问题不是对这些知识的学习，而是从程序的角度上的复习。通过此练习，读者应当能对量化程序的最基石的电子积分的求取有大概的印象；在以后调用电子积分时，心里可以有些底气。

通过 mol._basis 的信息，指出双氧水分子在 6-31G 基组下的 primitive 函数数量、以及原子轨道数量。这两个值事实上可以分别通过 mol.npgto_nr() 与 mol.nao 分别给出；请通过这两个结果验证你的答案。
请设计函数
```
def integral_ovlp_gaussian_s(alpha, beta, rAB):
"""Implement your function here"""
```
该函数代入 alpha \(\alpha\), beta \(\beta\) 作为被积分的 Gaussian 函数的指数系数，rAB \(\Vert \boldsymbol{A} - \boldsymbol{B} \Vert_2\) 作为两原子的距离，得到重叠积分

\[\int \tilde g_{s} (\alpha, \boldsymbol{r} - \boldsymbol{A}) \tilde g_{s} (\alpha, \boldsymbol{r} - \boldsymbol{B}) \, \mathrm{d} \boldsymbol{r} = \left( \frac{\pi}{\alpha + \beta} \right)^{3 / 2} \exp \left( - \frac{\alpha \beta}{\alpha + \beta} \Vert \boldsymbol{A} - \boldsymbol{B} \Vert_2 \right)\]
请设计函数
```
def integral_ovlp_primitive_s(alpha, beta, rAB):
"""Implement your function here"""
```
该函数的输入参数与上一小题一致，但输出的结果是 primitive 基组的积分，或者对于现在的问题，等价于归一化后的 Gaussian 函数的积分：

\[\int \phi_{s}^\mathrm{GTO} (\alpha, \boldsymbol{r} - \boldsymbol{A}) \phi_{s}^\mathrm{GTO} (\alpha, \boldsymbol{r} - \boldsymbol{B}) \, \mathrm{d} \boldsymbol{r} = \left( \frac{2 \alpha_{p \mu}}{\pi} \right)^{3/4} \left( \frac{2 \beta_{p \mu}}{\pi} \right)^{3/4} \int \tilde g_{s} (\alpha, \boldsymbol{r} - \boldsymbol{A}) \tilde g_{s} (\alpha, \boldsymbol{r} - \boldsymbol{B}) \, \mathrm{d} \boldsymbol{r}\]
根据下述的氢原子的原子轨道的表述

\[\begin{split}\begin{align} \phi_\mu^\mathrm{CGTO} (\boldsymbol{r}) = & + 0.0334946 \times \phi_s^\mathrm{GTO} (18.731137, \boldsymbol{r}) \\ & + 0.23472695 \times \phi_s^\mathrm{GTO} (2.8253937, \boldsymbol{r}) \\ & + 0.81375733 \times \phi_s^\mathrm{GTO} (0.6401217, \boldsymbol{r}) \\ \phi_\nu^\mathrm{CGTO} (\boldsymbol{r}) = & + 1.0 \times \phi_s^\mathrm{GTO} (0.1612778, \boldsymbol{r}) \end{align}\end{split}\]

请尝试利用上述函数，计算下述积分：\((\mu | \nu)\)、\((\mu | \mu)\)、\((\nu | \nu)\)。其中，\(\mu, \nu\) 都属于相同的氢原子。你也可以尝试构建一个用于 \(s\) 的原子轨道的积分函数来解决该问题。对于 \((\mu | \mu)\)、\((\nu | \nu)\)，你可以使用原子轨道的归一性，推断其积分的值为 1。
如果现在 \((\mu | \nu)\)、\((\mu | \mu)\)、\((\nu | \nu)\) 竖线左右的原子轨道分属于上面定义的双氧水的两个氢原子上 (相距大约 1.578 Angstrom，但读者应当能用程序给出距离的精确值)；那么积分的值是多少？

上述题目的结果应当能通过后文会作说明的 mol.intor("int1e_ovlp")[18:, 18:] 的结果来验证。

注意上述的积分对非 \(s\) 轨道譬如 \(p, d, f, \cdots\) 轨道，是不适用的。因此，上述的结论不能推广到所有的氧原子的原子轨道。任意角量子数的电子积分的学习与实践将是非常考验代数、程序实践与调试、程序效率提升等能力，一般的量化方法开发者与发展者应当不需要对此有很深的认识。

任务 (3)

请前往 Basis Set Exchange，找到氢与氧原子的 6-31G 基组，并与 mol._basis 的输出进行比对。以后也许会遇到需要手动设置基组的情况；Basis Set Exchange 可以满足通常基组的需求。
简单回顾一下 6-31G 基组的意义。6-31G 基组可以看作 Double-zeta 基组；zeta 指 Slater 函数的指数系数。以氧原子为例，其 \(1s\) 轨道仍然是 Single-zeta 的，并且这个 Single-zeta 使用 6 个 primitive 组合而成；\(2s\) 与 \(2p\) 轨道是 Double-zeta 的，其中一个 zeta 使用 3 个 primitive 组合，而另 1 个 zeta 则用 1 个 primitive 表示。尽管 6-31G 被称为 Double-zeta，但有两处有些名不副实。一来，Double-zeta 只针对价层轨道，内层轨道仍然是 Single-zeta 的。二来，zeta 在 STO-3G 有明确的与 Slater 函数的对应关系，但 6-31G 基组的参数却是从拟合原子能量获得的，不是真正的对 Slater 函数的拟合；换言之不是真正的 Slater 函数的参数 \(\zeta\)。

现在问，对于氧原子的每个 6-31G 的 \(s\) 轨道基组 (一共有 3 个)，有几个原子轨道在全空间没有节面 (类似于 \(1s\) 轨道)，有几个原子轨道有一个球形节面 (类似于 \(2s\) 轨道)？

2.2.2. 电子积分

在上面的任务 (2) 中，我们通过现成的公式计算了氢原子之间重叠积分。在 PySCF 中，包括重叠积分的各种电子积分，在不考虑效率的情况下，可以使用 mol.intor 函数给出。

[13]:

S = mol.intor("int1e_ovlp")

该重叠积分的两个维度大小均是双氧水分子的原子轨道数量：

[14]:

S.shape

[14]:

(22, 22)

除了重叠积分，我们还可以获得各种各样的积分；譬如动能积分

[49]:

T = mol.intor("int1e_kin")
T.shape

[49]:

(22, 22)

不仅是矩阵类型的积分，张量的积分也是可以导出的，譬如双电子积分

[15]:

eri = mol.intor("int2e")

它是四维原子轨道数量的张量：

[16]:

eri.shape

[16]:

(22, 22, 22, 22)

一般来说，有上述积分后，若不考虑效率与收敛问题，一个粗略的 HF 自洽场计算就是可以实现的了。

但不是所有积分的维度都一定与原子轨道数量相关。对于负值的偶极积分，它是三维张量，第一维度代表笛卡尔坐标的三个分量：

[17]:

mdip = mol.intor("int1e_r")
mdip.shape

[17]:

(3, 22, 22)

所有种类的电子积分调用方式可以参考 PySCF 文档。显然在梯度求取过程中，不会仅仅使用上面叙述的积分；所需要的积分、包括这些积分的公式符号，会在以后的文档中作说明。

2.2.3. 壳层分割与原子轨道分割

获得各种电子积分后，尽管我们确实可以做各种量化计算，但这些积分的每个元素代表何种意义？这就需要使用 PySCF 中的分割函数 mol.aoslice_by_atom 来辅助我们理解。

[21]:

mol.aoslice_by_atom()

[21]:

array([[ 0,  5,  0,  9],
       [ 5, 10,  9, 18],
       [10, 12, 18, 20],
       [12, 14, 20, 22]])

上述列表的值作如下解释：

列表 4 行、4 列。其中，4 行是指原子数量为 4；而列数是固定的。
第一行 [0, 5, 0, 9] 是第一个氧原子的分割；我们知道，氧原子的 6-31G 基组包含 5 种不同的基函数，其中 3 个为 \(s\) 轨道，2 个为 \(p\) 轨道；因此，5 称为氧原子的壳层数 (Shell)。列表的前两个元素 0, 5 表示第一个氧原子壳层的起始与终止位置。
我们也知道，氧原子 6-31G 基组包含 9 个原子轨道，其中 3 个为 \(s\) 轨道，6 个分别为 2 个 \(p\) 壳层所派生出来的 \(p_x, p_y, p_z\) 轨道。列表后两个元素 0, 9 表示第一个氧原子的原子轨道的起始与终止位置。
第二行 [5, 10, 9, 18] 是第二个氧原子的分割。我们能发现其中的起始位置 5, 9 与第一个氧原子分割的终止位置相同。注意这里使用 Python 尾后指针或索引的习惯，因此第二个氧原子的壳层包含的数值是 5, 6, 7, 8, 9，不包含 10。对于原子轨道分割也作相同的理解。
最后一行 [12, 14, 20, 22] 中，我们能知道分子总壳层数是 14，而原子轨道总数是 22。

我们方才求取了两个氢原子轨道的重叠积分。根据上述列表，我们应当知道，这两个氢原子的的原子轨道序号分别是 18, 19 与 20, 21；因此任务 (2) 中有关重叠积分计算的答案应当在下述矩阵中：

[22]:

mol.intor("int1e_ovlp")[18:22, 18:22]

[22]:

array([[1.        , 0.65829205, 0.0410283 , 0.20486272],
       [0.65829205, 1.        , 0.20486272, 0.48819655],
       [0.0410283 , 0.20486272, 1.        , 0.65829205],
       [0.20486272, 0.48819655, 0.65829205, 1.        ]])

壳层的分割以后会用得比较少，但原子轨道的分割会经常使用到。我们以后会经常使用以下的函数来调取原子轨道分割：

[28]:

def mol_slice(atm, mol=mol):
    _, _, p0, p1 = mol.aoslice_by_atom()[atm]
    return slice(p0, p1)

那么，第二个氧原子的原子轨道与分子内其它原子轨道的重叠积分就可以通过下述代码实现：

[45]:

print(S[mol_slice(1), :].shape)
S[mol_slice(1), :]

(9, 22)

[45]:

array([[ 0.00000,  0.00033,  0.02046,  0.00000,  0.00000,  0.00106,  0.00000,  0.00000,  0.05835,  1.00000,  0.23369,
         0.16728,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.00025,  0.01847,  0.05133,  0.07625],
       [ 0.00033,  0.02130,  0.14113,  0.00000,  0.00000,  0.04048,  0.00000,  0.00000,  0.32436,  0.23369,  1.00000,
         0.76364,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.00920,  0.12048,  0.32292,  0.41501],
       [ 0.02046,  0.14113,  0.33799,  0.00000,  0.00000,  0.12857,  0.00000,  0.00000,  0.49783,  0.16728,  0.76364,
         1.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.00000,  0.08429,  0.29491,  0.48399,  0.72901],
       [ 0.00000,  0.00000,  0.00000,  0.00955,  0.00000,  0.00000,  0.08729,  0.00000,  0.00000,  0.00000,  0.00000,
         0.00000,  1.00000,  0.00000,  0.00000,  0.50152,  0.00000,  0.00000,  0.00923,  0.04778,  0.00000,  0.00000],
       [ 0.00000,  0.00000,  0.00000,  0.00000,  0.00955,  0.00000,  0.00000,  0.08729,  0.00000,  0.00000,  0.00000,
         0.00000,  0.00000,  1.00000,  0.00000,  0.00000,  0.50152,  0.00000,  0.00000,  0.00000,  0.28199,  0.11808],
       [-0.00106, -0.04048, -0.12857,  0.00000,  0.00000, -0.07073,  0.00000,  0.00000, -0.21713,  0.00000,  0.00000,
         0.00000,  0.00000,  0.00000,  1.00000,  0.00000,  0.00000,  0.50152, -0.01384, -0.07167, -0.20142, -0.08434],
       [ 0.00000,  0.00000,  0.00000,  0.08729,  0.00000,  0.00000,  0.33799,  0.00000,  0.00000,  0.00000,  0.00000,
         0.00000,  0.50152,  0.00000,  0.00000,  1.00000,  0.00000,  0.00000,  0.11887,  0.21658,  0.00000,  0.00000],
       [ 0.00000,  0.00000,  0.00000,  0.00000,  0.08729,  0.00000,  0.00000,  0.33799,  0.00000,  0.00000,  0.00000,
         0.00000,  0.00000,  0.50152,  0.00000,  0.00000,  1.00000,  0.00000,  0.00000,  0.00000,  0.48364,  0.37477],
       [-0.05835, -0.32436, -0.49783,  0.00000,  0.00000, -0.21713,  0.00000,  0.00000, -0.39527,  0.00000,  0.00000,
         0.00000,  0.00000,  0.00000,  0.50152,  0.00000,  0.00000,  1.00000, -0.17830, -0.32487, -0.34546, -0.26769]])

2.3. 任务参考解答

2.3.1. 任务 (1)

get_r_AB_withForLoop 是原文中描述的方法：

[12]:

def get_r_AB_withForLoop(mol):
    A_t = mol.atom_coords()
    r_AB = np.empty((mol.natm, mol.natm))
    for A in range(mol.natm):
        for B in range(mol.natm):
            if A != B:
                r_AB[A, B] = np.linalg.norm(A_t[A] - A_t[B])
            else:
                r_AB[A, B] = np.infty
    return r_AB

get_r_AB_withVectorization 是不使用显示 for，或称使用向量化的方法：

[13]:

def get_r_AB_withVectorization(mol):
    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)
    return r_AB

两者的结果是相同的：

[14]:

np.allclose(get_r_AB_withForLoop(mol), get_r_AB_withVectorization(mol))

[14]:

True

显然，不使用显式 for 的语句执行速度会快不少：

[15]:

%%timeit -r 10 -n 1000
get_r_AB_withForLoop(mol)

86.5 µs ± 2.83 µs per loop (mean ± std. dev. of 10 runs, 1000 loops each)

[16]:

%%timeit -r 10 -n 1000
get_r_AB_withVectorization(mol)

35.8 µs ± 4.03 µs per loop (mean ± std. dev. of 10 runs, 1000 loops each)

2.3.2. 任务 (2)

2.3.2.1. 任务 (2.1)

这个任务的本意实际上是希望读者能从 mol._basis 的输出，心算判断 primitive 与原子轨道数量。但如果希望有一种准确无误的计算方法，还是需要写一个小程序。这样一个小程序尽管不难，但需要耗费一些代码。读者可以向其中插入一些打印语句来获取一些被迭代对象是什么。

[17]:

# construct dict that stores primitive and ao number for atoms
bas_number_dict = {}
for atom, bas_list in mol._basis.items():
    primitive_number = 0
    ao_number = 0
    for ao in bas_list:
        ao_shape = ao[0]  # 0: s; 1: p; 2: d; 3: f ...
        # since we use spherical atomic orbital basis,
        # one defined basis contribute to atomic orbital numbers are
        # s: 1, p: 3, d: 5; f: 7, ...
        ao_number += 2 * ao_shape + 1
        for primitive in ao[1:]:
            primitive_number += 2 * ao_shape + 1
    bas_number_dict[atom] = {
        "primitive_number": primitive_number,
        "ao_number": ao_number
    }
bas_number_dict

[17]:

{'O': {'primitive_number': 22, 'ao_number': 9},
 'H': {'primitive_number': 4, 'ao_number': 2}}

[18]:

# calculate total primitive and ao number for H2O2 molecule
primitive_number = 0
ao_number = 0
for atm_id in range(mol.natm):  # mol.natm = 4 for H2O2
    atm_symbol = mol.atom_symbol(atm_id)  # return "O" or "H" for H2O2
    primitive_number += bas_number_dict[atm_symbol]["primitive_number"]
    ao_number += bas_number_dict[atm_symbol]["ao_number"]
print("H2O2 primitive number: ", primitive_number)
print("H2O2 ao        number: ", ao_number)

H2O2 primitive number:  52
H2O2 ao        number:  22

在 PySCF 中，他们可以通过下述函数或属性直接生成：

[19]:

print("H2O2 primitive number: ", mol.npgto_nr())
print("H2O2 ao        number: ", mol.nao)

H2O2 primitive number:  52
H2O2 ao        number:  22

2.3.2.2. 任务 (2.2)

[20]:

def integral_ovlp_gaussian_s(alpha, beta, rAB):
    return (np.pi / (alpha + beta))**(3/2) * np.exp(- alpha * beta / (alpha + beta) * rAB**2)

2.3.2.3. 任务 (2.3)

[21]:

def integral_ovlp_primitive_s(alpha, beta, rAB):
    return (2 * alpha / np.pi)**(3/4) * (2 * beta / np.pi)**(3/4) * integral_ovlp_gaussian_s(alpha, beta, rAB)

2.3.2.4. 任务 (2.4)

我们首先看正确答案应该是多少。如果用原子轨道的归一化来判断结果，那么 \((\mu | \mu)\) 与 \((\nu | \nu)\) 都应当为 1。另外，如果已经看完这一份笔记，那么应该知道 mol.intor("int1e_ovlp")[18:, 18:] 记录的是与这两个氢原子有关的重叠积分，并且应当能判断 mol.intor("int1e_ovlp")[18, 19] 即 0.65829205 是 \((\mu | \nu)\) 的值。

[22]:

mol.intor("int1e_ovlp")[18:, 18:]

[22]:

array([[1.        , 0.65829205, 0.0410283 , 0.20486272],
       [0.65829205, 1.        , 0.20486272, 0.48819655],
       [0.0410283 , 0.20486272, 1.        , 0.65829205],
       [0.20486272, 0.48819655, 0.65829205, 1.        ]])

下面我们写一个计算 CGTO 的原子轨道积分的函数：

[23]:

def integral_ovlp_ao_s(alpha_orb_list, beta_orb_list, rAB):
    integral = 0.
    for alpha, coef_a in alpha_orb_list:
        for beta, coef_b in beta_orb_list:
            integral += coef_a * coef_b * integral_ovlp_primitive_s(alpha, beta, rAB)
    return integral

其中，对于 \(\mu\) 轨道，alpha_orb_list 即是

[24]:

mol._basis["H"][0][1:]

[24]:

[[18.731137, 0.0334946], [2.8253937, 0.23472695], [0.6401217, 0.81375733]]

上面的列表的构造在前文叙述过：每一个子列表代表一个 primitive 基组；每一个 primitive 基组的第一个值是 \(\alpha\)，而第二个值是 CGTO 的系数 \(d_{p \mu}\)。由于这一小题中，我们指定了两个原子轨道属于同一原子，因此两原子轨道的原子中心距离 \(\Vert \boldsymbol{A} - \boldsymbol{B} \Vert_2 = 0\)。可能有一些精度上的误差。

[25]:

integral_ovlp_ao_s(mol._basis["H"][0][1:], mol._basis["H"][0][1:], 0)  # (\mu | \mu) = mol.intor("int1e_ovlp")[18, 18]

[25]:

1.0000000284768251

[26]:

integral_ovlp_ao_s(mol._basis["H"][0][1:], mol._basis["H"][1][1:], 0)  # (\mu | \nu) = mol.intor("int1e_ovlp")[18, 19]

[26]:

0.6582920587123634

[27]:

integral_ovlp_ao_s(mol._basis["H"][1][1:], mol._basis["H"][1][1:], 0)  # (\nu | \nu) = mol.intor("int1e_ovlp")[19, 19]

[27]:

0.9999999999999999

2.3.2.5. 任务 (2.5)

这一小题与上一小题的唯一区别是，现在被积分的两个原子轨道分属于两个不同的原子；我们只要指定这两个原子距离 rHH 即可。注意这个距离的单位是 a.u.，或 Bohr 半径，不是 Angstrom。

[28]:

rHH = np.linalg.norm(mol.atom_coords()[2] - mol.atom_coords()[3])

[29]:

integral_ovlp_ao_s(mol._basis["H"][0][1:], mol._basis["H"][0][1:], rHH)  # (\mu | \mu) = mol.intor("int1e_ovlp")[18, 20]

[29]:

0.04102829995317654

[30]:

integral_ovlp_ao_s(mol._basis["H"][0][1:], mol._basis["H"][1][1:], rHH)  # (\mu | \nu) = mol.intor("int1e_ovlp")[18, 21]

[30]:

0.20486272719145937

[31]:

integral_ovlp_ao_s(mol._basis["H"][1][1:], mol._basis["H"][1][1:], rHH)  # (\nu | \nu) = mol.intor("int1e_ovlp")[19, 21]

[31]:

0.488196553296388

2.3.3. 任务 (3)

2.3.3.1. 任务 (3.2)

所有三个基组都是没有空间节面的，类似于 \(1s\) 轨道。

对于氧原子，6-31G 的 6, 3, 1 分别对应下述的基组输出：

[32]:

print("`6` primitive number:", len(mol._basis["O"][0][1:]))
print(mol._basis["O"][0])
print("`3` primitive number:", len(mol._basis["O"][1][1:]))
print(mol._basis["O"][1])
print("`1` primitive number:", len(mol._basis["O"][2][1:]))
print(mol._basis["O"][2])

`6` primitive number: 6
[0, [5484.6717, 0.0018311], [825.23495, 0.0139501], [188.04696, 0.0684451], [52.9645, 0.2327143], [16.89757, 0.470193], [5.7996353, 0.3585209]]
`3` primitive number: 3
[0, [15.539616, -0.1107775], [3.5999336, -0.1480263], [1.0137618, 1.130767]]
`1` primitive number: 1
[0, [0.2700058, 1.0]]

经常地，我们会说，6-31G 基组的 6 模拟内层 \(1s\) 轨道，而 31 通过两个原子轨道模拟价层 \(2s\) 与 \(2p\) 轨道。此话不假，但用来模拟 \(2s\) 轨道的函数事实上没有空间的节面。我们实际上会指望用 631 三个基组的某种线性组合，达到产生有节面的函数的效果；而这种线性组合就交给量化方法通过变分法来给出。这也在 (Szabo, 3.301-302) 公式中有所反映。

事实上，任务 (3.2) 使用了误导性的问法。我们需要对“原子轨道”一词作说明。

用语解释

“原子轨道”(Atomic Orbital) 一词在通篇文档中表示的是 CGTO 基组。尽管这个概念来源于物理，但它现在是一个纯粹的量化计算上的用语。这与通常意义下的原子轨道 (氢原子本征函数) 并不相同。