VMD实现动态选择蛋白附近的磷脂分子计算序参量的tcl代码实现问题 - 计算化学公社

proc orderparam-c3 { result { seltext all } } {
upvar $result arr
#获得要计算的帧数
set n [molinfo top get numframes]
#以POPC为例。计算18个C原子
for { set i 2 } { $i <= 20 } { incr i } {
puts $i
set sum 0.0
set nh 0
#每个C原子的所有帧循环
for { set frame 0 } { $frame < $n } { incr frame } {
#选择距离蛋白5A 以内的磷脂
set dynamic_lipid [atomselect top "($seltext) and lipid and same resid as within 5 of protein"]
$dynamic_lipid frame $frame
$dynamic_lipid update
set cp [atomselect top "($dynamic_lipid) and name C3$i"]
if { [$cp num] == 0 } {
puts "skipping $i"
continue
}
#选择与C相连的HX、HY
set hx [atomselect top "($dynamic_lipid) and name H${i}X"]
set hy [atomselect top "($dynamic_lipid) and name H${i}Y"]
#set hz [atomselect top "($dynamic_lipid) and name H${i}Z"]
$cp frame $frame
$hx frame $frame
$hy frame $frame
$cp update
$hx update
$hy update
#有num个磷脂分子被选中
set nres [$cp num]
#获得C原子坐标
set cpx [$cp get x]
set cpy [$cp get y]
set cpz [$cp get z]
#C-H坐标差
set hxx [vecsub $cpx [$hx get x]]
set hxy [vecsub $cpy [$hx get y]]
set hxz [vecsub $cpz [$hx get z]]
#计算与Z轴的夹角余弦平方
foreach dx $hxx dy $hxy dz $hxz {
set norm2 [expr {$dx*$dx + $dy*$dy + $dz*$dz}]
set sum [expr {$sum + $dz*$dz/$norm2}]
}
incr nh $nres
#如果C上有HY
if { [$hy num] != 0 } {
set hyx [vecsub $cpx [$hy get x]]
set hyy [vecsub $cpy [$hy get y]]
set hyz [vecsub $cpz [$hy get z]]
foreach dx $hyx dy $hyy dz $hyz {
set norm2 [expr {$dx*$dx + $dy*$dy + $dz*$dz}]
set sum [expr {$sum + $dz*$dz/$norm2}]
}
incr nh $nres
}
}
#序参量计算
set arr($i) [expr {-1.5*$sum/$nh + 0.5}]
}
}

复制代码