计算径向分布函数

一、项目背景与科学问题

径向分布函数 g(r) 是描述无序体系（液体、非晶、熔体）微观结构的核心统计量，它给出距离参考粒子 r 处出现另一个粒子的概率密度相对于均匀分布的比值。g(r) 直接与 X 射线/中子散射的结构因子 S(q) 互为傅里叶变换对，因此成为连接原子模拟与实验衍射数据的桥梁。从 g(r) 中可提取最近邻距离、配位壳层结构与配位数等关键结构信息，是验证力场准确性、研究液-固转变与熔体短程有序的标准工具。

本项目以液态铜（Cu）在 1500 K 的分子动力学模拟为对象，计算并分析其径向分布函数与累积配位数，讨论液态金属的结构特征与熔化前后结构变化。

二、计算方法与理论背景

2.1 径向分布函数定义

g(r) 定义为 g(r) = ρ(r)/ρ₀，其中 ρ(r) 为距离参考原子 r 处球壳内的局域数密度，ρ₀ = N/V 为平均数密度。计算时对轨迹中每一对原子统计距离直方图并归一化：g(r) 在无穷远处趋于 1（均匀分布）。对液态金属，g(r) 呈现”短程有序、长程无序”的特征：第一峰尖锐，随后振荡快速衰减至 1。

2.2 配位数

由 g(r) 可计算累积配位数 N(r) = 4πρ₀∫₀ʳ g(r′)r′²dr′，其第一壳层配位数（积分至第一谷）对液体给出平均最近邻原子数。对液态 Cu，第一壳层配位数约 11–12，略低于面心立方晶体的 12，反映液体中配位壳层的松散化。配位数是理解液体堆垛密度与输运性质的重要参数。

2.3 与散射实验的关系

g(r) 与结构因子 S(q) 满足 S(q) = 1 + 4πρ₀/q ∫₀^∞ [g(r)−1] r sin(qr) dr。因此模拟计算的 g(r) 可与中子/ X 射线散射实验测量的 S(q) 互验，是模拟-实验对照的核心桥梁。

三、计算设置与模型构建

• 软件：LAMMPS；势函数：EAM 嵌入原子势（描述 Cu 金属键合）。
• 体系：4096 个 Cu 原子，立方盒，周期性边界。
• 系综：NPT，温度 1500 K（高于熔点约 1358 K，液态）、压力 1 bar。
• 模拟：2 ns 平衡 + 2 ns 采样，步长 1 fs。
• 分析：对采样轨迹统计原子对距离直方图，归一化得 g(r)，并数值积分得 N(r)。

四、结果分析与案例图解读

图 4 的左面板给出液态 Cu 在 1500 K 的径向分布函数 g(r)（蓝色曲线）。g(r) 在短程区域呈现尖锐的第一峰，峰值位置 2.58 Å（红色虚线标注），峰高约 2.7，对应液态 Cu 的最近邻原子间距；第一谷位于约 3.3 Å，定义第一配位壳层边界；第二峰位于约 4.85 Å，强度明显减弱；随后 g(r) 围绕 1 振荡且振幅快速衰减，在约 6 Å 后趋近于 1。这一”第一峰尖锐、长程衰减”的形态是液态金属的典型特征，与 X 射线衍射实验得到的液态 Cu 结构因子转换结果一致。第一峰位 2.58 Å 略大于固态 Cu 的最近邻距离（约 2.55 Å），反映熔化后原子间距略微膨胀。

图 4 的右面板给出累积配位数 N(r) 曲线（橙色曲线）。N(r) 随 r 单调上升，在跨越第一配位壳层（r ≈ 3.3 Å）后进入平台区，第一壳层配位数约 11.5（红色虚线标注），与液态 Cu 的经典值一致。配位数接近 12 说明液态 Cu 的局域堆垛仍保留面心立方晶体的特征，但平台较固态更平缓，体现液体的结构松散化。N(r) 曲线的平台高度为后续计算扩散系数、粘度等输运性质的模型验证提供了结构基准。

两面板结合，左面板给出原子间距与壳层结构（g(r)），右面板给出定量配位数（N(r)），共同完成液态 Cu 微观结构的系统表征，为验证 EAM 势函数、对比实验散射数据与理解液态金属的短程有序提供依据。

五、结论与展望

本项目通过 LAMMPS 分子动力学计算了液态 Cu 在 1500 K 的径向分布函数，第一峰位于 2.58 Å、峰高约 2.7，第一壳层配位数约 11.5，结果与实验衍射数据一致。案例图以 g(r) 曲线与 N(r) 配位数曲线双面板呈现。

后续可拓展：(1) 计算液态 Cu 的自扩散系数与粘度，建立结构与输运性质关联；(2) 降温跨越熔点，原位观察 g(r) 峰分裂与晶化过程；(3) 对合金熔体（Cu-Ag、Cu-Ni）计算偏径向分布函数，研究化学短程有序。