做第一性原理弹性常数的同行都有一个共同体验：算个能带态密度可以一晚上搞定，但弹性常数一旦较真起来，两周出不了门。
弹性常数计算的麻烦不在于软件本身，而在于对精度的极致要求。弹性常数是能量的二阶导数——任何一个在能量和力计算中可以被忽略的微小误差，到了二阶导数层面都会被放大十倍。所以第一步不是打开VASP，是先把精度相关的参数做一次全方位的重新审视。
ENCUT值的选取是第一个容易栽坑的地方。很多人习惯用POTCAR里ENMAX的1.3倍，这个经验值对普通结构优化没问题，但对弹性刚度来说偏保守。实际经验是：对于含有过渡金属的体系，ENCUT至少取到ENMAX的1.5倍，如果有轻元素（C、N、O），还要再往上加。曾经算过一个含B的硼化物体系，ENCUT从500提到650eV之后，C₁₁从一个明显偏低的数值跳了将近8%。这个差距不是小修小补。
K点密度同样不能偷懒。做弹性常数时K点间距建议取结构优化时的70%——如果结构优化用了0.03的间距，那弹性矩阵阶段最好压到0.02左右。这个经验来自多次翻车之后的反省：Gibbs2的弹性代码在K点不够密的时候，应变-应力曲线的线性度会退化，最终C₄₄这种剪切模量被明显低估。
然后是应变幅度的选择。有限差分法的原理是用小应变扰动算出应力响应，但”小”到什么程度是个问题。经验值是正负1%为基线，但需要做一个收敛测试：从0.5%到2%，每隔0.25%算一组弹性常数，看哪个区间结果稳定。太小的应变会让数值噪声淹没信号，太大的应变会超出线性响应区。曾经有一组六方结构的数据，0.5%应变算出来的C₃₃连续三次都不收敛，放大到1.2%瞬间稳定。
这里有一个很多人忽略的细节：应变模式的对称性。六方晶系独立弹性常数是5个，立方是3个——如果施加的应变张量不对应这些独立模式，你就会在应力响应里混入多余的线性组合项。VASP的应变模式通过IBRION=6加NFREE来控制，但具体到每种晶系施加哪几组独立的应变张量，需要在OUTCAR里手动核对应力张量的Voigt分量是不是正确解锁了目标弹性常数。曾经帮一个课题组排查一个非对角线弹性常数C₁₂异常偏高的问题，最终发现是正交晶系的应变模式漏了一个剪切分量，导致C₁₂把C₆₆的信息混在了里面。
最后说一句不那么技术性的话：即便上述参数都收敛了，弹性常数的绝对值跟实验比差10%~15%是常态。这不是谁算错了——DFT算的是0K零压条件下的理论弹性张量，而实验数据通常是在室温有限压力下测的，中间隔了一个热膨胀修正。真正该对比的不是绝对值，是各向异性的趋势：Zener比、Pugh比、Poisson比的可靠性远比单个Cij的数值高。很多论文把C₁₁跟超声实验比对得严丝合缝，但C₄₄差了20%——这种选择性比对是审稿中最容易被挑出来的问题。
关于力学稳定性判据的补充：算完弹性常数之后，有一套Born稳定性判据需要过一遍。对立方晶系，C₁₁-C₁₂>0、C₁₁+2C₁₂>0、C₄₄>0三条必须同时满足。如果有一条违反，说明这个结构在当前条件下是力学不稳定的——不是计算错了，是材料本身的力学稳定性确实不支持这个结构存在。这一点在做高压相变预测时尤其重要，计算出的弹性常数先过力学稳定性检查，确认结构在给定压力下是力学上能独立存在的，再往下分析体模量和剪切模量。
https://www.keyanxueshu.com/
