8. 非自洽 GGA 密度泛函

这一节我们讨论非自洽 GGA 密度泛函；它也将是 XYG3 型泛函的铺垫。

非自洽密度泛函将泛函 A 的密度代入泛函 B 的能量表达式中，得到体系的能量。

非自洽泛函在密度泛函领域中，并不是主流。在 19 世纪 90 年代，曾经为了互补 HF 自洽场的相关效应与 GGA 的 SIE (Self-Interaction Error) 效应两者的不足，以 HF-DFT 的名义有所发展；但普遍来说，能量或分子结构的表现上，非自洽 GGA 泛函仍然没有显著地优于 GGA 泛函。不仅如此，非自洽泛函的梯度性质的公式形式与自洽泛函的形式有不少区别，增加了一些额外的计算量。

相比于 B2PLYP 泛函是基于自洽的 GGA 密度泛函构建而言，XYG3 泛函正是基于非自洽 GGA 密度泛函构建的。从程序的角度上，由于双杂化本身引入的计算量比较大，因此非自洽泛函本身的额外的计算量相比而言近乎于是可以忽略的了。这在以后的文档中会有实际的体会。

尽管从应用的角度来说，非自洽泛函本身的意义并不大；但如果我们相信密度泛函近似的误差可以分为密度误差与能量误差，那么非自洽泛函可以是密度与能量误差分析的有力工具。譬如说，若我们获得某些分子体系的 Full-CI 密度，并将它代入近似泛函中，就可以对泛函能量表达式的误差进行分析了。

[1]:

from pyscf import scf, gto, dft
import numpy as np
from functools import partial
import warnings

from pyxdh.Utilities.test_molecules import Mol_H2O2
from pyxdh.DerivOnce import GradNCDFT

warnings.filterwarnings("ignore")
np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
np.set_printoptions(5, linewidth=150, suppress=True)

[2]:

mol = Mol_H2O2().mol
grids = Mol_H2O2().gen_grids()

8.1. 量化软件的 HF-B3LYP 计算

8.1.1. PySCF 实现

我们拿 HF-B3LYP 作为非自洽泛函作为范例来考察。HF-B3LYP 的计算过程是，首先计算 HF 自洽场 scf_eng：

[3]:

scf_eng = scf.RHF(mol)
scf_eng.conv_tol = 1e-11
scf_eng.conv_tol_grad = 1e-9
scf_eng.kernel()

[3]:

-150.58503378083847

随后，我们构建用于计算 B3LYP 能量的类 nc_eng；注意 不要执行实际的计算，即不要使用 kernel 或使用 run 成员函数。

[4]:

nc_eng = dft.RKS(mol)
nc_eng.grids = grids

我们将 HF 自洽场 scf_eng 的密度代入 B3LYP nc_eng 的能量中：

\[E_\mathrm{HF-B3LYP} = E_\mathrm{B3LYP} [D_{\mu \nu}^\mathrm{HF}]\]

[5]:

nc_eng.energy_tot(dm=scf_eng.make_rdm1())

[5]:

-150.27716895192074

上述能量便是 HF-B3LYP 的总能量了。比较遗憾的是，几乎没有其它程序提供结果验证 HF-B3LYP 是否正确。

8.1.2. pyxdh 实现

pyxdh 使用 GradNCDFT 或 DipoleNCDFT 实现非自洽泛函能量的计算。我们不妨使用 GradNCDFT：

[6]:

from pyxdh.DerivOnce import GradNCDFT

GradNCDFT 是 GradSCF 的子类，因此它们的初始化、调用方式基本上是相同的；只是在初始化时，需要引入自洽泛函 scf_eng 和未经过 kernel 或 run 方法进行计算的非自洽泛函 nc_eng：

[7]:

config = {
    "scf_eng": scf_eng,
    "nc_eng": nc_eng
}
nch = GradNCDFT(config)

与普通的自洽场一样，通过 eng 属性就可以获得 HF-B3LYP 能量了。

[8]:

nch.eng

[8]:

-150.27716895192074

8.2. 实现参考

8.2.1. 非自洽 Fock 矩阵 \(F_{\mu \nu}^\mathrm{n}\), \(F_{pq}^\mathrm{n}\)

记号说明

从当前部分开始，

上下标 \(\mathrm{s}\) 代表自洽泛函；但通常省略不写；
上下标 \(\mathrm{n}\) 代表非自洽泛函。

在以后，我们经常会需要使用到非自洽泛函所对应的 Fock 矩阵。这里指出，自洽泛函 (对于 HF-B3LYP 而言自洽泛函是 HF) 的 MO 基组 Fock 矩阵 \(F_{pq}\) 应当是对角化的矩阵，且对角元与自洽泛函 (HF) 轨道能满足 \(F_{pp} = \varepsilon_p\) 的关系。对于 GradNCDFT 的实例，调取 MO 基组的 Fock 矩阵的方式与 GradSCF 实例的调取方式相同，即 nch.F_0_mo。我们验证其是否是对称矩阵，且对角元是否就是 HF 轨道能：

[9]:

np.allclose(nch.F_0_mo, np.diag(nch.e))

[9]:

True

但非自洽 Fock 矩阵是通过 B3LYP 泛函生成的 Fock 矩阵，其在 MO 基组下的表示一般地不是对称矩阵。非自洽 Fock 矩阵的 PySCF 生成方式是对非自洽泛函类 nc_eng 代入自洽泛函 scf_eng 的密度矩阵 \(D_{\mu \nu}\) 得到；其中 F_0_ao_nc 表示原子轨道基组的 \(F_{\mu \nu}^\mathrm{n}\)，而 F_0_mo_nc 表示分子轨道基组的 \(F_{pq}^\mathrm{n}\)；

\[F_{pq}^\mathrm{n} = C_{\mu p} F_{\mu \nu}^\mathrm{n} C_{\nu p}\]

[10]:

F_0_ao_nc = nc_eng.get_fock(dm=scf_eng.make_rdm1())
C = scf_eng.mo_coeff
F_0_mo_nc = C.T @ F_0_ao_nc @ C

我们可以简单观察一下非自洽的 MO 基组 Fock 矩阵。首先我们通过 pyxdh 给出占据与非占轨道的分割 so, sv：

[11]:

so, sv = nch.so, nch.sv

占据-占据部分的 \(F_{ij}^\mathrm{n}\) 表示为

[12]:

F_0_mo_nc[so, so]

[12]:

array([[-18.70465,   0.00009,  -0.04849,  -0.06047,  -0.00156,   0.02533,   0.01121,  -0.00991,   0.00898],
       [  0.00009, -18.66858,  -0.04947,   0.05645,   0.03169,   0.00482,   0.01978,  -0.00393,  -0.00309],
       [ -0.04849,  -0.04947,  -1.13497,  -0.0001 ,  -0.02707,  -0.02778,  -0.04806,   0.01615,  -0.01083],
       [ -0.06047,   0.05645,  -0.0001 ,  -0.82973,   0.03003,  -0.02548,   0.00621,   0.00707,  -0.02126],
       [ -0.00156,   0.03169,  -0.02707,   0.03003,  -0.47708,   0.01432,  -0.01259,  -0.00464,   0.00115],
       [  0.02533,   0.00482,  -0.02778,  -0.02548,   0.01432,  -0.45708,   0.02529,   0.01502,   0.0033 ],
       [  0.01121,   0.01978,  -0.04806,   0.00621,  -0.01259,   0.02529,  -0.41389,   0.00197,  -0.03418],
       [ -0.00991,  -0.00393,   0.01615,   0.00707,  -0.00464,   0.01502,   0.00197,  -0.26219,  -0.00388],
       [  0.00898,  -0.00309,  -0.01083,  -0.02126,   0.00115,   0.0033 ,  -0.03418,  -0.00388,  -0.26544]])

我们发现这是一个对称但并不是对角化的矩阵。有意思的是，其对角元与轨道能的值大体相近。我们列举出占据轨道能的值：

[13]:

scf_eng.mo_energy[so]

[13]:

array([-20.67103, -20.63571,  -1.58083,  -1.23933,  -0.76882,  -0.72551,  -0.58961,  -0.52943,  -0.5192 ])

对于非占-非占部分 \(F_{ab}^\mathrm{n}\) 来说也应当类似。而非占-占据部分 \(F_{ai}^\mathrm{n}\) 的值为

[14]:

F_0_mo_nc[sv, so]

[14]:

array([[ 0.03385, -0.01491, -0.00629,  0.01663,  0.00275, -0.00696, -0.00754,  0.00298, -0.00071],
       [ 0.01201,  0.0227 ,  0.00014, -0.00452, -0.00243, -0.00666,  0.00388, -0.00234, -0.00188],
       [ 0.00741,  0.01734, -0.00733,  0.00218, -0.0175 , -0.00175,  0.00629, -0.00331, -0.00352],
       [ 0.01849, -0.01636,  0.00826,  0.00928,  0.00255,  0.00243,  0.00097,  0.00174,  0.0006 ],
       [-0.00415,  0.02197,  0.00169,  0.01539, -0.00569, -0.00636,  0.00983, -0.00407,  0.00402],
       [-0.00323, -0.02753,  0.00781, -0.01028, -0.0118 , -0.00144, -0.00516, -0.00309, -0.0014 ],
       [-0.03135, -0.00558,  0.00415,  0.00697,  0.00632, -0.00457, -0.01114, -0.00306, -0.00469],
       [ 0.03239,  0.00018,  0.00073, -0.00973,  0.0027 , -0.00659,  0.01685,  0.00647,  0.00443],
       [-0.00862, -0.00485, -0.0134 ,  0.00327,  0.00597, -0.00452, -0.01746,  0.01562,  0.00135],
       [ 0.03555,  0.02291, -0.01811,  0.00674, -0.00828,  0.01476, -0.0427 , -0.00309, -0.01121],
       [-0.01729, -0.01734,  0.00822,  0.00745, -0.00033,  0.00609, -0.00038, -0.00131,  0.01795],
       [-0.01999,  0.09043,  0.00079, -0.00644,  0.00902, -0.00174,  0.01827, -0.00173,  0.00276],
       [ 0.0803 ,  0.01993,  0.0007 ,  0.00458,  0.00867,  0.00631,  0.03445, -0.0076 ,  0.01666]])

这部分的值都比较小，但并非是零。我们以后会知道 \(F_{ai}^\mathrm{n}\) 将会对非自洽泛函的梯度性质的计算有至关重要的贡献。

现在我们用 pyxdh 的实例 nch 生成非自洽 Fock 矩阵。nch 本身一般只会调用自洽泛函的结果与性质；但若要调用非自洽泛函，则先使用属性 nc_deriv，随后再调取相应的变量。我们验证使用 pyxdh 调用的 AO 基组非自洽 Fock 矩阵与 PySCF 生成的 F_0_ao_nc 结果相等：

[15]:

np.allclose(nch.nc_deriv.F_0_ao, F_0_ao_nc)

[15]:

True