LAMMPS三阶力常数实操

势能展开中的三阶项可以描述局部非谐性。本篇用一个可解析的Morse两原子模型，实际运行third_order、核对全部216个笛卡尔分量，并用独立总能差分检查其中一个分量。

下载三阶张量输入、原始稀疏文件、独立能量差分与解析核对脚本，解压进入run.sh目录。

1 提交小体系例子#

sbatch submit.slurm

核时P100模板1GPU1CPU，使用P100构建中的普通CPU路径，不启用GPU或KOKKOS后缀。套餐用户用yskit生成所属分区脚本，在已分配资源中执行bash run.sh；不要在登录节点模拟。已有work请重新解压。

版本LAMMPS30Sep2026、GNU C++13.3.0/C++17、Open MPI5.0.10；分析ASE3.29.0。三个张量运行分别是regular 0.001Å、regular 0.0001Å和eskm 0.0001Å；另有一次四点能量扫描。inputs是新输入，reference是本文原始实算，work是新运行。

2 三阶张量与模型#

采用能量定义Phi3(iα,jβ,kγ)=∂³E/(∂xiα∂xjβ∂xkγ)。两个原子各有3个坐标，因此本例稠密张量为6×6×6，共216个分量。

Morse参数D0=0.5kcal/mol、α=2Å⁻¹、r0=1.6Å，两个原子沿x轴放置在0与1.6Å。距离在势阱最低点附近，径向E‴(r0)=−6D0α³=−24kcal/mol/ų。改变1号与2号原子的位置会带来相反的距离导数符号，需要保留原子及坐标顺序。

两种质量设为1和2，仅用于检查输出定义；该能量三阶导数不按质量加权。它是解析教学模型，没有替代真实晶体的三阶力常数或热导率计算。

3 regular完整输入#

units real
atom_style atomic
atom_modify map array
boundary f f f
region box block -10 10 -10 10 -10 10 units box
create_box 2 box
mass 1 1.0
mass 2 2.0
create_atoms 1 single 0 0 0 units box
create_atoms 2 single 1.6 0 0 units box
pair_style morse 4.5
pair_coeff * * 0.5 2.0 1.6
neighbor 0.3 bin
thermo_style custom step atoms pe
thermo_modify norm no format float %.16g
run 0
third_order all regular 0.0001 file third-order.dat binary no
run 0
write_dump all custom final.dump id type x y z fx fy fz modify sort id format float %.16g

atom_modify map array在建盒前启用原子ID映射，这是此命令需要的输入前置设置。third_order all regular ...对all组求张量；file指定输出，binary no选择可直接阅读的文本。

命令完成后再次run 0，在恢复的原子位置上重算基准能量和力，再保存final.dump。原始张量文件与最后基准结构分开读取。

4 先固定输出的符号约定#

head -n 8 work/cases/regular-0p0001/third-order.dat

当前30Sep2026构建的原生文件来自原子力对两个坐标的二阶差分。由于F=−∂E/∂x，按本篇能量定义的Phi3为原生值的相反数。本例原生x2,x2,x2约+24，转换为能量三阶导数后约−24kcal/mol/ų。

本包analyze.py保留原生值并明确生成相反数，不直接覆盖原始文件。它用可解析模型检查这一关系；更换程序版本、输出格式或外部软件时，先核对定义，不能只凭“third order”文件名就忽略符号。

5 步长对照与真实误差#

两个regular步长的误差约随步长平方下降；0.0001Å在本例更准确，但不能据此给任意材料固定一个最小步长。实际任务还要权衡力的数值噪声、截断和位移大小。

6 独立总能核对#

cat work/cases/energy-control/energy.dat

energy-control在x2=r0±h、r0±2h求四个总能，h=0.001Å。中心差分为[E(+2h)−2E(+h)+2E(−h)−E(−2h)]/(2h³)，本次实际得到约−24.000120kcal/mol/ų，与解析−24和已转换张量的对应分量一致。

这项检查独立于third_order的力差分，可以用来确认能量定义的符号与单位。对更大体系可先检查一个明确分量，再扩展全部索引和目标软件映射。

提取、转换与重画#

module load ase/3.29.0
python analyze.py work
python plot_results.py work
python check_results.py work
head -n 12 work/analysis/tensor-components.csv
ls work/figures

CSV同时保留原生归一化值、按能量定义的Phi3和解析值；JSON保留完整6×6×6张量、两个步长误差及独立能量差分。查看本文数据时将work改成reference。

原生命令参数见third_order官方说明，具体文件读取和ESKM换算继续看张量格式与单位。