book/DFT/Density Functional Theory Modeling Mathematical Analysis.md

Density Functional Theory: Modeling, Mathematical Analysis, Computational Methods, and Applications

Preface

基本上来说，DFT 的思路是这样的：

• 使用一组基来表示单电子的波函数。可能使用平面波（对于周期性结构），也可能使用高斯函数（对于分子）。
  • 这个单电子的波函数并不是实际体系的波函数。 也就是说，在实际计算时，我们并不是把这些单电子波函数组合成 Slater 行列式就当成实际波函数来使用的（这样的话误差就太大了）。
• 根据波函数可以计算出电荷密度。
  • 这个电荷密度既是实际体系的电荷密度，也是假想体系的电荷密度。
  • 如果没有对于泛函的近似，那么这个电荷密度应该是精确的。
• 根据电荷密度可以计算出假想体系的哈密顿量。
  • 实际体系的哈密顿量不用算，已经是已知的了。
  • 假想体系中，一部分原本复杂的电子间相互作用被归结到外场中，这个外场的势（哈密顿量的表达式）由电荷密度决定。
• 根据哈密顿量可以再解出单电子波函数，或者计算出下一轮循环中的波函数。当然也可以计算出能量。
• 当两轮循环的能量差小于某个阈值时，就认为计算收敛了，也就是认为现在的电荷密度和波函数都是准确的了（如果不考虑交换关联泛函的近似）。

在实际的体系中，波函数的维度很高（无法被写为单电子波函数的乘积）:

\Psi = \Psi(x_1, x_2, \ldots, x_N)

（这里 x_i 表示包含了自旋的坐标。）

对应的定态哈密顿方程为：

-\frac12 \sum_i \nabla_i^2 + \int V_{\text{ext}}(r) \rho(r) + \sum_{i\ne j} \frac{1}{2|r_i - r_j|} = E \Psi

注意到，波函数对电子的描述有冗余：实际上电子是不可区分的，但波函数实际上描述的是互相之间可以区分的粒子。 这导致波函数需要更多的约束条件才能使得它处在一个合法的空间中，比如交换反对称性。 但无论如何，一旦给出波函数，就给每个电子标上了序号。 在后续计算有关单体和多体之间关系的问题时，可能需要注意，利用波函数的反对称性来简化计算（或者说，把不物理的冗余信息去掉一些）。 例如，显然可以根据波函数计算出“第一个电子”的电荷密度（把其它维度积分掉，然后取模方即可），即：

\rho_1(r_1) = \int |\Psi(x_1, x_2, \ldots, x_N)|^2 \, \dd{\sigma_1} \dd{x}_2 \ldots \dd{x}_N

但这显然是不物理的（什么叫做“第一个电子”呢？）。

因此一般来说，考虑的是所有电子的总的电荷密度，根据对称性，其实也就是上面那个值乘以 $N$； 或者，考虑两个电荷的共同的概率密度（因为库伦势中，每一项最多只涉及到两个电子）。 这个等会儿再说。

之前提到了实际计算时的过程。但在推导理论（而不关心具体的程序实现）时，使用的是另外一个思路：

• 首先，可以证明，由基态电荷密度，可以计算出：
  • 体系中有几个电子。
  • 每个原子核在什么位置，有多少电荷。
• 由这些信息，可以写出完整的薛定谔方程，进一步得到体系的所有信息。
• 也就是说，基态电荷密度本身就包含了体系的所有信息（不只是基态的信息，激发态也包含了）。
• 因此，求解第一性原理计算的问题，往往可以归结为求解基态（即能量最低态）的电荷密度的问题。
• 或者说：给定一个合法的电荷密度，它对应于多个合法的波函数，进一步对应于多个能量值。取最低能量值最低的那个电荷密度，就是我们要找的了。
• 可以想象：如果可以直接写出电荷密度到它对应的那些波函数的最低能量的泛函的表达式，求解这个问题就变得很容易了（求个导就行）。 实际上这个表达式写不出来，如果尽力写的话，就是：

  
    E[\rho] = \int V(r)\rho(r) \dd{r} + F[\rho]
  

  这里并没有多么高超的技术，只是把可以直接求出来（电荷密度可以直接决定，不取决于所取的波函数）的部分拿出来， 剩下的需要折腾波函数的部分归结到 F[\rho] 里面。

1. Review of Approximations for the Exchange-Correlation Energy in Density-Functional Theory

1.1 Basics of Density-Functional Theory

1.1.1 The Many-Body Problem

1.1.2 The Universal Density Functional

1.1.3 The Kohn–Sham Scheme

1.1.3.1 Decomposition of the Universal Functional

要寻找 F[\rho] 的表达式，按照 F[\rho] 的定义，就有这样的思路： 在所有合法的、对应的电荷密度一样的波函数中，找出一个动能和电子-电子相互作用能加起来最小的那个。 按照 Kohn-Sham 的思路，这个事情分为两步走： 第一步，在所有合法的、对应的电荷密度一样的、同时形式是单电子波函数的 Slater 行列式的波函数中， 找出能量最小的那个，并计算出从电荷密度到波函数的表达式、以及从波函数到能量的表达式，这样就可以得到一个从电荷密度到能量的表达式。 第二步，实际上的波函数并不是一个 Slater 行列式，我们再计算出实际上的波函数的能量比 Slater 行列式形式的波函数的能量要少多少， 把这个差值再补进去。

在接下来，使用 \Phi 来表示上述的那个 Slater 行列式形式的波函数，使用 \Psi 来表示实际的波函数。 这样就有：

  F[\rho] = T_{\text s}[\rho] + E_{\text H}[\rho] + E_{\text x}[\rho] + E_{\text c}[\rho]

其中 T_{\text s}[\rho] 是 Slater 波函数的动能：

  T_{\text s}[\rho] = \sum_i \langle \Phi | -\frac{1}{2} \nabla^2 | \Phi \rangle

它其实可以直接由 Slater 行列式中的单电子波函数写出来，它的表达式稍后再讨论。

E_{\text H}[\rho] 是经典的电子间相互作用势能（与具体波函数无关）：

  E_{\text H}[\rho] = \frac{1}{2} \int \int \frac{\rho(r)\rho(r')}{|r - r'|} \dd{r} \dd{r'}

E_{\text x}[\rho] 用来表示经典电子间相互作用势能与 Slater 行列式形式的波函数的电子间相互作用势能之间的差异， 它的表达式（竟然）是可以推导出来的：

  E_{\text x}[\rho] = -\frac{1}{2} \sum_\sigma \int \frac{|\gamma_\sigma(r, r')|^2}{|r - r'|} \dd{r} \dd{r'}

这里 \gamma_\sigma(r, r') 是 Slater 行列式的波函数的密度矩阵（将 N-1 个电子的部分积分去掉后得到的东西）：

  \gamma_\sigma(r, r') = N \int \Phi^*(r, \sigma, x_2, \ldots, x_N) \Phi(r', \sigma, x_2, \ldots, x_N) \dd{x}_2 \ldots \dd{x}_N

TODO: 这个密度矩阵是我理解的那个混态的密度矩阵吗？是不是存在什么简单的表达式，使得它可以直接从单电子波函数中计算出来，而不需要再对 Slater 行列式积分？

E_{\text c}[\rho] 用来表示实际波函数与 Slater 行列式形式的波函数的差异导致的能量变化，它包括动能和电子-电子相互作用两项：

  E_{\text c}[\rho] = \left(\langle \Psi | -\frac{1}{2} \sum_i \nabla_i^2 | \Psi \rangle - T_{\text s}[\rho]\right)
    + \left(\langle \Psi | -\frac{1}{2} \sum_{i\ne j} \frac{1}{|r_i - r_j|} | \Psi \rangle
    - E_{\text x}[\rho] - E_{\text H}[\rho]\right)

到现在为止，看起来好像没问题了，但仔细一想其实还有挺大的问题：

• 如何计算 \Psi 的表达式还没有给出。即使不给出一个严格的表达式，至少应该给出一个迭代求解的算法。
• E_{\text xc}[\rho] 的表达式也没有给出。
• 一旦给出的电荷密度到能量的表达式不是解析的（事实上确实如此），那么就需要一个迭代求解的算法来找到最低能量的电荷密度，而不能直接求导算出来。

如果解决了上面三个问题，那么原则上来说，整个问题就解决了。

1.1.3.2 The Kohn–Sham Equations

接下来考虑单电子波函数的情况，即尽可能将上面的东西都写成组成 \Phi 的单电子波函数的表达式。 一旦总的哈密顿量可以拆解成关于每个单电子的部分，那么就可以把多电子的薛定谔方程写成单电子的薛定谔方程。

我们先假定所有的单电子波函数的自旋都要么朝上要么朝下。在把大致的流程走完一遍之后，再将自旋纳入考虑。

在这之前，我们先重新审视“求基态波函数”这个事情。 它实际上可以看作是用拉格朗日法求极小值。 也就是，在 \langle\varphi|\varphi\rangle = 1 的约束下，求 E[\varphi] = \langle\varphi|H|\varphi\rangle 的极小值。 使用拉格朗日乘子法，那就是解下面这个方程：

  \delta\left(\langle\varphi|H|\varphi\rangle - \lambda(\langle\varphi|\varphi\rangle - 1)\right) = 0

对 \varphi 求偏微分，就可以得到：

  \Re\left(\langle\delta\varphi|H|\varphi\rangle - \lambda\langle\delta\varphi|\varphi\rangle\right) = 0

由 \delta\varphi 的任意性，就可以得到：

  H|\varphi\rangle = \lambda|\varphi\rangle

这样就得到了薛定谔方程，而基态能量就是拉格朗日乘子。

沿着同样的思路，我们来导出关于各个单电子波函数（称为 KS 轨道）的类似于薛定谔方程的那个东西（称为 KS 方程）。分为这样几步：

• 写出由 KS 轨道计算体系能量的表达式。
• 写出 KS 轨道受到的约束（波函数归一化）。
• 使用拉格朗日乘子法，得到要使得变分为零的式子。
• 对这个式子求偏导，令偏导等于零，得到 KS 方程。

第一步，由 KS 轨道计算体系能量的表达式是（把上面的结果抄下来）：

  E[\{\varphi_i\}] = \int V(r)\rho(r) \dd{r} + T_{\text s}[\{\varphi_i\}] + E_{\text H}[\{\varphi_i\}]
    + E_{\text x}[\{\varphi_i\}] + E_{\text c}[\{\varphi_i\}]

其中，电荷密度：

  \rho(r) = \sum_i |\varphi_i(r)|^2

动能项：

  T_{\text s}[\rho] = \sum_i \langle \varphi_i | -\frac{1}{2} \nabla^2 | \varphi_i \rangle

E_{\text H}[\rho] 直接写出来不方便，但一定是可以写出来的。我们等会儿只用到它的变分，所以等会儿再一起写。

E_{\text x}[\rho] 和 E_{\text c}[\rho] 仍然是不知道的。

第二步，KS 轨道的约束其实就是 $\langle\varphi_i|\varphi_i\rangle = 1$。这里有 N 个约束条件。把对应的拉格朗日乘子记作 $\varepsilon_i$。

第三步，使用拉格朗日乘子法，知道我们需要求解：

  \delta\left(\int V(r)\rho(r) \dd{r} + T_{\text s}[\{\varphi_i\}] + E_{\text H}[\{\varphi_i\}]
    + E_{\text x}[\{\varphi_i\}] + E_{\text c}[\{\varphi_i\}]
    - \sum_i \varepsilon_i(\langle\varphi_i|\varphi_i\rangle - 1)\right) = 0

第四步，要求括号里关于 \varphi_i 的变分为零。 对于其中可以直接由 \rho 表示的项， 按照 \frac{\delta F[\{\varphi_i\}]}{\delta \varphi_i} = \frac{\delta F[\rho]}{\delta \rho} \frac{\delta \rho}{\delta \varphi_i} 的道理， 我们可以先求括号内关于 \rho 的变分，再求 \rho 关于 \varphi_i 的变分，而后者很容易计算：

  \frac{\delta \rho}{\delta \varphi_i} = 2\varphi_i^*

仔细计算后可以得到：

  \left(V(r) -\frac12 \nabla^2 + \int\frac{\rho(r')}{|r - r'|}\dd{r'}
    + \frac{\delta V_{\text x}}{\delta \rho} + \frac{\delta V_{\text c}}{\delta \rho}\right) \varphi_i
    = \varepsilon_i \varphi_i

注意到，经典的电子间相互作用势能项前面的 1/2 被去掉了，这个咋一看有点奇怪，仔细算一下就知道是对的。

记 $v_{\text x} = \frac{\delta V_{\text x}}{\delta \rho}$，$v_{\text c} = \frac{\delta V_{\text c}}{\delta \rho}$， 那么就得到了 KS 方程：

  \left(-\frac12 \nabla^2 + V(r) + \int\frac{\rho(r')}{|r - r'|}\dd{r'} + v_{\text x}[\rho] + v_{\text c}[\rho]\right) \varphi_i
    = \varepsilon_i \varphi_i

这个方程给出了求解 KS 轨道（以及进一步得到 \Phi 和电荷密度等信息）的方法：求解前面那个算符的特征向量。

TODO: 为什么接下来会选择能量最低的几个轨道作为被占据的轨道？轨道能量与总能量之间的关系是什么？

1.1.3.3 Extension to Spin Density-Functional Theory

如果去掉之前“自旋要么朝上要么朝下”的约束，那么自旋维度上，单电子的波函数也可能会和空间维度上那样互相耦合， 也就是单个电子的自旋会处在一个混态或者叠加态上。 （有证明提到，这时其实仍然可以不考虑带自旋的电荷密度，但我们先不管它。）

TODO: 这时，不带自旋的基态电荷密度是否仍然可以决定一切？如果不可以的话，为什么分上下自旋的电荷密度就可以？

我们先按照这样一个假定来继续推导：分上下自旋的电荷密度就可以决定一切。 在这之后，所有原本的推导中，关于电荷密度的部分，就需要拆成两个部分。 因为总的电荷密度恰好是上下自旋的电荷密度之和，所以对于一些项不需要拆分。

仔细推导后可以发现，在 KS 方程中，需要改变的地方并不多（如果没有外加磁场的话），只有 v_{\text x} 和 v_{\text c} 的形式需要修改。 对于交换能，还可以推导得到一个简洁的结果：

  E_{\text x}[\rho_\uparrow, \rho_\downarrow] =  \frac12 (E_{\text x}[2\rho_\uparrow] + E_{\text x}[2\rho_\downarrow])

将电荷密度拆分为上下自旋的电荷密度的另外一个好处是，在之后拟合交换关联泛函的时候，可以得到更精确的结果。

1.1.4 The Generalized Kohn–Sham Scheme

TODO: 关于加入更多自由度来更好地拟合交换关联泛函的问题，用到了再看吧。

接下来要解决的问题：

• 如何近似表示交换关联能？
• 程序实现上，通常来说是怎样做的？