仓库源文，站点原文

layout: post title: GROMACS计算距离的方法及注意点 categories:

科 tags:
gmx

2018-01-14 20:27:30

低版本的GROMACS中(具体哪个版本就没有考古了, 猜想是3.x吧), g_distance计算距离时需要选择两个组, 然后程序会自动计算这两个组的原子两两配对之间的距离, 也可以计算这两个组的质心之间的距离. 可能是从4.x版本起吧, gmx distance使用的索引组的方法有了变化, 可以提供多个索引组, 每个组是独立的, 在其中列出要计算的原子对. 因此要计算距离的组中总原子个数必须是偶数, 否则会给出原子数目非偶数的错误. 例如, 我们要计算体系中3号原子和5, 9, 13, 19之间的距离, 那么可以在索引文件中定义如下组

<table class="highlighttable"><th colspan="2" style="text-align:left">index.ndx</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%">1 2</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%">[ dist ] 3 5 3 9 3 13 3 19 </pre></div> </td></tr></table>

编号对也可以放于多行中, 每行一个原子对编号. 但值得注意的是, 至少在我测试的 2016.4版本中, 在处理这些编号对时, GROMACS存在错误, 对1 2 2 3这样的编号对, 会先将中间重复的编号合并为一个, 变成1 2 3, 从而导致无法计算. 解决办法也很简单, 就是更换顺序, 保证相邻编号之间不存在重复. 如上面的例子改成1 2 3 2就可以正常计算了.

采用这种模式, 需要我们自己提供要计算的原子对. 当原子数较少时并无困难, 但原子数一多起来, 组合就多了, 手写不太现实, 所以我就写了一段代码来处理这个问题.

<table class="highlighttable"><th colspan="2" style="text-align:left">gmx_ndx.bsh</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34</pre></div></td><td class="code"><div class="highlight"><pre style="line-height:125%">usage="\ >>>>>>>>>>>>>>>> gmx_ndx <<<<<<<<<<<<<<<< >>>>>>>>>>>>>>>> Jicun LI <<<<<<<<<<<<<<<< >>>>>>>>>> 2016-09-15 09:49:36 <<<<<<<<< >> Usage: gmx_ndx <File.ndx> [Group 1] [Group 2] >> Default: gmx_ndx N/A Other SOL" [[ $# -lt 1 ]] && { echo "$usage"; exit; } File=index; [[ $# -ge 1 ]] && { File=${1%.ndx}; } Grp1=Other; [[ $# -ge 2 ]] && { Grp1=$2; } Grp2=SOL; [[ $# -ge 3 ]] && { Grp2=$3; } echo -e "\n生成 $Grp1 与 $Grp2 的组合新组 [ ${Grp1}_$Grp2 ] 并追加到 $File.ndx 文件" awk -v Grp1=$Grp1 -v Grp2=$Grp2 ' BEGIN{ RS="[" } $1==Grp1 { N1=NF-2 for(i=1; i<=N1; i++) G1[i]=$(i+2) } $1==Grp2 { N2=NF-2 for(i=1; i<=N2; i++) G2[i]=$(i+2) } END { print "[ "Grp1"_"Grp2" ]" for(i=1; i<=N1; i++) { for(j=1; j<=N2; j++) printf "%d %d ", G1[i], G2[j] print "" } } ' $File.ndx >> $File.ndx </pre></div> </td></tr></table>

运行脚本, 指定输出文件名称和两个组的名称, 这段代码会生成两个组的原子对组合, 并将其命名为一个新组追加到原来的索引文件中. 然后运行gmx distance -s -f -n -oav -oall, 选择这个新组, 就可以得到所要的距离文件了. 值得注意的是, 组中的原子对数目最好不要太多, 否则GROMACS运行时所需的内存太多, 可能导致失败.

使用这种模式计算原子之间的距离很简单, 但要想计算两个组的质心之间的距离就有点麻烦, 需要利用选区语法来完成. 例如我们要计算残基1-4与残基7-9质心之间的距离, 可以使用下面的方式:

gmx distance -s -f -select "com of resnr 1 to 4 plus com of resnr 7 to 9" -oav -oall

其中的选区语法请参考GROMACS选区(selection)语法及用法.