简介面向材料科学与计算模拟方向的用户聚焦VASP与Quantum EspressoQE的应力应变计算及Python后处理。压缩包共16个文件、约30KB包含8个.py脚本涵盖拉伸、剪切计算及绘图检查、4个输入文件.in、用于结构描述的POSCAR与POSCAR_rota以及说明文档README.md可支撑从DFT结构优化到应力应变曲线绘制的完整流程。已有931人浏览学习适合需要快速上手应力应变关系计算、理解VASP/QE输入输出及Python数据分析的科研人员。包内脚本区分VASP与QE两套路径分别提供无绘图和带绘图检查的版本并配套弛豫输入文件与旋转结构示例通过pymatgen、matplotlib等库用户可提取应力应变数据、绘制曲线并计算弹性模量、泊松比等参数是衔接第一性原理计算与力学性能分析的实用工具。1. 用 VASP 和 QE 算应力-应变关系为什么还要一个 Python 包做材料力学性能第一性原理计算时很少有比“拉一个方向看它怎么反抗”更直观的物理量了。我们想要的是应力-应变曲线横轴是施加的应变纵轴是计算得到的应力直线段的斜率就是杨氏模量再配合侧向收缩能拿到泊松比进而组装出弹性常数矩阵。VASP 和 Quantum ESPRESSOQE是当前最常用的两套平面波 DFT 代码都能输出应变下的应力张量但它们的输入格式、应力定义方式、提取路径完全不同。手动改 POSCAR 或 QE 输入做十几个变形构型再一个个打开输出文件抄应力这个流程不仅枯燥还容易因为正负号约定或单位换算出错。一个 Python 脚本包可以把“生成应变构型—提交计算—提取应力—拟合弹性常数”串成一条流水线。这篇文章面向已经跑通 VASP 或 QE 自洽计算的工程师讲清楚从应变施加到应力提取再到拟合的完整链路顺带给出一套可以直接改写的小脚本。2. 应力-应变计算的物理模型与 VASP/QE 实现差异2.1 弹性区间的胡克定律与应变矩阵晶体在弹性范围内的应力与应变满足广义胡克定律σ_ij C_ijkl ε_klC_ijkl 就是四阶弹性刚度张量对称性处理后可以用 6×6 的 Voigt 矩阵表示。实操上我们不需要从理论上推导太多对称性只要知道三斜晶系有 21 个独立弹性常数立方晶系只有 3 个C11、C12、C44。应变的具体定义是变形后的晶格矩阵除以原始晶格矩阵再减去单位矩阵ε (A·A0^(-1) (A·A0^(-1))^T) / 2 - I在实际计算里我们通常直接给晶格矩阵施加一个有限大小的应变比如沿某一轴拉伸 1%、2%、3%而不是用无穷小应变。这里有一个容易被忽略的点弹性常数是从应力-应变曲线在零应变附近的斜率拟合出来的有限应变量要足够小以保证在线性区内但又不能小到应力值被数值噪声淹没。我一般取 -3% 到 3%步长 1%针对不同方向分别做。2.2 VASP 与 QE 的应力输出约定差异VASP 的应力计算结果输出在 OUTCAR 中关键行是total stress格式如下total stress (kbars) -27.89031 -0.00000 -0.00000 -0.00000 -27.89031 -0.00000 -0.00000 -0.00000 -27.89031注意单位是 kbar也就是千巴换算成 GPa 要乘以 0.11 kbar 0.1 GPa。而且 VASP 输出的应力是“负值表示晶胞受压”它的符号约定与通常力学里的拉正压负有差异提取后要取相反数再用于作图。QE 的相对应输出在自洽计算的 stdout 或输出文件末尾类似total stress (Ry/bohr^3) (kbar) (GPa) -0.00289875 -0.00483530 -0.483530 -0.048353 -0.00000000 -0.00000000 0.000000 0.000000 -0.00000000 -0.00000000 0.000000 0.000000QE 会直接把 kbar 和 GPa 都打印出来对我们方便很多。但它输出的矩阵与 VASP 一样也是反号后的值也就是说 QE 的“应力”同样需要取负号才是通常意义上的应力。另外QE 的应力单位是 Ry/bohr^3手工解读时直接看 GPa 列即可。2.3 施加应变的两种常见做法2.3.1 直接修改晶格矩阵对 VASP把原始晶格矩阵乘以一个变形张量得到新的晶格后在 POSCAR 里替换。比如沿 x 轴施加 ε0.01 的单轴应变变形张量 F [[1.01, 0, 0], [0, 1, 0], [0, 0, 1]]。对 QE修改输入文件的 cell_parameters 块中对应行的坐标并同步修改 ibrav 或直接使用 ibrav0 配合 celldm。常见做法是 ibrav0 时把 A 矩阵整体替换避免手动算晶格常数。2.3.2 有限形变法在代码里的落地以 ASE 为例它可以直接读取 POSCAR 和 QE 输入然后施加应变。一个最简单的 VASP 应变构型生成脚本如下from ase.io import read, write from ase.constraints import ExpCellFilter from ase.filter import FrechetCellFilter import numpy as np atoms read(POSCAR_original) # 定义单轴应变 strain 0.01 F np.eye(3) F[0, 0] strain # 对晶格施加变形 atoms.set_cell(atoms.cell F, scale_atomsTrue) # 写出新的 POSCAR用于 VASP 计算 write(POSCAR_strain_001, atoms, formatvasp, vasp5True, directTrue)在 scale_atomsTrue 时所有原子坐标会随晶胞一起缩放。对于单轴应变如果希望原子能充分驰豫来对应变响应建议后续计算开启离子驰豫但晶格长度保持固定这样得到的应力才对应我们施加的应变。这个逻辑在 QE 里对应通过 calculationrelax 但设置 cell_dynamicsnone即固定晶胞只动原子。3. 用 Python 驱动 VASP 和 QE 批量产生应变-应力数据3.1 准备工作目录与计算状态控制批量计算的第一个常见坑是应变后的结构需要先做一次高精度静态计算还是直接做弛豫我的经验是如果应变幅度小于 3%且原始结构已经是充分优化的就直接做原子位置驰豫 固定晶格的静态计算也就是VASP设置 NSW50ISIF2固定晶胞形状和体积只驰豫原子位置ISMEAR-5 用于绝缘体或金属取整。QE用 calculationrelax但是把 cell_dynamics 关掉只让 ion 动。这样每个变形构型只需要跑一次不需要两步走。但对高各向异性或含氢键的体系建议先做纯静态再用沃罗诺伊约束或外部脚本读取应力否则应力的数值噪声会掩盖真正的线性关系。3.2 用 ASE 生成多组应变构型的完整脚本下面的脚本会遍历 x、y、z 三个轴每个轴生成从 -0.03 到 0.03 共 7 个应变并分别为 VASP 和 QE 写出输入文件from ase.io import read, write from ase.visualize import view import numpy as np import os atoms0 read(POSCAR_original) strains np.linspace(-0.03, 0.03, 7) axes [x, y, z] def apply_strain(atoms, axis, eps): F np.eye(3) if axis x: F[0, 0] eps elif axis y: F[1, 1] eps else: F[2, 2] eps atoms.set_cell(atoms.cell F, scale_atomsTrue) return atoms for axis in axes: for eps in strains: tag f{axis}{eps:.3f} atoms atoms0.copy() atoms apply_strain(atoms, axis, eps) os.makedirs(fVASP/{tag}, exist_okTrue) write(fVASP/{tag}/POSCAR, atoms, formatvasp, vasp5True, directTrue) os.makedirs(fQE/{tag}, exist_okTrue) write(fQE/{tag}/pw.in, atoms, formatespresso-in, pseudopotentials{Si: Si.pbe-n-rrkjus_psl.1.0.0.UPF}, kpts(4, 4, 4), ecutwfc40.0, input_data{calculation: relax})这里tag字符串里x0.010的形式很直观方便后面对应输出文件。write到 QE 格式时ASE 默认会给一个比较粗糙的输入实际生产环境我会复制一份标准的 PWscf 输入模板然后只替换 cell 和 atomic positions避免 ASE 默认参数不够收敛。3.2.1 QE 输入文件中固定晶格的关键参数ASE 生成的 pw.in 里如果使用input_data{calculation:relax}它内部会设置cell_dynamicsnone这能满足固定晶格的要求。如果你手动写模板需要确保包含这几行CONTROL calculationrelax outdir./tmp prefixstrain_calc / CELL cell_dynamicsnone / IONS ion_dynamicsbfgs /cell_dynamicsnone告诉 QE 只优化原子坐标不改变晶胞。如果用了calculationvc-relax就会改变晶胞导致应力对应变的关系失真。3.3 提交任务的 bash 循环写好后每个目录里都需要跑一次 VASP 或 pw.x。以下是一个简单的 bash 循环处理刚才生成的所有 VASP 目录#!/bin/bash for d in VASP/*/; do cd $d cp ../../../INCAR_strain ./ cp ../../../POTCAR ./ cp ../../../KPOINTS ./ mpirun -np 8 vasp_std run.out 21 cd ../../ done注意 INCAR 里必须包含ISIF2、NSW100、IBRION2并且ISMEAR根据体系选择。对于 QE 目录循环里写for d in QE/*/; do cd $d mpirun -np 8 pw.x -i pw.in pw.out 21 cd ../../ done这种批处理方式不需要 Python 参与跑完后统一提取输出。4. 提取应力应变数据并拟合弹性常数4.1 从 VASP 的 OUTCAR 提取应力张量VASP 的 OUTCAR 中每个离子步会打印一次total stress最好取最后一次离子步的值因为它对应收敛后的原子位置。提取脚本可以用正则表达式拆分三行六个数import re import numpy as np def read_vasp_stress(outcar_path): with open(outcar_path) as f: lines f.readlines() stress_lines [] for i, line in enumerate(lines): if total stress in line: vals lines[i1].split()[2:] lines[i2].split()[1:] lines[i3].split()[1:] stress_lines.append(np.array([float(v) for v in vals])) # 取最后一个离子步的输出 stress stress_lines[-1] # 单位 kbar - GPa并取负号得到正应力约定 stress_gpa -stress * 0.1 return stress_gpa同样地读取应变信息也很简单因为我们生成的目录名已经包含应变值直接从目录名字符串解析 epsfor d in strain_dirs: eps float(d.split()[1]) # 从 x0.010 解析 stress read_vasp_stress(f{d}/OUTCAR)4.2 从 QE 输出文件中提取应力QE 输出文件的total stress块在自洽循环结束后打印一次直接找这个关键字然后取 GPa 列的前三行即可def read_qe_stress(pw_out_path): with open(pw_out_path) as f: lines f.readlines() stress None for i, line in enumerate(lines): if total stress in line: base i 1 mat [] for row in range(3): parts lines[base row].split() # 文本中最后两列是 kbar 和 GPa取最后一列 GPa mat.append([float(parts[-1])]) stress np.array(mat).reshape(3, 3) # QE 输出同样需要取负号 return -stress这里直接把每行的最后一个数字当作 GPa。要注意 QE 输出文件中“total stress”块可能出现多次有些是在初始计算打印的有些是最后的。用pw_out文件的最后一次匹配结果通常最可靠上面代码用stress变量覆盖最后留的就是最后一次。4.3 拟合应力-应变关系提取到每个方向的应力分量后选择对应轴上的对角元应变与应力做线性拟合。比如说沿 x 轴应变时取 σ_xx 与 ε_xx。在弹性区非线性明显时我会限制拟合范围在 ±2% 以内并对原始数据做一次高阶项检查。from scipy.stats import linregress import numpy as np eps np.array([-0.03, -0.02, -0.01, 0.0, 0.01, 0.02, 0.03]) sigma np.array([...]) # 从 OUTCAR 或 QE 输出中提取的对应应力 mask np.abs(eps) 0.02 slope, intercept, rvalue, pvalue, stderr linregress(eps[mask], sigma[mask]) young_modulus slope # 单位 GPa这里的 slope 就是沿该轴的杨氏模量。如果是立方晶系可以用三个独立方向的结果结合应变矩阵反推 C11 和 C12但最稳妥的做法是使用多个应变模式如体积应变、剪切应变联合求解。4.3.1 用最小二乘同时拟合 C11、C12、C44下面给出一个可扩展的拟合思路通过构造不同应变模式把应力-应变关系写成线性方程组然后使用numpy.linalg.lstsq求解import numpy as np # 每个应变模式的行e11, e22, e33, 2e23, 2e13, 2e12 # 以立方晶系为例C11, C12, C44 为变量 A [] b [] for mode, (e111, e222, e333, e23, e13, e12) in enumerate(strains_modes): sigma1 stress_matrix[mode][0,0] sigma2 stress_matrix[mode][1,1] sigma3 stress_matrix[mode][2,2] sigma4 stress_matrix[mode][1,2] # 方程sigma1 C11*e111 C12*(e222e333) A.append([e111, e222 e333, 0.0]) b.append(sigma1) A.append([e222, e111 e333, 0.0]) b.append(sigma2) A.append([e333, e111 e222, 0.0]) b.append(sigma3) A.append([0.0, 0.0, e23]) b.append(sigma4) C np.linalg.lstsq(np.array(A), np.array(b), rcondNone)[0] print(fC11{C[0]:.2f} GPa, C12{C[1]:.2f} GPa, C44{C[2]:.2f} GPa)这段代码里 A 的每一行对应胡克定律在某个应变模式下的一个分量方程把多组数据堆叠成超定方程组最小二乘能得到一组平均意义上的弹性常数。要注意剪切应变和工程应变的换算通常在生成应变构型时需要在应变矩阵的非对角元上使用半系数。5. 验证计算结果的三个实用技巧5.1 检查能量-应变曲线的平滑度应力是能量对应变的一阶导如果应力-应变曲线噪声大直接看能量-应变曲线更明显。把每个应变构型的总能减去零应变总能画成抛物线形状如果某个点明显偏离抛物线说明该应变下的构型可能没有收敛需要检查电子步或原子坐标是否真的达到平衡。对 VASP 脚本可以这样提取能量grep without entropy OUTCAR | tail -1 | awk {print $4}对 QE 则是grep ! pw.out | tail -1 | awk {print $5}然后把对应能量值放入 Python 中做二次拟合二次项系数实际对应弹性常数的部分贡献可以用它作为 sanity check和应力线性拟合的斜率交叉验证。5.2 横向对比 VASP 与 QE 的结果同一套应变构型在两种代码中算出来的应力通常有 1-2% 的偏差主要来自平面波截断、赝势类型和布里渊区取样差异。如果某个方向两种代码结果偏差超过 5%优先检查 K 点密度和截断能是否有区别。建议两边都用收敛测试确定的同一套 K 点密度如 0.02 Å^-1 间隔。对比时输出一个小表应变方向VASP 斜率 (GPa)QE 斜率 (GPa)偏差x171.2170.50.4%y50.352.13.5%如果偏差集中在某个特定方向注意检查该方向的应变矩阵是否施加正确尤其是非对角元。5.3 用收敛测试确认应力值稳定应力是全局量对电子步收敛阈值非常敏感。VASP 中 EDIFF 保守设置为 1e-6 eV 时应力才足够稳定QE 中 conv_thr 建议取到 1e-8 或更低。改变这两个参数后应力变化应该在 0.1 GPa 以内否则说明收敛标准太松。更简单的验证技巧把零应变构型不加任何变形直接算一次应力应该非常接近零通常小于 0.5 kbar 是合理的如果零应变应力都偏大说明原始结构没有充分优化后面所有数据都不可信。最后所有脚本和输出文件组织好之后可以用一个 Python 脚本统一遍历所有数据目录生成一张总表每行是应变方向、应变值、应力张量六分量。这张表直接喂给拟合程序比每次手动复制输出要稳得多。对后续加入新应变模式或新材料的计算只需要改应变生成部分提取和拟合代码完全复用。本文还有配套的精品资源点击获取 Python脚本实现VASP与QE应力应变计算及弹性常数自动拟合
简介面向材料科学与计算模拟方向的用户聚焦VASP与Quantum EspressoQE的应力应变计算及Python后处理。压缩包共16个文件、约30KB包含8个.py脚本涵盖拉伸、剪切计算及绘图检查、4个输入文件.in、用于结构描述的POSCAR与POSCAR_rota以及说明文档README.md可支撑从DFT结构优化到应力应变曲线绘制的完整流程。已有931人浏览学习适合需要快速上手应力应变关系计算、理解VASP/QE输入输出及Python数据分析的科研人员。包内脚本区分VASP与QE两套路径分别提供无绘图和带绘图检查的版本并配套弛豫输入文件与旋转结构示例通过pymatgen、matplotlib等库用户可提取应力应变数据、绘制曲线并计算弹性模量、泊松比等参数是衔接第一性原理计算与力学性能分析的实用工具。1. 用 VASP 和 QE 算应力-应变关系为什么还要一个 Python 包做材料力学性能第一性原理计算时很少有比“拉一个方向看它怎么反抗”更直观的物理量了。我们想要的是应力-应变曲线横轴是施加的应变纵轴是计算得到的应力直线段的斜率就是杨氏模量再配合侧向收缩能拿到泊松比进而组装出弹性常数矩阵。VASP 和 Quantum ESPRESSOQE是当前最常用的两套平面波 DFT 代码都能输出应变下的应力张量但它们的输入格式、应力定义方式、提取路径完全不同。手动改 POSCAR 或 QE 输入做十几个变形构型再一个个打开输出文件抄应力这个流程不仅枯燥还容易因为正负号约定或单位换算出错。一个 Python 脚本包可以把“生成应变构型—提交计算—提取应力—拟合弹性常数”串成一条流水线。这篇文章面向已经跑通 VASP 或 QE 自洽计算的工程师讲清楚从应变施加到应力提取再到拟合的完整链路顺带给出一套可以直接改写的小脚本。2. 应力-应变计算的物理模型与 VASP/QE 实现差异2.1 弹性区间的胡克定律与应变矩阵晶体在弹性范围内的应力与应变满足广义胡克定律σ_ij C_ijkl ε_klC_ijkl 就是四阶弹性刚度张量对称性处理后可以用 6×6 的 Voigt 矩阵表示。实操上我们不需要从理论上推导太多对称性只要知道三斜晶系有 21 个独立弹性常数立方晶系只有 3 个C11、C12、C44。应变的具体定义是变形后的晶格矩阵除以原始晶格矩阵再减去单位矩阵ε (A·A0^(-1) (A·A0^(-1))^T) / 2 - I在实际计算里我们通常直接给晶格矩阵施加一个有限大小的应变比如沿某一轴拉伸 1%、2%、3%而不是用无穷小应变。这里有一个容易被忽略的点弹性常数是从应力-应变曲线在零应变附近的斜率拟合出来的有限应变量要足够小以保证在线性区内但又不能小到应力值被数值噪声淹没。我一般取 -3% 到 3%步长 1%针对不同方向分别做。2.2 VASP 与 QE 的应力输出约定差异VASP 的应力计算结果输出在 OUTCAR 中关键行是total stress格式如下total stress (kbars) -27.89031 -0.00000 -0.00000 -0.00000 -27.89031 -0.00000 -0.00000 -0.00000 -27.89031注意单位是 kbar也就是千巴换算成 GPa 要乘以 0.11 kbar 0.1 GPa。而且 VASP 输出的应力是“负值表示晶胞受压”它的符号约定与通常力学里的拉正压负有差异提取后要取相反数再用于作图。QE 的相对应输出在自洽计算的 stdout 或输出文件末尾类似total stress (Ry/bohr^3) (kbar) (GPa) -0.00289875 -0.00483530 -0.483530 -0.048353 -0.00000000 -0.00000000 0.000000 0.000000 -0.00000000 -0.00000000 0.000000 0.000000QE 会直接把 kbar 和 GPa 都打印出来对我们方便很多。但它输出的矩阵与 VASP 一样也是反号后的值也就是说 QE 的“应力”同样需要取负号才是通常意义上的应力。另外QE 的应力单位是 Ry/bohr^3手工解读时直接看 GPa 列即可。2.3 施加应变的两种常见做法2.3.1 直接修改晶格矩阵对 VASP把原始晶格矩阵乘以一个变形张量得到新的晶格后在 POSCAR 里替换。比如沿 x 轴施加 ε0.01 的单轴应变变形张量 F [[1.01, 0, 0], [0, 1, 0], [0, 0, 1]]。对 QE修改输入文件的 cell_parameters 块中对应行的坐标并同步修改 ibrav 或直接使用 ibrav0 配合 celldm。常见做法是 ibrav0 时把 A 矩阵整体替换避免手动算晶格常数。2.3.2 有限形变法在代码里的落地以 ASE 为例它可以直接读取 POSCAR 和 QE 输入然后施加应变。一个最简单的 VASP 应变构型生成脚本如下from ase.io import read, write from ase.constraints import ExpCellFilter from ase.filter import FrechetCellFilter import numpy as np atoms read(POSCAR_original) # 定义单轴应变 strain 0.01 F np.eye(3) F[0, 0] strain # 对晶格施加变形 atoms.set_cell(atoms.cell F, scale_atomsTrue) # 写出新的 POSCAR用于 VASP 计算 write(POSCAR_strain_001, atoms, formatvasp, vasp5True, directTrue)在 scale_atomsTrue 时所有原子坐标会随晶胞一起缩放。对于单轴应变如果希望原子能充分驰豫来对应变响应建议后续计算开启离子驰豫但晶格长度保持固定这样得到的应力才对应我们施加的应变。这个逻辑在 QE 里对应通过 calculationrelax 但设置 cell_dynamicsnone即固定晶胞只动原子。3. 用 Python 驱动 VASP 和 QE 批量产生应变-应力数据3.1 准备工作目录与计算状态控制批量计算的第一个常见坑是应变后的结构需要先做一次高精度静态计算还是直接做弛豫我的经验是如果应变幅度小于 3%且原始结构已经是充分优化的就直接做原子位置驰豫 固定晶格的静态计算也就是VASP设置 NSW50ISIF2固定晶胞形状和体积只驰豫原子位置ISMEAR-5 用于绝缘体或金属取整。QE用 calculationrelax但是把 cell_dynamics 关掉只让 ion 动。这样每个变形构型只需要跑一次不需要两步走。但对高各向异性或含氢键的体系建议先做纯静态再用沃罗诺伊约束或外部脚本读取应力否则应力的数值噪声会掩盖真正的线性关系。3.2 用 ASE 生成多组应变构型的完整脚本下面的脚本会遍历 x、y、z 三个轴每个轴生成从 -0.03 到 0.03 共 7 个应变并分别为 VASP 和 QE 写出输入文件from ase.io import read, write from ase.visualize import view import numpy as np import os atoms0 read(POSCAR_original) strains np.linspace(-0.03, 0.03, 7) axes [x, y, z] def apply_strain(atoms, axis, eps): F np.eye(3) if axis x: F[0, 0] eps elif axis y: F[1, 1] eps else: F[2, 2] eps atoms.set_cell(atoms.cell F, scale_atomsTrue) return atoms for axis in axes: for eps in strains: tag f{axis}{eps:.3f} atoms atoms0.copy() atoms apply_strain(atoms, axis, eps) os.makedirs(fVASP/{tag}, exist_okTrue) write(fVASP/{tag}/POSCAR, atoms, formatvasp, vasp5True, directTrue) os.makedirs(fQE/{tag}, exist_okTrue) write(fQE/{tag}/pw.in, atoms, formatespresso-in, pseudopotentials{Si: Si.pbe-n-rrkjus_psl.1.0.0.UPF}, kpts(4, 4, 4), ecutwfc40.0, input_data{calculation: relax})这里tag字符串里x0.010的形式很直观方便后面对应输出文件。write到 QE 格式时ASE 默认会给一个比较粗糙的输入实际生产环境我会复制一份标准的 PWscf 输入模板然后只替换 cell 和 atomic positions避免 ASE 默认参数不够收敛。3.2.1 QE 输入文件中固定晶格的关键参数ASE 生成的 pw.in 里如果使用input_data{calculation:relax}它内部会设置cell_dynamicsnone这能满足固定晶格的要求。如果你手动写模板需要确保包含这几行CONTROL calculationrelax outdir./tmp prefixstrain_calc / CELL cell_dynamicsnone / IONS ion_dynamicsbfgs /cell_dynamicsnone告诉 QE 只优化原子坐标不改变晶胞。如果用了calculationvc-relax就会改变晶胞导致应力对应变的关系失真。3.3 提交任务的 bash 循环写好后每个目录里都需要跑一次 VASP 或 pw.x。以下是一个简单的 bash 循环处理刚才生成的所有 VASP 目录#!/bin/bash for d in VASP/*/; do cd $d cp ../../../INCAR_strain ./ cp ../../../POTCAR ./ cp ../../../KPOINTS ./ mpirun -np 8 vasp_std run.out 21 cd ../../ done注意 INCAR 里必须包含ISIF2、NSW100、IBRION2并且ISMEAR根据体系选择。对于 QE 目录循环里写for d in QE/*/; do cd $d mpirun -np 8 pw.x -i pw.in pw.out 21 cd ../../ done这种批处理方式不需要 Python 参与跑完后统一提取输出。4. 提取应力应变数据并拟合弹性常数4.1 从 VASP 的 OUTCAR 提取应力张量VASP 的 OUTCAR 中每个离子步会打印一次total stress最好取最后一次离子步的值因为它对应收敛后的原子位置。提取脚本可以用正则表达式拆分三行六个数import re import numpy as np def read_vasp_stress(outcar_path): with open(outcar_path) as f: lines f.readlines() stress_lines [] for i, line in enumerate(lines): if total stress in line: vals lines[i1].split()[2:] lines[i2].split()[1:] lines[i3].split()[1:] stress_lines.append(np.array([float(v) for v in vals])) # 取最后一个离子步的输出 stress stress_lines[-1] # 单位 kbar - GPa并取负号得到正应力约定 stress_gpa -stress * 0.1 return stress_gpa同样地读取应变信息也很简单因为我们生成的目录名已经包含应变值直接从目录名字符串解析 epsfor d in strain_dirs: eps float(d.split()[1]) # 从 x0.010 解析 stress read_vasp_stress(f{d}/OUTCAR)4.2 从 QE 输出文件中提取应力QE 输出文件的total stress块在自洽循环结束后打印一次直接找这个关键字然后取 GPa 列的前三行即可def read_qe_stress(pw_out_path): with open(pw_out_path) as f: lines f.readlines() stress None for i, line in enumerate(lines): if total stress in line: base i 1 mat [] for row in range(3): parts lines[base row].split() # 文本中最后两列是 kbar 和 GPa取最后一列 GPa mat.append([float(parts[-1])]) stress np.array(mat).reshape(3, 3) # QE 输出同样需要取负号 return -stress这里直接把每行的最后一个数字当作 GPa。要注意 QE 输出文件中“total stress”块可能出现多次有些是在初始计算打印的有些是最后的。用pw_out文件的最后一次匹配结果通常最可靠上面代码用stress变量覆盖最后留的就是最后一次。4.3 拟合应力-应变关系提取到每个方向的应力分量后选择对应轴上的对角元应变与应力做线性拟合。比如说沿 x 轴应变时取 σ_xx 与 ε_xx。在弹性区非线性明显时我会限制拟合范围在 ±2% 以内并对原始数据做一次高阶项检查。from scipy.stats import linregress import numpy as np eps np.array([-0.03, -0.02, -0.01, 0.0, 0.01, 0.02, 0.03]) sigma np.array([...]) # 从 OUTCAR 或 QE 输出中提取的对应应力 mask np.abs(eps) 0.02 slope, intercept, rvalue, pvalue, stderr linregress(eps[mask], sigma[mask]) young_modulus slope # 单位 GPa这里的 slope 就是沿该轴的杨氏模量。如果是立方晶系可以用三个独立方向的结果结合应变矩阵反推 C11 和 C12但最稳妥的做法是使用多个应变模式如体积应变、剪切应变联合求解。4.3.1 用最小二乘同时拟合 C11、C12、C44下面给出一个可扩展的拟合思路通过构造不同应变模式把应力-应变关系写成线性方程组然后使用numpy.linalg.lstsq求解import numpy as np # 每个应变模式的行e11, e22, e33, 2e23, 2e13, 2e12 # 以立方晶系为例C11, C12, C44 为变量 A [] b [] for mode, (e111, e222, e333, e23, e13, e12) in enumerate(strains_modes): sigma1 stress_matrix[mode][0,0] sigma2 stress_matrix[mode][1,1] sigma3 stress_matrix[mode][2,2] sigma4 stress_matrix[mode][1,2] # 方程sigma1 C11*e111 C12*(e222e333) A.append([e111, e222 e333, 0.0]) b.append(sigma1) A.append([e222, e111 e333, 0.0]) b.append(sigma2) A.append([e333, e111 e222, 0.0]) b.append(sigma3) A.append([0.0, 0.0, e23]) b.append(sigma4) C np.linalg.lstsq(np.array(A), np.array(b), rcondNone)[0] print(fC11{C[0]:.2f} GPa, C12{C[1]:.2f} GPa, C44{C[2]:.2f} GPa)这段代码里 A 的每一行对应胡克定律在某个应变模式下的一个分量方程把多组数据堆叠成超定方程组最小二乘能得到一组平均意义上的弹性常数。要注意剪切应变和工程应变的换算通常在生成应变构型时需要在应变矩阵的非对角元上使用半系数。5. 验证计算结果的三个实用技巧5.1 检查能量-应变曲线的平滑度应力是能量对应变的一阶导如果应力-应变曲线噪声大直接看能量-应变曲线更明显。把每个应变构型的总能减去零应变总能画成抛物线形状如果某个点明显偏离抛物线说明该应变下的构型可能没有收敛需要检查电子步或原子坐标是否真的达到平衡。对 VASP 脚本可以这样提取能量grep without entropy OUTCAR | tail -1 | awk {print $4}对 QE 则是grep ! pw.out | tail -1 | awk {print $5}然后把对应能量值放入 Python 中做二次拟合二次项系数实际对应弹性常数的部分贡献可以用它作为 sanity check和应力线性拟合的斜率交叉验证。5.2 横向对比 VASP 与 QE 的结果同一套应变构型在两种代码中算出来的应力通常有 1-2% 的偏差主要来自平面波截断、赝势类型和布里渊区取样差异。如果某个方向两种代码结果偏差超过 5%优先检查 K 点密度和截断能是否有区别。建议两边都用收敛测试确定的同一套 K 点密度如 0.02 Å^-1 间隔。对比时输出一个小表应变方向VASP 斜率 (GPa)QE 斜率 (GPa)偏差x171.2170.50.4%y50.352.13.5%如果偏差集中在某个特定方向注意检查该方向的应变矩阵是否施加正确尤其是非对角元。5.3 用收敛测试确认应力值稳定应力是全局量对电子步收敛阈值非常敏感。VASP 中 EDIFF 保守设置为 1e-6 eV 时应力才足够稳定QE 中 conv_thr 建议取到 1e-8 或更低。改变这两个参数后应力变化应该在 0.1 GPa 以内否则说明收敛标准太松。更简单的验证技巧把零应变构型不加任何变形直接算一次应力应该非常接近零通常小于 0.5 kbar 是合理的如果零应变应力都偏大说明原始结构没有充分优化后面所有数据都不可信。最后所有脚本和输出文件组织好之后可以用一个 Python 脚本统一遍历所有数据目录生成一张总表每行是应变方向、应变值、应力张量六分量。这张表直接喂给拟合程序比每次手动复制输出要稳得多。对后续加入新应变模式或新材料的计算只需要改应变生成部分提取和拟合代码完全复用。本文还有配套的精品资源点击获取