7. GGA 核坐标二阶 U 矩阵

我们曾经已经讨论过 RHF 下二阶 U 矩阵的计算了;但对于 GGA 而言,它还要引入一部分 GGA 的贡献。这种 GGA 的贡献并非是相当容易求取的,其最为困难之处在于求得 \(A_{pq, rs}^\mathbb{A}\)\(F_{\mu \nu}^{A_t B_s}\) 的过程。尽管该量仅仅包含了一阶导数,但确实只在二阶导数的计算过程中才会用上。

我们这一节的的重点问题就会是讨论 \(A_{pq, rs}^\mathbb{A}\)\(F_{\mu \nu}^{A_t B_s}\) 的计算,进而以 GGA 的二阶“未旋转”的 U 矩阵来作验证。我们在 MP2 的推导过程中知道,我们并不需要在二阶解析梯度中用到二阶 U 矩阵,因此,二阶 U 矩阵只是验证中间变量是否计算正确的手段。

程序变量名变更

出于程序简便的考量,这一节中特别地,所有 U_1 代表的不是 \(\mathscr{U}_{pq}^\mathbb{A}\),而是原先用 U_1_nr 所指代的未经轨道旋转的 \(U_{pq}^\mathbb{A}\)

类似地,gradh_nr 会被 gradhhessh_nr 会被 hessh 替代。

7.1. 准备工作

[31]:
%matplotlib notebook

from pyscf import gto, scf, dft, lib, grad, hessian
import numpy as np
from functools import partial
import warnings
from matplotlib import pyplot as plt
from pyxdh.Utilities import NucCoordDerivGenerator, DipoleDerivGenerator, NumericDiff, GridHelper, KernelHelper, GridIterator
from pyxdh.DerivOnce import GradSCF
from pyxdh.DerivTwice import HessSCF

np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
np.allclose = partial(np.allclose, atol=1e-6, rtol=1e-4)
np.set_printoptions(5, linewidth=150, suppress=True)
warnings.filterwarnings("ignore")
[2]:
mol = gto.Mole()
mol.atom = """
O  0.0  0.0  0.0
O  0.0  0.0  1.5
H  1.0  0.0  0.0
H  0.0  0.7  1.0
"""
mol.basis = "6-31G"
mol.verbose = 0
mol.build()
[2]:
<pyscf.gto.mole.Mole at 0x7f8f545a6340>
[3]:
def mol_to_grids(mol, atom_grid=(75, 302)):
    grids = dft.Grids(mol)
    grids.atom_grid = atom_grid
    grids.becke_scheme = dft.gen_grid.stratmann
    grids.prune = None
    grids.build()
    return grids
grids = mol_to_grids(mol)
[4]:
def mol_to_scf(mol):
    scf_eng = dft.RKS(mol)
    scf_eng.grids = mol_to_grids(mol)
    scf_eng.xc = "B3LYPg"
    scf_eng.conv_tol = 1e-10
    return scf_eng.run()
[5]:
gradh = GradSCF({"scf_eng": mol_to_scf(mol), "cphf_tol": 1e-12, "rotation": False})
hessh = HessSCF({"deriv_A": gradh, "deriv_B": gradh, "rotation": False})
[6]:
nmo, nao, natm, nocc, nvir = gradh.nao, gradh.nao, gradh.natm, gradh.nocc, gradh.nvir
so, sv, sa = gradh.so, gradh.sv, gradh.sa
mol_slice = gradh.mol_slice
C, Co, Cv, e, eo, ev, D = gradh.C, gradh.Co, gradh.Cv, gradh.e, gradh.eo, gradh.ev, gradh.D
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
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
H_1_ao, S_1_ao, eri1_ao, F_1_ao = gradh.H_1_ao, gradh.S_1_ao, gradh.eri1_ao, gradh.F_1_ao
H_1_mo, S_1_mo, eri1_mo, F_1_mo = gradh.H_1_mo, gradh.S_1_mo, gradh.eri1_mo, gradh.F_1_mo
Ax0_Core, B_1, U_1, U_1_vo = gradh.Ax0_Core, gradh.B_1, gradh.U_1, gradh.U_1_vo
H_2_ao, S_2_ao, eri2_ao = hessh.H_2_ao, hessh.S_2_ao, hessh.eri2_ao
H_2_mo, S_2_mo, eri2_mo = hessh.H_2_mo, hessh.S_2_mo, hessh.eri2_mo
[38]:
grdh = GridHelper(mol, grids, D)
kerh = KernelHelper(grdh, "B3LYPg", deriv=3)
cx, xc = gradh.cx, kerh.xc
[20]:
ngrid = grdh.ngrid
ao_0, ao_1, ao_2, ao_3 = grdh.ao_0, grdh.ao_1, grdh.ao_2, grdh.ao_3
rho_1, rho_2 = grdh.rho_1, grdh.rho_2
A_rho_1 = grdh.A_rho_1.reshape(natm * 3, ngrid)
A_gamma_1 = grdh.A_gamma_1.reshape(natm * 3, ngrid)
fr, fg, frr, frg, fgg = kerh.fr, kerh.fg, kerh.frr, kerh.frg, kerh.fgg
frrr, frrg, frgg, fggg = kerh.frrr, kerh.frrg, kerh.frgg, kerh.fggg
[9]:
def grad_generator(mol):
    scf_eng = mol_to_scf(mol)
    config = {"scf_eng": scf_eng, "cphf_tol": 1e-12, "rotation": False}
    return GradSCF(config)
gradn = NucCoordDerivGenerator(mol, grad_generator)

与 RHF 相同地,我们也可以求取数值二阶 U 矩阵 n_U_2 \(U_{pq}^{\mathbb{AB}}\)

\[U_{pq}^\mathbb{AB} = \frac{\partial U_{pq}^\mathbb{A}}{\partial \mathbb{B}} + U_{pm}^\mathbb{B} U_{mq}^\mathbb{A}\]
[10]:
nd_U_1 = NumericDiff(gradn, lambda gradh_nr: gradh_nr.U_1).derivative.swapaxes(0, 1)
n_U_2 = nd_U_1 + np.einsum("Bpm, Amq -> ABpq", U_1, U_1)
n_U_2.shape
[10]:
(12, 12, 22, 22)

我们下面补充定义一个函数 two_dim_to_one,其作用是使类似于 \((A, t, \mu, \nu)\) 的维度变为 \((A_t, \mu, \nu)\),即将前两个维度合并。

[28]:
def two_dim_to_one(mat):
    dim = mat.shape
    return mat.reshape([dim[0] * dim[1]] + list(dim[2:]))

7.2. 一阶 A 张量计算

我们知道,二阶 U 矩阵的计算过程中,摆在我们面前的问题会有 Fock 矩阵的二阶 Skeleton 导数 \(F_{\mu \nu}^\mathbb{AB}\),以及一阶 A 张量 \(A_{\mu \nu, \kappa \lambda}^\mathbb{A}\) 或等价地 \(A_{pq, rs}^\mathbb{A}\) 的计算过程。看起来,\(F_{\mu \nu}^\mathbb{AB}\) 作为 Fock 矩阵的二阶 Skeleton 导数相对比较容易推导与求取;当然事实也确实如此。但从定义过程上,应当是一阶 A 张量的计算在前。我们回顾一下推导过程。

\[\frac{\partial F_{\mu \nu}^\mathbb{A}}{\partial \mathbb{B}} = F_{\mu \nu}^\mathbb{AB} + \frac{\partial A_{\mu \nu, \kappa \lambda}}{\partial \mathbb{A}} C_{\kappa m} U_{mk}^\mathbb{B} C_{\lambda k}\]

其中,等式左边与右边的第二项是可以确切定义的项,因此它们可以用以定义 \(F_{\mu \nu}^\mathbb{AB}\)。而一阶 A 张量的定义是根据等式右边的第二项而来:

\[A_{pq, rs}^\mathbb{A} = C_{\mu p} C_{\nu q} \frac{\partial A_{\mu \nu, \kappa \lambda}}{\partial \mathbb{A}} C_{\kappa r} C_{\lambda s} = C_{\mu p} C_{\nu q} A_{\mu \nu, \kappa \lambda}^\mathbb{A} C_{\kappa r} C_{\lambda s}\]

由于在实际运算过程中,使用以分子轨道的 \(A_{pq, rs}^\mathbb{A}\) 进行缩并的效率,比原始的 \(A_{\mu \nu, \kappa \lambda}^\mathbb{A}\) 要好不少,因此我们下面会实际介绍的是缩并函数 \(A_{pq, rs}^\mathbb{A} X_{rs}^\mathbb{B}\) Ax1_Core 的编写。

7.2.1. HF 贡献部分

我们先以 HF 贡献部分的程序编写,来熟悉 Ax1_Core 函数的编写原则。

Ax1_Core 的输入会是 \(p, q, r, s\) 所需要的轨道分割;输出则是另一个函数,暂且记为 fxfx 的输入则是 \(X_{rs}^\mathbb{B}\),输出的是维度以 \((\mathbb{A}, \mathbb{B}, p, q)\) 储存的,对 \(r, s\) 角标求和的 \(A_{pq, rs}^\mathbb{A} X_{rs}^\mathbb{B}\) 的张量。

对于 \(A_{pq, rs}^\mathbb{A} X_{rs}^\mathbb{B}\) 而言,我们先计算 Ax1_Core_HF \(A_{pq, rs}^{\mathrm{HF}, \mathbb{A}} X_{rs}^\mathbb{B}\),并且我们假定 \(\mathbb{B}\) 是 5 维度的,\((p, q, r, s)\) 分别代表非占-占据-非占-占据的分割 (在这里 \(p, q, r, s\) 暂且当作维度不确定的分子轨道):

[11]:
X = np.random.randn(5, nvir, nocc)
[12]:
def Ax1_Core_HF(sp, sq, sr, ss):                                                # Line  1
    size_A = U_1.shape[:-2]                                                     # Line  2
    def fx(X):                                                                  # Line  3
        dmX = C[:, sr] @ X @ C[:, ss].T                                         # Line  4
        dmX += dmX.swapaxes(-1, -2)                                             # Line  5
        size_B = dmX.shape[:-2]                                                 # Line  6
        dmX.shape = (np.prod(size_B, dtype=int), dmX.shape[-2], dmX.shape[-1])  # Line  7
        ax_ao = np.zeros([np.prod(size_A, dtype=int)] + list(dmX.shape))        # Line  8
        ax_ao += np.einsum("Auvkl, Bkl -> ABuv", eri1_ao, dmX)                  # Line  9
        ax_ao -= 0.5 * cx * np.einsum("Aukvl, Bkl -> ABuv", eri1_ao, dmX)       # Line 10
        ax_ao += ax_ao.swapaxes(-1, -2)                                         # Line 11
        Ax = C[:, sp].T @ ax_ao @ C[:, sq]                                      # Line 12
        Ax.shape = list(size_A) + list(size_B) + list(Ax.shape[-2:])            # Line 13
        return Ax                                                               # Line 14
    return fx                                                                   # Line 15

我们可以用下述代码验证结果:

\[A_{pq, rs}^{\mathrm{HF}, \mathbb{A}} X_{rs}^\mathbb{B} = \big( 4 (pq | rs)^\mathbb{AB} - (pr | qs)^\mathbb{AB} - (ps | qr)^\mathbb{AB} \big) X_{rs}^\mathbb{B}\]
[13]:
np.allclose(
    Ax1_Core_HF(sv, so, sv, so)(X),
    + 4 * np.einsum("Aaibj, Bbj -> ABai", eri1_mo[:, sv, so, sv, so], X)
    - cx * np.einsum("Aabij, Bbj -> ABai", eri1_mo[:, sv, sv, so, so], X)
    - cx * np.einsum("Aajib, Bbj -> ABai", eri1_mo[:, sv, so, so, sv], X)
)
[13]:
True

现在,我们对上述代码作细致的说明:

  • Line 2:记录 \(\mathbb{A}\) 的原始维度。由于 \(\mathbb{A}\) 在核坐标中,可能代表维度 \((A_t, )\),也可能代表 \((A, t)\)size_A 则记录下该维度。

  • Line 4-5:求出对 \(r, s\) 角标求和的 dmX \(C_{\kappa r} X_{rs}^\mathbb{B} C_{\lambda s} + C_{\lambda r} X_{rs}^\mathbb{B} C_{\kappa s}\),但保留 \(\mathbb{B}\) 的原始维度。

  • Line 6:记录 \(\mathbb{B}\) 的原始维度。同 Line 2。

  • Line 7:压平 \(\mathbb{B}\) 所指代的维度,并让 dmX 成为三维张量,维度为 \((\mathbb{B}, \kappa, \lambda)\)

  • Line 8:初始化 ax_ao,该张量在最后储存的会是对 \(\kappa, \lambda, r, s\) 角标求和的 \(A_{\mu \nu, \kappa \lambda}^{\mathrm{HF}, \mathbb{A}} C_{\kappa r} C_{\lambda s} X_{rs}^\mathbb{B}\),维度为四维度的 \((\mathbb{A}, \mathbb{B}, r, s)\)

  • Line 9-11:具体计算了 \(A_{\mu \nu, \kappa \lambda}^{\mathrm{HF}, \mathbb{A}} C_{\kappa r} X_{rs}^\mathbb{B} C_{\lambda s}\)

  • Line 12:对分子轨道张量进行进一步缩并,得到 \(A_{pq, rs}^{\mathrm{HF}, \mathbb{A}} X_{rs}^\mathbb{B}\)

  • Line 13:最终对 A 张量的维度作一定的改变回到应有的情况。

我们注意到在整个流程中,我们使用到了

\[X_{\kappa \lambda}^\mathbb{B} = C_{\kappa b} X_{bj}^\mathbb{B} C_{\lambda j} + \mathrm{swap} (\kappa, \lambda)\]

以及最后的计算中,用到 \(\mathrm{swap} (\mu, \nu)\)。因此即使是相当简单的 HF 贡献的一阶 A 张量缩并 (从公式表达上来看仅有三项),但为了效率上的考量,我们会用代码的复杂化作牺牲。在 GGA 的一阶 A 张量计算过程中,我们的整体过程也大致如此。

7.2.2. GGA 贡献部分:代码

这里的内容推导需要先参考 pyxdh A 张量的张量缩并方式 或 XYG3 CheatSheet。由于该张量计算方式的特殊性,我们几乎需要一次性地用代码生成;并且不能方便地对代码的正确性作直接的验证。因此,我们先对代码作说明,随后依着代码的流程进行公式推导。

[62]:
def Ax1_Core(sp, sq, sr, ss):
    size_A = U_1.shape[:-2]
    # Block 1: Generate `dmU`
    dmU = C @ U_1[:, :, so] @ Co.T
    dmU += dmU.swapaxes(-1, -2)
    def fx(X):
        dmX = C[:, sr] @ X @ C[:, ss].T
        dmX += dmX.swapaxes(-1, -2)
        size_B = dmX.shape[:-2]
        dmX.shape = (np.prod(size_B, dtype=int), dmX.shape[-2], dmX.shape[-1])
        ax_ao = np.zeros([np.prod(size_A, dtype=int)] + list(dmX.shape))
        # Block 2: Generate `grdit`
        grdit = GridIterator(mol, grids, D, deriv=3)
        for grdh in grdit:
            kerh = KernelHelper(grdh, xc, deriv=3)
            # Block 3: Define some kernel and density skeleton derivative
            pd_frr = kerh.frrr * two_dim_to_one(grdh.A_rho_1) + kerh.frrg * two_dim_to_one(grdh.A_gamma_1)
            pd_frg = kerh.frrg * two_dim_to_one(grdh.A_rho_1) + kerh.frgg * two_dim_to_one(grdh.A_gamma_1)
            pd_fgg = kerh.frgg * two_dim_to_one(grdh.A_rho_1) + kerh.fggg * two_dim_to_one(grdh.A_gamma_1)
            pd_fg = kerh.frg * two_dim_to_one(grdh.A_rho_1) + kerh.fgg * two_dim_to_one(grdh.A_gamma_1)
            pd_rho_1 = two_dim_to_one(grdh.A_rho_2)
            # Block 4: Form `dmX` density grid
            rho_X_0 = np.array([grdh.get_rho_0(dm) for dm in dmX])
            rho_X_1 = np.array([grdh.get_rho_1(dm) for dm in dmX])
            pd_rho_X_0 = np.array([two_dim_to_one(grdh.get_A_rho_1(dm)) for dm in dmX]).swapaxes(0, 1)
            pd_rho_X_1 = np.array([two_dim_to_one(grdh.get_A_rho_2(dm)) for dm in dmX]).swapaxes(0, 1)
            # Block 5: Define temporary M intermediates (Original and Skeleton derivative)
            M_0 = (
                    + np.einsum("g, Bg -> Bg", kerh.frr, rho_X_0)
                    + 2 * np.einsum("g, wg, Bwg -> Bg", kerh.frg, grdh.rho_1, rho_X_1)
            )
            M_1 = (
                    + 4 * np.einsum("g, Bg, rg -> Brg", kerh.frg, rho_X_0, grdh.rho_1)
                    + 8 * np.einsum("g, wg, Bwg, rg -> Brg", kerh.fgg, grdh.rho_1, rho_X_1, grdh.rho_1)
                    + 4 * np.einsum("g, Brg -> Brg", kerh.fg, rho_X_1)
            )
            pd_M_0 = (
                    + np.einsum("Ag, Bg -> ABg", pd_frr, rho_X_0)
                    + np.einsum("g, ABg -> ABg", kerh.frr, pd_rho_X_0)
                    + 2 * np.einsum("Ag, wg, Bwg -> ABg", pd_frg, grdh.rho_1, rho_X_1)
                    + 2 * np.einsum("g, Awg, Bwg -> ABg", kerh.frg, pd_rho_1, rho_X_1)
                    + 2 * np.einsum("g, wg, ABwg -> ABg", kerh.frg, grdh.rho_1, pd_rho_X_1)
            )
            pd_M_1 = (
                    + 4 * np.einsum("Ag, Bg, rg -> ABrg", pd_frg, rho_X_0, grdh.rho_1)
                    + 4 * np.einsum("g, Bg, Arg -> ABrg", kerh.frg, rho_X_0, pd_rho_1)
                    + 4 * np.einsum("g, ABg, rg -> ABrg", kerh.frg, pd_rho_X_0, grdh.rho_1)
                    + 8 * np.einsum("Ag, wg, Bwg, rg -> ABrg", pd_fgg, grdh.rho_1, rho_X_1, grdh.rho_1)
                    + 8 * np.einsum("g, Awg, Bwg, rg -> ABrg", kerh.fgg, pd_rho_1, rho_X_1, grdh.rho_1)
                    + 8 * np.einsum("g, wg, Bwg, Arg -> ABrg", kerh.fgg, grdh.rho_1, rho_X_1, pd_rho_1)
                    + 8 * np.einsum("g, wg, ABwg, rg -> ABrg", kerh.fgg, grdh.rho_1, pd_rho_X_1, grdh.rho_1)
                    + 4 * np.einsum("Ag, Brg -> ABrg", pd_fg, rho_X_1)
                    + 4 * np.einsum("g, ABrg -> ABrg", kerh.fg, pd_rho_X_1)
            )
            # Contribution 1: pdSkeleton_M * ao_grid
            contrib1 = np.zeros((natm * 3, dmX.shape[0], nao, nao))
            contrib1 += np.einsum("ABg, gu, gv -> ABuv", pd_M_0, grdh.ao_0, grdh.ao_0)
            contrib1 += np.einsum("ABrg, rgu, gv -> ABuv", pd_M_1, grdh.ao_1, grdh.ao_0)
            # Contribution 2: M * pdSkeleton_ao_grid
            tmp_contrib = (
                    - 2 * np.einsum("Bg, tgu, gv -> tBuv", M_0, grdh.ao_1, grdh.ao_0)
                    - np.einsum("Brg, trgu, gv -> tBuv", M_1, grdh.ao_2, grdh.ao_0)
                    - np.einsum("Brg, tgu, rgv -> tBuv", M_1, grdh.ao_1, grdh.ao_1)
            )
            contrib2 = np.zeros((natm, 3, dmX.shape[0], nao, nao))
            for A in range(natm):
                sA = mol_slice(A)
                contrib2[A, :, :, sA] += tmp_contrib[:, :, sA]
            contrib2.shape = (natm * 3, dmX.shape[0], nao, nao)
            # Block 6: U contribution to pdU_M
            rho_U_0 = np.einsum("Auv, gu, gv -> Ag", dmU, grdh.ao_0, grdh.ao_0)
            rho_U_1 = 2 * np.einsum("Auv, rgu, gv -> Arg", dmU, grdh.ao_1, grdh.ao_0)
            gamma_U_0 = 2 * np.einsum("rg, Arg -> Ag", grdh.rho_1, rho_U_1)
            pdU_frr = kerh.frrr * rho_U_0 + kerh.frrg * gamma_U_0
            pdU_frg = kerh.frrg * rho_U_0 + kerh.frgg * gamma_U_0
            pdU_fgg = kerh.frgg * rho_U_0 + kerh.fggg * gamma_U_0
            pdU_fg = kerh.frg * rho_U_0 + kerh.fgg * gamma_U_0
            pdU_rho_1 = rho_U_1
            # Block 7: Define temporary M intermediates (U derivative)
            pdU_M_0 = (
                    + np.einsum("Ag, Bg -> ABg", pdU_frr, rho_X_0)
                    + 2 * np.einsum("Ag, wg, Bwg -> ABg", pdU_frg, grdh.rho_1, rho_X_1)
                    + 2 * np.einsum("g, Awg, Bwg -> ABg", kerh.frg, pdU_rho_1, rho_X_1)
            )
            pdU_M_1 = (
                    + 4 * np.einsum("Ag, Bg, rg -> ABrg", pdU_frg, rho_X_0, grdh.rho_1)
                    + 4 * np.einsum("g, Bg, Arg -> ABrg", kerh.frg, rho_X_0, pdU_rho_1)
                    + 8 * np.einsum("Ag, wg, Bwg, rg -> ABrg", pdU_fgg, grdh.rho_1, rho_X_1, grdh.rho_1)
                    + 8 * np.einsum("g, Awg, Bwg, rg -> ABrg", kerh.fgg, pdU_rho_1, rho_X_1, grdh.rho_1)
                    + 8 * np.einsum("g, wg, Bwg, Arg -> ABrg", kerh.fgg, grdh.rho_1, rho_X_1, pdU_rho_1)
                    + 4 * np.einsum("Ag, Brg -> ABrg", pdU_fg, rho_X_1)
            )
            # Contribution 3: pdU_M * ao_grid
            contrib3 = np.zeros((natm * 3, dmX.shape[0], nao, nao))
            contrib3 += np.einsum("ABg, gu, gv -> ABuv", pdU_M_0, grdh.ao_0, grdh.ao_0)
            contrib3 += np.einsum("ABrg, rgu, gv -> ABuv", pdU_M_1, grdh.ao_1, grdh.ao_0)
            # GGA Contribution Summation
            ax_ao += contrib1 + contrib2 + contrib3
        # HF contribution
        ax_ao += np.einsum("Auvkl, Bkl -> ABuv", eri1_ao, dmX)
        ax_ao -= 0.5 * cx * np.einsum("Aukvl, Bkl -> ABuv", eri1_ao, dmX)
        # Swap mu, nu and final contraction
        ax_ao += ax_ao.swapaxes(-1, -2)
        Ax = C[:, sp].T @ ax_ao @ C[:, sq]
        Ax.shape = list(size_A) + list(size_B) + list(Ax.shape[-2:])
        return Ax
    return fx
[65]:
np.allclose(
    Ax1_Core(sv, so, sv, so)(X),
    gradh.Ax1_Core(sv, so, sv, so)(X)  # pyxdh approach
)
[65]:
True

我们下面对其中一部分的代码块作说明。由于记号表达的不顺畅,很多公式项难以确切地表示出来,可能需要读者自行理解了。

至于要如何验证,上述生成的或者 pyxdh 所生成的 Ax1_Core 是正确的,这还需要我们用它来生成二阶 B 矩阵,从而得到二阶 U 矩阵来验证。

7.2.3. Block 1:生成 dmU \(U_{\mu \nu}^\mathbb{A}\)

由于我们会经常地使用到以 \(C_{\mu m} U_{mi}^\mathbb{A} C_{\nu i} + \mathrm{swap} (\mu, \nu)\) 作为密度矩阵所生成的格点。我们定义 dmU \(U_{\mu \nu}^\mathbb{A}\)

\[U_{\mu \nu}^\mathbb{A} = C_{\mu m} U_{mi}^\mathbb{A} C_{\nu i} + \mathrm{swap} (\mu, \nu)\]

我们以后会用 \(\rho[U_{\mu \nu}^\mathbb{A}]\) 来表示使用了 \(U_{\mu \nu}^\mathbb{A}\) 的密度格点。相应地,我们也会用 \(\rho[X_{\mu \nu}^\mathbb{A}]\) 表示使用了 \(X_{\mu \nu}^\mathbb{A}\) 的密度格点。但若使用的是 \(D_{\mu \nu}\),我们可能会不明确写出其所使用的原子轨道密度。

7.2.4. Block 2:生成实例 grdit

pyxdh 中,实际用于格点积分的量是 GridIterator 的实例,而非以前文档中出现的 GridHelper 实例;但这两者是相当相似的。GridIterator 的实例 grdit 是一个迭代器,用于分批产生 DFT 积分格点。之所以这么做,是因为当格点数量太大时,就可以在较少的内存下分批次地作格点积分。它有些像 PySCF 中的 NumInt.block_loop,一定程度上就是它的外封装。

我们注意到下述语句:

for grdh in grdit:

我们说,每次 GridIterator 实例所迭代的内容可以看作是 GridHelper 的实例,同样地都能生成原子轨道格点、密度格点、以及各种梯度量,拥有极其相似的 API。但一方面,GridIterator 事实上与 GridHelper 毫无关系,只是从程序的调用方式上非常类似;二来,GridIterator 类还允许我们用函数,譬如 get_rho_0,从密度矩阵直接生成密度格点。

7.2.5. Block 3:生成密度与格点 Skeleton 导数偏导量

尽管我们以前不断声明,对解析梯度以 Skeleton 导数与 U 导数为分类进行拆分是没有明确的物理意义或依据的——特别是对于与 Fock 矩阵相关的量,由于 \(F_{\mu \nu}^\mathbb{A}\) 中同时存在着原子轨道与密度矩阵,因此作者认为,如果抛开约定俗成,似乎多少不太适合说 \(F_{\mu \nu}^\mathbb{A}\) 是 Fock 矩阵的 Skeleton 导数,因为它并不等于 \(\partial_\mathbb{A} F_{\mu \nu}\)

但将导数分离的做法而言,分成 Skeleton 导数与 U 导数还是对公式的整理有相当的好处。我们先回顾到,在零阶 A 张量计算过程中,我们使用到了

\[A_{\mu \nu, \kappa \lambda} X_{\kappa \lambda}^\mathbb{A} = (\mu \nu | \kappa \lambda) X_{\kappa \lambda}^\mathbb{A} - \frac{c_\mathrm{x}}{2} (\mu \kappa | \nu \lambda) X_{\kappa \lambda}^\mathbb{A} + M^\mathbb{A} \phi_\mu \phi_\nu + M_r^\mathbb{A} \phi_{r \mu} \phi_\nu + \mathrm{swap} (\mu, \nu)\]

其中,上式的前两项是 HF 贡献的部分,后面的两项则是 GGA 的贡献。那么,GGA 的贡献大致地拆分为三类贡献:一类是 \(M^\mathbb{A}\) 或类似项的 Skeleton 导数,一类是 \(\phi_\mu\) 或类似项的 Skeleton 导数,最后则是 \(M^\mathbb{A}\) 或类似项的 U 导数。

我们回顾到 M_0 \(M^\mathbb{B}\) 表达为

\[M^\mathbb{B} = f_{\rho \rho} \rho[X_{\kappa \lambda}^\mathbb{B}] + 2 f_{\rho \gamma} \rho_w \rho_w[X_{\kappa \lambda}^\mathbb{B}]\]

M_1 \(M_r^\mathbb{B}\) 表达为

\[M_r^\mathbb{B} = 4 f_{\rho \gamma} \rho_r \rho[X_{\kappa \lambda}^\mathbb{B}] + 8 f_{\rho \gamma} \rho_w \rho_r \rho_w[X_{\kappa \lambda}^\mathbb{B}] + 4 f_\gamma \rho_r[X_{\kappa \lambda}^\mathbb{B}]\]

那么,为了求取第一类导数,我们需要给出上式计算中出现项所对应的 Skeleton 导数。我们先求与 \(\rho[X_{\kappa \lambda}^\mathbb{B}]\) 无关的所有 Skeleton 导数。

  • pd_frr \(\partial_\mathbb{A} f_{\rho \rho} \xleftarrow{\text{Skeleton Derivative}} f_{\rho \rho \rho} \rho^\mathbb{A} + f_{\rho \rho \gamma} \gamma^\mathbb{A}\)

  • pd_frg \(\partial_\mathbb{A} f_{\rho \gamma} \xleftarrow{\text{Skeleton Derivative}} f_{\rho \rho \gamma} \rho^\mathbb{A} + f_{\rho \gamma \gamma} \gamma^\mathbb{A}\)

  • pd_fgg \(\partial_\mathbb{A} f_{\gamma \gamma} \xleftarrow{\text{Skeleton Derivative}} f_{\rho \gamma \gamma} \rho^\mathbb{A} + f_{\gamma \gamma \gamma} \gamma^\mathbb{A}\)

  • pd_fg \(\partial_\mathbb{A} f_{\gamma} \xleftarrow{\text{Skeleton Derivative}} f_{\rho \gamma} \rho^\mathbb{A} + f_{\gamma \gamma} \gamma^\mathbb{A}\)

  • pd_rho_1 \(\partial_\mathbb{A} \rho_r \xleftarrow{\text{Skeleton Derivative}} \rho_r^\mathbb{A}\)

7.2.6. Block 4:生成与 \(\rho[X_{\kappa \lambda}^\mathbb{B}]\) 有关的 Skeleton 导数

  • rho_X_0 \(\rho[X_{\kappa \lambda}^\mathbb{B}] = X_{\kappa \lambda}^\mathbb{B} \phi_\kappa \phi_\lambda\)

  • rho_X_1 \(\rho_r[X_{\kappa \lambda}^\mathbb{B}] = 2 X_{\kappa \lambda}^\mathbb{B} \phi_{r \kappa} \phi_\lambda\)

  • pd_rho_X_0 \(\rho^{A_t}[X_{\kappa \lambda}^\mathbb{B}] = - 2 X_{\kappa \lambda}^\mathbb{B} \phi_{t \kappa_A} \phi_\lambda\)

  • pd_rho_X_1 \(\rho_r^\mathbb{A}[X_{\kappa \lambda}^\mathbb{B}] = - 2 X_{\kappa \lambda}^\mathbb{B} \phi_{tr \kappa_A} \phi_\lambda - 2 X_{\kappa \lambda}^\mathbb{B} \phi_{r \kappa} \phi_{t \lambda_A}\)

由于 GridIterator 类提供了通过密度直接生成格点的函数,因此我们不需要额外地编写生成代码了。

7.2.7. Block 5:M 矩阵及其一阶 Skeleton 导数

  • M_0 \(M^\mathbb{B} = f_{\rho \rho} \rho[X_{\kappa \lambda}^\mathbb{B}] + f_{\rho \gamma} \rho_w \rho_w[X_{\kappa \lambda}^\mathbb{B}]\)

  • M_1 \(M_r^\mathbb{B} = 4 f_{\rho \gamma} \rho_r \rho[X_{\kappa \lambda}^\mathbb{B}] + 8 f_{\rho \gamma} \rho_w \rho_r \rho_w[X_{\kappa \lambda}^\mathbb{B}] + 4 f_\gamma \rho_r[X_{\kappa \lambda}^\mathbb{B}]\)

  • pd_M_0\(M^\mathbb{B}\) 的每一项依链式法则,分别求取 Skeleton 导数

  • pd_M_1\(M_r^\mathbb{B}\) 的每一项依链式法则,分别求取 Skeleton 导数

7.2.8. Block 6:U 矩阵作为密度矩阵所生成的格点

上面仅仅讨论了 Skeleton 导数的问题,但所有的 U 矩阵导数也是我们所需要求取的。我们注意到,在 GGA 对 A 张量缩并的贡献项中,\(\phi_{\mu}\)\(\phi_{r \mu}\) 是不会产生 U 导数贡献的。因此,我们只需要考虑与 \(\partial_\mathbb{A} M^\mathbb{B}\)\(\partial_\mathbb{A} M_r^\mathbb{B}\) 的项。

但我们还需要注意到,不是所有的项都需要被求到 U 导数。譬如,对于 \(\rho[X_{\kappa \lambda}^\mathbb{B}]\) 而言,我们会注意到

\[\rho[X_{\kappa \lambda}^\mathbb{B}] = X_{rs}^\mathbb{B} \phi_{\kappa} \phi_\lambda (C_{\kappa r} C_{\lambda s} + C_{\kappa s} C_{\lambda r})\]

上式中,\(X_{rs}^\mathbb{B}\) 是不适合求导数的;因此能得到 U 导数的项只有两个轨道系数矩阵 \(C_{\kappa r}\) 等项。但我们留意到一阶 A 张量的定义:

\[A_{pq, rs}^\mathbb{A} = C_{\mu p} C_{\nu q} \frac{\partial A_{\mu \nu, \kappa \lambda}}{\partial \mathbb{A}} C_{\kappa r} C_{\lambda s}\]

因此,事实上类似于 \(C_{\kappa r}\) 等项并不应当被求导。因此,我们之后在求取 \(\partial_\mathbb{A} M^\mathbb{B}\)\(\partial_\mathbb{A} M_r^\mathbb{B}\) 的 U 导数贡献时,就无需要考虑与 \(X_{\kappa \lambda}^\mathbb{B}\) 有关的密度格点了。

文档未完成

这份文档可能在短期之内不会再更新了。这是应由于作者的编写热情降低、以及公式符号的复杂化引起的。

作者希望读者在读到这里之后,能自行地推导 Fock 矩阵的二阶 Skeleton 导数,并且以此生成 GGA 的二阶解析 U 矩阵,并与数值 U 矩阵作核验。

这之后恐怕就没有更困难的计算了。对于 XYG3 型泛函的 Hessian 而言,只需要在现在的基础上,在 MP2 生成的部分注意 \(L_{ai}^\mathrm{PT2+}\) 的项,补上非自洽的 \(F_{ai}^\mathrm{n}\) 即可。对于其它导数,需要留意被求导量 \(\mathbb{A}, \mathbb{B}\) 的顺序是否正确。

如果要拿代码进行核验,不妨 Hack 一下 pyxdh 的代码。你或许会发现 HessNCDFTHessMP2HessXDH 的代码意外地很简单,因此真正需要加深理解的,只有 RHF, MP2 与 GGA 的 Hessian 求取,与二阶 U 矩阵的计算。有这些基础后,其它的解析导数譬如 B2PLYP 型、XYG3 型求取就会轻松许多。