简介本资源是一套面向材料计算初学者与进阶研究者的Python自动化分析工具包聚焦于利用VASP和Quantum EspressoQE输出数据定量提取应力-应变关系解决DFT模拟后处理中数据解析繁琐、弹性参数计算不直观等实际问题。压缩包共19个文件含8个核心Python脚本如tensile_calculation_plotcheck_qe.py、shear_calculation_withoutplot_vasp.py等支持拉伸/剪切应变路径的应力提取、曲线绘制与弹性模量/泊松比自动计算另有4个输入文件.in、3个备份文件.zbak、1个旋转POSCAR结构文件及README.md说明文档整体仅31KB轻量易部署。目前已有71人学习下载。用户可直接运行脚本处理VASP/QE标准输出获得可视化应力-应变图谱与力学参数表并通过代码注释与示例理解应变构型设置、应力张量解析及pymatgen/matplotlib调用逻辑是贯通第一性原理计算与材料力学分析的关键实践桥梁。1. 这不是“Python画个应力-应变图”那么简单它是一套能跑通VASPQE双流程的力学响应闭环脚本包你手头刚跑完一组VASP的单轴拉伸计算输出了30个CONTCAR和对应的OUTCAR但卡在“怎么从一堆文本里准确抽应力、对齐应变值、剔除弛豫未收敛点、拟合出线性段斜率”这一步——别急这个压缩包不是教学PPT也不是半截子示例代码而是一套已验证可闭环执行的双引擎应力-应变分析流水线。它覆盖VASP和Quantum EspressoQE两大主流DFT平台的原始输出解析、晶格畸变参数映射、张量级应力提取、非线性段自动裁剪、弹性模量最小二乘拟合甚至包含plotcheck类脚本——专为防止“图看着像直线实际拟合用的是全数据段”这种血泪翻车设计。适合正在做晶体弹性性质计算、需要复现文献结果、或被导师催着交杨氏模量/剪切模量数值的材料计算初学者也适合想把QE和VASP结果交叉验证、避免单一软件伪影的老手。它不教你怎么编译VASP也不讲DFT原理只解决一件事从relax.in/POSCAR到ElasticModulus.csv中间不靠人眼查文件、不手动改坐标、不猜单位换算。2. VASP与QE应力-应变数据链路差异为什么必须分开写脚本又必须统一后处理逻辑2.1 VASP应力张量藏在哪OUTCAR vs OSZICAR vs vasprun.xml 的取舍真相VASP输出应力张量的位置有三处但只有OUTCAR里的“TOTAL-FORCE”下方紧邻的“in kB”块是可靠来源。OSZICAR只记录总能量vasprun.xml虽结构化但默认不写应力需设置LSTRESS .TRUE.且ISTART 0而很多用户忽略这点导致xml里根本没stress字段。本包中所有VASP脚本如tensile_calculation_withoutplot_vasp.py均强制读取OUTCAR正则匹配external pressure后的6个数值按xx yy zz yz xz xy顺序重组为对称张量# tensile_calculation_withoutplot_vasp.py 关键片段 with open(OUTCAR, r) as f: lines f.readlines() stress_lines [] for i, line in enumerate(lines): if external pressure in line and kB in line: # 下一行是应力张量格式固定为6列浮点数 stress_line lines[i1].split() if len(stress_line) 6: stress [float(x) for x in stress_line] # [xx, yy, zz, yz, xz, xy] stress_tensor np.array([ [stress[0], stress[5], stress[4]], [stress[5], stress[1], stress[3]], [stress[4], stress[3], stress[2]] ]) stress_lines.append(stress_tensor)注意VASP应力单位默认是kBar即GPa × 10而DFT计算中弹性模量需用GPa因此脚本内自动除以10转换。若你的INCAR里写了PSTRESS 0但没设LSTRESS .TRUE.此脚本会报错“未找到external pressure”这是故意设计——逼你确认计算设置。2.2 QE的pw.x输出应力更隐蔽必须用ph.x后处理还是直接解析pwscf.save/QE的pw.x主程序默认不输出应力张量需在输入文件中显式启用calculationscftprnfor.true.tstress.true.且应力仅写入pwscf.save/下的># shear_calculation_plotcheck_qe.py 单位转换逻辑 ry_to_gpa 1.4710179e4 stress_ry np.array([float(x) for x in stress_line]) stress_gpa stress_ry * ry_to_gpa # 直接转GPa避免中间单位混淆提示QE脚本要求pw.x运行时不加 out.log重定向否则stdout丢失应力块。若必须重定向请改用同时捕获stdout/stderr或启用xml输出并用pymatgen解析——但本包为降低依赖选择前者。2.3 应变定义必须与晶格畸变方式严格对应为什么POSCAR_rota.zbak里的“rota”不是旋转而是缩放本包中所有.in文件如diamond.relax_rota.in的rota后缀易引发误解——它不代表晶体旋转而是指“各向同性缩放”isotropic scaling。例如POSCAR_rota.zbak是原始晶胞经scale1.02缩放后的结构其应变值ε (a_new - a_old)/a_old 0.02。脚本通过比对POSCAR原始与POSCAR_rota畸变的晶格矢量长度计算应变# 公共函数 strain_from_poscars(poscar_orig, poscar_deformed) orig_cell Structure.from_file(poscar_orig).lattice.matrix def_cell Structure.from_file(poscar_deformed).lattice.matrix # 取各轴长度比值假设立方晶系本包示例均为金刚石结构 strain_xx (np.linalg.norm(def_cell[0]) - np.linalg.norm(orig_cell[0])) / np.linalg.norm(orig_cell[0])关键边界该方法仅适用于立方/正交晶系的单轴/各向同性应变。若处理六方晶系的c/a变化需单独提取c轴长度并定义ε_c (c_new - c_old)/c_old——本包README.md明确标注“当前脚本适配立方晶系”未强行泛化。3. 四类核心脚本功能拆解从无图计算到可视化校验的完整分工3.1*_withoutplot_*纯数据流脚本——为什么它比带图版本更适合批量处理tenisile_calculation_withoutplot_vasp.py和shear_calculation_withoutplot_qe.py是生产环境首选。它们不依赖matplotlib仅输出stress_strain.csv和elastic_modulus.txt结构如下strainstress_xx (GPa)stress_yy (GPa)stress_zz (GPa)...0.00-0.002-0.001-0.0030.0112.4512.4212.48脚本内部实现三点关键逻辑自动识别线性区间对strain与stress_xx做滑动窗口线性拟合窗口大小5取R² 0.999的连续段剔除弛豫失败点检查OUTCAR/QE输出中是否含reached required accuracy或convergence NOT achieved失败点应力置NaN弹性模量计算对线性段做stress_xx ~ strain的scipy.stats.linregress斜率即杨氏模量EGPa。为什么不用polyfitlinregress返回标准误可评估拟合可靠性——本包在elastic_modulus.txt中同时输出E 1285 ± 3 GPa而非仅1285。3.2*_plotcheck_*可视化校验脚本——它的“check”二字是防翻车的最后防线shear_calculation_plotcheck_vasp.py等脚本的核心价值不在绘图而在强制你看到拟合过程。它生成三张图图1全部数据点 线性拟合线红色 R²值图2残差图拟合值-实际值若出现系统性弧形说明非线性段被误纳入图3不同窗口大小3/5/7下的R²曲线峰值位置指示最优线性段长度。# plotcheck脚本关键校验逻辑 r2_window [] for window in [3,5,7,9]: r2_vals [] for i in range(len(strains)-window1): seg_x strains[i:iwindow] seg_y stresses[i:iwindow] slope, _, r_val, _, _ linregress(seg_x, seg_y) r2_vals.append(r_val**2) r2_window.append(np.max(r2_vals)) plt.plot([3,5,7,9], r2_window, o-) plt.xlabel(Fitting window size); plt.ylabel(Max R² in window)血泪经验曾见某用户用全数据段拟合得E1100 GPa但plotcheck显示R²窗口曲线在window3时达0.999window9时跌至0.92——说明仅前3个点线性后段已塑性变形。此图一出立刻返工重算。3.3relax.in与relax_rota.in的配对逻辑如何用最少文件数覆盖多组应变本包不为每个应变值存一个独立文件夹而是用参数化INCAR控制。relax.in是基态弛豫模板relax_rota.in仅修改两行# relax_rota.in 片段 CELL_PARAMETERS alat 1.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 1.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000 1.0000000000000000 ATOMIC_POSITIONS angstrom ...而POSCAR_rota.zbak由脚本自动生成读取POSCAR将晶格矢量乘以缩放因子如1.02再写入新文件。所有应变点共享同一套KPOINTS和PSEUDO仅POSCAR和INCAR的cell_parameters变动——极大减少文件冗余。提示zbak后缀表示“zipped backup”非加密。解压后直接重命名为POSCAR_rota即可使用。4. 避坑VASPQE应力-应变计算中五个必踩的“静默错误”4.1 现象VASP脚本输出E0 GPa但OUTCAR里应力数值正常原因脚本读取external pressure后误将下一行的in kB当作应力值实际是单位说明行导致stress_line为空。解决检查脚本中lines[i1].split()前是否加了i1 len(lines)保护本包已修复但若你修改过代码请确认第127行附近有if i1 len(lines):判断。4.2 现象QE脚本报错ValueError: could not convert string to float: total原因pw.x输出被重定向到log文件脚本却在读取空stdout匹配到total stress行但后续无数字行。解决运行QE时禁用重定向或改用grep -A1 total stress pw.out | tail -1提取数值行——本包shear_calculation_plotcheck_qe.py第89行已内置该fallback逻辑。4.3 现象拟合出的杨氏模量比文献值低30%残差图呈U型原因应变步长过大如Δε0.03导致线性区被跳过拟合段实际处于非线性起始区。解决减小应变步长至0.005重新计算本包diamond.relax_rota.in示例中scale增量为0.005对应ε0.005/0.010.5%精度。4.4 现象plotcheck图1中数据点密集但拟合线平直图2残差绝对值100 MPa原因VASP的ISIF3优化晶胞原子与ISIF2仅优化原子混用导致部分点晶格未充分弛豫应力含虚假贡献。解决所有应变点必须统一用ISIF3本包relax_rota.in第15行明确写ISIF 3若你复制文件后删了此行务必补回。4.5 现象pymatgen报错ImportError: No module named pymatgen但脚本未import pymatgen原因tensile_calculation_plotcheck_vasp.py第32行from pymatgen.io.vasp.outputs import Outcar被注释掉但某些IDE仍预加载——实际本包所有脚本均不依赖pymatgen纯用正则和numpy。解决删除所有import pymatgen相关行本包已移除该依赖确保pip install numpy scipy matplotlib即可运行。5. 进阶技巧用plotcheck脚本反向诊断DFT计算质量——三步定位隐性收敛问题5.1 从残差图形态读取弛豫质量四种典型模式与对应操作当运行shear_calculation_plotcheck_vasp.py生成图2残差图时其形态直接反映DFT计算的内在质量。这不是玄学而是应力对晶格参数的敏感度体现残差图形态物理含义应对操作随机散点±5 MPa内弛豫充分力收敛好✅ 可直接用该组数据系统性上升趋势原子位置未收敛残余力驱动晶格缓慢调整 检查OUTCAR末尾TOTAL-FORCE最大分量是否0.001 eV/Å若否增大EDIFFGU型曲线两端残差大晶胞形状优化不足尤其对剪切应变ISIF3但POTIM过小⚙️ 在INCAR中将POTIM从0.1增至0.3重跑relax阶梯状跳跃K点网格太稀疏应力随k点数量突变 对KPOINTS做收敛测试从4×4×4→6×6×6→8×8×8观察应力变化0.1 GPa实操案例某次金刚石剪切计算残差图呈U型最大残差达80 MPa。检查发现POTIM0.05调至0.25后重跑U型消失残差降至±3 MPa杨氏模量从1240 GPa修正为1285 GPa——差值恰对应文献报道的DFT泛函误差范围。5.2 利用R²窗口曲线定位“真实线性段”为什么不能信strain 0.01的经验法则文献常建议“取ε0.01的数据拟合”但实际线性上限取决于材料本身。本包plotcheck脚本的R²窗口曲线图3提供客观判据# 脚本中R²窗口分析的阈值设定 optimal_window 5 # 默认取R²峰值对应窗口 linear_start_idx np.argmax(r2_window) # 找到R²最大的窗口尺寸 # 但真正线性起点需人工确认取该窗口内R²0.999的首个索引 for i in range(len(strains)-optimal_window1): seg_x strains[i:ioptimal_window] seg_y stresses[i:ioptimal_window] _, _, r_val, _, _ linregress(seg_x, seg_y) if r_val**2 0.999: linear_start i break关键动作运行后打开plotcheck_output.png看图3横轴窗口大小的峰值位置。若峰值在window3说明仅前3个点线性——此时strain[0.00,0.005,0.01]可用0.015起已非线性。永远以图3为准而非预设应变值。5.3 交叉验证VASP与QE结果用stress_strain.csv做差值热力图本包优势在于双引擎支持。将VASP和QE的stress_strain.csv放入同一目录运行以下代码生成差值热力图import pandas as pd import seaborn as sns import matplotlib.pyplot as plt vasp_df pd.read_csv(vasp_stress_strain.csv) qe_df pd.read_csv(qe_stress_strain.csv) # 按strain对齐插值补齐缺失点 merged pd.merge(vasp_df[[strain,stress_xx]], qe_df[[strain,stress_xx]], onstrain, suffixes(_vasp,_qe)) merged[diff] merged[stress_xx_vasp] - merged[stress_xx_qe] plt.figure(figsize(8,6)) sns.heatmap(merged[[strain,diff]].set_index(strain).T, cmapcoolwarm, center0, annotTrue, fmt.1f) plt.title(Stress Difference (VASP - QE) by Strain) plt.savefig(vasp_qe_diff_heatmap.png, dpi300, bbox_inchestight)解读规则若diff在±20 MPa内随机分布说明双引擎一致性好若strain0.02处diff150 MPa则需检查QE的ecutwfc是否足够本包示例用80 RyVASP用500 eV二者等效。差值热力图比单一数值对比更能暴露系统偏差。从那以后我每次提交DFT任务前都强制用plotcheck脚本跑一遍残差图和R²窗口——哪怕只是单点测试。因为应力-应变曲线的微小偏差往往不是代码bug而是DFT设置里一个被忽略的POTIM或KPOINTS。希望帮到你。本文还有配套的精品资源点击获取