2. LJ流体径向分布函数计算
从平衡轨迹到 g(r)：后处理流程与关键参数
本例对已平衡的LJ流体轨迹做后处理，计算径向分布函数 g(r)。核心思路是：以每个粒子为参考中心，统计距离 r 到 r+Δr 球壳内的粒子数，再对全部帧、全部参考粒子做平均，最后用理想气体均匀密度归一化。LJ约化单位下，温度 T*=1.0、密度 ρ*=0.80 的液体，g(r) 第一峰通常出现在 r*≈1.08–1.12，第一配位壳层配位数约 12。
计算式：g(r) = (1/(ρ·N·V_shell)) · Σ_frames Σ_i n_i(r, r+Δr)
其中 V_shell = (4/3)π[(r+Δr)³ − r³]，ρ = N/V 为体系数密度。
后处理需要三个输入：平衡后的轨迹文件（含坐标）、拓扑/原子数信息、以及积分参数（r_max、Δr）。LAMMPS 的 compute rdf 可直接完成统计，也可用 Python 读取 dump 自行实现。下面给出可直接运行的 LAMMPS 后处理脚本，复用平衡阶段的数据文件。
# in.rdf —— 对平衡后轨迹计算 g(r)
units lj
atom_style atomic
boundary p p p
read_data lj_fluid.data # 平衡后的构型
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
# 读取平衡轨迹（NVT 平衡段已保存的 dump）
# 若已有 data 文件，可直接从该构型继续采样
velocity all create 1.0 4928459 dist gaussian
fix 1 all nve
fix 2 all langevin 1.0 1.0 1.0 48279
thermo 100
timestep 0.005
# 关键：定义 RDF 计算，50 个 bin，截断 2.5σ
compute myrdf all rdf 50 1 1 cutoff 2.5
fix rdfout all ave/time 100 1 1000 c_myrdf[*] file rdf.dat mode vector
run 20000 # 采样 100 帧
运行后 rdf.dat 每行包含 r、g(r)、配位数。判断结果是否可信，先看 g(r) 在 r 小于约 0.9σ 时是否严格为 0（排除原子重叠），再看大 r 处是否收敛到 1.0。若大 r 处偏离 1 明显，说明采样帧数不足或体系未真正平衡。
从图中可读出三个判读要点：第一，r 小于约 0.9σ 区间 g(r) 必须为 0，否则说明初始构型存在原子重叠或积分步长过大；第二，第一峰位置 r*≈1.1 对应最近邻距离，峰高反映短程有序程度，温度升高时峰高下降、峰宽增大；第三，第一峰与第一极小之间的积分面积即第一配位壳层配位数，LJ液体在 ρ*=0.8、T*=1.0 附近约为 12，与面心立方密堆的 12 个最近邻一致。
若 g(r) 在大 r 处不收敛到 1，优先检查：采样帧数是否足够（建议 ≥100 帧）、r_max 是否超过盒子半边长（应 ≤ L/2）、体系是否真正平衡。r_max 超过 L/2 会因周期性边界导致统计失真。
结果验证与常见问题
| 观察项 | 正常表现 | 异常与原因 |
|---|---|---|
| r < 0.9σ 的 g(r) | 严格为 0 | 非零：原子重叠，检查初始构型或时间步长 |
| 第一峰位置 | r*≈1.08–1.12 | 偏移过大：密度或温度设置不符 |
| 大 r 收敛 | g(r)→1.0 | 偏离：采样不足或体系未平衡 |
| 第一配位壳层配位数 | ≈12（ρ*=0.8, T*=1.0） | 偏差大：积分区间选取不当 |
若需自行用 Python 实现，读取 dump 后对每帧、每个粒子统计球壳内邻居数，再除以 ρ·V_shell 并做帧平均即可。关键是把周期性边界的最小镜像距离算对，否则 g(r) 在小 r 处会出现虚假峰。