11. 单元课题:XYG3 能量计算

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

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

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

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

11.1. 程序流程导引

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

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

程序分为以下几个模块:

  1. 初始化
    • 引入库 (PySCF)
    • 定义分子 (PySCF)
    • 定义格点 (PySCF)
  2. 无需自洽场密度或轨道系数就能计算的变量 (自洽场无关)
    • 原子轨道积分 (PySCF)
    • 轨道格点与格点权重 (PySCF)
    • 原子核排斥能
  3. 需要代入自洽场密度或轨道系数的变量 (自洽场相关)
    • 库伦、交换积分
    • 密度格点与泛函格点 (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_Jgen_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_01rho_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

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