在完成分子动力学模拟之后，需要对模拟轨迹进行分析，计算径向分布函数（Radial Distribution Function, RDF）即是一种常见的分析。大多分子动力学模拟软件提供了计算RDF的功能，但往往不够灵活，也不能适应不同的轨迹格式。

MDtraj是一个用于处理分子动力学轨迹的Python包，它能够接受多种轨迹格式，将不同格式的轨迹作为统一的轨迹对象处理；同时，它也提供了许多性质计算的功能，其中就包括RDF；最后，通过使用Python代码，使得性质计算灵活了许多。

本文将先介绍RDF的计算原理，再以一段分子动力学轨迹为例，展示如何通过MDtraj计算RDF，并通过RDF计算配位数。已对RDF熟悉的读者，可跳过第一部分。

径向分布函数的计算原理

均匀分布

首先在1*1 的平面内均匀的生成若干个点，并选择(0.5, 0.5) 作为中心。现在我们想要知道，距离中心点为r 处点的数密度（后文简称为密度）。由于离散情形，实际上我们想知道的是r 到r + $\Delta$ r 这样的圆环内点的密度。可以想象，由于点的分布是均匀的，因此各处密度应该相等，并且与平均密度一致。由于面积为1，因此平均密度在数值上与点的数量相等。

首先将各个圆环画出来

为了计算圆环内点的密度，我们只需统计各圆环内点的数量，然后将其除以圆环的面积。注意到，一个点落到哪个圆环内，只与其距中心的距离R 有关，而与方向无关，这也正是径向的含义。因此，只需计算每个点距中心的距离，并且当r $\le$ R $<$ r+ $\Delta$ r 时，对应r 的圆环内点的计数加1。numpy已经提供了类似的功能，因此这里直接使用。

可以看见，在r < 0.5 的时候，密度接近于1000，也就是点的平均密度。而当r > 0.5 时，由于生成的点只在(0,1)的范围内，此时圆环越来越“空”，密度也快速下降。

但是，这和RDF有什么关系呢？事实上，RDF只不过是用密度再除以一下平均密度而已，$$g(r)=\frac{\rho (r)}{\rho }$$

这里提供一段简短的代码，读者可自行修改点的数量与圆环的数量，观察一下它们对RDF曲线的影响。

引入额外的点

均匀分布的点，其RDF平平无奇，因为任何位置的密度都与平均密度相同。让我们增加一些额外的点，让局部的密度上升，看看会发生什么。

可以看见，RDF曲线上出现了两处峰，这意味着在对应的r 处，点的密度突然增大了。让我们把散点画出来，看看发生了什么。

可以看见，在中心点的左下和右上方，分别有两块点密度很大的区域，RDF上的两个峰，正是对应了这两个区域的出现。这就是说，RDF可以告诉我们沿着径向密度的变化，或说包含了这部分信息，但请注意，RDF不包含方向信息，它只能告诉我们各个方向密度的平均值。

MDtraj 计算径向分布函数

基本上，计算分子动力学轨迹的RDF与前述计算RDF的方式是十分类似的，唯一不同的是，中心点通常不止一个，轨迹不止一帧。放轻松，这并不会带来更多本质上的复杂度，只不过前面的圆环在这里变成了球壳，剩下的只是按照前述的思路分别计算一帧一个中心点的情况，接着求平均值即可。

下面直接使用Mdtraj读取一段水分子的轨迹并计算径向分布函数

计算RDF

在约0.1 nm处有一个很强的RDF峰，这个峰对应这每个水分子上的氢原子（氢氧键的键长约为0.1 nm）

计算配位数

接下来我们计算一下配位数，所谓配位数，也就是距离中心原子距离$r^{\prime }$ 以内的原子的数量。由以下公式给出 $$n(r^{\prime })={\int }_0^{r^{\prime }}g(r)\ast \rho \ast 4\pi r^{2}dr$$

由于RDF已经告诉了我们每一个球壳内原子的密度信息，将RDF乘以平均的原子密度$\rho$，就得到了球壳内原子的密度；接着乘以球壳的体积$4\pi r^{2}dr$，就是球壳内原子的数量；最后累加一下(积分)，就可以知道距离中心原子距离$r^{\prime }$ 以内的原子的数量。

由于我们的RDF是离散的，因此将积分改为求和，并且使用以下公式计算球壳体积即可 $$\begin{array}{rl}Vn(r)& =\frac{4}{3}\pi (r+\Delta r{)}^3-\frac{4}{3}\pi r^3\\ & =\frac{4\pi }{3}(r^3+3r^2\Delta r+3r\Delta r^2+\Delta r^3)-\frac{4}{3}\pi r^3\\ & =\frac{4}{3}\pi (3r^2\Delta r+3r\Delta r^2+\Delta r^3)\\ & \approx 4\pi r^2\Delta r\end{array}$$ 由于$\Delta r$ 很小，因此可以忽略其二阶和三阶项

$$n(r^{\prime })=\Sigma g(r)\ast \rho \ast 4\pi r^2\Delta r$$

实际上也就是是将积分改为求和，$dr$替换为了$\Delta r$

可以看见，曲线在0.1 nm处出现了明显的转折，对应的数值大约为2，也即对应了水分子中两个氢原子。

寻峰与绘图

在实际应用中，我们往往想知道的是离特定原子最近的一层原子的数量，这对应着到RDF第一个峰的峰谷距离以内的原子数量。总而言之，我们不仅要计算RDF，还需要找到它峰谷的位置。寻峰实际上是一个很具有一般性的问题，下面通过一个简单的例子说明如何寻峰。首先计算一段RDF。

这里我们直接调用scipy中信号处理的接口

1. 先将RDF曲线平滑，便于寻峰。
2. 将曲线反转，通过寻找峰的方式找到峰谷

也可以不反转峰，直接寻峰，然后再读取props中的'right_bases'，不过亲测这样做有时候并不能正确地找到峰谷，请读者谨慎尝试

可以看见，平衡之后曲线峰的位置不变，峰谷也找得较为准确。最后，可以将前面的过程简单的封装一下，便于复用。

封装完成后，只需使用MDtraj读入轨迹，并指定中心原子和参数原子，即可迅速计算出RDF与配位数，并绘制好图片，十分地便利。对于有其他需求的读者，可以自行修改函数。

关于RDF与配位数计算的讨论就到这里，感谢读者愿意花费时间阅读，如有问题或者建议，欢迎在评论区讨论。