1. Lennard-Jones流体平衡态模拟

案例目标与物理模型

本案例使用 LAMMPS 对三维 Lennard-Jones 流体进行 NVT 系综模拟。LJ 势能形式为：

U(r) = 4ε [ (σ/r)^12 - (σ/r)^6 ]

模拟采用约化单位（ε=1, σ=1, m=1），截断半径取 r_c = 2.5σ，并在截断处施加长程能量与压力修正。体系为 500 个粒子，密度 ρ* = 0.80，目标温度 T* = 1.20，时间步长 Δt* = 0.002。初始构型为简单立方晶格，速度按目标温度随机赋值并去除质心动量。

LAMMPS 输入脚本

# ---------- 初始化 ----------
units           lj
atom_style      atomic
boundary        p p p

# ---------- 创建体系 ----------
lattice         sc 1.0772          # ρ* = 0.80 对应晶格常数
region          box block 0 5 0 5 0 5
create_box      1 box
create_atoms    1 box
mass            1 1.0

# ---------- 势函数 ----------
pair_style      lj/cut 2.5
pair_coeff      1 1 1.0 1.0 2.5
neighbor        0.3 bin
neigh_modify    every 1 delay 0 check yes

# ---------- 初始速度 ----------
velocity        all create 1.20 4928459 mom yes rot yes dist gaussian

# ---------- 平衡阶段：NVT ----------
fix             1 all nvt temp 1.20 1.20 1.0
timestep        0.002
thermo          100
thermo_style    custom step temp press pe ke etotal

# 先跑 5000 步让体系初步松弛
run             5000

# ---------- 采样阶段：继续 NVT，记录涨落 ----------
reset_timestep  0
fix             2 all ave/time 10 100 1000 temp press file thermo_avg.dat
run             20000

# ---------- 输出最终构型 ----------
write_data      lj_nvt_final.data

温度与压力涨落观察

NVT 系综下，温度由 Nosé-Hoover 热浴维持在 1.20 附近，但瞬时温度仍会围绕目标值涨落。压力在 LJ 流体中涨落幅度更大，因为压力包含位力项，对粒子局部构型非常敏感。判断平衡的核心依据是：

• 温度时间序列不再有单调漂移，围绕 1.20 做有界涨落；
• 压力时间序列的滑动平均趋于稳定，不再持续上升或下降；
• 势能与总能量在采样阶段无系统性趋势；
• 径向分布函数在后续分析中不再随时间明显变化（本案例不展开）。

平衡判据与结果判读

将 thermo_avg.dat 中的温度、压力对时间步作图，观察最后 10000 步的滑动平均。以下为典型平衡判读逻辑：

结果验证与常见问题

运行完成后，检查 thermo_avg.dat 中最后 10000 步的统计量。若温度均值在 1.20±0.02 范围内、压力滑动平均变化小于 0.1，可认为体系达到热力学平衡。常见问题包括：

• 温度失控：检查 fix nvt 的阻尼参数是否合理，LJ 单位下通常取 1.0 左右。
• 压力异常大：确认截断半径与长程修正是否一致，pair_style lj/cut 2.5 已包含能量修正，但压力修正需在 thermo_style 中确认。
• 体系未平衡：延长平衡阶段至 10000 步以上，或检查初始构型是否过于规则导致局部应力过大。