计算化学公社

标题: 在VMD里面计算水合物中每个水分子的氢键的数量，数量完全不符合实际 [打印本页]

作者
Author: 牧生 时间: 2021-11-4 15:30
标题: 在VMD里面计算水合物中每个水分子的氢键的数量，数量完全不符合实际
本帖最后由牧生于 2021-11-4 22:07 编辑

我按照教程中，金-水界面的模拟，z方向上，单个水分子形成氢键数量，可以重复出来相同的结果。但我自己做了个水合物（用的四点水模型），结果算出来的氢键与实际不符合。
步骤如下；
①启动VMD，载入gro，删除仅有一帧，再载入xtc，播放完毕
②在命令窗口里粘贴
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]]
}

③再在命令窗口粘贴
set znow 5; set zend 50; set step 0.5 //z方向从5到50 埃
while {$znow<$zend} {
set ztmp [expr $znow+$step]
set val [numhbavg "same resid as resname SOL and z>=$znow and z<$ztmp and x>5 and x<15 and y>5 and y<15" 1500 2000] //x，y均在5~15埃之间，选择1500~2000帧
puts [format "%4.2f %5.3f" [expr ($znow+$ztmp)/2.0] $val]
set znow $ztmp
}

得到结果：

5.25 11.903
5.75 11.912
6.25 11.916
6.75 11.893
7.25 11.920
7.75 11.913
8.25 11.922
8.75 11.940
9.25 11.927
9.75 11.927
10.25 11.935
10.75 11.938
11.25 11.931
11.75 11.928
12.25 11.923
12.75 11.893
13.25 11.929
13.75 11.916
14.25 11.917
14.75 11.927
15.25 11.938
15.75 11.934
16.25 11.924
16.75 11.928

平均一个水分子有12个氢键。按理，四点水模型也不会造成这个影响，请问一下该如何处理

md.gro和md.xtc放在了蓝奏云
https://www.lanzouw.com/iXFo6w532cj

作者
Author: sobereva 时间: 2021-11-5 17:11
我没时间细看。你在VMD里用Hbond的drawing method，恰当设置好几何判断阈值后，仔细从图上看看显示出的氢键是否合乎自己的期望、虚拟点有没有产生奇怪的额外影响。

作者
Author: 牧生 时间: 2021-11-5 21:27
本帖最后由牧生于 2021-11-5 21:30 编辑

(, 下载次数 Times of downloads: 37) (, 下载次数 Times of downloads: 34)

看了一会，觉得问题出在四点水模型。每个氧原子形成了四个氢键，每个氢原子有两个氢键。所以一个水分子就有8个氢键了。

为了得到正确结果，脚本该如何修改呢。

作者
Author: xjw 时间: 2022-4-14 17:17

牧生发表于 2021-11-5 21:27
看了一会，觉得问题出在四点水模型。每个氧原子形成了四个氢键，每个氢原子有两个氢键。所以一 ...

请问，正常情况下，四点水每个水分应该有几个氢键啊

作者
Author: 牧生 时间: 2022-4-14 17:55
本帖最后由牧生于 2022-4-14 17:57 编辑

xjw 发表于 2022-4-14 17:17
请问，正常情况下，四点水每个水分应该有几个氢键啊

虚原子只是为了使得电荷分布更合理，所以它本身不形成氢键。所以，不应该是四点水显示的样子作为氢键的数量，应该是三点水为准。
一个水分子，会与周围其他水分子形成四个氢键。

作者
Author: xjw 时间: 2022-4-14 17:56
本帖最后由 xjw 于 2022-4-14 17:59 编辑

牧生发表于 2022-4-14 17:55
无论三点水，或者四点水，都是四个氢键，只是氢键的强弱不等而已。

哦哦，那请问怎么算每个水的氢键啊，您给的哪个算法可以嘛？

作者
Author: 牧生 时间: 2022-4-14 18:04
本帖最后由牧生于 2022-4-14 18:15 编辑

xjw 发表于 2022-4-14 17:56
哦哦，那请问怎么算每个水的氢键啊，您给的哪个算法可以嘛？

应该是以实际情况为准。因为虚原子是人为引进的，它并不存在的，所以应该去掉虚原子以后，计算三点水模型。但是氢键的数量，也不一定程序算出多少就是多少，比如常温常压的水，一个水分子形成3个氢键，但是以IGMH分析，可以看到四个归属于氢键的等值面，但其中有一个氢键的强度很弱，虽然其分布的距离和角度满足氢键的特征，但是强度不如一般我们理解上的氢键。
仔细看这个
http://sobereva.com/54

作者
Author: xjw 时间: 2022-4-14 18:10

牧生发表于 2022-4-14 18:04
应该是以实际情况为准。因为虚原子是人为引进的，它并不存在的，所以应该去掉虚原子以后，计算三点水模型 ...

好的，谢谢哈

作者
Author: xjw 时间: 2022-4-15 10:38

牧生发表于 2022-4-14 18:04
应该是以实际情况为准。因为虚原子是人为引进的，它并不存在的，所以应该去掉虚原子以后，计算三点水模型 ...

您好，我想问下如果我想计算水的氢键可以用tip3p这种三点水嘛？是不是后期算氢键也方便一点啊

作者
Author: Frozen-Penguin 时间: 2022-4-15 12:05

xjw 发表于 2022-4-15 10:38
您好，我想问下如果我想计算水的氢键可以用tip3p这种三点水嘛？是不是后期算氢键也方便一点啊

选择力场和水模型的主要依据是是否适用于研究的体系，与分析是否方便无关。
如果用四点水计算氢键存在问题，可以考虑删去虚拟原子后再计算。

作者
Author: xjw 时间: 2022-4-15 16:32

Frozen-Penguin 发表于 2022-4-15 12:05
选择力场和水模型的主要依据是是否适用于研究的体系，与分析是否方便无关。
如果用四点水计算氢键存在问 ...

嗯嗯，好的，谢谢啦

作者
Author: chen0201 时间: 2023-9-6 11:10
本帖最后由 chen0201 于 2023-9-6 11:14 编辑

请问楼主，我的输入文件是.xyz，加了晶格以后，我也是按照sob老师和您这个操作一样执行阿，为什么我的第二步会报错阿？

欢迎光临计算化学公社 (http://bbs.keinsci.com/)