企业官网建设流程全解析

简介：面向材料科学领域研究者和计算模拟爱好者，提供一套基于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 矩阵：下标 1=xx、2=yy、3=zz、4=yz、5=xz、6=xy。例如 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。

立方晶系习惯上用三个模式，见下表：

第一个模式单独给 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(f"C11 = {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(path="POSCAR"): """读取 POSCAR 的缩放因子和前 3 行格矢, 返回实际格矢矩阵""" with open(path, errors="ignore") 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, mode="uniaxial_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", errors="ignore") 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(f"POSCAR_{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(outcar="OUTCAR"): """提取最后一个总应力块, 返回 6 分量 (xx,yy,zz,yz,xz,xy), 单位 GPa""" with open(outcar, errors="ignore") as f: lines = f.readlines() stress = None for i, line in enumerate(lines): if "FORCES: max atom, stress" in line: try: block = lines[i+1:i+8] 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(outcar="OUTCAR"): """取最后一个 free energy TOTEN 行的总能, 单位 eV""" with open(outcar, errors="ignore") 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/zz，m[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. QE（Quantum 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(pwo="si.pwo"): """从 pw.x 输出摘应力, 返回 6 分量 (xx,yy,zz,yz,xz,xy), 单位 GPa""" with open(pwo, errors="ignore") 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(pwo="si.pwo"): """取最后一个 ! 行的总能, 单位换算成 eV""" with open(pwo, errors="ignore") as f: txt = f.read() matches = re.findall(r"!+\s+total 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 的输出项放在同一张表里，写解析脚本时不容易搞混：

同一结构的同一种应变模式，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 GPa，QE 算出 1032，刚好差 10 倍；或者能量法出来的模量是 1.6，压根不像 GPa 量级。 原因：OUTCAR 里应力默认 kBar（1 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 算的弹性常数基本都能互相印证。希望这套方法帮到你，少走几步弯路。

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