计算化学公社

标题: 将GROMACS的原子电荷信息读入VMD的方法 [打印本页]

作者
Author: sobereva 时间: 2017-3-27 07:58
标题: 将GROMACS的原子电荷信息读入VMD的方法
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)下越红代表原子电荷越负，越蓝代表原子电荷越正。

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

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

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

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

作者
Author: 我本是个娃娃 时间: 2017-3-27 08:11
沙发是我的！

作者
Author: ruanyang 时间: 2017-3-27 08:29
非常不错的脚本，显示方式也很nice顶

作者
Author: Mirror 时间: 2017-3-27 08:44
不错，谢谢分享

作者
Author: guoy14iccas 时间: 2017-3-27 14:20
感谢sob老师，亲测GROMACS-4.6.7 可用。

作者
Author: diaok 时间: 2017-3-28 17:11
以前我也用的这种方法

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

作者
Author: sobereva 时间: 2017-3-28 21:41

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

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

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

作者
Author: diaok 时间: 2017-3-30 11:09

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

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

作者
Author: lao7 时间: 2019-12-23 11:07
这个文件是不是能总IR振动分析？请问哪里有Gromacs做IR分析的案例？谢谢！

作者
Author: sobereva 时间: 2020-1-3 20:05
由于发现文中的gmxoutchg与gromacs 2019不兼容，今日更新了gmxoutchg 1.1版，对gmx 2019也支持了

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

作者
Author: sobereva 时间: 2020-6-11 23:37

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

笔误，已改

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

作者
Author: lishine 时间: 2022-8-27 11:36

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

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

作者
Author: uenh1998 时间: 2022-8-29 20:11

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

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

作者
Author: mtlin 时间: 2024-2-19 14:14
Extract atomic charges from formatted .tpr file
Version 1.1, release date: 2020-Jan-3
Programmed by Sobereva (sobereva@sina.com)

Number of molecules:       2
Maximum number of atoms in a molecule:       15
Number of molecular type:       2
Number of atoms in molecular type 0: 15
Number of atoms in molecular type 1:21960
free(): invalid pointer
已放弃 (核心已转储)
各位老师，使用gmxoutchg1.1版本出现了这个报错，请问是什么原因呢？我是GROMACS version: 2022.6
这个报错是意味着不适用吗？

作者
Author: sobereva 时间: 2024-2-19 14:39

mtlin 发表于 2024-2-19 14:14
Extract atomic charges from formatted .tpr file
Version 1.1, release date: 2020-Jan-3
Programme ...

建议用http://bbs.keinsci.com/thread-37839-1-1.html提供的读取tpr的插件，更方便，直接载入tpr就有原子电荷信息了

作者
Author: mtlin 时间: 2024-2-19 14:47

sobereva 发表于 2024-2-19 14:39
建议用http://bbs.keinsci.com/thread-37839-1-1.html提供的读取tpr的插件，更方便，直接载入tpr就有原子 ...

老师，在linux系统上除了您开发的这个gmxoutchg，还有别的方法吗？

作者
Author: may123456 时间: 2024-5-7 15:00

uenh1998 发表于 2022-8-29 20:11
十分感谢您的指点！我试下！

请问有解决吗？ txt的格式也是UTF-8，赋值代码也加了，但是一运行就发生闪退的现象

作者
Author: kazuha 时间: 2025-4-24 11:08

may123456 发表于 2024-5-7 15:00
请问有解决吗？ txt的格式也是UTF-8，赋值代码也加了，但是一运行就发生闪退的现象

我也还是不行，请问有解决嘛

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