LinLin（二）续

2.6 Extension to finite temperature

---

2.6扩展到有限温度

到目前为止，我们一直关注零温度，从而关注多体系统的基态能量。对于有限的温度，量子系统不是最小化能量，而是最小化自由能。因此，具有更高能量的状态也将做出贡献，因此使用（多体）密度矩阵P（N）来描述系统变得更方便，其中在这里和下文中，我们使用上标（N）来强调N体量。密度矩阵是Hilbert空间上具有单位迹的非负自伴算子。
根据密度矩阵，自由能由
$$F[P^{(N)}]=\mathrm{Tr}{P}^{(N)}(H+{\beta }^{-1}\ln P^{(N)}),\text{\hspace{1em}(2.6.1)}$$
其中β是逆温度，${\beta }^{-1}Tr{P}^{(N)}ln{P}^{(N)}$是冯·诺依曼熵，而$Tr{P}^{(N)}H$是能量。注意，在零温度的极限下，β→∞，上述自由能减少到能量$Tr{P}^{(N)}H$。有限温度系统的状态由最小化自由能$F[P^{(N)}]$来确定，其约束条件是P（N）是自伴的、半正定的，并且它的迹是1。
由于P（N）是自伴的，我们可以把它的特征分解写成
$$P^{(N)}=\sum _i{f}_i|{\Psi }_i⟩⟨{\Psi }_i|,\text{\hspace{1em}(2.6.2)}$$
其中{fi}称为职业号码。由于P(N)是非负的并且Tr P(N)=1，我们有$$0\le f_i\le 1\text{and}\sum _i{f}_i=1.\text{\hspace{1em}(2.6.3)}$$
下面我们给出能量和熵：
从下面两个公式可以发现自由能等于能量和熵的和。
$$\text{\hspace{1em}(2.6.4)}$$
$$\text{\hspace{1em}(2.6.5)}$$
注意，对于P（N）是与波函数$|\Phi ⟩$相关联的投影算子的特殊情况，即$P^{(N)}=|\Phi ⟩⟨\Phi |$，我们有$Tr{P}^{(N)}H=<\Phi |H|\Phi >$和$Tr{P}^{(N)}ln{P}^{(N)}=0$，因此自由能减少到波函数$|\Phi ⟩$的能量。然而，熵项偏好$f_i$严格地介于0和1之间，以便使 $f_{i}ln{f}_i$为负，因此有限温度下的最小化密度矩阵通常不是投影仪（即，$f_i$为0或1）。
（2.6.1）极小化的欧拉——拉格朗日方程是
$$H+{\beta }^{-1}(\ln P^{(N)}+I)-\lambda I=0,\text{\hspace{1em}(2.6.6)}$$
这里λ是约束TRP（N）=1的拉格朗日乘子。因此，最小化器由下式给出
$$P_{\beta }^{(N)}=\frac{1}{Z}\exp (-\beta H),\text{\hspace{1em}(2.6.7)}$$
其中Z是归一化常数Z=Tr exp(-β H)，使得$Tr{P}^{(N)}\beta =1$；Z在物理学术语中被称为配分函数。

---

与零温度情况类似，有限温度DFT背后的思想是将自由能视为单体电子密度的泛函，单体电子密度由密度矩阵给出为
$$\rho (r)=\sum _{\sigma }N\int P^{(N)}((r,\sigma ),x_2,\dots ,x_N;(r,\sigma ),x_2,\dots ,x_N)\mathrm{d}{x}_2\cdots \mathrm{d}{x}_N.\text{\hspace{1em}(2.6.8)}$$
遵循约束最小化的过程，我们有
$$\begin{array}{rl}\text{F}& =\underset{P^{(N)}\in {\mathcal{D}}_N}{\inf }\mathrm{Tr}(P^{(N)}H)+{\beta }^{-1}\mathrm{Tr}(P^{(N)}\ln P^{(N)})\\ & =\underset{\rho \in {\mathcal{J}}_N}{\inf }\left\{\underset{P^{(N)}\in {\mathcal{D}}_N}{\inf }\mathrm{Tr}(P^{(N)}H)+{\beta }^{-1}\mathrm{Tr}(P^{(N)}\ln P^{(N)})\right\},\end{array}\text{\hspace{1em}(2.6.9)}$$
其中$D_N$表示N体密度矩阵的集合。注意，对于任何ρ∈J_N，至少存在一个P（N）给出密度，因为我们可以用$P^{(N)}=|\Psi ⟩⟨\Psi |$和$\Psi \in {\mathcal{A}}_N$给出密度。因此，约束极小化是很好的定义，我们可以进一步得出
$$F=\underset{\rho \in {\mathcal{J}}_N}{\inf }\left\{{\mathcal{F}}_{\beta }[\rho ]+\int \rho V_{\mathrm{ext}}\mathrm{d}r\right\},\text{\hspace{1em}(2.6.10)}$$
其中${\mathcal{F}}_{\beta }[\rho ]$是仅取决于动能、电子-电子排斥和熵（因此取决于温度）的泛函：
$${\mathcal{F}}_{\beta }[\rho ]=\underset{\begin{array}{c}{P}^{(N)}\in {\mathcal{D}}_N\\ P^{(N)}\mapsto \rho \end{array}}{\inf }\left(\mathrm{Tr}{P}^{(N)}(T+V_{\mathrm{ee}})+{\beta }^{-1}\mathrm{Tr}(P^{(N)}\ln P^{(N)})\right).\text{\hspace{1em}(2.6.11)}$$

2.8 周期系统的布里渊区采样

Kohn-Sham DFT也可以应用于具有周期电位的晶体系统。在本节中，我们讨论与布里渊区采样（称为kpoint采样）相关的问题。类似于到目前为止使用的符号，我们使用r∈R3作为实空间中的变量，$k\in {\mathbb{R}}^3$作为倒数空间中的变量。为了简单起见，我们在本节中只考虑无自旋粒子。根据1.4节的讨论，我们用$\mathbb{L}$表示Bravais格，用${\mathbb{L}}^{\ast }$表示它在互易空间中的对偶格。设$\Omega$，${\Omega }^{\ast }$分别为Bravais晶格和第一布里渊区中的晶胞。$\Omega ,\Omega \ast$的体积由$|\Omega |，|{\Omega }^{\ast }|$表示，它们根据
$$|{\Omega }^{\ast }|=\frac{(2\pi {)}^3}{|\Omega |}.$$

Supercell formulation

---

超晶胞模型构建
尽管周期系统的大小是无限的，但事实上，认为周期系统只不过是一个“巨型分子”是非常有帮助的，至少在启发式上是这样。换句话说，我们可以用一个有限大小的系统来近似周期系统，该系统分别沿着Bravais晶格向量$a_1、a_2、a_3$具有$N_1^l、N_2^l、N_3^l$个晶胞。这个虚构的系统通常被称为用${\Omega }^l$表示的超胞。因此，超晶胞中晶胞的总数为
$$N^{\ell }=N_1^{\ell }\times N_2^{\ell }\times N_3^{\ell }.$$
假设每个晶胞的电子数为N，系统中的电子总数为$N{N}^l$。那么Kohn-Sham DFT可以用满足正交性条件的$N{N}^l$单粒子轨道 { ψ i } i = 1 N N l \{ψ_i\}^{NN^l}_{i=1} {ψi​}i=1NNl​来表示。每个轨道应满足超胞上的周期性边界条件，即，
$${\psi }_i(\mathbf{r}+{\mathbf{a}}_{\alpha }{N}_{\alpha }^{\ell })={\psi }_i(\mathbf{r}),\forall \mathbf{r}\in {\Omega }^{\ell },\alpha =1,2,3.$$
这种特殊的周期性边界条件称为玻恩-冯-卡门边界条件。相应地，电子密度$\rho (r)$也可以周期性地扩展到$R^3$。因此，超胞上的Kohn-Sham能量泛函可以写成如前所述： E ℓ [ { ψ i } ] = 1 2 ∫ Ω ℓ ∑ i ∣ ∇ ψ i ( r ) ∣ 2 d r + ∫ Ω ℓ V e x t ( r ) ρ ( r ) d r (2.8.1) E^\ell[\{\psi_i\}]=\frac12\int_{\Omega^\ell}\sum_i\lvert\nabla\psi_i(r)\rvert^2\mathrm{d}r+\int_{\Omega^\ell}V_{\mathrm{ext}}(r)\rho(r)\mathrm{d}r\tag{2.8.1} Eℓ[{ψi​}]=21​∫Ωℓ​i∑​∣∇ψi​(r)∣2dr+∫Ωℓ​Vext​(r)ρ(r)dr(2.8.1) + 1 2 ∫ Ω ℓ ∫ R 3 ρ ( r ) ρ ( r ′ ) ∣ r − r ′ ∣ d r ′ d r + E x c ℓ [ ρ ] + E I I ℓ [ { R I } ] . +\frac12\int_{\Omega^\ell}\int_{\mathbb{R}^3}\frac{\rho(\boldsymbol{r})\rho(\boldsymbol{r}^{\prime})}{|\boldsymbol{r}-\boldsymbol{r}^{\prime}|}\mathrm{d}\boldsymbol{r}^{\prime}\mathrm{d}\boldsymbol{r}+E_{\mathrm{xc}}^\ell[\rho]+E_{\mathrm{II}}^\ell[\{\boldsymbol{R}_I\}]. +21​∫Ωℓ​∫R3​∣r−r′∣ρ(r)ρ(r′)​dr′dr+Excℓ​[ρ]+EIIℓ​[{RI​}].

这里，交换关联能泛函和原子核排斥能也带有上标$l$，以反映它们是相对于超胞定义的事实。超胞中的基态能量可以变分地获得。
最小化问题（2.8.2）被很好地定义，我们可以预期，当我们接近热力学极限（${\Omega }^l$→$R^3$）时，每晶胞的总能量（定义为$E^l$/$N^l$）将对周期系统有一个很好地定义的极限。然而，从这个公式中并不清楚如何简洁地写出热力学极限。
      min ⁡ { ψ i } i = 1 N N ℓ E ℓ [ { ψ i } ] s . t . ∫ Ω ℓ ψ i ∗ ( r ) ψ j ( r ) d r = δ i , j ,   ρ ( r ) = ∑ i ∣ ψ i ( r ) ∣ 2 (2.8.2) \begin{aligned}\ \ \ \ \ &\min_{{\{\psi_{i}\}_{i=1}^{NN}\ell}}\quad E^{\ell}[\{\psi_{i}\}]\\&\mathrm{s.t.}\int_{\Omega^\ell}\psi_i^*(\boldsymbol{r})\psi_j(\boldsymbol{r})\mathrm{d}\boldsymbol{r}=\delta_{i,j},\\&\ \rho(\boldsymbol{r})=\sum_i\lvert\psi_i(\boldsymbol{r})\rvert^2\end{aligned}\tag{2.8.2}      ​{ψi​}i=1NN​ℓmin​Eℓ[{ψi​}]s.t.∫Ωℓ​ψi∗​(r)ψj​(r)dr=δi,j​, ρ(r)=i∑​∣ψi​(r)∣2​(2.8.2)
此外，该公式具有两个主要的计算缺点:
首先，Kohn-Sham轨道的数量是$N{N}^l$。基于对角化的求解器的计算复杂度按O(($N{N}^l$)${}^3$)缩放。
第二，积分的定义域是${\Omega }^l$，随着$N^l$变大，这也变得非常昂贵。

Unit cell formulation

---

晶胞配方

幸运的是，所有这些问题都可以使用Bloch-Floquet分解来解决。由于外部势相对于晶胞仍然是周期性的，1.4节中的Bloch分解将单轨道指数i分解为两个指数(n，k），如下
$${\psi }_i(\mathbf{r}):={\psi }_{n,\mathbf{k}}(\mathbf{r})=\frac{1}{\sqrt{N^{\ell }}}{e}^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}}{u}_{n,\mathbf{k}}(\mathbf{r}),k\in {\Omega }^{\ast }.\text{\hspace{1em}(2.8.3)}$$
这里周期部分$u_{n，k}(r)$满足
$$u_{n,\mathbf{k}}(r+R)=u_{n,\mathbf{k}}(r),\forall R\in \mathbb{L}\text{\hspace{1em}(2.8.4)}$$
它是晶胞Ω中的周期性边界条件。为了便于后面的讨论，引入了归一化因子$\frac{1}{\sqrt{N^{\ell }}}$。现在我们假设每个k点具有相同数量的轨道，即n=1，…，N.如稍后将看到的，这对应于零温度下的绝缘系统的设置。

对于固定晶胞，Born-von Karman边界条件对相位因子施加了额外的条件，如
$$e^{\mathrm{i}\mathbf{k}\cdot (\mathbf{r}+{\mathbf{a}}_{\alpha }{N}_{\alpha }^{\ell })}=e^{\mathrm{i}\mathbf{k}\cdot \mathbf{r}},\forall \mathbf{r}\in {\Omega }^{\ell },\alpha =1,2,3,$$
或等同地，
$$e^{\mathrm{i}{N}_{\alpha }^{\ell }\mathbf{k}\cdot {\mathbf{a}}_{\alpha }}=1,\alpha =1,2,3.$$
因此，$\mathbf{k}$不可能是${\Omega }^{\ast }$中的任意点。在不失一般性的情况下，我们假设$N_1^l、N_2^l、N_3^l$都是偶数，并定义集合
$${\mathcal{K}}^{\ell }=\left\{\sum _{\alpha =1}^3\frac{m_{\alpha }}{N_{\alpha }^{\ell }}{\mathbf{b}}_{\alpha }|m_{\alpha }=-\frac{N_{\alpha }^{\ell }}{2}+1,\dots ,\frac{N_{\alpha }^{\ell }}{2},\alpha =1,2,3\right\}\subset {\Omega }^{\ast }.$$
则玻恩-冯卡门边界条件要求$k\in K^l$。注意，基数$|{\mathcal{K}}^l|=N^l$，因此允许的k个点的总数等于超胞${\Omega }^l$中的晶胞数量。
由于$K^l$是一个离散集，从${\psi }_{n,\mathbf{k}}$的正交性条件，我们有
$$\begin{array}{rl}{\delta }_{n^{\prime },n}{\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}}& ={\int }_{{\Omega }^{\ell }}{\psi }_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }(r){\psi }_{n,\mathbf{k}}(r)\mathrm{d}r\\ & =\frac{1}{N^{\ell }}{\int }_{{\Omega }^{\ell }}{e}^{\mathrm{i}(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{r}}{u}_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }(r)u_{n,\mathbf{k}}(r)\mathrm{d}r\\ & =\frac{1}{N^{\ell }}\sum _{\mathbf{R}\in \mathcal{L}}{\int }_{\Omega }{e}^{\mathrm{i}(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot (\mathbf{r}+\mathbf{R})}{u}_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }(\mathbf{r}+\mathbf{R})u_{n,\mathbf{k}}(\mathbf{r}+\mathbf{R})\mathrm{d}\mathbf{r}\\ & =\left(\frac{1}{N^{\ell }}\sum _{\mathbf{R}\in \mathbb{L}}{e}^{\mathrm{i}(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{R}}\right){\int }_{\Omega }{e}^{\mathrm{i}(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot \mathbf{r}}{u}_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }(\mathbf{r})u_{n,\mathbf{k}}(\mathbf{r})\mathrm{d}\mathbf{r}\\ & ={\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}}{\int }_{\Omega }{u}_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }(\mathbf{r})u_{n,\mathbf{k}}(\mathbf{r})\mathrm{d}\mathbf{r}.\end{array}\text{\hspace{1em}(2.8.5)}$$
在这里，我们使用了$u_{n，k}$的周期性以及$k-k’$的差异在互易空间中的任何方向上不能相差一个以上的晶格向量的事实。因此，超胞中ψ n，k的正交性条件可以等价地写成晶胞中$u_{n，k}$的正交性条件，即，
$$\text{\hspace{1em}(2.8.6)}$$
我们强调，当$k=k’$时，$u_{n，k}$ 和 $u_{n’，k’}$通常彼此不正交。
现在让我们用 { u n ， k } \{u_{n，k}\} {un，k​}重写（2.8.1）中每晶胞的能量。首先，每个晶胞的动能部分是（注意等式两边的归一化因子）

$$\frac{1}{N^{\ell }}{\int }_{{\Omega }^{\ell }}\sum _{n,\mathbf{k}}\frac{1}{2}|\nabla {\psi }_{n,\mathbf{k}}(\mathbf{r}){|}^2\mathrm{d}\mathbf{r}=\frac{1}{N^{\ell }}\sum _{\mathbf{k}}{\int }_{\Omega }\sum _n\frac{1}{2}|({\nabla }_r+\mathrm{i}\mathbf{k})u_{n,\mathbf{k}}(\mathbf{r}){|}^2\mathrm{d}\mathbf{r}.$$

由于集合KaTeX parse error: Got function '\boldsymbol' with no arguments as superscript at position 3: K^\̲b̲o̲l̲d̲s̲y̲m̲b̲o̲l̲{\ell}是离散布里渊区的均匀网格点的集合，$\mathbf{k}$上的和可以看作${\Omega }^{\ast }$中积分的黎曼和，即存在极限

带入上述从而将动能部分写为：
$$\frac{1}{|{\Omega }^{\ast }|}{\int }_{{\Omega }^{\ast }}{\int }_{\Omega }\sum _n\frac{1}{2}|({\nabla }_{\mathbf{r}}+\mathrm{i}\mathbf{k})u_{n,\mathbf{k}}(\mathbf{r}){|}^2\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{k}.\text{\hspace{1em}(2.8.8)}$$

（热力学下的）电子密度：

$$\rho (\mathbf{r})=\frac{1}{|{\Omega }^{\ast }|}{\int }_{{\Omega }^{\ast }}\sum _n|u_{n,\mathbf{k}}(r){|}^2\mathrm{d}k,\text{\hspace{1em}(2.8.9)}$$

它是晶胞Ω中的周期函数。由于能量泛函中的所有其他项仅取决于电子密度，因此每个晶胞的能量泛函可以使 { u n ， k } \{u_{n，k}\} {un，k​}在热力学极限下简洁地写成如下：
E [ { u n , k } ] = 1 ∣ Ω ∗ ∣ ∫ Ω ∗ ∫ Ω ∑ n = 1 N 1 2 ∣ ( ∇ r + i k ) u n , k ( r ) ∣ 2 d r d k + ∫ Ω V e x t ( r ) ρ ( r ) d r + 1 2 ∫ Ω ∫ R 3 ρ ( r ) ρ ( r ′ ) ∣ r − r ′ ∣ d r ′ d r + E x c [ ρ ] + E I I [ { R I } ] . \begin{gathered} \mathcal{E}[\{u_{n,\mathbf{k}}\}]= \frac1{|\Omega^*|}\int_{\Omega^*}\int_{\Omega}\sum_{n=1}^{N}\frac12|(\nabla_{\boldsymbol{r}}+\mathrm{i}\boldsymbol{k})u_{n,\boldsymbol{k}}(\boldsymbol{r})|^2 \mathrm{d}\boldsymbol{r} \mathrm{d}\boldsymbol{k}+\int_{\Omega}V_{\mathrm{ext}}(\boldsymbol{r})\rho(\boldsymbol{r}) \mathrm{d}\boldsymbol{r} \\ +\frac12\int_\Omega\int_{\mathbb{R}^3}\frac{\rho(\boldsymbol{r})\rho(\boldsymbol{r}^{\prime})}{|r-r^{\prime}|} \mathrm{d}\boldsymbol{r}^{\prime} \mathrm{d}r+E_{\mathrm{xc}}[\rho]+E_{\mathrm{II}}[\{\boldsymbol{R}_I\}]. \end{gathered} E[{un,k​}]=∣Ω∗∣1​∫Ω∗​∫Ω​n=1∑N​21​∣(∇r​+ik)un,k​(r)∣2drdk+∫Ω​Vext​(r)ρ(r)dr+21​∫Ω​∫R3​∣r−r′∣ρ(r)ρ(r′)​dr′dr+Exc​[ρ]+EII​[{RI​}].​
注意，实空间中的所有积分范围都限于晶胞Ω，除了长程库仑相互作用，其中r’变量包括晶胞Ω内的电子密度对静电相互作用的贡献以及其在R3中的所有周期图像。注意，ρ是周期函数，并且在R3上的积分导致发散的Hartree势。然而，这种发散将被Vext中的电子-原子核相互作用和$E_{II}[R_I]$中的原子核-原子核排斥所抵消。来自原子核的总静电能量和电子密度是明确定义的。这至少对于电荷中性系统是正确的[50]。交换关联泛函Exc[ρ]和原子核排斥能EII也应理解为每晶胞的能量。那么，热力学极限下周期系统的Kohn-Sham DFT可以变分地表述为
min ⁡ { u n , k } E [ { u n , k } ] s . t . ∫ Ω u n ′ , k ∗ ( r ) u n , k ( r ) d r = δ n ′ , n , ρ ( r ) = 1 ∣ Ω ∗ ∣ ∫ Ω ∗ ∑ n ∣ u n , k ( r ) ∣ 2 d k . (2.8.11) \begin{aligned} \operatorname*{min}_{\{u_{n,k}\}}& \mathcal{E}[\{u_{n,\mathbf{k}}\}] \\ \mathbf{s.t.}& \int_\Omega u_{n^{\prime},\boldsymbol{k}}^*(r)u_{n,\boldsymbol{k}}(r) \mathrm{d}r=\delta_{n^{\prime},n}, \\ &\rho(\boldsymbol{r})=\frac{1}{|\Omega^{*}|}\int_{\Omega^{*}}\sum_{n}\lvert u_{n,\boldsymbol{k}}(\boldsymbol{r})\rvert^{2} \mathrm{d}\boldsymbol{k}. \end{aligned}\tag{2.8.11} {un,k​}min​s.t.​E[{un,k​}]∫Ω​un′,k∗​(r)un,k​(r)dr=δn′,n​,ρ(r)=∣Ω∗∣1​∫Ω∗​n∑​∣un,k​(r)∣2dk.​(2.8.11)
上式的欧拉-拉格朗日方程为：
$$\begin{array}{rl}(H_{\mathbf{k}}{u}_{n,\mathbf{k}})(\mathbf{r})& :=\left(-\frac{1}{2}{\left(\nabla +\mathrm{i}\mathbf{k}\right)}^2+V_{\mathrm{eff}}[\rho ](\mathbf{r})\right)u_{n,\mathbf{k}}(\mathbf{r})\\ & ={\epsilon }_{n,\mathbf{k}}{u}_{n,\mathbf{k}}(\mathbf{r}),\mathbf{r}\in \Omega ,\mathbf{k}\in {\Omega }^{\ast },\end{array}\text{\hspace{1em}(2.8.12)}$$
服从电子密度的正交性条件和自洽条件。有效电位定义为
$$V_{\mathrm{eff}}[\rho ](r)=V_{\mathrm{ext}}(r)+{\int }_{{\mathbb{R}}^3}\frac{\rho (r^{\prime })}{|r-r^{\prime }|}\mathrm{d}{r}^{\prime }+V_{\mathrm{xc}}[\rho ](r).$$
特别地，如果Vext和ρ是关于Bravais晶格的周期函数，那么Veff[ρ]（r）也是周期函数。这证明了前面讨论中的布洛赫分解。

与超胞公式相比，（2.8.11）和（2.8.12）中给出的公式具有主要的计算优势。在布里渊区离散成$N^l$点后，每个k点对应于仅在晶胞上定义的解耦特征值问题。因此，关于$N^l$的计算成本从立方缩放减少到线性缩放。这使得即使使用适度的计算资源也能够使用大量的k点进行模拟。式$(\mathrm{2.8.12})$中的所有特征值{ε_{n，k}}的集合形成由Kohn-Sham DFT建模的能带结构。

---

现在我们回到Bloch分解之前的轨道{ ψn，k }。虽然每个ψ n，k都可以在任意有限大小的超胞内进行归一化，但这样的归一化ψ n，k在热力学极限下会弱收敛到零。这是因为ψ n，k不属于L2 ( R3 )，而只是H在R3中的一个广义特征函数.这类似于1.2节中平面波可以看成是动量算符的广义本征函数。为此，我们可以利用Bloch分解来定义广义特征函数，仍然记为{ ψn，k }为(对所有的r∈R3)
$$\text{\hspace{1em}(2.8.15)}$$
与( 2.8.3 )相比，我们降低了归一化因子。由此得到的{ψn，k}在分布意义下满足正交条件：
$${\int }_{{\mathbb{R}}^3}{\psi }_{n^{\prime },k^{\prime }}^{\ast }(r){\psi }_{n,k}(r)\mathrm{d}r=|{\Omega }^{\ast }|{\delta }_{n^{\prime },n}\delta (k^{\prime }-k).\text{\hspace{1em}(2.8.15)}$$
这暗示了在布里渊区积分时的归一化条件：
$$\frac{1}{|{\Omega }^{\ast }|}{\int }_{{\Omega }^{\ast }}{\int }_{{\mathbb{R}}^3}{\psi }_{n^{\prime },{\mathbf{k}}^{\prime }}^{\ast }(\mathbf{r}){\psi }_{n,\mathbf{k}}(\mathbf{r})\mathrm{d}\mathbf{r}\mathrm{d}\mathbf{k}={\delta }_{n^{\prime },\mathbf{n}}.\text{\hspace{1em}(2.8.16)}$$
选择( 2.8.14 )式的另一个好处是电子密度也可以用{ψ_{n，k}}来写：

Kohn-Sham DFT要求能量较低的轨道必须先于能量较高的轨道占据。
到目前为止，我们假设占据轨道的数目N对所有的k一致成立。只有在满足下面的带外隔离度条件时，情况才是如此：
$${\epsilon }_{N,\mathbf{k}}<{\epsilon }_{N+1,{\mathbf{k}}^{\prime }}\forall \mathbf{k},{\mathbf{k}}^{\prime }\in {\Omega }^{\ast }.\text{\hspace{1em}(2.8.17)}$$
满足隔离条件的系统称为绝缘系统，其具有的正带隙定义为
$$E_{\mathbf{g}}=\left(\underset{\mathbf{k}\in {\Omega }^{\ast }}{\inf }{\epsilon }_{N+1,\mathbf{k}}\right)-\left(\underset{\mathbf{k}\in {\Omega }^{\ast }}{\sup }{\epsilon }_{N,\mathbf{k}}\right).\text{\hspace{1em}(2.8.18)}$$
另一方面，当系统的隔离度被破坏时，称该系统为金属系统，并定义带隙Eg为0。在材料科学中，另一个常用术语是半导体体系，它是一种绝缘体系，但具有相对较小的带隙(通常< 1 eV)。因此，这个术语只在数量意义上成立。

对于金属体系，能量较低的轨道必须首先被占据，并且每k个点所占据的轨道数目在整个布里渊区是不均匀的。绝缘和金属系统的处理可以通过引入化学势$\mu$来统一，$\mu$定义了一组占位数
$$f_{n,\mathbf{k}}:=f_{\infty }({\epsilon }_{n,\mathbf{k}}-\mu )=\{\begin{array}{ll}1,& {\epsilon }_{\mathbf{n},\mathbf{k}}\le \mu ,\\ 0,& {\epsilon }_{\mathbf{n},\mathbf{k}}>\mu .\end{array}$$
然后，在求解Kohn - Sham方程( 2.8.12 )后，电子密度被定义为占据数
$$\rho (\mathbf{r})=\frac{1}{|{\Omega }^{\ast }|}{\int }_{{\Omega }^{\ast }}\sum _n{f}_{n,\mathbf{k}}|u_{n,\mathbf{k}}(\mathbf{r}){|}^2\mathrm{d}\mathbf{k}.\text{\hspace{1em}(2.8.19)}$$
化学势μ应自洽地调节，以满足每个晶胞有N个电子的条件，即：
$${\int }_{\Omega }\rho (r)\mathrm{d}r=N.\text{\hspace{1em}(2.8.20)}$$
因此，化学势μ也是与N相关的拉格朗日乘子.从这个角度出发，在给定带外隔离度的条件下，绝缘系统占位数的选择就变成了
$$f_{n,\mathbf{k}}=\{\begin{array}{ll}1,& 1\le n\le N,\\ 0,& n>N.\end{array}$$
最后，通过引入占位数$f_{n,k}$的熵项，周期系统的能量泛函也可以推广到有限温度情形.我们将在此省略细节。

2.10 几何优化和从头算分子动力学

在玻恩——奥本海默近似下，原子核由经典力学描述，电子根据与每个原子构型相关的基态波函数从属于原子核。总能量E({RI})提供了原子间势能的第一原理描述。
几何优化(或结构优化)问题旨在寻找势能面E({RI})上的(全局或局部)极小值，相应的极小值{R∫I}是(全局或局部)稳定构型。大多数自然存在的分子和固体都处于或接近这样的极小值。
定义第i个原子核的原子力为
F I ( { R I } ) = − ∂ E ( { R I } ) ∂ R I ; F_I(\{R_I\})=-\frac{\partial E(\{R_I\})}{\partial R_I}; FI​({RI​})=−∂RI​∂E({RI​})​;

那么，确定极小值的一个必要条件是，对于所有原子，原子力$F_I$都应该为零。
对于一个有M个原子核的封闭系统，原子核自由度的总能量为
E t o t = ∑ I = 1 M M I R ˙ I 2 2 + E ( { R I } I = 1 M ) , E_{\mathrm{tot}}=\sum_{I=1}^M\frac{M_I\dot{\boldsymbol{R}}_I^2}2+E(\{\boldsymbol{R}_I\}_{I=1}^M), Etot​=I=1∑M​2MI​R˙I2​​+E({RI​}I=1M​),
其中MI是第i个原子核的质量。
原子核的动力学性质可以用牛顿方程来研究
M I R ¨ I ( t ) = − ∂ E ∂ R I ( { R I ( t ) } ) = F I ( { R I ( t ) } ) . M_I\ddot{\boldsymbol{R}}_I(t)=-\frac{\partial E}{\partial\boldsymbol{R}_I}(\{\boldsymbol{R}_I(t)\})=\boldsymbol{F}_I(\{\boldsymbol{R}_I(t)\}). MI​R¨I​(t)=−∂RI​∂E​({RI​(t)})=FI​({RI​(t)}).

这被称为从头算分子动力学(AIMD)。几何优化和分子动力学消除了构建经验原子间势的需要，并使化学和材料科学中静态和动态性质研究的广泛应用成为可能。

在几何优化和AIMD中，关键因素是评估原子力。在三维空间中，这种导数信息正式需要3M独立电子结构计算，这将是极其昂贵的。幸运的是，赫尔曼——费曼定理允许我们评估原子力，与单个电子结构计算相比，几乎没有增加成本。

原子力的评估

为了引入赫尔曼——费曼定理，让我们首先考虑一个参数相关的哈密顿量H(λ)。对于每个给定的λ，我们假设基态是非简并的。它的特征值ε(λ)和本征函数$\psi (r;\lambda )$可以计算为
$$H(\lambda )\psi (\lambda )=\epsilon (\lambda )\psi (\lambda )\text{\hspace{1em}(2.10.4)}$$
特征值对λ的导数是
$$\begin{array}{rl}\frac{\mathrm{d}\epsilon (\lambda )}{\mathrm{d}\lambda }& =\frac{\mathrm{d}}{\mathrm{d}\lambda }⟨\psi (\lambda )|H(\lambda )|\psi (\lambda )⟩\\ & =⟨\psi (\lambda )|\frac{\mathrm{d}H(\lambda )}{\mathrm{d}\lambda }|\psi (\lambda )⟩+\epsilon (\lambda )\frac{\mathrm{d}}{\mathrm{d}\lambda }⟨\psi (\lambda )|\psi (\lambda )⟩\\ & =⟨\psi (\lambda )|\frac{\mathrm{d}H(\lambda )}{\mathrm{d}\lambda }|\psi (\lambda )⟩.\end{array}$$
这里我们对所有λ使用了归一化条件<φ(λ)|φ(λ)>=1。因此，赫尔曼——费曼定理表明，特征值相对于外部参数扰动的一阶导数简单地由哈密顿量相对于λ的变化给出。
赫尔曼——费曼定理的另一个观点可以从变分原理中理解
$$\epsilon (\lambda )=\underset{|\psi )\ne |0⟩}{\inf }E[\psi ,\lambda ]:=\underset{|\psi ⟩\ne |0⟩}{\inf }\frac{⟨\psi |H(\lambda )|\psi ⟩}{⟨\psi |\psi ⟩}.$$
因此，
$$\frac{\mathrm{d}\epsilon (\lambda )}{\mathrm{d}\lambda }={\frac{\partial E[\psi ,\lambda ]}{\partial \lambda }|}_{\psi =\psi (\lambda )}+\int {\frac{\delta E[\psi ,\lambda ]}{\delta \psi (r)}|}_{\psi =\psi (\lambda )}\frac{\partial \psi (r;\lambda )}{\partial \lambda }\mathrm{d}r=⟨\psi (\lambda )|\frac{\mathrm{d}H(\lambda )}{\mathrm{d}\lambda }|\psi (\lambda )⟩.$$
请注意，虽然导数$\frac{\partial \psi (r;\lambda )}{\partial \lambda }$一般不会消失，但由于欧拉——拉格朗日方程，它对$\frac{\mathrm{d}\epsilon (\lambda )}{\mathrm{d}\lambda }$没有贡献。

现在，对于Kohn-Sham DFT中的总能量，只有外部势Vext和离子——离子相互作用EII明确依赖于原子构型{RI}。此外，外部势一般可以分解为每个原子的贡献，如
V e x t ( r ; { R I } ) = ∑ I V e x t I ( r − R I ) . V_\mathrm{ext}(r;\{R_I\})=\sum_IV_\mathrm{ext}^I(\boldsymbol{r}-\boldsymbol{R}_I). Vext​(r;{RI​})=I∑​VextI​(r−RI​).
然后，应用赫尔曼——费曼定理，我们发现
F I = − ∂ E ( { R I } ) ∂ R I = − ∫ ρ ( r ) ∂ V e x t I ( r − R I ) ∂ R I   d r − ∂ E I I ( { R I } ) ∂ R I . \boldsymbol{F}_I=-\frac{\partial E(\{\boldsymbol{R}_I\})}{\partial\boldsymbol{R}_I}=-\int\rho(\boldsymbol{r})\frac{\partial V_\mathrm{ext}^I(\boldsymbol{r}-\boldsymbol{R}_I)}{\partial\boldsymbol{R}_I}\mathrm{~d}\boldsymbol{r}-\frac{\partial E_\mathrm{II}(\{\boldsymbol{R}_I\})}{\partial\boldsymbol{R}_I}. FI​=−∂RI​∂E({RI​})​=−∫ρ(r)∂RI​∂VextI​(r−RI​)​ dr−∂RI​∂EII​({RI​})​.

Car-Parrinello分子动力学

在AIMD中，电子结构问题必须完全自洽地解决，以便达到每个原子构型的基态。因此，这种类型的AIMD也被称为玻恩——奥本海默分子动力学(BOMD)。很长一段时间，尽管BOMD有明显的潜力，但它被认为过于昂贵。AIMD是由开创性的Car-Parrinello分子动力学(CPMD)[14]实现的，它引入了一个扩展的拉格朗日函数，包括原子核和电子的自由度：

2.11 Time-dependent density functional theory

---

2.11含时密度泛函理论

密度泛函理论指出，多体基态能量和多体基态波函数仅由电子密度ρ决定。

关键的工具是变分原理，它产生了约束极小化公式。

对于含时多体Schr ö dinger方程，Runge和Gross于1984年提出了含时密度泛函理论( TDDFT )。然而，TDDFT的严格基础与DFT相比要清晰得多，因此我们并没有从物理角度给出原创性的论证，而只是给出了结论。多体薛定谔方程为
$$\mathrm{i}{\partial }_t\Psi (t)=H(t)\Psi (t),\text{\hspace{1em}(2.11.1)}$$
其中Ψ ( t ) = Ψ ( x1 , … , xN , t)是带有初始条件Ψ ( t0 ) = Ψ 0的含时多体波函数。随时间变化的电子密度是
$$\rho (\mathbf{r},t)=N\sum _{\sigma }\int {\left|\Psi ((\mathbf{r},\sigma ),{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N,t)\right|}^2\mathrm{d}{\mathbf{x}}_2\dots \mathrm{d}{\mathbf{x}}_N.$$
与时间相关的哈密顿量：
$$H(t)=\sum _{i=1}^N(-\frac{1}{2}{\Delta }_{{\mathbf{r}}_i}+V_{\mathrm{ext}}({\mathbf{r}}_i,t))+\sum _{i<j}^N\frac{1}{|{\mathbf{r}}_i-{\mathbf{r}}_j|},$$

其中显式的时间依赖性来自外部势Vext。Runge和Gross的陈述是任意时刻t≥t0的多体波函数Ψ（t）由初始态Ψ 0和电子密度${\rho (r，s)}_{t_0\le s\le t}$的历史唯一确定。
此外，根据Hohenberg-Kohn定理，如果系统从多体基态开始，那么ψ 0也由ρ（r，t0）唯一确定。因此，多体系统的演化完全由密度${\rho (r，s)}_{t_0\le s\le t}$的演化决定。此外，类似于Kohn-Sham DFT的构造，Runge-Gross TDDFT也可以由R中的以下耦合含时薛定谔方程正式给出
i ∂ t ψ j ( x , t ) = ( − 1 2 Δ + V e x t ( r , t ) + ∫ ρ ( r ′ , t ) ∣ r − r ′ ∣ d r ′ + V x c [ { ρ } t 0 ≤ s ≤ t ] ( r ) ) ψ j ( x , t ) , ρ ( r , t ) = ∑ σ ∑ j = 1 N ∣ ψ j ( x , t ) ∣ 2 . (2.11.4) \begin{aligned} \mathrm{i}\partial_{t}\psi_{j}(\boldsymbol{x},t)& =\left(-\frac12\Delta+V_{\mathrm{ext}}(r,t)+\int\frac{\rho(r^{\prime},t)}{|r-r^{\prime}|}\mathrm{d}r^{\prime}+V_{\mathrm{xc}}\left[\{\rho\}_{t_0\leq s\leq t}\right](r)\right)\psi_j(x,t), \\ \rho(\boldsymbol{r},t)& =\sum_{\sigma}\sum_{j=1}^{N}\left|\psi_{j}(x,t)\right|^{2}. \end{aligned}\tag{2.11.4} i∂t​ψj​(x,t)ρ(r,t)​=(−21​Δ+Vext​(r,t)+∫∣r−r′∣ρ(r′,t)​dr′+Vxc​[{ρ}t0​≤s≤t​](r))ψj​(x,t),=σ∑​j=1∑N​∣ψj​(x,t)∣2.​(2.11.4)
与基态Kohn-Sham DFT相比，重要的区别在于交换相关势Vxc原则上不仅在空间变量r中是非局域的，而且在时间变量t中也是非局域的。
记忆效应可以跨越整个历史${\rho (r,s)}_{t_0\le s\le t}$。鉴于在基态Kohn-Sham DFT中找到精确的交换相关泛函已经很困难，进一步识别TDDFT中的记忆效应变得更加困难。在实际的TDDFT计算中，几乎所有的交换相关泛函都假设
V x c [ { ρ } t 0 ≤ s ≤ t ] ( r ) ≈ V x c K S [ ρ ( t ) ] ( r ) , (2.11.5) V_{\mathrm{xc}}\left[\{\rho\}_{t_0\leq s\leq t}\right](r)\approx V_{\mathrm{xc}}^{\mathrm{KS}}\left[\rho(t)\right](r),\tag{2.11.5} Vxc​[{ρ}t0​≤s≤t​](r)≈VxcKS​[ρ(t)](r),(2.11.5)
其中$V_{xc}^{KS}$是基态计算的交换相关泛函。方程（2.11.5）也称为TDDFT中的绝热近似。虽然这种近似在计算广泛的电子和光学性质时工作良好，但它也可能在定性上失败，即使精确的基态交换相关泛函可用。

从实际角度来看，对于给定的交换相关泛函，TDDFT计算的一个主要挑战是小时间步长，其通常在阿秒量级（1阿秒=10-18秒）。与原子运动的典型时间尺度（1飞秒=10-15秒）相比，TDDFT的时间尺度要快得多，因此研究电子和原子核的自由度同时传播的非绝热动力学可能具有挑战性。

最后，让我们根据时间相关密度矩阵重写（2.11.4）：

其中，$P(t)$满足：

这就是（自洽）量子刘维尔方程（也称为冯·诺依曼方程），它可以被视为TDDFT的内在表示。

END…