尧图网站建设 
简介面向材料科学领域的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.逻辑说明ISYM0是必须加的因为应变会降低体系对称性如果让 VASP 自动检测对称性它可能把原本独立的应变分量合并掉导致你无法区分某个剪切应力分量。ISIF3表示在弹性计算中让晶胞体积和形状一起优化但注意这是用在IBRION2的有限差分弹性常数计算中不是用在你要做应力-应变扫点的那一步。如果真的要做扫点ISIF应该改为 2 或直接固定晶格只放开原子位置。参数表如下参数值作用常见错误ISYM0关闭对称性不关闭时剪切应变分量被混淆ENMAX1.3×默认提高应力收敛度用默认截断能时应力误差约 0.5-1 kBarSIGMA0.05控制部分占据金属体系用 0.2绝缘体可以更低EDIFF1e-7电子步收敛应力对电子步收敛极敏感默认 1e-5 不够NSW60离子步上限应变后原子弛豫会慢用这套参数跑完后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参数说明smearingcold对金属和过渡金属氧化物的应力计算往往比默认的 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] ], dtypefloat) 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( rCELL_PARAMETERS.*?(?ATOMIC_SPECIES), fCELL_PARAMETERS (angstrom)\n \n.join( f{v[0]:.6f} {v[1]:.6f} {v[2]:.6f} for v in new_cell) \n, inp, flagsre.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_outputTrue, textTrue) with open(pw.strain_run.out, r) as f: out f.read() stress_kbar re.findall(rtotal\sstress\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,, headerstrain_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 矩阵单位 kBarQE 的输出里total stress也是一个 3×3 矩阵但默认单位是 Ry/Bohr³只有加上stress .true.并打开verbosityhigh才会在末尾输出换算后的 kBar 版本。字段映射如下物理量VASP 字段QE 字段单位应变POSCAR 晶格矩阵CELL_PARAMETERSÅ应力张量STRESStotal stresskBar外力FORCESForceseV/Å平衡体积VOLUMEcell volumeų一个实用的做法是写一个解析函数输入文件路径和代码类型输出统一单位下的应力张量。下面这段代码处理两种格式提取应力并把非对角分量排除掉因为对大多数正交晶系你关心的只有σxx, σyy, σzz。import re def read_stress_from_output(output_path, codevasp): with open(output_path, r) as f: text f.read() if code vasp: match re.search(rSTRESS\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(rstress 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(rtotal\sstress\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, codeqe) 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, indexFalse)逻辑说明把应力和应变关系做成一个 DataFrame 后后续的线性拟合就完全不需要再碰文本文件了。stress_yy在单轴应变下也不为零它是对横向响应的探测如果横向应力有异常大幅变化说明晶格方向设置错了——单轴应变应当只影响纵向分量横向分量保持接近零如果用的是固定横向尺寸的模式。4.3 Python 数据分析与可视化中容易踩到的小数精度问题对 0.5% 量级的应变VASP 的应力输出只有三位有效数字并且应力张量的输出精度由IO格式决定。OUTCAR 里的应力通常保留到小数点后 1 位kBar对弹性模量几个 GPa 量级的材料1 kBar 的误差会带来 10% 左右的模量误差。解决办法是不要在文本层面提高精度而是用 VASP 的DFTU或大截断能让应力在电子步上收敛到更小的值然后直接在 OUTCAR 里用程序读原始数值不要用人工复制粘贴。有一种常见错误用 Python 的float()转换1.23456E02这类字符串没问题但如果 QE 输出的负号前有空格比如-1.234E01split()会把它正确切分可如果是 Fortran 的write格式产生-1.234E01后面没有多余空格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(fC11 {linreg.slope:.2f} kBar {linreg.slope * 0.1:.2f} GPa) print(fR^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, fR² {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 比是否在应变后仍保持优化值——这是单元晶胞计算中最常被忽略的系统误差。本文还有配套的精品资源点击获取 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.逻辑说明ISYM0是必须加的因为应变会降低体系对称性如果让 VASP 自动检测对称性它可能把原本独立的应变分量合并掉导致你无法区分某个剪切应力分量。ISIF3表示在弹性计算中让晶胞体积和形状一起优化但注意这是用在IBRION2的有限差分弹性常数计算中不是用在你要做应力-应变扫点的那一步。如果真的要做扫点ISIF应该改为 2 或直接固定晶格只放开原子位置。参数表如下参数值作用常见错误ISYM0关闭对称性不关闭时剪切应变分量被混淆ENMAX1.3×默认提高应力收敛度用默认截断能时应力误差约 0.5-1 kBarSIGMA0.05控制部分占据金属体系用 0.2绝缘体可以更低EDIFF1e-7电子步收敛应力对电子步收敛极敏感默认 1e-5 不够NSW60离子步上限应变后原子弛豫会慢用这套参数跑完后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参数说明smearingcold对金属和过渡金属氧化物的应力计算往往比默认的 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] ], dtypefloat) 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( rCELL_PARAMETERS.*?(?ATOMIC_SPECIES), fCELL_PARAMETERS (angstrom)\n \n.join( f{v[0]:.6f} {v[1]:.6f} {v[2]:.6f} for v in new_cell) \n, inp, flagsre.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_outputTrue, textTrue) with open(pw.strain_run.out, r) as f: out f.read() stress_kbar re.findall(rtotal\sstress\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,, headerstrain_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 矩阵单位 kBarQE 的输出里total stress也是一个 3×3 矩阵但默认单位是 Ry/Bohr³只有加上stress .true.并打开verbosityhigh才会在末尾输出换算后的 kBar 版本。字段映射如下物理量VASP 字段QE 字段单位应变POSCAR 晶格矩阵CELL_PARAMETERSÅ应力张量STRESStotal stresskBar外力FORCESForceseV/Å平衡体积VOLUMEcell volumeų一个实用的做法是写一个解析函数输入文件路径和代码类型输出统一单位下的应力张量。下面这段代码处理两种格式提取应力并把非对角分量排除掉因为对大多数正交晶系你关心的只有σxx, σyy, σzz。import re def read_stress_from_output(output_path, codevasp): with open(output_path, r) as f: text f.read() if code vasp: match re.search(rSTRESS\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(rstress 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(rtotal\sstress\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, codeqe) 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, indexFalse)逻辑说明把应力和应变关系做成一个 DataFrame 后后续的线性拟合就完全不需要再碰文本文件了。stress_yy在单轴应变下也不为零它是对横向响应的探测如果横向应力有异常大幅变化说明晶格方向设置错了——单轴应变应当只影响纵向分量横向分量保持接近零如果用的是固定横向尺寸的模式。4.3 Python 数据分析与可视化中容易踩到的小数精度问题对 0.5% 量级的应变VASP 的应力输出只有三位有效数字并且应力张量的输出精度由IO格式决定。OUTCAR 里的应力通常保留到小数点后 1 位kBar对弹性模量几个 GPa 量级的材料1 kBar 的误差会带来 10% 左右的模量误差。解决办法是不要在文本层面提高精度而是用 VASP 的DFTU或大截断能让应力在电子步上收敛到更小的值然后直接在 OUTCAR 里用程序读原始数值不要用人工复制粘贴。有一种常见错误用 Python 的float()转换1.23456E02这类字符串没问题但如果 QE 输出的负号前有空格比如-1.234E01split()会把它正确切分可如果是 Fortran 的write格式产生-1.234E01后面没有多余空格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(fC11 {linreg.slope:.2f} kBar {linreg.slope * 0.1:.2f} GPa) print(fR^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, fR² {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 比是否在应变后仍保持优化值——这是单元晶胞计算中最常被忽略的系统误差。本文还有配套的精品资源点击获取