使用Multiwfn作电子密度差图

文/Sobereva   First release: 2011-Dec-24

1 前言

在Multiwfn里可以十分方便地做密度差图，在程序手册里也已经给出几个例子很清楚地介绍了绘制密度差图的方法，但是仍然总有人问绘制方法，所以在此写个专门帖子通过实例来说说绘制过程，会把每一步都做出详细解释。本文同时也会谈到密度差图的意义、用处。另外，对于作密度差图常常被忽视的原子密度球对称化问题也会着重说一下。本文所用Multiwfn版本为2.2.1 rev1，Gaussian版本为09W A02。

做密度差平面图还有一种老方法，是通过cubegen/Multiwfn生成格点文件，然后利用GsGrid求差值获取截面数据，经过Excel处理后导入Sigmaplot绘制的。Zhou2009近期在小木虫量化版上发表的《氢键二聚水量子化学图形制作过程》（http://muchong.com/bbs/viewthread.php?tid=3977525&fpage=1）对这种老方法做了完整、系统的总结。我相信，读者如果亲自用老办法和本文介绍的方法都做一遍密度差图，就会深感直接使用Multiwfn来绘制密度差图在操作上容易得多，也快捷得多。

本文所说的密度都是指电子密度。利用和本文完全相同的方法，还可以做其它函数的差值图，比如静电势差值图、动能密度差值图、电子定域化函数(ELF)差值图等等，只要在选择实空间函数时选择对应的实空间函数就可以了。


2 密度差图的种类

“密度差”就是指一些体系在各自特定的状态下彼此间的密度的差值，它是三维实空间函数。有三种情况常会涉及：
1 某个体系的某个状态的密度减去这个体系在另外一个状态的密度。例如，某分子的激发态的密度减去它在基态时的密度、某分子在有外加电场时的密度减去它在自由状态下的密度。注意在计算这两个状态的密度时分子的笛卡尔坐标必须一致，不能在各自状态下分别优化，而只能共用其中一个状态的几何构型。同时也不能允许量化程序为了能利用分子对称性而自动旋转、平移坐标。否则都导致两个状态下原子坐标没法对应，密度差也就没什么意义。

2 某个体系的密度减去组成它的各个片段的密度。例如，水分子二聚体的密度减去这两个水分子都在孤立时的密度。同样要注意，求各个片段的密度时，它们的笛卡尔坐标必须与在整个体系中的坐标一致。

3 某个体系的密度减去组成它的各个原子在自由状态下的密度，这也被称为变形密度(Deformation density)。这种密度差对于研究分子形成过程中电子密度的变化、探究化学键的本质很有用。

密度差必须做成图才便于被直观地考察，图分为三类：

1 曲线图：展现某条直线上密度差的数值，这种图对于密度差的研究用得较少。
2 平面图：能够清楚地展现一个截面上各处的密度差数值，十分有用，其中最适合展现密度差的是等值线图。平面图的缺点是难以从一个截面上了解整个空间中密度差的数值分布。
3 等值面图：在等值面上每个点的数值都等于人为指定的数值(这个值被称为isovalue)，通过它可以很直观地了解密度差在三维空间中的分布。但是等值面里面和外面的密度差数值没法一下子看清楚，必须反复调节isovalue考察等值面的变化才能获知。等值面图和平面图的用处是互补的。


3 Multiwfn中的自定义运算功能

Multiwfn可以直接计算密度并做曲线、平面、等值面图，Multiwfn也提供了自定义运算(custom operation)功能，用于对多个波函数的密度、静电势等实空间函数进行相互运算，求密度差实际上只是Multiwfn的自定义运算功能的一个具体应用而已。

比如，启动Multiwfn后载入的波函数是A.wfn，在自定义运算中设定有3个波函数将要依次对它进行运算，且算符和文件名分别写
+,madoka.wfn
x,TMA.wfn
-,homura.wfn
然后，假设要计算的属性选择的是电子密度，那么Multiwfn会依次计算这四个波函数文件的电子密度格点数据，并且对这些格点数据做如下运算：(A.wfn的密度+madoka.wfn的密度) * TMA.wfn的密度 - homora.wfn的密度
在计算后三个波函数的格点数据时Multiwfn都会自动使用与计算A.wfn格点数据时相同的格点设定，无需用户考虑格点位置的对应性问题。
Multiwfn可以对无限多个波函数进行自定义运算，支持的算符是+ - x /，对应加减乘除（乘号写为字母x）。

下面将通过一些实例，说明怎么用自定义运算功能绘制各种类型密度差图。


4 某体系某状态减去此体系另一个状态的密度差图

4.1 施加电场后COCl2电子密度的变化的等值面图

在外场作用下，电子密度分布会被极化。这里考察一下COCl2分子在C-O轴线方向上加上0.03a.u.均匀电场（1a.u.=51.421 V/Angstrom）后电子密度的变化。注意电场作用也会诱导体系几何结构发生改变，但绘制密度差要求几何结构保持不变，所以这里不考虑电场对分子结构的影响，也因此不在电场下重新优化结构。

首先生成COCl2在自由状态下的波函数，将以下内容保存为Gaussian输入文件并执行即可，wfn文件输出路径和文件名可以自定。这个结构是事先在b3lyp/6-31+G*下优化好的。
#p b3lyp/6-31+G* out=wfn nosymm

b3lyp/6-31+G* opted

0 1
C                  0.00000000    0.00000000    0.49676300
Cl                 0.00000000    1.46464400   -0.48297200
Cl                 0.00000000   -1.46464400   -0.48297200
O                  0.00000000    0.00000000    1.68005900

c:\COCl2.wfn
算完后将这个输入文件的route section部分末尾加上field=z+300关键词，并将输出的wfn文件名改为c:\COCl2_z+300.wfn，然后重新算一次以得到在Z方向上（C-O轴线方向）上加了0.03a.u.均匀电场后的波函数文件。

如果你对Gaussian缺乏经验的话，在计算用于做密度差的波函数文件时总要写上nosymm，并且使用笛卡尔坐标，这样可确保分子三维空间中的位置、朝向是确定的。nosymm关键词避免了Gaussian自动令分子旋转、平移以调整到标准朝向下。

我们先来做等值面图。启动Multiwfn，依次输入
c:\COCl2_z+300.wfn
5 //生成格点文件
0 //设定自定义运算
1 //只有1个文件要对已读入的波函数进行运算
-,c:\COCl2.wfn //运算符（减号）和文件名
1 //选电子密度函数。因此接下来得到的就是COCl2_z+300.wfn的密度减去COCl2.wfn的密度
2 //中等质量格点精度
现在Multiwfn开始依次计算这两个波函数文件的密度的格点数据，之后自动相减。接下来有诸多选项可以后处理生成的格点数据，用选项2可以将格点数据输出为cube文件，之后可以载入进gaussview、VMD、chemcraft、Molekel等软件观看等值面。Multiwfn也自带了预览等值面的功能，这里选择-1来激活之，如下图所示。由于默认的isovalue为0.05，而密度差数据通常很小，所以要将isovalue值降低一些，比如降到0.005，才能比较清楚地看到密度差图。isovalue的大小的设定没绝对标准，只要设定一个值，能让密度差特征充分清楚地展现出来就行了。设得太大时等值面会很小而看不清楚甚至消失不见，而设得太小会使等值面涵盖的空间范围太大，混在一起，看不出到底哪些区域的电子密度变化较大。



密度差数值与isovalue的符号相同和相反的等值面分别以绿色和蓝色显示。图中isovalue被设为0.005，也就是说绿色等值面上密度差值都为+0.005，体现了加了电场后电子密度增加的区域；蓝色的面上密度差都为-0.005，体现了加了电场后电子密度降低的区域。本例中外加电场矢量是由Z正方向指向负方向的，受其影响电子密度应趋向于向Z轴正方向转移，图中也验证了这一点，可以看到主要发生极化的是C-O键的pi电子云。下个例子将绘制穿过C-O键的垂直于分子平面的截面上的密度差图，这个截面对应于Y=0的XZ平面。

4.2 施加电场后COCl2电子密度的变化的等值线图

如果你刚做完上例子的等值面图，可以在格点数据后处理界面中选择0从而直接退回到主菜单。如果直接做这个例子，那么还是先启动Multiwfn，然后输入c:\COCl2_z+300.wfn来载入之。之后依次输入
4 //绘制平面图
0 //设定自定义运算，后面几步和4.1完全相同
1
-,c:\COCl2.wfn
1
2 //等值线图
直接按回车使用推荐的格点设定，即平面图中两个方向都是200个点。
2 //XZ平面
0 //表明这是Y=0的XZ平面
瞬间就弹出了下面的图



实线代表电子密度增加的区域，虚线代表减小的区域。这幅图比等值面图展现的密度差更为清楚具体。在图上点击右键关闭它之后可以看到一个菜单，其中有丰富的选项，用于保存图片、调节等值线设定、调节坐标轴刻度、选择是否将等值线的数值标在图中等等，请大家自行尝试，在手册对应的章节也都有细致的说明。从这个例子可见使用Multiwfn做密度差平面图是极其容易的，而且图像精细漂亮，可直接作为文献的插图。既不需要费时费事地去自行计算格点数据，也完全不需要在excel、sigmaplot里面折腾数据和调节设定。

4.3 电子激发前后的密度差图

还是以COCl2为例子，这次来看看第一单重激发态S1相对于基态的电子密度的变化。还是用4.1节的COCl2的输入文件，但这次将route section改为
#p b3lyp/6-31+G* out=wfn nosymm td(root=1,singlet) density
并且将末尾定义的wfn文件名改为c:\COCl2_S1.wfn。用Gaussian计算之，COCl2_S1.wfn就将包含S1态的波函数信息。density关键词一定要加，否则wfn文件包含的是基态波函数，而不是你指定的第root态的。

做等值面图和做截面图的方法和前两节的过程完全一致，将c:\COCl2_z+300.wfn替换为c:\COCl2_S1.wfn就行了。得到的等值面图如下所示，isovalue=0.03。可以看到在体系从S0垂直激发到S1后，碳和氧的Px轨道上电子密度都有所增加，这是由于电子主要跃迁到了pi反键轨道所致。另外，氧的Py轨道上电子密度有所减小。



4.4 计算Fukui函数预测反应位点

Fukui函数是常用的预测反应位点的实空间函数，尤其适合软酸软碱之间的反应（静电势的分析方法更适合用于分析硬酸硬碱反应）。Fukui函数有三类，分别用于预测亲核、亲电、自由基反应位点。预测亲电反应的函数f- = ρ(N)-ρ(N-1)，其中ρ(N)是分子正常状态下的密度函数，ρ(N-1)是分子少了一个电子，或者说分子的净电荷为+1状态下的密度。f-在哪个原子附近数值越大，说明亲电反应越容易发生在哪个原子上。注意计算ρ(N-1)时几何结构必须使用计算ρ(N)时的，这不仅是绘制密度差图的要求，这也是Fukui函数的定义本身所要求的。

这里计算呋喃的f-函数的等值面图。首先计算正常状态下的波函数文件，Gaussian输入文件如下，结构已在b3lyp/6-31G*下优化好。
#p b3lyp/6-31g* out=wfn nosymm

opted b3lyp/6-31G*

0 1
C                 -0.66222300   -0.93877200    0.00000000
C                  0.69784300   -0.97448500    0.00000000
C                  1.13199200    0.39408700    0.00000000
C                  0.00000000    1.14881300    0.00000000
O                 -1.10670000    0.35104000    0.00000000
H                 -1.43048300   -1.69690900    0.00000000
H                  1.31938500   -1.85911500    0.00000000
H                  2.14981500    0.75868100    0.00000000
H                 -0.19078800    2.21116400    0.00000000

c:\furan.wfn
算完后把净电荷和自旋多重度从0 1改成1 2，c:\furan.wfn改成c:\furan+1.wfn，重算一遍得到N-1电子时的波函数文件。
之后过程和4.1节完全相同，即在Multiwfn里依次输入
c:\furan.wfn
5
0
1
-,c:\furan+1.wfn
1
2
算完后选-1预览等值面，适当调节isovalue，以使得等值面区域只集中在个别原子上。当前体系isovalue设为0.02时比较清楚，如下图所示，等值面只在氧的邻位上出现。这显示出氧邻位碳附近的f-函数值大于间位碳附近的，预测了呋喃的亲电反应会发生在氧的邻位碳上，这与实验结论是一致的。




5 绘制某个体系与组成它的各个片段的密度差图

这一节以水二聚体为例介绍怎么绘制某个体系与组成它的各个片段的密度差图。二聚体的几何结构直接从S22弱相互作用数据库获得，如下图所示



首先获得二聚体的波函数文件，输入文件如下（注1）
#P B3LYP/6-31+G** out=wfn nosymm

From S22

0 1
O  -1.551007  -0.114520   0.000000
H  -1.934259   0.762503   0.000000
H  -0.599677   0.040712   0.000000
O   1.350625   0.111469   0.000000
H   1.680398  -0.373741  -0.758561
H   1.680398  -0.373741   0.758561

c:\waterdimer.wfn
二聚体算完之后，将后三个原子从输入文件中删掉，并将输出波函数文件名改为c:\water1.wfn，然后计算。算完后再将后三个原子恢复回来，而将前三个原子从输入文件中删掉，将输出波函数文件名改为c:\water2.wfn，再次计算。此时water1.wfn和water2.wfn里面就分别是第一个水分子和第二个水分子的波函数信息。

启动Multiwfn，依次输入
c:\waterdimer.wfn
4 //绘制平面图
0 //设定自定义运算
2 //有两个波函数将从当前载入的波函数中减去
-,c:\water1.wfn
-,c:\water2.wfn
1 //电子密度函数。最终得到的密度就是waterdimer.wfn的密度减去water1.wfn的密度再减去water2.wfn的密度
2 //等值线图
0 //设定绘图平面的延展范围，默认的值4.5bohr对此体系偏小，一部分感兴趣的区域会被边框截断，故自行设定
6 //延展距离设定为6bohr，对于此体系可以保证感兴趣的区域都能出现在图中
直接按回车使用推荐的格点设定
4 //通过输入三个原子的序号定义作图的平面
1,2,3 //1,2,3号分别是第一个水分子的三个原子序号
马上密度差图就会蹦出来，如下图所示



从这幅图上能看出两个水分子形成复合物后电子密度是如何转移的，实线和虚线代表结合后密度增加和减少的区域。1区域所示的密度增加是由于氧的丰富的电子受到氢原子核的静电吸引而产生的明显极化。氢附近的电子由于受到另一个水的氧附近的电子密度的静电互斥和交换互斥，导致2区域电子密度被挤走，向着水的其余部分转移。


6 绘制体系与组成它的每个原子的密度差图

做完上个例子后，一定已经知道怎么做体系与组成它的各原子的密度差图了（下面称为变形密度），只需要在填入每个片段波函数文件名的地方改为填入每个原子的波函数文件名就行了。然而，原子的数目通常很多，依次计算原子的波函数文件并加入自定义运算列表里是很费事的，而且这样也忽略了原子密度球对称化问题（第8节将谈到）。为了使变形密度图能够方便地绘制出来，Multiwfn可以自动调用Gaussian生成原子波函数文件并作球对称化，并自动将其中心挪到各个对应的原子位置上，求密度并作差值，完全不用人工干预。

这里来绘制镁卟吩平面的变形密度等值线图。此体系的波函数文件MN.wfn还是以前Zhou2009传给我的，是在B3lyp/6-31G*下生成的，结构如下



首先确保Gaussian环境变量已经设好，否则Multiwfn在调用Gaussian时会出错。对于Win7，设定方法是进入“系统”-“高级系统设置”-“高级”标签页，点“环境变量”，在用户变量中点击“新建”，变量名填GAUSS_EXEDIR，变量值填Gaussian的安装路径，比如d:\study\g09w。设过一次之后以后就不用再重新设了。

启动Multiwfn依次输入
MN.wfn
4 //绘制平面图
-2 //获取变形属性（注2）。这里的“属性”是指接下来你所选的实空间函数，因为接下来要选电子密度函数，所以将得到变形密度。
b3lyp/6-31G* //设定计算原子波函数所用的方法和基组。
d:\study\g09w\g09.exe //Gaussian可执行文件的路径。注意可执行文件名是g09.exe而不是g09w.exe。如果你懒得每次做变形密度图时都重新输入一遍路径，那么就在settings.ini文件里将gaupath=后面空一格写上这个路径，并且这路径两边要用双引号扩住，比如gaupath= "D:\study\G09W\g09.exe"。当Multiwfn发现这个路径正确，就不会再让你手动输入路径了。
现在Multiwfn自动调用Gaussian计算每个原子的波函数文件，由于此体系只含C H N Mg这四种元素，所以会调用Gaussian四次。待算完后接着输入
1 //电子密度函数
2 //等值线图
敲回车，用默认格点设定
4 //由于MN.wfn中镁卟吩不是恰好处在某个XY或YZ或XZ平面上，所以这次还是通过三个原子来定义平面。
19,22,23 //用19、22、24号原子定义绘图平面。这三个原子位于镁卟吩的三个顶角，因此做出来的图会覆盖整个镁卟吩。（也可以像第5节那样通过设定延展范围来调节作图区域大小）
很快就得到下图



可以看到，在每个共价键之间都有明显的正值区域，表明共价键的形成伴随着电子密度在成键原子间聚集。N和Mg之间是离子键，如预期的，并没有电子密度以两原子之间为中心聚集。由N向着Mg方向凸起来的等值线区域表现的是N与C成键后，在那片区域凝聚起来的孤对电子。


7 Multiwfn处理原子波函数的内部机制

本节介绍Multiwfn处理原子波函数的内部机制，以便用户更灵活地、更方便地运用Multiwfn做变形密度图。下文所谓的“当前目录”是调用Multiwfn时所在目录，如果通过双击图标来启动Multiwfn，则当前目录就是指Multiwfn可执行文件所在目录。

在用户启动Multiwfn后选择了作图的类型，并选了-2以获取变形属性时，Multiwfn会执行以下步骤：
1 删除当前目录下wfntmp文件夹然后重建此文件夹。此文件夹用于储存计算原子波函数的临时文件。
2 将当前目录下atomwfn目录下所有波函数文件拷贝到wfntmp文件夹内。
3 检查wfntmp文件夹是否已经有了当前体系中所有元素的波函数文件。例如，如果发现此文件夹里有N .wfn文件（元素符号为单字符时文件名必须留一个空格），它就被认为是氮元素的波函数文件。元素的波函数文件里原子中心必须位于原点。
4 如果体系中所有元素的波函数文件都能在wfntmp文件夹内找到，就跳到第7步。如果少了任何一个元素的波函数文件，就让用户输入计算它们的理论方法和基组。然后检查settings.ini下gaupath的路径是否正确，如果不正确就让用户输入Gaussian的可执行文件路径。
5 在wfntmp目录下生成缺失的元素的Gaussian输入文件，调用Gaussian执行之，以得到它们的波函数文件。Gaussian任务的输出文件和波函数文件也都生成在wfntmp目录下。
6 对上一步生成的元素波函数文件做球对称化处理，这将会改写元素的波函数文件。
7 将每个元素的波函数文件根据体系中相应元素原子的位置复制成原子的波函数文件。比如体系中3、6、7号原子都是氮原子，则程序会将N .wfn里的原子坐标分别移动到3、6、7号原子位置上并复制一份，分别命名为N    3.wfn、N    6.wfn和N    7.wfn。
8 将所有原子波函数加入自定义运算列表里，运算符都是负号。表明它们的属性将依次从体系中减去。

如果觉得每次调用Gaussian重新计算各元素的波函数太麻烦，可以自己先手动算好，并做球对称化处理，改成 [元素名].wfn 格式后放到atomwfn目录下。由上述过程描述可知，Gaussian发现已经有了相应元素的波函数文件就不会再通过Gaussian计算了。在Multiwfn的examples文件夹里就有一个atomwfn文件夹，里面存的是前四周期在ROHF（主族）或UB3LYP（过渡金属）下结合6-31G*事先计算好的元素波函数文件（已球对称化）。将这个atomwfn文件夹移动到当前目录下，就可以省得再调用Gaussian了。

如果想直接得到特定方法和基组下已球对称化好的元素波函数文件，目的是将之放到atomwfn目录下以避免每次都调用Gaussian计算它们，那么可以随便找个包含那些元素的波函数文件载入进Multiwfn，选完绘图类型后选-2，使元素波函数文件在wfntmp下生成，然后将它们拷到atomwfn目录下即可。另外，如果想一次就把前四个周期的元素波函数文件都生成出来，那么可以载入examples目录下的genatmwfn.pdb文件，它包含了前4个周期所有元素，因此用-2选项时就会生成所有这些元素的波函数文件。

Multiwfn只能对前四个周期元素自动通过Gaussian计算波函数文件并球对称化。在提示输入理论方法和波函数时，如果只输入基组，比如6-31G*，那么就会默认用ROHF计算主族元素，用UB3LYP计算过渡金属元素，这也是推荐的方法。或许这样的理论方法和你的体系波函数文件所用方法不符，但实际上没什么问题，因为原子密度对理论方法依赖性并不大，对基组依赖性也不很大，所以examples\atomwfn目录下预置的那些元素波函数文件就足以用于做大部分体系的密度差图了，起码定性是正确的。

如果同时输入方法名和基组名，那么必须严格按照诸如B3LYP/6-31G**这样的格式，不要在前面加上RO、R或者U，因为Multiwfn对于过渡金属会自动加上U，对主族自动加上RO。方法名只能用HF或DFT（双杂化泛函除外），Multiwfn的球对称方法无法处理其它情况。如果在自动调用Gaussian计算波函数文件时出错（现象是突然退出），那么应检查wfntmp目录下自动生成的Gaussian输入文件的route section语法是否正确，以及Gaussian产生的out文件中是否出现错误，比如没能收敛（建议换个泛函）。

如果你的体系是混合基组，尤其是常见的使用赝势的情况，比如计算二茂铁（铁用lanl2TZ、其它原子用6-31G*），那么建议你自行在lanl2TZ下计算出铁的波函数文件（最好也自行对之作球对称化处理），改名为Fe.wfn后放到atomwfn文件下，然后按照正常步骤做变形密度图，填入基组时填6-31G*。这样Multiwfn就会调用Gaussian在6-31G*下生成Fe以外原子的波函数文件，而对Fe直接使用你提供的lanl2TZ下的波函数文件。


8 原子密度的球对称化问题

原子在自由状态下的密度是球对称的。变形密度是实际体系与组成它的各个原子在自由状态下密度的差值，这也就要求原子波函数文件对应的密度必须是球对称的。然而，对一些元素的一些状态，在一般量子化学计算中得到的密度却不是球对称的。比如基态为s2p2组态的碳原子，有两个p轨道上占了单电子，还有一个p轨道是空的，所以空着的p轨道方向上电子密度低，碳原子的密度成了扁椭圆的（而在现实中，我们无法观测到空p轨道朝向哪里，只能观测到平均状态--球对称的密度）。原子密度不满足球对称会造成变形密度出现不合理性。比方说CH4的变形密度，本应该是有Td对称性的，如果用了扁椭圆的碳的密度，就会导致变形密度的对称性被破坏，缺失了物理意义。再比如，基态氧分子是s2p4组态，有一个双占据p轨道和两个单占据p轨道，因此氧的电子密度也不是球形的，而是长椭圆形。假设要算O2分子的变形密度，那么计算氧原子的密度时应当让双占据的p朝向什么方向？让双占据p轨道顺着分子轴没什么道理，毕竟两个双占据轨道在一起不会成键；然而如果将它冲着垂直于分子轴的方向，那么变形密度就不是本应有的D∞h对称性了。

因此，生成变形密度时必须要求原子波函数对应的密度是球对称的。对于第一、二、五主族元素，惰性气体元素，半满或满占据的过渡金属元素，一般量化算得的基态波函数的密度本身就是球对称的，所以不必去做球对称化处理。而其它元素就必须考虑如何获得球对称的密度。一般有这么几种办法：
1 用密度是球对称的组态代替基态。比如C，可以用sp3代替sp2p2，由于C在分子中往往是sp3杂化，所以这么做是有道理的。然而对于基态组态为s2p1的B，就没有合适的生成球对称密度的组态了，如果把p轨道的那个电子激发到3s轨道上，显然密度会太过于弥散，与真实很不符。所以这种球对称方法并不普适。
2 用分数占据的DFT计算。在这种形式中，不要求每个轨道电子占据数都是整数，比如可以让B的三个p轨道都占1/3个电子，这保证了密度是球对称的。这种方法很好，但是分数占据的DFT计算并没有被广泛支持，Gaussian也不支持（听CCL上有人说也支持了，但我没找到用法，我也没看到最新修订版介绍里面提到过）
3 做MCSCF计算。最后得到的其中每个组态波函数的密度不是球对称的，但MCSCF波函数的密度会是球对称的，产生的自然轨道及占据数和分数占据的DFT的情况是相似的。然而做这种MCSCF会把问题复杂化。
4 用一些s型GTF或STO函数拟合非球对称的原子密度，之后由这些s型函数重现出的密度就是球对称的了。这种方法也相对麻烦，费时费事。

在Multiwfn里，对原子密度球对称化的方式简单高效，对不同类型元素用不同处理方法。对于第三主族，将限制性开壳层计算中那个单占据p轨道旋转向其它两个方向并且都复制一份，然后将这三个轨道的占据数都设成1/3；对于第四主族，提供了两种球对称化方法，默认的是使用sp3组态，另一种做法和处理第三主族方式类似，可以通过将settings.ini里的SpherIVgroup设为1来启用；对于六、七主族，由于每个p轨道都是占据轨道（只是占据数不同），而占据轨道间形状差异很小（远小于占据轨道和空轨道的差异），因此直接将p轨道上全部电子在三个p轨道上平摊；对于Sc、Ti、V用的方法和第三主族类似，是将单占据的d轨道向其它方向旋转并复制（Sc的4s轨道被极化得较厉害，所以也通过旋转复制的策略削弱它造成的各向异性）；对于Fe、Co、Ni，将d轨道上beta电子都均摊到5个alpha的d轨道上。

注意在过渡金属原子波函数计算时，所用的理论方法必须能给出正确的4s和3d能级顺序，否则Multiwfn没法正确判断哪些是d轨而无法正确地球对称化。一般的泛函如B3LYP没问题，不能用HF，给出的顺序是错的。Multiwfn对第六、七主族和过渡金属的处理实际上并没有使之密度精确球对称化，但是密度的各向已经已经被削弱到可以忽略的程度。如果不希望Multiwfn自动对元素的波函数文件进行球对称化处理，可以将settings.ini里的ispheratm设为0来关闭之。

==============================
注1：关于应当用什么泛函生成弱相互作用体系的密度（假设都使用相同的可靠的几何结构），我想在这里表达一下我的观点。目前的泛函全都是近似的，而且不少都没有清楚的物理意义。在完全真实的密度下用这些泛函得不到真实的能量，而能够得到较好的能量/几何结构的泛函产生的密度未必比那些能量/几何结构预测得不好的泛函更好。究竟用什么泛函产生密度更合适，这没有确切结论，如果想获得在原理上很可靠的密度，那就去做高精度post-HF。但所幸，不同泛函之间产生的密度的差异远不像能量/几何结构差异那么大，尤其是对于弱相互作用体系来说。比如，就连公认对氢键等弱相互作用体系描述不好的B3LYP也能获得和MP2很相近的密度（见约化密度梯度分析弱相互作用的原文JACS,132,6498）。所以，第五节就用B3LYP来获得水二聚体的密度。M06-2X对弱相互作用的能量/结构预测上表现得比B3LYP好很多，会有不少人认为使用M06-2X产生密度更合理，表面上来看用它也的确不会有不妥之处，只不过并不能直接论证它的密度会比B3LYP更接近真实。而对于DFT-D类的泛函，对弱相互作用也普遍表现得较好，这是由于它们在原有泛函能量表达式中额外引入了经验校正项的结果，而单电子有效势算符并没有改变。故产生的密度实际上和原有泛函是一致的。（有一些DFT-D泛函在引入经验校正项的时候还重新拟合了原有泛函中的参数，这类泛函在原理上我不建议使用。因为经验校正项不光描述远程相关作用以改进弱相互作用预测能力，也会描述一定程度中程的相关作用，而原有泛函也描述了中程相关。因此重新拟合出的参数会弱化原有泛函展现的中程相关，这将可能导致产生的密度并没有充分体现出电子相关效应的影响，还不如原有泛函的密度更可靠。）

注2：如果选了-1，就会得到promolecule属性，它是体系中每个原子在自由状态下属性的叠加。也就是说，变性属性=体系实际属性-promolecule属性。如果实空间函数选的是电子密度，就将得到promolecule密度，它代表了分子结构刚刚形成，但是电子密度还没来得及弛豫时的密度。生成promolecule属性也需要生成原子波函数文件，和生成变性属性经历的流程是完全一致的。

[ Last edited by sobereva on 2011-12-24 at 17:05 ] 返回小木虫查看更多