LAMMPS计算径向分布函数：参数设置与物理含义的深度剖析

LAMMPS计算径向分布函数是分子动力学后处理中最基础也最常用的分析手段——无论是判断体系是否达到平衡态、识别晶体与液体的结构差异，还是表征合金的局部原子环境，RDF都是第一步要看的指标。但做过LAMMPS计算径向分布函数的人都知道，一条看起来”还行”的RDF曲线背后，参数选择的每一步都可能把物理信号掩盖在人为噪声里。

项目背景：金属玻璃Cu-Zr的结构表征

团队去年处理了一个Cu₆₄Zr₃₆金属玻璃项目，核心任务是用分子动力学验证急冷制备条件下非晶结构的稳定性，并与实验XRD数据做定量对比。LAMMPS计算径向分布函数在这个项目中的角色不只是”画一条曲线看看”——需要从RDF中提取配位数、第一壳层距离和结构因子，再与实验衍射数据交叉验证。

模拟参数：LAMMPS 29Oct2020版本，EAM势函数（Mendelev型Cu-Zr势），NPT系综，时间步1 fs，体系包含8000个原子。急冷方案从2000 K以10¹³ K/s的速率冷却至300 K，然后做100 ns的等温弛豫确保结构稳定。

Bin宽度：RDF形态的第一道关卡

LAMMPS计算径向分布函数时，`compute rdf`命令的第一个关键参数是bin数量（`nbins`），它决定了RDF曲线的空间分辨率。bin宽度Δr = rmax/nbins，其中rmax是截断距离。对于金属体系，第一近邻距离通常在2.5-3.0 Å附近，RDF第一峰的半高宽约0.15-0.20 Å——这意味着bin宽度必须小于0.05 Å才能充分分辨峰形。

团队最初用nbins=200配合rmax=10 Å，bin宽度为0.05 Å——这个参数对第一峰足够，但第二壳层（约4.5 Å）和第三壳层（约6.5 Å）的峰在0.05 Å分辨率下开始出现锯齿状噪声，原因是这些壳层的原子数密度远低于第一壳层，统计样本量不足。LAMMPS计算径向分布函数的统计精度在每个bin上的采样次数与该bin对应的体积成正比（体积随r²增长），远处bin的统计质量天然低于近处——这是RDF的物理特性，不是参数错误。

团队最终选择了折衷方案：用两套参数分别计算短程和长程RDF。短程（rmax=5 Å，nbins=200，Δr=0.025 Å）用于精确分辨第一、第二壳层峰形；长程（rmax=15 Å，nbins=300，Δr=0.05 Å）用于观察中程有序信号。两套RDF在5 Å附近拼接，短程部分在3 Å以内具有可靠统计精度。

截断距离：长程关联的边界

LAMMPS计算径向分布函数的截断距离rmax需要大于体系特征结构的最大关联长度。对于晶体体系，RDF在若干壳层后趋于恒定值1.0，截断距离取6-8 Å即可。但金属玻璃的中程有序特征可能延伸到10 Å以上——如果rmax过短，RDF在截断处会出现人为跳变，而长程振荡信号被完全丢失。

Cu₆₄Zr₃₆的RDF在10 Å附近仍有明显的第二宽峰残余，对应于中程有序结构。团队把rmax设为15 Å后，RDF在12-15 Å区间归一化至1.0附近，但在13 Å处出现了一个异常小峰——这个峰不是物理信号，而是截断边界处的统计噪声放大。LAMMPS计算径向分布函数在远距离的统计质量取决于模拟时长和原子数，8000原子+100 ns的模拟在r>12 Å处的采样统计量已经接近噪声边界。

团队的判断：金属玻璃的RDF分析中，rmax=10 Å是物理信息的可靠上限，超过10 Å的部分需要用结构因子S(q)间接获取——S(q)是RDF的傅里叶变换，在q空间的统计精度远优于r空间的远距离区域。

采样间隔：时间平均的统计平衡

LAMMPS计算径向分布函数的另一个关键参数是时间采样间隔——`fix ave/time`或`fix ave/correlate`中的`Nevery`参数决定了多少步取一次RDF样本。采样间隔过密（如每10步一次）会把相邻帧的原子构型高度关联的样本送入平均，降低统计独立性；采样间隔过疏（如每10000步一次）又可能让总样本数不足以支撑远距离bin的统计精度。

对于Cu₆₄Zr₃₆的100 ns弛豫过程，团队计算了速度自相关函数的衰减时间约0.5 ps（500步），据此设定Nevery=500步——确保相邻采样帧之间原子位置已充分 decorrelated。Nrepeat=20（连续20次采样做一次局部平均），Nfreq=10000步（每10000步输出一次时间平均RDF）。100 ns总模拟时长给出2000个独立RDF样本，统计精度在第一壳层处误差小于0.5%，在10 Å处误差约3-5%——可接受范围。

多组RDF：拆解混合体系的结构信号

Cu₆₄Zr₃₆是二元合金，总RDF是Cu-Cu、Cu-Zr、Zr-Zr三种部分RDF的加权叠加。LAMMPS计算径向分布函数支持分组计算（`compute rdf myRDF all type1 type2`），但这个分组功能在多组分体系中的使用方式需要仔细设计。

团队分别计算了Cu-Cu、Cu-Zr、Zr-Zr三种部分RDF。结果揭示了一个总RDF无法呈现的关键信息：总RDF的第一峰位置在2.75 Å，看起来像是单一特征距离；但部分RDF显示Cu-Zr第一峰在2.68 Å、Cu-Cu第一峰在2.55 Å、Zr-Zr第一峰在3.15 Å——总RDF的第一峰实际上是Cu-Cu和Cu-Zr两个峰的叠加，Zr-Zr的贡献在3.15 Å处独立成峰。如果只看总RDF做结构解读，Cu-Zr最近邻距离会被误判为2.75 Å而非真实的2.68 Å，偏差0.07 Å在合金结构分析中已不可忽略。

配位数分析同样需要部分RDF支撑。从总RDF第一壳层积分得到的配位数是11.2——这个数字无法区分Cu周围有多少Zr、Zr周围有多少Cu。部分RDF积分结果：每个Cu原子周围有5.3个Zr和5.0个Cu，每个Zr原子周围有9.4个Cu和1.2个Zr。Cu原子周围的化学短程有序参数远偏离随机混合值（6.4），证明Cu-Zr非晶中存在强烈的化学有序倾向——这个结论只能从部分RDF推导，总RDF无能为力。

从RDF到结构因子：与实验的定量对接

RDF在r空间的物理直觉清晰，但与XRD实验数据的对比需要经过傅里叶变换进入q空间。LAMMPS计算径向分布函数得到的g(r)通过S(q)=1+4πρ∫[g(r)-1]r²sin(qr)/(qr)dr转换为结构因子，其中ρ是平均数密度。团队用Python脚本实现了这个转换，截断积分上限取10 Å（RDF可靠范围），与实验XRD数据的S(q)在主峰位置（q≈2.8 Å⁻¹）和峰宽上吻合良好，峰值偏差约5%——这个偏差主要来源于EAM势函数对Cu-Zr混合势的描述精度，而非RDF参数设置问题。

从bin宽度到截断距离，从采样间隔到多组RDF拆解——LAMMPS计算径向分布函数的每个参数选择都需要基于体系的物理特征做针对性调整，而不是沿用”默认参数”。参数背后的判断力，才是结构分析从”画曲线”升级到”读信号”的关键。

──────────────────────────────────────────────────

更多内容请访问 https://www.keyanxueshu.com/