24小时热门版块排行榜    

查看: 1700  |  回复: 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 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 材料专硕306英一数二 +7 z1z2z3879 2026-03-16 9/450 2026-03-17 17:31 by ccjequ
[硕博家园] 湖北工业大学 生命科学与健康学院-课题组招收2026级食品/生物方向硕士 +3 1喜春8 2026-03-17 5/250 2026-03-17 17:18 by ber川cool子
[考研] 化学工程321分求调剂 +11 大米饭! 2026-03-15 14/700 2026-03-17 17:11 by ruiyingmiao
[考研] 290求调剂 +3 p asserby. 2026-03-15 4/200 2026-03-17 16:35 by wangkm
[考研] 材料与化工304求B区调剂 +7 邱gl 2026-03-11 8/400 2026-03-17 09:36 by 努力学习赚彩礼
[文学芳草园] 伙伴们,祝我生日快乐吧 +17 myrtle 2026-03-10 26/1300 2026-03-16 18:32 by 青橙Ln
[考研] 333求调剂 +3 文思客 2026-03-16 7/350 2026-03-16 18:21 by 文思客
[考研] 304求调剂 +3 曼殊2266 2026-03-14 3/150 2026-03-16 16:39 by houyaoxu
[考研] 283求调剂 +10 小楼。 2026-03-12 14/700 2026-03-16 16:08 by 13811244083
[教师之家] 焦虑 +7 水冰月月野兔 2026-03-13 9/450 2026-03-16 10:00 by Quakerbird
[考研] 一志愿哈工大材料324分求调剂 +5 闫旭东 2026-03-14 5/250 2026-03-14 14:53 by 木瓜膏
[考研] 招收0805(材料)调剂 +3 18595523086 2026-03-13 3/150 2026-03-14 00:33 by 123%、
[考研] 337一志愿华南理工0805材料求调剂 +7 mysdl 2026-03-11 9/450 2026-03-13 22:43 by JourneyLucky
[考研] 材料工程调剂 +9 咪咪空空 2026-03-12 9/450 2026-03-13 22:05 by 星空星月
[考研] 315求调剂 +9 小羊小羊_ 2026-03-11 10/500 2026-03-13 21:13 by SXNU李老师
[考研] 材料与化工085600调剂求老师收留 +9 jiaanl 2026-03-11 9/450 2026-03-13 20:22 by JourneyLucky
[考研] 材料工程调剂 +4 咪咪空空 2026-03-11 4/200 2026-03-13 19:57 by JourneyLucky
[考研] 328化工专硕求调剂 +4 。,。,。,。i 2026-03-12 4/200 2026-03-13 14:44 by JourneyLucky
[考研] 308求调剂 +3 是Lupa啊 2026-03-12 3/150 2026-03-13 14:30 by 求调剂zz
[考研] 大连大学化学专业研究生调剂 +3 琪久. 2026-03-10 8/400 2026-03-11 10:02 by 琪久.
信息提示
请填处理意见