LAMMPS温度与NPT副本交换

本组用108个Lennard-Jones粒子演示四副本交换：temper 做NVT温度交换，temper/npt 做NPT交换，temper/grem 做广义系综λ交换。每种方法都实际运行了4000步交换过程，并从四份重启文件续算1000步。单位采用约化 units lj，没有把这些粒子标成某种真实材料。

下载完整输入、重启实例、原始轨迹与三个分析工具，进入个人工作目录：

mkdir replica-demo
cd replica-demo
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-replicas-wave80-2d54dd6ee2e4.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

提交脚本申请4个CPU核和一个P100资源，顺序完成三种方法及续算。本例是CPU/MPI路径，没有GPU力计算或加速比。模块为 lammps/30Sep2026-p100-gpu（30Sep2026、GCC13.3、OpenMPI5.0.10、CUDA12.8），后处理用 ase/3.29.0。

计算结束后，在整包目录执行：

module purge
module load ase/3.29.0
python analyze.py work
python demultiplex.py work
python plot_results.py work
head -8 work/analysis/replica-state-snapshots.csv

已有结果位于 reference/，可把三个命令中的 work 换成 reference 先看数据。新运行请换一个解压目录，不覆盖既有work。analyze.py 核对交换映射、独立LJ能量/维里压力、动能温度和续算；demultiplex.py 按温度或λ状态整理快照；plot_results.py 出四张原图。

1 副本、进程和温度编号分别是什么#

四个副本分别在四个MPI分区中运行，每个分区1个MPI进程。原子坐标留在各自副本；交换发生时，副本所对应的温度状态会改变，速度随温度比进行缩放。副本编号0–3与温度编号0–3因此不能长期当成同一件事。

variable target world .7 .8 .9 1.0 给四个初始温度。每次交换都在相邻温度状态间提议，而不是固定让副本0永远与副本1交换。

2 NVT完整输入#

variable rep world 0 1 2 3
variable target world .7 .8 .9 1.0
variable seed world 10291 29381 38479 48571
units lj
atom_style atomic
boundary p p p
lattice fcc .84
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 1 2.5
pair_modify shift yes
neighbor .3 bin
neigh_modify every 1 delay 0 check yes
velocity all create ${target} ${seed} mom yes rot no dist gaussian
timestep .002
thermo 100
thermo_style custom step temp pe press vol
thermo_modify norm no format float %.16g
fix bath all nvt temp ${target} ${target} .2
run 1000
reset_timestep 0
variable measured equal temp
variable E equal pe
variable P equal press
variable V equal vol
fix stats all ave/time 100 1 100 v_measured v_E v_P v_V file stats-${rep}.dat format " %.16g"
dump frames all custom 100 trajectory-${rep}.dump id type x y z vx vy vz
dump_modify frames sort id format float %.16g
temper 4000 20 ${target} bath 0 12791
write_dump all custom post-${rep}.dump id type x y z vx vy vz modify sort id format float %.16g
run 0
print "POST ${rep} $(temp:%.16g) $(pe:%.16g) $(press:%.16g) $(vol:%.16g)" file post-${rep}.dat screen no
write_restart restart-${rep}.bin

先以不同随机种子初始化速度，NVT平衡1000步，再把计步器设为0，运行4000步副本交换。temper 4000 20 每20步尝试交换；seed1=0 轮流采用两种奇偶配对，第二个种子用于Metropolis判断。不是每次提议都会接受。

原生四分区启动关系如下；整包提交脚本已经处理，不必重复手动启动：

mpirun --bind-to none -np 4 lmp -partition 4x1 -in in.lammps

图像/副本数与MPI资源保持一致，world 变量也要有对应数量的取值。每个副本的轨迹、统计与重启文件都带自己的编号，避免覆盖。

3 NPT完整输入#

variable rep world 0 1 2 3
variable target world .7 .8 .9 1.0
variable seed world 10291 29381 38479 48571
units lj
atom_style atomic
boundary p p p
lattice fcc .84
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 1 2.5
pair_modify shift yes
neighbor .3 bin
neigh_modify every 1 delay 0 check yes
velocity all create ${target} ${seed} mom yes rot no dist gaussian
timestep .002
thermo 100
thermo_style custom step temp pe press vol
thermo_modify norm no format float %.16g
fix bath all npt temp ${target} ${target} .2 iso 1 1 2
run 1000
reset_timestep 0
variable measured equal temp
variable E equal pe
variable P equal press
variable V equal vol
fix stats all ave/time 100 1 100 v_measured v_E v_P v_V file stats-${rep}.dat format " %.16g"
dump frames all custom 100 trajectory-${rep}.dump id type x y z vx vy vz
dump_modify frames sort id format float %.16g
temper/npt 4000 20 ${target} bath 0 12791 1
write_dump all custom post-${rep}.dump id type x y z vx vy vz modify sort id format float %.16g
run 0
print "POST ${rep} $(temp:%.16g) $(pe:%.16g) $(press:%.16g) $(vol:%.16g)" file post-${rep}.dat screen no
write_restart restart-${rep}.bin

与NVT相比，这里使用 fix npt 和 temper/npt。四个副本都以压强1为目标，温度不同，盒子分别调整；本例没有交换不同压强，不能把它写成已经覆盖全部温压组合。NPT的接受判据与NVT不同，不要简单把 temper 搭在任意变体积流程上。

4 读取实际交换记录#

head -10 work/cases/nvt/log.lammps
tail -5 work/cases/nvt/log.lammps
head -8 work/cases/npt/stats-0.dat

全局日志 Step T0 T1 T2 T3 的列是副本，值是这个副本当前对应的温度编号，不是瞬时温度本身。分区日志 log.lammps.0 等和 stats-0.dat 记录实际温度、总势能、压强和体积。

每种颜色是一个状态编号。NVT例中副本会访问不同温度；这组NPT参数的交换较少，不能把程序能够交换等同于采样已经充分。

左图是各副本的连续轨迹采样温度，叉号是该输出时刻对应的目标状态。右图按温度状态整理，数据可能来自不同副本；瞬时动能温度会围绕目标波动，不会每一步等于设定值。

5 怎样判断状态间隔是否合适#

4000步、每20步一次，共200轮；奇数配对轮有两个相邻状态对，另一轮有一个，因此生产段有300次状态对提议。表中不是“接受轮数/200”。

本例NPT仅接受22/300个状态对，短计算不足以给出收敛热力学统计。实际应用应检查相邻温度的势能/焓分布重叠、各副本的温度游走和往返次数，再调整状态间隔、增加副本或延长采样；不要只追求某个统一的接受率数值。

本例盒子和轨迹都很小，只作命令与数据流程示范。小体系、短轨迹中的能量或体积变化不能用来宣称熔点、相变性质或采样增益已经收敛。

继续：gREM有效温度与λ交换、按状态整理与正确续算。官方说明：temper、temper/npt。