LAMMPS计算RDF：从轨迹到结构信息的完整分析链条

LAMMPS计算RDF是项目组用得最频繁的后处理分析——几乎每次MD模拟完成后第一件事就是画RDF。但用对了和用错了天差地别。项目组在液态Cu的RDF分析中，最初用compute rdf命令直接输出，结果第一峰高度明显偏低——与文献值有显著差距。排查后发现是采样时间不够（只用了10 ps轨迹），增加到200 ps后第一峰高度恢复正常，与文献吻合。RDF不只是”跑个命令”，采样策略和后处理直接影响结果可信度。

RDF的物理含义

径向分布函数g(r)描述的是：在距离参考原子r处找到另一个原子的概率密度，归一化到理想气体的均匀分布。

g(r) = V/N × <Σ δ(r - |r_i - r_j|)> / (4πr²Δr)

物理意义：

• g(r)=1 → 该距离处的分布与均匀气体相同（无结构）
• g(r)>1 → 该距离处原子富集（有序排列）
• g(r)<1 → 该距离处原子贫乏（排斥区）

对于液体，RDF通常显示几个振荡峰——第一峰对应最近邻配位壳层，第二峰对应次近邻。随r增大，振荡逐渐衰减到g(r)=1（长程无序）。

LAMMPS中的RDF计算

LAMMPS使用compute rdf命令计算RDF：

compute 1 all rdf 100

这会计算全体系的RDF，分100个bin，截断半径等于力场截断半径。但这个默认设置有几个问题：

问题1：截断半径不够大

默认截断等于力场截断（如12 Å），但RDF的结构信息可能延伸到15-20 Å。需要显式指定截断：

compute 1 all rdf 100 cutoff 15.0

问题2：bin数量影响分辨率

100个bin在15 Å截断下，每个bin=0.15 Å。对于液态金属的第一峰（半高宽约0.3 Å），0.15 Å的分辨率刚好够。如果峰更窄（如晶体），需要增加bin数：

compute 1 all rdf 200 cutoff 15.0

200个bin在15 Å下，每个bin=0.075 Å，分辨率足够分辨尖锐峰。

问题3：多组分体系的偏RDF

对于二元体系（如Cu-Ag合金），需要计算偏RDF：g_CuCu(r)、g_CuAg(r)、g_AgAg(r)。LAMMPS支持多类型RDF：

compute 1 all rdf 100 1 1 1 2 2 2

这里1 1计算type 1-type 1的RDF，1 2计算type 1-type 2，2 2计算type 2-type 2。输出文件会有6列（每对类型：r、g(r)、coordination(r)）。

从RDF提取结构信息

配位数：

第一配位数（最近邻原子数）通过对RDF第一峰积分获得：

N_1 = 4π ρ ∫₀^{r_min} r² g(r) dr

其中r_min是第一峰和第二峰之间的极小值位置，ρ是数密度。

LAMMPS的compute rdf输出第三列就是累积配位数，直接读取r_min处的值即可。液态Cu的第一配位数约11.5（固态FCC是12），反映了液态中最近邻数的略微减少。

结构因子：

RDF通过傅里叶变换与X射线/中子散射的结构因子S(Q)直接联系：

S(Q) = 1 + 4π ρ ∫₀^∞ [g(r) - 1] × sin(Qr)/(Qr) × r² dr

项目组用Python做这个变换：

import numpy as np

# 读取LAMMPS输出的RDF数据
data = np.loadtxt("rdf.txt", skiprows=4)
r = data[:, 0]
g = data[:, 1]

rho = 0.0755  # 液态Cu数密度 (atoms/ų)

# 结构因子计算
Q = np.linspace(0.5, 15, 200)
S = np.zeros_like(Q)
for i, q in enumerate(Q):
    integrand = (g - 1) * np.sin(q * r) / (q * r) * r**2
    S[i] = 1 + 4 * np.pi * rho * np.trapz(integrand, r)

# 输出
np.savetxt("structure_factor.txt", np.column_stack([Q, S]))

计算出的S(Q)第一峰位置在Q≈3.0 Å⁻¹附近——与液态Cu的X射线散射实验一致，验证了模拟结构的正确性。

有限尺寸校正

小体系（<1000原子）的RDF会受到周期性边界条件的影响——在r>L/2（L是盒子边长）处，RDF会因为镜像原子而失真。项目组的处理方法：

1. 只分析r<L/2的部分：最简单的做法，但损失了长程信息
2. 外推校正：假设g(r> L/2)=1，用解析方法校正有限尺寸效应。对于液体，这个假设是合理的（长程无序）
3. 增大体系：用更大的模拟盒子，确保L/2>15 Å（覆盖RDF的所有有意义结构）

项目组在液态Cu的模拟中用4096个原子（盒子约35 Å），L/2=17.5 Å，足够覆盖RDF的所有振荡结构。

采样策略：时间窗口的影响

RDF的统计精度取决于采样量。项目组测试了不同采样时长对液态Cu RDF第一峰高度的影响：

10 ps的采样高度偏低，因为短时间采样无法覆盖构型空间的充分采样。200 ps后基本收敛。标准差随√(1/t)减小——理论上要减半标准差需要4倍采样时间。

项目组的标准做法：先跑100 ps平衡，再做200 ps以上的生产模拟，每100步输出一次构型用于RDF计算。这样有2000个构型做平均，RDF的统计误差<2%。

多体系RDF对比分析

项目组在LiPF6/EC-DMC电解液体系中，需要同时分析多个偏RDF来理解溶剂化结构：

从偏RDF可以判断溶剂化壳层的组成：如果Li-O(EC)的第一峰很高而Li-F(PF6⁻)很低，说明Li⁺主要被溶剂分子包围（溶剂分离离子对, SSIP）；反之说明PF6⁻进入了第一溶剂化壳层（接触离子对, CIP）。这个信息对理解电解液的离子电导率和解离行为至关重要。

反思：RDF是必要但不充分的结构描述

RDF是一维投影——它把三维的结构信息压缩到了径向距离上。两个不同的三维结构可能给出相同的RDF。例如，BCC和HCP的RDF在某些距离上非常相似。要全面描述结构，还需要结合角分布函数（ADF）、键序参数（Steinhardt参数）和Voronoi分析。

RDF最大的价值在于与实验散射数据的直接对比——它是模拟与实验之间的桥梁。更多MD后处理分析的实战经验，可以参考分子动力学栏目，或返回科研学术网首页。