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\) 的公式是
其中,\(A_{ij}^{ab}\) 的形式比 Hartree-Fock 下多了响应张量 \(G_{ij}^{ab}\),并且要在 Exchange 积分贡献上乘上 B3LYP 泛函所指定的 Hartree-Fock 型 Exchange 贡献的系数 \(c_\mathrm{x}\) (对于 B3LYP 来说 \(c_\mathrm{x} = 0.2\)):
在迭代求解的过程中,我们不需要真正地使用四脚标的张量,因此得到下述的表达式:
其中,上述的广义密度定义为
上述的广义密度在实际应用中,可以由 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 拉格朗日量:
为了程序书写上的简便,并且区分 DFT 和与 MP2 拉格朗日矩阵类似的项区分开,我们定义张量
同时定义弛豫密度的部分贡献
那么我们可以重新写拉格朗日量为
我们首先构造张量 \(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 方程即可:
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
最后重构弛豫密度矩阵
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])