计算化学公社
标题: 绘制有填色效果的用于弱相互作用分析的RDG散点图的方法 [打印本页]
作者Author: sobereva 时间: 2017-12-17 20:41
标题: 绘制有填色效果的用于弱相互作用分析的RDG散点图的方法
绘制有填色效果的用于弱相互作用分析的RDG散点图的方法
Method to draw RDG scatter map with coloring effect for weak interaction analysis
文/Sobereva @北京科音 First release: 2017-Dec-17 Last update: 2022-Aug-11
2010年提出的用于分析弱相互作用的RDG方法（文献中也普遍叫NCI方法）已被广泛用于考察各种分子间/分子内弱相互作用了，笔者也写过不少相关文章，不了解此方法的者看《使用Multiwfn图形化研究弱相互作用》（http://sobereva.com/68）以及里面提及的相关文章，特别是推荐看《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（http://sobereva.com/667）里介绍的笔者的综述。能做RDG分析的程序不少，Multiwfn（http://sobereva.com/multiwfn）是其中最为流行、强大、好用的。在上面的博文里已经介绍了怎么在Multiwfn中直接绘制RDG vs sign(lambda2)rho的散点图来考察弱相互作用。一些文章里的这种散点图还加上了填色效果，可以使得对应不同横坐标的spike的颜色一目了然，便于与Multiwfn+VMD绘制的RDG填色等值面图相对应来讨论问题。其实这种图稍有photoshop使用技能的人都可以不太困难地作出来，就是把VMD的色彩刻度条在ps里拉伸成与散点图作图范围相同的大小，垫在Multiwfn给出的散点图下方的图层，然后再把散点图的图层当中作图区域的黑色部分以色彩范围选择方式选中，删除，透出来下层的色彩刻度层即可。不过肯定有不少人嫌这种做法麻烦，此文介绍一种利用gnuplot程序的简单快捷的方法绘制这种填色RDG散点图。
读者请使用2019-Aug-24及之后更新的Multiwfn。这里用通过苯酚二聚体来示例，相应的波函数文件是Multiwfn文件包里examples目录下的PhenolDimer.wfn。本文的操作在《使用Multiwfn做NCI分析展现分子内和分子间弱相互作用》（https://www.bilibili.com/video/av71561024）里也有视频演示。
启动Multiwfn，依次输入以下命令，让Multiwfn把此体系的RDG vs sign(lambda2)rho的散点数据导出到当前目录下的output.txt中。
examples\PhenolDimer.wfn
20 //弱相互作用图形化分析
1 //NCI分析
3 //高质量格点
2 //导出散点数据
gnuplot是个免费的轻量级的基于命令行的数据作图程序，各种系统都支持，可以在这里下载：http://www.gnuplot.info。本文用的是gnuplot 5.4 Windows版。将output.txt放到gnuplot目录下的bin子目录下，然后将Multiwfn目录下的examples\scripts\RDGscatter.gnu这个绘图脚本也拷到此目录下。开启操作系统的命令行模式（例如Windows下的cmd环境）并进入此目录，运行命令gnuplot RDGscatter.gnu（对于Windows用户，这一步不知道怎么弄的话直接把RDGscatter.gnu拖到gnuplot.exe图标上也行），就会在当前目录下产生RDGscatter.ps，这就是填色散点图的postscript格式的文件了。这是一种矢量图形格式，可无损缩放，很多程序都可以查看。比如可以直接用acrobat打开，打开后可以无损缩放。也可以用photoshop打开，打开的时候可以选择产生像素为多大的图片。如果机子里装了ghostscript程序，也可以用小巧且强大的看图程序irfanview观看。如果你懒得装单机程序，也可以用免费的在线程序https://cloudconvert.com/image-converter把ps格式转成常见图像格式。此例效果如下：
(, 下载次数 Times of downloads: 124)
RDGscatter.gnu脚本里有很多参数可以设定，比较关键的参数就是X,Y轴以及色彩刻度轴的上下限（x/y/cbrange后面的值）、标签的数值范围和步长（x/y/cbtic后面的值）、散点的大小（pointsize后面的值），以及色彩刻度的定义。笔者习惯在VMD中用-0.035~0.02来对RDG等值面着色，色彩刻度是默认的蓝-绿-红，因此脚本中可以看到这样的设定
set palette defined (-0.035 "blue",-0.0075 "green", 0.02 "red")
如果要把填色的散点图与VMD绘制的填色的RDG等值面图相对照，则二者色彩刻度设定必须严格一致。比如在Multiwfn目录下的examples\RDGfill.vmd文件就是VMD里绘制填色等值面图的脚本，这里面mol scaleminmax top 1那一行后面的值应该设为-0.035 0.02才能与上图来对照（默认就是如此）。在这种色彩刻度下绘制的苯酚二聚体的RDG填色图如下所示，很明显散点图上各个spike位置和RDG填色图上的等值面通过颜色很容易进行一一对应。
(, 下载次数 Times of downloads: 122)
笔者在《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（http://sobereva.com/598）中介绍的我提出的IRI方法比RDG方法明显更强大，不仅可以展现弱相互作用区域，还可以展现化学键作用区域，因此强烈建议用IRI取代RDG。Multiwfn里IRI分析的界面和RDG分析如出一辙，可以导出IRI vs sign(lambda2)rho的散点数据到output.txt。之后用examples\scripts\IRIscatter.gnu（2022-Jul-16及以后更新的Multiwfn才有）代替上文的RDGscatter.gnu就可以绘制出填色的散点图。如果把脚本里的横坐标范围设大，比如设到-0.05到0.03范围，还可以使散点图把化学键作用区域的spike展现出来。例子看Multiwfn手册4.20.4节。
笔者提出的IGMH方法现在也特别流行，远比RDG更适合专门考察片段间的弱相互作用，在《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用》（http://sobereva.com/621）博文中有详细介绍。IGMH方法也可以绘制填色散点图，这在此博文里提到的IGMH官方教程里有很具体的例子。
作者Author: teller3531 时间: 2017-12-19 11:02
这个真漂亮
作者Author: youyno 时间: 2017-12-19 11:20
我都是手工用PS通过渐变画出来的
作者Author: w33 时间: 2017-12-19 21:22
gnuplot网站为什么一直500错误啊
作者Author: sobereva 时间: 2017-12-19 21:53
你的网有问题
作者Author: hljiang531 时间: 2017-12-24 20:39
你好，请问有填充颜色的散点图会画不？
作者Author: sobereva 时间: 2017-12-24 20:53
病句，意义不明
作者Author: indec 时间: 2019-3-28 04:10
社长您好，我使用文中方法产生一系列化合物的 RDG 散点 output文件，但是其中只有一部分能够产生完整的散点图，另一部分用 gnuplot 画图时只能产生 坐标轴和色彩刻度，散点都没有了，cmd 也没提示任何错误。照理说我的处理方法都一样，不应该一部分成功一部分不成功呀。这是什么情况，该怎么处理，谢谢
作者Author: sobereva 时间: 2019-3-28 04:26
我从没遇到过，不好说
作者Author: indec 时间: 2019-3-28 04:49
因为我的体系较大，分析弱作用的区域也较大，格点精度较高，有没有可能是因为点太多，程序处理不过来？
作者Author: sobereva 时间: 2019-3-28 07:06
你可以先试试点数较少的情况，如果点数多了就那样，也许是gnuplot的问题。
作者Author: lijiayisjtu 时间: 2019-4-8 19:44
老师，对于填色图绿色区域，是不是对应于结构中的范德华作用的总和？比如是否能知道粉色框中范德华作用对应于填色图的那一块区域呢，谢谢！！
作者Author: sobereva 时间: 2019-4-9 07:38
怎么考察散点图不同位置对应什么区域这里提了
用Multiwfn+VMD做RDG分析时的一些要点和常见问题
http://sobereva.com/291（http://bbs.keinsci.com/thread-1206-1-1.html）
把某些三维空间区域的点屏蔽掉再绘制散点图进行对比便知
作者Author: xxzj 时间: 2019-6-19 14:40
老师，您好，我想问一下，通常设计分子内非共价键相互作用时，通常指的O....S,H.......S等原子间的作用，那么氧原子是否与氮原子之间也存在弱相互作用呢？辛苦老师了
作者Author: sobereva 时间: 2019-6-20 04:04
弱相互作用无处不在，只要离得近，显然N和O也会有不可忽视的弱相互作用，只不过基本上肯定是以互斥作用为主（静电互斥）
作者Author: xxzj 时间: 2019-6-20 08:17
老师，Multiwfn中有没有方法可以证明N和O之间是互斥作用，还是这个是这个常识，那么我研究N和O之间相互作用与O和S间相互作用的差异是不是就没有意义啦，辛苦老师
作者Author: sobereva 时间: 2019-6-21 02:13
必然是互斥作用，因为通常二者肯定都带负电荷，而且也不会形成sigma-hole之类形式的作用。有化学常识的人都明白这一点
而且在实际中，由于N和O之间的静电互斥，优化后的结构里二者必然离得较远，肉眼上也能判断它们之间不会有强烈的相互作用
作者Author: sobereva 时间: 2019-8-24 04:01
从2019-Aug-24及之后更新的Multiwfn开始，根据手册3.23.1节的说明使用RDGscatter.gnu代替本文的RDGmap.gnu，就不需要在作图前手动对output.txt进行处理了，明显更省事了。
作者Author: snljty 时间: 2020-2-29 14:15
本帖最后由 snljty 于 2021-5-6 09:59 编辑
考虑到熟悉python的筒子们应该比熟悉gnuplot的多一些，写了一个作图的python脚本分享。
主要是matplotlib自带的一大堆色带方式居然没有BGR，显得很不方便，只好自定义了一个。试了几种手册上和网上推荐的方法，这种比较简洁。
需要先安装python3，Numpy库和Matplotlib库。
这个脚本用的是含有5列数据的原始的output.txt。
下面是脚本内容。
- #! /usr/bin/env python3
- # -*- Coding: UTF-8 -*-
- r"""
- Draw RDG map scatter.
- """
- import time
- import numpy as np
- import matplotlib.pyplot as plt
- from matplotlib.colors import LinearSegmentedColormap
- from sys import argv
- time_s = time.time()
- cmap_bgr = LinearSegmentedColormap.from_list("bgr", ["#0000FF", "#00FF00", "#FF0000"])
- sl2r, RDG = np.loadtxt(argv[1] if len(argv) > 1 else "output.txt", unpack = True, usecols = (3, 4))
- time_m = time.time()
- print("Time used for reading data: %5.1lf s" % (time_m - time_s))
- fig, ax = plt.subplots(figsize = (9.6, 7.2))
- sc = ax.scatter(sl2r, RDG, c = sl2r, s = 1, vmin = -0.035, vmax = 0.02, cmap = cmap_bgr)
- ax.set_xlim(-0.05, 0.05)
- ax.set_ylim(0, 2)
- ax.set_xticks(np.linspace(-0.05, 0.05, 11))
- ax.set_yticks(np.linspace(0, 2, 11))
- ax.set_xlabel("$\\mathrm{sign}\\left(\\lambda_2\\right)\\rho$ (a.u.)")
- ax.set_ylabel("RDG")
- cbar = plt.colorbar(sc)
- cbar.set_ticks(np.linspace(-0.035, 0.02, 12))
- # fig.savefig("RDGmap.ps")
- fig.savefig("RDGmap.png")
- time_e = time.time()
- # plt.show()
- print("Time used for drawing: %5.1lf s" % (time_e - time_m))
- print("Total time elapsed: %5.1lf s" % (time_e - time_s))
python还是慢了一些，50万个格点的体系差不多要十几秒，大体系可能要数分钟。请耐心等待。
下面是苯酚二聚体的结果。
(, 下载次数 Times of downloads: 115)
作者Author: 南窗月 时间: 2020-7-1 00:37
老师好，我按照文中方法得到的ps文件为什么零字符，用PS打开也是空白
作者Author: sobereva 时间: 2020-7-1 01:53
从你的描述无从判断。先把视频里的例子完整重复一遍。还不行换个机子再试
作者Author: 南窗月 时间: 2020-7-1 09:34
Terminal type is now 'wxt'
gnuplot> load 'RDGmap.gnu'
"RDGmap.gnu", line 19: warning: Cannot find or open file "output.txt"
"RDGmap.gnu", line 19: No data in plot
老师好，我又试了几次，会显示这个，我已经把output.txt文件放在bin文件夹里了，为什么说读不了呢
作者Author: sobereva 时间: 2020-7-1 10:38
南窗月 发表于 2020-7-1 09:34
Terminal type is now 'wxt'
gnuplot> load 'RDGmap.gnu'
"RDGmap.gnu", line 19: warning: C ...
直接在操作系统的命令行里调用gnuplot执行脚本多好，干嘛非要先进入gnuplot再载入脚本
作者Author: 南窗月 时间: 2020-7-1 15:37
好的，谢谢老师，问题已解决
作者Author: 不可理喻 时间: 2021-5-5 21:33
本帖最后由 不可理喻 于 2021-5-5 21:38 编辑
老师你好，我想请问一下为什么我做出的散点图的横坐标是sign(l_2)r (a.u.)，而不是跟视频里一样是sign(λ_2)ρ (a.u.)
作者Author: sobereva 时间: 2021-5-6 02:10
截图，并且说清楚用的是Multiwfn哪个版本，什么系统的
作者Author: 不可理喻 时间: 2021-5-6 08:47
本帖最后由 不可理喻 于 2021-5-6 08:48 编辑
电脑系统是win10的，Multiwfn用的是3.7版本，gnuplot用的是5.4.1版本，图是我用PS2020打开ps文件后导出的图片
(, 下载次数 Times of downloads: 131)
D:\gnuplot\gnuplot\bin\RDGscatter.png
作者Author: sobereva 时间: 2021-5-7 02:58
photoshop的事
用acrobat或者irfanview+ghostscript都能正常看到字符
作者Author: 不可理喻 时间: 2021-5-7 13:49
问题已解决，谢谢老师！
作者Author: Manyu 时间: 2022-6-29 11:38
请问为何我下载了ghostscript，还是无法用irfanview打开RDGmap.ps
作者Author: sobereva 时间: 2022-6-30 00:48
ghostscript没恰当安装，或者没被irfanview识别，或者你的ps文件本身就有毛病
作者Author: 20mllu1 时间: 2022-7-18 10:02
老师您好，这是我算IRI散点图的输出文件（图片1）
根据老师您放B站的教程视频，我也是用UE删了前面三列（图片2）
但是画出的图是这样子的（图片3），横坐标是从-0.05到0.05，而Multiwfn直接绘图是这样的（图4），横坐标从-0.4到0.1.
然后我用记事本改了RDGscatter.gnu中x轴范围图（5），画出的散点图依旧没有任何变化。
请问这种情况该如何解决呢？
作者Author: sobereva 时间: 2022-7-18 11:44
目前版本Multiwfn根本就不需要自己改txt文件，仔细看1L
最新版Multiwfn的examples\scripts目录里直接给了IRIscatter.gnu，直接用就完了
作者Author: 好运来 时间: 2022-7-23 09:55
您好，我也出现了这个问题，请问您解决了吗
作者Author: 好运来 时间: 2022-7-23 11:34
老师您好，我完全跟着您的方法用gnuplot生成的.ps文件，ghostscipt和irfanview装好之后里面自带的例子可以打开，但我自己生成的ps文件就打不开，我也觉得是文件的问题，该怎么解决呢
作者Author: sobereva 时间: 2022-7-23 14:33
没具体信息没法说
至少得提供你的output.txt，能让别人重复你的情况
作者Author: jingetiema6112 时间: 2024-7-19 20:46
卢老师您好，请教您2个问题，（1）用DFT-B3LYP/6-311G**优化邻硝基苯胺的结构，并且基于这个优化结构算了一系列性质。但是当研究邻硝基苯胺中分子内氢键这种弱相互作用时，看您的博文说应该考虑色散。像这种情况还能不能用DFT-B3LYP/6-311G**优化的这个结构啊？还是必须考虑色散，比如在B3LYP-D3(BJ)/6-311G**水平上重新优化，再基于这个重新优化的结构研究弱相互作用？
（2）仍然是用DFT-B3LYP/6-311G**优化邻硝基苯胺的结构，但是预测某些性质时，用到的公式是基于DFT-B3PW91/6-31G*得到的经验公式。像这种情况是不是也需要把邻硝基苯胺在DFT-B3PW91/6-31G*下重新优化，然后再用这个公式进行计算？
作者Author: sobereva 时间: 2024-7-23 03:36
1 拿不准就一律带着色散校正，下文明确说了
谈谈“计算时是否需要加DFT-D3色散校正？”
http://sobereva.com/413
2 严格来说是
作者Author: jingetiema6112 时间: 2024-7-23 15:37
好的，感谢卢老师指导
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |