LAMMPS副本数据整理与续算

本组用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 先分清输出是在交换前还是交换后#

本版本的 temper 系列先完成这一段MD并写出dump/统计，再作该步交换，最后写全局状态映射。于是同为Step100：原始轨迹帧属于这次交换之前，全局日志Step100一行属于这次交换之后。

本例每20步交换、每100步存帧，所以Step100的原始帧使用Step80之后的映射；起始Step0使用初始映射。续算起点Step4000使用恢复后的初始行，其后的帧同样按这一规则处理。

工具没有直接把每帧配到同一步的新映射。它还单独保存最后一个交换前轨迹帧与显式交换后dump，实际核对：普通NVT/NPT的速度按 sqrt(Tnew/Told) 缩放；gREM的速度不因λ交换而改变。

先查看最后几行映射及交换后的快照：

tail -3 work/cases/nvt/log.lammps
head -12 work/cases/nvt/post-2.dump
head -10 work/analysis/replica-state-snapshots.csv

2 用工具按温度或λ整理快照#

module load ase/3.29.0
python analyze.py work
python demultiplex.py work
head -8 work/analysis/by-state/nvt/state-0.xyz

工具为六段计算、每段四个状态生成24个扩展XYZ文件。每帧保留原子ID、步数、源副本、状态编号及盒子尺寸；用于显示时将盒子原点平移至0并把坐标折回盒内，原始原点记录在 DEMULTIPLEX-VALIDATION.json，原始dump保留实际盒界。

按状态整理后的连续帧可能来自不同副本，不能直接用于MSD、VACF或扩散系数。 单个副本虽然原子轨迹连续，但温度或λ在改变，也不能自动当成固定温度的输运轨迹。XYZ在这里主要用于结构展示与构型快照整理。

3 从四份重启文件恢复#

每个副本分别有 restart-0.bin 至 restart-3.bin。全局日志最后一行是“副本→温度/λ状态”的映射，不能用最初0、1、2、3直接替代。程序自动生成以下实际输入并完成1000步续算。

NVT续算#

variable rep world 0 1 2 3
variable target world 0.7 0.8 0.9 1.0
variable active world 0 3 2 1
variable current world 0.7 1.0 0.9 0.8
read_restart ../nvt/restart-${rep}.bin
thermo 100
thermo_style custom step temp pe press vol
thermo_modify norm no format float %.16g
fix bath all nvt temp ${current} ${current} .2
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 1000 20 ${target} bath 0 12791 ${active}
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

NPT续算#

variable rep world 0 1 2 3
variable target world 0.7 0.8 0.9 1.0
variable active world 0 1 2 3
variable current world 0.7 0.8 0.9 1.0
read_restart ../npt/restart-${rep}.bin
thermo 100
thermo_style custom step temp pe press vol
thermo_modify norm no format float %.16g
fix bath all npt temp ${current} ${current} .2 iso 1 1 2
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 1000 20 ${target} bath 0 12791 1 ${active}
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

gREM续算#

variable rep world 0 1 2 3
variable target world 0.7 0.8 0.9 1.0
variable active world 0 2 1 3
variable current world 0.7 0.9 0.8 1.0
read_restart ../grem/restart-${rep}.bin
thermo 100
thermo_style custom step temp pe press vol
thermo_modify norm no format float %.16g
fix bath all nvt temp .8 .8 .2
fix ensemble all grem ${target} -.0005 -700 bath
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 f_ensemble 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/grem 1000 20 ${target} ensemble bath 0 12791 ${active}
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) $(f_ensemble:%.16g)" file post-${rep}.dat screen no
write_restart restart-${rep}.bin

target 是原有有序状态表，active 是最后映射索引，current 是各副本当前目标值。普通NVT/NPT用current重建目标；gREM按同一原始λ表初始化fix，再用active恢复正确λ分配。二进制重启应使用对应LAMMPS版本与兼容构建，不当作跨版本通用数据格式。

实际续算的初始映射与生产段最后映射一致，速度保持一致；有些原始坐标会因周期折回相差一个整数盒长，最小镜像位置保持一致。这不是原子在读重启时发生了真实跃迁。

本例续算重新初始化了交换随机数发生器，展示的是正确状态和动力学数据的恢复，不宣称与不中断计算的后续随机序列逐位一致。

4 独立检查势能、温度和压力#

analyze.py 从每帧位置和盒子独立重算截断平移LJ势能，再从速度重算动能温度，并根据两体维里重算压力。NPT小盒子变小时会纳入截断范围内额外的周期映像，避免只用简单最近镜像而漏项。该检查器假设最短盒长大于截断半径2.5，不能用于任意压缩体系。

生产段和续算分别核对；图中短轨迹的接受率只是本输入的观察，不是所有系统都应达到的标准值。

python check_results.py work
python plot_results.py work

分析与按状态整理共6490项核对通过，保留原生输入、轨迹、映射、统计和重启文件，便于追溯每个快照来自何处。短例子的数值与命令流程可以参考，材料热力学结论仍需合适体系、状态设计和更长采样。

继续：NVT/NPT副本交换、gREM广义系综。原生交换及输出说明见官方temper手册。