标题: 一个每一帧统计水分子数的程序，希望懂得人帮个忙 [打印本页] 作者Author: ruanyang 时间: 2014-11-19 13:52 标题: 一个每一帧统计水分子数的程序，希望懂得人帮个忙 程序： set time [clock format [clock seconds] -format {%b. %d, %Y %I:%M:%S %p}]
#selects the top open molecule's carbon
set nt [atomselect top "resname Gr1"]
# returns z min of the cnt
set NT_min [lindex [lindex [measure minmax $nt] 0] 2]
# returns z max of the cnt
set NT_max [lindex [lindex [measure minmax $nt] 1] 2]
set x_0 [lindex [measure center $nt] 0]
# returns y coord
set y_0 [lindex [measure center $nt] 1]
set NT_radius 6.780
set num_frames [molinfo top get numframes]
set x_list {}
set y_list {}
set wat_sel [ atomselect top resname SOL ]
set out [open "Gr1_y_list.txt" w]
puts $out "This file contains the output of the number of water molecules in the Gr1"
puts $out "by measuring the number of oxygen in the Gr1"
puts $out""
puts $out"This script was created by;"
puts $out "Christopher Stiles:cs145331@albany.edu or chris@cs86.com"
puts $out ""
puts $out "There time of creation was:"
puts $out "$time"
puts $out "Number of frames uesd:$num_frames"
puts $out "User defined radius of the Gr1:$NT_radius"
puts $out "Box dimensions:($Lx, $Ly,$Lz)"
puts $out ""
close out
set out [open "Gr1_y_list.txt" a]
for {set i 0} {$i < $num_frames} {incr i} {
$wat_sel frame $i
$wat_sel num
$wat_sel update
#$wat_sel set user 1.0
lappend x_list $i
lappend y_list [$wat_sel num]
puts $out "[$wat_sel num]"
}
close $out
Gr1:代表的是中间的那根碳管。谢谢！
作者Author: sobereva 时间: 2014-11-19 14:14
没明白你的意图。
贴代码的时候建议选择“代码”按钮然后再贴，这样不会被转义成表情，代码看着也整齐。作者Author: sobereva 时间: 2014-11-19 14:17
几行代码就能完成的事，搞得太复杂了。比如下面的代码显示每一帧resid 15的5埃内的水的数量，你把选择范围改成纳米管的区域即可（如waters and x> 10 and x<20）
for {set i 0} {$i<200} {incr i 1} {
atomselect top {same residue as{waters within 5 of resid 15}} frame $i
puts "Frame: $i number of waters: [expr [atomselect$i num]/3]"