超快动力学计算（二）：误差来源与相应的对策

      在此我们先回顾一下上期的几个主要知识点。首先，在绝热近似下，原子核与电子的自由度被分开考虑，因此衍生出了势能面的概念。对于一个确定的电子态，体系的势能是原子核位置的函数。对于基态直到第N激发态，都有属于自己的势能面，而且不同势能面之间不存在跃迁过程。超快动力学涉及到绝热近似失效的情况，不同势能面之间会出现非绝热的跃迁过程。两个势能面之间的非绝热耦合的强度正比于原子核的运动速度的大小，反比于发生跃迁的势能面对应的电子态的能量之差。

      在经典-量子框架下，通常使用每个构型下的绝热波函数作为基组来展开当前构型下含时演化的波函数（实际上也能使用透热态作为基组但是本期先不讨论）。原子核的运动方程仍然遵循经典力学，不同势能面之间的非绝热跃迁即为不同电子态之间的跃迁。主流的方法有两种：平均势能面方法和势能面跳跃方法。

      平均势能面方法又称为含时的自洽场（TDSCF方法），这种方法中电子和原子核互相在对方产生的平均场中演化，相对直观简便且容易引入动力学过程。

      势能面跳跃方法则依靠非绝热耦合矢量和原子速度等来定义体系在不同电子态（此处也能对应势能面）之间的跃迁概率，通过计算多组轨迹取统计平均。轨迹数目的分布满足约束条件：每个电子态上对应的轨迹数目与密度算符在该电子态对应的对角项矩阵元成正比。由于满足上述约束的跳跃方法理论上有无穷种，因此在实际操作中，研究者们更倾向于使用Tully提出的最少势能面跳跃方法（FSSH）。

      值得注意的是，上述两种算法本身并不依赖于电子结构计算的理论级别，换言之，上述的电子波函数可以是解析的理论模型，可以是半经验方法，可以是密度泛函理论（DFT）的单电子波函数，也可以是量子化学计算中的多组态波函数。理论模型一般对应一些简单二能级系统等，从事方法开发的工作者会用这些理论模型测试算法的数值稳定性和鲁棒性。半经验方法一般用于中大型体系的较为粗犷的研究。DFT是平衡了精度和计算量的选择，在凝聚态材料计算中最常用。多组态波函数是精度最高的选择，当研究对象是小分子体系时，如果计算资源允许，最好使用多组态波函数。不过相对的，多组态的使用门槛也很高，对基础理论不了解的使用者很容易白费算力，最后什么有意义的结果都得不到。

      在多组态波函数中，每一个势能面都可以由一个或几个不同电子占据态下的Slater行列式描述。而当波函数变为DFT对应的单电子波函数时，每个势能面对应的则是Kohn-Sham轨道不同的占据方式。在这里我们简单介绍一下当电子结构基于DFT计算的波函数时可能引入的误差及相应的对策：

      细致平衡（detailed balance）是一个非常重要但是又经常被忽视的概念，它代表体系最终应该达到的微观可逆的状态。合理的MQC框架下动力学的演化过程最终应该达到正确的平衡状态，而“正确的平衡状态”的判据即体系是否满足细致平衡，或者说体系中的电子态的占据数是否能满足Boltzmann分布等平衡态分布函数。

      当一个量子系统和一个经典热浴耦合时，平衡状态下的电子向高处和向低处跃迁的概率是相同的，最终导致每个量子态在平衡态的占据数都相同，此时对应的电子温度是无穷。MQC框架下的原子运动遵循经典动力学，因此理论上也存在上述风险。只不过在平均势能面和势能面跳跃中，都有能量守恒这个约束条件，因此可以在一定程度上限制电子态之间的跃迁过程。但是这并不能证明上述两种算法下 体系最终能满足细致平衡。

      Tully曾在2005年发表的工作中，通过一个简单的数学模型检验了电子运动对原子核完全没有反馈，平均势能面方法和最少势能面跳跃这三种方法是否能描述细致平衡。对于第一种情况，电子占据数对应的温度毫不意外的趋于了无穷。平均势能面方法虽然考虑了电子结构演化对原子核运动的反馈，但是却与玻尔兹曼分布存在一个显式的误差，该误差的大小与二能级系统能级之差的大小成反比。对于FSSH，只要能根据势能面跳跃后的能量变化沿着非绝热耦合矢量的方向修正原子的速度，则体系自然可以包含细致平衡；如果体系沿着非绝热耦合矢量方向的动能无法满足从低到高的跳跃的两个势能面对应的能量之差，则这个跳跃会被禁止，即“阻挫跳跃”。关于更详细的数学推导可以查看Tully的原文J. Chem. Phys. 122, 094102 (2005)。

      笔者在材料大会做报告时曾被问道了FSSH是否天然的包含细致平衡，当时笔者由于没有从经典路径近似下的FSSH调转过来（经典路径近似下显然需要对跳跃概率再乘以一个因子），所以想当然的说了没有。实际上原版的能更新原子速度的FSSH是天然包含细致平衡的，笔者在此对当时被误导的听众以及提问的何力新老师致以诚挚的歉意。

      要解决平均势能面方法不满足细致平衡条件的缺陷，最直接的方法是引入核量子效应。孟胜课题组曾经以环聚物分子动力学（RPMD）的方式近似引入了核量子效应，结合TDAP的平均势能面方法计算了单层石墨烯的载流子动力学过程并成功再现了与实验复合良好的空穴弛豫时间。

      在平均势能面方法中引入核量子效应会导致计算量倍增，而实时修正速度的FSSH对于凝聚态材料体系也非常昂贵。因此研究者们考虑如何在MQC框架下也能给平均势能面方法引入细致平衡，其中最具代表性的工作是汪林望老师在PWmat中实现的Boltzmann TDDFT。该方法定义了含时演化的波函数在两个态之间跃迁的载流子流，当载流子从高能量流向低能量时，不会受到任何的阻碍；而从低能量流向高能量时需要经过一个与退相干时间有关的修正。这种修正最终会影响体系的密度矩阵并改变含时波函数的演化趋势。同时，修正后的能量差也会通过合适的方法反馈给原子核，这样就维持了能量守恒的约束条件。在光致相变，非热熔化，电子束辐照分解中，这种方法带来的载流子随时间演化趋势都能与实验结果符合良好。

      一束包含了所有电子态的波包，随着时间的演化，可能会出现若干分支，而这些分支随着时间的演化会逐渐失去彼此的相干性。例如，体系的波函数原本在一个势能面上演化，当经过势能面的交点时，由于势能面交点附近的非绝热耦合很强，体系的量子态理论上会分裂到两个势能面上，随着时间的演化，当两个势能面隔得足够远时，两个势能面上的波函数分支不再相干并沿着各自的轨迹独立运动。显然，平均势能面方法中根本没有这种物理图像。在FSSH中轨迹分支看起来很自然，但是轨迹分支不等于退相干，FSSH中分支后的波函数存在过相干的过程，当波函数反复经过非绝热耦合很强的区域时，误差还会积累。

      最直接的引入退相干过程的方法仍然是引入核量子效应，当原子核的波函数交叠趋于0时，退相干自然发生了。孟胜课题组通过RP-TDAP对平均势能面方法引入核量子效应并成功模拟了臭氧分子的环态和基态的退相干过程。

      Prezhdo提出了退相干诱导的势能面跳跃方法（DISH）。该方法需要计算出展开含时演化波函数的所有绝热轨道各自的退相干时间，然后根据退相干时间产生一个泊松分布的随机数。通过该随机数与所有绝热轨道距离上次退相干事件发生的时间的相对大小，判断体系的含时波函数是否会发生退相干。Truhlar则对轨迹势能面跳跃（TSH）分别发展了基于波包和基于能量的退相干修正方法。

      汪林望老师在PWmat中实现了基于平均势能面方法的自然轨道分支（NOB）方法。该方法用密度矩阵的形式重写了含时的薛定谔方程，此时退相干过程可以表现为密度矩阵的对角项衰减为0。对密度矩阵作对角化则可以得到自然轨道和对应的本征值，自然轨道也可以用绝热的Kohn-Sham轨道展开。根据密度矩阵的本征值，可以定义一个熵。通过熵与给定阈值的相对大小来决定退相干过程是否发生。由于体系的波函数坍缩到某个自然轨道时，仍然可以写作多个Kohn-Sham轨道的线性组合，相比于DISH直接将一个Kohn-Sham轨道从含时波函数中剥离出来，NOB显得更自然。目前该方法已经被用于模拟冷冻电镜下分子的辐照分解过程，并证明了冷冻电镜的低温并不能抑制分子被辐照后的电离过程（详细内容可参考往期链接🔗）。

      上期提到过，MQC有一种特殊形式，即原子核的运动轨迹一直处在基态势能面上演化，电子的含时演化不会对原子核的运动产生反馈。在这种形式下做FSSH只需要预先用BOMD跑出一个基态的轨迹，计算轨迹上所有构型的波函数和本征值，然后直接构造非绝热耦合项并进行模拟。在读过本期前3部分的内容之后，读者可能会产生疑问，这种不修正原子核运动的FSSH不可能满足细致平衡条件，甚至可能都不满足能量守恒。确实如此，所以针对这个问题，Prezhdo将阻挫跳跃核速度标定近似转换为势能面跳跃的概率公式乘以一个玻尔兹曼因子，以此让体系在数学上依然满足细致平衡，至于物理上，仁者见仁智者见智吧。该方法也是PYXAID乃至Hefei-NAMD的基石。

      PWmat也能实现经典路径近似下的载流子动力学演化。注意到实际操作中非绝热耦合项可以直接由当前时刻和下一个时刻的波函数内积除以时间步长得到，PWmat在做BOMD产生原子核轨迹的过程中，直接输出了两个时刻间波函数的内积，这样就避免了再对轨迹上的构型做自洽计算求解波函数的过程，既节约了计算时间又节约了存储空间。在处理退相干和细致平衡的方法上，PWmat并没有采用FSSH，而是使用了由汪林望老师发展的P-matrix方法，该方法也是将含时薛定谔方程采用密度矩阵形式进行重新表述，通过密度矩阵非对角项衰减来描述退相干过程。此外，这些非对角项被分为两部分，分别对应能量增加和能量降低的电子跃迁。通过将玻尔兹曼因子应用于能量增加的跃迁来实现细致平衡。在P-matrix的加持下，只需要进行一次非绝热模拟就能得到和FSSH多组轨迹的统计平均一样的结果。

      正如本期第2和第3部分提到的，引入核量子效应就能解决MQC框架下的细致平衡和退相干问题。但是引入核量子效应之后的计算量会显著增大，因此研究者需要分清楚自己研究的体系是否必须引进核量子效应。

      首先，引入核量子效应并不意味着体系的电子结构需要用含时薛定谔方程描述；核量子效应也不一定要通过路径积分等昂贵计算才能引入。比如钟志诚课题组在研究SrTiO₃的量子顺电的问题中，就是只使用了基态的势能面训练的机器学习势函数，结合朗之万的量子热浴，就取得了非常好的效果。在这类问题中，核量子效应主要是用来描述重核被忽略的零点振动。

      与质子有关的过程则会复杂很多，质子的运动速度很快，往往不能用绝热近似描述。质子带来的量子隧穿效应往往会伴随着超快的电荷转移，此时必须使用非绝热的电子演化和核量子效应相结合的计算方式。

      对于常规体系，如果核量子效应只是用来补足MQC在描述细致平衡和退相干过程中可能的缺点，则完全可以用PWmat的Boltzmann TDDFT和NOB等方法平替。