简介这套压缩包面向材料科学领域需要使用VASP与QE开展力学性质计算的研究者和学生聚焦应力-应变关系的模拟与数据处理兼顾第一性原理仿真与Python后处理。包内共16个文件以8个Python脚本和4个输入文件为主另含POSCAR与POSCAR_rota结构文件、README说明等整体仅30KB轻量清晰。脚本覆盖拉伸与剪切两种典型变形模式并分别提供VASP和QE版本其中又分为带绘图与不带绘图两类可直接读取输出文件提取应力应变数据绘制曲线并计算弹性模量、泊松比等参数。已有931人学习下载。通过阅读脚本、输入文件和说明可理解通过修改晶格参数施加应变并求应力响应的DFT流程也可借鉴作者用Python组织、拟合和可视化力学数据的代码思路适合正在学习第一性原理计算或希望提升材料计算分析能力的读者参考。1. 从应力-应变关系出发VASP 和 QE 的弹性计算怎么分工材料的弹性常数决定了结构在外力下怎么变形也决定了声子热导、热膨胀系数和力学稳定性这些衍生性质。第一性原理算弹性常数实质是给晶体施加已知应变计算体系总能量对应变幅度的二阶导数得到刚度张量 Cij。VASP 和 Quantum ESPRESSO 都能做这件事差别在输入格式、应力输出方式和自动化流程。做材料计算的人通常两套代码都装VASP 用于高精度 PAW 势计算QE 用模守恒或超软赝势跑同样体系作为交叉验证。Python 在这里不是计算引擎而是粘合剂——批量生成变形晶格、扫描应变幅度、从输出文件里抽能量和应力、做二次多项式拟合。本文把这条路径完整走一遍命令和脚本可以直接改体系参数复用。先说明一个常见误区很多人觉得“算应力-应变关系”就要用弛豫后读取应力张量实际上从能量-应变曲线拟合弹性常数更稳健。直接读应力在某些代码里依赖应力定理的数值实现且需要多个非线性应变分量分解而能量是变分量对收敛标准和小应变扰动更敏感。所以下面的实操主线是固定应变值 → 单点能计算 → 收集能量-应变数据 → 二次多项式拟合而不是逐点读取应力再作商。2. 应变施加策略与能量-应变拟合两套代码共同的计算基础2.1 为什么用能量-应变曲线而不是直接读应力单位晶胞在微小变形下应变能可以展开为E(ε) E₀ V₀·σᵢεᵢ (V₀/2)·Cᵢⱼεᵢεⱼ …其中 V₀ 是平衡体积Cᵢⱼ 是二阶弹性常数σᵢ 是初始应力平衡结构下应近似为零。如果结构已经优化到力收敛线性项消失弹性常数直接来自二次项系数。这是能量法的核心假设。直接读应力的做法在某些时候不太可靠VASP 的 OUTCAR 里应力张量单位是 kBQE 的 stress 张量输出有自己的约定Ry/au³转换时容易出错更重要是应力对 k 点采样和截断能的敏感度高于总能量。能量法把问题归结为“用足够精度算很多个单点能”然后做数值微分。代价是需要多算 2 到 6 个应变方向收益是结果可复现、可反查拟合曲线。所以我的建议是用能量法求弹性常数用应力张量做收敛性检查和结构优化后的验证。两者结合比单看任何一者都可靠。2.2 Voigt 记号与应变矩阵构造二阶弹性刚度张量 C 是 4 阶张量有 81 个分量。由于对称性Voigt 记号把它压缩成 6×6 矩阵应力 σ₁, σ₂, σ₃ 对应 xx, yy, zzσ₄, σ₅, σ₆ 对应 yz, xz, xy。应变张量同理但剪切应变附带因子 2εxx e₁εyy e₂εzz e₃γyz e₄γxz e₅γxy e₆注意 VASP 和 QE 内部用的都是完整二阶张量但描述“施加什么应变”要用 Voigt 向量然后再还原成晶格变形矩阵。下面这个 Python 函数做的事就是这个还原import numpy as np def voigt_to_matrix(voigt): 输入 Voigt 应变向量 (e1,e2,e3,e4,e5,e6)返回对称应变张量矩阵 e1, e2, e3, e4, e5, e6 voigt return np.array([ [e1, e6/2, e5/2], [e6/2, e2, e4/2], [e5/2, e4/2, e3 ] ])注意主对角线就是正应变非对角元要除以 2因为矩阵对称化后张量分量 εxy ½·γxy。很多人第一次从这个函数转到 VASP 的 POSCAR 变形矩阵时少除了这个 2导致施加的剪切应变翻倍拟合出来的 C44 直接两倍偏差。2.3 应变模式与独立弹性常数数量不同晶系需要计算的应变模式数量不一样。立方晶系只有 3 个独立弹性常数 C11, C12, C44但施加同一个等体积应变的一个矩阵并不总能同时分离出所有分量。常用做法是施加三组不同模式应变模式Voigt 向量拟合目标单轴拉伸(δ, 0, 0, 0, 0, 0)C11双轴拉伸(δ, δ, 0, 0, 0, 0)C11 C12剪切(0, 0, 0, δ, 0, 0) 或 (0,0,0,0,0,δ)C44三组应变得到三个能量-应变曲线联立求解就得到三个独立常数。对于六角晶系需要五个独立常数对应至少 5 组应变模式正交晶系是 9 个三斜晶系则要 21 个。常见自动化工具会自动枚举这些模式手动做的时候先用晶系对称性简化。实际操作里每组应变取 59 个幅度从 -2% 到 2%步长 0.5%。少于 5 个点时二次拟合的自由度不够多于 9 个点则计算量翻倍而精度提升有限。2.4 线性应变与工程应变的选择对同一个变形用“线性应变” ε (a-a₀)/a₀ 和“工程应变” η (a²-a₀²)/(2a₀²) 会得到不同的能量解析形式。第一性原理代码里改变晶格常数时POSCAR 或 QE 的 celldm 参数直接按比例缩放所以默认对应的是线性应变。小应变范围内两者只差高阶项但拟合时如果用了工程应变而没用对应的体积变化公式会引入可察觉的系统误差。推荐统一用线性应变因为它直接对应代码里晶格常数缩放系数。设平衡晶格矢为 h₀3×3 矩阵施加应变 ε 后新晶格矢是h h₀ · (I ε)Python 实现def apply_strain(lattice_matrix, voigt): lattice_matrix: 3x3 晶格矢量行为 a,b,c voigt: Voigt 应变向量 strain voigt_to_matrix(voigt) deformed lattice_matrix (np.eye(3) strain) return deformed注意矩阵乘法方向POSCAR 里晶格矢量是行向量还是列向量取决于代码约定。VASP 的 POSCAR 第二行是缩放系数第三到第五行是晶格矢量每行一个矢量上面对应的变形矩阵乘法是 h_new h_old · (Iε)。有些脚本用列向量写法结果相同但要小心 numpy 的 方向颠倒带来的转置问题。3. VASP 实操从 POSCAR 变形到 OUTCAR 应力张量3.1 完整流程概览VASP 下计算弹性常数的标准流程是高精度结构优化得到平衡晶格常数和原子坐标固定晶格对每个应变模式生成多个变形 POSCAR对每个变形结构做静态自洽计算IALGO 默认的电子步收敛从 OUTCAR 或 OSZICAR 提取总能拟合能量-应变曲线得到弹性常数其中第一步和第三步要分开不能直接拿优化未收敛的结构施加应变否则初始应力项不为零能量展开式里的线性项会污染二次拟合。3.2 用 Python 生成 VASP 变形 POSCARimport numpy as np def read_poscar(filename): 读取 POSCAR 晶格矢量部分返回 3x3 数组 with open(filename) 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 write_poscar(filename, lattice, base_lines): 根据原 POSCAR 的元素信息写出新晶格 with open(filename, w) as f: for i, line in enumerate(base_lines[:2]): f.write(line) f.write(f {lattice[0][0]:.8f} {lattice[0][1]:.8f} {lattice[0][2]:.8f}\n) f.write(f {lattice[1][0]:.8f} {lattice[1][1]:.8f} {lattice[1][2]:.8f}\n) f.write(f {lattice[2][0]:.8f} {lattice[2][1]:.8f} {lattice[2][2]:.8f}\n) for line in base_lines[5:]: f.write(line) # 主逻辑读取平衡结构施加单轴应变批量写出 base_poscar open(POSCAR_opt).readlines() h0 read_poscar(POSCAR_opt) strains [x/1000 for x in range(-20, 21, 5)] # -2% 到 2% for i, delta in enumerate(strains): deformed apply_strain(h0, (delta, 0, 0, 0, 0, 0)) write_poscar(fPOSCAR_e1_{i:02d}, deformed, base_poscar)读入平衡结构后把缩放系数并入射格矢量里避免后续处理把缩放系数弄丢。输出文件名建议带应变方向和序号方便后续脚本正则匹配。3.3 INCAR 与 KPOINTS 参数选择弹性常数计算对收敛精度要求高于普通结构优化。推荐一组可复用的 INCARENCUT 1.3 * ENMAX # 平面波截断能VASP 推荐值的 1.3 倍 EDIFF 1E-8 # 电子步自洽收敛标准 EDIFFG -0.01 # 只做静态计算时此项不影响结果 ISMEAR -5 # 绝缘体/半导体用四面体方法 SIGMA 0.2 # ISMEAR0 时才有意义-5 时忽略 PREC Accurate IBRION -1 # 静态计算 NSW 0 LREAL .FALSE. # 精确计算禁用实空间投影对于金属体系ISMEAR 用 0 或 1 配合 SIGMA0.1四面体方法对部分占据不适用。ENCUT 检查办法是用平衡结构分别跑 1.0×ENMAX 和 1.3×ENMAX能量差小于 1 meV/atom 就够用。KPOINTS 至少要用 Monkhorst-Pack 自动网格密度参考cubic 晶系 k-spacing ≤ 0.03 Å⁻¹非立方晶系相应减少。三斜晶系或者超胞大的体系用 Gamma-centered 网格往往更稳。关键检查是增加 k 点数后 Cij 变化小于 1 GPa。3.4 批量提交计算并从 OUTCAR 提取数据生成完一组 POSCAR 后批量复制 INCAR、KPOINTS、POSCAR 和 POTCAR 到子目录然后提交 VASPfor d in strain_*; do cp INCAR KPOINTS POTCAR $d/ cd $d mpirun -np 16 vasp_std run.log cd .. done运行完后提取总能量for d in strain_*; do E$(grep free energy $d/OUTCAR | tail -1 | awk {print $5}) echo $d $E done注意 VASP 输出的能量有多种标签free energy是有限温度自由能energy(sigma-0)是外推的 0K 能量。静态计算用后者更准确因为 ISMEAR-5 时free energy和energy(sigma-0)数值不同。如果发现某个应变后结构对外加应变有明显应力响应再看 OUTCAR 的total stress块确认应变后残余应力与外加应变的对应关系。这能帮助判断 POSCAR 变形矩阵是否方向搞反。def parse_outcar_energy(filename): with open(filename) as f: for line in f: if energy(sigma-0) in line: return float(line.split()[1].split()[0])提取结果后直接进入第五章的拟合脚本。但先看 QE 那边怎么做两者组合才能互为对照。4. QE 实操用 pw.x 的 celldm 参数扫描应变4.1 从 VASP 到 QE 的输入结构转换QE 的 pw.x 输入文件以SYSTEMnamelist 里的晶格参数描述晶体。常用ibrav指定 Bravais 晶格类型celldm或a, b, c, cosAB等指定晶格常量。和 VASP 直接给 3×3 晶格矢量不同QE 用ibrav0时可以直接给CELL_PARAMETERS三行矢量这和 VASP 的 POSCAR 几乎对应。CONTROL calculation scf prefix si outdir ./tmp pseudo_dir ./pseudo/ verbosity high / SYSTEM ibrav 0 nat 2 ntyp 1 ecutwfc 40 ecutrho 320 occupationssmearing smearingmp degauss0.01 / ELECTRONS conv_thr 1.0d-9 mixing_beta 0.4 / ATOMIC_SPECIES Si 28.0855 Si.pz-vbc.UPF CELL_PARAMETERS {angstrom} 3.8400 0.0000 0.0000 1.9200 3.3260 0.0000 1.9200 1.1087 3.1356 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 0ecutwfc对应 VASP 的 ENCUTecutrho默认取 812 倍ecutwfc超软赝势需要更高。QE 里用degauss做 smearing此时平面波截断能对弹性常数的影响比 VASP 更明显务必做截断能收敛测试。4.2 用 Python 修改 CELL_PARAMETERS 施加应变等价于 VASP 的 POSCAR 变形直接操作 QE 的 CELL_PARAMETERS 块import re def apply_strain_qe(input_text, voigt): 读入 QE 输入文件字符串修改 CELL_PARAMETERS 晶格矢量 lines input_text.splitlines() for i, line in enumerate(lines): name line.split()[0].upper() if name CELL_PARAMETERS: cells [] for j in range(i1, i4): cells.append([float(x) for x in lines[j].split()]) h0 np.array(cells) deformed apply_strain(h0, voigt) for j in range(3): lines[i1j] f{deformed[j][0]:.8f} {deformed[j][1]:.8f} {deformed[j][2]:.8f} break return \n.join(lines)然后用循环生成多个输入文件for i, delta in enumerate(strains): text open(scf.in).read() deformed_text apply_strain_qe(text, (delta, 0, 0, 0, 0, 0)) open(fscf_e1_{i:02d}.in, w).write(deformed_text)注意 QE 的CELL_PARAMETERS默认单位是 alat如果在SYSTEM里写了ibrav0且没有celldm(1)单位默认为 Å规范起见在CELL_PARAMETERS那行显式写{angstrom}避免单位歧义。4.3 批量运行与数据提取for f in scf_e1_*.in; do mpirun -np 8 pw.x -in $f ${f%.in}.out done从输出文件提取能量和压力grep ! scf_e1_*.out grep P scf_e1_*.out | tail -1QE 的P 行给出外加流体静压力单位是 kbar。这个值在等体积应变下应为 0如果参考结构是平衡态。如果P的数值明显漂移说明参考结构没有充分优化。QE 没有像 VASP OUTCAR 那样方便的总能量提取energy(sigma-0)但!行给出的total energy在固定占据数下就是最终能量。金属体系用抹布时能量对 degauss 有依赖所以每组应变都用同一套 smearing 参数避免系统性偏差。还可以从PW的total stress (GPa)输出块读取完整应力张量单位已经转换成 GPa方便和 VASP 对比。def parse_qe_energy(filename): with open(filename) as f: lines f.readlines() for line in reversed(lines): if line.startswith(!): return float(line.split()[-2])这个函数从文件尾部反向搜索!行得到总能。注意 QE 输出里!行格式是! total energy ... Ry用负索引提取能正确处理多行输出。4.4 VASP 和 QE 数据格式对比项目VASPQE晶格描述POSCAR 三行矢量CELL_PARAMETERS 三行矢量截断能控制ENCUT单位 eVecutwfc/ecutrho单位 Ry能量单位eVRy1 Ry 13.6057 eV应力单位OUTCAR 里 kB输出里 GPa 或 kbar施加载荷后的计算必须重新做自洽单次 scf 即可处理剪切应变直接变形直接变形关键坑位是两个代码能量单位差了 13.6 倍拟合弹性常数时要统一成 eV 或者 J/m³不然 Cij 数值天然差两个数量级。5. Python 后处理拟合弹性常数并画出应力-应变线性段5.1 二次多项式拟合与 Cij 提取把 VASP 和 QE 提取的能量放到一起按应变值排序然后用 numpy 做最小二乘二次拟合import numpy as np def fit_elastic(data, strain_name, volume): data: [(strain, energy_eV), ...] volume: 平衡体积单位 Å^3 返回 (C_value_GPa, 拟合误差) strain np.array([d[0] for d in data]) energy np.array([d[1] for d in data]) coefs np.polyfit(strain, energy, 2) # 能量单位 eV体积单位 Å^3转 GPa 的系数 eV_per_Ang3_to_GPa 160.21766208 C 2 * coefs[0] / volume * eV_per_Ang3_to_GPa # 手动算不确定度 fit_e energy - np.polyval(coefs, strain) sigma np.sqrt(np.mean(fit_e**2)) error sigma / np.sqrt(len(strain)) * 2 * eV_per_Ang3_to_GPa / volume return C, error二次项系数乘 2 再除以体积才是弹性常数这个 2 来自能量展开的 ½ 因子。注意体积单位换算VASP 的 CONTCAR 里体积单位是 ųQE 的 CELL_PARAMETERS 如果用 Å 则一致如果用 alat 则需要乘 (celldm(1)*0.529177)³。np.polyfit默认返回最高次项在前所以coefs[0]是二次项系数。5.2 用应力-应变曲线的线性拟合交叉验证能量法算出弹性常数后可以再从应力-应变数据验证。VASP 的 OUTCAR 里total stress是应力张量用输出的应变对应的应力分量 σ 对 ε 做线性拟合斜率应该等于对应 Cijimport numpy as np stress_eps [] for i, delta in enumerate(strains): stress read_outcar_stress(fstrain_e1_{i:02d}/OUTCAR) stress_eps.append((delta, stress[0, 0])) # xx 分量 # stress 单位 kB1 kB 0.1 GPa stress_eps np.array(stress_eps) stress_eps[:, 1] * 0.1 # 转 GPa coef_slope np.polyfit(stress_eps[:, 0], stress_eps[:, 1], 1) print(fC11 from stress-strain slope: {coef_slope[0]:.1f} GPa)这个数值和能量法拟合出来的一致性应落在数值噪声范围内一般是 12 GPa 的偏差。如果偏差达到 5 GPa 以上首先检查应变定义应力是共轭于有限应变还是工程应变VASP 的应力输出对应的是无限小应变张量所以和线性应变一致。5.3 批量脚本与结果汇总下面把整个流程的可复用脚本框架画出来实际使用时按代码包里的目录结构调整路径import json results {} for code in [vasp, qe]: results[code] {} for mode in [uniaxial, biaxial, shear]: data collect_data(code, mode) C, err fit_elastic(data, volume, mode) results[code][mode] C with open(Cij_results.json, w) as f: json.dump(results, f, indent2)格式化成 JSON 的好处是后续不管喂给 origin 画图还是喂给机器学习模型都不用重新解析文本。5.4 画能量-应变曲线把拟合曲线和原始数据画在一张图里一目了然import matplotlib.pyplot as plt strain_points [d[0] for d in data] energy_points [d[1] for d in data] sm np.linspace(min(strain_points), max(strain_points), 50) plt.scatter(strain_points, energy_points, labelDFT data) plt.plot(sm, np.polyval(coefs, sm), r-, labelquadratic fit) plt.xlabel(Strain) plt.ylabel(Energy (eV)) plt.legend() plt.savefig(f{code}_{mode}_fit.png, dpi150)拟合曲线如果偏离数据点超过 1 meV说明应变范围取得太大非线性效应进入了三次项范围这时把应变范围收窄到 ±1%或者拟合时加三阶项看系数是否显著。正常材料在 2% 应变内三次项贡献小于 0.1%但如果体系本身有相变比如铁电材料要缩短范围。6. 验证弹性常数结果的四个收敛性检查弹性常数的精确值对计算参数极其敏感跑完一组数据不做检查直接发表很容易出问题。以下四个检查按顺序做任何一环不过都要返工。正式计算前先跑一组单点收敛测试固定 k 点从 1.0×ENMAX 扫到 1.3×ENMAX记录能量差再固定 ENCUT从 ×1 到 ×2 加密 k 点网格。如果截断能每提升 5% 能量变化超过 1 meV/atom弹性常数变化会超过 2 GPa。QE 同理ecutwfc 从 30 扫到 60 Ry记录 Cij 变化趋势。合格标准是相邻两个截断能下 Cij 差小于 0.5 GPa。应变幅度检查看二次拟合残差把残差对应变画出来如果呈 U 型分布说明二阶模型不够加个三次项试试如果残差呈现周期波动多半是 SCF 收敛精度不够EDIFF 或 conv_thr 还需要收紧一两个量级。参考结构必须是弛豫到压力为 0 的平衡晶格。判断方法是应变0 的那个点VASP 的静水压力小于 0.5 kbarQE 的 P 小于 0.5 kbar。如果初始压力太大拟合斜率和二次项都会失真弹性常数可能偏小约 510%。最终要把算出的 Cij 组成对称矩阵检查特征值是否全部为正对应 Born 力学稳定性判据。立方晶系的判据是 C44 0C11 |C12|C112C12 0。如果 C44 算出负值优先怀疑剪切应变定义里的 γ/2 因子没处理好其次检查 k 点密度是否太低。交叉验证用 VASP 和 QE 两套结果对比偏差大于 3 GPa 时排查赝势和截断能设置。用 0.5% 步长 5 个点扫一遍如果 Cij 随步长变化超过 1 GPa说明应变区间选得不合适重新选步长。本文还有配套的精品资源点击获取
第一性原理弹性常数计算：VASP与QE能量-应变拟合实战
简介这套压缩包面向材料科学领域需要使用VASP与QE开展力学性质计算的研究者和学生聚焦应力-应变关系的模拟与数据处理兼顾第一性原理仿真与Python后处理。包内共16个文件以8个Python脚本和4个输入文件为主另含POSCAR与POSCAR_rota结构文件、README说明等整体仅30KB轻量清晰。脚本覆盖拉伸与剪切两种典型变形模式并分别提供VASP和QE版本其中又分为带绘图与不带绘图两类可直接读取输出文件提取应力应变数据绘制曲线并计算弹性模量、泊松比等参数。已有931人学习下载。通过阅读脚本、输入文件和说明可理解通过修改晶格参数施加应变并求应力响应的DFT流程也可借鉴作者用Python组织、拟合和可视化力学数据的代码思路适合正在学习第一性原理计算或希望提升材料计算分析能力的读者参考。1. 从应力-应变关系出发VASP 和 QE 的弹性计算怎么分工材料的弹性常数决定了结构在外力下怎么变形也决定了声子热导、热膨胀系数和力学稳定性这些衍生性质。第一性原理算弹性常数实质是给晶体施加已知应变计算体系总能量对应变幅度的二阶导数得到刚度张量 Cij。VASP 和 Quantum ESPRESSO 都能做这件事差别在输入格式、应力输出方式和自动化流程。做材料计算的人通常两套代码都装VASP 用于高精度 PAW 势计算QE 用模守恒或超软赝势跑同样体系作为交叉验证。Python 在这里不是计算引擎而是粘合剂——批量生成变形晶格、扫描应变幅度、从输出文件里抽能量和应力、做二次多项式拟合。本文把这条路径完整走一遍命令和脚本可以直接改体系参数复用。先说明一个常见误区很多人觉得“算应力-应变关系”就要用弛豫后读取应力张量实际上从能量-应变曲线拟合弹性常数更稳健。直接读应力在某些代码里依赖应力定理的数值实现且需要多个非线性应变分量分解而能量是变分量对收敛标准和小应变扰动更敏感。所以下面的实操主线是固定应变值 → 单点能计算 → 收集能量-应变数据 → 二次多项式拟合而不是逐点读取应力再作商。2. 应变施加策略与能量-应变拟合两套代码共同的计算基础2.1 为什么用能量-应变曲线而不是直接读应力单位晶胞在微小变形下应变能可以展开为E(ε) E₀ V₀·σᵢεᵢ (V₀/2)·Cᵢⱼεᵢεⱼ …其中 V₀ 是平衡体积Cᵢⱼ 是二阶弹性常数σᵢ 是初始应力平衡结构下应近似为零。如果结构已经优化到力收敛线性项消失弹性常数直接来自二次项系数。这是能量法的核心假设。直接读应力的做法在某些时候不太可靠VASP 的 OUTCAR 里应力张量单位是 kBQE 的 stress 张量输出有自己的约定Ry/au³转换时容易出错更重要是应力对 k 点采样和截断能的敏感度高于总能量。能量法把问题归结为“用足够精度算很多个单点能”然后做数值微分。代价是需要多算 2 到 6 个应变方向收益是结果可复现、可反查拟合曲线。所以我的建议是用能量法求弹性常数用应力张量做收敛性检查和结构优化后的验证。两者结合比单看任何一者都可靠。2.2 Voigt 记号与应变矩阵构造二阶弹性刚度张量 C 是 4 阶张量有 81 个分量。由于对称性Voigt 记号把它压缩成 6×6 矩阵应力 σ₁, σ₂, σ₃ 对应 xx, yy, zzσ₄, σ₅, σ₆ 对应 yz, xz, xy。应变张量同理但剪切应变附带因子 2εxx e₁εyy e₂εzz e₃γyz e₄γxz e₅γxy e₆注意 VASP 和 QE 内部用的都是完整二阶张量但描述“施加什么应变”要用 Voigt 向量然后再还原成晶格变形矩阵。下面这个 Python 函数做的事就是这个还原import numpy as np def voigt_to_matrix(voigt): 输入 Voigt 应变向量 (e1,e2,e3,e4,e5,e6)返回对称应变张量矩阵 e1, e2, e3, e4, e5, e6 voigt return np.array([ [e1, e6/2, e5/2], [e6/2, e2, e4/2], [e5/2, e4/2, e3 ] ])注意主对角线就是正应变非对角元要除以 2因为矩阵对称化后张量分量 εxy ½·γxy。很多人第一次从这个函数转到 VASP 的 POSCAR 变形矩阵时少除了这个 2导致施加的剪切应变翻倍拟合出来的 C44 直接两倍偏差。2.3 应变模式与独立弹性常数数量不同晶系需要计算的应变模式数量不一样。立方晶系只有 3 个独立弹性常数 C11, C12, C44但施加同一个等体积应变的一个矩阵并不总能同时分离出所有分量。常用做法是施加三组不同模式应变模式Voigt 向量拟合目标单轴拉伸(δ, 0, 0, 0, 0, 0)C11双轴拉伸(δ, δ, 0, 0, 0, 0)C11 C12剪切(0, 0, 0, δ, 0, 0) 或 (0,0,0,0,0,δ)C44三组应变得到三个能量-应变曲线联立求解就得到三个独立常数。对于六角晶系需要五个独立常数对应至少 5 组应变模式正交晶系是 9 个三斜晶系则要 21 个。常见自动化工具会自动枚举这些模式手动做的时候先用晶系对称性简化。实际操作里每组应变取 59 个幅度从 -2% 到 2%步长 0.5%。少于 5 个点时二次拟合的自由度不够多于 9 个点则计算量翻倍而精度提升有限。2.4 线性应变与工程应变的选择对同一个变形用“线性应变” ε (a-a₀)/a₀ 和“工程应变” η (a²-a₀²)/(2a₀²) 会得到不同的能量解析形式。第一性原理代码里改变晶格常数时POSCAR 或 QE 的 celldm 参数直接按比例缩放所以默认对应的是线性应变。小应变范围内两者只差高阶项但拟合时如果用了工程应变而没用对应的体积变化公式会引入可察觉的系统误差。推荐统一用线性应变因为它直接对应代码里晶格常数缩放系数。设平衡晶格矢为 h₀3×3 矩阵施加应变 ε 后新晶格矢是h h₀ · (I ε)Python 实现def apply_strain(lattice_matrix, voigt): lattice_matrix: 3x3 晶格矢量行为 a,b,c voigt: Voigt 应变向量 strain voigt_to_matrix(voigt) deformed lattice_matrix (np.eye(3) strain) return deformed注意矩阵乘法方向POSCAR 里晶格矢量是行向量还是列向量取决于代码约定。VASP 的 POSCAR 第二行是缩放系数第三到第五行是晶格矢量每行一个矢量上面对应的变形矩阵乘法是 h_new h_old · (Iε)。有些脚本用列向量写法结果相同但要小心 numpy 的 方向颠倒带来的转置问题。3. VASP 实操从 POSCAR 变形到 OUTCAR 应力张量3.1 完整流程概览VASP 下计算弹性常数的标准流程是高精度结构优化得到平衡晶格常数和原子坐标固定晶格对每个应变模式生成多个变形 POSCAR对每个变形结构做静态自洽计算IALGO 默认的电子步收敛从 OUTCAR 或 OSZICAR 提取总能拟合能量-应变曲线得到弹性常数其中第一步和第三步要分开不能直接拿优化未收敛的结构施加应变否则初始应力项不为零能量展开式里的线性项会污染二次拟合。3.2 用 Python 生成 VASP 变形 POSCARimport numpy as np def read_poscar(filename): 读取 POSCAR 晶格矢量部分返回 3x3 数组 with open(filename) 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 write_poscar(filename, lattice, base_lines): 根据原 POSCAR 的元素信息写出新晶格 with open(filename, w) as f: for i, line in enumerate(base_lines[:2]): f.write(line) f.write(f {lattice[0][0]:.8f} {lattice[0][1]:.8f} {lattice[0][2]:.8f}\n) f.write(f {lattice[1][0]:.8f} {lattice[1][1]:.8f} {lattice[1][2]:.8f}\n) f.write(f {lattice[2][0]:.8f} {lattice[2][1]:.8f} {lattice[2][2]:.8f}\n) for line in base_lines[5:]: f.write(line) # 主逻辑读取平衡结构施加单轴应变批量写出 base_poscar open(POSCAR_opt).readlines() h0 read_poscar(POSCAR_opt) strains [x/1000 for x in range(-20, 21, 5)] # -2% 到 2% for i, delta in enumerate(strains): deformed apply_strain(h0, (delta, 0, 0, 0, 0, 0)) write_poscar(fPOSCAR_e1_{i:02d}, deformed, base_poscar)读入平衡结构后把缩放系数并入射格矢量里避免后续处理把缩放系数弄丢。输出文件名建议带应变方向和序号方便后续脚本正则匹配。3.3 INCAR 与 KPOINTS 参数选择弹性常数计算对收敛精度要求高于普通结构优化。推荐一组可复用的 INCARENCUT 1.3 * ENMAX # 平面波截断能VASP 推荐值的 1.3 倍 EDIFF 1E-8 # 电子步自洽收敛标准 EDIFFG -0.01 # 只做静态计算时此项不影响结果 ISMEAR -5 # 绝缘体/半导体用四面体方法 SIGMA 0.2 # ISMEAR0 时才有意义-5 时忽略 PREC Accurate IBRION -1 # 静态计算 NSW 0 LREAL .FALSE. # 精确计算禁用实空间投影对于金属体系ISMEAR 用 0 或 1 配合 SIGMA0.1四面体方法对部分占据不适用。ENCUT 检查办法是用平衡结构分别跑 1.0×ENMAX 和 1.3×ENMAX能量差小于 1 meV/atom 就够用。KPOINTS 至少要用 Monkhorst-Pack 自动网格密度参考cubic 晶系 k-spacing ≤ 0.03 Å⁻¹非立方晶系相应减少。三斜晶系或者超胞大的体系用 Gamma-centered 网格往往更稳。关键检查是增加 k 点数后 Cij 变化小于 1 GPa。3.4 批量提交计算并从 OUTCAR 提取数据生成完一组 POSCAR 后批量复制 INCAR、KPOINTS、POSCAR 和 POTCAR 到子目录然后提交 VASPfor d in strain_*; do cp INCAR KPOINTS POTCAR $d/ cd $d mpirun -np 16 vasp_std run.log cd .. done运行完后提取总能量for d in strain_*; do E$(grep free energy $d/OUTCAR | tail -1 | awk {print $5}) echo $d $E done注意 VASP 输出的能量有多种标签free energy是有限温度自由能energy(sigma-0)是外推的 0K 能量。静态计算用后者更准确因为 ISMEAR-5 时free energy和energy(sigma-0)数值不同。如果发现某个应变后结构对外加应变有明显应力响应再看 OUTCAR 的total stress块确认应变后残余应力与外加应变的对应关系。这能帮助判断 POSCAR 变形矩阵是否方向搞反。def parse_outcar_energy(filename): with open(filename) as f: for line in f: if energy(sigma-0) in line: return float(line.split()[1].split()[0])提取结果后直接进入第五章的拟合脚本。但先看 QE 那边怎么做两者组合才能互为对照。4. QE 实操用 pw.x 的 celldm 参数扫描应变4.1 从 VASP 到 QE 的输入结构转换QE 的 pw.x 输入文件以SYSTEMnamelist 里的晶格参数描述晶体。常用ibrav指定 Bravais 晶格类型celldm或a, b, c, cosAB等指定晶格常量。和 VASP 直接给 3×3 晶格矢量不同QE 用ibrav0时可以直接给CELL_PARAMETERS三行矢量这和 VASP 的 POSCAR 几乎对应。CONTROL calculation scf prefix si outdir ./tmp pseudo_dir ./pseudo/ verbosity high / SYSTEM ibrav 0 nat 2 ntyp 1 ecutwfc 40 ecutrho 320 occupationssmearing smearingmp degauss0.01 / ELECTRONS conv_thr 1.0d-9 mixing_beta 0.4 / ATOMIC_SPECIES Si 28.0855 Si.pz-vbc.UPF CELL_PARAMETERS {angstrom} 3.8400 0.0000 0.0000 1.9200 3.3260 0.0000 1.9200 1.1087 3.1356 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 0ecutwfc对应 VASP 的 ENCUTecutrho默认取 812 倍ecutwfc超软赝势需要更高。QE 里用degauss做 smearing此时平面波截断能对弹性常数的影响比 VASP 更明显务必做截断能收敛测试。4.2 用 Python 修改 CELL_PARAMETERS 施加应变等价于 VASP 的 POSCAR 变形直接操作 QE 的 CELL_PARAMETERS 块import re def apply_strain_qe(input_text, voigt): 读入 QE 输入文件字符串修改 CELL_PARAMETERS 晶格矢量 lines input_text.splitlines() for i, line in enumerate(lines): name line.split()[0].upper() if name CELL_PARAMETERS: cells [] for j in range(i1, i4): cells.append([float(x) for x in lines[j].split()]) h0 np.array(cells) deformed apply_strain(h0, voigt) for j in range(3): lines[i1j] f{deformed[j][0]:.8f} {deformed[j][1]:.8f} {deformed[j][2]:.8f} break return \n.join(lines)然后用循环生成多个输入文件for i, delta in enumerate(strains): text open(scf.in).read() deformed_text apply_strain_qe(text, (delta, 0, 0, 0, 0, 0)) open(fscf_e1_{i:02d}.in, w).write(deformed_text)注意 QE 的CELL_PARAMETERS默认单位是 alat如果在SYSTEM里写了ibrav0且没有celldm(1)单位默认为 Å规范起见在CELL_PARAMETERS那行显式写{angstrom}避免单位歧义。4.3 批量运行与数据提取for f in scf_e1_*.in; do mpirun -np 8 pw.x -in $f ${f%.in}.out done从输出文件提取能量和压力grep ! scf_e1_*.out grep P scf_e1_*.out | tail -1QE 的P 行给出外加流体静压力单位是 kbar。这个值在等体积应变下应为 0如果参考结构是平衡态。如果P的数值明显漂移说明参考结构没有充分优化。QE 没有像 VASP OUTCAR 那样方便的总能量提取energy(sigma-0)但!行给出的total energy在固定占据数下就是最终能量。金属体系用抹布时能量对 degauss 有依赖所以每组应变都用同一套 smearing 参数避免系统性偏差。还可以从PW的total stress (GPa)输出块读取完整应力张量单位已经转换成 GPa方便和 VASP 对比。def parse_qe_energy(filename): with open(filename) as f: lines f.readlines() for line in reversed(lines): if line.startswith(!): return float(line.split()[-2])这个函数从文件尾部反向搜索!行得到总能。注意 QE 输出里!行格式是! total energy ... Ry用负索引提取能正确处理多行输出。4.4 VASP 和 QE 数据格式对比项目VASPQE晶格描述POSCAR 三行矢量CELL_PARAMETERS 三行矢量截断能控制ENCUT单位 eVecutwfc/ecutrho单位 Ry能量单位eVRy1 Ry 13.6057 eV应力单位OUTCAR 里 kB输出里 GPa 或 kbar施加载荷后的计算必须重新做自洽单次 scf 即可处理剪切应变直接变形直接变形关键坑位是两个代码能量单位差了 13.6 倍拟合弹性常数时要统一成 eV 或者 J/m³不然 Cij 数值天然差两个数量级。5. Python 后处理拟合弹性常数并画出应力-应变线性段5.1 二次多项式拟合与 Cij 提取把 VASP 和 QE 提取的能量放到一起按应变值排序然后用 numpy 做最小二乘二次拟合import numpy as np def fit_elastic(data, strain_name, volume): data: [(strain, energy_eV), ...] volume: 平衡体积单位 Å^3 返回 (C_value_GPa, 拟合误差) strain np.array([d[0] for d in data]) energy np.array([d[1] for d in data]) coefs np.polyfit(strain, energy, 2) # 能量单位 eV体积单位 Å^3转 GPa 的系数 eV_per_Ang3_to_GPa 160.21766208 C 2 * coefs[0] / volume * eV_per_Ang3_to_GPa # 手动算不确定度 fit_e energy - np.polyval(coefs, strain) sigma np.sqrt(np.mean(fit_e**2)) error sigma / np.sqrt(len(strain)) * 2 * eV_per_Ang3_to_GPa / volume return C, error二次项系数乘 2 再除以体积才是弹性常数这个 2 来自能量展开的 ½ 因子。注意体积单位换算VASP 的 CONTCAR 里体积单位是 ųQE 的 CELL_PARAMETERS 如果用 Å 则一致如果用 alat 则需要乘 (celldm(1)*0.529177)³。np.polyfit默认返回最高次项在前所以coefs[0]是二次项系数。5.2 用应力-应变曲线的线性拟合交叉验证能量法算出弹性常数后可以再从应力-应变数据验证。VASP 的 OUTCAR 里total stress是应力张量用输出的应变对应的应力分量 σ 对 ε 做线性拟合斜率应该等于对应 Cijimport numpy as np stress_eps [] for i, delta in enumerate(strains): stress read_outcar_stress(fstrain_e1_{i:02d}/OUTCAR) stress_eps.append((delta, stress[0, 0])) # xx 分量 # stress 单位 kB1 kB 0.1 GPa stress_eps np.array(stress_eps) stress_eps[:, 1] * 0.1 # 转 GPa coef_slope np.polyfit(stress_eps[:, 0], stress_eps[:, 1], 1) print(fC11 from stress-strain slope: {coef_slope[0]:.1f} GPa)这个数值和能量法拟合出来的一致性应落在数值噪声范围内一般是 12 GPa 的偏差。如果偏差达到 5 GPa 以上首先检查应变定义应力是共轭于有限应变还是工程应变VASP 的应力输出对应的是无限小应变张量所以和线性应变一致。5.3 批量脚本与结果汇总下面把整个流程的可复用脚本框架画出来实际使用时按代码包里的目录结构调整路径import json results {} for code in [vasp, qe]: results[code] {} for mode in [uniaxial, biaxial, shear]: data collect_data(code, mode) C, err fit_elastic(data, volume, mode) results[code][mode] C with open(Cij_results.json, w) as f: json.dump(results, f, indent2)格式化成 JSON 的好处是后续不管喂给 origin 画图还是喂给机器学习模型都不用重新解析文本。5.4 画能量-应变曲线把拟合曲线和原始数据画在一张图里一目了然import matplotlib.pyplot as plt strain_points [d[0] for d in data] energy_points [d[1] for d in data] sm np.linspace(min(strain_points), max(strain_points), 50) plt.scatter(strain_points, energy_points, labelDFT data) plt.plot(sm, np.polyval(coefs, sm), r-, labelquadratic fit) plt.xlabel(Strain) plt.ylabel(Energy (eV)) plt.legend() plt.savefig(f{code}_{mode}_fit.png, dpi150)拟合曲线如果偏离数据点超过 1 meV说明应变范围取得太大非线性效应进入了三次项范围这时把应变范围收窄到 ±1%或者拟合时加三阶项看系数是否显著。正常材料在 2% 应变内三次项贡献小于 0.1%但如果体系本身有相变比如铁电材料要缩短范围。6. 验证弹性常数结果的四个收敛性检查弹性常数的精确值对计算参数极其敏感跑完一组数据不做检查直接发表很容易出问题。以下四个检查按顺序做任何一环不过都要返工。正式计算前先跑一组单点收敛测试固定 k 点从 1.0×ENMAX 扫到 1.3×ENMAX记录能量差再固定 ENCUT从 ×1 到 ×2 加密 k 点网格。如果截断能每提升 5% 能量变化超过 1 meV/atom弹性常数变化会超过 2 GPa。QE 同理ecutwfc 从 30 扫到 60 Ry记录 Cij 变化趋势。合格标准是相邻两个截断能下 Cij 差小于 0.5 GPa。应变幅度检查看二次拟合残差把残差对应变画出来如果呈 U 型分布说明二阶模型不够加个三次项试试如果残差呈现周期波动多半是 SCF 收敛精度不够EDIFF 或 conv_thr 还需要收紧一两个量级。参考结构必须是弛豫到压力为 0 的平衡晶格。判断方法是应变0 的那个点VASP 的静水压力小于 0.5 kbarQE 的 P 小于 0.5 kbar。如果初始压力太大拟合斜率和二次项都会失真弹性常数可能偏小约 510%。最终要把算出的 Cij 组成对称矩阵检查特征值是否全部为正对应 Born 力学稳定性判据。立方晶系的判据是 C44 0C11 |C12|C112C12 0。如果 C44 算出负值优先怀疑剪切应变定义里的 γ/2 因子没处理好其次检查 k 点密度是否太低。交叉验证用 VASP 和 QE 两套结果对比偏差大于 3 GPa 时排查赝势和截断能设置。用 0.5% 步长 5 个点扫一遍如果 Cij 随步长变化超过 1 GPa说明应变区间选得不合适重新选步长。本文还有配套的精品资源点击获取