谈谈VMD可视化程序的连接关系的判断和设置问题

sobereva

谈谈VMD可视化程序的连接关系的判断和设置问题

On the determination and setting of the connection relationship in VMD visualization program

文/Sobereva@北京科音

First release: 2020-Feb-25   Last update: 2023-Dec-29

化学体系可视化程序VMD的用户越来越多，经常有人在网上问我怎么在VMD里修改键连关系，这里就专门说一下这个问题。本文内容对应于撰文时最新的VMD正式版1.9.3，对其它版本可能适用也可能不适用。


当VMD载入pdb、xyz、gro等格式的结构文件时，程序就会自动判断连接关系。判断规则和Multiwfn、GaussView等程序类似，都是根据两个原子间距离以及这两个原子半径之和来判断的，参看此文的相关说明：《谈谈原子间是否成键的判断问题》（http://sobereva.com/414）。VMD里默认的line方式，以及常用的CPK、Licorice等显示风格都是根据连接关系显示的。

具体来说，若原子间距离小于两个原子半径之和的0.6倍就会被VMD视为成键。原子半径是VMD内置的，大部分是Bondi范德华半径。元素的默认半径没法在VMD图形界面里或配置文件里改，但载入结构后可以通过命令行改特定原子的半径。

如果两个原子间没有被判断为成键，但你希望让它们之间显示键，就选择Mouse - Add/Remove Bonds，然后点击这两个原子的正中心，它们就连上了。如果你点击两个已经成键的原子，它们之间的键就会被去掉。

在你点击原子时，VMD会自动给你点击的原子加上标签并显示在图上，如果不想显示，就按键盘上的1键，再点击有标签的原子，标签就被隐藏了。如果你想把一批标签都隐去或者彻底删掉，就进入Graphics - Labels，按住鼠标左键选中一批，然后点Hide按钮就可以隐藏，点Delete按钮就可以删掉。如果想让批量删除标签的操作尽可能省事，对于Windows版VMD，你可以在VMD目录下的vmd.rc里加入如下内容
proc da {} {
label delete Atoms all
}
启动VMD后，随时都可以在文本窗口里输入自定义命令da来删除所有原子标签。

如果有很多同类的键，它们的连接关系判断得都和你要求的不符，一个一个靠鼠标点击来修改连接关系很麻烦。此时可以在Graphics - Representation界面里把Drawing Method改为DynamicBonds，这代表基于当前结构实时判断成键并显示，判断阈值可以用旁边的Distance Cutoff控件来修改，恰当改成令显示的连接方式满足你要求的情况即可。但这种显示方式有个缺点就是不显示原子球，而只显示键，而且键是空心的，很不好看。因此为了让原子球也显示出来，你可以点Add Rep增加一个Representation，把Drawing Method设成VDW，并且把Sphere Scale改小到合适。

有的时候在用DynamicBonds方式显示时，通过调节Distance Cutoff虽然让A-B之间键连方式满足要求了，但同时C-D之间的连接方式变得却与期望的不符了，这种时候就需要灵活变通，多设几个Representation，每个都恰当用不同的选择语句选定绘制的区域（参考《VMD里原子选择语句的语法和例子》http://sobereva.com/504），并用恰当的Drawing Method和参数设定，从而让显示出的效果叠加起来正好满足你的要求。

有的文件格式包含了连接关系信息。比如AMBER、NAMD的拓扑文件都可以被VMD载入，载入的时候就相当于提供了连接关系了，之后再往这个ID里载入结构文件或轨迹时也因此都是按照拓扑文件里的连接关系来显示。VMD虽然不能载入GROMACS的拓扑文件，但是如果安装了http://bbs.keinsci.com/thread-37839-1-1.html中提供的tpr插件的话，载入GROMACS的tpr时也会从中读取正确的连接关系。还有的结构文件格式自己也带了连接关系段落，典型的就是mol2格式（mol格式也有，但VMD不能载入），相减《谈谈记录化学体系结构的mol2文件》（http://sobereva.com/655）。你用文本编辑器打开它的话，会看到@<TRIPOS>BOND字段，这部分记录的是键的形式，诸如5 3 27 2就代表这是第5号键，对应3与27原子间的双键。当VMD载入mol2格式文件的时候，就不会自动根据结构来判断连接关系，而是直接用文件里的连接关系。pdb格式里有个CONECT字段，也是记录连接关系用的，但是VMD并不会采用这里的连接关系。如果你想让VMD直接利用这里的连接关系，方法看《使VMD根据pdb文件中的CONECT字段设定原子连接关系》（http://sobereva.com/121）。

在GaussView里可以编辑键连关系，而且从GaussView 6开始，可以自定义不同元素之间成键的判断阈值，这在前述的《谈谈原子间是否成键的判断问题》文中已经说过了。改过之后只点击Edit - Rebond就可以按照用户自定义的阈值重新判断成键，因此在批量修改特定类型元素间的成键方式的时候比较方便。GaussView可以保存mol2格式的文件，里面的连接关系信息正对应于你在GaussView图形界面里当前看到的。因此，可以先让GaussView正确显示出你期望的连接关系，保存成mol2文件，再用VMD读取，此时VMD里显示的连接关系就和GaussView里精确一致了。

注意VMD里没法以GaussView的风格来显示双键、三键等成键形式。比如载入VMD的mol2文件离记录了某两个原子间形成的是二重键，但在VMD的Line、CPK、Licorice显示风格显示时看起来和单键一样。VMD的Bonds显示风格倒是能区分形式键级，但双键和三键分别是以两个、三个套筒形式显示的，非常难看，没实际意义。

VMD自身也可以保存mol2文件，因此你在VMD里通过Add/Remove bonds好不容易把连接关系改合适后，不妨通过File - Save coordinate功能保存成mol2文件，以后再次显示时就免得再修改一遍了。

常有人问我怎么载入动力学轨迹后，虽然第一帧的成键方式正确，但播放轨迹时成键方式有时连得乱七八糟，甚至跨越盒子。这是因为如前所述，连接关系要么是从拓扑文件里读取的（比如看Amber轨迹前要载入拓扑文件），要么是载入结构文件时自动判断的（比如看GROMACS轨迹前通常要载入与模拟体系原子信息对应的gro或pdb文件），或者是根据轨迹文件的第一帧结构判断的（比如载入molclus、xtb程序产生的多帧xyz文件），之后再载入轨迹或播放轨迹时都不会改变连接关系。由于动力学过程中可能发生分子解离、异构化，或者由于周期边界条件的原因分子时而完整时而被盒子截断，因此内存里记录的连接关系列表可能对于轨迹的某些帧并不适用，极度影响正常观看。像这种情况，可以用DynamicBonds显示方式，这样就会按照当前帧的坐标实时判断成键并显示（但并不会改写内存里的连接关系列表）；但如果你想用其它显示方式，就得让VMD根据当前这一帧的结构重新判断连接关系并改写内存里的连接关系列表，做法是在文本窗口输入mol bondsrecalc all; topo retypebonds。此时连接关系就是满足当前这一帧的情况了，而对其它帧可能就不适合了。为了方便，你可以定义重新判断连接关系的快捷键，甚至可以让VMD自动对每一帧实时判断连接关系，做法见《VMD初始化文件(vmd.rc)我的推荐设置》（http://sobereva.com/545）。

如果你想把当前内存里的连接关系列表以文本方式显示出来，就在文本窗口输入[atomselect top all] getbonds，输出的信息比如
{1 2 3 4 7} 0 0 0 {5 6 7 0} 4 4 {4 0}
这是个列表，里面共有8个成员（有的成员自身就是个列表），对应体系一共8个原子的成键情况，比如第0号原子与1、2、3、4、7成键，第2号原子与第0号原子成键，第3号原子也与第0号原子成键...注意这体现的是连接关系，而非键级。

连接关系也可以从自行提供的列表里读取。比如运行
[atomselect top all] setbonds {{1 2 3 4 7} 0 0 0 {5 6 7 0} 4 4 {4 0}}
就代表把当前体系的连接关系改成{1 2 3 4 7} 0 0 0 {5 6 7 0} 4 4 {4 0}所定义的。

Uus/pMeC6H4-/キ

lyj714 发表于 2023-8-29 09:31

这一段代码是出自哪个文件的呢？虽然看到命令行更改原子半径的语句了，但仍想试试看能不能动源文件；除了改掉一堆默认2.0埃的半径，还有点想补齐118个元素（

顺便说一句，隔壁ovito的原子半径定义可以在源码里面找到，倒数第二列的范德华半径与VMD一致，有些倒数第三列的共价半径则出自OpenBabel。

hanlu

sobereva 发表于 2023-11-10 07:07
合适

谢谢老师!

sobereva

hanlu 发表于 2023-11-9 20:25
哦哦，我看错了老师，谢谢老师，那老师我这个命令合适么？

合适

hanlu

sobereva 发表于 2023-11-9 17:32
本来就不生成md1.gro，你-o指定的明明是.xtc

哦哦，我看错了老师，谢谢老师，那老师我这个命令合适么？

sobereva

hanlu 发表于 2023-11-9 17:21
老师，我用的这个命令 ，gmx trjconv -f md.xtc -o md1.xtc -pbc mol -s md.tpr ，没有生成新的md1.gro,不 ...

本来就不生成md1.gro，你-o指定的明明是.xtc

hanlu

老师，我用的这个命令 ，gmx trjconv -f md.xtc -o md1.xtc -pbc mol -s md.tpr ，没有生成新的md1.gro,不过有md1.xtc文件，想问一下为啥没产生md1.gro文件，不过载入轨迹之后，没有乱建现象。

Gzh_NJ

lyj714 发表于 2023-8-29 09:31

感谢老师

Gzh_NJ

sobereva 发表于 2023-8-28 19:58
看源代码，或者载入个体系，然后查询原子的radius属性

谢谢卢老师

lyj714

Gzh_NJ 发表于 2023-8-28 19:06
谢谢卢老师，请教下老师VMD中内置的半径数值从哪可以查看？然后准备写到文章SI中作为判断成键断键的标准 ...

1. /*
   
2. * corresponding table of VDW radii.
   
3. * van der Waals radii are taken from A. Bondi,
   
4. * J. Phys. Chem., 68, 441 - 452, 1964,
   
5. * except the value for H, which is taken from R.S. Rowland & R. Taylor,
   
6. * J.Phys.Chem., 100, 7384 - 7391, 1996. Radii that are not available in
   
7. * either of these publications have RvdW = 2.00 Å.
   
8. * The radii for Ions (Na, K, Cl, Ca, Mg, and Cs are based on the CHARMM27
   
9. * Rmin/2 parameters for (SOD, POT, CLA, CAL, MG, CES) by default.
   
10. */
    
11. static const float pte_vdw_radius[] = {
    
12.     /* X  */ 1.50f, 1.20f, 1.40f, 1.82f, 2.00f, 2.00f,
    
13.     /* C  */ 1.70f, 1.55f, 1.52f, 1.47f, 1.54f,
    
14.     /* Na */ 1.36f, 1.18f, 2.00f, 2.10f, 1.80f,
    
15.     /* S  */ 1.80f, 2.27f, 1.88f, 1.76f, 1.37f, 2.00f,
    
16.     /* Ti */ 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f,
    
17.     /* Ni */ 1.63f, 1.40f, 1.39f, 1.07f, 2.00f, 1.85f,
    
18.     /* Se */ 1.90f, 1.85f, 2.02f, 2.00f, 2.00f, 2.00f,
    
19.     /* Zr */ 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f,
    
20.     /* Pd */ 1.63f, 1.72f, 1.58f, 1.93f, 2.17f, 2.00f,
    
21.     /* Te */ 2.06f, 1.98f, 2.16f, 2.10f, 2.00f,
    
22.     /* La */ 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f,
    
23.     /* Eu */ 2.00f, 2.00f, 2.00f, 2.00f, 2.00f,
    
24.     /* Er */ 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f,
    
25.     /* W  */ 2.00f, 2.00f, 2.00f, 2.00f, 1.72f, 1.66f,
    
26.     /* Hg */ 1.55f, 1.96f, 2.02f, 2.00f, 2.00f, 2.00f, 2.00f,
    
27.     /* Fr */ 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 1.86f,
    
28.     /* Np */ 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f,
    
29.     /* Md */ 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f, 2.00f,
    
30.     /* Mt */ 2.00f, 2.00f, 2.00f
    
31. };

复制代码

sobereva

Gzh_NJ 发表于 2023-8-28 19:06
谢谢卢老师，请教下老师VMD中内置的半径数值从哪可以查看？然后准备写到文章SI中作为判断成键断键的标准 ...

看源代码，或者载入个体系，然后查询原子的radius属性

Gzh_NJ

sobereva 发表于 2023-8-28 15:49
每个原子都有radius属性（范德华半径，大多来自Bondi范德华半径，但有些来路不明），原子间距离小于radiu ...

谢谢卢老师，请教下老师VMD中内置的半径数值从哪可以查看？然后准备写到文章SI中作为判断成键断键的标准。

sobereva

Gzh_NJ 发表于 2023-8-28 15:29
卢老师好，请教下VMD中观察轨迹时候具体的成键阈值是从哪里可以看到？还是需要按照1.6倍手动计算，这是我在 ...

每个原子都有radius属性（范德华半径，大多来自Bondi范德华半径，但有些来路不明），原子间距离小于radius之和的0.6倍才会被视为成键。内置的半径没法改，但可以载入后对原子的radius属性进行设置，然后用mol bondsrecalc all重新计算

Gzh_NJ

卢老师好，请教下VMD中观察轨迹时候具体的成键阈值是从哪里可以看到？还是需要按照1.6倍手动计算，这是我在发表文章中SI读到的成键阈值。

sobereva

Frank 发表于 2022-10-17 15:54
社长，请教个问题，如果是VMD可视化电化学界面PBC体系(XY有PBC，Z没有PBC)，应该如何使用wrap命令？此外， ...

看你载入VMD的文件的盒子信息怎么记录的。如果Z方向没有PBC，不想在Z方向被某种方式wrap的话，就把Z方向尺寸设得巨大

pbc wrap接-compound res/chain/fragment使得残基/链/片段保持完整

Frank

社长，请教个问题，如果是VMD可视化电化学界面PBC体系(XY有PBC，Z没有PBC)，应该如何使用wrap命令？此外，怎么保证盒子边缘的分子不会因为wrap而断开?