使用Multiwfn结合VMD分析和绘制分子表面静电势分布

:笔者后来又写了《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》(http://sobereva.com/443),是绘制分子表面静电势图最完美的解决方案。其中通过脚本,把此文的绘制标准静电势填色的范德华表面图的步骤简化到了极致,只需要用不到本文1/10的步骤和耗时就可以达到比本文更好的效果。建议读者先阅读443这篇文章后再读本文,两篇文章有互补性。


使用Multiwfn结合VMD分析和绘制分子表面静电势分布

文/Sobereva @北京科音
First release: 2013-Jul-28   Last update: 2018-Aug-19

 

1 前言

分子表面静电势图经常在文献中出现,不同表面区域静电势大小通过不同颜色展现,使分子表面上静电势的分布一目了然。分子表面一般都用Bader定义的范德华表面,即电子密度为0.001 e/Bohr^3的等值面。只要提供所需的输入数据,这种图在许多程序中都可以作。本文介绍的通过Multiwfn的定量分子表面分析功能结合VMD和Photoshop作分子表面静电势图步骤相对步骤较多,但是优点十分显著,就是可以在分子表面上显示出静电势极值点位置,可以十分灵活地调节显示效果,而且还可以通过相同的方法绘制出分子表面上静电势以外的实空间函数的分布,比如平均局部离子化能、局部电子亲和能、Fukui函数等等。另外还可以顺带着获得许多其它信息,比如原核与分子表面的距离,实空间函数在分子表面上分布的统计数据等等。只要了解每一步操作的意义,就会觉得其实根本不复杂,而且思想上获得极大的解放。

本文所用Multiwfn版本为3.3.5,可在http://sobereva.com/multiwfn上下载。VMD为1.9版,可在http://www.ks.uiuc.edu/Research/vmd/上免费下载。Photoshop为CS2版。本文将以一个简单体系呋喃为例进行说明,波函数文件在B3LYP/6-31G**下由Gaussian产生。虽然本文主要是介绍作图方法,但也会顺带着说一下定量分子表面分析的一些功能。阅读本文前建议先看看《使用Multiwfn的定量分子表面分析功能预测反应位点、分析分子间相互作用》(http://sobereva.com/159

计算静电势是较为耗时的过程。如果你想显著降低Multiwfn计算静电势步骤上的耗时,特别是对于大体系、大基组而言,强烈建议通过借用Gaussian里的cubegen工具来实现,只需把settings.ini里的cubegenpath参数设为cubegen的实际路径即可。详见http://sobereva.com/435


2 在Multiwfn中计算、统计、导出数据

启动Multiwfn后依次输入
furan.wfn
12  //定量分子表面分析
0  //开始计算

默认就是在电子密度为0.001的表面上计算静电势,并且格点间距为0.25Bohr。格点间距可以通过选项3来设定,格点间距越小,分子表面就会用越多的顶点来描述,定量统计值、极值点位置也会越精确,之后作出的分子表面填色图的色彩过渡也越光滑,但是计算耗时将越长。如果你的目的仅仅是绘制分子表面静电势图,可以忽略以下内容而直接跳到本节最后一段。

计算完毕后会看到分子表面上静电势极大、极小点的坐标和数值,并且输出大量统计数据,比如分子表面上静电势的最大/最小值、平均值、方差、电荷平衡度、表面积等等,它们对于了解分子特征、建立QSPR/QSAR方程预测分子理化性质和生物活性等问题都十分有用,见Multiwfn手册3.15.1节以及IJQC,85,676、JPCA,110,1005等文章的介绍。

此时会看到后处理菜单。点击0进入图形界面,并且将Ratio of atomic size设为4.0就可以清楚看到分子表面上静电势极大点(红点)和极小点(蓝点)的位置。如果点击Minimal/Maximum label,就会显示出极值点的编号,可以和命令行窗口显示的极值点信息相对照,得知极值点上的具体数值。点击右上角Return可以关闭图形窗口。

此例中各个极值点静电势数值如下(由于格点精度有限,所以等价的极值点的数值不完全一致,这里取了平均)
极小点1(最小点):-20.60 kcal/mol   体现氧的孤对电子对静电势的负贡献
极小点2、3:-15.23kcal/mol   在两个beta碳(邻位碳)正上方,体现pi电子对静电势的负贡献
极大点1、4(最大点):18.41kcal/mol   体现alpha位的氢原子所带的正电
极大点5、6:15.21kcal/mol   体现beta位的氢原子所带的正电
极大点2、3:-9.63kcal/mol   虽然是极大点,但是静电势数值为负,所以化学意义不大,可无视之

我们可以查看一下不同静电势区间内分子表面积,这对于了解分子表面静电势的定量分布很有益。做法是依次输入
9
all
-25,22  //统计范围。当前体系分子表面静电势范围为-20.60~18.41kcal/mol,这里将范围稍微扩大并取整数来定义范围
15  //将-25~22kcal/mol均匀分为15个区间获得表面积
3  //输入的单位为kcal/mol
立刻屏幕上就输出了不同静电势区间内的表面积,我们将Center这一列和Area这一列的数据分别从屏幕上拷下来(不会拷的读者见手册5.4节),并粘贴到诸如origin等作图工具里并作成条形图,就可以得到诸如以下图像

静电势的定量分布在这张图上一目了然。此分子不同静电势区间内的表面积分布还是相对比较均匀的。

如果选择选项11,则程序会输出每个原子附近的分子表面上的静电势统计数据,对于了解原子在此分子中的特征很有用。例如静电势平均值
Atom#    All/Positive/Negative average
    1   -10.08997        NaN  -10.08997
    2   -12.59141        NaN  -12.59141
    3   -12.64250        NaN  -12.64250
    4   -10.04843        NaN  -10.04843
    5   -13.54925        NaN  -13.54925
    6     6.87357    9.80155   -3.13475
    7     4.84667    8.19922   -4.56994
    8     4.84069    8.21673   -4.50535
    9     6.88875    9.78376   -3.16678
氧原子(5号)附近的分子表面静电势平均值最负(-13.55kcal/mol),这也容易理解,毕竟其孤对电子对静电势有很大的负贡献。由于其附近分子表面上没有正值区域,所以正值区域的静电势平均值显示的是NaN。呋喃中碳原子裸露在分子表面上的区域主要就是体现pi电子特征的区域,由于pi电子云使这部分区域静电势为负,所以碳原子附近静电势平均值也都为明显负值,并且无正值区域。同时从静电势平均值上也体现出beta碳(2、3号)比alpha碳(1、4号)的pi电子云更富集因而静电势更负。这种讨论分子表面上对应不同原子区域的定量数据的方法和分子表面极值点分析往往会得到共通的结论,但此方法可以得到更多的信息,比如alpha碳上没有出现极值点(因此对它没有任何描述),而通过分析它对应的分子表面上的局部区域的静电势平均值等数据就可以定量考察它的特征。这种分析方法是笔者独创的,也仅有Multiwfn支持。另外,通过选项12还可以查看分子表面上对应于指定分子片段区域的定量性质。

选过选项9之后会问你是否输出locsurf.pdb文件,通过此文件可以利用VMD查看不同原子对应的分子表面区域。由于这不是本文的重点,所以输入n不让程序输出。

利用选项10可以查看分子表面相对于指定坐标点或指定原子核的最远和最近距离,这对于讨论分子间相互作用导致分子表面的穿透距离很有用。此功能也可以用于计算分子半径或直径,见《谈谈分子半径的计算和分子形状的描述》(http://sobereva.com/190)的讨论。

现在选择2将分子表面极值点导出到当前目录下surfanalysis.pdb文件中,此文件中碳和氧原子分别对应分子表面静电势极大点和极小点,pdb文件的B因子那一列的数据是静电势数值(kcal/mol)。然后再选择6将所有分子表面顶点导出到当前目录下vtx.pdb文件中,B因子是每个顶点的静电势数值。另外,如果你还没有当前分子的几何结构文件,应再选择5然后输入furan.pdb把当前体系的坐标导出到当前目录下furan.pdb中。


3 在VMD中作图

由于Multiwfn所用的图形库的限制,Multiwfn自身无法直接产生不同颜色填充的分子表面图,所以需要借助VMD来实现。下面的操作步骤对于当前体系比较适合,对于其它体系,应举一反三,根据实际情况进行略微调整。

启动VMD,然后把furan.pdb、surfanalysis.pdb和vtx.pdb按顺序依次拖动到VMD主窗口里(VMD Main),它们的ID将分别是0、1、2。选择Display-Depth Cueing将它关掉,否则图像会有些朦胧。进入Graphics-Colors,选Display-Background-White将背景改为白色,并且在此界面的Color-Scale标签页里选择BWR,使分子表面的色彩根据数值范围由小到大以 蓝-白-红 的方式变化。选Display-Axes-Off不让坐标轴显示出来。

注:如果希望每次启动VMD时都自动做如上操作,可以在VMD目录下vmd.rc文件末尾添加以下4
display depthcue off
color scale method BWR
color Display Background white
axes location Off

进入Graphics-Representations,然后执行以下3步来分别设定分子结构、分子表面极值点、分子表面顶点的显示方式。(如果只是想简单地看一下表面静电势分布,只要做第3步即可)
(1)在Selected Molecule一栏里选择furan.pdb,Drawing Method选Licorice,Bond Radius减小到0.1。

(2)将Selected Molecule一栏切换到surfanalysis.pdb,在控制台输入mol modstyle 0 1 VDW 0.06 (0和1是显示方式编号和体系的ID,此命令代表用大小为0.06的VDW球显示。由于GUI中VDW球最小只能设到0.1,故这里用命令行来实现)。然后在Selected Atoms里输入carbon并回车,然后将Coloring Method选为ColorID,并且在右边新出现的框里选Orange2。此时分子表面极大点就通过橙色小圆球显示出来了。点击Create Rep按钮创建新显示方式,在Selected Atoms里输入oxygen并回车,然后将ColorID右边的框设为Cyan,此时分子表面极小点就通过青色圆球显示出来了。

(3)在Selected Molecule一栏里选择vtx.pdb,Drawing Method选Points,Size设为25(设多大合决于视角的远近,在当前视角下应当让size恰好足够大,使分子表面上的顶点紧密相连,不留明显空隙),Coloring Method选Beta(根据pdb文件里B因子那一列的数据,此例即静电势数值进行填色),在Trajectory标签页里将Color Scale Data Range填上-22和22并点击Set,代表色彩刻度设为-22~20kcal/mol(其实此例用默认的色彩刻度范围就可以,这里只是为了取个整)。现在分子表面填色图就出现了。越蓝的区域静电势越负,越红的区域越正,白色区域的静电势数值在0附近。

之后给图上加上色彩刻度轴。选Extensions-Visualization-Color Scale Bar,Color bar width设为0.08,Display title选on并且将Color bar title里写上ESP (kcal/mol),Minimum和Maximum scale value分别填-22和22,Number of axis labels输入10,Color labels选Black,Label format选Decimal。然后点Draw Color Scale Bar按钮,色彩刻度就出现在画面中了,并且VMD Main窗口中多出了一个名为Color Scale Bar的一项。然后调整它的大小和位置,即双击VMD Main窗口中Color Scale Bar那一项当中的F标签使之变为红色(即不让色彩刻度轴在画面中的位置冻结),而双击其它项目的F标签使它们的F变为黑色(让它们的位置冻结住)。然后激活VMD图形窗口,按t键进入平移模式,然后拖动鼠标将色彩刻度轴放置到合适位置,并且用鼠标滚轮调整它的大小。调合适之后再按r键恢复旋转视角模式,并且在VMD Main里将Color Scale Bar那一项的F重新双击成黑色,而其它三项的F重新双击为红色。

当前的显示效果如下

此图还有许多地方值得进一步调整,也就是让分子结构显示出来(现在都被表面顶点掩盖了)、给分子表面极值点标上具体数值、让表面极值点更清楚地显示,特别是藏在分子表面后方的极值点也显示出来。为达到这些目的,需要利用photoshop。由于步骤比较细碎,笔者不把所有具体操作都叙述一遍,否则太罗嗦。只要会基本的photoshop操作的人都应该能理解应该怎么实现。


4 通过Photoshop改进作图效果