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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 分子模拟 » Gromacs/Amber/NAMD » 【转载】计算不同z位置水能形成氢键数的VMD Tcl脚本
41/1返回列表
查看: 1960  |  回复: 3
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

childsliu

捐助贵宾 (正式写手)

二当家的

[交流] 【转载】计算不同z位置水能形成氢键数的VMD Tcl脚本 已有2人参与

From http://hi.baidu.com/sobereva/blo ... 595ef7fc037fbf.html
by Sobereva

有人请我写一个计算垂直于溶液界面的即不同z值处的水能形成氢键数的脚本，下面是我的VMD Tcl脚本，或许有人也用得着。
氢键判据用的是常用的35度3.5埃的几何判据，当然也可以直接在脚本里改。计算方法是，只要有水分子有一个原子在某一层里，
则这个水分子就认为属于这一层的水。对于每一帧，计算属于每一层的水selin与其它物质selbig之间的氢键数，氢键包括了这一层
中的水作为氢键受体和供体两种情况，其数目分别为代码中的变量a和b。并且加上这一层水内部之间的氢键数（变量c）的2倍。
a+b+2c除以这一层的水数，作为这一帧这一层的每个水的平均氢键数。脚本中循环轨迹中的每一帧，最终得到这一层平均氢键数。
nonum变量记录有多少帧在所设定的范围里没有水，这些帧不计算。#后面那行用于调试目的，要考察每帧结果就去掉开头的#。
首先运行下面的脚本，来加载实现这个功能的子程序
CODE:
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]]
set b [llength [lindex [measure hbonds 3.5 35 $selin $selbig] 0]]
set c [llength [lindex [measure hbonds 3.5 35 $selin] 0]]
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 不要恰顶着体系的边界，最好至少留出5埃的距离。选择范围中resname SOL代表了水。
CODE:
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

[ Last edited by lei0736 on 2009-11-25 at 11:23 ]
回复此楼

» 猜你喜欢

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

» 本主题相关价值贴推荐，对您同样有帮助:
下苦功，三个字，一个叫下，一个叫苦，一个叫功，一定要振作精神，下苦功。
1楼 2009-11-09 19:12:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

Wendee

至尊木虫 (著名写手)
★
小木虫(金币+0.5):给个红包，谢谢回帖交流
不错不错，感谢你的分享！
赞一下(1人) 回复此楼
2楼2011-03-17 22:39:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

长草团子

新虫 (初入文坛)
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
请问楼主，我试着运行了一下这个tcl脚本，也将其改成了我的参数，但是在运行的时候却出现了问题，在调用numhbavg的时候，显示错误：无效的命令名“”，请问这是为什么啊？
赞一下 回复此楼
3楼2018-03-30 09:11:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

长草团子

新虫 (初入文坛)
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
请问楼主，我试着运行了一下这个tcl脚本，也将其改成了我的参数，但是在运行的时候却出现了问题，在调用numhbavg的时候，显示错误：无效的命令名“”，请问这是为什么啊？
赞一下 回复此楼
4楼2018-03-30 14:37:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 childsliu 的主题更新
41/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 326求调剂 +4 诺贝尔化学奖觊� 2026-03-15 7/350 2026-03-16 17:11 by 诺贝尔化学奖觊�
[考研] 070303一志愿西北大学学硕310找调剂 +5 d如愿上岸 2026-03-12 8/400 2026-03-16 15:19 by peike
[考研] 一志愿华中师范071000，325求调剂 +6 RuitingC 2026-03-12 6/300 2026-03-16 14:50 by 可淡不可忘
[考研] 070303 总分349求调剂 +3 LJY9966 2026-03-15 5/250 2026-03-16 14:24 by xwxstudy
[硕博家园] 深圳大学硕士招生（2026秋，传感器方向，仅录取第一志愿） +4 xujiaoszu 2026-03-11 8/400 2026-03-16 09:45 by xujiaoszu
[考研] 070305求调剂 +3 mlpqaz03 2026-03-14 4/200 2026-03-15 11:04 by peike
[考研] 材料工程327求调剂 +3 xiaohe12w 2026-03-11 3/150 2026-03-14 20:20 by ms629
[考研] 307求调剂 +7 超级伊昂大王 2026-03-10 7/350 2026-03-14 00:49 by JourneyLucky
[考研] 材料与化工（0856）304求B区调剂 +6 邱gl 2026-03-12 7/350 2026-03-13 23:24 by 邱gl
[考研] 材料专硕288分求调剂 一志愿211 +4 在家想你 2026-03-11 4/200 2026-03-13 22:49 by JourneyLucky
[考研] 求材料调剂 085600英一数二总分302 前三科235 精通机器学习 一志愿哈工大 +4 林yaxin 2026-03-12 4/200 2026-03-13 22:04 by 星空星月
[考研] 304求调剂 +7 7712b 2026-03-13 7/350 2026-03-13 21:42 by peike
[考研] 【考研调剂求收留】 +3 Ceciilia 2026-03-11 3/150 2026-03-13 20:18 by JourneyLucky
[硕博家园] 085600 260分求调剂 +3 天空还下雨么 2026-03-13 5/250 2026-03-13 18:46 by 天空还下雨么
[考研] 求b区学校调剂 +3 周56 2026-03-11 3/150 2026-03-13 16:20 by JourneyLucky
[考研] 求调剂 +7 18880831720 2026-03-11 7/350 2026-03-13 16:10 by JourneyLucky
[考研] 工科调剂 +4 Jiang191123！ 2026-03-11 4/200 2026-03-13 15:15 by Miko19
[考研] 08食品或轻工求调剂，本科发表3篇sci一区top论文，一志愿南师大食品科学与工程 +3 我是一个兵， 2026-03-10 3/150 2026-03-13 10:21 by Yuyi.
[考研] 279求调剂 +3 莫xiao 2026-03-10 4/200 2026-03-11 08:06 by 斩魂滴兔子！
[考研] 收调剂 +7 调剂的考研学生 2026-03-10 7/350 2026-03-10 17:57 by 麦茶汤圆
信息提示
关闭
请填处理意见
关闭