python态密度曲线绘制_使用Multiwfn绘制态密度(DOS)图考察电子结构

使用Multiwfn绘制态密度(DOS)图考察电子结构

文/Sobereva@北京科音

First release: 2019-May-14  Last update: 2020-Nov-17

0 前言

基于主流量子化学程序产生的波函数信息，使用Multiwfn可以对孤立体系非常简单、快速、灵活地产生质量非常理想的态密度图，对于分析化学体系的电子结构特征很有帮助。本文将介绍各类态密度图(Total DOS、Partial DOS、Overlap population DOS、Local DOS)的原理，以及它们在Multiwfn程序中的绘制方法。DOS虽然是第一性原理计算范畴的重要概念，但本文对DOS的绘制只涉及量子化学计算范畴，Multiwfn也不支持基于第一性原理计算产生的数据绘制DOS。

Multiwfn从很早期的版本开始就支持DOS的绘制，后来又不断进行了改进和完善，如今已非常强大，显著强于其它任何能绘制DOS的程序。Multiwfn的DOS绘制功能已经被大量发表的文章所使用，比如J. Comput. Chem., 38, 1574 (2017)、RSC Adv., 5, 78192 (2015)、Nature Communication, 8, 14551 (2017)、Phys. Chem. Chem. Phys., 19, 23373 (2017)、J. Phys. Chem. C, 119, 8349 (2015)、J. Phys. Chem. A, 121, 4009 (2017)、Chem. Eur. J., 24, 17046 (2018)等等等等，这些文章也可以当成例子。

本文的内容非常深入细致，如果你很急着作出TDOS/PDOS/OPDOS图，可以先直接参考本文第3节的例子，之后若想全面了解Multiwfn的绘制DOS的功能，再看本文其它部分。本文内容对应Multiwfn官网上的最新版本，不要用太老的版本，否则情况和本文可能明显不符！Multiwfn可在其官网http://sobereva.com/multiwfn免费下载，相关基础知识看《Multiwfn入门tips》(http://sobereva.com/167)和《Multiwfn FAQ》(http://sobereva.com/452)。

1 态密度的物理意义与种类

态密度(Density-of-states，DOS)是一个以能量为变量的函数。DOS(E)代表在能量为E的位置，单位能量区间内态的数目。“态”具体指什么，看这个词具体用在什么地方。

一般说DOS这个词较多的场合是第一性原理领域，此时的态通常指的是解KS-DFT方程得到的单电子轨道。第一性原理研究的主要是周期性体系，有k点的概念，不同k点处轨道能级也不同。第一性原理计算可以得到物理意义严格的DOS，即通过对k空间进行积分得到。

量子化学主要研究的是孤立体系(分子、团簇等)，当使用DOS考察体系电子结构时，态指的是分子轨道(可以由HF、半经验、KS-DFT等各种基于单电子近似的理论方法计算产生)。孤立体系计算的时候没有k点的概念，算出来的轨道能级都是离散分布的，因此DOS这个词对于孤立体系来说本身没什么意义，画出图来的话，就是下面这样，对应的公式也一起给出了(为了避免和后文说的其它类型的DOS混淆，这里用了Total DOS (TDOS)这个词)。

上图中虚线位置是HOMO的位置，公式里ε是轨道能量，δ是狄拉克δ函数(不懂的话Google)。从上图可见，孤立体系的TDOS图本身只不过相当于在分子轨道能级位置画一个竖线而已，这种图一点意思也没有，没法传递出比轨道能级分布本身更多的信息。

如果通过一些函数将上图中这些离散的竖线依次进行人为地展宽然后再累加，从图中将可以一目了然地看到不同能量区域内轨道分布的密集程度，如下图所示。注意这种展宽纯粹是人为的，仅仅是为了便于直观分析而做的，这样得到的DOS曲线并不是像第一性原理计算那样是以严格方式得到的。

不同的展宽函数把离散的竖线展宽出的峰的形状不同。绘制DOS时展宽函数一般用Gaussian函数，在Multiwfn里也允许用Lorentzian函数，或者二者混合产生的Pseudo-Voigt函数。函数表达式这里就不列出了，读者可以在Multiwfn手册里的介绍DOS绘制功能的3.12节中看到。不管用的是哪个展宽函数，都会牵扯到半高全宽(FWHM)这个参数，即竖线展宽出来的峰在高度一半的位置峰的宽度。显然FWHM设得越大峰越宽，曲线看起来就越平滑、越连续，而FWHM越小则对峰的分辨率越高，绘制上图用的FWHM为0.05 a.u.。FWHM这个值没有什么物理意义，设多大是具有主观任意性的，一般的设定标准是让得到的图适合、便于说明自己想说明的问题。下图分别是FWHM改为0.1 a.u.和0.02 a.u.时候得到的图。

在《使用Multiwfn绘制光电子谱》(http://sobereva.com/478)一文中，得到光电子谱(PES)的方式在本质上和上述绘制TDOS其实如出一辙。

还可以定义Partial DOS (PDOS)，可以叫做“局部DOS”或“分数DOS”，它体现的是特定的片段对TDOS贡献的曲线。若恰当定义片段，通过PDOS图就可以较好地把握不同能量区间的轨道的本质、主体构成。显然，如果定义的所有片段的并集等于整个体系的话，那么每个片段的PDOS曲线累加在一起正好就是TDOS。对于某个片段A，其PDOS定义如下。其中F是展宽函数，Θ是轨道成份。

计算轨道成份的方法很多，详见《谈谈轨道成份的计算方法》(http://sobereva.com/131)和《分子轨道成分的计算》(http://sioc-journal.cn/Jwk_hxxb/CN/abstract/abstract340458.shtml)里的介绍。在绘制DOS的功能中，Multiwfn支持Mulliken、SCPA、Hirshfeld和Becke方法计算轨道成份。

重叠布居DOS(Overlap population DOS，OPDOS)对于考察片段间的相互作用比较有用。笔者在《Multiwfn支持的分析化学键的方法一览》(http://sobereva.com/471)里专门介绍过什么叫Mulliken键级，这和重叠布居(Overlap population)是一码事，可以体现原子间是成键作用还是反键作用，没看过此文的话一定要仔细看一下。文中还提到，重叠布居可以精确分解为各个轨道的贡献，将这个思想和DOS图搞到一起，就是OPDOS的定义了。下式是片段A与B之间的OPDOS的表达式

式中的那个Χ项，也等价于假设i轨道是占据的，然后这个轨道对片段A和B之间的重叠布居的贡献量除以轨道占据数。在某个能量范围，若OPDOS数值明显为正，就说明这个能量范围里的轨道如果被电子占据了的话，会对两个片段间的结合起到积极贡献，换句话说，这些轨道对这两个片段而言是成键轨道；若数值为负，说明这些轨道如果被电子占据了的话，将对片段间的结合产生不利贡献，相当于起到反键轨道作用。

虽然OPDOS确实能展现片段间相互作用的情况，但是由于Mulliken键级本身对键的强度反映不好，所以也别太把OPDOS当回事。而且，有弥散函数的情况下，OPDOS结果会很离谱，可能数据完全没有意义。所以，若计算波函数时有弥散函数则不要讨论OPDOS，若必须讨论OPDOS且之前的计算用了弥散函数，应把弥散函数去掉算个单点重新得到波函数后再做分析。Multiwfn是分析成键问题的百宝箱，里面有大把方法可以很好地讨论化学键特征，所以千万别拘泥于拿OPDOS说事，强烈建议仔细看《Multiwfn支持的分析化学键的方法一览》(http://sobereva.com/471)。另外，之前笔者还见到有人试图拿两个片段各自的PDOS曲线讨论它们的成键，以为在某些能量区间它们同时较大就说明二者间有成键作用，这明显毫无意义，且不说它们是相位相同方式重叠还是相位相反重叠，首先这种情况都没法体现这两个片段间存在轨道重叠，因为PDOS图根本都体现不出片段的轨道分布在三维空间中的什么区域。

还有一种DOS叫Local DOS，这个留到本文第4节专门介绍。

我们来看一个典型的DOS例子，是Multiwfn的原文J. Comput. Chem., 33, 580 (2012)里的图，将二茂铁的TDOS、PDOS、OPDOS绘制在了一起。为了便于说清楚轨道特征，一部分轨道的等值面图是用Multiwfn的主功能0绘制出来后手动ps到恰当位置的。

图中左侧是对应TDOS、PDOS曲线的轴，右侧是对应OPDOS曲线的轴。DOS或OPDOS的数值的绝对大小是我们不需要关心的，刻度轴的标签都抹掉也无妨，需要关心的仅仅是不同能量区间之间曲线的相对高度。为了便于直观考察不同能量区间内轨道的本质，绘制时定义了三个片段：

• 片段1：Fe的所有轨道

• 片段2：碳的垂直于茂环的p轨道。由于茂环在XY平面上，因此它们也即碳的pz轨道，合在一起相当于茂环的pi轨道

• 片段3：茂环的碳的其它轨道，即s、px、py

由图可见，对于能量比较低区间的轨道，诸如-0.53 a.u.的那一批，从图上立马就可以判断出它们基本都是由碳的s、px、py轨道贡献的，结合轨道等值面图也看出确实如此。而对于能量为-0.25 a.u.附近的轨道，蓝色曲线和红色曲线都较高，说明这些轨道应当是由茂环的pi轨道与铁的轨道混合产生的，确实从轨道等值面上也可以证实是由铁的d轨道与茂环的pi轨道以相位相同叠加组合出的。在HOMO轨道位置，Fe片段的PDOS与TDOS非常接近，说明HOMO主要是由Fe的轨道构成的，茂环的参与较弱。我们再看绿色的OPDOS曲线。在比如-0.25 a.u.附近，绿色曲线明显为正，说明这些轨道的存在对于茂环与铁的结合起到积极作用，这也容易理解，毕竟是相位相同叠加，所以对于此二者而言起到的是显著的成键轨道作用。对于非占据轨道区域，可见OPDOS曲线都为显著的负值，这也就是说，如果它们被电子占了，倾向于让茂环与铁的结合散架。确实从LUMO的轨