- 2021-06-19 08:19:44
前一篇文章介绍了VMD氢键插件的使用, 这篇文章就再接再厉, 说明下GROMACS氢键分析工具hbonds的使用, 顺便发布一个gmx_hbdat脚本, 使其能够整理氢键数据.
无论准备使用GMX的哪个工具, 首先要做的是查看其文档, 以及GMX手册中的相关说明, 确保能理解工具给出的结果. 对有些工具, 文档和手册中的说明太多简略, 有时候还必须查看相关的源代码. 这是使用GMX工具的必修课, 切记.
GMX的氢键工具hbond选项众多, 细究起来还是比较复杂的, 我们这里只介绍最基本的功能, 也就是分析一组原子内, 或两组原子间的氢键. 这需要指定两个索引组. 如果两个索引组完全相同, 那得到的就是这个索引组内部的氢键; 如果两个索引组完全不同（即两组不会同时包含相同的原子）, 那得到的就是两个索引组之间的氢键. 这是目前仅支持的两种情况. 如果两个索引组既不完全相同, 却又部分相同, 那就无法计算, 结果当然也不可靠.
分析结果主要涉及三个文件,
hbnum.xvg: 每个时刻的氢键数目, 只满足距离标准的氢键数目,hbond.ndx: 列出整个分析轨迹中所有可能氢键涉及的原子编号. 注意编号从1开始, 不同于VMD的从0开始.hbmap.xpm: 所有可能氢键的时间演化, 给出每一分析帧中每个氢键是否存在.
如果查阅MD文章中的氢键分析部分, 可以发现很多文章给出的结果表格都是AMBER格式的, 其中会给出每个氢键所涉及的残基/原子, 每个氢键的占有率. 我们利用GMX的拓扑文件以及后两个文件, 容易整理得到这样的氢键数据. hbond.ndx中给出了原子编号, 根据编号查看gro文件中的信息, 或者拓扑中的数据, 就可以得到相关原子的残基, 原子类型, 名称等数据. 处理下hbmap.xpm文件就可以计算出每个氢键的占有率.
简单示例
继续以上篇文章的体系做示例吧. 这是一个包合物体系, 我们计算两个分子间的氢键. 先使用make_ndx做个索引文件, 其中包含两个分子各自的索引组. 然后运行命令
gmx hbond -f traj.xtc -s topol.tpr -n index.ndx -num -hbn -hbm提示时, 选两个分子各自的索引组, 我用的是MOLA, MOLB.
运行完, 就得到了三个文件输出文件.
| hbnum.xvg | |
|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | @ title "Hydrogen Bonds"
@ xaxis label "Time (ps)"
@ yaxis label "Number"
@TYPE xy
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend "Hydrogen bonds" # 同时满足距离标准, 角度标准的氢键数
@ s1 legend "Pairs within 0.35 nm" # 只满足距离标准的氢键数
0 2 9
0.02 2 8
0.04 2 6
0.06 2 6
0.08 2 6
0.1 2 6
0.12 2 6
... ... |
| hbond.ndx | |
|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | [ MOLA ] # MOLA的所有原子
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
... ...
[ donors_hydrogens_MOLA ] # MOLA的供体原子, 氢原子
8 9
12 13
22 23
... ...
[ acceptors_MOLA ] # MOLA的受体原子
3 8 12 16 22 24 29 33 37 43 45 50 54 58 64
... ...
[ MOLB ] # MOLB的所有原子
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
... ...
[ donors_hydrogens_MOLB ] # MOLB的供体原子, 氢原子
148 201
[ acceptors_MOLB ] # MOLB的受体原子
148 149
[ hbonds_MOLA-MOLB ] # 整段分析轨迹中MOLA和MOLB之间所有可能的氢键
64 65 149
146 147 148
148 201 43 |
| hbmap.xpm | |
|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | /* XPM */
/* This file can be converted to EPS by the GROMACS program xpm2ps */
/* title: "Hydrogen Bond Existence Map" */
/* legend: "Hydrogen Bonds" */
/* x-label: "Time (ps)" */
/* y-label: "Hydrogen Bond Index" */
/* type: "Discrete" */
static char *gromacs_xpm[] = {
"101 3 2 1",
" c #FFFFFF " /* "None" */, # 无氢键, 空白
"o c #FF0000 " /* "Present" */, # 有氢键, 红色o
/* x-axis: 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 0.66 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 1.08 1.1 1.12 1.14 1.16 1.18 1.2 1.22 1.24 1.26 1.28 1.3 1.32 1.34 1.36 1.38 1.4 1.42 1.44 1.46 1.48 1.5 1.52 1.54 1.56 1.58 */
/* x-axis: 1.6 1.62 1.64 1.66 1.68 1.7 1.72 1.74 1.76 1.78 1.8 1.82 1.84 1.86 1.88 1.9 1.92 1.94 1.96 1.98 2 */
/* y-axis: 0 1 2 */
# 三条氢键, 与前面`hbond.ndx`中的 [ hbonds_MOLA-MOLB ] 对应
"ooooooooooooooooo o ooooo oo oooooooo oo ooooooooooooooooooooooooooooooooooooooooooooooooooooo",
" oooooo ooo oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo",
"ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo" |
使用bash脚本hbdat.bsh获取类似AMBER的氢键数据
bash hbdat.bsh得到hbdat.dat
| hbdat.dat | |
|---|---|
1 2 3 4 | #Donor Hydrogen Acceptor Occupancy%
64#"MOLA"_MOL@OH 65#"MOLA"_MOL@HO 149#"MOLB"_MOL@O_3 100.0
146#"MOLA"_MOL@OH 147#"MOLA"_MOL@HO 148#"MOLB"_MOL@OH 70.3
148#"MOLB"_MOL@OH 201#"MOLB"_MOL@HO 43#"MOLA"_MOL@OH 87.1 |
其中原子的标识为 总体编号#"分子名称"_残基名称@原子名称. 此数据可以与前一篇文章中的数据对比, 可以看到氢键所涉及的原子完全一致, 但由于VMD和GMX的氢键判断标准无法精确地相互转换, 所以占有率只是基本一致, 无法完全一致.
其他
分析氢键时的另一个常见问题是, 如何分析GMX不支持原子间的氢键, 如涉及卤素原子的氢键. 这可以通过修改拓扑文件, 将卤素原子改为O或N, 骗过GMX来达到目的. 因为GMX不会检查实际原子是什么元素, 只看其名称. 修改了原子名称GMX就会以为原子真的变了, 即便名不符实. 但注意不要把拓扑搞混了, 这种修改过的拓扑只能用来计算氢键, 不要用它来做其他计算.
