计算化学公社

标题: 使用tcl脚本计算氢键数结果有误 [打印本页]

作者
Author: zhangfaxue 时间: 2024-12-7 13:58
标题: 使用tcl脚本计算氢键数结果有误
脚本如下，想要计算两个组分氢键数目随着时间的变化，之前测试了表面活性剂和水分子之间氢键数目发现计算出来结果全是3，为了验证脚本是否存在问题就修改了两个组分都为水分子，结果计算还是氢键数目全部是3，请教各位哪里存在问题？

pbc wrap -center com -centersel "resname WATE" -all

proc hb {fps1 fps2 dt interval} {
set cum 0
set outputFile "hb-WATE-vmd.txt"
set outFile [open $outputFile "w"]
set targetTime 0

for {set i $fps1} {$i <= $fps2} {incr i 1} {
set currentTime [expr $i * $dt]

set dd [atomselect top "resname WATE" frame $i]
set aa [atomselect top "resname WATE" frame $i]

set temp [measure hbonds 3.5 35 $dd $aa]

set hbondsCount [llength $temp]

if {($currentTime - $targetTime) >= $interval} {
set targetTime [expr $targetTime + $interval]

puts $outFile "$currentTime $hbondsCount"

}
}

close $outFile
puts "Output written to $outputFile"

}
(, 下载次数 Times of downloads: 8)

作者
Author: Uus/pMeC6H4-/キ 时间: 2024-12-7 14:17
measure hbonds返回的是由给体原子index、受体原子index、氢原子index共3个元素数量相同的子列表构成的列表，而llength直接返回的是列表中元素数量而不关心元素是不是还有列表，所以直接把measure hbonds返回结果存到变量$temp里再用llength $temp可不就取成3了。统计氢键数量应该用其中一个子列表来数，比如用第一个子列表的命令写作[llength [lindex $temp 0]]。

有个小小的请求，如果像上帖这样提问得到回答，可否提供个反馈呢？回个帖或者评个分都行，无论回答是很有帮助还是完全不对，提问者的回应可以给回答者一个确认，也可以帮助有相同问题的后来者解惑……

作者
Author: fhh2626 时间: 2024-12-7 20:39
VMD的extension里面有统计氢键的插件

作者
Author: zhangfaxue 时间: 2024-12-8 09:49

Uus/pMeC6H4-/キ发表于 2024-12-7 14:17
measure hbonds返回的是由给体原子index、受体原子index、氢原子index共3个元素数量相同的子列表构成的列表 ...

谢谢回复！问题已解决。回复是最基本的，很抱歉上一条我疏忽忘记了没有反馈。非常感谢您的帮助。

作者
Author: zhangfaxue 时间: 2024-12-8 09:49

fhh2626 发表于 2024-12-7 20:39
VMD的extension里面有统计氢键的插件

好的好的谢谢

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