分子与固体的计算

本文是「分子与固体中的化学成键」系列的第 1 篇（共 10 篇），整理自 Richard Dronskowski 所著 Chemical Bonding: From Plane Waves via Atomic Orbitals 第 2 章的中文译本（DOI）。文中引用的文献依据原书参考文献清单（DOI）整理，列于文末。

分子的电子结构

假设我们希望计算一个由原子组成的微小系统的电子结构（因此，一切都由原子核和电子构成）。化学家通常把这样的实体称为分子，例如 H$_2$O、NH$_3$ 等。我们不考虑光谱学问题，只想正确求得分子的结构和总能量。在这种情况下，表征这个刚性分子的、与时间无关的波函数 $\Psi(\boldsymbol r)$ 就是我们的目标；我们需要求解定态薛定谔方程 (2.1) 来寻找它。这个方程提出于近一个世纪之前 (1926a; 1926b)\cite{schrodinger1926aquantisierung,schrodinger1926bquantisierung}：

$$ H\,\Psi(\boldsymbol r)=\left[-\frac{\hbar^2}{2m}\nabla^2+V(\boldsymbol r)\right]\Psi(\boldsymbol r)=E\,\Psi(\boldsymbol r)\tag{2.1} $$

哈密顿算符 $H$ 描述分子的全部能量，包括电子和原子核的动能与势能。对于 H$_2$O，有 3 个原子核，以及 $2\times1$ (H) $+8$ (O) $=10$ 个电子。用更堂皇的说法来表述，薛定谔方程对应于一个在 $3N$ 维空间中、具有奇异势的线性椭圆微分算子的本征值问题（其中 $N$ 为电子数，H$_2$O 中为 10），而自旋和泡利原理 (1925)\cite{pauli1925den} 又使它更加复杂。¹事实上，这个问题如此困难，以至于会让任何计算机都爆掉（对于具有合理规模的分子而言）。当然，我们可以先冻结原子核坐标，从而大幅简化问题；对于 H$_2$O，就是把三个原子核固定在空间中，让它们静止不动。这种所谓的玻恩–奥本海默近似 (Szabo & Ostlund, 1989)\cite{szabo1989modern} 是相当不错的选择，除非光谱学非常重要；但如前所述，我们并不关心这一点。接着，还可以把原子核之间的库仑排斥（一个常数项）设为零，待求得解之后再把它加回去。通过这样做之后，哈密顿量就只涉及电子了，这也是我们讨论所谓电子结构理论的原因：

$$ H=-\sum_{i=1}^{N}\frac{\hbar^2}{2m}\nabla_i^2-\frac{1}{4\pi\varepsilon_0}\sum_{i=1}^{N}\sum_A^M\frac{Z_Ae^2}{r_{iA} }{+\frac{1}{4\pi\varepsilon_0}\sum_{i=1}^{N}\sum_{j\gt i}^{N}\frac{e^2}{r_{ij} }}\tag{2.2} $$

请仔细看看式 (2.2)：等号之后只剩下电子的动能（左边一项）、电子与原子核之间的库仑势（中间一项），以及电子之间的排斥（右边一项）。电子用小写字母 $(i,j)$ 表示，原子核则用大写字母 $(A)$ 表示。真正的专家接下来会采用原子单位²，省去若干自然常数（$\hbar$、$m$、$\varepsilon_0$ 和 $e$）；不过，我们不再为这些细节费心，相关内容可以在分子或固体量子化学的书籍中找到 (Szabo & Ostlund, 1989; Mayer, 2003; Dronskowski, 2005)\cite{szabo1989modern,mayer2003simple,dronskowski2005computational}。为了着手处理一个复杂分子，先考虑一个简单的原子，把波函数 $\Psi(\boldsymbol r)$ 近似为各个电子占据的单电子波函数（也称为原子轨道）$\varphi_a(\boldsymbol r_1)$ 的所谓哈特里乘积，如下式所示：

$$ \Psi(\boldsymbol r)=\varphi_a(\boldsymbol r_1)\varphi_b(\boldsymbol r_2)\varphi_c(\boldsymbol r_3)\cdots\tag{2.3} $$

这其实是一个不错的试探形式，因为它使我们接近化学家所采用的准电子概念 (Primas & Müller-Herold, 1984)\cite{primas1984elementare}。还需要补充的是，上述原子轨道同样是薛定谔方程的解，只是其核电荷为某一数值（可能经过屏蔽），而且只有一个电子，即所谓的单电子问题。这些原子轨道集合，再加上构造原理（必须用电子正确地“填充”原子轨道），就能让我们理解元素周期表。³如果所考虑的是氧原子，这些轨道就叫作 $1s$、$2s$ 和三个 $2p$ 轨道，因此我们已经认识它们了。不过，回到前面提到的分子，我们会用原子轨道线性组合 (LCAO) 来近似分子单电子轨道 $\psi_i(\boldsymbol r)$；此时，原子轨道通常记作 $\varphi_\mu(\boldsymbol r)$：

$$ \psi_i(\boldsymbol r)={\sum_A\sum_{\mu\in A} }c_{\mu i}\varphi_\mu(\boldsymbol r)\tag{2.4} $$

符号 $\mu$ 遍历所有原子 $A$ 上的全部轨道。一些理论家有充分的理由强调，上面的 LCAO 表达式所用的原子轨道 $\varphi_\mu$ 最好已经针对分子中的情况作过适配或变形，因为分子中的势与原子中的势不同。因此，在 LCAO 的语境中，通常不再使用“原子轨道”这一术语，而以更一般的以原子为中心的基函数来代替，用它代表系统（分子）中的一部分（原子）。无论如何，它都是属于某个特定原子、或以该原子为中心的单电子函数。

电子是费米子，即自旋为 $1/2$ 的粒子，因此所得波函数必须具有反对称性；以全部单电子波函数构成的所谓斯莱特行列式可以保证这一点。泡利原理在这里起着支配作用。这就对应于哈特里–福克 (HF) 方法 (Fock, 1930a; Fock, 1930b)\cite{fock1930anaherungsmethode,fock1930bselfconsistent}，它正确地包含了 $N$ 个电子之间的全部交换相互作用：⁴

$$ \Psi(\boldsymbol r)=\frac{1}{\sqrt{N!} }\begin{vmatrix}\psi_1(\boldsymbol r_1)&\psi_2(\boldsymbol r_1)&\cdots&\psi_N(\boldsymbol r_1)\\\psi_1(\boldsymbol r_2)&\psi_2(\boldsymbol r_2)&\cdots&\psi_N(\boldsymbol r_2)\\\vdots&\vdots&\ddots&\vdots\\\psi_1(\boldsymbol r_N)&\psi_2(\boldsymbol r_N)&\cdots&\psi_N(\boldsymbol r_N)\end{vmatrix}\tag{2.5} $$

这种著名且极为成功的 HF 策略着眼于全电子波函数，对于不太大的分子，被认为是一个很好的选择。寻找这样的波函数，就相当于求解一个福克式的本征值方程：

$$ F\psi_i(\boldsymbol r)=\varepsilon_i\psi_i(\boldsymbol r)\tag{2.6} $$

其中，福克算符 $F$ 由单电子哈密顿量 $h$、直接库仑项 $J_j$，以及纯粹源于量子力学的交换库仑项 $K_j$ 构成。交换项反映了波函数的反对称性与双体算符相结合所产生的结果 (Szabo & Ostlund, 1989; Helgaker et al., 2000; Mayer, 2003)\cite{szabo1989modern,helgaker2000molecular,mayer2003simple}：

$$ F=h+\sum_{j=1}^{N}\left(J_j-K_j\right)\tag{2.7} $$

虽然式 (2.5) 中的单行列式方法非常灵活，但通常采用的（最简单的）求解方式会给出正交、正则且离域的轨道。如果需要从波函数中提取更多“化学”信息，也可以把这些正则轨道变换成更加局域的轨道。

电子相关与半经验方法

多电子 HF 波函数仍然缺少电子相关能 $E^{\mathrm{corr} }$（尽管斯莱特行列式已处理了自旋，无论自旋如何，电子仍会彼此排斥），因此，可以用 HF 能量与精确能量 $E^{\mathrm{exact} }$ 的偏差来衡量相关能的大小：

$$ E^{\mathrm{corr} }=E^{\mathrm{exact} }-E^{\mathrm{HF} }\tag{2.8} $$

这是 Löwdin 所给出的定义。改进 HF 结果有多种途径：可以引入多个斯莱特行列式，它们对应于电子在分子轨道中的不同占据情况，这就是所谓的（受限）组态相互作用 (CI)，它很容易从化学角度得到直观说明 (Coffey & Jug, 1974)\cite{coffey1974pedagogic}；也可以采用 Møller–Plesset (MP) 微扰理论方法，或者耦合簇理论（可直观理解为把 CI 改造成具有大小一致性的形式，但通常的截断耦合簇方法不具有变分性）。分子量子化学领域的文献极为丰富。如有需要，感兴趣的读者可以参考同样优秀的、专门介绍这些内容的教科书 (Szabo & Ostlund, 1989; Helgaker et al., 2000; Mayer, 2003)\cite{szabo1989modern,helgaker2000molecular,mayer2003simple}。

对于大分子，上述方法仅就计算量而言就变得难以处理（HF 理论的计算量按四阶增长，涉及相关轨道的方法则按五阶或更高阶增长）。因此，聪明的量子化学家不得不筛选 HF 理论中的积分，确定哪些足够重要而必须计算，哪些可以简化或直接舍弃，由此诞生了极为成功的半经验分子轨道理论 (Bredow & Jug, 2017)\cite{bredow2017semiempirical}，例如完全忽略微分重叠 (CNDO)。把研究对象扩展到晶体物质时，情况就更加严峻；不过，对于绝缘体和半导体，仍然可以进行包括周期性 MP 微扰修正在内的 HF 计算 (Usvyat et al., 2017)\cite{usvyat2017periodic}。对于其他固态材料，尤其是金属，HF 理论根本不是一个好的起点，因为从一开始，电子相关就同样重要。固体中的电子密度更高，这是由于其堆积更加致密（这里采用的是一个不够严格的论证），所以不能像处理许多简单分子那样忽略电子相关。正因为如此，只要涉及金属性，用周期性边界条件（见下文）运行 HF 理论来模拟晶体，可能就毫无帮助，HF 理论甚至会出现惊人的失败。⁵这也是分子量子化学与固体量子化学之间出现令人遗憾的分裂的主要原因。以固态形式出现的元素，大多数都具有金属性，这一点不应被忘记。

密度泛函理论

因此，固态理论家必须寻找另一种计算多体电子结构的方法。在物理学家的帮助下，他们确实找到了 (Inkson, 1986)\cite{inkson1986many}。结果表明，任何具有相互作用电子的晶体（也包括分子和原子，这一点同样值得感谢）的基态总能量，都可以表示为电子密度的泛函⁶：

$$ E\{\rho(\boldsymbol r)\}=T_0\{\rho(\boldsymbol r)\}+\int d\boldsymbol r\,V_{\mathrm{ext} }(\boldsymbol r)\rho(\boldsymbol r)+E_{\mathrm H}+E_{\mathrm{XC} }\{\rho(\boldsymbol r)\}\tag{2.9} $$

等号之后只有四个组成部分。首先是非相互作用电子的动能——所以记作 $T_0$——这些电子具有与相互作用电子相同的密度；

第二项是原子核产生的库仑势；第三项是简单的哈特里项⁷；第四项是交换–相关项，用于修正把电子描述成彼此不相互作用这一严重失实的假设。密度泛函理论 (DFT) 就是通过这种极为巧妙的方式发挥作用的。回过头来看，它取得成功的原因很容易理解。首先，$T_0\{\rho(\boldsymbol r)\}$ 相当大，但只要给定一组单电子波函数（轨道），就可以直接计算出来（即求二阶导数）；其次，修正项 $E_{\mathrm{XC} }\{\rho(\boldsymbol r)\}$ 往往相当小，除非材料的“相关性”太强——物理学家在这里可能会想到过渡金属氧化物等材料。正因如此，使用 $\Psi(\boldsymbol r)$ 的薛定谔多体方程，被改写成了使用 $\psi_i(\boldsymbol r)$ 的一组单电子 Kohn–Sham 方程。这些方程提出于 20 世纪 60 年代中期 (Hohenberg & Kohn, 1964; Kohn & Sham, 1965)\cite{hohenberg1964inhomogeneous,kohn1965self}：

$$ \left[-\frac{\hbar^2}{2m}\nabla^2+V_{\mathrm{eff} }(\boldsymbol r)\right]\psi_i(\boldsymbol r){=\varepsilon_i\,\psi_i(\boldsymbol r)}\tag{2.10} $$

在这里，只要有一个可靠的交换与相关泛函，有效势就可以处理多体效应。于是，我们又回到了一种便于解释的单电子理论，多体效应则被巧妙地“夹带”进这个框架中，成为一种经过强化的哈特里理论。事实上，这个想法多少有些顺理成章，因为 Slater (1951)\cite{slater1951simplification} 很早就认识到，HF 理论中的多体非局域交换势可以用电子密度的立方根来近似，这正是 DFT 的关键。

在进一步从交换–相关项的角度考察电子结构理论中这个真正伟大的思想 (Parr & Yang, 1989)\cite{parr1989density} 之前，先为感兴趣的读者提供一些文献，大概是合适的。尤其是，DFT 已经彻底改变了计算物理和计算化学，我们完全可以把它视为当代理论工作者最主要的工作工具。关于 DFT 的起源 (Jones & Gunnarsson, 1989)\cite{jones1989density}，以及对其未来前景的初步展望 (Jones, 2015; Springborg & Dong, 2017)\cite{jones2015density,springborg2017density}，都已有出色的概述。事实上，DFT 在日常研究中已变得如此普遍，以至于许多实际从事研究的科学家都可以讨论它的优缺点 (Teale et al., 2022)\cite{teale2022dft}，并列出大量最新文献。

回到至关重要的交换–相关项：当 $E_{\mathrm{XC} }\{\rho(\boldsymbol r)\}=0$ 时，我们仍然困在“哈特里地狱”里。当然，精确泛函是未知的，但局域密度近似 (LDA) 采用

$$ E_{\mathrm{XC} }^{\mathrm{LDA} }\{\rho(\boldsymbol r)\}=\int d\boldsymbol r\,\rho(\boldsymbol r)\epsilon_{\mathrm{XC} }\{\rho(\boldsymbol r)\}\tag{2.11} $$

对于许多密度缓慢变化的情况，这就已经给出了出乎意料的准确结果，虽然最初人们并没有预料到这一点。因此，在空间中的每一点，都加入一定的交换与相关贡献，其数值来自相同密度的电子气计算，并已制成表格。为了进一步改进，采用如下形式的近似广义梯度修正泛函 (GGA) 会很有帮助：

$$ E_{\mathrm{XC} }^{\mathrm{GGA} }\{\rho(\boldsymbol r)\}=\int d\boldsymbol r\,\rho(\boldsymbol r)\epsilon_{\mathrm{XC} }\{\rho(\boldsymbol r),\nabla\rho(\boldsymbol r)\}\tag{2.12} $$

这样，密度的梯度也进入了泛函表达式。接着可以采用 meta-GGA 泛函（通常还依赖 Kohn–Sham 轨道的动能密度；某些形式也使用密度的拉普拉斯算符），再进一步则是包含精确交换的杂化泛函（HF 理论又一次登场了），等等。如今，可供使用的密度泛函数量繁多，使 DFT 几乎带上了半经验的色彩；不过，实际使用的泛函种类相当有限，而且不同研究群体的选择也不同。分子量子化学家 (Koch & Holthausen, 2001)\cite{koch2001chemists} 可能更喜欢看起来略带经验性、但表现很好的 B3LYP 泛函，而不少固态理论物理学家仍然使用 LDA：简单而透明。不过，要正确得到（磁性）过渡金属各同素异形体的相对稳定性，物理学研究就必须采用 GGA (Dronskowski, 2005)\cite{dronskowski2005computational}。几乎不可能提出一个好的推荐，而又不被某些研究群体“干掉”，因此，提供一份合适的参考文献 (Teale et al., 2022)\cite{teale2022dft} 或许就足够了。

倒易空间与布洛赫定理

具体到固态，情况又如何？前面已经提到，DFT 诞生于固态研究（以及金属物理）。由于大多数固态材料以晶体物质的形式存在，我们可以利用晶体学家处理晶体时一直采用的、基本相同的数学方法，因此需要引入倒易空间。在倒易空间中，可以引入一个新的量子数，记作 $\boldsymbol k$，它既具有分数性，也具有方向性；还可以用晶体结构中的平移矢量 $\boldsymbol T$，对任何满足薛定谔方程或 Kohn–Sham 方程的波函数进行平移。换一种说法，具有平移不变性的固体，其全部电子结构都可以归回（或“反折叠”）到倒易空间的晶胞，即布里渊区中。这个固态科学中极为重要的定理最初由 Bloch 提出 (Bloch, 1928)\cite{bloch1928quantenmechanik}，形式看起来清爽而简单，但意义极为深远：

$$ \Psi_{\boldsymbol k}(\boldsymbol r+\boldsymbol T)=\Psi_{\boldsymbol k}(\boldsymbol r)e^{i\boldsymbol k\boldsymbol T}\tag{2.13} $$

借助平移对称性，我们不必计算一个由例如 1 mol（约 $6\times10^{23}$ 个）原子组成的体系，只需处理布里渊区即可，多么了不起的成就！因此，任何针对晶体物质的固态电子结构计算，都要在 $\boldsymbol k$ 空间中的不同点上对结果进行采样，这是使一个无限体系在某种意义上变得“有限”所付出的代价，而指数前因子起着决定性作用。它的作用类似于分子情形下 LCAO 的混合系数 $c_{\mu i}$（式 (2.4)）。图 2.1 示意了这一思想 (Deringer & Dronskowski, 2013)\cite{deringer2013computational}。

为求完整，再补充一点：如今，倒易空间或 $\boldsymbol k$ 空间的采样已是常规工作，实施方法也有多种。例如，可以选择“最佳”的采样点或最高效的网格 (Chadi & Cohen, 1973; Monkhorst & Pack, 1976)\cite{chadi1973special,monkhorst1976special}，随后采用最高效的积分方法 (Blöchl et al., 1994)\cite{blochl1994improved}。更多细节可参见其他文献 (Schwarz & Blaha, 2017)\cite{schwarz2017dft}；在我看来，这些内容相当技术化，虽然确实重要。

平面波基组

前面已经说过，平移一个波函数，只相当于把它乘以相位因子 $e^{i\boldsymbol k\boldsymbol T}$，仅此而已，这确实大大简化了问题。Bloch 本人也颇为欣喜地发现，通过简单的傅里叶分析就可以看出，这个波函数与自由电子的平面波之间，唯一的差别只是一个周期性调制。反映晶体平移不变性的，是周期函数 $u_{\boldsymbol k}(\boldsymbol r)$，满足 $u_{\boldsymbol k}(\boldsymbol r+\boldsymbol T)=u_{\boldsymbol k}(\boldsymbol r)$；$e^{i\boldsymbol k\cdot\boldsymbol T}$ 则是平移相位因子。这一区别，对基函数的选择具有重要影响。这种周期性调制意味着，平面波是适应晶体边界条件的对称性匹配波函数。因此，任何周期性波函数都可以表示为一系列指数函数（或者正弦、余弦函数）之和；这些函数从一开始就是完全离域的，求和延伸至某个晶格矢量 $\boldsymbol G$：

$$ \Psi_{\boldsymbol k}(\boldsymbol r)=\frac{e^{i\boldsymbol k\boldsymbol r} }{\sqrt{\Omega} }{\sum_{\boldsymbol G} }c_{\boldsymbol k}(\boldsymbol G)e^{i\boldsymbol G\boldsymbol r}\tag{2.14} $$

这带来了巨大的数学优势。例如，可以用一个简单的控制参数来调节基组的质量（以平面波的最大能量衡量），还可以精确计算原子间的 Hellmann–Feynman 力（也称为核梯度）。最大截断能可表示为

$$ \frac{\hbar^2}{2m}|\boldsymbol k+\boldsymbol G|^2\leq E_{\mathrm{cut} }.\tag{2.15} $$

另一方面，只有价电子才能得到高效处理，因为位于价轨道之下的类芯轨道具有强烈的节点特征，需要极其庞大的平面波基组。用一个简单的模型体系最容易说明这一点：由 Na 原子组成的一维晶体。⁸在布里渊区边缘，即特殊点 X 处，延展波函数如图 2.2 所示。

当我们从一个晶胞走到另一个晶胞（或从一个原子走到另一个原子）时，波函数像正弦波、余弦波或平面波一样，上下起伏，正负交替。这是平移不变性以及区边界——即倒易空间的 X 点——共同导致的结果（$k=\pi/T$，使波函数改变符号）。不过，在原子附近，波函数的振荡更加剧烈，呈现出节点，原因很简单：$3s$ 轨道与能量更低的芯轨道 $2p$、$2s$ 和 $1s$ 正交。如果要描述这种节点结构，就需要数量极其庞大、延伸至很高能量的平面波，这会很快把平面波理论“逼死” (Dronskowski, 2005)\cite{dronskowski2005computational}。

赝势与 PAW 方法

因此，需要赝势理论（或有效芯势理论）：用作用于价电子的软势代替真实的核势，并舍弃整个原子芯。为什么可以这样做？任何价电子不仅受到原子核的吸引，还会被位于原子核与价电子之间、能量更低的芯电子排斥。因此，剩下的只是一个较弱的赝势 (Hellmann, 1935; Hellmann, 1937a; Hellmann, 1937b)\cite{hellmann1935new,hellmann1937akvantovaya,hellmann1937beinfuhrung}。这一策略最初由 Hellmann 想到⁹，相当于把元素周期表的思想引入量子力学。它既有合理依据，又存在需要审慎看待的问题，因为实现它的方法有很多种，事实上有无限多种。实际使用的赝势可以是模守恒的（在指定参考态中保持截断半径以内的径向积分范数，并匹配截断半径以外的波函数）、能量守恒的（给出相同的本征值）、形状守恒的（距原子核一定距离之后给出相同的波函数），等等。这种任意性是一个问题，至少在支持全电子计算的人看来如此。

目前，所谓的投影缀加波 (PAW) 理论很可能是最好的赝势理论，因为其中的价电子部分在某种程度上“依托于”先前的全电子计算 (Blöchl, 1994)\cite{blochl1994projector}。实际上，PAW 理论是赝势理论与线性化缀加平面波 (LAPW) 理论的一次美好联姻；后者可以说是迄今发明的最精确的全电子能带结构方法 (Andersen, 1975; Schwarz et al., 2002)\cite{andersen1975linear,schwarz2002electronic}。对于位于某个原子球内的每个原子，PAW 赝势理论把价能级真实的、具有节点的波函数 $|\psi_i\rangle$ 定义为

$$ |\psi_i\rangle=|\widetilde\psi_i\rangle+\sum_{\mu,R}c_{\mu R}\left(|\varphi_{\mu R}\rangle-|\widetilde\varphi_{\mu R}\rangle\right)\tag{2.16} $$

它由同一个球内完全用平面波表示的赝波函数（等号后的第一项），加上一个修正项（第二项）构成；修正项包含全电子“分波”和赝“分波”（这是物理学家对数值轨道或基函数的称呼）。因此，如前所述，赝波函数与全电子轨道相衔接；PAW 赝势理论由此消除了前面提到的大部分任意性。具体细节见其他文献，其中的处理颇为巧妙复杂 (Goedecker & Saha, 2017)\cite{goedecker2017eliminating}。

平面波 DFT 计算流程

事情就是这样继续展开的。对于给定的固态材料，使用完全离域的平面波基组，在倒易空间中求解 Kohn–Sham 方程。为了做到这一点，只使用价轨道，并通过赝势理论（或 PAW 理论）巧妙地重新定义它们所感受到的（核）势。电子之间的全部相互作用，则由适用于具体问题的各种 DFT“方言”（即交换–相关泛函）来处理。对于真正的延展材料，采用的晶胞可以是晶体学晶胞，也可以是由它导出的、更小的原胞，或更大的超胞。研究分子对象时，也可以使用超胞，只需在分子周围留出足够的空白空间。随后，计算要反复进行，不仅为了达到自洽，还要在倒易空间的不同 $\boldsymbol k$ 点进行，以采样其性质。这样可以得到质量良好的密度（用于改进势）和态密度。另一种做法是沿倒易空间中的某些路径移动，考察自洽能级如何变化；正是这些能带结构计算，赋予了该方法它的名称。¹⁰在充分收敛的电子结构基础上——也就是说，能量相对于 $\boldsymbol k$ 空间采样和变分自由度（以平面波的一定截断能衡量的基组大小）已经收敛——可以进一步计算原子间作用力，求得（准谐）声子，这是获得有限温度下自由能所需的最重要的组成部分。例如，可以获得吉布斯自由能，并由此巧妙地导出热力学性质 (Stoffel & Dronskowski, 2017)\cite{stoffel2017lattice}，前提是不存在太多作为不稳定征兆的虚频声子。不过，这是另一个理论问题，这里不再进一步说明。

关于如何对完美的（以及有缺陷的）固态材料进行这类电子结构模拟，论述几乎数不胜数，其中也包括从深刻的化学视角出发的研究 (Bredow et al., 2009)\cite{bredow2009theory}。甚至还存在完全不同的途径，它们不从单电子图景出发，而是采用双电子图景 (Plekhanov & Tchougréeff, 2017)\cite{plekhanov2017resonating}。同样，这里也没有篇幅介绍那些突破 DFT 计算极限的超大规模模拟。为此，化学家必须找到把量子力学（用于体系中有意思的、“具有反应性”的部分）与分子力学（用于周围较少引人关注的外层部分）耦合起来的方法，而这类方法也确实存在 (Catlow et al., 2017)\cite{catlow2017quantum}。

结构模型与价的问题

到目前为止，我们已经讨论了电子结构理论，至少是其基本原理。剩下的“唯一”问题，就是提出一个好的化学组成（包含哪些原子？）和一个好的结构模型（这些原子在哪里？），之后就可以开始计算了。结构模型尤其重要。在处理烃类等简单分子时，构建模型往往相当容易，因为即使分子尚未被合成出来，我们也知道原子之间的那些“棒”：C 有四根，N 有三根，O 有两根，氢有一根。仅凭这一点，就很容易为这样的分子提出合理的结构模型。¹¹对于未知的固态材料，例如含有 3 个 Pu、11 个 Hf 和 7 个 Cr 原子的金属间化合物（我不知道是否存在这种化合物），其结构可能完全未知，因为没人知道 Pu 与 Hf、Hf 与 Cr，以及 Pu 与 Cr 之间分别有多少根“棒”。这个问题——（固态）无机化学中尚未解决的价的问题——使这门科学比（分子）有机化学复杂几个数量级，但没有人愿意公开承认这一点。

因此，可以尝试（或者说不得不先尝试）大量不同的结构模型，再通过进化算法、粒子群算法或相关算法，对它们进行筛选、杂交和增殖 (Oganov, 2010; Wang et al., 2010; Lonie & Zurek, 2011; Yu et al., 2017)\cite{oganov2010modern,wang2010crystal,lonie2011xtalopt,yu2017predicting}，以充分表征构型空间，随后按照能量对所有候选结构排序。希望我们寻找的材料就是能量最低的那一个，即基态结构，但其他结构也可能很重要。这一点无论怎样强调都不为过。

稳定性、亚稳态与动力学

虽然本章讨论的是纯粹的理论，但仍需提醒一句。由于这件事极其重要，确实应该放在这里，而不是放进附录。仅通过理论发现的稳定新化合物，往往被惯常地定义为总能量低于竞争相（起始材料或分解产物）的化合物；用简化的化学语言来说，它们就是放热化合物。¹²如果声子谱中没有虚频模，理论家会格外高兴，因为他们把这些材料称为谐近似下动力学稳定的材料，这又是一个稳定性判据。不过，那些并非放热、而是吸热的化合物（它们的能量更高而不是更低），则被认为不大可能、甚至不可能存在。

这种说法在理论上（或者物理上）听起来虽然正确，但从化学角度，也就是在真实世界中，却完全错了。化学家特别善于借助动力学“绕过”热力学，制备吸热（不稳定）的化合物。几乎整个有机化学研究的都是不稳定的（吸热的）化合物，它们本应自发分解成热力学终点 CO$_2$ 和 H$_2$O；然而，谢天谢地，这些分子（还有你和我）确实存在，因为阻止分解的活化势垒足够高。化学家把这类化合物称为亚稳态化合物。同样，无机分子肼 N$_2$H$_4$ 的吸热量为 $+159$ kJ mol$^{-1}$，即它是不稳定的，但它可以被制备出来，而且完全可以操作使用（处于亚稳态）。甚至固态炸药叠氮化铅 Pb(N$_3$)$_2$，只要足够谨慎，也可以定量制备——直到你对它的触碰过于猛烈，它便伴随着一声巨响释放出 556 kJ mol$^{-1}$，生成 Pb 和 N$_2$。上述能量判据会把这些有趣的化合物全部漏掉；它们比那些老掉牙、索然无味的 NaCl 和 MgO——热力学终点——有意思得多。

此外，许多亚稳态化合物确实具有虚频声子模。这虽然让理论家感到不安，大自然却很平静，让这些亚稳态化合物分解，有时很快，有时则缓慢地持续许多年。在后一种情况下，人们有充裕的时间，利用亚稳态化合物继续开展化学研究。不要忘记，自 Ostwald 的时代起，人们就已知道亚稳态相的存在，以及它们比稳定相更快结晶的倾向，这就是使熵产生最小化的 Ostwald 分步规则 (van Santen, 1984)\cite{santen1984ostwald}。有时，一种曾经被认为“稳定”的化合物，直到几十年之后，才因为偶然发现了能量更低的多晶型而显露出其亚稳性。具有类异质石墨烯层状结构的叠氮化铜 $\beta$-CuN$_3$ (Liu et al., 2015)\cite{liu2015cun} 就是一个很好的例子。

因此，未来的理论工作者最好不要过分看重上述能量（或热力学）判据，还应考虑动力学，尤其是活化势垒。实际存在的化合物未必是放热的；它们之所以存在，是因为阻止分解的活化势垒足够高。理论上还有许多工作要做！

系列导航

1. 分子与固体的计算（本文）
2. 分子与固体的分析
3. “七武士”：七类材料及其化学成键
4. 混合阴离子、复杂阴离子与复杂阳离子
5. 电池材料中的共价性与离子性
6. 分子晶体、氢键与其他次级成键
7. Zintl 相、金属、巡游磁体与钢
8. 高压下的奇特物质
9. 固体中的多中心成键
10. 普通化学与量子化学的几则趣味知识