玄铸计算平台案例7：FCC Al 的径向分布函数计算-CSDN博客

本文链接：https://blog.csdn.net/Piaen/article/details/148567386

径向分布函数的定义及物理意义

定义与公式

径向分布函数（Radial Distribution Function, g®描述了在一个粒子周围距离为 r 处找到另一个粒子的相对概率，相对于理想气体中粒子均匀分布的情况。
常见的表达式为：
$\frac{1}{4\pi r^2 \rho} \frac{dN(r)}{dr}$
其中

$d N (r)$ 是粒子在距离 $r$ 到 $r + d r$ 之间的平均数目，
$\rho$ 是体系的平均粒子密度，
$4\pi r^2 dr$ 是以粒子为中心、半径为 $r$ 的球壳体积。

物理意义

当 $g (r) = 1$ ：粒子分布如同理想气体，没有相关性。
当 $g (r) > 1$ ：在距离 $r$ 处找到粒子的概率高于随机分布，通常对应有结构或配位关系（如液体、固体的局部有序）。
当 $g (r) < 1$ ：在距离 $r$ 处找到粒子的概率低于随机分布，表示排斥效应或空洞区域。
$g (r)$ 的第一个峰通常对应最近邻粒子的平均位置。

在电池领域，RDF可分析锂离子溶剂化壳层；在催化反应和生物医药领域，RDF则帮助阐明分子间相互作用机制。在材料科学中，它用于表征合金原子排列

文献计算结果

J. Mater. Environ. Sci 9.1 (2018): 93-99
https://www.researchgate.net/profile/Hanae-Chabba/publication/322195076_Modeling_aluminum_using_molecular_dynamics_simulation/links/5e64b8ab299bf1744f67ca1f/Modeling-aluminum-using-molecular-dynamics-simulation.pdf
关注该文献中通过lammps计算,采用EAM势，得到的FCC (面心立方) Al的晶胞的径向分布函数(红色)

在这里， $g (r)$ 体现了原子层层排列的规则性，它揭示了晶体中每个原子周围的空间有多大概率分布着其他原子，是描述晶体有序性、热振动、结构缺陷最基础的函数。

每个清晰的峰代表一个壳层，对应晶体中某类近邻原子的平均距离。
峰的位置 = 特定晶格距离
峰的高度 = 排列有序程度（越高越有序）

结果表明，由EAM势得到的Al_fcc晶格最近邻原子间距离为3 Å

计算软件和参数设置	计算软件：LAMMPS
计算软件和参数设置	势函数	温度
Al (fcc）	EAM 势	300K

传统的建模及计算过程

以传统方法LAMMPS计算Al_fcc的RDF为例,流程包括建模、势函数下载、输入文件编写、数据处理及作图等步骤

建模过程

先从Materials Project数据库中读取Al-bcc晶格的晶格常数为：4.04Å
编辑Al_fcc.in文件来建模超胞大小为101010的Al-fcc晶胞
运行后，文件夹下会有一个命名为‘Al_fcc.data’的结构文件

势函数的获得

需要去网站上获取计算所需的势函数，另存至本地，再上传到集群用于lammps计算
Al_Zhou04.eam.alloy

计算过程

首先，需要对结构进行弛豫，使原子系统在目标温度（300 K）下达到稳定振动状态；随后才能进行 RDF 计算。
因此，需要分别编写用于弛豫过程和 RDF 计算的 LAMMPS 输入文件 rdf_300K.in

计算结果

计算结束后，需要先对结果进行数据处理，再下载到本地，通过代码或一些作图软件绘制RDF图

数据处理

计算结束后，rdf的数据输出结果在文件rdf_300K.txt文件中，此时文件有4列，下面给出示例

bin_index     r     g_ij1     g_ij2

这通常是内部结构安排，为了多原子对预留空间,因为作图只需要两列，所以需要提取中间两列到rdf_1.txt文件

之前设置的每 100 步取一次数据,每次写入时会在每组数据后输出一个“时间步数标签”行,所以rdf_1.txt文件每隔100步还会出现 100 ，因此需要去掉这一行。

具体数据处理指令如下：

awk '{print $2, $3}' rdf_300K.txt > rdf_1.txt
awk 'NF==2' rdf_1.txt > rdf_clean.txt

第一步数据处理后的rdf_1.txt文件

径向分布函数RDF作图

首先需要将rdf_clean.txt文件下载到本地
编写代码(或用其他作图软件)作图
从传统方法Lammps计算得到Al_fcc晶胞的径向分布函数来看，第一个尖峰出现在2.84Å 左右，与文献值（3.00Å）接近，差距可能来自不同计算软件，参数设置，以及晶胞大小不同，在合理范围内。

使用 MatModeler 进行建模计算

只需要建模，提交任务(这一步设置参数)，查看结果(直接输出了RDF图)三步即可

建模过程

只需要直接从内置数据库中导入Al晶胞（FCC对应的空间群为Fm-3m),再扩胞即可.具体操作如下：

导入：新建文件 - 从数据库中导出 - 从晶体材料数据库导入-仅包含输入元素-选择Al
扩胞：直接选’扩胞’，设置 $10 * 10 * 10$ 即可

计算过程

分子动力学模拟选择M-Mat模块
不需要单独上传力场文件，MatModeler已支持gaff小分子力场，M3GNET和多种eam力场
不需要单独编写输入文件，直接在任务提交时，设置参数即可

建模完成后点击右侧边栏的请选择APP，选择M-mat，进行任务参数配置，选择从可视化窗口导入，选择当前文件（ $A l - 10 * 10 * 10. c i f$ ）,点击确认
直接在性质计算栏选择‘径向分布函数’(可以直接设置参考中心和周围粒子选择，如果不设置，则会计算以结构文件中第一个分子的分子类型为中心，所有分子类型对其的RDF)
力场选择eam,直接在eam力场中选择Al_zhou.eam.alloy势函数用作计算
预设模拟参数选择‘custom’(Matmodeler还支持两种预设的流程beginner和product)，在这里直接设置模拟的温度为300K，时间步长dt为0.002ps，模拟步数nsteps为10000步
最后提交任务