研究自旋禁阻反应的新方法——多态自旋混合（MSSM）模型

大多数化学反应发生在基态势能面上（如图1(1)所示；注意：本文不关心光化学反应），但是，在有些反应特别是过渡金属参与的反应中，会涉及到多个自旋态。Shaik等人最早提出了“两态反应性”（two-state reactivity）的概念，这自然而然地又引申出“多态反应性”（multiple-state reactivity），见综述[1]。另一方面，根据多态反应过程中的反应物、过渡态、中间体、产物的自旋变化，又出现了自旋禁阻反应（见图1(2)）和自旋加速反应（图1(3)）的说法[2]。在文献中，讲自旋禁阻反应的比较多，其次是讲两态反应，其实只是关注问题的角度不同，都属于多态反应的特例。在本文中一律讲多态反应，但是为了照顾一些人的习惯，标题仍保留了自旋禁阻反应的说法。

图1

多态反应能够发生，其本质是相对论效应中的自旋轨道耦合 （SOC）在起作用。严格考虑SOC就要用到各种二分量或四分量相对论的计算方法（如ADF的spin-orbit ZORA DFT、NWChem的SODFT、DIRAC的各种二分量和四分量从头算方法）。二分量和四分量方法在理论上很完美，但在技术上存在一些困难：

• 二分量和四分量方法的计算量比标量方法大得多，难以用于稍微大一些的分子体系
• 在公开的量子化学程序中，二分量和四分量方法没有解析Hessian（补充：仅CFOUR程序的SO-CCSD有解析Hessian，但是目前只能算闭壳层小分子），虽然可以通过数值Hessian去优化过渡态、IRC以及计算热化学量，但仅适用于小分子
• 开壳层体系难收敛，或者收敛到激发态上去
• 开壳层体系的初猜占据难以控制，很容易布居到虚轨道旋量（二分量和四分量已经没有轨道概念，准确的叫法是虚旋量），导致收敛更慢
• 开壳层体系的波函数存在数值不稳定问题，结果可重复性差

鉴于此，二分量和四分量方法在化学反应的实际计算中用的很少。

目前，研究多态反应的首选方法还是最低能量交点（MECP）方法。即，首先用标量方法寻找各种可能的自旋态势能面上的驻点。第二步，根据各种自旋态驻点能量的高低，对那些存在势能面交叉的自旋态优化MECP。由于Harvey的MECP优化程序以及卢天修改的sobMECP程序的加持，优化MECP的难度大大降低。最后一步是在MECP结构上，用二分量或四分量方法，或者更便宜的“标量+SOC微扰”方法（如TDDFT+SOC、CASPT2+SOC）进行单点能校正。由于MECP在大多数情况下不会导致新的过渡态结构（即图2(b)的Landau-Zener情形），最后一步往往就省了。MECP方法虽然已被广泛使用，但是也存在以下问题：

• 如果MECP 导致新的过渡态（即图2(a)的情形），由于它不是驻点，因此无法计算该点的振动频率和热化学量
• Harvey与合作者的研究发现，在某些多态反应中，MECP 导致的“过渡态”与真实过渡态结构之间存在较大差异，导致过渡态能量被低估[3]
• MECP只能考虑两个自旋态的交点，而对于5d元素参与的多态反应，经常有两个以上自旋态发生强烈耦合（甚至不需要交叉）

图2

2018年，Yang、Gagliardi、Truhlar提出了两态自旋混合（TSSM）模型[4]。在该模型中，SO矩阵被近似为一个2阶的有效SO模型哈密顿，其中两个自旋态之间的SOC效应用一个不依赖结构的经验常数表示。SO模型哈密顿可以展开成一元二次方程，于是混合自旋基态的能量就是一元二次方程的第一个根。相应地，混合自旋基态的梯度、Hessian公式也可以很容易地从能量表达式导出。随后，日本埼玉大学的高柳敏幸（T. Takayanagi）团队用TSSM模型结合DFT，系统地研究了许多涉及3d和4d元素的多态反应（见文献[5]以及作者其它引用该文献的论文）。他们还发现，3d和4d元素体系的结果对TSSM模型的SO经验常数的选取并不敏感，因此可以在50-400 cm−1的范围随意取一个值，免去了优化参数的麻烦。

作为MECP的替代，TSSM不存在MECP的前两个问题，但是第三个问题仍然没有解决，即不能用于多个自旋态发生强耦合的情况。为此，我们把TSSM模型推广到任意多个自旋态的情况，提出了MSSM模型，并用Fortran 90编写了相应的程序。程序有两种运行方式：

1. 内置在量子化学程序BDF中，通过BDFOPT 模块的关键词Multistate，以及选项nsoc（n=2,3,...,9）进行调用。（不过，因为上半年忙于上课和修改学生毕业论文，一直没来得及测试和更新，在已经上线的BDF中只能做结构优化。年底前会提交具有完整功能的程序。）

2. 以独立程序MultiState@GWEV（简称MS@GWEV）的形式发布，包括源代码、运行脚本、手册、算例，在Gaussian 16 中通过External 关键词调用。有人问GWEV是什么缩写。答案：GWEV=Gaussian’s pancake Wraps EVerything。

需要强调的是，MSSM并不是TSSM的简单重复，因为按照TSSM的做法，需要先写出一元N次方程的能量表达式，再进行求导。了解一点代数的朋友都知道N>2时是根本不可能的，因为一元三次方程、一元四次方程通解的表达式极其复杂，而四次以上无解析解。其实MSSM的能量和求导公式的推导并不难，参见手册，如果以前看过Hartree-Fock的求导公式，马上就会想到了。

MS@GWEV的使用手册和下载链接：

https://github.com/zorkzou/MultiState

应用案例

Takayanagi等人已经用TSSM模型研究了大量3d、4d元素的多态反应，MSSM的结果应该也差不多。Yang、Gagliardi、Truhlar用TSSM模型研究了成键饱和的5d金属钨某种络合物的脱氢反应[4]，我们的MSSM计算可以得到完全相同的结果，说明TSSM与MSSM-2是等价的。通过与二分量相对论方法进行对比，TSSM/MSSM的反应能量误差约3 kcal/mol，Harvey等人此前对3d元素的多态反应做过MECP和二分量方法的对比测试，误差也在这个范围，说明误差是可以接受的。我们最关心的是MSSM用于不饱和5d金属参与的多态反应的情况，于是对两个反应进行了标量DFT，MSSM DFT和二分量DFT计算：

• OsO+ + CH4 → OOs(CH2)+ + H2
• W + CH4 → WCH2 + H2

其中前者包含三种自旋态，后者包含四种自旋态。详细结果及讨论见论文：J. Chem. Phys. 158, 224110 (2023) (https://doi.org/10.1063/5.0151630)

关于MSSM适用性的总结：

1. 对于Xe 之前元素体系的多态反应，结果对SOC 参数不敏感，建议取值50-400 cm-1。
2. 对于5d 元素体系的多态反应，建议取值500-2500 cm-1。但是其中成键不饱和的5d元素体系，MSSM的反应能量不可靠，还需要用二分量、四分量相对论方法在MSSM 优化的结构上进行单点能校正，不过MSSM 计算的其它结果（结构、振动频率、热化学量）和二分量计算结果差不多，还是可以用的。
3. 在进行MSSM计算前，需要先用常规的标量方法优化各种自旋态势能面上的驻点，从中可以判断有哪些自旋态参与了反应，以及势能面交叉是否导致了新的过渡态。

参考文献：

[1] S. Shaik, Isr. J. Chem. 60, 938 (2020).

[2] J. N. Harvey, WIREs Comput. Mol. Sci.4, 1 (2014).

[3] D. Ricciarelli, L. Belpassi, J. N.Harvey, and P. Belanzoni, Chem. Eur. J. 26, 3080 (2020).

[4] Yang, Gagliardi, Truhlar, Transitionstates of spin-forbidden reactions, PCCP, 20, 4129 (2018).

[5] T. Takayanagi and T. Nakatomi, J.Comput. Chem. 39, 1319 (2018).