手机版

LAMMPS计算径向分布函数 — 从命令设置到结构解读的实战指南

发布时间：2026-07-15 来源：科研学术网

小中大

字号：

LAMMPS计算径向分布函数是分子动力学分析中最基础也最常用的操作之一。径向分布函数g(r)描述了以某个原子为中心、距离r处找到另一个原子的概率密度相对于理想气体的比值。LAMMPS计算径向分布函数的结果直接反映体系的局部结构特征——液态、非晶态、晶态的g(r)形态截然不同。

一、径向分布函数的物理意义

g(r)的定义：以A类原子为中心，在距离r到r+dr的球壳内找到B类原子的数目，除以相同体积内理想气体的平均数目。公式为g_AB(r) = V/(N_A × N_B × 4πr²dr) × ⟨Σ_i∈A Σ_j∈B δ(r – |r_ij|)⟩。

g(r)的典型特征：

液态金属：第一峰尖锐（g≈3-4），位于最近邻距离处；第二峰较弱，之后迅速趋于1.0
非晶固体：第一峰与液态类似但更宽，第二峰分裂为两个子峰（这是非晶结构的标志性特征）
晶态：尖锐的峰对应晶格位置，峰间g(r)=0；温度升高后峰展宽

通过g(r)可以提取的关键信息：配位数（第一峰积分）、键长（第一峰位置）、结构有序度（峰的锐度）、中间程有序（第二峰特征）。

二、LAMMPS中compute RDF命令详解

LAMMPS计算径向分布函数使用compute rdf命令。基本语法：

compute 1 all rdf 100

这会计算all-all类型的RDF，分100个bin，截断半径等于力场截断半径。输出4列：r、g(r)、累积配位数coord(r)。

计算特定类型间的RDF：

compute 1 all rdf 100 1 1 1 2 2 2

这会计算三组RDF：1-1、1-2、2-2类型间的。输出列数 = 2 + 2×N_pairs。

关键参数设置：

bin数量：100-200个bin足够。太少g(r)曲线粗糙，太多每个bin内统计样本不足导致噪声增大。bin宽度 = 截断半径/bin数，通常0.02-0.05 Å。

截断半径：默认等于力场的LJ截断半径（通常10-12 Å）。如果需要看长程结构（如中间程有序），需要增大截断半径。方法：在run之前设置pair_style lj/cut 15.0（扩大到15 Å），同时compute rdf的截断也会自动跟随。

输出频率：通过fix ave/time控制。典型设置：每100步采样一次，每1000步输出一次平均值。这样每个输出数据点是10个采样的平均，降低统计噪声。

fix 1 all ave/time 100 10 1000 c_1[*] file rdf.txt mode vector

三、不同体系的RDF特征与分析

液态金属（液态Fe为例）：Fe-Fe RDF在2.48 Å处有第一峰（g≈3.2），对应最近邻配位数约12（fcc/bcc液态的典型值）。第二峰在4.7 Å处，之后g(r)快速趋于1.0。第一峰的宽度反映热运动幅度——在1811 K（熔点）时半高宽约0.3 Å，在2000 K时展宽到0.4 Å。

聚合物（PE熔体为例）：C-C RDF在1.53 Å处有尖锐的第一峰（化学键），3.8 Å和5.2 Å处有较宽的峰（次近邻和第三近邻）。与液态金属不同，聚合物的g(r)在长程不趋于1.0而是趋于0.8左右——这是因为聚合物链的连通性约束使得长程密度略低于平均。

水溶液体系：O-O RDF在2.75 Å处有第一峰（氢键距离），3.4 Å处有第二峰（第二配位层）。第一谷在3.3 Å处，积分到第一谷得到配位数约4.2（实验值4.0，偏差5%）。第一峰的强度g≈4.5反映水中强氢键网络。

四、实战复盘：液态Cu的RDF计算

项目目标：计算液态Cu在1500 K（高于熔点1358 K）的RDF，与中子散射实验对比。

体系：8000个Cu原子，EAM势函数（Mishin势），NPT系综，1500 K，时间步长1 fs。

第一步：体系平衡。先在2000 K下NPT平衡5 ns（让体系完全熔化），然后降温到1500 K再平衡5 ns。确认温度波动<2%，密度稳定在8.02 g/cm³（实验值8.04 g/cm³）。

第二步：RDF计算设置。

compute myRDF all rdf 200 1 1
fix outRDF all ave/time 100 10 1000 c_myRDF[*] file cu_rdf.txt mode vector

截断半径12 Å（跟随EAM势的截断），200个bin，bin宽度0.06 Å。每100步采样，每1000步（1 ps）输出一次平均值。

第三步：产出模拟50 ns，获得50个RDF数据块。最后对所有数据块做时间平均，得到统计误差<2%的g(r)曲线。

第四步：结果分析。Cu-Cu RDF第一峰在2.55 Å处，g_max=2.85。第一谷在3.45 Å处，g_min=0.42。积分第一峰到第一谷得到配位数11.8（实验值12.0，偏差1.7%）。

与中子散射实验对比：第一峰位置偏差0.03 Å（实验2.52 Å），第一峰高度偏差6%（实验g_max=3.04）。考虑到EAM势是拟合固态数据的外推，液态RDF的精度已经相当好。

五、常见问题与解决方法

统计噪声过大：如果g(r)曲线在长程剧烈震荡不收敛到1.0，通常是采样不足。解决方法：延长模拟时间，或增大bin宽度（减少bin数量）。经验法则：每个bin内的原子对数目应>1000才能得到<5%的统计误差。

截断半径不够大：需要看第二、第三配位层时，力场截断半径可能不够。解决方法：在分析阶段单独设置一个更大的截断半径。在LAMMPS中可以用rerun命令重新读取轨迹文件，用更大的截断半径计算RDF，不影响原始模拟。

多组分体系的RDF组合：含N种原子的体系有N×(N+1)/2种RDF对。例如含Cu、Zr、Ti三元素的金属玻璃，需要计算Cu-Cu、Cu-Zr、Cu-Ti、Zr-Zr、Zr-Ti、Ti-Ti共6组RDF。LAMMPS的compute rdf命令支持一次计算多组RDF，但输出文件会有很多列，需要仔细对应。