企业官网建设流程全解析

简介：面向材料科学与计算模拟方向的用户，聚焦VASP与Quantum Espresso（QE）的应力应变计算及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 ESPRESSO（QE）是当前最常用的两套平面波 DFT 代码，都能输出应变下的应力张量，但它们的输入格式、应力定义方式、提取路径完全不同。手动改 POSCAR 或 QE 输入做十几个变形构型，再一个个打开输出文件抄应力，这个流程不仅枯燥，还容易因为正负号约定或单位换算出错。一个 Python 脚本包可以把“生成应变构型—提交计算—提取应力—拟合弹性常数”串成一条流水线。这篇文章面向已经跑通 VASP 或 QE 自洽计算的工程师，讲清楚从应变施加到应力提取再到拟合的完整链路，顺带给出一套可以直接改写的小脚本。

2. 应力-应变计算的物理模型与 VASP/QE 实现差异

2.1 弹性区间的胡克定律与应变矩阵

晶体在弹性范围内的应力与应变满足广义胡克定律：

σ_ij = C_ijkl ε_kl

C_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.1（1 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.000000

QE 会直接把 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 或直接使用 ibrav=0 配合 celldm。常见做法是 ibrav=0 时把 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_atoms=True) # 写出新的 POSCAR，用于 VASP 计算 write('POSCAR_strain_001', atoms, format='vasp', vasp5=True, direct=True)

在 scale_atoms=True 时，所有原子坐标会随晶胞一起缩放。对于单轴应变，如果希望原子能充分驰豫来对应变响应，建议后续计算开启离子驰豫，但晶格长度保持固定，这样得到的应力才对应我们施加的应变。这个逻辑在 QE 里对应通过 calculation='relax' 但设置 cell_dynamics='none'，即固定晶胞，只动原子。

3. 用 Python 驱动 VASP 和 QE 批量产生应变-应力数据

3.1 准备工作目录与计算状态控制

批量计算的第一个常见坑是：应变后的结构需要先做一次高精度静态计算，还是直接做弛豫？我的经验是，如果应变幅度小于 3%，且原始结构已经是充分优化的，就直接做原子位置驰豫 + 固定晶格的静态计算，也就是：

• VASP：设置 NSW=50，ISIF=2（固定晶胞形状和体积，只驰豫原子位置），ISMEAR=-5 用于绝缘体或金属取整。
• QE：用 calculation='relax'，但是把 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_atoms=True) 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(f"VASP/{tag}", exist_ok=True) write(f"VASP/{tag}/POSCAR", atoms, format='vasp', vasp5=True, direct=True) os.makedirs(f"QE/{tag}", exist_ok=True) write(f"QE/{tag}/pw.in", atoms, format='espresso-in', pseudopotentials={'Si': 'Si.pbe-n-rrkjus_psl.1.0.0.UPF'}, kpts=(4, 4, 4), ecutwfc=40.0, input_data={'calculation': 'relax'})

这里tag字符串里x+0.010的形式很直观，方便后面对应输出文件。write到 QE 格式时，ASE 默认会给一个比较粗糙的输入，实际生产环境我会复制一份标准的 PWscf 输入模板，然后只替换 cell 和 atomic positions，避免 ASE 默认参数不够收敛。

3.2.1 QE 输入文件中固定晶格的关键参数

ASE 生成的 pw.in 里如果使用input_data={'calculation':'relax'}，它内部会设置cell_dynamics='none'，这能满足固定晶格的要求。如果你手动写模板，需要确保包含这几行：

&CONTROL calculation='relax' outdir='./tmp' prefix='strain_calc' / &CELL cell_dynamics='none' / &IONS ion_dynamics='bfgs' /

cell_dynamics='none'告诉 QE 只优化原子坐标，不改变晶胞。如果用了calculation='vc-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 2>&1 cd ../../ done

注意 INCAR 里必须包含ISIF=2、NSW=100、IBRION=2，并且ISMEAR根据体系选择。对于 QE 目录，循环里写：

for d in QE/*/; do cd $d mpirun -np 8 pw.x -i pw.in > pw.out 2>&1 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[i+1].split()[2:] + lines[i+2].split()[1:] + lines[i+3].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

同样地，读取应变信息也很简单，因为我们生成的目录名已经包含应变值，直接从目录名字符串解析 eps：

for d in strain_dirs: eps = float(d.split('+')[1]) # 从 'x+0.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*(e222+e333) 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), rcond=None)[0] print(f"C11={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 间隔）。对比时输出一个小表：

如果偏差集中在某个特定方向，注意检查该方向的应变矩阵是否施加正确，尤其是非对角元。

5.3 用收敛测试确认应力值稳定

应力是全局量，对电子步收敛阈值非常敏感。VASP 中 EDIFF 保守设置为 1e-6 eV 时应力才足够稳定；QE 中 conv_thr 建议取到 1e-8 或更低。改变这两个参数后应力变化应该在 0.1 GPa 以内，否则说明收敛标准太松。更简单的验证技巧：把零应变构型不加任何变形直接算一次，应力应该非常接近零（通常小于 0.5 kbar 是合理的），如果零应变应力都偏大，说明原始结构没有充分优化，后面所有数据都不可信。

最后，所有脚本和输出文件组织好之后，可以用一个 Python 脚本统一遍历所有数据目录，生成一张总表：每行是应变方向、应变值、应力张量六分量。这张表直接喂给拟合程序，比每次手动复制输出要稳得多。对后续加入新应变模式或新材料的计算，只需要改应变生成部分，提取和拟合代码完全复用。

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