刘建川・GROMACS 实战 | 任意点周围指定半径内分子数密度分布计算

大家好，我是刘建川，14 年深耕 GROMACS 分子模拟领域，始终专注分子模拟实战。在分子动力学模拟结果分析中，统计两类分子间的空间分布规律，可通过径向分布函数（RDF）实现，使用 GROMACS 自带的 gmx rdf 工具即可完成。但如果需要以空间中任意一点为球心，统计指定半径范围内目标分子的数密度分布，GROMACS 原生工具暂不支持。我们提供了自定义脚本实现该功能，可加入模拟之家 QQ 群（709020941）获取脚本文件。

一、单帧结构数密度计算

对于单帧结构文件，直接调用 Python 脚本即可完成计算，执行命令如下：

python num_density.py nvt.gro SOL 6.2 6.2 6.2 6

参数说明：

• nvt.gro：待计算的模拟体系结构文件
• SOL：待统计的目标分子残基名，需与 gro 文件内的残基名称保持一致
• 6.2 6.2 6.2：自定义球体中心的 x、y、z 坐标
• 6：统计球体的半径，单位为 nm

命令执行完成后，会生成 SOL_densi_num.txt文件，记录目标分子在指定球心、指定半径范围内的数密度分布数据。

二、多帧轨迹平均统计

单帧结构的计算结果不具备统计意义，分子动力学模拟通常需要对多帧轨迹取平均值，才能得到可靠的统计结果。可通过拆分轨迹 + Shell 循环调用脚本的方式实现。

2.1 拆分轨迹为单帧结构

首先使用 gmx trjconv 命令，将完整轨迹拆分为逐帧保存的 gro 结构文件：

gmx trjconv -f md.xtc -s md.tpr -sep -o N.gro

其中 -sep参数表示将每一帧结构单独保存，输出文件名为 N0.gro、N1.gro……依次按帧编号命名。

2.2 批量计算并取平均值

编写 Shell 循环脚本，批量调用 Python 脚本处理所有单帧文件，最终输出平均后的数密度分布。脚本内容如下，使用前需将 Num=50修改为实际生成的 gro 文件总数：

Name=SOL
Num=50
rm tmp1
touch tmp1
for((i=0;i<=$Num;i++))
do
python3 num_density.py N${i}.gro $Name 6.2 6.2 6.2 6
sed -e "1,$ p" -n ${Name}_densi_num.txt | awk '{print $2}' > tmp2
paste tmp2 tmp1 > tmp3
mv tmp3 tmp1
done
awk 'sum=0;{for(i=1;i<=NF;i++)sum+=$i;print sum/NF;}' tmp1 > tmp4
sed -e "1,$ p" -n ${Name}_densi_num.txt | awk '{print $1}' > tmp5
paste tmp5 tmp4 > ${Name}_avg.txt
rm tmp*

脚本运行完成后，会生成 SOL_avg.txt文件，里面记录了多帧平均后的目标分子数密度分布结果。

完成以上操作后，得到的统计结果具备充分的统计意义，可用于分析特定位点、界面区域、限域空间内的分子聚集与分布特征，为相关机制研究提供量化的数据支撑。