简介：这份资源面向材料科学计算方向的研究生与科研人员，聚焦如何用Python驱动VASP与Quantum Espresso完成应力—应变关系计算，解决第一性原理力学性质模拟中数据提取、处理与可视化的问题。压缩包共16个文件，约30KB，以8个Python脚本为核心，配合4个输入文件、POSCAR结构文件及README说明，覆盖拉伸与剪切两类计算场景，并区分VASP与QE两套流程，脚本还提供是否绘图的可选版本，便于按需调用。已有931人学习下载，说明其在同类计算任务中具备一定参考价值。读者可从中获得读取输出文件、提取应力应变数据、绘制曲线并进一步拟合弹性模量与泊松比等参数的完整脚本框架，同时借助示例输入与结构文件快速复现计算流程，理解两种DFT软件在应变模拟中的衔接方式，适合作为力学性质计算的入门模板与二次开发起点。
1. 从一份能跑的应力应变脚本包说起：VASP 与 QE 双引擎的 Python 落地
很多人第一次算弹性常数，卡住的不是 DFT 本身，而是「应变怎么加、应力从哪读、数据怎么对齐」。VASP 和 Quantum Espresso 各自输出格式不同，一个走 OUTCAR 里的应力张量，一个走 pwscf 的 output 或 XML，手动抄数据基本等于自找麻烦。这个压缩包 StrengthCalculation-strain-stress-master 把拉伸和剪切两条线都拆成了独立脚本，VASP 和 QE 各一套，还配了 diamond 的 POSCAR、relax 输入和旋转版本，等于把「改晶格 → 跑静态 → 提应力 → 拟合」整条链路摊开给你看。适合已经装好 VASP 或 QE、能跑通单点能计算、但还没把弹性常数流程串起来的材料计算从业者。Python 在这里不是主角，是胶水，负责把两个引擎的输出粘成一条应力应变曲线。
2. 脚本包结构拆解：拉伸与剪切两条线怎么分
2.1 文件命名里的信息量
拿到压缩包先别急着跑，把文件名过一遍就能看出作者的意图。tensile_calculation_withoutplot_vasp.py和tensile_calculation_plotcheck_vasp.py是一对，前者只算不画，后者带绘图检查；QE 侧同样有tensile_calculation_withoutplot_qe.py和tensile_calculation_plotcheck_qe.py。剪切线是shear_calculation_plotcheck_vasp.py、shear_calculation_withoutplot_vasp.py、shear_calculation_plotcheck_qe.py、shear_calculation_withoutplot_qe.py。这种「withoutplot / plotcheck」的拆分很实用：批量跑的时候用 withoutplot 省时间，调参阶段用 plotcheck 当场看曲线是否线性。
结构文件有POSCAR和POSCAR_rota，对应relax.in、relax_rota.in、diamond.relax.in、diamond.relax_rota.in。rota后缀说明作者考虑了晶格取向问题——金刚石结构在不同晶向下弹性响应不同，旋转后的 POSCAR 用来验证各向异性。README.md 是唯一的文档入口，.gitattributes说明这包是从 git 仓库导出的，版本管理痕迹还在。
2.2 拉伸与剪切的物理区别在脚本里怎么体现
拉伸计算改的是晶格常数，沿某个方向拉长或压缩，其他方向可能固定也可能按泊松比松弛。剪切计算改的是晶格矢量之间的夹角，POSCAR 里表现为基矢的非对角项变化。脚本里对这两种形变的处理逻辑不同：拉伸通常只动一个晶格参数，剪切要构造完整的形变矩阵。
常见做法是定义一个形变梯度矩阵，对原始晶格矢量做线性变换。拉伸对应对角矩阵，剪切对应非对角元非零的矩阵。脚本里如果直接改 POSCAR 的缩放系数，那只适合各向同性拉伸；要算完整的弹性常数矩阵，得按 Voigt 记号逐个施加应变模式。
提示：先确认脚本里应变的定义是工程应变还是真应变。小应变下两者差别不大，但应变加到 2% 以上时，拟合出的弹性常数会有可观测的偏差。
2.3 输入文件与脚本的对应关系
relax.in是 QE 的输入模板，diamond.relax.in是金刚石的具体算例。VASP 侧没有单独的 INCAR 模板，说明脚本可能直接生成 INCAR 或者依赖你手动准备。POSCAR 是 VASP 的结构文件，QE 用CELL_PARAMETERS和ATOMIC_POSITIONS卡片，两者格式不通用，脚本里应该有转换逻辑或者分别读取的分支。
跑之前建议先手动跑一个应变点，确认 VASP 或 QE 能正常输出应力。VASP 看 OUTCAR 里in kB那行的应力张量，QE 看 output 里total stress段落。如果这一步就报错，后面脚本再对也没用。
3. 环境准备与单点验证：跑脚本前必须过的三道关
3.1 VASP 与 QE 的最小可跑配置
VASP 需要 POTCAR、POSCAR、INCAR、KPOINTS 四个文件。POTCAR 按元素顺序拼接，POSCAR 用包里的，INCAR 至少设IBRION = -1（静态计算）、ISIF = 2（算应力但不改结构）、NSW = 0。KPOINTS 用 Gamma 点或 Monkhorst-Pack 网格，金刚石结构用 8×8×8 起步比较稳。
QE 的输入文件是单个.in，里面包含&CONTROL、&SYSTEM、&ELECTRONS三个 namelist，后面跟ATOMIC_SPECIES、CELL_PARAMETERS、ATOMIC_POSITIONS、K_POINTS卡片。relax.in里calculation = 'scf'还是'relax'决定了跑的是单点还是结构优化。算应力必须用tprnfor = .true.和tprnstr = .true.，否则 output 里没有应力张量。
# VASP 单点验证：确认 OUTCAR 里有应力输出 mpirun -np 4 vasp_std > vasp.log grep "in kB" OUTCAR # 应该看到类似：in kB -1.234 0.567 0.567 0.000 0.000 0.000# QE 单点验证：确认 output 里有 total stress mpirun -np 4 pw.x -in diamond.relax.in > qe.log grep -A 3 "total stress" qe.log # 应该看到 3x3 应力张量，单位是 Ry/bohr^3 或 GPaVASP 的应力单位默认是 kB，1 kB = 0.1 GPa，换算时别搞错。QE 的应力单位取决于&CONTROL里的设置，常见做法是在后处理时统一转成 GPa。
3.2 Python 依赖与 pymatgen 的取舍
脚本大概率依赖 numpy 和 matplotlib，可能还用了 pymatgen 读 OUTCAR。pymatgen 的Outcar类能直接提取应力张量，省去手写正则的麻烦。但 pymatgen 安装有时会卡在依赖上，如果只是读几个数，用正则匹配 OUTCAR 更轻量。
# 轻量读取 VASP OUTCAR 应力张量的常见写法 import re import numpy as np def read_stress_outcar(filename): """从 OUTCAR 提取最后一个应力张量，返回 3x3 numpy 数组""" with open(filename, 'r') as f: lines = f.readlines() # 从后往前找，取最后一个 "in kB" 后面的数据 for i in range(len(lines) - 1, -1, -1): if 'in kB' in lines[i]: # 应力值在下一行，6 个分量：xx yy zz xy yz zx vals = [float(x) for x in lines[i + 1].split()] # 组装成 3x3 对称矩阵 stress = np.array([ [vals[0], vals[3], vals[5]], [vals[3], vals[1], vals[4]], [vals[5], vals[4], vals[2]] ]) return stress * 0.1 # kB 转 GPa raise ValueError("OUTCAR 里没找到应力数据")这段代码的关键在应力分量的顺序。VASP 输出的是 Voigt 记号：xx, yy, zz, xy, yz, zx。组装成 3×3 矩阵时，xy 对应 (0,1) 和 (1,0)，yz 对应 (1,2) 和 (2,1)，zx 对应 (2,0) 和 (0,2)。顺序搞反了，剪切分量会错位，拟合出的 C44 会偏。
3.3 应变步长的选择与收敛测试
应变步长不是越小越好。步长太小（比如 0.1%），应力响应可能被数值噪声淹没；步长太大（比如 5%），超出弹性范围，曲线弯曲，拟合出的弹性常数偏小。常见做法是在 ±2% 范围内取 5 到 7 个点，步长 0.5% 或 1%。
每个应变点都要重新跑静态计算，K 点网格和截断能要和结构优化时一致。如果 relax 用的 8×8×8，静态也得用 8×8×8，否则能量和应力不可比。收敛标准EDIFF建议设到 1E-6 甚至 1E-7，应力对电子步收敛更敏感。
注意：QE 的
ecutwfc和ecutrho在应变计算中要保持不变。有人为了省时间在静态计算里降截断能，结果应力张量整体偏移，弹性常数全错。
4. 拉伸与剪切脚本实操：从改 POSCAR 到拟合弹性常数
4.1 拉伸计算：改哪个晶格参数、怎么改
拉伸计算的核心是构造一系列形变后的结构。以金刚石为例，原始 POSCAR 的晶格矢量是三行，缩放系数在第一行。沿 z 方向拉伸，就是把第三行矢量乘以 (1+ε)，其他两行不变。脚本里如果直接改缩放系数，那是各向同性拉伸，只能算体弹模量。
# 沿指定方向施加单轴应变，生成一系列 POSCAR import numpy as np def generate_strained_poscar(base_lattice, strain_list, direction=2): """ base_lattice: 3x3 晶格矢量矩阵，每行一个矢量 strain_list: 应变值列表，如 [-0.02, -0.01, 0, 0.01, 0.02] direction: 0/1/2 对应 x/y/z 方向 返回：每个应变对应的晶格矩阵列表 """ strained = [] for eps in strain_list: new_lat = base_lattice.copy() # 只改指定方向的晶格矢量长度 new_lat[direction] = base_lattice[direction] * (1 + eps) strained.append(new_lat) return strained # 读取原始 POSCAR 的晶格部分 def read_poscar_lattice(filename): with open(filename) as f: lines = f.readlines() scale = float(lines[1].strip()) lattice = np.array([[float(x) for x in lines[i].split()] for i in range(2, 5)]) * scale return lattice这段代码只做了单轴拉伸。要算完整的弹性常数矩阵，需要分别沿 x、y、z 拉伸，再施加 xy、yz、zx 三个剪切模式。每个模式对应一列弹性常数。脚本包里的tensile_calculation_*应该覆盖了单轴拉伸，shear_calculation_*覆盖剪切。
参数说明：strain_list的范围建议 ±0.02，点数 5 到 7 个。direction参数在单轴拉伸时只改一个方向，但实际计算中其他方向可能因为泊松效应产生应力，所以 ISIF 要设成 2（允许原子弛豫但不改晶胞），或者设成 4（允许改晶胞形状但体积不变）——具体看你要算的是哪个弹性常数。
4.2 剪切计算：形变矩阵的构造
剪切比拉伸麻烦，因为要改的是晶格矢量之间的夹角。以 xy 剪切为例，形变矩阵是单位矩阵加上一个非对角元：
# 构造剪切形变矩阵并应用到晶格 def apply_shear(lattice, eps, plane='xy'): """对晶格施加剪切应变，eps 是剪切量""" deform = np.eye(3) if plane == 'xy': deform[0, 1] = eps # x 方向矢量在 y 方向的分量 elif plane == 'yz': deform[1, 2] = eps elif plane == 'zx': deform[2, 0] = eps # 形变后的晶格 = 原始晶格 @ 形变矩阵的转置 return lattice @ deform.T剪切应变下，晶胞体积基本不变，但对称性降低。VASP 的 ISIF 要设成 2，让原子在固定晶胞内弛豫。QE 里calculation = 'scf'配合tprnstr = .true.即可。
剪切计算最容易翻车的地方是应变方向。xy 剪切和 yx 剪切在弹性常数矩阵里是同一个分量，但形变矩阵的写法不同。如果脚本里deform[0,1]和deform[1,0]混用，算出的 C44 和 C55 会对调。建议先跑一个已知材料验证，比如金刚石的 C44 约 580 GPa，算出来差太多就是方向搞错了。
4.3 应力提取与弹性常数拟合
拿到一系列应变和对应的应力后，拟合就是线性回归。弹性常数 Cij = dσi / dεj，在弹性范围内应力应变是线性的。用 numpy 的 polyfit 一次多项式即可。
# 从应力-应变数据拟合弹性常数 import numpy as np def fit_elastic_constant(strain_list, stress_list): """ strain_list: 应变值列表 stress_list: 对应应力分量列表（GPa） 返回：弹性常数（GPa）和拟合优度 R^2 """ coeffs = np.polyfit(strain_list, stress_list, 1) slope = coeffs[0] # 弹性常数 # 计算 R^2 p = np.poly1d(coeffs) yhat = p(strain_list) ybar = np.mean(stress_list) ss_res = np.sum((stress_list - yhat) ** 2) ss_tot = np.sum((stress_list - ybar) ** 2) r2 = 1 - ss_res / ss_tot return slope, r2 # 示例：单轴拉伸沿 z 方向，取 σ_zz 对 ε_zz 的斜率 strains = [-0.02, -0.01, 0.0, 0.01, 0.02] stresses = [-12.5, -6.3, 0.1, 6.2, 12.4] # 单位 GPa C33, r2 = fit_elastic_constant(strains, stresses) print(f"C33 = {C33:.1f} GPa, R^2 = {r2:.4f}")R² 低于 0.99 就要检查：应变范围是否太大、某个应变点的计算是否没收敛、应力提取是否取错了行。金刚石这类高对称材料，线性应该非常好，R² 接近 1。
提示：拟合时截距应该接近零。如果截距明显偏离零，说明零应变点的结构没有完全弛豫，或者应力提取有系统误差。
5. 避坑与排查：应力应变计算里最容易翻车的五件事
5.1 现象：应力张量全是零
原因：VASP 的 INCAR 里没设ISIF = 2或IBRION = -1，或者 QE 的tprnstr没打开。VASP 默认 ISIF=2 会算应力，但如果 NSW=0 且 IBRION=-1，应力应该正常输出。QE 的tprnstr默认是.false.，必须显式设成.true.。
解决：检查 INCAR 和 QE 输入文件，确认应力输出开关打开。VASP 还可以看 OUTCAR 里有没有FORCES acting on ions和Stress tensor段落。
5.2 现象：不同应变点的能量不连续
原因：K 点网格或截断能在不同应变点之间变了。有人为了省时间，在小应变时用低截断，大应变时用高截断，导致能量基准不一致。
解决：所有应变点用完全相同的计算参数。写脚本时把 KPOINTS 和 INCAR 模板固定，只改 POSCAR 或 CELL_PARAMETERS。
5.3 现象：拟合出的弹性常数偏小
原因：应变范围太大，超出了线性弹性区。金刚石的弹性线性范围大概在 ±1% 以内，加到 ±3% 曲线就弯了。
解决：把应变范围缩到 ±1%，或者用二次多项式拟合后取一次项系数。更稳妥的做法是先跑一个大范围看曲线拐点，再在拐点内取点。
5.4 现象：QE 和 VASP 算出的弹性常数差很多
原因：两者的赝势不同、交换关联泛函可能不同、应力单位换算有误。VASP 的 kB 转 GPa 是乘 0.1，QE 的 Ry/bohr³ 转 GPa 要乘 14710.5。
解决：先统一泛函（都用 PBE）和赝势类型（都用 PAW 或都用超软），再核对单位换算。如果还差，检查 QE 的ecutwfc是否足够，QE 对截断能比 VASP 敏感。
5.5 现象：脚本跑完没有输出文件
原因：VASP 或 QE 的可执行文件路径没设对，或者 mpirun 的进程数和 K 点并行不匹配。脚本里如果硬编码了vasp_std或pw.x，环境变量 PATH 里没有就会静默失败。
解决：在脚本里加错误检查，跑完检查 OUTCAR 或 output 文件是否存在且非空。常见做法是用subprocess.run的check=True让异常抛出，而不是默默跳过。
6. 进阶技巧：用 plotcheck 脚本做快速验证与批量扫描
plotcheck版本的脚本价值在于边算边看。我一般会先用它跑一个应变点，确认应力提取和绘图链路通了，再切到withoutplot批量跑。批量跑的时候，把应变列表写成循环，每个应变生成一个目录，跑完统一收集数据。
# 批量扫描应变并收集应力，适合 withoutplot 模式 import os import subprocess import numpy as np def batch_strain_scan(base_dir, strain_list, direction=2): """在每个应变下生成结构、跑 VASP、收集应力""" results = [] for eps in strain_list: work_dir = os.path.join(base_dir, f"strain_{eps:+.4f}") os.makedirs(work_dir, exist_ok=True) # 生成 POSCAR（省略具体实现，参考 4.1 节） # 复制 INCAR、KPOINTS、POTCAR 到 work_dir # 跑 VASP subprocess.run("mpirun -np 4 vasp_std > vasp.log", shell=True, cwd=work_dir, check=True) # 提取应力 stress = read_stress_outcar(os.path.join(work_dir, "OUTCAR")) results.append((eps, stress[2, 2])) # 取 σ_zz return results # 跑完直接拟合 strains, stresses = zip(*batch_strain_scan("./scan", [-0.01, -0.005, 0, 0.005, 0.01])) C33, r2 = fit_elastic_constant(list(strains), list(stresses)) print(f"C33 = {C33:.1f} GPa, R^2 = {r2:.4f}")这个批量脚本的关键是每个应变独立目录，避免文件覆盖。check=True保证 VASP 报错时脚本停下来，而不是继续跑下一个应变点。收集完数据直接拟合，R² 不合格就回去查哪个应变点出了问题。
验证方法上，我习惯拿金刚石或硅做基准。金刚石的 C11 约 1076 GPa、C12 约 125 GPa、C44 约 577 GPa，算出来在这个量级附近就说明流程通了。如果差一个数量级，多半是单位换算错了；如果差 20% 到 30%，检查赝势和截断能。
注意：QE 的应力输出在 XML 文件里更规整，
pw.x的 text output 有时会截断。如果脚本读 text output 不稳定，改用style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;" />
