12.3分析脚本

在蛋白质附近找到水

他的例子为每一帧轨迹找到蛋白质附近的水，并写出包含这些水的PDB文件：

set sel [atomselect top "water and same residue as (within 2 of protein)"]

set n [molinfo top get numframes]

for { set i 0 } { $i < $n } { incr i } {

$sel frame $i

$sel update

$sel writepdb water_$i.pdb

}

frame选项设置选择的帧，update告诉原子选择重新计算哪些水靠近蛋白质，writepdb将选择的水写入文件。

一个选择的总质量

proc total_mass {selection} {

set sum 0

foreach mass [$selection get mass] {

set sum [expr {$sum + $mass}]

}

return $sum

}

注意上面示例中的expr命令后面的花括号。省略这些大括号会导致脚本运行速度慢三倍！这个故事的寓意是：总是给传递给expr的表达式加上大括号。

下面是做同样事情的另一种方法（稍微慢一点）。这是有效的，因为从选择返回的质量是一个列表的列表。将其放在eval的引号内使其成为一个向量序列，因此vecadd命令将对其起作用。

proc total_mass1 {selection} {

set mass [$selection get mass]

eval "vecadd $mass"

}

坐标最小和最大

找到给定分子在x， y和z方向上的最小和最大坐标值（参见测量命令‘ minmax ’）。函数取分子id并返回两个向量；第一个包含最小值，第二个包含最大值。

proc minmax {molid} {

set sel [atomselect top all]

set sx [$sel get x]

set sy [$sel get y]

set sz [$sel get z]

set minx [lindex $sx 0]

set miny [lindex $sy 0]

set minz [lindex $sz 0]

set maxx $minx

set maxy $miny

set maxz $minz

foreach x $sx y $sy z $sz {

if {$x < $minx} {set minx $x} else {if {$x > $maxx} {set maxx $x}}

if {$y < $miny} {set miny $y} else {if {$y > $maxy} {set maxy $y}}

if {$z < $minz} {set minz $z} else {if {$z > $maxz} {set maxz $z}}

}

return [list [list $minx $miny $minz] [list $maxx $maxy $maxz]]

}

回转半径

计算一个选择的旋转半径（参见测量rgyr）。旋转半径的平方定义为pi mi(~ ri−~ rc)2/ pi mi。这使用了本章前面定义的质心函数；一个更快的版本将用测量中心取代它。注意，measure rgyr命令所做的事情与这个脚本相同，只是速度快得多。

proc gyr_radius {sel} {

# make sure this is a proper selection and has atoms

if {[$sel num] <= 0} {

error "gyr_radius: must have at least one atom in selection"

}

# gyration is sqrt( sum((r(i) - r(center_of_mass))^2) / N)

set com [center_of_mass $sel]

set sum 0

foreach coord [$sel get {x y z}] {

set sum [vecadd $sum [veclength2 [vecsub $coord $com]]]

}

return [expr sqrt($sum / ([$sel num] + 0.0))]

}

把这个应用到alanin.pdb 坐标文件

vmd > mol new alanin.pdb

vmd > set sel [atomselect top all]

vmd > gyr_radius $sel

Info) 5.45443

均方根偏差RMSD

计算轨迹两帧间选择的均方根差。这需要一个选择和两个帧的值进行比较。

proc frame_rmsd {selection frame1 frame2} {

set mol [$selection molindex]

# check the range

set num [molinfo $mol get numframes]

if {$frame1 < 0 || $frame1 >= $num || $frame2 < 0 || $frame2 >= $num} {

error "frame_rmsd: frame number out of range"

}

# get the first coordinate set

set sel1 [atomselect $mol [$selection text] frame $frame1]

set coords1 [$sel1 get {x y z}]

# get the second coordinate set

set sel2 [atomselect $mol [$selection text] frame $frame2]

set coords2 [$sel2 get {x y z}]

# and compute the rmsd values

set rmsd 0

foreach coord1 $coords1 coord2 $coords2 {

set rmsd [expr $rmsd + [veclength2 [vecsub $coord2 $coord1]]]

}

# divide by the number of atoms and return the result

return [expr $rmsd / ([$selection num] + 0.0)]

}

下面使用帧rmsd函数列出分子在整个轨迹上的rmsd，与第一个帧相比。

vmd > mol new alanin.psf

vmd > mol addfile alanin.dcd

vmd > set sel [atomselect top all]

vmd > for {set i 0} {$i < [molinfo top get numframes]} {incr i} {

? puts [list $i [frame_rmsd $sel $i 0]]

? }

0 0.0

1 0.100078

2 0.291405

3 0.523673

....

97 20.0095

98 21.0495

99 21.5747

最后一个示例展示了如何设置beta字段。这很有用，因为其中一种上色方法是“Beta”，它使用Beta值根据当前的色阶给分子上色。（这也可以用占用字段来完成。）因此，重新定义beta值允许您根据自己的定义为分子上色。一个有用的例子是根据到特定点的距离给分子上色（在这种情况下，根据到病毒中心的距离（0,0,0）给脊髓灰质炎病毒原聚体上色有助于显示表面特征）。

proc betacolor_distance {sel point} {

# get the coordinates

foreach coord [$sel get {x y z}] {

# get the distance and put it in the "newbeta" list

set dist [veclength2 [vecsub $coord $point]]

lappend newbeta $dist

}

# set the beta term

$sel set beta $newbeta

}

这里有一种用法：

# load pdb2plv.ent using anonymous ftp to the PDB

vmd > mol new 2plv

vmd > set sel [atomselect top all]

vmd > betacolor_distance $sel {0 0 0}

然后转到图形菜单，将“着色方法”设置为“Beta”。