使用Multiwfn图形化研究弱相互作用
杨伟涛课题组在名为Revealing Noncovalent Interactions (JACS,132,6498-6506)一文中提出了一种新的可视化研究弱相互作用的方法，概念简单清晰，具有广泛意义，很有实用价值，本文目的之一在于介绍、推广这一方法，但与原文的角度有一些不同。使用这一方法需要计算空间上各点的RDG函数和sign(λ_2)ρ函数的值，虽然可以使用杨伟涛课题组开发的免费的NCIplot程序(http://www.chem.duke.edu/~yang/Software/softwareNCI.html)进行计算，但是使用起来有诸多不便。在Multiwfn 1.5程序中不仅可以实现NCIplot的所有功能，还做了很大改进，使这一分析方法更容易投入实际应用，本文目的之二就是结合实例介绍Multiwfn做这种分析时的操作方法。
这种分析方法的研究对象从距离上讲是中、长程相互作用，从类型上讲主要包括氢键、静电、范德华作用，也能展现位阻作用。虽然原文作者称这种研究方法的对象是非共价相互作用，笔者认为称为弱相互作用更为妥当，文中将使用这种称呼。这种分析方法也存在一些局限、弊端、随意性，这将在日后另文讨论，本文只以正面角度来讨论。本文所用的Multiwfn是一个免费、开源、方便、灵活的波函数分析程序，可以从http://multiwfn.codeplex.com下载最新版本，文中例子中的操作步骤是针对1.5版而言的，若以后版本中选项有所改动，根据程序中选项的含义就能明白怎么操作。
1. 用RDG函数等值面显示弱相互作用区域
此方法的主要目的在于凸显出体系中涉及弱相互作用的区域，以便直观地了解到分子中哪些区域与弱相互作用有关。也就是定义一个实空间函数，使其数值能够区分开体系中具有不同特征的区域。此方法使用的是约化密度梯度函数(Reduced density gradient, 下文称RDG)，其表达式为1/(2*(3*π^2)^(1/3))*|▽ρ(r)|/ρ(r)^(4/3)，其中▽是梯度算符，|▽ρ(r)|即电子密度梯度的模。实际上前面的常数项（其值为0.16162）可以忽略掉，只考察|▽ρ(r)|/ρ(r)^(4/3)部分。一个分子体系主要由以下四部分构成： [表1]
表格中的大、小的标准比较模糊，只是定性的。如果我们想要找弱相互作用区域，利用RDG函数的数值大小差异就可以将“原子核附近”和“分子边缘”区域去掉，但“弱相互作用区域”和“化学键附近”的RDG函数值、|▽ρ(r)|值都比较小，区分不开，但ρ(r)存在一定差异。所以，结合使用RDG函数和ρ(r)函数，就可以确定分子中哪些区域涉及弱相互作用。
如果设定立方网格，使网格中的点能够覆盖整个体系，做这些点上ρ(r) vs. RDG的散点图，就可以把上述概念图形化且定量地表述出来。下面来做甲烷二聚体的这种图。
首先要利用Gaussian生成体系的波函数文件(.wfn)。写一个甲烷二聚体的输入文件，route section写上#p opt mp2/aug-cc-pVDZ density out=wfn，坐标后面空一行写上.wfn文件的输出路径，附件里的methanedimer.gjf是已写好
的。用Gaussian03执行后得到methanedimer.wfn。
我建议在Gaussian输入文件中的Title Section部分写上与route section一致的内容（当然不要把#也写进去），这些内容将会自动传递到.wfn文件中的第一行，以免以后忘记.wfn是在什么条件下生成的。也可以不生成.wfn文件，因为Multiwfn可以读入.fch文件，所得结果将是一样的。要注意.wfn文件记录的波函数的基组角动量最高只支持到f函数（因为更高角动量函数在.wfn文件中没有标准的定义），如果用涉及到g角动量的基组，比如cc-pVQZ，就只能通过.fch文件把波函数信息传递给Multiwfn。
先把multiwfn.exe目录下的Settings.ini里的RDG_maxrho设为0.0（注意等号后面要留空格），其原因后文会解释。然后启动Multiwfn，依次输入： c:methanedimer.wfn //输入文件的路径
100 //功能100，其中包含Multiwfn中比较杂的功能 1 //绘制“函数1 vs. 函数2”散点图并生成相应格点文件
1,13 //输入函数1和函数2的序号，分别作为散点图的横轴和纵轴。在Multiwfn支持的函数中ρ(r)是第1号，RDG函数是第13号。
2 //用中等质量的网格，总共约512000个点，x,y,z方向的具体点数通过使x,y,z方向格点间距相等来自动确定。网格的区域自动往分子外延展6 bohr。
现在开始计算格点数据。格点数越多、体系所含Gauss函数越多，计算速度越慢，计算时间与二者都成正比。计算完毕后，输入
4 //设定散点图X轴 0,0.35 //X轴上下限的值 5 //设定散点图Y轴 0,2 //Y轴上下限的值 -1 //绘制散点图
很快ρ(r) vs. RDG的散点图就弹出来了，如图2左图所示。在图上点右键可以关闭，然后选功能1可以将图保存到当前目录下（即Multiwfn.exe所在目录下，后同）。如果对图的效果不满意，可以选功能2将数据导出到当前目录下output.txt，然后用origin、sigmaplot等程序做散点图。其中前三列是数据点的坐标，后两列是两个函数的值。
[图2]
按上述步骤绘制甲烷孤立状态的散点图得到图2右图（设定网格时选3，用高质量网格）。从图2可以看出，体系中存在与不存在弱相互作用时散点图最主要区别在于图中最左侧是否有一个竖条，在原文中这被称为spike。这个竖条上的点正是弱相互作用区域“RDG数值为0~中，ρ(r)^(4/3)数值为小”所对应的点；在右侧也有一个区域RDG 值接近0，这是C-H键区域“RDG数值为0~较小，ρ(r)^(4/3)数值为中”所对应的格点；图中坐标轴范围的更右侧就是原子核附近区域的点；图中左上角的尖峰再往上继续延伸就是分子边缘的区域，虽然离分子越远的地方|▽ρ(r)|和ρ(r)^(4/3)都越小，但后者比前者减小得更快，所以离分子越远RDG值越大，并直至无限大，可自行调整坐标轴观看。不同体系的散点图的成键、原子核附近、分子边缘区域都是类似的，一个体系中是否含有弱相互作用，就是看在ρ(r)较小区域是否有spike出现，这是此分析方法的要点。当然，网格不能太稀疏，如果在弱相互作用区域恰好没有点，spike也不会出现。
接下来，要用等值面确定这些对应于弱相互作用区域的点在实空间上的位置。在计算完格点后的那个菜单中，输入7，然后输入想看的等值面就可以观看第2个函数（即RDG函数）的等值面。其0.5的等值面如图3左图所示
[图3]
图上在两个甲烷中间出现了封闭的等值面，描绘的正是二者间的范德华作用，但是在分子上也出现了三角形封闭的等值面。这是因为成键区域和弱相互作用区域的RDG函数数值范围有很大程度的重叠，如果在图2左图上做一个y=0.5的直线，会发现这条直线不仅贯穿spike，还贯穿成键区域，所以相应的RDG=0.5的封闭等值面不止一个，而在成键区域附近也会出现。这也是前面所说，必须再通过ρ(r)函数区分开成键和弱相互作用区域。
屏蔽掉成键区域，也就是将ρ(r)值稍微大一些的区域，比如ρ(r)=0.05以上的RDG函数的数值设得很大，比如设成100，这样在散点图上y=0.5的直线就不会经过那个区域了，等值面也就只剩下弱相互作用区域。具体做法是在之前的菜单中选择-3，然后输入0,0.05，再输入100，这就表明将ρ(r)的范围在[0,0.05]以外区域的点的RDG函数数值设为100。重新绘制散点图，得到了图4的结果，等值面也变为了所期望的图3右图的情况。
[图4]
一般来讲，观看RDG函数一般观看的是0.5的等值面，这没有什么严格的物理意义，只是等值面大小比较适中。由于弱相互作用区域的ρ(r)一般不会越过0.05，散点图上y=0.5的直线在ρ(r)<=0.05的区域内也只与spike相交，所以每次作弱相互作用等值面图时没必要再考察散点图，只需直接将ρ(r)>0.05的点的RDG函数值设为较大数值就行了。由于这个步骤经常要做，为了方便，Multiwfn在settings.ini里面有一个RDG_maxrho参数，凡是涉及到计算RDG函数的功能，只要某个点的ρ(r)大于这个参数，这个点的RDG值就自动被设为100，这个参数默认被设定为0.05。所以用户就不需要再考虑屏蔽掉成键区域了，这已由Multiwfn自动完成。当然，在后文中作完整的散点图时、或者就是想通过RDG函数研究成键区域时，应当关闭这个功能，将RDG_maxrho设为0.0就代表关闭此功能。
2. 合理地设定网格