GROMACS分析教程:氢键分析

类别:    标签: gmx vmd 氢键   阅读次数:   版权: (CC) BY-NC-SA

卢天给过一个VMD的tcl脚本, 用于计算不同z位置水能形成的平均氢键数. 但轨迹大了以后, VMD的tcl脚本分析起来有点吃力, 且得到的结果与GROMACS的默认氢键标准存在差距. 这里我以此为例, 展示下如何组合GROMACS自带的工具进行复杂一点的分析, 并且利用gmx trjconv的外挂脚本功能让分析自动化. 示例所用GROMACS版本为5.1.4.

运行模拟

示例文件用的是GROMACS自带的spc216.gro, 只不过沿z方向将盒子扩大了一些, 留出一定的真空层. 运行NVT模拟.

gmx grompp
gmx mdrun

预处理轨迹

得到轨迹后, 先对轨迹进行预处理

gmx trjconv -f traj_comp.xtc -pbc whole -o traj_pbc.xtc

选择0, 对整个体系进行PBC校正, 保证分子完整. 如果需要, 根据情况进行其他特殊处理.

分析单帧构型

在对整个轨迹进行分析前, 先选择一帧构型测试命令和脚本的正确性, 这样可以节约时间.

gmx trjconv -f traj_pbc.xtc -o traj.gro -dump 0 -sep

选择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" same mol as (name OW and (3<z and z<3.2))
"Wexc" not same mol as (name OW and (3<z and z<3.2))

如果利用质心条件的话, 可以使用下面的语句

resname SOL and (3<res_com z and res_com z<3.2)
或 resname SOL and (res_com z 3 to 3.2)

如果要同时选择Wsel和Wexc的话, 最好先定义一个变量, 这样写起来简洁, 执行起来效率也更高

Wsel = resname SOL and (res_com z 3 to 3.2); "Wsel" Wsel; "Wexc" not Wsel

最终, 我们可以使用下面的命令获得gmx hbond所需要的分组及其原子数目

gmx select -f traj0.gro -s -select 'Wsel = resname SOL and (res_com z 3 to 3.2); "Wsel" Wsel; "Wexc" not Wsel' -os traj0.xvg -on traj0.ndx

traj0.xvg最后一行的第二列是Wsel中的原子数目, 其1/3就是我们所需要的N(water); traj0.ndx中保存了我们需要的两个分组.

对原子组进行分析

获得了原子组的索引文件, 就可以用它进行氢键数目分析了.

gmx hbond -f traj0.gro -n traj0.ndx -num traj0.xvg

运行上面的命令时会提示选择两个分组, 手动选择不利于自动化, 我们可以利用echo命令和管道来替代手动选择

echo 0 0 | gmx hbond -f traj0.gro -n traj0.ndx -num traj0_00.xvg
echo 0 1 | gmx hbond -f traj0.gro -n traj0.ndx -num traj0_01.xvg

traj0_00.xvg最后一行的第二列是我们需要的N(Wsel-Wsel), traj0_01.xvg最后一行的第二列则是我们需要的N(Wsel-Wexc).

整理报告分析结果

有了这两个数据, 以及前面的水分子数目, 我们就可以计算出每个水分子所成的平均氢键数目了. 当然还是需要用脚本来自动计算, 这样后面才能实现自动化.

Nsel=$(tail -n 1 traj0.xvg)
N00=$( tail -n 1 traj0_00.xvg)
N01=$( tail -n 1 traj0_01.xvg)
echo $Nsel $N00 $N01 | awk '{print; print $1, (2*$5+$8)/($2/3) >>"HB.xvg"}'

上面的脚本先获取每个文件的最后一行, 然后使用管道将其传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

echo 0 0 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_00.xvg
echo 0 1 | 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}; do
	file=traj$i
	gmx select -f $file.gro <略> -os $file.xvg -on $file.ndx

	echo 0 0 | gmx hbond -f $file.gro -n $file.ndx -num ${file}_00.xvg
	echo 0 1 | 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=$1
file=traj$i
<下同>

这样我们直接使用

gmx trjconv -f traj_pbc.xtc -o traj.gro -sep -exec "bash hb.bsh"

就可以自动对整条轨迹进行分析, 得到所需要的数据了. 对我们的测试轨迹, 得到的平均数目为3.5左右, 符合预期.

总结说明

附: 类似功能的VMD tcl脚本

使用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.xtc

set 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 100
set Rcut  3.5
set Angle 40.

puts "#Frm  #Wsel #Wexc  #sel-sel #sel-exc #exc-sel   #HB/Wat"

for {set i 0} {$i<=$Nfrm} {incr i} {
	$Wsel frame $i
	$Wsel update
	$Wexc frame $i
	$Wexc update

	set Nsel [expr [$Wsel num]/3.]
	set Nexc [expr [$Wexc num]/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]"
}
随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: GROMACS和VMD中的氢键判定标准
后一篇: 我使用的VMD初始化脚本

访问人次(2015年7月 9日起): | 最后更新: 2017-06-23 03:12:17 UTC | 版权所有 © 2008 - 2017 Jerkwin