在GaussView里绘制分子表面静电势填色图的过程

在GaussView里绘制分子表面静电势填色图的过程

文/Sobereva(2)

First release: 2014-Sep-3   Last update: 2016-Jun-24


分子表面静电势填色图是计算化学研究中经常要绘制的一种图,在gview里的绘制方法非常简单,虽然操作步骤是最最最最最初级的问题,但是还是有初学者经常因此犯难并反复提问,此贴就简单说一下绘制流程。

PS:虽然gview绘制方便效果也好,但是如果想展现更多定量信息,更灵活,显得更高大上,就必须用Multiwfn结合VMD绘制,参见《使用Multiwfn结合VMD分析和绘制分子表面静电势分布》(http://sobereva.com/196)。另外,如果想让分子表面上静电势的色彩有明确的色阶以使不同区域静电势的定量差异显得更清晰,需要利用Molekel,做法见《谈谈Molekel做静电势填色等值面图的方法及误区》(http://sobereva.com/98)。


此文以多巴胺为例,这是输入文件
%chk=c:\gtest\dopamine.chk
#P B3LYP/6-31G* opt

dopamine

 0 1
O     -2.33158827    1.90324395    0.05204669
O     -3.45288968   -0.47124460    0.37321012
N      4.33709192    0.16551438    0.29362167
C      2.10661405   -0.10993743   -0.67966314
C      0.62079598   -0.23561182   -0.42292440
C      2.91795648    0.08557132    0.60357302
C     -0.18237422    0.90195216   -0.31692881
C      0.02151120   -1.46774074   -0.25587568
C     -1.53226977    0.81838896   -0.05200952
C     -1.34075114   -1.56264519    0.01270596
C     -2.12492987   -0.43572743    0.11734436
H      2.29355697    0.72769638   -1.34826040
H      2.47698171   -0.99578880   -1.18305280
H      2.54052562    0.95777703    1.14088437
H      2.76403205   -0.77011970    1.25210418
H      0.25802965    1.87739170   -0.45049059
H      0.60483520   -2.36786844   -0.33618849
H     -1.79443393   -2.53198365    0.13749855
H      4.88321312    0.19170936    1.13290303
H      4.54272968    1.00577810   -0.21284385
H     -1.82929975    2.69323340   -0.08852494
H     -3.74330646   -1.36791983    0.46123396

 

算完之后,用gview打开dopamine.chk,选Results-Surfaces/Contours。点Cube Actions-New Cube,然后设成下面这样

对于HF、DFT计算,Density Matrix都是选SCF。假设你用的是MP2方法,就得选MP2。Grid一栏控制格点精度,直接影响图像质量。如果体系较大,为了节省时间应当选Coarse。而选择Fine则意义不大,因为Medium时给出的图像质量已经足够好了。默认是不计算内层电子的贡献的,如果选上Use full density matrix则也会计算内层电子贡献。对一般分子选不选Use full density matrix对结果没什么影响,若为了保险和结果精确可以选上,而对于一些含碱、碱土金属的分子如果不选上可能结果定性错误。

点Ok后,gview开始调用Gaussian自带的cubegen程序产生电子密度的格点数据,因此屏幕上出现一个命令行窗口。等这个窗口自动关闭了,说明算完了,并在Cube Available一栏里出现了电子密度格点数据信息。

范德华表面有不同定义,最常用的是Bader对范德华表面的定义,也就是电子密度为0.001的等值面。因此把窗口中的Density右边的值改为0.001。然后点击Surface Actions按钮,出现两个选项:New Surface代表产生当前格点数据的等值面,此时也就是绘制出分子范德华表面。New Mapped Surface则不仅产生前格点数据的等值面,还把指定的函数在这个等值面上的值以不同颜色映射到这个表面上。对于当前目的,我们应当选New Mapped Surface,然后看到新窗口

什么都不用改,默认的映射到等值面的函数就是ESP(静电势),直接选Ok即可(对于MP2计算,此时Density Matrix也得选MP2,否则是基于HF密度计算静电势,没那么准确)。

之后会看到图形窗口左下角显示Generating mesh...。这个计算过程会先产生电子密度0.001等值面的顶点,这些顶点一起描述了分子范德华表面,然后会计算静电势在每个顶点上的函数值,之后图形窗口里就出现了静电势填色等值面。

这样的图形不太好考察,每个原子位置都看不清楚,因此做一些调整。在图上点右键,选择View-Display Format,在General标签页里点击Background Color右边的那个按钮,选白色。然后进入Surface标签页,Format选择Transparent。然后会看到一个Transparent滑条,越往右边滑会看到静电势的填色效果越浓。

色彩刻度轴在图中上方,默认情况下,静电势从-0.08054变化到0.08054对应于颜色从红色变化到蓝色(遗憾的是gview里色彩过渡方式没法自己选择,和笔者的习惯相反,笔者喜欢蓝色代表负静电势而红色代表正静电势)。这个色彩刻度范围是程序自己设的,建议设成一个比较整的数以便于标注,如-0.08到0.08。如果把色彩刻度上下限调小,比如改为-0.06到0.06,那么静电势很正或很负的区域的颜色就会变得更为浓郁。

最终效果就是本文开头的图,颜色鲜明,同时不挡着原子。

已有 14 条评论

  1. mr

    其实Gview中有一选项,叫transparency,选50%就不错了

  2. merlot

    Sob老师,由于我的笔记本电脑打不开chk文件,我尝试了使用fchk按照上述步骤想要得到静电势填色图,但是得到的确实所谓的静电势等值面图,两头显示的数值是e^(-3),我想知道错在什么地方?我仔细对比的操作,唯一的不同应该就是输入文件不同,但是fchk不是包含了chk的信息么?

    1. sobereva

      gview载入chk的话,会自动调用formchk转换成fchk之后再载入,和你直接载入fchk是一码事。
      如果得到了静电势等值面图,肯定是你在gv里操作不对。第一步你得选电子密度,这样得到的表面肯定是分子表面。第二步你再选静电势,静电势就会投影上去。

      1. merlot

        找到问题了,从Cube Actions-New Cube后,我没选择total density 而是ESP,才导致后续问题,真是一步错步步错。谢谢老师分享和指正。

      2. Yoyo

        请问老师,new cube里的选项都是灰色的,不能选择total density,是哪里出现问题了吗

  3. 狒狒

    我算了一个水分子的例子,结构优化完成后,打开chk文件,按照这个教程,最后的时候提示Warning, Surface map failed! Opening temporary array file,这是怎么回事呢

    最后显示的好象是表面电子密度,不是静电势,谢谢

  4. 学者

    您好,我用GaussView,根本打不开chk文件,是不是应该fchk文件,但我用fchk文件,计算出就一种颜色,怎么回事?

    1. sobereva

      就一种颜色要么是作图过程不对(比如显示等值面时你选的是静电势,而不是电子密度),要么就是色彩刻度设得不合理,没法将不同静电势区域以颜色区分开。
      gview打不开chk,要么是你机子里没在安装Gaussian(因为gview会先调用Gaussian里的formchk工具自动转换为fch再载入),要么就是虽然装了Gaussian,但是和你产生这个chk文件用的Gaussian不是同一个版本。

  5. Justin

    请问分子表面某一点的坐标和静电势的对应关系可以输出吗

    1. sobereva

      这要用Multiwfn,可以得到分子表面上各个顶点的具体静电势值(Multiwfn输出的vtx.pdb里就有),见
      使用Multiwfn结合VMD分析和绘制分子表面静电势分布
      http://sobereva.com/196

  6. 北落师门

    请问一下,在已知晶体结构,有cif文件的情况下,是不是可以不用计算chk文件直接计算静电势呢?应该怎么做?非常感谢!

    1. sobereva

      如果你就是想研究分子在晶体结构中的构象的静电势,把分子坐标从cif里提出来,冻结重原子优化氢原子(晶体的氢原子位置严重不准),然后再产生.fch或.wfn文件用Multiwfn或gview绘制。

  7. lairylil

    您好,sob老师,我们计算多用WFA程序做静电势分析,有没有关于这个程序的绘制静电势的方法呀

    1. sobereva

      用WFA干嘛,WFA能做的用我开发的Multiwfn都能做,而且更灵活、功能更多、做得更好。

添加新评论