vasp计算晶体弹性常数

我做过一个立方晶体弹性常数计算的项目，客户手里有一种新型高温合金的候选相，想快速判断它在室温下是不是力学稳定。我这个项目最核心的交付物就是左侧的能量-应变曲线和右侧的弹性常数矩阵热力图。

弹性常数 C_ij 描述的是晶体在小应变下的应力响应。对于立方晶系，独立的弹性常数只有三个：C11、C12、C44。我的做法是在优化后的平衡结构上加一系列小应变 ε，然后算体系总能量变化 ΔE/V，再用抛物线拟合得到二阶导数，进而求出 C_ij。左侧图里三条曲线分别对应三种应变模式：蓝色 ε1 对应单轴应变，用来提取 C11；绿色 ε1+ε2 对应双轴应变，和 C11+C12 相关；红色 ε4 对应剪切应变，直接给出 C44。每条曲线都是开口向上的抛物线，极小点在 ε=0 处，曲率越大说明对应刚度越大。

右侧图是拟合出来的弹性常数矩阵。C11=240 GPa，C12=135 GPa，C44=120 GPa。立方晶系的力学稳定性有三个 Born 判据：C11 > |C12|、C11 + 2C12 > 0、C44 > 0。我的结果三个都满足，而且 C44 比较大，说明这个相对剪切变形有较强的抵抗力。客户原来担心这种新相会不会像某些高温相一样 C44 接近零导致软模，结果 120 GPa 打消了这个顾虑。

具体操作上，我用 VASP 的 IBRION=2 做离子弛豫，先在零应变下优化晶胞；然后对每种应变模式取 7 个应变点（−3% 到 +3%），步长 1%，分别算总能。注意这里不是算应力，而是用能量-应变关系拟合，这样对应变模式的覆盖更完整，也避免了应力计算中 Pulay 应力的影响。拟合我用二次函数 E(ε) = E0 + 0.5·V0·C·ε²，其中 V0 是平衡体积。

踩坑经验有两个。第一个是应变范围不能太大：超过 4% 时高阶弹性项会混进来，抛物线拟合会低估或高估 C_ij；我试过 5% 范围，C11 偏低约 8%。第二个是 K 点和截断能：弹性常数对基组不敏感，但对 K 点收敛敏感，我最终用 600 eV 截断、12×12×12 K 点。第三个是应力法 vs 能量法：应力法算得快，但要求 INCAR 里 ISIF=2 且应变很小；能量法更稳，适合期刊发表。

工程经验上，我建议拿到 C_ij 后顺手算几个工程常数：体积模量 B = (C11+2C12)/3，剪切模量 G 用 Voigt/Reuss/Hill 平均，杨氏模量 E = 9BG/(3B+G)。我这个材料 B≈170 GPa、G≈95 GPa、E≈240 GPa，和 C11 接近，说明材料偏脆硬。最后我也提醒客户，0 K 的弹性常数比室温略高 5-10%，如果做结构强度评估，建议用有限温度分子动力学再修正。

更多 VASP 弹性常数和力学性质计算的案例，我汇总在 [vasp计算晶体弹性常数](https://www.keyanxueshu.com/category/dft/) 栏目里。[vasp计算晶体弹性常数](https://www.keyanxueshu.com/)

右侧弹性常数矩阵热力图还能读出更多力学稳定性信息。立方晶系的三个独立常数 C11、C12、C44 满足 Born 判据后，还可以算弹性各向异性因子 A = 2C44/(C11−C12)。我们这个案例 A = 2×120/(240−135) ≈ 2.29，大于 1 说明材料在 <100> 和 <111> 方向的剪切响应有差异。A=1 表示完全各向同性。对于多晶材料，工程上更关心 Hill 平均模量；对于单晶或取向生长的薄膜，A 的值会直接影响裂纹扩展方向。我建议在交付报告中增加一个“弹性各向异性”小节，把 A、B、G、E 一起列出。这样客户如果要做晶体取向优化，就有据可依。