简介本资源是一套面向材料计算初学者与科研实践者的Python自动化分析工具包聚焦于VASP与Quantum EspressoQE第一性原理模拟后的应力-应变关系建模与力学参数提取。针对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 ESPRESSOQE完成多构型弛豫→单点应力提取→数据对齐→张量拟合→弹性常数导出。适合正在做金属/陶瓷/二维材料力学性质预测的计算材料方向从业者尤其当你需要向实验组提供可比的杨氏模量、泊松比或 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应力单位是 kB1 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, rcondNone)[0].reshape(6,6) # 最小二乘拟合参数说明rcondNone避免病态矩阵警告reshape(6,6)得到 Voigt 格式的弹性矩阵 Cᵢⱼ最终 C 是对称矩阵需强制C (C C.T)/2。4. 避坑应力-应变拟合中 4 个血泪经验换来的高频翻车点4.1 现象拟合 R² 0.95C₁₁ 出现负值原因应变步长过大0.015激活了非线性项或某组构型未完全弛豫IBRION2但EDIFFG过松。VASP 中EDIFFG −0.01仅保证力收敛到 0.01 eV/Å对弹性常数而言太粗糙。解决将EDIFFG设为−0.001并检查OSZICAR中最后 3 步的F值是否全部 0.001对 QEconv_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₁₁ 0C₃₃ 0C₄₄ 0C₆₆ 0C₁₁ − 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 (C11C12)*C33 - 2*C13^2 -0.8 GPa 0。5.2 各向异性用 Zener 比率量化“有多不均匀”Zener 各向异性比率 A 2C₄₄ / (C₁₁ − C₁₂) 是立方材料的黄金指标。A 1 表示完全各向同性A 1 表示剪切各向异性更强如 Cu 的 A ≈ 3.2A 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_VASP6×6QE 拟合得 C_QE6×6然后计算相对差异矩阵 Dᵢⱼ |C_VASPᵢⱼ − C_QEᵢⱼ| / max(|C_VASPᵢⱼ|, |C_QEᵢⱼ|)若所有 Dᵢⱼ 0.1212%则标记✅ 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。希望帮到你。本文还有配套的精品资源点击获取