`_) 的时间消耗；因此一般来说 :code:`timeit` 所给出的平均时间比起 :code:`time` 所给出的挂墙时间要少一些．不过 :code:`timeit` 命令会尝试多次执行，因此时间会跑得长一些．该命令也是通常评测代码效率所更推荐的方法．
+::
+
+ >>> %timeit d = {i for i in range(10000000)}
+ 1.56 s ± 42.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
+
+如果需要在一个 Cell 而非一行代码中中评测时间消耗，则需要使用 :code:`%%time` 与 :code:`%%timeit` 分别代替 :code:`%time` 与 :code:`%timeit`．
+
+.. note ::
+ 在 Windows 下，执行 :code:`%time` 后不会出现 CPU 时间．这是作为操作系统的 Windows 所给予的限制．在非 Windows 系统，包括 WSL，则会显示 CPU 时间．
+
+多矩阵连乘
+::::::::::
+
+对于矩阵连乘 :math:`R_{im} = r_{ij} r_{jk} r_{kl} r_{lm}`，至少有三种做法；若 :math:`\mathbf{r}` 是由 NumPy 生成的随机 50 维矩阵，则
+::
+
+ >>> r = np.random.rand(50, 50)
+ >>> %timeit R = r @ r @ r @ r
+ 26.1 µs ± 1.66 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
+
+ >>> %timeit R = np.einsum("ij, jk, kl, lm -> im", r, r, r, r)
+ 1.72 s ± 6.94 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
+
+ >>> %timeit R = np.einsum("ij, jk, kl, lm -> im", r, r, r, r, optimize=True)
+ 286 µs ± 5.79 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
+
+因此，完成上述命令的最快方式显然是传统的矩阵乘积．对于多矩阵的乘积，:code:`numpy.einsum` 会使用未优化计算复杂度的方式进行计算 (就本例而言，计算复杂度是 :math:`O (N^5)`；但通常我们都会认为上述运算的复杂度在 :math:`O (N^3)` 至 :math:`O (N^2 \log N)` 之间)．而经过优化的 :code:`numpy.einsum` 则可以正确地处理上述计算为不高于 :math:`O (N^3)` 的复杂度，在 50 维下其计算效率比未优化的 :code:`numpy.einsum` 要高效一些，但为此有不小的效率损耗．
+
+不过，如果矩阵维度变小，未优化过的 :code:`numpy.einsum` 反而会快一些．我们现在看看三维矩阵的情况：
+::
+
+ >>> r = np.random.rand(3, 3)
+ >>> %timeit R = r @ r @ r @ r
+ 2.5 µs ± 101 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
+
+ >>> %timeit R = np.einsum("ij, jk, kl, lm -> im", r, r, r, r)
+ 11.9 µs ± 655 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
+
+ >>> %timeit R = np.einsum("ij, jk, kl, lm -> im", r, r, r, r, optimize=True)
+ 217 µs ± 2.25 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
+
+因此，论效率上，公式表达式与程序代码关系不友好的矩阵相乘记号是最快的；而使用 :code:`numpy.einsum` 不是最效率的；同时，如果处理的问题维度较小，或不优化与优化的计算复杂度没有改变时，使用未优化的 :code:`numpy.einsum` 有时比优化的版本还快一些．
+
+当然，作为开发方法的工作者，自然会对效率上的要求有所降低，因此，通常情况下直接使用优化的 :code:`numpy.einsum` 未尝不可，因为它的代码本身与公式的对应关系非常显然．很多时候，教程中就会使用这种可能偏低效的方法了．
+
+矩阵构建
+--------
+
+创建一个新的全零矩阵可以通过两种途径：
+::
+
+ >>> # 通过向 numpy.zeros 传入 tuple 型数组
+ >>> np.zeros((2, 3))
+ >>> # 也可以通过已有矩阵所导出的 tuple 作为变量
+ >>> np.zeros(a.shape)
+ >>> # 或者使用 numpy.zeros_like 来构建与传入矩阵相同维度的全零矩阵
+ >>> np.zeros_like(a)
+ array([[0, 0, 0],
+ [0, 0, 0]])
+
+创建对角阵则可以使用
+::
+
+ >>> np.eye(3)
+ array([[1., 0., 0.],
+ [0., 1., 0.],
+ [0., 0., 1.]])
+
+而经常地，我们会从本征值向量 :math:`\boldsymbol{e}` 展开成二维分子轨道 Fock 矩阵 :math:`\mathbf{F}`，这个过程通常可以由下述技巧完成：
+::
+
+ >>> dim = 4
+ >>> e = np.arange(dim)
+ >>> e * np.eye(dim)
+ array([[0., 0., 0., 0.],
+ [0., 1., 0., 0.],
+ [0., 0., 2., 0.],
+ [0., 0., 0., 3.]])
+
+而在处理 MP2 计算时，其分母项中会出现张量 :math:`\mathcal{E}_{ab}^{ij} = \varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b`；在这里我们以比较简单的矩阵 :math:`\mathcal{E}_{c}^{k} = \varepsilon_k - \varepsilon_c` 来举例子．我们可以通过改变矩阵的维度的技巧获得：
+::
+
+ >>> # 定义变量
+ >>> dim = 4
+ >>> k = np.arange(-1, -dim - 1, -1)
+ >>> c = np.arange(2, 2 * dim + 2, 2)
+ >>> k
+ array([-1, -2, -3, -4])
+ >>> c
+ array([2, 4, 6, 8])
+ >>> # 计算矩阵
+ >>> k.reshape(-1, 1) # 或 k.reshape(4, 1)
+ array([[-1],
+ [-2],
+ [-3],
+ [-4]])
+ >>> k.reshape(-1, 1) - c # 即 E_c^k 矩阵
+ array([[ -3, -5, -7, -9],
+ [ -4, -6, -8, -10],
+ [ -5, -7, -9, -11],
+ [ -6, -8, -10, -12]])
+
+其中用到了矩阵或向量的大小重新定义的函数 :code:`numpy.reshape`．该函数输入为新矩阵大小的 tuple 型变量；也支持用 -1 让程序推断该维度的值：
+::
+
+ >>> a.reshape(3, 2)
+ >>> a.reshape(-1, 2)
+ >>> a.T
+ array([[0, 3],
+ [1, 4],
+ [2, 5]])
+
+如果只是将矩阵压平成为向量，还可以使用 :code:`numpy.ravel` 函数：
+::
+
+ >>> a.reshape(-1)
+ >>> a.ravel()
+ array([0, 1, 2, 3, 4, 5])
+
+向量视图
+--------
+
+在这次教程中，出现了少数代码，这些代码的理解必须要基于简单的 NumPy 向量的向量视图 (View) 的概念．这些概念不存在于 Fortran 与 C，它与 Python 本身不具有明确指针多少有些关系．我们知道，Fortran 的向量通常就可以当做指针来看待；而 C 或 C++ 的向量还多一种引用的描述方式．对于 Python，它一般不太容易写出其引用与指针，因此我们不太容易把握在完成向量操作时，是否真的对原来的向量作了操作，导致了原始数据的破坏；或者是否我们复制出一个新的向量，造成了内存空间的浪费．
+
+NumPy 的向量类可以简单地看作由底层数据和表面形状 (shape) 构成．NumPy 很少采用真正的深层复制 (Deep Copy)，即很少将底层数据复制到另一个变量中．深层复制的通常做法是
+::
+
+ >>> d = a.copy()
+
+以后对 :code:`d` 的任何数据、形状的改动，都不会影响 :code:`a`．反之亦然．
+
+而更多时候是引用．它不将数据复制出来，但包含表面形状的信息．在最为简单的情况下，可以直接理解为一种引用．例如，向量的索引相当于对其对应的原始数据的引用：
+::
+
+ >>> d = np.arange(4)
+ >>> d[2] = 10
+ >>> d
+ array([ 0, 1, 10, 3])
+
+但还有一些更为特殊的操作，这些不能简单地看作通常的引用．例如我们可以令 :code:`v` 是 :code:`d` 的一种视窗：:code:`v` 是 :code:`d` 若干个元素的引用；对 :code:`v` 的形状的改变不会对 :code:`d` 产生影响，但对其数据的改动则会直接改动 :code:`d` 的数据：
+::
+
+ >>> d = np.arange(4)
+ >>> # v 是 d 的视图，并非将数据复制给了 v，数据还是从 d 读出来
+ >>> v = d[0:3:2]
+ >>> v
+ array([0, 2])
+ >>> # 更改 v 的形状对 d 没有影响
+ >>> v.shape = 2, 1
+ >>> v
+ array([[0],
+ [2]])
+ >>> d.shape
+ (4,)
+ >>> # 更改 v 的数据对 d 有影响，这类似于引用关系
+ >>> v[:] = np.array([[-2], [-6]])
+ >>> d
+ array([-2, 1, -6, 3])
+ >>> # 但下面这句语句并非是给视图更改数据
+ >>> # 创建了新的向量赋值给 v，自此 v 与 d 不存在相互关系
+ >>> v = np.array([[-3], [-9]])
+ >>> d
+ array([-2, 1, -6, 3])
+
+其它函数
+--------
+
+:code:`numpy.linalg.norm` 模长函数
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+对于向量模长，可以简单地调用它来计算；对于矩阵，它等同于化为向量：
+::
+
+ >>> np.linalg.norm(a)
+ >>> np.linalg.norm(a.ravel())
+ 7.416198487095663
+
+:code:`numpy.linalg.eigh` 本征系统
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+若现在有对称矩阵 :math:`\textbf{f}`，则其本征值与本征向量可以借助该函数获得：
+::
+
+ >>> f = np.array(
+ ... [[2., 3., 3.],
+ ... [3., 2.33, 3.],
+ ... [3., 3., 3.]])
+ >>> eig, vec = np.linalg.eigh(f)
+ >>> # 本征值
+ >>> eig
+ array([-0.85515066, -0.27773769, 8.46288834])
+ >>> # 本征向量
+ >>> vec
+ array([[-0.7806939 , -0.29943621, -0.54850251],
+ [ 0.61076205, -0.55134403, -0.56832163],
+ [ 0.13223751, 0.77868974, -0.61331519]])
+
+:code:`numpy.allclose` 判断矩阵相同
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+现在我们对本征方才的本征系统的简单性质作验证．首先，本征向量所构成的矩阵是正交矩阵，即其逆应当等于转置：
+::
+
+ >>> vec.T
+ array([[-0.7806939 , 0.61076205, 0.13223751],
+ [-0.29943621, -0.55134403, 0.77868974],
+ [-0.54850251, -0.56832163, -0.61331519]])
+ >>> np.linalg.inv(vec)
+ array([[-0.7806939 , 0.61076205, 0.13223751],
+ [-0.29943621, -0.55134403, 0.77868974],
+ [-0.54850251, -0.56832163, -0.61331519]])
+ >>> np.allclose(vec.T, np.linalg.inv(vec))
+ True
+
+本征系统本质上可以看作是一种矩阵对角化．我们验证一下对角化前后的矩阵是否一致：
+::
+
+ >>> f_after_diag = vec @ (eig * np.eye(eig.shape[0])) @ vec.T
+ >>> f_after_diag
+ array([[2. , 3. , 3. ],
+ [3. , 2.33, 3. ],
+ [3. , 3. , 3. ]])
+ >>> np.allclose(f_after_diag, f)
+ True
diff --git a/source/xyg3_energy.ipynb b/source/xyg3_energy.ipynb
index 0ca0663..5850989 100644
--- a/source/xyg3_energy.ipynb
+++ b/source/xyg3_energy.ipynb
@@ -1,1289 +1,1289 @@
-{
- "cells": [
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# B3LYP 能量与 XYG3 能量"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "在这一节中，我们会在闭壳层下给出 B3LYP 与 XYG3 能量；同时验证与 DFT 计算有关的矩阵，以及了解较为基础的 DFT 格点积分有关的代码书写．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "from pyscf import gto, scf, dft, mp, ao2mo, lib\n",
- "import numpy as np"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "这一节我们还会对格点作简单的可视化．因此需要引入下面的工具．其中，`%matplotlib notebook` 是在 Jupyter Notebook 中嵌入交互式的图像的工具．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "# %matplotlib notebook\n",
- "\n",
- "from mpl_toolkits.mplot3d import Axes3D\n",
- "from matplotlib import pyplot as plt\n",
- "from matplotlib.colors import LogNorm"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "**提示**\n",
- "\n",
- "这里把 `%matplotlib notebook` 注释掉了，但这只是因为生成该网页的 `nbsphinx` 似乎未必允许 GUI 形式的输出．如果你是用的是货真价实的 Jupyter Notebook，这个 Magic Command 对于 3D 图像的交互可视化确实是很有用的．\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## XYG3 型双杂化泛函参考文献"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "XYG3 型双杂化泛函 (xDH, XYG3 type of Double Hybrid functional) 是一系列引入精确交换能与 PT2 形式交换能的泛函家族，其最初的泛函是 XYG3．其它典型的泛函有 XYGJ-OS、xDH-PBE0 等．\n",
- "\n",
- "* XYG3\n",
- " * Zhang, Y.; Xu, X.; Goddard, W. A. Doubly Hybrid Density Functional for Accurate Descriptions of Nonbond Interactions, Thermochemistry, and Thermochemical Kinetics. *Proc. Natl. Acad. Sci. U.S.A.* **2009**, *106* (13), 4963–4968.\n",
- " * doi: [10.1073/pnas.0901093106](https://dx.doi.org/10.1073/pnas.0901093106)\n",
- "* XYGJ-OS\n",
- " * Zhang, I. Y.; Xu, X.; Jung, Y.; Goddard, W. A. A Fast Doubly Hybrid Density Functional Method close to Chemical Accuracy Using a Local Opposite Spin Ansatz. *Proc. Natl. Acad. Sci. U.S.A.* **2011**, *108* (50), 19896–19900.\n",
- " * doi: [10.1073/pnas.1115123108](https://doi.org/10.1073/pnas.1115123108)\n",
- "* xDH-PBE0\n",
- " * Zhang, I. Y.; Su, N. Q.; Brémond, É. A. G.; Adamo, C.; Xu, X. Doubly Hybrid Density Functional xDH-PBE0 from a Parameter-Free Global Hybrid Model PBE0. *J. Chem. Phys.* **2012**, *136* (17), 174103.\n",
- " * doi: [10.1063/1.3703893](https://doi.org/10.1063/1.3703893)\n",
- "\n",
- "对 XYG3 泛函的一些测评、原理与展望等说明，有且不限于以下的文献：\n",
- "\n",
- "* The XYG3 Type of Doubly Hybrid Density Functionals\n",
- " * Su, N. Q.; Xu, X. *WIREs Comput. Mol. Sci.* **2016**, *6* (6), 721–747.\n",
- " * doi: [10.1002/wcms.1274](https://doi.org/10.1002/wcms.1274)\n",
- "* Development of New Density Functional Approximations\n",
- " * Su, N. Q.; Xu, X. *Annu. Rev. Phys. Chem.* **2017**, *68* (1), 155–182.\n",
- " * doi: [10.1146/annurev-physchem-052516-044835](https://doi.org/10.1146/annurev-physchem-052516-044835)\n",
- "* Doubly Hybrid Density Functionals That Correctly Describe Both Density and Energy for Atoms\n",
- " * Su, N. Q.; Zhu, Z.; Xu, X. *Proc. Natl. Acad. Sci. U.S.A.* **2018**, *115* (10), 2287–2292.\n",
- " * doi: [10.1073/pnas.1713047115](https://doi.org/10.1073/pnas.1713047115)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## 预备工作"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 顶层函数计算 B3LYP 能量"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "由于我们准备与 Gaussian 核对能量，因此这里使用 VWN3 型的 LDA 相关泛函作为 B3LYP 中 LDA 相关泛函．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mol = gto.Mole()\n",
- "mol.atom = \"\"\"\n",
- "O 1.0 0.0 0.0\n",
- "H 1.0 1.0 0.0\n",
- "H 1.0 0.0 1.0\n",
- "\"\"\"\n",
- "mol.basis = \"6-31G\"\n",
- "mol.build()"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "scf_eng = dft.RKS(mol)\n",
- "scf_eng.xc = \"b3lypg\" # compare that to gaussian\n",
- "scf_eng.grids.atom_grid = (99, 590)\n",
- "scf_eng.grids.build()\n",
- "scf_eng.conv_tol = 1e-13\n",
- "energy_scf = scf_eng.kernel()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### B3LYP 轨道构造的 MP2 相关能"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "我们可以直接将上述的 B3LYP 轨道代入 MP2 相关能公式中．求得的 MP2 相关能将会是 XYG3 能量构成的一部分．在这里我们可以预先生成之．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "mp2_eng = mp.MP2(scf_eng)\n",
- "energy_mp2_corr, _ = mp2_eng.kernel()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### B3LYP 重要中间矩阵"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "在上一份文档中，我们已经对众多 HF 中间矩阵作了较为充分的说明．在这里，我们就简单地将这些重要矩阵写在一起．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "nao = mol.nao\n",
- "nmo = scf_eng.mo_energy.shape[0]\n",
- "nelec = mol.nelectron\n",
- "nocc = mol.nelec[0]\n",
- "nvir = nmo - nocc\n",
- "\n",
- "S = mol.intor('int1e_ovlp_sph')\n",
- "T = mol.intor('int1e_kin_sph')\n",
- "Vnuc = mol.intor('int1e_nuc_sph')\n",
- "eri = mol.intor('int2e_sph')\n",
- "\n",
- "C = scf_eng.mo_coeff\n",
- "Co = C[:, :nocc]\n",
- "Cv = C[:, nocc:]\n",
- "e = scf_eng.mo_energy\n",
- "eo = e[:nocc]\n",
- "ev = e[nocc:]\n",
- "\n",
- "D = scf_eng.make_rdm1()\n",
- "F = scf_eng.get_fock()\n",
- "\n",
- "J = scf_eng.get_j()\n",
- "K = scf_eng.get_k()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "通常来说，由于 HF 与 B3LYP 的计算方式几乎一样，因此上一份文档中提及的大多数性质都能保证；但关于 Fock 矩阵的验证上，需要注意到只有一部分性质仍然存在："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "print(\"Eigenvalues of Fock :\", np.allclose(C.T @ F @ C, np.eye(nmo) * e))\n",
- "print(\"Fock matrix decompose :\", np.allclose(T + Vnuc + J - 0.5 * K, F))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "这也就意味着对于 DFT 方法，我们仍然需要了解 Fock 矩阵的构建方式．对于能量亦是如此："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose((D * (T + Vnuc + 0.5 * J - 0.25 * K)).sum() + mol.energy_nuc(), \\\n",
- " scf_eng.energy_tot())"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "下面一节中，我们就将解决如何求取 Fock 矩阵与 B3LYP 能量的问题．"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## B3LYP 交换相关势与交换相关能"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 通过 PySCF 高级函数生成"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "我们首先通过 PySCF 较为顶层的函数来生成交换相关势 $V^\\mathrm{xc}_{\\mu \\nu} [\\mathbf{D}]$ 与交换相关能 $E^\\mathrm{xc} [\\mathbf{D}]$．这两者分别可以通过 Fock 矩阵与 B3LYP 能量来验证．下面的代码会一次性地生成这两者，以及对密度矩阵 $D_{\\mu \\nu}$ 的电子数的数值求和："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "xc_n, xc_e, xc_v = \\\n",
- " scf_eng._numint.nr_rks(mol, scf_eng.grids, scf_eng.xc, D)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "这个函数的主体是 `scf_eng._numint`，它是 `dft.numint.Numint` 类型；我们可以直接调用这个类型的函数进行计算，其结果是相同的："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "xc_n_ni, xc_e_ni, xc_v_ni = \\\n",
- " dft.numint.nr_rks(scf_eng._numint, mol, scf_eng.grids, scf_eng.xc, D)\n",
- "print(np.allclose(xc_n_ni, xc_n))\n",
- "print(np.allclose(xc_e_ni, xc_e))\n",
- "print(np.allclose(xc_v_ni, xc_v))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "该函数的传入参数分别是分子构型、格点信息、泛函信息以及 AO 密度矩阵．"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 输出 1：电子数和"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "我们知道，水分子的电子数是 10 个．该函数的第一个输出即可以验证电子数是否合理．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "print(xc_n)\n",
- "np.allclose(xc_n, nelec)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 输出 2：交换相关能"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "第二个输出是交换相关能．获得该能量后，我们可以验证 B3LYP 总能量了．但在此之前，我们需要先知道 B3LYP 作为杂化泛函，其杂化 HF 型交换能，即交换积分的比例系数 $c_\\mathrm{x}$．这里我们使用下面的方法导出系数．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "c_x = scf_eng._numint.hybrid_coeff(scf_eng.xc)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "随后我们可以通过杂化泛函的计算通式给出 B3LYP 的总能量．这里同时列出杂化泛函的计算式与 HF 的计算式，相以比对：\n",
- "\n",
- "\\begin{align}\n",
- "E^\\textrm{Hyb} &= D_{\\mu \\nu} (T_{\\mu} + V_{\\mu}^\\mathrm{Nuc} + \\frac{1}{2} J_{\\mu \\nu} [\\mathbf{D}] - \\frac{1}{4} c_x K_{\\mu \\nu} [\\mathbf{D}]) + E^\\mathrm{xc} [\\mathbf{D}] + E^\\mathrm{Nuc} \\\\\n",
- "E^\\textsf{HF} &= D_{\\mu \\nu} (T_{\\mu} + V_{\\mu}^\\mathrm{Nuc} + \\frac{1}{2} J_{\\mu \\nu} [\\mathbf{D}] - \\frac{1}{4} K_{\\mu \\nu} [\\mathbf{D}]) + E^\\mathrm{Nuc}\n",
- "\\end{align}"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose((D * (T + Vnuc + 0.5 * J - 0.25 * c_x * K)).sum() + xc_e + mol.energy_nuc(), \\\n",
- " scf_eng.energy_tot())"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "**提示**\n",
- "\n",
- "一般地，在程序中对 DFT 部分的计算与对交换积分的计算是分开的；因此，从实现的角度上讲，$E^\\mathrm{xc} [\\mathbf{D}]$ 不包含交换积分．\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 输出 3：交换相关势"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "类似于交换相关能，当获得交换相关势后，我们就可以验证 Fock 矩阵．我们也对杂化泛函与 HF 的 Fock 矩阵作比对：\n",
- "\n",
- "\\begin{align}\n",
- "F_{\\mu \\nu}^\\textrm{Hyb} &= T_{\\mu \\nu} + V_{\\mu \\nu}^\\mathrm{Nuc} + J_{\\mu \\nu} [\\mathbf{D}] - \\frac{1}{2} c_x K_{\\mu \\nu} [\\mathbf{D}] + V_{\\mu \\nu}^\\mathrm{xc} [\\mathbf{D}] \\\\\n",
- "F_{\\mu \\nu}^\\textsf{HF} &= T_{\\mu \\nu} + V_{\\mu \\nu}^\\mathrm{Nuc} + J_{\\mu \\nu} [\\mathbf{D}] - \\frac{1}{2} K_{\\mu \\nu} [\\mathbf{D}]\n",
- "\\end{align}"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose(T + Vnuc + J - 0.5 * c_x * K + xc_v, F)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "在 PySCF 中，杂化泛函的 `get_veff` 函数与 HF 中的功能略微不同；它同时包含交换相关势 $V_{\\mu \\nu}^\\mathrm{xc} [\\mathbf{D}]$ 的贡献："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose(scf_eng.get_veff(), J - 0.5 * c_x * K + xc_v)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## 格点积分方法"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### DFT 计算使用的格点"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "上面使用了顶层的 `dft.numint.NumInt.nr_rks` 函数生成了计算 B3LYP 能量时所需要的交换相关势 $V^\\mathrm{xc}_{\\mu \\nu} [\\mathbf{D}]$ 与交换相关能 $E^\\mathrm{xc} [\\mathbf{D}]$；但该函数比较高级．如果我们需要作一些底层的修改，该函数就不合适了．为此，我们会初步地对其中的底层实现作基础的了解．\n",
- "\n",
- "DFT 的矩阵与能量的计算关键是格点积分部分．之后的几小节将会在假定格点信息已经具备的情况下，通过对格点数据的求和来得到我们期望的交换相关势与能量．\n",
- "\n",
- "在这一小节中，我们仅仅是研究格点信息，因此会使用非常小的格点数据进行讨论．我们所选取的格点大小是径向 4 个点，角向 14 个点的格点积分．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "grid_test = dft.gen_grid.Grids(mol)\n",
- "grid_test.atom_grid = (4, 14)\n",
- "grid_test.build()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "我们知道，得到 DFT 的交换相关势相当于需要泛函对密度或动能密度的一阶梯度，因此这里设置梯度量为 1．通过 PySCF 的中层函数，我们可以给出目前格点下所需要的众多信息："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "ao_deriv = 1\n",
- "grid_ao, grid_mask, grid_weight, grid_coords = \\\n",
- " next(scf_eng._numint.block_loop(mol, grid_test, nao, ao_deriv, 2000))\n",
- "print(\"Shape of arrays:\")\n",
- "print(\"ao :\", grid_ao.shape)\n",
- "print(\"mask :\", grid_mask.shape)\n",
- "print(\"weight :\", grid_weight.shape)\n",
- "print(\"coords :\", grid_coords.shape)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "其中，`grid_coords` 变量代表格点的坐标，`grid_weight` 变量代表对应坐标的格点权重．我们可以使用下面的代码进行可视化．注意到因为一部分格点的权重为零，因此在作对数图可视化时需要加一个小量 (在实际进行积分时则不需要引入)；这里假的小量为 $10^{-5}$．这一般来说不影响我们对格点数据的定性观察．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "weight_shifted = grid_weight + 1e-5\n",
- "fig = plt.figure()\n",
- "ax = fig.add_subplot(111, projection='3d')\n",
- "sc = ax.scatter(grid_coords.T[0], grid_coords.T[1], grid_coords.T[2], \\\n",
- " c=weight_shifted, \\\n",
- " norm=LogNorm(vmin=weight_shifted.min(), vmax=weight_shifted.max()))\n",
- "fig.colorbar(sc)\n",
- "# fig.show()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "关于变量 `grid_ao` 的意义，以及 `grid_weight` 的意义，将会在下面两小节中进行说明．\n",
- "\n",
- "刚才的格点 `grid_test` 数量显然太少，只有 $4 \\times 14 = 168$ 个点．在实际的量化计算中，显然这是不合适的．我们现在回到真实环境下所使用的格点 (99, 590)，来进行后续的计算．\n",
- "\n",
- "对于变量 `grid_mask`，猜测为因为在一些原子轨道积分 `grid_ao` 下，值总是为零或非常小，于是为了加速格点求和，忽略这些格点积分的数值．在这份笔记中，我们暂时不需要关心这些为加速运算速度而引入的数组．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "ao_deriv = 1\n",
- "grid_ao, grid_mask, grid_weight, grid_coords = \\\n",
- " next(scf_eng._numint.block_loop(mol, scf_eng.grids, nao, ao_deriv, 2000))\n",
- "print(\"Shape of arrays:\")\n",
- "print(\"ao :\", grid_ao.shape)\n",
- "print(\"mask :\", grid_mask.shape)\n",
- "print(\"weight :\", grid_weight.shape)\n",
- "print(\"coords :\", grid_coords.shape)\n",
- "ngrids = grid_weight.shape[0]"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "**注意**\n",
- "\n",
- "上面的代码中，我们对 `scf_eng._numint.block_loop` 函数使用的操作是将迭代过程中的一组数据取出．这也就意味着，实际上该函数是一个迭代函数，它完全可能生成不只一次格点信息．该函数传入的第 5 个参数是当前的内存信息；如果内存不足以存下足够的格点，那么格点将会多次生成．\n",
- "\n",
- "因此，在真正计算大分子或大格点时，不可以写成上面的样子；而必须要通过迭代器实现：\n",
- "\n",
- "```python\n",
- "grid_ao, grid_mask, grid_weight, grid_coords = \\\n",
- " next(scf_eng._numint.block_loop(mol, scf_eng.grids, nao, ao_deriv, 2000))\n",
- " # Your manuplation on DFT grid comes here ...\n",
- "```\n",
- "\n",
- "之所以在这里，我们可以只取一组数据就能进行 DFT 计算，是因为当前体系的分子足够小、格点足够小、基组也足够小，因此能在 2GB 左右的内存内进行给单积分计算．\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 输出 1：电子数和格点求取，格点积分实现"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "首先，我们会使用格点积分处理最简单的电子数的和．这里的记号会混用 Einstein Convention 与积分记号，因此只会在引入格点积分时出现求和记号，其它时候不引入；符号上多少会复杂一些，但应当不会产生歧义．\n",
- "\n",
- "我们知道，电子密度的和可以用下述的方式给出：($D_{ij}$ 与 $D_{\\mu \\nu}$ 分别是 MO 基组与 AO 基组密度矩阵，它们两者不同但有下述的关系)\n",
- "\n",
- "\\begin{align}\n",
- "n &= D_{ij} = 2 \\delta_{ij} = \\mathrm{tr} (\\mathbf{D}^\\mathrm{MO}) \\\\\n",
- "&= C_{\\mu i} S_{\\mu \\nu} C_{\\nu j} = \\mathrm{tr} (\\mathbf{C}^\\mathrm{occ, T} \\mathbf{S} \\mathbf{C}^\\mathrm{occ}) = \\mathrm{tr} (\\mathbf{C}^\\mathrm{occ, T} \\mathbf{C}^\\mathrm{occ} \\mathbf{S}) = \\mathrm{tr} (\\mathbf{D}^\\mathrm{AO} \\mathbf{S}) \\\\\n",
- "&= (C_{\\mu i} C_{\\nu j} \\delta_{ij}) S_{\\mu \\nu} = D_{\\mu \\nu} S_{\\mu \\nu}\n",
- "\\end{align}\n",
- "\n",
- "下面我们只验证 $n = D_{\\mu \\nu} S_{\\mu \\nu}$："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose((D * S).sum(), nelec)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "上式完全是矩阵的运算．但如果我们希望使用格点积分，对空间变量 $\\boldsymbol{r}$ 显式地积分，那么上式可以根据重叠矩阵的定义写为\n",
- "\n",
- "\\begin{equation}\n",
- "n = D_{\\mu \\nu} \\int \\phi_{\\mu} (\\boldsymbol{r}) \\phi_{\\nu} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
- "\\end{equation}\n",
- "\n",
- "随后我们就可以进行格点积分计算了．格点积分的实现方式是对空间坐标积分量离散化并继而求和．具体地来说，如果现在的被积函数是 $f(\\boldsymbol{r})$，那么对于一系列已知坐标的、角标记为 $g$ 的格点 (即坐标集合) $\\{ \\boldsymbol{r}_g \\}_g$ 与格点权重 (即实数集合) $\\{ w_g \\}_g$，有\n",
- "\n",
- "\\begin{equation}\n",
- "\\int f(\\boldsymbol{r}) \\mathrm{d} \\boldsymbol{r} \\simeq \\sum_{g} f(\\boldsymbol{r}_g) w_g\n",
- "\\end{equation}\n",
- "\n",
- "如果现在将 $f(\\boldsymbol{r}_g)$ 记为 $f_g$，即将函数 $f(\\boldsymbol{r})$ 映射到格点集合 $\\{ f_g \\}_g$，那么我们可以使用 Einstein Convention 简记上述积分为\n",
- "\n",
- "\\begin{equation}\n",
- "\\int f(\\boldsymbol{r}) \\mathrm{d} \\boldsymbol{r} \\simeq f_g w_g\n",
- "\\end{equation}\n",
- "\n",
- "回到电子数的格点积分．如果我们定义对原子轨道基组的格点映射为 $\\phi_{\\mu} (\\boldsymbol{r}) \\mapsto \\{ \\phi_{g \\mu} \\}_g$，那么电子数积分就可以写为\n",
- "\n",
- "\\begin{equation}\n",
- "n = D_{\\mu \\nu} \\phi_{g \\mu} \\phi_{g \\nu} w_g\n",
- "\\end{equation}\n",
- "\n",
- "下面我们就通过格点积分验证电子数的和．这里指出，[上文中](#DFT-格点) `block_loop` 函数所声称的 `grid_ao` 是四部分构成的，其中的第一部分 `grid_ao[0]` 是原子轨道基组的格点 $\\phi_{g \\mu}$，而剩下三部分 `grid_ao[1:4]` 则是原子轨道基组的 $(x, y, z)$ 梯度分量 $\\nabla_t \\phi_{g \\mu}$，其中 $t \\in (x, y, z)$．在电子数求和时只需要 $\\phi_{g \\mu}$ 即可，而 $\\nabla_t \\phi_{g \\mu}$ 在交换相关势与能量求取是会使用．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "xc_n_sum = float(np.einsum(\"uv, gu, gv, g\", D, grid_ao[0], grid_ao[0], \\\n",
- " grid_weight, optimize=True))\n",
- "xc_n_sum"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose(xc_n_sum, xc_n)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "我们甚至可以顺便验证 $S_{\\mu \\nu} = \\phi_{g \\mu} \\phi_{g \\nu} w_g$，不过因为格点积分的精度不算太高，因此在使用 `np.allclose` 验证矩阵相同时，需要把默认精度 $10^{-8}$ 降低到 $10^{-7}$."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose(np.einsum(\"gu, gv, g\", grid_ao[0], grid_ao[0], \\\n",
- " grid_weight, optimize=True), \\\n",
- " S, atol=1e-7)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 输出 2：交换相关能格点求取"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "交换相关能从定义上，其积分表达式为\n",
- "\n",
- "\\begin{equation}\n",
- "E^\\textrm{xc} [\\rho] = \\int f^\\textrm{xc} [\\rho, \\nabla \\rho] \\rho(\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "显然地，上式经过格点化后可以写为\n",
- "\n",
- "\\begin{equation}\n",
- "E^\\textrm{xc} [\\rho] = f^\\textrm{xc}_g \\rho_g w_g\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "看起来非常简洁．但是上式的 $f^\\textrm{xc}_g$ 与 $\\rho_g$ 都还不是已知量．$\\rho_g$ 的构造方法比较容易：\n",
- "\n",
- "\\begin{equation}\n",
- "\\rho_g = \\rho (\\boldsymbol{r}_g) = D_{\\mu \\nu} \\phi_{\\mu} (\\boldsymbol{r}_g) \\phi_{\\nu} (\\boldsymbol{r}_g) = D_{\\mu \\nu} \\phi_{g \\mu} \\phi_{g \\nu}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "为了构造 $f^\\textrm{xc}_g$，我们需要传入关于 $\\rho, \\boldsymbol{\\nabla} \\rho$ 的信息，并且了解作为 DFT 近似核心的 $f^\\textrm{xc}_g$．PySCF 中，这些计算已经通过接口函数 `eval_xc` 完成，其底层一般来说是 LibXC 库函数．我们需要准备的信息不多，只需要下面三行代码即足够："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "grid_rho = np.einsum(\"uv, tgu, gv -> tg\", D, grid_ao, grid_ao[0], optimize=True)\n",
- "grid_rho[1:4] *= 2\n",
- "grid_exc, grid_vxc = scf_eng._numint.eval_xc(scf_eng.xc, grid_rho, deriv=1)[:2]"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "关于上面三行代码，我们还有不少需要说明的地方．首先，关于第一行，实际上是合并了两个操作，同时导出了 $\\rho_g$ (`grid_rho[0]`) 与 $\\nabla_t \\rho_g$ (`grid_rho[1:4]`)．$\\rho_g$ 的导出方式上面已经叙述了；而 $\\nabla_t \\rho_g$ 的导出方式如下："
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\\begin{align}\n",
- "\\nabla_t \\rho_g &= \\nabla_t \\rho (\\boldsymbol{r}_g) = D_{\\mu \\nu} \\nabla_t \\big( \\phi_{\\mu} (\\boldsymbol{r}_g) \\phi_{\\nu} (\\boldsymbol{r}_g) \\big) \\\\\n",
- "&= D_{\\mu \\nu} \\nabla_t \\phi_{\\mu} (\\boldsymbol{r}_g) \\phi_{\\nu} (\\boldsymbol{r}_g) + D_{\\mu \\nu} \\phi_{\\mu} (\\boldsymbol{r}_g) \\nabla_t \\phi_{\\nu} (\\boldsymbol{r}_g) \\\\\n",
- "&= D_{\\mu \\nu} (\\nabla_t \\phi_{g \\mu}) \\phi_{g \\nu} + D_{\\mu \\nu} \\phi_{g \\mu} (\\nabla_t \\phi_{g \\nu}) \\\\\n",
- "&= 2 D_{\\mu \\nu} (\\nabla_t \\phi_{g \\mu}) \\phi_{g \\nu}\n",
- "\\end{align}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "联系到 `grid_ao[1:4]` 储存了 $\\nabla_t \\phi_{g \\mu}$，而 `grid_ao[0]` 储存的是 $\\phi_{g \\mu}$，应当不难理解上述代码的第一行为何可以同时导出 $\\rho_g$ 和 $\\nabla_t \\rho_g$ 了．但是我们看到上面公式中有两倍，但 $\\nabla_t \\rho_g$ 公式中只有一倍，因此我们有必要对 `grid_ao[1:4]` 部分单独乘上 2 倍．"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "第三行代码则是具体执行计算 $f^\\mathrm{xc}_g$ 的代码．根据该代码的 [API 文档](https://sunqm.github.io/pyscf/dft.html#pyscf.dft.libxc.eval_xc)，由于我们不仅需要计算泛函的能量本身，还要计算其梯度以构造 Fock 矩阵，因此设置 `deriv=1`；由此，`eval_xc` 函数的输出尽管是一个四元素 tuple，但只有前两个元素非 `None`，且分别为泛函能量格点 `grid_exc` 与梯度格点 `grid_vxc`；后两个元素本来应当是泛函二阶梯度与三阶梯度量，但在 GGA 的 SCF 中不需要计算它们．关于梯度格点，将在下一小节中说明．"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "我们在构造 `grid_rho[0]` 后，我们还可以先简单验证其可以导出电子数之和对于水分子应当为 10：\n",
- "\n",
- "\\begin{equation}\n",
- "n = \\rho_g w_g\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "(grid_rho[0] * grid_weight).sum()"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "能量值的导出也相当容易：我们可以发现，`grid_exc` 的维度就是格点大小："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "grid_exc.shape"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "那么只要简单地套到能量格点公式 $E^\\textrm{xc} [\\rho] = f^\\textrm{xc}_g \\rho_g w_g$ 中，就可以给出交换相关能．我们验证这一结果．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "xc_e_sum = (grid_exc * grid_rho[0] * grid_weight).sum()\n",
- "xc_e_sum"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose(xc_e_sum, xc_e)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 输出 3：交换相关势格点求取"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "交换相关势的求取公式仍然不算复杂．在这一小节，我们先列出公式与代码，再对它们作较为详细的解释；下一小节中，我们会简单地推导交换相关势的计算公式．\n",
- "\n",
- "交换相关势的计算公式为\n",
- "\n",
- "\\begin{align}\n",
- "V^\\textrm{xc}_{\\mu \\nu} [\\rho] \n",
- "&= w_g (\\partial_\\rho f^\\mathrm{xc}_{g}) \\phi_{g \\mu} \\phi_{g \\nu} \\\\\n",
- "&\\quad + 2 w_g (\\partial_\\gamma f^\\mathrm{xc}_{g}) (\\nabla_t \\rho_g) (\\nabla_t \\phi_{g \\mu}) \\phi_{g \\nu} \\\\\n",
- "&\\quad + 2 w_g (\\partial_\\gamma f^\\mathrm{xc}_{g}) (\\nabla_t \\rho_g) (\\nabla_t \\phi_{g \\nu}) \\phi_{g \\mu}\n",
- "\\end{align}"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "grid_frho, grid_fgamma = grid_vxc[:2]\n",
- "\n",
- "xc_v_sum = np.einsum(\"g, g, gu, gv -> uv\", grid_weight, grid_frho, \\\n",
- " grid_ao[0], grid_ao[0], optimize=True)\n",
- "xc_v_sum += 2 * np.einsum(\"g, g, tg, tgu, gv -> uv\", grid_weight, grid_fgamma, \\\n",
- " grid_rho[1:4], grid_ao[1:4], grid_ao[0], optimize=True)\n",
- "xc_v_sum += 2 * np.einsum(\"g, g, tg, tgv, gu -> uv\", grid_weight, grid_fgamma, \\\n",
- " grid_rho[1:4], grid_ao[1:4], grid_ao[0], optimize=True)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "我们可以验证上面导出的交换相关势 `xc_v_sum` 与使用高级函数 `nr_rks` 所导出的 `xc_v` 是相同的："
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "np.allclose(xc_v_sum, xc_v)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "上述代码中需要解释的是第一行．`grid_vxc` 是 `eval_xc` 函数所导出的第二个变量．我们知道 `eval_xc` 会返回四个变量，分别是关于能量格点、一阶导数格点、二阶导数格点、三阶导数格点．那么 `grid_vxc` 所表征的是一阶导数格点．\n",
- "\n",
- "尽管我们现在接触的是以 GGA 为底层的 B3LYP，但这里有必要简单了解 meta-GGA．meta-GGA 的泛函变量一般至多是下述四个数量值的几种：\n",
- "\n",
- "\\begin{equation}\n",
- "\\rho, \\quad \n",
- "\\gamma = \\boldsymbol{\\nabla} \\rho \\cdot \\boldsymbol{\\nabla} \\rho, \\quad\n",
- "\\boldsymbol{\\nabla}^2 \\rho, \\quad\n",
- "\\tau = | \\boldsymbol{\\nabla} \\phi_\\mu |^2\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "因此，在程序导出泛函的梯度时，也会相应地导出\n",
- "\n",
- "\\begin{equation}\n",
- "\\left(\n",
- "\\frac{\\partial f^\\mathrm{xc}}{\\partial \\rho}, \n",
- "\\frac{\\partial f^\\mathrm{xc}}{\\partial \\gamma}, \n",
- "\\frac{\\partial f^\\mathrm{xc}}{\\partial (\\boldsymbol{\\nabla}^2 \\rho)}, \n",
- "\\frac{\\partial f^\\mathrm{xc}}{\\partial \\tau}\n",
- "\\right)\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "上面的写法不适合得到文档的行内表示，因此后面会简记上面的记号为\n",
- "\n",
- "\\begin{equation}\n",
- "(\n",
- "\\partial_\\rho f^\\mathrm{xc},\n",
- "\\partial_\\gamma f^\\mathrm{xc},\n",
- "\\partial_{\\boldsymbol{\\nabla}^2 \\rho} f^\\mathrm{xc},\n",
- "\\partial_\\tau f^\\mathrm{xc}\n",
- ")\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "上面写出的函数偏导实际上是关于空间变量 $\\boldsymbol{r}$ 的函数，只是因为记号复杂而将其省略．如果我们写得详细一些，譬如 GGA Kernel 下的 $\\partial_\\rho f^\\mathrm{xc}$，应该表示为"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\\begin{equation}\n",
- "(\\partial_\\rho f^\\mathrm{xc}) (\\boldsymbol{r}) = \\frac{\\partial f^\\mathrm{xc} (\\rho, \\gamma)}{\\partial \\rho (\\boldsymbol{r})}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "当要化为格点时，即取 $\\boldsymbol{r}$ 为确定的坐标 $\\boldsymbol{r}_g$ 时，简记下述记号\n",
- "\n",
- "\\begin{equation}\n",
- "\\partial_\\rho f^\\mathrm{xc}_g = (\\partial_\\rho f^\\mathrm{xc}) (\\boldsymbol{r}_g)\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "对于 $\\partial_\\gamma f^\\mathrm{xc}_g$ 亦如此．这样，我们就解释了上面代码所对应的公式的记号了．"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "对于 RKS 而言，`grid_vxc` 变量确实就是四维 tuple 变量；因此，第一行代码\n",
- "\n",
- "```\n",
- "grid_frho, grid_fgamma = grid_vxc[:2]\n",
- "```\n",
- "\n",
- "所作的就是将 $\\partial_\\rho f^\\mathrm{xc}_g$ 赋值给 `grid_frho`；$\\partial_\\gamma f^\\mathrm{xc}_g$ 赋值给 `grid_fgamma`；这些变量的维度都是格点大小．由于 B3LYP 是 GGA，因此 $\\partial_{\\boldsymbol{\\nabla}^2 \\rho} f^\\mathrm{xc}, \\partial_\\tau f^\\mathrm{xc}$ 没有意义，即 `grid_vxc[2:4]` 不应是有意义的量．下面的代码可以简单说明这些．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "print(grid_vxc.__len__())\n",
- "print(grid_vxc[2:4])\n",
- "print(grid_frho.shape)\n",
- "print(grid_fgamma.shape)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "**提示**\n",
- "\n",
- "但对于 UKS 而言，由于 $\\alpha$ 与 $\\beta$ 密度不同，因此 `eval_xc` 给出的一阶导数格点尽管仍然是四维 tuple 变量，但矩阵的数量总共有 9 个．更多信息需要参考 [API 文档](https://sunqm.github.io/pyscf/dft.html#pyscf.dft.libxc.eval_xc)．\n",
- "\n",
- "
"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 交换相关势的推导"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "这里的推导仅仅是表明交换相关势中每一项、以及它们的系数是如何导出的；更为详细与严谨的推导还需要参考其它文献．"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "泛函的变分一般来说只有放在积分的环境下才具有意义．根据 [Wikipedia](https://en.wikipedia.org/wiki/Functional_derivative#Formula) 的说明，如果能量泛函记为"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\\begin{equation}\n",
- "E^\\textrm{xc} [\\rho] = \\int f^\\textrm{xc} (\\rho, \\gamma) \\rho(\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "我们可以对 GGA 型的泛函一次变分记为"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\\begin{align}\n",
- "\\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\rho (\\boldsymbol{r})} \\phi (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r} \n",
- "&= \\int \\left[ \\frac{\\partial f^\\textrm{xc} (\\rho, \\gamma)}{\\partial \\rho (\\boldsymbol{r})} \\phi (\\boldsymbol{r})\n",
- "+ \\frac{\\partial f^\\textrm{xc} (\\rho, \\gamma)}{\\partial \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r})} \\cdot \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\phi (\\boldsymbol{r}) \\right] \n",
- "\\, \\mathrm{d} \\boldsymbol{r} \\\\\n",
- "&= \\int \\left[ \\frac{\\partial f^\\textrm{xc} (\\rho, \\gamma)}{\\partial \\rho (\\boldsymbol{r})} \\phi (\\boldsymbol{r})\n",
- "+ \\frac{\\partial f^\\textrm{xc} (\\rho, \\gamma)}{\\partial \\gamma (\\boldsymbol{r})} 2 \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r}) \\cdot \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\phi (\\boldsymbol{r}) \\right] \n",
- "\\, \\mathrm{d} \\boldsymbol{r}\n",
- "\\end{align}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "第一个等式是变分本身的定义 (需要利用广义 Stocks 定理)；第二个等式利用到连续偏导以及关于 $\\gamma$ 与 $\\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho$ 关系的一个性质："
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\\begin{equation}\n",
- "\\frac{\\partial \\gamma (\\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho)}{\\partial \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r})}\n",
- "= \\frac{\\partial (\\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho \\cdot \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho)}{\\partial \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r})}\n",
- "= 2 \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r})\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "上面出现记号 $\\boldsymbol{\\nabla}_{\\boldsymbol{r}'}$；这仅仅是普通的梯度记号，一般会记为 $\\nabla$；这里记得较为复杂，只是强调梯度的坐标变量 $\\boldsymbol{r}'$ 与被积元变量 $\\boldsymbol{r}$ 不是相同的．如果我们使用 Einstein Convention，对 $\\boldsymbol{r}'$ 更换为 $t \\in (x, y, z)$ 并在实际计算过程中对 $t$ 求和，来记 GGA 的一次变分；同时按照 [上一小节](#输出-3：交换相关势格点求取) 简化偏导数记号；那么可以写为"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\\begin{equation}\n",
- "\\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\rho (\\boldsymbol{r})} \\phi (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
- "= \\int \\left[ (\\partial_{\\rho} f^\\mathrm{xc}) \\phi + 2 (\\partial_{\\gamma} f^\\mathrm{xc}) (\\nabla_t \\rho) (\\nabla_t \\phi) \\right] \\, \\mathrm{d} \\boldsymbol{r}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "现在我们考察上述结论是如何应用到 $V^\\textrm{xc}_{\\mu \\nu} [\\rho]$ 的导出过程中的．我们知道，$V^\\textrm{xc}_{\\mu \\nu} [\\rho]$ 是 Fock 矩阵构成的一部分；Fock 矩阵的构成方式是将 Hartree-Fock 能量对分子轨道进行变分而产生的．能量可以分为单电子、库伦积分、交换积分与泛函积分的贡献，那么 Fock 矩阵也会分为这四部分；关于这点我们已经在 [上文](#输出-3：交换相关势) 中验证过了．那么泛函积分对 Fock 矩阵的贡献可以表示为\n",
- "\n",
- "\\begin{equation}\n",
- "F_{\\mu \\nu} [\\rho] \\leftarrow V^\\textrm{xc}_{\\mu \\nu} [\\rho] = \\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\phi_{\\mu} (\\boldsymbol{r})} \\phi_{\\nu} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "由于密度的表达式可以写为 (该公式不使用 Einstein Convention)\n",
- "\n",
- "\\begin{equation}\n",
- "\\rho(\\boldsymbol{r}) = \\rho [\\{ \\phi_{\\mu} \\}_\\mu] (\\boldsymbol{r}) = \\int \\sum_{\\nu} \\phi_{\\nu} (\\boldsymbol{r}) \\phi_{\\nu} (\\boldsymbol{r}') \\delta(\\boldsymbol{r} - \\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r}'\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "因此密度对轨道的变分可以写作 (该公式不使用 Einstein Convention)\n",
- "\n",
- "\\begin{align}\n",
- "\\frac{\\delta \\rho}{\\delta \\phi_\\mu (\\boldsymbol{r})} \n",
- "&= \\int \\sum_{\\nu} \\frac{\\partial \\phi_{\\nu}}{\\partial \\phi_{\\mu} (\\boldsymbol{r})} \\phi_{\\nu} (\\boldsymbol{r}') \\delta(\\boldsymbol{r} - \\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r}' \\\\\n",
- "&= \\int \\phi_{\\mu} (\\boldsymbol{r}') \\delta(\\boldsymbol{r} - \\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r}' = \\phi_\\mu (\\boldsymbol{r})\n",
- "\\end{align}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "因此，根据连续变分规则，有\n",
- "\\begin{equation}\n",
- "V^\\textrm{xc}_{\\mu \\nu} [\\rho] = \\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\rho (\\boldsymbol{r})} \\frac{\\delta \\rho}{\\delta \\phi_\\mu (\\boldsymbol{r})} \\phi_{\\nu} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
- "= \\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\rho (\\boldsymbol{r})} \\phi_\\mu (\\boldsymbol{r}) \\phi_{\\nu} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "那么我们用上面已经有的结论，令 $\\phi (\\boldsymbol{r}) = \\phi_\\mu (\\boldsymbol{r}) \\phi_{\\nu} (\\boldsymbol{r})$，并注意到"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\\begin{equation}\n",
- "\\nabla_t \\phi = \\nabla_t (\\phi_\\mu \\phi_\\nu) = (\\nabla_t \\phi_\\mu) \\phi_\\nu + (\\nabla_t \\phi_\\nu) \\phi_\\mu\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "那么\n",
- "\n",
- "\\begin{equation}\n",
- "V^\\textrm{xc}_{\\mu \\nu} [\\rho] = \\int \\left[ (\\partial_{\\rho} f^\\mathrm{xc}) \\phi_\\mu \\phi_\\nu + 2 (\\partial_{\\gamma} f^\\mathrm{xc}) (\\nabla_t \\rho) (\\nabla_t \\phi_\\mu) \\phi_\\nu + 2 (\\partial_{\\gamma} f^\\mathrm{xc}) (\\nabla_t \\rho) (\\nabla_t \\phi_\\nu) \\phi_\\mu \\right] \\, \\mathrm{d} \\boldsymbol{r}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "将上式格点化，就立即得到 [上一小节](#输出-3：交换相关势格点求取) 所使用的格点公式了．"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## XYG3 能量"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### 手动设置 B3LYP 泛函形式"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "XYG3 作为非自洽泛函，其密度与轨道从 B3LYP 获得，但其非自洽能量泛函则使用自己的泛函形式．如果要在尚未实现 XYG3 的软件中获得 XYG3 能量，就必须手动设置泛函参数．在实际使用 XYG3 参数前，我们先尝试对 B3LYP 进行手动的参数设置，并验证结果是否与上面的计算相同．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "b3try_eng = dft.RKS(mol)\n",
- "b3try_eng.xc = \"HF*0.2 + .08*LDA + .72*B88, .81*LYP + .19*VWN3\" # compare that to gaussian\n",
- "b3try_eng.grids.atom_grid = (99, 590)\n",
- "b3try_eng.grids.build()\n",
- "\n",
- "b3try_eng.kernel()\n",
- "print(\"Compare My B3LYPG: \", np.allclose(b3try_eng.e_tot, scf_eng.e_tot))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "### XYG3 泛函形式与能量"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "XYG3 型双杂化泛函属于第五阶泛函 (按 Predew 的 [Jacob 阶梯](https://dx.doi.org/10.1063/1.1390175) 说法)，其 [定义](https://dx.doi.org/10.1073/pnas.0901093106) 是\n",
- "\n",
- "\\begin{equation}\n",
- "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}\n",
- "\\end{equation}\n",
- "\n",
- "对应到程序中，每一项的系数则展开为\n",
- "\n",
- "\\begin{equation}\n",
- "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}\n",
- "\\end{equation}\n",
- "\n",
- "对于 XYG3，其系数的确定是\n",
- "\n",
- "\\begin{equation}\n",
- "c_1 = 0.8033, \\quad c_2 = 0.2107, \\quad c_3 = 0.3211, \\quad c_4 = 1 - c_3\n",
- "\\end{equation}\n",
- "\n",
- "因此，XYG3 的泛函形式为\n",
- "\n",
- "\\begin{equation}\n",
- "E_\\mathrm{xc}^\\textsf{XYG3} = 0.8033 E_\\mathrm{x}^\\textsf{exact} - 0.0140 E_\\mathrm{x}^\\textsf{LDA} + 0.2107 E_\\mathrm{x}^\\textsf{GGA} + 0.6789 E_\\mathrm{c}^\\textsf{GGA} + 0.3211 E_\\mathrm{c}^\\textsf{PT2}\n",
- "\\end{equation}"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "在 PySCF 中，泛函形式的确定通过字符串进行定义．该字符串的定义方式可以很灵活，下面的代码只是其中一种定义方式．其中，逗号前后分别分割了交换能与相关能的部分．由于 PySCF 中不对 PT2 型相关能在 DFT 部分进行定义，因此在定义泛函的时候，不将 $0.3211 E_\\mathrm{c}^\\textsf{PT2}$ 写入字符串．\n",
- "\n",
- "当泛函的形式确定后，我们只要简单地把自洽场过程中收敛的 B3LYP 密度代入 XYG3 型泛函能量表达式中，即可得到最终的 XYG3 能量了．"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "nscf_eng = dft.RKS(mol)\n",
- "nscf_eng.xc = \"0.8033*HF - 0.0140*LDA + 0.2107*B88, 0.6789*LYP\"\n",
- "nscf_eng.grids.atom_grid = (99, 590)\n",
- "nscf_eng.grids.build()\n",
- "\n",
- "xyg3_energy = nscf_eng.energy_tot(dm=D) + mp2_eng.e_corr * 0.3211\n",
- "print(\"XYG3 energy allclose: \", np.allclose(xyg3_energy, -0.76282393305943E+02))"
- ]
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Python 3",
- "language": "python",
- "name": "python3"
- },
- "language_info": {
- "codemirror_mode": {
- "name": "ipython",
- "version": 3
- },
- "file_extension": ".py",
- "mimetype": "text/x-python",
- "name": "python",
- "nbconvert_exporter": "python",
- "pygments_lexer": "ipython3",
- "version": "3.6.6"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# B3LYP 能量与 XYG3 能量"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "在这一节中，我们会在闭壳层下给出 B3LYP 与 XYG3 能量；同时验证与 DFT 计算有关的矩阵，以及了解较为基础的 DFT 格点积分有关的代码书写．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from pyscf import gto, scf, dft, mp, ao2mo, lib\n",
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "这一节我们还会对格点作简单的可视化．因此需要引入下面的工具．其中，`%matplotlib notebook` 是在 Jupyter Notebook 中嵌入交互式的图像的工具．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# %matplotlib notebook\n",
+ "\n",
+ "from mpl_toolkits.mplot3d import Axes3D\n",
+ "from matplotlib import pyplot as plt\n",
+ "from matplotlib.colors import LogNorm"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "**提示**\n",
+ "\n",
+ "这里把 `%matplotlib notebook` 注释掉了，但这只是因为生成该网页的 `nbsphinx` 似乎未必允许 GUI 形式的输出．如果你是用的是货真价实的 Jupyter Notebook，这个 Magic Command 对于 3D 图像的交互可视化确实是很有用的．\n",
+ "\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## XYG3 型双杂化泛函参考文献"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "XYG3 型双杂化泛函 (xDH, XYG3 type of Double Hybrid functional) 是一系列引入精确交换能与 PT2 形式交换能的泛函家族，其最初的泛函是 XYG3．其它典型的泛函有 XYGJ-OS、xDH-PBE0 等．\n",
+ "\n",
+ "* XYG3\n",
+ " * Zhang, Y.; Xu, X.; Goddard, W. A. Doubly Hybrid Density Functional for Accurate Descriptions of Nonbond Interactions, Thermochemistry, and Thermochemical Kinetics. *Proc. Natl. Acad. Sci. U.S.A.* **2009**, *106* (13), 4963–4968.\n",
+ " * doi: [10.1073/pnas.0901093106](https://dx.doi.org/10.1073/pnas.0901093106)\n",
+ "* XYGJ-OS\n",
+ " * Zhang, I. Y.; Xu, X.; Jung, Y.; Goddard, W. A. A Fast Doubly Hybrid Density Functional Method close to Chemical Accuracy Using a Local Opposite Spin Ansatz. *Proc. Natl. Acad. Sci. U.S.A.* **2011**, *108* (50), 19896–19900.\n",
+ " * doi: [10.1073/pnas.1115123108](https://doi.org/10.1073/pnas.1115123108)\n",
+ "* xDH-PBE0\n",
+ " * Zhang, I. Y.; Su, N. Q.; Brémond, É. A. G.; Adamo, C.; Xu, X. Doubly Hybrid Density Functional xDH-PBE0 from a Parameter-Free Global Hybrid Model PBE0. *J. Chem. Phys.* **2012**, *136* (17), 174103.\n",
+ " * doi: [10.1063/1.3703893](https://doi.org/10.1063/1.3703893)\n",
+ "\n",
+ "对 XYG3 泛函的一些测评、原理与展望等说明，有且不限于以下的文献：\n",
+ "\n",
+ "* The XYG3 Type of Doubly Hybrid Density Functionals\n",
+ " * Su, N. Q.; Xu, X. *WIREs Comput. Mol. Sci.* **2016**, *6* (6), 721–747.\n",
+ " * doi: [10.1002/wcms.1274](https://doi.org/10.1002/wcms.1274)\n",
+ "* Development of New Density Functional Approximations\n",
+ " * Su, N. Q.; Xu, X. *Annu. Rev. Phys. Chem.* **2017**, *68* (1), 155–182.\n",
+ " * doi: [10.1146/annurev-physchem-052516-044835](https://doi.org/10.1146/annurev-physchem-052516-044835)\n",
+ "* Doubly Hybrid Density Functionals That Correctly Describe Both Density and Energy for Atoms\n",
+ " * Su, N. Q.; Zhu, Z.; Xu, X. *Proc. Natl. Acad. Sci. U.S.A.* **2018**, *115* (10), 2287–2292.\n",
+ " * doi: [10.1073/pnas.1713047115](https://doi.org/10.1073/pnas.1713047115)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 预备工作"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 顶层函数计算 B3LYP 能量"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "由于我们准备与 Gaussian 核对能量，因此这里使用 VWN3 型的 LDA 相关泛函作为 B3LYP 中 LDA 相关泛函．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mol = gto.Mole()\n",
+ "mol.atom = \"\"\"\n",
+ "O 1.0 0.0 0.0\n",
+ "H 1.0 1.0 0.0\n",
+ "H 1.0 0.0 1.0\n",
+ "\"\"\"\n",
+ "mol.basis = \"6-31G\"\n",
+ "mol.build()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "scf_eng = dft.RKS(mol)\n",
+ "scf_eng.xc = \"b3lypg\" # compare that to gaussian\n",
+ "scf_eng.grids.atom_grid = (99, 590)\n",
+ "scf_eng.grids.build()\n",
+ "scf_eng.conv_tol = 1e-13\n",
+ "energy_scf = scf_eng.kernel()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### B3LYP 轨道构造的 MP2 相关能"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "我们可以直接将上述的 B3LYP 轨道代入 MP2 相关能公式中．求得的 MP2 相关能将会是 XYG3 能量构成的一部分．在这里我们可以预先生成之．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mp2_eng = mp.MP2(scf_eng)\n",
+ "energy_mp2_corr, _ = mp2_eng.kernel()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### B3LYP 重要中间矩阵"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "在上一份文档中，我们已经对众多 HF 中间矩阵作了较为充分的说明．在这里，我们就简单地将这些重要矩阵写在一起．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "nao = mol.nao\n",
+ "nmo = scf_eng.mo_energy.shape[0]\n",
+ "nelec = mol.nelectron\n",
+ "nocc = mol.nelec[0]\n",
+ "nvir = nmo - nocc\n",
+ "\n",
+ "S = mol.intor('int1e_ovlp_sph')\n",
+ "T = mol.intor('int1e_kin_sph')\n",
+ "Vnuc = mol.intor('int1e_nuc_sph')\n",
+ "eri = mol.intor('int2e_sph')\n",
+ "\n",
+ "C = scf_eng.mo_coeff\n",
+ "Co = C[:, :nocc]\n",
+ "Cv = C[:, nocc:]\n",
+ "e = scf_eng.mo_energy\n",
+ "eo = e[:nocc]\n",
+ "ev = e[nocc:]\n",
+ "\n",
+ "D = scf_eng.make_rdm1()\n",
+ "F = scf_eng.get_fock()\n",
+ "\n",
+ "J = scf_eng.get_j()\n",
+ "K = scf_eng.get_k()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "通常来说，由于 HF 与 B3LYP 的计算方式几乎一样，因此上一份文档中提及的大多数性质都能保证；但关于 Fock 矩阵的验证上，需要注意到只有一部分性质仍然存在："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(\"Eigenvalues of Fock :\", np.allclose(C.T @ F @ C, np.eye(nmo) * e))\n",
+ "print(\"Fock matrix decompose :\", np.allclose(T + Vnuc + J - 0.5 * K, F))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "这也就意味着对于 DFT 方法，我们仍然需要了解 Fock 矩阵的构建方式．对于能量亦是如此："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose((D * (T + Vnuc + 0.5 * J - 0.25 * K)).sum() + mol.energy_nuc(), \\\n",
+ " scf_eng.energy_tot())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "下面一节中，我们就将解决如何求取 Fock 矩阵与 B3LYP 能量的问题．"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## B3LYP 交换相关势与交换相关能"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 通过 PySCF 高级函数生成"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "我们首先通过 PySCF 较为顶层的函数来生成交换相关势 $V^\\mathrm{xc}_{\\mu \\nu} [\\mathbf{D}]$ 与交换相关能 $E^\\mathrm{xc} [\\mathbf{D}]$．这两者分别可以通过 Fock 矩阵与 B3LYP 能量来验证．下面的代码会一次性地生成这两者，以及对密度矩阵 $D_{\\mu \\nu}$ 的电子数的数值求和："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "xc_n, xc_e, xc_v = \\\n",
+ " scf_eng._numint.nr_rks(mol, scf_eng.grids, scf_eng.xc, D)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "这个函数的主体是 `scf_eng._numint`，它是 `dft.numint.Numint` 类型；我们可以直接调用这个类型的函数进行计算，其结果是相同的："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "xc_n_ni, xc_e_ni, xc_v_ni = \\\n",
+ " dft.numint.nr_rks(scf_eng._numint, mol, scf_eng.grids, scf_eng.xc, D)\n",
+ "print(np.allclose(xc_n_ni, xc_n))\n",
+ "print(np.allclose(xc_e_ni, xc_e))\n",
+ "print(np.allclose(xc_v_ni, xc_v))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "该函数的传入参数分别是分子构型、格点信息、泛函信息以及 AO 密度矩阵．"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 输出 1：电子数和"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "我们知道，水分子的电子数是 10 个．该函数的第一个输出即可以验证电子数是否合理．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(xc_n)\n",
+ "np.allclose(xc_n, nelec)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 输出 2：交换相关能"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "第二个输出是交换相关能．获得该能量后，我们可以验证 B3LYP 总能量了．但在此之前，我们需要先知道 B3LYP 作为杂化泛函，其杂化 HF 型交换能，即交换积分的比例系数 $c_\\mathrm{x}$．这里我们使用下面的方法导出系数．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "c_x = scf_eng._numint.hybrid_coeff(scf_eng.xc)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "随后我们可以通过杂化泛函的计算通式给出 B3LYP 的总能量．这里同时列出杂化泛函的计算式与 HF 的计算式，相以比对：\n",
+ "\n",
+ "\\begin{align}\n",
+ "E^\\textrm{Hyb} &= D_{\\mu \\nu} (T_{\\mu} + V_{\\mu}^\\mathrm{Nuc} + \\frac{1}{2} J_{\\mu \\nu} [\\mathbf{D}] - \\frac{1}{4} c_x K_{\\mu \\nu} [\\mathbf{D}]) + E^\\mathrm{xc} [\\mathbf{D}] + E^\\mathrm{Nuc} \\\\\n",
+ "E^\\textsf{HF} &= D_{\\mu \\nu} (T_{\\mu} + V_{\\mu}^\\mathrm{Nuc} + \\frac{1}{2} J_{\\mu \\nu} [\\mathbf{D}] - \\frac{1}{4} K_{\\mu \\nu} [\\mathbf{D}]) + E^\\mathrm{Nuc}\n",
+ "\\end{align}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose((D * (T + Vnuc + 0.5 * J - 0.25 * c_x * K)).sum() + xc_e + mol.energy_nuc(), \\\n",
+ " scf_eng.energy_tot())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "**提示**\n",
+ "\n",
+ "一般地，在程序中对 DFT 部分的计算与对交换积分的计算是分开的；因此，从实现的角度上讲，$E^\\mathrm{xc} [\\mathbf{D}]$ 不包含交换积分．\n",
+ "\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 输出 3：交换相关势"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "类似于交换相关能，当获得交换相关势后，我们就可以验证 Fock 矩阵．我们也对杂化泛函与 HF 的 Fock 矩阵作比对：\n",
+ "\n",
+ "\\begin{align}\n",
+ "F_{\\mu \\nu}^\\textrm{Hyb} &= T_{\\mu \\nu} + V_{\\mu \\nu}^\\mathrm{Nuc} + J_{\\mu \\nu} [\\mathbf{D}] - \\frac{1}{2} c_x K_{\\mu \\nu} [\\mathbf{D}] + V_{\\mu \\nu}^\\mathrm{xc} [\\mathbf{D}] \\\\\n",
+ "F_{\\mu \\nu}^\\textsf{HF} &= T_{\\mu \\nu} + V_{\\mu \\nu}^\\mathrm{Nuc} + J_{\\mu \\nu} [\\mathbf{D}] - \\frac{1}{2} K_{\\mu \\nu} [\\mathbf{D}]\n",
+ "\\end{align}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose(T + Vnuc + J - 0.5 * c_x * K + xc_v, F)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "在 PySCF 中，杂化泛函的 `get_veff` 函数与 HF 中的功能略微不同；它同时包含交换相关势 $V_{\\mu \\nu}^\\mathrm{xc} [\\mathbf{D}]$ 的贡献："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose(scf_eng.get_veff(), J - 0.5 * c_x * K + xc_v)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## 格点积分方法"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### DFT 计算使用的格点"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "上面使用了顶层的 `dft.numint.NumInt.nr_rks` 函数生成了计算 B3LYP 能量时所需要的交换相关势 $V^\\mathrm{xc}_{\\mu \\nu} [\\mathbf{D}]$ 与交换相关能 $E^\\mathrm{xc} [\\mathbf{D}]$；但该函数比较高级．如果我们需要作一些底层的修改，该函数就不合适了．为此，我们会初步地对其中的底层实现作基础的了解．\n",
+ "\n",
+ "DFT 的矩阵与能量的计算关键是格点积分部分．之后的几小节将会在假定格点信息已经具备的情况下，通过对格点数据的求和来得到我们期望的交换相关势与能量．\n",
+ "\n",
+ "在这一小节中，我们仅仅是研究格点信息，因此会使用非常小的格点数据进行讨论．我们所选取的格点大小是径向 4 个点，角向 14 个点的格点积分．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "grid_test = dft.gen_grid.Grids(mol)\n",
+ "grid_test.atom_grid = (4, 14)\n",
+ "grid_test.build()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "我们知道，得到 DFT 的交换相关势相当于需要泛函对密度或动能密度的一阶梯度，因此这里设置梯度量为 1．通过 PySCF 的中层函数，我们可以给出目前格点下所需要的众多信息："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ao_deriv = 1\n",
+ "grid_ao, grid_mask, grid_weight, grid_coords = \\\n",
+ " next(scf_eng._numint.block_loop(mol, grid_test, nao, ao_deriv, 2000))\n",
+ "print(\"Shape of arrays:\")\n",
+ "print(\"ao :\", grid_ao.shape)\n",
+ "print(\"mask :\", grid_mask.shape)\n",
+ "print(\"weight :\", grid_weight.shape)\n",
+ "print(\"coords :\", grid_coords.shape)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "其中，`grid_coords` 变量代表格点的坐标，`grid_weight` 变量代表对应坐标的格点权重．我们可以使用下面的代码进行可视化．注意到因为一部分格点的权重为零，因此在作对数图可视化时需要加一个小量 (在实际进行积分时则不需要引入)；这里假的小量为 $10^{-5}$．这一般来说不影响我们对格点数据的定性观察．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "weight_shifted = grid_weight + 1e-5\n",
+ "fig = plt.figure()\n",
+ "ax = fig.add_subplot(111, projection='3d')\n",
+ "sc = ax.scatter(grid_coords.T[0], grid_coords.T[1], grid_coords.T[2], \\\n",
+ " c=weight_shifted, \\\n",
+ " norm=LogNorm(vmin=weight_shifted.min(), vmax=weight_shifted.max()))\n",
+ "fig.colorbar(sc)\n",
+ "# fig.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "关于变量 `grid_ao` 的意义，以及 `grid_weight` 的意义，将会在下面两小节中进行说明．\n",
+ "\n",
+ "刚才的格点 `grid_test` 数量显然太少，只有 $4 \\times 14 = 168$ 个点．在实际的量化计算中，显然这是不合适的．我们现在回到真实环境下所使用的格点 (99, 590)，来进行后续的计算．\n",
+ "\n",
+ "对于变量 `grid_mask`，猜测为因为在一些原子轨道积分 `grid_ao` 下，值总是为零或非常小，于是为了加速格点求和，忽略这些格点积分的数值．在这份笔记中，我们暂时不需要关心这些为加速运算速度而引入的数组．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ao_deriv = 1\n",
+ "grid_ao, grid_mask, grid_weight, grid_coords = \\\n",
+ " next(scf_eng._numint.block_loop(mol, scf_eng.grids, nao, ao_deriv, 2000))\n",
+ "print(\"Shape of arrays:\")\n",
+ "print(\"ao :\", grid_ao.shape)\n",
+ "print(\"mask :\", grid_mask.shape)\n",
+ "print(\"weight :\", grid_weight.shape)\n",
+ "print(\"coords :\", grid_coords.shape)\n",
+ "ngrids = grid_weight.shape[0]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "**注意**\n",
+ "\n",
+ "上面的代码中，我们对 `scf_eng._numint.block_loop` 函数使用的操作是将迭代过程中的一组数据取出．这也就意味着，实际上该函数是一个迭代函数，它完全可能生成不只一次格点信息．该函数传入的第 5 个参数是当前的内存信息；如果内存不足以存下足够的格点，那么格点将会多次生成．\n",
+ "\n",
+ "因此，在真正计算大分子或大格点时，不可以写成上面的样子；而必须要通过迭代器实现：\n",
+ "\n",
+ "```python\n",
+ "grid_ao, grid_mask, grid_weight, grid_coords = \\\n",
+ " next(scf_eng._numint.block_loop(mol, scf_eng.grids, nao, ao_deriv, 2000))\n",
+ " # Your manuplation on DFT grid comes here ...\n",
+ "```\n",
+ "\n",
+ "之所以在这里，我们可以只取一组数据就能进行 DFT 计算，是因为当前体系的分子足够小、格点足够小、基组也足够小，因此能在 2GB 左右的内存内进行给单积分计算．\n",
+ "\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 输出 1：电子数和格点求取，格点积分实现"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "首先，我们会使用格点积分处理最简单的电子数的和．这里的记号会混用 Einstein Convention 与积分记号，因此只会在引入格点积分时出现求和记号，其它时候不引入；符号上多少会复杂一些，但应当不会产生歧义．\n",
+ "\n",
+ "我们知道，电子密度的和可以用下述的方式给出：($D_{ij}$ 与 $D_{\\mu \\nu}$ 分别是 MO 基组与 AO 基组密度矩阵，它们两者不同但有下述的关系)\n",
+ "\n",
+ "\\begin{align}\n",
+ "n &= D_{ij} = 2 \\delta_{ij} = \\mathrm{tr} (\\mathbf{D}^\\mathrm{MO}) \\\\\n",
+ "&= C_{\\mu i} S_{\\mu \\nu} C_{\\nu j} = \\mathrm{tr} (\\mathbf{C}^\\mathrm{occ, T} \\mathbf{S} \\mathbf{C}^\\mathrm{occ}) = \\mathrm{tr} (\\mathbf{C}^\\mathrm{occ, T} \\mathbf{C}^\\mathrm{occ} \\mathbf{S}) = \\mathrm{tr} (\\mathbf{D}^\\mathrm{AO} \\mathbf{S}) \\\\\n",
+ "&= (C_{\\mu i} C_{\\nu j} \\delta_{ij}) S_{\\mu \\nu} = D_{\\mu \\nu} S_{\\mu \\nu}\n",
+ "\\end{align}\n",
+ "\n",
+ "下面我们只验证 $n = D_{\\mu \\nu} S_{\\mu \\nu}$："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose((D * S).sum(), nelec)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "上式完全是矩阵的运算．但如果我们希望使用格点积分，对空间变量 $\\boldsymbol{r}$ 显式地积分，那么上式可以根据重叠矩阵的定义写为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "n = D_{\\mu \\nu} \\int \\phi_{\\mu} (\\boldsymbol{r}) \\phi_{\\nu} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "\\end{equation}\n",
+ "\n",
+ "随后我们就可以进行格点积分计算了．格点积分的实现方式是对空间坐标积分量离散化并继而求和．具体地来说，如果现在的被积函数是 $f(\\boldsymbol{r})$，那么对于一系列已知坐标的、角标记为 $g$ 的格点 (即坐标集合) $\\{ \\boldsymbol{r}_g \\}_g$ 与格点权重 (即实数集合) $\\{ w_g \\}_g$，有\n",
+ "\n",
+ "\\begin{equation}\n",
+ "\\int f(\\boldsymbol{r}) \\mathrm{d} \\boldsymbol{r} \\simeq \\sum_{g} f(\\boldsymbol{r}_g) w_g\n",
+ "\\end{equation}\n",
+ "\n",
+ "如果现在将 $f(\\boldsymbol{r}_g)$ 记为 $f_g$，即将函数 $f(\\boldsymbol{r})$ 映射到格点集合 $\\{ f_g \\}_g$，那么我们可以使用 Einstein Convention 简记上述积分为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "\\int f(\\boldsymbol{r}) \\mathrm{d} \\boldsymbol{r} \\simeq f_g w_g\n",
+ "\\end{equation}\n",
+ "\n",
+ "回到电子数的格点积分．如果我们定义对原子轨道基组的格点映射为 $\\phi_{\\mu} (\\boldsymbol{r}) \\mapsto \\{ \\phi_{g \\mu} \\}_g$，那么电子数积分就可以写为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "n = D_{\\mu \\nu} \\phi_{g \\mu} \\phi_{g \\nu} w_g\n",
+ "\\end{equation}\n",
+ "\n",
+ "下面我们就通过格点积分验证电子数的和．这里指出，[上文中](#DFT-格点) `block_loop` 函数所声称的 `grid_ao` 是四部分构成的，其中的第一部分 `grid_ao[0]` 是原子轨道基组的格点 $\\phi_{g \\mu}$，而剩下三部分 `grid_ao[1:4]` 则是原子轨道基组的 $(x, y, z)$ 梯度分量 $\\nabla_t \\phi_{g \\mu}$，其中 $t \\in (x, y, z)$．在电子数求和时只需要 $\\phi_{g \\mu}$ 即可，而 $\\nabla_t \\phi_{g \\mu}$ 在交换相关势与能量求取是会使用．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "xc_n_sum = float(np.einsum(\"uv, gu, gv, g\", D, grid_ao[0], grid_ao[0], \\\n",
+ " grid_weight, optimize=True))\n",
+ "xc_n_sum"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose(xc_n_sum, xc_n)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "我们甚至可以顺便验证 $S_{\\mu \\nu} = \\phi_{g \\mu} \\phi_{g \\nu} w_g$，不过因为格点积分的精度不算太高，因此在使用 `np.allclose` 验证矩阵相同时，需要把默认精度 $10^{-8}$ 降低到 $10^{-7}$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose(np.einsum(\"gu, gv, g\", grid_ao[0], grid_ao[0], \\\n",
+ " grid_weight, optimize=True), \\\n",
+ " S, atol=1e-7)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 输出 2：交换相关能格点求取"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "交换相关能从定义上，其积分表达式为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "E^\\textrm{xc} [\\rho] = \\int f^\\textrm{xc} [\\rho, \\nabla \\rho] \\rho(\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "显然地，上式经过格点化后可以写为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "E^\\textrm{xc} [\\rho] = f^\\textrm{xc}_g \\rho_g w_g\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "看起来非常简洁．但是上式的 $f^\\textrm{xc}_g$ 与 $\\rho_g$ 都还不是已知量．$\\rho_g$ 的构造方法比较容易：\n",
+ "\n",
+ "\\begin{equation}\n",
+ "\\rho_g = \\rho (\\boldsymbol{r}_g) = D_{\\mu \\nu} \\phi_{\\mu} (\\boldsymbol{r}_g) \\phi_{\\nu} (\\boldsymbol{r}_g) = D_{\\mu \\nu} \\phi_{g \\mu} \\phi_{g \\nu}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "为了构造 $f^\\textrm{xc}_g$，我们需要传入关于 $\\rho, \\boldsymbol{\\nabla} \\rho$ 的信息，并且了解作为 DFT 近似核心的 $f^\\textrm{xc}_g$．PySCF 中，这些计算已经通过接口函数 `eval_xc` 完成，其底层一般来说是 LibXC 库函数．我们需要准备的信息不多，只需要下面三行代码即足够："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "grid_rho = np.einsum(\"uv, tgu, gv -> tg\", D, grid_ao, grid_ao[0], optimize=True)\n",
+ "grid_rho[1:4] *= 2\n",
+ "grid_exc, grid_vxc = scf_eng._numint.eval_xc(scf_eng.xc, grid_rho, deriv=1)[:2]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "关于上面三行代码，我们还有不少需要说明的地方．首先，关于第一行，实际上是合并了两个操作，同时导出了 $\\rho_g$ (`grid_rho[0]`) 与 $\\nabla_t \\rho_g$ (`grid_rho[1:4]`)．$\\rho_g$ 的导出方式上面已经叙述了；而 $\\nabla_t \\rho_g$ 的导出方式如下："
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\\begin{align}\n",
+ "\\nabla_t \\rho_g &= \\nabla_t \\rho (\\boldsymbol{r}_g) = D_{\\mu \\nu} \\nabla_t \\big( \\phi_{\\mu} (\\boldsymbol{r}_g) \\phi_{\\nu} (\\boldsymbol{r}_g) \\big) \\\\\n",
+ "&= D_{\\mu \\nu} \\nabla_t \\phi_{\\mu} (\\boldsymbol{r}_g) \\phi_{\\nu} (\\boldsymbol{r}_g) + D_{\\mu \\nu} \\phi_{\\mu} (\\boldsymbol{r}_g) \\nabla_t \\phi_{\\nu} (\\boldsymbol{r}_g) \\\\\n",
+ "&= D_{\\mu \\nu} (\\nabla_t \\phi_{g \\mu}) \\phi_{g \\nu} + D_{\\mu \\nu} \\phi_{g \\mu} (\\nabla_t \\phi_{g \\nu}) \\\\\n",
+ "&= 2 D_{\\mu \\nu} (\\nabla_t \\phi_{g \\mu}) \\phi_{g \\nu}\n",
+ "\\end{align}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "联系到 `grid_ao[1:4]` 储存了 $\\nabla_t \\phi_{g \\mu}$，而 `grid_ao[0]` 储存的是 $\\phi_{g \\mu}$，应当不难理解上述代码的第一行为何可以同时导出 $\\rho_g$ 和 $\\nabla_t \\rho_g$ 了．但是我们看到上面公式中有两倍，但 $\\nabla_t \\rho_g$ 公式中只有一倍，因此我们有必要对 `grid_ao[1:4]` 部分单独乘上 2 倍．"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "第三行代码则是具体执行计算 $f^\\mathrm{xc}_g$ 的代码．根据该代码的 [API 文档](https://sunqm.github.io/pyscf/dft.html#pyscf.dft.libxc.eval_xc)，由于我们不仅需要计算泛函的能量本身，还要计算其梯度以构造 Fock 矩阵，因此设置 `deriv=1`；由此，`eval_xc` 函数的输出尽管是一个四元素 tuple，但只有前两个元素非 `None`，且分别为泛函能量格点 `grid_exc` 与梯度格点 `grid_vxc`；后两个元素本来应当是泛函二阶梯度与三阶梯度量，但在 GGA 的 SCF 中不需要计算它们．关于梯度格点，将在下一小节中说明．"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "我们在构造 `grid_rho[0]` 后，我们还可以先简单验证其可以导出电子数之和对于水分子应当为 10：\n",
+ "\n",
+ "\\begin{equation}\n",
+ "n = \\rho_g w_g\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "(grid_rho[0] * grid_weight).sum()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "能量值的导出也相当容易：我们可以发现，`grid_exc` 的维度就是格点大小："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "grid_exc.shape"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "那么只要简单地套到能量格点公式 $E^\\textrm{xc} [\\rho] = f^\\textrm{xc}_g \\rho_g w_g$ 中，就可以给出交换相关能．我们验证这一结果．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "xc_e_sum = (grid_exc * grid_rho[0] * grid_weight).sum()\n",
+ "xc_e_sum"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose(xc_e_sum, xc_e)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 输出 3：交换相关势格点求取"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "交换相关势的求取公式仍然不算复杂．在这一小节，我们先列出公式与代码，再对它们作较为详细的解释；下一小节中，我们会简单地推导交换相关势的计算公式．\n",
+ "\n",
+ "交换相关势的计算公式为\n",
+ "\n",
+ "\\begin{align}\n",
+ "V^\\textrm{xc}_{\\mu \\nu} [\\rho] \n",
+ "&= w_g (\\partial_\\rho f^\\mathrm{xc}_{g}) \\phi_{g \\mu} \\phi_{g \\nu} \\\\\n",
+ "&\\quad + 2 w_g (\\partial_\\gamma f^\\mathrm{xc}_{g}) (\\nabla_t \\rho_g) (\\nabla_t \\phi_{g \\mu}) \\phi_{g \\nu} \\\\\n",
+ "&\\quad + 2 w_g (\\partial_\\gamma f^\\mathrm{xc}_{g}) (\\nabla_t \\rho_g) (\\nabla_t \\phi_{g \\nu}) \\phi_{g \\mu}\n",
+ "\\end{align}"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "grid_frho, grid_fgamma = grid_vxc[:2]\n",
+ "\n",
+ "xc_v_sum = np.einsum(\"g, g, gu, gv -> uv\", grid_weight, grid_frho, \\\n",
+ " grid_ao[0], grid_ao[0], optimize=True)\n",
+ "xc_v_sum += 2 * np.einsum(\"g, g, tg, tgu, gv -> uv\", grid_weight, grid_fgamma, \\\n",
+ " grid_rho[1:4], grid_ao[1:4], grid_ao[0], optimize=True)\n",
+ "xc_v_sum += 2 * np.einsum(\"g, g, tg, tgv, gu -> uv\", grid_weight, grid_fgamma, \\\n",
+ " grid_rho[1:4], grid_ao[1:4], grid_ao[0], optimize=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "我们可以验证上面导出的交换相关势 `xc_v_sum` 与使用高级函数 `nr_rks` 所导出的 `xc_v` 是相同的："
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "np.allclose(xc_v_sum, xc_v)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "上述代码中需要解释的是第一行．`grid_vxc` 是 `eval_xc` 函数所导出的第二个变量．我们知道 `eval_xc` 会返回四个变量，分别是关于能量格点、一阶导数格点、二阶导数格点、三阶导数格点．那么 `grid_vxc` 所表征的是一阶导数格点．\n",
+ "\n",
+ "尽管我们现在接触的是以 GGA 为底层的 B3LYP，但这里有必要简单了解 meta-GGA．meta-GGA 的泛函变量一般至多是下述四个数量值的几种：\n",
+ "\n",
+ "\\begin{equation}\n",
+ "\\rho, \\quad \n",
+ "\\gamma = \\boldsymbol{\\nabla} \\rho \\cdot \\boldsymbol{\\nabla} \\rho, \\quad\n",
+ "\\boldsymbol{\\nabla}^2 \\rho, \\quad\n",
+ "\\tau = | \\boldsymbol{\\nabla} \\phi_\\mu |^2\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "因此，在程序导出泛函的梯度时，也会相应地导出\n",
+ "\n",
+ "\\begin{equation}\n",
+ "\\left(\n",
+ "\\frac{\\partial f^\\mathrm{xc}}{\\partial \\rho}, \n",
+ "\\frac{\\partial f^\\mathrm{xc}}{\\partial \\gamma}, \n",
+ "\\frac{\\partial f^\\mathrm{xc}}{\\partial (\\boldsymbol{\\nabla}^2 \\rho)}, \n",
+ "\\frac{\\partial f^\\mathrm{xc}}{\\partial \\tau}\n",
+ "\\right)\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "上面的写法不适合得到文档的行内表示，因此后面会简记上面的记号为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "(\n",
+ "\\partial_\\rho f^\\mathrm{xc},\n",
+ "\\partial_\\gamma f^\\mathrm{xc},\n",
+ "\\partial_{\\boldsymbol{\\nabla}^2 \\rho} f^\\mathrm{xc},\n",
+ "\\partial_\\tau f^\\mathrm{xc}\n",
+ ")\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "上面写出的函数偏导实际上是关于空间变量 $\\boldsymbol{r}$ 的函数，只是因为记号复杂而将其省略．如果我们写得详细一些，譬如 GGA Kernel 下的 $\\partial_\\rho f^\\mathrm{xc}$，应该表示为"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "(\\partial_\\rho f^\\mathrm{xc}) (\\boldsymbol{r}) = \\frac{\\partial f^\\mathrm{xc} (\\rho, \\gamma)}{\\partial \\rho (\\boldsymbol{r})}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "当要化为格点时，即取 $\\boldsymbol{r}$ 为确定的坐标 $\\boldsymbol{r}_g$ 时，简记下述记号\n",
+ "\n",
+ "\\begin{equation}\n",
+ "\\partial_\\rho f^\\mathrm{xc}_g = (\\partial_\\rho f^\\mathrm{xc}) (\\boldsymbol{r}_g)\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "对于 $\\partial_\\gamma f^\\mathrm{xc}_g$ 亦如此．这样，我们就解释了上面代码所对应的公式的记号了．"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "对于 RKS 而言，`grid_vxc` 变量确实就是四维 tuple 变量；因此，第一行代码\n",
+ "\n",
+ "```\n",
+ "grid_frho, grid_fgamma = grid_vxc[:2]\n",
+ "```\n",
+ "\n",
+ "所作的就是将 $\\partial_\\rho f^\\mathrm{xc}_g$ 赋值给 `grid_frho`；$\\partial_\\gamma f^\\mathrm{xc}_g$ 赋值给 `grid_fgamma`；这些变量的维度都是格点大小．由于 B3LYP 是 GGA，因此 $\\partial_{\\boldsymbol{\\nabla}^2 \\rho} f^\\mathrm{xc}, \\partial_\\tau f^\\mathrm{xc}$ 没有意义，即 `grid_vxc[2:4]` 不应是有意义的量．下面的代码可以简单说明这些．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(grid_vxc.__len__())\n",
+ "print(grid_vxc[2:4])\n",
+ "print(grid_frho.shape)\n",
+ "print(grid_fgamma.shape)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "\n",
+ "**提示**\n",
+ "\n",
+ "但对于 UKS 而言，由于 $\\alpha$ 与 $\\beta$ 密度不同，因此 `eval_xc` 给出的一阶导数格点尽管仍然是四维 tuple 变量，但矩阵的数量总共有 9 个．更多信息需要参考 [API 文档](https://sunqm.github.io/pyscf/dft.html#pyscf.dft.libxc.eval_xc)．\n",
+ "\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 交换相关势的推导"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "这里的推导仅仅是表明交换相关势中每一项、以及它们的系数是如何导出的；更为详细与严谨的推导还需要参考其它文献．"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "泛函的变分一般来说只有放在积分的环境下才具有意义．根据 [Wikipedia](https://en.wikipedia.org/wiki/Functional_derivative#Formula) 的说明，如果能量泛函记为"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "E^\\textrm{xc} [\\rho] = \\int f^\\textrm{xc} (\\rho, \\gamma) \\rho(\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "我们可以对 GGA 型的泛函一次变分记为"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\\begin{align}\n",
+ "\\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\rho (\\boldsymbol{r})} \\phi (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r} \n",
+ "&= \\int \\left[ \\frac{\\partial f^\\textrm{xc} (\\rho, \\gamma)}{\\partial \\rho (\\boldsymbol{r})} \\phi (\\boldsymbol{r})\n",
+ "+ \\frac{\\partial f^\\textrm{xc} (\\rho, \\gamma)}{\\partial \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r})} \\cdot \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\phi (\\boldsymbol{r}) \\right] \n",
+ "\\, \\mathrm{d} \\boldsymbol{r} \\\\\n",
+ "&= \\int \\left[ \\frac{\\partial f^\\textrm{xc} (\\rho, \\gamma)}{\\partial \\rho (\\boldsymbol{r})} \\phi (\\boldsymbol{r})\n",
+ "+ \\frac{\\partial f^\\textrm{xc} (\\rho, \\gamma)}{\\partial \\gamma (\\boldsymbol{r})} 2 \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r}) \\cdot \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\phi (\\boldsymbol{r}) \\right] \n",
+ "\\, \\mathrm{d} \\boldsymbol{r}\n",
+ "\\end{align}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "第一个等式是变分本身的定义 (需要利用广义 Stocks 定理)；第二个等式利用到连续偏导以及关于 $\\gamma$ 与 $\\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho$ 关系的一个性质："
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\frac{\\partial \\gamma (\\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho)}{\\partial \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r})}\n",
+ "= \\frac{\\partial (\\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho \\cdot \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho)}{\\partial \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r})}\n",
+ "= 2 \\boldsymbol{\\nabla}_{\\boldsymbol{r}'} \\rho (\\boldsymbol{r})\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "上面出现记号 $\\boldsymbol{\\nabla}_{\\boldsymbol{r}'}$；这仅仅是普通的梯度记号，一般会记为 $\\nabla$；这里记得较为复杂，只是强调梯度的坐标变量 $\\boldsymbol{r}'$ 与被积元变量 $\\boldsymbol{r}$ 不是相同的．如果我们使用 Einstein Convention，对 $\\boldsymbol{r}'$ 更换为 $t \\in (x, y, z)$ 并在实际计算过程中对 $t$ 求和，来记 GGA 的一次变分；同时按照 [上一小节](#输出-3：交换相关势格点求取) 简化偏导数记号；那么可以写为"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\rho (\\boldsymbol{r})} \\phi (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "= \\int \\left[ (\\partial_{\\rho} f^\\mathrm{xc}) \\phi + 2 (\\partial_{\\gamma} f^\\mathrm{xc}) (\\nabla_t \\rho) (\\nabla_t \\phi) \\right] \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "现在我们考察上述结论是如何应用到 $V^\\textrm{xc}_{\\mu \\nu} [\\rho]$ 的导出过程中的．我们知道，$V^\\textrm{xc}_{\\mu \\nu} [\\rho]$ 是 Fock 矩阵构成的一部分；Fock 矩阵的构成方式是将 Hartree-Fock 能量对分子轨道进行变分而产生的．能量可以分为单电子、库伦积分、交换积分与泛函积分的贡献，那么 Fock 矩阵也会分为这四部分；关于这点我们已经在 [上文](#输出-3：交换相关势) 中验证过了．那么泛函积分对 Fock 矩阵的贡献可以表示为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "F_{\\mu \\nu} [\\rho] \\leftarrow V^\\textrm{xc}_{\\mu \\nu} [\\rho] = \\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\phi_{\\mu} (\\boldsymbol{r})} \\phi_{\\nu} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "由于密度的表达式可以写为 (该公式不使用 Einstein Convention)\n",
+ "\n",
+ "\\begin{equation}\n",
+ "\\rho(\\boldsymbol{r}) = \\rho [\\{ \\phi_{\\mu} \\}_\\mu] (\\boldsymbol{r}) = \\int \\sum_{\\nu} \\phi_{\\nu} (\\boldsymbol{r}) \\phi_{\\nu} (\\boldsymbol{r}') \\delta(\\boldsymbol{r} - \\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r}'\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "因此密度对轨道的变分可以写作 (该公式不使用 Einstein Convention)\n",
+ "\n",
+ "\\begin{align}\n",
+ "\\frac{\\delta \\rho}{\\delta \\phi_\\mu (\\boldsymbol{r})} \n",
+ "&= \\int \\sum_{\\nu} \\frac{\\partial \\phi_{\\nu}}{\\partial \\phi_{\\mu} (\\boldsymbol{r})} \\phi_{\\nu} (\\boldsymbol{r}') \\delta(\\boldsymbol{r} - \\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r}' \\\\\n",
+ "&= \\int \\phi_{\\mu} (\\boldsymbol{r}') \\delta(\\boldsymbol{r} - \\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r}' = \\phi_\\mu (\\boldsymbol{r})\n",
+ "\\end{align}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "因此，根据连续变分规则，有\n",
+ "\\begin{equation}\n",
+ "V^\\textrm{xc}_{\\mu \\nu} [\\rho] = \\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\rho (\\boldsymbol{r})} \\frac{\\delta \\rho}{\\delta \\phi_\\mu (\\boldsymbol{r})} \\phi_{\\nu} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "= \\int \\frac{\\delta E^\\textrm{xc} [\\rho]}{\\delta \\rho (\\boldsymbol{r})} \\phi_\\mu (\\boldsymbol{r}) \\phi_{\\nu} (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "那么我们用上面已经有的结论，令 $\\phi (\\boldsymbol{r}) = \\phi_\\mu (\\boldsymbol{r}) \\phi_{\\nu} (\\boldsymbol{r})$，并注意到"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\\begin{equation}\n",
+ "\\nabla_t \\phi = \\nabla_t (\\phi_\\mu \\phi_\\nu) = (\\nabla_t \\phi_\\mu) \\phi_\\nu + (\\nabla_t \\phi_\\nu) \\phi_\\mu\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "那么\n",
+ "\n",
+ "\\begin{equation}\n",
+ "V^\\textrm{xc}_{\\mu \\nu} [\\rho] = \\int \\left[ (\\partial_{\\rho} f^\\mathrm{xc}) \\phi_\\mu \\phi_\\nu + 2 (\\partial_{\\gamma} f^\\mathrm{xc}) (\\nabla_t \\rho) (\\nabla_t \\phi_\\mu) \\phi_\\nu + 2 (\\partial_{\\gamma} f^\\mathrm{xc}) (\\nabla_t \\rho) (\\nabla_t \\phi_\\nu) \\phi_\\mu \\right] \\, \\mathrm{d} \\boldsymbol{r}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "将上式格点化，就立即得到 [上一小节](#输出-3：交换相关势格点求取) 所使用的格点公式了．"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## XYG3 能量"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### 手动设置 B3LYP 泛函形式"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "XYG3 作为非自洽泛函，其密度与轨道从 B3LYP 获得，但其非自洽能量泛函则使用自己的泛函形式．如果要在尚未实现 XYG3 的软件中获得 XYG3 能量，就必须手动设置泛函参数．在实际使用 XYG3 参数前，我们先尝试对 B3LYP 进行手动的参数设置，并验证结果是否与上面的计算相同．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "b3try_eng = dft.RKS(mol)\n",
+ "b3try_eng.xc = \"HF*0.2 + .08*LDA + .72*B88, .81*LYP + .19*VWN3\" # compare that to gaussian\n",
+ "b3try_eng.grids.atom_grid = (99, 590)\n",
+ "b3try_eng.grids.build()\n",
+ "\n",
+ "b3try_eng.kernel()\n",
+ "print(\"Compare My B3LYPG: \", np.allclose(b3try_eng.e_tot, scf_eng.e_tot))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### XYG3 泛函形式与能量"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "XYG3 型双杂化泛函属于第五阶泛函 (按 Predew 的 [Jacob 阶梯](https://dx.doi.org/10.1063/1.1390175) 说法)，其 [定义](https://dx.doi.org/10.1073/pnas.0901093106) 是\n",
+ "\n",
+ "\\begin{equation}\n",
+ "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}\n",
+ "\\end{equation}\n",
+ "\n",
+ "对应到程序中，每一项的系数则展开为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "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}\n",
+ "\\end{equation}\n",
+ "\n",
+ "对于 XYG3，其系数的确定是\n",
+ "\n",
+ "\\begin{equation}\n",
+ "c_1 = 0.8033, \\quad c_2 = 0.2107, \\quad c_3 = 0.3211, \\quad c_4 = 1 - c_3\n",
+ "\\end{equation}\n",
+ "\n",
+ "因此，XYG3 的泛函形式为\n",
+ "\n",
+ "\\begin{equation}\n",
+ "E_\\mathrm{xc}^\\textsf{XYG3} = 0.8033 E_\\mathrm{x}^\\textsf{exact} - 0.0140 E_\\mathrm{x}^\\textsf{LDA} + 0.2107 E_\\mathrm{x}^\\textsf{GGA} + 0.6789 E_\\mathrm{c}^\\textsf{GGA} + 0.3211 E_\\mathrm{c}^\\textsf{PT2}\n",
+ "\\end{equation}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "在 PySCF 中，泛函形式的确定通过字符串进行定义．该字符串的定义方式可以很灵活，下面的代码只是其中一种定义方式．其中，逗号前后分别分割了交换能与相关能的部分．由于 PySCF 中不对 PT2 型相关能在 DFT 部分进行定义，因此在定义泛函的时候，不将 $0.3211 E_\\mathrm{c}^\\textsf{PT2}$ 写入字符串．\n",
+ "\n",
+ "当泛函的形式确定后，我们只要简单地把自洽场过程中收敛的 B3LYP 密度代入 XYG3 型泛函能量表达式中，即可得到最终的 XYG3 能量了．"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "nscf_eng = dft.RKS(mol)\n",
+ "nscf_eng.xc = \"0.8033*HF - 0.0140*LDA + 0.2107*B88, 0.6789*LYP\"\n",
+ "nscf_eng.grids.atom_grid = (99, 590)\n",
+ "nscf_eng.grids.build()\n",
+ "\n",
+ "xyg3_energy = nscf_eng.energy_tot(dm=D) + mp2_eng.e_corr * 0.3211\n",
+ "print(\"XYG3 energy allclose: \", np.allclose(xyg3_energy, -0.76282393305943E+02))"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.6.6"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}