12.2使用atomselect命令

原子选择是获取分子中原子信息的主要方法。它分为两个步骤。第一步是在给定选择文本、分子id和可选框架号的情况下创建一个选择。这是由一个名为atomselect的函数完成的，该函数返回新原子选择的名称。第二步是使用创建的选择来访问关于选择中的原子的信息。

Atom选择是作为Tcl函数实现的。从atomselect返回的数据是要使用的函数的名称。名称的形式为atomselect%d，其中“%d”是非负数（例如“atomselect0”，“atomselect26”，…）。

使用atomselect命令创建的函数的方法是将名称存储到一个变量中，然后在需要时使用该变量获取名称。

vmd> set sel [atomselect top "water"]

atomselect3

vmd> $sel text

water

这相当于说

vmd> atomselect3 text

考虑这个问题的最简单方法是atomselect命令创建一个对象。

要从对象获取信息，必须向它发送命令。因此，在上面的示例（atomselect1 num）中，对象“atomselect1”被发送命令“num”，该命令要求对象返回选择中的原子数。这些派生的对象函数（名字类似atomselect3的）有很多选项，如9.3.2节所述，例如给定选择

vmd> set sel [atomselect top "resid 4"]

atomselect4

可以使用命令获取选择中每个原子的原子名

vmd> $sel get name

N H CA CB C O

记住，这和vmd> atomselect4 get name一样

通过提交一个列表可以请求多个属性，所以如果你想知道哪些原子在主链上，

vmd> $sel get {name backbone}

{N 1} {H 0} {CA 1} {CB 0} {C 1} {O 1}

原子的坐标是

vmd> $sel get {x y z}

{0.710000 4.211000 1.093000} {-0.026000 3.700000 0.697000} {0.541000

4.841000 2.388000} {-0.809000 4.462000 2.976000} {1.591000 4.371000

3.381000} {2.212000 5.167000 4.085000}

注意，从get命令返回的数据格式取决于请求的属性数量。如果您只请求一个属性，就像上面的get name示例一样，那么您将返回一个简单的元素列表。另一方面，如果您请求两个或更多属性，您将返回一个子列表。具体来说，它是一个大小为n的列表，其中每个元素本身是一个大小为i的列表，其中n是选择中的原子数，i是请求的属性数。

如果一次只检索一个属性，您的脚本将运行得更快，因为VMD不必为每个属性构造子列表。记住，在Tcl中，你可以使用foreach命令一次遍历多个列表：

foreach resid [$sel get resid] resname [$sel get resname] {

# process each resid and resname here

}

一个你可以用坐标快速构建的函数是计算几何中心的方法(不完全是质量中心；这有点难)。这也使用了一些在关于向量和矩阵的章节中讨论的向量命令，但你应该能够从上下文中找出它们。

proc geom_center {selection} {

# set the geometrical center to 0

set gc [veczero]

# [$selection get {x y z}] returns a list of {x y z}

# values (one per atoms) so get each term one by one

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

# sum up the coordinates

set gc [vecadd $gc $coord]

}

# and scale by the inverse of the number of atoms

return [vecscale [expr 1.0 /[$selection num]] $gc]

}

有了这个定义，您可以这样说（假设使用前面的atomselection示例创建了$sel）

vmd> geom_center $sel

0.703168 4.45868 2.43667

我将逐行浏览这个示例。该函数命名为geom center，并接受一个参数，即选择的名称。第一行将变量“gc”设置为零向量，即0 0 0。在第二行代码中，发生了两件事。首先，命令

$selection get {x y z}

执行，字符串被替换为结果，即

{0.710000 4.211000 1.093000} {-0.026000 3.700000 0.697000} {0.541000

4.841000 2.388000} {-0.809000 4.462000 2.976000} {1.591000 4.371000

3.381000} {2.212000 5.167000 4.085000}

这是一个包含6个项的列表（每个项代表选择中的每个原子），每个项是一个包含三个元素的列表，即x、y和z坐标，按此顺序。

“foreach”命令将列表拆分为六个项，并逐项向下查看列表，将变量“coord”设置为每个连续的项。在循环中，$coord的值被加到total sum中。

最后一行返回所选原子的几何中心。由于几何中心被定义为坐标向量的总和除以元素的数量，到目前为止我只计算了向量的总和，我需要元素数量的倒数，这是用表达式完成的

expr 1.0 / [$selection num]

“1.0”中的小数很重要，否则Tcl会进行整数除法。最后，这个值用于缩放坐标向量的总和（使用vecscale），它返回新值，该值本身作为过程的结果返回。质心函数稍微难一点因为你需要得到质量以及x， y， z的值，然后把它们分解成分量。质心的公式是

proc center_of_mass {selection} {

# some error checking

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

error "center_of_mass: needs a selection with atoms"

}

# set the center of mass to 0

set com [veczero]

# set the total mass to 0

set mass 0

# [$selection get {x y z}] returns the coordinates {x y z}

# [$selection get {mass}] returns the masses

# so the following says "for each pair of {coordinates} and masses,

# do the computation ..."

foreach coord [$selection get {x y z}] m [$selection get mass] {

# sum of the masses

set mass [expr $mass + $m]

# sum up the product of mass and coordinate

set com [vecadd $com [vecscale $m $coord]]

}

# and scale by the inverse of the number of atoms

if {$mass == 0} {

error "center_of_mass: total mass is zero"

}

# The "1.0" can’t be "1", since otherwise integer division is done

return [vecscale [expr 1.0/$mass] $com]

}

vmd> center_of_mass $sel

Info) 0.912778 4.61792 2.78021

“get”的反义词是“set”。许多关键字（最明显的是“x”、“y”和“z”）可以设置为新值。例如，这允许更改原子坐标、更新占用值或添加用户力。您还可以更改resname、segid等，这在VMD中可能比手工编辑PDB文件更容易。

set sel [atomselect top "index 5"]

$sel get {x y z}

{1.450000 0.000000 0.000000}

$set set {x y z} {{1.6 0 0}}

注意，就像get选项返回一个列表列表一样，set选项也需要一个列表列表，这就是为什么需要额外的花括号。同样，这必须是一个大小为n的列表，包含大小为i的元素，例外情况是，如果n = 1，列表被复制了足够多的次数，所以每个原子有一个元素。

# get two atoms and set their coordinates

set sel [atomselect top "index 6 7"]

$sel set {x y z} { {5 0 0} {7.6 5.4 3.2} }

在本例中，索引为6的原子将其（x, y, z）值设置为5 0 0，索引为7的原子将其坐标更改为7.6 5.4 3.2。

通过获取坐标、改变它们（比如添加一些偏移量）和替换它，可以以这种方式移动原子。下面是一个这样做的函数：

proc moveby {sel offset} {

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

lappend newcoords [vecadd $coord $offset]

}

$sel set {x y z} $newcoords

}

要使用这个函数(在本例中，对选择“$movesel”应用（x y z) =(0.1 -2.8 9)的偏移量）：

moveby $movesel {0.1 -2.8 9}

但是，为了简化问题，在选项中添加了一些选项来处理移动（这些命令也在c++中实现，并且比Tcl版本快得多）。这些函数是moveby、moveto和move。前两个是位置向量，最后一个是变换矩阵。

第一个命令moveby将选择中的每个原子按给定的矢量偏移量移动。

$sel moveby {1 -1 3.4}

第二个是moveto，它将选择中的所有原子移动到给定的坐标（对于一个以上原子的选择使用这个方法会很奇怪，但这是允许的）。例子:

$sel moveto {-1 1 4.3}

最后一个是move，它将给定的变换矩阵应用于每个原子坐标。它最适合用于围绕给定轴旋转一组原子，如

$sel move [trans x 90]

它将所选内容围绕x轴旋转90度。当然，可以使用任何变换矩阵。

下面是一个更有用的例子，它将CA-CB键周围的侧链原子旋转10度。

# get the sidechain atoms (CB and onwards)

set sidechain [atomselect top "sidechain residue 22"]

# get the CA coordinates -- could do next two on one line ...

set CA [atomselect top "name CA and residue 22"]

set CAcoord [lindex [$CA get {x y z}] 0]

# and get the CB coordinates

set CB [atomselect top "name CB and residue 22"]

set CBcoord [lindex [$CB get {x y z}] 0]

# apply a transform of 10 degrees about the given bond axis

$sidechain move [trans bond $CAcoord $CBcoord 10 deg]