LAMMPS gREM广义副本交换

本组用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 λ不是动能温度#

本例动能恒温浴固定为0.8，四个广义状态的λ分别为0.7、0.8、0.9、1.0。fix grem 定义 Teff=λ+η(H−H0)，η取−0.0005、H0取−700。这里使用NVT，所以对应的H项在这个实现设置中是势能E；不要额外把动能再加进去。

不同λ改变构型采样所用的广义系综。Teff、λ和动能温度是三个不同的量，不能把广义态的结构简单当作处于温度λ的正则样本。

2 完整输入与输出#

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 .8 .8 .2
fix ensemble all grem ${target} -.0005 -700 bath
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 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 4000 20 ${target} ensemble 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) $(f_ensemble:%.16g)" file post-${rep}.dat screen no
write_restart restart-${rep}.bin

fix bath ... temp .8 .8 控制动能浴；fix ensemble ... grem 施加广义系综；temper/grem 提议交换λ状态。原子坐标仍留在各自副本。不同于普通温度交换，这里交换λ时不按λ比重新缩放速度，因为动能浴没有从0.8改成另一个温度。

程序在每100步输出实际温度、总势能、压强、体积及 f_ensemble 有效温度。format " %.16g" 保留精度并显式保留列间空格，便于直接用NumPy或文本工具读取。

head -8 work/cases/grem/stats-0.dat
tail -5 work/cases/grem/log.lammps
cat work/cases/grem/post-0.dat

3 独立核对有效温度#

分析脚本使用这个快照在交换前对应的λ，与原生势能代入公式核对。下面直接查看实际CSV中的值：

module load ase/3.29.0
python - <<'PY'
import csv
errors=[]
for row in csv.DictReader(open('work/analysis/replica-state-snapshots.csv')):
    if row['stage']!='grem':
        continue
    predicted=float(row['target'])-0.0005*(float(row['energy'])+700)
    errors.append(abs(predicted-float(row['effective_temperature'])))
print('Maximum effective-temperature discrepancy:',max(errors))
PY

实际生产段Teff约0.6353–0.9275，始终为正；动能温度则围绕0.8波动。图中也对照NPT副本各自的体积变化：

4 观察实际交换与参数限制#

这组gREM输入演示λ交换，但最低λ状态在短生产段中没有充分游走，不能据此声称相变采样已经有效。实际体系的η、H0、λ间隔要根据能量分布与目标区间选择，并检查相邻分布重叠和状态往返；本例参数不是可直接迁移到任意材料的通用设置。

广义系综数据要用于正则物理量时，需要适当重加权和充分统计。本篇仅整理实际输入、有效温度和交换记录，没有计算重加权自由能或相变点。

5 gREM续算时需要重建什么#

fix grem 不把自身λ、η、H0等信息写入二进制重启，因此续算要重建恒温与grem设置，并恢复最后的状态映射。示例工具自动读取全局日志最后一行。

fix grem 的初始λ与 temper/grem 入参使用同一有序λ表，真正恢复的分配通过末尾 active 索引指定；不要把两处参数都预先重复重排。当前版本的正确输入见续算教程。

官方说明：temper/grem、fix grem。