VSEPR模型^[1]是本科无机化学教学的重点之一，目的在于教会学生原子之间经典的键合规则，从而了解化合物的结构。张祖德所著的《综合化学：无机化学∙分析化学∙有机化学》中有这样一道习题^[2]：

在方程式：6P₂H₄ + 10Cs → Cs₂P₄ + 8CsPH₂ + 4H₂中，生成的P₄²⁻离子团的几何构型为()

(A) 平面四边形  (B) 正四面体

(C) 变形四面体  (D) 直线型

学生可以运用VSEPR模型解决这一道题目，在溶液中，P₄²⁻离子团能被阳离子稳定，它的构型符合传统价键理论的推断，即平面正方形。事实上，如引入原子团簇这一概念，可深入探究得：在真空中，P₄²⁻团簇的构型是多样的，且具有特殊的物理化学性质。

原子团簇具备确定的结构与组成，又处于纳米尺度，可能兼具纳米材料与分子材料的特性，一直以来备受关注^[3]。原子团簇易制备出纯净产物，已知的纳米团簇化合物，如金纳米团簇，在吸附、催化、导电率等方面有着不俗的表现^[4]。

Alexandrova等曾经在《Advances in Chemical Physics》一书中“Sub-nano Clusters Is the Last Frontier of Inorganic Chemistry”这一节中提到：纳米团簇通常具有不符合我们直觉的不寻常的结构。人们发现纳米团簇的键合规则比传统分子中的更复杂，而且往往是前所未有的。而纳米团簇在材料中的应用开发最终需要对团簇中的化学键有深刻的理解^[5]。

目前，随着计算机技术的飞速发展和普及，计算化学逐步被非物理化学专业研究生以及高年级本科生所接受并使用^[6–8]。计算化学是研究物质结构与性质的有力工具，不仅能够帮助学生更好地学习物质结构与性质，又可以拓展前沿知识，与将来的科学研究接轨。

P₄²⁻团簇首先由Korber等^[9]报道，他们基于碳和磷在元素周期表中的对角线关系，推测出芳香烃C₄H₄²⁻的同系物P₄²⁻也具有芳香性，随后基于磷盐Cs₂P₄中的平面正方形的P₄²⁻离子单元对其6π芳香体系进一步研究。液氨中Cs₂P₄的合成方法如下：

以二磷烷作为新的磷源，在液氨中成功地合成了Cs₂P₄。Korber等将P₂H₄在−78 ℃下用铯处理，在停止放出气体之后，将干燥的溶液冷凝到反应容器中。反应混合物在−40 ℃下静置几天后，其主要产物Cs₂P₄以透明的黄色立方体状晶体析出。反应方程式为：

$6{{\rm{P}}_2}{{\rm{H}}_4} + 10{\rm{Cs}} \to {\rm{C}}{{\rm{s}}_2}{{\rm{P}}_4} + 8{\rm{CsP}}{{\rm{H}}_2} + 4{{\rm{H}}_2}$

X射线晶体结构分析显示，近似正方形的平面P₄²⁻单元存在于这种金属有机化合物之中。随后，Korber等在四氢呋喃中用铯处理白磷，使之溶解在液氨中，反应得到P₄²⁻以副产物形式存在于上清液中，从而得到溶液环境中的P₄²⁻离子。

Korber等又利用不同的理论水平(HF/6-311++G**，HF/6-311++G(3df, 3pd)，HF/aug-cc-pVTZ)计算出P₄²⁻的π分子轨道，与C₆H₆、C₅H₅⁻、C₄H₄²⁻、P₅⁻、S₄²⁺这些已知芳香族物种的π分子轨道相同。在HF/6-311++G**级别下用AIM方法计算了P₄²⁻的键级，与Forst-Musulin模型吻合。在HF/AUG-CC-pVQZ级别下用GIAO方法计算得到P₄²⁻的化学位移为359.6，与实验值的偏差仅为11。这一系列计算证实了属于6π芳香体系P₄²⁻的存在^[9]。

2005年，Korber等通过ELF方法分析P₄²⁻与其同系物的化学键，结合盆布居数方差研究电子离域和芳香性，对P₄²⁻的芳香性的研究表明：较大的布居数和电子波动方差体现在P₄²⁻的孤对电子上，而不像典型的芳香烃一样体现在化学键的盆布居中。实际上，P₄²⁻中的孤对电子正是其芳香性特征的体现，Korber等因此提出了“孤对电子芳香性”的概念^[10]。

然而，Pn₄²⁻ (Pn = P，As，Sb)双电荷阴离子不是电子稳定的物种，在既往的研究中，研究者们多选用碱金属阳离子来稳定Pn₄²⁻离子，以保证不会显著改变它们的电子结构和几何结构，但其研究也或多或少受到阳离子的影响。同时，HF方法对于Pn₄²⁻离子计算精度明显不够，因为HF方法本身缺少对电子相关作用的考虑。本文期望通过更精确的计算化学方法对前人的研究进行验证，进一步详细研究P₄²⁻的结构特征与理化性质。

原子团簇可能具有多种稳定构型，而传统的VESPR理论并不能很好地比较一些构型的稳定性。本文利用molclus软件对团簇进行结构搜索，从1000个初始团簇构型中筛选出能够稳定存在的结构。在PBE0-D3(BJ)/def2-SVP水平下^{[11, 12]}对其结构进行优化和振动分析，同时考察部分团簇结构之间相互转化的过渡态，并通过内禀反应坐标(IRC)验证了过渡态结构的准确性。此外，在高精度的CCSD(T)/jun-cc-pVQZ水平^[13–15]下进行了单点计算，并获取各结构的Gibbs自由能，给出了各团簇结构的相对能量和相互转换的势能面。团簇结构的搜索流程如图 1所示。

首先利用molclus程序^[16]插件genmer来生成最初的P₄²⁻团簇结构：修改参数文件genmer.ini中的参数为“ngen = 1000; iradtype = 0; radscale = 0.8; 4 rmax 3.0”，即规定期望生成的团簇数目为1000个，P原子之间接触采用共价半径，设定的P―P之间的距离小于3 Å (1 Å = 0.1 nm)。此时程序给出1000个可能的初猜结构。

接着，在MOPAC ^[17]中对这些结构使用半经验方法PM6-DH+ ^[18]进行初步优化。运用molclus自带的isostat工具，通过比对结构的绝对能量大小剔除掉优化后相同的结构，此时剩下4个结构，P₄²⁻结构的初筛已完成。

最后，在Gaussian软件^[19]中，在PBE0-D3(BJ)/def2-SVP水平^[9]下优化这四个结构，并将收敛标准设定为tight，通过振动分析确保优化收敛到极小点。PBE0杂化泛函适合对主族元素形成的原子团簇进行优化计算^[20]，相对较小的def2-SVP基组被证明兼顾较低的计算成本和较高的准确性^[21]。

在Gaussian软件中，利用QST2方法，两两寻找这四个团簇转化的过渡态，并对每个过渡态进行内禀反应坐标(IRC)计算，每个过渡态IRC计算结果显示能量曲线最低点分别指向其对应的反应物和产物，说明计算所得的过渡态是反应通道上的真实过渡态。研究结果得到3个过渡态。

将优化后的4个团簇结构与3个过渡态结构在CCSD(T)/jun-cc-pVQZ级别下进行单点计算获得单点能，并将振动分析获得的Gibbs自由能校正量与之相加，即得这些物种的Gibbs自由能。

计算所得的P₄²⁻团簇构型及相对能量如图 2所示。共有4种较稳定存在的P₄²⁻团簇，按相对能量从小到大分别编号为结构1、2、3和4。结构1 (D_4h)呈平面四边形，与Korber等实验中^[9]存在于Cs₂P₄晶体的“近似正方形的平面四边形”相吻合。4个P原子位于正方形的4个顶点，P－P键键长2.17 Å。结构2 (C_2v)由两个非平面的共边等边三元环构成，P－P键键长为2.26 Å，两根非桥头P与桥头P所成键的夹角为96.1°。结构3 (C_2h)是一个平面结构，键长为2.03 Å，P－P键夹角为109.6°。结构4 (C_2v)中，一个等腰三元环的顶端P水平地连接着一个P原子，4个P原子位于同一平面内。它的3种P－P键长分别为2.00、2.10、2.47 Å，其中键长为2.10、2.00 Å的两根键形成的夹角为143.9°。我们得到了4种较稳定存在的P₄²⁻团簇。观察结构分析，这些团簇之间结构上有一定相似之处。例如结构3仿佛是结构2中两根P－P键被“拉长”所得。分子构型的转变中常常会有过渡态存在，过渡态理论也是化学动力学研究的热点之一。在PBE0-D3(BJ)/def2-SVP级别下进行过渡态搜索共找到3个过渡态，结构1与3之间、结构2与3之间、结构4与3之间的过渡态分别用“TS1-3”“TS3-2”“TS3-4”表示。

过渡态搜索结果见图 3，从图中可以看出，结构1是热力学产物，而结构2是动力学产物。结构1并不能直接转化为结构2和结构4，而需要经过结构3方可转化。首先，结构1中的一根P－P键逐渐拉长，其Mayer键级^[22]从1.227降低至0.418。从结构1到结构3转化的能垒约为160 kJ∙mol⁻¹，根据过渡态理论，这一过程在0 K下的半衰期约为1.20 $ \text{×}\text{} $10¹⁵ s ^[23]，所以是极难进行的。但从结构3向结构1转化则是非常容易，半衰期仅有0.0332 s。因此，即便有结构3生成，也会非常快地向结构1转化。而结构3向结构4或结构2的转化则十分容易。讨论反应是否容易发生应当取决速步的自由能垒，以向结构2转化为例，反应能垒仅26 kJ∙mol⁻¹。用单分子过渡态理论计算速率，反应速率达1.7 $ \text{×}\text{} $10⁸ s⁻¹。结构3向结构4转化的反应能垒为63 kJ∙mol⁻¹，反应速率达59.384 s⁻¹，结构3向结构1转化的反应速率达20.805 s⁻¹这两个转化途径的反应速率远远小于结构3向结构2转化的速率。结构3转化得到的稳定结构的结构2的相对能量也低于结构3、结构4的值，过渡态TS3-2属于C2点群，其两根被“压缩”的P－P的Mayer键级为1.480。

Hammond假说^[24]指出，放热反应中，过渡态结构与反应物类似；吸热反应中，过渡态结构与产物类似。由图 2所示的相对能量和构型不难看出，P₄²⁻团簇的结构转化这一过程符合哈蒙德假说：TS1-3的结构不类似于呈环状的结构1，其两条等长的P－P键键长为2.05 Å，已经非常接近于结构3中2.03 Å的键长，与结构3近似；TS3-4与结构4有着相似的键长和键角，有结构4中三元环的相似结构；TS3-2中新键还未生成，还不具备结构2中三配位P的结构，与结构3更接近。

电子定域是常用于分析电子结构的工具，ELF函数^[25]的投影图能够清晰地呈现电子的ELF函数数值，这一值处于0到1之间，值越高则代表着该区域的定域性越强，电子不能随意地在此区域离域，相反，在值较低的区域，电子定域性弱，很容易离­­域。

图 4中，结构1的投影图显示，4个P原子构成的平面正方形的结构中间的ELF值很低，证明在这个区域内电子的离域性很强。结构2、3、4内部则没有这样的区域。结构2作为一个类似蝶形的结构，选取其中一个三元P环进行ELF分析得到如图 4-2所示的结果。结构3的4个磷原子整体呈现线状连接，1号P与4号P较远，即使它与结构1、结构4同为平面结构，其离域能力也大大降低。结构4独有的三元P环结构使环中间的ELF值降低至0.5–0.6，证明其有一定的离域能力，但离域仅仅发生在三元P环之中，与结构1的四元环离域相比，离域程度要小。

体系的芳香性是因为π电子多中心离域而产生的，在计算ELF的时候若只考虑π电子，这对于考察芳香性将十分有价值。于是我们使用ELF-π^[26]来更直观地揭示电子离域路径。ELF-π的二分点，在数学上被定义为该函数的(3, −1)型临界点，它对应于ELF-π的等值面刚开始分裂成两个的位置，其值已经被证明可以用来量化共轭的强度。人们认为，分叉点的ELF-π值越高，电子越容易在该点连接的两个区域之间离域。我们对P₄²⁻团簇稳定结构中的三个平面结构进行ELF-π值计算：平面正方形的结构1的ELF-π值为0.63，结构3的ELF-π值为0.22，结构4的ELF-π值为0.58。这些数值与投影图所示的离域程度的大小是吻合的。

我们运用Multiwfn软件^[27]对P₄²⁻团簇进行ELF盆分析旨在得到更精确的盆布居数并进一步得出盆内的定域化指数(localization index，LI)和盆间的离域化指数(delocalization index，DI)。LI值衡量平均来说有多少电子定域在盆里；DI值衡量的是平均有多少电子离域在两个盆之间，或者说盆之间共享了多少电子。结合LI和DI的数值，我们能够对结构1不同部位的离域情况有更清晰的把握，表 1列举了P₄²⁻团簇结构1的盆布居数。盆C(P)中的电子是磷的内核电子，V(P)指的是P原子的非成键区域的电子，这里应当指P上的孤对电子，V(P, P)代表这个盆的电子来自于两个P的价层电子，即相互作用的电子。C(P)的积分值为9.845，接近10.000，与P原子内核电子为10是相符的。在传统价键理论中，P－P键属于单键，两个P原子共用一对电子，那么积分值应当为2.000，这里实际的积分值为1.613，与2还是有一定差距的，说明纳米团簇的键合规则与传统价键理论相比是有区别的。孤对电子的布居数积分值高达3.801，这也是和经典化学理论中的解释不同。事实上，团簇的键合方式不能仅仅依靠经典化学理论来解释，对于团簇的计算与分析是基于实空间的，经典化学中的许多概念描述仅便于讨论，实际上并不准确。

图 5展示了结构1的ELF盆分析，已经被标注的区域是吸引子聚集的位置，也就是盆的位置。橙色球体代表着磷原子，相应地9、10、11、12代表的是C(P)，2、4、5、8代表的是V(P)，1、3、6、7代表的是V(P, P)，使用数字编号可以使接下来的DI值和LI值的分析更加直观。

DI(9, 10)代表相邻P的C(P)之间的离域化指数，由表 2离域化指数矩阵可知仅为0.004，而对角P的DI(9, 11)值则更低，仅为0.001，这些数值非常接近于0，这表明在4个磷原子的内核区域之间电子很难相互离域。DI(9, 1)的值为0.136，这代表着C(P)和V(P, P)之间的电子有一定的离域。DI(2, 9)和DI(2, 1)的值分别为0.521和0.616，这一值较大，证明体系中的离域多在V(P)与C(P)，V(P)与V(P, P)之间，且离域较为容易。表 3中，从结构1的定域化指数来看，P的内核的盆的LI为9.192，与盆C(P)的电子布居数9.845相比差值为0.653，说明P原子大部分内核电子都待在内核区域，只有小部分电子会离域。而V(P, P)的LI为0.638，和它的电子布居数1.613有0.975的差。

V(P)的LI值为2.333，而它的电子布居数为3.801，差值为1.468，这一差值是三组比较中最大的，说明盆V(P)中的电子容易与外界电子发生交换。一般来说价层ELF盆的差值普遍较大，但盆V(P)的差值比经典理论中认为形成P－P键的盆V(P, P)的值还要大，再加之DI值也证明了体系中的离域多体现在V(P)与C(P)，V(P)与V(P, P)之间，从而证明V(P)是体系离域的关键所在。这也证实Korber等提出的“孤对芳香性”中的观点，P₄²⁻中的孤对电子是其芳香性特征体现，这一点传统芳香族化合物并不具备。

本文从量子化学的角度对一道本科无机化学习题进行深入研究，对P₄²⁻团簇结构进行分析，介绍了团簇在真空中结构的多样性，而非从传统视角判断构型问题。回顾这一道题目，在溶液中，P₄²⁻离子团被阳离子稳定，它的构型符合传统价键理论的推断，即平面正方形；但在真空中，P₄²⁻团簇的构型是多样的，且具有特殊的物理化学性质。

本文基于Korber等对于P₄²⁻芳香性的研究与P₄²⁻团簇的计算，进一步利用molclus软件，进行P₄²⁻团簇结构搜索，删繁去重，筛选出稳定的团簇结构。在PBE0-D3(BJ)/def2-SVP水平下对搜索到的结构进行几何优化、过渡态搜索与IRC计算，呈现出P₄²⁻团簇相互转化的历程。在高精度的CCSD(T)/jun-cc-pVQZ水平下对各结构进行单点计算，详细地从量子化学的角度介绍P₄²⁻团簇的物理化学性质。运用Multiwfn软件对P₄²⁻团簇的几种稳定结构做ELF函数分析，并从多个角度探究体系的离域情况，运用ELF盆分析阐述“孤对芳香性”的特点。

本科的无机化学教学更多基于经典化学理论，本文实施的内容可以视为教学中的拓展与补充。团簇中的许多概念都是违背化学直觉(chemical intuition)的，量子化学是基于实空间的分析，其结果与经典化学概念并不契合。为学生引入原子团簇的概念，可以引导学生从全新的角度看待物质的微观结构，激发学生自主探索的兴趣。本文所使用的多个计算化学软件功能强大，组合使用可以完成几何优化、振动分析、势能面扫描、搜索过渡态等多种任务，进行波函数分析各类性质，且几乎不需要计算机专业知识，易于实践。同时，本文所介绍的团簇研究方法和流程可普遍适用于各类团簇研究，学生亦可自行运用软件探索，了解感兴趣的物质真实的微观结构。