VASP与QE应力应变计算全解析：从DFT参数到Python拟合

简介：面向材料科学领域的DFT计算学习者，这份资料将第一性原理软件VASP与Quantum Espresso中的力学计算流程，封装成可直接运行的Python脚本，适合已有一定计算基础、希望自动化处理应力应变数据的用户。压缩包共16个文件，以八个Python脚本为核心，另配有QE和VASP的输入文件、POSCAR结构文件以及说明文档，整体大小只有三十KB，轻量且便于按需修改。脚本分别对应拉伸与剪切两类形变场景，同时提供VASP与QE两种版本，可以读取输出文件中的应力应变信息，调用绘图模块生成曲线，并进一步计算弹性模量、泊松比等参数，有助于对比两种软件在不同体系下的计算表现。目前已有九百三十一人学习下载，适合需要快速搭建应变计算流程、深入理解DFT力学分析的科研工作者与高年级学生。

1. 从弹性常数到应力应变关系，DFT 计算的起点不是脚本而是力学模型

很多人拿到“使用VASP和QE计算应力和应变关系”这类任务，第一反应是去搜 VASP 的 INCAR 参数、QE 的输入文件模板，或者直接找 Python 拟合脚本。但实际在组里待过三年以上的人都知道，这两个程序在应力应变计算上的最大差异，根本不在于输入格式，而在于它们对“应力”和“应变”这两个物理量的定义方式不同。VASP 默认输出的是 virial 应力张量，单位是 kBar，而 QE 以 Ry/Bohr³ 为应力单位，两者换算因子记错了，后面拟合出来的弹性常数直接差一个数量级。

这个标题里的“和”字很关键，它暗示的不是二选一，而是要在两套代码之间做交叉验证。常见做法是，用 VASP 做高精度的结构优化并计算完整的弹性常数张量，再用 QE 做应力-应变扫点，验证同一材料在线性区间内的响应是否一致。对 5 年以上从业者来说，真正容易翻车的地方反而不是 DFT 参数，而是应变矩阵的施加方式、参考结构的对称性保持，以及 Python 拟合时线性区间切片的选择。这篇文章就按照从原理到拟合的路径，把这条链路上每个环节的参数含义和坑位讲清楚。

2. 应变施加的力学定义与 VASP 的弹性常数计算路径

2.1 从有限应变理论到晶格矩阵的雅可比行列式

应变张量 ε 的施加方式直接决定计算结果的物理意义。在 DFT 代码里，通常通过修改晶格矩阵 L 来实现应变：L' = L · (I + ε)，这里的 ε 是对称张量，包含 6 个独立分量。对正交晶系，这 6 个分量直接对应三个轴向拉伸和三个剪切；对非正交晶系，必须注意应变矩阵是在分数坐标还是笛卡尔坐标下作用。VASP 在处理 POSCAR 时，默认使用笛卡尔坐标下的晶格矢量，所以你在加应变时，一定要先确认原胞的基矢是否正交化过。

提示：对六方晶系和三角晶系，千万不要把 POSCAR 里的晶格常数直接乘上 (1 + ε)，因为晶格矢量存在非对角分量。正确做法是用原胞基矢构造 3×3 矩阵，左乘应变矩阵，再重新把矩阵写成 POSCAR 格式。

计算弹性常数 Cij 的标准做法是对晶格施加若干组有限应变，将总能量对应变展开：

E(ε) = E₀ + V₀ · Σ σᵢεᵢ + (V₀/2) · Σ Cᵢⱼεᵢεⱼ + O(ε³)

其中 V₀ 是平衡体积，σᵢ 是应力。如果你只需要应力-应变关系，就不需要做二阶差分拟合能量，直接读取 VASP 输出的 stress 矩阵，然后画出应力的某个分量与应变分量的线性段，斜率就是相应的弹性常数。这个思路在下载的压缩包里通常对应一个strain_calc目录，里面有一组预先设定好形变量的 POSCAR。

2.2 VASP 的 INCAR 参数怎么设才能同时拿应力和能量

要在 VASP 里既拿到准确应力，又不把结构优化带偏，关键是关掉对称性对力的混合。下面这套参数是我在 fcc 铜、bcc 铁、六方钛和钙钛矿氧化物上反复调过的起点：

System = strain_stress_test ISTART = 0 ICHARG = 2 PREC = Accurate ENMAX = 1.3 * ENCUT_default EDIFF = 1E-7 EDIFFG = -0.002 ISMEAR = 0 SIGMA = 0.05 IBRION = 2 ISIF = 3 NSW = 60 POTIM = 0.2 ISYM = 0 LREAL = .FALSE.

逻辑说明：ISYM=0是必须加的，因为应变会降低体系对称性，如果让 VASP 自动检测对称性，它可能把原本独立的应变分量合并掉，导致你无法区分某个剪切应力分量。ISIF=3表示在弹性计算中让晶胞体积和形状一起优化，但注意这是用在IBRION=2的有限差分弹性常数计算中，不是用在你要做应力-应变扫点的那一步。如果真的要做扫点，ISIF应该改为 2 或直接固定晶格，只放开原子位置。

参数表如下：

用这套参数跑完后，OUTCAR 里找到TOTAL-FORCE (eV/Angst)上方有FORCES: max atom, min atom，而应力在STRESS标签下，单位是 kBar。注意 VASP 输出的应力是直接量，不是部分占据修正后的量，这与其他程序输出名义应力不同。

2.3 VASP linux 怎么测试：一个最小算例验证应力输出

很多刚接触 VASP linux 环境的人不知道如何验证安装是否正常，更不知道如何验证应力计算是否可靠。最快的方法是用 fcc 铝做一个单轴拉伸测试，因为铝的弹性常数是各向异性极小且 C11、C12、C44 已知值明确。

# 从官网下载 Al 的 POTCAR（或从已有赝势库复制） cp ~/potpaw/POTCAR.Al/POTCAR . # 构造 1x1x1 晶胞，晶格常数 4.05 Å，方向沿 x cat > POSCAR << EOF Al_fcc_strain_test 4.05 0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0 1 Direct 0.0 0.0 0.0 EOF

然后创建 INCAR 和 KPOINTS。KPOINTS 用 12×12×12 的 Monkhorst-Pack，位移 0.0。跑完后用grep "TOTAL-S" OUTCAR看应力张量。对未加应变的平衡结构，应力应当接近零；如果输出应力大于 0.5 kBar，通常说明 KPOINTS 或 ENCUT 不够密集，原因是对金属铝，应力张量的收敛速度显著慢于总能。测试时如果发现应力矩阵对角元都在同一个数量级、但非对角元不接近零，那就先查 POSCAR 的对称性和 ISYM 是否真的生效了。

3. 用 QE 做应力-应变扫点时的输入文件设计与参数匹配

3.1 QE 的 stress 计算单元与单位换算

QE 计算应力的方式与 VASP 不同，它通过 Pulay 修正后的密度矩阵直接计算应力张量，在pw.x输入中只要添加stress = .true.就可以了。这个开关会触发程序在自洽循环结束后额外计算一次应力张量，输出单位是 kBar 除以一个换算因子——准确来说是 Ry/Bohr³，换算到 kBar 需要乘以 147.105。

如果是让 PWscf 在每一个 SCF 步之后都输出应力，需要用testress = .true.配合etot_conv_thr控制；否则只在最后输出一次。这里有一个常见陷阱：在结构优化vc-relax过程中，stress = .true.才会在每个优化步计算应力，而在普通relax中它默认不重复计算，只输出起始结构的那一次。

说到单位，还有个细节：QE 输出的应力张量多数版本的cell和atom单位都是 Ry/Bohr³，但pw.x的output里会有一行把它转成 kBar。你在写 Python 解析脚本时，直接读total stress那一行，它已经换算过。如果用pw_export.x或第三方库读内部张量，就要自己乘 147.105。

3.2 QE 输入文件中决定应力准确度的三个关键参数

QE 中对应力收敛影响最大的参数不是 k 点，而是ecutwfc和ecutrho。ecutrho对应电荷密度截止，对压强和应力的影响比ecutwfc更明显，因为应力里的动能部分和高波矢分量耦合更强。经验值是让ecutrho取ecutwfc的 8 到 12 倍（对应 norm-conserving 赝势与 ultra-soft 赝势的差别）。如果你使用 PAW 型设置，QE 会用ecutrho = 4 * ecutwfc左右就能收敛，但应力响应要求高时建议拉到 8 倍再做收敛测试。

&CONTROL calculation = 'scf' prefix = 'strain' outdir = './tmp' pseudo_dir = '../pseudo' verbosity = 'high' stress = .true. / &SYSTEM ibrav = 0 nat = 2 ntyp = 1 ecutwfc = 60 ecutrho = 480 occupations = 'smearing' smearing = 'cold' degauss = 0.01 / &ELECTRONS conv_thr = 1.0e-8 mixing_beta = 0.3 / CELL_PARAMETERS (angstrom) 2.86 0.00 0.00 1.43 2.48 0.00 0.00 0.00 3.40 ATOMIC_SPECIES Ti 47.867 Ti.pbe-spfn-rrkjus_psl.0.1.UPF ATOMIC_POSITIONS (crystal) Ti 0.0 0.0 0.0 Ti 0.5 0.5 0.5 K_POINTS {automatic} 11 11 7 0 0 0

参数说明：smearing='cold'对金属和过渡金属氧化物的应力计算往往比默认的 gaussian 收敛得更平滑，尤其是对 d 电子体系的应力张量。degauss设得太大（超过 0.02）会产生明显的人工应力展宽，但设得太小会让 k 点离散化引入噪声。对 hcp Ti 这种结构，k 点网格要按倒空间对应关系取，所以用了 11×11×7 而不是均匀立方网格。

3.3 用 python 操控 QE 生成系列应变结构

在下载的 zip 里，一般会有一个run_qe_strain.py脚本，它做的事情是：对平衡结构施加不同幅度的应变，改写CELL_PARAMETERS，重复提交pw.x，并把输出文件里的应力和应变提取到 CSV。这里给出一个最小可用的生成器，不用依赖 pymatgen，只用 numpy 就能改晶格矩阵：

import numpy as np import subprocess, re base_cell = np.array([ [2.86, 0.00, 0.00], [1.43, 2.48, 0.00], [0.00, 0.00, 3.40] ], dtype=float) strain_values = np.linspace(-0.02, 0.02, 9) results = [] for eps in strain_values: # 默认沿 x 轴单轴应变 strain_mat = np.eye(3) strain_mat[0, 0] += eps new_cell = base_cell @ strain_mat with open('pw.strain.in', 'r') as f: inp = f.read() new_inp = re.sub( r'CELL_PARAMETERS.*?(?=ATOMIC_SPECIES)', f'CELL_PARAMETERS (angstrom)\n' + '\n'.join( f'{v[0]:.6f} {v[1]:.6f} {v[2]:.6f}' for v in new_cell) + '\n', inp, flags=re.S ) with open('pw.strain_run.in', 'w') as f: f.write(new_inp) subprocess.run(['mpirun', '-np', '4', 'pw.x', '-in', 'pw.strain_run.in'], capture_output=True, text=True) with open('pw.strain_run.out', 'r') as f: out = f.read() stress_kbar = re.findall(r'total\s+stress\s+=\s+([-?\d\.E+]+)', out) stress_ry = float(stress_kbar[0]) * 147.105 # 这一行期间要确认单位 results.append((eps, stress_ry)) np.savetxt('strain_stress_qe.csv', np.array(results), delimiter=',', header='strain_x,stress_kbar')

代码里的@是矩阵乘法，不是逐元素乘，很多新手在这里直接把矩阵的每个元素都加了eps，导致剪切分量被带入。re.sub只替换从CELL_PARAMETERS到ATOMIC_SPECIES之间的文本段，避免误伤其他标签。单位换算那行要特别注意：有些版本的 QE 的输出已经是 kBar，再乘 147.105 就错了；稳妥做法是检查stress行的注释，或者对平衡结构跑一次，看应力是否接近零来判断是否需要换算。

4. 用 Python 解析 VASP 和 QE 输出文件的异同

4.1 VASP OUTCAR 与 QE 输出文件的字段映射

VASP 的 OUTCAR 里STRESS出现在电子步收敛之后，格式是 3×3 矩阵，单位 kBar；QE 的输出里total stress也是一个 3×3 矩阵，但默认单位是 Ry/Bohr³，只有加上stress = .true.并打开verbosity='high'才会在末尾输出换算后的 kBar 版本。字段映射如下：

一个实用的做法是写一个解析函数，输入文件路径和代码类型，输出统一单位下的应力张量。下面这段代码处理两种格式，提取应力，并把非对角分量排除掉，因为对大多数正交晶系，你关心的只有σxx, σyy, σzz。

import re def read_stress_from_output(output_path, code='vasp'): with open(output_path, 'r') as f: text = f.read() if code == 'vasp': match = re.search(r'STRESS\s*\(kBar\)\s*(.*?)(?=\n\s*\n)', text, re.S) if match: block = match.group(1).strip().split('\n') else: # 新版本 vasp.6.x 的格式有差异 match = re.search(r'stress matrix\s*\(kBar\)\s*(.*?)\n\s*\n', text, re.S) block = match.group(1).strip().split('\n') stress_arr = [] for line in block: parts = line.replace('-', ' -').split() stress_arr.append([float(x) for x in parts if re.match(r'^-?\d', x)]) return np.array(stress_arr) elif code == 'qe': match = re.search(r'total\s+stress\s+=\s*([-\d.E+]+)\s+([-\d.E+]+)\s+([-\d.E+]+)', text) if match: return np.array([[float(match.group(1)), 0, 0], [0, float(match.group(2)), 0], [0, 0, float(match.group(3))]]) return None

这个函数的问题在于正则表达式对 VASP 的输出依赖空行位置，如果 ISIF 或参数设置不同导致输出格式变化，就会匹配失败。更好的方式是按行扫描STRESS之后连续 3 行非空文本；但这里给出简化版本作为起点已经足够。

4.2 从 QE 输出中提取应变对应力并用 pandas 组装数据表

实际处理中不只是单个结构，而是一系列应变的扫描结果。把 QE 输出和 VASP 输出统一组装成 DataFrame，是后续拟合最省心的路径。这里建议用 pandas，虽然标题里只是说 Python，但 pymatgen 和 ASE 都不适合做数据清洗时的人为介入——它们封装的太深，容易掩盖单位错误。

import pandas as pd data_list = [] for eps, fpath in strain_output_files: stress_matrix = read_stress_from_output(fpath, code='qe') sigma_xx = stress_matrix[0, 0] # 已经是 kBar data_list.append({'strain_xx': eps, 'stress_xx': sigma_xx, 'stress_yy': stress_matrix[1, 1]}) df = pd.DataFrame(data_list) df['strain_xx'] = df['strain_xx'].astype(float) df.to_csv('combined_strain_stress.csv', index=False)

逻辑说明：把应力和应变关系做成一个 DataFrame 后，后续的线性拟合就完全不需要再碰文本文件了。stress_yy在单轴应变下也不为零，它是对横向响应的探测，如果横向应力有异常大幅变化，说明晶格方向设置错了——单轴应变应当只影响纵向分量，横向分量保持接近零（如果用的是固定横向尺寸的模式）。

4.3 Python 数据分析与可视化中容易踩到的小数精度问题

对 0.5% 量级的应变，VASP 的应力输出只有三位有效数字，并且应力张量的输出精度由IO格式决定。OUTCAR 里的应力通常保留到小数点后 1 位（kBar），对弹性模量几个 GPa 量级的材料，1 kBar 的误差会带来 10% 左右的模量误差。解决办法是不要在文本层面提高精度，而是用 VASP 的DFT+U或大截断能让应力在电子步上收敛到更小的值，然后直接在 OUTCAR 里用程序读原始数值，不要用人工复制粘贴。

有一种常见错误：用 Python 的float()转换1.23456E+02这类字符串没问题，但如果 QE 输出的负号前有空格，比如-1.234E+01，split()会把它正确切分；可如果是 Fortran 的write格式产生-1.234E+01后面没有多余空格，re.findall就会匹配出错误的分组。稳妥的解析方式是先按正则抓出所有科学计数法的 token，再重新排列成 3×3。

5. 用 Python 拟合线性区间，验证应力应变关系的可靠性

弹性常数拟合本质上是对一组 (ε, σ) 数据做线性回归，但难点在于“线性区间”的界定。对金属材料，应变超过 ±3% 后应力-应变曲线迅速偏离线性；对陶瓷和半导体，线性区间更窄，约 ±1%。所以我一般不在全区间做拟合，而是先画散点图，再用简单的误差条选择区间。

import numpy as np from scipy import stats df = pd.read_csv('combined_strain_stress.csv') # 只取应变绝对值在 0.02 以内的点做初始拟合 mask = np.abs(df['strain_xx']) < 0.02 x = df.loc[mask, 'strain_xx'].values y = df.loc[mask, 'stress_xx'].values linreg = stats.linregress(x, y) sigma = linreg.stderr print(f"C11 = {linreg.slope:.2f} kBar = {linreg.slope * 0.1:.2f} GPa") print(f"R^2 = {linreg.rvalue**2:.4f}, stderr = {sigma:.2f} kBar")

线性回归得到的斜率就是弹性常数 C11（如果施加的是单轴 x 应变）。将 kBar 转成 GPa 是乘 0.1，因为 1 kBar = 0.1 GPa。很多人在这一步把单位搞混，因为 VASP 里 1 kBar 等于 1 GPa 的十分之一，从 OUTCAR 里直接读出的 C11 如果 600 左右，单位是 kBar，换成 GPa 就变成 60，显然不对，要对上预期量级。

更稳妥的区间选择是用残差分析：把全区间数据做二次拟合，看哪一段二次项系数几乎为零。也可以用交叉验证的方式——逐步扩大中间区间，观察斜率变化率小于 1% 的区间就是可用区间。这种做法比直接用固定阈值可靠得多。

# 从中心点向两侧扩展，观察斜率漂移 sorted_df = df.sort_values('strain_xx') for half_window in np.linspace(0.005, 0.03, 6): mask = (sorted_df['strain_xx'] > -half_window) & \ (sorted_df['strain_xx'] < half_window) sub = sorted_df[mask] if len(sub) >= 3: slope, intercept, r, p, se = stats.linregress( sub['strain_xx'], sub['stress_xx']) print(f"±{half_window*100:.2f}%: C11 = {slope*0.1:.2f} GPa, " f"R² = {r**2:.5f}, se = {se:.2f} kBar")

这段代码的用意不是挑选一个数字，而是观察 C11 随着窗口变化是否稳定。如果你看到 C11 随窗口扩大单调下降，说明已经进入塑性或非线性区间，必须缩小范围。如果 C11 上下波动大于 3%，大概率是计算噪声或 k 点密度不够，而不是拟合问题。

关于拟合法还有一个进阶想法：不要把应力和应变直接做最小二乘，而是用对称弹性常数约束。对各向同性材料，C12 = C11 - 2C44，如果你同时拟合了三组应变方向的数据，可以联立方程组解出 C11、C12、C44，并用线性代数验证矩阵的正定性。这种约束拟合在用 Python 的scipy.optimize做的时候不会增加太多复杂度，但能显著减少单轴拟合的方向偏差。

应力和应变关系计算的最后一环不是拟合本身，而是对比 VASP 和 QE 的结果。用前面构造的脚本分别计算同一应变序列下的应力，画在同一张图上，如果两条曲线在各应变点处的应力差超过 2-3 kBar，优先检查赝势版本和ecutrho的收敛性，而不是急着改 Python 拟合代码。两个代码对同一个物理结果的差异应控制在 1-2% 以内才算计算可信。对六方晶系，还要额外检查 c/a 比是否在应变后仍保持优化值——这是单元晶胞计算中最常被忽略的系统误差。

本文还有配套的精品资源，点击获取