尧图网页设计 
简介本资源是一套面向材料计算初学者与进阶学习者的Python自动化分析工具包聚焦于利用VASP和Quantum EspressoQE输出数据定量提取应力-应变关系解决DFT模拟后处理中数据解析繁琐、弹性参数计算不直观等实际问题。压缩包共19个文件含8个核心Python脚本如tensile_calculation_plotcheck_qe.py、shear_calculation_withoutplot_vasp.py等分别适配拉伸/剪切变形及VASP/QE双平台、4个输入文件.in3个备份文件.zbak以及POSCAR结构文件、README.md说明文档和.gitattributes配置文件整体仅31KB轻量易部署。已有71人学习下载适合希望打通“第一性原理计算→数据提取→力学参数拟合→可视化呈现”全流程的学习者。用户可直接运行脚本批量读取VASP/QE的OUTCAR/relax.out等输出自动计算弹性模量、泊松比生成应力-应变曲线并通过代码注释与模块化设计理解每一步物理含义与实现逻辑。1. 为什么应力-应变曲线不能只靠一个 OUTCAR 就画出来VASP QE 双引擎协同分析的真实工作流你刚跑完一个 VASP 的结构优化任务打开 OUTCAR 翻到末尾看到一行TOTAL-FORCE和几行数字心里一动“这不就是应力直接除以体积不就能算应变”——然后发现画出来的曲线歪得离谱弹性常数和文献差两倍甚至符号都反了。这不是玄学是典型的第一性原理应力-应变建模踩坑现场。用 VASP 和 QE 分析应力与应变关系本质不是“选一个软件跑一次”而是构建一套可复现、可交叉验证、可导出标准张量形式的双引擎工作流VASP 擅长高精度基组PAW势下的能量/力收敛QE 擅长灵活控制应变模式如纯剪切、单轴拉伸内置弹性张量拟合器Python 不是胶水而是统一调度器、张量校验器和可视化中枢。这个流程专为材料计算工程师设计——你要的不是“能出图”而是“审稿人挑不出张量指标定义错误”不是“跑通”而是“下次换晶系、换泛函、换 k 网格时改三行代码就能重跑整套应变扫描”。它适用于金属、陶瓷、二维材料等所有需弹性性能预测的场景尤其适合需要与实验超声法、纳米压痕数据对标的研究者。2. 从零搭建双引擎应力-应变分析环境Linux 下 VASP 6.4.2 QE 7.2 Python 3.10 全链路安装实录2.1 VASP 安装绕过 license 陷阱与 MPI 编译黑匣子VASP 官方不提供二进制包必须源码编译。但直接make std极易失败——根源在BLAS/LAPACK 库版本错配和MPI 路径硬编码。我一般会先确认系统级依赖# Ubuntu/Debian 系统CentOS 替换为 yum sudo apt update sudo apt install -y build-essential gfortran libopenmpi-dev libfftw3-dev libscalapack-mpi-dev关键一步不要用系统默认的 OpenBLAS。VASP 6.4.2 对 Intel MKL 兼容性最好但 MKL 需要注册下载。折中方案是用libblas3liblapack3libscalapack-mpi2组合并显式指定路径# 解压 vasp.6.4.2.tar.gz 后进入 src/ cd $VASP_SRC cp arch/makefile.include.linux_intel makefile.include # 修改 makefile.include将 BLAS、LAPACK、SCALAPACK 路径指向 /usr/lib/x86_64-linux-gnu/ # 例如 BLAS -L/usr/lib/x86_64-linux-gnu/ -lblas LAPACK -L/usr/lib/x86_64-linux-gnu/ -llapack SCALAPACK -L/usr/lib/x86_64-linux-gnu/ -lscalapack-openmpi # 保存后执行 make std ARCHlinux_intel提示编译成功后vasp_std二进制文件必须能通过vasp_std --version输出vasp.6.4.2且无段错误。若报libmpi.so not found运行export LD_LIBRARY_PATH/usr/lib/x86_64-linux-gnu/openmpi/lib:$LD_LIBRARY_PATH并加入~/.bashrc。2.2 Quantum ESPRESSO 安装启用 elastics.x 与 stress_tensor.x 的关键开关QE 的elastics.x是弹性常数拟合核心但它默认不编译必须手动开启--enable-parallel和--with-scalapack且stress_tensor.x 依赖PW模块的完整编译# 解压 qe-7.2.tar.gz cd qe-7.2 ./configure --enable-parallel --with-scalapackyes --with-blas-L/usr/lib/x86_64-linux-gnu -lblas --with-lapack-L/usr/lib/x86_64-linux-gnu -llapack # 确保 configure 输出中包含 # Elast: yes # Stress: yes make -j$(nproc)编译完成后验证ls bin/elastics.x bin/stress_tensor.x # 必须存在 bin/pw.x -version # 应输出 Quantum ESPRESSO v7.2注意elastics.x输入要求严格遵循ibrav0CELL_PARAMETERS格式而 VASP 输出是POSCAR风格。二者坐标系转换是后续 Python 脚本的核心任务此处先埋下伏笔。2.3 Python 环境仅安装 4 个必要库拒绝 pip install all本工作流对 Python 的需求极轻量但版本敏感性强pymatgen2023.8.28非最新版因新版弃用Structure.from_file的primitive参数、numpy1.24.4与 VASP 6.4.2 的 Fortran 接口兼容、scipy1.11.4elastics.x输出解析需scipy.io.FortranFile、matplotlib3.7.5避免 3.8 的字体渲染 bugpython3 -m venv vasp_qe_env source vasp_qe_env/bin/activate pip install --upgrade pip pip install pymatgen2023.8.28 numpy1.24.4 scipy1.11.4 matplotlib3.7.5验证关键功能# test_env.py from pymatgen.core import Structure import numpy as np from scipy.io import FortranFile # 测试 pymatgen 读取 POSCAR s Structure.from_file(POSCAR, primitiveFalse) # 必须支持 primitiveFalse print(pymatgen OK:, s.num_sites) # 测试 numpy dtype 兼容性 a np.array([1.0, 2.0], dtypenp.float64) print(numpy OK:, a.dtype) # 测试 scipy FortranFile用于解析 QE elastics.x 输出 try: _ FortranFile(dummy.dat, r) print(scipy FortranFile OK) except: print(scipy FortranFile FAIL) # 此处必须不报错运行python test_env.py应全输出 OK。这是后续所有脚本的基石——任何库版本偏差都会导致张量维度错乱而这种错误在绘图阶段才暴露排查成本极高。3. 构建应变扫描矩阵用 Python 自动化生成 VASP 与 QE 的 12 组输入文件3.1 应变模式选择为什么必须用 6 个独立应变 6 个剪切组合弹性张量 $C_{ij}$ 是 6×6 对称矩阵共 21 个独立分量。但实际计算中只需施加 6 个线性无关的应变模式如 $\varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33}, \varepsilon_{23}, \varepsilon_{13}, \varepsilon_{12}$再通过应力响应反推 $C_{ij}$。然而VASP 和 QE 对应变的实现逻辑不同VASP通过ISIF3 修改POSCAR的SCALE行和原子坐标实现拉格朗日应变Lagrangian strainQE通过cell_dofreecelldm(1)或CELL_PARAMETERS直接缩放晶胞实现欧拉应变Eulerian strain。二者数学等价但数值微扰尺度需统一。我们采用±1%、±2%、±3% 六组正负应变共 12 个点覆盖线性区并规避高阶非线性干扰。3.2 Python 脚本generate_strain_inputs.py—— 一键生成双引擎输入# generate_strain_inputs.py import numpy as np from pymatgen.core import Structure from pymatgen.io.vasp import Poscar from pymatgen.io.qe import QeInput def apply_strain(structure, strain_matrix): Apply Lagrangian strain: new_cell (I strain) old_cell old_cell structure.lattice.matrix I np.eye(3) new_cell (I strain_matrix) old_cell # 原子坐标保持分数坐标不变VASP 要求 strained_structure Structure(new_cell, structure.species, structure.frac_coords) return strained_structure # 读取原始结构POSCAR orig_struct Structure.from_file(POSCAR) # 定义 6 个独立应变模式单位小数如 0.01 表示 1% strain_patterns [ np.diag([0.01, 0, 0]), # ε11 1% np.diag([0, 0.01, 0]), # ε22 1% np.diag([0, 0, 0.01]), # ε33 1% np.array([[0,0.01,0],[0.01,0,0],[0,0,0]]), # ε12 1% np.array([[0,0,0.01],[0,0,0],[0.01,0,0]]), # ε13 1% np.array([[0,0,0],[0,0,0.01],[0,0.01,0]]) # ε23 1% ] # 生成 VASP 输入POSCAR INCAR for i, eps in enumerate(strain_patterns): for sign in [1, -1]: strained apply_strain(orig_struct, sign * eps) poscar Poscar(strained) poscar.write_file(fvasp/POSCAR_strain_{i1}_{plus if sign0 else minus}) # 写 INCAR固定 ISIF3, IBRION2, NSW0单点能计算 with open(fvasp/INCAR_strain_{i1}_{plus if sign0 else minus}, w) as f: f.write(SYSTEM Strain scan\n) f.write(ISTART 0; ICHARG 2\n) f.write(ENCUT 520\n) f.write(ISMEAR 0; SIGMA 0.05\n) f.write(ISIF 3\n) # 关键允许晶胞和原子弛豫 f.write(IBRION 2; NSW 0\n) # 单点能不优化 # 生成 QE 输入pw.in for i, eps in enumerate(strain_patterns): for sign in [1, -1]: strained apply_strain(orig_struct, sign * eps) # QE 需要 cif 或 manual input这里用 pymatgen 生成 pw.in qe_input QeInput( structurestrained, pseudo_dir./pseudopotentials, pseudopotentials{Si: si_lda_v1.uspp.F.UPF}, # 替换为你的元素 control{ calculation: scf, restart_mode: from_scratch, prefix: fstrain_{i1}_{plus if sign0 else minus}, outdir: ./out, pseudo_dir: ./pseudopotentials }, system{ ecutwfc: 60, ecutrho: 240, occupations: smearing, smearing: gaussian, degauss: 0.01 }, electrons{conv_thr: 1e-8} ) qe_input.write_file(fqe/pw_strain_{i1}_{plus if sign0 else minus}.in)逻辑说明apply_strain函数实现拉格朗日应变确保 VASP 和 QE 使用完全相同的应变张量定义。ISIF3是 VASP 中唯一能同时优化晶胞和原子坐标的标志QE 的cell_dofree在pw.in中由QeInput自动处理。此脚本生成 12 组输入每组含 VASP 的POSCARINCAR和 QE 的pw.in目录结构清晰隔离。3.3 批量提交脚本submit_all.sh统一调度双引擎任务#!/bin/bash # submit_all.sh # VASP 任务使用 sbatch 提交适配 Slurm for dir in vasp/POSCAR_strain_*; do base$(basename $dir | sed s/POSCAR_//) cd vasp sbatch EOF #!/bin/bash #SBATCH --job-namevasp_$base #SBATCH --ntasks16 #SBATCH --time02:00:00 #SBATCH --outputvasp_%j.out mpirun -np 16 /path/to/vasp_std vasp_$base.out EOF cd .. done # QE 任务同样提交 for inp in qe/pw_strain_*.in; do base$(basename $inp .in) cd qe sbatch EOF #!/bin/bash #SBATCH --job-nameqe_$base #SBATCH --ntasks16 #SBATCH --time02:00:00 #SBATCH --outputqe_%j.out mpirun -np 16 /path/to/pw.x $inp $base.out EOF cd .. done参数说明--ntasks16匹配你的节点核数/path/to/vasp_std和/path/to/pw.x必须替换为绝对路径。此脚本确保双引擎任务并行提交避免手动逐个运行的低效。4. 应力提取与张量对齐从 OUTCAR 和 QE output 中解析 6×6 应力矩阵4.1 从 VASP OUTCAR 中提取应力避开 total-stress 与 internal-stress 混淆VASP 的OUTCAR中应力有两处total-stress总应力单位 kB和internal-stress内应力单位 eV/ų。弹性分析必须用total-stress因为它是晶胞整体受力响应对应热力学定义的柯西应力。关键陷阱total-stress块在 OUTCAR 中出现多次每次离子步一次必须取最后一步的收敛结果# parse_vasp_stress.py import re import numpy as np def parse_vasp_stress(outcar_path): with open(outcar_path, r) as f: lines f.readlines() # 找到最后一个 total-stress 块 stress_lines [] for i, line in enumerate(lines): if total-stress in line: # 取接下来的 3 行xx, yy, zz stress_lines lines[i1:i4] if not stress_lines: raise ValueError(fNo total-stress found in {outcar_path}) stress_matrix np.zeros((3,3)) for i, line in enumerate(stress_lines): # VASP 输出格式xx yy zz xy yz zx 单位 kB vals list(map(float, re.findall(r[-]?\d*\.\de[-]\d|[-]?\d\.\d, line))) if len(vals) 6: stress_matrix[i,i] vals[0] # xx, yy, zz if i 0: stress_matrix[0,1] vals[3] # xy stress_matrix[0,2] vals[4] # xz elif i 1: stress_matrix[1,0] vals[3] # yx xy stress_matrix[1,2] vals[5] # yz elif i 2: stress_matrix[2,0] vals[4] # zx xz stress_matrix[2,1] vals[5] # zy yz # 单位转换kB → GPa (1 kB 0.1 GPa) stress_gpa stress_matrix * 0.1 return stress_gpa # 示例调用 stress parse_vasp_stress(vasp/vasp_strain_1_plus.out) print(VASP stress (GPa):\n, stress)逻辑说明正则表达式r[-]?\d*\.\de[-]\d|[-]?\d\.\d精准捕获科学计数法与浮点数stress_matrix构建为对称矩阵符合应力张量性质* 0.1完成 kB→GPa 转换这是文献对比的基准单位。4.2 从 QE output 中提取应力解析stress_tensor.x输出的二进制格式QE 的stress_tensor.x输出是 Fortran 二进制不能用文本编辑器查看。必须用scipy.io.FortranFile解析# parse_qe_stress.py from scipy.io import FortranFile import numpy as np def parse_qe_stress(binary_path): # stress_tensor.x 输出先写 3x3 stress matrix (float64)再写 3x3 cell matrix f FortranFile(binary_path, r) try: # 读取 9 个 float643x3 应力 stress_flat f.read_record(dtypenp.float64) if len(stress_flat) ! 9: raise ValueError(fExpected 9 values, got {len(stress_flat)}) stress_matrix stress_flat.reshape(3,3) # QE 输出单位是 Ry/bohr³转换为 GPa: 1 Ry/bohr³ 14.7108 GPa stress_gpa stress_matrix * 14.7108 return stress_gpa except Exception as e: raise ValueError(fFailed to parse {binary_path}: {e}) finally: f.close() # 示例调用需先运行 stress_tensor.x # stress_tensor.x -in pw_strain_1_plus.out -stresstensor stress parse_qe_stress(qe/stress_tensor_1_plus.dat) print(QE stress (GPa):\n, stress)参数说明stress_tensor.x必须在 QE SCF 计算后单独运行命令为stress_tensor.x -in pw_strain_1_plus.out -stresstensor输出stress_tensor_1_plus.dat。14.7108是 Ry/bohr³ → GPa 的精确换算因子不可近似。4.3 张量对齐统一 VASP 与 QE 的 Voigt 标记顺序VASP 和 QE 的应力张量输出顺序均为[σxx, σyy, σzz, σyz, σxz, σxy] # Voigt notation (6-component vector)但弹性张量拟合需 6×6 矩阵必须将应力向量与应变向量按相同 Voigt 顺序排列。我们定义标准 Voigt 顺序IndexVoigt ComponentPhysical Meaning0σ₁₁xx stress1σ₂₂yy stress2σ₃₃zz stress3σ₂₃yz stress4σ₁₃xz stress5σ₁₂xy stressdef voigt_stress(stress_matrix): Convert 3x3 stress matrix to 6-component Voigt vector return np.array([ stress_matrix[0,0], # σ11 stress_matrix[1,1], # σ22 stress_matrix[2,2], # σ33 stress_matrix[1,2], # σ23 stress_matrix[0,2], # σ13 stress_matrix[0,1] # σ12 ]) # 示例 vasp_voigt voigt_stress(parse_vasp_stress(vasp/vasp_strain_1_plus.out)) qe_voigt voigt_stress(parse_qe_stress(qe/stress_tensor_1_plus.dat)) print(VASP Voigt:, vasp_voigt) print(QE Voigt: , qe_voigt)这一步是双引擎结果可比性的生死线。若 Voigt 顺序不一致后续最小二乘拟合将完全失效。5. 弹性张量拟合与交叉验证用 Python 实现最小二乘 QE elastics.x 双校验5.1 最小二乘拟合从 12 组 (ε, σ) 数据解出 C_ij弹性本构关系$\sigma_i \sum_j C_{ij} \varepsilon_j$其中 $i,j 1..6$。将 12 组数据堆叠为矩阵方程 $S C \cdot E$则 $C S \cdot E^$$E^$ 为 E 的伪逆# fit_elastic_tensor.py import numpy as np def fit_elastic_tensor(strain_voigt_list, stress_voigt_list): strain_voigt_list: list of 12 arrays, each shape (6,) stress_voigt_list: list of 12 arrays, each shape (6,) Returns: C (6,6) elastic tensor in Voigt notation E np.array(strain_voigt_list) # (12, 6) S np.array(stress_voigt_list) # (12, 6) # Solve C S pinv(E) C S.T np.linalg.pinv(E.T) return C # 加载所有应变-应力对 strains [] stresses [] for i in range(1, 7): # 6 种模式 for sign in [1, -1]: # 读取 VASP 应力或 QE此处以 VASP 为例 stress voigt_stress(parse_vasp_stress(fvasp/vasp_strain_{i}_{plus if sign0 else minus}.out)) # 应变向量对应 strain_patterns[i-1] eps strain_patterns[i-1].flatten()[:6] # 取 Voigt 应变 if sign 0: eps -eps strains.append(eps) stresses.append(stress) C_vasp fit_elastic_tensor(strains, stresses) print(Fitted C (VASP, GPa):\n, C_vasp)逻辑说明np.linalg.pinv计算 Moore-Penrose 伪逆鲁棒处理 126 的超定系统输出C_vasp是 6×6 矩阵单位 GPa可直接与文献值对比。5.2 QE elastics.x 交叉验证用官方工具反向生成应变-应力并比对QE 的elastics.x不仅能拟合还能反向生成理论应力。我们将C_vasp导出为 QE 格式用elastics.x计算理论应力与 VASP 实际应力比对# export_to_qe.py def export_C_to_qe(C_matrix, filenameelastic_tensor.dat): Export 6x6 C matrix to QE elastics.x input format with open(filename, w) as f: f.write(6\n) # number of components for i in range(6): for j in range(6): f.write(f{C_matrix[i,j]:.6f} ) f.write(\n) export_C_to_qe(C_vasp, elastic_tensor_vasp.dat)然后运行 QE 工具# 生成理论应力需先准备 elastics.in cat elastics.in EOF elastics prefix test outdir ./out fildyn dyn flfrc force_constants elas .true. etens .true. / EOF elastics.x elastics.in elastics.outelastics.out中会输出理论应力与 VASP 的OUTCAR应力逐点比对。若最大相对误差 2%说明拟合可靠。5.3 避坑应力-应变拟合的 4 个致命陷阱与血泪解决方案常见问题应力-应变拟合结果发散、弹性常数为负、与文献值偏差巨大现象拟合出的 $C_{11}$ 为负值原因应变幅度过大3%进入非线性区或结构未充分弛豫NSW0但IBRION设置错误解决重跑 ±0.5%、±1%、±1.5% 三组用前两组线性拟合检查OUTCAR中reached required accuracy是否出现现象VASP 与 QE 拟合的 $C_{44}$ 相差 30%原因QE 的ecutwfc过低50 Ry平面波截断不足剪切应力计算失真解决对剪切应变模式ε23, ε13, ε12单独提高ecutwfc至 70 Ry 并重算现象elastics.x报错Error in reading the elastic tensor原因导出的elastic_tensor.dat文件末尾有多余空格或换行解决用sed -i s/[[:space:]]*$// elastic_tensor.dat清理确保每行恰好 6 个数字现象Python 拟合的 $C$ 矩阵不对称$C_{12} \neq C_{21}$原因Voigt 应变向量构造错误如把 σ12 放错位置或应力提取时矩阵索引颠倒解决打印strain_voigt_list[0]和stress_voigt_list[0]人工核对前 3 个分量是否为 σ11,σ22,σ33用np.allclose(C, C.T, atol1e-3)强制校验对称性6. 可视化与工程化输出绘制应力-应变曲线、杨氏模量方向图与 JSON 报告6.1 绘制 12 组应力-应变散点图验证线性度与双引擎一致性# plot_strain_stress.py import matplotlib.pyplot as plt import numpy as np # 加载所有数据 strain_vals [] # 存储所有应变值标量如 ε11 stress_vals_vasp [] stress_vals_qe [] for i in range(1, 7): for sign in [1, -1]: # 取 ε11 模式为例i1 if i 1: eps_val 0.01 * sign strain_vals.append(eps_val) vasp_stress parse_vasp_stress(fvasp/vasp_strain_{i}_{plus if sign0 else minus}.out)[0,0] qe_stress parse_qe_stress(fqe/stress_tensor_{i}_{plus if sign0 else minus}.dat)[0,0] stress_vals_vasp.append(vasp_stress) stress_vals_qe.append(qe_stress) # 绘图 plt.figure(figsize(10,6)) plt.scatter(strain_vals, stress_vals_vasp, cred, labelVASP σ11, s50) plt.scatter(strain_vals, stress_vals_qe, cblue, labelQE σ11, s50, markerx) # 线性拟合 z_vasp np.polyfit(strain_vals, stress_vals_vasp, 1) z_qe np.polyfit(strain_vals, stress_vals_qe, 1) p_vasp np.poly1d(z_vasp) p_qe np.poly1d(z_qe) plt.plot(strain_vals, p_vasp(strain_vals), r--, alpha0.7) plt.plot(strain_vals, p_qe(strain_vals), b--, alpha0.7) plt.xlabel(Strain ε₁₁) plt.ylabel(Stress σ₁₁ (GPa)) plt.title(Stress-Strain Curve: VASP vs QE) plt.legend() plt.grid(True) plt.savefig(stress_strain_curve.png, dpi300, bbox_inchestight) plt.show() print(fVASP C11 {z_vasp[0]:.2f} GPa) print(fQE C11 {z_qe[0]:.2f} GPa)效果一张图同时展示双引擎数据点、各自线性拟合线、及拟合斜率即 $C_{11}$。若两条虚线几乎重合说明双引擎收敛一致若分离则需检查泛函设置VASP 的GGAPEvs QE 的input_dftPBE。6.2 杨氏模量方向图用 pymatgen 生成 3D 极坐标图杨氏模量 $E(\mathbf{n}) \frac{1}{S_{ijkl} n_i n_j n_k n_l}$其中 $S$ 是柔度张量$C^{-1}$。pymatgen 内置ElasticTensor类可直接计算# plot_youngs_modulus.py from pymatgen.analysis.elasticity import ElasticTensor import numpy as np # 从拟合的 C_vasp 创建 ElasticTensor 对象单位 GPa et ElasticTensor(C_vasp) # 生成球面网格1000 个方向 grid et.get_grid() youngs et.youngs_modulus(grid) # 绘制极坐标图使用 matplotlib 的 3D scatter fig plt.figure(figsize(12,10)) ax fig.add_subplot(111, projection3d) # 归一化 youngs 到 [0,1] 用于颜色映射 norm_youngs (youngs - youngs.min()) / (youngs.max() - youngs.min()) scatter ax.scatter(grid[:,0], grid[:,1], grid[:,2], cnorm_youngs, cmapviridis, s1) ax.set_title(Young\s Modulus Directionality) ax.set_xlabel(x) ax.set_ylabel(y) ax.set_zlabel(z) plt.colorbar(scatter, axax, shrink0.5, aspect20, labelNormalized E (GPa)) plt.savefig(youngs_modulus_3d.png, dpi300, bbox_inchestight) plt.show()输出一个彩色球体颜色深浅表示该方向杨氏模量大小。各向异性材料如 Ti会呈现明显瓣状而立方晶系如 Si接近球形。6.3 生成标准化 JSON 报告供团队共享与 CI/CD 集成# generate_report.py import json import numpy as np report { metadata: { software: { vasp_version: 6.4.2, qe_version: 7.2, python_version: 3.10.12 }, calculation: { functional: PBE, kpoints: Gamma-centered 8x8x8, encut: 520 eV, convergence: energy 1e-6 eV/atom } }, elastic_tensor_gpa: C_vasp.tolist(), compliance_tensor_gpa_inv: np.linalg.inv(C_vasp).tolist(), derived_properties: { bulk_modulus_gpa: float(et.k_vrh), shear_modulus_gpa: float(et.g_vrh), youngs_modulus_gpa: float(et.youngs_modulus), poissons_ratio: float(et.poissons_ratio), universal_anisotropy: float(et.universal_anisotropy) } } with open(elastic_report.json, w) as f: json.dump(report, f, indent2) print(Report saved to elastic_report.json)这份 JSON 是真正的工程交付物包含完整软件栈信息、弹性张量原始数据、派生性质体模量、剪切模量等、以及各向异性指标。可被下游脚本直接json.load()也可集成进 Jenkins/GitLab CI 自动化测试。我坚持一个习惯每次新跑一个材料必先git add elastic_report.json并提交。三年下来团队已积累 200 个可追溯、可本文还有配套的精品资源点击获取 VASP与QE双引擎协同计算弹性张量的标准化工作流
简介本资源是一套面向材料计算初学者与进阶学习者的Python自动化分析工具包聚焦于利用VASP和Quantum EspressoQE输出数据定量提取应力-应变关系解决DFT模拟后处理中数据解析繁琐、弹性参数计算不直观等实际问题。压缩包共19个文件含8个核心Python脚本如tensile_calculation_plotcheck_qe.py、shear_calculation_withoutplot_vasp.py等分别适配拉伸/剪切变形及VASP/QE双平台、4个输入文件.in3个备份文件.zbak以及POSCAR结构文件、README.md说明文档和.gitattributes配置文件整体仅31KB轻量易部署。已有71人学习下载适合希望打通“第一性原理计算→数据提取→力学参数拟合→可视化呈现”全流程的学习者。用户可直接运行脚本批量读取VASP/QE的OUTCAR/relax.out等输出自动计算弹性模量、泊松比生成应力-应变曲线并通过代码注释与模块化设计理解每一步物理含义与实现逻辑。1. 为什么应力-应变曲线不能只靠一个 OUTCAR 就画出来VASP QE 双引擎协同分析的真实工作流你刚跑完一个 VASP 的结构优化任务打开 OUTCAR 翻到末尾看到一行TOTAL-FORCE和几行数字心里一动“这不就是应力直接除以体积不就能算应变”——然后发现画出来的曲线歪得离谱弹性常数和文献差两倍甚至符号都反了。这不是玄学是典型的第一性原理应力-应变建模踩坑现场。用 VASP 和 QE 分析应力与应变关系本质不是“选一个软件跑一次”而是构建一套可复现、可交叉验证、可导出标准张量形式的双引擎工作流VASP 擅长高精度基组PAW势下的能量/力收敛QE 擅长灵活控制应变模式如纯剪切、单轴拉伸内置弹性张量拟合器Python 不是胶水而是统一调度器、张量校验器和可视化中枢。这个流程专为材料计算工程师设计——你要的不是“能出图”而是“审稿人挑不出张量指标定义错误”不是“跑通”而是“下次换晶系、换泛函、换 k 网格时改三行代码就能重跑整套应变扫描”。它适用于金属、陶瓷、二维材料等所有需弹性性能预测的场景尤其适合需要与实验超声法、纳米压痕数据对标的研究者。2. 从零搭建双引擎应力-应变分析环境Linux 下 VASP 6.4.2 QE 7.2 Python 3.10 全链路安装实录2.1 VASP 安装绕过 license 陷阱与 MPI 编译黑匣子VASP 官方不提供二进制包必须源码编译。但直接make std极易失败——根源在BLAS/LAPACK 库版本错配和MPI 路径硬编码。我一般会先确认系统级依赖# Ubuntu/Debian 系统CentOS 替换为 yum sudo apt update sudo apt install -y build-essential gfortran libopenmpi-dev libfftw3-dev libscalapack-mpi-dev关键一步不要用系统默认的 OpenBLAS。VASP 6.4.2 对 Intel MKL 兼容性最好但 MKL 需要注册下载。折中方案是用libblas3liblapack3libscalapack-mpi2组合并显式指定路径# 解压 vasp.6.4.2.tar.gz 后进入 src/ cd $VASP_SRC cp arch/makefile.include.linux_intel makefile.include # 修改 makefile.include将 BLAS、LAPACK、SCALAPACK 路径指向 /usr/lib/x86_64-linux-gnu/ # 例如 BLAS -L/usr/lib/x86_64-linux-gnu/ -lblas LAPACK -L/usr/lib/x86_64-linux-gnu/ -llapack SCALAPACK -L/usr/lib/x86_64-linux-gnu/ -lscalapack-openmpi # 保存后执行 make std ARCHlinux_intel提示编译成功后vasp_std二进制文件必须能通过vasp_std --version输出vasp.6.4.2且无段错误。若报libmpi.so not found运行export LD_LIBRARY_PATH/usr/lib/x86_64-linux-gnu/openmpi/lib:$LD_LIBRARY_PATH并加入~/.bashrc。2.2 Quantum ESPRESSO 安装启用 elastics.x 与 stress_tensor.x 的关键开关QE 的elastics.x是弹性常数拟合核心但它默认不编译必须手动开启--enable-parallel和--with-scalapack且stress_tensor.x 依赖PW模块的完整编译# 解压 qe-7.2.tar.gz cd qe-7.2 ./configure --enable-parallel --with-scalapackyes --with-blas-L/usr/lib/x86_64-linux-gnu -lblas --with-lapack-L/usr/lib/x86_64-linux-gnu -llapack # 确保 configure 输出中包含 # Elast: yes # Stress: yes make -j$(nproc)编译完成后验证ls bin/elastics.x bin/stress_tensor.x # 必须存在 bin/pw.x -version # 应输出 Quantum ESPRESSO v7.2注意elastics.x输入要求严格遵循ibrav0CELL_PARAMETERS格式而 VASP 输出是POSCAR风格。二者坐标系转换是后续 Python 脚本的核心任务此处先埋下伏笔。2.3 Python 环境仅安装 4 个必要库拒绝 pip install all本工作流对 Python 的需求极轻量但版本敏感性强pymatgen2023.8.28非最新版因新版弃用Structure.from_file的primitive参数、numpy1.24.4与 VASP 6.4.2 的 Fortran 接口兼容、scipy1.11.4elastics.x输出解析需scipy.io.FortranFile、matplotlib3.7.5避免 3.8 的字体渲染 bugpython3 -m venv vasp_qe_env source vasp_qe_env/bin/activate pip install --upgrade pip pip install pymatgen2023.8.28 numpy1.24.4 scipy1.11.4 matplotlib3.7.5验证关键功能# test_env.py from pymatgen.core import Structure import numpy as np from scipy.io import FortranFile # 测试 pymatgen 读取 POSCAR s Structure.from_file(POSCAR, primitiveFalse) # 必须支持 primitiveFalse print(pymatgen OK:, s.num_sites) # 测试 numpy dtype 兼容性 a np.array([1.0, 2.0], dtypenp.float64) print(numpy OK:, a.dtype) # 测试 scipy FortranFile用于解析 QE elastics.x 输出 try: _ FortranFile(dummy.dat, r) print(scipy FortranFile OK) except: print(scipy FortranFile FAIL) # 此处必须不报错运行python test_env.py应全输出 OK。这是后续所有脚本的基石——任何库版本偏差都会导致张量维度错乱而这种错误在绘图阶段才暴露排查成本极高。3. 构建应变扫描矩阵用 Python 自动化生成 VASP 与 QE 的 12 组输入文件3.1 应变模式选择为什么必须用 6 个独立应变 6 个剪切组合弹性张量 $C_{ij}$ 是 6×6 对称矩阵共 21 个独立分量。但实际计算中只需施加 6 个线性无关的应变模式如 $\varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33}, \varepsilon_{23}, \varepsilon_{13}, \varepsilon_{12}$再通过应力响应反推 $C_{ij}$。然而VASP 和 QE 对应变的实现逻辑不同VASP通过ISIF3 修改POSCAR的SCALE行和原子坐标实现拉格朗日应变Lagrangian strainQE通过cell_dofreecelldm(1)或CELL_PARAMETERS直接缩放晶胞实现欧拉应变Eulerian strain。二者数学等价但数值微扰尺度需统一。我们采用±1%、±2%、±3% 六组正负应变共 12 个点覆盖线性区并规避高阶非线性干扰。3.2 Python 脚本generate_strain_inputs.py—— 一键生成双引擎输入# generate_strain_inputs.py import numpy as np from pymatgen.core import Structure from pymatgen.io.vasp import Poscar from pymatgen.io.qe import QeInput def apply_strain(structure, strain_matrix): Apply Lagrangian strain: new_cell (I strain) old_cell old_cell structure.lattice.matrix I np.eye(3) new_cell (I strain_matrix) old_cell # 原子坐标保持分数坐标不变VASP 要求 strained_structure Structure(new_cell, structure.species, structure.frac_coords) return strained_structure # 读取原始结构POSCAR orig_struct Structure.from_file(POSCAR) # 定义 6 个独立应变模式单位小数如 0.01 表示 1% strain_patterns [ np.diag([0.01, 0, 0]), # ε11 1% np.diag([0, 0.01, 0]), # ε22 1% np.diag([0, 0, 0.01]), # ε33 1% np.array([[0,0.01,0],[0.01,0,0],[0,0,0]]), # ε12 1% np.array([[0,0,0.01],[0,0,0],[0.01,0,0]]), # ε13 1% np.array([[0,0,0],[0,0,0.01],[0,0.01,0]]) # ε23 1% ] # 生成 VASP 输入POSCAR INCAR for i, eps in enumerate(strain_patterns): for sign in [1, -1]: strained apply_strain(orig_struct, sign * eps) poscar Poscar(strained) poscar.write_file(fvasp/POSCAR_strain_{i1}_{plus if sign0 else minus}) # 写 INCAR固定 ISIF3, IBRION2, NSW0单点能计算 with open(fvasp/INCAR_strain_{i1}_{plus if sign0 else minus}, w) as f: f.write(SYSTEM Strain scan\n) f.write(ISTART 0; ICHARG 2\n) f.write(ENCUT 520\n) f.write(ISMEAR 0; SIGMA 0.05\n) f.write(ISIF 3\n) # 关键允许晶胞和原子弛豫 f.write(IBRION 2; NSW 0\n) # 单点能不优化 # 生成 QE 输入pw.in for i, eps in enumerate(strain_patterns): for sign in [1, -1]: strained apply_strain(orig_struct, sign * eps) # QE 需要 cif 或 manual input这里用 pymatgen 生成 pw.in qe_input QeInput( structurestrained, pseudo_dir./pseudopotentials, pseudopotentials{Si: si_lda_v1.uspp.F.UPF}, # 替换为你的元素 control{ calculation: scf, restart_mode: from_scratch, prefix: fstrain_{i1}_{plus if sign0 else minus}, outdir: ./out, pseudo_dir: ./pseudopotentials }, system{ ecutwfc: 60, ecutrho: 240, occupations: smearing, smearing: gaussian, degauss: 0.01 }, electrons{conv_thr: 1e-8} ) qe_input.write_file(fqe/pw_strain_{i1}_{plus if sign0 else minus}.in)逻辑说明apply_strain函数实现拉格朗日应变确保 VASP 和 QE 使用完全相同的应变张量定义。ISIF3是 VASP 中唯一能同时优化晶胞和原子坐标的标志QE 的cell_dofree在pw.in中由QeInput自动处理。此脚本生成 12 组输入每组含 VASP 的POSCARINCAR和 QE 的pw.in目录结构清晰隔离。3.3 批量提交脚本submit_all.sh统一调度双引擎任务#!/bin/bash # submit_all.sh # VASP 任务使用 sbatch 提交适配 Slurm for dir in vasp/POSCAR_strain_*; do base$(basename $dir | sed s/POSCAR_//) cd vasp sbatch EOF #!/bin/bash #SBATCH --job-namevasp_$base #SBATCH --ntasks16 #SBATCH --time02:00:00 #SBATCH --outputvasp_%j.out mpirun -np 16 /path/to/vasp_std vasp_$base.out EOF cd .. done # QE 任务同样提交 for inp in qe/pw_strain_*.in; do base$(basename $inp .in) cd qe sbatch EOF #!/bin/bash #SBATCH --job-nameqe_$base #SBATCH --ntasks16 #SBATCH --time02:00:00 #SBATCH --outputqe_%j.out mpirun -np 16 /path/to/pw.x $inp $base.out EOF cd .. done参数说明--ntasks16匹配你的节点核数/path/to/vasp_std和/path/to/pw.x必须替换为绝对路径。此脚本确保双引擎任务并行提交避免手动逐个运行的低效。4. 应力提取与张量对齐从 OUTCAR 和 QE output 中解析 6×6 应力矩阵4.1 从 VASP OUTCAR 中提取应力避开 total-stress 与 internal-stress 混淆VASP 的OUTCAR中应力有两处total-stress总应力单位 kB和internal-stress内应力单位 eV/ų。弹性分析必须用total-stress因为它是晶胞整体受力响应对应热力学定义的柯西应力。关键陷阱total-stress块在 OUTCAR 中出现多次每次离子步一次必须取最后一步的收敛结果# parse_vasp_stress.py import re import numpy as np def parse_vasp_stress(outcar_path): with open(outcar_path, r) as f: lines f.readlines() # 找到最后一个 total-stress 块 stress_lines [] for i, line in enumerate(lines): if total-stress in line: # 取接下来的 3 行xx, yy, zz stress_lines lines[i1:i4] if not stress_lines: raise ValueError(fNo total-stress found in {outcar_path}) stress_matrix np.zeros((3,3)) for i, line in enumerate(stress_lines): # VASP 输出格式xx yy zz xy yz zx 单位 kB vals list(map(float, re.findall(r[-]?\d*\.\de[-]\d|[-]?\d\.\d, line))) if len(vals) 6: stress_matrix[i,i] vals[0] # xx, yy, zz if i 0: stress_matrix[0,1] vals[3] # xy stress_matrix[0,2] vals[4] # xz elif i 1: stress_matrix[1,0] vals[3] # yx xy stress_matrix[1,2] vals[5] # yz elif i 2: stress_matrix[2,0] vals[4] # zx xz stress_matrix[2,1] vals[5] # zy yz # 单位转换kB → GPa (1 kB 0.1 GPa) stress_gpa stress_matrix * 0.1 return stress_gpa # 示例调用 stress parse_vasp_stress(vasp/vasp_strain_1_plus.out) print(VASP stress (GPa):\n, stress)逻辑说明正则表达式r[-]?\d*\.\de[-]\d|[-]?\d\.\d精准捕获科学计数法与浮点数stress_matrix构建为对称矩阵符合应力张量性质* 0.1完成 kB→GPa 转换这是文献对比的基准单位。4.2 从 QE output 中提取应力解析stress_tensor.x输出的二进制格式QE 的stress_tensor.x输出是 Fortran 二进制不能用文本编辑器查看。必须用scipy.io.FortranFile解析# parse_qe_stress.py from scipy.io import FortranFile import numpy as np def parse_qe_stress(binary_path): # stress_tensor.x 输出先写 3x3 stress matrix (float64)再写 3x3 cell matrix f FortranFile(binary_path, r) try: # 读取 9 个 float643x3 应力 stress_flat f.read_record(dtypenp.float64) if len(stress_flat) ! 9: raise ValueError(fExpected 9 values, got {len(stress_flat)}) stress_matrix stress_flat.reshape(3,3) # QE 输出单位是 Ry/bohr³转换为 GPa: 1 Ry/bohr³ 14.7108 GPa stress_gpa stress_matrix * 14.7108 return stress_gpa except Exception as e: raise ValueError(fFailed to parse {binary_path}: {e}) finally: f.close() # 示例调用需先运行 stress_tensor.x # stress_tensor.x -in pw_strain_1_plus.out -stresstensor stress parse_qe_stress(qe/stress_tensor_1_plus.dat) print(QE stress (GPa):\n, stress)参数说明stress_tensor.x必须在 QE SCF 计算后单独运行命令为stress_tensor.x -in pw_strain_1_plus.out -stresstensor输出stress_tensor_1_plus.dat。14.7108是 Ry/bohr³ → GPa 的精确换算因子不可近似。4.3 张量对齐统一 VASP 与 QE 的 Voigt 标记顺序VASP 和 QE 的应力张量输出顺序均为[σxx, σyy, σzz, σyz, σxz, σxy] # Voigt notation (6-component vector)但弹性张量拟合需 6×6 矩阵必须将应力向量与应变向量按相同 Voigt 顺序排列。我们定义标准 Voigt 顺序IndexVoigt ComponentPhysical Meaning0σ₁₁xx stress1σ₂₂yy stress2σ₃₃zz stress3σ₂₃yz stress4σ₁₃xz stress5σ₁₂xy stressdef voigt_stress(stress_matrix): Convert 3x3 stress matrix to 6-component Voigt vector return np.array([ stress_matrix[0,0], # σ11 stress_matrix[1,1], # σ22 stress_matrix[2,2], # σ33 stress_matrix[1,2], # σ23 stress_matrix[0,2], # σ13 stress_matrix[0,1] # σ12 ]) # 示例 vasp_voigt voigt_stress(parse_vasp_stress(vasp/vasp_strain_1_plus.out)) qe_voigt voigt_stress(parse_qe_stress(qe/stress_tensor_1_plus.dat)) print(VASP Voigt:, vasp_voigt) print(QE Voigt: , qe_voigt)这一步是双引擎结果可比性的生死线。若 Voigt 顺序不一致后续最小二乘拟合将完全失效。5. 弹性张量拟合与交叉验证用 Python 实现最小二乘 QE elastics.x 双校验5.1 最小二乘拟合从 12 组 (ε, σ) 数据解出 C_ij弹性本构关系$\sigma_i \sum_j C_{ij} \varepsilon_j$其中 $i,j 1..6$。将 12 组数据堆叠为矩阵方程 $S C \cdot E$则 $C S \cdot E^$$E^$ 为 E 的伪逆# fit_elastic_tensor.py import numpy as np def fit_elastic_tensor(strain_voigt_list, stress_voigt_list): strain_voigt_list: list of 12 arrays, each shape (6,) stress_voigt_list: list of 12 arrays, each shape (6,) Returns: C (6,6) elastic tensor in Voigt notation E np.array(strain_voigt_list) # (12, 6) S np.array(stress_voigt_list) # (12, 6) # Solve C S pinv(E) C S.T np.linalg.pinv(E.T) return C # 加载所有应变-应力对 strains [] stresses [] for i in range(1, 7): # 6 种模式 for sign in [1, -1]: # 读取 VASP 应力或 QE此处以 VASP 为例 stress voigt_stress(parse_vasp_stress(fvasp/vasp_strain_{i}_{plus if sign0 else minus}.out)) # 应变向量对应 strain_patterns[i-1] eps strain_patterns[i-1].flatten()[:6] # 取 Voigt 应变 if sign 0: eps -eps strains.append(eps) stresses.append(stress) C_vasp fit_elastic_tensor(strains, stresses) print(Fitted C (VASP, GPa):\n, C_vasp)逻辑说明np.linalg.pinv计算 Moore-Penrose 伪逆鲁棒处理 126 的超定系统输出C_vasp是 6×6 矩阵单位 GPa可直接与文献值对比。5.2 QE elastics.x 交叉验证用官方工具反向生成应变-应力并比对QE 的elastics.x不仅能拟合还能反向生成理论应力。我们将C_vasp导出为 QE 格式用elastics.x计算理论应力与 VASP 实际应力比对# export_to_qe.py def export_C_to_qe(C_matrix, filenameelastic_tensor.dat): Export 6x6 C matrix to QE elastics.x input format with open(filename, w) as f: f.write(6\n) # number of components for i in range(6): for j in range(6): f.write(f{C_matrix[i,j]:.6f} ) f.write(\n) export_C_to_qe(C_vasp, elastic_tensor_vasp.dat)然后运行 QE 工具# 生成理论应力需先准备 elastics.in cat elastics.in EOF elastics prefix test outdir ./out fildyn dyn flfrc force_constants elas .true. etens .true. / EOF elastics.x elastics.in elastics.outelastics.out中会输出理论应力与 VASP 的OUTCAR应力逐点比对。若最大相对误差 2%说明拟合可靠。5.3 避坑应力-应变拟合的 4 个致命陷阱与血泪解决方案常见问题应力-应变拟合结果发散、弹性常数为负、与文献值偏差巨大现象拟合出的 $C_{11}$ 为负值原因应变幅度过大3%进入非线性区或结构未充分弛豫NSW0但IBRION设置错误解决重跑 ±0.5%、±1%、±1.5% 三组用前两组线性拟合检查OUTCAR中reached required accuracy是否出现现象VASP 与 QE 拟合的 $C_{44}$ 相差 30%原因QE 的ecutwfc过低50 Ry平面波截断不足剪切应力计算失真解决对剪切应变模式ε23, ε13, ε12单独提高ecutwfc至 70 Ry 并重算现象elastics.x报错Error in reading the elastic tensor原因导出的elastic_tensor.dat文件末尾有多余空格或换行解决用sed -i s/[[:space:]]*$// elastic_tensor.dat清理确保每行恰好 6 个数字现象Python 拟合的 $C$ 矩阵不对称$C_{12} \neq C_{21}$原因Voigt 应变向量构造错误如把 σ12 放错位置或应力提取时矩阵索引颠倒解决打印strain_voigt_list[0]和stress_voigt_list[0]人工核对前 3 个分量是否为 σ11,σ22,σ33用np.allclose(C, C.T, atol1e-3)强制校验对称性6. 可视化与工程化输出绘制应力-应变曲线、杨氏模量方向图与 JSON 报告6.1 绘制 12 组应力-应变散点图验证线性度与双引擎一致性# plot_strain_stress.py import matplotlib.pyplot as plt import numpy as np # 加载所有数据 strain_vals [] # 存储所有应变值标量如 ε11 stress_vals_vasp [] stress_vals_qe [] for i in range(1, 7): for sign in [1, -1]: # 取 ε11 模式为例i1 if i 1: eps_val 0.01 * sign strain_vals.append(eps_val) vasp_stress parse_vasp_stress(fvasp/vasp_strain_{i}_{plus if sign0 else minus}.out)[0,0] qe_stress parse_qe_stress(fqe/stress_tensor_{i}_{plus if sign0 else minus}.dat)[0,0] stress_vals_vasp.append(vasp_stress) stress_vals_qe.append(qe_stress) # 绘图 plt.figure(figsize(10,6)) plt.scatter(strain_vals, stress_vals_vasp, cred, labelVASP σ11, s50) plt.scatter(strain_vals, stress_vals_qe, cblue, labelQE σ11, s50, markerx) # 线性拟合 z_vasp np.polyfit(strain_vals, stress_vals_vasp, 1) z_qe np.polyfit(strain_vals, stress_vals_qe, 1) p_vasp np.poly1d(z_vasp) p_qe np.poly1d(z_qe) plt.plot(strain_vals, p_vasp(strain_vals), r--, alpha0.7) plt.plot(strain_vals, p_qe(strain_vals), b--, alpha0.7) plt.xlabel(Strain ε₁₁) plt.ylabel(Stress σ₁₁ (GPa)) plt.title(Stress-Strain Curve: VASP vs QE) plt.legend() plt.grid(True) plt.savefig(stress_strain_curve.png, dpi300, bbox_inchestight) plt.show() print(fVASP C11 {z_vasp[0]:.2f} GPa) print(fQE C11 {z_qe[0]:.2f} GPa)效果一张图同时展示双引擎数据点、各自线性拟合线、及拟合斜率即 $C_{11}$。若两条虚线几乎重合说明双引擎收敛一致若分离则需检查泛函设置VASP 的GGAPEvs QE 的input_dftPBE。6.2 杨氏模量方向图用 pymatgen 生成 3D 极坐标图杨氏模量 $E(\mathbf{n}) \frac{1}{S_{ijkl} n_i n_j n_k n_l}$其中 $S$ 是柔度张量$C^{-1}$。pymatgen 内置ElasticTensor类可直接计算# plot_youngs_modulus.py from pymatgen.analysis.elasticity import ElasticTensor import numpy as np # 从拟合的 C_vasp 创建 ElasticTensor 对象单位 GPa et ElasticTensor(C_vasp) # 生成球面网格1000 个方向 grid et.get_grid() youngs et.youngs_modulus(grid) # 绘制极坐标图使用 matplotlib 的 3D scatter fig plt.figure(figsize(12,10)) ax fig.add_subplot(111, projection3d) # 归一化 youngs 到 [0,1] 用于颜色映射 norm_youngs (youngs - youngs.min()) / (youngs.max() - youngs.min()) scatter ax.scatter(grid[:,0], grid[:,1], grid[:,2], cnorm_youngs, cmapviridis, s1) ax.set_title(Young\s Modulus Directionality) ax.set_xlabel(x) ax.set_ylabel(y) ax.set_zlabel(z) plt.colorbar(scatter, axax, shrink0.5, aspect20, labelNormalized E (GPa)) plt.savefig(youngs_modulus_3d.png, dpi300, bbox_inchestight) plt.show()输出一个彩色球体颜色深浅表示该方向杨氏模量大小。各向异性材料如 Ti会呈现明显瓣状而立方晶系如 Si接近球形。6.3 生成标准化 JSON 报告供团队共享与 CI/CD 集成# generate_report.py import json import numpy as np report { metadata: { software: { vasp_version: 6.4.2, qe_version: 7.2, python_version: 3.10.12 }, calculation: { functional: PBE, kpoints: Gamma-centered 8x8x8, encut: 520 eV, convergence: energy 1e-6 eV/atom } }, elastic_tensor_gpa: C_vasp.tolist(), compliance_tensor_gpa_inv: np.linalg.inv(C_vasp).tolist(), derived_properties: { bulk_modulus_gpa: float(et.k_vrh), shear_modulus_gpa: float(et.g_vrh), youngs_modulus_gpa: float(et.youngs_modulus), poissons_ratio: float(et.poissons_ratio), universal_anisotropy: float(et.universal_anisotropy) } } with open(elastic_report.json, w) as f: json.dump(report, f, indent2) print(Report saved to elastic_report.json)这份 JSON 是真正的工程交付物包含完整软件栈信息、弹性张量原始数据、派生性质体模量、剪切模量等、以及各向异性指标。可被下游脚本直接json.load()也可集成进 Jenkins/GitLab CI 自动化测试。我坚持一个习惯每次新跑一个材料必先git add elastic_report.json并提交。三年下来团队已积累 200 个可追溯、可本文还有配套的精品资源点击获取