LAMMPS 应变应力与弹性常数
应变应力与弹性常数,包含可提交输入、实际终端命令、数据和配图。
阅读本篇目录
本文使用鸢算科技集群 LAMMPS 30 Sep 2026 的真实计算。下载输入、参考轨迹、分析脚本与数据图(约 72.5 MB),解压后进入 lammps-examples。cases 保存参考结果,inputs 保存干净输入,使用 new-case.py 创建自己的计算副本。
要完成什么#
先把 FCC 镍优化到零压力,再从同一个 relaxed.restart 分别施加 x 方向 −0.4%、−0.2%、0、0.2%、0.4% 的应变,其他方向尺寸固定。每个应变点优化内部坐标,记录 pxx、pyy 和 pzz。
LAMMPS 压力为压缩正号,绘制拉伸正应力需取负号;metal 的压力单位是 bar,乘 1e−4 转成 GPa。对 σxx–εxx 和 σyy–εxx 的斜率,得到本模型 C11 与 C12 的教学估计。
这里只扰动一个正应变分量,没有测得完整 6×6 弹性矩阵,也不能据此给出 C44。各向异性、有限温度和大应变问题,需要相应的全部扰动、热平均或非线性分析。
1 准备并提交#
python new-case.py ni-elastic
cd work/ni-elastic
sbatch run.sh
run.sh 是本算例匹配的核时 P100 或 V10016G 模板。套餐用户按自己的权限使用 yskit 生成提交脚本,保留下面的计算命令。完整模拟在计算节点执行。其他例子回到资料包根目录,替换目录名后分别提交。
| 目录 | 内容 |
|---|---|
ni-elastic |
独立输入、参考输出和可重跑命令 |
2 实际计算命令与输入#
ni-elastic#
mpirun -np 1 lmp -sf gpu -pk gpu 1 -in in.lammps > stdout.log 2>&1
下面是 ni-elastic/in.lammps 的完整输入;其他目录的完整输入也在资料包中。
units metal
atom_style atomic
boundary p p p
lattice fcc 3.52
region box block 0 6 0 6 0 6
create_box 1 box
create_atoms 1 box
mass 1 58.6934
pair_style eam
pair_coeff * * Ni_u3.eam
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
thermo 100
thermo_style custom step time atoms temp pe ke etotal press vol density
thermo_modify norm no format float %.12g
fix relax all box/relax iso 0 vmax .001
minimize 1e-12 1e-8 3000 20000
unfix relax
write_restart relaxed.restart
clear
units metal
atom_style atomic
read_restart relaxed.restart
pair_style eam
pair_coeff * * Ni_u3.eam
change_box all x scale 0.996 remap
minimize 1e-12 1e-8 2000 10000
run 0
print "-0.004 $(pxx:%.12g) $(pyy:%.12g) $(pzz:%.12g) $(pe:%.12g)" append elastic.dat screen no
clear
units metal
atom_style atomic
read_restart relaxed.restart
pair_style eam
pair_coeff * * Ni_u3.eam
change_box all x scale 0.998 remap
minimize 1e-12 1e-8 2000 10000
run 0
print "-0.002 $(pxx:%.12g) $(pyy:%.12g) $(pzz:%.12g) $(pe:%.12g)" append elastic.dat screen no
clear
units metal
atom_style atomic
read_restart relaxed.restart
pair_style eam
pair_coeff * * Ni_u3.eam
change_box all x scale 1 remap
minimize 1e-12 1e-8 2000 10000
run 0
print "0 $(pxx:%.12g) $(pyy:%.12g) $(pzz:%.12g) $(pe:%.12g)" append elastic.dat screen no
clear
units metal
atom_style atomic
read_restart relaxed.restart
pair_style eam
pair_coeff * * Ni_u3.eam
change_box all x scale 1.002 remap
minimize 1e-12 1e-8 2000 10000
run 0
print "0.002 $(pxx:%.12g) $(pyy:%.12g) $(pzz:%.12g) $(pe:%.12g)" append elastic.dat screen no
clear
units metal
atom_style atomic
read_restart relaxed.restart
pair_style eam
pair_coeff * * Ni_u3.eam
change_box all x scale 1.004 remap
minimize 1e-12 1e-8 2000 10000
run 0
print "0.004 $(pxx:%.12g) $(pyy:%.12g) $(pzz:%.12g) $(pe:%.12g)" append elastic.dat screen no
3 查看输出与实际数据#
tail -n 30 stdout.log
ls -lh
stdout.log 是计算输出;轨迹保存在 trajectory.lammpstrj(适用的算例),结构保存在 final.data。静态扫描和 run 0 对照不一定生成动力学轨迹。以下动态均值取最后一个热力学输出块的后半段;静态结果取对应输出,单位随输入确定。
| 算例 | 原子数 | 平均温度 | 平均密度 | 最后势能 |
|---|---|---|---|---|
ni-elastic |
864 | 0 | 8.94114 | -3844.80011 |
LJ/珠簧例子使用约化单位;metal 为 K、g/cm³、eV;水的 real 单位为 K、g/cm³、kcal/mol。势能列是整个体系的总量,跨不同原子数不可直接相比。
本次拟合得到 C11 = 233.27 GPa、C12 = 154.30 GPa。输入 metal 中应力的 bar→GPa 转换和负号已在分析脚本中实现。

4 提取数据与重画#
回到附件根目录,使用已经安装的 ASE Python 环境执行真实数据分析脚本:
module load ase/3.29.0
python analyze.py --root .
ls analysis
ls figures
脚本使用 NumPy/Matplotlib,读取参考日志、轨迹和相关函数,生成 CSV 与 PNG。查看 analyze.py 中本篇对应目录的处理段,可以把路径换成你的新结果。完整脚本按本专题的全部参考目录组织,单个新算例不需要假造其他目录来运行;提取对应段即可。