4. 第一性原理计算框架

弹性常数描述的是晶体在微小应变下能量对形变的二阶响应，因此它天然属于“基态能量—应变关系”的问题。第一性原理计算之所以能给出弹性常数，是因为密度泛函理论（DFT）提供了从原子序数和晶格结构出发、不依赖经验拟合势即可获得体系基态能量与应力的途径。理解这条链路，比记住某个软件命令更重要：给定应变 → 自洽求解电子基态 → 得到能量或应力 → 拟合二阶导数 → 弹性常数。本节围绕这条链路，说明DFT基本思想、Kohn-Sham方程、交换关联泛函、平面波基组与赝势方法各自在弹性常数计算中扮演什么角色，以及哪些参数会直接影响结果的可靠性。

一、DFT基本思想：把多电子问题化为可解的密度问题

真实晶体是“原子核 + 大量电子”的多体系统，直接求解多电子薛定谔方程不可行。Hohenberg-Kohn定理给出两个关键结论：体系基态电子密度唯一决定外势，从而决定体系所有基态性质；存在一个关于电子密度的能量泛函，其最小值对应基态能量。于是，原本依赖3N个电子坐标的波函数问题，被转化为只依赖三维空间电子密度的问题。对弹性常数而言，这意味着我们不需要追踪每个电子的运动细节，只需在给定晶格形变下，求出体系的最低总能量E(ε)或应力张量σ(ε)。

但密度泛函本身并未给出能量泛函的显式形式，尤其是动能和交换关联部分。Kohn-Sham方案的作用，就是把最难处理的多体相互作用“打包”进交换关联泛函，而把其余部分写成可计算的形式。

二、Kohn-Sham方程：自洽求解电子基态

Kohn-Sham（KS）方程将相互作用电子系统映射为一个无相互作用电子在有效势中运动的辅助系统。其形式为：

[ -ħ²/2m ∇² + V_eff(r) ] ψ_i(r) = ε_i ψ_i(r)
V_eff(r) = V_ext(r) + V_H(r) + V_xc(r)

其中V_ext是原子核外势，V_H是电子间库仑排斥的Hartree势，V_xc是交换关联势。电子密度由占据态波函数模方求和得到：ρ(r)=Σ_i |ψ_i(r)|²。由于V_eff依赖ρ，而ρ又依赖ψ_i，KS方程必须自洽迭代求解：先猜一个初始密度，构造有效势，解方程得到新波函数和新密度，再判断是否收敛。收敛后，总能量可写为：

E_total = T_s[ρ] + E_ext[ρ] + E_H[ρ] + E_xc[ρ] + E_ion-ion

在弹性常数计算中，KS自洽求解是内层循环。外层循环则是对晶格施加一系列应变，例如对立方晶系施加正应变ε₁、ε₂、ε₃和剪切应变ε₄、ε₅、ε₆，每个应变下都重新做一次自洽计算，得到对应的E(ε)或应力。弹性常数由能量对strain的二阶导数定义：

C_ij = (1/V₀) ∂²E/∂ε_i∂ε_j |_{ε=0}

实际操作中，常用有限差分法：取一组小应变（如±0.5%、±1.0%、±1.5%），对每个应变计算能量，再用二次多项式拟合E(ε)，二次项系数乘以1/V₀即得C_ij。也可以直接用应力—应变关系：σ_i = Σ_j C_ij ε_j，通过计算不同应变下的应力张量，用线性拟合得到C_ij。两种方法在数值上应互相验证。

三、交换关联泛函：精度与适用性的核心选择

交换关联泛函是DFT中唯一需要近似的地方，也是弹性常数计算误差的主要来源之一。常用泛函包括：

• LDA（局域密度近似）：只依赖局域电子密度。对简单金属和离子晶体尚可，但通常高估结合能、低估晶格常数，导致弹性常数偏硬。
• GGA（广义梯度近似）：引入密度梯度，如PBE、PW91。对大多数金属和半导体，晶格常数和弹性常数比LDA更接近实验，是弹性常数计算的首选起点。
• meta-GGA与杂化泛函：如SCAN、HSE06。对强关联或需要更高精度的体系，可改善电子结构描述，但计算成本显著增加。弹性常数计算中，若GGA结果与实验偏差较大，可尝试SCAN或HSE06验证。

选择泛函时，应优先考虑与实验晶格常数和弹性常数已有数据对比。例如，对α-Fe，PBE给出的C₁₁、C₁₂、C₄₄与实验偏差通常在5%～10%以内；对金刚石，LDA会明显高估C₄₄。若研究的是合金或缺陷体系，还需注意泛函对磁性和电子关联的敏感性。

四、平面波基组与赝势：把无限问题变成有限可算问题

晶体具有周期性，Bloch定理允许将波函数写成平面波与周期函数的乘积。平面波基组天然适合周期性体系，且可通过截断能E_cut统一控制基组完备性。截断能越高，基组越完备，总能量和应力越收敛，但计算量也越大。弹性常数对截断能比总能量更敏感，因为它是能量的二阶导数。通常需要做收敛性测试：固定其他参数，逐步提高E_cut，观察C_ij变化小于1%～2%时，认为收敛。

然而，全电子平面波方法需要处理芯电子的剧烈振荡波函数，成本极高。赝势方法将芯电子与原子核一起冻结为有效离子实，只显式处理价电子。这样波函数在芯区变得平滑，所需平面波数量大幅减少。常用赝势包括模守恒赝势、超软赝势和投影缀加波（PAW）方法。PAW在弹性常数计算中尤其常用，因为它兼顾了精度和效率，对过渡金属和稀土体系也能给出可靠结果。

赝势的另一个关键参数是价电子配置。例如，计算Ti的弹性常数时，若只把3d和4s作为价电子，可能忽略3p半芯态的影响，导致C₄₄偏差。此时应使用包含3p半芯态的赝势。类似地，对碱金属和碱土金属，需注意是否包含p半芯态。

五、在弹性常数计算中的适用性与常见陷阱

DFT计算弹性常数的适用性取决于体系类型和参数选择。对金属、半导体和离子晶体，PBE+PAW+平面波方法通常能在5%～15%内重现实验弹性常数。但以下情况需要特别注意：

• 强关联体系：如过渡金属氧化物，GGA可能给出错误的磁基态和弹性常数，需引入DFT+U或杂化泛函。
• 范德华体系：层状材料如石墨，GGA无法描述层间弱相互作用，需加色散修正（DFT-D3、vdW-DF）。
• 有限温度：DFT给出的是0 K静态弹性常数。若需室温值，需考虑声子贡献和热膨胀，可通过准谐近似或分子动力学实现。
• 应变幅度：应变过大（>2%）会引入非线性，拟合二阶导数时需用更小应变；应变过小则能量差可能被数值噪声淹没。建议先做应变幅度收敛测试。

总结来说，DFT为弹性常数计算提供了从电子结构出发的可靠框架。Kohn-Sham方程负责自洽求解基态，交换关联泛函决定精度上限，平面波基组和赝势方法决定计算可行性。理解每个环节的物理意义和参数影响，才能在实际计算中做出合理选择，并对结果的可靠性做出正确判断。