简介面向材料科学领域研究者和计算模拟爱好者提供一套基于VASP与Quantum Espresso计算应力应变关系的脚本资源。内容围绕密度泛函理论模拟后的数据后处理展开帮助用户从结构优化与应变计算的输出文件中提取应力与应变数据绘制并拟合应力应变曲线进而获得弹性模量、泊松比等关键力学参数可显著减少手工处理大量输出数据的重复劳动。压缩包共16个文件主体为8个Python脚本分别覆盖拉伸与剪切两类变形工况并针对VASP和QE两条计算路线分别编写便于按实际使用的软件选取另含4个计算输入文件、POSCAR及旋转结构文件以及说明文档整体仅30KB轻量易部署。脚本内置了“仅计算”和“计算并绘图”两种运行模式满足不同操作习惯配合说明文档和示例输入既能引导初学者快速掌握流程也方便有一定基础的用户修改参数迁移到其他晶体结构或力学性质计算中。目前已有932人学习下载适合材料、物理、化学等方向需要开展第一性原理力学模拟的科研人员与研究生参考使用。1. 为什么用 VASP 和 QE 算应力-应变关系从胡克定律到弹性矩阵材料计算里常见的一个需求是晶体结构已经拿到了下一步想知道这材料“硬不硬、脆不脆”也就是杨氏模量、剪切模量、泊松比。这些量背后全是应力-应变关系曲线——一组从多个小应变计算里拟合出来的数据。VASP 和 QE 是目前第一性原理领域最常用的两套 DFT 代码前者收敛稳定、流程成熟后者开源透明、输入文件好改两者配合 Python 的 numpy/scipy 做应变生成、输出解析和曲线拟合就能在本地把弹性常数计算流程完整跑起来。适用人群是材料计算方向的学生或工程师手里已经有 POSCAR 或 QE 输入文件缺的是一套成体系的应力和应变计算脚本以及能交叉验证的做法。2. 胡克定律、应变模式和拟合方法应力-应变计算背后的物理2.1 胡克定律和 Voigt 排序把 6 个应力分量读成一条直线弹性区内应力与应变满足广义胡克定律 σ_ij C_ijkl ε_kl其中 C 是四阶刚度张量。四阶张量共 81 个分量看着复杂但真实材料里由对称性约束后独立分量很少。工程和 DFT 后处理都不直接操作四阶张量而是用 Voigt 排序压成 6×6 矩阵下标 1xx、2yy、3zz、4yz、5xz、6xy。例如 C11 对应 x 方向单轴拉伸的弹性模量C44 对应剪切方向。晶体对称性决定独立常数的数量立方晶系只有 3 个C11、C12、C44六方晶系 5 个C11、C12、C13、C33、C44正交晶系 9 个。QE 和 VASP 输出的应力张量分量顺序基本都是 (xx, yy, zz, yz, xz, xy)解析时对齐这个顺序是 Python 脚本里最容易出差异的一步。获取弹性常数的常见做法不是直接拟合“应力-应变直线”而是用能量-应变关系。把应变 ε 代入总能公式展开到二阶得到 ΔE(ε) V0/2 × Σ_i Σ_j C_ij ε_i ε_j加高阶项。为什么用能量而不是应力两条路线信息等价但第一性原理的总能量是变分极值收敛后能量噪声小应力张量虽然也打印在 OUTCAR 或 pwo 文件里但数值噪声相对大原子位置不完美收敛时应力会明显漂移。能量-应变法对新人更友好拟合残差能直接反映数据质量这也是后面三章统一采用的主路线。2.2 应变模式的选择三类典型模式覆盖立方晶系直接拿任意 3×3 应变张量改晶格是最容易翻车的做法。一旦打破了晶体对称性不等价原子数量增加k 点归约的对称性降级能量很容易出现无规律的毛刺。常见做法是先选定体系里“对称等价的应变模式”再通过几个模式叠加反解出独立的 C_ij。立方晶系习惯上用三个模式见下表应变模式 (εxx, εyy, εzz)拟合关系 ΔE/V0可解出的组合(δ, 0, 0)(1/2) C11 δ²C11(δ, δ, 0)(C11 C12) δ²C11 C12(δ, δ, -2δ)3(C11 - C12) δ²C11 - C12第一个模式单独给 C11第二个和第三个联立给 C12。C44 则需要剪切型应变模式比如把应变张量设成分量 εxy 非零的组合此时保证晶格行列式的一阶变化为零对应纯剪切。六方和正交晶系就是把上表按轴向和剪切方向扩展解析式更长但原理一致。2.3 应变步长和拟合阶数先看信噪比再定参数应变范围习惯取 -2% 到 2%。这个跨度既不会进入非线性区又让能量差远大于数值噪声。单点 SCF 能量在正常收敛下的噪声约 1 meV应变 1% 时典型能量变化几十到几百 meV信噪比足够。步长建议取 5~9 个点均匀分布在区间内不要只取三个点——三点二次拟合对线性项和噪声完全没有防御能力。用 numpy 的 polyfit 做二次拟合即可import numpy as np # deltas: 应变量, 如 [-0.02, -0.015, -0.01, -0.005, 0.0, 0.005, 0.01, 0.015, 0.02] # energies: 对应总能, 单位 eV, 从 VASP 或 QE 输出中解析 deltas np.array([-0.02, -0.015, -0.01, -0.005, 0.0, 0.005, 0.01, 0.015, 0.02]) energies np.array([...]) # 由实际计算填充 volume_A3 159.99 # 应变前原胞体积, 单位 ų coef np.polyfit(deltas, energies, 2) c11_gpa 2.0 * coef[0] / volume_A3 * 160.217 # eV/ų 转 GPa residual energies - np.polyval(coef, deltas) rms_mev np.sqrt(np.mean(residual**2)) * 1000.0 print(fC11 {c11_gpa:.1f} GPa, 拟合残差 RMS {rms_mev:.2f} meV)残差 RMS 超过 2 meV 时先回头查收敛精度和 k 点密度而不是加三次方项硬凑。二次多项式的线性项反映零应变的微小应力残留这是正常的保留即可。3. VASP修改 POSCAR 施加应变、INCAR 参数和 OUTCAR 应力解析3.1 用 Python 生成应变后的 POSCAR核心是格矢乘以 (Iε)VASP 的 POSCAR 前 3 行是晶格向量。施加应变时把应变矩阵 ε 加到单位阵 I 上得到变形梯度 F I ε新格矢 旧格矢 × F。原子分数坐标保持不动笛卡尔坐标自动跟随晶格一起变形。写一个通用 Python 函数来生成整个应变系列import numpy as np def read_poscar(pathPOSCAR): 读取 POSCAR 的缩放因子和前 3 行格矢, 返回实际格矢矩阵 with open(path, errorsignore) as f: lines f.readlines() scale float(lines[1]) lattice np.array([[float(x) for x in line.split()] for line in lines[2:5]]) return lattice * scale def strain_lattice(lattice, eps, modeuniaxial_x): 给格矢施加应变, mode 支持 uniaxial_x / biaxial_xy / tri_axial_iso I np.eye(3) S np.zeros((3, 3)) if mode uniaxial_x: S[0, 0] eps elif mode biaxial_xy: S[0, 0] S[1, 1] eps elif mode tri_axial_iso: S[0, 0] S[1, 1] eps S[2, 2] -2.0 * eps # 剪切模式在这里按需扩展 return lattice (I S) def write_poscar(path, lattice, rest_lines): 写出带应变的新 POSCAR, rest_lines 是从第 6 行起的原子信息 with open(path, w) as f: f.write(Strained by python\n) f.write(1.0\n) for row in lattice: f.write(f{row[0]:.10f} {row[1]:.10f} {row[2]:.10f}\n) f.write(rest_lines) if __name__ __main__: old read_poscar(POSCAR) with open(POSCAR, errorsignore) as f: rest \n.join(f.readlines()[5:]) for delta in np.linspace(-0.02, 0.02, 9): new_lat strain_lattice(old, delta, uniaxial_x) write_poscar(fPOSCAR_{delta:.3f}, new_lat, rest)注意两点。第一数值格式用%.10f而不是默认的简洁格式避免格矢和分数坐标之间出现微小不一致这对高度收敛的计算会产生 1% 量级的弹性常数差异。第二rest_lines 保存了 POSCAR 从第 6 行开始的元素名、离子数、坐标类型和分数坐标这些内容在应变前后完全不变。生成的 POSCAR 文件名带上应变量后续解析时可以直接对应。3.2 INCAR 参数设置单点计算配小的 EDIFF应变扫描每步只做一次单点自洽不做原子弛豫。INCAR 最小集合如下PREC Accurate ENCUT 1.3 * ENMAX EDIFF 1E-6 ISMEAR -5 SIGMA 0.1 IBRION -1 NSW 0 ISIF 3 LWAVE .FALSE. LCHARG .FALSE.几个参数的选择理由要讲清楚。PREC Accurate 和 ENCUT 1.3 倍 ENMAX 是为了让 PAW 投影在应变后依然稳定截断不足会让应力张量出现明显漂移。EDIFF 取 1E-6 而非默认的 1E-4保证总能噪声压在亚 meV 量级默认值在小应变扫描下会让能量差淹没在噪声里。ISMEAR -5 适用于绝缘体和半导体金属体系换成 ISMEAR 1 并配合 SIGMA 0.1再用 SIGMA → 0 外推总能量否则金属半占据会在应变下产生虚假能量漂移。ISIF 3 在这里是让 VASP 输出完整应力张量IBRION -1 和 NSW 0 明确禁止离子移动。LWAVE 和 LCHARG 关闭后节省磁盘空间单点计算不需要后续续算。3.3 解析 OUTCAR 的应力kB 和 GPa 的换算别弄反VASP 的 OUTCAR 里总应力张量出现在文本FORCES: max atom, stress之后以 3×3 矩阵形式打印单位是 kBar。1 kBar 0.1 GPa。OUTCAR 在离子步里会多次出现这个块单点计算只有一次但保险起见解析时取最后一次出现的块import numpy as np def read_vasp_stress_gpa(outcarOUTCAR): 提取最后一个总应力块, 返回 6 分量 (xx,yy,zz,yz,xz,xy), 单位 GPa with open(outcar, errorsignore) as f: lines f.readlines() stress None for i, line in enumerate(lines): if FORCES: max atom, stress in line: try: block lines[i1:i8] matrix [] for bl in block: parts bl.split() if len(parts) 3: matrix.append([float(x) for x in parts[:3]]) if len(matrix) 3: m np.array(matrix[:3]) stress np.array([m[0,0], m[1,1], m[2,2], m[1,2], m[0,2], m[0,1]]) * 0.1 except (ValueError, IndexError): continue if stress is None: raise RuntimeError(没读到应力, 检查 OUTCAR 是否正常算完) return stress def read_vasp_energy(outcarOUTCAR): 取最后一个 free energy TOTEN 行的总能, 单位 eV with open(outcar, errorsignore) as f: lines f.readlines() energies [float(l.split()[4]) for l in lines if l.startswith( free energy TOTEN)] if not energies: raise RuntimeError(OUTCAR 里没有 TOTEN 行) return energies[-1]这里把 3×3 矩阵按行主序还原成 6 分量对角元是 xx/yy/zzm[1,2]、m[0,2]、m[0,1] 对应 yz/xz/xy。乘 0.1 完成 kBar 到 GPa 的换算。能量读取取 TOTEN 行的第 5 列这是 OUTCAR 里的自由能也就是常规单点能量。批量跑多个应变时可以写个简单 shell 循环for d in -0.020 -0.015 -0.010 -0.005 0.000 0.005 0.010 0.015 0.020; do cp POSCAR_${d} POSCAR mpirun -np 8 vasp_std run_${d}.log python3 -c import parse_vasp; print($d, parse_vasp.read_vasp_energy()) done每个应变目录独立运行避免 OUTCAR 互相覆盖。4. QEQuantum ESPRESSO用 CELL_PARAMETERS 施加应变和应力读取4.1 pw.x 输入文件里怎么改晶格QE 施加应变通常借助 ibrav 0 时的 CELL_PARAMETERS 段。把笛卡尔格矢按 3×3 矩阵写在里面单位在括号里指定。一个典型的硅输入文件如下CONTROL calculation scf prefix si pseudo_dir ./ outdir ./tmp / SYSTEM ibrav 0 nat 2 ntyp 1 ecutwfc 40 ecutrho 320 occupations fixed / ATOMIC_SPECIES Si 28.085 Si.pz-vbc.UPF CELL_PARAMETERS (angstrom) 5.43000000 0.00000000 0.00000000 0.00000000 5.43000000 0.00000000 0.00000000 0.00000000 5.43000000 ATOMIC_POSITIONS (crystal) Si 0.00 0.00 0.00 Si 0.25 0.25 0.25 K_POINTS (automatic) 8 8 8 0 0 0注意原子坐标用 crystal 类型。应变后格矢变化但分数坐标不变原子在实空间的位置自动跟随晶格变形这和 VASP 里保留分数坐标是同一个道理。生成一系列应变输入时只需要替换 CELL_PARAMETERS 的三行ATOMIC_POSITIONS 整块不动。这里额外提一句用 ibrav 0 时 celldm(1) 最好设成 1.0让 CELL_PARAMETERS 的单位由括号里的angstrom明确控制避免和默认 Bohr 单位混淆。4.2 从 pw.x 输出里读应力和总能抓最后一组数字pw.x 在自洽计算结束后打印total stress后面按 Ry/bohr³、kbar、GPa 三组顺序给出 6 个分量。稳妥的做法是取标识行后面一行里的全部数字分成三组各 6 个直接取第三组作为 GPa 值不要自己再从 Ry/bohr³ 手动换算省一次单位坑import re import numpy as np def _isfloat(s): try: float(s) return True except ValueError: return False def read_qe_stress(pwosi.pwo): 从 pw.x 输出摘应力, 返回 6 分量 (xx,yy,zz,yz,xz,xy), 单位 GPa with open(pwo, errorsignore) as f: txt f.read() idx txt.rfind(total stress) if idx -1: raise RuntimeError(pwo 输出里没有 total stress) seg txt[idx: idx 200].split() nums [float(x) for x in seg if _isfloat(x)] if len(nums) 18: raise RuntimeError(stress 行数字不足, 检查输出完整性) gpa nums[12:18] return np.array(gpa) def read_qe_energy(pwosi.pwo): 取最后一个 ! 行的总能, 单位换算成 eV with open(pwo, errorsignore) as f: txt f.read() matches re.findall(r!\stotal energy\s*\s*([-\d.]), txt) if not matches: raise RuntimeError(没找到总能量行) return float(matches[-1]) * 13.6057 # Ry - eV符号约定方面pw.x 打印的应力数值正值表示体系倾向于压缩对应“外部压力为正”的习惯和材料力学里“拉应力为正”相反。用能量-应变法拟合弹性常数时不涉及应力符号因为能量和 ε² 挂钩但如果想用直接应力 σ Cε 做交叉核对就需要把 QE 打印值取反再把符号统一成拉正压负。4.3 两个代码的应力和能量读取对照把 VASP 和 QE 的输出项放在同一张表里写解析脚本时不容易搞混数据来源VASPQE总能量行free energy TOTEN! total energy应力标记FORCES: max atom, stresstotal stress应力原始单位kBarRy/bohr³ kbar GPa晶格修改位置POSCAR 前 3 行CELL_PARAMETERS原子坐标类型保留分数坐标crystal 坐标同一结构的同一种应变模式VASP 和 QE 拟合出的 C11 应在 5% 以内符合这本身就是交叉验证的第一步。偏差超过这个范围时回第 5 章查单位、符号和收敛。5. 应力-应变计算的 5 个常见坑踩坑记录5.1 应变步长太小或太大能量差被噪声吃掉现象拟合出的弹性常数随步长明显漂移和文献值对不上。 原因步长小于 0.3% 时应变能量差只有几个 meV与 SCF 数值噪声和 k 点不完全收敛处于同一量级步长大于 3% 时三次以上非线性项开始混进来二次拟合系数被系统性高估或低估。 解决先在一个应变模式下用 0.5%、1%、2% 三组步长试算二次项系数稳定在 1% 以内再用。默认取 -2% 到 2% 共 9 个点是经验上最稳的组合。5.2 单位换算出错kB、GPa、eV/ų 三个坑现象VASP 算的 C11 是 103.2 GPaQE 算出 1032刚好差 10 倍或者能量法出来的模量是 1.6压根不像 GPa 量级。 原因OUTCAR 里应力默认 kBar1 kBar 0.1 GPa而能量-应变拟合直接用 eV 和 ų 时会得到 eV/ų1 eV/ų 约等于 160.217 GPa。很多人处理好应力单位却忘了能量也需要换算或者反过来。 解决在脚本顶部把换算系数写成常量比如KBAR_TO_GPA 0.1、EV_A3_TO_GPA 160.217所有解析函数统一输出 GPa最终汇总表里只出现一种单位。5.3 剪切应变破坏了晶格对称性现象加剪切应变后原本一个原胞能算的体系突然变慢或能量-应变曲线出现不连续毛刺。 原因任意应变矩阵把格矢扭转后晶格对称性下降不等价原子变多k 点归约失效计算量成倍上涨。 解决选应变模式前先确认晶系点群。剪切应变要用对称等价的组合而不是裸的工程剪切或者直接改用原胞而非超胞。跑正式扫描前先做一次零应变单点计算确认应变的 KPOINTS 设置和 ISMEAR 仍合理。5.4 EDIFF 太松或 ISMEAR 不合适能量曲线出现毛刺现象ΔE 对 δ 的散点明显偏离抛物线拟合残差大于 2 meV。 原因EDIFF 1E-4 时总能误差在 meV 量级小应变能量差被噪声盖住金属体系用 ISMEAR -5 强制绝缘体处理半占据误差随应变乱跳。 解决EDIFF 调到 1E-6 甚至 1E-7金属用 ISMEAR 1 加适当 SIGMA再做 SIGMA → 0 外推。同一个应变集合的 k 点和截断设置保持一致中途不要改动。5.5 符号约定没对齐直接应力法会给出反向弹性常数现象VASP 和 QE 算出的 C11 一个为正一个为负或者和实验符号相反。 原因VASP 和 QE 打印应力张量的正负方向约定不同差一个整体负号。能量-应变法对正负应变对称不受影响但一旦想用直接应力 σ Cε 做交叉核对符号就出问题。 解决统一以“拉应力为正”。读取 QE 应力后取反VASP 的 OUTCAR 应力按同样方向约定处理后再放回胡克定律。符号归一的方法是在解析函数里加一行注释标明当前约定避免两个月后回来看脚本时重新踩一遍。6. 自动应变扫描脚本VASP 和 QE 交叉验证一条龙6.1 用一份 Python 脚本把两套流程串起来把前三章的步骤合并输入是基础 POSCAR 和基础 QE 输入输出是弹性常数对比表。脚本骨架如下import numpy as np CONFIG { deltas: [-0.02, -0.015, -0.01, -0.005, 0.0, 0.005, 0.01, 0.015, 0.02], mode: uniaxial_x, volume_A3: 159.99, } def fit_c11(deltas, energies, volume_A3): coef np.polyfit(deltas, energies, 2) return 2.0 * coef[0] / volume_A3 * 160.217 # GPa # 主流程: # 1) 对每个 delta 生成 POSCAR / QE 输入 # 2) 分别调用 vasp_std 和 pw.x 运行 SCF # 3) 调用 read_vasp_energy / read_qe_energy 收集能量 # 4) fit_c11 得到 VASP 和 QE 各自的 C11, 并排打印实际运行时每个应变最好独立工作目录避免 OUTCAR 和 pwo 互相覆盖。脚本里保留一份 CSV 输出把 delta、VASP 能量、QE 能量、两种代码拟合的 C11 都记录下来。6.2 收尾的自检清单验证计算是否做对的几个检查点零应变能量与文献值对比偏差大于 0.1 eV/atom 就先查赝势和截断正负应变点的能量曲线对称性抛物线带一点三次项正常但严重偏离说明应变模式选错两种代码结果差超过 5% 时先回 5.2 查单位再回 5.4 查收敛。我的习惯是每个项目把应变扫描脚本和结果 CSV 存在同一目录算完顺手把零应变能量、拟合残差、两种代码对比写进一个 README。这套流程跑过不同晶系后会发现真正影响结果的大多是细节单位、符号、拟合阶数。这三件事确认了VASP 和 QE 算的弹性常数基本都能互相印证。希望这套方法帮到你少走几步弯路。本文还有配套的精品资源点击获取