24小时热门版块排行榜    

查看: 1910  |  回复: 4

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 ]
回复此楼

» 猜你喜欢

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

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

coolrainbow

木虫 (著名写手)

未来国家冻凉


小木虫(金币+0.5):给个红包,谢谢回帖交流
这是sob写的那个脚本?
技术博客: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
氢键的角度有两种定义方式, 如图所示
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三章一定要打印下来仔仔细细看
5楼2010-06-29 15:05:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 jianchaoyv 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[教师之家] 售T0P一区SCI文章,我:8O5.51.O.54,科目齐全,可+急 +3 68gx6bh3ei 2026-07-09 3/150 2026-07-10 06:53 by pbuc7it0i8
[硕博家园] 售T0P一区SCI文章,我:8O5.51.O.54,科目齐全,可+急 +3 68gx6bh3ei 2026-07-09 3/150 2026-07-10 06:33 by pbuc7it0i8
[硕博家园] 售T0P一区SCI文章,我:8O5.51.O.54,科目齐全,可+急 +3 ph7bxbho8f 2026-07-09 3/150 2026-07-10 05:53 by j4frgf6zcw
[教师之家] 售T0P一区SCI文章,我:8O5.51.O.54,科目齐全,可+急 +3 al9t2f50x3 2026-07-09 4/200 2026-07-10 05:33 by j4frgf6zcw
[考研] 售T0P一区SCI文章,我:8O5.51.O.54,科目齐全,可+急 +4 al9t2f50x3 2026-07-09 4/200 2026-07-10 05:33 by j4frgf6zcw
[基金申请] b口会评 +6 当空照 2026-07-08 6/300 2026-07-09 20:49 by tonyliu0401
[基金申请] 信息学部会评时间 +3 xuel2011 2026-07-06 3/150 2026-07-09 19:05 by mollyzhang_2003
[基金申请] 祈祷青基必中 50+3 hadengry 2026-07-09 9/450 2026-07-09 17:37 by 文杰猪猪
[考研] 在职考研 +6 wu5d 2026-07-08 7/350 2026-07-09 16:40 by 化工小霸王呀
[论文投稿] MSER送审了还被拒稿 +3 S778950906 2026-07-08 5/250 2026-07-09 15:46 by S778950906
[基金申请] 关于如何从代码看上不上会 +11 且听虎啸 2026-07-09 14/700 2026-07-09 15:27 by libeny
[有机交流] 硝基还原 20+3 暮年飞雪 2026-07-07 5/250 2026-07-09 13:30 by tfang
[基金申请] 今年E04面上 +15 布布和一二 2026-07-05 21/1050 2026-07-09 10:25 by 学员fOzRO9
[硕博家园] 大龄残疾硕士的一点执念 +14 3776737374 2026-07-06 22/1100 2026-07-09 06:07 by lovechizi
[基金申请] 生命口会评 +4 wuke100666 2026-07-07 4/200 2026-07-08 16:38 by 逍遥不够自在
[教师之家] 网络举报学术不端成为了男女之间情杀的手段了,真荒唐。。。 +3 zju2000 2026-07-04 3/150 2026-07-07 21:00 by hyx552200
[有机交流] Suzuki偶联中硼酸酯原料掉下的硼酸酯是否可以自身偶联,生成B2Pin2杂质 500+3 天下大同 2026-07-06 3/150 2026-07-07 10:31 by martyr_oyx
[论文投稿] 刑法学论文投稿求助 +3 高源678 2026-07-04 5/250 2026-07-06 12:20 by 平天尺
[论文投稿] Journal of Environmental Chemical Engineering +3 雪花天涯 2026-07-05 3/150 2026-07-06 09:00 by bobvan
[基金申请] 最后一年,祈求好运 +3 pyw1992 2026-07-06 3/150 2026-07-06 08:53 by 苏_frp
信息提示
请填处理意见