LAMMPS DPD流体与扩散取数

DPD粒子通常代表粗粒化组分，不能默认一个粒子就是一个真实原子。DPD相互作用包括保守力、耗散力和随机力；本例先关闭随机项核对两粒子的力，再运行有热涨落的216粒子流体。

1. 运行例子#

下载完整输入、运行脚本和绘图包，解压进入 lammps-mesoscale，在已分配的CPU作业内执行：

bash run.sh

也可 sbatch job.slurm 提交。CPU分区按56核整节点分配和计费；本例只用一个线程，使用 lammps/30Sep2026-8173m 的普通CPU样式、oneAPI2024.2构建，绘图用 sisl/0.16.4。本例不是GPU加速测试。CPU队列处于维护时请待恢复后提交，不在登录节点做生产模拟。

脚本按输入逐项运行，在 raw/ 保留原生输出，在 results/ 输出全数组和汇总，在 plots/ 生成本次图片。三个教程共用一个包，完整运行一次后只读结果和重画即可。

2. 先把两粒子的力看清楚#

cat inputs/dpd-approach.lmp

units lj
atom_style atomic
boundary f f f
region box block -2 2 -2 2 -2 2
create_box 1 box
create_atoms 1 single 0 0 0 units box
create_atoms 1 single 0.6 0 0 units box
mass 1 1
group left id 1
group right id 2
velocity left set 0.1 0 0 units box
velocity right set -0.1 0 0 units box
comm_modify vel yes
pair_style dpd 0 1 92873
pair_coeff * * 25 4.5
neighbor 0.2 bin
thermo_style custom step atoms pe ke
thermo_modify norm no
run 0
write_dump all custom state.dump id type x y z vx vy vz fx fy fz modify sort id format float %.17g
variable E equal pe
print "${E}" file energy.txt

pair_style dpd 0 1 92873 中温度设为0，使随机项为零；这只是分离各力项的静态检查，不是常温流体。两粒子距0.6、截断1、A=25，保守力大小为25×(1−0.6)=10。两粒子相向速度各0.1时，耗散项再增加0.144，总力大小 10.144；相背则为 9.856。三种状态的保守势能均为 2.0，耗散力不作为可积势能加进去。

3. 完整的短流体输入#

cat inputs/dpd-fluid.lmp

units lj
atom_style atomic
boundary p p p
lattice sc 3.0
region box block 0 6 0 6 0 6
create_box 1 box
create_atoms 1 box
mass 1 1
velocity all create 1.0 71943 mom yes rot no dist gaussian
comm_modify vel yes
pair_style dpd 1 1 92873
pair_coeff * * 25 4.5
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes
fix integrate all nve
compute msd all msd com yes
compute rdf all rdf 30
fix rdfout all ave/time 10 10 100 c_rdf[*] file rdf.dat mode vector
timestep 0.01
thermo 100
thermo_style custom step temp pe ke etotal c_msd[4]
thermo_modify norm no
dump tr all custom 10 trajectory.dump id type xu yu zu vx vy vz fx fy fz
dump_modify tr sort id format float %.17g
run 1000

体系有216个单位质量粒子，数密度3、温度参数1、时间步0.01，共1000步，即10个约化时间单位。所有量都用LJ约化单位；未指定映射前不能写成ps、Å或cm²/s。DPD本身包含配对恒温作用，此例再用NVE积分位置和速度，不额外套NVT恒温器。

compute msd ... com yes 去掉整体平移。轨迹输出 xu yu zu 是解缠坐标，适合计算位移；将这些坐标画出来会有粒子跨越原模拟盒，不能理解成粒子真的逃出周期体系。

4. 从终端读取MSD与RDF#

tail -n 12 raw/dpd-fluid/rdf.dat
cat results/SUMMARY.json

本次末时刻的去质心MSD为 17.6157496。对末半窗口的MSD斜率除以6，得到约 0.36299，仅作为“如何取数”的约化单位示例；10时间单位、单条轨迹不足以宣称扩散系数已收敛。

RDF图取最后100步窗口，最大距离为相互作用截断1；不是完整长程结构。程序用周期最小像独立统计粒子对，核对全部30个RDF和配位数分箱；MSD则从原生解缠坐标独立重算。

5. 重画而不重跑#

module load sisl/0.16.4
python -B verify.py --reference
python -B plot.py

固定随机种子和单进程便于复跑这个小例。改用多进程、不同粒子排序或邻居遍历时，随机序列对应关系可能改变；科学比较应看统计和守恒关系，不能要求任意并行配置都逐帧相同。