密度泛函理论平面波基组展开

密度泛函理论：平面波基组展开

基本理论与概念

密度泛函理论 (Density Functional Theory, DFT) 是一种研究多电子体系电子结构的量子力学方法。密度泛函理论在物理、 化学、材料上都有广泛的应用, 特别是用来研究分子和凝聚态 的性质, 是凝聚态物理和计算化学领域最常用的方法之一。

误差的几个方面：

• 物理模型的建立（模型误差, 测量误差: 更妙的物理）
• 数值分析（截断误差：更好的数值方法）
• 计算机程序 (舍入误差: 更长的字长）

原胞

原胞的概念, 就是在完全定义一个无限扩展的周期性材料时, 含有必要的最少原子个数的超晶胞。

绝热近似、 HF 方法

我们一般用用简谐近似描述核的运动。玻恩–奥本海默近似考虑核不动。

HF 近似，哈密顿量拆成单电子求和，波函数拆成乘积的形式。某个(与其它电子相互作用的)电子的运动可以近似地看成是它在其他电子的平均密度所产生的电磁场中的运动 (平均场理论)。

$✓$ 多电子的薛定谔方程可以化简为 $\mathrm{N}$ 个单电子薛定谔方程
$✓$ 多电子问题化为单电子问题！
$✓$ Coulomb项或Hartree项需要自洽求解。

DFT 与 KS 方程

DFT：总能量是密度的泛函！！！总能量对密度的变分极小就是体系的基态。动能用密度不能准确云表达，但用波函数可以。
考虑薛定谔方程，
$$\left[-\frac{h^2}{2m}\sum _{i=1}^N{\nabla }_i^2+\sum _{i=1}^{N}V\left({\mathbf{r}}_i\right)+\sum _{i=1}^N\sum _{j<i}U\left({\mathbf{r}}_i,{\mathbf{r}}_j\right)\right]\psi =E\psi$$
其中 m 为电子质量。括号中的三项依次表示电子的动能、电子和原子核之间的作用能、不同电子之间的作用能。$\psi$ 是电子波函数，它依赖于所有电子的空间坐标，有多少电子，全波函数就有三倍函数的自变量。$E$ 是基态能量。

当电子数量很多的时候（N 个），维数如此之高，计算量如此之大，这时候，我们可以考虑把电子波函数表达成 N 个单个电子波的乘积，即 Hartree 乘积，
$$\psi ={\psi }_1(r){\psi }_2(r)\cdots {\psi }_N(r)$$
对于波函数而言，N 个电子在特定坐标系下出现的概率为，
$${\psi }^{\ast }\left({\mathbf{r}}_1,\cdots ,{\mathbf{r}}_N\right)\psi \left({\mathbf{r}}_1,\cdots ,{\mathbf{r}}_N\right)$$
空间中某个具体位置上的电荷密度 $n(\mathbf{r})$ 可以用单电子波函数的形式写出来，
$$n(\mathbf{r})=2\sum _i{\psi }_i^{\ast }(\mathbf{r}){\psi }_i(\mathbf{r})$$
此式中, 针对电子系统所占据的全部单电子波函数进行求和。因此, 求 和项中的表达式就是在单电子波函数 ${\psi }_i(\mathbf{r})$ 中一个电子位于 $\mathbf{r}$ 处的概率值。表达式中出现因子 2 是因为电子具有自旋, 且泡利不相容原理 表明:每个单电子波函数能够被不同自旋的两个电子所占据。

Schrödinger 方程得到的基态能量是电荷密度的唯一函数（泛函）。使整体泛函最小化的电荷密度就是对应于 Schrödinger 方程完全解的真实电荷密度。

能量泛函可以写为，
E [ { ψ i } ] = E known  [ { ψ i } ] + E x c [ { ψ i } ] E\left[\left\{\psi_{i}\right\}\right]=E_{\text {known }}\left[\left\{\psi_{i}\right\}\right]+E_{\mathrm{xc}}\left[\left\{\psi_{i}\right\}\right] E[{ψi​}]=Eknown ​[{ψi​}]+Exc​[{ψi​}]
其中, 将泛函分开为能够写成简单解析形式的一项 E k n o w n [ { ψ i } ] E_{\mathrm{known}}\left[\left\{\psi_{i}\right\}\right] Eknown​[{ψi​}]，
E known  [ { ψ i } ] = − h 2 m ∑ i ∫ ψ i ∗ ∇ 2 ψ i   d 3 r + ∫ V ( r ) n ( r ) d 3 r + e 2 2 ∬ n ( r ) n ( r ′ ) ∣ r − r ′ ∣ d 3 r   d 3 r ′ + E i o n \begin{aligned} E_{\text {known }}\left[\left\{\psi_{i}\right\}\right]=&-\frac{h^{2}}{m} \sum_{i} \int \psi_{i}^{*} \nabla^{2} \psi_{i} \mathrm{~d}^{3} r+\int V(\boldsymbol{r}) n(\boldsymbol{r}) \mathrm{d}^{3} r \\ &+\frac{\mathrm{e}^{2}}{2} \iint \frac{n(\boldsymbol{r}) n\left(\boldsymbol{r}^{\prime}\right)}{\left|\boldsymbol{r}-\boldsymbol{r}^{\prime}\right|} \mathrm{d}^{3} r \mathrm{~d}^{3} r^{\prime}+E_{\mathrm{ion}} \end{aligned} Eknown ​[{ψi​}]=​−mh2​i∑​∫ψi∗​∇2ψi​ d3r+∫V(r)n(r)d3r+2e2​∬∣r−r′∣n(r)n(r′)​d3r d3r′+Eion​​
以及其他所有部分 $E_{\mathrm{X}\mathrm{C}}$。
在上式右侧, 依次分别是: 电子的动能, 电子和原子核之间的库伦作用, 电子之间的库伦作用, 原子核之间的库伦作用。对于 $V_{\mathrm{x}\mathrm{c}}$ 能量泛函完整表达式中的另一项 E X C [ { ψ i } ] E_{\mathrm{XC}}\left[\left\{\psi_{i}\right\}\right] EXC​[{ψi​}], 是交换关联泛函, 它所定义的是没有 包括在 “ known” 这一项中所有其他的量子力学效应。

事实上，求解正确的电荷密度，可以表示为求解一套方程，而其中每个方程都只与一个电子相关。
Kohn - Sham 方程的表达式为
$$\left[-\frac{h^2}{m}{\nabla }^2+V(\mathbf{r})+V_{\mathrm{H}}(\mathbf{r})+V_{\mathrm{X}\mathrm{C}}(\mathbf{r})\right]{\psi }_i(\mathbf{r})={\epsilon }_i{\psi }_i(\mathbf{r})$$

第一项表示一个电子和所有原子核之间的相互作用。第二个也称为 Hatree 势能,可以写为
$$V_{\mathrm{H}}(\mathbf{r})={\mathrm{e}}^2\int \frac{n\left({\mathbf{r}}^{\prime }\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}{\mathrm{d}}^3{r}^{\prime }$$

这个势能描述的是一个 Kohn - Sham 方程所考虑的单个电子, 与该问题中全部电子所产生的总电荷密度之间的库伦排斥作用。$V_{XC}$ 为交换关联能的 “泛函导数”, 即
$$V_{\mathrm{X}\mathrm{C}}(\mathbf{r})=\frac{\delta E_{\mathrm{X}\mathrm{C}}(\mathbf{r})}{\delta n(\mathbf{r})}$$

自洽迭代

在以上有关 Kohn - Sham 方程的讨论中, 如果你感到似乎是陷人 了一个循环, 那么这个感觉是对的。为了求解 Kohn-Sham 方程, 需要 确定 Hatree 势能; 而为了得到 Hatree 势能, 又需要知道电荷密度; 但为 了找到电荷密度, 又必须知道单电子波函数方程; 为了知道这些波函 数, 又必须求解 Kohn-Sham 方程。为了打破这一循环, 这个问题通常 用迭代算法来处理, 其过程简述如下:
（1）定义一个初始的、尝试性的电荷密度 $n(\mathbf{r})$ 。
(2) 求解由尝试性的电荷密度所确定的 Kohn - Sham 方程, 得到 单电子波函数 ${\psi }_i(\mathbf{r})$ 。
（3）计算由第(2)步 Kohn - Sham 单粒子波函数所确定的电荷密 度, 即 $n_{\mathrm{K}\mathrm{S}}(\mathbf{r})=2\sum {\sum }_i^{\ast }(\mathbf{r}){\psi }_i(\mathbf{r})$ 。
（4）比较计算得到的电荷密度 $n_{\mathrm{K}\mathrm{S}}(\mathbf{r})$ 和在求解 Kohn-Sham 方程 时使用的电荷密度 $n(\mathbf{r})$ 。如果两个电荷密度相同, 则这就是基态电荷 密度, 并可将其用于计算总能。如果两个电荷密度不同, 则用某种方式 对尝试性电荷密度进行修正, 然后再从第(2)步重新开始。

这就是自洽求解的过程。

交换关联势

交换关联能 $E_{xc}[\rho ]$ 包含了总能量表达式前几项不能精确表达的部分, 包含了:

• 电子交换能 $K_{ij}$
• 电子关联能
• 电子动能的不准确部分 $(T_0[\rho ]$ 不同于真 实动能 $T_e[\rho ]$ )
• 库伦势导致的自相互作用修正(在HF中,仅 当 $K_{ii}=J_{ii}$, 自相互作用抵消)。

交换关联势是交换关联能的变分导数。从自洽迭代过程中，计算看起来很容易，但是交换关联势的确定是不容易的。实际上,我们并不清楚交换关联泛函的真实形式, 尽管 Hohenberg Kohn 定理肯定它是确定存在的。幸运的是, 对于均匀电子气这种情形, 该泛函可以直接导出。在此情形下, 电荷密度在空间所有点上都是常 数, 即 $n(\mathbf{r})=$ 常数。对于任何真实材料而言, 这种情形的意义可能不 大, 其原因在于: 正是由于电荷密度的变化才确定了化学键, 也才使材料更有意义。但均匀电子气给出了实际使用 Kohn - Sham 方程的可行方法。为了做到这一点, 我们把每个位置的交换关联势能都设定为已知的交换关联势能, 这个已知的交换关联能是根据该位置所观测到的电荷密度, 由均匀电子气得到的, 即
$$V_{\mathrm{X}\mathrm{C}}(\mathbf{r})=V_{\mathrm{X}\mathrm{C}}^{\text{electron gas }}[n(\mathbf{r})]$$
这一近似仅仅使用了局域密度来确定近似的交换关联泛函, 因此称作局域密度近似(Local Density Approximation,LDA)。LDA 使我们能够完 全确定地写出 Kohn -Sham 方程, 但需要注意的是: 用这些方程所得的 结果并不能严格求解真实的 Schrödinger 方程, 因为我们并没有使用真 实的交换关联泛函。

$$E_{xc}[\rho ]=E_x[\rho ]+E_c[\rho ]$$

• 一般形式:
  $E_{xc}[\rho ]=\int d\vec{r}{\epsilon }_{xc}\left[\rho ,\nabla \rho ,{\nabla }^2\rho ,\cdots \right]\rho (\vec{r})$

• 广义梯度近似（GGA）:
  $E_{xc}[\rho ]=\int d\vec{r}{\epsilon }_{xc}[\rho ,\nabla \rho ]\rho (\vec{r})$

• 局域密度近似（LDA）:
  $E_{xc}[\rho ]=\int d\vec{r}{\epsilon }_{xc}[\rho ]\rho (\vec{r})$

  • 每一个小区域内的电子密度是常数,可以用均匀电子气描述
  • 均匀电子气的交换关联能可以算出
  • 整个体系的交换关联能是各个小区域的交换关联能的代数和

现今最广泛使用的局域势是: B88 (Becke, 1988)、PW91(Perdew & Wang, 1991)、P86 (Perdew, 1986), LYP (Lee at al, 1988)。

交换关联泛函关于导数是不连续的。

DFT本质上是一个面向基态的理论。 因此, 采用LDA及GGA的DFT框架对于半导体 的带隙的预测非常差。采用GW方法后, 才得到了合理的结果(最近采用metaGGA泛函的DFT, 对于绝缘体和半导体的带隙预言已用GGA的计算量达到了GW的精度)。在激发态的计算中, 采用LDA的DFT还不足以与半经验的、及修正的 $\mathrm{H}\mathrm{F}$ 方法竞争(新近发展的TDDFT, 及采用求解 Bethe-Salpter Equation的方法可以很好地计算 光吸收谱)。

赝势

赝势是:

• 针对原子设计的或构造的, 替换原子的全 电子势的一种近似势, 可减少基底的数目, 从而减少计算量
• 这种近似势把内壳层电子隐去了, 只考虑 价电子且价电子的應波函数没有节点
• 可以部分考虑相对论效应

赝势改变的是外势项。

DFT 应用

纳米金属催化剂的催化效果取决于纳米颗粒的表面反应活性，而表面反应活性整体上可以看作是纳米颗粒形状和每种原子反应活性的的复杂函数。那么，能否在金属纳米颗粒的形状尺寸，与作为合成氨反应催化剂的活性之间简历直接联系？事实上，这个过程总的化学反应由金属催化剂表面上的 12 个单独步骤完成的。这些步骤的速率很大程度上上取决于所涉及金属原子的区域配位关系。
在金属铜当中加入一些比它大的原子，容易使它脆化。原因在于，铜的晶粒和晶粒之间有个晶界，杂质原子更容易聚集在晶界处，改变了晶界的电子结构。通过电子结构计算，我们可以得出脆化的原因。
不同温度、压强下的物质化学结构会发生改变，包括原子间的组合关系，分子的结构等等，通过这，我们可以预测行星的一些材料性质。我们关注的是不同相对空间位置下原子的能量。原子是由原子核和电子构成，而原子核的体量远比电子大得多，所以电子对于环境变化的相应是远远大于原子核的。我们定义一个表示能量最低的基态能量函数，它的自变量是各个原子核的位置，我们也称之为绝热势能面。知道了绝热势能面，我们也就知道材料的能量将如何变化。

赝势平面波 DFT 方法

各项能量表达

单粒子的薛定谔方程:
$$H{\Psi }_{\mathbf{k}}(r)=E(\mathbf{k}){\Psi }_{\mathbf{k}}(r)$$
波函数的展开式为
$${\Psi }_{\mathbf{k}}(r)=\sum _{\mathbf{G}}{c}_G{\Phi }_G^{\mathbf{k}}(\mathbf{r})$$
总能量为(vasp 的实际计算种还包含核-核库伦相互作用)
$$E_{\text{tot }}=T+E_{ext}+E^{coul}+E_{xc}(\mathbf{r})+E_{nn},$$
其中动能为
$$T=\sum _j\int d\mathbf{r}{\Phi }_j^{\ast }(\mathbf{r})\left(-\frac{1}{2}{\nabla }^2\right){\Phi }_j(\mathbf{r}),$$
电子库伦项为
$$E^{\text{coul }}=\frac{1}{2}\iint d\mathbf{r}d{\mathbf{r}}^{\prime }\frac{\rho (\mathbf{r})\rho \left({\mathbf{r}}^{\prime }\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|},$$
外势项写为
$$E_{ext}=\sum _{j,\mu ,l}\iint d{\mathbf{r}}^{\prime }d\mathbf{r}({\Phi }_j^f(\mathbf{r})\sum _{m!}{Y}_{lm}(\mathbf{r})V_l(\left(\mathbf{r}-{\mathbf{R}}_{\mu l}\mid \right)Y_{lm}^{\ast }\left({\mathbf{r}}^{\prime }\right){\Phi }_j\left({\mathbf{r}}^{\prime }\right)$$
交换关联能为
$$\begin{array}{rl}{E}_{xc}& =\int C_{xc}\rho (\mathbf{r})d\mathbf{r},\\ {\hat{V}}_{xc}(\mathbf{r})& =\frac{\delta E_{xc}}{\delta \rho (\mathbf{r})}=C_{xc}+\rho \frac{\partial C_{xc}}{\partial \rho }\end{array}$$
核-核库伦项为
$$E_{nn}=\frac{1}{2}\sum _{\mu \ne v}\frac{Z_{\mu }{Z}_v}{\left|{\mathbf{R}}_{\mu }-{\mathbf{R}}_v\right|},$$
其中密度为 $\rho (\mathbf{r})={\sum }_j{\Phi }_j^{\ast }(\mathbf{r}){\Phi }_j(\mathbf{r})$.

Fourier 展开

波函数可以找一个基组展开，
$$H\Psi =E\Psi$$

$$\Psi =\sum C_{\alpha I}^k{\Phi }_{\beta I}^k$$

$$\sum \left[<{\Phi }_{\alpha I}^k|H|{\Phi }_{\beta J}^k>-E(k)<{\Phi }_{\alpha I}^k\mid {\Phi }_{\beta J}^k>\right]C_{\beta J}^k=0$$

$$\sum \left[H_{\alpha I,\beta J}^k-E(k)S_{\alpha I,\beta J}^k\right]C_{\beta J}^k=0$$

周期函数对于描述体相材料非常有用, 因为至少对于无缺陷的材料而言, 电荷密度和波函数实际上都是空间周期性函数。傅里叶展开正是利用了这种周期性质。
$$f(\mathbf{r})=\sum _{\mathbf{G}}f(\mathbf{G})e^{i\mathbf{G}\cdot \mathbf{r}}$$
$$\begin{array}{rl}f(\mathbf{G})& =\frac{1}{\Omega }\int f(\mathbf{r})e^{-i\mathrm{G}\cdot \mathbf{r}}d\mathbf{r}\\ f(-\mathbf{G})& =\frac{1}{\Omega }\int f(\mathbf{r})e^{i\mathrm{G}\cdot \mathbf{r}}d\mathbf{r}\end{array}$$
$$\frac{1}{\Omega }\int e^{i\left(\mathbf{G}-{\mathbf{G}}^{\prime }\right)\cdot \mathbf{r}}d\mathbf{l}={\delta }_{\mathrm{G},{\mathrm{G}}^{\prime }}$$

波函数的动量空间表达式:
$${\Phi }_j^{\mathbf{k}}(\mathbf{r})=\sum _{\mathbf{G}}{\Phi }_j(\mathbf{k}+\mathbf{G})e^{i(\mathbf{k}+\mathbf{G})\mathbf{r}},$$
密度的动量空间表达式:
$$\rho (\mathbf{r})=\sum _{\mathbf{G}}\rho (\mathbf{G})e^{i\mathbf{G}\mathbf{r}},$$

动能项
$$\begin{array}{rl}T& =\sum _j\int d\mathbf{r}{\Phi }_j^{\ast }(\mathbf{r})\left(-\frac{1}{2}{\nabla }^2\right){\Phi }_j(\mathbf{r})\\ & =\frac{1}{2}\Omega \sum _{j,\mathbf{k},\mathbf{G}}{\left|{\Phi }_j(\mathbf{k}+\mathbf{G})\right|}^2|\mathbf{k}+\mathbf{G}{|}^2\end{array}$$
其中，利用到
$$\begin{array}{rl}& {\Phi }_j=\sum _G{\Phi }_j(\vec{k}+\vec{G})e^{i(\vec{k}+\vec{G})\cdot \vec{r}}\\ & \int e^{i\left(\vec{G}-{\vec{G}}^{\prime }\right)\cdot \vec{r}}d\vec{r}=\Omega {\delta }_{\vec{G}{\vec{G}}^{\prime }}\end{array}$$
交换关联项
$${ϵ}_{xc}(\mathbf{r})=\sum _{\mathbf{G}}{ϵ}_{xc}(\mathbf{G})e^{i\mathbf{G}\mathbf{r}}$$
$${\hat{V}}_{xc}(\mathbf{r})=\sum _{\mathbf{G}}{V}_{xc}(\mathbf{G})e^{i\mathbf{G}\mathbf{r}}$$
$$E_{xc}=\int {ϵ}_{xc}\rho (\mathbf{r})d\mathbf{r}=\Omega \sum _{\mathbf{G}}{ϵ}_{xc}(\mathbf{G}){\rho }^{\ast }(\mathbf{G}),$$
库伦项
利用 $\int \rho (r)e^{iGr}dr=\Omega {\rho }^{\ast }(G)$,
$$\rho (\mathbf{r})=\sum _{\mathbf{G}}\rho (\mathbf{G})e^{i\mathrm{G}\mathrm{r}},$$
这里的 $\rho (\mathbf{G})={\rho }^{\ast }(-\mathbf{G})$。
$$\begin{array}{rl}{V}_{coul}(\mathbf{r})& =\int \frac{\rho \left({\mathbf{r}}^{\prime }\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}d{\mathbf{r}}^{\prime }\\ & \equiv \sum _{\mathbf{G}}{V}_{coul}(\mathbf{G})e^{i\mathrm{G}\mathrm{r}}\\ & =\sum _{\mathbf{G}}\frac{4\pi \rho (\mathbf{G})}{|\mathbf{G}{|}^2}{e}^{i\mathrm{G}\mathrm{r}}\end{array}$$
$$\begin{array}{rl}\int \frac{\rho \left({\mathbf{r}}^{\prime }\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}d{\mathbf{r}}^{\prime }& =\sum \int \frac{\rho (\mathbf{G})}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}{e}^{i\mathrm{G}\cdot r^{\prime }}d{\mathbf{r}}^{\prime }\\ & =-\sum \rho (\mathbf{G})\int \frac{1}{|\mathbf{x}|}{e}^{i\mathrm{G}\cdot (\mathbf{r}-\mathbf{x})}d\mathbf{x}\\ & =-\sum \rho (\mathbf{G})e^{i\mathrm{G}\cdot \mathrm{r}}\int \frac{1}{|\mathbf{x}|}{e}^{-i\mathrm{G}\cdot \mathbf{x}}d\mathbf{x},\\ & =\sum \rho (\mathbf{G})\frac{4\pi }{\mid {\mathbf{G}}^2}{e}^{i\mathrm{G}\cdot \mathrm{r}}(\text{ if }G\ne 0).\end{array}$$

外势场项
$$\begin{array}{rl}{E}_{\text{ext }}& =\sum _{j,\mu ,l}\iint d{\mathbf{r}}^{\prime }d\mathbf{r}{\Phi }_j^{\ast }(\mathbf{r})\sum _m{Y}_{lm}(\mathbf{r})V_l\left(\left|\mathbf{r}-{\mathbf{R}}_{\mu }\right|\right)Y_{lm}^{\prime }\left({\mathbf{r}}^{\prime }\right){\Phi }_j\left({\mathbf{r}}^{\prime }\right)\\ & =\Omega \sum _{j,\mathbf{k},,,{\mathrm{G}}_{,}{\mathrm{G}}^{\prime }}{\Phi }_j^{\ast }(\mathbf{G}){\Phi }_j\left({\mathbf{G}}^{\prime }\right)S\left({\mathbf{G}}^{\prime }-\mathbf{G}\right)V_{l,\mathbf{k},\mathrm{G},{\mathrm{G}}^{\prime }}^{ps}\end{array}$$
其中 $S$ 为结构因子 (structure factor)
$$S\left({\mathbf{G}}^{\prime }-\mathbf{G}\right)=\frac{1}{N}\sum _{\mu }{e}^{i\left({\mathrm{G}}^{\prime }-\mathrm{G}\right)\cdot {\mathrm{R}}_{\mu }}$$

如果体系含有多种类原子, 没有简单的结构因子,相应项应变为
$$S\left({\mathbf{G}}^{\prime }-\mathbf{G}\right)V_{l,\mathrm{k},\mathrm{G},{\mathrm{G}}^{\prime }}^{p^{\prime s}}\to \frac{1}{N}\sum _{\mu }{e}^{i\left({\mathrm{G}}^{\prime }-\mathrm{G}\right)\cdot {\mathrm{R}}_{\mu }{V}_{l,\mathrm{k},\mathrm{G},{\mathrm{G}}^{\prime }}^{\mu ,ps}}$$
其中
$$\begin{array}{rl}{V}_{l,\mathbf{k},\mathbf{G},{\mathrm{G}}^{\prime }}^{ps}=& \frac{1}{{\Omega }_{at}}(2l+1)4\pi P_l(\cos \gamma )\\ & \int V_l(r)j_l(|\mathbf{k}+\mathbf{G}|r)j_l\left(\left|\mathbf{k}+{\mathbf{G}}^{\prime }\right|r\right)r^{2}dr\end{array}$$

$$\cos \gamma =\frac{(\mathbf{k}+\mathbf{G})\cdot \left(\mathbf{k}+{\mathbf{G}}^{\prime }\right)}{\left|\mathbf{k}+\mathbf{G}\Vert \mathbf{k}+{\mathbf{G}}^{\prime }\right|},$$
其中 $j_1$ 和 $P_l$ 是spherical Bessel 及Legendre polynomials。

如果 $V_l(r)=V^{\text{core }}(r)+V_l^{\mathrm{N}\mathrm{L}}(r)$
$$\begin{array}{rl}{E}_{ext}=\Omega \sum _{\mathbf{G}}S(\mathbf{G})V^{\text{core }}(\mathbf{G}){\rho }^{\ast }(\mathbf{G})& \\ & +\Omega \sum _{i,l,\mathbf{G},{\mathbf{G}}^{\prime }}{\Phi }_i^{\ast }(\mathbf{G}){\Phi }_i\left({\mathbf{G}}^{\prime }\right)S\left({\mathbf{G}}^{\prime }-\mathbf{G}\right)V_{l,\mathbf{k},\mathbf{G},{\mathbf{G}}^{\prime }}^{NL}\end{array}$$

动量空间总能量及发散问题

考虑动量空间总能量及发散问题
$$\begin{array}{rl}{E}_{\text{tot }}=& \Omega \{\sum _{j,\mathbf{G}}\frac{1}{2}{\left|{\Phi }_j(\mathbf{k}+\mathbf{G})\right|}^2|\mathbf{k}+\mathbf{G}{|}^2\\ & +\sum _{\mathbf{G}}\frac{1}{2}{V}_{\text{coul }}(\mathbf{G}){\rho }^{\ast }(\mathbf{G})\\ & +\sum _{G}S(\mathbf{G})V^{\text{core }}(\mathbf{G}){\rho }^{\ast }(\mathbf{G})\\ & +\sum _{j,l,{\mathbf{G}}^{\prime }}{\Phi }_j^{\ast }(\mathbf{k}+\mathbf{G}){\Phi }_j\left(\mathbf{k}+{\mathbf{G}}^{\prime }\right)S\left({\mathbf{G}}^{\prime }-\mathbf{G}\right)V_{l,{\mathbf{k}}_j,\mathbf{G},{\mathbf{G}}^{\prime }}^{NL}\\ & +\sum _G{ϵ}_{xc}(\mathbf{G}){\rho }^{\ast }(\mathbf{G})\}+\frac{1}{2}\sum _{\mu \ne v}\frac{Z_{\mu }{Z}_v}{\left|{\mathbf{R}}_{\mu }-{\mathbf{R}}_v\right|}.\end{array}$$

考虑极限作用下的库仑项和外势项
密度展开 $\rho =\frac{1}{\Omega }\sum \left(Z_{\mu }+{\beta }_{\mu }|\mathbf{G}{|}^2\Omega \right)+O^3(\mathbf{G})$
$$\begin{array}{rl}{V}_{\text{coul }}(\mathbf{G})& =\frac{4\pi \rho (\mathbf{G})}{|\mathbf{G}{|}^2}\\ & =\uparrow \frac{4\pi }{\Omega }\sum _{\mu }\left(\frac{Z_{\mu }}{|\mathbf{G}{|}^2}+{\beta }_{\mu }\Omega \right)+O^1(\mathbf{G})\\ V^{\text{core }}(\mathbf{G})& =\frac{1}{\Omega }\sum _{\mu }\left(-\frac{4\pi Z_{\mu }}{|\mathbf{G}{|}^2}+{\alpha }_{\mu }\Omega \right)+O^1(\mathbf{G})\end{array}$$
电子电子相互作用，
$$E_{nn}=\frac{1}{2}\sum _{\mu \ne v}\frac{Z_{\mu }{Z}_v}{\left|{\mathbf{R}}_{\mu }-{\mathbf{R}}_v\right|}.$$
求和在实空间中趋向于无穷，因此当 $k\to 0$ 在 $\mathrm{k}$ 空间中是发散的。因为整个系统是电中性的，所以这些发散必须被不同项抵消。我们把 $E_{nn}$ 分成两项，在 $k$ 空间中计算的发散项，在实空间中计算的其他项。

离子相互作用项
$E_{nn}=\frac{1}{2}{\sum }_{\mu ,v}\frac{4\pi }{\Omega }{\sum }_{\mathrm{G}\ne 0}\frac{Z_{\mu }{Z}_V}{|\mathbf{G}{|}^2}{e}^{-\frac{|G{|}^2}{4x^2}}{e}^{i\mathrm{G}\cdot \left({\mathbf{R}}_{\mu }-{\mathbf{R}}_{\nu }\right)}+\frac{1}{2}{\sum }_{\mu \ne v}{Z}_{\mu }{Z}_V\frac{\mathrm{erfc}\left(|{\mathbf{R}}_{\mu }-{\mathbf{R}}_v|\kappa \right)}{|{\mathbf{R}}_{\mu }-{\mathbf{R}}_V|}$
$-\frac{1}{2}{\sum }_{\mu ,v}\frac{4\pi }{\Omega }{Z}_{\mu }{Z}_v\frac{1}{4{\kappa }^2}-\sum Z_{\mu }^2\frac{\kappa }{\sqrt{\pi }}+{\lim }_{\mathrm{G}\to 0}\frac{1}{2}{\sum }_{\mu ,v}\frac{4\pi }{\Omega }{Z}_{\mu }{Z}_v\frac{1}{|\mathbf{G}{|}^2}$
$E_{nn}\equiv {\gamma }_{\text{Ewald }}+{\lim }_{\mathbf{G}\to 0}\frac{1}{2}{\sum }_{\mu ,v}\frac{4\pi }{\Omega }{Z}_{\mu }{Z}_v\frac{1}{|\mathbf{G}{|}^2}$
其中的 ${\gamma }_{\text{Ewald }}$ 是所有的 $G\ne 0$ 项.

三个发散项求和
$$\begin{array}{rl}& \underset{\mathbf{G}\to 0}{\lim }\Omega \{\sum _{\mathrm{G}}\frac{1}{2}{V}_{\text{coul }}(\mathbf{G}){\rho }^{\ast }(\mathbf{G})\\ & +\sum _{G}S(\mathbf{G})V^{\text{rore }}(\mathbf{G}){\rho }^{\ast }(\mathbf{G})\}\\ & +\frac{1}{2}\sum _{\mu \ne v}\frac{Z_{\mu }{Z}_v}{\left|{\mathbf{R}}_{\mu }-{\mathbf{R}}_v\right|}\end{array}$$
$$\begin{array}{rl}=& \underset{\mathbf{G}\to 0}{\lim }\{\frac{1}{2}\sum _{\mathbf{G}}\sum _{\mu }\left(\frac{4\pi Z_{\mu }}{|\mathbf{G}{|}^2}+4\pi {\beta }_{\mu }\Omega \right){\rho }^{\ast }(\mathbf{G})\\ & +\sum _{\mathbf{G}}S(G)\sum _{\mu }\left(-\frac{4\pi Z_{\mu }}{|\mathbf{G}{|}^2}+{\alpha }_{\mu }\Omega \right){\rho }^{\ast }(\mathbf{G})\}\\ & +{\gamma }_{\text{Ewald }}+\frac{1}{2}\underset{\mathbf{G}\to 0}{\lim }\sum _{\mu ,v}\frac{4\pi Z_{\mu }{Z}_V}{\Omega |\mathbf{G}{|}^2}\end{array}$$
注意, ${\lim }_{G\to 0}S(G)\to 1$, 及系数的 $1/2$ 问题。

$$\begin{array}{rl}=& 2\pi \sum _v{\beta }_v\sum _{\mu }{Z}_{\mu }-\frac{1}{2}\underset{\mathbf{G}\to 0}{\lim }\sum _{\mu ,v}\frac{4\pi Z_V{Z}_{\mu }}{\Omega |\mathbf{G}{|}^2}\\ & -\frac{1}{2}4\pi \sum _v{Z}_{\mu }+\sum {\alpha }_v\sum Z_{\mu }\\ & +{\gamma }_{\text{Ewald }}+\frac{1}{2}\underset{\mathbf{G}\to 0}{\lim }\sum _{\mu ,v}\frac{4\pi Z_{\mu }{Z}_v}{\Omega |\mathbf{G}{|}^2}\\ =& \sum _v{\alpha }_v\sum Z_{\mu }+{\gamma }_{\text{Ewald }}.\end{array}$$

动量空间 KS 方程

Kohn-Sham 方程
$$\begin{array}{c}\left[-\frac{1}{2}{\nabla }^2+{\hat{V}}_{eff}(\mathbf{r})\right]{\Phi }_j(\mathbf{r})={ϵ}_j{\Phi }_j(\mathbf{r}),\\ {\hat{V}}_{eff}(\mathbf{r})=\sum _{ju}{V}_l\left(\left|\mathbf{r}-{\mathbf{R}}_{il}\right|\right)F_l+V_{\text{coul }}(\mathbf{r})+{\hat{V}}_{xc}(\mathbf{r})\end{array}$$
先把波函数用平面波展开，
$$\begin{array}{rl}& \sum _{\mathbf{G}}[-\frac{1}{2}{\nabla }^2+\sum _{\mu ,l}{V}_l\left(\left|\mathbf{r}-{\mathbf{R}}_{\mu }\right|\right){\mathcal{P}}_l\\ & +V_{\text{coul }}(\mathbf{r})+{\hat{V}}_{xc}(\mathbf{r})]{\Phi }_j(\mathbf{k}+\mathbf{G})e^{i(\mathbf{k}+\mathbf{G})\mathbf{r}}\\ =& \sum _{\mathbf{G}}{ϵ}_j{\Phi }_j(\mathbf{k}+\mathbf{G})e^{i(\mathbf{k}+\mathbf{G})\mathbf{r}}\end{array}$$

动能项，
$$\begin{array}{rl}& \int e^{-i\left(\mathbf{k}+{\mathbf{G}}^{\prime }\right)\mathbf{r}}\sum _{\mathbf{G}}\left[-\frac{1}{2}{\nabla }^2\right]{\Phi }_j(\mathbf{k}+\mathbf{G})e^{i(\mathbf{k}+\mathbf{G})\mathbf{r}}d\mathbf{r}\\ =& \sum _{\mathbf{G}}\frac{1}{2}|\mathbf{k}+\mathbf{G}{|}^2{\Phi }_j(\mathbf{k}+\mathbf{G})\int e^{-i\left(\mathbf{k}+{\mathbf{G}}^{\prime }\right)\mathbf{r}}{e}^{i(\mathbf{k}+\mathbf{G})\mathbf{r}}d\mathbf{r}\\ =& \Omega \sum _{\mathbf{G}}\frac{1}{2}{\left|\mathbf{k}+{\mathbf{G}}^2\right|}_j(\mathbf{k}+\mathbf{G}){\delta }_{\mathbf{G},{\mathbf{G}}^{\prime }}\end{array}$$
库伦项和交换关联项
$$\begin{array}{rl}& V(r)=\sum _{G^{\prime \prime }}V\left(G^{\prime \prime }\right)e^{i{G}^{\prime \prime }r}\\ & V(r){\Phi }_j(r)=\sum _{G,G^{\prime \prime }}V\left(G^{\prime \prime }\right)e^{i{G}^{\prime \prime }r}{\Phi }_j(k+G)e^{i(k+G)r}\\ & \int V(r){\Phi }_j(r)e^{-i(k+G^{\prime })r}dr=\\ & \sum _{G,G^{\prime \prime }}V\left(G^{\prime \prime }\right){\Phi }_j(k+G)\delta \left(G^{\prime \prime }+G-G^{\prime }\right)\Omega =\\ & =\sum _{G}V\left(G^{\prime }-G\right){\Phi }_j(k+G)\Omega \end{array}$$
外势场项
$$\begin{array}{rl}& \sum _{l,\mu ,G}\int \left(V_l\left(\left|\mathbf{r}-{\mathbf{R}}_{\mu }\right|\right){\mathcal{P}}_l{e}^{i(\mathbf{k}+\mathbf{G})\mathbf{r}}\right)e^{-i\left(\mathbf{k}+{\mathbf{G}}^{\prime }\right)\mathbf{r}}d\mathbf{r}{\Phi }_j(\mathbf{k}+\mathbf{G})\\ =& \Omega \sum _{l,\mathbf{G}}{\Phi }_j(\mathbf{k}+\mathbf{G})S\left(\mathbf{G}-{\mathbf{G}}^{\prime }\right)V_{l,\mathbf{k},\mathbf{G},{\mathbf{G}}^{\prime }}^{ps}\end{array}$$
对于半局域㕍势, 可分为两项
$$\begin{array}{rl}& \sum _{l,G^{\prime }}S\left(G^{\prime }-G\right)V_{lkG{G}^{\prime }}^{PS}{\Phi }_j\left(k+G^{\prime }\right)=\\ & =\sum _{G^{\prime }}S\left(G^{\prime }-G\right)V^{\text{core }}\left(G-G^{\prime }\right){\Phi }_j\left(k+G^{\prime }\right)+\sum _{l,G^{\prime }}S\left(G^{\prime }-G\right)V_{lkG{G}^{\prime }}^{NL}{\Phi }_j\left(k+G^{\prime }\right)\end{array}$$

动量空间的 Kohn-Sham 方程
$$\sum _{{\mathrm{G}}^{\prime }}\left[{\frac{1}{2}\left|\mathbf{k}+{\mathbf{G}}^{\prime }\right||}^2{\delta }_{{\mathrm{G}\mathrm{G}}^{\prime }}+V_{{\mathrm{k}\mathrm{G}\mathrm{G}}^{\prime }}^{eff}\right]{\Phi }_j\left({\mathrm{G}}^{\prime }\right)={ϵ}_j{\Phi }_j(\mathbf{G})$$
其中动量空间的有效势:
$$\begin{array}{rl}{V}_{\mathbf{k}\mathbf{G}\mathbf{G}}^{eff}=& V_{\text{coul }}\left(\mathbf{G}-{\mathbf{G}}^{\prime }\right)\\ & +{\hat{V}}_{xc}\left(\mathbf{G}-{\mathbf{G}}^{\prime }\right)+S\left({\mathbf{G}}^{\prime }-\mathbf{G}\right)(V^{\mathrm{core}}\left(\mathbf{G}-{\mathbf{G}}^{\prime }\right)\\ & +\sum _l{V}_{l,\mathbf{k},\mathbf{G},{\mathbf{G}}^{\prime }}^{\mathrm{N}\mathrm{L}}).\end{array}$$

截断能

平面波本身是一组正交完备的周期性函数组。对于周期性系统选择平面波基组是非常便利的。平面波本身是正交的, 所以交叠矩阵为单位矩阵。 势能矩阵就是势能的傅里叶变换。 整个计算过程非常简单。
但是平面波本身有无限多, 实际计算过程中只需要选取能量较低的平面波, 所以基于平面波程序中都 有平面波截断能量。由于平面波本身是完备的, 所以可以通过增加截断能量来增加计算精度。

平面波能量截断
$${\Phi }_{j\mathrm{k}}(\mathrm{r})=\sum _{\mathrm{G}}{C}_{j\mathrm{k}}(\mathrm{k}+\mathrm{G})e^{i(\mathrm{k}+\mathrm{G})\mathrm{r}},H{\Phi }_{j\mathrm{k}}=E_{j\mathrm{k}}{\Phi }_{j\mathrm{k}}$$

$$\frac{{\hslash }^2}{2m}|\vec{k}+\vec{G}{|}^2\le E_{\text{cutoff }}$$

倒空间波函数
当求解 $k$-空间Kohn-Sham方程时, 对给定 $k$, 我们得到一系列的能级和 本征函数, 使得我们可得到能带或 色散关系。然而, Kohn-Sham方程是 需要自洽求解的方程, 自洽的有效 势只与坐标有关。
$${\Phi }_{j\mathrm{k}}(\mathrm{r})=\sum _{\mathrm{G}}{C}_{j\mathrm{k}}(\mathrm{k}+\mathrm{G})e^{i(\mathrm{k}+\mathrm{G})\mathrm{r}},H{\Phi }_{j\mathrm{k}}=E_{j\mathrm{k}}{\Phi }_{j\mathrm{k}}$$

倒空间求和
看看如何计算电子密度,
$$\rho (\mathbf{r})=\sum _j{\Phi }_j^{\ast }(\mathbf{r}){\Phi }_j(\mathbf{r}),$$
实际上是需要对不同的 $k$ 求平均。
$$\rho (\mathbf{r})=\sum _{j\mathbf{k}}{\mathbf{w}}_j{\Phi }_{j\mathbf{k}}^{\prime }(\mathbf{r}){\Phi }_{j\mathbf{k}}(\mathbf{r})$$
对 $k$ 空间求平均的离散化会引起误差，不同的k点越多, 计算精度越高, 计算量越大。
$$n(\vec{r})=\frac{1}{N_k}\sum _{\vec{k},i}^{occ}{n}_{i,\vec{k}}(\vec{r})$$
需要在 $k$ 空间中积分:
$$\frac{1}{{\Omega }_{BZ}}{\int }_{BZ}d\vec{k}{f}_i(\vec{k})\stackrel{\text{ discretize }}{⟶}\frac{1}{N_k}\sum _{\vec{k}}{f}_i(\vec{k})$$

平面波 DFT 方法概述

我们由变分原理可推导出 Kohn-Sham 方程 这是一个非相互作用电子的Schrödinger方程
$$\left[-\frac{1}{2}{\nabla }_i^2+{\hat{V}}_{eff}(\mathbf{r})\right]{\varphi }_i^{KS}(\mathbf{r})={ϵ}_i{\varphi }_i^{KS}(\mathbf{r})$$
更具体地，
$$[-\frac{1}{2}{\nabla }_i^2+\underset{{\hat{V}}_{eff}(\mathbf{r})}{\underbrace{{\hat{v}}_{ext}(\mathbf{r})+{\hat{U}}_{cl}(\mathbf{r})+{\hat{V}}_{xc}(\mathbf{r})}}]{\varphi }_i^{KS}(\mathbf{r})={ϵ}_i{\varphi }_i^{KS}(\mathbf{r})$$
能量表达式
$$\begin{array}{c}E[\rho ]=\sum _{i=1}^N{ϵ}_i-\frac{1}{2}\iint \frac{\rho (\mathbf{r})\rho \left({\mathbf{r}}^{\prime }\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}d\mathbf{r}d{\mathbf{r}}^{\prime }\\ -\int {\hat{V}}_{xc}(\mathbf{r})\rho (\mathbf{r})d\mathbf{r}+E_{xc}[\rho ]\\ E=\sum _i^N{\epsilon }_i-E_H[\rho ]+E_{\mathrm{x}\mathrm{c}}[\rho ]-\int \frac{\delta E_{\mathrm{x}\mathrm{c}}[\rho ]}{\delta \rho (\mathbf{r})}\rho (\mathbf{r})d\mathbf{r}\end{array}$$
$${\hat{V}}_{xc}(\mathbf{r})=\frac{\delta E_{xc}[\rho (\mathbf{r})]}{\delta \rho (\mathbf{r})}$$

Kohn-Sham轨道 ${\varphi }_i^{KS}(\mathbf{r})$, 可用来计算总的电子密度:
$$\rho (\mathbf{r})=\sum _{i=1}^N{\left|{\varphi }_i^{KS}(\mathbf{r})\right|}^2.$$

电子密度 $\rho (\mathbf{r})$ 又可用于计算改进的有效 势 ${\hat{V}}_{eff}(\mathbf{r})$, 这样导致了一个自洽计算:
$$\begin{array}{rl}{\rho }_0(\mathbf{r}))⟶& {\hat{V}}_{eff}(\mathbf{r})\to {\varphi }_i^{KS}(\mathbf{r})\\ & \downarrow \\ & \rho (\mathbf{r})\end{array}$$
动量空间的 Kohn-Sham 方程
$$\sum _{{\mathrm{G}}^{\prime }}\left[\frac{1}{2}\left|\mathbf{k}+{\mathbf{G}}^{\prime }\right|{}^2{\delta }_{{\mathrm{G}\mathrm{G}}^{\prime }}+V_{\mathrm{k}G{G}^{\prime }}^{eff}\right]{\Phi }_j\left({\mathbf{G}}^{\prime }\right)={ϵ}_j{\Phi }_j(\mathbf{G}),$$
Ewald 求和
$$\begin{array}{rl}{E}_{nn}=& \frac{1}{2}{\sum }_{\mu ,v}\frac{4\pi }{\Omega }{\sum }_{\mathrm{G}\ne 0}\frac{Z_{\mu }{Z}_V}{|\mathbf{G}{|}^2}{e}^{-\frac{|\mathrm{G}{|}^2}{4x^2}}{e}^{i\mathrm{G}\cdot \left({\mathbf{R}}_{\mu }-{\mathbf{R}}_v\right)}\\ & +\frac{1}{2}\sum Z_{\mu }{Z}_V\frac{\mathrm{erfc}\left(|{\mathbf{R}}_{\mu }-{\mathbf{R}}_v|\kappa \right)}{|{\mathbf{R}}_{\mu }-{\mathbf{R}}_v|}\\ & -\frac{1}{2}{\sum }_{\mu \ne v}\frac{4\pi }{\Omega }{Z}_{\mu }{Z}_V\frac{1}{4{\kappa }^2}-\sum Z_{\mu }^2\frac{\kappa }{\sqrt{\pi }}\\ & +{\lim }_{\mathbf{G}\to 0}\frac{1}{2}\sum \frac{4\pi }{\Omega }{Z}_{\mu }{Z}_V\frac{1}{|\mathbf{G}{|}^2}\end{array}$$

$$\begin{array}{rl}{E}_{\text{tot }}=& \sum _j{ϵ}_j-\Omega \frac{1}{2}\sum _{\mathbf{G}\ne 0}{V}_{\text{coul }}(\mathbf{G}){\rho }^{\ast }(\mathbf{G})\\ & -\Omega \sum _{\mathbf{G}\ne 0}\left[V_{xc}(\mathbf{G})-{ϵ}_{xc}(\mathbf{G})\right]{\rho }^{\ast }(\mathbf{G})\\ & +\sum _v{\alpha }_v\sum _{\mu }{Z}_{\mu }+{\gamma }_{\text{Ewald }}.\end{array}$$
能量的截断，
$${\Phi }_{j\mathrm{k}}(\mathrm{r})=\sum _{\mathrm{G}}{C}_{j\mathrm{k}}(\mathrm{k}+\mathrm{G})e^{i(\mathrm{k}+\mathrm{G})\mathrm{r}},H{\Phi }_{j\mathrm{k}}=E_{j\mathrm{k}}{\Phi }_{j\mathrm{k}}$$

• $k$ 空间KS方程是无穷维矩阵问题, 必须截断
• 对给定波矢 $k$, 平面波的数目由动能截断决定
  $$\frac{{\hslash }^2}{2m}|\vec{k}+\vec{G}{|}^2\le E_{\text{cutoff }}$$
• 动能截断与赝势选取有关
• 动能截断导致基底的不完备性
• 增加动能截断有总能量的一致收敛性

Monkhorst-Pack 方法选取 k 点，
$${\vec{k}}_{n_1,n_2,n_3}=\sum _1^3\frac{2n_i-N_i-1}{2N_i}{\vec{G}}_i$$
$$\left[-\frac{1}{2}{\nabla }_i^2+{\hat{V}}_{eff}(\mathbf{r})\right]{\varphi }_i^{KS}(\mathbf{r})={ϵ}_i{\varphi }_i^{KS}(\mathbf{r})$$

求能量极小是一个计算数学问题，有很多方法。
$$\begin{array}{ll}\text{ 全局总能量极小化方法 }& \text{ 局部总能量极小化方法 }\\ \text{ 单形法 (simplex) }& \text{ 牛顿法(Newton-Raphson) }\\ \text{ 模拟退火法 (simulated }& \text{ 割线法 (Secant) }\\ \text{ annealing) }& \text{ 最陡下降法(steepest) }\\ \text{ 遗传算法(genetic algorithm) }& \text{ 共轭梯度法 (conjugate }\\ \text{ 蚁群算法 (ant colony algorithm) }& \text{ gradient) }\\ \text{ 图论算法 (Graph Cut) }& \text{ 。。。。。 }\\ \text{ 粒子群算法(Particle Swam }& \\ \text{ Optimization)。。。。。。 }& \end{array}$$

参考书

• Sholl D, Steckel J A. Density functional theory: a practical introduction[M]. John Wiley & Sons, 2011.（有中译本）

• Martin R M. Electronic structure: basic theory and practical methods[M]. Cambridge university press, 2020.