计算化学公社

标题: VMD里原子选择语句的语法和例子 [打印本页]

作者
Author: sobereva 时间: 2019-8-19 22:02
标题: VMD里原子选择语句的语法和例子

VMD里原子选择语句的语法和例子

Syntax and examples of atomic selection in VMD

文/Sobereva@北京科音 2019-Aug-19

1 前言

VMD（http://www.ks.uiuc.edu/Research/vmd/）是极其强大、灵活的化学体系可视化程序，笔者之前也写过不少相关文章，见http://sobereva.com/category/VMD/。VMD的选择语句(selection)用来选择满足特定要求的体系中的原子，其用法极度简单灵活，对于VMD的使用至关重要。

选择语句在VMD里用的地方非常多，无处不在。比如如果想在图形窗口只显示指定的区域，那么可以在Graphics - Representation里在selected atoms的地方写上选择语句。如果想把指定的区域保存成新文件，那么可以在File - Save coordinates里在Selected atoms写上选择语句。很多自带的插件也需要选择语句，比如VMD自带的径向分布函数计算插件，里面selection 1、selection 2就是让你输入选择语句的地方。选择语句在VMD的分析脚本的编写中起到极为关键性角色，用于创建atomselect对象、作为自带命令传入的参数。我在“北京科音分子动力学与GROMACS培训班”（http://www.keinsci.com/workshop/KGMX_content.html）里非常深入系统讲解VMD分析脚本的编写并给出巨量例子，其中大量用到选择语句的知识。顺带一提，如今GROMACS里也可以用selections语句，和VMD很大程度一致，但不完全一致，在培训里我也有专门讲解。

鉴于经常有人问VMD的选择语句怎么用，每次回复很麻烦，笔者遂专门写个小文说一下。本文内容对应VMD 1.9.3。下文从简单到复杂进行讲解。

2 单关键词（Single words）

有一些关键词可以直接选择特定原子，以下举例一部分：
all：所有原子
none：不选择任何原子
noh：氢以外的原子（即重原子）
ion：离子
water：水
backbone：生物大分子骨架
sidechain：生物大分子侧链
protein：蛋白
nucleic：核酸
helix：螺旋
alpha_helix：alpha螺旋（是helix中的子集，较长一段螺旋才算）
sheet：折叠
turn：转角
coil：盘绕
alpha：蛋白质的alpha碳
acidic：PH=7时带负电氨基酸
basic：PH=7时带正电氨基酸
charged：acidic和basic的并集
neutral：电中性氨基酸
polar：极性残基
hydrophobic：疏水性残基
bonded：成键的原子
hetero：非蛋白质和核酸的部分
carbon、hydrogen、oxygen、nitrogen、sulfur：相应元素。对于其它元素没法这么输入元素名来选择

这些单关键词实际上可以在Representation界面里的Selections标签页里的Singlewords直接看到，可见可以用的单关键词远不止上述这些。有些单关键词其实是复合选择语句，比如你选中hetero，就会看到其定义其实是not (protein or nucleic)。

(, 下载次数 Times of downloads: 199)

注意有些情况下，单关键词未必能如实选择相应的区域。比如你载入的结构里有水，如果输入文件里水的残基名很特殊，比如叫FFF，那么VMD就不会把这个残基识别成水分子，用water关键词的时候也没法选中这些水。

3 一般关键词

用下面这些关键词可以通过属性选取原子，都是后面要接参数的
name：原子名。例：name OW选择原子名叫OW的原子
index：原子序号（从0开始！）。例：index 4
serial：原子序号（从1开始）
type：原子类型。例：type CA选择CA原则类型
element：元素名。例：element P选择磷原子
resname：残基名。例：resname ALA代表选择丙氨酸
residue：残基编号，从0开始。例resid 372代表选择372号残基
resid：残基编号，从1开始。若结构文件里有残基号则与之一致
chain：链名。例：chain B代表选择B链
fragment：片段编号。VMD对每个键连的片段自动设定一个编号。例：fragment 4代表选择片段4
numbonds：成键数目。例：numbonds=2或numbonds 2代表选形成了两个键的原子
structure：二级结构。例：structure H代表选择螺旋(helix)区域
x,y,z：X/Y/Z笛卡尔坐标
vx,vy,vz：X/Y/Z方向速度
beta：pdb文件中的beta值
occupancy：pdb文件中的原子占有率
mass：原子质量
charge：原子电荷
phi、psi：蛋白质骨架角度
radius：原子半径
...等等

每个属性后面能接什么值，在Selections标签页里都能看到，不确定的话看一眼便知：

(, 下载次数 Times of downloads: 194)

许多属性并非对于任何输入文件都能用。比如：
·使用charge属性，必须输入的文件里体现了原子电荷才行，比如可以用mol2或pqr，后者详见《使用Multiwfn+VMD以原子着色方式表现原子电荷、自旋布居、电荷转移、简缩福井函数》（http://sobereva.com/425）。
·使用beta属性，通常需要用pdb文件作为输入，因为里面专门有一列记录B因子信息。
·用type的话必须载入拓扑文件才行。
·用vx、vy、vz的话，对于GROMACS用户，参看《使VMD读入Gromacs产生的trr轨迹中速度信息的方法》（http://sobereva.com/117）。
·element信息是很多文件里没有的，比如GROMACS的.gro文件里就没体现

4 选择语句中可利用的规则

在选择语句中有以下规则可以利用，通过组合、嵌套，使得选择语句无比强大
·可以写多个参数一次选择一批，彼此间用空格分隔
·可以用... to ...选择特定范围
·可以用与、或、非这些逻辑关系：and、or、not
·可以用( )或{ }指定语句处理的优先顺序
·双引号内的字符会被视为整体，并且可以使用正则表达式
·用单引号扩住则里面的字符可以避免被转义
·可以用判断语句：<, <=, =, >=, >, !=
·可以用函数：sqr（平方）, sqrt（开根号）, abs（绝对值）, sin, cos, tan, atan, asin, acos, sinh, cosh, tanh, exp, log, log10
·支持运算符：+ - * /。可以用^或**来表示多少次方
·特殊选择方式：
within 5 of AAA：距离AAA 5埃以内的原子。选取时不考虑周期边界条件，用pbwithin则考虑
exwithin 5 of AAA ：同上，但不包含AAA自身
withinbonds 2 of AAA：距离AAA不超过两个键的原子
same p as AAA：与AAA选区的p属性相同的部分
ringsize 5 from AAA：处于AAA中五元环上的原子
maxringsize 6 from AAA：处于AAA中<=六元环的原子

下面来看一些具体例子
index 5 to 200 210：序号在5~200内的原子以及210号原子
protein or nucleic：蛋白质与核酸的原子
resname ALA CYS ARG：丙氨酸、半胱氨酸、精氨酸原子
backbone not helix：除了螺旋区域以外的骨架原子
name CA CB 或 name "CA|CB" 或 name "C[AB]" 或 name "C(A|B)"：名为CA和CB的原子
name "C."：名字为两个字符且第一个字符为C的原子
name "CE[1-3]"：名字为CE1、CE2、CE3的原子
name 'O5*'：叫O5*的原子。注意原子名带星号的在选取时要用单引号括住以免转义
resname 'CA2+'：残基名是CA2+的原子（二价钙离子）。名字带正负号的也要用单引号括住以免转义
mass > 5：质量大于5的原子
abs(charge)>1：电荷大小超过1的原子
x<6 and x>3：选择x在3~6埃区域内的一层原子
x>1 and x<8 and y>24 and y<35 and z>1 and z<5：一个矩形区域内的原子
sqr(x-5)+sqr(y+4)+sqr(z) < sqr(5) ：以(5,-4,0)点为中心半径5埃以内的原子
((x-33)^2+(y-14.5)^2)<12^2 and z<40 and z>10：选择以x=33、y=14.5埃为中心，半径为12埃，z范围在10~40埃的柱形区域
x+y+z<80：斜切面内侧的原子（回忆平面方程）
not {oxygen and numbonds=0}：扣除孤立的氧原子（可以用于去除X光衍射pdb文件里的结晶水）
within 6 of protein：距离蛋白质6埃以内的原子
not within 5 of resname ADP：距离名为ADP的分子5埃以外的原子
water within 5 of residue 8 to 44：距离8~44号残基5埃以内的水
withinbonds 2 of index 31：距离编号为31原子的两个键及以内的原子
maxringsize 6 from protein：蛋白当中所有六元及六元以下环上的原子
same resname as resid 33：所有与33号残基相同名称的残基
same residue as {protein within 5 of nucleic}：与核酸的原子相距5埃以内的蛋白的原子，并且把被截断的残基保留完整
x > 15 and not same fragment as {exwithin 8 of protein}：蛋白质以及蛋白质8埃范围以外的原子，保留完整片段，同时x坐标得大于15埃

以上例子中，涉及到坐标、速度变量的，属于动态选区，即随着帧号变化被选择的原子会可能发生变化。观看这些选区的时候，注意在Representation界面的Trajectory标签页里要把Update Selection Every Frame选上，否则选中的原子是对刚选中时那一帧而言的，不会随着轨迹播放被动态更新。在一些VMD的插件中，比如计算rdf的Radial Pair Distribution Function g(r)插件里，当Selection文本框里用了动态选区时，应当把Update Selections复选框选上，否则也由于不会被动态更新而和期望的不符。

作者
Author: 函数与激情 时间: 2021-9-2 19:52
请问如何输出选中原子的index呢

作者
Author: sobereva 时间: 2021-9-2 21:23

函数与激情发表于 2021-9-2 19:52
请问如何输出选中原子的index呢

[atomselect top "选择语句"] list

作者
Author: 函数与激情 时间: 2021-9-2 21:49

sobereva 发表于 2021-9-2 21:23
[atomselect top "选择语句"] list

太感谢老师啦

作者
Author: wuyanzhou 时间: 2021-10-19 14:49
老师您好。若是想要选取目标原子周围2.6埃--5.2埃壳层内的的第32号原子，请问选择语句如何书写呢？
我用以下语句进行提取，也有结果输出，但是不知道是否正确。
set seltte [atomselect top "type 32 and same residue as {exwithin 2.6 of index $i} and same residue as {pbwithin 5.2 of index $i}"]

作者
Author: sobereva 时间: 2021-10-20 02:28

wuyanzhou 发表于 2021-10-19 14:49
老师您好。若是想要选取目标原子周围2.6埃--5.2埃壳层内的的第32号原子，请问选择语句如何书写呢？
我用以 ...

type后面是跟着原子类型，怎么能接原子序号
原子序号要么用index要么serial选定

作者
Author: wuyanzhou 时间: 2021-10-22 15:28

sobereva 发表于 2021-10-20 02:28
type后面是跟着原子类型，怎么能接原子序号
原子序号要么用index要么serial选定

老师，我想重新描述一下我的问题。
在我的溶液团簇中有A,B,C三种分子，现在我想找到在A分子周围，径向距离上的B分子的分布。整个思路是，先找到A,B两个分子的质心，然后通过质心距离进行判断筛选。
问题卡在了：如何通过原子序号来确定分子序号，然后得到分子的质心呢？
例如A分子上代表性的原子是氟原子和氧原子。atomselect top 后面怎么写呢？

作者
Author: sobereva 时间: 2021-10-23 04:01

wuyanzhou 发表于 2021-10-22 15:28
老师，我想重新描述一下我的问题。
在我的溶液团簇中有A,B,C三种分子，现在我想找到在A分子周围，径向距 ...

比如same fragment as element O就能把带氧原子的分子里所有原子选中
之后怎么干看你想具体怎么写脚本。如果你要得到这些分子的残基号，可以用$sel get resid获得$sel里所有原子的resid号，然后去重，就得到这些分子的resid号列表

作者
Author: wuyanzhou 时间: 2021-10-24 14:59

sobereva 发表于 2021-10-23 04:01
比如same fragment as element O就能把带氧原子的分子里所有原子选中
之后怎么干看你想具体怎么写脚本。 ...

好的，老师。非常感谢您的回复。

作者
Author: wuyanzhou 时间: 2021-10-25 15:48
本帖最后由 wuyanzhou 于 2021-10-25 15:51 编辑

老师，想麻烦您给看下这段程序。运行之后在第四行的位置显示invalid command name. (, 下载次数 Times of downloads: 154)

作者
Author: sobereva 时间: 2021-10-26 14:52

wuyanzhou 发表于 2021-10-25 15:48
老师，想麻烦您给看下这段程序。运行之后在第四行的位置显示invalid command name.

一行一行运行，把嵌套的命令拆开来一条一条执行，自然就会搞清楚

作者
Author: wuyanzhou 时间: 2021-10-26 15:09

sobereva 发表于 2021-10-26 14:52
一行一行运行，把嵌套的命令拆开来一条一条执行，自然就会搞清楚

老师，我单独运行$sel get resid也是没有问题的。但是如何把运行后的结果赋给新的链表是我遇到的问题。
第四行： set resid [$sel get resid]这条语句，运行后不能将$sel get resid的残基序号置换到名为resid链表中。
。

作者
Author: sobereva 时间: 2021-10-27 13:05

wuyanzhou 发表于 2021-10-26 15:09
老师，我单独运行$sel get resid也是没有问题的。但是如何把运行后的结果赋给新的链表是我遇到的问题。
...

试试改成set resid "[$sel get resid]"

作者
Author: tuzhidingdong 时间: 2021-11-1 22:51
请问如何在vmd中显示指定实空间范围内的VOLUME数据？

作者
Author: sobereva 时间: 2021-11-2 11:00

tuzhidingdong 发表于 2021-11-1 22:51
请问如何在vmd中显示指定实空间范围内的VOLUME数据？

看不懂你的意思
如果是只显示某个区域的格点数据的等值面，VMD没法直接实现。需要通过Multiwfn对格点数据进行屏蔽（参考Multiwfn手册4.13.4节里的例子）之后导出格点数据再在VMD里看，或者在Chimera里显示等值面的时候把不想要的地方擦除。

作者
Author: tuzhidingdong 时间: 2021-11-19 15:16

sobereva 发表于 2021-11-2 11:00
看不懂你的意思
如果是只显示某个区域的格点数据的等值面，VMD没法直接实现。需要通过Multiwfn对格点数 ...

嗯，就是这个意思，感谢社长

作者
Author: DZW 时间: 2021-12-17 20:45
请问VMD中有显示/隐藏极性氢原子和非极性氢原子的语法吗? 还是说这个得自己写脚本？

作者
Author: sobereva 时间: 2021-12-17 21:07

DZW 发表于 2021-12-17 20:45
请问VMD中有显示/隐藏极性氢原子和非极性氢原子的语法吗? 还是说这个得自己写脚本？

没有现成语法区分极性和非极性氢。除非你载入的文件有原子电荷信息，可以根据原子电荷大小区分

作者
Author: Aristotler 时间: 2022-3-8 21:23
老师请问：在gromacs 里如何选区某残基的中的全部某类元素呢？比如如何选区DME这个残基中全部的O原子（包含O1 和O2）呢？我写resname DME and name O1 and O2 感觉不对，请指教

作者
Author: sobereva 时间: 2022-3-9 05:14

Aristotler 发表于 2022-3-8 21:23
老师请问：在gromacs 里如何选区某残基的中的全部某类元素呢？比如如何选区DME这个残基中全部的O原子（包 ...

resname DME and name O1 O2

作者
Author: Aristotler 时间: 2022-3-9 07:08

sobereva 发表于 2022-3-9 05:14
resname DME and name O1 O2

非常感谢老师

作者
Author: 牧生 时间: 2022-5-7 20:15
本帖最后由牧生于 2022-5-9 09:29 编辑

请教一个问题，我做水合物模拟后，一个氯离子跑到了盒子的边缘，我想选择这个氯离子附近5 A的所有分子做分析，我使用same residue as pbwithin 5 of residue 1209, 结果得到穿过盒子边缘的两个半球形 (, 下载次数 Times of downloads: 107)

我该使用什么样的语句才能顺利取出边界上的一个完整的圆球形呢。？

作者
Author: sobereva 时间: 2022-5-9 21:56

牧生发表于 2022-5-7 20:15
请教一个问题，我做水合物模拟后，一个氯离子跑到了盒子的边缘，我想选择这个氯离子附近5 A的所有分子做分 ...

你得用trjconv处理轨迹，让氯离子在盒子里居中，并且把此时盒子外的分子wrap进盒子，然后再用VMD选择

作者
Author: 牧生 时间: 2022-5-10 11:27
本帖最后由牧生于 2022-5-11 19:30 编辑

sobereva 发表于 2022-5-9 21:56
你得用trjconv处理轨迹，让氯离子在盒子里居中，并且把此时盒子外的分子wrap进盒子，然后再用VMD选择

我使用 gmx make_ndx -f md.gro -o index.ndx 选择需要分析的氯离子

再用gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -pbc mol -center -n index.ndx，

尽管氯离子放到盒子的中间位置，但实际上氯离子仍然处在盒子边缘上，并不在盒子内部的中间。
(, 下载次数 Times of downloads: 111) (, 下载次数 Times of downloads: 124)
请帮忙再看一下。

2022.5.11 问题已经解决。

解决方法：①直接将MD结束后得到的gro转为pdb格式。②用gv打开这个pdb，会自动去掉虚原子，再存为pdb，③文本软件打开这个pdb，确认前几行有盒子的尺寸信息，如果没有盒子信息，自己手动加上。④使用Multiwfn扩展为2X2X2，存为pdb格式。⑤用VMD打开这个盒子，用语法去选择想要的原子及周围的分子，根本不用想命令，也不用处理轨迹，也不用wrap，非常简单。
看起来似乎步骤多，但都只是在打开和保存，转化格式，且不用费脑子。

有一点缺陷，但不知道是哪一步有缺陷。最终用VMD选择的小球形中，CL原子可能会变成C，但这个问题不大，自己用文本编辑软件打开最后一步的pdb，一看就能找到是哪个。

作者
Author: 一条君 时间: 2022-11-13 11:26
本帖最后由一条君于 2022-11-13 11:27 编辑

【within 5 of AAA：距离AAA 5埃以内的原子。选取时不考虑周期边界条件，用pbwithin则考虑】这点强烈推荐给周期性体系用VMD中tcl命令作结构统计，配位数统计，vmd边界原子选取等等的同学，自己用了一年才发现这个问题！！！
类似帖子如http://bbs.keinsci.com/thread-17972-1-1.html【求助如何正确统计离子的数目】；http://bbs.keinsci.com/thread-16153-1-1.html【请教：周期性体系统计分子数目的问题】

作者
Author: ginlpein 时间: 2023-1-11 22:11
社长，想问一下“一般关键词”这一批里面，哪些是直接从原文件中读取的？哪些是读取原子坐标后VMD根据程序判定自动生成的？
比如，在PDB文件没有键接信息的情况下，residue似乎是根据VMD对成键判定再对各个分子标号，而resid似乎是纯粹看原文件中的情况。
请问这方面是否有情况的总结？
不同软件里生成出来的PDB经常有很多不同的情况，比如不登记残基名称/序号、原子序号生成时候重新打乱等等情况，所以需要tcl脚本重新搞一下才能后续用的比较顺畅。

作者
Author: sobereva 时间: 2023-1-12 06:36

ginlpein 发表于 2023-1-11 22:11
社长，想问一下“一般关键词”这一批里面，哪些是直接从原文件中读取的？哪些是读取原子坐标后VMD根据程序 ...

输入文件里没有明确记录的信息肯定是VMD自己判断的，不清楚的话测试几次就知道了

作者
Author: Jus 时间: 2023-1-29 22:20

sobereva 发表于 2023-1-12 06:36
输入文件里没有明确记录的信息肯定是VMD自己判断的，不清楚的话测试几次就知道了

sob老师您好，我是想用angle.tcl脚本计算夹角，这个脚本里是用四个变量来确定两个向量，我想计算的夹角，是两个原子构成的向量与Z轴形成的夹角。我在前两个select里输入的是resid 691 and type O58和resid 691 and type O119确定的两个原子（我也不确定是不是这样写的），然后剩下的两个select不知道写啥来描述Z轴。
万分感谢老师指点！

作者
Author: sobereva 时间: 2023-1-30 23:43

Jus 发表于 2023-1-29 22:20
sob老师您好，我是想用angle.tcl脚本计算夹角，这个脚本里是用四个变量来确定两个向量，我想计算的夹角， ...

我不知道你说的angle.tcl是什么东西

作者
Author: Jus 时间: 2023-1-31 08:22
本帖最后由 Jus 于 2023-1-31 09:08 编辑

sobereva 发表于 2023-1-30 23:43
我不知道你说的angle.tcl是什么东西

老师您好，是在 http://bbs.keinsci.com/thread-14821-1-1.html 这个帖子里的angle.tcl脚本，内容如下#---------------------------------------------------
set outfile [open angle.dat w]
set select1 "protein and resid 1"
set select2 "protein and resid 2"
set select3 "protein and resid 2"
set select4 "protein and resid 3"
#---------------------------------------------------
set nf [molinfo top get numframes]
set sel1 [atomselect top "$select1"]
set sel2 [atomselect top "$select2"]
set sel3 [atomselect top "$select3"]
set sel4 [atomselect top "$select4"]
for { set i 1 } { $i <= $nf } { incr i } {
      $sel1 frame $i
      set V1 [measure center "$sel1"]
      $sel2 frame $i
      set V2 [measure center "$sel2"]
      $sel3 frame $i
      set V3 [measure center "$sel3"]
      $sel4 frame $i
      set V4 [measure center "$sel4"]
      set VA [vecsub $V1 $V2]
      set VB [vecsub $V3 $V4]
      set COSAB [expr [vecdot $VA $VB]/([veclength $VA]*[veclength $VB])]
      set ANGLE [expr acos($COSAB)*180/3.1415926]
      puts $outfile "[expr $ANGLE]"
}
close $outfile
puts "All Done!"

作者
Author: sobereva 时间: 2023-2-2 01:32

Jus 发表于 2023-1-31 08:22
老师您好，是在 http://bbs.keinsci.com/thread-14821-1-1.html 这个帖子里的angle.tcl脚本，内容如下#-- ...

我没用过那个脚本，应当问原作者

计算与Z轴夹角自己写单独的脚本就完了

北京科音分子动力学与GROMACS培训班里详细讲怎么编写脚本，对于与Z轴夹角问题也给了相关的脚本

(, 下载次数 Times of downloads: 90)

(, 下载次数 Times of downloads: 76)

作者
Author: Gzh_NJ 时间: 2024-3-13 15:37
卢老师好，请教下老师我想选择体系中的CO2和N2该如何使用选择语法，谢谢卢老师。

作者
Author: sobereva 时间: 2024-3-14 10:27

Gzh_NJ 发表于 2024-3-13 15:37
卢老师好，请教下老师我想选择体系中的CO2和N2该如何使用选择语法，谢谢卢老师。

没有结构文件没法回答
这取决于体系里都有什么、原子名、残基名、原子序号范围
根据这些信息自行判断怎么用合适的选择语句

作者
Author: Gzh_NJ 时间: 2024-3-15 11:23

sobereva 发表于 2024-3-14 10:27
没有结构文件没法回答
这取决于体系里都有什么、原子名、残基名、原子序号范围
根据这些信息自行判断怎 ...

谢谢卢老师，我再上传结构文件在另一个帖子进行求助。老师，我是想用CP2K培训班统计产物脚本时想统计产物CO2和H2O随着时间变化，尝试一段时间后没有成功，麻烦老师有空看看，谢谢老师。

作者
Author: yihanxu 时间: 2024-3-28 18:00
老师好，我想请教一下，为什么同一个残基它的residue 165、但resid 195？resid看起来对同一个蛋白的同一个残基是固定的，不同文献中这个蛋白质上的这个残基都是SER195。
谢谢老师！

作者
Author: sobereva 时间: 2024-3-28 20:47

yihanxu 发表于 2024-3-28 18:00
老师好，我想请教一下，为什么同一个残基它的residue 165、但resid 195？resid看起来对同一个蛋白的同一个 ...

residue：残基编号，从0开始
resid：残基编号，从1开始。若结构文件里有残基号则与之一致

作者
Author: yihanxu 时间: 2024-3-28 21:38

sobereva 发表于 2024-3-28 20:47
residue：残基编号，从0开始
resid：残基编号，从1开始。若结构文件里有残基号则与之一致

明白了，谢谢老师！

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