5. 简单实现 CEPA 方法#
创建时间：2021-09-04
本文档会简单地实现闭壳层 CEPA(\(n\)) 方法。
CEPA(\(n\)) 不太完整地说，是基于 CCSD 与 CISD 的近似。一般认为这类方法是严格的 CCSD 近似，对 MP2 甚至 MP3 有所提升；但依据推导方式不同，可以是 CISD 的近似，也可以是 CISD 的补充。CEPA(\(n\)) 方法的相对完善的综述可以参考 [1]。
import numpy as np
from pyscf import gto, scf, ci, lib
from pyscf.cc import ccsd
from opt_einsum import contract as einsum
np.set_printoptions(6, suppress=True, linewidth=120)
5.1. CISD 能量#
5.1.1. 准备工作#
我们以键长 0.96 Angstrom、键角 104.5° 的水分子作为研究对象。基组是 cc-pVDZ。
mol = gto.Mole(atom="O; H 1 0.96; H 1 0.96 2 104.5", basis="cc-pVDZ", verbose=0).build()
mf_scf = scf.RHF(mol).run()
mf_ci = ci.CISD(mf_scf).run()
我们首先考察 CISD 能量的计算。回顾到 CISD 计算方程组是
其中，我们规定 \(|\Psi\rangle\) 是半归一化 (Intermediate Normalized) 的波函数：
\(|\Phi_0\rangle\) 是 HF 基态波函数，\(|\Phi_i^a\rangle\) 是单激发波函数、\(|\Phi_{ij}^{ab}\rangle\) 是双激发波函数。需要强调的是，由于是空间轨道，因此这些波函数并不反映物理实在，特别是我们无法考察自旋算符 \(\hat S{}^2\) 的本征态。它们只是用来计算能量的方便的记号。
若要针对一次激发、二次激发分别给出 CISD 方程，那么
5.1.2. 激发系数#
通过 PySCF 的 ci attribute，进而通过 cisdvec_to_amplitudes 函数可以得到 \(|\Psi\rangle\) 的激发系数，但该激发系数并不是归一化的。我们手动将其进行半归一化，得到系数 t1 \(t_i^a\) 与 t2 \(t_{ij}^{ab}\)。
civec = mf_ci.ci / mf_ci.ci[0]
_, t1, t2 = mf_ci.cisdvec_to_amplitudes(civec)
t1.shape, t2.shape
((5, 19), (5, 5, 19, 19))
5.1.3. 能量表达式验证#
首先我们要给出电子积分。在有足够内存的情况下，可以用 ccsd 的 _make_eris_incore 实现：
eris = ccsd._make_eris_incore(mf_ci, mf_scf.mo_coeff)
我们通篇文档不关心涉及 \(\hat H\) 的计算过程。我们暂时只需要其中的
该张量可以通过 eris.ovov 调出，维度是 \((i, a, j, b)\)：
eris.ovov.shape
(5, 19, 5, 19)
对于半归一化方法，不论是 MP2, CEPA(\(n\)), CI(S)D 或 CC(S)D，下式在闭壳层下总是成立的：
2 * einsum("iajb, ijab ->", eris.ovov, t2) - einsum("iajb, ijba ->", eris.ovov, t2)
-0.2053384394211093
我们再次回顾 PySCF 计算得到的 CISD 相关能是下述非常接近的结果：
mf_ci.e_corr
-0.20533844297533488
5.2. CEPA(\(n\)) 计算#
5.2.1. CEPA(\(n\)) 原理#
CEPA(\(n\)) 方程组是
其中，
对电子能定义为
def pair_energy(eris_ovov, t2):
return 2 * einsum("iajb, ijab -> ij", eris_ovov, t2) - einsum("iajb, ijba -> ij", eris_ovov, t2)
CEPA(2) 由于不确定其 \(B_i\) 的定义，因此这里不作实现。当然，CEPA(3) 的定义也可能存在疑问。
def cepa_shift(cepa_n, t1, t2, e_ij):
e_i, e_j = e_ij.sum(axis=1), e_ij.sum(axis=0)
A_ij, B_i = np.zeros_like(e_ij), np.zeros_like(e_i)
if cepa_n == 0:
pass
elif cepa_n == 1:
A_ij = 0.5 * (e_i[:, None] + e_j[None, :])
B_i = e_i
elif cepa_n == 3:
A_ij = e_i[:, None] + e_j[None, :] - e_ij
B_i = - e_ij.diagonal() + 2 * e_i
else:
raise ValueError("cepa_n value error for " + cepa_n)
return A_ij, B_i
5.2.2. 迭代法求取 CISD 系数#
我们首先考虑二次激发系数 \(t_{ij}^{ab}\) 的求取。
如果我们记等式左为 \(\mathscr{v}_{ij}^{ab}\)，那么该式整理为
但需要注意，这不是一个好的激发系数 \(t_{ij}^{ab}\) 的更新策略。这里指出，由于 \(\langle \Phi_{ij}^{ab} | \hat H - E_\textsf{HF} | \Phi_0 \rangle = -D_{ij}^{ab}\) 是 \(\mathscr{v}_{ij}^{ab}\) 的重要贡献项，因此不妨将上式写为
上式是我们实际会使用到的激发系数 \(t_{ij}^{ab}\) 的更新策略。
在迭代计算激发系数时，我们需要给定初始激发系数。一般来说使用 MP2 激发系数即可：
其中，d1 \(D_i^a = \varepsilon_i - \varepsilon_a\)，d2 \(D_{ij}^{ab} = \varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b\)。
def corr_cepa(cepa_n):
# Prepare D_i^a, D_ij^ab
d0, d1, d2 = mf_ci.cisdvec_to_amplitudes(mf_ci.make_diagonal(eris))
d1 = d0 - d1; d2 = d0 - d2
# Prepare initial t_i^a, t_ij^ab
t0, t1, t2 = 1, np.zeros_like(d1), eris.ovov.swapaxes(1, 2) / d2
civec = mf_ci.amplitudes_to_cisdvec(t0, t1, t2)
# Iteration
for it in range(100):
t1_old, t2_old = t1, t2
_, v1, v2 = mf_ci.cisdvec_to_amplitudes(mf_ci.contract(civec, eris))
e_ij = pair_energy(eris.ovov, t2)
A_ij, B_i = cepa_shift(cepa_n, t1, t2, e_ij)
t1 = t1_old + (v1 - B_i[:, None] * t1_old) / d1
t2 = t2_old + (v2 - A_ij[:, :, None, None] * t2_old) / d2
e_corr = pair_energy(eris.ovov, t2).sum()
civec = mf_ci.amplitudes_to_cisdvec(t0, t1, t2)
if np.linalg.norm(t1 - t1_old) + np.linalg.norm(t2 - t2_old) < 1e-5:
print("Total Iterations: ", it)
print("Correlation Energy: ", e_corr)
break
corr_cepa(0)
Total Iterations: 26
Correlation Energy: -0.2167752177602909
corr_cepa(1)
Total Iterations: 38
Correlation Energy: -0.21352333911480398
corr_cepa(3)
Total Iterations: 53
Correlation Energy: -0.21129784897010107