请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 1024|回复: 10

[Multiwfn资源与经验] 巨大体系的范德华表面静电势图的快速绘制方法

[复制链接]

1万

帖子

25

威望

2万

eV
积分
43663

管理员

公社社长

发表于 2019-5-7 16:11:11 | 显示全部楼层 |阅读模式
巨大体系的范德华表面静电势图的快速绘制方法

文/Sobereva@北京科音
First release: 2019-May-7  Last update: 2019-May-13


有人问我,怎么绘制他的含有多达336个原子的分子的分子表面静电势图。对这么大体系,绘制这种图看似极其困难,很难吃得消,但其实只要方法得当,哪怕只有台普通四核机子,也完全没有压力,而且绘制的图像效果极佳,本文就来介绍一下流程。绘制用的步骤和《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》(http://sobereva.com/443)里介绍的基本一样,但牵扯到一些其它过程和技巧,没读过此文的读者一定要先阅读一下。本文用到的Multiwfn可以在其官网http://sobereva.com/multiwfn上免费下载,读者用的Multiwfn版本应是2019年4月或以后更新的。

本文笔者的条件是:CPU为Intel i7-2630QM(4核),操作系统是Win7-64bit,Multiwfn是2019-May-5版,xtb是2019-Apr-16版,Gaussian是G16 A.03版,VMD是1.9.3版,这几个程序里只有Gaussian是收费的。笔者跑xtb时是在VMware装的CentOS 7.4虚拟机下跑的,因为此程序目前只有Linux版,不会装Linux的话看《在VMware 15中安装CentOS 7.6的完整过程视频演示》(http://sobereva.com/454)。

要绘制分子表面静电势图,不管通过什么途径,肯定得先产生波函数,至少得算一次单点任务。但如果你只有一台四核机子,遇见300多个原子的体系,用Gaussian算DFT就甭想了,哪怕在ORCA里用纯泛函开RI,照样很难算得动。虽然用半经验方法如PM6算几百个原子很快,但是一般的程序都不支持基于半经验方法得到的波函数去算静电势。碰到这种情况,最佳的出路就是用Grimme的xtb程序。此程序的安装和基本用法在《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》(http://sobereva.com/421)已经介绍了,不会用xtb的人一定要看一下。

先把那个336原子的大分子的结构文件载入Multiwfn(支持mol、mol2、pdb、xyz、gjf等众多格式),进入主功能100的子功能2,选择相应选项导出为xyz文件,比如叫336.xyz。然后把此文件随便放到一个空目录下,在此目录下运行此命令通过xtb做单点计算:xtb 336.xyz --molden(如果这个结构是未经优化过的,可再加上--opt让xtb优化此结构)。在笔者的4核机子上只花了20秒钟就算完了,在当前目录下出现了molden.input,这是Molden输入文件,可以被Multiwfn读取并做波函数分析。

把Multiwfn目录下的settings.ini里的nthreads设为机子里实际的CPU物理核心数,使得CPU运算能力能在接下来的计算中充分发挥。

下面产生此体系的电子密度格点数据。启动Multiwfn,载入molden.input,然后依次输入以下命令:
5  //计算格点数据
1  //计算电子密度
3  //高质量格点
计算过程花费仅不到1分钟,然后选2导出格点数据,此时当前目录下出现了density.cub,将之改名为density1.cub备用。

然后我们要计算静电势格点数据,用Multiwfn调用Gaussian目录下的cubegen速度会很快,但是这要求必须以fch作为输入文件,所以我们要把当前的波函数导出为fch文件。接着在Multiwfn里输入
0  //返回主菜单
100  //其它功能part 1
2  //导出文件
7  //导出fch文件
336.fch  //输出文件名
马上当前目录下就有了336.fch。现在退出Multiwfn。

我们修改Multiwfn目录下的settings.ini,将其中的cubegenpath路径设为你机子里实际的Gaussian目录下的cubegen可执行文件的路径,比如D:\study\G16W\cubegen.exe。关于这点,在http://sobereva.com/435里有更多说明。

当前这个体系很大,当Multiwfn调用cubegen去处理这个文件时,会由于G16的cubegen默认可用内存上限不够而报错。为解决这个问题,进入操作系统的控制面板,选择“系统”-“高级系统设置”-“高级”-“环境变量”,然后在用户变量里点击“新建”,变量名填GAUSS_MEMDEF,变量值填2GB,然后点“确定”。之后cubegen处理fch文件时最大可用内存就为2GB了,对当前体系完全够了。

启动Multiwfn,依次输入
336.fch
5  //计算格点数据
12  //计算静电势
3  //高质量格点
现在会看到Multiwfn正在调用cubegen计算静电势。这一步是整个过程最耗时的一步,花了50分钟(而在笔者的2*E5-2696v3 36核服务器上三分多钟就能算完)。然后选2导出格点数据,此时当前目录下出现了totesp.cub,将之改名为ESP1.cub。

现在把density1.cub和ESP1.cub都拷到VMD目录下,把Multiwfn目录下的examples\drawESP目录下的ESPiso.vmd也拷到VMD目录下。启动VMD,在文本窗口输入source ESPiso.vmd并回车,马上就可以看到下图:

1.jpg

这图虽然效果已经很好了,但看着还有点乱,因此我们进入Graphics - Materials,选里面的EdgyGlass材质(这是用于显示当前分子表面用的材质),然后适当调节各个参数令效果尽可能好,最后改成下面这样:

2.png

然后进入Graphics - Representation,点击对应显示分子结构的rep,把CPK显示方式切换为Licorice,并让键的粗细从默认的0.3改为0.2。然后点击对应Isosurface显示方式的那个rep,选择Trajectory标签页,把Color Scale Data Range范围设成比默认-0.03~0.03更宽的-0.04~0.04使得颜色变化更柔和些。最后图像效果如下,两个视角都给了,可见效果相当令人满意!(可能有人觉得这个图还没一开始看到的酷炫,但此时的图更容易分辨分子表面上不同区域静电势的差异)

3.jpg

图中红色和蓝色分别对应于静电势为正和为负的区域。如果你想把色彩刻度轴加上去也很容易,做法在http://sobereva.com/443文中已经介绍了。

本文的做法绘制甚至四五百个原子的体系的静电势图也是完全没问题的,哪怕只有一台四核机子,也完全可以算得动,只不过是时间长一些而已。本文给出的静电势分布其实算不上很准确,毕竟xtb程序做的GFN-xTB计算只相当于半经验层面的DFT方法,用的是极小基,不过对于几百个原子的静电势图的绘制,谁也不会那么关注定量大小,只要定性正确就足够了,当前绘制的图完全可以满足这个需要。本文的例子也充分体现出,对于好几个百原子的精度要求不很高的波函数分析,xtb+Multiwfn可谓是黄金组合。

下面是xtb产生的波函数与精度理想的B3LYP/def2-TZVP波函数绘制的分子表面静电势图的对比,可见差异不是很大,基于xtb的结果已经完全可以满足需要了。
compare.jpg

不过,如果你想基于xtb的波函数按照《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》(http://sobereva.com/159)去定量考察分子表面静电势特征,那么xtb的波函数质量还是差点意思。比如对于上面这个多巴胺的分子表面静电势最大、最小点数值,基于xtb算的结果与基于B3LYP/def2-TZVP波函数算的结果甚至能差出20 kcal/mol的程度。而对于计算拟合静电势电荷,比如CHELPG,基于xtb波函数算的结果与基于B3LYP/def2-TZVP算的结果差异最大的可达0.23,这也不是个小数目。所以,xtb波函数对应的静电势可以满足肉眼观看、定性分析需要,但定量层面上还明显达不到准确的标准。

评分

参与人数 7eV +33 收起 理由
wangyanforward + 3 好物!已收藏!
kulaomega + 5 牛!
Novice + 5 我很赞同
zsu007 + 5 牛!
lerel + 5 好物!
exity + 5 精品内容
alonewolfyang + 5 好物!已收藏

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

58

帖子

0

威望

533

eV
积分
591

Level 4 (黑子)

发表于 2019-5-7 16:20:05 | 显示全部楼层
...",花了50分钟"...

58

帖子

0

威望

920

eV
积分
978

Level 4 (黑子)

发表于 2019-5-7 19:37:33 | 显示全部楼层
点赞!我以前是用orca开RI算一个单点分析过上百个原子,,贼慢,这个方法快的太多了

115

帖子

0

威望

601

eV
积分
716

Level 4 (黑子)

发表于 2019-5-7 21:58:27 | 显示全部楼层
其实我已经用过这个方法了,哈哈,就是没有对比过xTB的结果和Gaussian算的差异是否显著,不知道社长知道不

1万

帖子

25

威望

2万

eV
积分
43663

管理员

公社社长

 楼主| 发表于 2019-5-8 09:35:49 | 显示全部楼层
exity 发表于 2019-5-7 16:20
...",花了50分钟"...

若有个像样的服务器,整个绘制过程从头到尾不会超过10分钟。文中刻意用4核机子举例,主要目的是说明,哪怕是不是专门做计算的组,没有服务器,用个普通PC机按此文做法也能算得动。
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

1万

帖子

25

威望

2万

eV
积分
43663

管理员

公社社长

 楼主| 发表于 2019-5-8 10:00:30 | 显示全部楼层
Novice 发表于 2019-5-7 21:58
其实我已经用过这个方法了,哈哈,就是没有对比过xTB的结果和Gaussian算的差异是否显著,不知道社长知道不

文末加入了对比图
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

1

帖子

0

威望

13

eV
积分
14

Level 1 能力者

发表于 2019-5-10 21:25:34 | 显示全部楼层
请问老师,本文所使用的xtb需要与gaussian联用吗

1万

帖子

25

威望

2万

eV
积分
43663

管理员

公社社长

 楼主| 发表于 2019-5-10 23:09:01 | 显示全部楼层
haowei 发表于 2019-5-10 21:25
请问老师,本文所使用的xtb需要与gaussian联用吗

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

43

帖子

0

威望

268

eV
积分
311

Level 3 能力者

发表于 2019-5-12 18:41:03 | 显示全部楼层
xtb有个--esp calculate electrostatic potential on VdW-grid选项,要是xtb可以计算cube格式电势静数据,再结合xtb生成的电子密度cube数据,是不是可以绘制更大分子的表面静电势图?


1万

帖子

25

威望

2万

eV
积分
43663

管理员

公社社长

 楼主| 发表于 2019-5-13 00:17:12 | 显示全部楼层
hzamber 发表于 2019-5-12 18:41
xtb有个--esp calculate electrostatic potential on VdW-grid选项,要是xtb可以计算cube格式电势静数据, ...

这个选项只是生成分子表面上的点的静电势,没法实现本文的效果
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

1万

帖子

25

威望

2万

eV
积分
43663

管理员

公社社长

 楼主| 发表于 2019-5-13 00:29:36 | 显示全部楼层
在本文末尾添加了一段,强调了xtb给出的波函数用于定量分析还精度不够
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。达5000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2019-7-21 07:04 , Processed in 0.159944 second(s), 28 queries .

快速回复 返回顶部 返回列表