使用Multiwfn绘制原子轨道图形、研究原子壳层结构及相对论效应的影响

使用Multiwfn绘制原子轨道图形、研究原子壳层结构及相对论效应的影响文/Sobereva2012-Jul-91前言这个帖子主要介绍怎么用Multiwfn程序（）结合Gaussian绘制各种类型的原子轨道图形，包括角度和径向部分，在绘制过程中能加深一些对原子轨道的理解，如原子轨道间的正交性和钻穿效应。本文绘制轨道并不是像一般教材中通过原子轨道波函数的解析形式来绘制的，解析的方式可以用matlab、mathematica等程序绘制，本文是通过Multiwfn靠Gaussian输出的单原子体系波函数信息绘制的。在绘制过程中可以使没用过Multiwfn的人熟悉Multiwfn的基本绘图操作，对于有一定经验的用户也能学到一些特殊技巧。文中还将利用Multiwfn简要讨论相对论效应对轨道径向分布产生的影响，读者可以同时了解到在Gaussian中使用全电子标量相对论计算的基本方法。最后还将通过绘制各种实空间函数展现原子各个主层特征。本文介绍的方法和作出来的图对于讲授结构化学课程的老师也我想比较有用，很适合向学生们展示一些基本概念。本文用的Multiwfn为2.4版，Gaussian为G09A02。实际上，原子轨道只有对于类氢原子体系（一个核+单个电子）才是物理意义严格的，对于多电子原子体系，原子轨道模型只是近似的描述，但还是很合用的。类氢原子轨道波函数是径向部分波函数与角度部分（球谐函数）的乘积。比s角动量更高的原子轨道有的角度部分是复数，复数型原子轨道难以图形表示，用起来也不方便，因此一般都是将复数型原子轨道线性组合成实数型来用（它们将不再是Lz算符的本征函数而没法讨论磁量子数）。本文说的原子轨道都是指实数型原子轨道，教科书上的原子轨道图形也一般是实型的。而本文所谓的真实原子轨道，则是指实数型的类氢原子轨道。2绘制s,p,d,f,g角动量原子轨道的角度部分图形这里我们先不考虑径向部分，假定是个任意的常数，这里先来通过绘图将s,p,d,f,g角动量原子轨道的角度部分表现出来。s,p,d角动量的原子轨道图形想必大家已很熟悉，但是f、g的图形可能不少读者还不怎么印象深刻，此节将绘制它们。我们先建立一个Gaussian输入文件，内容如下。使Gaussian运行它后每个“分子轨道”都对应一个原子轨道，这样用一般方法观看分子轨道就等于观看原子轨道了。%chk=c:\gtest\atom.chk#phf/genpop=fullguess=(cards,only,save)Atom11HH0S11.00.11.P11.00.11.D11.00.11.F11.00.11.G11.00.11.****25(f2.0)-11.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.0.1.0这个输入文件看起来可能觉得比较古怪，这里进行解释。Gaussian用的是高斯函数作为基函数，高斯函数又细分为笛卡尔型高斯函数和球谐型高斯函数，后者的角度部分和真实原子轨道的角度部分是一一对应的，所以我们应当用球谐型高斯函数。对于此例自定义基组的情况，默认就是用球谐型高斯函数，不必手动加5d7f关键词。对于这个问题的更细致讨论见《谈谈5d，6d型d壳层波函数与它们在Gaussian程序中的标识》（）和《球谐型与笛卡尔型Gauss函数的转换关系》（）这个体系电子数为0，即质子体系。通过自定义基组方式，给这个质子加上从s到g角动量基函数各一个，比如D11.00.11.就代表加上指数为0.1（因为我们这一节忽略原子轨道径向部分，所以此值是随意取的）的d壳层的共5个基函数，且每个基函数都只含一个高斯函数。当前体系总基函数数目是1+3+5+7+9=25个。但是我们并不能在自定义基组后直接就这么计算，否则从pop=full输出的信息会看到“分子轨道”中存在基函数的混合，绘制“分子轨道”图形就不能对应于原子轨道图形了，所以我们必须写上guess=(cards,only,save)并自行设定初猜。其中cards代表自行从输入文件后面读取初猜的分子轨道系数，而不用程序自动给的初猜；only代表不做迭代，否则又会引起基函数的混合；save代表将初猜波函数写入chk文件中（默认对于only任务是不写入）。25(f2.0)代表自行写的初猜信息是每行25个值（恰好每行代表一个分子轨道的25个基函数的系数，比较清晰易读），而且每个值用fortran语言的f2.0浮点格式（占两个位置，即一位整数和一个小数点符号自身的占位，比如3.1415就会表示成3.，这样虽然精度极低但对于目前问题最适合，十分紧凑）。后面的-1代表用自己设的初猜替换所有轨道的初猜。接下来就是自己设的初猜轨道系数了，我们要让i号“分子轨道”正好对应i号基函数，因此只有对角线系数为1其余为0。输入文件末尾的一个0代表自定义轨道初猜信息已经写完了。pop=full其实没有意义，只是这将便于从输出文件中检查每个“分子轨道”的组成，看看是否如预期的每个“分子轨道”只在一个基函数上有值且系数为1，即基函数没有发生混合。用Gaussian计算这个任务，然后将chk用formchk转换为fch文件，然后启动Multiwfn，输入fch文件的路径，然后选0，通过点击弹出的图形界面右下角的分子轨道标号就可以看相应分子轨道，在此例中即原子轨道的图形了。在点击一个轨道标签后，也可以用键盘的上下键切换，在浏览一批轨道时比用鼠标点更为方便。建议在启动Multiwfn之前将Multiwfn目录下的settings.ini里的aug3D参数从默认的6调大到10，然后保存，这样可以避免以较小isovalue显示外层原子轨道等值面时在边缘被默认的比较窄的格点数据空间范围所截断。在Multiwfn中观看原子轨道时建议将isovalue从默认的0.05减小到0.03，否则原子轨道的一些特征表现不出来。这里我们随便选一个轨道，比如15号轨道，从pop=full给出的信息中看到这个轨道对应于F+3，注意这个绝非代表磁量子数为3的f轨道。根据《球谐型与笛卡尔型Gauss函数的转换关系》一文提供的信息，我们知道这个轨道用笛卡尔形高斯函数表示为√(5/8)*XXX-3/√8*XYY，可想而知这个轨道是处在XY平面上的。下图将f原子轨道和g原子轨道等值面图形汇总：7个f（10~16号轨道）9个g（17~25号轨道）下面我们用Multiwfn作轨道的平面图，这里以第15号分子轨道为例，在Z=0的XY平面上作图最能充分表现它的特征。如果已经打开了Multiwfn，先将它关闭，然后在settings.ini里将idelvirorb值从默认的1设为0，然后保存。默认的1代表在一些可能比较耗时的任务中删掉fch文件中的前10个虚轨道以外的虚轨道以节约计算时间，但是删了虚轨道就达不到本文目的了，所以将此参数改为0并保存，以避免Multiwfn这么做（以后的Multiwfn版本中这个设定有可能还会改变，注意参见手册和此参数在settings.ini里的注释）。启动Multiwfn，依次输入o//直因为之前已经输入过一次当前体系的fch文件的路径了，所以这次直接写字母o就可以打开上一次载入的文件4//作平面图4//要作的函数为轨道波函数值15//作第15号轨道1//填色图直接敲回车用默认的格点设定200,2000//设定作图延展距离，默认的延展距离对于此例偏小，会看到轨道外部被截断8//设延展距离为8Bohr1//作XY平面0//Z=0图像立刻弹出来，但是基本是绿色，看不出什么特征，这是因为默认的色彩刻度范围太大，不适合当前情况。遂在图上点击右键关闭之，选1，输入-0.09,0.09修改彩刻度上下限。如果想让图像上同时出现等值线，可以再选选项2。最后选-1重新作图，得到下图（关闭图像后选0可以将图保存为当前目录下的文件名以DISLIN开头图形文件）我们接下来绘制这条轨道的电子密度的地形图+投影图。虽然目前它是空轨道（其它轨道也都是空着的），但是我们可以通过修改波函数信息来设它为双占据。选-5从刚才的后处理界面退回到主界面，依次输入6//修改波函数26//设定轨道占据数15//选15号轨道2//占据数设为2（双占据轨道）q//退回-1//退到主菜单。由于其它轨道都是空轨道，目前只有15号轨道有电子占据，因此照常作密度图时，就等于只作15号轨道的密度图了4//作平面图1//电子密度5//地形图+投影图直接敲回车用默认格点设定100,1001//XY平面（由于之前已经调整了延展距离，所以这次不用再设一遍）0//Z=0图上几乎一片空白，显然默认设定不适合当前体系，所以我们还要调整。在Multiwfn目前版本绘制地形图时，地形图的Z轴范围不允许改变。明显是因为当前的Z轴范围过大（-4到