手机版

LAMMPS扩散系数计算 — 从均方位移到输运性质的实战方法

发布时间：2026-07-15 来源：科研学术网

小中大

字号：

LAMMPS扩散系数计算是分子动力学后处理分析中最常见的操作之一。扩散系数表征了粒子在体系中的迁移能力，对液体输运、离子传导、药物释放等过程的研究至关重要。LAMMPS扩散系数计算的核心原理是Einstein关系：在长时极限下，均方位移（MSD）与时间呈线性关系，斜率与扩散系数成正比。但”长时极限”到底多长、线性区间怎么选、周期性边界怎么处理，这些细节决定了结果的可靠性。

一、扩散系数的物理基础

Einstein关系：MSD(t) = ⟨|r(t) – r(0)|²⟩ = 6Dt（三维扩散），其中D是扩散系数。这个关系只在扩散运动进入正常扩散（normal diffusion，MSD ∝ t）阶段后才成立。在短时间内，粒子处于弹道运动区（MSD ∝ t²），还没开始扩散；在中间过渡区，MSD与时间的关系是非线性的。

从MSD计算扩散系数：对MSD-t曲线的线性区域做最小二乘拟合，斜率k = 6D，所以D = k/6。二维扩散用D = k/4，一维用D = k/2。

典型数值范围：液态金属的扩散系数约10⁻⁹ m²/s（如液态Fe在熔点附近约4.5×10⁻⁹ m²/s）；水分子自扩散系数约2.3×10⁻⁹ m²/s；聚合物中的小分子扩散约10⁻¹⁰ m²/s；离子在固体电解质中约10⁻¹²到10⁻¹⁰ m²/s。

二、LAMMPS中compute msd命令

LAMMPS计算MSD使用compute msd命令：

compute 1 all msd com yes

com=yes表示在计算MSD前去除质心运动。如果不去除质心漂移，整体平动会污染MSD，导致扩散系数严重高估。

输出4列：MSD_x、MSD_y、MSD_z、MSD_total。

通过fix ave/time输出：

fix 1 all ave/time 100 1 100 c_1[4] file msd.txt

每100步输出一次总MSD。

关键参数设置：

时间步长：扩散系数计算需要足够长的时间窗口来观察正常扩散区。对于液态金属，通常需要10-50 ns的轨迹才能进入线性区。时间步长1 fs时，需要1000万-5000万步。

原子unwrapping：粒子穿越周期性边界时坐标会”跳变”，导致MSD突然下降。必须使用unwrapped坐标：dump 1 all custom 100 dump.lammpstrj id type xu yu zu。xu/yu/zu是unwrapped坐标，记录粒子的真实连续位置。

参考起点：compute msd以模拟的第一步为参考起点（t=0）。这意味着MSD会随着模拟时间不断增大，长模拟的MSD可能很大。替代方法：用时间原点平均（time-origin average），即对不同的t=0起点分别计算MSD再平均。这需要在后处理中实现，LAMMPS内置的compute msd不做时间原点平均。

三、线性区选择与拟合方法

MSD-t曲线通常有三个区域：

弹道区（0-1 ps）：MSD ∝ t²，粒子在势阱中振荡
过渡区（1-100 ps）：MSD与t的关系非线性
正常扩散区（>100 ps）：MSD ∝ t，可提取扩散系数

线性区的选择是扩散系数计算中最主观的步骤。推荐方法：

绘制log(MSD) vs log(t)图，正常扩散区的斜率应接近1.0
在斜率为9-1.1的区间内做线性拟合
拟合的R²应大于0.99

误差估计：扩散系数的统计误差来自轨迹长度的有限性。经验法则：模拟时间应至少是线性区起始时间的10倍。如果线性区从100 ps开始，模拟至少需要1 ns。更好的做法是分块计算——把轨迹分成5段，每段独立算D，取平均值和标准差。

四、实战复盘：液态Cu的扩散系数

体系：8000个Cu原子，EAM势（Mishin势），NPT 1500 K，时间步长1 fs，总模拟50 ns。

第一步：MSD计算。compute msd com yes，每100步输出一次。50 ns轨迹共50万个数据点。

第二步：log-log图分析。在0-0.5 ps区间斜率≈2.0（弹道区），0.5-20 ps斜率从2.0降到1.0（过渡区），20 ps以后斜率稳定在1.0±0.05（正常扩散区）。

第三步：线性拟合。在50 ps到5000 ps区间拟合MSD = 6Dt。拟合得到斜率k = 3.42×10⁻⁵ Å²/fs = 3.42×10⁻⁹ m²/s。D = k/6 = 5.7×10⁻¹⁰ m²/s。

但实验值（中子散射）给出液态Cu在1500 K的D约4.2×10⁻⁹ m²/s。模拟值偏低约15%。

第四步：误差分析。分5段计算：D值在4.8-6.5×10⁻¹⁰之间，平均值5.6±0.7×10⁻¹⁰ m²/s。统计误差约12%。

偏差来源分析：EAM势拟合的是固态性质，对液态的描述可能有系统偏差。另外，NPT系综下的密度波动也影响扩散——如果密度比实验高2%，扩散系数会偏低约8%。

五、常见问题与解决方法

MSD不线性：如果MSD-t曲线在长时间仍不呈线性，可能是以下原因：(1) 体系处于玻璃态或凝胶态，粒子被束缚，呈现亚扩散行为（MSD ∝ t^α, α<1）；(2) 模拟时间太短，还没进入正常扩散区；(3) 体系太小，周期性边界效应显著——粒子在短时间内就”看到”了自己的镜像。解决方法：增大体系（至少让盒子边长 > 2×扩散距离）或延长模拟时间。

各向异性扩散：在层状材料或受限体系中，不同方向的扩散系数不同。LAMMPS的compute msd输出x/y/z三个方向的MSD，应分别拟合。例如石墨烯层间水的扩散：面内D约2.0×10⁻⁹ m²/s，面外D约0.3×10⁻⁹ m²/s，差了近7倍。

质心漂移：即使设置了com=yes，如果体系有净流动（如剪切流），质心去除可能不够。检查方法：单独计算质心的MSD，如果它不为零，说明有系统漂移。