计算纯水体系下的氢键数，结果不符合实际

airbang

老师您好，我按照您计算不同z位置水能氢键数的程序，建立了一个纯体相水（5nm*5nm*5nm）的体系，用NAMD进行了模拟5000步  模型图、运行文件和参数文件都放在了附件
步骤如下：按照您的要求我计算时我在xy轴空出了5埃 氢键判据用的是常用的30度3.5埃的几何判据  计算出的结果不符合实际
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 30 $selbig $selin] 0]]
set b [llength [lindex [measure hbonds 3.5 30 $selin $selbig] 0]]
set c [llength [lindex [measure hbonds 3.5 30 $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]]
}
for {set i 10} {$i<=30} {incr i} {
set k [expr $i*0.1]
set now [numhbavg "same resid as water  and x<45 and x>-45 and y<45 and y>-45 and z<[expr $k+0.1] and z>=$k" 4900 5000]
puts [format "%4.2f %4.2f %5.3f" $k [expr $k+0.1] $now]
}
结果如下：1.00 1.10 2.476
1.10 1.20 2.276
1.20 1.30 2.130
1.30 1.40 2.423
1.40 1.50 2.182
1.50 1.60 2.381
1.60 1.70 2.467
1.70 1.80 2.320
1.80 1.90 2.368
1.90 2.00 2.214
2.00 2.10 2.345
2.10 2.20 2.692
2.20 2.30 2.267
2.30 2.40 2.526
2.40 2.50 2.333
2.50 2.60 2.533
2.60 2.70 2.409
2.70 2.80 2.389
2.80 2.90 2.333
2.90 3.00 2.375
3.00 3.10 2.391

sobereva

注意VMD里氢键键角的定义和习俗不同，一般说的30度氢键键角判据对应VMD的氢键键角判据大约43度

另外，VMD里氢键计算命令没法考虑周期边界条件，会一定程度低估氢键数

airbang

sobereva 发表于 2022-5-19 22:10
注意VMD里氢键键角的定义和习俗不同，一般说的30度氢键键角判据对应VMD的氢键键角判据大约43度

另外，VM ...

好的 谢谢sob老师