VASP计算晶体弹性常数：Cij矩阵构建、应力-应变法与收敛陷阱解析

C11被低估40%的警告

弹性常数矩阵（Cij）是预测材料热力学稳定性、声子谱及力学行为的基础输入。本项目在计算MgO（岩盐结构）的弹性常数时，初始采用应力-应变法（IBRION=6），应变步长取0.005，结果C11=198 GPa，与实验值297 GPa偏差达33%。排查后发现，问题出在KPOINTS网格过稀：初始采用6×6×6的Monkhorst-Pack网格，而弹性常数计算对K点收敛的要求远高于静态优化。将K点加密至12×12×12后，C11回升至286 GPa，误差缩至3.7%。这个经历揭示了一个被普遍忽略的事实：弹性常数计算不是IBRION=6换一行代码那么简单，ENCUT和KPOINTS的收敛标准必须比静态优化再上一个台阶。

困境累积：应变步长的”Goldilocks”区域

应力-应变法的精度直接受应变步长影响。步长太小（<0.001），应力响应淹没在数值噪声中；步长太大（0.02），进入非线性区，胡克定律失效。本项目在Ti3AlC2（MAX相）的Cij计算中测试了5个步长：0.001、0.003、0.005、0.01、0.02。结果发现：0.001步长下C11=348 GPa，0.003步长下C11=346 GPa，两者相差0.6%，均在误差允许范围内；但0.02步长下C11骤降至312 GPa，偏差达10%。这是因为MAX相的层状结构在2%应变下已触发层间滑移，应力-应变关系偏离线性。最终项目组选定0.003作为标准步长，兼顾精度与稳定性。对于能量-应变法（IBRION=2+手动应变），步长可放宽至0.01，因为能量对应变的二阶导数更鲁棒，但计算代价是6个应变模式各做一次优化，总耗时增加5倍。

关键抉择：应力-应变法 vs 能量-应变法

两种方法各有适用边界。应力-应变法（IBRION=6）自动施加6个独立应变模式，适合立方、六方等低对称体系，但要求应力收敛极为严格（EDIFFG=-1e-3 eV/Å）。本项目在计算Al（FCC）时，应力-应变法耗时仅12分钟（8核），但FCC只有3个独立弹性常数（C11、C12、C44），IBRION=6仍按6个模式计算，效率浪费。改用能量-应变法后，只计算3个关键应变模式（εxx、εyy+εzz、εyz），耗时降至6分钟，结果一致。对于正交晶系（如BaTiO3），有9个独立Cij，IBRION=6的优势才显现。因此项目组的规则是：立方/六方用能量-应变法，正交/三斜用应力-应变法，单斜体系两种方法都不理想，建议改用有限差分法或DFPT（IBBRION=7/8）。

解决验证：三类典型体系的参数模板

经过6个体系的验证，本项目建立了弹性常数计算的参数矩阵：

验证结果：

– MgO：C11=286 GPa（实验297），C12=92 GPa（实验95），C44=152 GPa（实验156），误差均<4%

– Al：C11=108 GPa（实验114），C12=62 GPa（实验62），C44=28 GPa（实验32），偏差最大的是C44（12.5%），原因是Al的低频声子对K点极度敏感，需进一步用DFPT精修

– Ti3AlC2：C11=346 GPa，C12=67 GPa，C13=102 GPa，C33=288 GPa，C44=121 GPa，与Wang等人（PRB 75, 174111）的DFT结果偏差<2%

反思边界：温度效应与零温近似

VASP计算的是0 K弹性常数，而实验通常在室温测量。对于MgO这种热膨胀系数低的体系，0 K与300 K差异约5%，可以接受；但对于Ti3AlC2这种层状MAX相，层间剪切模量C44在300 K时下降约15%，因为热振动显著削弱了层间范德华耦合。项目组在Al的C44计算中额外做了准谐近似（QHA）修正：通过计算0 K、100 K、200 K、300 K的弹性常数，拟合得到dC44/dT=-0.04 GPa/K，300 K预测值为28-0.04×300=16 GPa，与实验32 GPa仍有差距。这说明QHA对于自由电子金属（Al）的剪切模量低估严重，需引入非谐效应（MD模拟）才能收敛。当前结果适用于0 K力学筛选和高压相稳定性预测，温度相关的弹性行为需补充有限温度MD或QHA计算。详细信息请参考科研学术网首页，或浏览VASP/DFT第一性原理栏目获取赝势与泛函选择策略。