LAMMPS并行副本动力学实操

本组使用LAMMPS30Sep2026的原生REPLICA功能：四个粒子中，1号粒子沿x跨越解析双阱，另外三个局域谐振子作为陪伴自由度。所有粒子还在各自的y/z势阱内运动。使用LJ约化单位，图中X是占位标签，不代表某种真实元素；能垒1不是1 eV。

模型为1号粒子的 U=(x²−1)²+5(y+3)²+5z²，其他粒子为 U=5(x+1)²+5[y−(2id−5)]²+5z²。初态x=−1，另一稳态x=+1，解析能垒为1。这里用来学习输入、事件和时钟；材料扩散率还需要适合材料的势函数、事件定义、去相关和统计收敛。

下载原生输入、实际输出、CSV和图片，解压到个人目录：

mkdir prd-tad-demo
cd prd-tad-demo
curl -fL https://yuansuan.top/assets/downloads/full-function/lammps-rare-event-wave83-877009b7a3cb.zip -o examples.zip
unzip examples.zip
sbatch submit.slurm

提交脚本会依次运行各个小算例，然后自动提取数据、导出结构和画图。使用 lammps/30Sep2026-p100-gpu 基础模块，原生GCC13.3、OpenMPI5.0.10、CUDA12.8构建；后处理用 ase/3.29.0。这些命令走CPU/MPI路径，并不因模块名字含gpu而自动使用GPU力计算。脚本申请5个MPI进程和一个P100配额；PRD/TAD各自使用所需的分区数。

完成后，在解压目录运行：

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

已有实算结果位于 reference/，想先读图时可把后处理参数work换为reference。新计算请重新解压到新目录，保留自己的既有结果。

1 先认识四个粒子和两个势阱#

左图是本模型已知的解析能量，右图的实点来自原生初态坐标，空心点仅标出1号粒子的另一稳态。三个谐振子不会发生另一种材料相变；模型只有一个激活坐标。

2 四副本完整输入#

包中 inputs/prd-4/in.lammps 如下。pair_style zero 提供邻域框架，实际能量和力由addforce中的解析表达式给出；fix_modify field energy yes 将对应势能计入pe。

variable replica uloop 4
units lj
atom_style atomic
atom_modify map array
boundary f f f
region box block -5 5 -8 8 -5 5 units box
create_box 1 box
mass 1 1
create_atoms 1 single -1 -3 0 units box
create_atoms 1 single -1 -1 0 units box
create_atoms 1 single -1 1 0 units box
create_atoms 1 single -1 3 0 units box
pair_style zero 2.5
pair_coeff * *
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
variable dy atom y-(2*id-5)
variable active atom id==1
variable U atom v_active*(x*x-1)*(x*x-1)+(1-v_active)*5*(x+1)*(x+1)+5*v_dy*v_dy+5*z*z
variable Fx atom -v_active*4*x*(x*x-1)-(1-v_active)*10*(x+1)
variable Fy atom -10*v_dy
variable Fz atom -10*z
fix field all addforce v_Fx v_Fy v_Fz energy v_U
fix_modify field energy yes
compute physicalT all temp
compute_modify physicalT extra/dof 0
thermo 1000
thermo_style custom step c_physicalT pe ke
thermo_modify norm no format float %.16g
velocity all create 0.3 32613 mom no rot no dist gaussian temp physicalT
fix integrate all nve
variable bathseed equal 12341+17*v_replica
fix bath all langevin 0.3 0.3 0.2 ${bathseed} zero no
compute event all event/displace 1.0
timestep .005
min_style cg
dump trajectory all custom 100 trajectory-${replica}.dump id type x y z vx vy vz fx fy fz
dump_modify trajectory sort id format float %.16g
prd 50000 20 2 200 40 event 596203 min 0.0 1.0e-8 500 1000 temp .3 vel all gaussian time steps
unfix bath
run 0
write_dump all custom final-${replica}.dump id type x y z vx vy vz fx fy fz modify sort id format float %.16g
print "COMPLETE ${replica} $(step:%.0f) $(pe:%.16g)" file complete-${replica}.dat screen no

关键部分是：

prd 50000 20 2 200 40 event 596203 min 0.0 1.0e-8 500 1000 temp .3 vel all gaussian time steps

• 50000：每个副本的步数预算，采用 time steps；不包含去相位和淬火最小化的额外开销。
• 20：每隔20步淬火，比较前后稳态，检测事件。
• 2和200：每轮去相位进行两次速度随机化，每次运行200步。实际材料要检查去相关是否足够。
• 40：识别连续相关事件的时间窗口，必须是事件检查间隔的整数倍。
• event：前面 compute event all event/displace 1.0 的ID；这里按淬火后坐标判断，阈值1小于两稳态间距2。
• temp .3：去相位速度目标温度；langevin也使用0.3。种子控制随机序列，不能靠挑种子来代替统计收敛。

3 MPI分区数必须一致#

整包已自动启动。若要单独研究四副本输入，应先分配计算资源，再在解压目录建立独立运行文件夹：

module purge
module load lammps/30Sep2026-p100-gpu
mkdir manual-prd4
cp inputs/prd-4/in.lammps manual-prd4/
cd manual-prd4
mpirun -np 4 lmp -partition 4x1 -in in.lammps
cd ..

variable replica uloop 4 与4x1匹配；每个副本一个MPI进程。两副本输入在 inputs/prd-2/，对应2x1。不要只改mpirun进程数而保留不匹配的分区和变量。

4 看实际事件和时钟#

cat work/cases/prd-2/launcher.log
cat work/cases/prd-4/launcher.log
head -8 work/analysis/prd-4-events.csv

Clock是累计方法时钟，原生列以步为单位；乘dt=0.005换成约化时间。表中Clock属于最后一条事件记录，并不是50000步结束时的总Clock。 每个副本的Step与聚合Clock不是同一列。

原生输出有时重复打印同一个事件编号。配套脚本只在展示时合并完全相同的事件记录，原始日志保留。不要按日志行数统计事件，也不要把12/5当成四副本相对两副本的加速比；这只是两段有限随机样本。

5 导出真实稳态结构#

python export_events.py work
head -8 work/analysis/prd-4-events.xyz

导出的XYZ依次包含初态和12个实际事件稳态，共13帧，带particle_id、active及对应Clock信息。X只是占位物种，坐标使用LJ约化长度。它不是等时间间隔MD轨迹，不能直接据此算MSD/VACF。

原理、参数和适用条件见 LAMMPS PRD说明 和 event/displace。