7. XYG3 的 CP-KS 方程计算相关¶

通过上一节，我们基本了解了 CP-HF 方程的实际代码过程．这一节，我们会引入 DFT 的内容，首先对 B3LYP 的极化率作说明，以了解 CP-KS 方程的具体解法；随后对 XYG3 的极化率作叙述，以了解 XYG3 电场下的梯度计算流程．

7.1. 准备工作¶

7.1.1. 分子计算¶

In [1]:

# 环境搭建
import psi4
import numpy as np

# 引入 DIIS 模块
import sys
sys.path.append("include")
from DIIS_helper import DIIS_helper

# 简化矩阵输出
np.set_printoptions(8, 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,
})

# 波函数信息
scf_e, scf_wfn = psi4.energy("B3LYP", molecule=mol, return_wfn=True)

# 验证结果
psi4.compare_values(scf_e, -76.3771828949, 6, 'B3LYP Energy')

        B3LYP Energy......................................................PASSED

Out[1]:

True

7.1.2. 中间变量¶

In [2]:

# 电子积分引擎
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 基组双电子排斥积分

# DFT 积分引擎
V_pot = scf_wfn.V_potential()  # DFT 积分引擎

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

In [3]:

# 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'))  # 未占轨道能级
dip_psi4 = mints.ao_dipole()
dip = np.array([np.asarray(mat) for mat in dip_psi4])  # AO 基组偶极矩积分
dip_vo = np.einsum("guv, ua, vi -> gai", dip, Cv, Co, optimize=True)  # 分子轨道 ai 的偶极矩积分

In [4]:

# PT2 导出变量
# 轨道能之差的张量 D_{ij}^{ab} -> d[i, a, j, b]
d_ovov = (eo.reshape(-1, 1, 1, 1) - ev.reshape(-1, 1, 1) + eo.reshape(-1, 1) - ev)
# 轨道能之差 D_i^a -> d_vo[a, i]
d_vo = - ev.reshape(-1, 1) + eo
# 全 MO 基组 ERI <pq|rs> = (pr|qs) = g_{pq}^{rs} = g_{rs}^{pq} -> g[p, r, q, s]
g_pqrs = np.einsum("up, vr, uvkl, kq, ls -> prqs", C, C, eri, C, C, optimize=True)
g_ovov = g_pqrs[:nocc, nocc:, :nocc, nocc:]  # <ij|ab> = g_{ij}^{ab}
t_ovov = g_ovov / d_ovov  # 轨道对激发振幅

7.1.3. XYG3 非自洽泛函有关准备¶

In [5]:

# 定义 XYG3 非自洽泛函
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')

    lda_x = psi4.core.LibXCFunctional("XC_LDA_X", restricted)
    lda_x.set_alpha(-0.0140)
    sup.add_x_functional(lda_x)

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

    lyp_c = psi4.core.LibXCFunctional("XC_GGA_C_LYP", restricted)
    lyp_c.set_alpha(0.6789)
    sup.add_c_functional(lyp_c)

    sup.set_x_alpha(0.8033)

    return sup

# 计算 XYG3 能量
nscf_e, nscf_wfn = psi4.energy("SCF", dft_functional=build_xyg3_nc_superfunctional, return_wfn=True)
# 非自洽部分 DFT 积分引擎
Vn_pot = nscf_wfn.V_potential()

In [6]:

# XYG3 非自洽 Fock 矩阵
Vn = psi4.core.Matrix(nbf, nbf)
Vn_pot.set_D([scf_wfn.Da()])
Vn_pot.compute_V([Vn])
Fn = T + V
Fn += 2 * np.einsum("uvkl, kl -> uv", eri, D)
Fn -= 0.8033 * np.einsum("ukvl, kl -> uv", eri, D)
Fn += Vn
Fn_pq = C.T @ Fn @ C

7.1.4. CP-KS 方程变量¶

In [7]:

diis = None
USE_DIIS = True
if USE_DIIS:
    diis = DIIS_helper()
MAX_ITER = 100
CONV = 1.e-8

7.2. B3LYP 极化率¶

B3LYP 所涉及到的 CP-KS 方程与 CP-HF 唯一的不同是在四脚标张量 \(\textbf{A}\) 中存在 DFT 所特有的响应张量 \(\textbf{G}\)；\(\textbf{G}\) 相当于 DFT 势函数的导数，或者说 DFT Kernel 的二阶导数．我们通过构建张量 \(\textbf{A}\) 的过程来了解张量 \(\textbf{G}\) 的实际应用方式．

我们知道，在 HF 方法中，可以使用直接 CP-HF 方法，也可以使用迭代法求解．在 CP-KS 方程中，一般不直接导出 \(\mathbf{G}\) 张量，因此这里始终会使用迭代法求解．迭代求解 B3LYP 极化率所需要的旋转矩阵 \(U_{ai}^g\) 的公式是

\begin{equation} U_{ai}^g = \frac{1}{\varepsilon_i - \varepsilon_a} \left( A_{ij}^{ab} U_{bj}^g - \mu_{ai}^g \right) \end{equation}

其中，\(A_{ij}^{ab}\) 的形式比 Hartree-Fock 下多了响应张量 \(G_{ij}^{ab}\)，并且要在 Exchange 积分贡献上乘上 B3LYP 泛函所指定的 Hartree-Fock 型 Exchange 贡献的系数 \(c_\mathrm{x}\) (对于 B3LYP 来说 \(c_\mathrm{x} = 0.2\))：

\begin{equation} A_{pr}^{qs} = 4 (pq|rs) - c_\mathrm{x} (pr|qs) - c_\mathrm{x} (ps|qr) + 4 G_{pr}^{qs} \end{equation}

在迭代求解的过程中，我们不需要真正地使用四脚标的张量，因此得到下述的表达式：

\begin{equation} A_{ij}^{ab} U_{bj} = C_{\mu i} C_{\nu a} (4 J_{\mu \nu} [D_{\mu \nu}] - c_\mathrm{x} K_{\nu \mu} [D_{\mu \nu}] - c_\mathrm{x} K_{\mu \nu} [D_{\mu \nu}] + 4 G_{\mu \nu}[D_{\mu \nu}]) \end{equation}

其中，上述的广义密度定义为

\begin{equation} D_{\mu \nu} = C_{\mu b} U_{bj}^g C_{\nu j} \end{equation}

上述的广义密度在实际应用中，可以由 psi4.core.JK.D 方法给出；而 \(G_{\mu \nu}[D_{\mu \nu}]\) 则可以由 DFT 积分引擎的 psi4.core.VBase.compute_Vx 在代入广义密度的情况下给出．下面我们就进行 CP-KS 方程的求解．似乎在 Psi4 中，DFT 积分的耗时比较多，所以下述代码的执行需要花一些时间．

In [8]:

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

# 定义左矢为 Cv，右矢为空矩阵
C_right_list = []
for g in range(3):
    jk.C_left_add(psi4.core.Matrix.from_array(Cv))
    mat = psi4.core.Matrix(nbf, nvir)
    C_right_list.append(np.asarray(mat))
    jk.C_right_add(mat)

In [9]:

# 初猜
U = - dip_vo / d_vo
# 旧矩阵
U_old = np.copy(U)
# DIIS 初始化
diis = []
if USE_DIIS:
    for _ in range(3):
        diis.append(DIIS_helper())

In [10]:

%%time
for it in range(1, MAX_ITER + 1):
    # 更新 JK 引擎的右矢为 U_bj C_lj
    for g in range(3):
        C_right_list[g][:] = Co @ U[g].T

    # 计算 J[D]、K[D] 积分
    jk.compute()
    for g in range(3):
        J = np.asarray(jk.J()[g])
        K = np.asarray(jk.K()[g])

        # 构建 DFT 的 G 响应
        Vx = psi4.core.Matrix(nbf, nbf)  # 空矩阵储存 G 响应
        D_jk = jk.D()[g]  # 相当于 Cv @ U[g] @ Co.T
        V_pot.compute_Vx([D_jk], [Vx])  # 计算 G 响应到 Vx 矩阵中
        Vx = np.asarray(Vx)  # 转换为 NumPy 矩阵

        # 可以更新 U 矩阵了
        # Unew_ai^g = [ C_ui C_va (4 * JD_uv - cx KD_vu - cx KD_uv + 4 GD_uv) - mu_ai^g ] / (e_i - e_a)
        U[g] = Cv.T @ (4 * J - 0.2 * K.T - 0.2 * K + 4 * Vx) @ Co - dip_vo[g]
        U[g] /= d_vo

    # 若使用了 DIIS 加速，则执行下述插值
    if USE_DIIS:
        for g in range(3):
            diis[g].add(U[g], U[g] - U_old[g])
            U[g] = diis[g].extrapolate()

    # 检查收敛情况
    rms = np.linalg.norm(U - U_old)
    # print('CPHF Iteration {:3d}: RMS = {:14.10f}'.format(it, rms))

    # 判断是否收敛
    if (rms < CONV):
        print("CPHF Converged in {:3d} iterations!".format(it))
        break
    else:
        U_old = np.copy(U)

CPHF Converged in  12 iterations!
CPU times: user 16.7 s, sys: 504 ms, total: 17.2 s
Wall time: 17.2 s

由此，我们可以给出 B3LYP 的极化率张量：

In [11]:

alpha_b3lyp = np.einsum("gai, fai -> fg", U, dip_vo) * 4
alpha_b3lyp.round(decimals=6)

Out[11]:

array([[ 1.414654, -0.      ,  0.      ],
       [-0.      ,  7.259427, -0.      ],
       [-0.      , -0.      ,  6.452491]])

不过上述的计算结果可能无法很好地与 Guassian 的结果对上，并且差距相当可观．Gaussian 输入卡：B3LYP

In [12]:

alpha_b3lyp_gaussian = np.array([1.4146668, 7.2595695, 6.4526498])
(np.diag(alpha_b3lyp) - alpha_b3lyp_gaussian).round(decimals=7)

Out[12]:

array([-0.0000132, -0.0001426, -0.0001585])

这可能与格点积分的选取有很大的关系．下面是在 Gaussian 中更改格点的精度所产生的 B3LYP 的极化率，不同格点精度的极化率的差基本上与上面的误差的数量级相同．因此，我们认为我们还是基本重复出了 B3LYP 的极化率．

In [13]:

print("Grid UltraFine -      Fine: ",
      np.array([1.4146668, 7.2595695, 6.4526498]) - np.array([1.4146945, 7.2598606, 6.4524672]))
print("Grid SuperFine - UltraFine: ",
      np.array([1.414654, 7.2595522, 6.4525927]) - np.array([1.4146668, 7.2595695, 6.4526498]))

Grid UltraFine -      Fine:  [-0.0000277 -0.0002911  0.0001826]
Grid SuperFine - UltraFine:  [-0.0000128 -0.0000173 -0.0000571]

7.3. XYG3 偶极矩与自然轨道¶

提示

之后会约定，在符号右上标 \(\mathrm{s}\) 的为自洽场泛函 (XYG3 中则表示 B3LYP)，而上标 \(\mathrm{n}\) 的则为非自洽场泛函 (XYG3 中则表示其自身)．

在后文中，会不加说明地代入 XYG3 的 Hartree-Fock 型 Exchange 积分缩放系数 \(c_\mathrm{x}^\mathrm{n} = 0.8033\) 与 B3LYP 型的 \(c_\mathrm{x}^\mathrm{s} = 0.2\)，以及 XYG3 的 MP2 型 Correlation 积分缩放系数 \(c_\mathrm{c}^\mathrm{n} = 0.3211\)．

XYG3 的偶极矩求取方式与 MP2 非常类似，区别在于 CP-HF 方程更变为 CP-KS 方程，以及 MP2 拉格朗日矩阵更变为 XYG3 拉格朗日矩阵．由于 CP-KS 方程的导出前提是 \(\partial_\xi F_{pq}^\mathrm{s} = 0\)，因此 XYG3 的 CP-KS 方程形式没有变化；但拉格朗日矩阵的导出则是从能量的梯度所产生的，因此其形式与 MP2 会有所不同．

7.3.1. XYG3 拉格朗日矩阵的构建¶

首先，我们写出 XYG3 拉格朗日量：

\begin{equation} L_{ai}^\mathrm{n} = \frac{1}{2} P_{jk}^{(2)} A^\mathrm{s,}{}_{aj}^{ik} + \frac{1}{2} P_{bc}^{(2)} A^\mathrm{s,}{}_{ab}^{ic} - c_\mathrm{c}^\mathrm{n} T_{jk}^{ab} g_{ib}^{jk} + c_\mathrm{c}^\mathrm{n} T_{ij}^{bc} g_{aj}^{bc} + F_{ai}^\mathrm{n} \end{equation}

为了程序书写上的简便，并且区分 DFT 和与 MP2 拉格朗日矩阵类似的项区分开，我们定义张量

\begin{equation} a_{pr}^{qs} = 4 (pq|rs) - c_\mathrm{x} (pr|qs) - c_\mathrm{x} (ps|qr) = A^\mathrm{s,}{}_{pr}^{qs} - 4 G^\mathrm{s,}{}_{pr}^{qs} \end{equation}

同时定义弛豫密度的部分贡献

\begin{equation} p_{pq}^{(2)} = \left\{ \begin{matrix} P_{ij}^{(2)} = - c_\mathrm{c}^\mathrm{n} T_{ik}^{ab} t_{jk}^{ab}, & p = i \text{ and } q = j \\ P_{ab}^{(2)} = c_\mathrm{c}^\mathrm{n} T_{ij}^{ac} t_{ij}^{bc}, & p = a \text{ and } q = b \\ 0, & \text{otherwise} \end{matrix} \right. \end{equation}

那么我们可以重新写拉格朗日量为

\begin{equation} L_{ai}^\mathrm{n} = \frac{1}{2} p_{pq}^{(2)} a_{ap}^{iq} - c_\mathrm{c}^\mathrm{n} T_{jk}^{ab} g_{ib}^{jk} + c_\mathrm{c}^\mathrm{n} T_{ij}^{bc} g_{aj}^{bc} + 2 G_{ai}^\mathrm{s} [p_{pq}^{(2)}] + F_{ai}^\mathrm{n} \end{equation}

我们首先构造张量 \(a_{pr}^{qs}\) 与弛豫密度部分贡献 \(p_{pq}^{(2)}\)：

In [14]:

T_ovov = 2 * t_ovov - t_ovov.swapaxes(1,3)  # 闭壳层的 T 张量
P2_oo = - 0.3211 * np.einsum("iakb, jakb -> ij", T_ovov, t_ovov)  # 占据-占据弛豫密度
P2_vv = 0.3211 * np.einsum("iajc, ibjc -> ab", T_ovov, t_ovov)  # 非占-非占弛豫密度
p2_pq = np.block([
    [P2_oo, np.zeros((nocc, nvir))],
    [np.zeros((nvir, nocc)), P2_vv],
])  # 弛豫密度部分贡献
# A 张量的非 DFT 贡献
a = 4 * g_pqrs - 0.2 * np.einsum("prqs -> pqrs", g_pqrs) - 0.2 * np.einsum("psqr -> pqrs", g_pqrs)

随后就可以构建 \(L_{ai}^\mathrm{n}\) 了：

In [15]:

# 第一段：非 DFT 贡献
L_vo = 0.5 * np.einsum("pq, aipq -> ai", p2_pq, a[nocc:, :nocc, :, :])
L_vo += - 0.3211 * np.einsum("jakb, ijbk -> ai", T_ovov, g_pqrs[:nocc, :nocc, nocc:, :nocc])
L_vo += 0.3211 * np.einsum("ibjc, abjc -> ai", T_ovov, g_pqrs[nocc:, nocc:, :nocc, nocc:])
# 第二段：DFT 贡献
Gs_vo = psi4.core.Matrix(nbf, nbf)
V_pot.compute_Vx([psi4.core.Matrix.from_array(C @ p2_pq @ C.T)], [Gs_vo])
L_vo += 2 * Cv.T @ Gs_vo @ Co
# 第三段：非自洽 Fock 矩阵
L_vo += Fn_pq[nocc:, :nocc]

7.3.2. XYG3 CP-KS 方程¶

随后解下述 CP-KS 方程即可：

\begin{equation} P_{ai}^{(2)} = \frac{1}{\varepsilon_i - \varepsilon_a} \big( A^\mathrm{s,}{}_{ij}^{ab} P_{bj}^{(2)} + L_{ai} \big) \end{equation}

In [16]:

# JK 引擎
jk = psi4.core.JK.build(scf_wfn.basisset())
jk.initialize()
jk.C_left_add(psi4.core.Matrix.from_array(Cv))
mat = psi4.core.Matrix(nbf, nvir)
C_right = np.asarray(mat)
jk.C_right_add(mat)

In [17]:

# 初猜与收敛设定
P2_vo = L_vo / d_vo
P2_vo_old = np.copy(P2_vo)
diis = None
if USE_DIIS:
    diis = DIIS_helper()

In [18]:

%%time
# CP-HF 方程
for it in range(1, MAX_ITER + 1):
    C_right[:] = Co @ P2_vo.T
    jk.compute()
    J = np.asarray(jk.J()[0])
    K = np.asarray(jk.K()[0])
    Vx = psi4.core.Matrix(nbf, nbf)
    D_jk = jk.D()[0]
    V_pot.compute_Vx([D_jk], [Vx])
    Vx = np.asarray(Vx)
    P2_vo = Cv.T @ (4 * J - 0.2 * K.T - 0.2 * K + 4 * Vx) @ Co + L_vo
    P2_vo /= d_vo

    # 若使用了 DIIS 加速，则执行下述插值
    if USE_DIIS:
        diis.add(P2_vo, P2_vo - P2_vo_old)
        P2_vo = diis.extrapolate()

    # 检查收敛情况
    rms = np.linalg.norm(P2_vo - P2_vo_old)
    # print('CPHF Iteration {:3d}: RMS = {:14.10f}'.format(it, rms))

    # 判断是否收敛
    if (rms < CONV):
        print("CPHF Converged in {:3d} iterations!".format(it))
        break
    else:
        P2_vo_old = np.copy(P2_vo)

CPHF Converged in   7 iterations!
CPU times: user 3.2 s, sys: 104 ms, total: 3.31 s
Wall time: 3.3 s

最后重构弛豫密度矩阵

\begin{equation} P_{pq}^{(2)} = \begin{pmatrix} P_{ij}^{(2)} & P_{ia}^{(2)} \\ P_{ai}^{(2)} & P_{ab}^{(2)} \end{pmatrix} \end{equation}

In [19]:

P2_pq = np.block([
    [P2_oo, P2_vo.T],
    [P2_vo, P2_vv]
])

In [20]:

P2 = C @ P2_pq @ C.T

7.3.3. XYG3 自然轨道占据数与偶极矩¶

有了弛豫密度后，许多自然轨道占据数与偶极矩就非常容易给出．自然轨道表示如下：

In [21]:

# 定义自洽场密度
D_pq = np.zeros_like(D)
for i in range(5):
    D_pq[i, i] = 1
# 对自洽场密度与弛豫密度的和进行对角化
no_occ, no_coef = np.linalg.eigh(P2_pq + D_pq)
# 从而得到自然轨道的占据数
no_occ[::-1] * 2

Out[21]:

array([1.9999934 , 1.99409371, 1.98761029, 1.9803091 , 1.97868861, 0.02074072, 0.01844683,
       0.0116838 , 0.00571251, 0.00147259, 0.00073034, 0.00028439, 0.00023373])

B3LYP 贡献的电子云偶极矩与原子核偶极矩表示如下．Gaussian 输入卡：B3LYP

In [22]:

# 原子核偶极矩
mol_geo = np.asarray(mol.geometry())
neu_charge = []
for neu in range(mol_geo.shape[0]):
    neu_charge.append(mol.charge(neu))
dip_mol = np.zeros(3)
for neu in range(mol_geo.shape[0]):
    for ind in range(3):
        dip_mol[ind] += neu_charge[neu] * mol_geo[neu][ind]

# B3LYP 偶极矩
dip_b3lyp = 2 * np.einsum("fuv, uv -> f", dip, D)
dip_b3lyp = dip_b3lyp.round(decimals=6)
# 验证结果
psi4.compare_values((dip_mol + dip_b3lyp)[2], 1.031112, 6, 'B3LYP Energy')

        B3LYP Energy......................................................PASSED

Out[22]:

True

XYG3 泛函贡献的偶极矩部分则可以写为弛豫部分与 B3LYP 贡献部分的加和：

In [23]:

dip_xyg3 = 2 * np.einsum("fuv, uv -> f", dip, P2)
dip_xyg3 = dip_xyg3.round(decimals=6)
dip_tot_xyg3 = dip_mol + dip_b3lyp + dip_xyg3
dip_tot_xyg3

Out[23]:

array([0.        , 0.        , 1.07524207])