内容发布更新时间 : 2026/6/25 22:36:50星期一 下面是文章的全部内容请认真阅读。
[图13]
虽然能显示高斯类型格点文件的等值面的程序很多，但支持将一个函数数值用不同颜色填到另外一个函数的等值面上的
可视化程序比较有限，常用的GaussView虽然支持，但操作不便而且不够强大。VMD是观看、分析分子动力学结果最重要的软件之一，它在映射颜色和显示等值面方面也很好用，等值面又光滑又有光泽，填色的色彩变化细腻，调整等值面也比较容易，而且运行流畅。VMD可以免费在http://www.ks.uiuc.edu/Developme ... cgi?PackageName=VMD下载。
首先安装VMD，然后将func1.cub和func2.cub复制到VMD安装后的目录下，即vmd.exe所在路径。然后在此目录下编写一个文本文件，后缀为.vmd，比如ltwd.vmd。在里面填上如下内容： mol new func1.cub mol addfile func2.cub mol delrep 0 top
mol representation CPK 1.0 0.3 18.0 16.0 mol addrep top
mol representation Isosurface 0.50000 1 0 0 1 1 mol color Volume 0 mol addrep top
mol scaleminmax top 1 -0.04 0.02 color scale method BGR
保存文件后，启动VMD，选File-Load State，选择ltwd.vmd，图13下方的图就显示来了。若想把背景改成黑色，选Graphics-Colors-Display-Background-8 White。图中的弯曲的片状等值面边缘略有锯齿，可以通过增加此处格点密度来改善。从颜色上可看出，弯曲的片状等值面描述的是二聚体之间范德华作用，但部分区域也有微弱的位阻效应。圆片等值面只有中间呈蓝色，说明H与O之间形成了氢键，但并不像甲酸二聚体的氢键那么强。
VMD里面每个操作对应一条语句，载入.vmd本质上就是让.vmd文件中的语句全部执行，这就免得每次都手动执行载入文件、设定参数的一系列繁琐步骤。比如mol new func1.cub的含义就是读入当前路径下func1.cub文件（默认的当前路径一般就是vmd.exe所在路径），mol scaleminmax top 1 -0.04 0.02代表将1号representation（对应于显示填色等值面的那个图层）色彩刻度的下上限分别设为-0.04和0.02，color scale method BGR代表将色彩刻度由小到大设为蓝->绿->红。将背景设为白色的命令是color Display Background white，若将此命令添加到.vmd里面就能在载入.vmd文件时顺便执行，使背景自动设为白色。这些命令在VMD手册上都有解释。这些命令也可以在VMD的控制台下直接运行，控制台通过Extensions-Tk Console进入，比如想把色彩刻度改为从-0.05到0.05，就在控制台执行mol scaleminmax top 1 -0.05 0.05。有很多命令在VMD的GUI上有相应的选项，其中有些通过GUI操作会容易得多，比如调整等值面数值，可以在Graphics-Representations里面的列表中选定第二个显示模式（即Style为isosurface的那个），然后将下方Range的下上限分别设为比如0和1，点回车，之后拉动滑条就能在0~1范围内改变等值面。
另外，Chemcraft也可以实现相同功能，对初学者来说使用更简单，不过效果不如VMD，而且收费。使用方法： 先打开func2.cub，再点左下角Add cube，选择func1.cub。将Contour value设为0.5，敲回车，然后点Show isosurface。然后把Map other:选为2，再把Values range分别填-0.04和0.02，敲回车。
这种分析方法也可以用于分析晶体间内部弱相互作用。虽然Gaussian的PBC功能并不给力，但是简单的PBC计算还是可以胜任的。由于Gaussian的PBC计算是以高斯型函数作为基函数，所以波函数信息可以直接被Multiwfn读入并进行分析。这里将以金刚石晶体作为例子。
还是先获得波函数文件。Gaussian的输入文件就是压缩包里的pbc_diamond.gjf。运行之后，用formchk将.chk转化为.fch文件。当然，用.wfn作为Multiwfn的波函数输入也可以。由于计算的是素晶胞，波函数信息也只含有两个碳原子的，获得的等值面显然体现不出周期性，然而计算复晶胞又太费时。在Multiwfn里提供了一个晶胞波函数平移复制
功能，可以将这素晶胞波函数信息扩展为足够大的复晶胞波函数。复制次数越多，体系中高斯函数越多，计算格点时越慢，为避免计算格点时间太长，这里只向各方向平移复制两次。
启动Multiwfn后输入： c:pbc_diamond.fch 6 //修改波函数功能 14 //平移复制体系
4.7523,0,0 //第一个平移向量。平移向量在输出文件中的PBC vector段落，或者查看.fch中的Translation vectors字段。不要用输入文件中的平移向量，因为Gaussian可能自动改动坐标，平移向量也就变了。 1 //单位为bohr
2 //复制两次。接下来再对另外两个方向做相同的平移复制。 2.3762,4.1156,0 1 2
2.3762,1.3719,3.8802 1 2
0 //平移复制已完毕，保存当前波函数到当前目录下new.wfn，也就是压缩包中的pbc_diamond_dup.wfn。
由于得到的体系是斜着的，把所有原子都纳入立方网格中必将造成很多格点落在体系外而浪费掉。实际上只要取体系中间一个原子（28号）作为网格中心，然后向四周延展一些距离，所得等值面就足够体现周期性了。现在生成格点文件，依次输入
pbc_diamond_dup.wfn 100 1 15,13
7 //用两个原子的中点作为网格中心
28,28 //当输入的两个原子序号相同，则以此原子为中心向四周延展。Multiwfn默认延展距离是6 bohr，对此例子比较合适，不用修改。 80,80,80 //各方向格点数
算完后，用功能3保存格点文件，按照上例方法，用VMD显示结果如下。可见，由于碳原子间距离较近，原子空穴间充满红色等值面，体现很强的位阻效应。
[图14]
5 波函数质量产生的影响
在前面的例子中多数用的是B3LYP/6-31G*波函数，必定有人会认为用B3LYP/6-31G*计算以范德华作用为主的弱相互作用体系很不合理，但实际上相互作用能与理论方法关系更大，其值很差不代表相应的ρ(r)就会差，原文作者认为使用这种分析方法时用B3LYP/6-31G*就够了，这样计算量也小，而且改为MP2/6-311++G**后等值面并没有什么改变。甚至ρ(r)粗糙到只用自由电子密度叠加都能得到定性一致的结果，所以这种分析方法对电子密度的质量很不敏感，使之用于大体系成为了可能。
6 Promolecule近似及对结果产生的影响
使原子坐标保持在形成分子时的状态，将自由原子密度叠加得到的密度称为Promolecule密度，可以视为在形成分子前，电子密度尚未驰豫的电子密度。构建这种密度不需要量子化学计算，只要有分子结构文件和自由原子密度就能十分容易地构建。由于在弱相互作用区域Promolecule密度和实际密度差异并不像在成键区域那么显著，Promolecule密度在此分析方法中可以近似代替实际电子密度，等值面的基本形状不会有太大差异，但是定量上，也即等值面细部特征会有一定差异，因为电子密度驰豫后会在吸引作用区域聚集，尤其是会在位阻效应区域疏散以减小斥力，位阻效应越强，改变量越大。
原子在自由状态的真实电子密度是球对称的，但是对于比如基态氧原子，p轨道一个是双占据两个是单占据，量化算出
来的电子密度不是球对称的，这将导致分析结果出现不应有的取向性，所以需要将之球对称化。球对称化方法不是唯一的，比如可以简单地对三个p轨道占据数取平均，也有人利用GVB方法解决，原文的方法是用s型STO函数拟合B3LYP/6-31G*原子密度，第一、二、三周期的原子分别用1、2、3个STO拟合，以对应壳层数。这不仅解决了球对称问题，还有另一个好处，就是描述原子密度的函数大大减少了，6-31G*描述氧原子用28个GTO，而拟合后只用2个STO，这使得计算RDG函数、sign(λ_2)ρ的速度也大大加快。实际上STO用的少主要在于双电子积分时的困难，由于STO能正确地表现原子轨道波函数随r增大的收敛行为和原点处Cusp的特征，如果研究内容不涉及到双电子积分，比如这种只依赖ρ(r)的分析方法，用STO又便宜又好。
在Multiwfn中使用Promolecule密度计算RDG与sign(λ_2)ρ函数，只需在选择要计算的函数的界面中选择后面带着\的相应函数即可。注意此时体系中只能存在前三周期的原子。由于不涉及到波函数信息，只要将分子结构传递给Multiwfn就可以，所以既可以通过.wfn、.fch文件将分子结构传入Multiwfn，也可以通过Multiwfn支持的.pdb文件导入分子结构。在Promolecule近似下苯酚二聚体的散点图和等值面图如下所示