应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » Gromacs/Amber/NAMD » 【讨论】关于不同z值处的水能形成氢键数的脚本讨论
51/1返回列表
查看: 1561  |  回复: 4
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

jianchaoyv

金虫 (小有名气)

[交流] 【讨论】关于不同z值处的水能形成氢键数的脚本讨论已有3人参与

《转载》
氢键判据用的是常用的35度3.5埃的几何判据，当然也可以直接在脚本里改。计算方法是，只要有水分子有一个原子在某一层里，则这个水分子就认为属于这一层的水。对于每一帧，计算属于每一层的水selin与其它物质selbig之间的氢键数，氢键包括了这一层中的水作为氢键受体和供体两种情况，其数目分别为代码中的变量a和b。并且加上这一层水内部之间的氢键数（变量c）的2倍。a+b+2c除以这一层的水数，作为这一帧这一层的每个水的平均氢键数。脚本中循环轨迹中的每一帧，最终得到这一层平均氢键数。nonum变量记录有多少帧在所设定的范围里没有水，这些帧不计算。#后面那行用于调试目的，要考察每帧结果就去掉开头的#。
首先运行下面的脚本，来加载实现这个功能的子程序
proc numhbavg {sel fps1 fps2} {
set selin [atomselect top $sel]
set selbig [atomselect top "same resid as exwithin 3.5 of $sel"]
set k 0.0
set nonum 0
for {set i $fps1} {$i<=$fps2} {incr i} {
$selin frame $i
$selin update
$selbig frame $i
$selbig update
if {[$selin num]!=0} {
set a [llength [lindex [measure hbonds 3.5 35 $selbig $selin] 0]] #measure hbonds 3.5 35 $selbig $selin这个命令返回什么样的值？怎么又先用了lindex、llength命令？
set b [llength [lindex [measure hbonds 3.5 35 $selin $selbig] 0]] #这句跟上句的区别？
set c [llength [lindex [measure hbonds 3.5 35 $selin] 0]]  #这句为何只有一个$selin？一般都是像上2句那样有2个$selin $selbig或$selbig $selin？
set k [expr $k+($a+$b+2*$c)*3.0/[$selin num]]
#puts "fps:$i $a+$b+[expr 2*$c] num_water:[expr [$selin num]/3.0] avg:[expr $k+($a+$b+2*$c)*3.0/[$selin num]]"
} else {incr nonum}
}
if {[expr $fps2-$fps1+1]==$nonum} {return "no result"}
return [expr $k/[expr $fps2-$fps1+1-$nonum]]
}
然后下面的循环会调用这个子程序来输出每一层的平均氢键数，这里假设要计算z=4.0~5.6埃的数据，间隔为0.1埃，且限定20 for {set i 40} {$i<=55} {incr i} {
set k [expr $i*0.1]
set now [numhbavg "same resid as resname SOL and x<40 and x>20 and y<40 and y>20 and z<[expr $k+0.1] and z>=$k" 100 150]
puts [format "%4.2f %4.2f %5.3f" $k [expr $k+0.1] $now]
}
我这里随便算一个主要由水构成的普通的体系，水形成的氢键数大概在3.1左右，标准放宽到4.0埃，40度，则可形成氢键数约为3.6。在冰中由于结构十分有序，可形成4个氢键，在液态情况下分子的动能必然造成氢键的破坏，所以结果是很合理的。
帧数范围越大、xy平面越大计算越慢，这个脚本计算速度比较慢，不要一下将范围设得太大。
输出结果如下，前两列代表统计的z值范围，第三列是水的平均氢键数
4.00 4.10 3.015
4.10 4.20 3.161
4.20 4.30 3.201
4.30 4.40 3.159
4.40 4.50 3.237
4.50 4.60 3.130
4.60 4.70 3.201
4.70 4.80 3.338
4.80 4.90 3.201
4.90 5.00 3.041
5.00 5.10 3.122
5.10 5.20 3.182
5.20 5.30 3.160
5.30 5.40 3.309
5.40 5.50 3.189
5.50 5.60 3.327

请各位指点一下，我把自己看不懂的写在上面了，不过可以对程序解释的越详细越好，更便于本人的正确理解！！期待各位高手指点！尤其是bay__gulf ！
另：氢键几何判据中角度到底指的哪个？vmd中的measure hbonds 的命令中的角度指的是：“the angle formed by the donor, hydrogen, and acceptor must be less than angle from 180 degrees. “   而我看到的文献Phys. Chem. Chem. Phys. , 2004, 6, 829–835 中氢键几何判据：“The bond angle a between the O–O direction and the molecular O–H direction of the donor, where H is the hydrogen which forms the bond, is lower than a threshold angle aC ”这两个好像不一致啊？

[ Last edited by jianchaoyv on 2010-6-29 at 10:54 ]
回复此楼

» 猜你喜欢

» 本主题相关商家推荐: (我也要在这里推广)

» 本主题相关价值贴推荐，对您同样有帮助:
1楼2010-06-29 10:41:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

coolrainbow

木虫 (著名写手)

未来国家冻凉
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
这是sob写的那个脚本？
赞一下(1人) 回复此楼
技术博客：http://hi.baidu.com/coolrainbow/blog
2楼2010-06-29 11:03:21
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bay__gulf

金虫 (著名写手)

刘苏州
★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
lei0736(金币+2):谢谢 2010-06-29 11:31:54
氢键的角度有两种定义方式, 如图所示
赞一下(2人) 回复此楼
3楼2010-06-29 11:20:23
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)
程序来源于http://hi.baidu.com/sobereva/blo ... 595ef7fc037fbf.html
麻烦最上面程序中还有问题请教！请给解释一下
另外，bay__gulf 那个vmd图画的与vmd的ug中描述好像不符吧？？

[ Last edited by jianchaoyv on 2010-6-29 at 14:43 ]
赞一下 回复此楼
4楼2010-06-29 14:15:05
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bay__gulf

金虫 (著名写手)

刘苏州
★ ★ ★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
ghcacj(金币+4):谢谢 2010-06-29 15:07:14
hbonds cutoff angle selection1 [selection2]:  Find all hydrogen bonds in the given selection(s), using simple geometric criteria. Donor and acceptor must be within the cutoff distance, and the angle formed by the donor, hydrogen, and acceptor must be less than angle from 180 degrees. Only non-hydrogen atoms are considered in either selection. If both selection1 and selection2 are given, the selection1  is considered the donor and selection2 is considered the acceptor. If only one selection is given, all non-hydrogen atoms in the selection are considered as both donors and acceptors. The two selections must be from the same molecule. The function returns three lists; each element in each list corresponds to one hydrogen bond. The first list contains the indices of the donors, the second contains the indices of the acceptors, and the third contains the index of the hydrogen atom in the hydrogen bond.

Known Issue: The output of hbonds cannot be considered 100% accurate if the donor and acceptor selection share a common set of atoms.
=============

$selbig $selin这个命令返回什么样的值？怎么又先用了lindex、llength命令？
请在上面找答案

#这句跟上句的区别？
逐字比对

这句为何只有一个$selin？一般都是像上2句那样有2个$selin $selbig或$selbig $selin？
请在上面找答案

氢键几何判据中角度到底指的哪个？
请在上面找答案

vmd 的ug 较长, 全部打印不现实
但8，10，11三章一定要打印下来仔仔细细看
赞一下(2人) 回复此楼
5楼2010-06-29 15:05:34
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 jianchaoyv 的主题更新
51/1返回列表
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭