VASP与QE弹性常数计算：Python统一处理应力应变与拟合

简介这份资源面向材料科学计算方向的研究生与科研人员聚焦如何用Python串联VASP与Quantum Espresso两大DFT软件完成应力—应变关系的提取、处理与可视化适合已具备第一性原理计算基础、希望系统掌握力学性质后处理流程的中高级用户。压缩包共16个文件约30KB以8个Python脚本为核心配合4个输入文件及POSCAR、POSCAR_rota等结构文件另有README说明与gitattributes配置脚本按拉伸与剪切两类工况分别提供VASP和QE版本并区分是否绘图便于按需选用。目前已有931人学习下载。读者可从中获得读取输出文件、提取应力应变数据、绘制曲线并进一步拟合弹性模量与泊松比的完整脚本范例同时借助示例输入与结构文件快速复现计算流程理解两种软件在应变模拟中的差异与衔接方式。1. 从一条应力应变曲线说起VASP 与 QE 的弹性常数怎么算、Python 怎么串起来做第一性原理计算的人迟早会碰到同一个需求给一个晶体施加一系列微小形变算出对应的应力再拟合出弹性常数。VASP 和 Quantum ESPRESSOQE都能做这件事但两者的输入输出格式、应力单位、形变施加方式完全不同手工来回倒数据非常容易出错。标题里的「使用 VASP 和 QE 计算应力和应变关系」说的就是这条链路用两种主流 DFT 引擎分别跑形变-应力数据再用 Python 统一后处理拟合出 C11、C12、C44 这类弹性常数或者直接得到应力-应变斜率。这套流程适合做力学性质、相稳定性、材料硬度的从业者也适合刚接触弹性常数计算、想找一条能复现路径的新手。核心难点不在 DFT 本身而在三件事形变怎么加才不破坏对称性、应力单位怎么统一、拟合时哪些点该丢。下面按「先立住原理再动手复现最后讲坑」的顺序拆开讲。2. 应力应变与弹性常数先把物理量和单位对齐2.1 弹性常数到底在算什么弹性常数描述的是晶体在小形变下的线性响应。广义胡克定律写成矩阵形式就是 σ_i C_ij · ε_j其中 σ 是应力Voigt 记法 6 分量ε 是应变C_ij 是 6×6 的弹性刚度矩阵。对立方晶系独立分量只有 C11、C12、C44 三个六方晶系有 5 个正交晶系有 9 个。你要算几个独立分量取决于晶体对称性而不是想算几个就算几个。实际计算里我们不去直接解这个矩阵而是对晶体施加一组特定的应变模式让某个应变分量非零、其余为零然后读取 DFT 输出的应力张量。比如施加 ε_1 δ其余为 0那么 σ_1 C11·δσ_2 C12·δσ_3 C12·δ。这样一次形变就能同时拿到 C11 和 C12。C44 则需要施加剪切应变 ε_4即 ε_yz。提示应变幅度一般取 ±0.001 到 ±0.01 之间。太小则应力数值噪声占主导太大则进入非线性区拟合斜率会偏。2.2 VASP 和 QE 的应力输出差异VASP 的应力在 OUTCAR 里单位是 kB千巴1 kB 0.1 GPa。QE 的应力在输出文件里单位是 kbar同样 1 kbar 0.1 GPa但 QE 默认输出的是负的应力即压强符号约定和 VASP 相反。这是最容易翻车的地方如果你把两边的数直接拼在一起拟合斜率符号会乱。项目VASPQE应力位置OUTCAR 中 in kB 行pwscf.out 中 total stress 行单位kBkbar符号约定正为拉伸正为压缩需取反形变施加修改 POSCAR 晶格修改 CELL_PARAMETERS推荐截断能1.3×ENMAXecutwfc 收敛测试统一到 GPa 并修正符号后两边数据才能放进同一个拟合脚本。Python 在这里的角色就是解析、单位换算、符号修正、拟合一条龙。2.3 用 Python 解析 VASP 应力最小可跑脚本先看 VASP 侧。OUTCAR 里应力行格式固定用正则抓取即可。下面这段脚本读取 OUTCAR提取最后一步的应力张量转成 GPa。import re import numpy as np def read_vasp_stress(outcar_path): 从 VASP OUTCAR 读取最后一步应力返回 3x3 数组单位 GPa with open(outcar_path, r) as f: lines f.readlines() # 从后往前找取最后一次出现的应力块 stress [] for i in range(len(lines) - 1, -1, -1): if in kB in lines[i]: # 应力块共 3 行每行 6 个数取前 3 个为 xx yy zz for j in range(3): nums re.findall(r[-]?\d\.\d, lines[i j]) stress.append([float(nums[0]), float(nums[1]), float(nums[2])]) break stress np.array(stress) # kB - GPa: 1 kB 0.1 GPa return stress * 0.1 if __name__ __main__: s read_vasp_stress(OUTCAR) print(VASP stress (GPa):\n, s)逻辑说明正则[-]?\d\.\d匹配带符号浮点数避免把行首的 in kB 误抓。从后往前找保证拿到收敛后的最后一步。参数上stress * 0.1是单位换算如果你用的 VASP 版本输出单位不同需要核对 OUTCAR 表头。返回的是 3×3 矩阵对角元是正应力非对角元是剪切应力。2.4 用 Python 解析 QE 应力并修正符号QE 的输出格式和 VASP 不同应力块在 total stress 之后单位 kbar且符号需要取反。import re import numpy as np def read_qe_stress(out_path): 从 QE pwscf.out 读取总应力返回 3x3 数组单位 GPa with open(out_path, r) as f: text f.read() # 定位 total stress 块 match re.search(rtotal\sstress\s\(kbar\)\s*\n\s*([-\d\.\s])\n, text) if not match: raise ValueError(未找到 total stress 块) nums [float(x) for x in match.group(1).split()] # QE 输出 6 个分量: xx yy zz xy xz yz s np.array(nums[:6]) # kbar - GPa: 1 kbar 0.1 GPa且 QE 符号为压缩正需取反 s_gpa -s * 0.1 # 组装成 3x3 mat np.array([[s_gpa[0], s_gpa[3], s_gpa[4]], [s_gpa[3], s_gpa[1], s_gpa[5]], [s_gpa[4], s_gpa[5], s_gpa[2]]]) return mat if __name__ __main__: s read_qe_stress(pwscf.out) print(QE stress (GPa):\n, s)逻辑说明QE 的应力块是 6 个数一行顺序为 xx yy zz xy xz yz需要映射到 3×3 矩阵。-s * 0.1同时完成符号修正和单位换算。参数上如果你的 QE 版本输出的是 total stress 带多个空格正则里的\s能兼容。注意 QE 的剪切分量顺序和 VASP 不同VASP 是 xx yy zz xy yz zx映射时要小心。3. 形变怎么加VASP 和 QE 的应变施加实操3.1 应变模式与晶格矩阵的对应关系施加应变本质上是修改晶格矢量。给定原始晶格矩阵 A3×3每行是一个晶格矢量施加应变 ε 后新晶格 A A · (I ε)其中 I 是单位矩阵ε 是对称应变张量。对立方晶系施加 ε_1 δ 时ε 矩阵只有 (1,1) 元为 δ其余为 0。施加剪切 ε_4 时ε 矩阵的 (2,3) 和 (3,2) 元为 δ/2。这一步用 Python 生成新晶格最方便避免手算出错。import numpy as np def apply_strain(lattice, strain_vec): lattice: 3x3 原始晶格; strain_vec: 6 分量 Voigt 应变 返回形变后的晶格矩阵 exx, eyy, ezz, gyz, gxz, gxy strain_vec eps np.array([[exx, gxy/2, gxz/2], [gxy/2, eyy, gyz/2], [gxz/2, gyz/2, ezz]]) return lattice (np.eye(3) eps) # 示例立方晶格施加 0.005 的 xx 应变 A np.array([[5.43, 0, 0], [0, 5.43, 0], [0, 0, 5.43]]) A_new apply_strain(A, [0.005, 0, 0, 0, 0, 0]) print(A_new)逻辑说明Voigt 记法的剪切分量 γ 对应张量里的 ε_ij γ/2i≠j这是工程应变和张量应变的区别搞错会让 C44 偏大一倍。参数上strain_vec的 6 个分量顺序是 xx yy zz yz xz xy和 QE 输出顺序一致。生成的 A_new 直接替换 POSCAR 或 CELL_PARAMETERS 即可。3.2 VASP 批量形变计算的组织方式实际做的时候不会手动改 POSCAR 几十次。常见做法是写一个 Python 脚本对每个应变幅度生成一个子目录复制 INCAR、KPOINTS、POTCAR只改 POSCAR 的晶格行然后批量提交。#!/bin/bash # 批量生成 VASP 形变目录并提交 for strain in -0.01 -0.005 0.0 0.005 0.01; do dirstrain_${strain} mkdir -p $dir cp INCAR KPOINTS POTCAR $dir/ # 用 Python 生成该应变下的 POSCAR python3 gen_poscar.py --strain $strain --out $dir/POSCAR cd $dir mpirun -np 16 vasp_std vasp.log cd .. done wait逻辑说明每个应变一个目录互不干扰。gen_poscar.py内部调用上面的apply_strain读取原始 POSCAR 的晶格和原子坐标只改晶格行原子分数坐标保持不变因为应变是均匀施加在晶格上的。参数上-np 16按你的核数调整应变点建议至少 5 个含 0覆盖正负两侧以便线性拟合。3.3 QE 的 CELL_PARAMETERS 修改与输入模板QE 的形变通过CELL_PARAMETERS卡片实现单位可以是 angstrom 或 alat。推荐用 angstrom避免和 celldm 混淆。CONTROL calculation scf prefix strain outdir ./tmp / SYSTEM ibrav 0 nat 2 ntyp 1 ecutwfc 60 ecutrho 480 / ELECTRONS conv_thr 1.0d-8 / ATOMIC_SPECIES Si 28.0855 Si.pbe-n-rrkjus_psl.1.0.0.UPF CELL_PARAMETERS angstrom 5.4567 0.0000 0.0000 0.0000 5.4567 0.0000 0.0000 0.0000 5.4567 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逻辑说明ibrav 0表示手动指定晶格CELL_PARAMETERS angstrom后面三行就是形变后的晶格。ecutwfc和ecutrho需要做收敛测试一般 ecutrho 取 ecutwfc 的 8 倍 ultrasoft或 4 倍norm-conserving。参数上conv_thr建议 1e-8 以下应力对收敛更敏感。K 点密度对弹性常数影响很大8×8×8 是硅的常用值你的体系需要自己测。3.4 用 Python 驱动 QE 批量计算QE 的批量提交和 VASP 类似但输入文件是单个 pwscf.in可以用 Python 的 subprocess 或直接 shell 循环。import subprocess import os def run_qe_strain(strain_list, base_inputpwscf.in): 对每个应变生成输入并运行 QE with open(base_input) as f: template f.read() for s in strain_list: workdir fqe_strain_{s} os.makedirs(workdir, exist_okTrue) # 生成形变后的晶格行 lat apply_strain(np.array([[5.4567,0,0],[0,5.4567,0],[0,0,5.4567]]), [s, 0, 0, 0, 0, 0]) cell_lines \n.join( %.6f %.6f %.6f % tuple(row) for row in lat) new_input template.replace(CELL_PARAMETERS angstrom\n 5.4567 0.0000 0.0000\n 0.0000 5.4567 0.0000\n 0.0000 0.0000 5.4567, fCELL_PARAMETERS angstrom\n{cell_lines}) with open(f{workdir}/pwscf.in, w) as f: f.write(new_input) # 运行 QE subprocess.run([mpirun, -np, 16, pw.x, -in, pwscf.in], cwdworkdir, stdoutopen(f{workdir}/pwscf.out, w))逻辑说明template.replace把原始晶格行替换成形变后的注意替换字符串要和模板完全一致否则替换失败。subprocess.run的cwd参数让 QE 在各自目录运行避免临时文件冲突。参数上-np 16按机器调整QE 的并行效率对 k 点和平面波都有依赖建议先做小体系测试。4. 拟合与验证从应力数据到弹性常数4.1 线性拟合的 Python 实现拿到每个应变下的应力后对每个独立分量做线性拟合。以立方晶系为例σ_1 对 ε_1 的斜率就是 C11σ_2 对 ε_1 的斜率是 C12。import numpy as np from numpy.polynomial import polynomial as P def fit_elastic(strain_list, stress_list): strain_list: 每个形变的应变向量; stress_list: 对应应力矩阵 返回 C11, C12, C44 strains np.array(strain_list) stresses np.array(stress_list) # shape (N, 3, 3) # 提取正应力分量 s_xx stresses[:, 0, 0] s_yy stresses[:, 1, 1] # 对 xx 应变拟合 e_xx strains[:, 0] c11 np.polyfit(e_xx, s_xx, 1)[0] c12 np.polyfit(e_xx, s_yy, 1)[0] # 剪切需要单独施加 yz 应变的计算 # 这里假设 strains 里已有剪切数据 return c11, c12 # 示例数据 strains [[-0.01,0,0,0,0,0], [-0.005,0,0,0,0,0], [0,0,0,0,0,0], [0.005,0,0,0,0,0], [0.01,0,0,0,0,0]] stresses [[[-1.2,0,0],[0.3,0,0],[0.3,0,0]], # 示例值单位 GPa [[-0.6,0,0],[0.15,0,0],[0.15,0,0]], [[0,0,0],[0,0,0],[0,0,0]], [[0.6,0,0],[-0.15,0,0],[-0.15,0,0]], [[1.2,0,0],[-0.3,0,0],[-0.3,0,0]]] c11, c12 fit_elastic(strains, stresses) print(fC11 {c11:.1f} GPa, C12 {c12:.1f} GPa)逻辑说明np.polyfit(x, y, 1)返回一次多项式系数第一个是斜率。对 C11 用 σ_xx 对 ε_xx 拟合对 C12 用 σ_yy 对 ε_xx 拟合。参数上应变点要包含 0 且正负对称这样截距应该接近 0如果截距明显不为 0说明初始结构没弛豫好。剪切分量 C44 需要单独施加 yz 应变用 σ_yz 对 ε_yz 拟合。4.2 体积模量和泊松比的交叉验证拟合出 C11、C12 后立方晶系的体积模量 B (C11 2C12)/3这个值可以和状态方程拟合得到的 B0 对比。如果两者差超过 10%说明你的弹性常数有问题。这是最实用的自检手段。验证项公式合理偏差体积模量 B(C112C12)/3与 EOS 拟合 B0 差 10%剪切模量 GC44与文献差 15%泊松比 νC12/(C11C12)0.2~0.4 常见力学稳定性C11-C120, C112C120, C440必须满足注意力学稳定性判据不满足时说明你的结构在该应变下不稳定或者计算参数太差不要强行拟合。4.3 收敛测试怎么做才不浪费时间弹性常数对截断能和 K 点都敏感。常见做法是固定 K 点扫截断能看 C11 变化小于 1 GPa 就停再固定截断能扫 K 点。不要一上来就做全收敛测试先用一个中等参数跑一遍完整流程确认脚本没问题再回头做收敛。# 收敛测试的简单框架 for ecut in [40, 50, 60, 70, 80]: c11 run_one_strain(ecutecut, kpoints[8,8,8]) print(fecut{ecut}, C11{c11:.1f})逻辑说明run_one_strain是你封装好的单点计算函数返回拟合后的 C11。参数上截断能步长取 10 eVK 点步长取 2。如果 C11 在连续两个参数下变化小于 1 GPa可以认为收敛。5. 避坑与排查那些让弹性常数偏一倍的操作5.1 应力符号搞反C12 直接变负现象拟合出的 C12 是负值或者 C11 和 C12 的比值明显不对。原因QE 的应力符号和 VASP 相反QE 输出的是压缩为正直接拿来用会让斜率符号翻转。解决在解析 QE 应力时统一取反或者在拟合前对 QE 数据乘 -1。验证方法施加正应变时正应力应该为正拉伸如果为负就是符号错了。5.2 剪切应变用了工程应变而非张量应变现象C44 比文献值大一倍左右。原因Voigt 记法里的剪切分量 γ 对应张量应变 ε_ij γ/2如果施加形变时直接用了 γ 而没有除以 2实际应变就是两倍。解决在apply_strain里对非对角元除以 2或者施加应变时直接用张量形式。检查方法施加 ε_yz 0.005 时晶格矩阵的 (2,3) 和 (3,2) 元应该是 0.0025 倍的晶格常数而不是 0.005。5.3 原子位置没弛豫应力里混入内力现象应变为 0 时应力不为 0或者拟合截距明显偏离 0。原因初始结构的原子位置不是平衡位置内部有残余应力。解决先做一次完整弛豫ISIF3 或 QE 的 relax拿到平衡晶格和原子坐标再以此为起点施加应变。施加应变后如果只改晶格不弛豫原子对于某些结构如金刚石结构原子位置会偏离平衡需要做固定晶格的离子弛豫ISIF2。5.4 K 点太少导致应力噪声大现象不同应变下的应力值跳动拟合的 R² 很低。原因应力对 K 点采样比能量更敏感K 点太少时布里渊区积分误差大。解决把 K 点密度提高一倍或者用更密的 Monkhorst-Pack 网格。对于金属体系还需要加展宽。检查方法同一应变下用两套 K 点算应力差应该小于 0.1 GPa。5.5 应变幅度选得太大进入非线性区现象应力-应变曲线明显弯曲线性拟合残差大。原因应变超过 1% 后高阶弹性项贡献不可忽略。解决把应变幅度降到 ±0.005 以内或者用二次多项式拟合并取一次项系数。对于软材料如层状结构应变幅度要更小。检查方法分别用 ±0.005 和 ±0.01 拟合如果 C11 差超过 5%说明应变太大。6. 把两种引擎的数据合起来用一个可复用的后处理技巧实际项目里你可能一部分体系用 VASP 算一部分用 QE 算最后要放在一张图里对比。这时候最省事的做法不是改计算而是写一个统一的读取层把两边的输出都转成同一个数据结构一个包含strain6 分量和stress3×3GPa的记录列表。下面这个类可以直接抄。import numpy as np import re class StressStrainSet: def __init__(self): self.records [] # list of (strain_vec, stress_mat) def add_vasp(self, outcar, strain_vec): s read_vasp_stress(outcar) self.records.append((np.array(strain_vec), s)) def add_qe(self, out_path, strain_vec): s read_qe_stress(out_path) self.records.append((np.array(strain_vec), s)) def fit_component(self, i, j, strain_idx): 拟合 stress[i,j] 对 strain[strain_idx] 的斜率 e np.array([r[0][strain_idx] for r in self.records]) s np.array([r[1][i, j] for r in self.records]) # 去掉应变为 0 的点做拟合避免截距干扰 mask np.abs(e) 1e-8 if mask.sum() 2: raise ValueError(有效应变点不足) slope np.polyfit(e[mask], s[mask], 1)[0] return slope # 使用示例 ss StressStrainSet() ss.add_vasp(strain_0.005/OUTCAR, [0.005,0,0,0,0,0]) ss.add_vasp(strain_-0.005/OUTCAR, [-0.005,0,0,0,0,0]) ss.add_qe(qe_strain_0.005/pwscf.out, [0.005,0,0,0,0,0]) ss.add_qe(qe_strain_-0.005/pwscf.out, [-0.005,0,0,0,0,0]) c11 ss.fit_component(0, 0, 0) print(f合并拟合 C11 {c11:.1f} GPa)逻辑说明StressStrainSet把不同来源的数据统一成(strain_vec, stress_mat)元组fit_component按分量索引拟合。参数上strain_idx指定用哪个应变分量做自变量比如 0 对应 ε_xx3 对应 ε_yz。mask去掉应变为 0 的点是为了避免截距影响斜率如果你的数据截距确实为 0不去掉也可以。这个技巧的价值在于你不需要关心数据来自 VASP 还是 QE只要读取层正确后面的拟合和验证完全复用。我一般会把这个类存成一个elastic_utils.py每个项目直接 import。踩过的坑是 QE 的剪切分量顺序和 VASP 不同add_qe里已经做了映射但如果你自己写解析一定要核对输出文件的表头顺序。最后一个习惯每次拟合完先把 C11、C12、C44 和体积模量打印出来和文献值或 EOS 拟合值对一眼。如果对不上先查符号和单位再查应变幅度和 K 点。这套流程我跑了十几个体系九成的问题都出在这四个地方。希望帮到你。本文还有配套的精品资源点击获取