{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# GGA 方法 U 矩阵与 B2PLYP 型泛函核坐标梯度" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们进一步讨论 GGA 方法的 U 矩阵求取方式。在求取 U 矩阵的过程中，我们利用到 GGA 方法下的 CP-HF 方程，进而可以用来求取含有 GGA 的 B2PLYP 型泛函的弛豫密度矩阵 $D_{pq}^\\mathrm{MP2}$，进而得到 B2PLYP 型泛函核坐标梯度。\n", "\n", "在 pyxdh 中，B2PLYP 型泛函与 MP2 方法使用相同的流程进行处理；这两者相当相似。这一节的大多数公式和做法都可以参考 MP2 梯度部分；但在一些公式细节上，需要留意细微的差别。\n", "\n", "这一节中提到的 U 矩阵求取仅仅包括非占-占据部分 $U_{ai}^{A_t}$。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 准备工作" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "计算 B2PLYP 型泛函的梯度，我们用与 MP2 方法一样的 `GradMP2` 类。由于涉及到大批格点的计算，我们这里暂时假定供给 np.einsum 的内存大小是足够大的。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook\n", "\n", "from pyscf import gto, scf, dft, dft, lib, hessian\n", "from pyscf.scf import cphf\n", "import numpy as np\n", "from functools import partial\n", "import warnings\n", "from matplotlib import pyplot as plt\n", "from pyxdh.Utilities import NucCoordDerivGenerator, DipoleDerivGenerator, NumericDiff, GridHelper, KernelHelper\n", "from pyxdh.DerivOnce import GradMP2, GradSCF\n", "\n", "np.einsum = partial(np.einsum, optimize=[\"greedy\", 1024 ** 3 * 1000 / 8])\n", "np.einsum_path = partial(np.einsum_path, optimize=[\"greedy\", 1024 ** 3 * 1000 / 8])\n", "np.allclose = partial(np.allclose, atol=1e-6, rtol=1e-4)\n", "np.set_printoptions(5, linewidth=150, suppress=True)\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mol = gto.Mole()\n", "mol.atom = \"\"\"\n", "O 0.0 0.0 0.0\n", "O 0.0 0.0 1.5\n", "H 1.0 0.0 0.0\n", "H 0.0 0.7 1.0\n", "\"\"\"\n", "mol.basis = \"6-31G\"\n", "mol.verbose = 0\n", "mol.build()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def mol_to_grids(mol, atom_grid=(75, 302)):\n", " grids = dft.Grids(mol)\n", " grids.atom_grid = atom_grid\n", " grids.becke_scheme = dft.gen_grid.stratmann\n", " grids.prune = None\n", " grids.build()\n", " return grids\n", "grids = mol_to_grids(mol)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "但需要注意到，既然我们计算的是 B2PLYP，那么其泛函也需要更改为 B2PLYP 的自洽场过程中使用的泛函：" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def mol_to_scf(mol):\n", " scf_eng = dft.RKS(mol)\n", " scf_eng.grids = mol_to_grids(mol)\n", " scf_eng.xc = \"0.53*HF + 0.47*B88, 0.73*LYP\"\n", " scf_eng.conv_tol = 1e-10\n", " return scf_eng.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "与 MP2 方法不同的是，我们需要额外定义 PT2 相关能系数 `cc` $c_\\mathrm{c}$：" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "gradh = GradMP2({\"scf_eng\": mol_to_scf(mol), \"cc\": 0.27})\n", "cc = gradh.cc" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "nmo, nao, natm, nocc, nvir, cx = gradh.nao, gradh.nao, gradh.natm, gradh.nocc, gradh.nvir, gradh.cx\n", "mol_slice = gradh.mol_slice\n", "so, sv, sa = gradh.so, gradh.sv, gradh.sa\n", "C, Co, Cv, e, eo, ev, D = gradh.C, gradh.Co, gradh.Cv, gradh.e, gradh.eo, gradh.ev, gradh.D\n", "H_0_ao, S_0_ao, eri0_ao, F_0_ao = gradh.H_0_ao, gradh.S_0_ao, gradh.eri0_ao, gradh.F_0_ao\n", "H_0_mo, S_0_mo, eri0_mo, F_0_mo = gradh.H_0_mo, gradh.S_0_mo, gradh.eri0_mo, gradh.F_0_mo" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def to_natm_3(mat: np.ndarray):\n", " shape = list(mat.shape)\n", " shape = [int(shape[0] / 3), 3] + shape[1:]\n", " return mat.reshape(shape)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "H_1_ao, S_1_ao, eri1_ao = to_natm_3(gradh.H_1_ao), to_natm_3(gradh.S_1_ao), to_natm_3(gradh.eri1_ao)\n", "H_1_mo, S_1_mo, eri1_mo = to_natm_3(gradh.H_1_mo), to_natm_3(gradh.S_1_mo), to_natm_3(gradh.eri1_mo)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "grdh = GridHelper(mol, grids, D)\n", "kerh = KernelHelper(grdh, \"0.53*HF + 0.47*B88, 0.73*LYP\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "为了简化后面程序代码，我们将一部分与格点相关的 (基于 $D_{\\mu \\nu}$ 自洽场密度给出的) 变量定义如下：" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "ao_0, ao_1, ao_2 = grdh.ao_0, grdh.ao_1, grdh.ao_2\n", "rho_1, rho_2 = grdh.rho_1, grdh.rho_2\n", "A_rho_1, A_rho_2 = grdh.A_rho_1, grdh.A_rho_2\n", "fr, fg, frr, frg, fgg = kerh.fr, kerh.fg, kerh.frr, kerh.frg, kerh.fgg" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def grad_generator(mol):\n", " scf_eng = mol_to_scf(mol)\n", " config = {\"scf_eng\": mol_to_scf(mol), \"cc\": 0.27}\n", " return GradMP2(config)\n", "gradn = NucCoordDerivGenerator(mol, grad_generator)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## GGA 方法 U 矩阵的求取" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fock Skeleton 导数" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "求取 U 矩阵 (非占-占据部分) 的推导前提是 Fock 矩阵 $F_{ai} = 0$。我们首先回顾原子轨道下，GGA 的 Fock 矩阵定义：\n", "\n", "$$\n", "F_{\\mu \\nu} = h_{\\mu \\nu} + (\\mu \\nu | \\kappa \\lambda) D_{\\kappa \\lambda} - \\frac{c_\\mathrm{x}}{2} (\\mu \\kappa | \\nu \\lambda) D_{\\kappa \\lambda} + v_{\\mu \\nu}^\\mathrm{xc}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "前三项我们已经在 RHF 的 U 矩阵计算时知道了求取方法；现在我们要考虑第四项 $v_{\\mu \\nu}^\\mathrm{xc}$ 导数的计算。我们再回顾 $v_{\\mu \\nu}^{\\mathrm{xc}}$ 的定义：\n", "\n", "$$\n", "v_{\\mu \\nu}^\\mathrm{xc} = f_\\rho \\phi_\\mu \\phi_\\nu + 2 f_\\gamma \\rho_r (\\phi_{r \\mu} \\phi_{\\nu} + \\phi_{\\mu} \\phi_{r \\nu})\n", "$$\n", "\n", "我们对 Skeleton 导数的定义是不包含 $C_{\\mu p}$ 的导数；剩下的部分则归为 U 矩阵导数。我们将 $v_{\\mu \\nu}^{\\mathrm{xc}}$ 的 Skeleton 导数写为 $v_{\\mu \\nu}^{\\mathrm{xc}, A_t}$。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们先依据链式法则，考察 $f_\\rho$ 的全导数：\n", "\n", "$$\n", "\\frac{\\partial f_\\rho}{\\partial A_t} = f_{\\rho \\rho} \\frac{\\partial \\rho}{\\partial A_t} + 2 f_{\\rho \\gamma} \\rho_r \\frac{\\partial \\rho_r}{\\partial A_t}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们已经对 $\\partial_{A_t} \\rho$ 与 $\\partial_{A_t} \\rho_r$ 作过讨论：\n", "\n", "$$\n", "\\begin{align}\n", "\\frac{\\partial \\rho}{\\partial A_t}\n", "&= - \\phi_{t \\mu_A} \\phi_\\nu D_{\\mu \\nu} - \\phi_\\mu \\phi_{t \\nu_A} D_{\\mu \\nu} + 4 \\phi_\\mu \\phi_\\nu U_{mi}^{A_t} C_{\\mu m} C_{\\nu i} \\\\\n", "&= \\rho^{A_t} + 4 \\phi_\\mu \\phi_\\nu U_{mi}^{A_t} C_{\\mu m} C_{\\nu i} \\\\\n", "\\frac{\\partial \\rho_r}{\\partial A_t}\n", "&= - 2 \\phi_{tr \\mu_A} \\phi_\\nu D_{\\mu \\nu} - 2 \\phi_{r \\mu} \\phi_{t \\nu_A} D_{\\mu \\nu} + 4 (\\phi_{r \\mu} \\phi_\\nu + \\phi_\\mu \\phi_{r \\nu}) U_{mi}^{A_t} C_{\\mu m} C_{\\nu i} \\\\\n", "&= \\rho_r^{A_t} + 4 (\\phi_{r \\mu} \\phi_\\nu + \\phi_\\mu \\phi_{r \\nu}) U_{mi}^{A_t} C_{\\mu m} C_{\\nu i}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "对于 $f_\\gamma$ 的全导数也是相同的：\n", "\n", "$$\n", "\\frac{\\partial f_\\gamma}{\\partial A_t} = f_{\\rho \\gamma} \\frac{\\partial \\rho}{\\partial A_t} + 2 f_{\\gamma \\gamma} \\rho_r \\frac{\\partial \\rho_r}{\\partial A_t}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "上述出现了 U 矩阵的部分看作是 U 导数，否则为 Skeleton 导数。那么，$v_{\\mu \\nu}^\\mathrm{xc}$ 的所有贡献项就容易地用链式法则表示如下：\n", "\n", "$$\n", "\\begin{align}\n", "\\partial_{A_t} v_{\\mu \\nu}^\\mathrm{xc} \\xleftarrow{\\text{Skeleton derivative}} v_{\\mu \\nu}^{\\mathrm{xc}, A_t}\n", "&= (f_{\\rho \\rho} \\rho^{A_t} + 2 f_{\\rho \\gamma} \\rho_w \\rho_w^{A_t}) \\phi_\\mu \\phi_\\nu + 2 (f_{\\rho \\gamma} \\rho^{A_t} + 2 f_{\\gamma \\gamma} \\rho_w \\rho_w^{A_t}) \\rho_r \\rho_r^{A_t} (\\phi_{r \\mu} \\phi_{\\nu} + \\phi_{\\mu} \\phi_{r \\nu}) \\\\\n", "&\\quad + f_\\rho \\phi_\\mu^{A_t} \\phi_\\nu + f_\\rho \\phi_\\mu \\phi_\\nu^{A_t} \\\\\n", "&\\quad + 2 f_\\gamma \\rho_r^{A_t} (\\phi_{r \\mu} \\phi_{\\nu} + \\phi_{\\mu} \\phi_{r \\nu})\n", "+ 2 f_\\gamma \\rho_r (\\phi_{r \\mu}^{A_t} \\phi_{\\nu} + \\phi_{\\mu}^{A_t} \\phi_{r \\nu} + \\phi_{r \\mu} \\phi_{\\nu}^{A_t} + \\phi_{\\mu} \\phi_{r \\nu}^{A_t})\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "上述出现了 U 矩阵的部分看作是 U 导数，否则为 Skeleton 导数。那么，$v_{\\mu \\nu}^\\mathrm{xc}$ 的所有贡献项就容易地用链式法则表示如下：\n", "\n", "$$\n", "\\begin{align}\n", "\\partial_{A_t} v_{\\mu \\nu}^\\mathrm{xc} \\xleftarrow{\\text{Skeleton derivative}} v_{\\mu \\nu}^{\\mathrm{xc}, A_t}\n", "&= (f_{\\rho \\rho} \\rho^{A_t} + 2 f_{\\rho \\gamma} \\rho_w \\rho_w^{A_t}) \\phi_\\mu \\phi_\\nu \\\\\n", "&\\quad + 2 (f_{\\rho \\gamma} \\rho^{A_t} + 2 f_{\\gamma \\gamma} \\rho_w \\rho_w^{A_t}) \\rho_r (\\phi_{r \\mu} \\phi_{\\nu} + \\phi_{\\mu} \\phi_{r \\nu}) \\\\\n", "&\\quad + f_\\rho \\phi_\\mu^{A_t} \\phi_\\nu + f_\\rho \\phi_\\mu \\phi_\\nu^{A_t} \\\\\n", "&\\quad + 2 f_\\gamma \\rho_r^{A_t} (\\phi_{r \\mu} \\phi_{\\nu} + \\phi_{\\mu} \\phi_{r \\nu})\n", "+ 2 f_\\gamma \\rho_r (\\phi_{r \\mu}^{A_t} \\phi_{\\nu} + \\phi_{\\mu}^{A_t} \\phi_{r \\nu} + \\phi_{r \\mu} \\phi_{\\nu}^{A_t} + \\phi_{\\mu} \\phi_{r \\nu}^{A_t}) \\\\\n", "&= \\frac{1}{2} f_{\\rho \\rho} \\rho^{A_t} \\phi_\\mu \\phi_\\nu + f_{\\rho \\gamma} \\rho_w \\rho_w^{A_t} \\phi_\\mu \\phi_\\nu \\\\\n", "&\\quad + 2 f_{\\rho \\gamma} \\rho^{A_t} \\rho_r \\phi_{r \\mu} \\phi_{\\nu} + 4 f_{\\gamma \\gamma} \\rho_w \\rho_w^{A_t} \\rho_r \\phi_{r \\mu} \\phi_{\\nu} \\\\\n", "&\\quad + f_\\rho \\phi_\\mu^{A_t} \\phi_\\nu \\\\\n", "&\\quad + 2 f_\\gamma \\rho_r^{A_t} \\phi_{r \\mu} \\phi_{\\nu} + 2 f_\\gamma \\rho_r \\phi_{r \\mu}^{A_t} \\phi_{\\nu} + 2 f_\\gamma \\rho_r \\phi_{\\mu}^{A_t} \\phi_{r \\nu} \\\\\n", "&\\quad + \\mathrm{swap} (\\mu, \\nu)\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "利用 $\\partial_{A_t} \\phi_\\mu = - \\phi_{t \\mu_A}$ 等结论，我们会将上式写为\n", "\n", "$$\n", "\\begin{align}\n", "\\partial_{A_t} v_{\\mu \\nu}^\\mathrm{xc} \\xleftarrow{\\text{Skeleton derivative}} v_{\\mu \\nu}^{\\mathrm{xc}, A_t}\n", "&= \\frac{1}{2} f_{\\rho \\rho} \\rho^{A_t} \\phi_\\mu \\phi_\\nu + f_{\\rho \\gamma} \\rho_w \\rho_w^{A_t} \\phi_\\mu \\phi_\\nu \\\\\n", "&\\quad + 2 f_{\\rho \\gamma} \\rho^{A_t} \\rho_r \\phi_{r \\mu} \\phi_{\\nu} + 4 f_{\\gamma \\gamma} \\rho_w \\rho_w^{A_t} \\rho_r \\phi_{r \\mu} \\phi_{\\nu} \\\\\n", "&\\quad + 2 f_\\gamma \\rho_r^{A_t} \\phi_{r \\mu} \\phi_{\\nu} \\\\\n", "&\\quad - f_\\rho \\phi_{t \\mu_A} \\phi_\\nu - 2 f_\\gamma \\rho_r \\phi_{tr \\mu_A} \\phi_{\\nu} - 2 f_\\gamma \\rho_r \\phi_{t \\mu_A} \\phi_{r \\nu} \\\\\n", "&\\quad + \\mathrm{swap} (\\mu, \\nu)\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们定义 `F_1_ao_GGA` $v_{\\mu \\nu}^{\\mathrm{xc}, A_t}$ 为 (维度为 $(A, t, \\mu, \\nu)$)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "(4, 3, 22, 22)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_1_ao_GGA = (\n", " + 0.5 * np.einsum(\"g, Atg, gu, gv -> Atuv\", frr, A_rho_1, ao_0, ao_0)\n", " + np.einsum(\"g, wg, Atwg, gu, gv -> Atuv\", frg, rho_1, A_rho_2, ao_0, ao_0)\n", " + 2 * np.einsum(\"g, Atg, rg, rgu, gv -> Atuv\", frg, A_rho_1, rho_1, ao_1, ao_0)\n", " + 4 * np.einsum(\"g, wg, Atwg, rg, rgu, gv -> Atuv\", fgg, rho_1, A_rho_2, rho_1, ao_1, ao_0)\n", " + 2 * np.einsum(\"g, Atrg, rgu, gv -> Atuv\", fg, A_rho_2, ao_1, ao_0)\n", ")\n", "for A in range(natm):\n", " sA = mol_slice(A)\n", " F_1_ao_GGA[A, :, sA, :] -= np.einsum(\"g, tgu, gv -> tuv\", fr, ao_1[:, :, sA], ao_0)\n", " F_1_ao_GGA[A, :, sA, :] -= 2 * np.einsum(\"g, rg, trgu, gv -> tuv\", fg, rho_1, ao_2[:, :, :, sA], ao_0)\n", " F_1_ao_GGA[A, :, sA, :] -= 2 * np.einsum(\"g, rg, tgu, rgv -> tuv\", fg, rho_1, ao_1[:, :, sA], ao_1)\n", "F_1_ao_GGA += F_1_ao_GGA.swapaxes(-1, -2) # swap (u, v)\n", "F_1_ao_GGA.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们再考虑 Fock 矩阵的 Skeleton 导数 `F_1_ao`：\n", "\n", "$$\n", "F_{\\mu \\nu}^{A_t} = h_{\\mu \\nu}^{A_t} + (\\mu \\nu | \\kappa \\lambda)^{A_t} D_{\\kappa \\lambda} - \\frac{c_\\mathrm{x}}{2} (\\mu \\kappa | \\nu \\lambda)^{A_t} D_{\\kappa \\lambda} + v_{\\mu \\nu}^{\\mathrm{xc}, A_t}\n", "$$" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(4, 3, 22, 22)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F_1_ao = (\n", " + H_1_ao\n", " + np.einsum(\"Atuvkl, kl -> Atuv\", eri1_ao, D)\n", " - 0.5 * cx * np.einsum(\"Atukvl, kl -> Atuv\", eri1_ao, D)\n", " + F_1_ao_GGA\n", ")\n", "F_1_ao.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "pyxdh 中的 `F_1_ao` 也可以用来验证上述的生成结果：" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(F_1_ao, to_natm_3(gradh.F_1_ao))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "在 PySCF 中，也可以用 `make_h1` 来验证上述结果：" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(F_1_ao, hessian.rks.Hessian(gradh.scf_eng).make_h1(C, gradh.mo_occ))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "那么分子轨道下的 $F_{pq}^{A_t}$ 可以表示为：\n", "\n", "$$\n", "F_{pq}^{A_t} = C_{\\mu p} F_{\\mu \\nu}^{A_t} C_{\\nu q}\n", "$$" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "F_1_mo = np.einsum(\"up, Atuv, vq -> Atpq\", C, F_1_ao, C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### A 张量" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们已经对 $\\partial_{A_t} \\rho$ 与 $\\partial_{A_t} \\rho_r$ 作过讨论：\n", "\n", "$$\n", "\\begin{align}\n", "\\frac{\\partial \\rho}{\\partial A_t}\n", "&= \\rho^{A_t} + 4 \\phi_\\mu \\phi_\\nu U_{mi}^{A_t} C_{\\mu m} C_{\\nu i} \\\\\n", "\\frac{\\partial \\rho_r}{\\partial A_t}\n", "&= \\rho_r^{A_t} + 4 (\\phi_{r \\mu} \\phi_\\nu + \\phi_\\mu \\phi_{r \\nu}) U_{mi}^{A_t} C_{\\mu m} C_{\\nu i}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "从 RHF 推导 A 张量的经验，我们知道，A 张量的定义应是 Fock 矩阵的 U 导数所产生的贡献：\n", "\n", "$$\n", "\\frac{\\partial F_{\\mu \\nu}}{\\partial A_t} \\xleftarrow{\\text{U derivative}} A_{\\mu \\nu, \\kappa \\lambda} C_{\\kappa m} C_{\\lambda i} U_{mi}^{A_t}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们仍然根据定义\n", "\n", "$$\n", "v_{\\mu \\nu}^\\mathrm{xc} = f_\\rho \\phi_\\mu \\phi_\\nu + 2 f_\\gamma \\rho_r (\\phi_{r \\mu} \\phi_{\\nu} + \\phi_{\\mu} \\phi_{r \\nu})\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "可以推知其 U 矩阵的贡献是\n", "\n", "$$\n", "\\begin{align}\n", "\\partial_{A_t} v_{\\mu \\nu}^\\mathrm{xc} \\xleftarrow{\\text{U derivative}}\n", "&\\quad \\big[ 4 (f_{\\rho \\rho} \\phi_\\kappa \\phi_\\lambda + 2 f_{\\rho \\gamma} \\rho_r \\phi_{r \\kappa} \\phi_\\lambda + 2 f_{\\rho \\gamma} \\rho_r \\phi_\\kappa \\phi_{r \\lambda}) \\phi_\\mu \\phi_\\nu \\\\\n", "&\\quad + 8 (f_{\\rho \\gamma} \\phi_\\kappa \\phi_\\lambda + 2 f_{\\gamma \\gamma} \\rho_w \\phi_{w \\kappa} \\phi_\\lambda + 2 f_{\\gamma \\gamma} \\rho_r \\phi_\\kappa \\phi_{r \\lambda}) \\rho_r (\\phi_{r \\mu} \\phi_{\\nu} + \\phi_{\\mu} \\phi_{r \\nu}) \\\\\n", "&\\quad + 8 f_\\gamma (\\phi_{r \\kappa} \\phi_\\lambda + \\phi_\\kappa \\phi_{r \\lambda}) (\\phi_{r \\mu} \\phi_{\\nu} + \\phi_{\\mu} \\phi_{r \\nu}) \\big] U_{mi}^{A_t} C_{\\mu m} C_{\\nu i}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "上式中，方括号所包含的部分就是 GGA 对 A 张量 $A_{\\mu \\nu, \\kappa \\lambda}$ 的贡献；但这个表达式实在是相当长，我们可以依靠对称性作一些简化：\n", "\n", "$$\n", "\\begin{align}\n", "A_{\\mu \\nu, \\kappa \\lambda} \\xleftarrow{\\text{GGA contrib}} A_{\\mu \\nu, \\kappa \\lambda}^\\text{GGA}\n", "&= f_{\\rho \\rho} \\phi_\\kappa \\phi_\\lambda \\phi_\\mu \\phi_\\nu + 4 f_{\\rho \\gamma} \\rho_r \\phi_{r \\kappa} \\phi_\\lambda \\phi_\\mu \\phi_\\nu \\\\\n", "& + 4 f_{\\rho \\gamma} \\phi_\\kappa \\phi_\\lambda \\rho_r \\phi_{r \\mu} \\phi_{\\nu} + 16 f_{\\gamma \\gamma} \\rho_w \\phi_{w \\kappa} \\phi_\\lambda \\rho_r \\phi_{r \\mu} \\phi_{\\nu} \\\\\n", "& + 8 f_\\gamma \\phi_{r \\kappa} \\phi_\\lambda \\phi_{r \\mu} \\phi_{\\nu} \\\\\n", "& + \\mathrm{swap} (\\mu, \\nu) + \\mathrm{swap} (\\kappa, \\lambda)\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(22, 22, 22, 22)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A_0_ao_GGA = (\n", " + np.einsum(\"g, gk, gl, gu, gv -> uvkl\", frr, ao_0, ao_0, ao_0, ao_0)\n", " + 4 * np.einsum(\"g, rg, rgk, gl, gu, gv -> uvkl\", frg, rho_1, ao_1, ao_0, ao_0, ao_0)\n", " + 4 * np.einsum(\"g, gk, gl, rg, rgu, gv -> uvkl\", frg, ao_0, ao_0, rho_1, ao_1, ao_0)\n", " + 16 * np.einsum(\"g, wg, wgk, gl, rg, rgu, gv -> uvkl\", fgg, rho_1, ao_1, ao_0, rho_1, ao_1, ao_0)\n", " + 8 * np.einsum(\"g, rgk, gl, rgu, gv -> uvkl\", fg, ao_1, ao_0, ao_1, ao_0)\n", ")\n", "A_0_ao_GGA += A_0_ao_GGA.swapaxes(-1, -2)\n", "A_0_ao_GGA += A_0_ao_GGA.swapaxes(-3, -4)\n", "A_0_ao_GGA.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "随后我们就可以给出 GGA 下的 A 张量：\n", "\n", "$$\n", "A_{\\mu \\nu, \\kappa \\lambda} = 4 (\\mu \\nu | \\kappa \\lambda) - c_\\mathrm{x} (\\mu \\kappa | \\nu \\lambda) - c_\\mathrm{x} (\\mu \\lambda | \\nu \\kappa) + A_{\\mu \\nu, \\kappa \\lambda}^\\mathrm{GGA}\n", "$$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "A_0_ao = 4 * eri0_ao - cx * eri0_ao.swapaxes(-2, -3) - cx * eri0_ao.swapaxes(-1, -3) + A_0_ao_GGA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们可以拿现成的程序来验证上述 A 张量是否生成正确。pyxdh 中，它仍然可以用 `Ax0_Core` 与张量作乘积得到：(下式对除了 $p, q$ 之外的角标求和)\n", "\n", "$$\n", "A_{pq, mi} X_{mi} = C_{\\mu p} C_{\\nu q} A_{\\mu \\nu, \\kappa \\lambda} C_{\\kappa m} X_{mi} C_{\\lambda i}\n", "$$" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "X = np.random.randn(nmo, nocc)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(\n", " np.einsum(\"up, vq, uvkl, km, mi, li -> pq\", C, C, A_0_ao, C, X, Co),\n", " gradh.Ax0_Core(sa, sa, sa, so)(X)\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "**演示 A 张量 GGA 贡献部分的缩并效率**\n", "\n", "你可能会意识到上述生成 `A_0_ao_GGA` 的耗时相当长，但 `Ax0_Core` 的耗时却很短。请尝试用 `np.einsum_path` 分析时间复杂度，并指出一种策略，使得计算 $A_{pq, mi} X_{mi}$ 的耗时与 `Ax0_Core` 相对接近。\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": { "heading_collapsed": true }, "source": [ "### CP-HF 方程" ] }, { "cell_type": "markdown", "metadata": { "hidden": true, "scrolled": true }, "source": [ "既然我们已经知道 A 张量的计算方式，那么就可以通过 CP-HF 方程对 U 矩阵进行求解。首先我们仍然可以通过下式给出 B 矩阵：\n", "\n", "$$\n", "B_{pq}^\\mathbb{A} = F_{pq}^\\mathbb{A} - S_{pq}^\\mathbb{A} \\varepsilon_q - \\frac{1}{2} A_{pq, kl} S_{kl}^\\mathbb{A}\n", "$$" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "hidden": true }, "outputs": [], "source": [ "B_1 = (\n", " + F_1_mo\n", " - np.einsum(\"Atpq, q -> Atpq\", S_1_mo, e)\n", " - 0.5 * gradh.Ax0_Core(sa, sa, so, so)(S_1_mo[:, :, so, so])\n", ")" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "随后就可以通过 CP-HF 方程求解 `U_1_vo` $U_{ai}^{A_t}$ (在这份文档中，我们一般统一称 CP-KS 方程为 CP-HF 方程，因为解法与意义都几乎一致)：\n", "\n", "$$\n", "- (\\varepsilon_a - \\varepsilon_i) U_{ai}^\\mathbb{A} - A_{ai, bj} U_{bj}^\\mathbb{A} = B_{ai}^\\mathbb{A}\n", "$$" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "hidden": true }, "outputs": [], "source": [ "U_1_vo = cphf.solve(\n", " gradh.Ax0_Core(sv, so, sv, so, in_cphf=True),\n", " e,\n", " gradh.scf_eng.mo_occ,\n", " B_1[:, :, sv, so].reshape(natm * 3, nvir, nocc),\n", " tol=1e-6,\n", ")[0].reshape(natm, 3, nvir, nocc)" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "我们可以将其与 pyxdh 所给出的 U 矩阵作验证：" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "hidden": true, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(U_1_vo.ravel(), gradh.U_1_vo.ravel())" ] }, { "cell_type": "markdown", "metadata": { "hidden": true }, "source": [ "我们也可以用数值的 U 矩阵中非占-占据的分块来验证上述解析 U 矩阵：" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "hidden": true }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('