博文
GROMACS分析教程：氢键分析
|
2017-01-09 16:56:34
卢天给过一个VMD的tcl脚本, 用于计算不同z位置水能形成的平均氢键数. 但轨迹大了以后, VMD的tcl脚本分析起来有点吃力, 且得到的结果与GROMACS的默认氢键标准存在差距. 这里我以此为例, 展示下如何组合GROMACS自带的工具进行复杂一点的分析, 并且利用gmx trjconv的外挂脚本功能让分析自动化. 示例所用GROMACS版本为5.1.4.
运行模拟示例文件用的是GROMACS自带的spc216.gro, 只不过沿z方向将盒子扩大了一些, 留出一定的真空层. 运行NVT模拟.
得到轨迹后, 先对轨迹进行预处理
选择0, 对整个体系进行PBC校正, 保证分子完整. 如果需要, 根据情况进行其他特殊处理.
分析单帧构型在对整个轨迹进行分析前, 先选择一帧构型测试命令和脚本的正确性, 这样可以节约时间.
选择0, 输出整个体系, 得到traj0.gro.
选择需要的原子组我们的最终目的是分析不同高度层内水分子所成的平均氢键数目, 要使用的命令是gmx hbond. 使用这个命令分析氢键时, 需要指定两个组, 它们可以相同也可以不同. 如果相同, 分析的是本组分子之间形成的氢键; 如果不同, 分析的则是两组分子之间形成的氢键.
对于我们的目的而言, 我们需要两个组, 组1 Wsel是某高度层内的水分子, 组2 Wexc是排除Wsel后的其余水分子. 这样我们就可以计算组1之间的氢键数目N(Wsel-Wsel), 组1和组2之间的氢键数目N(Wsel-Wesc), 而Wsel中的水分子所成的平均氢键数目 Nhb = ( 2*N(Wsel-Wsel)+N(Wsel-Wexc) )/N(water). N(Wsel-Wsel)之前的因子2是因为gmx hbond报告同组分子所成的氢键数目时会排除重复, 而我们不需要这样做.
选择不同高度层内的分子, 选择方式有几种策略, 或者根据氧原子的位置, 或者根据分子的质心位置, 几何中心位置, 或者根据任意原子位置等, 具体的语法看参考gmx select的语法及用法.
选中一定高度层内的水分子, 最简单的策略是根据氧原子的位置, 先选中符合条件的氧原子, 然后扩展到与氧原子所属分子相同的原子. 下面是的语句选择z坐标位于3到3.2之间的水分子以及其余水分子:
如果利用质心条件的话, 可以使用下面的语句
如果要同时选择Wsel和Wexc的话, 最好先定义一个变量, 这样写起来简洁, 执行起来效率也更高
最终, 我们可以使用下面的命令获得gmx hbond所需要的分组及其原子数目
traj0.xvg最后一行的第二列是Wsel中的原子数目, 其1/3就是我们所需要的N(water); traj0.ndx中保存了我们需要的两个分组.
对原子组进行分析获得了原子组的索引文件, 就可以用它进行氢键数目分析了.
运行上面的命令时会提示选择两个分组, 手动选择不利于自动化, 我们可以利用echo命令和管道来替代手动选择
traj0_00.xvg最后一行的第二列是我们需要的N(Wsel-Wsel), traj0_01.xvg最后一行的第二列则是我们需要的N(Wsel-Wexc).
整理报告分析结果有了这两个数据, 以及前面的水分子数目, 我们就可以计算出每个水分子所成的平均氢键数目了. 当然还是需要用脚本来自动计算, 这样后面才能实现自动化.
上面的脚本先获取每个文件的最后一行, 然后使用管道将其传awk, awk计算出所需要的值并保存到HB.xvg中.
通用单帧构型分析脚本将前面的命令写到一个bash脚本中, 就可以自动执行上面的分析了
| hb.bsh | |
|---|---|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 | i=0
file=traj$i
gmx select -f $file.gro -s -select 'Wsel= resname SOL and (res_com z 3 to 3.2); "Wsel" Wsel; "Wexc" not Wsel' -os $file.xvg -on $file.ndx
echo00 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_00.xvg
echo01 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_01.xvg
Nsel=$(tail -n 1$file.xvg)N00=$( tail -n 1${file}_00.xvg)N01=$( tail -n 1${file}_01.xvg)echo$Nsel$N00$N01 | awk '{print; print $1, (2*$5+$8)/($2/3) >>"HB.xvg"}'
rm -rf $file.gro $file.xvg $file.ndx ${file}_00.xvg ${file}_01.xvg
|
我们在最前面定义了帧号, 这样脚本就很容易用于其他帧了, 只要改变帧号对应的变量i就可以了. 脚本的最后我们删除了用到的中间文件, 这样既可以避免重复运行GROMACS工具时因备份文件太多而导致的错误, 也可以让我们的目录更清爽.
分析整条轨迹既然已经可以对一帧构型进行分析, 并完成了一个通用的分析脚本, 那么将脚本用于多个构型就比较简单了. 最直接的方式就是先使用gmx trjconv输出所有的构型, 然后循环处理每帧构型. 比如我们有帧号为0到100的多个构型文件, 那么在bash中可以使用
| hb.bsh | |
|---|---|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 | for i in {0..100}; dofile=traj$i
gmx select -f $file.gro <略> -os $file.xvg -on $file.ndx
echo00 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_00.xvg
echo01 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_01.xvg
Nsel=$(tail -n 1$file.xvg)N00=$( tail -n 1${file}_00.xvg)N01=$( tail -n 1${file}_01.xvg)echo$Nsel$N00$N01 | awk '{print; print $1, (2*$5+$8)/($2/3) >>"HB.xvg"}'
rm -rf $file.gro $file.xvg $file.ndx ${file}_00.xvg ${file}_01.xvg
done |
这种方式的缺点在于需要先输出一大堆构型文件, 当处理的帧数很多时, 就凌乱了. 好在 gmx trjconv支持一个外挂脚本的选项, -exec "命令", 可以在输出每一帧后对此帧构型执行指定的命令, 脚本或程序, 且以帧号作为命令行参数. 举例来说, 如果选项为-exec "cmd", 那么对输出每帧构型后对其执行的命令就是cmd 帧号. 利用这一功能, 我们就不再需要自己写循环了, 只要将我们上面的脚本改为以帧号为输入参数就可以了
| hb.bsh | |
|---|---|
1
2
3 | i=$1file=traj$i
<下同>
|
这样我们直接使用
就可以自动对整条轨迹进行分析, 得到所需要的数据了. 对我们的测试轨迹, 得到的平均数目为3.5左右, 符合预期.
总结说明gmx trjconv的-sep选项只支持gro或pdb格式, 不支持二进制格式构型文件.
这种组合GROMACS已有工具, 以及简单数据处理小脚本的分析方式不需要自己编写代码处理轨迹, PBC, 分析等问题, 只需要一些代码处理输出输入文件, 流程具有通用性, 但难度远小于自己从头写代码处理轨迹.
这种处理方法在分析时每次只处理一帧构型, 对于需要多帧构型才能计算的物理量就不适用了. 在那种情况下, 只能先输出所有的信息, 再统一处理了.
这种方式每次只处理一帧构型, 效率不是太高, 但对机器的要求低, 不像VMD那样要先载入整条轨迹, 需要很大内存.
如果你对于效率很在意, 那么可以采用并行的方式进行处理: 先输出所有的构型文件, 使用脚本并行处理左右构型文件, 最后再将这些文件的结果合并. 至于如何并行处理, 这里就不再细说了, 可参考我以前的两篇博文Bash脚本实现批量作业并行化, GNU Parallel.
如果你想试验一下这种处理方式, 那可以试着完成我前面提出的一个问题, 蛋白口袋或纳米通道内水分子的个数统计. 当然, 这些口袋或通道都是柔性的, 否则的话, 就没有必要这么处理了.
使用VMD的tcl脚本进行分析时, 思路也是类似的. 但值得注意的是, VMD的默认氢键标准与GROMACS不同, 这是因为目前存在多种判断氢键的标准. 简言之, VMD的 3.5埃-40度 标准大致和GROMACS的标准一致. 此外, VMD还没有考虑PBC的问题. 详情可见前一篇博文GROMACS的默认氢键标准.
| HB.tcl | |
|---|---|
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 | #mol delete all#mol new conf.gro traj_pbc.xtcset sel "same resid as resname SOL and z 30 to 32 "set Zset "z 30 to 32 and x 4 to 14 and y 4 to 14"set Wsel [atomselect top "same resid as name OW and $Zset"]set Wexc [atomselect top "same resid as name OW and not($Zset)"]set Nfrm 100set Rcut 3.5set Angle 40.
puts"#Frm #Wsel #Wexc #sel-sel #sel-exc #exc-sel #HB/Wat"for{set i 0}{$i<=$Nfrm}{incr i}{$Wselframe$i$Wsel update
$Wexc frame $i$Wexc update
set Nsel [expr[$Wselnum]/3.]set Nexc [expr[$Wexcnum]/3.]set Nss [llength[lindex[measure hbonds $Rcut$Angle$Wsel]0]]set Nse [llength[lindex[measure hbonds $Rcut$Angle$Wsel$Wexc]0]]set Nes [llength[lindex[measure hbonds $Rcut$Angle$Wexc$Wsel]0]]puts"$i $Nsel $Nexc $Nss $Nse $Nes [expr (2.*$Nss+$Nse+$Nes)/$Nsel]"} |
◆本文地址: http://jerkwin.github.io/2017/01/09/GROMACS分析教程-氢键分析/, 转载请注明◆
https://blog.sciencenet.cn/blog-548663-1026522.html
上一篇：GROMACS和VMD中的氢键判定标准
下一篇：我使用的VMD初始化脚本
