Using GaussView element coloring scheme in VMD

文/Sobereva@北京科音   2022-Sep-16

GaussView的元素配色总体上来说不错，而且对所有元素都定义了唯一的颜色，在File - Preference - Colors - Element Colors中可以看到所有元素的颜色，如下所示

VMD是免费、灵活、强大而且超级流行的化学体系可视化工具。在VMD里也可以按照元素着色，见《在VMD程序里对不同元素的原子用不同颜色显示的方法》（//www.umsyar.com/624）。然而，至少对于笔者撰文时的最新的正式版1.9.3版来说，VMD自带的元素着色定义很少，对绝大多数元素都是统一用的褐色，导致很多情况下没法区分不同元素，而且颜色也不美观。这里介绍在VMD里使用GaussView里元素颜色定义的做法。

从//www.umsyar.com/attach/652/gview_color.tcl下载gview_color.tcl文件，将之放到VMD目录下（即VMD启动后在文本窗口里输入pwd命令后显示的目录），然后在VMD的vmd.rc文件末尾加入一行proc gview {} {source gview_color.tcl}。如果你不了解vmd.rc文件的话，看《VMD初始化文件(vmd.rc)我的推荐设置》（//www.umsyar.com/545）。

启动VMD后，在文本窗口里输入gview，就可以把默认的元素着色方案替换成和GaussView相同的了。之后载入个含有元素信息的结构文件（比如pdb、xyz等），然后在Graphics - Representation里把Coloring Method设为Element即可看到按元素着色的效果。如果想和GaussView显示的图像特征尽可能接近，将Drawing Method设为CPK，Bond Radius设为0.2，Material设为AOShiny。对顺铂体系显示效果如下，左边是VMD的，右边是GaussView的

上图还可以看到在键的着色方面有差异，VMD的键的两边的颜色分别对应两个原子的，而GaussView都是白色的。如果想让VMD在这点和GaussView统一，可以把当前的Representation里的Bond Radius设为0使得键不显示，然后再建立一个Representation，用CPK，把Sphere Scale设0，Bond Radius设0.15，Coloring Method选ColorID并指定为白色，之后效果如下所示，可见和GaussView显示的很接近了

再把原子半径、光源方向、材质微调一下，就更像GaussView的效果了。不过，GaussView里对多重键的显示效果是VMD里怎么设都模仿不来的。

The way to display atoms of different elements in different colors in VMD program

文/Sobereva@北京科音

First release: 2021-Nov-11  Last update: 2022-Sep-16

1 前言

之前笔者写了大量的将Multiwfn与VMD相结合的文章，比如《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（//www.umsyar.com/598）、《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》（//www.umsyar.com/443）、《用VMD绘制艺术级轨道等值面图的方法》（//www.umsyar.com/449）等等，并且做法已被大量同行所使用，令VMD被 的研究者们用得也越来普遍，并有越来越多的人在思想家公社QQ群和计算化学公社论坛里问我VMD的使用问题。其中有一个问题被问得越来越频繁，已经成了半周经问题，也就是怎么对不同元素的原子用不同颜色显示。为了免得老得重复回答，特此写一个小文章专门说一下这事。本文的情况是针对VMD 1.9.3版而言的，其它版本的情况可能相同也可能不同。

笔者还有一篇文章介绍了怎么让VMD对元素的着色和另一个流行的可视化程序GaussView相同，建议阅读，见《在VMD中使用GaussView的元素着色的方法》（//www.umsyar.com/652）。

2 name和element属性

VMD里每个原子都有各种属性。name和element在VMD里是每个原子的两种不同的属性，name是原子的名字，element是元素周期表里的元素符号（第一个字母大写，第二个小写）。在VMD的文本窗口中，输入[atomselect top all] get name命令可以列出当前体系里所有原子的name，输入[atomselect top all] get element命令可以列出当前体系里所有原子的element，都是按原子序号顺序输出。

3 不同格式的文件记录的信息

不同输入文件提供给VMD的信息是不一样的，这里对几种典型的文件来说一下。

• xyz文件：介绍见《谈谈记录化学体系结构的xyz文件》（//www.umsyar.com/477）。在标准的xyz文件中，第一列是每个原子的元素符号。将xyz文件载入VMD后，name和element属性值都是元素符号。

注意有些xyz文件不标准，里面记录的是原子名而非元素符号，比如VMD载入gro文件后直接保存出来的xyz文件里就是原子名。这样的xyz文件载入VMD后name属性值是文件里记录的原子名，element会根据原子名里面的字母部分去猜。比如原子名是HG3的话会被认作汞，因此element被设为Hg。如果猜不出来，比如原子名是CG，则element属性将为X。注意VMD这么猜元素很容易猜错，比如alpha碳在氨基酸里的标准原子名是CA，如果xyz文件里记录的是原子名，VMD就会把它视为是钙，从而令element属性值为Ca。

• gro文件：即GROMACS的结构文件。里面对每个原子记录的是原子名，比如氨基酸里面的碳可以有CA、CB、CG、CG1等等原子名。载入VMD后，name就是这些原子名，而element皆为X。

• mol2文件：同上。

• pdb文件：在标准的pdb文件中，第三列是原子名，最后一列是元素符号。因此，pdb文件载入VMD后可以给VMD分别提供name和element信息，因此element不用根据原子名猜了。但是有些pdb文件不规矩，没有记录元素符号的最后一列，这样的pdb文件载入后element就都为X。

• cub文件：cub文件经常通过VMD程序绘制成等值面图，见比如《在VMD里将cube文件瞬间绘制成效果极佳的等值面图的方法》（//www.umsyar.com/483）。如《Gaussian型cube文件简介及读、写方法和简单应用》（//www.umsyar.com/125）所述，此文件里没记录原子名，但记录了各个原子的元素在周期表里的序号。此文件载入VMD后，name和element属性都为元素符号。

还值得一提的是VMD里每个原子还有个type属性，记录的是原子类型，可以用[atomselect top all] get type命令查询。诸如amber的拓扑文件prmtop里记录了原子类型，mol2格式也记录了原子类型（但很多程序输出的mol2文件里原子类型和元素符号相同），因此载入后type属性就是文件里记录的原子类型。而xyz、gro、pdb格式里没记录原子类型，载入后type属性值会和name相同。cub格式里也没记录原子类型，载入后type、name和element都相同。

4 对原子的着色设置

这里以氯仿为例进行说明，此分子的pdb文件可以在//www.umsyar.com/attach/624/CHCl3.pdb下载。

将它载入VMD后，在Graphics - Representation里把Drawing method改为CPK，目前看到的是下图

可见，碳和氯的颜色完全一样，都是青色。为什么没有区分开？原因很简单，如上图所见，Coloring method目前是默认的Name，即根据原子的name属性决定颜色。我们来看一下当前着色用的颜色定义。进入Graphics - Colors，在Categories里选Name，点击C，然后看到下图

如图可见，当前根据name进行着色时并不区分C和Cl，只要首字母是C就都按用cyan颜色着色。

如果想按照元素来着色，就在Graphics - Representation里把Coloring method改为Element，此时图像如下，可见氯和碳的颜色区分开了

氯用棕色明显不好看。我们要把它改为常用的绿色，就进入Graphics - Colors，在Categories里选Element，点击Cl，再选green，此时如下图所示

如果想微调green颜色的具体色彩定义，可以在上图的Color Definitions标签页下方修改红、绿、蓝三种颜色分量的大小。

记住，只有载入的是pdb、xyz这样能直接给VMD提供元素信息的输入文件，才能像上面这样按照element属性来着色。而用比如gro文件的时候，由于不能提供元素信息，就只能按照name来着色了。如果你是GROMACS用户又想按元素着色，就得把gro转成比如pdb格式，按照pdb格式规范在合适的列上补上元素信息后再载入；或者用比如gmx editconf -f md.tpr -o md.pdb把tpr文件转成pdb文件，由于tpr文件里本身有元素信息，所以这样得到的pdb文件里最后一列直接记录了元素信息（顺带一提，之后你可以删除当前仅有的一帧，然后再往这个体系中载入gro文件或轨迹文件，此时载入的只有坐标，而元素信息还是载入pdb时提供给VMD的）。

5 设置默认着色方式

如果你想将element作为默认着色方式，并且Cl元素默认用绿色，免得每次在VMD里还得如上操作一遍，可以在vmd.rc文件末尾加入以下两行
mol default color Element
color Element Cl green
每次启动VMD后就会自动执行这两条命令修改默认设置。如果不了解vmd.rc文件的话，看《VMD初始化文件(vmd.rc)我的推荐设置》（//www.umsyar.com/545）。

另外，如果你想默认用CPK方式显示体系结构，在vmd.rc末尾还可以再加上mol default style CPK。

还值得一提的是VMD有个colordefs.dat文件，记录了默认的着色用的规则。对于Windows版，此文件在VMD目录下的scripts\vmd目录下；对于Linux版，如果安装到了默认目录，此文件在/usr/local/lib/vmd/scripts/vmd目录下。用文本编辑器打开此文件可以看到

Name   H   white
Name   O   red
Name   N   blue
Name   C   cyan
Name   S   yellow
Name   P   tan
Name   Z   silver
...略
Element C  cyan
Element Ca ochre
Element Cd ochre
Element Ce ochre
Element Cf ochre
Element Cl ochre
Element Cm ochre
...略

你可以直接在colordefs.dat里修改按照name着色、按照element着色时默认的色彩，比如可以把Element Cl ochre改为Element Cl green使得程序按照element着色时对Cl元素默认用绿色。这里设置的优先级低于自己在vmd.rc末尾添加诸如color Element Cl green这样的设置命令。

On the image clarity issues in graphing related to computational chemistry

文/Sobereva@北京科音  2021-Jul-16

0 前言

寡人在网上答疑的时候，总是看到有人问“得到的图像不清晰怎么办”、“怎么让图像更清晰”这种问题，语义不明，也不知道对方说的“清晰”到底指什么，对方也不说他觉得不清晰是怎么个不清晰法，也不给出个截图来说明，令人感觉很莫名其妙。笔者忍不住写个小文，把各种初学者所谓的“不清晰”的可能情况都罗列一下，便于他们自己对照判断如何解决。此文也有助于一些缺乏相关常识的计算化学工作者了解怎么改进图像质量。此文内容很多地方用流行的可视化程序VMD（http://www.ks.uiuc.edu/Research/vmd/）来举例说明，这里是针对VMD 1.9.3版而言的。

1 图片分辨率太低

关于图片分辨率的事我在《谈谈怎么正确认识论文投稿时对图像分辨率的要求》（//www.umsyar.com/511）里专门做过科普。有些人所说的不清晰可能是指图片的分辨率（当前语境下等于图片的像素）太低。显然，让产生图片的程序输出更高分辨率的图像就完了。比如：
(1)Multiwfn程序（//www.umsyar.com/multiwfn）保存一维、二维、三维图像的分辨率分别通过settings.ini文件里的graph1Dsize、graph2Dsize、graph3Dsize参数来设置。
(2)GaussView通过File - Save Image方式保存图片的像素由界面上的Enlarge Width and Height by控制。比如如果这里设的是3x，而窗口的大小是400*300，那么保存出的图像就是1200*900像素。
(3)VMD可以通过File - Render方式保存图片，如果渲染器选的是Snapshot或者Tachyon (internal, in-memory rendering)，那么渲染得到的图片和OpenGL窗口的大小是相同的。如果想要更高分辨率的图，要么把窗口拉大再渲染，要么在用第三方渲染器POV-Ray或VMD自带的独立的Tachyon渲染器渲染的时候直接指定渲染出的图片像素。用Tachyon渲染器的时候如何通过命令行指定渲染出的图片分辨率在《用Multiwfn+VMD做RDG分析时的一些要点和常见问题》（//www.umsyar.com/291）里我专门说过，这里不再累述。

有些可视化程序没有单独的保存图片的功能，你只能靠截屏来得到图片文件。此时你若要想得到高像素的图像，就应当把窗口最大化，并且把窗口里的体系也放大到尽可能大，然后再截图。如果还嫌图片像素不够高，去找个高分辨率屏幕的机子，在那上面截图。

2 过度有损压缩图像

十分常用的jpg是一种有损压缩的图像格式。Photoshop、IrfanView等像样的图片编辑或者图片观看程序在保存图像成jpg格式时，都可以选择jpg文件的质量；设得越高，则图像质量相对于原图降低得越少，但图像文件就越大。显然，如果保存时图片质量设得太低，必然得到的图像会比较模糊，并且还可能伴有色彩失真。

对于科研方面的图，特别是可能用到文章里的图，我平时都存成png格式，这是非常好的无损压缩格式，不会像jpg格式那样降低图片质量（代价是尺寸更大），又比bmp那种未经压缩的无损格式小得多得多，也比常见的tif那种无损压缩格式尺寸更小。

3 图像被拉大了

有些人在保存/生成图片的时候图像的像素并不是很高，然而他把图插入到文章里之后，还手动把图的尺寸给拉大了，这时候图片能不模糊么？只有ps、eps、wmf等格式储存的矢量图才能无损缩放，而png、bmp、tif、jpg等位图格式只要放大了图片注定会看起来模糊，这属于图像处理的最基本常识。

4 一些程序自动压缩了图片

有些人问，怎么图片原本看起来是清晰的，插入XXX程序之后就不清晰了。这看具体情况。比如powerpoint（笔者这里用的是2016），在“选项”-“高级”的“图像大小和质量”一栏中，有个选项是“不压缩文件中的图像”，这个选项默认是没有选中的。这种情况下，在保存pptx文件时，对于ppi（和《谈谈怎么正确认识论文投稿时对图像分辨率的要求》//www.umsyar.com/511里说的dpi本质上相同）过高的图像，其像素就会自动压缩到特定的ppi。如果你把一个图片插入到powerpoint里，并且把图片尺寸缩小，然后保存pptx文件，此时由于图片的ppi往往会高于阈值，遂被powerpoint自动压缩，即当前pptx里的这个图片的像素已经低于原图了。之后重新打开这个pptx文件，若再把这张图拉回原先的大小，就会看到图片已经模糊了。只有把“不压缩文件中的图像”选上才可以避免powerpoint自动压缩图像。

还有其它的程序可能也有根据ppi自动压缩图片的设置，应仔细去选项菜单里仔细看看，弄不明白的话可以Google搜搜看有无别人遇到过并提供了解决方案。

5 直接截图，而不是保存图像

有些可视化程序为了显示速度比较快，直接在屏幕上显示的时候用略低的图像质量，而在保存图片的时候才自动设高图像质量。因此，可视化程序凡是有专门的保存图片的功能的话就不要直接截图。

对于Multiwfn显示的三维分子结构图、等值面图，特别是使用mesh或者solid face+mesh风格显示等值面的时候，尤其建议保存成高像素的图片文件，然后再缩成想要的尺寸，这比起直接在三维窗口里看到的图像质量好太多了，还顺带起到抗锯齿效果（见后文）。

6 保存图像的格式不合适

要分清楚矢量图和位图的区别。如果你的图片没有平滑的色彩过渡，而是以线条、文字和离散的颜色构成的，那么保存成矢量图比较好，可以无损缩放，而且线条看起来很平滑。比如在Multiwfn里按照《使用Multiwfn绘制红外、拉曼、UV-Vis、ECD、VCD和ROA光谱图》（//www.umsyar.com/224）和《使用Multiwfn绘制NMR谱》
（//www.umsyar.com/565）绘制光谱图、按照Multiwfn手册4.4节的例子绘制等值线图，我都建议以矢量图形式保存为pdf格式的文件。如果图片有连续的色彩过渡，比如Multiwfn、VMD、ChimeraX等程序显示的等值面图，以及Multiwfn主功能4绘制的填色图，则建议用png等位图格式。

如果违背上面的建议，把线条、文字和离散的颜色构成的图保存成png等格式，线条看起来可能就不如保存成矢量格式时看着那么平滑。如果把色彩连续变化的图保存成矢量图，图片看起来就会有很难看的色阶。

7 被投稿系统降低了图像质量

在大多数期刊投稿时，在最后确认提交之前，一般投稿系统会生成一个用于预览的pdf文件。转换成pdf文件的时候会对图片自动进行处理，例如把过高分辨率的图像自动降低分辨率、适当对图片进行有损压缩从而避免pdf尺寸太大，这可能导致在预览pdf文件时发现图片清晰度变低了，有点糊了。这一般没办法避免，也不用太在意，一般也不至于影响审稿。这时候的图片清晰度和最终出版时候的图片清晰度没直接关系。

8 没用抗锯齿

显示3D物体时，如果没在可视化程序里开启抗锯齿，而且显示器的dpi又不是太高，就会明显看到3D物体边缘的线条毛毛糙糙的，可能有的人说的不清晰是指这个。很多可视化程序里都有抗锯齿（antialiasing）选项，打开之后就能让边缘锯齿感减轻很多，变得平滑不少。有些程序里抗锯齿有不同级别，越高级别抗锯齿效果越好，但代价就是对GPU性能要求越高。有些可视化程序里可能没有抗锯齿选项，但也没关系，你可以进入显卡驱动，往往在里面可以设置强行对某个程序用某种抗锯齿设定。比如VMD的抗锯齿选项（Display - Antialiasing）往往没法直接选，因此你可以在显卡驱动里直接设，具体怎么操作和你的显卡的显示芯片制造商以及驱动版本有关，自行摸索或Google。如果你的显卡驱动是Windows自带的，那么是没有显卡驱动面板的，应当自行去显示芯片的厂家官网上下载并装上。

还有一种人工对图片实现抗锯齿的做法是保存很高像素的图像，然后用Photoshop、IrfanView等程序把图片缩小（降低像素），此时程序会做重新采样，可以等效实现抗锯齿效果。

9 没有关闭雾化效果

很多程序默认都开启了雾化效果来体现物体距离屏幕的位置。比如在VMD里，默认是黑背景，并且开着雾化效果（即Display - Depth cueing是打开的）。此时距离屏幕越远的物体就显得越黑。但由于VMD对雾化的默认衰减设置不太理想，所以即便是距离屏幕较近的物体看起来也是雾蒙蒙的，在《用Multiwfn+VMD做RDG分析时的一些要点和常见问题》（//www.umsyar.com/291）的一开始就有个示意图。有些初学者由于缺乏常识，就把雾化效果理解成了不清晰。我建议用VMD的人要么关闭Depth cueing，要么在Display - Display Settings里恰当设置Cue Mode，使雾化效果至少不会影响到距离屏幕较近的物体，而只对中、远距离的物体才生效。

10 绘制等值面图时格点间距太大

不管什么可视化程序，在绘制某个函数的等值面图之前都需要获得这个函数在特定一块三维空间中均匀分布的各个格点上的函数值，这叫格点数据。格点间距越大，得到的等值面就越粗糙、越有棱有角，格点间距很大时甚至有的地方的等值面可能还不连续，也许有些人说的不清晰就是指这种情况。格点间距如何影响RDG函数的等值面质量在《用Multiwfn+VMD做RDG分析时的一些要点和常见问题》（//www.umsyar.com/291）的2.1节做了展现。显然，等值面效果不好的问题只要把格点间距设小一些即可缓解，在诸如Multiwfn程序里提供了非常丰富灵活的方式设定格点。但无疑减小格点间距的代价是计算的总格点数变多了，耗时增高了。

11 光源设置问题

显示3D物体时，如果光源设置不当，可能导致有些地方太暗，令有的人觉得不清晰。解决方法就是调节光源设置。比如在VMD 1.9.3里有4个光源，默认开启前两个，可以在Display里选择相应的光源来切换打开或关闭的状态。另外，在Mouse - Move Light里选择一个光源后，可以在图形窗口里拖动鼠标来调节光源位置。而在Multiwfn里，可以在窗口上方的菜单里选择Other settings - Set lightings，然后选择哪些光源打开或关闭。

12 显示风格设得不合理

大多数可视化程序里都可以设置不同的显示效果，如果设置不当可能会造成看起来不清晰。比如VMD 1.9.3程序默认用的是line方式显示，而且粗细只有1.0，因此直接显示出的体系看起来太细，如果是没有成键的单个原子，图形窗口里会显示成只有一个像素大小的点，视力不好的同志不仅根本看不清楚，还甚至会以为是VMD显示有问题，漏掉了原子，笔者在网上答疑时就见过这样的人。显然把显示方式改成诸如CPK、Licorice之类的就能看清楚了。笔者在《VMD初始化文件(vmd.rc)我的推荐设置》（//www.umsyar.com/545）中推荐的默认设置中把默认的line风格的线条粗细设大到了2.0也是为了让默认的显示效果看着更清楚。

A simple example of performing ab initio dynamics (AIMD) using ORCA

Sobereva@北京科音

First release: 2020-Dec-13  Last update: 2021-Jul-7

1 前言

基于 方法的动力学一般称为从头算动力学（Ab initio molecular dynamics, AIMD），相比于基于一般的经典力场的动力学，其关键优势在于精度高、普适性强、能够描述化学反应，代价是耗时相差N个数量级。ORCA 程序有不错的做BOMD形式的从头算动力学的功能，使用很方便，而且本身ORCA做DFT的效率又高，是做孤立体系AIMD的首选程序之一。虽然有些特性不支持，比如没法像Gaussian的BOMD那样直接做准经典动力学，不能根据原子距离等标准判断什么时候自动结束任务等等，但都不是大问题。对于跑跑普通的AIMD来说，笔者感觉ORCA明显比Gaussian的ADMP或BOMD更好用（Gaussian的AIMD输入文件较为抽象，手册相应部分写得很烂，而且连个像样的热浴都没有），而且速度也明显更快。

笔者发表的18碳环（cyclo[18]carbon）的研究论文中，其单体和二聚体做的分子动力学就是用ORCA跑的，分别见《揭示各种新奇的碳环体系的振动特征》（//www.umsyar.com/578）对中Chem. Asian J., 16, 56 (2021)和《全面探究18碳环独特的分子间相互作用与pi-pi堆积特征》（//www.umsyar.com/572）中对Carbon, 171, 514 (2021)一文的介绍。在《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（//www.umsyar.com/696）中介绍的笔者的Inorg. Chem., 62, 19986 (2023)文章中还用ORCA跑了B6C6N6的高温动力学以证明其稳定性。在《18个氮原子组成的环状分子长什么样？一篇文章全面揭示18氮环的特征！》（//www.umsyar.com/696）中介绍的笔者的ChemPhysChem, 25, e202400377 (2024)文章中还用ORCA跑了18氮环的动力学以考察其热分解行为。这些文章可以作为ORCA跑AIMD的典型范例，很推荐阅读和引用。

笔者在北京科音高级 培训班（http://www.keinsci.com/workshop/KAQC_content.html）中会用多达两百多页的ppt专门深入详细讲AIMD的模拟，其中也包括ORCA做AIMD的各种细节、大量技巧以及诸多实例，并且此培训中还会系统讲授ORCA程序的使用，因此是使用ORCA专门做AIMD的读者一定不可错过的培训。而本文只是举一个简单的例子，帮助读者快速了解ORCA如何做最简单的AIMD。注意ORCA只适合跑孤立体系的从头算动力学，如果是做周期性体系的从头算动力学，CP2K是最佳的选择，在笔者讲授的北京科音CP2K第一性原理计算培训班（http://www.keinsci.com/workshop/KFP_content.html）中有极为全面、系统、详细的讲解并给了大量例子。

本文内容适用于Multiwfn最新版本、VMD 1.9.3、ORCA 4.2.x及以后版本的情况。

2021-Jul-7 针对ORCA 5.0的补充说明：对于2021-Jul-7及以后更新的Multiwfn，进入本文所述的Multiwfn产生输入文件的界面后，可以通过选项-11选择适合的ORCA版本，默认为ORCA 5.0及以后版本，而下文内容对应的是4.2.x版的情况。对于>=5.0版，Multiwfn自动设的热浴是CSVR，比之前版本唯一支持的Berendsen热浴更好，而且同样普适。而且在run后面多加了CenterCOM，这是从5.0版本开始支持的消除整体质心运动的选项，而下文里提到的constraint add center这一行就没有了。

2 实例：[Al(H2O)6]3+与NH3之间的质子转移

ORCA的MD模块的开发者网站上有个真空中[Al(H2O)6]3+与NH3之间的质子转移的动画，如下所示

可见NH3向带正电的[Al(H2O)6]3+逐渐靠近，水上的一个质子转移到了氨气分子上，然后由于静电互斥，NH4+就逐渐远离[Al(H2O)5(OH)]2+了。这里我们试图重现这个过程。下面提到的文件都可以在//www.umsyar.com/attach/576/file.zip中获得。

我们首先获得[Al(H2O)6]3+的基本合理的结构。当然用ORCA优化也可以，这里笔者习惯性地用Gaussian来优化。在GaussView里搭建Al(H2O)6，保存为Al(H2O)6_optfreq.gjf，将关键词改为# B3LYP/6-31G* opt freq，将电荷改为3，然后用Gaussian运行之，就得到了优化后的[Al(H2O)6]3+结构。再用GaussView打开输出文件Al(H2O)6_optfreq.out，把一个氨气分子画在一个水的旁边，如下所示，然后保存为complex.gjf。

Multiwfn程序（//www.umsyar.com/multiwfn）有很便利的生成ORCA常见任务的输入文件的功能，见《详谈Multiwfn产生ORCA 程序的输入文件的功能》（//www.umsyar.com/490），这里我们用Multiwfn生成ORCA的AIMD任务的标准输入文件。

启动Multiwfn，载入complex.gjf，然后输入
oi  // 生成ORCA输入文件
0  // 选择任务类型
6  // 分子动力学
1  // 计算级别用B97-3c
在当前目录下就得到了ORCA的AIMD任务的标准输入文件complex.inp。

在complex.inp里面将相应几行改成下面这样：
%maxcore  10000
%pal nprocs   10 end
timestep 1.0_fs
initvel 298.15_K
* xyz   3   1

现在的complex.inp的完整内容如下。以#为开头的行代表后面的设置被注释掉了，不会生效，想启用可以去掉#。此文件里也有大量Multiwfn自动添加的注释，只要一看注释马上就明白相应的行是什么含义，巨贴心，都省得查手册了。

! B97-3c noautostart miniprint nopop
%maxcore  10000
%pal nprocs   10 end
%md
#restart ifexists  # Continue MD by reading [basename].mdrestart if it exists. In this case "initvel" should be commented
#minimize  # Do minimization prior to MD simulation
 timestep 1.0_fs  # This stepsize is safe at several hundreds of Kelvin
 initvel 298.15_K no_overwrite # Assign velocity according to temperature for atoms whose velocities are not available
 thermostat berendsen 298.15_K timecon 30.0_fs  # Target temperature and coupling time constant
 dump position stride 1 format xyz filename "pos.xyz"  # Dump position every "stride" steps
#dump force stride 1 format xyz filename "force.xyz"  # Dump force every "stride" steps
#dump velocity stride 1 format xyz filename "vel.xyz"  # Dump velocity every "stride" steps
#dump gbw stride 20 filename "wfn"  # Dump wavefunction to "wfn[step].gbw" files every "stride" steps
 constraint add center 0..22  # Fix center of mass at initial position
 run 2000  # Number of MD steps
end
* xyz   3   1
[坐标部分]
*

下面简单说一下complex.inp里这些设置的含义、为什么这么设。
• B97-3c是一个又便宜又不错的计算级别，在ORCA里还自动会启用RIJ加速，速度很快，因此很适合跑AIMD，描述当前体系没有问题。B97-3c的介绍见《详谈Multiwfn产生ORCA 程序的输入文件的功能》（//www.umsyar.com/490）的相应部分。但B97-3c也不是什么时候都能用，比如18碳环用纯泛函描述皆失败，见笔者在Carbon, 165, 468 (2020) 里的讨论，显然就不能用这个了，笔者跑涉及18碳环的AIMD的时候都用的是ωB97X-D3。
• %pal nprocs   10 end代表用10核。笔者当前任务实际上是在双路E5-2696v3共36个物理核心的机子上跑的，但却故意用了10核。这是因为根据笔者以前的测试，发现ORCA做DFT的AIMD的并行效率不理想，尤其是对于小体系，用核数太多甚至反倒速度更慢（大家可以对实际情况短暂跑比如5步实测一下设多少核速度最快）。一般来说就设10核就行了，当机子里有明显更多核的时候，可以跑多个AIMD任务来充分利用计算资源，但应当对CPU内核进行绑定，否则AIMD计算速度可能显著降低，见《通过设置CPU内核绑定降低ORCA同时做多任务的耗时》（//www.umsyar.com/553）。
• %maxcore  10000代表每个ORCA进程最多用10000MB。其实完全没必要这么大，普通泛函下的AIMD不怎么耗内存，给1000都绝对够了。
• restart ifexists这句被注释掉了。如果你的AIMD想续跑，且当前目录下有之前跑出来的与当前任务同名的后缀为.mdrestart的文件，可以取消注释，任务就会延续之前AIMD最后的状态续跑，新轨迹会在原有轨迹文件后面续写。
• minimize这句被注释掉了。如果取消注释的话，动力学开始前会自动在当前级别下做几何优化。
• timestep 1.0_fs代表动力学步长用1 fs。对于此例常温下的模拟，1 fs步长是合适的，改大的话可能造成动力学不稳定，而改小的话跑同样的时间长度需要更多步数导致需要更多耗时。如果追求绝对稳妥，或者是在明显更高温度下模拟，可以用0.5 fs。
• initvel 298.15_K代表根据298.15 K温度通过Maxwell分布随机初始化原子速度。no_overwrite代表如果之前已经有初速度信息了（比如可以是续跑时从之前的mdrestart文件里继承来的），就不再产生新的初速度。
• thermostat berendsen 298.15_K timecon 30.0_fs代表用Berendsen热浴控温在298.15 K（此例笔者用的ORCA 4.2.1只能用这个热浴，据悉从ORCA 5.0开始可以用更好的velocity rescale热浴），时间常数为30 fs，通常这个时间常数是合适的。
• dump position stride 1 format xyz filename "pos.xyz"代表每一步都往当前目录下的pos.xyz文件里写入当前的坐标，因此pos.xyz是多帧xyz格式的轨迹文件。此格式的介绍见《谈谈记录化学体系结构的xyz文件》（//www.umsyar.com/477）。
• dump force和dump velocity开头的行都被注释掉了，这两行用于把模拟过程中的原子受力和原子速度分别保存到相应xyz文件里。dump gbw开头的行也被注释掉了，它可以在MD过程中每隔指定步数保存一次gbw文件，之后可以通过写脚本调用Multiwfn对它们进行批量分析，从而考察动力学过程中电子结构变化，得到丰富的有化学意义的信息（如动力学中的电荷转移情况、成键变化情况、电子定域性变化等等，在北京科音高级 培训班里笔者会给出很多这种例子）
• constraint add center 0..22代表将当前整个体系（0号到22号原子）的质心约束在初始位置。之所以这么做，是因为尽管ORCA在产生初速度时已经去掉了整体平移的速度分量，但实际模拟过程中由于数值问题，仍然往往会看到有明显的整体漂移的现象，因而有碍观察（在VMD里观看轨迹时还得做一下align才能消掉）。因此直接把质心位置约束住就没有这个问题了。
• run 2000代表总共跑2000步，当前步长是1 fs，即最多跑2 ps。当前这个模拟的目的是观察到质子转移，跑多长时间合适并不好预估，所以可以一次先跑2000步看看，若不够到时候再续跑。值得一提的是，至少对于笔者现在用的ORCA 4.2.1，我发现如果一次跑的步数很多，达到2000步左右之后，之后每一步的耗时会有逐渐的上升趋势，原因不明。因此我建议每次跑最多不宜超过3000步，如果需要跑更长的话，最好分多段来跑。

现在用ORCA运行complex.inp，模拟过程中可以看到每一步的实时情况：

Step是当前的步数，Time是当前的时间，t_SCF和t_Grad分别是算这一步的SCF过程和计算受力的耗时，二者加和就是这一步的总耗时。可见每一步耗时约6~7秒，乘以要跑的步数就可以估计跑完整个任务的耗时。后面还可以看到每一步的温度、动能、势能等信息。由于当前体系原子数很少，温度相对于热浴的参考温度波动大是很正常的事。而且由于用了热浴，所以总能量E_tot也有明显波动。

VMD是观看动力学轨迹最好的程序，可以在http://www.ks.uiuc.edu/Research/vmd/免费下载。建议在模拟过程中，隔一阵子就用VMD把pos.xyz载入进去，看看当前的动态行为如何、跑成什么结构了。跑到289步的时候，笔者看了一眼pos.xyz，发现质子不仅已经完全转移，而且NH4+都已经跑走一定距离了，所以就没必要继续跑了，故把ORCA直接杀掉了。

模拟过程中ORCA还产生其它一些文件。complex.md.log是ORCA的MD模块输出的信息，相当于整个输出文件中的子集。complex.mdrestart是用来续跑的文件，每一步都会往里面写入当前步的时间、坐标、速度、受力等信息。complex-md-ener.csv是把每一步的时间、耗时、温度、能量等信息以csv格式保存的文件。还有其它一些零碎的文件，通常不是一般用户关心的，这里就不说了，都可以放心删掉，留着也没用。

现在我们来看模拟轨迹。建议大家根据《VMD初始化文件(vmd.rc)我的推荐设置》（//www.umsyar.com/545）里说的修改VMD的初始化文件，从而添加自定义命令bt，这样在播放轨迹的时候对每一帧会自动重新判断成键关系。

启动VMD，载入pos.xyz（在本文文件包里已提供，共290帧）。然后在文本窗口输入bt，按回车，使得每帧都更新连接关系。然后再输入bw，按回车使得背景变为白色。选Graphics - Representation，将Drawing Method改为CPK。然后点击VMD界面右下角的三角播放动画，会看到如下结果。如果想在VMD图形窗口中显示出帧号或时间，看《使VMD播放轨迹时同步显示帧号》（//www.umsyar.com/13）。

可见，我们跑出来的动力学现象和本文一开始的那个官方动画几乎完全一致。都是两个反应物先接近，然后形成复合状态时质子在二者之间震荡，最后NH4+跑掉。

如果想把质子转移情况通过距离随时间的变化曲线方式呈献给读者，可以在点击VMD的三维图形窗口后，按键盘上的2键（之后按r可以恢复为默认的旋转模式），然后点击两个原子正中心，二者之间就会增加Bond label（默认是以白色字显示距离，在黑背景下才看得清楚），这里笔者把N和转移过去的质子之间增加了Bond label。然后进入Graphics - Labels，切换到Bonds，选Graph，点击Show preview复选框，然后点击1:H 1:N这项，就会看到距离变化已经显示在预览窗口了，如下所示，其中红色和蓝色竖条标记的分别是最小距离和最大距离位置和数值。

如果点击Graph按钮，就会把曲线显示在大窗口中，如下所示。可见，在大约60帧，也即60 fs左右，N-H键就算是基本形成了，之后N-H键不断振动。

以类似方式，我们可以标记Al与N的距离，随时间变化如下所示。可见Al与N先接近，质子转移完毕后，二者就逐渐远离了。

在Labels窗口里还可以点击save，把距离变化数据保存到文本文件里，之后可以导入Origin等程序里绘制曲线。

用VMD还可以测量角度、二面角的变化，分别是在图形窗口里按3然后点三个原子、按4然后点四个原子进行标记，之后在Graphics - Labels里观看。

3 总结&其它

通过上面的例子，可以看到ORCA做AIMD是相当容易的，只要把Multiwfn支持的含有结构信息的文件（如pdb/gjf/xyz/mol/mol2/fch等等，见Multiwfn手册2.5节）载入Multiwfn，敲几下键盘产生标准AIMD任务的输入文件，然后根据实际情况稍微改几个设定就可以跑了。

以几十核的一般双路服务器的运算能力，ORCA里用B97-3c跑几十原子有机体系的几十ps的动力学不是特别困难的事。不过，能跑的时间尺度仍远远比不上xtb跑半经验层面DFT的GFN-xTB方法的动力学，xtb跑动力学的粗浅介绍和例子看此文的相应部分：《使用Molclus结合xtb做的动力学模拟对瑞德西韦(Remdesivir)做构象搜索》（http://bbs.keinsci.com/thread-16255-1-1.html）。因此，拿ORCA跑DFT的动力学之前，先拿xtb初步跑跑，找找感觉，大体摸索出自己期望的现象能出现的条件（如温度、初始结构、反应物相对位置和碰撞方式等），然后再用DFT跑通常是比较好的做法，免得做昂贵的DFT的MD试来试去把时间都耽误了。

本文只涉及了VMD一丁点皮毛，VMD对于做动力学的人是必须玩得非常溜的。笔者在北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/workshop/KGMX_content.html）里对VMD有非常深入全面的介绍，包括tcl脚本的编写。利用VMD的tcl脚本可以对轨迹做千变万化的分析，有些分析诸如质心距离变化、平面间夹角变化、某几何变量分布统计、不同结构出现比率等，是必须靠写脚本才能实现的。

光是分析分析动力学过程的能量、结构变化是很肤浅的。利用Multiwfn，可以对ORCA跑的动力学的全过程的电子结构做深入透彻的分析，从而考察化学键、弱相互作用、电荷分布等等在动力学过程中的变化，由此能够从提供深入的视角，使研究文章信息更丰富、明显更上档次。非常建议详细阅读《详谈Multiwfn的命令行方式运行和批量运行的方法》（//www.umsyar.com/612），里面第4节专门讲了怎么做这样的分析，你会发现特别容易也特别有价值。

如果Multiwfn创建ORCA做动力学输入文件的功能对你的实际研究产生了帮助，希望在写文章的时候提及诸如The input file of ab-initio molecular dynamics was prepared with the help of Multiwfn program并引用Multiwfn原文，这是对Multiwfn此功能开发最好的支持。

Displaying atomic forces calculated by Gaussian in VMD

文/Sobereva@北京科音  2020-Sep-13

在 研究中，原子受力往往是非常值得关注的量。Gaussian程序可以通过force关键词计算原子的受力，然后在GaussView中可以根据原子受力大小进行着色，但是没法绘制出受力矢量。本文介绍如何在免费、灵活、强大的VMD程序中，通过笔者自写的tcl脚本，将Gaussian计算的原子受力非常直观地通过箭头展现出来，这对于一些 问题的研究很有益处。VMD可以在http://www.ks.uiuc.edu/Research/vmd/免费下载，本文使用VMD 1.9.3，Gaussian使用G16 A.03。

在此处下载绘图脚本和示例文件：//www.umsyar.com/attach/568/file.zip。将其中的forcevec.tcl和drawarrow.tcl复制到VMD目录下。

用文本编辑器打开forcevec.tcl，可看到开头有几个设置：
set filename：载入的Gaussian输出文件路径
set sclfac：以a.u.为单位的原子受力矢量乘上这个系数就是绘制出的箭头的长度。应根据实际图像效果反复尝试找到最合适的
set rad：箭头的粗细
set color：箭头的颜色名
set crit：如果某个原子的受力大小小于所有原子最大受力大小乘以这个系数，则这个原子的受力就不会用箭头绘制出来，由此可避免在受力相对非常小的原子上出现非常短的、没什么实际意义的箭头。默认为0.1

例1：联苯的基态极小点结构下在S1势能面上的原子受力

本文文件包里的biphenyl_S1force.out是在优化过的联苯的基态极小点结构下，用TDDFT算的S1激发态的原子受力的输出文件，用到的关键词在此文件里已经体现了。

把biphenyl_S1force.out载入GaussView，保存成pdb文件，然后把pdb文件载入到VMD里。然后将biphenyl_S1force.out挪到VMD目录下，用文本编辑器打开forcevec.tcl，把set filename后面的内容改为biphenyl_S1force.out，把set sclfac后面的内容改为20，然后保存文件。之后在VMD的文本窗口里输入source forcevec.tcl，脚本就从biphenyl_S1force.out里读取原子受力，然后绘制出了箭头。把背景改成白色，在Graphics - Representation里把显示方式改为Licorice后看到下图

根据原子受力可见，在体系所处的这个Franck-condon点位置，联苯中央的C-C键倾向于缩短。实际上在PBE0/6-311G*级别下，联苯的基态极小点结构下这个C-C键键长为1.478埃，而同级别下用TDDFT优化的S1态极小点下这个键长为1.413埃，明显短了很多。另外，从上图可见其它原子也有显著受力，这体现出苯环结构将要自发进行调整，确实S1态极小点下苯环的各个键长相对于S0态都有显著变化。

例2：晕苯与富勒烯复合物

笔者之前在PM6-D3下对晕苯与富勒烯形成的复合物进行了优化。这里看一下如果在优化的复合物结构的基础上，若把富勒烯与晕苯距离稍微拉远一点，原子受力是什么样的。本文文件包里complex_force.out就是在人为拉远的结构下用PM6-D3做force任务的输出文件，complex_force.pdb是相应的结构文件。将complex_force.pdb载入VMD，把complex_force.out复制到VMD目录下，将forcevec.tcl里的set filename后面写上complex_force.out。由于弱相互作用对应的原子间相互作用很弱，所以在set sclfac后面写上一个比较大的值，这里用1500。在VMD文本窗口里运行source forcevec.tcl，进一步调节下作图设置后看到下图

可见，由于富勒烯与晕苯之间的色散吸引作用，再加上当前结构下二者间距比平衡距离更远，两个分子彼此间挨得较近的原子上都出现了令二者距离进一步拉近的力。

例3：外电场下的18碳环的受力

笔者对18碳环及类似体系做过大量的理论研究，汇总见//www.umsyar.com/carbon_ring.html。其中，笔者在名为Ultrastrong Regulation Effect of the Electric Field on the All-Carboatomic Ring Cyclo[18]Carbon（https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/cphc.202000903）的论文中全面、深入研究了外电场对18碳环的几何结构、电子结构、电子吸收光谱等特征的影响，此文的介绍见《一篇文章深入揭示外电场对18碳环的超强调控作用》（//www.umsyar.com/570）。在此文中通过类似于上文的做法绘制了下面的图（需要利用专门的tcl脚本）

此图体现的信息的详细解释请看上面提到的论文，这里只是简单提一下。这个图绘制了各个原子的受力，还用绿色箭头展现了左、右两侧原子的总体受力，用紫色箭头展现了上、下两侧原子的总体受力，划分方式在图(a)中通过虚线体现了。原子的颜色体现的是ADCH方法算的原子电荷，实现方法在《使用Multiwfn+VMD以原子着色方式表现原子电荷、自旋布居、电荷转移、简缩福井函数》（//www.umsyar.com/425）中介绍了。图(a)体现出18碳环在原始结构下，如果刚加上外电场，此刻的受力并不会立刻造成碳环的整体变形，但会明显造成键长发生改变的倾向。图(b)体现出，在外电场下当碳环的结构已稍微经过弛豫、出现了一些变形后，此刻的原子受力就会明显倾向于让18碳环在顺着电场方向拉长，而在垂直于电场方向压扁，从而使碳环变得更加像椭圆。图(c)是在电场下优化后的结构，但是此刻撤掉了外电场，由绿色箭头可见此时18碳环有明显的恢复成原本圆形的倾向，体现出弹性特征。此图清晰地展现出电场下18碳环这个体系的独特的力学特性，如果不把受力绘制成箭头，很难予以这样深入的探讨。

若要绘制一批原子的总受力，可以在source forcevec.tcl之后，把下面的代码拷到VMD文本窗口里。这部分代码原本是笔者用来绘制上图中左侧5个原子总受力用的。4 5 6 7 8是被计算总受力的原子序号，总受力矢量及其模也会被输出到屏幕上。在绘制箭头时，为了让位置比较舒服，箭头中心位置取的是$atmleft 3 9 2 10的几何中心，即4 5 6 7 8加上相邻的3 9 2 10四个原子。

set atmleft "4 5 6 7 8"
draw color green
set fxleft 0
set fyleft 0
set fzleft 0
foreach i $atmleft {
set fxleft [expr $fxleft+$fx($i)]
set fyleft [expr $fyleft+$fy($i)]
set fzleft [expr $fzleft+$fz($i)]
}
puts "Left side force: $fxleft $fyleft $fzleft"
set normval [expr sqrt($fxleft**2+$fyleft**2+$fzleft**2)]
puts "Norm: $normval"
drawarrow "serial $atmleft 3 9 2 10" $fxleft $fyleft $fzleft $sclfac 0.2

Method for plotting molecular vibrational vectors calculated by Gaussian in VMD

文/Sobereva@北京科音

 First release: 2020-Sep-8    Last update: 2024-Feb-4

Gaussian用户观看freq任务产生的振动矢量一般都是通过GaussView看（虽然也有ChemCraft等其它一些程序也可以看）。然而，起码对于GaussView 6来说，GaussView显示振动矢量的一个很大不足是箭头太细，而且头部不够粗，导致有时候都看不清楚，放在文章里不够美观。另外，GaussView绘制分子结构的作图选项不够灵活，而且还收费。VMD是极其流行的化学体系可视化程序，免费、灵活、图像效果好，本文介绍如何通过笔者写的VMD作图脚本非常方便地绘制Gaussian的振动分析任务产生的振动矢量。VMD可以在http://www.ks.uiuc.edu/Research/vmd/免费下载。

在这里下载笔者编写的绘图脚本和示例文件：//www.umsyar.com/attach/567/file.zip。此脚本至少对于目前撰文时的VMD正式版中最新的1.9.3、Gaussian 09和16是完全适用的。

这里以绘制多巴胺的振动矢量为例进行演示。把文件包里的dopamine.out放到VMD目录下，这是多巴胺的Gaussian的freq任务的输出文件。然后我们得把这个.out文件转化成一个VMD可以认的结构文件的格式，比如可以把此文件载入GaussView，然后另存为.pdb或.mol2文件。也可以下载Multiwfn（//www.umsyar.com/multiwfn），启动Multiwfn后载入此文件，然后选主功能100的子功能2，通过相应选项导出为.pdb或.xyz文件。

把文件包里的drawarrow.tcl和GauNorm.tcl都放到VMD目录下，然后用文本编辑器打开GauNorm.tcl，把开头的set filename后面的文件名改为dopamine.out。之后启动VMD，把多巴胺的结构文件载入VMD，然后在文本窗口输入source GauNorm.tcl执行此脚本，此时振动矢量信息就被读入了，与此同时定义了名为norm的绘制振动矢量的命令。之后在VMD的文本窗口输入比如norm 4，就可以把4号振动模式通过箭头画出来。

norm后面还可以接第2个参数，用来设置箭头长度是正则矢量的几倍，数值越大箭头越长，默认是3。norm后面还可以接第3个参数，用来设置箭头的半径，默认为0.05。比如norm 5 6 0.07就代表用6倍长度、0.07的半径绘制第5个振动矢量。默认是用黄色绘制箭头，如果想用别的颜色，把GauNorm.tcl中的draw color后面的yellow改成其它颜色名，比如cyan。

此例输入norm 17 5，然后令分子以CPK方式显示（在Graphics - Representation里把Drawing method改为CPK，再把Sphere Scale设为0.6），效果如下，可见非常理想！和GaussView显示的相对比，可见展现的信息是相同的，而GaussView画的箭头相比之下明显太小气了。

注意GauNorm.tcl开头还有个set ilinear语句，如果当前体系是线型体系，必须把后面的值改为1。

如果你在Gaussian做freq或opt freq任务中按照《在Gaussian中做限制性优化的方法》（//www.umsyar.com/404）中的做法将N个原子的笛卡尔坐标冻结了，运行source GauNorm.tcl之前必须把里面set nfreeze后面的值设为N（默认为0，没有原子被冻结）。

如果你希望让箭头的始端位于各个原子上（和GaussView的风格一致），就把本文文件包里的drawarrow2.tcl放到VMD目录下，把GauNorm.tcl里的两处drawarrow都改为drawarrow2并保存。之后再按上文绘图即可。用norm 17 4命令，效果如下

Calculation of dipole moment of molecular fragments and dipole moment of monomers in complexes using Multiwfn

文/Sobereva@北京科音  2020-Jun-23

1 原理

时不时有人问怎么计算一个分子中某些片段的偶极矩，或者问怎么计算一个分子复合物中各个分子的偶极矩，确实这对于分析局部的极性很有意义，在本文就通过例子对做法进行说明。做这种分析必须使用2020-Jun-22之后更新的Multiwfn程序，可以在//www.umsyar.com/multiwfn免费下载，相关基本知识看《Multiwfn FAQ》（//www.umsyar.com/452）和《Multiwfn入门tips》（//www.umsyar.com/167）。此功能需要使用含有波函数信息的文件作为输入文件，比如wfn、wfx、mwfn、fch、molden等，产生的方式见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（//www.umsyar.com/379）。

简单来说，计算一个体系中由某些原子构成的片段F的偶极矩矢量就是计算

其中R是原子核位置，Z是原子核电荷，ρ是电子密度，w是原子的权重函数。计算片段偶极矩矢量需要用Multiwfn的主功能15（模糊空间分析），目前支持Hirshfeld、Becke、Hirshfeld-I形式的原子权重函数。本文就用比较常用，而且物理意义比较清楚的Hirshfeld权重函数来算。

需要注意的是，如果某个片段的净电荷不为零，则这个片段的偶极矩是有原点依赖性的，即平移体系会令它发生改变。多数情况下被计算的片段的净电荷不为0，却也没有办法完全避免原点依赖性，只能是在讨论的时候留意这个问题，予以恰当讨论。

笔者之前还写过一篇与本文有一定关系的文章《使用Multiwfn+VMD绘制片段贡献的跃迁偶极矩矢量》（//www.umsyar.com/396），感兴趣者推荐看看，这对于考察体系的不同区域对电子激发的振子强度的影响很有益。

2 计算分子片段的偶极矩实例

本节计算下面这个单重态[Ru(bpy)3]2+阳离子体系的各个配体的偶极矩，其wfn文件可以在这里下载：//www.umsyar.com/attach/558/Ru_bpy_3.zip。此文件的结构使用BP86泛函结合SDD和6-31G*进行优化，波函数用B3LYP结合SDD和6-311G*在IEFPCM模型描述的水环境下产生。此体系分子结构如下：

启动Multiwfn后载入Ru_bpy_3.wfn，然后依次输入
15  //模糊空间分析
-1  //选择原子权重函数
3   //基于内置原子密度的Hirshfeld权重函数
-5  //定义功能1、2计算时考虑的原子范围
2-21  //体系中一个配体里的原子序号
2  //计算原子和分子多极矩

此时Multiwfn依次计算各个原子的多极矩，最后给出这些原子对应的片段的电子数、偶极矩和不同形式的多极矩：

             *****  Molecular dipole and multipole moments  *****
Total number of electrons:     81.427849
Molecular dipole moment (a.u.):      -0.000000      4.348703      0.000000
Molecular dipole moment (Debye):     -0.000000     11.053300      0.000000
Magnitude of molecular dipole moment (a.u.&Debye):      4.348703     11.053300
Molecular quadrupole moments (Standard Cartesian form):
XX=  -41.258381  XY=   -0.000000  XZ=  -15.828189
YX=   -0.000000  YY=  -13.049129  YZ=   -0.000000
ZX=  -15.828189  ZY=   -0.000000  ZZ=  -33.042376
Molecular quadrupole moments (Traceless Cartesian form):
XX=  -18.212629  XY=   -0.000000  XZ=  -23.742283
YX=   -0.000000  YY=   24.101250  YZ=   -0.000000
ZX=  -23.742283  ZY=   -0.000000  ZZ=   -5.888621
Magnitude of the traceless quadrupole moment tensor:   25.129610
Molecular quadrupole moments (Spherical harmonic form):
Q_2,0 =  -5.888621   Q_2,-1=  -0.000000   Q_2,1= -27.415227
Q_2,-2=  -0.000000   Q_2,2 = -24.429930
Magnitude: |Q_2|=   37.189945
Molecular octopole moments (Cartesian form):
XXX=    0.0000  YYY= -455.8273  ZZZ=   -0.0000  XYY=   -0.0000  XXY= -217.0765
XXZ=    0.0000  XZZ=    0.0000  YZZ= -176.7460  YYZ=    0.0000  XYZ=  -85.0737
Molecular octopole moments (Spherical harmonic form):
Q_3,0 =    -0.0000  Q_3,-1=   -20.8698  Q_3,1 =     0.0000
Q_3,-2=  -329.4892  Q_3,2 =    -0.0000  Q_3,-3=  -154.4790  Q_3,3 =     0.0000
Magnitude: |Q_3|=    364.5030

这些量的具体定义看Multiwfn手册3.18.3节，这里不细说。我们当前关心的只有
Molecular dipole moment (a.u.):      -0.000000      4.348703      0.000000
这一行，即偶极矩矢量为(0.0,4.3487,0.0) a.u.。

然后我们计算其它片段的偶极矩，依次输入
-5  //重新定义片段
42-61
2  //计算原子和分子多极矩
-5  //重新定义片段
22-41
2  //计算原子和分子多极矩

汇总一下，三个配体的偶极矩分别为：
2-21片段： ( 0.000000  4.348703  0.000000) a.u.
42-61片段：( 3.766463 -2.173911  0.001788) a.u.
22-41片段：(-3.766463 -2.173911 -0.001788) a.u.

下面，我们用VMD程序将偶极矩的箭头画出来以便于观看。VMD可以在http://www.ks.uiuc.edu/Research/vmd/免费下载，本文用的是VMD 1.9.3。wfn格式没法载入VMD，所以我们先把当前体系转换为VMD可以载入的xyz格式。在当前Multiwfn窗口中输入
0  //返回主菜单
100  //其它功能(Part 1)
2   //导出文件
2   //导出为xyz格式
直接按回车用默认路径，然后Ru_bpy_3.xyz就产生在当前目录下了。

将Ru_bpy_3.xyz载入VMD，然后将Multiwfn的examples\scripts目录下的drawarrow.tcl里的全部内容复制到VMD文本窗口里，这样就定义了一个新的名为drawarrow的命令，我们用这个命令来绘制片段偶极矩。

在VMD的文本窗口里运行以下命令
draw color green
drawarrow "serial 2 to 21" 0.000000 4.348703 0.000000 0.7
draw color red
drawarrow "serial 42 to 61" 3.766463 -2.173911  0.001788 0.7
draw color yellow
drawarrow "serial 22 to 41" -3.766463 -2.173911 -0.001788 0.7

这些命令代表分别用绿、红、黄色箭头把三个片段偶极矩绘制出来。双引号里面是选择相应片段里的原子的选择语句，语法详见《VMD里原子选择语句的语法和例子》（//www.umsyar.com/504）。之后的三个参数是偶极矩的X/Y/Z分量，最后的参数0.7是在绘图时把偶极矩数值乘上这个系数再绘制，起到控制箭头长度、令图像美观的目的。绘制出来的图像如下。分子的显示方式可以在Graphics - Representation界面调节。

偶极矩是从负电中心指向正电中心。由图可见每个片段的始端都是两个氮的方向，这是因为氮带负电，而且有丰富的孤对电子。显然当前的结果完全合理。

如果之后想删掉箭头，在文本窗口运行draw delete all命令。

值得一提的是，当前这个体系里三个片段的偶极矩是对称分布的，这是因为中心金属正好在坐标原点上。由于每个配体的净电荷不为零，因此如果在做产生波函数的任务之前把体系平移一下，最终得到的三个箭头就不以Ru为中心对称了。所以要注意原点位置问题。

有读者问如果片段里的原子序号是不连续的怎么输入。假设片段里原子序号是4,9-15,89-100,111，上面例子里的双引号里就写serial 4 9 to 15 89 to 100 111，即把-替换为[空格]to[空格]，把逗号替换为[空格]，即可。

3 计算分子复合物中单体偶极矩实例

下面这个例子我们计算一下苯酚二聚体中每个苯酚单体的偶极矩，并且连同二聚体的总偶极矩一起用VMD画出来。

启动Multiwfn，然后输入
examples\phenoldimer.wfn  //自带的苯酚二聚体波函数文件
15  //模糊空间分析
-1  //选择原子权重函数
3   //Hirshfeld
2  //计算原子和分子多极矩
由于我们当前还没定义片段，所以Multiwfn最后给出的下面的信息直接就是二聚体的偶极矩了
Molecular dipole moment (a.u.):       1.227306     -0.128087      0.650833

之后和上一节的例子一样，对第一个苯酚（1~13号原子）和第二个苯酚（14~26号原子）分别计算片段偶极矩，结果分别为(0.570356 -0.356257 0.284025) a.u.和(0.656950 0.228171 0.366808) a.u.。这两个苯酚的偶极矩的矢量和与前面输出的二聚体的偶极矩是相同的。

和前例一样，把当前结构导出成xyz文件并载入VMD，并且把drawarrow.tcl里的内容拷到VMD的文本窗口里运行，之后输入以下命令进行绘制：
draw color green
drawarrow all 1.227306 -0.128087 0.650833 2
draw color red
drawarrow "serial 1 to 13" 0.570356 -0.356257 0.284025 2
draw color yellow
drawarrow "serial 14 to 26" 0.656950 0.228171 0.366808 2

这些语句代表二聚体的偶极矩用绿色箭头绘制，all代表整个体系。两个苯酚的偶极矩分别用红色和黄色绘制。由于苯酚的偶极矩较小，为了让图中箭头足够长，所以命令末尾用了2来让箭头长度变成偶极矩大小的两倍。

得到的图像如下

可见由于两个苯酚的偶极矩的方向基本相同，二者的矢量和所对应的二聚体的偶极矩比单体偶极矩大得多。

My recommended settings of VMD initialization file (vmd.rc)

文/Sobereva@北京科音

 First release: 2020-Apr-1   Last update: 2023-Apr-4

VMD（http://www.ks.uiuc.edu/Research/vmd/）是用得极为广泛的化学体系可视化程序，由于其极度灵活，有很多技巧可以使其用起来更方便。

VMD启动时会先用初始化文件对一些设置进行初始化，即执行里面的各种命令，用户也可以往里添加额外的命令。对于Windows版来说，这个文件就是VMD目录下的vmd.rc。对于Linux版来说，这个文件叫.vmdrc，VMD会先在当前目录下搜索，没有的话就去找~/.vmdrc，还没有的话去找$VMDDIR/.vmdrc（这里$VMDDIR环境变量是没有预先定义的），如果还找不到此文件，就会用默认设置。
注：Linux下的.vmdrc文件默认出现在安装目录下，比如以默认路径安装会出现为/usr/local/lib/vmd/.vmdrc，但如果不配置VMDDIR环境变量的话这个文件并不会被VMD在启动时自动载入。Linux下.vmdrc一般都是自行在用户主目录下创建。

在此我将我自己的初始化文件里的设置进行分享，其中额外添加的内容如下（放到原有内容后头即可）。下文的叙述是对撰文时最新版本VMD 1.9.3而言的。

mol default style {Lines 2.0}
display depthcue off
#color Display Background white
#axes location Off
display rendermode GLSL
display distance -8.0

proc bw {} {color Display Background white}
proc bb {} {color Display Background black}

user add key Right {animate next}
user add key Left {animate prev}
user add key Up {animate goto [expr $vmd_frame([molinfo top])+10]}
user add key Down {animate goto [expr $vmd_frame([molinfo top])-10]}
user add key b {mol bondsrecalc all; topo retypebonds}

proc bt {} {
global vmd_frame
trace variable vmd_frame([molinfo top]) w updatebond
}
proc updatebond {args} {
mol bondsrecalc all
topo retypebonds
}
proc bn {} {
global vmd_frame
trace vdelete vmd_frame([molinfo top]) w updatebond
}

proc fog {} {
display depthcue on
display cuemode Linear
display cuestart 1.75
display cueend 2.5
}

下面解释一下做这些设置有什么好处。

程序默认的显示方式是Lines，但是线的粗细太细，往往看不清楚，所以用mol default style {Lines 2.0}将默认的显示方式改为两倍粗细的Lines。

程序默认开着雾化，即让距离镜头越远的物体的颜色混入越多的背景色。这会导致在黑色背景下物体的颜色显得不够鲜艳，而在白色背景下物体又显得有点雾蒙蒙，因此用display depthcue off将雾化效果关掉。

#color Display Background white这行是被注释掉的。如果你想让VMD启动后默认就用白背景，就把#去掉。

#axes location Off这行也是被注释掉的，如果你想让VMD默认不显示坐标轴，就把#去掉。

VMD默认用称作Normal的Rendermode，但此时有些材质的显示效果很差，甚至Transparent材质根本没法正确显示出透明效果。因此通过display rendermode GLSL将默认的Rendermode设为效果好得多的GLSL。

有很多人肯定早已发现画面边缘的物体畸变得特别厉害，很难看。通过display distance -8.0语句可以充分避免。但导致一个问题就是原本在窗口左下方的坐标轴看不到了，需要坐标轴的时候可以选Display - Axes - Origin让坐标轴显示在窗口中央。

下面这两行是自定义命令。在VMD的文本窗口里输入bw（意为background white）就可以令背景立刻变为白的，输入bb就可以令背景立刻变为黑的，非常方便。
proc bw {} {color Display Background white}
proc bb {} {color Display Background black}

下面的内容是设置用户自定义快捷键。载入轨迹后，在图形窗口处于被激活的状态时（激活窗口就是鼠标点击这个窗口的意思），按左、右键就可以分别后退1帧、向前1帧，按上、下键就可以分别增加10帧、后退10帧。这使得观看轨迹方便很多。
user add key Right {animate next}
user add key Left {animate prev}
user add key Up {animate goto [expr $vmd_frame([molinfo top])+10]}
user add key Down {animate goto [expr $vmd_frame([molinfo top])-10]}

关于VMD判断原子间连接关系的问题我在《谈谈VMD可视化程序的连接关系的判断和设置问题》（//www.umsyar.com/534）里有非常详细的说明。为了能很方便地让VMD对当前帧根据当前结构重新判断连接关系，增加了下面这句。图形窗口处于激活状态时按b键（意为bond）就能重新判断连接关系。
user add key b {mol bondsrecalc all; topo retypebonds}

如果想播放的时候实时自动更新连接关系，而不需要每次都按b键，靠以下语句可以实现。也就是定义了一个bt命令，如果在命令行窗口输入了bt，那么每当当前top体系的帧号发生了变化，就会调用updatebond命令自动来更新连接关系。这样做的代价就是对较大体系，每播放一帧都要根据距离重新判断连接关系，播放时会比较卡。
proc bt {} {
global vmd_frame
trace variable vmd_frame([molinfo top]) w updatebond
}
proc updatebond {args} {
mol bondsrecalc all
topo retypebonds
}

下面的内容定义bn命令。如果不想自动更新成键方式了，可以输入bn命令取消掉对帧号的跟踪即可。
proc bn {} {
global vmd_frame
trace vdelete vmd_frame([molinfo top]) w updatebond
}

雾化效果绝非毫无意义。只要恰当使用，就可以让近距离的物体完全不受影响，而让偏远的原子恰当地雾化，避免扰乱视觉、妨碍清楚地观看近距离的物体。下面的语句是自定义命令，只要在文本窗口输入fog，就可以打开雾化并且使用在我来看比较合适的雾化设置。如果觉得实际效果不好，需要进一步调节，可以用Display - Display Settings，修改里面的Cue Start和Cue End。
proc fog {} {
display depthcue on
display cuemode Linear
display cuestart 1.75
display cueend 2.5
}

On the determination and setting of the connection relationship in VMD visualization program

文/Sobereva@北京科音

First release: 2020-Feb-25   Last update: 2023-Dec-29

化学体系可视化程序VMD的用户越来越多，经常有人在网上问我怎么在VMD里修改键连关系，这里就专门说一下这个问题。本文内容对应于撰文时最新的VMD正式版1.9.3，对其它版本可能适用也可能不适用。

当VMD载入pdb、xyz、gro等格式的结构文件时，程序就会自动判断连接关系。判断规则和Multiwfn、GaussView等程序类似，都是根据两个原子间距离以及这两个原子半径之和来判断的，参看此文的相关说明：《谈谈原子间是否成键的判断问题》（//www.umsyar.com/414）。VMD里默认的line方式，以及常用的CPK、Licorice等显示风格都是根据连接关系显示的。

具体来说，若原子间距离小于两个原子半径之和的0.6倍就会被VMD视为成键。原子半径是VMD内置的，大部分是Bondi范德华半径。元素的默认半径没法在VMD图形界面里或配置文件里改，但载入结构后可以通过命令行改特定原子的半径。

如果两个原子间没有被判断为成键，但你希望让它们之间显示键，就选择Mouse - Add/Remove Bonds，然后点击这两个原子的正中心，它们就连上了。如果你点击两个已经成键的原子，它们之间的键就会被去掉。

在你点击原子时，VMD会自动给你点击的原子加上标签并显示在图上，如果不想显示，就按键盘上的1键，再点击有标签的原子，标签就被隐藏了。如果你想把一批标签都隐去或者彻底删掉，就进入Graphics - Labels，按住鼠标左键选中一批，然后点Hide按钮就可以隐藏，点Delete按钮就可以删掉。如果想让批量删除标签的操作尽可能省事，对于Windows版VMD，你可以在VMD目录下的vmd.rc里加入如下内容
proc da {} {
label delete Atoms all
}
启动VMD后，随时都可以在文本窗口里输入自定义命令da来删除所有原子标签。

如果有很多同类的键，它们的连接关系判断得都和你要求的不符，一个一个靠鼠标点击来修改连接关系很麻烦。此时可以在Graphics - Representation界面里把Drawing Method改为DynamicBonds，这代表基于当前结构实时判断成键并显示，判断阈值可以用旁边的Distance Cutoff控件来修改，恰当改成令显示的连接方式满足你要求的情况即可。但这种显示方式有个缺点就是不显示原子球，而只显示键，而且键是空心的，很不好看。因此为了让原子球也显示出来，你可以点Add Rep增加一个Representation，把Drawing Method设成VDW，并且把Sphere Scale改小到合适。

有的时候在用DynamicBonds方式显示时，通过调节Distance Cutoff虽然让A-B之间键连方式满足要求了，但同时C-D之间的连接方式变得却与期望的不符了，这种时候就需要灵活变通，多设几个Representation，每个都恰当用不同的选择语句选定绘制的区域（参考《VMD里原子选择语句的语法和例子》//www.umsyar.com/504），并用恰当的Drawing Method和参数设定，从而让显示出的效果叠加起来正好满足你的要求。

有的文件格式包含了连接关系信息。比如AMBER、NAMD的拓扑文件都可以被VMD载入，载入的时候就相当于提供了连接关系了，之后再往这个ID里载入结构文件或轨迹时也因此都是按照拓扑文件里的连接关系来显示。VMD虽然不能载入GROMACS的拓扑文件，但是如果安装了http://bbs.keinsci.com/thread-37839-1-1.html中提供的tpr插件的话，载入GROMACS的tpr时也会从中读取正确的连接关系。还有的结构文件格式自己也带了连接关系段落，典型的就是mol2格式（mol格式也有，但VMD不能载入），相减《谈谈记录化学体系结构的mol2文件》（//www.umsyar.com/655）。你用文本编辑器打开它的话，会看到@<TRIPOS>BOND字段，这部分记录的是键的形式，诸如5 3 27 2就代表这是第5号键，对应3与27原子间的双键。当VMD载入mol2格式文件的时候，就不会自动根据结构来判断连接关系，而是直接用文件里的连接关系。pdb格式里有个CONECT字段，也是记录连接关系用的，但是VMD并不会采用这里的连接关系。如果你想让VMD直接利用这里的连接关系，方法看《使VMD根据pdb文件中的CONECT字段设定原子连接关系》（//www.umsyar.com/121）。

在GaussView里可以编辑键连关系，而且从GaussView 6开始，可以自定义不同元素之间成键的判断阈值，这在前述的《谈谈原子间是否成键的判断问题》文中已经说过了。改过之后只点击Edit - Rebond就可以按照用户自定义的阈值重新判断成键，因此在批量修改特定类型元素间的成键方式的时候比较方便。GaussView可以保存mol2格式的文件，里面的连接关系信息正对应于你在GaussView图形界面里当前看到的。因此，可以先让GaussView正确显示出你期望的连接关系，保存成mol2文件，再用VMD读取，此时VMD里显示的连接关系就和GaussView里精确一致了。

注意VMD里没法以GaussView的风格来显示双键、三键等成键形式。比如载入VMD的mol2文件离记录了某两个原子间形成的是二重键，但在VMD的Line、CPK、Licorice显示风格显示时看起来和单键一样。VMD的Bonds显示风格倒是能区分形式键级，但双键和三键分别是以两个、三个套筒形式显示的，非常难看，没实际意义。

VMD自身也可以保存mol2文件，因此你在VMD里通过Add/Remove bonds好不容易把连接关系改合适后，不妨通过File - Save coordinate功能保存成mol2文件，以后再次显示时就免得再修改一遍了。

常有人问我怎么载入动力学轨迹后，虽然第一帧的成键方式正确，但播放轨迹时成键方式有时连得乱七八糟，甚至跨越盒子。这是因为如前所述，连接关系要么是从拓扑文件里读取的（比如看Amber轨迹前要载入拓扑文件），要么是载入结构文件时自动判断的（比如看GROMACS轨迹前通常要载入与模拟体系原子信息对应的gro或pdb文件），或者是根据轨迹文件的第一帧结构判断的（比如载入molclus、xtb程序产生的多帧xyz文件），之后再载入轨迹或播放轨迹时都不会改变连接关系。由于动力学过程中可能发生分子解离、异构化，或者由于周期边界条件的原因分子时而完整时而被盒子截断，因此内存里记录的连接关系列表可能对于轨迹的某些帧并不适用，极度影响正常观看。像这种情况，可以用DynamicBonds显示方式，这样就会按照当前帧的坐标实时判断成键并显示（但并不会改写内存里的连接关系列表）；但如果你想用其它显示方式，就得让VMD根据当前这一帧的结构重新判断连接关系并改写内存里的连接关系列表，做法是在文本窗口输入mol bondsrecalc all; topo retypebonds。此时连接关系就是满足当前这一帧的情况了，而对其它帧可能就不适合了。为了方便，你可以定义重新判断连接关系的快捷键，甚至可以让VMD自动对每一帧实时判断连接关系，做法见《VMD初始化文件(vmd.rc)我的推荐设置》（//www.umsyar.com/545）。

如果你想把当前内存里的连接关系列表以文本方式显示出来，就在文本窗口输入[atomselect top all] getbonds，输出的信息比如
{1 2 3 4 7} 0 0 0 {5 6 7 0} 4 4 {4 0}
这是个列表，里面共有8个成员（有的成员自身就是个列表），对应体系一共8个原子的成键情况，比如第0号原子与1、2、3、4、7成键，第2号原子与第0号原子成键，第3号原子也与第0号原子成键...注意这体现的是连接关系，而非键级。

连接关系也可以从自行提供的列表里读取。比如运行
[atomselect top all] setbonds {{1 2 3 4 7} 0 0 0 {5 6 7 0} 4 4 {4 0}}
就代表把当前体系的连接关系改成{1 2 3 4 7} 0 0 0 {5 6 7 0} 4 4 {4 0}所定义的。

使用Multiwfn和VMD绘制平均局部离子化能(ALIE)着色的分子表面图（含视频演示）

Using Multiwfn and VMD to plot map of molecular surface colored by averaged local ionization energy (ALIE)

文/Sobereva@北京科音

First release: 2019-Sep-17  Last update: 2020-Oct-10

平均局部离子化能(ALIE)是一个对于考察分子亲电反应位点、活性非常有用的量，在《静电势与平均局部离子化能综述合集》（http://bbs.keinsci.com/thread-219-1-1.html）我给出了很多综述、介绍性文章以及笔者写过的相关博文，不熟悉ALIE者强烈建议阅读。Multiwfn是考察ALIE最方便、灵活、强大的程序。通常分析ALIE时都是考察它在分子表面的分布。之前笔者有个文章《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》（//www.umsyar.com/443）介绍了如何将Multiwfn与VMD程序联用非常方便、快速地绘制高质量的静电势着色的分子范德华表面图。由于后来又有不少人问我怎么类似地也绘制ALIE着色的分子表面图，因此这里笔者就专门说一下。本文对应于Windows环境，若想在Linux下实现同样绘制效果，请读者灵活变通。

本文的过程配有全程演示视频，强烈建议一看：https://www.bilibili.com/video/av68116168/

首先大家去//www.umsyar.com/multiwfn下载Multiwfn程序，必须是2019-Sep-17及之后更新的版本，之前的版本没有下面提到的脚本。如果还没装VMD，去http://www.ks.uiuc.edu/Research/vmd/下载。本文用的是VMD 1.9.3。

之后进行以下操作：
·将examples\scripts里的ALIE.vmd复制到VMD目录下
·将examples\scripts里的ALIE_isoext.bat和ALIE_isoext.txt复制到Multiwfn.exe所在目录下
·编辑ALIE_isoext.bat批处理文件，将里面的文件路径改成实际要考察的，把VMD路径也改成实际的。输入文件可以用.wfn、.fch、.molden等
·双击ALIE_isoext.bat，Multiwfn将被自动调用进行计算，算完后会在VMD目录下产生avglocion.cub、density.cub、surfanalysis.pdb
·启动VMD，在命令行窗口输入source ALIE.vmd，就可以得到想要的图像了

下图是绘制出来的菲的ALIE着色的分子表面图像。

图中的表面对应的是电子密度为0.0005 a.u.的等值面。之所以不是常用的0.001 a.u.等值面（通常被用来作为分子范德华表面的定义），是因为此时ALIE着色效果不好，不容易通过颜色区分不同区域ALIE分布差异。上图中青色圆球是分子表面上ALIE的极小点，体现电子被束缚得极弱的位置，也因此体现出容易发生亲电反应的位点。ALIE是按照蓝-白-红的色彩过渡方式着色的，因此上图中蓝色区域都是ALIE比较小的区域，电子比较容易参与反应。

为了效果更好，可以用Tachyon渲染器渲染，在视频里已经演示了。

若你研究的体系不容易通过色彩区分ALIE在不同区域的相对大小，或者整个分子表面颜色全都是红色或蓝色，说明当前用的ALIE的色彩刻度范围不合适，可以在VMD命令行窗口输入诸如mol scaleminmax 0 1 0.31 0.37，这里0.31和0.37分别是色彩刻度下限和上限（亦可以在VMD的图形界面里修改，做法是进入Graphics - Representations，在窗口上方选择当前是Isosurface的那项，点击Trajectory标签页，在文本框里输入色彩刻度的下限和上限之后按回车）。色彩刻度范围怎么设合适和具体体系有关，可以反复试试。如果心里没谱的话，可以先把下限和上限分别设成当前分子表面上的ALIE实际的最小值和最大值，以让色彩变化能完整表现分子表面上的ALIE分布。这两个值在Multiwfn计算过程的输出信息里就可以看到（将ALIE_isoext.bat第一行末尾加上[空格]> out.txt，脚本调用Multiwfn期间输出的信息就会导出到当前目录下的out.txt，其中搜Minimal value:和Maximal value:），然后手动将这两个值除以27.2114换算为以Hartree为单位的值，作为色彩刻度下限和上限即可。

对于平面体系，如果想准确考察ALIE极值点与原子或者键的相对位置，最好改成正交视角，即选择Display - Orthographic，否则会由于近大远小的缘故可能导致误判位置。以正交视角显示时的菲分子的图像如下，位于蓝色区域的原子或键容易参与亲电反应。

值得一提的是，笔者在《一篇最全面、系统的研究新颖独特的18碳环的理论文章》（//www.umsyar.com/524）介绍的论文中使用本文的方法考察了电子结构非常特殊的18碳环体系，得到的图像很漂亮（如下所示），体现出的信息又非常有价值，很建议读者看看此文中关于反应性的相关讨论。

如果想像上图一样把色彩刻度轴画出来，按照《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图（含视频演示）》（//www.umsyar.com/443）里面介绍的做法就可以实现。如果你想查询分子表面ALIE极值点具体的值，也效仿此文里面说的做法。