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 结果目录怎么读#

2 用命令看事件数和低温时钟#

head -8 work/analysis/prd-2-events.csv
head -8 work/analysis/prd-4-events.csv
cat work/analysis/tad-5-events.csv

也可以直接读提取后的汇总：

python - <<'PYDATA'
import json
from pathlib import Path
for r in json.loads(Path('work/analysis/RESULTS.json').read_text())['results']:
    print(r['name'], r.get('events', r.get('executed_events')))
PYDATA

本次输出依次是 prd-2 5、prd-4 12、tad-5 1。这是有限样本中的事件计数，不是速率常数。

3 从原生能量核对力#

本包另有一个真实原生控制输入，将1号粒子置于 (-0.7,-2.95,0.03)，分别沿x/y/z移动±1e-5，用六个能量单点核对中心位置的力。

cat work/cases/control/cartesian-energy-control.dat
cat work/analysis/native-energy-force-control.csv

这里做的是笛卡尔坐标位移差分。中心原生力为 (-1.428,-0.5,-0.3)，分别对应负坐标能量导数。可直接计算：

python - <<'PYDATA'
import csv
from pathlib import Path
p=Path('work/cases/control/cartesian-energy-control.dat')
a={r[0]:list(map(float,r[1:])) for r in (x.split() for x in p.read_text().splitlines())}
for i,k in enumerate('xyz'):
    numeric=-(a[k+'plus'][0]-a[k+'minus'][0])/(2e-5)
    print(k, 'native=', a['BASE'][i+1], 'energy-FD=', numeric)
PYDATA

此例核对输入的力/能量一致性；它不能替代材料势函数的训练和验证。

4 导出并读取实际稳态XYZ#

python export_events.py work
python - <<'PYDATA'
from ase.io import read
for name in ['prd-2','prd-4','tad-5']:
    images=read('work/analysis/'+name+'-events.xyz', index=':')
    print(name, len(images), images[0].positions[0])
PYDATA

本次分别得到6、13、2帧，包含初态。坐标来自原生事件快照，保留粒子ID与激活自由度标记；最后的热态终点不混进淬火稳态序列。坐标使用约化单位，不能直接按Å解释。

图像和事件CSV的间隔不均匀；PRD还涉及副本同步，TAD还涉及反射。它们适合查看跃迁顺序和方法时钟，不可直接当作固定时间间隔轨迹去算MSD、VACF或振动频谱。

5 重新出图#

python plot_results.py work
ls work/figures

新运行生成四张PNG；已有参考图在reference/figures。绘图保留实际事件点和原生NEB点，不对事件数拟合或把曲线外推为材料扩散率。

继续看 PRD输入 和 TAD输入。