Post

Study Note：Kohn–Sham DFT

Posted
By Bane Dysta
39 min read
Study Note：Kohn–Sham DFT

在HF 推导笔记中，我们从分子的能量组成出发，用单个 Slater 行列式表示多电子波函数，再通过能量变分求出轨道。接下来想了解密度泛函理论（density functional theory，DFT），一个自然的问题是：既然它以电子密度为基本变量，为什么实际计算时仍然需要轨道？

Kohn–Sham（KS）方法正是为了解决这个问题。密度让我们可以换一种方式描述基态能量；轨道则让其中很难直接用密度计算的动能变得可以处理。下面先从电子密度能够告诉我们什么讲起，再逐步写出能量表达式，最后推导决定这些轨道的 KS 方程。

本文讨论固定核构型下、非相对论、零温的基态问题。轨道推导采用无自旋极化的闭壳层体系与实空间轨道，所有量均使用原子单位。

1 电子密度与能量

在 Born–Oppenheimer（BO）近似下，可以先固定原子核的位置，求电子在这套核构型下的状态。电子能量依旧有三种来源：电子运动的动能、电子与核之间的吸引能，以及电子之间的排斥能。核间排斥能由核构型决定，在这次电子计算中是一个常数，最后再加上即可。

多电子波函数需要同时描述所有电子的坐标。电子密度提出的问题更简单：空间中某个位置附近，平均有多少电子？无论体系有多少电子，这个量都只需要指定一个三维位置。

把电子数密度对全部空间积分，应当得到电子总数。因此，电子密度满足：

\[\rho(\mathbf r)\geq 0,\qquad\int\rho(\mathbf r)\,d\mathbf r=N\]

式中：

  • $\rho(\mathbf r)$：位置 $\mathbf r$ 处的电子数密度，表示该位置附近单位体积内的平均电子数。
  • $\mathbf r$：电子的三维空间位置。
  • $N$：体系的电子总数。
  • $\int\cdots d\mathbf r$：把相应的量对整个三维空间积分。

知道电子密度以后，有一部分能量已经可以直接计算。原子核的位置与电荷都已给定，它们在每个位置产生的核吸引势能为：

\[v_{\mathrm{ext}}(\mathbf r)=-\sum_A\frac{Z_A}{\lvert\mathbf r-\mathbf R_A\rvert}\]

式中：

  • $v_{\mathrm{ext}}(\mathbf r)$：一个电子在位置 $\mathbf r$ 处受到的外势，这里是所有原子核产生的吸引势能；下标 $\mathrm{ext}$ 表示来自电子体系外部的势。
  • $A$：原子核的编号。
  • $Z_A$：第 $A$ 个原子核的电荷数，例如氢核为 $1$。
  • $\mathbf R_A$：第 $A$ 个原子核的位置。
  • $\lvert\mathbf r-\mathbf R_A\rvert$：电子位置与这个原子核之间的距离。
  • $\sum_A$：将所有原子核的贡献相加。

某处的外势乘以该处的电子密度，再对空间积分，就得到全部电子的核吸引能：

\[E_{\mathrm{ext}}[\rho]=\int v_{\mathrm{ext}}(\mathbf r)\rho(\mathbf r)\,d\mathbf r\]

式中：

  • $E_{\mathrm{ext}}[\rho]$：电子在给定外势中的总势能，这里就是电子与核的吸引能。
  • 方括号 $[\rho]$：表示这个量依赖整份电子密度分布。给它一份密度，就可以算出一个能量数值，这种关系称为泛函（functional）。

外势能已经有了明确的密度表达式。接下来的问题是：电子动能和电子之间的排斥能，也能放到一个以密度为变量的理论中吗？

2 密度变分原理

Hohenberg–Kohn（HK）定理给出了基态 DFT 的出发点。先按非简并基态理解：对于固定电子数、固定电子间相互作用形式的体系，基态密度确定外势，允许相差一个只改变能量零点的常数。外势确定以后，电子哈密顿量及其基态也随之确定。因此，基态密度中包含了确定基态性质所需的信息。HK 原始论文建立了这一关系及相应的密度变分原理。

为了看清能量泛函怎样定义，可以用约束搜索（constrained search）的方式来想：先给定一份密度，在所有能够产生这份密度的合法多电子波函数中，选出动能与电子排斥能之和最低的那个。

\[F[\rho]=\min_{\Psi\to\rho}\langle\Psi\vert\hat T+\hat V_{\mathrm{ee}}\vert\Psi\rangle\]

式中：

  • $F[\rho]$：在给定密度下能够达到的最低内部电子能量，包含动能与电子间排斥能。
  • $\Psi$：描述全部电子的多电子波函数，要求对所有电子位置与自旋取值计算的总概率为 $1$，并满足交换电子坐标时变号的反对称性。
  • $\Psi\to\rho$：只在能够产生指定密度 $\rho$ 的波函数中寻找。
  • $\min$：取所有允许选择中的最低值。
  • $\hat T$：全部电子的动能算符。
  • $\hat V_{\mathrm{ee}}$：全部电子对之间的库仑排斥算符；下标 $\mathrm{ee}$ 表示电子与电子之间的作用。
  • 字母上方的帽号：表示算符，即作用于波函数的数学操作。
  • $\langle\Psi\vert\hat T+\hat V_{\mathrm{ee}}\vert\Psi\rangle$：在波函数 $\Psi$ 下计算动能与电子排斥能的期望值，也就是这个态对应的平均内部能量。

同一份密度对应的外势能已经由上一节的积分确定，因此加上它，就得到了在这个外势下的电子能量泛函：

\[E[\rho]=F[\rho]+E_{\mathrm{ext}}[\rho]\]

式中：

  • $E[\rho]$：在已给定的外势下，这份密度对应的电子能量泛函。

上面的内部能量泛函对不同原子、分子具有同样的定义，因为电子的动能算符和库仑排斥形式相同；不同核构型的差别进入外势能。这就是这里所谓的普适泛函（universal functional）。约束搜索的具体表述可参看 Levy 的论文。

现在让密度也发生变化。对每一份允许的密度，内部都已经选取了最低能量波函数；再在这些密度之间取最低值，就相当于在全部允许波函数中寻找基态：

\[E_0=\min_{\rho}E[\rho]=E[\rho_0]\]

式中：

  • $E_0$：给定外势下的精确基态电子能量。
  • $\rho_0$：达到这一最低能量的基态电子密度。
  • 此处的 $\min_{\rho}$：在能够由合法、归一化的 $N$ 电子波函数产生、且能量有限的密度中寻找，电子数保持不变。

这就建立了密度的变分问题。不过，定义出一个泛函，还需要找到具体的计算办法。约束搜索的定义中仍然包含多电子波函数，直接按它计算并不轻松。KS 方法接下来要做的，是把一大部分能量改写成可以用单电子轨道计算的形式。

3 Kohn–Sham 参考体系

3.1 同密度的轨道表示

先想一个计算起来简单的体系：让电子之间没有显式的两两相互作用，每个电子都在同一个有效势中运动。这样的无相互作用体系可以用单电子轨道描述，动能也容易逐轨道计算。

KS 方法寻找这样的参考体系，使它在基态给出的电子密度与真实相互作用体系相同。以下采用常见的双占据、闭壳层参考体系，并假定所讨论的目标密度能由相应的无相互作用局域势基态表示。

每个占据空间轨道中有两个电子，因此参考体系的电子密度为：

\[\rho(\mathbf r)=2\sum_{i=1}^{n_{\mathrm{occ}}}\phi_i(\mathbf r)^2,\qquad n_{\mathrm{occ}}=N/2\]

式中：

  • $\phi_i(\mathbf r)$：第 $i$ 个占据 KS 空间轨道，用来构造参考体系的电子密度；本文采用实轨道。
  • $i,j$：占据空间轨道的编号，取值为 $1,\ldots,n_{\mathrm{occ}}$。
  • $n_{\mathrm{occ}}$：占据空间轨道数；下标 $\mathrm{occ}$ 是 occupied 的缩写，表示被电子占据。
  • $\phi_i(\mathbf r)^2$：一个电子处于这个实轨道时的概率密度。
  • 前面的 $2$：每个占据空间轨道贡献两个电子。

这些轨道满足正交归一条件：

\[\int\phi_i(\mathbf r)\phi_j(\mathbf r)\,d\mathbf r=\delta_{ij}\]

式中：

  • $\delta_{ij}$：Kronecker delta，两个编号相同时为 $1$，不同时为 $0$。这里表示每个轨道自身的重叠为 $1$，不同轨道的重叠为 $0$。

KS 轨道在这里承担的是表示密度和计算参考动能的任务。真实体系的相互作用与相关效应，将通过完整的能量泛函进入后面的计算。

3.2 无相互作用动能

在能够产生给定密度的参考行列式中，取动能最低的一个，将相应动能记作无相互作用动能。用它的占据轨道计算：

\[T_s[\rho]=2\sum_i^{\mathrm{occ}}\int\phi_i(\mathbf r)\hat t\phi_i(\mathbf r)\,d\mathbf r,\qquad\hat t=-\frac{1}{2}\nabla^2\]

式中：

  • $T_s[\rho]$：给定密度对应的无相互作用参考动能；下标 $s$ 表示这个单电子参考体系。
  • $\hat t$：一个电子的动能算符。
  • $\nabla^2$：对三个空间坐标分别求二阶导数后相加，称为拉普拉斯算符。
  • $\sum_i^{\mathrm{occ}}$：对全部占据空间轨道求和，后文的双重占据轨道求和也使用这一范围。

到这里，轨道重新出现的原因就清楚了：它们为参考动能提供了一个可计算的表达式。接下来，将真实内部能量与这份参考动能之间的差别单独收集起来，就能形成 KS 的能量分解。

4 Kohn–Sham 能量分解

4.1 Hartree 能

给定密度以后，还可以先计算一个直观的电子排斥能：将电子密度看成连续分布的电荷，按每两小块电荷之间的距离计算库仑排斥，再将所有位置对的贡献相加。

每一对位置会被交换顺序计算两次，因此要乘上二分之一：

\[E_H[\rho]=\frac{1}{2}\iint\frac{\rho(\mathbf r)\rho(\mathbf r')}{\lvert\mathbf r-\mathbf r'\rvert}\,d\mathbf r\,d\mathbf r'\]

式中：

  • $E_H[\rho]$：Hartree 能，即电子密度按连续电荷分布计算的经典库仑排斥能。
  • $\mathbf r’$：另一个独立的空间积分位置，撇号用于与 $\mathbf r$ 区分。
  • $\lvert\mathbf r-\mathbf r’\rvert$：两个空间位置之间的距离。
  • $\iint\cdots d\mathbf r\,d\mathbf r’$：分别对两个三维位置积分。

这部分可以直接由密度计算。真实电子间作用中的交换、相关，以及真实动能与参考动能之间的差别，还需要在能量表达式中补齐。

4.2 交换关联能

将内部能量中参考动能和 Hartree 能以外的部分统一定义为交换关联能：

\[E_{\mathrm{xc}}[\rho]=F[\rho]-T_s[\rho]-E_H[\rho]\]

式中：

  • $E_{\mathrm{xc}}[\rho]$：交换关联能，即完整内部能量减去无相互作用动能与经典库仑能后剩下的部分。
  • 下标 $\mathrm{xc}$：由 exchange 与 correlation 的缩写组成，分别表示交换和相关。

将这个定义按动能与电子排斥能重新展开，可以看得更具体：

\[E_{\mathrm{xc}}[\rho]=\{T[\rho]-T_s[\rho]\}+\{V_{\mathrm{ee}}[\rho]-E_H[\rho]\}\]

式中：

  • $T[\rho]$：在定义 $F[\rho]$ 的约束搜索中，最低内部能量波函数所具有的真实动能。
  • $V_{\mathrm{ee}}[\rho]$：同一个波函数所具有的电子间排斥能。
  • 第一对花括号：真实动能与参考动能之间的差别。
  • 第二对花括号：真实电子排斥能与经典密度库仑能之间的差别。

所以，选择无相互作用参考体系之后，真实体系的动能差也进入了交换关联能。通过这个定义，内部能量的各项能够完整地拼回去。

将它们与外势能相加，得到 KS 的电子能量表达式：

\[\boxed{E[\rho]=T_s[\rho]+E_{\mathrm{ext}}[\rho]+E_H[\rho]+E_{\mathrm{xc}}[\rho]}\]

这四项的来源分别是参考轨道的动能、核对电子的吸引、电子密度的经典库仑排斥，以及补齐剩余内部能量的交换关联项。采用精确交换关联泛函时，这是一种精确的能量分解；实际计算需要为交换关联泛函选择具体的近似。

Kohn 与 Sham 的原始论文从密度变分原理构造了相应的自洽单电子方程。下面把从能量到方程的求导过程展开。

5 Kohn–Sham 方程的推导

5.1 有效势

要找到能量最低的密度，需要先知道：密度在某个位置改变一点，会让能量怎样变化？普通函数用导数描述一个变量变化时的影响；泛函则用一整个位置函数记录不同位置的密度变化所带来的影响，这就是泛函导数（functional derivative）。

以 Hartree 能为例，将它的一阶变化写成：

\[\delta E_H=\int v_H(\mathbf r)\delta\rho(\mathbf r)\,d\mathbf r,\qquad v_H(\mathbf r)=\frac{\delta E_H}{\delta\rho(\mathbf r)}\]

式中：

  • $\delta\rho(\mathbf r)$：位置 $\mathbf r$ 处电子密度的一点微小变化。
  • $\delta E_H$：由这些密度变化引起的 Hartree 能的一阶变化，只保留与变化量成正比的部分。
  • $v_H(\mathbf r)$：Hartree 势，是密度在这个位置变化时，对 Hartree 能一阶变化的贡献系数。
  • $\delta E_H/\delta\rho(\mathbf r)$：Hartree 能对密度的泛函导数，整个空间中的这些变化系数组成 Hartree 势。
  • 写在一个量前的 $\delta$：表示取一阶变化；带两个编号的 $\delta_{ij}$ 仍表示前面定义的 Kronecker delta。

Hartree 能包含两个密度因子，它们都要参与变化：

\[\delta E_H=\frac{1}{2}\iint\frac{\delta\rho(\mathbf r)\rho(\mathbf r')+\rho(\mathbf r)\delta\rho(\mathbf r')}{\lvert\mathbf r-\mathbf r'\rvert}\,d\mathbf r\,d\mathbf r'\]

在第二项中互换两个积分位置，两项就完全相同。它们相加后消去前面的二分之一：

\[\delta E_H=\int\left[\int\frac{\rho(\mathbf r')}{\lvert\mathbf r-\mathbf r'\rvert}\,d\mathbf r'\right]\delta\rho(\mathbf r)\,d\mathbf r\]

对照一阶变化的定义，直接读出：

\[v_H(\mathbf r)=\int\frac{\rho(\mathbf r')}{\lvert\mathbf r-\mathbf r'\rvert}\,d\mathbf r'\]

这正是整份电子密度在指定位置产生的经典排斥势能。外势能对密度是线性的，它的一阶变化系数就是外势本身。对交换关联能也作同样的定义：

\[v_{\mathrm{xc}}(\mathbf r)=\frac{\delta E_{\mathrm{xc}}}{\delta\rho(\mathbf r)},\qquad\delta E_{\mathrm{xc}}=\int v_{\mathrm{xc}}(\mathbf r)\delta\rho(\mathbf r)\,d\mathbf r\]

式中：

  • $v_{\mathrm{xc}}(\mathbf r)$：交换关联势，描述密度在不同位置发生微小变化时，交换关联能随之变化的程度。
  • $\delta E_{\mathrm{xc}}$：交换关联能的一阶变化。

将这三个势收集起来，定义 KS 有效势：

\[v_s(\mathbf r)=v_{\mathrm{ext}}(\mathbf r)+v_H(\mathbf r)+v_{\mathrm{xc}}(\mathbf r)\]

式中：

  • $v_s(\mathbf r)$：KS 参考体系的有效势，由外势、Hartree 势与交换关联势共同组成。

电子能量的一阶变化便可以写成：

\[\delta E=\delta T_s+\int v_s(\mathbf r)\delta\rho(\mathbf r)\,d\mathbf r\]

式中：

  • $\delta E$：电子总能量的一阶变化。
  • $\delta T_s$：参考动能的一阶变化。

到这里，三个势的来由都已经明确：它们分别是三个势能项对密度变化的响应。下一步把密度的变化换成轨道的变化，就能得到决定轨道的方程。

5.2 轨道变分

密度由占据轨道的平方和组成，因此改变轨道时，密度的一阶变化为：

\[\delta\rho(\mathbf r)=4\sum_i^{\mathrm{occ}}\phi_i(\mathbf r)\delta\phi_i(\mathbf r)\]

式中：

  • $\delta\phi_i(\mathbf r)$：第 $i$ 个实空间轨道的微小变化。
  • 系数 $4$：双占据贡献一个 $2$，对轨道平方取一阶变化再贡献一个 $2$。

参考动能也是轨道的二次型。假定轨道及其变化在远处衰减，使分部积分的边界项消失，动能算符作用于左右轨道所产生的两项可以合并为：

\[\delta T_s=4\sum_i^{\mathrm{occ}}\int\delta\phi_i(\mathbf r)\hat t\phi_i(\mathbf r)\,d\mathbf r\]

把这两式代入电子能量的一阶变化，得到：

\[\delta E=4\sum_i^{\mathrm{occ}}\int\delta\phi_i(\mathbf r)[\hat t+v_s(\mathbf r)]\phi_i(\mathbf r)\,d\mathbf r\]

轨道可以变化，但每个轨道都要保持归一化，不同轨道之间还要正交。因此，引入拉格朗日乘子，将能量和这些约束一起处理：

\[\mathcal L=E[\rho]-2\sum_{ij}^{\mathrm{occ}}\varepsilon_{ij}\left[\int\phi_i(\mathbf r)\phi_j(\mathbf r)\,d\mathbf r-\delta_{ij}\right]\]

式中：

  • $\mathcal L$：将电子能量与轨道正交归一约束合在一起的拉格朗日函数。
  • $\varepsilon_{ij}$：为每一项约束引入的待定系数，组成一个实对称的拉格朗日乘子矩阵。
  • $\sum_{ij}^{\mathrm{occ}}$：让两个编号各自遍历所有占据空间轨道后求和。
  • 方括号中的差：实际轨道重叠偏离正交归一要求的程度；满足约束时为零。
  • 约束项前的 $2$：乘子的定义约定，使后面的乘子本征值直接对应通常的轨道能量。

轨道重叠中左右两个轨道也都要取变化。利用乘子矩阵的对称性，将它们合并后：

\[\delta\mathcal L=4\sum_i^{\mathrm{occ}}\int\delta\phi_i(\mathbf r)\left\{[\hat t+v_s(\mathbf r)]\phi_i(\mathbf r)-\sum_j^{\mathrm{occ}}\varepsilon_{ji}\phi_j(\mathbf r)\right\}\,d\mathbf r\]

式中：

  • $\delta\mathcal L$：拉格朗日函数的一阶变化。

在驻点处，这个变化对任意微小轨道变化都应为零，因此每个轨道对应的花括号必须为零：

\[[\hat t+v_s(\mathbf r)]\phi_i(\mathbf r)=\sum_j^{\mathrm{occ}}\varepsilon_{ji}\phi_j(\mathbf r)\]

此时乘子矩阵还可能包含非对角元素。可以在占据轨道之间作正交线性组合，使这个实对称矩阵对角化。这样的组合保持占据空间以及轨道平方和不变，因此密度、能量与有效势也保持不变。

在这组正则轨道下，每个轨道就分别满足：

\[\boxed{\left[-\frac{1}{2}\nabla^2+v_{\mathrm{ext}}(\mathbf r)+\int\frac{\rho(\mathbf r')}{\lvert\mathbf r-\mathbf r'\rvert}\,d\mathbf r'+v_{\mathrm{xc}}(\mathbf r)\right]\phi_i(\mathbf r)=\varepsilon_i\phi_i(\mathbf r)}\]

式中：

  • $\varepsilon_i$：正则轨道 $i$ 对应的乘子本征值，称为 KS 轨道能量。

这就是闭壳层形式的 Kohn–Sham 方程。左侧依次包含一个电子的动能、核吸引、由电子密度产生的 Hartree 势，以及交换关联势。我们从四项能量出发，对后三项取密度变化、对动能取轨道变化，最终把它们组合成了这个单电子方程。

6 交换关联泛函近似

6.1 局域密度近似

KS 方程已经写出，实际计算中还需要给交换关联能选择一个具体表达式。局域密度近似（local density approximation，LDA）的做法是：在每个位置，用与该处密度相同的均匀电子气，估计这里每个电子的交换关联能。

将每个电子的这份能量乘以当地电子密度，再对空间积分：

\[E_{\mathrm{xc}}^{\mathrm{LDA}}[\rho]=\int\rho(\mathbf r)\epsilon_{\mathrm{xc}}^{\mathrm{unif}}(\rho(\mathbf r))\,d\mathbf r\]

式中：

  • $E_{\mathrm{xc}}^{\mathrm{LDA}}$：采用 LDA 计算的交换关联能。
  • $\epsilon_{\mathrm{xc}}^{\mathrm{unif}}(\rho)$：密度为 $\rho$ 的均匀电子气中，每个电子的交换关联能。
  • 上标 $\mathrm{unif}$：uniform 的缩写，表示空间中密度处处相同的均匀电子气。
  • $\epsilon$：这里采用的每电子能量符号，与前面表示轨道乘子的 $\varepsilon$ 区分。

因为积分中的能量密度只取决于同一位置的密度，对它求导只需使用普通函数的乘积求导法则：

\[v_{\mathrm{xc}}^{\mathrm{LDA}}(\mathbf r)=\epsilon_{\mathrm{xc}}^{\mathrm{unif}}(\rho(\mathbf r))+\rho(\mathbf r)\left.\frac{d\epsilon_{\mathrm{xc}}^{\mathrm{unif}}}{d\rho}\right|_{\rho=\rho(\mathbf r)}\]

式中：

  • $v_{\mathrm{xc}}^{\mathrm{LDA}}$：LDA 交换关联能对应的交换关联势。
  • $d\epsilon_{\mathrm{xc}}^{\mathrm{unif}}/d\rho$：每电子交换关联能随均匀电子气密度变化的普通导数。
  • 竖线及其下标：先对均匀电子气的密度变量求导，再代入当前空间位置的密度值。

这给出了一个具体例子：选定能量泛函以后，交换关联势通过求导得到，两者需要配套使用。

6.2 广义梯度近似

实际分子的密度在空间中并不均匀。广义梯度近似（generalized gradient approximation，GGA）进一步让能量密度依赖电子密度在附近变化的快慢，即密度梯度：

\[E_{\mathrm{xc}}^{\mathrm{GGA}}[\rho]=\int f(\rho(\mathbf r),\nabla\rho(\mathbf r))\,d\mathbf r\]

式中：

  • $E_{\mathrm{xc}}^{\mathrm{GGA}}$：采用 GGA 计算的交换关联能。
  • $f$：给定密度及其梯度以后，算出单位体积交换关联能的函数。
  • $\nabla\rho$：密度沿三个空间方向的变化率组成的向量，称为密度梯度。

密度改变时，它的梯度也会改变，因此一阶变化有两部分：

\[\delta E_{\mathrm{xc}}^{\mathrm{GGA}}=\int\left[\frac{\partial f}{\partial\rho}\delta\rho+\frac{\partial f}{\partial(\nabla\rho)}\cdot\nabla\delta\rho\right]\,d\mathbf r\]

式中：

  • $\partial f/\partial\rho$：固定梯度时，能量密度随电子密度变化的偏导数。
  • $\partial f/\partial(\nabla\rho)$：分别对梯度的三个分量求偏导后组成的向量。
  • $\nabla\delta\rho$：微小密度变化本身的空间梯度。
  • 两个向量之间的中点：对应分量相乘后相加，即向量的点积。

将第二项分部积分，并令密度变化在边界处消失，就能把梯度从密度变化上移走，再读出泛函导数：

\[v_{\mathrm{xc}}^{\mathrm{GGA}}(\mathbf r)=\frac{\partial f}{\partial\rho}-\nabla\cdot\frac{\partial f}{\partial(\nabla\rho)}\]

式中：

  • $v_{\mathrm{xc}}^{\mathrm{GGA}}$：GGA 交换关联能对应的交换关联势，右侧各量在位置 $\mathbf r$ 处计算。
  • $\nabla\cdot$：散度操作，即对向量的三个分量分别沿相应空间方向求导后相加。

这个势仍然是 KS 方程中乘在轨道上的一个位置函数。对于交换关联能的更一般形式，还需要按它依赖的变量继续求导。例如，常见杂化泛函含有依赖占据轨道的精确交换项，其通常的实现会出现非局域交换算符，进入广义 Kohn–Sham（generalized Kohn–Sham，GKS）形式。Seidl 等人的论文讨论了这种推广。

7 有限基组表示

到这里，KS 方程仍然是关于轨道函数的微分方程。与 HF 计算一样，可以选定一组已知基函数，用待定系数将轨道展开，把求函数的问题转成求矩阵本征向量的问题。

\[\phi_p(\mathbf r)=\sum_{\mu=1}^{M}C_{\mu p}\chi_\mu(\mathbf r)\]

式中：

  • $p$：任意空间轨道的编号，可以对应占据轨道或未被占据的虚轨道。
  • $M$：基函数总数。
  • $\chi_\mu$：第 $\mu$ 个预先选定的空间函数，用来通过线性组合近似待求轨道。这些函数常以原子核为中心，因此习惯上称为 AO 基函数；AO 是 atomic orbital，即原子轨道的缩写。
  • $\mu,\nu$：基函数的编号，取值为 $1,\ldots,M$。
  • $C_{\mu p}$：基函数 $\mu$ 在线性组合成轨道 $p$ 时的展开系数。

将 KS 方程投影到各个基函数上，定义重叠矩阵与 KS 算符矩阵：

\[S_{\mu\nu}=\int\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)\,d\mathbf r,\qquad H_{\mu\nu}^{\mathrm{KS}}=\int\chi_\mu(\mathbf r)[\hat t+v_s(\mathbf r)]\chi_\nu(\mathbf r)\,d\mathbf r\]

式中：

  • $S_{\mu\nu}$：两个 AO 基函数的重叠积分，组成重叠矩阵 $S$。
  • $H_{\mu\nu}^{\mathrm{KS}}$：KS 单电子算符在这两个 AO 基函数之间的矩阵元，组成矩阵 $H^{\mathrm{KS}}$。

有限基组中的轨道方程便写成：

\[H^{\mathrm{KS}}C=SC\varepsilon,\qquad C^TSC=I\]

式中：

  • $C$：轨道展开系数矩阵，每一列对应一个轨道。
  • $\varepsilon$：在这里表示全部正则轨道能量排成的对角矩阵，维数为 $M\times M$。
  • 上标 $T$：矩阵转置，即把行、列互换。
  • $I$：单位矩阵，对角元素为 $1$，其余元素为 $0$。
  • 第二个等式：在有重叠的 AO 基组中，MO 轨道仍满足正交归一条件；MO 是 molecular orbital，即分子轨道的缩写。

由占据轨道更新密度矩阵与电子密度：

\[P=2C_{\mathrm{occ}}C_{\mathrm{occ}}^T,\qquad\rho(\mathbf r)=\sum_{\mu\nu}P_{\mu\nu}\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)\]

式中：

  • $C_{\mathrm{occ}}$：只保留占据轨道各列的展开系数矩阵。
  • $P$：包含双占据因子 $2$ 的 AO 密度矩阵。
  • $P_{\mu\nu}$：电子密度展开后，基函数乘积 $\chi_\mu\chi_\nu$ 前面的系数。
  • $\sum_{\mu\nu}$：两个编号分别遍历全部基函数后求和。

因此，KS 与 HF 可以采用相似的矩阵求解框架。它们各自的能量表达式决定了所需构造的有效单电子算符。对 LDA、GGA，交换关联积分通常还需要在数值格点上计算。

8 自洽求解与总能量

8.1 SCF 迭代

Hartree 势和交换关联势依赖密度，密度又由待求的占据轨道产生。于是，求解轨道方程以后还要检查：新轨道给出的密度，是否与构造这次有效势所用的密度一致？

从一个密度初猜出发，可以反复进行：

\[\rho^{(k)}\;\longrightarrow\;v_s[\rho^{(k)}]\;\longrightarrow\;\{\phi_i^{(k+1)}\}\;\longrightarrow\;\rho^{(k+1)}\]

式中：

  • 上标 $(k)$、$(k+1)$：当前轮与下一轮迭代的编号。
  • $\rho^{(k)}$：当前轮作为输入的电子密度。
  • $v_s[\rho^{(k)}]$：用这份密度构造的 KS 有效势。
  • ${\phi_i^{(k+1)}}$：在这个势下求出的新占据轨道集合。
  • $\rho^{(k+1)}$：由新占据轨道重新计算的密度。

输入密度与输出密度相互一致时，就达到自洽。这就是自洽场（self-consistent field，SCF）求解。闭壳层基态计算通常占据最低的若干个空间轨道，直到电子总数得到满足；实际迭代还常配合密度混合、DIIS 等收敛技术。

8.2 能量计算

求出轨道能量以后，电子总能量仍按前面的四项能量表达式计算。为了看清轨道能量求和与总能量之间的关系，将每个正则轨道的 KS 方程乘以该轨道并积分，再对双占据轨道求和：

\[2\sum_i^{\mathrm{occ}}\varepsilon_i=T_s+E_{\mathrm{ext}}+2E_H+\int\rho(\mathbf r)v_{\mathrm{xc}}(\mathbf r)\,d\mathbf r\]

这里用到了密度乘以 Hartree 势的积分等于两倍 Hartree 能。与四项能量表达式对照，就得到：

\[\boxed{E[\rho]=2\sum_i^{\mathrm{occ}}\varepsilon_i-E_H[\rho]+E_{\mathrm{xc}}[\rho]-\int\rho(\mathbf r)v_{\mathrm{xc}}(\mathbf r)\,d\mathbf r}\]

轨道能量之和中，Hartree 相互作用计入了两次，需要减去一次。交换关联部分在轨道能量之和中表现为密度乘势的积分，因此还需要用真正的交换关联能将它替换。这条求和关系对应本文推导的局域乘法势 KS 方程，在自洽的正则轨道下使用。

最后加上固定核构型的核间排斥能，就得到这套构型的总能量：

\[E_{\mathrm{tot}}=E[\rho]+E_{\mathrm{nuc}},\qquad E_{\mathrm{nuc}}=\sum_{A<B}\frac{Z_AZ_B}{R_{AB}}\]

式中：

  • $E_{\mathrm{tot}}$：固定核构型下的总能量。
  • $E_{\mathrm{nuc}}$：核间库仑排斥能，在电子 SCF 迭代中保持不变。
  • $B$：另一个原子核的编号，$A<B$ 保证每一对核只计算一次。
  • $Z_B$：第 $B$ 个原子核的电荷数。
  • $R_{AB}$：第 $A$、$B$ 个原子核之间的距离。

9 推导回顾

最初的问题是怎样用密度来计算基态能量。HK 变分原理给出了以密度寻找基态的依据，KS 参考体系则提供了逐轨道计算动能的办法。在这个基础上，把剩余内部能量定义进交换关联项，就得到了完整的四项能量表达式。

对 Hartree 能和交换关联能取密度变化，分别得到相应的势；再把密度变化换成轨道变化，并加入正交归一约束，就得到 KS 方程。整个求解最终由下面三条关系连接起来：

\[E[\rho]=T_s[\rho]+E_{\mathrm{ext}}[\rho]+E_H[\rho]+E_{\mathrm{xc}}[\rho]\] \[\left[-\frac{1}{2}\nabla^2+v_{\mathrm{ext}}(\mathbf r)+v_H(\mathbf r)+v_{\mathrm{xc}}(\mathbf r)\right]\phi_i(\mathbf r)=\varepsilon_i\phi_i(\mathbf r)\] \[\rho(\mathbf r)=2\sum_i^{\mathrm{occ}}\phi_i(\mathbf r)^2\]

第一条给出要极小化的能量；第二条决定产生这个密度的参考轨道；第三条由占据轨道重新构造密度，使前后可以自洽。这也回答了开头的问题：密度是基态理论的基本变量，KS 轨道让参考动能和密度更新成为可以实际求解的计算步骤，而交换关联泛函承接了参考体系与真实相互作用体系之间剩余的能量差别。

参考资料

Quantum Chemistry
DFT Kohn-Sham Electronic Structure SCF Theory
This post is licensed under CC BY 4.0 by the author.
Total hits!