5. XYG3 能量与非自洽 Fock 矩阵¶

在这一节中，我们会以 XYG3 能量计算与非自洽 Fock 矩阵的构建为例，讲述如何重复 B3LYP 能量、MP2 能量、XYG3 能量，以及其中更多的细节．尽管一般来说，学习量化程序的第一步是 SCF 程序的实现；但 Post-HF 方法的程序实现实际上未必需要这一基础．作者希望这一节可以为后面的 CP 方程提供比较充分的实践基础．

这一节统一假定 Restricted 参考态．

5.1. 准备工作¶

5.1.1. 分子计算¶

我们现在做一些必要的前期准备．下面定义的众多量会频繁使用．

In [1]:

# 环境搭建
import psi4
import numpy as np

# 简化矩阵输出
np.set_printoptions(5, linewidth=100, suppress=True)

# 输出文件
psi4.set_output_file("output.dat", False)

# 设置内存 0.5 GB
psi4.set_memory(int(5e8))

# 设置分子坐标
mol = psi4.geometry("""
    O  0.000000000000    -0.000000000000    -0.079135765807
    H  0.000000000000     0.707106781187     0.627971015380
    H  0.000000000000    -0.707106781187     0.627971015380
    symmetry c1
""")

# 设置计算选项
psi4.set_options({
    'basis':               '6-31g',
    'scf_type':            'pk',
    'mp2_type':            'conv',
    'e_convergence':        1e-8,
    'd_convergence':        1e-8,
    'dft_spherical_points': 590,
    'dft_radial_points':    99,
})

In [2]:

# 波函数信息，包括 HF 与 B3LYP
# B3LYP 在 Gaussian 中需要设置 Int=UltraFine，并且只能对到小数点第 7 位
hf_e, hf_wfn = psi4.energy("HF", molecule=mol, return_wfn=True)
scf_e, scf_wfn = psi4.energy("B3LYP", molecule=mol, return_wfn=True)
mp2_e = psi4.energy("MP2", molecule=mol, return_wfn=False)

我们首先可以验证一下 Gaussian 的结果与上述的计算结果是否能对到第六位小数．Gaussian 输入卡：HF, MP2, B3LYP

In [3]:

psi4.compare_values(hf_e, -75.9697009555, 6, 'HF Energy')
psi4.compare_values(scf_e, -76.3771828949, 6, 'B3LYP Energy')
psi4.compare_values(mp2_e, -0.76104035644031E+02, 6, 'MP2 Energy')

        HF Energy.........................................................PASSED
        B3LYP Energy......................................................PASSED
        MP2 Energy........................................................PASSED

Out[3]:

True

5.1.2. 中间变量¶

下面的变量是从波函数信息中导出的量．由于我们主要使用 B3LYP 的计算结果，因此不作额外标记的量一般都是 B3LYP 的计算结果．一些电子积分的结果对于 HF 与 B3LYP 是相同的．同时，由于是 Restricted 参考态，因此 \(\alpha\) 与 \(\beta\) 自旋的导出矩阵的结果应当是相同的．

In [4]:

# 电子积分引擎
# 由于我们现在所使用的电子积分引擎是与方法无关的，因此这里无所谓使用哪个方法的波函数．
mints = psi4.core.MintsHelper(scf_wfn.basisset())
# 积分
T = np.asarray(mints.ao_kinetic())  # AO 基组动能积分
V = np.asarray(mints.ao_potential())  # AO 基组电子-核势能积分
eri = np.asarray(mints.ao_eri())  # AO 基组双电子排斥积分

In [5]:

# 数值量
nocc = scf_wfn.nalpha()
nbf = mints.nbf()
nmo = scf_wfn.nmo()
nvir = nmo - nocc
print("nmo = nbf? ", nmo == nbf)

nmo = nbf?  True

In [6]:

# HF 导出变量
F_hf = np.asarray(hf_wfn.Fa())  # AO 基组 Fock 矩阵
D_hf = np.asarray(hf_wfn.Da())  # AO 基组单电子密度
C_hf = np.asarray(hf_wfn.Ca())  # AO <-> MO 轨道系数 C_{up}
Co_hf = np.asarray(hf_wfn.Ca_subset('AO', 'OCC'))  # 占据轨道系数
Cv_hf = np.asarray(hf_wfn.Ca_subset('AO', 'VIR'))  # 未占轨道系数
e_hf = np.asarray(hf_wfn.epsilon_a())  # 轨道能级
eo_hf = np.asarray(hf_wfn.epsilon_a_subset('AO', 'OCC'))  # 占据轨道能级
ev_hf = np.asarray(hf_wfn.epsilon_a_subset('AO', 'VIR'))  # 未占轨道能级

In [7]:

# B3LYP 导出变量
F = np.asarray(scf_wfn.Fa())  # AO 基组 Fock 矩阵
D = np.asarray(scf_wfn.Da())  # AO 基组单电子密度
C = np.asarray(scf_wfn.Ca())  # AO <-> MO 轨道系数 C_{up}
Co = np.asarray(scf_wfn.Ca_subset('AO', 'OCC'))  # 占据轨道系数
Cv = np.asarray(scf_wfn.Ca_subset('AO', 'VIR'))  # 未占轨道系数
e = np.asarray(scf_wfn.epsilon_a())  # 轨道能级
eo = np.asarray(scf_wfn.epsilon_a_subset('AO', 'OCC'))  # 占据轨道能级
ev = np.asarray(scf_wfn.epsilon_a_subset('AO', 'VIR'))  # 未占轨道能级
# DFT 相关
Vxc = np.asarray(scf_wfn.Va())  # AO 基组 V_xc 密度泛函一次导数矩阵
V_pot = scf_wfn.V_potential()  # DFT 积分引擎

5.2. HF Fock 矩阵重复¶

5.2.1. 普通做法¶

我们首先验证一下最简单的 HF 方法下的 Fock 矩阵的构造．我们知道，Restricted HF 的 AO 基组 Fock 矩阵可以写作

\begin{align} F_{\mu \nu} [P_{\mu \nu}] = T_{\mu \nu} + V_{\mu \nu}^\mathrm{ext} + 2 J_{\mu \nu} [P_{\mu \nu}] - K_{\mu \nu} [P_{\mu \nu}] \end{align}

其中，

\begin{align} T_{\mu \nu} &= \langle \mu | - \frac{1}{2} \nabla^2 | \nu \rangle \\ V_{\mu \nu}^\mathrm{ext} &= \langle \mu | v^\mathrm{ext} | \nu \rangle \\ J_{\mu \nu} [P_{\mu \nu}] &= (\mu \nu | i i) = (\mu \nu | \kappa \lambda) C_{\kappa i} C_{\lambda i} = (\mu \nu | \kappa \lambda) P_{\mu \nu} \\ K_{\mu \nu} [P_{\mu \nu}] &= (\mu i | \nu i) = C_{\kappa i} (\mu \kappa | \nu \lambda) C_{\lambda i} = (\mu \kappa | \nu \lambda) P_{\mu \nu} \\ \end{align}

同时，AO 基组 Fock 矩阵到 MO 基组 Fock 矩阵的变换可以通过下式确定：

\begin{align} F_{pq} [P_{\mu \nu}] = C_{\mu p} F_{\mu \nu} [P_{\mu \nu}] C_{\nu q} \end{align}

对于 Hartree-Fock，MO 基组的 Fock 矩阵应当是对角阵，且对角值为轨道能量．我们就验证这些性质．

In [8]:

# 构建 AO Fock 矩阵
F_hf_my = np.zeros_like(F_hf)  # 构建零矩阵
F_hf_my += T  # 动能部分
F_hf_my += V  # 势能部分
F_hf_my += 2 * np.einsum("uvkl, ki, li -> uv", eri, Co_hf, Co_hf, optimize=True)  # Coulomb 积分部分
F_hf_my -= np.einsum("ukvl, ki, li -> uv", eri, Co_hf, Co_hf, optimize=True)  # Exchange 积分部分
# 构造 J, K 矩阵的两行代码还可以写作
#   F_hf_my += 2 * np.einsum("uvkl, kl -> uv", eri, D_hf)  # Coulomb 积分部分
#   F_hf_my -= np.einsum("ukvl, kl -> uv", eri, D_hf)  # Exchange 积分部分

# 检验 AO Fock 矩阵是否构建正确
print("AO Fock matrix allclose: ", np.allclose(F_hf_my, F_hf))
# 检验 MO Fock 矩阵对角元是否是本征值
print("MO Fock equal to eigenvalues: ", np.allclose(np.diag(C_hf.T @ F_hf_my @ C_hf), e_hf))

AO Fock matrix allclose:  True
MO Fock equal to eigenvalues:  True

由此，Hartree-Fock 能量也可以直接获得：

\begin{align} E^\textsf{HF} [P_{\mu \nu}] = P_{\mu \nu} (T_{\mu \nu} + V_{\mu \nu}^\mathrm{ext} + F_{\mu \nu} [P_{\mu \nu}]) + E^\mathrm{Nuc} \end{align}

In [9]:

hf_e_my = (D_hf * (T + V + F_hf_my)).sum() + mol.nuclear_repulsion_energy()
psi4.compare_values(hf_e, hf_e_my, 6, 'Constructed HF Energy')

        Constructed HF Energy.............................................PASSED

Out[9]:

True

提示

在这里，密度矩阵的定义是 \(P_{\mu \nu} = C_{\mu i} C_{\nu i}\)，与 Szabo (3.145) 的定义少一倍．因此，Hartree-Fock 能量中电子所贡献的部分与 (3.184) 的公式也差一倍．

提示

很多变量都会是单电子密度矩阵的相关变量，但这会把公式弄得比较糟糕．在之后，我们可能会较少地写 \([P_{\mu \nu}]\) 的记号，除非比较有必要．

5.2.2. 使用 JK 引擎积分¶

在后面叙述 CP 方程时，我们会经常地使用 JK 引擎．尽管这一节不会大量使用该引擎，但对它的了解将会方便后面对 CP 方程程序化的理解．其使用方式如下．

JK 引擎解决的问题是，当传入两个系数矩阵 \(C_{\kappa i}^\mathrm{left}\) 与 \(C_{\lambda i}^\mathrm{right}\) 时，给出其对应原子轨道基组的 Coulomb 与 Exchange 积分

\begin{align} J_{\mu \nu} &= (\mu \nu | \kappa \lambda) C_{\kappa i}^\mathrm{left} C_{\lambda i}^\mathrm{right} \\ K_{\mu \nu} &= (\mu \kappa | \nu \lambda) C_{\kappa i}^\mathrm{left} C_{\lambda i}^\mathrm{right} \end{align}

显然这是一个更广泛的解决方案．若局限于 Hartree-Fock 的 Coulomb 与 Exchange 矩阵，我们也可以使用 JK 引擎解决：只要我们传入的 \(C_{\kappa i}^\mathrm{left}\) 与 \(C_{\lambda i}^\mathrm{right}\) 都是 \(C_{\mu i}\) 即可．

In [10]:

# JK 引擎初始化
jk_hf = psi4.core.JK.build(hf_wfn.basisset())
jk_hf.initialize()

# 设定左矢与右矢轨道系数 (实际上相当于设定传入单电子密度)
jk_hf.C_left_add(psi4.core.Matrix.from_array(Co_hf))
jk_hf.C_right_add(psi4.core.Matrix.from_array(Co_hf))

# JK 引擎计算与结果输出
jk_hf.compute()
J_hf_my = np.asarray(jk_hf.J()[0])
K_hf_my = np.asarray(jk_hf.K()[0])
D_hf_my = np.asarray(jk_hf.D()[0])

# 比较结果
print(" Coulomb matrix allclose: ",
      np.allclose(np.einsum("uvkl, ki, li -> uv", eri, Co_hf, Co_hf, optimize=True), J_hf_my))
print("Exchange matrix allclose: ",
      np.allclose(np.einsum("ukvl, ki, li -> uv", eri, Co_hf, Co_hf, optimize=True), K_hf_my))
print(" Density matrix allclose: ",
      np.allclose(Co_hf @ Co_hf.T, D_hf_my))

 Coulomb matrix allclose:  True
Exchange matrix allclose:  True
 Density matrix allclose:  True

5.3. HF 偶极矩重复¶

这里简单介绍 HF 偶极矩的计算的方式．偶极矩的投影在 Hartree-Fock 下非常容易计算，只需要单电子密度、偶极积分与原子核位置就可以获得：

\begin{equation} \mu_f = P_{\mu \nu} ( \nu | -f | \mu ) + Z_A f_A \end{equation}

其中 \(f\) 代表坐标取向 \(x, y, z\)，\(Z_A\) 与 \(f_A\) 分别代表原子 \(A\) 的电荷数值与其坐标在取向 \(f\) 上的投影．偶极矩本身为矢量：

\begin{equation} \boldsymbol{\mu} = (\mu_x, \mu_y, \mu_z) \end{equation}

5.3.1. 电子云贡献¶

我们只需要拿出偶极矩积分即可．由于偶极矩分三个方向，因此偶极矩的积分会储存在一个列表中．我们将这个列表中的内容替换为 NumPy 向量对象．需要注意这些积分已经将电子的负电荷所产生的负号包含进去了．

In [11]:

# Hartree-Fock 偶极矩 AO 基组积分
dip_mat_hf = mints.ao_dipole()
for ind in range(3):
    dip_mat_hf[ind] = np.asarray(dip_mat_hf[ind])

随后我们就将这些矩阵与密度矩阵作元素相乘，即得到电子云贡献的偶极矩．注意到因为 RHF 下密度矩阵的定义，因此需要对结果乘上两倍．该偶极矩的物理量是原子单位，即电子电荷乘上玻尔半径．

In [12]:

dip_hf = []
for ind in range(3):
    dip_hf.append((dip_mat_hf[ind] * D_hf).sum() * 2)
dip_hf

Out[12]:

[-9.32554394266127e-17, 1.2847475302363298e-08, -0.04940294102679377]

5.3.2. 原子核贡献¶

我们知道，由于空间坐标涉及到坐标系的选取，因此在计算过程中，电子云相对于原点的位置选取上的改变会直接对偶极矩的值产生改变．不同的原点选取不会影响作为偶极矩对电场导数的极化率，或者对偶极矩空间坐标的导数，但会影响偶极矩本身的数值．

在 Psi4 中，分子在构建过程中就会进行坐标转换或平移，因此我们需要知道更新之后的原子坐标．这可以通过下述代码实现 (单位为玻尔半径)：

In [13]:

mol_geo = np.asarray(mol.geometry())
mol_geo

Out[13]:

array([[ 0.     ,  0.     , -0.14954],
       [ 0.     ,  1.33624,  1.18669],
       [ 0.     , -1.33624,  1.18669]])

乍看之下似乎变化不大，但实际上有微小区别；如果我们在输入分子坐标的地方对分子进行平移，则会看到比较明显的区别．同时，每个原子的打包在列表的原子核电荷可以通过下述代码给出 (单位为电子电荷)：

In [14]:

neu_charge = []
for neu in range(mol_geo.shape[0]):
    neu_charge.append(mol.charge(neu))
neu_charge

Out[14]:

[8.0, 1.0, 1.0]

最后，我们就给出原子核对偶极矩的贡献．

In [15]:

dip_mol = [0., 0., 0.]
for neu in range(mol_geo.shape[0]):
    for ind in range(3):
        dip_mol[ind] += neu_charge[neu] * mol_geo[neu][ind]
dip_mol

Out[15]:

[0.0, 0.0, 1.1770270745569982]

5.3.3. 总偶极矩¶

总偶极矩就是电子云贡献与原子核贡献的简单加和．

In [16]:

dip_tot_hf = [v1 + v2 for (v1, v2) in zip(dip_hf, dip_mol)]
dip_tot_hf

Out[16]:

[-9.32554394266127e-17, 1.2847475302363298e-08, 1.1276241335302044]

偶极矩在 SCF 方法下，包括一般的杂化 DFT 近似方法，都是比较容易实现的．但对于 MP2 等方法，则需要构建其一阶约化密度，并且这种密度的构建方式是通过对能量求导数获得的．因此，为获得 MP2 方法的偶极矩，我们需要求解 Z-Vector 方程．这在之后的教程中会提及．

5.4.2. MP2 能量重复¶

5.4.1. Kutzelnigg-Mukherjee (KM) 记号¶

在重复 MP2 能量前，我们需要确定一些 Restricted 参考态下常用的记号．这里采用 KM 记号 (或见 Psi4NumPy 教程 8a) 书写 (但我尚不习惯这些记号，因此可能有写得不太对的地方)．与 MP2 能量或 CP 方程有关的变量通常只有双电子积分 (ERI)，以及双轨道激发振幅 (Amplitude)；因此只给出与这两者有关的记号定义．

\begin{align} g_{ij}^{ab} &= \langle i j | a b \rangle \\ \bar g_{ij}^{ab} &= \langle i j \Vert a b \rangle = \langle i j | a b \rangle - \langle i j | b a \rangle \\ \mathscr{E}_{ab}^{ij} &= \varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b \\ t_{ij}^{ab} &= g_{ij}^{ab} (\mathscr{E}_{ab}^{ij})^{-1} \\ \bar t_{ij}^{ab} &= \bar g_{ij}^{ab} (\mathscr{E}_{ab}^{ij})^{-1} = t_{ij}^{ab} - t_{ij}^{ba} \end{align}

其中，在字母上不带横的记号通常可以表示相反自旋 (OS, Opposite Spin) 的贡献，而带一横的记号则可以表示相同自旋 (SS, Same Spin) 的贡献．

5.4.2. MP2 能量重复¶

我们可以将 MP2 能量写作下述形式：

\begin{equation} E_\mathrm{corr}^\mathsf{MP2} = (t_{ij}^{ab})^2 \mathcal{E}_{ab}^{ij} + \frac{1}{2} (\bar t_{ij}^{ab})^2 \mathcal{E}_{ab}^{ij} \end{equation}

下面就来验证 MP2 能量的构造．

In [17]:

# 构建轨道能之差的张量
d_hf = (eo_hf.reshape(-1, 1, 1, 1) - ev_hf.reshape(-1, 1, 1) + eo_hf.reshape(-1, 1) - ev_hf)
# 构建分子轨道积分 <ij|ab> = (ia|jb) = g_{ij}^{ab} -> g[i, a, j, b]
g_hf = np.einsum("ui, va, uvkl, kj, lb -> iajb", Co_hf, Cv_hf, eri, Co_hf, Cv_hf, optimize=True)
# 构建轨道对激发振幅
t_hf = g_hf / d_hf
# 计算 MP2 能量
mp2_corr_my = (t_hf ** 2 * d_hf + 0.5 * (t_hf - t_hf.swapaxes(1,3)) ** 2 * d_hf).sum()
# 对比方才计算的 MP2 能量
psi4.compare_values(mp2_corr_my, mp2_e - hf_e, 6, 'Constructed MP2 Correlation')

        Constructed MP2 Correlation.......................................PASSED

Out[17]:

True

事实上，ERI 积分可以直接从 Psi4 的积分引擎给出；结果是一样的．

In [18]:

# 这里传入的变量必须是 psi4.core.Matrix 型
Co_hf_psi4 = hf_wfn.Ca_subset('AO', 'OCC')
Cv_hf_psi4 = hf_wfn.Ca_subset('AO', 'VIR')
# 对比 psi4 的积分引擎与手算分别给出的 MO ERI 张量
print("mints.mo_eri vs. hand-made mo-eri tensor: ",
      np.allclose(mints.mo_eri(Co_hf_psi4, Cv_hf_psi4, Co_hf_psi4, Cv_hf_psi4), g_hf))

mints.mo_eri vs. hand-made mo-eri tensor:  True

5.5. B3LYP Fock 矩阵重复¶

现在我们将重复 B3LYP 的 Fock 矩阵．相对于 Hartree-Fock 方法，B3LYP 的 Fock 矩阵多出一项交换相关势能，同时其 K 积分的比例由杂化泛函的精确交换比例所确定：

\begin{equation} F_{\mu \nu} = T_{\mu \nu} + V_{\mu \nu}^\mathrm{ext} + 2 J_{\mu \nu} - c_\mathrm{x} K_{\mu \nu} + V_{\mu \nu}^\mathrm{xc} \end{equation}

\(V_{\mu \nu}^\mathrm{xc}\) 可以由 DFT 积分引擎给出，该积分引擎已经在中间矩阵处定义．事实上，B3LYP 波函数 scf_wfn 已经储存了交换相关势，其值从 scf_wfn.Va() 调出．表面上，这好像是该波函数所对应的势能；但这个势能不包含外势，或者说在普通分子体系下，核与电子的相互作用．因此，如果我们取看 hf_wfn.Va()，就会发现这是个零矩阵．我们可以用 DFT 积分引擎算出的交换相关势与波函数的交换相关势进行比对来验证结果．

In [19]:

# 创建一个零矩阵，它将储存交换相关势
Vxc_my = psi4.core.Matrix(nbf, nbf)
# 在计算交换相关势前，需要指定该势能所处的密度环境
V_pot.set_D([scf_wfn.Da()])
# 计算交换相关势，零矩阵 Vxc_my 在运行后将会被替换为非零矩阵
V_pot.compute_V([Vxc_my])
# 比对结果
print("wfn.V_potential vs. wfn.Va: ",
      np.allclose(Vxc_my, Vxc))

wfn.V_potential vs. wfn.Va:  True

那么仿照 Hartree-Fock，B3LYP 的 Fock 矩阵可以容易地给出．我们这里利用了 B3LYP 的 \(c_\mathrm{x} = 0.2\)．

In [20]:

# 构建 AO Fock 矩阵
F_my = np.zeros_like(F)  # 构建零矩阵
F_my += T  # 动能部分
F_my += V  # 势能部分
F_my += 2 * np.einsum("uvkl, kl -> uv", eri, D)  # Coulomb 积分部分
F_my -= 0.2 * np.einsum("ukvl, kl -> uv", eri, D)  # Exchange 积分部分
F_my += Vxc_my.np  # DFT 交换相关势

# 检验 AO Fock 矩阵是否构建正确
print("AO Fock matrix allclose: ", np.allclose(F_my, F))
# 检验 MO Fock 矩阵对角元是否是本征值
print("MO Fock equal to eigenvalues: ", np.allclose(np.diag(C.T @ F_my @ C), e))

AO Fock matrix allclose:  True
MO Fock equal to eigenvalues:  True

有了 Fock 矩阵后，我们似乎也可以试着给出 B3LYP 的能量；但需要注意到，DFT 的交换相关能的计算不是靠其 AO 基组的交换相关势给出的，而需要单独地算交换相关能才行．幸运的是，这部分能量已经在调用 DFT 积分引擎时已经计算完毕．我们只需要 V_pot.quadrature_values() 调用其结果即可．利用的公式是

\begin{align} E^\textsf{B3LYP} = P_{\mu \nu} (2 T_{\mu \nu} + 2 V_{\mu \nu}^\mathrm{ext} + 2 J_{\mu \nu} - c_\mathrm{x} K_{\mu \nu}) + E^\mathrm{xc} [P_{\mu \nu}] + E^\mathrm{Nuc} \end{align}

In [21]:

# 动能、势能、Coulomb 积分、Exchange 积分贡献
scf_e_my = (D * (2 *T + 2 * V + 2 * np.einsum("uvkl, kl -> uv", eri, D)
                 - 0.2 * np.einsum("ukvl, kl -> uv", eri, D))).sum()
# 交换相关能
scf_e_my += V_pot.quadrature_values()["FUNCTIONAL"]
# 核排斥
scf_e_my += mol.nuclear_repulsion_energy()

# 与 Psi4 计算结果比对
psi4.compare_values(scf_e_my, scf_e, 6, 'Constructed B3LYP Correlation')

        Constructed B3LYP Correlation.....................................PASSED

Out[21]:

True

5.6. XYG3 非自洽 Fock 矩阵与能量¶

5.6.1. 非自洽 Fock 矩阵与能量的定义¶

我们知道，XYG3 基于 B3LYP 自洽的密度与轨道，使用 XYG3 密度泛函进行能量计算．在后面，我们会记 \(\mathrm{s}\) 为自洽泛函记号 (B3LYP)，而 \(\mathrm{n}\) 为非自洽泛函记号 (XYG3)．其能量形式可以表达为

\begin{align} E^\textsf{XYG3} = P_{\mu \nu}^\mathrm{s} (2 T_{\mu \nu} + 2 V_{\mu \nu}^\mathrm{ext} + 2 J_{\mu \nu} [P_{\mu \nu}^\mathrm{s}] - c_\mathrm{x}^\mathrm{n} K_{\mu \nu} [P_{\mu \nu}^\mathrm{s}]) + E^\mathrm{xc, n} [P_{\mu \nu}^\mathrm{s}] + c_\mathrm{c}^\mathrm{n} E^\textsf{PT2} + E^\mathrm{Nuc} \end{align}

因此，其能量表达式的大部分项都非常容易获得；其中的 \(E^\textsf{PT2}\) 一般是通过在 MP2 公式中代入 B3LYP 轨道系数计算而得；比较麻烦的项是 \(E^\mathrm{xc, n} [P_{\mu \nu}^\mathrm{s}]\)．同时，XYG3 的非自洽 Fock 矩阵的形式为

\begin{equation} F_{\mu \nu}^\mathrm{n} = T_{\mu \nu} + V_{\mu \nu}^\mathrm{ext} + 2 J_{\mu \nu} [P_{\mu \nu}^\mathrm{s}] - c_\mathrm{x}^\mathrm{n} K_{\mu \nu} [P_{\mu \nu}^\mathrm{s}] + V_{\mu \nu}^\mathrm{xc, n} [P_{\mu \nu}^\mathrm{s}] \end{equation}

比较麻烦的是 \(V_{\mu \nu}^\mathrm{xc, n} [P_{\mu \nu}^\mathrm{s}]\)．

5.6.2. XYG3 泛函定义¶

为了计算上述两个比较麻烦的量，第一步是在程序中对 XYG3 的非自洽泛函进行程序中的定义．XYG3 泛函，或第五阶泛函 (按 Predew 的 Jacob 阶梯说法) 的定义是

\begin{equation} E_\mathrm{xc}^\textsf{R5} = E_\mathrm{xc}^\textsf{LDA} + c_1 (E_\mathrm{x}^\textsf{exact} - E_\mathrm{x}^\textsf{LDA}) + c_2 \Delta E_\mathrm{x}^\textsf{GGA} + c_3 (E_\mathrm{c}^\textsf{PT2} - E_\mathrm{c}^\textsf{LDA}) + c_4 \Delta E_\mathrm{c}^\textsf{GGA} \end{equation}

对应到程序中，每一项的系数则展开为

\begin{equation} E_\mathrm{xc}^\textsf{R5} = (1 - c_1 - c_2) E_\mathrm{x}^\textsf{LDA} + c_2 E_\mathrm{x}^\textsf{GGA} + (1 - c_3 - c_4) E_\mathrm{c}^\textsf{LDA} + c_4 E_\mathrm{c}^\textsf{GGA} + c_1 E_\mathrm{x}^\textsf{exact} + c_3 E_\mathrm{c}^\textsf{PT2} \end{equation}

对于 XYG3，其系数的确定是

\begin{equation} c_1 = 0.8033, \quad c_2 = 0.2107, \quad c_3 = 0.3211, \quad c_4 = 1 - c_3 \end{equation}

因此，XYG3 的泛函形式为

\begin{equation} E_\mathrm{xc}^\textsf{XYG3} = -0.0140 E_\mathrm{x}^\textsf{LDA} + 0.2107 E_\mathrm{x}^\textsf{GGA} + 0.6789 E_\mathrm{c}^\textsf{GGA} + 0.8033 E_\mathrm{x}^\textsf{exact} + 0.3211 E_\mathrm{c}^\textsf{PT2} \end{equation}

在程序中，真正为 DFT 积分引擎所关心的部分是上式的前三项；而后两项应该在 HF 或 Post-HF 框架下完成．这里，我们定义

\begin{equation} E^\mathrm{xc, n} = -0.0140 E_\mathrm{x}^\textsf{LDA} + 0.2107 E_\mathrm{x}^\textsf{GGA} + 0.6789 E_\mathrm{c}^\textsf{GGA} \end{equation}

在 Psi4 中，泛函形式的确定是通过新建一个函数来进行．其确定过程如下：

In [22]:

def build_xyg3_nc_superfunctional(name, npoints, deriv, restricted):

    # 新建一个空白泛函，最后返回该泛函
    sup = psi4.core.SuperFunctional.blank()

    # 对泛函进行基本描述
    sup.set_name('XYG3NC')
    sup.set_description('    XYG3 Non-Consistent Functional without MP2 Part\n')

    # -0.0140 * LDA Exchange
    # -0.0140 = 1 - 0.8033 - 0.2107
    lda_x = psi4.core.LibXCFunctional("XC_LDA_X", restricted)
    lda_x.set_alpha(-0.0140)
    sup.add_x_functional(lda_x)

    # 0.2107 * B88 Exchange
    gga_x = psi4.core.LibXCFunctional("XC_GGA_X_B88", restricted)
    gga_x.set_alpha(0.2107)
    sup.add_x_functional(gga_x)

    # 0.6789 * LYP Correlation
    # 0.6789 = 1 - 0.3211
    lyp_c = psi4.core.LibXCFunctional("XC_GGA_C_LYP", restricted)
    lyp_c.set_alpha(0.6789)
    sup.add_c_functional(lyp_c)

    # 0.8033 Exact Exchange
    sup.set_x_alpha(0.8033)
    # 我想上述的精确交换在这里不设好像也可以

    return sup

5.6.3. XYG3 交换相关能与交换相关势¶

随后我们就可以对上面两个最麻烦的部分进行计算了．我们的任务是拿到 DFT 积分引擎；尽管原则上不需要先进行自洽场计算，但我还没找到一个可以直接可以使用的 DFT 积分引擎的方案 (应当是交换相关势函数设得不好)，因此就先进行一次 DFT 计算，从非自洽泛函的波函数对象中拿出 DFT 引擎，再代入自洽密度，得到我们想要的结果．

提示

这里计算非自洽泛函的自洽场纯粹是权宜之计！只为获得 DFT 积分引擎．

In [23]:

nscf_e, nscf_wfn = psi4.energy("SCF", dft_functional=build_xyg3_nc_superfunctional, return_wfn=True)

In [24]:

# 非自洽势能定为 Vn，其代入的密度是自洽场泛函的 AO 基组密度
Vn = psi4.core.Matrix(nbf, nbf)
Vn_pot = nscf_wfn.V_potential()
Vn_pot.set_D([scf_wfn.Da()])
Vn_pot.compute_V([Vn])
# 在获得非自洽交换相关势后，同时可以得到交换相关能
Vn_pot.quadrature_values()["FUNCTIONAL"]

Out[24]:

-2.004883653728054

5.6.4. XYG3 能量¶

由此，我们可以获得 XYG3 能量．其计算过程如下：

In [25]:

#** XYG3 能量 (不包含 PT2)
# 一般 SCF 贡献
xyg3_e = (D * (2 * T + 2 * V + 2 * np.einsum("uvkl, kl -> uv", eri, D)
               - 0.8033 * np.einsum("ukvl, kl -> uv", eri, D))).sum()
# 交换相关贡献
xyg3_e += Vn_pot.quadrature_values()["FUNCTIONAL"]
# 核排斥贡献
xyg3_e += mol.nuclear_repulsion_energy()

#** PT2 部分
# 构建 MO 基组 ERI
Co_psi4 = scf_wfn.Ca_subset('AO', 'OCC')
Cv_psi4 = scf_wfn.Ca_subset('AO', 'VIR')
g_psi4 = mints.mo_eri(Co_psi4, Cv_psi4, Co_psi4, Cv_psi4)
g = np.asarray(g_psi4)
# 构建轨道能之差的张量
d = (eo.reshape(-1, 1, 1, 1) - ev.reshape(-1, 1, 1) + eo.reshape(-1, 1) - ev)
# 构建轨道对激发振幅
t = g / d
# 计算 XYG3 PT2 能量
xyg3_pt2 = 0.3211 * (t ** 2 * d + 0.5 * (t - t.swapaxes(1,3)) ** 2 * d).sum()

# 最终 XYG3 能量
xyg3_e += xyg3_pt2
xyg3_e

Out[25]:

-76.28239363067564

我们与 XYG3 在该分子下的参考值进行比对，在 6 位小数内正确．

In [26]:

psi4.compare_values(xyg3_e, -0.76282393305943E+02, 6, 'XYG3 Energy')

        XYG3 Energy.......................................................PASSED

Out[26]:

True

5.6.5. XYG3 非自洽 Fock 矩阵¶

我们按照与构造 B3LYP 的 Fock 矩阵类似的方法，就可以构造 XYG3 非自洽 Fock 矩阵．这里不提供参考值比对的情况了．我们给出下述代码作为生成 AO 基组的 XYG3 非自洽 Fock 矩阵：

In [27]:

Fn = T + V
Fn += 2 * np.einsum("uvkl, kl -> uv", eri, D)
Fn -= 0.8033 * np.einsum("ukvl, kl -> uv", eri, D)
Fn += Vn

以 B3LYP 自洽轨道为组成的 MO 基组的 XYG3 非自洽 Fock 矩阵可以以下述代码给出．我们可以发现，该 Fock 矩阵并不是对角矩阵，但其对角元不大；同时其对角元的值与轨道能也基本相近．该矩阵的占据-非占 (occupied-virtual) 部分将会是 XYG3 的 CP 方程求解过程中非常重要的项．

In [28]:

Fn_mo = C.T @ Fn @ C
Fn_mo

Out[28]:

array([[-20.23753,  -0.02323,   0.     ,  -0.01456,   0.     ,  -0.01381,  -0.     ,   0.     ,
         -0.00786,   0.     ,  -0.02656,  -0.     ,  -0.03727],
       [ -0.02323,  -1.27763,  -0.     ,  -0.03914,  -0.     ,   0.00513,  -0.     ,   0.     ,
          0.01026,   0.     ,  -0.01151,  -0.     ,   0.01092],
       [  0.     ,  -0.     ,  -0.6303 ,   0.     ,  -0.     ,   0.     ,  -0.00611,  -0.00518,
          0.     ,   0.     ,   0.     ,   0.02084,  -0.     ],
       [ -0.01456,  -0.03914,   0.     ,  -0.54857,   0.     ,   0.00842,  -0.     ,  -0.     ,
         -0.00336,  -0.     ,   0.00734,   0.     ,  -0.01383],
       [  0.     ,  -0.     ,  -0.     ,   0.     ,  -0.468  ,   0.     ,   0.     ,   0.     ,
         -0.     ,   0.00389,   0.     ,   0.     ,   0.     ],
       [ -0.01381,   0.00513,   0.     ,   0.00842,   0.     ,   0.14582,  -0.     ,  -0.     ,
          0.05169,   0.     ,   0.00862,  -0.     ,   0.02059],
       [ -0.     ,  -0.     ,  -0.00611,  -0.     ,   0.     ,  -0.     ,   0.23946,   0.03564,
          0.     ,  -0.     ,  -0.     ,   0.0502 ,   0.     ],
       [  0.     ,   0.     ,  -0.00518,  -0.     ,   0.     ,  -0.     ,   0.03564,   0.91335,
         -0.     ,   0.     ,  -0.     ,   0.01797,   0.     ],
       [ -0.00786,   0.01026,   0.     ,  -0.00336,  -0.     ,   0.05169,   0.     ,  -0.     ,
          1.07281,   0.     ,  -0.00273,  -0.     ,   0.00324],
       [  0.     ,   0.     ,   0.     ,  -0.     ,   0.00389,   0.     ,  -0.     ,   0.     ,
          0.     ,   1.0787 ,   0.     ,   0.     ,   0.     ],
       [ -0.02656,  -0.01151,   0.     ,   0.00734,   0.     ,   0.00862,  -0.     ,  -0.     ,
         -0.00273,   0.     ,   1.17073,  -0.     ,  -0.00232],
       [ -0.     ,  -0.     ,   0.02084,   0.     ,   0.     ,  -0.     ,   0.0502 ,   0.01797,
         -0.     ,   0.     ,  -0.     ,   1.25907,  -0.     ],
       [ -0.03727,   0.01092,  -0.     ,  -0.01383,   0.     ,   0.02059,   0.     ,   0.     ,
          0.00324,   0.     ,  -0.00232,  -0.     ,   1.66487]])