VASP 弹性常数与力学性质入门

想知道材料“硬不硬”“沿哪个方向更容易拉伸”，可以从弹性常数开始。最直观的做法是：在优化结构上施加几个小应变，计算应力，再拟合应力对应变的斜率。 本例使用金刚石结构硅，展示普通用户可以照着完成的计算流程。

1. 先决定原子要不要跟着放松#

两者都有用途，但不能混在同一个张量里。应变点不能使用完全放开的晶胞优化，否则外加变形可能被优化掉。 常用弛豫离子设置是 ISIF=2，只优化原子。

本例沿用结构优化教程的硅结构，使用VASP 6.6.1、PBE、550 eV截断、7×7×7 Gamma网格、EDIFF=1E-6；原子优化使用EDIFFG=-0.02 eV/Å。GPU分支采用NVHPC 25.3、CUDA 12.8与配套MPI编译。

2. 实际操作顺序#

1. 先优化未变形结构，检查残余压力与原子力。
2. 复制出−1%、−0.5%、0、+0.5%、+1%五个应变目录。改变晶格时保持分数坐标，让原子先随晶胞一起变形。
3. 分别施加x方向正应变和xy剪切。立方晶体只有C11、C12、C44三个独立常数；其他晶系需要更多独立变形。
4. 固定离子组设置 NSW=0；弛豫离子组设置下方参数。在每个目录准备匹配POTCAR，执行 sbatch run.sh。
5. 读取OUTCAR应力，把应变作横轴、应力作纵轴；用正负应变共同拟合斜率。
6. 用张量计算模量和方向性质，检查是否满足相应稳定性条件。

IBRION = 2
NSW = 80
ISIF = 2
EDIFF = 1E-6
EDIFFG = -0.02

应力符号与剪切定义必须一致。 VASP的OUTCAR采用压缩为正的应力约定；本例脚本经ASE转为拉伸为正，再换算为GPa。xy工程剪切γ对应对称应变矩阵中的εxy=εyx=γ/2；不要把二者混用，否则C44会差一倍。OUTCAR的kB表示kbar，10 kbar=1 GPa。

3. 硅的实际结果#

图中每个点均来自对应变形后的VASP计算。固定离子与弛豫离子的区别由实际结果展示，不能预先假设所有Cij都相同。这里的0.02 eV/Å适合入门流程；如果目标是很小的模量差别，尤其是内部弛豫很弱的小应变点，应检查更小力阈值、应变范围与k点是否改变结论。

4. 从Cij还能得到什么#

对于立方晶体，体积模量为 B=(C11+2C12)/3。剪切模量可分别做Voigt、Reuss平均，再取Hill平均。用B、G得到多晶近似的杨氏模量E与泊松比ν。这些平均值不等于某个单晶方向的杨氏模量。

立方晶体在零外压附近的常用稳定性条件是 C11−C12>0、C11+2C12>0、C44>0。满足条件只说明这一结构对小均匀变形稳定；是否存在不稳定声子、是否会分解成其他相，还要分别判断。非零外压时使用与外压一致的稳定性判据。

利用弹性模量和质量密度，还能在连续介质近似下估计纵波、横波声速及弹性Debye温度。这些是由弹性模型推导的量，不能替代完整声子色散。由B/G、泊松比或经验硬度公式讨论脆韧与硬度时，应写明是经验趋势，不能直接把它当成实验压痕硬度或断裂韧性。

5. 二维材料和理想强度怎样区别处理#

二维体系的三维应力含有真空高度，直接比较GPa会随真空变化。面内弹性常数通常转换为N/m：C₂D=C₃D×晶胞垂直高度，注意GPa与Å的单位换算，1 GPa·Å=0.1 N/m。

大应变拉伸/剪切需要沿指定加载方向逐点优化，并检查结构变化。固定横向晶格与横向应力为零是两种不同边界；曲线最大应力也只是当前理想模型的峰值，不是含缺陷材料的工程屈服强度。

下载输入、实算数据与图片（ZIP） · 下载44项结果表（CSV）。公开包不含POTCAR和WAVECAR。

官方参考：应力与晶胞自由度、有限位移与响应、VASPKIT力学分析。