2. RMP3 与 RMP4 能量#
创建时间：2019-11-01
这一节我们讨论 RMP3 与 RMP4 的能量计算。
读过 Szabo 第六章的人相信对 MPn 方法有所了解。不过我们这里不关注 MPn 方法的公式推导，并且将眼光局限于 Restricted 方法。事实上写这篇文档时，尽管曾经学习过一般的 MPn 与 CCPT 的推导方式，但并没有尝试推导过 Restricted 情况下的推导，仅仅是将书上出现的公式程序化而已。
这一篇文档的主要参考是 Helgaker et al. 2013 教材 [1] 的 section 14.4。关于 Spin-Orbital MP3 与 RMP3 的另一种实现方式，可以参考以 Szabo [2] 公式为蓝本的 Psi4NumPy 简要代码 [3]。
import numpy as np
from pyscf import scf, gto
from functools import partial
np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
import warnings
warnings.filterwarnings("ignore")
import re
np.set_printoptions(5, linewidth=150, suppress=True)
2.1. 分子体系与标准结果#
我们所使用的分子是非对称的双氧水分子，基组为 6-31G。比较常用但名称不太不常见的变量有
so,sv表示占据轨道、非占轨道的分割 (split)C,e分别表示分子轨道系数 \(C_{\mu p}\) 与轨道能 \(\varepsilon_p\)
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()
natm = mol.natm
nao = nmo = mol.nao
nocc = mol.nelec[0]
nvir = nmo - nocc
so, sv = slice(0, nocc), slice(nocc, nmo)
scf_eng = scf.RHF(mol)
scf_eng.conv_tol = 1e-12
scf_eng.max_cycle = 128
scf_eng.kernel()
-150.5850337808384
C, e = scf_eng.mo_coeff, scf_eng.mo_energy
Co, Cv = C[:, so], C[:, sv]
eo, ev = e[so], e[sv]
我们的参考值来自于 Gaussian 的计算 (输入卡 H2O2-MP4.gjf、输出文件 H2O2-MP4.out)。
with open("H2O2-MP4.out", "r") as f:
gaussian_output = f.readlines()
def gaussian_line(string):
for line in gaussian_output:
if string in line:
return line[:-1].replace("D+", "E+").replace("D-", "E-")
其中一些结论有：
ref_mp2_corr = float(gaussian_line("E2 =").split()[2])
ref_mp3_corr = float(gaussian_line("E3=").split()[1])
ref_mp4_S = float(gaussian_line("MP4(S)=").split()[1])
ref_mp4_D = float(gaussian_line("MP4(D)=").split()[1])
ref_mp4_Q = float(gaussian_line("MP4(R+Q)=").split()[1])
ref_mp4_T = float(gaussian_line("E4(SDTQ)=").split()[1]) - float(gaussian_line("E4(SDQ)=").split()[1])
print("MP2 Corr: {:16.10f}".format(ref_mp2_corr))
print("MP3 Corr: {:16.10f}".format(ref_mp3_corr))
print("MP4 S : {:16.10f}".format(ref_mp4_S))
print("MP4 D : {:16.10f}".format(ref_mp4_D))
print("MP4 T : {:16.10f}".format(ref_mp4_T))
print("MP4 Q : {:16.10f}".format(ref_mp4_Q))
MP2 Corr: -0.2690117593
MP3 Corr: 0.0074517392
MP4 S : -0.0043821604
MP4 D : -0.0088190269
MP4 T : -0.0067126030
MP4 Q : 0.0011355775
2.2. RMP2 相关能量#
这里基本参照了 Helgaker 书的叙述思路，因此可能与我曾经写过的文档记号和变量名有微妙的区别。
变量名 |
公式表达式 |
意义 |
程序的角标顺序 |
出处 |
|---|---|---|---|---|
|
\(\varepsilon_{ij}^{ab}\) |
轨道能差 |
\((i, a, j, b)\) |
eq (14.4.18) |
|
\(g_{pqrs}\) |
分子轨道 ERI 积分 |
\((p, q, r, s)\) |
- |
|
\({t_{ij}^{ab}}^{(1)}\) |
一阶激发系数 |
\((i, a, j, b)\) |
eq (14.4.41) |
|
\(L_{pqrs}\) |
定义量 |
\((p, q, r, s)\) |
eq (13.7.15) |
|
\({\bar t_{ij}^{ab}}^{(1)}\) |
定义量 |
\((i, a, j, b)\) |
eq (14.4.42) |
上述定义的一些实际表达式为
eri_ao = mol.intor("int2e")
eri_mo = np.einsum("uvkl, up, vq, kr, ls -> pqrs", eri_ao, C, C, C, C)
E_iajb = - e[so, None, None, None] + e[None, sv, None, None] - e[None, None, so, None] + e[None, None, None, sv]
g_mo = eri_mo
t_iajb = - g_mo[so, sv, so, sv] / E_iajb
L_mo = 2 * g_mo - g_mo.swapaxes(-1, -3)
T_iajb = 4 * t_iajb - 2 * t_iajb.swapaxes(-1, -3)
在这些定义下，我们可以相当容易地写出 RMP2 的能量 (Helgaker, eq (14.4.55))
energy_mp2_corr = (t_iajb * L_mo[so, sv, so, sv]).sum()
energy_mp2_corr
-0.26901177170795454
2.3. RMP3 相关能#
RMP3 能量可以根据 Helgaker, eq (14.4.61-62) 构建：
mp3_X_iajb = (
+ 0.5 * np.einsum("icjd, acbd -> iajb", t_iajb, g_mo[sv, sv, sv, sv])
+ 0.5 * np.einsum("kalb, kilj -> iajb", t_iajb, g_mo[so, so, so, so])
+ np.einsum("iakc, bjkc -> iajb", t_iajb, L_mo[sv, so, so, sv])
- np.einsum("kajc, bcki -> iajb", t_iajb, g_mo[sv, sv, so, so])
- np.einsum("kaic, bjkc -> iajb", t_iajb, g_mo[sv, so, so, sv])
)
关于 Helgaker 书上的 eq (14.4.59)，可能是由于我们的实现与之实现方式不同，因此 Helgaker 书中的 \({\tilde t_{ij}^{ab}}^{(1)}\) 可以当作 \({\bar t_{ij}^{ab}}^{(1)}\) 使用。因此，
energy_mp3_corr = (T_iajb * mp3_X_iajb).sum()
energy_mp3_corr
0.007451743819516383
2.4. RMP4(SDQ) 相关能#
在 RMP4 相关能计算过程中，由于其中涉及到的三激发项从实现上较为复杂，因此我们将会对其分开讨论。
2.4.1. RMP4(S) 相关能#
定义以下变量：
变量名 |
公式表达式 |
意义 |
程序的角标顺序 |
出处 |
|---|---|---|---|---|
|
\(\varepsilon_{i}^{a}\) |
轨道能差 |
\((i, a)\) |
eq (14.4.18) |
|
\({t_{i}^{a}}^{(2)}\) |
二阶激发系数 |
\((i, a)\) |
eq (14.4.50) |
\(\varepsilon_i^a\) 的定义是显然的；\({t_{i}^{a}}^{(2)}\) 的定义为
定义以下变量：
变量名 |
公式表达式 |
意义 |
程序的角标顺序 |
出处 |
|---|---|---|---|---|
|
\(\varepsilon_{i}^{a}\) |
轨道能差 |
\((i, a)\) |
eq (14.4.18) |
|
\({t_{i}^{a}}^{(2)}\) |
二阶激发系数 |
\((i, a)\) |
eq (14.4.50) |
E_ia = - e[so, None] + e[None, sv]
t_2_ia = (
+ np.einsum("kald, kild -> ia", t_iajb, L_mo[so, so, so, sv])
- np.einsum("kcid, adkc -> ia", t_iajb, L_mo[sv, sv, so, sv])
)
t_2_ia /= E_ia
由此，根据 Helgaker eq (14.4.79) 与 eq (14.4.83)，有
mp4_S_iajb = (
+ np.einsum("jc, aibc -> iajb", t_2_ia, g_mo[sv, so, sv, sv])
- np.einsum("kb, aikj -> iajb", t_2_ia, g_mo[sv, so, so, so])
)
以及
energy_mp4_S = (T_iajb * mp4_S_iajb).sum()
energy_mp4_S
-0.004382160888639593
2.4.2. RMP4(D) 相关能#
定义以下变量：
变量名 |
公式表达式 |
意义 |
程序的角标顺序 |
出处 |
|---|---|---|---|---|
|
\({t_{ij}^{ab}}^{(2)}\) |
二阶激发系数 |
\((i, a, j, b)\) |
eq (14.4.51) |
\({t_{ij}^{ab}}^{(2)}\) 的定义为
注意到这里的算符 \(\hat P_{ij}^{ab}\) 是对称化算符，其出处是 Helgaker eq (13.7.13)：
为此，我们将 \({t_{ij}^{ab}}^{(2)}\) 的后一个关于 \((i, a, j, b)\) 的张量 \(\left( {t_{ik}^{ac}}^{(1)} L_{bjkc} - {t_{kj}^{ac}}^{(1)} g_{bcki} - {t_{ki}^{ac}}^{(1)} g_{bjkc} \right)\) 先使用变量 tmp_symm 储存，随后使用 np.transpose 转置来执行对称化。
tmp_symm = (
+ np.einsum("iakc, bjkc -> iajb", t_iajb, L_mo[sv, so, so, sv])
- np.einsum("kajc, bcki -> iajb", t_iajb, g_mo[sv, sv, so, so])
- np.einsum("kaic, bjkc -> iajb", t_iajb, g_mo[sv, so, so, sv])
)
tmp_symm += tmp_symm.transpose((2, 3, 0, 1))
t_2_iajb = (
- np.einsum("icjd, acbd -> iajb", t_iajb, g_mo[sv, sv, sv, sv])
- np.einsum("kalb, kilj -> iajb", t_iajb, g_mo[so, so, so, so])
- tmp_symm
)
t_2_iajb /= E_iajb
由此，根据 Helgaker eq (14.4.80) 与 eq (14.4.83)，有
mp4_D_iajb = (
+ 0.5 * np.einsum("icjd, acbd -> iajb", t_2_iajb, g_mo[sv, sv, sv, sv])
+ 0.5 * np.einsum("kalb, kilj -> iajb", t_2_iajb, g_mo[so, so, so, so])
+ 1.0 * np.einsum("iakc, bjkc -> iajb", t_2_iajb, L_mo[sv, so, so, sv])
- 1.0 * np.einsum("kajc, bcki -> iajb", t_2_iajb, g_mo[sv, sv, so, so])
- 1.0 * np.einsum("kaic, bjkc -> iajb", t_2_iajb, g_mo[sv, so, so, sv])
)
以及
energy_mp4_D = (T_iajb * mp4_D_iajb).sum()
energy_mp4_D
-0.008819028049601586
2.4.3. RMP4(Q) 相关能#
在 Helgaker 的书中，RMP4(Q) 的相关能没有再作额外定义。根据 Helgaker eq (14.4.82) 与 eq (14.4.83)，有
mp4_Q_iajb = (
+ 0.5 * np.einsum("kalb, icjd, kcld -> iajb", t_iajb, t_iajb, g_mo[so, sv, so, sv])
+ 1.0 * np.einsum("iakc, jbld, kcld -> iajb", t_iajb, t_iajb, L_mo[so, sv, so, sv])
- 1.0 * np.einsum("iakc, lbjd, kcld -> iajb", t_iajb, t_iajb, L_mo[so, sv, so, sv])
+ 0.5 * np.einsum("kaic, lbjd, kcld -> iajb", t_iajb, t_iajb, g_mo[so, sv, so, sv])
+ 0.5 * np.einsum("kajd, lbic, kcld -> iajb", t_iajb, t_iajb, g_mo[so, sv, so, sv])
- 1.0 * np.einsum("iakb, lcjd, lckd -> iajb", t_iajb, t_iajb, L_mo[so, sv, so, sv])
- 1.0 * np.einsum("iajc, kbld, kcld -> iajb", t_iajb, t_iajb, L_mo[so, sv, so, sv])
)
以及
energy_mp4_Q = (T_iajb * mp4_Q_iajb).sum()
energy_mp4_Q
0.001135577514293328
由此，我们可以计算 RMP4(SDQ) 所贡献的相关能：
energy_mp4SDQ_corr = energy_mp4_S + energy_mp4_D + energy_mp4_Q
energy_mp4SDQ_corr
-0.01206561142394785
2.5. RMP4(T) 相关能#
RMP4(T) 从实现上来讲，最简单的做法需要消耗 \(O^3 V^3\) 的内存。我们先从简单的实现入手；随后会简要给出一个中间张量的内存消耗是 \(O^2 V^2\) 的算法。
2.5.1. RMP4(T) 相关能：简单实现#
定义以下变量：
变量名 |
公式表达式 |
意义 |
程序的角标顺序 |
出处 |
|---|---|---|---|---|
|
\(\varepsilon_{ijk}^{abc}\) |
轨道能差 |
\((i, a, j, b, k, c)\) |
eq (14.4.18) |
|
\({t_{ijk}^{abc}}^{(2)}\) |
二阶激发系数 |
\((i, a, j, b, k, c)\) |
eq (14.4.52) |
\(\varepsilon_{ijk}^{abc}\) 的定义是显然的；\({t_{ijk}^{abc}}^{(2)}\) 的定义为
这里的算符 \(\hat P_{ijk}^{abc}\) 仍然是对称算符，但其形式更复杂 (Helgaker, eq (13.7.14))：
E_iajbkc = (
- e[so, None, None, None, None, None]
+ e[None, sv, None, None, None, None]
- e[None, None, so, None, None, None]
+ e[None, None, None, sv, None, None]
- e[None, None, None, None, so, None]
+ e[None, None, None, None, None, sv]
)
t_2_iajbkc = (
+ np.einsum("iajd, ckbd -> iajbkc", t_iajb, g_mo[sv, so, sv, sv])
- np.einsum("ialb, cklj -> iajbkc", t_iajb, g_mo[sv, so, so, so])
)
t_2_iajbkc = (
+ t_2_iajbkc.transpose((0, 1, 2, 3, 4, 5))
+ t_2_iajbkc.transpose((0, 1, 4, 5, 2, 3))
+ t_2_iajbkc.transpose((2, 3, 0, 1, 4, 5))
+ t_2_iajbkc.transpose((2, 3, 4, 5, 0, 1))
+ t_2_iajbkc.transpose((4, 5, 0, 1, 2, 3))
+ t_2_iajbkc.transpose((4, 5, 2, 3, 0, 1))
)
t_2_iajbkc = - t_2_iajbkc / E_iajbkc
由此，根据 Helgaker eq (14.4.81) 与 eq (14.4.83)，有
mp4_T_iajb = (
+ np.einsum("iajckd, bckd -> iajb", t_2_iajbkc, L_mo[sv, sv, so, sv])
- np.einsum("kajcid, kdbc -> iajb", t_2_iajbkc, g_mo[so, sv, sv, sv])
- np.einsum("iakblc, kjlc -> iajb", t_2_iajbkc, L_mo[so, so, so, sv])
+ np.einsum("lakbic, kjlc -> iajb", t_2_iajbkc, g_mo[so, so, so, sv])
)
以及
energy_mp4_T = (T_iajb * mp4_T_iajb).sum()
energy_mp4_T
-0.006712603570763939
由此，我们可以给出完整的 RMP4 相关能矫正了：
energy_mp4_corr = energy_mp4_S + energy_mp4_D + energy_mp4_T + energy_mp4_Q
energy_mp4_corr
-0.018778214994711787
2.5.2. RMP4(T) 相关能：中间张量 \(O^2 V^2\) 内存实现#
这里我们通过限制 \({t_{ijk}^{abc}}^{(2)}\) 中的 \(b\) 与 \(k\) 并作求和的方式，将中间矩阵的储存大小降为 \(O^2 V^2\) 并给出 RMP4(T) 的能量。
这个严格来说还是至少消耗了 \(O^1 V^3\) 的内存，因为在计算过程中用到了张量 \(g_{aicd}\)。
首先指出，像上述代码中使用 np.transpose 的方式来处理 \(\hat P_{ijk}^{abc}\) 在这里可能不适用；因此需要手输所有的中间张量的缩并过程。下述函数可以用来辅助我们进行带有指标转换的张量缩并。
# https://stackoverflow.com/a/11122744/9647779
def substitute_string(string, rule):
str_lst = list(string)
rule1, rule2 = rule.replace(" ", "").split("->")
idx_list = [[i.start() for i in re.finditer(c, string)] for c in rule1]
for idx, pos_list in enumerate(idx_list):
for pos in pos_list:
str_lst[pos] = rule2[idx]
return "".join(str_lst)
譬如我们现在需要将缩并的目标 \({t_{ij}^{ad}}^{(1)} g_{ckbd}\) 的 \((i, a, j, b, k, c)\) 转换为 \((j, b, k, c, i, a)\)，那么我们执行下述代码就可以生成转换后的缩并字符串了：
substitute_string("iajd, ckbd", "iajbkc -> jbkcia")
'jbkd, aicd'
后文中出现的代码就是受这个小函数的帮助而写成的。
事实上，下述程序几乎与上面的程序等价，只是每次处理 \({t_{ijk}^{abc}}^{(2)}\) 张量时，总是先固定 \(b, k\) 指标，而对其它的张量进行正常的运算，求出固定了 \(b, k\) 的 RMP4(T) 相关能；随后我们再对 \(b, k\) 指标进行循环，把所有 \(b, k\) 相关能贡献总和起来，得到最终的 RMP4(T) 相关能。
energy_mp4_T_by_ckiter = 0
for k in range(nocc):
for b in range(nvir):
bn = b + nocc
t_2_bk_iajb = - (
+ np.einsum("iajd, cd -> iajc", t_iajb[:, :, :, :], g_mo[sv, k, bn, sv])
- np.einsum("ial, clj -> iajc", t_iajb[:, :, :, b], g_mo[sv, k, so, so])
+ np.einsum("iad, jcd -> iajc", t_iajb[:, :, k, :], g_mo[bn, so, sv, sv])
- np.einsum("ialc, jl -> iajc", t_iajb[:, :, :, :], g_mo[bn, so, so, k])
+ np.einsum("jid, cad -> iajc", t_iajb[:, b, :, :], g_mo[sv, k, sv, sv])
- np.einsum("jla, cli -> iajc", t_iajb[:, b, :, :], g_mo[sv, k, so, so])
+ np.einsum("jd, aicd -> iajc", t_iajb[:, b, k, :], g_mo[sv, so, sv, sv])
- np.einsum("jlc, ail -> iajc", t_iajb[:, b, :, :], g_mo[sv, so, so, k])
+ np.einsum("cid, jad -> iajc", t_iajb[k, :, :, :], g_mo[bn, so, sv, sv])
- np.einsum("cla, jli -> iajc", t_iajb[k, :, :, :], g_mo[bn, so, so, so])
+ np.einsum("cjd, aid -> iajc", t_iajb[k, :, :, :], g_mo[sv, so, bn, sv])
- np.einsum("cl, ailj -> iajc", t_iajb[k, :, :, b], g_mo[sv, so, so, so])
)
t_2_bk_iajb /= (
- e[so, None, None, None]
+ e[None, sv, None, None]
- e[None, None, so, None]
+ e[None, None, None, sv]
- e[k] + e[bn]
)
mp4_T_bk_iajb = (
+ np.einsum("iajd, bd -> iajb", t_2_bk_iajb, L_mo[sv, bn, k, sv])
- np.einsum("idja, db -> iajb", t_2_bk_iajb, g_mo[ k, sv, sv, bn])
- np.einsum("ialb, jl -> iajb", t_2_bk_iajb, L_mo[ k, so, so, bn])
+ np.einsum("laib, jl -> iajb", t_2_bk_iajb, g_mo[ k, so, so, bn])
)
energy_mp4_T_by_ckiter += (T_iajb * mp4_T_bk_iajb).sum()
energy_mp4_T_by_ckiter
-0.006712603570763938