3.6. GGA 自洽场计算¶
上一节我们已经讨论了 LDA 的自洽场计算与交换或相关能的计算。但一者，LDA 由于其精度不足以解决普遍的化学问题，因此并不是当前主流的泛函；二者，XYG3 等 xDH 型双杂化泛函所使用的基础泛函来自于 GGA；三者，在 PySCF 中 LDA 与 GGA 的实现方式比较不同。
这一节我们讨论 GGA 的自洽场计算，这会是理解后续文档的重要基础。这里所指的 GGA 允许杂化，但从程序实现的角度上，不允许 LDA 与 meta-GGA。尽管说 GGA 原则上是兼容 LDA 的，但 GGA 因为其泛函核形式中需要 \(\partial_r \rho = \rho_r\)，从而有许多 LDA 泛函中没有的表达式；这在以后推导梯度时会有实际体会。因此，从程序的角度上，GGA 退化到 LDA 尽管程序框架没有变化，但代码细节上的变动将会非常巨大。我们以后也仅仅关注 GGA 的代码实现。
GGA 的种类繁多，我们在这里就以 B3LYP 举例。
[1]:
import numpy as np
import warnings
from pyscf import scf, gto, dft
from pkg_resources import resource_filename
from pyxdh.Utilities import FormchkInterface, GridHelper, KernelHelper, GridIterator
from pyxdh.Utilities.test_molecules import Mol_H2O2
from pyxdh.DerivOnce import GradSCF
from functools import partial
np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
warnings.filterwarnings("ignore")
np.set_printoptions(5, linewidth=150, suppress=True)
dft.numint.libxc = dft.xcfun
ni = dft.numint.NumInt()
3.6.1. 量化软件的 GGA 计算¶
3.6.1.1. PySCF 计算¶
我们仍然拿双氧水分子进行计算，格点取 (99, 590)。关于分子、格点初始化的说明，以及如何用 PySCF 构建分子 mol，见 RHF 自洽场文档的说明；关于如何用 PySCF 构建格点 grids，见 格点生成部分。
[2]:
mol = Mol_H2O2().mol
grids = Mol_H2O2().gen_grids()
nao = mol.nao
ngrids = grids.weights.size
那么 PySCF 对 B3LYP 的计算可以通过下述的代码实现：
[3]:
scf_eng = dft.RKS(mol)
scf_eng.xc = "B3LYPg"
scf_eng.grids = grids
scf_eng.kernel()
[3]:
-151.3775435605392
需要注意，PySCF 的泛函命名规则遵循大多数除 Gaussian 以外的其它量化软件的规则，因此默认的 B3LYP 关键词下，其中具有的 LDA 相关能贡献是 VWN5 泛函；但 Gaussian 采用的 VWN3。为了取得与 Gaussian 计算较为相近的结果，这里的泛函名称使用 B3LYPg，以使其中的 LDA 相关泛函为 VWN3。
3.6.1.2. Gaussian 计算¶
Gaussian 的计算结果储存在 pyxdh 库的资源文件夹下，输入卡文件是 H2O2-B3LYP-freq.gjf，输出的 formchk 文件是 H2O2-B3LYP-freq.fchk。输入卡内容是
%chk=H2O2-B3LYP-freq
#p RB3LYP/6-31G nosymm SCF(VeryTight) Freq Int(Grid=99590)
H2O2 B3LYP Frequency
0 1
O 0.0 0.0 0.0
O 0.0 0.0 1.5
H 1.0 0.0 0.0
H 0.0 0.7 1.0
[4]:
ref_fchk = FormchkInterface(resource_filename("pyxdh", "Validation/gaussian/H2O2-B3LYP-freq.fchk"))
Gaussian 计算给出的 B3LYP 自洽场能量用下述方法调取：
[5]:
ref_fchk.total_energy()
[5]:
-151.3775436089372
尽管我们已经使 PySCF 的自洽场计算配置尽可能接近 Gaussian，但结果上仍然有一些差别；尽管在能量上看不出来，但在求 B3LYP 的梯度时，这种差别就会变得明显。
3.6.1.3. pyxdh 计算¶
pyxdh 对 B3LYP 的计算与 RHF 的计算方式是相同的：
[6]:
config = {"scf_eng": scf_eng}
scfh = GradSCF(config)
scfh.eng
[6]:
-151.3775435605392
绝大多数情况下，GGA 与 HF 在 pyxdh 下可以等同地处理。
3.6.2. 格点与泛函相关定义¶
我们单列一小节对格点与泛函相关的内容作定义。尽管在上一节中，我们已经对 LDA 型泛函的格点与泛函已经有所了解，但为了以后文档的统一，在这里我们系统地介绍格点与泛函记号，并且介绍 pyxdh 对格点与泛函的处理程序 GridHelper、KernelHelper。这一小节中有不少变量未必需要用在 GGA 的能量计算过程中，但为了避免以后零散的说明，我们仍然对它们进行讨论。
格点积分有关的量包括格点本身的性质、轨道或密度及其梯度格点，以及泛函核格点．轨道或密度在原子坐标下的梯度我们会在将来叙述；这一小节的梯度指的是电子坐标的导数．
3.6.2.1. 泛函核无关部分¶
记号说明
- \(\rho\) 代表电子态密度密度
- \(\rho_r = \partial_r \rho\)
- \(\rho_{rw} = \partial_r \partial_w \rho\)
- \(\gamma = \rho_r \rho_r\) 表示密度梯度量
- 轨道与轨道的电子坐标导数记号类同，其使用范例参见 动能积分。
注意
以后在公式中，将不会再出现格点记号 \(g\)；但程序中仍然需要显式地考虑格点的指标。一般来说，
- 除了轨道格点外的所有格点张量，格点指标总是在最后一维度；譬如 \(\rho_{rw} \rightarrow \rho_{rwg}\)；
- 包含 AO 轨道指标的张量，AO 轨道指标在最后维度，而格点指标在 AO 轨道指标之前；譬如 \(\phi_{rw \mu} \rightarrow \phi_{rw g \mu}\)。
3.6.2.1.1. GridHelper 生成格点¶
GridHelper 可以一次性生成 xDH 型泛函二阶梯度所需要的轨道与密度梯度格点。其中的电子坐标梯度格点和其它标量有：
ngrid格点数量weight格点权重ao各阶电子坐标偏导数的 AO 轨道格点ao_0轨道格点 \(\phi_\mu\)ao_1轨道格点一阶导数 \(\phi_{r \mu} = \partial_r \phi_\mu\)ao_2轨道格点二阶导数 \(\phi_{r w \mu} = \partial_r \partial_w \phi_\mu\)ao_3轨道格点三阶导数 \(\phi_{r w t \mu} = \partial_r \partial_w \partial_t \phi_\mu\)ao_2T轨道格点二阶导数，但两个坐标分量打包在一个维度中 \(\phi_{T \mu} = \partial_{T_1} \partial_{T_2} \phi_\mu\)ao_3T轨道格点三阶导数，但其中两个坐标分量打包在一个维度中 \(\phi_{T r \mu} = \partial_{T_1} \partial_{T_2} \partial_r \phi_\mu\)rho_0密度格点 \(\rho = D_{\mu \nu} \phi_\mu \phi_\nu\)rho_1密度格点一阶导数 \(\rho_r = \partial_r \rho = 2 D_{\mu \nu} \phi_{r \mu} \phi_\nu\)rho_2密度格点二阶导数 \(\rho_{rw} = \partial_r \partial_w \rho = 2 D_{\mu \nu} (\phi_{r w \mu} \phi_\nu + \phi_{r \mu} \phi_{w \nu})\)rho_01密度格点与其一阶导数的合并张量；只用于生成泛函核导数
GridHelper 的初始化可以通过下述语句实现，需要代入的量是分子结构、格点、以及电子态密度。下述 grdh 是 GridHelper 的一个实例：
[7]:
grdh = GridHelper(mol, grids, scfh.D)
若要调用上述的轨道格点，譬如要查看轨道格点的三阶导数 \(\phi_{rwt \mu}\) 的维度信息，只要用下述代码即可：
[8]:
grdh.ao_3.shape
[8]:
(3, 3, 3, 130776, 22)
注意到 3 代表空间中的三个坐标分量 \((x, y, z)\)，130,776 代表格点数量，22 代表 AO 基组数量。
任务 (1)
相比于以前我公式的推导，这里公式的记号因为去除了格点记号，因此不能单纯地从角标查看张量维度的信息了．请自行查看上述各个张量的维度信息。对于张量维度的把握是正确处理量化公式与 numpy 程序实现的第一步，也是至关重要的一步。
(演示) 我们在上一节的 轨道格点 一段中已经尝试生成轨道格点。请尝试仿照上一节的文档，生成上述格点，并用
np.allclose验证结果。提示：
ni.block_loop(mol, grids, nao, 1, g_mem)的1只会生成 AO 轨道格点的一阶电子坐标导数；试将1替换为3，就可以生成三阶导数了。生成密度格点一阶电子坐标导数的另一个看起来更合理的做法是 \(\rho_r = D_{\mu \nu} (\phi_{r \mu} \phi_\nu + \phi_\mu \phi_{r \nu})\)。试问为何 \(\rho_r = 2 D_{\mu \nu} \phi_{r \mu} \phi_\nu\) 也是正确的？
这是一个非常关键的问题。很多时候我们需要利用 \(\mu, \nu\) 的对称性，但不是所有 \(\mu, \nu\) 都具有对称性，也不是所有具有对称性的表达式都可以被简化。对这类问题的理解会极大帮助我们正确地推导公式并作公式与代码的简化。
3.6.2.1.2. GridIterator 生成格点¶
GridHelper 作为示范范例会在文档中经常使用。但在 上一节的习题 中，我们知道一般来说，对于稍大一些的分子或稍高一些的格点精度，一般电脑的内存就无法装下所有的 DFT 格点。因此，如果在计算实际体系时使用一次性生成所有格点的 GridHelper，很容易出现内存不够的情况。为此，分批地生成格点是非常有必要的。分批生成格点的好处在于缓解内存压力，但缺点是每次生成的格点都需要重新计算一次，即以计算时间作为代价换取内存空间。
在 pyxdh 中，分批生成格点的类为 GridIterator；这个类经常用于写实际计算用的 pyxdh 程序。其程序的使用方式非常接近 GridHelper，只是需要作为类迭代器来使用。
作为例子，我们看利用 GridIterator 生成 \(\rho_r\) 的过程。我们定义 GridIterator 的实例是 grdit 并分配 50 MB 内存用于在每个 batch 中生成轨道格点；GridIterator 生成的 \(\rho_r\) 的变量名为 rho_1。
[9]:
grdit = GridIterator(mol, grids, scfh.D, deriv=3, memory=100)
rho_1 = np.zeros((3, grids.weights.size))
g_start, g_end, g_count = 0, 0, 0
for it in grdit:
g_end += grdit.ngrid
rho_1[:, g_start:g_end] = it.rho_1
g_start = g_end
g_count += 1
print(np.allclose(rho_1, grdh.rho_1))
print("Iteration times:", g_count)
True
Iteration times: 10
上面代码中，关键的两行分别是第 4 行的迭代过程，和第 6 行调用 it.rho_1。变量 it 在迭代过程中，从调用方式上可以几乎与 GridHelper 实例相同。
由于 GridIterator 调用上多少有些方便，加之双氧水分子体系也足够小，因此在文档中通常不用 GridIterator。
3.6.2.2. 泛函核有关部分¶
记号说明
- \(f\) 代表泛函核；泛函核满足关系：在函数图景下 \(E_\mathrm{xc} = \int f[\rho] \rho(\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r}\)，或格点积分下，\(E_\mathrm{xc} = f \rho\)
- \(f_\rho = \partial_\rho (f \rho)\)；注意不是 \(\partial_\rho f\)，这种记号可能引起歧义但足够简洁
- \(f_\gamma = \partial_\gamma (f \rho)\)
- \(f_{\rho \gamma} = \partial_\rho \partial_\gamma (f \rho)\)，其它高阶导数同理
- \(c_\mathrm{x}\) 代表杂化泛函中的精确交换积分贡献．
注意
所有 DFT 格点权重将会归并到泛函核向量中。举例来说，使用前一节的带格点 \(g\) 与权重向量 \(w_g\) 的两个表达式
在以后的文档中则将会写为
KernelHelper 类通过读入必要的密度格点信息，返回泛函核及其相对于电子密度、密度梯度量 \(\rho, \gamma\) 的格点向量。该类有以下成员变量：
c_x泛函杂化系数 \(c_\mathrm{x}\)exc带权重的 \(f\)fr带权重的 \(f_\rho\)fg带权重的 \(f_\gamma\)frr带权重的 \(f_{\rho \rho}\)frg带权重的 \(f_{\rho \gamma}\)fgg带权重的 \(f_{\gamma \gamma}\)frrr带权重的 \(f_{\rho \rho \rho}\)frrg带权重的 \(f_{\rho \rho \gamma}\)frgg带权重的 \(f_{\rho \gamma \gamma}\)fggg带权重的 \(f_{\gamma \gamma \gamma}\)
KernelHelper 的初始化可以通过下述语句实现，需要代入的量是 GridHelper 或 GridIterator 的实例、泛函名称以及求导级别。下述 kerh 是 KernelHelper 的一个实例：
[10]:
kerh = KernelHelper(grdh, "B3LYPg", deriv=3)
若要调用上述的泛函核或其导数格点，譬如要查看带权重的 \(f_{\rho \gamma}\) 的维度信息，只要用下述代码即可：
[11]:
kerh.frg.shape
[11]:
(130776,)
需要注意所有泛函核格点都是一维向量；尽管看上去 \(f_{\rho \gamma}\) 有两个角标，连同格点指标一起似乎应该组成三维张量。这是因为 \(\rho, \gamma\) 对于当前体系来说是独一无无二的；不像 \(\phi_{rw \mu}\) 中，\(r, w\) 分别可能有三种取值，而 \(\mu\) 在当前的 6-31G 双氧水下有 22 中取值，因此 \(\phi_{rw \mu}\) 连同格点指标应是四维张量。
任务 (2)
(演示) 我们在上一节的 Slater 交换能计算 一段中已经尝试生成轨道格点。请尝试仿照上一节的文档，生成上述带权重泛函核格点，并用
np.allclose验证结果。提示：尝试使用
ni.eval_xc("B3LYPg", grdh.rho_01, deriv=3)，查看器返回值。(可选)
KernelHelper在初始化时，不只可以代入GridHelper的实例，还可以代入GridIterator的实例。试用GridIterator实例生成完整的带权重的泛函核格点。
3.6.3. GGA 自洽场实现参考¶
前三个问题是电子数、交换相关能与交换相关势。这些可以通过函数 ni.nr_rks 生成。
[12]:
ni.nr_rks.__func__
[12]:
<function pyscf.dft.numint.nr_rks(ni, mol, grids, xc_code, dms, relativity=0, hermi=0, max_memory=2000, verbose=None)>
[13]:
xc_n, xc_e, xc_v = ni.nr_rks(mol, grids, "B3LYPg", scfh.D)
3.6.3.1. 电子数 \(n_\mathrm{nelec}\)¶
电子数其实通过分子本身就已经被定义：
[14]:
mol.nelectron
[14]:
18
但我们可以借下述格点积分验证生成的密度是否合理：
[15]:
(grids.weights * grdh.rho_0).sum()
[15]:
18.000000186404467
我们可以验证上述结果与 PySCF 的结果相等：
[16]:
np.allclose((grids.weights * grdh.rho_0).sum(), xc_n)
[16]:
True
3.6.3.2. 交换相关能 \(E_\mathrm{xc}\)¶
[17]:
(kerh.exc * grdh.rho_0).sum()
[17]:
-14.506876761363431
我们可以验证上述结果与 PySCF 的结果相等：
[18]:
np.allclose((kerh.exc * grdh.rho_0).sum(), xc_e)
[18]:
True
之所以这里的代码不需要乘上 grids.weights，但计算电子数需要乘上 grids.weights，是因为 KernelHelper 已经将格点权重乘在泛函核的格点中了。在以后的文档中，一般不会出现需要显式地乘以 grids.weights 的情况。
3.6.3.3. 交换相关势 \(v_{\mu \nu}^\mathrm{xc} [\rho]\)¶
[19]:
np.allclose(
+ np.einsum("g, gu, gv -> uv", kerh.fr, grdh.ao_0, grdh.ao_0)
+ 2 * np.einsum("g, rg, rgu, gv -> uv", kerh.fg, grdh.rho_1, grdh.ao_1, grdh.ao_0)
+ 2 * np.einsum("g, rg, gu, rgv -> uv", kerh.fg, grdh.rho_1, grdh.ao_0, grdh.ao_1),
xc_v
)
[19]:
True
任务 (3)
以前我们生成密度的一阶梯度 \(\rho_r\) 时提到，那时的公式与代码中的 \(\phi_{r \mu} \phi_\nu + \phi_\mu \phi_{r \nu}\) 可以简化为 \(2 \phi_{r \mu} \phi_\nu\) 进行计算。试问现在生成交换相关势时，是否也可以这么简化？为什么？
(可选) 你可能已经理解不可以像生成 \(\rho_r\) 时那样简化 \(\phi_{r \mu} \phi_\nu + \phi_\mu \phi_{r \nu}\) 了，但你仍然可以依靠 \(\mu, \nu\) 角标的对称性质，对上面代码块的计算耗时优化到原先的 2/3 倍。请提出你的解决方案。
依靠 \(\mu, \nu\) 角标的对称性质将会在以后经常使用；这不仅会提高代码效率，同时也简化公式推导过程。
(可选) 我们没有给出交换相关势的推导过程．请尝试推导交换相关势，并将你的推导与上面的公式对应，以熟悉这份笔记的记号体系。
既然 \(v_{\mu \nu}^\mathrm{xc} [\rho]\) 是与密度有关的量，那么其构成中，哪些张量会具体地因与体系密度不同而变化，而那些则始终不变？你是否认为用 \(D_{\kappa \lambda}\) 替换掉方括号中的 \(\rho\) 是合理的行为？
3.6.3.4. 重叠积分 \(S_{\mu \nu}\)¶
DFT 格点积分的一个比较有意思的应用方式是求电子积分。但需要指出，DFT 格点积分远不如 Gaussian 基组的解析积分来得快，精度上也有瑕疵。
[20]:
np.allclose(
mol.intor("int1e_ovlp"),
np.einsum("g, gu, gv -> uv", grdh.weight, grdh.ao_0, grdh.ao_0),
atol=1e-7
)
[20]:
True
[21]:
%%timeit -r 7 -n 7
mol.intor("int1e_ovlp")
739 µs ± 380 µs per loop (mean ± std. dev. of 7 runs, 7 loops each)
[22]:
%%timeit -r 7 -n 7
np.einsum("g, gu, gv -> uv", grdh.weight, grdh.ao_0, grdh.ao_0)
19.1 ms ± 1.34 ms per loop (mean ± std. dev. of 7 runs, 7 loops each)
3.6.3.5. 动能积分 \(T_{\mu \nu}\)¶
[23]:
np.allclose(
mol.intor("int1e_kin"),
- 0.5 * np.einsum("g, gu, gvr -> uv", grdh.weight, grdh.ao_0, grdh.ao_2.diagonal(axis1=0, axis2=1)),
atol=1e-6
)
[23]:
True
[24]:
%%timeit -r 7 -n 7
mol.intor("int1e_kin")
786 µs ± 415 µs per loop (mean ± std. dev. of 7 runs, 7 loops each)
[25]:
%%timeit -r 7 -n 7
- 0.5 * np.einsum("g, gu, gvr -> uv", grdh.weight, grdh.ao_0, grdh.ao_2.diagonal(axis1=0, axis2=1))
28.2 ms ± 894 µs per loop (mean ± std. dev. of 7 runs, 7 loops each)
3.6.3.6. Fock 矩阵 (GGA) \(F_{\mu \nu} [R_{\kappa \lambda}]\)¶
手动实现
回顾 RHF 的 Fock 矩阵：
对于 GGA 而言，其 Hamiltonian Core、库伦积分仍然相同，但交换积分前需要乘以杂化系数 \(c_\mathrm{x}\)，并且要加上交换相关势。因此，
根据 B3LYP 的原始文献 [Bec93] 的式 (3)，B3LYP 的 Exact 交换系数应当是 \(c_\mathrm{x} = 0.2\)。我们可以通过下述方式给出交换系数 cx：
[26]:
cx = ni.hybrid_coeff("B3LYPg")
cx
[26]:
0.2
[27]:
np.random.seed(1)
R = np.random.random((nao, nao))
R += R.T
以此可以给出 F_0_ao_R \(F_{\mu \nu} [R_{\kappa \lambda}]\)：
[28]:
F_0_ao_R = scfh.H_0_ao + scf_eng.get_j(dm=R) - 0.5 * cx * scf_eng.get_k(dm=R) + ni.nr_rks(mol, grids, "B3LYPg", R)[2]
PySCF 实现
与 RHF 相同地，可以使用 get_fock 成员函数给出任意密度下的 Fock 矩阵：
[29]:
np.allclose(scf_eng.get_fock(dm=R), F_0_ao_R)
[29]:
True
pyxdh 实现
pyxdh 的实例 scfh 不提供任意密度下的 Fock 矩阵求取方法；因此只能给出自洽场密度下的 Fock 矩阵 \(F_{\mu \nu} [D_{\kappa \lambda}]\)：
[30]:
np.allclose(scfh.F_0_ao, scf_eng.get_fock())
[30]:
True
3.6.3.7. 电子态能量 (GGA) \(E_\mathrm{elec}[X_{\mu \nu}]\)¶
回顾 RHF 的 电子态能量：
类似于 Fock 矩阵，GGA 的电子态能量的形式是
其中 \(f^R\), \(\rho^R\) 表示这些密度格点是由 \(R_{\mu \nu}\) 所产生的。
PySCF 实现
[31]:
eng_elec_R = scf_eng.energy_elec(dm=R)[0]
eng_elec_R
[31]:
123.10771512394984
pyxdh 实现
pyxdh 的实例 scfh 也不提供任意密度下的电子态能量，只能给出自洽场密度下的总能量；若要获取电子态能量，需要减去核与核排斥能：
[32]:
np.allclose(scfh.eng - scf_eng.energy_nuc(), scf_eng.energy_elec()[0])
[32]:
True
任务 (4)
(可选) 请尝试不依靠 PySCF 的成员函数
get_j、get_k、ni_rks以及分子轨道系数 \(C_{\mu p}\)，但可以使用mol的成员函数、ni.block_loop,ni.eval_xc或GridHelper,KernelHelper、交换系数 \(c_\mathrm{x} = 0.2\)、以及 numpy，给出任意密度下双氧水分子 B3LYP 的体系能量 \(E_\mathrm{elec} [R_{\mu \nu}]\)，并与eng_elec_R的结果作对比。杂化 GGA 泛函的电子态能量，除去 \(h_{\mu \nu}\) 贡献项之外，是否可以由 \(\rho\) 或 \(D_{\mu \nu}\) 确定而不依赖 \(C_{\mu p}\)？
既然交换相关势 \(v_{\mu \nu}^\mathrm{xc} [\rho]\) 是由交换相关能 \(E_\mathrm{xc} [\rho]\) 对密度的变分 \(\rho\) 而来：
\[v_{\mu \nu}^\mathrm{xc} [\rho] = \langle \mu | \frac{\delta E_\mathrm{xc} [\rho]}{\delta \rho} | \nu \rangle\]那么交换相关能是否可以从交换相关势与电子密度的乘积得来？
\[\begin{split}\begin{align} E_\mathrm{xc} [\rho] \stackrel{?}{=} v_{\mu \nu}^\mathrm{xc} [\rho] D_{\mu \nu} &= \delta_{ij} \langle i | \frac{\delta E_\mathrm{xc} [\rho]}{\delta \rho} | j \rangle \\ &= \int \frac{\delta E_\mathrm{xc} [\rho]}{\delta \rho} \rho (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} \\ \end{align}\end{split}\]
注意
我们在 RHF 与 GGA 自洽场两节中都提到，pyxdh 的实例无法提供任意密度下的电子态能量或 Fock 矩阵，或者说，在泛函 A 的能量表达式中代入 任意 密度。但是 pyxdh 同时也可以计算 XYG3 型泛函的能量；而 XYG3 型泛函是非自洽型泛函，需要在泛函 A 的能量表达式中代入 指定的 泛函 B 的自洽场密度。
这一定程度上是处于代码安全的角度出发而产生的结果。关于上述看似矛盾的功能描述，请参考后一节关于 XYG3 型泛函的文档。
3.6.4. 参考任务解答¶
3.6.4.1. 任务 (1)¶
3.6.4.1.1. 任务 (1.1)¶
[33]:
print("{:8s}{}".format("Name", "Shape"))
print("{:8s}{}".format("------", "------"))
for name in ["weight", "ao", "ao_0", "ao_1", "ao_2", "ao_3", "ao_2T", "ao_3T", "rho_0", "rho_1", "rho_2", "rho_01"]:
print("{:8s}{}".format(name, grdh.__dict__[name].shape))
Name Shape
------ ------
weight (130776,)
ao (20, 130776, 22)
ao_0 (130776, 22)
ao_1 (3, 130776, 22)
ao_2 (3, 3, 130776, 22)
ao_3 (3, 3, 3, 130776, 22)
ao_2T (6, 130776, 22)
ao_3T (3, 6, 130776, 22)
rho_0 (130776,)
rho_1 (3, 130776)
rho_2 (3, 3, 130776)
rho_01 (4, 130776)
3.6.4.1.2. 任务 (1.2) 演示¶
下面我们生成张量，并逐一与 GridHelper 的实例 grdh 的成员变量作比较。
格点数量 ngrid
[34]:
ngrid = grids.weights.size
ngrid == grdh.ngrid
[34]:
True
格点权重 weight
[35]:
weight = grids.weights
np.allclose(weight, grdh.weight)
[35]:
True
各阶电子坐标偏导数的 AO 轨道格点 ao
[36]:
ao = np.zeros((20, ngrid, nao))
g_start, g_end, g_mem = 0, 0, 2000
for inner_ao, _, _, _ in ni.block_loop(mol, grids, nao, 3, g_mem):
g_end = g_start + ao.shape[1]
ao[:, g_start:g_end, :] = inner_ao
g_start = g_end
[37]:
np.allclose(ao, grdh.ao)
[37]:
True
轨道格点 ao_0 \(\phi_\mu\)
[38]:
ao_0 = ao[0]
np.allclose(ao_0, grdh.ao_0)
[38]:
True
轨道格点一阶导数 ao_1 \(\phi_{r \mu} = \partial_r \phi_\mu\)
[39]:
ao_1 = ao[1:4]
np.allclose(ao_1, grdh.ao_1)
[39]:
True
轨道格点二阶导数 ao_2 \(\phi_{r w \mu} = \partial_r \partial_w \phi_\mu\)
[40]:
XX, XY, XZ, YY, YZ, ZZ = range(4, 10)
ao_2 = ao[([XX, XY, XZ],
[XY, YY, YZ],
[XZ, YZ, ZZ],
), :, :]
np.allclose(ao_2, grdh.ao_2)
[40]:
True
轨道格点二阶导数，但两个坐标分量打包在一个维度中 ao_2T \(\phi_{T \mu} = \partial_{T_1} \partial_{T_2} \phi_\mu\)
[41]:
ao_2T = ao[4:10]
np.allclose(ao_2T, grdh.ao_2T)
[41]:
True
轨道格点三阶导数 ao_3 \(\phi_{r w t \mu} = \partial_r \partial_w \partial_t \phi_\mu\)
[42]:
XXX, XXY, XXZ, XYY, XYZ, XZZ, YYY, YYZ, YZZ, ZZZ = range(10, 20)
ao_3 = ao[([[XXX, XXY, XXZ], [XXY, XYY, XYZ], [XXZ, XYZ, XZZ]],
[[XXY, XYY, XYZ], [XYY, YYY, YYZ], [XYZ, YYZ, YZZ]],
[[XXZ, XYZ, XZZ], [XYZ, YYZ, YZZ], [XZZ, YZZ, ZZZ]],
), :, :]
np.allclose(ao_3, grdh.ao_3)
[42]:
True
轨道格点三阶导数，但其中两个坐标分量打包在一个维度中 ao_3T \(\phi_{T r \mu} = \partial_{T_1} \partial_{T_2} \partial_r \phi_\mu\)
[43]:
XXX, XXY, XXZ, XYY, XYZ, XZZ, YYY, YYZ, YZZ, ZZZ = range(10, 20)
ao_3T = ao[([XXX, XXY, XXZ, XYY, XYZ, XZZ],
[XXY, XYY, XYZ, YYY, YYZ, YZZ],
[XXZ, XYZ, XZZ, YYZ, YZZ, ZZZ],
), :, :]
np.allclose(ao_3T, grdh.ao_3T)
[43]:
True
密度格点 rho_0 \(\rho = D_{\mu \nu} \phi_\mu \phi_\nu\)
[44]:
rho_0 = np.einsum("uv, gu, gv -> g", scfh.D, ao_0, ao_0)
np.allclose(rho_0, grdh.rho_0)
[44]:
True
密度格点一阶导数 rho_1 \(\rho_r = 2 D_{\mu \nu} \phi_{r \mu} \phi_\nu\)
[45]:
rho_1 = 2 * np.einsum("uv, rgu, gv -> rg", scfh.D, ao_1, ao_0)
np.allclose(rho_1, grdh.rho_1)
[45]:
True
密度格点与其一阶导数的合并张量 rho_01
[46]:
rho_01 = np.vstack([rho_0[None, :], rho_1])
np.allclose(rho_1, grdh.rho_1)
[46]:
True
密度格点二阶导数 rho_2 \(\rho_{rw} = 2 D_{\mu \nu} (\phi_{r w \mu} \phi_\nu + \phi_{r \mu} \phi_{w \nu})\)
[47]:
rho_2 = 2 * (
+ np.einsum("uv, rwgu, gv -> rwg", scfh.D, ao_2, ao_0)
+ np.einsum("uv, rgu, wgv -> rwg", scfh.D, ao_1, ao_1)
)
np.allclose(rho_2, grdh.rho_2)
[47]:
True
3.6.4.1.3. 任务 (1.3)¶
首先我们验证 \(\rho_r = D_{\mu \nu} (\phi_{r \mu} \phi_\nu + \phi_\mu \phi_{r \nu})\) 是正确的。我们重新定义 rho_1 变量：
[48]:
rho_1 = (
+ np.einsum("uv, rgu, gv -> rg", scfh.D, grdh.ao_1, grdh.ao_0)
+ np.einsum("uv, gu, rgv -> rg", scfh.D, grdh.ao_0, grdh.ao_1)
)
np.allclose(rho_1, grdh.rho_1)
[48]:
True
但是我们会发现，如果定义临时张量 T \(T_{r \mu \nu} = D_{\mu \nu} \phi_{r \mu} \phi_\nu\) 与 R \(R_{r \mu \nu} = D_{\mu \nu} \phi_\mu \phi_{r \nu}\)，那么这两个张量应当是不相等的。(出于内存大小考虑，程序里的 T 与 R 已经对格点指标求和)
[49]:
T = np.einsum("uv, rgu, gv -> ruv", scfh.D, grdh.ao_1, grdh.ao_0)
R = np.einsum("uv, gu, rgv -> ruv", scfh.D, grdh.ao_0, grdh.ao_1)
np.allclose(T, R)
[49]:
False
尽管上述两个张量在最后两个维度上作转置之后结果便是一样的了，即 \(T_{r \mu \nu} = R_{r \nu \mu}\)：
[50]:
np.allclose(T, R.swapaxes(-1, -2))
[50]:
True
现在暂时不使用 Einstein Summation Convention。根据上述分析，如果对于表达式
不对任何指标求和而给出关于 \(r, \mu, \nu\) 或再多包含一个格点指标维度的张量，那么我们不可以将上式缩减为 \(2 D_{\mu \nu} \phi_{r \mu} \phi_\nu\)。
但若对于表达式
情况就不同了。由于处在求和过程中，因此维度相同的角标 \(\mu, \nu\) 在任意一个被求和项中可以相互交换位置。因此，
同时注意到，由于密度矩阵 \(D_{\mu \nu}\) 是对称矩阵，因此 \(D_{\mu \nu} = D_{\nu \mu}\)，因此
因此，我们可以利用 \(\mu, \nu\) 的对称性简化 \(\rho_r\) 的表达式，但需要注意到其前提是表达式对 \(\mu, \nu\) 同时求和，以及利用到了 \(D_{\mu \nu}\) 作为对称矩阵的性质。
3.6.4.2. 任务 (2)¶
3.6.4.2.1. 任务 (2.1) 演示¶
Exact 交换系数
根据 B3LYP 的原始文献 [Bec93] 的式 (3)，B3LYP 的 Exact 交换系数应当是 \(c_\mathrm{x} = 0.2\)。在 PySCF 中，交换系数可以通过 hybrid_coeff 给出：
[51]:
cx = ni.hybrid_coeff("B3LYPg")
cx
[51]:
0.2
带权重的泛函核及其电子坐标梯度
首先我们对比一下上一节调用 eval_xc 的代码与这一节的代码：
ni.eval_xc("Slater", rho_s_0, deriv=0)[0] # Last section
ni.eval_xc("B3LYPg", grdh.rho_01, deriv=3) # This section
首先泛函名称改变了。
其次，上一节的 rho_s_0 是电子态密度 \(\rho\)，但这一节的 grdh.rho_01 不仅有电子态密度 \(\rho\)，还有密度的电子坐标分量一阶梯度 \(\rho_r\) (回顾 rho_01 应该是 \(4 \times n_\mathrm{grid} \times n_\mathrm{AO}\) 大小的张量)。上一节的泛函是 LDA 泛函，因此只需要代入密度即可；但这一节是 GGA 泛函，因此还需要代入密度梯度量 \(\gamma\)，但同时 \(\gamma = \rho_r \rho_r\)，因此代入密度梯度也是等价的。
最后，上一节由于只通过密度计算泛函的能量，因此只需要泛函核的格点 \(f\) 即可；但这一节还需要导出关于密度 \(\rho\) 与密度梯度量 \(\gamma\) 的格点导数，因此需要设定梯度最高为三阶梯度。
eval_xc 共返回四个变量，分别对应 0-3 阶密度与密度梯度量的导数：
[52]:
grid_exc, grid_vxc, grid_fxc, grid_kxc = ni.eval_xc("B3LYPg", grdh.rho_01, deriv=3)
零阶导数 exc \(f\)
[53]:
exc = grid_exc * grids.weights
[54]:
print(np.allclose(exc, kerh.exc))
True
一阶导数 fr \(f_\rho\), fg \(f_\gamma\)
[55]:
fr = grid_vxc[0] * grids.weights
fg = grid_vxc[1] * grids.weights
[56]:
print(np.allclose(fr, kerh.fr))
print(np.allclose(fg, kerh.fg))
True
True
二阶导数 frr \(f_{\rho \rho}\), frg \(f_{\rho \gamma}\), fgg \(f_{\gamma \gamma}\)
[57]:
frr = grid_fxc[0] * grids.weights
frg = grid_fxc[1] * grids.weights
fgg = grid_fxc[2] * grids.weights
[58]:
print(np.allclose(frr, kerh.frr))
print(np.allclose(frg, kerh.frg))
print(np.allclose(fgg, kerh.fgg))
True
True
True
三阶导数 frrr \(f_{\rho \rho \rho}\), frrg \(f_{\rho \rho \gamma}\), frgg \(f_{\rho \gamma \gamma}\), fggg \(f_{\gamma \gamma \gamma}\)
[59]:
frrr = grid_kxc[0] * grids.weights
frrg = grid_kxc[1] * grids.weights
frgg = grid_kxc[2] * grids.weights
fggg = grid_kxc[3] * grids.weights
[60]:
print(np.allclose(frrr, kerh.frrr))
print(np.allclose(frrg, kerh.frrg))
print(np.allclose(frgg, kerh.frgg))
print(np.allclose(fggg, kerh.fggg))
True
True
True
True
3.6.4.2.2. 任务 (2.2) 可选¶
我们以生成当前密度下的 B3LYP 交换相关能作为例子。回顾交换相关能 eng_xc 的表达式是 \(E_\mathrm{xc} = f \rho\)。
[61]:
eng_xc = 0
grdit = GridIterator(mol, grids, scfh.D, deriv=3, memory=100)
g_count = 0
for it in grdit:
ker_it = KernelHelper(it, "B3LYPg", deriv=0)
eng_xc += (ker_it.exc * it.rho_0).sum()
g_count += 1
print("Iterated times:", g_count)
print("B3LYP xc energy:", eng_xc)
Iterated times: 10
B3LYP xc energy: -14.506876761363433
可以发现与后文中计算得到的交换相关能的结果几乎相等。
3.6.4.3. 任务 (3)¶
3.6.4.3.1. 任务 (3.1)¶
显然是不行的。如果我们按照下述错误的交换相关势
并与正确的交换相关势 xc_v \(v_{\mu \nu}^\mathrm{xc}\) 作对比，可以看到是无法对比成功的：
[62]:
np.allclose(
+ np.einsum("g, gu, gv -> uv", kerh.fr, grdh.ao_0, grdh.ao_0)
+ 4 * np.einsum("g, rg, rgu, gv -> uv", kerh.fg, grdh.rho_1, grdh.ao_1, grdh.ao_0),
xc_v
)
[62]:
False
这是因为等式两边并没有对 \(\mu, \nu\) 求和，因此不能像密度求取过程一样，简单地对于某一个分项对换 \(\mu, \nu\) 而保持其它项不对换，所得结果还能相等。
3.6.4.3.2. 任务 (3.2) 可选¶
尽管交换相关势
不能通过合并分项进行简化，但我们可以利用交换相关势关于 \(\mu, \nu\) 对称的性质，进行对称性的简化：
其中，上式的 \(\mathrm{swap} (\mu, \nu)\) 是指等式右边的所有项对换 \(\mu, \nu\) 的结果；特别地，对于上述表达式而言，指代的是 \(\frac{1}{2} f_\rho \phi_\nu \phi_\mu + 2 f_\gamma \rho_r \phi_{r \nu} \phi_{\mu}\)
现在我们定义两个函数，他们都返回交换相关势 \(v_{\mu \nu}^\mathrm{xc} [\rho]\)，只是 vxc_no_swap 不使用 \(\mathrm{swap} (\mu, \nu)\)，而 vxc_with_swap 使用了 \(\mathrm{swap} (\mu, \nu)\)。定义函数而非变量的目的只是为了方便时间效率测评。
[63]:
def vxc_no_swap():
return (
+ np.einsum("g, gu, gv -> uv", kerh.fr, grdh.ao_0, grdh.ao_0)
+ 2 * np.einsum("g, rg, rgu, gv -> uv", kerh.fg, grdh.rho_1, grdh.ao_1, grdh.ao_0)
+ 2 * np.einsum("g, rg, gu, rgv -> uv", kerh.fg, grdh.rho_1, grdh.ao_0, grdh.ao_1)
)
def vxc_with_swap():
ret = (
+ 0.5 * np.einsum("g, gu, gv -> uv", kerh.fr, grdh.ao_0, grdh.ao_0)
+ 2 * np.einsum("g, rg, rgu, gv -> uv", kerh.fg, grdh.rho_1, grdh.ao_1, grdh.ao_0)
)
return ret + ret.swapaxes(-1, -2)
上述两个做法都可以得到正确的交换相关势 xc_v \(v_{\mu \nu}^\mathrm{xc} [\rho]\)：
[64]:
print(np.allclose(vxc_no_swap(), xc_v))
print(np.allclose(vxc_with_swap(), xc_v))
True
True
但是从时间消耗上，后者的时间消耗大约是前者的 \(2/3\)：
[65]:
%%timeit -r 7 -n 7
vxc_no_swap()
83.2 ms ± 4.69 ms per loop (mean ± std. dev. of 7 runs, 7 loops each)
[66]:
%%timeit -r 7 -n 7
vxc_with_swap()
51.6 ms ± 4.05 ms per loop (mean ± std. dev. of 7 runs, 7 loops each)
这是因为后者少去一次张量缩并 \(f_\gamma \rho_r \phi_{\mu} \phi_{r \nu}\)。张量缩并的耗时比张量求和、标量与张量乘积的时间都要多得多。
3.6.4.3.3. 任务 (3.3) 可选¶
在这里我们暂时不使用 Einstein Summation Convention 以及格点积分，而纯粹地讨论数学推导。这个问题并不是非常显然或者易于回答的。
错误推导
下面我们介绍一种错误推导方法。蓝色高亮表示推导过程中发生变化的项，红色高亮表示推导错误的部分。为了与程序的部分使用的记号比较接近，这里定义 \(F[\rho, \gamma] = f[\rho, \gamma] \rho\)，因此交换相关能可以表示为 \(E_\mathrm{xc} = \int F[\rho, \gamma] \, \mathrm{d} \boldsymbol{r}\)。
如果现在我们让等式左边为交换相关势 \(v_{\mu \nu}^\mathrm{xc} [\rho]\)， \(\frac{\partial F}{\partial \rho}\) 为格点 \(f_\rho\)，\(\frac{\partial F}{\partial \gamma}\) 为格点 \(f_\gamma\)，\(\nabla \rho (\boldsymbol{r})\) 为格点 \(\rho_r\)，\(\phi_\mu (\boldsymbol{r})\) 为格点 \(\phi_\mu\)，\(\nabla \phi_\mu (\boldsymbol{r})\) 为格点 \(\phi_{r \mu}\)，那么我们就立即得到了
这里补充一点。对于张量的梯度点积 \(\nabla \rho (\boldsymbol{r}) \cdot \nabla \phi_\mu (\boldsymbol{r})\)，如果我们令电子坐标分量 \(r \in \{ x, y, z \}\)，同时注意到粗斜体的电子坐标 \(\boldsymbol{r} = (x, y, z)\) (粗斜体的 \(\boldsymbol{r}\) 与斜体的 \(r\) 是两个无关的记号)，那么根据点积的定义，
上面的推导过程其实很有意思，因为看上去都非常合理，并且结果也是正确的。但是，上面的推导在两个关键步骤上错了。第一个问题在前两处红色项上标出：泛函的变分应当在积分的环境下定义。第二个问题在第三处红色项标出：尽管从结果上是正确的，但将 \(\delta \nabla \rho / \delta \rho\) 等同于梯度符号 \(\nabla\) 是完全错误的。下面介绍更合理的两种推导方案，推导过程中的一些细节可以参考 维基百科页面。
为了简化记号，定义 \(\phi (\boldsymbol{r}) = \phi_\mu (\boldsymbol{r}) \phi_\nu (\boldsymbol{r})\)。
根据泛函的定义推导
利用泛函变分等式
另一种做法是利用以下表达式 (可以参考维基百科，也可以参考 Parr [PW89] 式 (A.13))：
之所以上式不用在积分环境中写出，应当是因为对于任何泛函变分的积分，上述表达式总是成立的。因此，
这很像是微积分中的分部积分 (乘法导数的逆运用)；而根据 Gauss 定理 (参考 维基百科 Divergence theorem)，
其中，上式中的 \(\boldsymbol{n}\) 是球面积分中，被积球面的单位法向量，且方向远离被积球体。上式的等式右是环路积分。我们知道等式左边的被积对象是全三维空间；如果我们将符号写得更为严谨一些，
其中 \(\Omega\) 是任意单连通的 (但应包含双氧水分子几乎所有电子云的) 区域。上述等式右边的项一般认为是零；这可以应当可以用波函数 (\(\phi = \phi_\mu \phi_nu\) 两个轨道波函数的乘积，密度量 \(\nabla \rho\) 也与波函数有关) 在无穷远处的渐进性质、以及交换相关泛函量 \(\partial_\gamma F\) 的渐进性质严格证明，但作者还没有能力进行证明。
回到 \(v_{\mu \nu}^\mathrm{xc} [\rho]\) 的推导，通过上面的讨论知道其中的蓝色项积分的值应当为零，因此：
3.6.4.3.4. 任务 (3.4)¶
\(f_\rho\), \(f_\gamma\) 和 \(\rho_r\) 会随着密度的变化而变化。作者认为用 \(D_{\kappa \lambda}\) 替换掉方括号中的 \(\rho\) 是合理的行为。
3.6.4.4. 任务 (4)¶
3.6.4.4.1. 任务 (4.1) 可选¶
电子密度 R \(R_{\mu \nu}\) 定义如下：
[67]:
np.random.seed(1)
R = np.random.random((nao, nao))
R += R.T
我们希望得到的结果是 eng_elec_R \(E_\mathrm{xc} [R_{\mu \nu}]\)：
[68]:
eng_elec_R = scf_eng.energy_elec(dm=R)[0]
eng_elec_R
[68]:
123.10771512394973
由于后面会有很多机会使用 GridHelper 和 KernelHelper，这里就不使用它们；我们顺便借此巩固 PySCF 的 DFT 格点程序。
首先，我们生成 H \(h_{\mu \nu} = t_{\mu \nu} + v_{\mu \nu}^\mathrm{nuc}\), J_R \(J_{\mu \nu} [R_{\kappa \lambda}] = (\mu \nu | \kappa \lambda) R_{\kappa \lambda}\), K_R \(K_{\mu \nu} [R_{\kappa \lambda}] = (\mu \kappa | \nu \lambda) R_{\kappa \lambda}\)，并定义 eri0 \((\mu \nu | \kappa \lambda)\)：
[69]:
eri0 = mol.intor("int2e")
H = mol.intor("int1e_kin") + mol.intor("int1e_nuc")
J_R = np.einsum("uvkl, kl -> uv", eri0, R)
K_R = np.einsum("ukvl, kl -> uv", eri0, R)
随后我们要构建原子轨道格点。由于我们只是要计算交换相关能，因此只需要 GGA 泛函所必需的 \(\rho_r^R\) 的组成部分，即 ao_0 \(\phi_\mu\) 与 ao_1 \(\phi_{r \mu}\)；因此只需要至多 deriv=1 的原子轨道导数。
[70]:
ao = np.zeros((4, grids.weights.size, mol.nao))
ni = dft.numint.NumInt()
g_start = 0
for inner_ao, _, _, _ in ni.block_loop(mol, grids, mol.nao, deriv=1, max_memory=50):
ao[:, g_start:g_start+inner_ao.shape[-2]] = inner_ao
g_start += inner_ao.shape[-2]
ao_0 = ao[0]
ao_1 = ao[1:4]
那么密度格点 rho_0_R \(\rho^R = R_{\mu \nu} \phi_\mu \phi_\nu\) 与 rho_1_R \(\rho_r^R = 2 R_{\mu \nu} \phi_{r \mu} \phi_\nu\) 就可以如下给出：
[71]:
rho_0_R = np.einsum("uv, gu, gv -> g", R, ao_0, ao_0)
rho_1_R = 2 * np.einsum("uv, rgu, gv -> rg", R, ao_1, ao_0)
随后我们计算格点 exc_R \(f^R\)。注意到这是零阶导数项，因此在 eval_xc 中设定 deriv=0 即可。代入 eval_xc 的密度是 \(4 \times n_\mathrm{ngrid} \times n_\mathrm{AO}\) 大小的 \(\rho^R\) 与 \(\rho_r^R\) 的合并张量 rho_01_R。同时需要注意，根据这篇文档的定义，格点 \(f^R\) 应当已经乘以格点权重。
[72]:
rho_01_R = np.vstack([rho_0_R, rho_1_R])
exc_R = ni.eval_xc("B3LYPg", rho_01_R, deriv=0)[0]
exc_R *= grids.weights
最后，我们通过下式给出密度 \(R_{\mu \nu}\) 下的电子态能量 eng_elec_task_R：
[73]:
cx = ni.hybrid_coeff("B3LYPg") # B3LYP cx = 0.2
eng_elec_task_R = (
+ ((H + 0.5 * J_R - 0.25 * cx * K_R) * R).sum()
+ (exc_R * rho_0_R).sum()
)
eng_elec_task_R
[73]:
123.1077151246414
我们可以与 eng_elec_R 的结果作比较：
[74]:
np.allclose(eng_elec_task_R, eng_elec_R)
[74]:
True
3.6.4.4.2. 任务 (4.2)¶
可以。杂化泛函确实是密度的泛函，而不需要是轨道的泛函。需要指出，轨道系数 \(C_{\mu p}\) 一般认为具有比 \(D_{\mu \nu}\) 更大的信息，因为 (对于闭壳层) \(D_{\mu \nu} = 2 C_{\mu i} C_{\nu i}\)，因此密度不包含任何非占据轨道的信息。
之所以提出这个问题，是因为作者很长时间都困扰于杂化泛函中的交换积分 \(K_{\mu \nu} [R_{\kappa \lambda}]\) 是否可以被密度量表示。作者现在的理解是，既然交换积分可以被写作密度矩阵 \(R_{\kappa \lambda}\) 的泛函，那么它应当也是 \(\rho\) 的泛函；只是具体形式如何，作者也不清楚。
上面的措辞十分模糊；为了对这个问题有更清晰的表达，我们先考虑使用格点积分生成在密度 \(R_{\mu \nu}\) 下库伦积分 \(J^R\) (不考虑复共轭)。这一部分解答仍然使用 Einstein Summation Convention。
库伦积分值的格点积分导出
我们先借用上一小题的 R \(R_{\mu \nu}\) 和库伦积分矩阵 J_R \(J_{\mu \nu} [R_{\kappa \lambda}]\)，得到正确的库伦积分值：
[75]:
0.5 * (J_R * R).sum()
[75]:
953.225897393552
现在我们从更为根本的方式出发：
注意到这里的积分不再是对单个的电子坐标 \(\boldsymbol{r}\)，而是要同时对 \(\boldsymbol{r}'\) 进行积分。这里破例使用格点下标 \(g\) 与权重记号 \(w_g\) 并与电子坐标 \(\boldsymbol{r}\) 关联；同时格点下标 \(G\) 与权重记号 \(w_G\) 与电子坐标 \(\boldsymbol{r}'\) 关联。我们同时定义 \(d_{gG} = \Vert \boldsymbol{r}_g - \boldsymbol{r}'_G \Vert_2\)，那么库伦积分值用格点积分表示为
现在我们着手计算上述库伦积分。由于生成 \(d_{gG}\) 的过程非常消耗内存，我们暂且将格点 grids_corse 的精度降到 (30, 110)。需要注意这是非常低精度的格点，只能作定性分析。
[76]:
grids_corse = Mol_H2O2().gen_grids(30, 110)
ngrid_corse = grids_corse.weights.size
ngrid_corse
[76]:
9864
下面定义格点坐标与格点权重 weights_corse \(w_g\)：
[77]:
coords_corse = grids_corse.coords
weights_corse = grids_corse.weights
通过格点坐标则可以定义格点距离 dist_corse \(d_{gG}\)；为了避免相同格点距离为零而导致 \(1 / d_{gG}\) 没有意义，定义当 \(g = G\) 时 \(d_{gG} = +\infty\)：
[78]:
dist_corse = np.linalg.norm(coords_corse[:, None] - coords_corse, axis=-1)
dist_corse += np.diag(np.ones(ngrid_corse) * np.inf)
随后我们生成原子轨道格点 ao_0_corse \(\phi_{g \mu}\) 与密度格点 rho_0_corse \(\rho^R_g\)：
[79]:
ao_0_corse = np.zeros((ngrid_corse, mol.nao))
g_start = 0
for inner_ao, _, _, _ in ni.block_loop(mol, grids_corse, mol.nao, deriv=0, max_memory=50):
ao_0_corse[g_start:g_start+inner_ao.shape[-2]] = inner_ao
g_start += inner_ao.shape[-2]
rho_0_corse = np.einsum("uv, gu, gv -> g", R, ao_0_corse, ao_0_corse)
最后我们可以求得库伦积分值 \(J^R = \frac{1}{2} w_g w_G \rho_g^R \rho_G^R / d_{gG}\)：
[80]:
0.5 * np.einsum("g, G, g, G, gG ->", weights_corse, weights_corse, rho_0_corse, rho_0_corse, 1 / dist_corse)
[80]:
945.167184682469
由于这是粗的格点，因此我们无法将上述值与解析积分的值作精确的对比 (\(\sim 953.23\))；但两者大致接近。
交换积分值的格点积分导出
使用格点积分进行交换积分值计算的策略是相同的。我们仍然先借用上一小题的 R \(R_{\mu \nu}\) 和交换积分矩阵 K_R \(K_{\mu \nu} [R_{\kappa \lambda}]\)，得到正确的交换积分值：
[81]:
0.5 * (K_R * R).sum()
[81]:
999.3581121411252
但若我们着手考虑格点积分的计算方式：
我们发现几乎很难再化简为 \(\rho^R (\boldsymbol{r})\) 的形式了。取而代之地，通常可以定义下述非定域密度量
那么交换积分值可以是
我们定义非定域密度量 rho_0_nlc_corse \(\rho^{1, R}_{gG} = D_{\mu \nu} \phi_{g \mu} \phi_{G \nu}\)：
[82]:
rho_0_nlc_corse = np.einsum("uv, gu, Gv -> gG", R, ao_0_corse, ao_0_corse)
rho_0_nlc_corse.shape
[82]:
(9864, 9864)
那么交换积分值 \(K^R\) 计算可得
[83]:
np.einsum("g, G, gG, gG, gG ->", weights_corse, weights_corse, rho_0_nlc_corse, rho_0_nlc_corse, 1 / dist_corse)
[83]:
array(1982.57724)
显然这个值与解析结果 (\(\sim 999.36\)) 有些差距，但定性上基本是正确的。
总结
在 DFT 教材中 (譬如 Parr [PW89] p.40 式 (2.5.26))，交换能确实没有写作 \(\rho(\boldsymbol{r})\) 而只能写作类似于 \(\rho^1(\boldsymbol{r}, \boldsymbol{r}')\) 的形式；但即使如此，生成交换积分也不需要引入多余的轨道信息，单纯的密度信息应当就足够了。从这个角度上，HF 或者杂化泛函可以看作是严格地在 Hohenberg-Kohn 框架下的近似，即电子态能量仅仅是密度的泛函。
但一方面，不少学者认为 HF 或者杂化 GGA 泛函不属于 Kohn-Sham 框架，为此提出 Generalized Kohn-Sham 框架以及多种符合 Kohn-Sham 轨道的优化 (Optimized Effective Potential) 框架。另一方面，对于双杂化泛函、RPA 型泛函或者一部分 meta-GGA 泛函而言，其能量泛函的表达式已经包括不仅仅密度、还包含轨道的信息；尽管通常也认为，轨道可以看作密度的泛函，但从形式上已经与最初的 Hohenberg-Kohn (也许) 期望的 DFT 近似相差甚远。
“密度的泛函”，可以是极为宽泛的概念，也可以是极为狭隘的概念。在回答题目最初的用词中，“轨道的泛函”指的是在程序中不需要明确写出 \(C_{\mu p}\) 的泛函。其实这句话一方面肯定 HF 和杂化 GGA 泛函是纯粹的密度泛函，另一方面否定了双杂化泛函是纯粹的密度泛函；即使从程序的角度上措辞是严谨的，但这才是真正不伦不类的说法。但不论怎样，作者认为，我们期待的是能更好地解释问题的、在可容忍的代价下给出更为精准结果的泛函；对于是否是“纯粹的密度泛函”这样的标签，尽管确实这会帮助我们更好地厘清问题与泛函发展方向，但作用也仅限于此。
限于作者的能力与知识面，这里不能再讨论更多了。
3.6.4.4.3. 任务 (4.3)¶
不能。我们拿自洽场密度来讨论这个问题。显然，xc_v \(v_{\mu \nu}^\mathrm{xc} [\rho]\) 与密度矩阵 D \(D_{\mu \nu}\) 的乘积求和
[84]:
(xc_v * scfh.D).sum()
[84]:
-18.597820201880936
与 xc_e \(E_\mathrm{xc} [\rho]\) 是不相等的：
[85]:
xc_e
[85]:
-14.506876761363406
这里指出，一般来说
的 \(\stackrel{?}{=}\) 会替换为 \(\neq\)。
这是作者作为初学者很长一段时间没有搞清楚的问题。现在定义积分泛函
其实一个很简单的例子是若 \(E[\rho]\) 是库伦积分 \(J[\rho]\)，即
在这种情况下，通常可以写库伦积分的变分为 (Parr, p.248, (A.11))
因此，
这是一个非常直观的例子，因为泛函变分 \(\delta_\rho J[\rho]\) 与密度 \(\rho\) 乘积的积分得到的是正好 2 倍的 \(J[\rho]\)。
作者认为，普遍地来说，若 \(\frac{F[\rho]}{\rho}\) 与 \(\rho\) 无关，那么泛函变分与密度乘积的积分会还原泛函的值本身。譬如，令 \(F[\rho]\) 是核排斥势能 \(V^\mathrm{nuc} [\rho]\)，那么应当有
因为恰好下式成立：
3.6.5. 参考文献¶
| [Bec93] | (1, 2) Axel D. Becke. Density-functional thermochemistry. III. the role of exact exchange. J. Chem. Phys., 98(7):5648–5652, apr 1993. doi:10.1063/1.464913. |
| [PW89] | (1, 2) Robert G. Parr and Yang Weitao. Density-Functional Theory of Atoms and Molecules (International Series of Monographs on Chemistry, No. 16). Oxford University Press, 1989. ISBN 9780195042795. URL: https://www.amazon.com/Density-Functional-Molecules-International-Monographs-Chemistry/dp/01950427942. |