自己的本研课题已经大致确立在KS-Inversion这个主题上。
本篇文章是2019年nc上关于KS-Inversion的一篇论文 [1] 的笔记，这应该代表这个领域当前的最前沿进展了。
由于这篇文章的一个核心优化点在基组函数上，在欣赏他们的工作之前，先补充一些自己之前遗漏的知识点。
Why Gaussian Basis Lose
基组是如何工作的
先前在介绍DFT时，我们稍微提了一下在计算机中，我们求解薛定谔的方式。例如KS中单电子的轨道方程为
\[\hat{H} \phi_i = \epsilon_i \phi_i\]用基组展开分子轨道：$\phi_i = \sum_\nu C_{\nu i} \chi_\nu$，代入上式得到
\[\hat{H} \sum_\nu C_{\nu i} \chi_\nu = \epsilon_i \sum_\nu C_{\nu i} \chi_\nu\]为了把算符方程变成计算机认识的矩阵方程，我们在等式左边同时左乘另一个基函数 $\chi_\mu$，并在全空间做积分
\[\sum_\nu C_{\nu i} \int \chi_\mu \hat{H} \chi_\nu d\mathbf{r} = \epsilon_i \sum_\nu C_{\nu i} \int \chi_\mu \chi_\nu d\mathbf{r}\]定义左边的积分为哈密顿矩阵元 $H_{\mu\nu}$，右边的积分为基组的重叠矩阵元 $S_{\mu\nu}$：\(\sum_\nu H_{\mu\nu} C_{\nu i} = \epsilon_i \sum_\nu S_{\mu\nu} C_{\nu i}\)
因此薛定谔方程在基组下可以写成矩阵形式
\[HC = ESC\]轨道系数矩阵 $C$ 的第 $i$ 列向量 $\mathbf{c}i = [C{1i}, C_{2i}, \dots, C_{Mi}]^T$ 完整决定了第 $i$ 个分子轨道的空间形态与数学轮廓。
本征能量矩阵是一个对角矩阵，对角线上的第 $i$ 个元素 $E_{ii}$（通常记为 $\varepsilon_i$）代表了第 $i$ 个分子轨道的本征能量。
\[E = \begin{pmatrix} \varepsilon_1 & 0 & \dots & 0 \\ 0 & \varepsilon_2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \varepsilon_M \end{pmatrix}\]可能你会问能量不tm是未知的吗，好问题，问就是自洽求解（x）。这个留给读者思考，等我想清楚了再补上。
在有限基组下，求出来的系数矩阵$C$本质上是一个近似解，或者说真实解在有限空间的投影。这也就是所谓基组截断误差。
高斯基组（GTO）
高斯基组函数在直角坐标下的形式为
\[\chi_{lmn}(\mathbf{r} - \mathbf{R}_A) = N \cdot (x - x_A)^l (y - y_A)^m (z - z_A)^n \cdot e^{-\alpha |\mathbf{r} - \mathbf{R}_A|^2}\]其中$N$是归一化常数，如果用于拟合波函数，那么前置多项式 $x^l y^m z^n$ 的内涵就是角动量与轨道形状，衰减项 $e^{-\alpha r^2}$则体现波函数的径向衰减。基组数学原点会选择与原子核的物理位置完全重合。
这个基组可以一眼看出不正交，之所以被广泛使用完全是因为它算得快，具有独特的高斯乘积定理，在计算多中心两电子库仑排斥积分时快得惊人，能把四中心积分直接变成双中心积分。
\[\langle AB | \frac{1}{r_{12}} | CD \rangle = \int_{\Omega_1} \int_{\Omega_2} \frac{\chi_A^*(\mathbf{r}_1) \chi_B(\mathbf{r}_1) \cdot \chi_C^*(\mathbf{r}_2) \chi_D(\mathbf{r}_2)}{|\mathbf{r}_1 - \mathbf{r}_2|} d\mathbf{r}_1 d\mathbf{r}_2\] \[\langle AB | \frac{1}{r_{12}} | CD \rangle = K_{AB} K_{CD} \int_{\Omega_1} \int_{\Omega_2} \frac{e^{-\gamma_1 (\mathbf{r}_1 - \mathbf{R}_P)^2} \cdot e^{-\gamma_2 (\mathbf{r}_2 - \mathbf{R}_Q)^2}}{|\mathbf{r}_1 - \mathbf{r}_2|} d\mathbf{r}_1 d\mathbf{r}_2\]但缺陷也是很明显的
-
原子核临界区（$r \to 0$），高斯函数 $e^{-\alpha r^2}$ 在 $r=0$ 处的一阶导数严格等于 0，在核中心平滑且圆秃的。但由于原子核对电子有极强的静电吸引力，电子波函数或密度在原子核正中心（$r=0$）会形成一个尖锐的尖峰（Cusp），该处的数学一阶导数不为零。用一组一阶导数为0的圆滑函数，去强行拼凑一个导数不为零的尖锐锥形，在数学上永远会留下一层微小的残差。
-
远离核的无穷远处（$r \to \infty$），不难知电子出现的概率（密度）应该呈指数级缓慢衰减（$\sim e^{-r}$），而高斯函数由于指数的平方项，其衰减速度比真实的指数衰减快得多，导致远区拟合很差。
这些缺陷似乎可以通过增加基组数目来一定程度上缓解。理论上没错，但除了算力上的问题，单纯靠增大基组数量，在两个关键物理边界上收效甚微。
而且，在极其有限的分子空间内塞入越来越多的原始高斯函数时意味着不同的高斯基函数之间会发生严重的重叠。
| 在数学上，这意味着它们的重叠矩阵 $S$（$S_{\mu\nu} = \langle \chi_\mu | \chi_\nu \rangle$）中，很多行向量或列向量会变得几乎一模一样，于是重叠矩阵 $S$ 的本征值中会出现大量极度接近于 0 的数值。这直接导致哈密顿矩阵演变为极端病态矩阵，在求解广义本征值方程 $HC=ESC$ 算 $S^{-1/2}$ 时，计算机的双精度浮点数（16位有效数字）会产生巨大的舍入误差。这会导致自洽场（SCF）迭代剧烈摆动、无法收敛甚至直接数值发散。 |
Why FE Basis Win
现在来正式看这篇论文，会发现KS-Inversion虽然针对的都是一个问题，但思路真的是五花八门。这一篇在数学上更是难啃。
Adjoint Method
为了防止被论文里的数学吓到，先写个背景介绍作为过渡。所谓伴随态法（Adjoint Method），针对的是一类受到偏微分方程或大型代数方程约束的优化问题。具体来说，我们的目标是最小化（或最大化）一个目标泛函 $\mathcal{J}$
\[\min_{\mathbf{m}} \mathcal{J}(\mathbf{u}(\mathbf{m}), \mathbf{m}) \quad \text{subject to} \quad \mathbf{R}(\mathbf{u}, \mathbf{m}) = \mathbf{0}\]其中
-
一堆控制变量 $\mathbf{m}$：这是我们可以主动调节的参数
-
一堆状态变量 $\mathbf{u}$：这是系统被动响应产生的状态
-
一堆约束方程 $\mathbf{R}(\mathbf{u}, \mathbf{m}) = \mathbf{0}$：状态变量和控制变量之间必须满足的物理定律。只要 $\mathbf{m}$ 确定了，通过求解这个方程就能唯一确定 $\mathbf{u}$，即 $\mathbf{u} = \mathbf{u}(\mathbf{m})$。
想要优化当然就要进行梯度下降之类的操作，也就是计算 $\mathcal{J}$ 关于控制变量 $\mathbf{m}$ 的梯度 $\frac{d \mathcal{J}}{d \mathbf{m}}$。
\[\frac{d \mathcal{J}}{d \mathbf{m}} = \frac{\partial \mathcal{J}}{\partial \mathbf{m}} + \frac{\partial \mathcal{J}}{\partial \mathbf{u}} \frac{d \mathbf{u}}{d \mathbf{m}}\]但会发现 $\frac{d \mathbf{u}}{d \mathbf{m}}$这玩意一般是没那么好求的，我们需要对约束方程 $\mathbf{R}(\mathbf{u}, \mathbf{m}) = \mathbf{0}$ 求全导数
\[\frac{\partial \mathbf{R}}{\partial \mathbf{u}} \frac{d \mathbf{u}}{d \mathbf{m}} + \frac{\partial \mathbf{R}}{\partial \mathbf{m}} = \mathbf{0} \implies \frac{d \mathbf{u}}{d \mathbf{m}} = - \left( \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \right)^{-1} \frac{\partial \mathbf{R}}{\partial \mathbf{m}}\]注意这玩意是个矩阵，每个元素都要到约束方程里求一遍梯度，这谁受得了。
这时候就可以用伴随态法取个巧，由于 $\mathbf{R}(\mathbf{u}, \mathbf{m}) = \mathbf{0}$ 严格成立，我们可以在目标泛函上加上这一项，其值保持不变。引入乘子向量（也就是伴随态变量） $\boldsymbol{\lambda}$
\[\mathcal{L}(\mathbf{u}, \mathbf{m}, \boldsymbol{\lambda}) = \mathcal{J}(\mathbf{u}, \mathbf{m}) - \boldsymbol{\lambda}^T \mathbf{R}(\mathbf{u}, \mathbf{m})\]虽然看着和拉格朗日乘子很像，但注意伴随向量 $\boldsymbol{\lambda}$ 是状态变量 $\mathbf{u}$ 和控制变量 $\mathbf{m}$ 的函数，不是一个常数向量。引入的目的不是简单的凑拉格朗日乘子法，且看我操作。
因为 $\mathcal{L} = \mathcal{J}$，所以 $\frac{d \mathcal{J}}{d \mathbf{m}} = \frac{d \mathcal{L}}{d \mathbf{m}}$。我们展开其全导数公式（这里的 $d\boldsymbol{\lambda}$ 项因为 $\mathbf{R}=0$ 而被消掉）
\[\frac{d \mathcal{L}}{d \mathbf{m}} = \frac{\partial \mathcal{J}}{\partial \mathbf{m}} + \frac{\partial \mathcal{J}}{\partial \mathbf{u}} \frac{d \mathbf{u}}{d \mathbf{m}} - \boldsymbol{\lambda}^T \left( \frac{\partial \mathbf{R}}{\partial \mathbf{m}} + \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \frac{d \mathbf{u}}{d \mathbf{m}} \right)\] \[\frac{d \mathcal{L}}{d \mathbf{m}} = \left( \frac{\partial \mathcal{J}}{\partial \mathbf{m}} - \boldsymbol{\lambda}^T \frac{\partial \mathbf{R}}{\partial \mathbf{m}} \right) + {\left( \frac{\partial \mathcal{J}}{\partial \mathbf{u}} - \boldsymbol{\lambda}^T \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \right)} \frac{d \mathbf{u}}{d \mathbf{m}}\]为了彻底消灭掉后面那一项无法计算的 $\frac{d \mathbf{u}}{d \mathbf{m}}$，我们强制让它前面的括号项等于 0，得到
\[\frac{\partial \mathcal{J}}{\partial \mathbf{u}} - \boldsymbol{\lambda}^T \frac{\partial \mathbf{R}}{\partial \mathbf{u}} = \mathbf{0}\]这也就是神奇的伴随方程（Adjoint Equation）：\(\left( \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \right)^T \boldsymbol{\lambda} = \left( \frac{\partial \mathcal{J}}{\partial \mathbf{u}} \right)^T\)
得到了伴随态解 $\boldsymbol{\lambda}(\mathbf{u},\mathbf{m})$，原本全导数公式中耦合项就严格归零。全空间的总梯度公式化简为极其精简的形式
\[\frac{d \mathcal{J}}{d \mathbf{m}} = \frac{\partial \mathcal{J}}{\partial \mathbf{m}} - \boldsymbol{\lambda}^T \frac{\partial \mathbf{R}}{\partial \mathbf{m}}\]由此梯度下降的计算效率得到飞跃。可喜可贺。
正片开始
即使你看明白伴随态法的基本想法，你会发现论文里的公式还是很难懂，，
还是老一套，我们把KS-Inversion的问题，转换为泛函优化问题。论文中采用的泛函是一个带有自定义空间权重的主拟合项，和一个带自定义超参数的正则项（防止过拟合，确保势能平滑）
\[\min_{v_{\text{xc}}} \mathcal{J}[v_{\text{xc}}, \{\psi_i\}] = \mathcal{F}[\{\psi_i\}] + \lambda \mathcal{R}[v_{\text{xc}}]\] \[\mathcal{F}[\{\psi_i\}] = \frac{1}{2} \int_{\Omega} w(\mathbf{r}) \left( \rho(\mathbf{r}) - \rho_{\text{data}}(\mathbf{r}) \right)^2 d\mathbf{r}\] \[\mathcal{R}[v_{\text{xc}}] = \frac{1}{2} \int_{\Omega} |\nabla v_{\text{xc}}(\mathbf{r})|^2 d\mathbf{r}\]对应到前面的理论，这里的交换关联势就是控制变量，而KS轨道波函数自然就是状态变量。
带空间权重 $w(\mathbf{r})$ 是为了平衡不同位置的密度误差的权重，我们知道由于 $\rho_{\text{data}}(\mathbf{r})$ 不可能是真实密度，会由于高斯基底的问题在核位置极其不可信，这部分的权重就会相应地调小，算是一个人工修正。
注意到我们有约束函数
\[\hat{H}[v_{\text{xc}}] \psi_i(\mathbf{r}) = \varepsilon_i \psi_i(\mathbf{r}), \quad i = 1, 2, \dots, N_{\text{occ}}\]其中，单电子哈密顿算符为
\[\hat{H}[v_{\text{xc}}] = -\frac{1}{2}\nabla^2 + v_{\text{ext}}(\mathbf{r}) + v_{\text{H}}(\mathbf{r}) + v_{\text{xc}}(\mathbf{r})\]同时，分子轨道必须满足正交归一化约束
\[\int_{\Omega} \psi_i(\mathbf{r}) \psi_j(\mathbf{r}) d\mathbf{r} = \delta_{ij}\]由于泛函中的状态变量${\psi_i}$是任意的，由此正交归一化约束不能省。
于是为了使用伴随态法进行梯度下降，我们构造广义泛函（吓哭了）
\[\mathcal{L} = \mathcal{J} - \sum_{i} \int_{\Omega} p_i(\mathbf{r}) \left( \hat{H}\psi_i(\mathbf{r}) - \varepsilon_i \psi_i(\mathbf{r}) \right) d\mathbf{r} - \sum_{i,j} \mu_{ij} \left( \int_{\Omega} \psi_i \psi_j d\mathbf{r} - \delta_{ij} \right)\]依葫芦画瓢写出伴随方程（这里进化成了偏微分方程）
\[\left( \hat{H}[v_{\text{xc}}] - \varepsilon_k \right) p_k(\mathbf{r}) = 4 w(\mathbf{r}) \left( \rho(\mathbf{r}) - \rho_{\text{data}}(\mathbf{r}) \right) \psi_k(\mathbf{r}) - 2 \sum_{j=1}^{N_{\text{occ}}} \mu_{kj} \psi_j(\mathbf{r})\] \[\mu_{kj} = 2 \int_{\Omega} w(\mathbf{r}) \left( \rho(\mathbf{r}) - \rho_{\text{data}}(\mathbf{r}) \right) \psi_k(\mathbf{r}) \psi_j(\mathbf{r}) d\mathbf{r}\]解出来后代入，可以得到全空间显式的梯度表达
\[\frac{\delta \mathcal{J}}{\delta v_{\text{xc}}(\mathbf{r})} = \frac{\delta \mathcal{L}}{\delta v_{\text{xc}}(\mathbf{r})} = \sum_{k=1}^{N_{\text{occ}}} \psi_k(\mathbf{r}) p_k(\mathbf{r}) - \lambda \nabla^2 v_{\text{xc}}(\mathbf{r})\]具体的推导细节对我来说还是略显吃力，很多地方还不理解。这里先把方程摆上来再说，关键是体会思路（）
wait，所以基组呢？
差点在数学里陷得太深，忘了这篇文章的主题。这篇文章的另一个亮点，就是交换相关势不再采用高斯基组展开，而是采用有限元基组 $N_A(\mathbf{r})$。
\[v_{\text{xc}}(\mathbf{r}) = \sum_{A=1}^{M} V_A N_A(\mathbf{r})\]所以关于交换相关势的梯度可以直接变成关于系数的梯度
\[G_A = \frac{\partial \mathcal{J}}{\partial V_A} = \sum_{k=1}^{N_{\text{occ}}} \int_{\Omega} \psi_k(\mathbf{r}) p_k(\mathbf{r}) N_A(\mathbf{r}) d\mathbf{r} + \lambda \sum_{B=1}^{M} K_{AB} V_B\]其中 $K_{AB} = \int_{\Omega} \nabla N_A(\mathbf{r}) \cdot \nabla N_B(\mathbf{r}) d\mathbf{r}$ 是有限元理论中经典的全局刚度矩阵（Stiffness Matrix）。
什么，你问我有限元基组 $N_A(\mathbf{r})$ 到底是什么东西？我只能说抽象极了，以后有空再补。
要不干脆把标题改了吧，，真不想写了
以上！
