11. 单元课题：XYG3 能量计算

在结束这一单元前，我们通过完成一个比较完整的项目，回顾公式记号与最为基础的程序实现。这个比较完整的项目就是计算 XYG3 能量。

这个课题要求，除了 电子积分、轨道与泛函格点 外，尽可能只使用 不超过 numpy 的工具。

我们以后可能会使用一些程序上的技巧、以及 pyxdh 中所提供的一些便利工具来缩短文档和代码的篇幅；但作者认为，若要成为程序开发者，需要对一些必要的底层方法进行了解。这是作者编写这份课题的初衷。

作者认为，这份课题的所有代码未必需要亲手写一遍；这篇文档的代码前都会有导引，若看到导引就能知道代码大致是怎么写的 (调用哪个函数、或者能查阅到以前阅读过的哪篇文档的哪一小节、或者能正确地查阅到 PySCF 的 API 文档)，并且能不通过程序验证、正确地说出每个变量的维度，我认为就达成了作者期望读者阅读文档的目的了。

11.1. 程序流程导引

在设计程序之前，我们要知道 XYG3 的能量是如何给出的。

首先，我们需要跑一次 B3LYP 得到其密度矩阵 \(D_{\mu \nu}\) 与轨道系数 \(C_{\mu p}\)；
其次，我们将密度矩阵代入 XYG3 的 GGA 分项进行计算，得到其能量的 GGA 部分；
最后，将 B3LYP 得到的轨道系数代入 PT2 计算，得到 XYG3 能量的 PT2 部分。

程序分为以下几个模块：

初始化
- 引入库 (PySCF)
- 定义分子 (PySCF)
- 定义格点 (PySCF)
无需自洽场密度或轨道系数就能计算的变量 (自洽场无关)
- 原子轨道积分 (PySCF)
- 轨道格点与格点权重 (PySCF)
- 原子核排斥能
需要代入自洽场密度或轨道系数的变量 (自洽场相关)
- 库伦、交换积分
- 密度格点与泛函格点 (PySCF)
- 交换相关势与 Fock 矩阵
- GGA 能量
- SCF 循环
- XYG3 的 GGA 分项能量
- XYG3 的 PT2 分项能量
- XYG3 总能量

我们假定内存空间总是足够的。上述标记 PySCF 的部分是指我们允许在这些代码中使用 PySCF 程序，其它部分一概不允许 (包括不允许使用 pyxdh)。依据这些提示，读者应当能大致构思出程序框架，并能在 3 天时间以内从头写一个 XYG3 能量计算程序。

我们下面给出参考程序。

11.2. 初始化部分

11.2.1. 引入库

在引入 Python 库时，我们要考虑到以下方面：

我们需要使用到 PySCF 中的分子定义与电子积分引擎 gto、DFT 计算与泛函格点引擎 dft
np.einsum 的优化选项 optimize 需要常开
numpy 的输出稍简洁一些，这里使用 5 位小数输出
为了减少输出，因此不输出 Python 的 warning 信息

[1]:

import numpy as np
import scipy
import warnings
from pyscf import gto, dft
from functools import partial

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

11.2.2. 定义分子

如以前一样，我们定义如下的双氧水分子，基组为 6-31G：

[2]:

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

[2]:

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

我们以后可能会非常经常地使用占据、非占轨道数量与分割 (slice)，以及原子数量：

natm 原子数
nao 原子轨道数，nmo 分子轨道数，一般来说在量化程序中，两者相等
nocc 占据轨道数，nvir 非占轨道数
so, sv, sa 分别代表占据、非占、全轨道的分割

[3]:

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

11.2.3. 定义格点

我们定义下述 (99, 590) 格点：

[4]:

grids = dft.Grids(mol)
grids.atom_grid = (99, 590)
grids.becke_scheme = dft.gen_grid.stratmann
grids.build()

[4]:

<pyscf.dft.gen_grid.Grids at 0x7ff62f314780>

我们也会经常使用格点数量 ngrids：

[5]:

ngrid = grids.weights.size

11.3. 自洽场无关部分

11.3.1. 原子轨道积分

我们定义下述经常使用的原子轨道积分：

T 动能积分 \(t_{\mu \nu} = \langle \mu | - \frac{1}{2} \nabla^2 | \nu \rangle\)
Vnuc 外势能积分 \(v^\mathrm{nuc}_{\mu \nu} = \langle \mu | - \frac{Z_M}{| \boldsymbol{r} |} | \nu \rangle_{\boldsymbol{r} \rightarrow \boldsymbol{M}}\)
H_0_ao Hamiltonian Core 积分 \(h_{\mu \nu} = t_{\mu \nu} + v^\mathrm{nuc}_{\mu \nu}\)
S_0_ao 重叠积分 \(S_{\mu \nu} = \langle \mu | \nu \rangle\)
eri0_ao 双电子排斥积分 (ERI) \((\mu \nu | \kappa \lambda)\)

[6]:

T = mol.intor("int1e_kin")
Vnuc = mol.intor("int1e_nuc")
H_0_ao = T + Vnuc
S_0_ao = mol.intor("int1e_ovlp")
eri0_ao = mol.intor("int2e")
X = np.linalg.inv(np.linalg.cholesky(S_0_ao).T)

11.3.2. 轨道格点与权重

格点积分过程中会经常使用 PySCF 的 NumInt，在此我们用 ni 来表示 NumInt 的一个实例：

[7]:

ni = dft.numint.NumInt()

我们定义权重格点 weight \(w\)；它仅用于与泛函格点乘积，在公式中不会出现：

[8]:

weight = grids.weights

在计算 DFT 能量的过程中，我们至多使用轨道对电子坐标的一阶梯度。我们将会生成下述轨道格点：

ao_0 \(\phi_{\mu}\)
ao_1 \(\phi_{r \mu}\)

[9]:

ao = np.zeros((4, ngrid, nao))  # 4 refers to (noderiv, x_deriv, y_deriv, z_deriv)
g_start = 0
for inner_ao, _, _, _ in ni.block_loop(mol, grids, nao, deriv=1, max_memory=2000):
    ao[:, g_start:g_start+inner_ao.shape[-2]] = inner_ao
    g_start += inner_ao.shape[-2]
ao_0 = ao[0]
ao_1 = ao[1:4]

11.3.3. 原子核排斥能

回顾到原子核排斥能 \(E_\mathrm{nuc}\) 的表达式为

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

其中，

\[\begin{split}\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}\end{split}\]

我们定义如下变量

Z_A \(Z_A\) 原子核电荷数
A_t \(A_t\) 原子坐标分量
r_AB \(r_{AB}\) 原子间距离矩阵 (对角元设定为正无穷)
E_nuc \(E_\mathrm{nuc}\) 原子核排斥能

[10]:

Z_A = mol.atom_charges()
A_t = mol.atom_coords()
r_AB = np.linalg.norm(A_t[:, None, :] - A_t[None, :, :], axis=-1)
r_AB += np.diag(np.ones(natm) * np.inf)
E_nuc = 0.5 * (Z_A[None, :] * Z_A[:, None] / r_AB).sum()
E_nuc

[10]:

37.884674408641274

11.4. 自洽场相关部分

11.4.1. 库伦、交换积分

回顾到对于任意对称的、原子轨道下的密度矩阵 R \(R_{\mu \nu}\)，库伦积分与交换积分可以表示为

\[\begin{split}\begin{align} J_{\mu \nu} [R_{\kappa \lambda}] &= (\mu \nu | \kappa \lambda) R_{\kappa \lambda} \\ K_{\mu \nu} [R_{\kappa \lambda}] &= (\mu \kappa | \nu \lambda) R_{\kappa \lambda} \end{align}\end{split}\]

尽管以前一般来说，在讨论库伦与交换积分、以及它们对 Fock 矩阵的贡献时，会使用自洽场给出的密度矩阵；但我们希望手写一个自洽场迭代过程，因此需要将库伦、交换积分写成如下定义的代入密度、输出矩阵的函数 gen_J、gen_K 的形式：

[11]:

def gen_J(R):
    return np.einsum("uvkl, kl -> uv", eri0_ao, R)

def gen_K(R):
    return np.einsum("ukvl, kl -> uv", eri0_ao, R)

尽管上述计算中，我们还用到了 ERI 积分 eri0_ao \((\mu \nu | \kappa \lambda)\)，但由于 ERI 积分是由分子决定而不依赖于密度或泛函形式，因此在文档中我们可以不将 ERI 积分作为传入参数；类似地还有原子轨道格点。

11.4.2. 密度与泛函格点

密度 \(\rho\) 与密度梯度格点 \(\rho_r\) 格点需要通过密度矩阵与轨道格点获得：

\[\begin{split}\begin{align} \rho &= \phi_\mu \phi_\nu D_{\mu \nu} \\ \rho_r &= 2 \phi_{r \mu} \phi_\nu D_{\mu \nu} \end{align}\end{split}\]

我们仍然用函数定义上述过程，函数名为 gen_rho_grid，输入密度矩阵，输出密度 rho_0 \(\rho\)、密度梯度 rho_1 \(\rho_r\) 与 rho_01；rho_01 是密度与密度梯度的合并张量，仅仅用于生成泛函格点。

[12]:

def gen_rho_grid(D):
    rho_0 = np.einsum("gu, gv, uv -> g", ao_0, ao_0, D)
    rho_1 = 2 * np.einsum("rgu, gv, uv -> rg", ao_1, ao_0, D)
    rho_01 = np.concatenate(([rho_0], rho_1), axis=0)
    return rho_0, rho_1, rho_01

泛函的格点需要通过 PySCF 的 numint.eval_xc 函数给出，它需要我们输入泛函名称 (这个例子是 "B3LYPg")、轨道格点 rho_01，依次输出杂化系数 cx \(c_\mathrm{x}\)、泛函核 exc \(f\)、两个泛函核一阶导数 fr \(f_\rho\) 与 fg \(f_\gamma\)。需要注意，

我们不追求代码优化，因此不论是否需要泛函核一阶导数，我们都计算之。因此，numint.eval_xc 的第三个参数 deriv 始终设为 1；
这里再强调，出于公式简化的目的，输出的 \(f\), \(f_\rho\), \(f_\gamma\) 是已经被乘以了格点加权过的值；因此，在设计生成泛函格点权重的函数时，还需要将看似无关的格点权重引入。

[13]:

def gen_kernel_grid(xc_code, rho_01):
    cx = ni.hybrid_coeff(xc_code)
    exc, (fr, fg, _, _), _, _ = ni.eval_xc(xc_code, rho_01, deriv=1)
    exc *= weight
    fr *= weight
    fg *= weight
    return cx, exc, fr, fg

11.4.3. 交换相关势与 Fock 矩阵

回顾交换相关势 \(v_{\mu \nu}^\mathrm{xc} [\rho]\) 的计算：

\[v_{\mu \nu}^\mathrm{xc} [\rho] = f_\rho \phi_\mu \phi_\nu + 2 f_\gamma \rho_r (\phi_{r \mu} \phi_{\nu} + \phi_{\mu} \phi_{r \nu})\]

为了计算交换相关势，我们需要利用到所有生成密度与泛函格点的工具；因此，在函数 gen_vxc 的函数中，我们需要代入泛函名称与密度矩阵：

[14]:

def gen_vxc(xc_code, D):
    rho_0, rho_1, rho_01 = gen_rho_grid(D)
    _, _, fr, fg = gen_kernel_grid(xc_code, rho_01)
    vxc = (
        + np.einsum("g, gu, gv -> uv", fr, ao_0, ao_0)
        + 2 * np.einsum("g, rg, rgu, gv -> uv", fg, rho_1, ao_1, ao_0)
        + 2 * np.einsum("g, rg, gu, rgv -> uv", fg, rho_1, ao_0 ,ao_1)
    )
    return vxc

我们又知道，Fock 矩阵可以通过如下方式构造 (若 \(\rho\) 对应的密度矩阵是 \(D_{\mu \nu}\))：

\[F_{\mu \nu} [D_{\kappa \lambda}] = h_{\mu \nu} + J_{\mu \nu} [D_{\kappa \lambda}] - \frac{c_\mathrm{x}}{2} K_{\mu \nu} [D_{\kappa \lambda}] + v_{\mu \nu}^\mathrm{xc} [\rho]\]

因此，Fock 矩阵的生成函数 gen_F_0_ao 可以表示为

[15]:

def gen_F_0_ao(xc_code, D):
    cx = ni.hybrid_coeff(xc_code)
    F_0_ao = H_0_ao + gen_J(D) - 0.5 * cx * gen_K(D) + gen_vxc(xc_code, D)
    return F_0_ao

11.4.4. GGA 能量

回顾 GGA 总能量表达式为 (利用到交换相关能 \(E_\mathrm{xc, GGA} = f \rho\))

\[E_\mathrm{GGA} [D_{\kappa \lambda}] = h_{\mu \nu} D_{\mu \nu} + \frac{1}{2} J_{\mu \nu} [D_{\kappa \lambda}] D_{\mu \nu} - \frac{c_\mathrm{x}}{4} K_{\mu \nu} [D_{\kappa \lambda}] D_{\mu \nu} + f \rho\]

我们定义交换相关能 \(E_\mathrm{GGA}\) 的计算过程为 gen_energy_elec：

[16]:

def gen_energy_elec(xc_code, D):
    cx = ni.hybrid_coeff(xc_code)
    rho_0, _, rho_01 = gen_rho_grid(D)
    cx, exc, _, _ = gen_kernel_grid(xc_code, rho_01)
    energy_elec = ((H_0_ao + 0.5 * gen_J(D) - 0.25 * cx * gen_K(D)) * D).sum()
    energy_elec += (exc * rho_0).sum()
    return energy_elec

11.4.5. SCF 循环

通过前面的准备，我们已经可以实现 B3LYP 的自洽场循环了。

我们在以前的笔记中定义过 \(X_{\mu \nu}\)；其满足 \(X_{\mu \kappa} S_{\mu \nu} X_{\nu \lambda} = \delta_{\kappa \lambda}\) 或 \(\mathbf{X}^\dagger \mathbf{S} \mathbf{X} = \mathbf{1}\)。但若使用 scipy.linalg.eigh (API 文档)，我们可以直接求解 \(\mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C} \mathbf{\varepsilon}\)，避免手动定义 \(X_{\mu \nu}\) 并简化代码。

下述的自洽场过程的密度初猜是零。如果我们按照 Szabo 书中的 SCF 循环过程，应当会遇到密度矩阵剧烈振荡从而无法收敛的情况。为此，我们引入下述代码

D = 0.3 * D + 0.7 * D_old

它意味着每次迭代只更新 30% 的密度，剩余的 70% 密度仍然保留上一次迭代的结果；或者说我们让密度的更新过程更为保守。依靠这个小技巧，我们可以成功地在 100 步循环以内收敛密度，并且不会让密度振荡。

通过 SCF 过程，我们可以获得 B3LYP 下的轨道系数 C \(C_{\mu p}\)、轨道能 e \(\varepsilon_p\)、密度矩阵 D \(D_{\mu \nu}\)。尽管我们也计算出了 B3LYP 收敛后的 Fock 矩阵与能量，但这些对 XYG3 能量的计算没有贡献。

[17]:

C = e = NotImplemented
D = np.zeros((nao, nao))
D_old = np.zeros((nao, nao)) + 1e-4
count = 0

while (not np.allclose(D, D_old, atol=1e-8, rtol=1e-6)):
    if count > 500:
        raise ValueError("SCF not converged!")
    count += 1
    D_old = D
    F_0_ao = gen_F_0_ao("B3LYPg", D)
    e, C = scipy.linalg.eigh(F_0_ao, S_0_ao)  # Solve FC = SCe
    D = 2 * C[:, so] @ C[:, so].T
    if count > 1:
        D = 0.3 * D + 0.7 * D_old             # For convergence

E_elec = gen_energy_elec("B3LYPg", D)
E_tot = E_elec + E_nuc

print("SCF Converged in          ", count, " loops")
print("Electronic energy (B3LYP) ", E_elec, " a.u.")
print("Total energy      (B3LYP) ", E_tot, " a.u.")

SCF Converged in           54  loops
Electronic energy (B3LYP)  -189.2622179191188  a.u.
Total energy      (B3LYP)  -151.37754351047752  a.u.

11.4.6. XYG3 的 GGA 分项能量

由此，我们可以通过代入 B3LYP 密度，计算 XYG3 泛函的 GGA 部分能量 E_xyg3_gga \(E_\mathrm{XYG3, GGA}\)：

[18]:

E_xyg3_gga = gen_energy_elec("0.8033*HF - 0.0140*LDA + 0.2107*B88, 0.6789*LYP", D)
E_xyg3_gga

[18]:

-188.94500780243624

11.4.7. XYG3 的 PT2 分项能量

PT2 的能量计算与 MP2 能量计算完全一致，除了需要乘以相关系数 \(c_\mathrm{c}\)；对于 XYG3 而言，cc \(c_\mathrm{c} = 0.3211\)：

\[E_\mathrm{XYG3, PT2} = T_{ij}^{ab} t_{ij}^{ab} D_{ij}^{ab}\]

其中，

\[\begin{split}\begin{align} (pq|rs) &= (\mu \nu | \kappa \lambda) C_{\mu p} C_{\nu q} C_{\kappa r} C_{\lambda s} \\ D_{ij}^{ab} &= \varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b \\ t_{ij}^{ab} &= (ia|jb) / D_{ij}^{ab} \\ T_{ij}^{ab} &= c_\mathrm{c} (2 t_{ij}^{ab} - t_{ij}^{ba}) \end{align}\end{split}\]

我们分别用 eri0_mo, D_iajb, t_iajb, T_iajb 表示 \((pq|rs)\), \(D_{ij}^{ab}\), \(t_{ij}^{ab}\), \(T_{ij}^{ab}\)。需要注意，以后的记号中，\(T_{ij}^{ab}\) 会表示经过 \(c_\mathrm{c}\) 相乘后的张量；这会在后面的文档中再强调。这么做的目的是让 PT2 的梯度推导表达式与 MP2 的表达式基本一致。

所有张量在程序中的角标顺序是 \(i, a, j, b\)；这可能与公式中让人直接联想到的顺序不太相同。

[19]:

cc = 0.3211
eri0_mo = np.einsum("uvkl, up, vq, kr, ls", eri0_ao, C, C, C, C)
D_iajb = e[so, None, None, None] - e[None, sv, None, None] + e[None, None, so, None] - e[None, None, None, sv]
t_iajb = eri0_mo[so, sv, so, sv] / D_iajb
T_iajb = cc * (2 * t_iajb - t_iajb.swapaxes(-1, -3))

由此，XYG3 的 PT2 分项能量 E_xyg3_pt2 \(E_\mathrm{XYG3, PT2}\) 则为

[20]:

E_xyg3_pt2 = (T_iajb * t_iajb * D_iajb).sum()
E_xyg3_pt2

[20]:

-0.13594842432740734

11.4.8. XYG3 总能量

至此，我们已经求得所有 XYG3 的能量分项。现在我们只要将下述三个分项：GGA、PT2、原子核排斥能的结果相加就可以了：

\[E_\mathrm{XYG3} = E_\mathrm{XYG3, GGA} + E_\mathrm{XYG3, PT2} + E_\mathrm{nuc}\]

[21]:

E_xyg3 = E_xyg3_gga + E_xyg3_pt2 + E_nuc
E_xyg3

[21]:

-151.19628181812237

关于这个能量是否正确，可以看看上一篇文档的内容。