基于位点基底的电荷与激子输运非绝热分子动力学中静电势能的高效计算方法

《Journal of Chemical Theory and Computation》：Efficient Calculation of Electrostatic Energies for Large-Scale Nonadiabatic Molecular Dynamics in a Site Basis

本文提出了一种结合阻尼位移力（DSF）实空间静电求和与加减法方案的高效计算策略，用于位点基底非绝热分子动力学（NAMD）模拟。该方案将计算Nmol个不同电荷态静电能的复杂度从O(Nmol)降至O(Nmol0)，实现了与系统规模无关的恒定计算开销。在蒽（ANT）晶体空穴输运模拟中，该方法证实了静电势能涨落（对角静电无序）会降低空穴离域度（IPR）和迁移率，使模拟结果与实验值更为吻合，同时维持了瞬态离域（transient delocalization）的输运机制。

　　

引言

在凝聚态体系的经典分子动力学（MD）或蒙特卡洛（MC）模拟中，长程静电相互作用是描述系统结构和功能的关键，但也是计算成本最高的能量项。Ewald求和及其近线性标度版本（如P3M、PME、SPME）是周期性边界条件下模拟的标准方法。然而，当系统存在多种电荷或激发态时，需要在每个MD时间步或MC移动中为每个状态计算静电能量，这极大地增加了计算负担。

本文作者所在课题组近期开发了一种用于分子材料中电荷输运的非绝热分子动力学模拟方法，称为基于片段轨道的面跳跃（FOB-SH）。该方法模拟了晶体或非晶分子组装体中电荷载流子（过剩电子或空穴）的量子力学演化。对于一个包含N_mol个分子的组装体，它需要计算N_mol个电荷态的经典力场势能（即“位点能”），这些态仅在携带电荷的分子上有所不同。直接计算N_mol个位点能需要计算N_mol个Ewald求和，这对于真正纳米尺度的系统（N_mol通常大于100）来说在计算上是不可行的。

本文旨在介绍一种高效的计算方案，用于计算N_mol个不同电荷态的静电位点能和力，我们称之为“加减法方案”。该方案首先计算所有分子处于中性电荷态时的静电位点能和力，然后计算校正项，以解释在每个态中有一个分子带电而非中性的情况。研究表明，当加减法方案与实空间静电求和技术（如阻尼位移力（DSF）方法）结合时，计算N_mol个电荷态相对于单个电荷态的计算开销是系统规模无关的常数，即O(N_mol⁰)，而非O(N_mol)。这一进展使得我们能够在电荷和激子输运的非绝热分子动力学模拟中包含完整的静电相互作用。

方法

Ewald求和

为了设定背景，我们给出了周期性复制点电荷集合静电能的Ewald求和标准表达式（原子单位制），假设为导电“锡箔”边界条件：

E^EW= E^real+ E^rec+ E^self

其中，E^real是实空间求和，E^rec是倒空间求和，E^self是位置无关的能量自项。

阻尼位移力（DSF）方法

DSF方法是一种Ewald求和的替代方法，其优点在于仅需要成对实空间求和，无需倒空间求和。DSF静电能量E^DSF由下式给出：

E^DSF= (1/2) Σ_i=1^N_atΣ_j=1^N_atΣ_n' q_iq_jR^DSF(|R_ij+ n|)

其中，R^DSF是DSF静电势函数，Θ是Heaviside函数，截断所有距离大于R_cut的原子对相互作用。

多电荷态：加减法方案

考虑一个包含N_mol个相同分子的大型周期性超晶胞，每个分子由N_apm个原子组成，总原子数为N_at= N_molN_apm。每个分子可以存在于两种电荷态：中性态和带电态。系统可以存在于N_mol个不同的电荷态中，其中在电荷态I中，分子I带电，所有其他分子J ≠ I为中性。

使用加减法方案，对于给定的电荷态I，实空间求和E_I^real可以写为：

E_I^real= E_neutral^real+ ΔE_intra^real+ ΔE_inter^real

其中，E_neutral^real是中性系统的实空间静电能量，ΔE_intra^real和ΔE_inter^real分别是分子内和分子间校正项，用于解释分子I带电而非中性的事实。计算E_neutral^real的计算量标度为O(N_atN_at^cut)，其中N_at^cut是给定原子在半径为R_cut^EW的球内的平均原子数。计算校正项的计算量标度为O(N_apmN_at^cut)。E_neutral^real项只需计算一次，而校正项需要为N_mol个不同电荷态计算，计算量为O(N_molN_apmN_at^cut) = O(N_atN_at^cut)。因此，N_mol个实空间求和的计算量标度为O(N_atN_at^cut)。同样的标度考虑也适用于实空间力。

对于倒空间部分，使用加减法方案，倒空间能量E_I^rec可以写为：

E_I^rec= (2π/V) Σ_k≠0(1/|k|²) exp(-|k|²/(4η²)) × |Σ_j=1^N_atq_jⁿexp(i k·R_j) + Σ_i∈I^N_apm(q_i^c- q_iⁿ) exp(i k·R_i)|²

等式右边的第一个求和（对N_at个原子）仅需计算一次，计算量为O(N_kN_at)，其中N_k是求和中包含的倒空间矢量数。第二个求和（对N_apm个原子）需要计算N_mol次，计算量为O(N_kN_apmN_mol) = O(N_kN_at)。因此，N_mol个倒空间求和的计算量标度为O(N_kN_at)。然而，倒空间力的情况则不同，其计算量标度不利，为O(N_kN_at²)，这使得计算N_mol个倒空间求和在计算上变得极其昂贵。

结合DSF的加减法方案

DSF中的相互作用项与Ewald求和的短程实空间部分形式相似，表明DSF方法在电荷态数量方面应表现出与实空间Ewald求和同样良好的标度行为。使用加减法方案，DSF静电能量和力可以以O(N_atN_at^cut)的计算量计算。因此，DSF方法允许以相对于单个电荷态计算而言很小的开销来计算N_mol个电荷态的静电相互作用。计算N_mol个电荷态静电能量的计算量仅仅是单个电荷态计算量的系统规模无关的常数倍，而不是N_mol倍。

FOB-SH非绝热分子动力学

在FOB-SH中，有机半导体（OSs）的价带（或导带）由以下哈密顿量描述：

H = Σ_k^N_mol?_k|?_k???_k| + Σ_k≠l^N_molH_kl|?_k???_l|

其中，?_k是给定分子k上的正交化HOMO（LUMO），用于空穴（电子）输运，?_k是位点能，H_kl是?_k和?_l之间的电子耦合。N_mol个位点能?_k是使用参数化的经典力场计算的，其中分子k处于带电态，所有其他分子l ≠ k处于中性态。当包含静电学时，这构成了一个计算瓶颈，我们旨在通过使用上述新的DSF静电学加减法方案来消除这个瓶颈。

模拟细节

非极化力场

蒽分子的中性态由标准GAFF力场参数建模。为了模拟蒽分子的带电阳离子态，首先通过q_i^c,FF= q_i^n,FF+ Δq_i^DFT获得一组新的力场电荷，其中Δq_i^DFT= q_i^c,DFT– q_i^n,DFT。然后，通过适当缩放相应的DFT键位移Δr^DFT= r^c,DFT– r^n,DFT，并将其添加到中性态的力场平衡键长中，来调整带电态的力场平衡键长，r^c,FF= r^n,FF+ β Δr^DFT。缩放常数β经过调整，直到从力场获得的真空中两个蒽分子之间空穴自交换的4点重组能（λ^4p）与DFT计算获得的值匹配，得出β = 0.88。

通过电荷缩放的隐式极化力场

具有固定原子点电荷的力场无法解释电子极化和介电屏蔽。这导致电子转移的外层或外部重组能和活化能被高估，以及位点能的热涨落被夸大。本文采用了一种简单的电荷缩放方案，其中通过缩放常数γ缩放上述获得的固定原子电荷，q_i^n,FF→ q?_i^n,FF= γ q_i^n,FF和q_i^c,FF→ q?_i^c,FF= γ q_i^n,FF+ Δq_i^DFT。缩放常数经过迭代优化，直到从MD模拟获得的蒽晶体中两个相邻分子之间空穴转移的总重组能λ^st与使用下面详述的极化力场从MD模拟获得的总重组能匹配。确定最佳值为γ = 0.80。

极化力场

极化力场采用了与非极化力场相同的键合相互作用、点电荷和范德华相互作用。电子极化性通过诱导原子点偶极子来解释，各向同性原子极化率取自Amoeba极化力场（α_H= 0.496 ?³, α_C= 1.334 ?³）。

结果

DSF与Ewald求和的比较

DSF能量和相应的力在CP2K包中实现。研究发现，对于所有α值，能量和力的误差最初随着截止半径R_cut的增加而急剧减小，然后饱和。在最大的截止半径R_cut= 35 ?和α = 0.05时，误差最小，小于1%。这表明，通过适当选择α，静电能量和力可以非常接近精确Ewald求和的结果。然而，35 ?的大截止半径在计算上是昂贵的。研究发现，在R_cut= 16 ?和α = 0.09 ?^–1时，可以在精度和计算成本之间达到合理的折衷，此时能量和力的平均误差小于5%，最大误差小于30%。这些值用于下面介绍的所有DSF计算。

接下来，研究了DSF方法在描述静电能量和力的热涨落方面的准确性。DSF和Ewald给出了非常相似的能量涨落，标准差σ分别为134和138 meV。DSF和Ewald能量之间的差异始终低于约20 meV，证实了在使用R_cut和α的最优参数时，DSF相对于Ewald具有良好的准确性。对于力的比较，DSF和Ewald力之间的差异在所有时间都保持在最大差异<0.025 mHa bohr^–1以下，这比几何优化中典型的力收敛阈值低一个数量级。

最后，分析了DSF实现的总能量守恒。在NVE系综中进行的经典MD模拟显示，DSF和Ewald模拟中的能量漂移均保持在10^–11Ha atom^–1ps^–1的数量级，总能量的相对涨落（σ(E_tot)/?E_tot?）对于DSF和Ewald分别为0.04%和0.03%，这意味着DSF算法被正确实现。

系统规模和加速的标度

为了研究DSF计算的计算成本相对于系统规模的标度，考虑了一系列规模递增的超晶胞，并计算了使用DSF、精确Ewald求和和SPME计算单个中性电荷态的静电能量和力。正如预期的那样，DSF计算随系统规模线性标度（在对数（时间）-对数（系统规模）图上斜率为1.0），与SPME类似（斜率为0.9）。相反，精确Ewald求和的标度不利（斜率为2.8）。

接下来，考虑计算系统N_mol个不同带电态的静电能量和力的时间作为超晶胞规模的函数。在不使用加减法方案的情况下，相对于单个电荷态计算的计算开销简单地增加了N_mol倍。将加减法方案仅应用于SPME和Ewald的实空间部分并不能改变这种情况。瓶颈在于倒空间求和，它不能使用加减法方案直接优化。相比之下，使用DSF结合加减法方案，所有N_mol个电荷态的能量和力计算仍然随系统规模线性标度（斜率=1.2）；仅仅预因子或在对数-对数图上的截距改变了约5倍。

为了收敛FOB-SH模拟中的电荷迁移率，通常需要数百个分子的系统规模。根据计时结果，在给定的MD时间步计算所有N_mol= 378个电荷态，使用蛮力Ewald求和将花费超过1000秒，使用SPME将花费约20秒，而使用DSF结合加减法方案仅需约0.5秒。显然，蛮力Ewald甚至SPME计算N_mol个电荷态很快变得不切实际，而DSF结合加减法方案使得这些计算以线性标度且具有合理小的预因子变得可行。

重组能

总重组能是针对ANT晶体中两个最近邻分子之间的空穴转移计算的，这两个分子具有“P”（平行）或“T”（T形）取向，分别与b和a晶向紧密对齐。它包括内部和外部贡献，并通过三种不同方式计算：在0 K时根据4点方案λ^4p，使用超晶胞的几何优化；在300 K时根据平均能隙λ^st；以及根据能隙涨落λ^var.，使用MD模拟。

研究发现，在所有三种重组能定义中，使用DSF或Ewald静电学的非极化MD非常接近。DSF和Ewald对于λ^4p的差异可以忽略不计，对于λ^st的差异小于7%。有限温度导致重组能小幅增加。此外，来自Stokes位移的重组能λ^st和来自涨落的重组能λ^var.相当相似，这意味着系统处于线性响应范围内，其中两种重组能定义完全等价。

总重组能的最佳估计由电子极化MD模拟提供，对于P和T取向，λ^st的值分别为190 ± 8和187 ± 8 meV。从这些值中减去内部贡献λ_i^4p= 142 meV，得到外部贡献λ_o^st= 48和45 meV。电子极化导致外部重组能减少，对于Ewald静电学，沿P和T方向的因子p^MD= λ_o^st(nonpol)/λ_o^st(pol)分别为1.61和1.46。

电子极化对外部重组能的影响也可以从Marcus的连续介质静电学表达式估计为Pekar因子的比率p^c，将光学介电常数设置为1（对应于电子非极化力场）或实验值（对应于（理想）极化力场），p^c= (1 – 1/?_s)/(1/?_op– 1/?_s)。插入蒽的光学和静态介电常数的实验值，?_op= 2.6和?_s= 3.2，得到p^c= 9.5。p^c与p^MD相比被严重高估，指出了连续介质模型的内在局限性，导致重组能对光学介电常数的依赖性过强。相反，观察到上述p^MD值在两个晶向上平均接近?_op^1/2= 1.61。这表明，如果电荷按因子γ = ?_op^–1/2= 0.62缩放，则可以从非极化MD模拟中获得极化MD的重组能估计值。

遵循这一观察结果，缩放固定点电荷，并在使用DSF静电学的隐式极化MD模拟中计算总重组能。研究发现，电荷缩放因子γ = 0.80给出的沿P和T方向的重组能λ^st分别为188 ± 11和191 ± 9 meV，与上述极化MD模拟的值非常接近。该电荷缩放值用于ANT中空穴输运的FOB-SH模拟。

电荷输运模拟

现在转向ANT导电a-b平面中空穴输运的模拟，并使用FOB-SH非绝热分子动力学结合新的DSF静电学加减法方案，研究静电学对空穴波函数离域和迁移率的影响。为此，进行了FOB-SH模拟，其中位点能?_k的静电贡献使用非极化力场和隐式极化力场（γ = 0.80）计算，并将结果与忽略所有静电相互作用（称为“无静电学”）的结果进行比较。

第一个重要观察结果是，在FOB-SH中包含静电学对ANT的堆积结构几乎没有影响。当静电相互作用被排除或分别使用未缩放或缩放的电荷包含时，分子的质心（COM）的RMSD保持不变。这意味着范德华相互作用足以准确再现ANT的堆积结构。

接下来，考虑IPR，它是空穴波函数离域分子数的度量。研究发现，如果排除静电学，所有轨迹平均的IPR约为5。包含具有隐式极化的静电学导致IPR降低至约3.5，对于非极化静电学则降低至约3.0。IPR降低的原因是静电相互作用导致位点能涨落增加（有时也称为对角静电无序），导致空穴载流子波函数在非绝热动力学过程中瞬态占据的能带态（本征态）更加局域化。位能涨落的增加与包含外部/静电贡献时重组能增加一致。

空穴载流子波函数的MSD显示，无论是否包含静电学，MSD都随时间线性增加，表明是爱因斯坦扩散。然而，当包含静电学时，斜率略小，并且这种效应沿a晶向比沿b晶向更明显。将MSD的斜率转换为迁移率，发现当包含具有隐式极化的静电学时，迁移率略有下降，沿b方向从μ_b= 3.5