5.2 有限应变求弹性张量
本课输入均为未执行的教学示例,数值是待测试的起点;图为原创示意,不是实际计算结果。使用授权 VASP 和 PAW 数据,替换所有占位符,记录程序版本并自行检验收敛。
5.2.1 模型、单位与记录
能量用 eV,长度用 Å,力用 eV/Å;1 kbar = 0.1 GPa。明确每原子、每分子、每原胞或每计算胞的归一化。记录 PAW 名称、版本、ZVAL、ENMAX 及允许的哈希;不分发 POTCAR。SCF 收敛只验证所选电子问题,目标物性的收敛仍需单独检查。
原创示意图;曲线与空白数据区用于说明或填写真实结果,不代表已运行计算。
5.2.2 完整案例:步骤、解释与检验
测量什么?
体积模量只描述均匀压缩;弹性张量则描述每一种小形变引起的各应力分量,包括拉伸与剪切。可以把晶体看成有方向性的弹簧网络,沿 x 难拉伸不代表难剪切。小应变下 σi=ΣjCijηj。本节使用工程 Voigt 应变 η=(εxx,εyy,εzz,2εyz,2εxz,2εxy),应力顺序为 (σxx,σyy,σzz,σyz,σxz,σxy)。剪切中的因子 2 若处理错误,会直接污染 C44。
“冻结离子”张量让原子随晶胞做仿射形变,不额外弛豫内部位置;“弛豫离子”张量允许原子在受应变晶胞内重新调整。内部调整通常降低能量并软化响应,二者不是可随意互换的结果。
实操顺序
起点晶胞应接近零压,其剩余应力明显小于施加应变产生的应力。记录参考应力。对六个工程应变分量分别施加正负应变,初步可测试 0.0025、0.005、0.01 三种幅度。用 F=I+ε 构造形变;工程剪切 ηj 对应对称非对角元 ηj/2。这是小应变处理,不能直接当作大应变本构关系。
每个受应变晶胞先可做冻结离子静态计算;若需要弛豫离子结果,用 ISIF=2 保持该晶胞不变而只优化原子,然后以统一电子参数做最终静态应力计算。对每个应力分量与应变拟合,或使用上面的中心差分。先查看完整矩阵中 Cij−Cji 的不对称程度,再按晶体点群施加约束。很大的不对称是误差线索,不要一开始就通过平均隐藏。
VASP 路线与符号
内置有限差分是对参考结构进行响应计算的另一条路线,不是把同一段输入再用于每个手工应变点。应核对所用版本支持的参数,并辨认 OUTCAR 中冻结离子、离子贡献及最终相关张量;不要随便复制出现的第一块矩阵。把 IBRION=6 改成 8 并不会自动得到完整的应变弹性张量。
VASP 的应力输出符号中,对角元为正表示压缩。若采用拉伸为正的连续介质约定,应统一变号并将 kbar 转为 GPa;还要把 VASP 打印的分量次序映射到本节的 Voigt 次序。可通过微小各向同性膨胀检查符号:相对于平衡态,恢复压力应向拉伸一侧变化。
正确解读张量
零外应力下,立方晶体的弹性稳定条件为 C11−C12>0、C11+2C12>0、C44>0。这只检验局域长波弹性稳定性,不能证明所有有限波矢声子都稳定。有限压力下需要使用相应的应力修正稳定条件。一般晶体不能机械地要求每个 Cij 都为正,因为非对角项可以是负值。
应变幅度与电子精度要联合收敛:δ 太小会将应力噪声放大,δ 太大则引入非线性弹性。比较不同幅度的斜率,提高截断能,并检查等价分量是否一致。二维材料要报告 N/m 的二维刚度或明确有效厚度;真空稀释后的 GPa 不是材料本征值。
练习:给立方晶体施加 xx 单轴小应变,判断哪些应力分别给出 C11、C12,再施加剪切核对因子 2。比较冻结离子与弛豫离子的差别,不要预设每个分量软化程度一样。
5.2.3 未执行输入与分析框架
本课输入均为未执行的教学示例,数值是待测试的起点;图为原创示意,不是实际计算结果。使用授权 VASP 和 PAW 数据,替换所有占位符,记录程序版本并自行检验收敛。
5.2.3.1 输入片段 1
# External finite-strain series: internal coordinates relax, strain is fixed
IBRION = 2
NSW = 150
ISIF = 2
EDIFFG = -0.002
# Save final stress using an additional NSW=0 static run
5.2.3.2 输入片段 2
# Alternative response calculation on the well-relaxed reference structure
IBRION = 6
ISIF = 3
NFREE = 2
POTIM = 0.015
NSW = 1