3. LJ流体均方位移与扩散系数

任务定位与物理量定义

本例在NVE系综下运行三维LJ流体，从平衡轨迹中计算均方位移（MSD），并用Einstein关系拟合扩散系数 D。核心公式为：

MSD(t) = <|r_i(t) - r_i(0)|²>
Einstein关系（三维）：D = (1/6) · lim_{t→∞} d[MSD(t)]/dt

关键约束：NVE下总能量守恒，温度由初始速度标定决定；MSD必须使用未折叠坐标（unwrapped），否则粒子跨越周期边界时位移会被错误截断。扩散系数拟合应取MSD进入线性区间的时段，避开弹道区（t ≲ 0.5）和统计噪声主导的长时尾部。

LAMMPS输入脚本与关键设置

以下脚本在NVE系综下运行LJ流体，输出未折叠坐标用于MSD计算。单位采用LJ约化单位（σ=ε=m=1）。

# LJ流体NVE模拟 + MSD输出
units           lj
atom_style      atomic
boundary        p p p
region          box block 0 10 0 10 0 10
create_box      1 box
create_atoms    1 random 500 12345 box
mass            1 1.0

pair_style      lj/cut 2.5
pair_coeff      1 1 1.0 1.0 2.5

velocity        all create 1.0 4928459 dist gaussian
fix             1 all nve
timestep        0.002

thermo          500
run             2000

# 重置时间原点，输出未折叠坐标
reset_timestep  0
dump            d1 all custom 100 traj_unwrap.lammpstrj id xu yu zu
dump_modify     d1 sort id
run             50000

MSD计算与Einstein拟合

使用Python对轨迹进行后处理，计算多时间原点平均的MSD，并在线性区间拟合斜率。

import numpy as np

def read_lammps_traj(filename):
    """读取LAMMPS自定义dump，返回坐标数组 (n_frames, n_atoms, 3)"""
    frames = []
    with open(filename) as f:
        while True:
            line = f.readline()
            if not line:
                break
            if 'ITEM: TIMESTEP' in line:
                f.readline()  # timestep
                f.readline()  # ITEM: NUMBER OF ATOMS
                n_atoms = int(f.readline())
                f.readline()  # ITEM: BOX BOUNDS
                for _ in range(3):
                    f.readline()
                f.readline()  # ITEM: ATOMS id xu yu zu
                coords = np.array([f.readline().split()[1:4] for _ in range(n_atoms)], dtype=float)
                frames.append(coords)
    return np.array(frames)

traj = read_lammps_traj('traj_unwrap.lammpstrj')
n_frames, n_atoms, _ = traj.shape
dt = 0.002 * 100  # 每100步输出一次，时间间隔

# 多时间原点平均MSD
msd = np.zeros(n_frames)
count = np.zeros(n_frames)
for lag in range(n_frames):
    disp = traj[lag:] - traj[:n_frames-lag]
    sq_disp = np.sum(disp**2, axis=2)  # (n_frames-lag, n_atoms)
    msd[lag] = np.mean(sq_disp)
    count[lag] = n_frames - lag

time = np.arange(n_frames) * dt

# 选择线性拟合窗口：t* ∈ [1, 5]
fit_mask = (time >= 1.0) & (time <= 5.0)
coeffs = np.polyfit(time[fit_mask], msd[fit_mask], 1)
slope = coeffs[0]
D = slope / 6.0

print(f"拟合斜率 = {slope:.4f}")
print(f"扩散系数 D = {D:.4f} (LJ单位)")

温度对扩散行为的影响

在三个温度下重复上述流程，观察D随T*的变化。以下为教学示例的预期趋势（非实测数据）：

结果验证与常见问题

• MSD曲线形状检查：短时应有弹道区（MSD∝t²），中时应进入线性区（MSD∝t），长时应因统计不足而发散。若全程线性，可能时间原点平均不足或输出间隔过大。
• 扩散系数合理性：LJ液体在T*=1.0、ρ*=0.8附近，D约在0.03–0.08（LJ单位）量级。若D接近0，检查是否误用了折叠坐标。
• 温度一致性：用 compute temp 输出平均温度，确认与目标T*一致。NVE下温度漂移过大说明步长或截断设置有问题。
• 拟合窗口选择：建议用MSD对t的log-log图确认线性区起点，避免将弹道区纳入拟合。