本文是「第一性原理的微观计算模拟」系列的第 11 篇（共 11 篇），内容整理自同名书稿第 3 章，公式编号与原书一致。文中引用的文献依据原书参考文献清单整理并列于文末，按原书清单顺序从 1 开始编号。\nocite{*}
与其他方法相比，第一性原理计算方法最大的优势在于：① 可以给出准确的电子结构；② 对于成分复杂的体系，不需要拟合多参数的原子间相互作用势就可以得到可靠的能量。随着第一性原理计算方法在应用学科以及工程技术学科的广泛应用，其第二个优势在实践中显得尤为重要。但是，通常情况下，第一性原理计算方法在体系的空间尺度上有着较大限制，因此在能量计算上经常需要通过特殊的处理来阐述体系的稳定性以及其他性质。本节主要介绍四种常见的能量讨论或者处理方法。
缺陷形成能
在超单胞体系中，考虑带电量为 $q$ 的结构缺陷 $D$，其形成能 $\Delta H_{D,q}^{\mathrm{f}}(E_{\mathrm{F}},\mu)$ 定义为\cite{lany2008assessment}
$$ \Delta H_{D,q}^{\mathrm{f}}(E_{\mathrm{F}},\mu)=E_{D,q}-E_{\mathrm{ref}}+q(E_{\mathrm{V}}+E_{\mathrm{F}})+\sum_{i}n_i\mu_i \tag{3.655} $$
式中：$E_{D,q}$ 与 $E_{\mathrm{ref}}$ 分别为包含缺陷 $D$ 且带电量为 $q$ 的体系与完美体系的总能；$E_{\mathrm{V}}$ 为含缺陷体系价带顶能量。式（3.655）清楚地表明 $\Delta H_{D,q}^{\mathrm{f}}(E_{\mathrm{F}},\mu)$ 取决于缺陷的价态（$q$）、载流子类型及密度（费米能级 $E_{\mathrm{F}}$ 的位置）、材料合成环境（化学势 $\mu_i$）这三个条件。据此可以画出 $\Delta H^{\mathrm{f}}$ 随各条件变化的关系曲线图，进而可以展开非常细致和详尽的讨论\cite{lany2008assessment,vandewalle2004first}。
因为计算量的限制，一般而言超单胞不可能取得太大。目前比较常见的超单胞包含约 150 个原子。在这种尺寸下的单胞内即使只引入一个缺陷，缺陷浓度也将达到或接近 1%，这个值比实际浓度高好几个数量级。因此，为了使计算结果符合实验观测结果，必须计入某些修正项。
第一个修正与 $E_{\mathrm{V}}$ 有关。缺陷的存在可以严重地改变其附近的能带结构，而我们希望使用的则是远离该缺陷处的 $E_{\mathrm{V}}$\cite{vandewalle2004first}。一种可行的解决方法是用参考体系的价带顶能量 $E_{\mathrm{V}}^{\mathrm{ref}}$ 代替 $E_{\mathrm{V}}$。这种方法需要额外考虑如何使二者对齐。周期性边界条件（PBC）使得带电缺陷与其映像产生静电相互作用，由此而产生的额外的能量项使得整个实空间中的位势产生变化，从而导致缺陷体系的能级相对于参考体系有一个平移。这个平移量无法直接得出。因此 Van de Walle 和 Laks 等人提出了对齐两能级的方法\cite{vandewalle2004first,laks1992native}：设超单胞沿 $x$ 方向边长最大，分别计算带缺陷超单胞和完美超单胞的静电势分布，然后对 $Oyz$ 平面进行平均，得到 $V(x)$。取所有平面中距离缺陷最远的那个平面 $\boldsymbol{x}_0$，定义位势校正项 $\Delta V$ 如下：
$$ \Delta V=V_D(\boldsymbol{x}_0)-V_{\mathrm{ref}}(\boldsymbol{x}_0) \tag{3.656} $$
因此，$\Delta H_{D,q}^{\mathrm{f}}$ 可重新写为
$$ \Delta H_{D,q}^{\mathrm{f}}=E_{D,q}-E_{\mathrm{ref}}+\sum_{i}n_i\mu_i+q(E_{\mathrm{V}}+E_{\mathrm{F}}+\Delta V) \tag{3.657} $$
第二个修正也与周期性边界条件的使用有关。在带电缺陷及其映像的相互作用下，会产生专属于缺陷轨道／悬挂键交叠的能带 $E_D$。这与孤立带电缺陷的情况有偏差。因为在孤立带电缺陷的情况下，因缺陷而产生的应该是位于 $\varGamma$ 点处的平坦的孤立能级。因此 Wei 建议加入色散修正\cite{wei2004overcoming}：
$$ E_{\mathrm{dis}}=q(E_D(\varGamma)-E_D(\boldsymbol{k})) \tag{3.658} $$
式中括号内第一项是单 $\varGamma$ 点的计算结果，第二项是标准的 Monkhorst-Pack 多 $\boldsymbol{k}$ 点的计算结果。最后将 $E_{\mathrm{dis}}$ 加入式（3.657）中。
第三个常用的修正是 Makov–Payne 有限尺寸修正；此处写出的 $1/L$ 项是带电缺陷的单极项\cite{makov1995periodic}：
$$ E_{\mathrm{MP}}=\dfrac{q^2\alpha}{2\epsilon L} \tag{3.659} $$
式中：$\alpha$ 为 Madelung 常数；$\epsilon$ 为介电常数；$L$ 为超单胞的边长。Makov 和 Payne 指出，这个修正可以有效地提高能量计算相对于超单胞大小的收敛速度。 未修正的 $E_{D,q}$ 含周期映像的有限尺寸误差；应依据所用 Madelung 常数和形成能约定，显式施加相应修正，必要时还要考虑 $1/L^3$ 项。
上述三个修正都与有限大小的超单胞以及周期性边界条件有关。也就是说，如果体系足够大，这几个修正都可以不要。Castleton、Höglund 和 Mirbt 针对这种论点做了非常细致的研究\cite{castleton2006managing}。现行计算能力下不可能采用足够大的体系，因此他们利用多项式拟合计算无限大体系下的缺陷形成能 $\Delta H_0^{\mathrm{f}}$：
$$ \Delta H^{\mathrm{f}}(L)=\Delta H_0^{\mathrm{f}}+\frac{a}{L}+\frac{b}{L^{3}} \tag{3.660} $$
$\Delta H^{\mathrm{f}}(L)$ 是采用不同大小的超单胞得出的形成能。而所有这些超单胞必须保证形状完全相同（也即三个方向重复次数相同）。通过研究 InP 中十一种不同缺陷不同价态，他们得出结论：① 用多项式拟合式（3.660）是最为可靠和准确的方法，当然计算量也极大；② 对于单个超单胞计算，位势校正（见式（3.656））在大多数情况下可以给出合理的答案；③ 色散修正作用在受主态上的效果要优于其在施主态上的效果，但是总的来说，$E_{\mathrm{dis}}$ 并不可靠，有时候甚至会修正到相反的方向；④ Makov-Payne 修正在很多情况下无法给出比未修正的结果更好的结果，最能发挥效用的情况是原子弛豫较小的体系。
最后来讨论一下原子的化学势问题。对化合物而言，化学势是一个重要的环境变量，描述了体系所处的外部环境，例如各组分的贫富程度。因此 $\mu_i$ 可以在一个范围内变化，上限一般取为该元素单质处于稳定状态时单个原子的能量，超出这个限值，该元素的原子将会形成单质沉积下来，而不会形成化合物。以 Fe2O3 为例，$\mu_{\mathrm{Fe}}^{\max}=\mu_{\mathrm{Fe[bulk]}}$，即单质晶体中每个 Fe 原子的平均能量（极端富铁环境）。而 $\mu_{\mathrm{O}}^{\max}=E_{\mathrm{O_2}}/2$，即 O2 分子中每个 O 原子的平均能量（极端富氧环境）。特别需要注意的是，这里讨论的化学势，包括极值，本质上都是自由能，因此其大小取决于环境的温度及压强，很多情况下不可以将 $\mu_i^{\max}$ 简化为 0 K 下单质体系的内能，虽然很多情况下这种简化是合理的。为了将各组分的化学势联系在一起，还需要假设各组分永远与它们的化合物处于相平衡状态。同样以 Fe2O3 为例，则有
$$ 2\mu_{\mathrm{Fe}}+3\mu_{\mathrm{O}}=E_{\mathrm{Fe_2O_3[bulk]}} $$
由此可以确定各元素化学势的下限：
$$ \mu_{\mathrm{Fe}}^{\min}=\frac{1}{2}\left(E_{\mathrm{Fe_2O_3[bulk]}}-\frac{3}{2}E_{\mathrm{O_2}}\right) \tag{3.661} $$
$$ \mu_{\mathrm{O}}^{\min}=\frac{1}{3}(E_{\mathrm{Fe_2O_3[bulk]}}-2\mu_{\mathrm{Fe[bulk]}}) \tag{3.662} $$
化学势还可以用来确定杂质在材料中的溶解度。对于杂质 C，化学势 $\mu_{\mathrm{C}}$ 的下限为负无穷大，此时杂质在环境中的浓度为 0；上限则取为该元素单质的平均原子能量。但是通常情况下还必须考虑更严格的限制条件。这是因为杂质可能与材料中的某种元素结合成稳定的化合物而作为沉积物析出。Van de Walle 等人举了 Mg 掺杂于 GaN 中的例子\cite{vandewalle2004first}：Mg 可以占据 Ga 位而与 N 形成 Mg3N2，因此有条件
$$ 3\mu_{\mathrm{Mg}}+2\mu_{\mathrm{N}}=E_{\mathrm{Mg_3N_2}} \tag{3.663} $$
这相当于规定了 $\mu_{\mathrm{Mg}}$ 的新上限，且将其与 $\mu_{\mathrm{N}}$ 联系了起来。同时考虑电中性 Mg 在 Ga 位的形成能 $\Delta H_{\mathrm{MgGa},0}^{\mathrm{f}}$，且
$$ \Delta H_{\mathrm{MgGa},0}^{\mathrm{f}}=E_{\mathrm{MgGa},0}-E_{\mathrm{ref}}-\mu_{\mathrm{Mg}}+\mu_{\mathrm{Ga}} \tag{3.664} $$
将方程（3.663）和方程（3.664）联立起来，可以求得不同条件下 Mg 在 Ga 位上的最大溶解度，即最低形成能。Van de Walle 等人发现，在极端富氮条件下，Mg 在 Ga 位上的溶解度将达到最大\cite{vandewalle2004first}。
表面能
对于单质，表面能的计算公式比较简单：
$$ E_{\mathrm{surf}}=\frac{1}{2S}(E_{\mathrm{slab}}(N)-N\cdot E_{\mathrm{coh}}) \tag{3.665} $$
式中：$S$ 为表面积；$N$ 为拥有两个表面的体系（slab）所包含的原子数；$E_{\mathrm{coh}}$ 为该物质的聚合能。
化合物的情况要更复杂一些，因为其表面能与晶体在选定方向上的原子层堆垛情况有关。在大多数情况下，从任意位置分离晶体时，分离面两侧的两个表面不一致。以钙钛矿结构的 LaCoO3 为例，沿[001]方向，其原子层的堆垛顺序为 LaO-CoO2-LaO。因此，每一个(001)截面（不考虑重构）必然包含一个 LaO 面和一个 CoO2 面。在这种情况下表面能实际上应该由这两种面所共享。一般来说，对于沿某一方向呈 ABAB 形式堆垛的情况，需要考虑两个体系，其中体系 1 的两端面均为 $A$，而体系 2 的两端面均为 $B$。将这两个体系彼此首尾相接，可以构造出一个符合化学式的周期性单胞，记为体系 3。因此，体系 1 与体系 2 的总能之和与体系 3 的能量差实际上源于这四个表面。由此可以定义平均表面能为\cite{lee2009initio}
$$ \bar{E}_{\mathrm{surf}}=\frac{1}{4S}(E_{\mathrm{sys1}}+E_{\mathrm{sys2}}-E_{\mathrm{sys3}}) \tag{3.666} $$
式中：$E_{\mathrm{sys1}}$ 和 $E_{\mathrm{sys2}}$ 均为弛豫后的体系能量。
Eglitis 与 Vanderbilt 改写了方程（3.666），将每个表面的表面能分为刚性断裂项与弛豫项\cite{eglitis2007initio,eglitis2008initio}，有
$$ E_{\mathrm{surf}}(A)=E^{\mathrm{unr}}+E^{\mathrm{ref}}(A)=\frac{1}{4S}(E_{\mathrm{sys1}}^{\mathrm{unr}}+E_{\mathrm{sys2}}^{\mathrm{unr}}-E_{\mathrm{sys3}})+\frac{1}{2S}(E_{\mathrm{sys1}}-E_{\mathrm{sys1}}^{\mathrm{unr}}) \tag{3.667} $$
式中：$E^{\mathrm{unr}}$ 代表不经弛豫、原子均处于完美晶体的格点位置时的体系能量。这个公式强调了两种端面弛豫的差异。
表面巨势
表面巨势 $\varOmega$（有时也称表面自由能）经常与表面能同时使用，有时甚至比表面能更为重要。因为它可以描述不同生长条件下体系不同晶面的热力学稳定性，从而确定晶体生长的形状。不考虑熵的贡献，单位面积的 $\varOmega$ 定义为\cite{piccinin2008first,piccinin2010alloy}
$$ \varOmega=\frac{1}{2S}\left(E_{\mathrm{tot}}-\sum_{i}N_i\mu_i\right) \tag{3.668} $$
式中：$S$ 为所模拟体系的表面积，有系数 2 是因为厚板模型（slab model）在周期性边界条件下有两个相同的表面；$E_{\mathrm{tot}}$ 是体系的总能；$N_i$ 和 $\mu_i$ 分别代表第 $i$ 种原子的个数和化学势，一般而言，化学势 $\mu_i$ 是温度和压强的函数。体系处于稳态时 $\varOmega$ 最小，因此由 $\varOmega$ 可以预测给定条件下体系的表面组分、形态、吸附构型，以及颗粒形状等。详细的讨论可参阅文献\cite{piccinin2010alloy,bottin2003stability,reuter2004oxide}。本节中，我们按照文献\cite{bottin2003stability}中的思路，讨论钙钛矿结构的 LaCoO3 的最稳定表面。钙钛矿结构的低指数面比较复杂。为简单起见，这里仅考虑非重构的六个低指数表面：LaO-(001)、CoO2-(001)、O2-(110)、LaCoO-(110)、Co-(111) 及 LaO3-(111)。设体系恒与体相 LaCoO3 热平衡，则有
$$ \mu_{\mathrm{La}}+\mu_{\mathrm{Co}}+3\mu_{\mathrm{O}}=E_{\mathrm{LaCoO_3}} \tag{3.669} $$
式中：$E_{\mathrm{LaCoO_3}}$ 为每个 LaCoO3 立方单胞的能量。为了保证 La、Co 不在表面上以单质形式析出且 O 元素不以分子态逃逸，要求
$$ \begin{cases} \mu_{\mathrm{La}}\leqslant\mu_{\mathrm{La}}^{0}\\ \mu_{\mathrm{Co}}\leqslant\mu_{\mathrm{Co}}^{0}\\ \mu_{\mathrm{O}}\leqslant\mu_{\mathrm{O}}^{0} \end{cases} \tag{3.670} $$
如果进一步要求表面上不允许存在二元的金属氧化物，如 La2O3 及 CoO 等，则应引入不等式
$$ \begin{cases} 2\mu_{\mathrm{La}}+3\mu_{\mathrm{O}}\leqslant\mu_{\mathrm{La_2O_3}}=E_{\mathrm{La_2O_3}}\\ \mu_{\mathrm{Co}}+\mu_{\mathrm{O}}\leqslant\mu_{\mathrm{CoO}}=E_{\mathrm{CoO}} \end{cases} \tag{3.671} $$
平衡条件式（3.669）可以用来消除变量 $\mu_{\mathrm{Co}}$。所以不等式组（3.670）和不等式组（3.671）可以约化为
$$ \begin{cases} \mu_{\mathrm{La}}\leqslant\mu_{\mathrm{La}}^{0}\\ \mu_{\mathrm{La}}+3\mu_{\mathrm{O}}\geqslant E_{\mathrm{LaCoO_3}}-\mu_{\mathrm{Co}}^{0}\\ \mu_{\mathrm{O}}\leqslant\mu_{\mathrm{O}}^{0}\\ \mu_{\mathrm{La}}+\mu_{\mathrm{O}}\leqslant E_{\mathrm{La_2O_3}}+E_{\mathrm{CoO}}-E_{\mathrm{LaCoO_3}} \end{cases} \tag{3.672} $$
用 $\mu_{\mathrm{La}}$ 和 $\mu_{\mathrm{O}}$ 分别定义参考点 $\mu_{\mathrm{La}}^{0}$ 和 $\mu_{\mathrm{O}}^{0}$，其中 $\mu_{\mathrm{La}}^{0}$ 为单质 La 理想晶体的聚合能，而 $\mu_{\mathrm{O}}^{0}=E_{\mathrm{O_2}}/2$，即 O2 分子能量的均分值。这样，可以定义
$$ \begin{cases} \Delta\mu_{\mathrm{La}}=\mu_{\mathrm{La}}-\mu_{\mathrm{La}}^{0}\\ \Delta\mu_{\mathrm{O}}=\mu_{\mathrm{O}}-\mu_{\mathrm{O}}^{0} \end{cases} \tag{3.673} $$
类似地，可以定义 $\mu_{\mathrm{Co}}^{0}$ 为单质 Co 理想晶体的聚合能。将式（3.668）、式（3.669）和式（3.673）联立，可得
$$ \begin{aligned} \varOmega={}&\frac{1}{2S}[E_{\mathrm{slab}}-N_{\mathrm{Co}}E_{\mathrm{LaCoO_3}}-\mu_{\mathrm{O}}^{0}(N_{\mathrm{O}}-3N_{\mathrm{Co}})-\mu_{\mathrm{La}}^{0}(N_{\mathrm{La}}-N_{\mathrm{Co}})]\\ &-\frac{1}{2S}[\Delta\mu_{\mathrm{O}}(N_{\mathrm{O}}-3N_{\mathrm{Co}})+\Delta\mu_{\mathrm{La}}(N_{\mathrm{La}}-N_{\mathrm{Co}})] \end{aligned} \tag{3.674} $$
对于表面构型 $A$，式（3.674）可写成
$$ \begin{aligned} \Omega(A)={}&\frac{E_{\mathrm{slab}}^{A}-N_{\mathrm{Co}}E_{\mathrm{LaCoO_3}} -\mu_{\mathrm O}^{0}(N_{\mathrm O}-3N_{\mathrm{Co}}) -\mu_{\mathrm{La}}^{0}(N_{\mathrm{La}}-N_{\mathrm{Co}})}{2S}\\ &-\frac{\Delta\mu_{\mathrm O}(N_{\mathrm O}-3N_{\mathrm{Co}}) +\Delta\mu_{\mathrm{La}}(N_{\mathrm{La}}-N_{\mathrm{Co}})}{2S}. \end{aligned} \tag{3.675} $$
基于上述讨论，可以将方程（3.674）作为二元一次函数，在 $\Delta\mu_{\mathrm{La}}$ 和 $\Delta\mu_{\mathrm{O}}$ 所展开的平面上计算各个面的表面巨势，在取值许可的范围内最低的 $\varOmega$ 即为该条件下 LaCoO3 最稳定的表面。将式（3.668）至式（3.675）中出现的所有参量都利用 DFT 方法（如采用 VASP）求出，可知 $\varOmega$ 的取值范围由下列三个边界条件确定：
$$ \begin{gathered} \Delta\mu_{\mathrm{La}}\leqslant0\ \mathrm{eV},\quad\Delta\mu_{\mathrm{O}}\leqslant0\ \mathrm{eV}\\ \Delta\mu_{\mathrm{La}}+3\Delta\mu_{\mathrm{O}}\geqslant-11.23\ \mathrm{eV} \end{gathered} $$
图 3.21 给出了 LaCoO3 最稳定表面的相图。可见，在大多数情况下，LaO-(001) 都是最稳定的面。但是在富氧、贫镧条件下，LaO-(111) 面将转变为基态表面。
图 3.21 LaCoO3 最稳定表面的相图
注：左上角代表富氧-贫镧环境，右下角代表富镧-贫氧环境。黑色区域中，二元金属氧化物将在表面析出。
我们可以进一步研究不同环境下 LaCoO3 小颗粒的构型。对于三维体系，最稳定的构型由表面自由能 $F$ 极小值确定。$F$ 可以表示为表面巨势 $\varOmega$ 对体系表面的面积分，即
$$ F=\unicode{x222F}_{A(V)}\varOmega(\hat{\boldsymbol{n}})\,\mathrm{d}A \tag{3.676} $$
式中：$\hat{\boldsymbol{n}}$ 代表表面 $A$ 的法向。对晶体的微观模型而言，$\hat{\boldsymbol{n}}$ 基本不可能连续变化，因此小颗粒的构型一般为由若干个低指数面包围的多面体。具体确定这个多面体的形状则要用到 Wulff 构建法。这里不做理论上的讨论，仅仅给出实际操作步骤：设第 $i$ 个表面的法向为 $\hat{\boldsymbol{n}}_i$，计算表面巨势 $\varOmega_i$，则该表面过点 $\varOmega_i\hat{\boldsymbol{n}}_i$。所有这些面所包围的最小的闭合多面体即为给定条件下 LaCoO3 小颗粒最稳定的构型。图 3.22 给出了富氧、富镧条件下的 LaCoO3 表面稳定性相图及小颗粒的构型。当 $\Delta\mu_{\mathrm{O}}=0$ 且 $\Delta\mu_{\mathrm{La}}=-7.6\ \mathrm{eV}$ 时，$\varOmega_{(001)}(\mathrm{LaO})$ 和 $\varOmega_{(111)}(\mathrm{LaO_3})$ 比较接近，所以此时 LaCoO3 小颗粒是由六个{001}面和八个{111}面组成的十四面体。而当 $\Delta\mu_{\mathrm{La}}=0$ 且 $\Delta\mu_{\mathrm{O}}=-8.0\ \mathrm{eV}$ 时，$\varOmega_{(001)}(\mathrm{LaO})$ 明显小于其他各面的巨势，这时 LaCoO3 小颗粒是由六个{001}面组成的正六面体。
$\varOmega$ 还可以用于研究暴露于气体氛围中的金属／化合物表面形貌。Reuter 和 Scheffler 发展了一套普遍的方法来讨论这个问题。该方法称为受限热力学平衡（constrained thermodynamic equilibrium）方法，可参阅文献\cite{reuter2001composition,reuter2003first,reuter2003composition}。
图 3.22 LaCoO3 表面稳定性相图及小颗粒的构型
（a）富氧（$\Delta\mu_{\mathrm{O}}=0$）环境下 LaCoO3 表面稳定性相图及小颗粒的构型；
（b）富镧（$\Delta\mu_{\mathrm{La}}=0$）环境下 LaCoO3 表面稳定性相图及小颗粒的构型
集团展开与二元合金相图
集团展开（cluster expansion，CE）是一种在统计力学中被广泛使用的重整化方法。近年来，这种方法已经发展成为描述二元体系相互作用以及吸附分子间相互作用的一种极其有效的工具\cite{vandewalle2002alloy,vandewalle2002automating,vandewalle2002self}。
集团展开是一种基于晶格理论的方法，旨在通过简化相互作用力，实现对体系中的相互作用能量的有效描述。这种方法的基础在于将相互作用能量作为体系中原子配置的函数，然后用一个展开式来表示，其中每一项对应于不同的相互作用集团（例如单个原子、原子对、原子三元组等）。对于二元体系，可以利用集团展开方法构建一个有效的哈密顿量，用来描述不同相互作用的能量影响。此外，这种方法也可以应用于更复杂的体系，例如描述吸附分子间的相互作用。
在集团展开理论中，体系的总能 $E$ 可以表示为组成该体系的原子组态 $\boldsymbol{\sigma}$ 的函数：
$$ E(\boldsymbol\sigma)=J_0+\sum_iJ_i\sigma_i +\sum_{i\lt j}J_{ij}\sigma_i\sigma_j +\sum_{i\lt j\lt k}J_{ijk}\sigma_i\sigma_j\sigma_k+\cdots \tag{3.677} $$
式中格点的组态 $\boldsymbol{\sigma}$ 可以利用 Ising 模型或 Lattice Gas 模型（见第 8 章）表示，而 $E_0$、$V_i$、$V_{ij}$ 分别为空项、格点项、二体项。更高阶的 $V_{ij}$ 称为 $N$ 体项。这些项均为待定系数，统称为等效集团相互作用（effective cluster interaction，ECI）。因此，设有一个由两种元素共 $N$ 个原子组成的体系，总可以给定一个包含 $m$ 项（例如包含格点项、第一近邻二体项、第二近邻二体项、共线的三体项、非共线且所围面积最小的三体项、非共线且所围面积最小的四体项等等）的划分，并用方程（3.677）表达其总能。对于二元体系，容易看到可能的组态总数是 $2^{N}$，每一种组态对应一个能量。假设已有随机选取的 $l$ 个组态的能量，则可以由此构建一个 $l\times m$ 的系数矩阵 $\boldsymbol{C}$，并满足线性方程组
$$ \boldsymbol{CV}=\boldsymbol{E} \tag{3.678} $$
式中：$\boldsymbol{V}$ 是所有的 ECI 组成的 $m$ 维列矢量；$\boldsymbol{E}=[E_1,E_2,\cdots,E_l]^{\mathrm{T}}$。用最小二乘法拟合式（3.678），即可得到方程（3.677）所需的所有 ECI，进而用这些参量遍历所有 $2^{N}$ 个组态，找出二元体系在给定成分下的能量最低态，由此即可确定二元合金的相图\cite{zarkevich2004reliable,drautz2004ordering}。
习题
- 请简要概述第一性原理计算的基本原理，并用自己的话解释什么是密度泛函理论。
- 考虑一个简单的一维无限深势阱中的粒子，请使用第一性原理计算方法计算其能量本征值和波函数。这个问题的解析解是什么？将计算结果与解析解进行比较。
- 对于一个简单的二维平面正方形晶格，使用第一性原理计算方法计算其能带结构并描述计算过程。请解释布里渊区域的概念及这个晶格中电子的行为。
- 请简要描述赝势的概念，以及它在第一性原理计算中的作用。请比较赝势方法与全电子计算方法的优缺点。
- 由方程（3.390）推导 Perdew-Wang 以及 Vosko-Wilk-Nusair 形式的关联势 $V_{\mathrm{c}}(r_{\mathrm{s}})$。
- 推导式（3.564）。
- 证明式（3.561）与式（3.570）等价。
- 由赝波函数的正交性及 $\varPsi_{\mathrm{ps}}^{lm}(\boldsymbol{r})$ 是哈密顿量 $-\boldsymbol{\nabla}^{2}/2+V_{\mathrm{ps}}^{\mathrm{loc}}+\delta V_l$ 的本征函数证明方程（3.271）。
- 证明式（3.284）。