GROMACS计算距离的方法及注意点

• 2018-01-14 20:27:30

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

1
2

[ dist ]
3 5 3 9 3 13 3 19

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

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

 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

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

运行脚本, 指定输出文件名称和两个组的名称, 这段代码会生成两个组的原子对组合, 并将其命名为一个新组追加到原来的索引文件中. 然后运行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)语法及用法.