尧图网页设计 
简介这份资源面向材料科学计算方向的研究生与科研人员聚焦如何用Python驱动VASP与Quantum Espresso完成应力—应变关系计算解决第一性原理力学性质模拟中数据提取、处理与可视化的问题。压缩包共16个文件约30KB以8个Python脚本为核心配合4个输入文件、POSCAR结构文件及README说明覆盖拉伸与剪切两类计算场景并区分VASP与QE两套流程脚本还提供是否绘图的可选版本便于按需调用。已有931人学习下载说明其在同类计算任务中具备一定参考价值。读者可从中获得读取输出文件、提取应力应变数据、绘制曲线并进一步拟合弹性模量与泊松比等参数的完整脚本框架同时借助示例输入与结构文件快速复现计算流程理解两种DFT软件在应变模拟中的衔接方式适合作为力学性质计算的入门模板与二次开发起点。1. 从一份能跑的应力应变脚本包说起VASP 与 QE 双引擎的 Python 落地很多人第一次算弹性常数卡住的不是 DFT 本身而是「应变怎么加、应力从哪读、数据怎么对齐」。VASP 和 Quantum Espresso 各自输出格式不同一个走 OUTCAR 里的应力张量一个走 pwscf 的 output 或 XML手动抄数据基本等于自找麻烦。这个压缩包 StrengthCalculation-strain-stress-master 把拉伸和剪切两条线都拆成了独立脚本VASP 和 QE 各一套还配了 diamond 的 POSCAR、relax 输入和旋转版本等于把「改晶格 → 跑静态 → 提应力 → 拟合」整条链路摊开给你看。适合已经装好 VASP 或 QE、能跑通单点能计算、但还没把弹性常数流程串起来的材料计算从业者。Python 在这里不是主角是胶水负责把两个引擎的输出粘成一条应力应变曲线。2. 脚本包结构拆解拉伸与剪切两条线怎么分2.1 文件命名里的信息量拿到压缩包先别急着跑把文件名过一遍就能看出作者的意图。tensile_calculation_withoutplot_vasp.py和tensile_calculation_plotcheck_vasp.py是一对前者只算不画后者带绘图检查QE 侧同样有tensile_calculation_withoutplot_qe.py和tensile_calculation_plotcheck_qe.py。剪切线是shear_calculation_plotcheck_vasp.py、shear_calculation_withoutplot_vasp.py、shear_calculation_plotcheck_qe.py、shear_calculation_withoutplot_qe.py。这种「withoutplot / plotcheck」的拆分很实用批量跑的时候用 withoutplot 省时间调参阶段用 plotcheck 当场看曲线是否线性。结构文件有POSCAR和POSCAR_rota对应relax.in、relax_rota.in、diamond.relax.in、diamond.relax_rota.in。rota后缀说明作者考虑了晶格取向问题——金刚石结构在不同晶向下弹性响应不同旋转后的 POSCAR 用来验证各向异性。README.md 是唯一的文档入口.gitattributes说明这包是从 git 仓库导出的版本管理痕迹还在。2.2 拉伸与剪切的物理区别在脚本里怎么体现拉伸计算改的是晶格常数沿某个方向拉长或压缩其他方向可能固定也可能按泊松比松弛。剪切计算改的是晶格矢量之间的夹角POSCAR 里表现为基矢的非对角项变化。脚本里对这两种形变的处理逻辑不同拉伸通常只动一个晶格参数剪切要构造完整的形变矩阵。常见做法是定义一个形变梯度矩阵对原始晶格矢量做线性变换。拉伸对应对角矩阵剪切对应非对角元非零的矩阵。脚本里如果直接改 POSCAR 的缩放系数那只适合各向同性拉伸要算完整的弹性常数矩阵得按 Voigt 记号逐个施加应变模式。提示先确认脚本里应变的定义是工程应变还是真应变。小应变下两者差别不大但应变加到 2% 以上时拟合出的弹性常数会有可观测的偏差。2.3 输入文件与脚本的对应关系relax.in是 QE 的输入模板diamond.relax.in是金刚石的具体算例。VASP 侧没有单独的 INCAR 模板说明脚本可能直接生成 INCAR 或者依赖你手动准备。POSCAR 是 VASP 的结构文件QE 用CELL_PARAMETERS和ATOMIC_POSITIONS卡片两者格式不通用脚本里应该有转换逻辑或者分别读取的分支。跑之前建议先手动跑一个应变点确认 VASP 或 QE 能正常输出应力。VASP 看 OUTCAR 里in kB那行的应力张量QE 看 output 里total stress段落。如果这一步就报错后面脚本再对也没用。3. 环境准备与单点验证跑脚本前必须过的三道关3.1 VASP 与 QE 的最小可跑配置VASP 需要 POTCAR、POSCAR、INCAR、KPOINTS 四个文件。POTCAR 按元素顺序拼接POSCAR 用包里的INCAR 至少设IBRION -1静态计算、ISIF 2算应力但不改结构、NSW 0。KPOINTS 用 Gamma 点或 Monkhorst-Pack 网格金刚石结构用 8×8×8 起步比较稳。QE 的输入文件是单个.in里面包含CONTROL、SYSTEM、ELECTRONS三个 namelist后面跟ATOMIC_SPECIES、CELL_PARAMETERS、ATOMIC_POSITIONS、K_POINTS卡片。relax.in里calculation scf还是relax决定了跑的是单点还是结构优化。算应力必须用tprnfor .true.和tprnstr .true.否则 output 里没有应力张量。# VASP 单点验证确认 OUTCAR 里有应力输出 mpirun -np 4 vasp_std vasp.log grep in kB OUTCAR # 应该看到类似in kB -1.234 0.567 0.567 0.000 0.000 0.000# QE 单点验证确认 output 里有 total stress mpirun -np 4 pw.x -in diamond.relax.in qe.log grep -A 3 total stress qe.log # 应该看到 3x3 应力张量单位是 Ry/bohr^3 或 GPaVASP 的应力单位默认是 kB1 kB 0.1 GPa换算时别搞错。QE 的应力单位取决于CONTROL里的设置常见做法是在后处理时统一转成 GPa。3.2 Python 依赖与 pymatgen 的取舍脚本大概率依赖 numpy 和 matplotlib可能还用了 pymatgen 读 OUTCAR。pymatgen 的Outcar类能直接提取应力张量省去手写正则的麻烦。但 pymatgen 安装有时会卡在依赖上如果只是读几个数用正则匹配 OUTCAR 更轻量。# 轻量读取 VASP OUTCAR 应力张量的常见写法 import re import numpy as np def read_stress_outcar(filename): 从 OUTCAR 提取最后一个应力张量返回 3x3 numpy 数组 with open(filename, r) as f: lines f.readlines() # 从后往前找取最后一个 in kB 后面的数据 for i in range(len(lines) - 1, -1, -1): if in kB in lines[i]: # 应力值在下一行6 个分量xx yy zz xy yz zx vals [float(x) for x in lines[i 1].split()] # 组装成 3x3 对称矩阵 stress np.array([ [vals[0], vals[3], vals[5]], [vals[3], vals[1], vals[4]], [vals[5], vals[4], vals[2]] ]) return stress * 0.1 # kB 转 GPa raise ValueError(OUTCAR 里没找到应力数据)这段代码的关键在应力分量的顺序。VASP 输出的是 Voigt 记号xx, yy, zz, xy, yz, zx。组装成 3×3 矩阵时xy 对应 (0,1) 和 (1,0)yz 对应 (1,2) 和 (2,1)zx 对应 (2,0) 和 (0,2)。顺序搞反了剪切分量会错位拟合出的 C44 会偏。3.3 应变步长的选择与收敛测试应变步长不是越小越好。步长太小比如 0.1%应力响应可能被数值噪声淹没步长太大比如 5%超出弹性范围曲线弯曲拟合出的弹性常数偏小。常见做法是在 ±2% 范围内取 5 到 7 个点步长 0.5% 或 1%。每个应变点都要重新跑静态计算K 点网格和截断能要和结构优化时一致。如果 relax 用的 8×8×8静态也得用 8×8×8否则能量和应力不可比。收敛标准EDIFF建议设到 1E-6 甚至 1E-7应力对电子步收敛更敏感。注意QE 的ecutwfc和ecutrho在应变计算中要保持不变。有人为了省时间在静态计算里降截断能结果应力张量整体偏移弹性常数全错。4. 拉伸与剪切脚本实操从改 POSCAR 到拟合弹性常数4.1 拉伸计算改哪个晶格参数、怎么改拉伸计算的核心是构造一系列形变后的结构。以金刚石为例原始 POSCAR 的晶格矢量是三行缩放系数在第一行。沿 z 方向拉伸就是把第三行矢量乘以 (1ε)其他两行不变。脚本里如果直接改缩放系数那是各向同性拉伸只能算体弹模量。# 沿指定方向施加单轴应变生成一系列 POSCAR import numpy as np def generate_strained_poscar(base_lattice, strain_list, direction2): base_lattice: 3x3 晶格矢量矩阵每行一个矢量 strain_list: 应变值列表如 [-0.02, -0.01, 0, 0.01, 0.02] direction: 0/1/2 对应 x/y/z 方向 返回每个应变对应的晶格矩阵列表 strained [] for eps in strain_list: new_lat base_lattice.copy() # 只改指定方向的晶格矢量长度 new_lat[direction] base_lattice[direction] * (1 eps) strained.append(new_lat) return strained # 读取原始 POSCAR 的晶格部分 def read_poscar_lattice(filename): with open(filename) as f: lines f.readlines() scale float(lines[1].strip()) lattice np.array([[float(x) for x in lines[i].split()] for i in range(2, 5)]) * scale return lattice这段代码只做了单轴拉伸。要算完整的弹性常数矩阵需要分别沿 x、y、z 拉伸再施加 xy、yz、zx 三个剪切模式。每个模式对应一列弹性常数。脚本包里的tensile_calculation_*应该覆盖了单轴拉伸shear_calculation_*覆盖剪切。参数说明strain_list的范围建议 ±0.02点数 5 到 7 个。direction参数在单轴拉伸时只改一个方向但实际计算中其他方向可能因为泊松效应产生应力所以 ISIF 要设成 2允许原子弛豫但不改晶胞或者设成 4允许改晶胞形状但体积不变——具体看你要算的是哪个弹性常数。4.2 剪切计算形变矩阵的构造剪切比拉伸麻烦因为要改的是晶格矢量之间的夹角。以 xy 剪切为例形变矩阵是单位矩阵加上一个非对角元# 构造剪切形变矩阵并应用到晶格 def apply_shear(lattice, eps, planexy): 对晶格施加剪切应变eps 是剪切量 deform np.eye(3) if plane xy: deform[0, 1] eps # x 方向矢量在 y 方向的分量 elif plane yz: deform[1, 2] eps elif plane zx: deform[2, 0] eps # 形变后的晶格 原始晶格 形变矩阵的转置 return lattice deform.T剪切应变下晶胞体积基本不变但对称性降低。VASP 的 ISIF 要设成 2让原子在固定晶胞内弛豫。QE 里calculation scf配合tprnstr .true.即可。剪切计算最容易翻车的地方是应变方向。xy 剪切和 yx 剪切在弹性常数矩阵里是同一个分量但形变矩阵的写法不同。如果脚本里deform[0,1]和deform[1,0]混用算出的 C44 和 C55 会对调。建议先跑一个已知材料验证比如金刚石的 C44 约 580 GPa算出来差太多就是方向搞错了。4.3 应力提取与弹性常数拟合拿到一系列应变和对应的应力后拟合就是线性回归。弹性常数 Cij dσi / dεj在弹性范围内应力应变是线性的。用 numpy 的 polyfit 一次多项式即可。# 从应力-应变数据拟合弹性常数 import numpy as np def fit_elastic_constant(strain_list, stress_list): strain_list: 应变值列表 stress_list: 对应应力分量列表GPa 返回弹性常数GPa和拟合优度 R^2 coeffs np.polyfit(strain_list, stress_list, 1) slope coeffs[0] # 弹性常数 # 计算 R^2 p np.poly1d(coeffs) yhat p(strain_list) ybar np.mean(stress_list) ss_res np.sum((stress_list - yhat) ** 2) ss_tot np.sum((stress_list - ybar) ** 2) r2 1 - ss_res / ss_tot return slope, r2 # 示例单轴拉伸沿 z 方向取 σ_zz 对 ε_zz 的斜率 strains [-0.02, -0.01, 0.0, 0.01, 0.02] stresses [-12.5, -6.3, 0.1, 6.2, 12.4] # 单位 GPa C33, r2 fit_elastic_constant(strains, stresses) print(fC33 {C33:.1f} GPa, R^2 {r2:.4f})R² 低于 0.99 就要检查应变范围是否太大、某个应变点的计算是否没收敛、应力提取是否取错了行。金刚石这类高对称材料线性应该非常好R² 接近 1。提示拟合时截距应该接近零。如果截距明显偏离零说明零应变点的结构没有完全弛豫或者应力提取有系统误差。5. 避坑与排查应力应变计算里最容易翻车的五件事5.1 现象应力张量全是零原因VASP 的 INCAR 里没设ISIF 2或IBRION -1或者 QE 的tprnstr没打开。VASP 默认 ISIF2 会算应力但如果 NSW0 且 IBRION-1应力应该正常输出。QE 的tprnstr默认是.false.必须显式设成.true.。解决检查 INCAR 和 QE 输入文件确认应力输出开关打开。VASP 还可以看 OUTCAR 里有没有FORCES acting on ions和Stress tensor段落。5.2 现象不同应变点的能量不连续原因K 点网格或截断能在不同应变点之间变了。有人为了省时间在小应变时用低截断大应变时用高截断导致能量基准不一致。解决所有应变点用完全相同的计算参数。写脚本时把 KPOINTS 和 INCAR 模板固定只改 POSCAR 或 CELL_PARAMETERS。5.3 现象拟合出的弹性常数偏小原因应变范围太大超出了线性弹性区。金刚石的弹性线性范围大概在 ±1% 以内加到 ±3% 曲线就弯了。解决把应变范围缩到 ±1%或者用二次多项式拟合后取一次项系数。更稳妥的做法是先跑一个大范围看曲线拐点再在拐点内取点。5.4 现象QE 和 VASP 算出的弹性常数差很多原因两者的赝势不同、交换关联泛函可能不同、应力单位换算有误。VASP 的 kB 转 GPa 是乘 0.1QE 的 Ry/bohr³ 转 GPa 要乘 14710.5。解决先统一泛函都用 PBE和赝势类型都用 PAW 或都用超软再核对单位换算。如果还差检查 QE 的ecutwfc是否足够QE 对截断能比 VASP 敏感。5.5 现象脚本跑完没有输出文件原因VASP 或 QE 的可执行文件路径没设对或者 mpirun 的进程数和 K 点并行不匹配。脚本里如果硬编码了vasp_std或pw.x环境变量 PATH 里没有就会静默失败。解决在脚本里加错误检查跑完检查 OUTCAR 或 output 文件是否存在且非空。常见做法是用subprocess.run的checkTrue让异常抛出而不是默默跳过。6. 进阶技巧用 plotcheck 脚本做快速验证与批量扫描plotcheck版本的脚本价值在于边算边看。我一般会先用它跑一个应变点确认应力提取和绘图链路通了再切到withoutplot批量跑。批量跑的时候把应变列表写成循环每个应变生成一个目录跑完统一收集数据。# 批量扫描应变并收集应力适合 withoutplot 模式 import os import subprocess import numpy as np def batch_strain_scan(base_dir, strain_list, direction2): 在每个应变下生成结构、跑 VASP、收集应力 results [] for eps in strain_list: work_dir os.path.join(base_dir, fstrain_{eps:.4f}) os.makedirs(work_dir, exist_okTrue) # 生成 POSCAR省略具体实现参考 4.1 节 # 复制 INCAR、KPOINTS、POTCAR 到 work_dir # 跑 VASP subprocess.run(mpirun -np 4 vasp_std vasp.log, shellTrue, cwdwork_dir, checkTrue) # 提取应力 stress read_stress_outcar(os.path.join(work_dir, OUTCAR)) results.append((eps, stress[2, 2])) # 取 σ_zz return results # 跑完直接拟合 strains, stresses zip(*batch_strain_scan(./scan, [-0.01, -0.005, 0, 0.005, 0.01])) C33, r2 fit_elastic_constant(list(strains), list(stresses)) print(fC33 {C33:.1f} GPa, R^2 {r2:.4f})这个批量脚本的关键是每个应变独立目录避免文件覆盖。checkTrue保证 VASP 报错时脚本停下来而不是继续跑下一个应变点。收集完数据直接拟合R² 不合格就回去查哪个应变点出了问题。验证方法上我习惯拿金刚石或硅做基准。金刚石的 C11 约 1076 GPa、C12 约 125 GPa、C44 约 577 GPa算出来在这个量级附近就说明流程通了。如果差一个数量级多半是单位换算错了如果差 20% 到 30%检查赝势和截断能。注意QE 的应力输出在 XML 文件里更规整pw.x的 text output 有时会截断。如果脚本读 text output 不稳定改用style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;" /> VASP与QE双引擎Python脚本：应力应变计算与弹性常数拟合实战
简介这份资源面向材料科学计算方向的研究生与科研人员聚焦如何用Python驱动VASP与Quantum Espresso完成应力—应变关系计算解决第一性原理力学性质模拟中数据提取、处理与可视化的问题。压缩包共16个文件约30KB以8个Python脚本为核心配合4个输入文件、POSCAR结构文件及README说明覆盖拉伸与剪切两类计算场景并区分VASP与QE两套流程脚本还提供是否绘图的可选版本便于按需调用。已有931人学习下载说明其在同类计算任务中具备一定参考价值。读者可从中获得读取输出文件、提取应力应变数据、绘制曲线并进一步拟合弹性模量与泊松比等参数的完整脚本框架同时借助示例输入与结构文件快速复现计算流程理解两种DFT软件在应变模拟中的衔接方式适合作为力学性质计算的入门模板与二次开发起点。1. 从一份能跑的应力应变脚本包说起VASP 与 QE 双引擎的 Python 落地很多人第一次算弹性常数卡住的不是 DFT 本身而是「应变怎么加、应力从哪读、数据怎么对齐」。VASP 和 Quantum Espresso 各自输出格式不同一个走 OUTCAR 里的应力张量一个走 pwscf 的 output 或 XML手动抄数据基本等于自找麻烦。这个压缩包 StrengthCalculation-strain-stress-master 把拉伸和剪切两条线都拆成了独立脚本VASP 和 QE 各一套还配了 diamond 的 POSCAR、relax 输入和旋转版本等于把「改晶格 → 跑静态 → 提应力 → 拟合」整条链路摊开给你看。适合已经装好 VASP 或 QE、能跑通单点能计算、但还没把弹性常数流程串起来的材料计算从业者。Python 在这里不是主角是胶水负责把两个引擎的输出粘成一条应力应变曲线。2. 脚本包结构拆解拉伸与剪切两条线怎么分2.1 文件命名里的信息量拿到压缩包先别急着跑把文件名过一遍就能看出作者的意图。tensile_calculation_withoutplot_vasp.py和tensile_calculation_plotcheck_vasp.py是一对前者只算不画后者带绘图检查QE 侧同样有tensile_calculation_withoutplot_qe.py和tensile_calculation_plotcheck_qe.py。剪切线是shear_calculation_plotcheck_vasp.py、shear_calculation_withoutplot_vasp.py、shear_calculation_plotcheck_qe.py、shear_calculation_withoutplot_qe.py。这种「withoutplot / plotcheck」的拆分很实用批量跑的时候用 withoutplot 省时间调参阶段用 plotcheck 当场看曲线是否线性。结构文件有POSCAR和POSCAR_rota对应relax.in、relax_rota.in、diamond.relax.in、diamond.relax_rota.in。rota后缀说明作者考虑了晶格取向问题——金刚石结构在不同晶向下弹性响应不同旋转后的 POSCAR 用来验证各向异性。README.md 是唯一的文档入口.gitattributes说明这包是从 git 仓库导出的版本管理痕迹还在。2.2 拉伸与剪切的物理区别在脚本里怎么体现拉伸计算改的是晶格常数沿某个方向拉长或压缩其他方向可能固定也可能按泊松比松弛。剪切计算改的是晶格矢量之间的夹角POSCAR 里表现为基矢的非对角项变化。脚本里对这两种形变的处理逻辑不同拉伸通常只动一个晶格参数剪切要构造完整的形变矩阵。常见做法是定义一个形变梯度矩阵对原始晶格矢量做线性变换。拉伸对应对角矩阵剪切对应非对角元非零的矩阵。脚本里如果直接改 POSCAR 的缩放系数那只适合各向同性拉伸要算完整的弹性常数矩阵得按 Voigt 记号逐个施加应变模式。提示先确认脚本里应变的定义是工程应变还是真应变。小应变下两者差别不大但应变加到 2% 以上时拟合出的弹性常数会有可观测的偏差。2.3 输入文件与脚本的对应关系relax.in是 QE 的输入模板diamond.relax.in是金刚石的具体算例。VASP 侧没有单独的 INCAR 模板说明脚本可能直接生成 INCAR 或者依赖你手动准备。POSCAR 是 VASP 的结构文件QE 用CELL_PARAMETERS和ATOMIC_POSITIONS卡片两者格式不通用脚本里应该有转换逻辑或者分别读取的分支。跑之前建议先手动跑一个应变点确认 VASP 或 QE 能正常输出应力。VASP 看 OUTCAR 里in kB那行的应力张量QE 看 output 里total stress段落。如果这一步就报错后面脚本再对也没用。3. 环境准备与单点验证跑脚本前必须过的三道关3.1 VASP 与 QE 的最小可跑配置VASP 需要 POTCAR、POSCAR、INCAR、KPOINTS 四个文件。POTCAR 按元素顺序拼接POSCAR 用包里的INCAR 至少设IBRION -1静态计算、ISIF 2算应力但不改结构、NSW 0。KPOINTS 用 Gamma 点或 Monkhorst-Pack 网格金刚石结构用 8×8×8 起步比较稳。QE 的输入文件是单个.in里面包含CONTROL、SYSTEM、ELECTRONS三个 namelist后面跟ATOMIC_SPECIES、CELL_PARAMETERS、ATOMIC_POSITIONS、K_POINTS卡片。relax.in里calculation scf还是relax决定了跑的是单点还是结构优化。算应力必须用tprnfor .true.和tprnstr .true.否则 output 里没有应力张量。# VASP 单点验证确认 OUTCAR 里有应力输出 mpirun -np 4 vasp_std vasp.log grep in kB OUTCAR # 应该看到类似in kB -1.234 0.567 0.567 0.000 0.000 0.000# QE 单点验证确认 output 里有 total stress mpirun -np 4 pw.x -in diamond.relax.in qe.log grep -A 3 total stress qe.log # 应该看到 3x3 应力张量单位是 Ry/bohr^3 或 GPaVASP 的应力单位默认是 kB1 kB 0.1 GPa换算时别搞错。QE 的应力单位取决于CONTROL里的设置常见做法是在后处理时统一转成 GPa。3.2 Python 依赖与 pymatgen 的取舍脚本大概率依赖 numpy 和 matplotlib可能还用了 pymatgen 读 OUTCAR。pymatgen 的Outcar类能直接提取应力张量省去手写正则的麻烦。但 pymatgen 安装有时会卡在依赖上如果只是读几个数用正则匹配 OUTCAR 更轻量。# 轻量读取 VASP OUTCAR 应力张量的常见写法 import re import numpy as np def read_stress_outcar(filename): 从 OUTCAR 提取最后一个应力张量返回 3x3 numpy 数组 with open(filename, r) as f: lines f.readlines() # 从后往前找取最后一个 in kB 后面的数据 for i in range(len(lines) - 1, -1, -1): if in kB in lines[i]: # 应力值在下一行6 个分量xx yy zz xy yz zx vals [float(x) for x in lines[i 1].split()] # 组装成 3x3 对称矩阵 stress np.array([ [vals[0], vals[3], vals[5]], [vals[3], vals[1], vals[4]], [vals[5], vals[4], vals[2]] ]) return stress * 0.1 # kB 转 GPa raise ValueError(OUTCAR 里没找到应力数据)这段代码的关键在应力分量的顺序。VASP 输出的是 Voigt 记号xx, yy, zz, xy, yz, zx。组装成 3×3 矩阵时xy 对应 (0,1) 和 (1,0)yz 对应 (1,2) 和 (2,1)zx 对应 (2,0) 和 (0,2)。顺序搞反了剪切分量会错位拟合出的 C44 会偏。3.3 应变步长的选择与收敛测试应变步长不是越小越好。步长太小比如 0.1%应力响应可能被数值噪声淹没步长太大比如 5%超出弹性范围曲线弯曲拟合出的弹性常数偏小。常见做法是在 ±2% 范围内取 5 到 7 个点步长 0.5% 或 1%。每个应变点都要重新跑静态计算K 点网格和截断能要和结构优化时一致。如果 relax 用的 8×8×8静态也得用 8×8×8否则能量和应力不可比。收敛标准EDIFF建议设到 1E-6 甚至 1E-7应力对电子步收敛更敏感。注意QE 的ecutwfc和ecutrho在应变计算中要保持不变。有人为了省时间在静态计算里降截断能结果应力张量整体偏移弹性常数全错。4. 拉伸与剪切脚本实操从改 POSCAR 到拟合弹性常数4.1 拉伸计算改哪个晶格参数、怎么改拉伸计算的核心是构造一系列形变后的结构。以金刚石为例原始 POSCAR 的晶格矢量是三行缩放系数在第一行。沿 z 方向拉伸就是把第三行矢量乘以 (1ε)其他两行不变。脚本里如果直接改缩放系数那是各向同性拉伸只能算体弹模量。# 沿指定方向施加单轴应变生成一系列 POSCAR import numpy as np def generate_strained_poscar(base_lattice, strain_list, direction2): base_lattice: 3x3 晶格矢量矩阵每行一个矢量 strain_list: 应变值列表如 [-0.02, -0.01, 0, 0.01, 0.02] direction: 0/1/2 对应 x/y/z 方向 返回每个应变对应的晶格矩阵列表 strained [] for eps in strain_list: new_lat base_lattice.copy() # 只改指定方向的晶格矢量长度 new_lat[direction] base_lattice[direction] * (1 eps) strained.append(new_lat) return strained # 读取原始 POSCAR 的晶格部分 def read_poscar_lattice(filename): with open(filename) as f: lines f.readlines() scale float(lines[1].strip()) lattice np.array([[float(x) for x in lines[i].split()] for i in range(2, 5)]) * scale return lattice这段代码只做了单轴拉伸。要算完整的弹性常数矩阵需要分别沿 x、y、z 拉伸再施加 xy、yz、zx 三个剪切模式。每个模式对应一列弹性常数。脚本包里的tensile_calculation_*应该覆盖了单轴拉伸shear_calculation_*覆盖剪切。参数说明strain_list的范围建议 ±0.02点数 5 到 7 个。direction参数在单轴拉伸时只改一个方向但实际计算中其他方向可能因为泊松效应产生应力所以 ISIF 要设成 2允许原子弛豫但不改晶胞或者设成 4允许改晶胞形状但体积不变——具体看你要算的是哪个弹性常数。4.2 剪切计算形变矩阵的构造剪切比拉伸麻烦因为要改的是晶格矢量之间的夹角。以 xy 剪切为例形变矩阵是单位矩阵加上一个非对角元# 构造剪切形变矩阵并应用到晶格 def apply_shear(lattice, eps, planexy): 对晶格施加剪切应变eps 是剪切量 deform np.eye(3) if plane xy: deform[0, 1] eps # x 方向矢量在 y 方向的分量 elif plane yz: deform[1, 2] eps elif plane zx: deform[2, 0] eps # 形变后的晶格 原始晶格 形变矩阵的转置 return lattice deform.T剪切应变下晶胞体积基本不变但对称性降低。VASP 的 ISIF 要设成 2让原子在固定晶胞内弛豫。QE 里calculation scf配合tprnstr .true.即可。剪切计算最容易翻车的地方是应变方向。xy 剪切和 yx 剪切在弹性常数矩阵里是同一个分量但形变矩阵的写法不同。如果脚本里deform[0,1]和deform[1,0]混用算出的 C44 和 C55 会对调。建议先跑一个已知材料验证比如金刚石的 C44 约 580 GPa算出来差太多就是方向搞错了。4.3 应力提取与弹性常数拟合拿到一系列应变和对应的应力后拟合就是线性回归。弹性常数 Cij dσi / dεj在弹性范围内应力应变是线性的。用 numpy 的 polyfit 一次多项式即可。# 从应力-应变数据拟合弹性常数 import numpy as np def fit_elastic_constant(strain_list, stress_list): strain_list: 应变值列表 stress_list: 对应应力分量列表GPa 返回弹性常数GPa和拟合优度 R^2 coeffs np.polyfit(strain_list, stress_list, 1) slope coeffs[0] # 弹性常数 # 计算 R^2 p np.poly1d(coeffs) yhat p(strain_list) ybar np.mean(stress_list) ss_res np.sum((stress_list - yhat) ** 2) ss_tot np.sum((stress_list - ybar) ** 2) r2 1 - ss_res / ss_tot return slope, r2 # 示例单轴拉伸沿 z 方向取 σ_zz 对 ε_zz 的斜率 strains [-0.02, -0.01, 0.0, 0.01, 0.02] stresses [-12.5, -6.3, 0.1, 6.2, 12.4] # 单位 GPa C33, r2 fit_elastic_constant(strains, stresses) print(fC33 {C33:.1f} GPa, R^2 {r2:.4f})R² 低于 0.99 就要检查应变范围是否太大、某个应变点的计算是否没收敛、应力提取是否取错了行。金刚石这类高对称材料线性应该非常好R² 接近 1。提示拟合时截距应该接近零。如果截距明显偏离零说明零应变点的结构没有完全弛豫或者应力提取有系统误差。5. 避坑与排查应力应变计算里最容易翻车的五件事5.1 现象应力张量全是零原因VASP 的 INCAR 里没设ISIF 2或IBRION -1或者 QE 的tprnstr没打开。VASP 默认 ISIF2 会算应力但如果 NSW0 且 IBRION-1应力应该正常输出。QE 的tprnstr默认是.false.必须显式设成.true.。解决检查 INCAR 和 QE 输入文件确认应力输出开关打开。VASP 还可以看 OUTCAR 里有没有FORCES acting on ions和Stress tensor段落。5.2 现象不同应变点的能量不连续原因K 点网格或截断能在不同应变点之间变了。有人为了省时间在小应变时用低截断大应变时用高截断导致能量基准不一致。解决所有应变点用完全相同的计算参数。写脚本时把 KPOINTS 和 INCAR 模板固定只改 POSCAR 或 CELL_PARAMETERS。5.3 现象拟合出的弹性常数偏小原因应变范围太大超出了线性弹性区。金刚石的弹性线性范围大概在 ±1% 以内加到 ±3% 曲线就弯了。解决把应变范围缩到 ±1%或者用二次多项式拟合后取一次项系数。更稳妥的做法是先跑一个大范围看曲线拐点再在拐点内取点。5.4 现象QE 和 VASP 算出的弹性常数差很多原因两者的赝势不同、交换关联泛函可能不同、应力单位换算有误。VASP 的 kB 转 GPa 是乘 0.1QE 的 Ry/bohr³ 转 GPa 要乘 14710.5。解决先统一泛函都用 PBE和赝势类型都用 PAW 或都用超软再核对单位换算。如果还差检查 QE 的ecutwfc是否足够QE 对截断能比 VASP 敏感。5.5 现象脚本跑完没有输出文件原因VASP 或 QE 的可执行文件路径没设对或者 mpirun 的进程数和 K 点并行不匹配。脚本里如果硬编码了vasp_std或pw.x环境变量 PATH 里没有就会静默失败。解决在脚本里加错误检查跑完检查 OUTCAR 或 output 文件是否存在且非空。常见做法是用subprocess.run的checkTrue让异常抛出而不是默默跳过。6. 进阶技巧用 plotcheck 脚本做快速验证与批量扫描plotcheck版本的脚本价值在于边算边看。我一般会先用它跑一个应变点确认应力提取和绘图链路通了再切到withoutplot批量跑。批量跑的时候把应变列表写成循环每个应变生成一个目录跑完统一收集数据。# 批量扫描应变并收集应力适合 withoutplot 模式 import os import subprocess import numpy as np def batch_strain_scan(base_dir, strain_list, direction2): 在每个应变下生成结构、跑 VASP、收集应力 results [] for eps in strain_list: work_dir os.path.join(base_dir, fstrain_{eps:.4f}) os.makedirs(work_dir, exist_okTrue) # 生成 POSCAR省略具体实现参考 4.1 节 # 复制 INCAR、KPOINTS、POTCAR 到 work_dir # 跑 VASP subprocess.run(mpirun -np 4 vasp_std vasp.log, shellTrue, cwdwork_dir, checkTrue) # 提取应力 stress read_stress_outcar(os.path.join(work_dir, OUTCAR)) results.append((eps, stress[2, 2])) # 取 σ_zz return results # 跑完直接拟合 strains, stresses zip(*batch_strain_scan(./scan, [-0.01, -0.005, 0, 0.005, 0.01])) C33, r2 fit_elastic_constant(list(strains), list(stresses)) print(fC33 {C33:.1f} GPa, R^2 {r2:.4f})这个批量脚本的关键是每个应变独立目录避免文件覆盖。checkTrue保证 VASP 报错时脚本停下来而不是继续跑下一个应变点。收集完数据直接拟合R² 不合格就回去查哪个应变点出了问题。验证方法上我习惯拿金刚石或硅做基准。金刚石的 C11 约 1076 GPa、C12 约 125 GPa、C44 约 577 GPa算出来在这个量级附近就说明流程通了。如果差一个数量级多半是单位换算错了如果差 20% 到 30%检查赝势和截断能。注意QE 的应力输出在 XML 文件里更规整pw.x的 text output 有时会截断。如果脚本读 text output 不稳定改用style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;" />