第一性原理计算应力-应变曲线的完整工作流与弹性常数拟合

简介：本资源是一套面向材料计算初学者与科研实践者的Python自动化分析工具包，聚焦于VASP与Quantum Espresso（QE）第一性原理模拟后的应力-应变关系建模与力学参数提取。针对DFT计算中反复修改晶格、批量读取输出、拟合弹性常数等繁琐环节，提供即用型脚本与结构化流程，显著降低后处理门槛。压缩包共16个文件，含8个核心Python脚本（如tensile_calculation_plotcheck_vasp.py、shear_calculation_withoutplot_qe.py等，分别适配拉伸/剪切变形及VASP/QE双平台）、4个典型输入文件（.in格式，含diamond.relax_rota.in等带旋转应变的弛豫模板）、2个POSCAR结构文件及README.md说明文档，整体仅30KB，轻量易部署。已有932人学习下载，用户可直接复用代码解析应力张量、绘制应力-应变曲线、计算弹性模量与泊松比，并通过对比VASP与QE双引擎输出，深入理解不同赝势与k点设置对力学响应的影响。

1. 为什么应力-应变曲线不能只靠一个 VASP 或 QE 单点计算就画出来？

你手头有一份《使用VASP和QE计算应力和应变关系_Python_下载.zip》，解压后看到几个 Python 脚本、INCAR/KPOINTS/POTCAR 模板，还有几组不同晶格常数的结构文件——但直接跑完发现：应力值跳变、应变拟合不线性、弹性常数算出来是负的。这不是你代码写错了，而是掉进了第一性原理计算应力-应变关系最典型的认知陷阱：应力不是材料固有属性，而是特定构型在特定计算设置下的响应函数；而应变本身没有绝对零点，必须定义参考态、控制对称性、约束自由度，并在足够小但非无穷小的位移步长下采样。这个 ZIP 包本质是一套面向工程复现的“应力-应变闭环工作流”：它用 Python 统一调度 VASP 和 Quantum ESPRESSO（QE）完成多构型弛豫→单点应力提取→数据对齐→张量拟合→弹性常数导出。适合正在做金属/陶瓷/二维材料力学性质预测的计算材料方向从业者，尤其当你需要向实验组提供可比的杨氏模量、泊松比或 C₁₁/C₁₂ 值，而不是仅仅贴一张 DOS 图时。它不解决“怎么装 VASP”，但能让你避开 90% 因输入参数漂移、晶格缩放方式错误、应力张量坐标系混淆导致的翻车。

---

2. 从晶体结构出发：构建等效应变路径的三种可靠做法

应力-应变关系建模的第一步，从来不是写 Python，而是决定“你到底想测哪一种应变”。常见误区是直接对原胞做各向同性缩放（isotropic scaling），结果发现立方晶系的 C₄₄ 算不准、六方结构的 c/a 比失真、甚至正交晶系因对称性破缺引入虚假剪切应力。我们必须回到弹性理论本源：广义胡克定律 σᵢⱼ = Cᵢⱼₖₗ εₖₗ 中，独立弹性常数个数由晶体点群决定（立方 3 个、六方 5 个、三斜 21 个）。因此，构建应变路径必须满足两个刚性条件：① 每组应变构型保持原始空间群对称性（或明确降群）；② 应变类型覆盖所有独立 εₖₗ 分量。以下是经某高校计算材料实验室三年实测验证的三种落地路径：

2.1 立方晶系：6 组正交应变 + 1 组静水压应变

这是最稳妥的起点。对原始晶胞施加如下六种微小应变（|ε| ≤ 0.015）：

• ε₁₁ = ε, ε₂₂ = ε₃₃ = 0 （沿 a 轴拉伸）
• ε₂₂ = ε, ε₁₁ = ε₃₃ = 0 （沿 b 轴拉伸）
• ε₃₃ = ε, ε₁₁ = ε₂₂ = 0 （沿 c 轴拉伸）
• ε₁₂ = ε, 其余为 0 （面内剪切）
• ε₁₃ = ε, 其余为 0 （a-c 面剪切）
• ε₂₃ = ε, 其余为 0 （b-c 面剪切）
  再加一组静水压应变：ε₁₁ = ε₂₂ = ε₃₃ = ε（各向同性压缩/膨胀）。

提示：不要用scale_lattice直接缩放晶格矢量！必须用应变张量 ε 作用于原始基矢aᵢ→aᵢ'= (δᵢⱼ + εᵢⱼ)aⱼ，否则会破坏晶胞正交性，导致 VASP 报错ERROR: the lattice vectors are not linearly independent。

2.2 六方晶系：5 组应变需同步调控 c/a 比

六方结构有 5 个独立弹性常数（C₁₁, C₁₂, C₁₃, C₃₃, C₄₄），其应变设计必须耦合 a 轴与 c 轴变化。推荐采用以下组合（ε = 0.008）：

• ε₁₁ = ε₂₂ = ε, ε₃₃ = 0 （面内均匀拉伸，c 不变）
• ε₃₃ = ε, ε₁₁ = ε₂₂ = 0 （c 轴拉伸，a 不变）
• ε₁₁ = −ε, ε₂₂ = ε, ε₃₃ = 0 （面内剪切，保持体积）
• ε₁₃ = ε, ε₂₃ = ε, 其余为 0 （c 轴耦合剪切）
• ε₁₁ = ε₂₂ = ε₃₃ = ε （静水压，此时 c/a 自动弛豫）
  关键点：第 5 组必须开启ISIF = 4（仅优化 c/a，固定体积），否则 QE 会把 c/a 当作自由变量乱调，导致 C₁₃ 计算失真。

2.3 正交及低对称晶系：用strain.py自动生成合规构型

对称性低于四方的结构，手动设计易漏独立分量。ZIP 包中strain.py内置了空间群识别模块（基于 spglib），可自动判断点群并生成最小完备应变集。使用方式：

python strain.py --structure POSCAR --symmetry hP --strain-step 0.007 --output-dir strained_structures/

该命令会：① 读取 POSCAR 并检测空间群；② 根据--symmetry强制指定（如hP=六方，cF=面心立方）；③ 在strained_structures/下生成 5~21 个子目录，每个含POSCAR_strained和strain_tensor.txt。
逻辑说明：脚本内部调用spglib.get_spacegroup()获取 Hall 符号，再查表匹配国际符号（如Pnma→ 正交，独立常数 9 个），最后按 Voigt 标记法枚举所有非零 εᵢⱼ 分量。--strain-step是最大应变幅值，超过 0.015 会导致高阶项干扰线性拟合，低于 0.005 则信噪比不足。

---

3. VASP 与 QE 的应力输出对齐：为什么你拿到的 stress.txt 总是“对不上”？

即使应变构型完全一致，VASP 和 QE 输出的应力张量也常出现 0.5 GPa 量级偏差——这不是精度问题，而是二者默认参考系、单位、符号约定根本不同。若直接拼接两套数据拟合弹性常数，结果必然发散。必须做三重对齐：

3.1 坐标系统一：强制采用“晶胞坐标系”而非“笛卡尔坐标系”

VASP 的OUTCAR中in kB单位的应力张量（TOTAL-FORCE后的in kB行）默认是晶胞坐标系（即以晶格矢量为基底），而 QE 的matdyn.out或prefix.stress默认输出笛卡尔坐标系（单位 GPa）。不对齐会导致 C₁₂/C₂₁ 混淆。解决方案：

• VASP：无需改动，OUTCAR中应力已为晶胞系；
• QE：在pw.x输入文件中添加：

  CONTROL calculation = 'scf' tprnfor = .true. tstress = .true. SYSTEM ibrav = 0 ELECTRONS ...

  并在&SYSTEM块中显式指定celldm(1)和celldm(3)（对六方），确保晶胞参数被正确解析；最关键的是，在运行后用post_processing/stress_cart2cell.py脚本转换：

  # stress_cart2cell.py import numpy as np cell = np.loadtxt("cell.dat") # 3x3 晶格矩阵，单位 Bohr stress_cart = np.loadtxt("prefix.stress") # 3x3，单位 GPa # 转换公式：σ_cell = (V / V₀) * M⁻¹ σ_cart M⁻ᵀ，其中 M 是晶格矩阵，V₀ 是原始体积 M_inv = np.linalg.inv(cell) V0 = np.abs(np.linalg.det(cell)) V_new = V0 # 假设体积未变（小应变近似） stress_cell = (V_new / V0) * M_inv @ stress_cart @ M_inv.T np.savetxt("stress_cell.txt", stress_cell, fmt="%.6f")

  参数说明：cell.dat必须与当前计算的celldm严格一致；V_new/V0近似为 1，但若做静水压应变则需代入实际体积比。

3.2 单位与符号：全部转为 GPa 并确认压应力为负

VASPOUTCAR应力单位是 kB（1 kB = 10 kBar = 1 GPa），但数值前缀带负号表示压应力（compression），这与力学惯例一致；QE 默认也是 GPa，但部分旧版本输出正值为压应力。必须校验：对静水压压缩构型（ε = −0.01），应力对角元应全为正值（如 5.2 GPa），表示材料受压。若 QE 输出为 −5.2，则需全局取反：

awk '{print -$1, -$2, -$3}' prefix.stress > stress_positive.txt

注意：此操作仅适用于静水压或各向同性应变；剪切应变的符号必须保留原始方向性。

3.3 零点校准：用未应变结构的应力作为 baseline

即使完美弛豫，DFT 计算的零应变结构应力也不严格为零（常有 ±0.1 GPa 漂移），直接减去会导致拟合截距偏移。正确做法是：对每组应变构型，计算 Δσᵢⱼ = σᵢⱼ(ε) − σᵢⱼ(ε=0)，再用 Δσ 对 ε 线性回归。ZIP 包中fit_elastic.py默认启用此模式：

# fit_elastic.py 关键段 ref_stress = np.loadtxt("unstrained/stress_cell.txt") # 3x3 for eps_dir in glob("strained_*/"): stress = np.loadtxt(f"{eps_dir}/stress_cell.txt") delta_stress = stress - ref_stress # 逐元素相减 eps_vec = parse_strain_from_dir(eps_dir) # 解析出 6 维 Voigt 应变向量 X.append(eps_vec) y.append(delta_stress.flatten()) # 展平为 6 维应力向量 X, y = np.array(X), np.array(y) C = np.linalg.lstsq(X, y, rcond=None)[0].reshape(6,6) # 最小二乘拟合

参数说明：rcond=None避免病态矩阵警告；reshape(6,6)得到 Voigt 格式的弹性矩阵 Cᵢⱼ；最终 C 是对称矩阵，需强制C = (C + C.T)/2。

---

4. 避坑：应力-应变拟合中 4 个血泪经验换来的高频翻车点

4.1 现象：拟合 R² < 0.95，C₁₁ 出现负值

原因：应变步长过大（>0.015）激活了非线性项，或某组构型未完全弛豫（IBRION=2但EDIFFG过松）。VASP 中EDIFFG = −0.01仅保证力收敛到 0.01 eV/Å，对弹性常数而言太粗糙。
解决：将EDIFFG设为−0.001，并检查OSZICAR中最后 3 步的F=值是否全部 < 0.001；对 QE，conv_thr = 1.D-8且forc_conv_thr = 1.D-4。

4.2 现象：C₄₄ 显著偏离文献值（如 Al 的 C₄₄ 应为 26–28 GPa，算出 15 GPa）

原因：剪切应变构型（如 ε₁₂）破坏了原始对称性，但ISIF设置错误。例如对立方 Al 施加 ε₁₂ = 0.01 后，若用ISIF = 3（全自由），晶胞会自发调整为单斜，引入额外应力分量。
解决：剪切应变必须用ISIF = 2（仅原子位置弛豫，晶胞固定），并在POSCAR中显式写入变形后的晶格矢量（非原始矢量）。

4.3 现象：QE 计算中stresses输出全为 0.0

原因：tstress = .true.未写入&CONTROL块，或calculation = 'scf'时未启用tprnfor = .true.（QE 要求力与应力同时输出才计算应力）。
解决：严格按 3.1 节的CONTROL块模板书写；运行后检查prefix.xml中是否有<stress>标签。

4.4 现象：Python 脚本报错ValueError: shapes (n,6) and (6,) not aligned

原因：fit_elastic.py期望输入X为(n_samples, 6)的应变矩阵，但某组strain_tensor.txt只写了上三角（如ε₁₁ ε₂₂ ε₃₃ ε₁₂ ε₁₃ ε₂₃），而另一组写了全张量（9 个数），导致维度错乱。
解决：所有strain_tensor.txt必须统一为 6 维 Voigt 向量（顺序：11,22,33,23,13,12）。用validate_strain.py批量校验：

for f in glob("*/strain_tensor.txt"): vec = np.loadtxt(f) if len(vec) == 9: # 提取 Voigt 顺序：xx,yy,zz,yz,xz,xy voigt = [vec[0], vec[4], vec[8], vec[5], vec[2], vec[1]] np.savetxt(f, voigt, fmt="%.6f")

---

5. 用 Python 实现弹性常数的交叉验证与物理合理性检验

拟合出 6×6 的 Cᵢⱼ 矩阵只是开始，真正决定结果能否交付的关键，在于三重物理检验：热力学稳定性判据、各向异性指标、以及与实验/高精度计算的交叉比对。ZIP 包中的validate_elastic.py就是为此而生——它不画图、不输出报告，只做一件事：用 0.1 秒告诉你这组 Cᵢⱼ 是否可信。

5.1 热力学稳定性：5 条判据必须全部通过

对不同晶系，稳定性条件差异极大。例如立方晶系只需 C₁₁ > |C₁₂| 且 C₄₄ > 0；而六方结构需同时满足：

• C₁₁ > 0
• C₃₃ > 0
• C₄₄ > 0
• C₆₆ > 0
• C₁₁ − C₁₂ > 0
• (C₁₁ + C₁₂)C₃₃ − 2C₁₃² > 0
  脚本自动识别晶系并执行：

def check_stability(C, crystal_system): if crystal_system == "cubic": return (C[0,0] > abs(C[0,1])) and (C[3,3] > 0) elif crystal_system == "hexagonal": # C[0,0]=C11, C[0,1]=C12, C[0,2]=C13, C[2,2]=C33, C[3,3]=C44, C[5,5]=C66 c11, c12, c13, c33, c44, c66 = C[0,0], C[0,1], C[0,2], C[2,2], C[3,3], C[5,5] return all([ c11 > 0, c33 > 0, c44 > 0, c66 > 0, c11 > c12, (c11 + c12) * c33 > 2 * c13**2 ]) # 其他晶系类似...

若任一判据失败，脚本立即终止并打印：❌ FAILED: Hexagonal stability criterion (C11+C12)*C33 - 2*C13^2 = -0.8 GPa < 0。

5.2 各向异性：用 Zener 比率量化“有多不均匀”

Zener 各向异性比率 A = 2C₄₄ / (C₁₁ − C₁₂) 是立方材料的黄金指标。A = 1 表示完全各向同性；A > 1 表示剪切各向异性更强（如 Cu 的 A ≈ 3.2）；A < 1 则反常（可能计算误差）。脚本输出：

✅ Zener anisotropy ratio A = 2.87 (Cu literature: 3.21 ± 0.15) → Deviation: -10.6%, acceptable for LDA/PBEsol

参数说明：deviation基于 PBEsol 泛函的典型误差范围（±12%），若用 PBE 则放宽至 ±15%。

5.3 交叉验证：用 VASP 与 QE 结果互为基准

最硬核的验证，是让两套计算自证清白。脚本要求：

• VASP 拟合得 C_VASP（6×6）
• QE 拟合得 C_QE（6×6）
  然后计算相对差异矩阵 Dᵢⱼ = |C_VASPᵢⱼ − C_QEᵢⱼ| / max(|C_VASPᵢⱼ|, |C_QEᵢⱼ|)
  若所有 Dᵢⱼ < 0.12（12%），则标记✅ Cross-validation PASSED；否则列出超标项：

⚠️ Cross-validation WARNING: C12 difference = 18.3% (VASP: 62.1, QE: 50.8 GPa) → Check ISIF setting in shear strain jobs

这比任何单点测试都可靠——因为 VASP 和 QE 的误差来源完全不同（赝势 vs PAW，平面波截断 vs 赝势投影），它们一致的地方，大概率就是物理真实。

我带过的某跨平台系统项目里，曾因忽略 Zener 比率检验，把一组因ENCUT不足导致的 C₄₄ 偏低数据当真，后续做分子动力学时剪切模量失配，浪费两周重算。现在我的习惯是：fit_elastic.py运行完，第一眼必看validate_elastic.py输出的 ✅/❌，第二眼扫 Zener 偏差，第三眼盯交叉验证的 C12。只要这三关过，这组弹性常数就能进论文 Supporting Information。希望帮到你。

本文还有配套的精品资源，点击获取