借助Multiwfn实现物理化学性质描述符的可视化：以单环B5N5为例

程学礼，赵燕云，李峰，杨永娟，韩建梅

1 泰山学院化学化工学院，山东 泰安 271000

2 泰山学院物理与电子工程学院，山东 泰安 271000

分子及其聚集体的物理化学性质是借助一系列描述符来实现的，精确定义并合理使用的描述符可以准确表达分子的结构、性质、反应特征及其内在联系。各种形式的分子轨道是为人熟知的描述符，人们热衷于把波函数以“轨道”这种图形化方式展示出来，实验化学家也习惯用HOMO和LUMO的概念解释实验事实[1,2]。随着计算化学的迅速崛起和普及，各种用于量子化学计算、动力学模拟和分子结构可视化的软件被广泛开发并不断更新换代，主流的量子化学计算软件和动力学软件本身带有强大的图形界面，但由于高昂的版权而限制它们在大学课堂的普及。近年来，一些免费软件受到广大科研工作者的青睐，特别是一些开源免费的图形软件如Multiwfn[3]、VMD[4]、CYLview[5]、Molekel[6]等。这些专业软件开源、免费、功能强大且用途广泛，方便绘制精美的科研论文插图，极大促进了理论与计算化学的发展。例如，Multiwfn是由北京科音自然科学研究中心卢天博士开发的功能强大的免费波函数分析软件，从2012年到目前的10年时间里，使用该软件发表的SCI论文已超过11000篇，用户更是遍布世界各地。结合免费软件VMD，Multiwfn可以实现几乎所有常见物理化学概念和描述符的可视化，模拟各种分子轨道[7]，以及包括扫描隧道显微镜STM图在内的实验数据模拟。

2018年9月10日召开的全国教育大会强调了学校的美育育人功能，要求提高学生的审美和人文素养，构建德智体美劳全面培养的教育体系。2019年6月《中共中央国务院关于深化教育教学改革 全面提高义务教育质量的意见》明确要求“坚持‘五育’并举，全面发展素质教育”。对化学相关专业学生展现化学中的美，在培养学生科学素养的同时进行美育教育，也应当成为大学化学课堂教学的重要内容。本文以B5N5单环为例，可视化展示多种物理化学描述符，同时展示化学的美。

随着C18碳环的成功合成，其独特的结构与性质引起广泛关注[8]。B5N5单环是氮化硼二维材料的潜在构筑单元[9]，是C10单环的等电子体，符合4n+ 2的芳香性规则。虽然在实验上还未证实，但预测B5N5将比C18和B9N9更容易维持平面环状结构，并且理论研究已经可以准确预测其物理化学性质，其氮和硼交替的平面结构很适合作为物理化学性质描述符可视化的理想模型。在结构优化的基础上，本文可视化了B5N5的分子结构及轨道、电子定域函数、分子表面静电势、平均局部离子化能和局部电子亲合能，用LOL-π分析展示了面上和面外π电子分布，模拟了STM模拟图，并用IRI和NCI分析了其成键情况。这些有益的讨论有助于学生理解结构化学课程中的Hückel分子轨道理论和共轭π键，以及有机化学中的亲电/亲核作用位点和弱相互作用方式，辅助教师课堂教学。借助免费软件实现物理化学性质描述符可视化是把科研成果成功融入课堂教学培养本科生创新能力的极好范例。

1 计算方法

本文用Gaussian 09程序包[10]在M06-2X-D3/def2-TZVP水平上完成B5N5单环的结构优化，同时进行频率验证以保证优化构型为无虚频的局部极小点。M06-2X在研究主族元素体系和弱相互作用方面表现优异，使用D3校正可确保准确描述环内B原子间的多中心键[11]。如无特别说明，本文主要以Multiwfn3.8 (dev)实现物理化学性质描述符的定量化并可视化。

2 结果与讨论

2.1 B5N5单环的结构

实验和理论证实C18的稳定构型为单-三键交替(polyyinic)结构[8,12]。近期我们[13]也报道了B9N9的结构并与C18比较，文中单-三键交替的C18结构及Cartesian坐标示于补充材料表S1。与B9N9类似，B5N5单环只存在所有B＝N键等长(0.1326 nm)的构型，属于D5h点群。优化后的坐标示于表S2。理论计算[12]证实M06-2X方法和带有长程校正的ωB97XD泛函能很好地区分C18环的单-三键交替的结构。因此，采用SCF = NoVarAcc关键词请求完全积分准确度，本文在ωB97XD/def2-TZVP水平上重新优化了B5N5的结构，进一步证实B5N5仅存在所有键等长的结构。图1显示，两种泛函得到的结构几乎完全相同。

图1 M06-2X-D3/def2-TZVP和ωB97XD/def2-TZVP水平上优化的B5N5结构参数

2.2 分子轨道

1) 分子轨道和能级。

分子轨道的形状和能级是讨论分子成键情况的基础方法，能够直观地揭示原子间的连接和键合情况。作为1,3-丁二烯和苯分子构型和离域π键的补充，B5N5的分子轨道和能级能够辅助学生深刻理解结构化学中的Hückel分子轨道方法和共轭分子结构。对B5N5分子，图2显示，第21-30号轨道为π轨道，其中第21号(HOMO-5)和第22号(HOMO-4)轨道遍及整个十元环，表明形成了离域大π键；第21和22号轨道为非简并，第23-30号轨道均为双重简并；能量更低的成键轨道主要是σ轨道，如第20号轨道反映了分子中形成了10个B—Nσ键。能量更高的轨道为未填充电子的空轨道，图2表明第31-33号轨道为π*反键轨道，其中第31号(LUMO)轨道均匀分散在整个分子体系。计算结果表明HOMO和LUMO间的能隙为9.435 eV。

图2 B5N5的部分分子轨道能级图(空轨道以紫色显示)及用GaussView 5.0.8绘制的分子轨道图。红线为费米能级，绿线为偏电压VB (以V为单位)。绿线和蓝线之间的4个轨道对STM图有贡献。能量以eV为单位

对共轭分子，π轨道是人们关注的焦点，能够帮助大学生理解有机化学的共轭分子和离域π键，并由此衍生出许多重要理论。本文介绍用定域化轨道定位函数(Localized orbital locator，LOL)[14,15]和电子定域化函数(ELF)[14,16,17]展示共轭分子的成键情况。

2) LOL-π分析。

仅由π电子贡献的定域化轨道定位函数(LOL-π)是揭示π电子特征的流行方法，并能以填色图的方式清晰地展示离域电子的分布。图2表明，B5N5的π成键轨道分在环上/下方离域(第22号，简并的25和26，以及简并的29和30号轨道)和沿着环平面离域(第21、23、24、27和28号轨道)两类。由π电子贡献的定域化轨道定位函数(LOL-π)填色图示于图3，面外和面上的LOL-π填色图均展示了B＝N双键及π电子离域状态。并且，N原子上的π电子明显多于B原子，说明N上的孤电子对离域到整个分子形成第2个离域大π键。图3由Multiwfn3.8 (dev)绘制，读者可根据个人喜好灵活设置。

图3 (a) B5N5平面上方1.2 Bohr处以彩虹色填色的面外π成键轨道LOL-π平面图；(b) 以黑-橙-黄渐进色填色的面上π成键轨道LOL-π平面图

3) ELF分析。

电子定域函数ELF是衡量电子定域性和离域性的三维函数，是研究化学键、离域π键和多中心键的重要方法，用Multiwfn3.8 (dev)绘制的B5N5分子平面(XY平面)填色图见图4。图4完美地展示了N上孤电子对的均匀分布情况，但缺陷是没有很好地展示π的离域性和多中心键。图片可根据个人爱好自由设置，带等高线的ELF填色图示于补充材料图S1。

图4 B5N5 XY平面的ELF表面填色图(Z = 0.6 Bohr)

4) STM模拟图。

根据Tersoff-Hamann (TF)模型[18,19]，STM图像是分子中特定区域局域态密度的反映，实际上与HOMO附近的某些成键轨道有关，目前Multiwfn支持模拟分子体系常高和常电流模式下的STM图像。分子体系的费米能级EF为EHOMO和ELUMO的平均值，对B5N5分子EF= -4.973 eV。如果要显示常高模式下第27-30号轨道(HOMO和HOMO-1)贡献的STM图像，偏压VB应设置在-7.287 - -5.499 V。当设置VB为-6.0 V时，轨道能高于-10.973 eV的占据轨道将被显示(图2中蓝线)，刚好对应第27-30号轨道。此时，常高模式下用Multiwfn模拟的STM图像如图5所示。

图5 偏电压为-6.0 V下的B5N5模拟STM图像

2.3 表面静电势及静电势极值点

原子电荷是反映分子中电荷分布的基本概念，是对位于原子中心点电荷分布的简单、直观的描述方式，即位于原子中心的点电荷[20]。1955年提出的Mulliken电荷是最古老且被普遍接受的原子电荷描述符，目前流行的还有NPA电荷、Hirshfeld电荷等形式。一般来说，中性分子中同时存在负静电势和正静电势区域，分别与亲电试剂和亲核试剂作用。范德华填色的表面静电势以更为直观的可视化方式考察分子表面电荷分布，通过不同颜色展现不同表面区域静电势的相对大小，在分析分子间静电相互作用、预测反应位点和分子性质等方面具有重要意义[21,22]。按照网络文章《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》[23]，采用Bader定义的范德华表面，单环B5N5静电势等值面示于图6。图6清楚地展示了B原子带正电，对应各自的局域静电势极大点130.4 kJ·mol-1，5个B原子共同作用，使得全局静电势极大点(147.9 kJ·mol-1)出现在环中心；N原子带负电，对应着静电势极小点-64.4 kJ·mol-1。

图6 由Multiwfn3.8 (dev)和VMD 1.9.3绘制的B5N5表面静电势(Isovalue = 0.001 e·Bohr-3)及静电势极值点(kJ·mol-1)

2.4 亲电、亲核反应位点

准确地预测分子的亲电、亲核反应位点是有机化学家梦寐以求的美好愿景。随着理论化学的发展，已有多种方法做到精确地定量描绘分子亲电和亲核位点，如计算福井函数、平均局部离子化能(ALIE)、局部电子亲合能(LEA)等[24,25]。ALIE和LEA指数可以借助Multiwfn和VMD很方便地可视化展示(Multiwfn3.8 (dev)手册4.12.2和4.12.13节)。ALIE和LEA着色的分子表面图示于图7。绘图时，默认的色度条显示范围为0.32-0.36 a.u.，分子被红色区域笼罩，表明B5N5是缺电子分子，很难给出电子；调节色度显示范围为0.42-0.46 a.u.可得图7(a)。青色圆球(N原子附近)是分子表面上ALIE的极小点，体现了电子被束缚得最弱的位置，也因是容易发生亲电反应的位点。有意思的是，LEA着色图说明，B5N5最容易结合电子的位置竟然是5个硼原子所围成五元环的中心和B—B键(相邻2个硼原子连线)的上下两侧，基本对应每个氮原子的斜上方，如图7(b)所示，这也意味着受到亲核试剂进攻时，B5N5环容易被破坏。增大色度范围至-1.00 - -0.20 a.u.，图7(c)显示，在分子平面上，B—N键的中心位置也较易受亲核试剂进攻而发生插入反应。

图7 B5N5的ALIE着色图(a)，以及色度范围为-0.80 - -0.30 a.u. (b)和-1.00 - -0.20 a.u. (c)时的LEA着色图

2.5 IRI分析和NCI分析

分子间和分子内的弱相互作用是化学家关注的焦点之一，近年来也发展出一些分析弱相互作用的新方法，如相互作用区域指示函数(IRI)[26,27]、非键作用函数(NCI)[28]、基于Hirshfeld划分的独立梯度模型(IGMH)[29]等。其中IRI函数不依赖于波函数，只要获得电子密度及电子密度梯度即可计算，且能够同时图形化展示化学键和弱相互作用，特别是环内原子间的相互作用。平面分子B5N5的IRI等值面填色图(图8)按网络文章《使用IRI方法图形化考察化学体系中的化学键和弱相互作用》3.2部分[30]绘制。图8不仅很理想地展示了B和N间的共价作用，还完美地描绘了B5N5的对称结构。有趣的是，在环内五角星空洞的中心存在明显的成键作用，它应该对应5个B原子间形成的多中心键。另一种IRI绘制方法(图S2)和NCI分析(图9)也证明了在环中心存在成键作用(绿色区域)，而不是库仑排斥。非键作用(NCI)分析通过约化密度梯度函数(Reduced density gradient，RDG)构建一个实空间函数，以凸显弱相互作用区域，因此也被称为RDG分析，是2010年发展起来的方法，与更新的IRI方法相比，仅能展示环中心的弱相互作用。

图8 B5N5的IRI填色平面图(Isovalue = 1.0 a.u.)

图9 由Multiwfn3.8 (dev)和VMD 1.9.3绘制的NCI分析等值面

3 结语

结合单环B5N5的结构特点，本文介绍了B5N5的分子结构及轨道、分子表面静电势、π电子分布和分子成键情况等物理化学性质描述符的可视化展示，而这仅仅是众多展示方法中的部分可视化案例。这些描述符的可视化不仅有益于本科生充分理解结构化学中化学键和弱相互作用的形成及特征，特别是共轭体系的离域键和芳香性概念，表面静电势、ALIE和LEA等描述符的图形化表达还可以帮助学生掌握有机化学中的分子表面特征和亲电/亲核位点。近年来发展起来的图形化方法还有很多，例如可以可视化方式描述分子内和分子间的弱相互作用，以及反应过程中的物理化学性质变化等。借助化学软件实现物理化学性质描述符的可视化，必将对课堂教学改革起到巨大促进作用：(1) 将抽象、枯燥的物理量和物理学及化学概念以直观、生动的方式展现出来，有利于提高学生的学习积极性和主动性，从而提高教学质量；(2) 把科研方法应用到课堂教学中，本身是对大学生科学素养和科研素质的培养和启迪；(3) 加强大学化学课堂的美育教育，让学生认识化学中的美，培养学生的审美能力并激发学生的创造性；(4) 把以前没法讲、讲不透的知识模块以简单直观的方式传授给学生，有助于学生形成逻辑严密且有机统一的化学基本概念和理论体系。

补充材料：可通过链接http://www.dxhx.pku.edu.cn免费下载。