谈谈怎么考察、计算、分析化学体系的电子密度

谈谈怎么考察、计算、分析化学体系的电子密度
How to investigate, calculate and analyze electron density of chemical systems

文/Sobereva@北京科音  2024-Jun-13

1 前言

笔者经常看到有人问“请问如何使用高斯计算电子云密度”之类的问题。之前我在《Multiwfn FAQ》（http://sobereva.com/452）的Q45里专门回答过这种问题，由于这个话题较大，现在觉得值得再多说说，所以写成个单独的短文。

首先要明确，电子密度（electron density）或者所谓的电子云密度，是一个三维空间函数，没有唯一的考察和分析方法。有人问“怎么计算原子的电子云密度”这种问题显得莫名其妙，波函数分析领域有不同方法如Hirshfeld、Hirshfeld-I、MBIS、Voronoi、AIM等都可以把化学体系整个三维空间划分为独立的原子空间，每个原子的空间里每个位置都有相应的电子密度值，你到底想计算哪里的？再比如“请问如何使用高斯计算电子云密度”这种提问也十分不明不白，本来Gaussian就不是专门干这个的，而且你想以什么形式展现？计算目的何在？明显基本情况都没说清楚，显然只有对电子密度的基本概念完全稀里糊涂的人才会这么问。

2 获得电子密度的方法

在说怎么考察电子密度之前先简谈一下电子密度数据能通过什么方式获得：
(1)实验方法：通过电子衍射、X光衍射、中子衍射，都可以测定出电子密度在三维空间中的分布，从而给出格点数据，描述三维空间里均匀分布的各个格点位置的电子密度值。在测定的电子密度解析度足够高的情况下还能进一步确定原子坐标，看《实验测定分子结构的方法以及将实验结构用于量子化学计算需要注意的问题》（http://sobereva.com/569）。
(2)理论计算方法：对孤立体系用量子化学程序（如Gaussian），以及对周期性体系用第一性原理程序（如CP2K），都可以计算电子的波函数，并进而得到电子密度。理论方法、基组等因素都直接影响得到的电子密度质量。理论方法对基态的电子密度分布描述的精度对比可参考J. Chem. Theory Comput., 13, 4753 (2017)。在常用的基组范畴中对于价层电子的描述质量，def2-SVP或6-31G**算是可以接受的底限，6-311G(2d,p)或def-TZVP算是不错，def2-TZVPP就算很好了。对于CP2K常见的基组来说，DZVP-MOLOPT-SR-GTH是底限，TZVP-MOLOPT-GTH算是不错，TZV2PX-MOLOPT-GTH就算非常好了。

除上述外，还有特别简单的构造电子密度分布的方法，称为准分子密度（promolecular density）。这就是把各个原子在孤立状态下的电子密度分布根据当前体系的核坐标进行叠加得到的，这样的电子密度相当于没有考虑由于原子间相互作用而造成的电子转移和极化情况的“零阶”的电子密度。这种电子密度不在下文讨论范畴中。

量子化学或第一性原理程序在计算后，往往可以导出电子密度的格点数据。例如CP2K可导出电子密度的cube文件，记录了均匀分布的格点上的电子密度，cube文件的介绍见《Gaussian型cube文件简介及读、写方法和简单应用》（http://sobereva.com/125），它可以直接用于VMD、VESTA等程序观看电子密度等值面等目的，也被一些程序支持用于定量分析。很多程序也可以给出记录波函数信息的文件，如Gaussian可以给出fch/wfn/wfx文件，ORCA和CP2K可以给出molden文件，都可以被Multiwfn等波函数分析程序用来做电子密度的相关分析，参考《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）和《详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件》（http://sobereva.com/651）。

3 考察电子密度的方法

波函数分析领域博大精深，电子密度的分析仅仅属于波函数分析的一小部分而已，下面会按类别把分析、表征电子密度的方法做简要的提及，这些分析都是Multiwfn完美支持的，不仅可以用于孤立体系也可以用于周期性体系的分析。Multiwfn（http://sobereva.com/multiwfn）是波函数分析领域功能最丰富、最易用、最流行的程序之一，不了解者建议阅读《Multiwfn入门tips》（http://sobereva.com/167）和《Multiwfn FAQ》（http://sobereva.com/452）。

(1)计算某个点的电子密度。用Multiwfn载入波函数文件后在主功能1中输入要考察的坐标即可得到这个位置上电子密度在内的一大堆属性。这有一些实际用处，比如考察一些芳环体系的各个原子在分子平面上方1埃位置的电子密度可以用来考察亲电反应的优先位点，看比如《亲电取代反应中活性位点预测方法的比较》（http://www.whxb.pku.edu.cn/CN/abstract/abstract28694.shtml）中的应用。

(2)图形化考察电子密度在三维空间的分布。这主要用于定性考察电子密度分布情况，包括绘制电子密度在某条路径上的曲线图、特定平面上的图（填色图、等值线图、地形图、梯度线图、向量场图等）、等值面图（适合展现立体分布情况），可以分别用Multiwfn的主功能3、4、5实现，分别把Multiwfn手册4.3、4.4、4.5节的例子看过一遍就懂了。

注意一般来说直接对电子密度作图没什么意义，因为电子密度分布是很“单调”和“乏味”的，即一般在原子核位置是个电子密度极大的峰，向四周以指数型下降，在这样的图上很难观察到化学上感兴趣的信息。但如果你考察的是价层电子密度，则可以讨论很多化学问题，十分建议参看笔者发表的《Revealing Molecular Electronic Structure via Analysis of Valence Electron Density》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252），里面有大量生动的例子。价层密度的绘制参考Multiwfn手册里4.6.2节的例子。

值得一提的是，电子密度的0.001 a.u.等值面有特殊的意义，通常被视为气相下分子的范德华表面，里面包围了体系大约98%的电子。电子密度的0.0015 a.u.等值面还可以用来定义动力学直径，见《使用Multiwfn计算分子的动力学直径》（http://sobereva.com/503）。

Multiwfn还能只绘制特定轨道贡献的电子密度，看《在Multiwfn中单独考察pi电子结构特征》（http://sobereva.com/432）。

(3)计算原子电荷（atomic charge）。原子电荷相当于原子的核电荷数减去原子带的电子数，这可以非常清晰、定量地展现实际化学环境中各个原子带的净电荷量。这种分析在文献中被普遍使用，非常有实用价值，比如《18碳环等电子体B6N6C6独特的芳香性：揭示碳原子桥联硼-氮对电子离域的关键影响》（http://sobereva.com/696）一文通过原子电荷考察了B9N9和B6N6C6的原子带电情况。我强烈建议阅读《一篇深入浅出、完整全面介绍原子电荷的综述文章已发表！》（http://sobereva.com/714）和《原子电荷计算方法的对比》（http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml）以全面了解原子电荷的各种相关知识和计算方法的特点及差异。可以通过Multiwfn的主功能7的各种选项以ADCH、Hirshfeld-I、RESP、EEM、AIM等许多不同方法来计算原子电荷，把手册4.7节的例子好好看一遍就能很快掌握了，对于周期性体系的原子电荷的计算我还专门写了文章《使用Multiwfn对周期性体系计算Hirshfeld(-I)、CM5和MBIS原子电荷》（http://sobereva.com/712）。

将某个片段中的所有原子的电荷加和就得到了片段电荷（fragment charge），可以确切了解片段的带电情况、分析电子转移量。在Multiwfn中计算片段电荷很方便，在主功能7里选择选项-1，输入片段里原子的序号来定义片段，之后再照常选择一种方法计算原子电荷，则原子电荷输出后就会输出片段电荷。之前还录过演示视频《使用Multiwfn计算分子的某个片段的电荷》（https://www.bilibili.com/video/av26312703/）。

还有个概念叫布居分析（population analysis），所谓的布居就是指某个对象带的电子数。除了计算原子带的电子数外，在Multiwfn的主功能7中基于波函数信息还可以用Mulliken、SCPA等方法给出各个基函数、壳层、原子轨道上带的电子数，参考比如《利用布居分析判断基函数与原子轨道的对应关系》（http://sobereva.com/418）。非常灵活的Multiwfn还可以给出各个原子上带的特定类型电子（如sigma、pi电子）的量，参考《使用Multiwfn基于Hirshfeld-I划分计算特定类型电子在各个原子上的分布量》（http://sobereva.com/697）。

(4)电子密度差（有个很难听、我很讨厌的别称叫差分电荷密度）。这是在不同状态或不同体系间将电子密度相互求差得到的，由此可以展现出单个体系电子密度分布无法展现出的信息，在定量分析电子转移、体系对外场的响应、考察共价键的形成情况等方面有重要价值。电子密度差有多种类型，见《使用Multiwfn作电子密度差图》（http://sobereva.com/113）中的说明，里面还有在Multiwfn中的绘制过程的介绍。电子密度差在周期性体系的研究中也很重要，而且对一些体系将密度差转化为电荷位移曲线后在分析讨论时可以更为清晰准确，看《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（http://sobereva.com/638）。

预测反应位点常用的福井函数和双描述符本质上也是电子密度差，看《概念密度泛函综述和重要文献合集》（http://bbs.keinsci.com/thread-384-1-1.html）和《使用Multiwfn计算双描述符势预测反应位点》（http://sobereva.com/708）里面提到的相关资料。分析（超）极化率本质常用的（超）极化率密度也是基于密度差所定义的，看《使用Multiwfn计算（超）极化率密度》（http://sobereva.com/305）。

(5)一些与电子密度分布密切相关的量，例如：
• 电偶极矩、电多极矩：体系的电荷分布可以做电多极展开成为单极矩+偶极矩+四极矩+八极矩...的描述形式，其中每个部分都可以分为核电荷与电子密度各自贡献的部分，Multiwfn手册3.18.3节有完整的计算公式和相关说明。偶极矩可以体现正电中心和负电中心间的分离程度，四极矩可以体现电荷分布偏离球对称的程度。在Multiwfn中载入波函数文件后，选择主功能300的子功能5即可计算它们。Multiwfn还可以计算各个原子和片段的偶极矩和多极矩，看《使用Multiwfn计算分子片段的偶极矩和复合物中单体的偶极矩》（http://sobereva.com/558）。

• <r^2>：这衡量电子分布的延展程度，见《电子空间范围<r^2>和电子径向分布函数的含义以及在Multiwfn中的计算》（http://sobereva.com/616）。

• 电子密度的径向分布函数：这可以考察相对于某个点不同距离的单位厚度壳层内的电子数目的相对差异，在http://sobereva.com/616里也有介绍。

(6)电子密度的拓扑分析。任何三维函数都可以在Multiwfn的主功能2中做拓扑分析，从而得到它们的各种临界点（梯度为0的点）的坐标，以及相应位置的各种实空间函数的值。做电子密度的拓扑分析是Bader提出的非常知名的Atom-in-molecules (AIM)分析中用到的关键分析手段，此分析得到的键临界点（BCP）通常是原子间相互作用路径上最有代表性的一个点，通过这个位置的各种函数值可以讨论很多问题。例如《透彻认识氢键本质、简单可靠地估计氢键强度：一篇2019年JCC上的重要研究文章介绍》（http://sobereva.com/513）指出利用BCP处的电子密度可以很好地预测氢键强度，Multiwfn手册4.2.1节和《计算分子内氢键键能的几种方法》（http://sobereva.com/522）都有具体例子。再比如，AIM分析得到的环临界点上的电子密度和这个环的芳香性有密切联系，看《衡量芳香性的方法以及在Multiwfn中的计算》（http://sobereva.com/176）中的相关部分。关于更多的AIM和临界点的相关知识参考《AIM学习资料和重要文献合集》（http://bbs.keinsci.com/thread-362-1-1.html）和《使用Multiwfn做拓扑分析以及计算孤对电子角度》（http://sobereva.com/108）。

(7)电子密度的盆分析。盆分析是指用特定实空间函数定义的零通量面将三维空间划分为一个个独立的子区域，这被成为盆，然后在盆中对特定函数进行积分，就可以考察很多问题。对电子密度就可以进行盆分析，得到的叫原子盆，在里面积分电子密度对应的就是AIM方法定义的原子带的电子数，用Multiwfn的主功能17可以实现，看《使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析》（http://sobereva.com/179）和手册4.17.1节。对价层电子密度或ELF函数进行盆分析，然后在盆里积分电子密度，能够获得成键区域电子数、孤对电子数等信息，看《Revealing Molecular Electronic Structure via Analysis of Valence Electron Density》（http://www.whxb.pku.edu.cn/EN/10.3866/PKU.WHXB201709252）中的讨论。Multiwfn甚至还可以对密度差做盆分析然后在盆内积分密度差，这可以定量考察特征区域电子密度的变化量，http://sobereva.com/179和手册4.17.4节都有例子。孤对电子的ELF盆、电子密度=0.001 a.u.等值面内和ELF=0.5等值面内交集区域内的电子数称为high ELF localization domain population (HELP)，是孤对电子数的一种特殊定义，与电离能、前线轨道能量等不少体系的属性有联系，参看手册4.17.8节。

(8)电子密度的衍生函数。有很多函数衍生自电子密度，它们的分析也是波函数分析的重要组成部分、有巨大实际意义。比如《使用Multiwfn做IGMH分析非常清晰直观地展现化学体系中的相互作用
http://sobereva.com/621》、《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》（http://sobereva.com/598）、《一篇最全面介绍各种弱相互作用可视化分析方法的文章已发表！》（http://sobereva.com/667）里介绍的IGMH和IRI分析就是基于电子密度及其导数做的，在分析化学体系中的相互作用方面已经非常流行。再比如电子密度的拉普拉斯函数及其椭率对于考察化学键的特征很有意义，参考《Multiwfn支持的分析化学键的方法一览》（http://sobereva.com/471）和《AIM键临界点处电子密度拉普拉斯值符号判断相互作用类型失败原因的图形分析》（http://sobereva.com/161）。静电势也是基于电子密度定义的，重要性极高，我写过大量相关文章，汇总看《静电势与平均局部离子化能相关资料合集》（http://bbs.keinsci.com/thread-219-1-1.html）。信息论在化学体系中的应用已被广泛探索，其中有很多函数和量是基于电子密度定义，参看《使用Multiwfn计算各种与信息论相关的量（information-theoretic quantities）》（http://sobereva.com/537）。

4 总结

上文对电子密度相关分析做了一个较为全面的说明。从本文可见，电子密度相关的分析超级丰富，所以提问时切勿简简单单问别人“怎么分析/计算/考察电子密度？”，必须说清楚你想怎么分析讨论。如果你对方法学不懂，就要告诉别人你的分析目的是什么，这样别人才能告诉你具体怎么做，而且如果对于你的目的不适合光靠分析电子密度来说事，别人也能及早指出。从上文也看到，电子密度相关分析主要都是Multiwfn等波函数分析程序的事，所以提问时也别问诸如“怎么用Gaussian计算电子密度？”这种问题。

限于本文篇幅，上面的各种分析方法只能是三言两语简单提及，所有这些方法在“量子化学波函数分析与Multiwfn程序培训班”（http://www.keinsci.com/workshop/WFN_content.html）里有最全面系统的讲解并给了大量例子，可以一次性完整学透，欢迎参加。