计算任意一点周围指定半径内的目标分子数密度分布

对于A、B两类分子来说，要统计A分子周围B分子的分布概率，可以计算RDF来表征，使用GROMACS的rdf命令工具即可计算。但是，如果想计算空间任意一点为球心，其指定半径内的目标分子数密度分布，GROMACS还不支持。我们写了一个代码，以实现这个功能。可以加入模拟之家QQ群获取：709020941。      

1.分析单个结构        

运行下面的脚本命令，即可计算：      

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

nvt.gro为需要计算的整个模拟盒子；SOL即为残基名，为需要统计分布的分子，可以改为目标分子的残基名（需要和gro文件里面对应）；6.2 6.2 6.2为所定义的球中心坐标xyz值；6为所需要计算的球体半径。执行代码后会生成SOL_densi_num.txt文件，里面记录了SOL分子在以6.2 6.2 6.2为中心，6nm半径内的数密度分布。      

2.分析多个结构        

上述命令只能计算一个结构的分布，但是对于分子动力学模拟来说不具有统计意义，需要计算多帧结构的平均值才能代表最终的数值。基于此，可以使用shell脚本循环调用该py程序来实现。      

调用之前，需要使用GROMACS的轨迹提取命令，把轨迹转成单帧的gro格式的结构，运行：      

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

上述命令中，-sep 即为把每一帧结构单独储存为N*.gro格式，这里的*为每一帧的ID。然后使用下面命令循环调用py脚本计算平均值：（这里N=50需要修改为自己生成的gro文件个数）      

Name=SOLNum=50rm tmp1touch tmp1for((i=0;i<=$Num;i++))dopython3 num_density.py N${i}.gro $Name 6.2 6.2 6.2 6sed -e "1,$ p" -n ${Name}_densi_num.txt | awk '{print $2}' > tmp2paste tmp2 tmp1 > tmp3mv tmp3 tmp1doneawk 'sum=0;{for(i=1;i<=NF;i++)sum+=$i;print sum/NF;}' tmp1 > tmp4sed -e "1,$ p" -n ${Name}_densi_num.txt | awk '{print $1}' > tmp5paste tmp5 tmp4 > ${Name}_avg.txtrm tmp*

上述运行完之后，得到SOL_avg.txt的文件，里面记录了SOL分子在以6.2 6.2 6.2为中心，6nm半径内的数密度分布。如下图所示：