利用基于直接动力学的轨线面跳跃方法处理非绝热过程

彭佳伟 ，谢宇 ，胡德平 ，杜利凯 ，兰峥岗 ，＊

1中国科学院青岛生物能源与过程研究所中科院生物基材料重点实验室，山东 青岛 266101

2华中农业大学信息学院湖北省农业生物信息学重点实验室，武汉 430070

3中国科学院大学，北京 100049

1 引言

19世纪以来，随着人们对化学问题研究的不断深入，在原子分子水平上理解化学现象变得越来越重要。由于原子核和电子在质量上存在很大的差异(一般3个数量级左右)，电子的运动比原子核快的多，所以从理论上来说，可以近似地将两者运动分开进行处理，即采用所谓的“玻恩-奥本海默”(Born-Oppenheimer)近似1，该近似能够有效地描述很多发生在电子基态上的化学过程。

有一类问题不能在BO近似的框架下去研究，即非绝热过程2–4。该类过程广泛存在于众多化学变化中，特别是当反应涉及电子激发态的时候。在该过程中，沿着原子核坐标变化的方向，不同电子态的势能面可能会出现交叉或避免交叉。当体系处于交叉点附近时，不同电子态之间存在较强的耦合，体系很容易从一个电子态跃迁到另一电子态上。这时，电子运动和核运动就耦合在一起，无法分离，因此基于 BO近似的理论方法变得不再适用2–4。

多年以来，研究者们试图发展各种理论工具，用于研究非绝热过程。一些严格的量子动力学方法，可以给出很好的描述，但往往由于计算量太大而难以处理非常大的体系。比如：标准的量子动力学方法、多组态含时Hartree(MCTDH)方法5以及多层MCTDH(ML-MCTDH)方法6–9等。一些基于随机波函数的方法10,11也被用来求解量子动力学。

对于处理复杂体系，一类常见的方法是使用量子耗散动力学，考察体系的约化密度矩阵随时间的演化，包括 Redfield方法12、级联运动方程(HEOM)方法13以及各种基于路径积分的方法14等；这些方法往往偏重于对物理模型的理解，需要使用特定的哈密顿形式。

处理非绝热动力学的另外一个思路是引入半经典近似。较严格的半经典方法15可以比较完善地处理非绝热动力学，但面临计算量太大的问题。因此一直以来，研究者们一直试图发展一些实用性强的半经典方法，比如：基于Mapping哈密顿2,15–22的准经典动力学方法、Ehrenfest方法和其变 种2,15,23–28以 及 轨 线 面 跳 跃 方 法等2,15,24,25,29–37。尽管这些方法在理论上，存在不严格性，但是由于实用性强而一直被研究者使用。在以上方法中，轨线面跳跃方法由于其易于程序化、较好的平衡了计算量和准确性，模拟的结果可以比较好的给出化学图像(如体系是处在激发态还是基态，是解离态还是束缚态)，因此受到了国内外学术界的广泛关注。特别是该方法易于和直接动力学结合37–46，能够在动力学演化中直接调用电子结构计算，所以可以用于研究真实的复杂多原子分子体系。如果将其和基于多尺度的QM/MM47–52方法结合，则可用于研究溶液、生物和材料中的非绝热动力学53–62，因此近年来获得了很快的发展。

将轨线面跳跃动力学和直接动力学结合，其数值实现涉及很多细节，包括初始抽样、动力学演化、电子结构计算的处理以及如何将动力学和量化计算进行较好的结合等诸多细节。为了比较全面的讨论这些内容，本文拟针对该问题，讨论如何将轨线面跳跃方法和直接动力学进行结合，对这一领域的发展与应用进行一定程度的梳理。

文章将从以下四个部分进行展开，基本理论部分主要对轨线面跳跃方法的基本理论进行回顾；方法实现部分主要对轨线面跳跃方法的数值实现进行介绍；方法应用部分主要对基于直接动力学的轨线面跳跃方法处理的一些典型问题进行介绍；挑战与展望部分主要就对轨线面跳跃方法下一步发展的一些方向进行展望。

需要说明的是，由于该研究课题涉及非常多的内容。已经有很多文献和综述曾经系统的讨论过其中一些内容2,63–67，包括轨线面跳跃动力学的基本理论以及发展，以及如何对其进行一些特殊近似用于处理复杂体系等。为了不和这些内容重复，本文的重点将关注于如何将面跳跃与直接动力学更好的进行结合，使其可以用于描述分子体系的非绝热动力学。同时简单讨论该类研究的一些例子。对于该理论的具体应用，我们也仅仅列举了一些分子体系的例子。由于作者知识有限，论述中难免会出现一些缺陷，内容也可能会有所遗漏，对此敬请读者原谅。

2 基本理论

2.1 重要概念

2.1.1 哈密顿量和绝热态、透热态

对于分子体系，含时的薛定谔方程可以写成

其中H包括原子核和电子两部分，称为总体哈密顿

而 He包括电子的动能(Te)、原子核之间的排斥相互作用(Vnn)、电子之间的排斥相互作用(Vee)以及电子与原子核之间的吸引相互作用(Vne)，称为电子哈密顿

整体的波函数 Ψ(R, r, t)可以用绝热态电子波函数严格展开，即以所谓的玻恩-奥本海默-黄(Born-Oppenheimer-Huang)2–4展开

其中 χi(R, t)为原子核的波函数，ϕi(r； R)为电子的波函数。R、r分别为原子核和电子的坐标矢量。

所谓的绝热态，即电子哈密顿的本征态，满足方程

其中Ei(R)为特定原子核构型下第i个绝热电子态所对应的本征能量。利用电子波函数的正交归一性，可得到原子核波函数的演化方程

其中 Fij(R)为一阶耦合矢量，Gij(R)为二阶耦合系数2–4：

对于一阶耦合矢量，很容易发现其具有反厄米性质。一阶耦合矢量可以表示为2–4

从上式可看出，两电子态间的耦合强度与其本征能量的差值成反比。当两个电子态间的能量差接近于零时，一阶耦合矢量趋向于无穷大，因此在势能面交叉点附近，很容易发生非绝热跃迁。

当势能面出现简并时，一阶耦合矢量趋向发散，因此很难将其用于量子动力学演化。为了规避这一问题，研究人员通常构建一组新的电子态，即所谓的透热态2–4。该态可以通过对绝热态的酉变换得到，即

在该表象下，电子的基函数不再是电子哈密顿的本征态，不同透热态的耦合变成了光滑的势能耦合。

尽管利用透热态可以避免一阶耦合矢量发散的问题，但是对于分子体系，透热态构建本身并不容易。严格意义上的透热态只有在双原子分子体系和电子态完备的多原子体系中才能够得到2,4,68，而通常所说的透热态都是近似的。对于多原子分子体系，曾有综述指出2，建立透热态常用的方法主要可以分为三类。第一类方法是基于透热态表象中一阶耦合项消失的特点，直接通过积分求解透热态3,4。该方法原则上是精确的方法，但是它需要计算积分路径上各个构型的非绝热耦合项，并不容易，并且结果依赖于路径的选择。第二类方法是通过使波函数本身或者特定的物理量变得平滑，但该方法依赖于体系的物理性质。第三种方法是基于势能面本身，假设透热态哈密顿矩阵的解析形式，通过将透热态的势能矩阵对角化来得到绝热态，结合从头算获得的绝热态能量，通过拟合方法，建立起绝热哈密顿。对于更为复杂的问题(如复杂体系内的激发态能量和电子转移等)，也有一些其它的透热态化方法，如库仑耦合方法69,70、Generalized Mulliken-Hush (GMH)方法71及其变体70、定域透热态化算法72和多态密度泛函方法73。

2.1.2 势能面交叉

在非绝热过程中，势能面的交叉具有不同形态。当两个电子态自旋相同时，出现的势能面相交称为锥形交叉2,65。该类交叉与很多分子的光化学性质息息相关，如DNA的光稳定性74,75、荧光蛋白等功能蛋白的光反应76,77以及有机共轭分子的光异构化77–80等。在交叉点上，由于两个电子态的严格简并，非绝热耦合矢量会发散。根据 noncrossing rule81,82，双原子分子体系不存在锥形交叉。对于多原子分子体系，根据电子态的对称性，锥形交叉可分为三类2,65。第一种是两个电子态具有不一样的对称性，两者沿着全对称的振动坐标方向，势能面发生交叉。由于两个电子态的波函数属于不同的不可约表示，即具有不一样的对称性，就会存在严格的能量简并。然后，如果特定的振动导致分子的对称性降低，两个电子态变成同样的不可约表示，具有同样的对称性后，则两者开始出现耦合，简并消失。此时，导致两个不同电子态的波函数成为同样不可约表示的振动模式与产生势能面交叉的全对称振动模式就构成了一个两维的坐标空间，在该空间内，两个势能面就变成了圆锥交叉的结构。第二种是对应于高对称性分子，如具有三重态旋转轴的分子，对于其属于简并的电子态，如果分子振动破坏了这种简并，那么就会导致出现锥形交叉。其实该现象就是化学中经常出现的姜-泰勒效应83。第三种情况，就是分子本身没有什么对称性，电子态也不存在对称性，但是由于核运动造成了电子态的偶然简并。过去该类交叉被认为很难出现，后来大家发现这种偶然的交叉其实广泛的存在于光化学中2,65。

不同自旋态的交叉，顾名思义，指的是自旋不同的电子态出现的交叉，在该类交叉附近，体系会发生自旋翻转84–88，该类情况广泛的存在于过渡金属的光反应89,90以及一些氧化还原反应中91,92。需要说明的是，在该类交叉中，电子态耦合可以使用自旋-透热态的框架描述，即电子态此时虽然是电子哈密顿的本征态，但是波函数是在限制自旋后获得的。不同自旋态的耦合是自旋-轨道耦合。对于自旋不同的电子态而言，实际上该处理方式隶属于透热态的框架。如果在相对论的框架下，求解狄拉克方程，重新考虑该问题，那么才能真正生成绝热电子态。

有时，具有同样自旋的不同势能面会发生靠近，但无法实现真正交叉，这种情况称为避免交叉，分子聚集体中激发态的电子和能量转移过程往往与这种交叉有关93,94。

2.2 面跳跃动力学的基本理论

所谓的面跳跃方法其实是一种半经典的动力学方法23,24,30,37,95，可以对很多非绝热过程进行研究。在该方法中，一般会先选择合适的方法进行取样，然后进行相应动力学的模拟。其中原子核自由度以经典的方式在单一势能面上进行演化。在整个演化过程中，非绝热跃迁通过轨线从一个电子态跳跃到另一电子态上来描述。最后，通过演化一系列轨线，进行统计平均获得最终的动力学特征。因此，一个比较核心的问题是如何引入跳跃几率来决定什么时候轨线发生跳跃。为了解决这一问题，研究人员发展了很多不同的跳跃算法，这也是不同轨线面跳跃方法的主要区别之一。

总体来说，面跳跃算法可以分为两大类，其中一类不需要计算电子的运动，而通过一些物理模型预先定义跳跃的几率，比如Landau-Zener算法和其变种37,95。另一类则是通过求解电子的运动方程，建立起电子运动与跳跃几率之间的联系，借此衡量跳跃发生的可能性，比如Tully发展的最少面跳跃算法23,24,30。基于以上的这些算法，轨线面跳跃方法产生了很多不同类型的变种67,96–98，包括Ring-Polymer Surface Hopping (RPSH)99, Phase-Corrected FSSH (PC-FSSH)100, Semiclassical Monte Carlo (SCMC)101, Augmented FSSH102, Coherence Penalty Functional (CPF)103, Min-Cost Algorithm104,Flexible Surface Hopping (FSH)105, Self-consistent FSSH (SC-FSSH)106, Markov State Surface Hopping(MSSH)107, Global Flux Surface Hopping (GFSH)108,Liouville Space FSSH109, GFSH (LS-GFSH)110,Decoherence induced Surface Hopping (DISH)111,Second-Quantized Surface Hopping (SQUASH)112,Fragment Orbital-based Surface Hopping (FOB-SH)113,Consensus Surface Hopping114等。

目前，基于Tully发展的最少面跳跃方法在从头算动力学中应用最多38–46。基于此，本文将对该方法进行详细的介绍。

2.2.1 电子态的演化

在轨线面跳跃方法中，电子具有量子特性，因此在特定核构型R(t)下，电子波函数(, ； )tφrR满足含时薛定谔方程

利用基组的正交归一性，可以得到展开系数cj(t)的演化方程

其中Hji为电子哈密顿矩阵，dij为一阶非绝热耦合矢量，v为原子核的速度矢量。如果我们使用密度矩阵语言，上面的表述可以写为

其中ρij(t) = ci(t)cj(t)*为电子运动的密度矩阵。

如果考虑电子波函数展开时所选的基组类型，可以得到不同表象下电子的演化方程。在绝热态下，由于Hji= Vjδji。因此电子的演化方程满足

在透热态下，一阶耦合矢量为零。相应的电子演化方程满足

2.2.2 原子核的演化

在轨线面跳跃方法中，原子核具有经典粒子特性，因此遵循牛顿力学运动方程。对于质量为M的原子核，它的演化方程满足

其中F为原子核在轨线当前态上所受的作用力，即

2.2.3 最少面跳跃算法

目前，最常用跳跃算法之一是耶鲁大学Tully发展的最少面跳跃算法30。假设体系有N条轨线，任意t时刻处于第i态上的轨线数为Ni(t)。当轨线足够多时，跳跃几率的平均值与轨线分布应该一致，满足

第i态发生跳跃的几率可以用轨线在第i态上变化率表示

因此跳跃几率可以展开成以下形式

由于第i态发生跳跃的几率等价于第i态跳到其它电子态几率的和，因此可以认为第i态跳到第j态上的几率满足方程

在最少面跳跃算法中，由于我们只关注轨线从当前态跳出的几率，因此最终的跳跃几率应该为非负值，即满足30

为了进一步确定当前态 i最终跳到哪个电子态 k上，需要在每一步中都产生一个[0,1]之间的随机数，如果它满足

则轨线可能发生跳跃；否则轨线就停留在当前态上。

2.2.4 速度调整

在轨线面跳跃中，当轨线跳跃发生时，由于势能会突然的变化，从而导致能量上的不守恒。为了解决这一问题，通常的做法是进行速度的调整。假定t时刻，轨线从电子态i跳到电子态j上，势能由Ei变为Ej。跳跃后原子核的速度满足以下方程

其中Mβ为原子核的质量，γij为缩放因子，vβ为跳跃前的速度，为跳跃后的速度，ω为沿着修正方向的N维向量。所以动能的改变量为

其中aij和bij满足以下关系，即

轨线发生跳跃，修正速度一般选取速度改变量较小的那个。由于从路径积分理论可以导出，在非绝热过程中，如果使用半经典方式处理原子核的运动，其所受的力23,24,30,31,115–119与非绝热耦合矢量相关，因此很多工作建议速度修正的方向采用非绝热耦合矢量方向23,24,30,31。对于非绝热耦合矢量不容易确定的情况，也可以采用一些其它的修正方向65,120,121。

3 方法实现

从数值上实现轨线面跳跃方法与直接动力学的结合，一般包括五个步骤：

(1) 根据所研究的问题，选择合适的取样方法，并获得一系列轨线的初始坐标、动量和电子态。对于每条轨线，在特定的分子结构下，通过量子化学计算，获得不同电子态下体系的能量、梯度以及电子态之间的非绝热耦合项。

(2) 对轨线进行动力学演化，其中核的演化遵从牛顿定理，电子的演化使用薛定谔方程。

(3) 计算跳跃几率，并与产生的随机数进行比较，判断是否满足跳跃条件。如若满足，轨线跳跃到相应电子态，并根据能量守恒，对速度进行调整。

(4) 重复(2)–(4)，直到演化结束。

(5) 最后对大量轨线的结果进行统计分析。

由于轨线面跳跃方法的实现涉及较多内容，因此关于每一步的实现细节我们将做如下展开。

3.1 采样的实施细节

在轨线面跳跃方法中，由于涉及原子核的经典化，因此需要确定轨线的初始条件，包括初始坐标、动量以及电子态。对于取样方法，其本质是将初始量子态投影到相空间中，得到分布函数，进而利用随机的方式产生初始的坐标和动量。目前常用的方法有 Wigner取样、Action-angle取样以及分子动力学取样。

3.1.1 Wigner取样方法

在量子力学中，我们往往采用波函数来描述体系的状态。但是本文中所涉及到的直接动力学方法都是对体系做了半经典近似，即通过轨线演化的方式来描述动力学过程。所以在做动力学之前，我们需要获得轨线的初始条件，即位置(q)和动量(p)。所以我们面临的首要问题是如何通过合理的方式将量子的初态对应到经典相空间里的状态，建立起经典相空间中的分布函数，然后通过分布函数，对q和p进行初始采样。

获得量子和经典对应关系的一个方法就是做Wigner变换。对于一个量子状态，我们都可以通过Wigner分布函数获得经典相空间中的q和p。相应的Wigner分布函数128,129满足

其中ρ是密度算符。

通常情况下，在初始采样过程中，我们认为体系势能面最小值附近区域能够用二次函数来近似，通过这种方法可以获得体系的振动模式，一般为3N - 6或3N - 5个(N为原子的个数)。因此在多维空间中，每个振动模式的基态波函数就可以推导出来。其Wigner分布函数可以简化成两维高斯分布，即

其中Qi、Pi分别为无量纲(dimensionless)正则坐标和动量。

利用Wigner取样，我们就可以得到一系列轨线的初始状态。对于体系处于振动激发态时，Wigner分布函数在一些区域会出现负值130,131，所以较少用于此种情况下的取样。

利用Wigner分布函数也可以进行特定温度T下的取样128。

其中α = tanh(ωћ/2kBT)，ω是频率，kB是玻尔兹曼常数，Qi、Pi分别为无量纲正则坐标和动量。

需要说明的是，Wigner分布函数的取样并不局限于特定的量子态形式。当密度矩阵已知，原则上就可以求出分布函数的对应形式。但问题在于对于任意的波函数形式，无法简单获得经典相空间分布函数的解析表达式，因此采样的时候，多考虑使用谐振子近似。

3.1.2 Action-angle取样方法

在轨线面跳跃方法中，原子核初始的坐标和动量也能够利用经典的 Action-angle取样2,123得到。对于正则振动体系，如果认为经典和量子的能量一致，则相应的无量纲正则坐标和动量分别满足以下表达

其中αi为[0,2π]之间随机的角度，ni为谐振子的主量子数。当ni= 0时，我们就得到了振动基态时坐标和动量的分布函数。同样我们可以对于特定振动模式，做特定振动态的激发，这就需要将对应振动模式下的 ni设定为相应值。最后对于特定温度T下的取样，可以通过玻尔兹曼分布函数对 ni进行取样。

其中ωi为第i个模式的振动频率。

3.1.3 分子动力学取样方法

对于轨线面跳跃方法中的初始参数，还可以利用分子动力学取样方法获得。所谓分子动力学取样方法指的是按照研究问题的类型，先选择一个有代表性的初始结构。以该结构作为起点，选择合适的动力学方法，进行动力学演化，使体系达到平衡状态。一般地判断标准是温度收敛到既定值，或出现微小的波动。接下来以此状态为出发点，继续进行动力学演化，每隔一段时间，收集一系列结构以及对应的速度。通常情况下，对于收集的样本，需要进行进一步的筛选与过滤。该类方法看似简单，但是依然需要小心的操作，避免体系的不平衡和动力学模拟尽量经历各种状态。此外，该方法最大的缺点是无法保证精细平衡原理132。

3.1.4 采样中考虑跳跃的几率

对于基态上的初始构型，我们可以很容易地利用上述三种取样方法获得。但是对于初始电子态的选择，有以下几种方法。一种常用的方法就是根据实际问题指定某个电子态，比如初始时将所有轨线都放到特定的激发态上，开始演化。也可以设定一个能量窗口，只在窗口之内进行激发。此外，有些工作建议通过考虑跃迁几率判断是否激发。对于每一构型从基态S0受光照激发到激发态Sk的几率满足133

其中μk0为跳跃偶极距，fk0为振子强度，ωk0为两电子态间的能级差。这样只要产生一系列[0,1]之间的随机数，通过与Pk0比较，就可以决定该构型处于哪个电子态上。

3.2 常见的量化方法选择

非绝热过程往往涉及多个电子态，因此需要选择合适的量子方法来获得各个电子态，尤其是电子激发态的信息。在量化方法选择上，一方面要求所选方法的计算速度应该尽可能的快；另外一方面计算结果尽可能准确，特别是考虑到从头算动力学中每一个时间都需要进行量化计算，所以平衡计算准确性和效率是非常关键的。

就现有的量化方法而言，没有一种方法能够做到对激发态的描述面面俱到。高精度从头算方法(如 CASSCF134,135、CASPT2136、MRCI137)，尽管计算的结果更加准确，但需要庞大的计算量，一般只能应用于小分子体系。此外，这些基于活性空间的方法，正确选择活化轨道对于产生正确的计算结果至关重要。半经验方法(例如 OM2/MRCI138–141、AM1/CIS142)，尽管计算效率很高，但是结果受参数的影响比较大。其它的方法143–147(如 EOM-CC、CC2、ADC(2))虽然可以处理电子激发态的计算，但是计算量依旧不小，而且也面临各类问题。目前，广泛应用的TDDFT方法147，相对而言具有易用的特点，也确实能用于处理一些问题。但是也严重依赖于泛函的选择，特别是需要小心电荷转移态的问题。该方法基于 DFT，所以对基态和第一激发态间的锥形交叉的描述存在问题148。近年来，TDDFT在原有理论的框架下，出现了一些新的处理方法，比如使用 Spin-Flip TDDFT149,150，已经可以用于处理基态和第一激发态间的锥形交叉。其外，基于Spin-Flip的各类方法150–155、PP-RPA/TDA 的方法156、Dual-Functional TDA方法157以及多态 DFT73等则可以正确处理该类交叉。

环境效应(如溶剂化效应)对于非绝热过程非常重要38,58。常用的方法包括量子力学、分子力学结合的方法(QM/MM)、郎之万方程的方法以及连续介质模型的方法等。其中QM/MM方法将研究的体系划分为两部分52。对于反应中心的部分使用量子化学的方法进行处理，对于环境则使用分子力学的方法进行处理，然后考虑两者之间的相互作用。通过结合该方法和面跳跃方法，研究者们已经可以处理溶液和生物体系中的非绝热过程53–62。

所谓的郎之万方程方法本质上是牛顿力学的延伸。考虑体系在有耗散和粘度的介质中运动，则可以通过郎之万方程以简明的方式处理环境效应。理论上来说，只要选用了合适的参数，该方法可以得到与显示溶剂模拟相似的结果158,159。该方法也曾被用于描述复杂体系的非绝热过程，如有机光伏体系的激发态动力学104,160–163。

在处理溶剂效应时，还存在一类常用的方法，即认为体系是可极化的连续介质，代表方法如PCM方法164，或COSMO165。该方法可以大幅降低处理环境的计算量。但是该方法中溶剂体系不是利用具体分子处理的，因此无法给出任何来自于溶剂分子构型相关的效应。

3.3 非绝热耦合项的计算

在绝热态下，演化电子的运动方程需要用到非绝热耦合项。对于该项的计算有解析法和数值法两种方式。原则上通过波函数的求导，可以获得非绝热耦合矢量的解析表达4,166–174。也有一些特定的量化方法在特定的量化程序实现了利用解析公式求解非绝热矢量的功能。对于这种情况，可以通过在从头算动力学中直接调用量化程序获得该非绝热矢量。因此对于有该功能的量化方法和程序，编写接口实现从头算层面上的面跳跃动力学并不困难。

但是对于一些量化方法，还无法通过简单调用量化程序获得非绝热耦合。此时，可以通过近似的数值差分方法来求解。

真正在动力学演化过程使用的非绝热耦合矩阵为利用数值差分，可以将上式转换为41,121,175–177

其中Δt为步长。

原则上，上式适用于任意波函数的理论框架。例如，在组态相互作用的框架下，电子态波函数可以写成不同组态的线性组合，然后原则上就可以求解了。

对于今天的激发态计算，TDDFT作为常用的方法之一147,178，被广泛使用。在此基础上求解非绝热耦合矢量也获得了很多关注166–170,179。但是至今，并不是所有的软件都可以实现该计算，所以下面以TDDFT为例子解释如何利用上式求解。对于基态的电子波函数，可以用占据的Kohn-Sham轨道构建的Slater行列式表示。利用Casida关系，激发态波函数能够近似的展开成CIS形式121,175,178,180–182，即

其中 φiCaSF(r；R( t ))表示从闭壳层的组态出发进行i→a单激发产生的波函数，考虑到电子自旋后，可以获得自旋匹配的组态波函数为

需要指出的是，在从头算动力学中，因为波函数相位的演化都直接存在于含时的展开系数中，见

2.2.1 章节。所以此处的电子波函数均是直接从量化计算获得，因此所有的波函数都是实数，即CI系数 ci,aK，分子轨道系数 Bjm( t)和所有积分最后均为实数。

在线性响应的TDDFT中，往往求解Casida方程获得激发能147,178

其中ω表示激发能，X和Y矢量表示单电子在不同轨道间激发和退激发的系数。对应的A矩阵包含单电子激发过程中占据与非占据轨道的能量差，以及单电子激发之间的耦合，B矩阵表征激发和退激发的耦合。对此方程的具体推导见文献147,178。

对于ci,aK可以通过TDDFT计算给出的X和Y矢量近似获得。需要指出的是，如果使用TDDFT/TDA近似，上式变为

因此，计算仅仅获得X矢量，只有激发没有退激发，这个时候，组态展开系数仅仅与X有关，CIS类型的波函数可以较为严格的构造出来。但是当使用 TDDFT/RPA(或 TDHF)时候，存在退激发项Y，这个时候，原则上无法建立其严格对应的波函数，只能近似的获得波函数。但是由于 Y往往很小，所以对于数值计算，近似波函数依然可以给出较为合理的非绝热耦合。

此外，数值方法可以用在其它理论中，如ADC(2)，此时，使用一阶响应的系数作为单激发的组态系数建立近似波函数。对于高精度的CASSCF等，波函数展开有时需要考虑双激发的贡献。

3.4 演化表象的选择

轨线面跳跃方法的应用中，一个关键是如何面对其表象依赖性。形式上来说，该方法的公式似乎可以在各种表象中都成立。常见的方法也有几种。

(1) 所有演化都在透热态下进行，即电子的演化和跳跃几率的计算都在透热态下进行；核的运动也由当前所处的透热态主导。

(2) 所有演化都在绝热态下进行，即电子的演化和跳跃几率的计算都在绝热态下进行；核的运动也由当前所处的绝热态主导。

(3) 此外还有一类方法，两种表象都涉及了。此时，核运动由当前所处的绝热态主导。电子运动先在透热态下进行演化，然后转换到绝热态下计算跳跃几率67,106,183,184。值得注意的是，混合情况下得到的结果原则上与绝热态的演化结果一致。此类方法包括 Granucci等人提出的局域透热态方法184，以及王林军等人在处理Trivial crossing问题时候发展的Self-Consistent最少面跳跃方法67,106。

尽管理论上来说这三种形式都可以选用，但就轨线面跳跃方法本身而言，第二种24和第三种方式一般能够获得比较有效地计算结果。也有工作证明似乎第三种比第二种方法更具有数值稳定性184,185。

3.5 运动方程的实现细节

在轨线面跳跃方法中，原子核的演化可以采用速度 Verlet方法和四阶及以上的龙格-库塔方法。其中速度Verlet方法，坐标和动量分别按以下方式进行更迭：

其中 M-1是原子核质量倒数向量，M-1= (,..., MNatoms-1)。

值得注意的是，在积分经典和量子运动方程时，选择合适的步长非常重要。对于原子核的运动，由于采用的是经典牛顿运动方程，原则上可以选择相对较大的时间步长，通常为0.1–0.5 fs左右。对于电子的运动，由于在求解的过程中有指数的存在，计算结果很容易出现震荡。为了避免这一问题，需要选用非常小的步长。

在解析势能面的模型体系中，原则上可以将电子和原子核的演化使用同一较小的步长，以获得数值稳定性和收敛性。但是在从头算动力学中，这样的做法会导致用量子计算处理大量类似的结构。所以在从头算动力学中，为了让计算效率达到最大，整个动力学过程经常对电子和核的演化采用不同的步长65。即原子核以一定的步长进行演化，每个步长又被分割成若干份来演化电子。至于电子演化时每个时刻所需要的能量和非绝热耦合等物理量可以通过线性插值方式获得。

3.6 退相干校正

在轨线面跳跃方法中，一个比较困难的问题是如何正确描述电子相干。对于这方面内容，文献有详细的探讨28,102,186–191。从电子运动的密度矩阵来看这个问题比较清晰。如果将电子和核放在一起处理，电子体系的约化密度矩阵的非对角元为

即两个电子态上对应核波函数的直积。如果两个电子态势能面差别较大，可以预计该积分会快速趋向于零。但是对于最少面跳跃动力学而言

这样即使离开相互作用区，该项依然不会衰减到零。

不少文献从一个更严格的理论角度看这个问题28,66,102,189,190。电子体系的演化过程会和核发生相互作用，即电子体系为开放体系而非封闭体系，因此电子体系不能使用纯态来刻画，其运动自然也不满足薛定谔方程。因此对其的描述应该是使用约化密度矩阵，但是至今，我们都无法获得该演化方程的准确形式。因此导致该问题缺乏良好的解决方案。

为了解决该问题，也有研究提出了一些修正方案。其中常用的方法之一是 Granucci和Persico186借助了 Zhu和 Truhlar等人28在发展Decay of Mixing动力学时候的方法，提出了对系数采取如下方法进行修正

其中τkl为衰减时间，l为当前所处的电子态，k为其它的电子态，α为参数(常用0.1 hartree)，Ekin为原子核动能。

除此之外，还有一些工作值得注意。Fang和Hammes-Schiffer等人在面跳跃算法中192，考虑当跳跃发生后，重置电子的波函数，让当前态的布局数为1，其余为零，来实现退相干修正。此外，早期的工作(如 Neria和 Nitzan等人193)还曾提出将不同电子态上核波函数利用高斯波包进行展开，通过近似估算这些高斯波包在其中心位置受力后发生的运动，来求解不同电子态上核波函数的重叠，进而获得退相干速率。后来，Schwartz、Bittner、Prezhdo和Rossky等人将这种方法的应用范围进一步拓展，使用高斯波包来估计混合-量子经典动力学中的退相干189,190。近年来，Granucci等人发展的 Overlap Decoherence Correction (ODC)方法187也属于该类方法的范畴。但是对于高斯波包宽度的选择并不是一个简单的问题。Shenvi、Subotnik和杨伟涛基于高斯波包发展了相位校正轨线面跳跃方法100，在此基础上，Shenvi和杨伟涛做了进一步的拓展，使该方法能够通过相位校正解决部分退相干194。另外，Subotnik等人借用了量子-经典刘维尔方程的一些内容，发展了A-FSSH方法102，可以不通过参数设定来直接实现退相干。此外，Prezhdo等人也发展基于退相干诱导面跳跃方法来考虑退相干修正111。近期，方维海和李新奇等人也发展了基于量子测量理论来估算退相干修正的理论方法195,196，并将其用于面跳跃动力学中195。

4 轨线面跳跃方法的应用

4.1 在DNA光诱导反应中的应用

在自然界中，DNA很容易吸收紫外光。但之后，它们往往能够保持相对稳定。这就意味着在该类体系中存在某种特定的机制，能够有效降低荧光产率，以及减少化学反应发生的可能性。这种机制即所谓的“光稳定性”74,75。该光稳定性是自然选择的结果，对于生命的存在和延续具有重要的意义。

在过去的二十多年里，随着实验技术和理论方法的发展，使我们对这一过程有了许多的认识74,75。简单来说，DNA受光激发后，在激发态上可能会运动到势能面交叉点附近，然后跳回基态。之后，体系过多的能量被环境吸收，转化为热能，这样有效防止了 DNA被光反应破坏，维持了“光稳定性”。

对于DNA光稳定性的研究，在实验上可以利用时间分辨光谱技术。由于该技术只能获得反应的时间，无法直接推断出反应的细节。因此，只进行实验研究是远远不够的。如果想要进一步理解实验的数据，就需要配合相关的理论模拟。近年来，出现了很多从各种角度关注该类体系激发态的理论工作。通过这些工作的开展，使我们对DNA光稳定性的本质有了更加深入的理解。其中，基于轨线面跳跃方法可以给出激发态寿命、反应通道和分支比，以及关键的分子运动等信息，因此被广泛的用于DNA研究工作42,57,85,180,197–210。

作为DNA的基本组成单元，单个碱基动力学的重要性不言而喻。光照条件下，单个碱基往往会表现出非常有趣和复杂的激发态衰减现象，因此吸引了很多研究。其中腺嘌呤是一种被大家广泛关注的分子。Babatti和Lischka等人在MR-CIS水平上模拟了气相中 9H-腺嘌呤激发态衰减197，发现该过程主要包括两个步骤，首先发生的是S3→S2→ S1的超快弛豫衰减，需要～22 fs；随后发生的是 S1→S0的弛豫衰减，需要～0.5 ps。其中第二步是通过2E锥形交叉(环上C2的扭曲)，来实现的。但是当使用了 OM2/MRCI后，Thiel等人发现6S1锥形交叉(环上C6的扭曲)是主要通道42。Barbatti等人利用不同的电子结构方法对腺嘌呤动力学重新进行了模拟198，发现从头算MRCIS、半经典 OM2/MRCI以及基于不同泛函的 TDDFT方法给出了不同的结果。尽管前两者都给出了快速的衰减动力学，但是具体的反应通道不同。MRCIS给出了2E是主要的衰减通道，但是活性空间的设置对于非主要通道的所占的比例有影响。OM2/MRCI计算中，不同的计算设置会影响激发态的寿命，但对整个衰减机理的影响很小，依然是6S1通道。在基于TDDFT的模拟中，没有一种泛函能够模拟出令人满意的结果，甚至是激发态寿命也完全不同。后来Barbatti等人通过在ADC(2)水平上进行的面跳跃动力学模拟发现，2E是依然是主要的衰减通道；但是当激发能变的更高时，6S1通道的重要性也开始上升180。值得一提的是，其实不管动力学模拟使用OM2/MRCI、MRCIS还是ADC(2)，都能较好的重复实验观测的寿命211–216。

对于水溶液中腺嘌呤的激发态衰减过程，Thiel等人在 QM/MM (QM = OM2/MRCI)水平上进行了模拟57，发现溶液中的反应和气相类似，超过90%的轨线借助于6S1锥形交叉回到基态，小于10%的轨线借助于2E锥形交叉回到基态。Barbatti在ADC(2)水平上，也研究了腺嘌呤和几个水分子形成的团簇的激发态衰减过程。研究发现9H-腺嘌呤和水形成团簇后，衰减机理基本保持不变。但是他也发现对于其异构体7H-腺嘌呤，其激发态动力学则不一致，溶剂和7H-腺嘌呤之间会发生电子转移，导致激发态衰减的发生199。

Thiel等人采用同样的 QM/MM (QM =OM2/MRCI)方法对腺嘌呤嵌入到 DNA的单链和双链200,201中的非绝热动力学过程进行了研究，发现衰减到基态的时间比真空中和水溶热中的10倍还要长。此外，当嵌入到(dA)10中时，反应和气相和溶液相同，6S1占主导。当嵌入到(dA)10·(dT)10中时，由于腺嘌呤和胸腺嘧啶之间存在氢键，6S1通道被完全抑制，2E变成主要的通道。

除了腺嘌呤外，其他的碱基和其衍生物也获得了很多关注，有很多面跳跃动力学方面的研究85,202,205–210。此外，在DNA体系中，人们已经知道光稳定性还与碱基的配对类型、堆积效应以及嘧啶二聚反应有关。对这些因素的研究，也出现了一些基于轨线面跳跃方法的研究工作。如Groenhof等人模拟鸟嘌呤-胞嘧啶碱基对在气相中和嵌入到 DNA链中的光激发过程203，发现光吸收可以促进 H从腺嘌呤转移到胞嘧啶。Doltsinis也曾研究过该碱基对在气相和水溶液中的动力学，也发现了类似H转移通道60。关于嘧啶的二聚体生成，González等人204也曾用面跳跃动力学描述该类过程中的一些关键反应步骤。

4.2 在光诱导异构化反应中的应用

在自然界中，光诱导异构化反应十分常见。无论在有机共轭分子中，还是在一些功能蛋白中，都有它的存在。在这类反应中往往涉及非绝热过程76–80，近年来，利用轨线面跳跃方法对该类反应的研究取得了一定的成果。研究表明，该类反应大部分是通过C＝C、C＝N或者N＝N双键的扭曲发生的，往往伴随着不同分子构型之间的变换。

Barbatti等人用高精度的CASSCF以及MRCI方法对 CH2NH2+的非绝热动力学过程进行了模拟40。当初始电子态为第二激发态时，激发态弛豫过程主要分为两步。首先是通过C―N键的拉伸，体系从第二激发态回到第一激发态，体系保持一个较高的平面性；然后才是通过 C―N键的拉伸或 C―N键的扭转等运动回到基态。类似的结果也被他人发现，如Tavernelli121,175、兰峥岗39,46等人用 TDDFT方法、Thiel等人217用 OM2/MRCI的方法等。值得一提的是，CH2NH2+是一种质子化的 Schiff碱，一般作为 Retinal的模型体系218来研究，当该类质子化 Schiff碱存在更多的共轭单元时，就可能出现 one-bond-flip以及 Hula-Twist等反应机理218。

Granucci等人184、Thiel等人217以及 Barbatti和Lischka等人219都利用轨线面跳跃方法对乙烯的光异构化过程进行了模拟。当体系被激发后，体系沿着C＝C双键进行异构化，然后两端的C之间会形成双自由基和电荷转移态，导致体系突然极化。双键一端的C由sp2杂化变为sp3杂化，在整个过程中，C发生金字塔化运动。这些结论与Martínez77,78利用AIMS方法研究的结果一致。

Marx和Doltsinis61等人、Granucci38,220等人和Barbatti等人221还对偶氮苯的顺反异构进行了模拟，发现该类分子在光激发后，两种构型E和Z之间存在转换，这一过程涉及旋转、反转、协同反转以及反转与旋转并存四种可能的途径。至于采取哪种途径，主要和体系的初始条件有关。

对于偶氮苯体系，除了有顺反异构以外，Weingart和Thiel等人222发现该分子在光异构化过程存在手性通道，经过此通道的轨线，旋光性会得到保持。在光化学中发现手性导致异构化通道的选择性也在其他体系中有所讨论，如方维海和崔刚龙等人对芳香偶氮吡唑化合物的工作223，以及兰峥岗等人对仿生光开关化合物异构化的研究工作224。

近年来，随着QM/MM技术的发展，对于分子体系光异构化过程，分子体系与外部环境的相互作用也越来越多的被考虑进来。比如，Persico等人225、Marx等人61,226也对偶氮苯体系体系进行了一系列详细的动力学研究，包括其在溶剂和体相材料中的非绝热动力学。另外，方维海和崔刚龙等人227利用基于QM/MM非绝热动力学方法，对偶氮苯在多肽体系的光异构化过程进行了模拟，发现多肽和偶氮苯之间的相互作用不仅能够加速蛋白的折叠和解折叠过程、调控二级结构比例，而且影响偶氮苯的异构化。

4.3 在光诱导金属配位络合物体系中的应用

过渡金属配位化合物在光物理和光化学领域受到广泛的关注。含铜、钌等配位化合物，特别是Cu(I)通常在可见光区域展现出非常有趣和复杂的激发态动力学现象，包括金属-配体电荷转移(MLCT)，较长的三重态MLCT寿命，以及光诱导的配体运动228–231。在光伏材料设计领域，由于Cu(I)具有更廉价和丰富的储量，作为昂贵的Ru(II)和Os(II)配位化合物的可能替代材料，有望在光催化剂，染料敏化太阳能电池(DSSCs)和有机发光二极管(OLED)等材料设计中发挥更大作用231–233。

杜利凯和兰峥岗等人利用基于 DFT/TDDFT的轨线面跳跃方法研究了典型含铜Cu(I)配位化合物([Cu(dmp)2]+, dmp = 2,9-dimethyl-1, 10-phenanthroline)的超快激发态光诱导动力学234。他们发现激发态金属-配体电荷迁移(1MLCT)能够引发[Cu(dmp)2]+复合物配体结构发生扁平化运动。通过全原子尺度非绝热动力学的模拟，可以清楚地分辨出不同分子运动模式以及配体取代基基团在Cu(I)配体扁平化运动过程中的作用。在激发态动力学的早期(30–50 fs)，体系主要发生超快的内转换弛豫到S1态。之后，在接近200–300 fs的时候，开始出现一些结构扁平化运动。在 600–750 fs的时候，配体间的扁平化角度趋于平稳。整个结构扁平化弛豫时间接近675 fs，这一结果与最近的一些实验结论基本一致235–237。不同的实验中，S2衰减到 S1的时间尺度为 47–80 fs，扁平化运动发生的时间尺度约为340–800 fs。因此激发态动力学模拟可以很好的重现实验的观测。此外，模拟也发现了一些分子运动的细节，如通过对 Cu―N键长和N―Cu―N键角的傅里叶变换，可以确定一些最重要的振动频率。另外，模拟发现整个扁平化运动中，两个dmp配体能够保持较好的平面性。对于同样的体系，Tavernelli等人利用QM/MM轨线面跳跃方法也进行了研究238，发现激发态过程受溶剂的影响很小。另外，他们还利用主成分分析的(PCA)方法确定了在激发衰减过程中，配体结构发生扁平化运动和其他一些关键模式起到了主导作用，并使用这些模式建立势能面，进行量子波包动力学的模拟。

对于较重的过渡金属，Tavernelli等人利用TDDFT-QM/MM方法研究了水溶热中含钌配位化合物([Ru(bpy)3]2+, bpy = 2,2′-bipyridine)的激发态光诱导动力学88。他们发现光激发生成单重态后，体系容易发生系间窜越变成三重态，这与实验上的观察十分吻合。另外，他们还发现溶剂能够在弛豫过程中起到很重要的作用。

4.4 在光伏材料中的应用

光伏材料对于能源利用过程至关重要，近年来有机太阳能电池、染料敏化太阳能电池等新型太阳能电池受到了广泛的关注。对于这些新型电池，在光电转换过程中，通常会伴随着超快的激发态能量和电子的转移。因此轨线面跳跃方法也被大家广泛的用于对该类体系的研究。

对于有机光伏体系，Fernandez-Alberti、Tretiak以及Roitberg等做了很多重要的工作104,160–162,239,240。例如，他们利用基于AM1/CIS的轨线面跳跃方法对一类典型枝杈状有机共轭长链分子的激发态衰减过程进行了模拟160,162。在该类分子中，短的共轭枝杈部分被激发后，会发生从短枝杈到长共轭枝杈单元的单向能量转移，导致长的枝杈部分被激发，最终导致体系中共轭最大的枝杈被激发。因而导致最后的荧光产生在具有最大共轭体系的单元上。在整个激发态衰减过程中，电子能量转移是由C≡C伸长导致的。兰峥岗等人后来在该类体系的高精度计算结果也支持该结论241,242。

帅志刚等人利用基于 TDKS-DFTB的轨线面跳跃方法对开环和闭环两种类型的 DPDBF激发态衰减过程进行了模拟243，他们发现在开环的DPDBF中，低频扭转运动能够强烈地耦合电子的运动，有效地耗散体系的能量。然而在闭环体系中，这种运动被化学键限制，导致体系的无辐射衰减速率比开环体系慢的多。他们推测在凝聚态中，开环体系的低频运动也会被限制，因而无辐射衰减的通道将会变慢，最终会出现明显的聚集诱导发光现象。

在有机光伏材料中，由于涉及的体系都很庞大，因此一般会采用 AM1/CIS或 TDDFTB等近似的量化方法与轨线面跳跃方法进行结合，处理该类体系的光诱导反应过程。近年来，随着计算机技术的发展，也出现了一些基于更高精度方法(如TDDFT)的直接动力学模拟。这些方法的引入，使我们能够从更加精确的角度去理解有机光伏材料中的非绝热过程，如激子动力学问题244–247等。但是值得注意的是，使用基于TDDFT的直接动力学处理这些问题的研究需要很大的计算量，因此这类工作依然不多。所以，发展有效的方法处理该类问题就非常重要。

对于一些更为复杂的光伏材料和光催化体系，也有研究人员尝试利用轨线面跳跃方法处理相关问题。Prezhdo等人248–250提出对于拓展体系，电子的激发本身对于体系的整体电子密度影响应该不是很大，因此可以使用经典路径近似，认为激发态受力和基态受力基本一致，即用基态的受力代替激发态的受力。这样，只需要在基态上来计算动力学轨迹，然后沿着该轨迹计算能量和非绝热耦合量，存储其随时间演化的信息，然后再求解面跳跃动力学。这种办法可以用于对非常庞大拓展体系非绝热动力学的研究，因此被他用来研究染料敏化太阳能电池中的电子转移问题，有机光伏体系中的电子转移问题以及光催化、钙钛矿等方面的问题248–255。但是值得注意的是这种处理方法建立的前提是电子激发对体系受力影响不大。此外，由于处理的是扩展体系，往往需要利用单电子图像，在轨道的框架下来计算轨道能量和非绝热耦合矢量。该方法没有建立起电子激发过程对核运动的反馈，因此无法用来描述化学键的生成和断裂作用。在光催化中，它仅仅能够用来理解空穴和电子的转移问题，无法理解涉及到化学变化的过程。

4.5 在其它领域中的应用

除了上述提到的应用以外，轨线面跳跃方法还被用来研究许多广泛的其他问题。例如分子的光解256,257，光诱导质子转移或质子耦合的电子转移过程(如方维海和崔刚龙等人258,259的工作)，化学和生物发光(如刘亚军等人的工作260,261)，一些蛋白(如视觉神经蛋白、荧光蛋白等)中的生色团的光反应62,218等。由于篇幅所限，此处就不一一讨论了。总之，面跳跃的方法极大的推动了我们对于光化学过程的理解。

4.6 发展趋势

近年来，除了发展面跳跃动力学方法本身外，理论化学家也一直在尝试利用基于直接动力学的轨线面跳跃方法去处理一些更加复杂地情况，包括进行多尺度模拟、将研究体系进一步扩大以及实现计算效率的进一步提升等。通过这些工作的开展，能够为轨线面跳跃方法的应用另辟蹊径。

在利用轨线面跳跃方法进行多尺度模拟方面，除了使用QM/MM等方法外，也有工作结合了不同的动力学模拟方法实现尺度跨越。比如将激发态的面跳跃动力学和基态的长时动力学模拟结合，理解光化学和基态反应过程的耦合。这方面的工作，如Doltsinis等人利用从头算分子动力学和经典动力学相结合的方法对二硫键的断裂过程进行了研究59，方维海和崔刚龙等人227对于对偶氮苯在多肽体系的光异构化过程的模拟。多尺度模拟手段也被用于理解分子光开关和多肽的连接体系262，其光反应和基态反应的耦合。此外，一个紧密相关的工作是Martínez等人结合 AIMS处理非绝热动力学和 BO动力学处理后续的基态反应，对H2CSO动力学进行了研究263。虽然没有使用面跳跃，但是思想类似。他们发现当硫化物被激发后，体系在锥形交叉点附近会出现不同的反应路径，生成各种不同类型的产物。这些耦合不同时间尺度的模拟，其实可以帮助我们更好的理解光化学和基态化学过程是怎么耦合的，是该领域发展的重要方向之一。

近年来，也有很多努力将轨线面跳跃方法和更多的从头算方法结合，用于理解各类问题。比如将面跳跃与TDDFTB的方法结合等264,265。其中Mitrić等人利用基于 TDDFTB的轨线面跳跃方法研究了腺嘌呤和水构成的团簇的激发态动力学264，该模拟中整个体系均使用了TDDFTB。Barbatti等人利用类似方法对两种cycloparaphenylene分子在气相中的激发态动力学进行了研究265。另外，他们还对TD-LC-DFTB、TD-DFTB、SD-LDDFTB的动力学表现进行了比较。这些工作对于我们突破现有的量化计算的尺度提出了一个研究思路。

此外针对 TDDFT无法正确描述基态—激发态交叉势能面的问题148，也有工作将Spin-Flip的TDDFT和面跳跃动力学结合261，改善了对内转换回到基态动力学过程的描述。因为DFT类型计算比高精度耗时少，因此该类工作可以使得我们开始讨论一些原来高精度无法讨论体系中出现的动力学。此外，也有工作将 CASPT2与面跳跃直接动力学结合，然后将其用于研究非绝热过程266。因为CASPT2的结果很多时候和实验对比较好，因此该发展对大家精确理解一些小体系的非绝热过程有所帮助。

Elstner等人利用DFTB的相关理论，在单电子图像下，结合QM/MM来模拟DNA体系中的电荷转移问题183。该课题组也曾经利用TDDFTB层面的面跳跃理解有机光伏体系激子扩散过程267。Akimov试图在片段分子轨道和拓展的休克尔理论基础上发展出透热态哈密顿的构建方法，然后将其用于非绝热动力学268。Blumberger等人结合分子力场和局域分子轨道，在单电子近似下，来发展透热态哈密顿，然后基于此发展了 FOB-SH方法113，用于理解电荷转移和输运问题。该类方法由于基于近似的电子结构理论，原则上可以处理包含大量分子构成的聚集体间的电荷输运等非绝热问题。

在提高计算效率方面，Martínez等人将面跳跃，构造透热态哈密顿的方法，QM/MM以及GPU技术相结合，可以用来处理超大体系激发态动力学过程。在最近的一项工作中，他们利用GPU加速在全原子水平上，直接研究了一个光合作用蛋白的激发态能量转移过程269，为理解复杂体系的光化学问题提供了新的研究思路。

在非绝热动力学中，另外一个挑战是处理系间蹿越的问题，即体系的自旋态发生变化。近年来，也有一些工作试图处理该类问题，比如Thiel和崔刚龙等人在 Spin-Diabatic表象下处理系间穿越问题270。González等人在 Spin-Adiabatic表象下处理该类问题44,87,185。

此外，如何分析面跳跃动力学的模拟结果也是一个重要的研究方向。比如，近期，研究者们利用机器学习算法中的的无监督学习方法，实现了通过降维方法理解非绝热动力学中主要的结构变化，自动建立主导反应坐标238,271,272。

近些年相干多维电子光谱在实验上取得了很大的进步273。利用相干多维光谱不仅能够使我们得到传统时间分辨光谱的信息，而且有助于我们理解电子态之间的耦合、电声子耦合等与非绝热动力学有关的重要信息。此外，超快X-ray光谱技术的发展274,275，使得我们可以直接指认很多电子激发以及激发态上结构变化特征，因此利用基于从头算的轨线面跳跃方法，来模拟这类光谱也必将是未来发展方向之一。

5 挑战与展望

轨线面跳跃方法从提出到现在已经发展了四十余年，逐渐成为处理非绝热问题中最受欢迎的半经典方法之一。近些年，在将其与从头算动力学结合上，以及一些实际应用方面，该类研究都取得了很大的进步。此外，值得注意的是近年来国内理论化学界在基于直接动力学的轨线面跳跃方法研究方面也取得了很多进展。除了文中提及的工作，也有一些其它的重要工作37,223,243,276–280，就不一一列举了。

尽管如此，利用轨线面跳跃方法去处理非绝热问题依然会面临很多的挑战。在理论研究方面。轨线面跳跃方法作为一种半经典的理论，本身并不是一个严格的理论。至少，对于该理论本身，现在无法通过简单的理论方法，建立起回到更为严格理论方法的道路。因此提出的各种校正方法往往都具有一定的随意性，很多时候也包含需要预设的参数。因此发展理论本身是还有很长的路要走。在实际应用方面。由于轨线面跳跃方法处理实际问题的过程中，存在很多细节，每一步出现问题都可能对计算结果产生重要的影响。比如量化方法的选择，不同的计算方法可能给出完全不同的结果。对于演化时间，目前大多数的模拟都限制在皮秒级以内，而一些复杂体系中的反应通常发生在更长的时间尺度，涉及更大的体系上，模拟此类跨越时间空间的多尺度问题依旧还比较困难。此外，因为核的运动是经典处理的，也无法描述隧穿等很多关键问题。

尽管轨线面跳跃方法受到了广泛关注。但我们应该清楚地认识到该方法在处理非绝热问题时依然会面临很多挑战。因此进一步发展更严格的半经典和全量子方法来处理非绝热问题也不失为一种好的选择。这也将是未来理论化学发展的重要方向之一。