LAMMPS压控与三斜晶胞实操

需要让晶胞随压力变化时，先决定是否控温，以及盒子各个方向、倾斜自由度如何变化。本文在同一个108粒子LJ模型上实际对照NPH、各向同性NPT、各向异性NPT和三斜NPT，给出完整输入与真实体积/压力数据。

1. 下载与提交#

下载全部原生输入、实际轨迹、CPU单点对照与绘图工具：

mkdir barostat-demo
cd barostat-demo
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-thermostats-wave91-df7f1ec315d1.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

脚本申请1个P100和1个CPU核，顺序完成12个CPU计算条件、12个GPU力对照和12个Kokkos对照；GPU末态另做相同坐标CPU单点。每次复算请使用新的解压目录。直接调用lmp的命令在已分配的计算资源中执行，登录节点可做后处理。

模块为 lammps/30Sep2026-p100-gpu、lammps/30Sep2026-p100-kokkos，编译使用GCC13.3、OpenMPI5.0.10、CUDA12.8；Python数据工具使用 ase/3.29.0。本组CPU执行路径用于数值对照，不是CPU节点性能数据。

模型、初始温度和势参数见温控实操。本组dt=0.002、10000步，目标压力1.0，Pdamp=2.0；NPT另设温度1.0、Tdamp=0.2。压力和阻尼时间都采用LJ约化单位。

2. 三斜NPT完整输入#

units lj
atom_style atomic
boundary p p p
lattice fcc .7
region box block 0 3 0 3 0 3
create_box 1 box
create_atoms 1 box
mass 1 1
pair_style lj/cut 2.5
pair_coeff * * 1.0 1.0
pair_modify shift yes
neighbor .3 bin
neigh_modify every 1 delay 0 check yes
velocity all create 1.8 20261029 mom yes rot no dist gaussian
change_box all triclinic
fix integrator all npt temp 1.0 1.0 .2 tri 1.0 1.0 2.0
timestep .002
thermo 100
thermo_style custom step time temp pe ke etotal press vol density ecouple econserve
thermo_modify norm no format float %.16g
dump trajectory all custom 1000 trajectory.dump id type xu yu zu vx vy vz fx fy fz
dump_modify trajectory sort id format float %.16g
run 10000
unfix integrator
run 0
write_dump all custom physical.dump id type x y z fx fy fz modify sort id format float %.16g
print "PHYSICAL $(pe:%.16g) $(vol:%.16g)" file physical-energy.dat screen no
print "COMPLETE" file complete.dat screen no

change_box all triclinic 先把盒子表示变成三斜形式，再用 tri 允许长度和倾斜自由度参与压控。不能只给正交盒子的NPT行加 tri 而漏掉盒子前提。

单独运行：

module purge
module load lammps/30Sep2026-p100-gpu
cd inputs/npt-tri-cpu
lmp -in in.lammps
cd ../..

3. 四种输入各自控制什么#

这些fix本身包含积分，不再给同组加NVE。NPH没有温度目标，不应要求T始终等于1；它也不要求粒子普通总能量 pe+ke 恒定，体系与压控自由度交换能量。ecouple 和 econserve 提供能量交换/扩展守恒量的相关诊断，不能把它们当作不含压控器的粒子内能。

4. 读取实际密度与压力变化#

压力目标为1，不表示每个时刻的瞬时压力都精确为1。小体系压力波动通常较明显；本文有限时间内的密度变化用于说明输入如何作用，不据此给出材料平衡密度、相稳定性或压缩率。

5. 三斜盒子的体积不能照抄边界宽度乘积#

原生dump的 BOX BOUNDS xy xz yz 行保存包围盒边界和倾斜量，不是直接把三个包围盒跨度相乘。配套解析器恢复真正的晶格矩阵H，再用 det(H) 得体积，与原生thermo逐帧核对；数密度应为108/V。

导出的XYZ保留完整3×3周期晶胞。未包裹坐标会按真实原点平移到同一参考，晶胞和粒子相对位置保持；不是把倾斜盒子硬改成正交盒子。

module purge
module load ase/3.29.0
python analyze.py work
python export_states.py work
python plot_results.py work
python check_results.py work

原生输出在 work/cases/，逐时CSV、汇总和含完整周期晶胞的XYZ在 work/analysis/，6张图在 work/figures/。reference/ 是本文实际数据；未重算时，可以把以上 work 改为 reference 重现图表。X为LJ粒子占位符，XYZ坐标、晶胞和 lj_velocity 仍是LJ约化单位，不能直接当成Å、eV或真实元素。

python - <<'PYCODE'
import csv
with open('work/analysis/npt-tri-cpu-thermo.csv') as f:
    rows = list(csv.DictReader(f))
for r in [rows[0], rows[-1]]:
    print(r['time'], r['pressure'], r['volume'], r['density'], 108/float(r['volume']))
PYCODE

真实体系要结合晶体对称性、界面方向、周期性、固定区域和外加应力选择压控自由度；允许全部倾斜并不总是合适。先核对平均体积、应力分量和时间步/阻尼，再做更长采样。GPU输入和精度范围见后端对照。

原生压力选项：fix npt/nph。