LAMMPS EOS与晶格体积实操
终端提取九点EOS、使用ASE拟合晶格与体积模量，独立检查能量导数和原生压力。
阅读本篇目录
获得一组体积与能量后，可以用状态方程拟合平衡体积和体积模量。本篇直接使用前两篇的五个真实势模型，说明从原生文件到ASE拟合、检查拟合范围，再用独立差分对照压力的步骤。
下载五种势模型输入、冻结势文件、实际EOS和压力数据、ASE分析脚本，解压进入run.sh目录。
1 提交整套流程#
sbatch submit.slurm
模板申请核时P100分区1GPU1CPU，顺序完成五种势的普通EOS、三种EAM的GPU对照及五组普通路径压力差分。共有8个EOS后端运行、每个9个体积点，另有5个压力检查、每个3点，共87次固定结构能量与压力求值；按五个势模型工作流组织。
套餐用户用yskit生成所属分区脚本，在已分配资源中执行bash run.sh；不要在登录节点计算。多卡仅使用分配的CUDA可见设备和CPU亲和范围，每卡一个势模型工作流，每程序1MPI1OMP。已有work请重新解压。
LAMMPS30Sep2026、GNU C++13.3.0/C++17、Open MPI5.0.10、CUDA GPU包混合精度；分析ASE3.29.0。普通路径在同一P100节点运行，不是CPU专用节点的性能数据。inputs保存新输入、DATA为本例势文件、reference为实际原始输出、work为新计算。
2 先看原始文件与单位#
cat work/cases/eam-funcfl-cpu/eos.dat
cat work/cases/eam-fs-cpu/eos.dat
每行EOS后依次为常规格常数a、整盒体积V、整盒势能E、原生压力P和原子数N。metal单位下a为Å、V为ų、E为eV、P为bar；thermo_modify norm no表示没有按原子数归一化。
拟合前将V、E分别除以N。FCC每原子体积为a³/4、BCC为a³/2、金刚石为a³/8，用相应关系把拟合V₀换回常规晶胞格常数，不能把每原子体积直接开立方就当a。
3 用ASE完成拟合#
module load ase/3.29.0
python - <<'PY'
import numpy as np
from ase.eos import EquationOfState
from ase.units import GPa
t = np.loadtxt('work/cases/eam-funcfl-cpu/eos.dat', usecols=(1,2,3,4,5))
n = t[0,4]
v = t[:,1] / n
e = t[:,2] / n
eos = EquationOfState(v[2:7], e[2:7], eos='birchmurnaghan')
v0, e0, bulk = eos.fit()
print('V0 per atom / A^3:', v0)
print('FCC lattice a0 / A:', (4*v0)**(1/3))
print('B / GPa:', bulk/GPa)
PY
脚本取九点中的中央五点，即基准格常数±2%范围；本包分析还对全部九点另做拟合，保留两种结果与残差。要用自己的体系，先确认采样包住最低点，再选择合适范围和方程，不直接截取任意五行。
4 实际局部拟合结果#
| 势模型 | 元素/结构 | 局部拟合a₀/Å | 局部拟合B/GPa | 九点拟合B/GPa |
|---|---|---|---|---|
| eam-funcfl | Ni/fcc | 3.520001 | 180.4987 | 180.1922 |
| eam-alloy | Ni/fcc | 3.519997 | 180.8931 | 181.7480 |
| eam-fs | Fe/bcc | 2.855236 | 176.5702 | 169.7258 |
| adp | Al/fcc | 4.049999 | 77.6742 | 77.9154 |
| meam | Si/diamond | 5.431001 | 97.5623 | 97.3951 |
这些是当前势、当前结构、静态体积路径和所选拟合方程得到的参数，不是实验值，也没有包含温度或零点振动。九点与五点结果不同，尤其Fe例子的B对范围较敏感；后续应按所需精度细化体积点、检查残差和拟合方式。
5 独立计算能量导数与压力#
压力应与固定原子缩放路径下的−dE/dV相符。本包在每种势基准格常数的0.99倍附近，再独立计算a−0.0001、a、a+0.0001Å三个点，使用两端整盒E与V构造差分，并和中心点原生virial压力比较。
一个完整压力控制输入为：
units metal
atom_style atomic
boundary p p p
lattice fcc 3.52
region box block 0 4 0 4 0 4
create_box 1 box
create_atoms 1 box
pair_style eam
pair_coeff * * ../DATA/Ni_u3.eam
neighbor 0.3 bin
thermo_style custom step atoms pe vol press
thermo_modify norm no format float %.16g
variable a index 3.4847 3.4848 3.4849
label scan
change_box all x final 0 $(4*v_a) y final 0 $(4*v_a) z final 0 $(4*v_a) remap units box
run 0
print "EOS ${a} $(vol:%.16g) $(pe:%.16g) $(press:%.16g) $(atoms:%.0f)" append eos.dat screen no
write_dump all custom atoms-${a}.dump id type x y z fx fy fz modify sort id format float %.16g
next a
jump SELF scan
cat work/cases/eam-funcfl-pressure/eos.dat
cat work/analysis/pressure-check.csv
| 势模型 | 原生压力/bar | 能量差分压力/bar | 差分−原生/bar |
|---|---|---|---|
| eam-funcfl | 58029.468034 | 58029.482656 | 0.014621 |
| eam-alloy | 58112.670800 | 58112.684132 | 0.013333 |
| eam-fs | 45704.925423 | 45704.927537 | 0.002114 |
| adp | 23863.274058 | 23863.275268 | 0.001210 |
| meam | 31373.163374 | 31373.167513 | 0.004140 |
本例五组差值均小于0.02bar。分析用ASE的bar常数将eV/ų转换成bar；若只把GPa和bar混用，数值会相差10000倍。
这个导数检查针对原生能量与原生压力的一致性；状态方程拟合本身仍是近似。拟合曲线的压力导数不必在宽体积区间逐点精确等于原生压力，因此还要保留拟合残差和范围对照。
6 换成自己的模型#
先确定势文件、元素映射与初始结构，再在合理体积范围固定形状扫描；逐点检查原子数、实际体积和能量。采样包住最低点后，用ASE拟合并检查范围依赖。合金、缺陷或低对称模型通常还需要在固定晶胞下松弛内部原子，任意晶胞形状的优化则需要相应自由度，不能只沿本例的一个立方格常数扫描。
提取与重画#
module load ase/3.29.0
python analyze.py work
python plot_results.py work
python check_results.py work
cat work/analysis/eos-fits.csv
cat work/analysis/pressure-check.csv
eos-fits.csv包含局部五点与完整九点两种拟合的晶格常数、每原子体积、能量、体积模量和拟合残差；pressure-check.csv保存原生压力与独立能量差分。JSON还保存所有原始扫描值和GPU对照。查看本文数据时将work改成reference。