径向分布函数模拟计算

我第一次做 径向分布函数模拟计算是分析一段 10 ns 水的轨迹，bin 宽度设了 0.1 Å，画出来的 g_OO(r) 锯齿状乱抖，我还以为体系结构异常。后来把 bin 调到 0.02 Å 并跑满 50 ns 平均，曲线立刻平滑成教科书形状。那次让我记住，径向分布函数模拟计算的曲线质量一半在统计，一半在参数。

一、为什么 径向分布函数模拟计算先做参数

g(r) 是直方图统计，bin 宽度、最大统计距离 r_max、参与帧数和粒子数共同决定曲线平滑度和可信范围。bin 太宽丢细节、太窄涨落大；r_max 超过盒子一半会因 PBC 自相关污染。我做 径向分布函数模拟计算时，默认 bin 0.02~0.05 Å、r_max < L/2、轨迹至少 20 ns 且多帧采样。曾经 bin 0.1 把水第二峰削平，误判”长程有序弱”。

二、核心原理：从轨迹到直方图

每帧统计所有粒子对落在 [r, r+dr] 壳层的数量，除以总壳层体积和平均密度得到 g(r)，跨帧平均降噪。关键认知：g(r) 的峰位由结构决定、峰高由温度和密度决定，所以对比不同体系必须同温同密度。我用 g(r) 校验力场时，先看第一峰位和配位数是否贴近实验散射数据，再谈别的。

三、关键技术要点：PBC 与混匀

径向分布函数模拟计算最致命的是 r_max 超过 L/2，远壳层统计会把自己盒子的镜像算进来，长程 g(r) 失真。我默认 r_max 取 0.45L 并对边界做最小镜像修正。另一个坑是体系未平衡就统计，早期弛豫的结构混进平均，第一峰被拉偏。还有人只用重原子算却忘了不同元素要分 g_AB 单独算，异类配位信息全混在同类里。

四、实操流程：从轨迹到 g(r)

我的标准动作：轨迹处理（PBC 修复、去水、居中）→ 定 bin 与 r_max → 多帧统计平均 → 积分配位数 → 和实验/X 射线对标。最容易被砍的是平衡判断，但 径向分布函数模拟计算 的峰值可信度全靠它，我默认先看 RMSD 确认平衡再统计。

补一个实算：Ca²⁺ 水溶液 g_CaO(r)，10 ns 平均第一峰 2.4 Å、配位数 7.5 还在抖；跑到 50 ns 收敛到 2.42 Å、配位数 8.1 稳定，和实验八配位一致。如果只交 10 ns 结果，你会误判配位数”约 7~8 不定”。这之后我做 g(r)，轨迹长度先按峰位收敛为准，不再固定 10 ns。

五、常见踩坑：分辨率与采样不足

径向分布函数模拟计算常见错是盒子太小（< 30 Å）导致统计粒子少，第三壳层纯噪声；或只取 100 帧间隔过大，自相关没打破。另一个坑是把 g(r) 当唯一结构判据，忽略它给不了角度信息（需配角分布函数），误以为”峰在就是对的”。

六、复盘：g(r) 是力场的试金石

回过头看，径向分布函数模拟计算最值钱的是用最低成本验证”你的力场画出来的液体像不像真的”，它是力场试金石。我接项目默认拿 g(r) 先对标再往下做性质。把 bin、r_max、平衡三件事做对，曲线才敢拿去和散射实验比。被证明有用的，是愿意为多跑 40 ns 让配位数收敛的耐心。下次有人给我一张 g(r)，我会先问——你的 bin 和轨迹长度够吗。