将GROMACS的原子电荷信息读入VMD的方法

sobereva · 发表于 Post on 2017-3-27 07:58:17

2023-Jun-8注：推荐用http://bbs.keinsci.com/thread-37839-1-1.html里介绍的VMD载入tpr文件的插件，比本文的做法更方便。VMD载入tpr文件后直接就有了tpr中的原子电荷信息

将GROMACS的原子电荷信息读入VMD的方法

Method to load GROMACS atomic charge information into VMD

文/Sobereva @北京科音

First release: 2017-Mar-27 Last update: 2022-Apr-12

VMD可以载入GROMACS的gro、trr、xtc文件，但不会载入速度、受力、电荷。以前写过怎么让VMD把GROMACS产生的速度和受力载入并绘制出来的博文：
《使VMD实时显示gromacs轨迹中原子的受力》（http://sobereva.com/36）
《使VMD读入Gromacs产生的trr轨迹中速度信息的方法》（http://sobereva.com/117）
有很多VMD里面的插件，比如显示偶极矩、计算静电势、绘制红外光谱的插件都依赖于原子电荷，有时候我们还需要编写依赖于原子电荷的脚本、使用依赖于原子电荷的选择语句，因此把GROMACS模拟时用的原子电荷载入VMD也是非常重要的，这里说说怎么做。

首先在这里下载笔者开发的gmxoutchg程序：http://sobereva.com/soft/gmxoutchg_1.1.rar
其中gmxoutchg.exe是Windows版可执行文件，没后缀的gmxoutchg是Linux下可执行文件，gmxoutchg.f90是源文件。经测试此程序兼容GROMACS 2016.1、2018.8和2019.3版，对其它版本兼容性未经测试。

tpr文件里蕴含了模拟所需的一切信息，包括原子电荷。用GROMACS里的dump命令将二进制的tpr文件转化成文本文件dump.txt，执行以下命令：
gmx dump -s test.tpr > dump.txt

然后将dump.txt放到gmxoutchg.exe所在目录下，运行gmxoutchg.exe，程序就会解析其中的数据，在当前目录下产生charge.txt。其中包含体系所有原子的原子电荷，顺序和结构文件里的原子顺序完全一致。其中第一列、第二列、第三列分别是原子电荷、分子序号、原子序号。

将charge.txt放到VMD目录后，在载入对应的结构/轨迹后，可以用以下tcl脚本将原子电荷从charge.txt中读入。

set sel [atomselect top all]
set natom [$sel num]
set rdchg [open "charge.txt" r]
set chglist {}
for {set iatm 0} {$iatm<=[expr $natom-1]} {incr iatm} {
gets $rdchg line
scan $line "%f" chg
lappend chglist $chg
}
$sel set charge $chglist
close $rdchg

如果想验证是否正确载入了，可以用[atomselect top all] get charge命令把所有原子的原子电荷输出出来（当然，也可以把all换成选择语句只输出指定部分的原子电荷）。

也可以将Coloring method设成Charge，直接从颜色上检验是否正确载入了。下面是乙醇盒子，在默认色彩刻度(RWB)下越红代表原子电荷越负，越蓝代表原子电荷越正。

顺带一提，对于一些体系通过恰当设定显示方式，可以令原子电荷分布一目了然。比如下图的分子，一个显示方式设为Licorice并且把键调细，另一个显示方式是CPK，把圆球调大，用透明材质，用charge来着色，效果挺不错。

有了原子电荷信息可以直接使用VMD的依赖于原子电荷的插件了，比如Extensions-Visualization-Dipole Moment Watcher可以观看基于原子电荷计算的偶极矩矢量，就是上图的红色箭头。

我本是个娃娃 · 发表于 Post on 2017-3-27 08:11:29

沙发是我的！

ruanyang · 发表于 Post on 2017-3-27 08:29:57

非常不错的脚本，显示方式也很nice顶

Mirror · 发表于 Post on 2017-3-27 08:44:47

不错，谢谢分享

guoy14iccas · 发表于 Post on 2017-3-27 14:20:36

感谢sob老师，亲测GROMACS-4.6.7 可用。

diaok · 发表于 Post on 2017-3-28 17:11:43

以前我也用的这种方法

不过前不久在别的群里有人提了个更方便的方法
如果只是想要体系的电荷 gmx可以输出pqr格式文件
gmx editconf -f a.tpr -mead a.pqr
直接载入vmd就补上了电荷和半径

sobereva · 发表于 Post on 2017-3-28 21:41:29

diaok 发表于 2017-3-28 17:11
以前我也用的这种方法

不过前不久在别的群里有人提了个更方便的方法

这样做的一个缺点是原子电荷保存精度太低，只有两位小数，整个分子所有原子电荷加和会往往能偏离整数不少（特别是进一步通过脚本计算静电作用能的话误差会较大）

diaok · 发表于 Post on 2017-3-30 11:09:56

sobereva 发表于 2017-3-28 21:41
这样做的一个缺点是原子电荷保存精度太低，只有两位小数，整个分子所有原子电荷加和会往往能偏离整数不 ...

的确是这样。。
这种方法只能用来初步可视化电荷的分布了

lao7 · 发表于 Post on 2019-12-23 11:07:04

这个文件是不是能总IR振动分析？请问哪里有Gromacs做IR分析的案例？谢谢！

sobereva · 发表于 Post on 2020-1-3 20:05:08

由于发现文中的gmxoutchg与gromacs 2019不兼容，今日更新了gmxoutchg 1.1版，对gmx 2019也支持了

mol · 发表于 Post on 2020-6-11 17:04:48

sob老师您好，在默认色彩刻度(BWR)下越红代表原子电荷越负，越蓝代表原子电荷越正。这句感觉有误，因该是“默认色彩刻度(RWB)下越红代表原子电荷越负，越蓝代表原子电荷越正吧”

sobereva · 发表于 Post on 2020-6-11 23:37:58

mol 发表于 2020-6-11 17:04
sob老师您好，在默认色彩刻度(BWR)下越红代表原子电荷越负，越蓝代表原子电荷越正。这句感觉有误，因该是“ ...

笔误，已改

uenh1998 · 发表于 Post on 2022-1-5 09:35:28

早上测试了一下1.1版与gmx2020版的兼容情况，发现不兼容。

lishine · 发表于 Post on 2022-8-27 11:36:40

uenh1998 发表于 2022-1-5 09:35
早上测试了一下1.1版与gmx2020版的兼容情况，发现不兼容。

我下载了1.1版本，也出现了您这种情况，以下是我调试结果和解决方法：
1. Error finding molblock的根本原因在于：gmx dump出来的文件编码格式不对，使得sob老师的代码读出来的都是乱码。所以我另存为了一下，使他变成UTF-8的格式，成功读取。
2. 代码有点小问题（或许是我误删除），源代码中，第一次使用maxtype时没有给他赋值，导致初始值时负数，因此需要加上一句赋值代码

uenh1998 · 发表于 Post on 2022-8-29 20:11:44

lishine 发表于 2022-8-27 11:36
我下载了1.1版本，也出现了您这种情况，以下是我调试结果和解决方法：
1. Error finding molblock的根本 ...

十分感谢您的指点！我试下！

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[GROMACS] 将GROMACS的原子电荷信息读入VMD的方法

评分 Rate

评分 Rate

评分 Rate

浏览过的版块