How to calculate charge mobility in molecular materials from surface hopping non-adiabatic molecular dynamics beyond the hopping/band paradigm

Charge transport in high mobility organic semiconductors is in an intermediate regime between small polaron hopping and band transport limits. We have recently shown that surface hopping non-adiabatic molecular dynamics is a powerful method for prediction of charge transport mechanisms in organic materials and for near-quantitative prediction of charge mobilities at room temperature where the effects of nuclear zero-point motion and tunneling are still relatively small [S. Giannini et al., Nat. Commun., 2019, 10, 3843]. Here we assess and critically discuss the extensions to Tully’s original method that have led to this success: (i) correction for missing electronic decoherence, (ii) detection of trivial crossings and (iii) removal of decoherence correction-induced spurious charge transfer. If any one of these corrections is not included, the charge mobility diverges with system size, each for different physical reasons. Yet if they are included, convergence with system size, detailed balance and good internal consistency are achieved.


Introduction
Organic semiconductors are promising materials for a large range of electronic applications. [1][2][3][4] Their flexibility and tunability are advantageous for organic photovoltaics or diodes, but their charge carrier mobility is still moderate compared to inorganic semiconductors. 5 The experimental and computational quest for OSs with larger mobility is hampered by insufficient knowledge of the charge transport (CT) mechanism in such materials. The parameters determining the CT usually range in a regime where the standard transport theories (band theory and the hopping model) fail. 6 Without an existing model at hand, a direct propagation of electron-nuclear dynamics is required. While improved theoretical models have recently been proposed, e.g., transient localization theory, 7-9 explicit propagation of coupled electronnuclear dynamics is arguably the most promising method to reveal the true nature of charge carriers in organic semiconductor materials. The explicit propagation is founded on rigorous quantum mechanical principles, is free of many of the assumptions that limit the predictitive power of theoretical models, and thus encompasses a wide range of possible transport mechanisms.
Among the numerous existing methods to propagate nuclei and electrons in a non-adiabatic framework at a molecular scale (e.g., ab initio multiple spawning, 10,11 exact factorization, [12][13][14] Ehrenfest dynamics 15 ), Tully's fewest switches surface hopping (FSSH) has become one of the most popular methods. 16,17 The FSSH algorithm relies on the usual molecular dynamics (MD) framework and additionally integrates electronic dynamics explicitly solving the time-dependent Schrödinger equation. By allowing instantaneous vertical transitions (hops) between potential energy surfaces, FSSH also includes the feedback between nuclei and electrons. The classical treatment of the nuclei and the ad hoc -but based on physical arguments -probability for hops between electronically excited states permits a fast, yet accurate, propagation of the dynamics in many situations. The properties required for FSSH (excited state energies, forces and non-adiabatic coupling vectors) can be calculated on-the-fly using time-dependent density functional theory 11,18,19 or semi-empirical Hamiltonians. [20][21][22][23][24][25] The latter approach renders FSSH attractive for study of condensed matter problems and large systems. 26,27 Moreover, FSSH has a number of desirable attributes: it conserves total energy and it was shown to obey detailed balance to a good approximation. 25,[28][29][30][31][32][33] However, this method also has a number of well-known shortcomings and a plethora of variants appeared in the literature since the original work of Tully to address them. [34][35][36][37][38] Among these various intrinsic issues, some shortcomings hinder particularly the simulation of charge transport: (i) the decay of the electronic coherences between adiabatic states (decoherence) is missing in the original formulation, 39 (ii) undetected trivial crossings may lead to unphysical long-range charge transfers, (iii) the common decoherence correction schemes induce spurious long-range charge transfers and (iv) nuclear quantum effects that are particularly important at low temperatures, such as zero-point energy and tunneling, are missing. In this paper we investigate in detail how issues (i)-(iii) affect the FSSH simulation of charge transport in organic materials and we analyse the performance of various correction schemes to mitigate these shortcomings. For inclusion of nuclear quantum effects in the surface hopping simulation of electron transfer, we refer to a recent publication. 40 The aim is to identify the best ''set-up'' or ''flavour'' of FSSH simulations for reliable calculation of charge mobilities in molecular materials. Fig. 1 illustrates issues (i)-(iii). The lack of an inherent decoherence mechanism within FSSH (see Fig. 1(A)) is a wellknown issue, often raised in the literature when studying excited states and relaxation processes. 41 Without a correction to enforce a decoherence, the electronic dynamics is strongly biased and that impacts the charge transfer rate. 42 Numerous correction schemes have been developed, 21,[43][44][45][46][47][48][49][50][51] yet several open questions remain, especially in relation to charge transport simulations. What is the impact of the different decoherence schemes on the equilibrium distribution of states? How important is the decoherence to calculate the electronic mobility using FSSH simulations? On the other hand, the correct detection of trivial crossings has been often an overlooked problem, though it limits substantially the accuracy of charge transport simulations via FSSH. When a trivial crossing occurs, it must be taken care of with an update of state indices, otherwise unphysical charge transfer will occur (see Fig. 1(B)). But such events are often undetected by the original FSSH algorithm, due to the finite MD timestep of the simulation. Different approaches have been developed to detect trivial crossings and to update the state indices accordingly. [52][53][54][55][56][57] Finally, Fig. 1(C) shows another source of spurious transfer recently pointed out by our group 26 and Wang and coworkers. 58 This transfer is induced by the common decoherence correction schemes, and, if not removed, will render any mobility calculation erroneous.
The goal of this paper is to establish the best practice for FSSH simulation of charge transport in real materials by determining the best set-up in terms of decoherence correction, elimination of spurious long-range charge transfers, detection of trivial crossings and appropriate definition of mean-squaredisplacement for mobility calculation. We present here a thorough study of the role of the decoherence correction schemes in both equilibrium and transport properties. We assess the necessity of a state-tracking algorithm to detect and to take care of the trivial crossings and of a correction to remove the spurious charge transfers induced by the decoherence correction schemes. We also discuss the definitions of electronic populations (somewhat ambiguous in FSSH due to the simultaneous propagation of quantum and surface states 59 ) and we compare different definitions for the mean-square displacement calculation. To test these various corrections and set-ups and explore the role of decoherence and state-tracking algorithms in large molecular The lack of decoherence biases the electronic wavefunction C(t) that becomes inconsistent with the active state wavefunction c a (t) after passing through an avoided crossing. We discuss different decoherence corrections in Section 2.2. (B) A trivial crossing of two potential energy surfaces leads to unphysical long-range charge transfer. We use a state-tracking algorithm to detect trivial crossings and reorder the states (see Section 2.3). (C) Decoherence correction-induced spurious charge transfer (DCICT). We develop a spurious charge transfer correction (SCTC) algorithm described in Section 2.4. systems, we needed an efficient method to propagate surface hopping trajectories. We relied on our recently developed fragment-orbital based surface hopping (FOB-SH), 25,60 a semiempirical approach that is designed to determine efficiently and yet accurately the electronic Hamiltonian and nuclear derivatives (forces, non-adiabatic coupling vectors) in large organic crystals. We conclude that a combination of a decoherence scheme, a trivial crossing detection and a correction of spurious long-range transfer permits converging the electronic mobility with system size and MD timestep. This paper is organized as follows. In Section 2, we provide a short summary of the FOB-SH method followed by a discussion of the different existing decoherence schemes, a detailed description of the trivial crossing issue and the state-tracking algorithm used in this work. We also outline the spurious charge transfer correction developed by us. We present the alternative electronic propagation proposed by Hammes-Schiffer to deal with the presence of forbidden transitions (elimination of classically forbidden hops, EFH) and we discuss the different electronic populations commonly used in the literature and the various MSD definitions. We then provide in Section 3 the details of the molecular systems (force field and parameters) we have used to study equilibrium and dynamical properties. Our results are discussed in Section 4: we first investigate the impact of decoherence on equilibrium properties for chains of ethylene-like molecules of different lengths and test the EFH propagation to improve the internal consistency. We then focus on the effect of the state-tracking algorithm and the decoherence correction scheme on electronic mobility and the inverse participation ratio (IPR, a measure of the size of the charge carrier) in an embedded chain of real anthracene molecules. We conclude our work in Section 5.

Fragment-orbital based surface hopping
To simulate the transport of excess charge carriers in molecular systems, we have recently developed an efficient fragment orbital based surface hopping (FOB-SH) framework. In this section, we will summarize the main elements of FOB-SH and its strengths. In Fig. 2 we report a simplified flowchart of the FOB-SH algorithm to help in discussing its fundamental features. For a more detailed consideration of FOB-SH, we refer to our previous publications. [24][25][26] The FOB-SH method is based on the following assumptions. (i) The full many-body electronic wavefunction is replaced by a one-particle wavefunction C(t) for the excess charge carrier. (ii) The excess charge wavefunction C(t) can be expanded in a localized, non-orthogonal basis set made up of fragment orbitals that mediate its transport in the system (usually these orbitals are SOMOs of the isolated molecules {j m }). C(t) takes the form: where R is the 3N vector of nuclear positions and M is the number of fragment orbitals mediating the charge transfer. (iii) The electronic Hamiltonian in the basis of the fragment orbitals is represented by a parametrized tight-binding Hamiltonian. To facilitate the forthcoming propagation equations, Löwdin orthogonalization of the basis set {j m } is applied to define the orthogonal localized basis set {f l }, where T ml = [S À1/2 ] ml , with S the overlap matrix of the fragment orbital basis set ( % S ml = hj m |j n i). The excess charge wavefunction is now: Inserting eqn (3) into the time-dependent Schrödinger equation, one obtains where H kl = hf k |H|f l i, with H the electronic Hamiltonian and d kl = hf k | _ f l i the non-adiabatic coupling elements (NACEs) of the localized orthogonal basis set. As those NACEs are generally close to zero, we label the orthogonal localized basis as a diabatic basis.
To carry out simulations on large systems and long time scales, we designed a parametrized approach to determine the electronic Hamiltonian H kl , thus avoiding explicit expensive electronic structure calculations. The diagonal elements H kk = hf k |H|f k i, which corresponds to the energy of a charge localized on molecule k, are calculated via a classical force field where molecule k is charged and all the other M À 1 molecules are neutral. The off-diagonal elements H kl = hf k |H|f l i, which correspond to the electronic coupling matrix elements or transfer integral, are calculated using our recently developed analytic overlap method (AOM). 61 This method relies on the assumption of a linear relationship between off-diagonal elements H kl and % S kl (i.e., the overlap between the fragment orbitals (j k ,j l ) projected into Slater-type functions), namely H kl = C % S kl . C is a fitting parameter and can be obtained by correlating the overlap % S kl with high quality DFT calculations. 61 This method allows the calculation of H kl for a cost several orders of magnitude lower than that of explicit electronic structure calculations. It was found that errors are less than a factor of 2 with respect to reference coupling values, obtained with approximate coupled cluster (SCS-CC2)/Generalized Mulliken Hush calculations, 62 which spanned 5 orders of magnitude. We refer to our previous paper 61 for a more detailed description of the AOM. It is worth noting the analogy between the calculation of the FOB-SH electronic Hamiltonian and the empirical valence bond approach of Warshel and co-workers, 63 where the electronic Hamiltonian is also built from the classical force field for the diagonal elements and different parametrizations for the off-diagonal elements. As indicated in Fig. 2 this Hamiltonian is a key feature of the FOB-SH method that allows fast computation.
Besides the electronic Hamiltonian matrix elements, the timedependent Schrödinger equation (eqn (3)) requires the determination of the NACEs d kl which can be related to the non-orthogonal Both D kl 0 and _ T are obtained from the finite difference between t and t + Dt. We have taken special care that the nuclear positions R t and R t+Dt are translated within the centerof-mass frame at each timestep. Otherwise, in the case of nonzero center-of-mass nuclear velocities, the overlap element hj k (R t )|j l (R t+Dt )i would have arbitrary values. Fig. 2 Scheme of the FOB-SH (fragment orbital-based surface hopping) algorithm. Different colors represent improvements of the algorithm necessary to fulfil: trivial crossing detection, detailed balance and energy conservation and internal consistency. In red are reported modifications and properties that have been analysed in the present work. Symbols are defined according to equations in the text. RK: Runge-Kutta algorithm, EFH: elimination of forbidden hops, AOM: analytic overlap method, SC-FSSH: self-consistent fewest switches surface hopping, FSSH: fewest switches surface hopping, NACV: non-adiabatic coupling vector, SCTC: spurious charge transfer correction, MSD: mean squared displacement (eqn (30)), IPR: inverse participation ratio (eqn (35)).
We now turn to the propagation of the nuclei according to the velocity-Verlet as shown in Fig. 2 and to the force calculation. In the FSSH algorithm, the nuclei evolve on one adiabatic energy surface E a (where E a = [H ad ] aa , with H ad = U † HU and U the unitary transformation matrix that diagonalizes H to H ad ). It is worth noticing that before starting the electronic integration, the phase of the eigenvectors forming U must be checked and made consistent along the trajectory for an accurate calculation of _ U at a later stage of the algorithm (eqn 8). Since H and U are real, this amounts to a check of the sign of the eigenvectors (''check sign'' in Fig. 2). 64 From U, we can define the adiabatic U ki f k , that form the adiabatic basis.
The nuclear force acting on nucleus I is F I,a = Àr I E a and can be obtained from the Hellmann-Feynman theorem: where [r I H] kl = r I hf k |H|f l i. We refer to our previous paper 25 for the derivation of eqn (6). The gradients of the electronic Hamiltonian matrix diagonal elements are obtained directly from the classical force field, whereas a finite difference approach is used for the matrix off-diagonal elements based on the AOM. 60 We note that the finite difference for the off-diagonal gradients requires an order of N atom M calculations of H kl elements that would make explicit electronic structure calculations unaffordable. The nuclear forces on a given adiabatic state a obtained in eqn (6) consist of a linear combination of the diagonal and off-diagonal forces on the diabats, with a weighting that is proportional to the projection of the adiabats on the diabatsthe weighting takes into account the effect of charge delocalization on the adiabatic forces. Finally, the core of the FSSH method is the choice of active surface E a on which the nuclei evolve and the feedback of the electronic dynamics onto the nuclear motion. In Tully's approach, 39 the active surface is decided in two steps: (i) a new state is chosen via a stochastic process and (ii) the energy conservation requirement is applied to determine whether the change in active state is energetically possible. The stochastic process (i) is based on the hopping probabilities calculated at each timestep t between the active surface and all the other states j: where d ad ja = hc j | _ c a i are the adiabatic NACEs, which are calculated from the diabatic NACEs (d kl [D] kl ), The adiabatic coefficients c j are the expansion coefficients of the electronic wavefunction in the adiabatic basis, CðtÞ ¼ P M i¼1 c i ðtÞc i ðRðtÞÞ. The probability to remain in state a is simply g aa ¼ 1 À P jaa g ja . After the calculation of the probability g ja , a random number is drawn to decide whether a hop can be attempted to a new state n. If so, the following condition should hold to ensure energy conservation, where E a and E n are the potential energies and T a and T n are the nuclear kinetic energies before and after the hop. When a hop is attempted from state a to state n, all quantities E a , T a and E n are already known. To ensure eqn (9) is satisfied, the nuclear kinetic energy (i.e., the nuclear velocities) must be adapted.
Based on the theoretical studies of Pechukas 65 and Herman, 66 Tully prescribed to adjust the velocity component in the direction of the non-adiabatic coupling vectors (NACVs) d ad I,an = hc a |r I c n i. 39 But if there is not enough kinetic energy along the NACVs to satisfy eqn (9), the hop is rejected, the active state remains in state a and the velocity components along the NACVs' direction are reversed. 67 To apply the NACV-oriented adjustment in the FOB-SH framework, we have derived an exact expression for the NACVs in terms of available nuclear gradients in the diabatic basis as well as an efficient approximation. 25 However, the rejection of hops causes an inconsistency in FSSH populations. Tully's hopping probability (eqn (7)) was designed to ensure for a model two-state system that -on average -the wavefunction C(t) is similar to the adiabatic active state. On the other hand, the energy conservation criterion leads to rejecting some classical nuclear hops along the dynamics. Without any correction, the electronic wavefunction will over-populate excited states that are high in energy and therefore unreachable for the classical nuclei. This yields the so-called FSSH internal inconsistency, i.e., a divergence between C(t) and c a (t) (as illustrated in Fig. 1(A)).
Due to this internal inconsistency, two different adiabatic populations coexist in the FSSH algorithm, the quantum amplitude averaged over the trajectory, h|c i (t)| 2 i trj and the surface population, where a n (t) is the index of the active state at time t of trajectory n. The internal inconsistency of FSSH leads to a divergence of those two adiabatic populations. In Sections 2.2 and 2.5, we will discuss different remedies to correct for this inconsistency of FSSH.

Decoherence corrections
The lack of an inherent decoherence mechanism is a longstanding issue of FSSH, already mentioned in Tully's original paper 39 and often advocated in the literature when studying excited state dynamics and relaxation processes. 41 After leaving an avoided crossing where the adiabatic states mix, the full electronic-nuclear wavefunction splits into two sub-wavepackets w i (R) and w j (R), which evolve on different adiabatic surfaces E i and E j . Immediately after the crossing, the center and/or the phase of each wavepacket may diverge in phase space, decreasing the wavefunction coherence between surfaces i and j, Ð dRw i Ã ðRÞw j ðRÞ. This effect is not taken into account in standard FSSH, where the coherence term (i.e., c i *c j ) remains finite. The lack of decoherence ruins the dynamics of the system, leading to the failures of FSSH for some important processes. Rossky and co-workers 41 found that, in the absence of decoherence, decay rates from excited states to the ground state are too fast, yielding incorrect excited state dynamics. Landry and Subotnik 42,68 have later shown that the decay in the charge transfer rate between two molecules obtained with FSSH does not follow the behaviour predicted by Marcus theory.
Since the pioneering work of Rossky and collaborators, 43 numerous correction schemes have been suggested in the literature to tackle the decoherence problem. The most common can be divided into three main categories. (i) Collapsing approaches, in which the electronic wavefunction is reset to the active state C(t) = c a (t) when a given criterion is fulfilled. Criteria suggested in the literature rely on collapsing events after each attempted or successful hop, after each successful hop 50 or when the adiabatic NACEs fall below a threshold. 44 (ii) Exponential damping approaches, in which all non-active adiabatic populations c i are damped at each time step c ic i exp(ÀDt/t ia ), while the active state population is scaled to ensure norm conservation. t ia is the decoherence time. 45,47 (iii) Stochastic damping approaches that rely on random numbers to determine whether the wavefunction is collapsed. 42,51,69 In the last category, each component of the wave-vector containing the expansion coefficients c i with i a a (where a is the active state index) is reset to zero whenever the collapsing probability is larger than a given random number (Z A [0,1]) drawn at each time step. The relative population is transferred to the active state in order to conserve the norm. Within this method, the probability of a collapsing event can be expressed as g collapse i = Dt/t ia in which Dt is the MD timestep. A longer decoherence time t ia results in a lower probability to collapse g collapse i .
As far as we know, no exact expression was derived in the literature to calculate the decoherence time t ia in the context of mixed-quantum classical approaches. However, different formulations were either proposed based on physically grounded justifications 45,46 or derived using approximations for the evolution of nuclear wavepackets. 41,43,48 More recently, using controlled approximations, a decoherence time has been derived from quantum classical Liouville equation (QCLE) formalisms. 36 Those expressions rely on the absence of decoherence when the potential energy surfaces are close to each other or when nuclei are fixed. The energy baseddecoherence time (EDC) proposed by Persico and Granucci (starting from an original expression suggested by Truhlar and co-workers 45,46 ) has the aforementioned characteristics and is widely used in the literature: Here T a is the nuclear kinetic energy and C 0 and E 0 are parameters to determine. We note that the system size implicitly enters into eqn (11) through the nuclear kinetic energy T a (an extensive quantity). Therefore, we suggest normalizing the nuclear kinetic energy by the number of degrees of freedom involved in the FSSH algorithm T a -T a /N dof . By taking the first term of eqn (11), we obtain i.e., the fastest decoherence time possible (due to the Heisenberg uncertainty principle) and free of any ad hoc parameters.
Other expressions for the decoherence time, derived for condensed phase systems and frozen Gaussians travelling on different potential energy surfaces, involve nuclear forces rather than the energy. For instance, Rossky and co-workers 41,43 derived: where the sum goes over the N nuclei of the system, F I i (t) and F I a (t) are the instantaneous forces on decoherent and active states, respectively, and a I is a parameter dependent on the frozen Gaussian width, which has a simple expression in the high temperature limit, a I = 6M I k B T/h 2 , where M I is the mass of the Ith nucleus.
Finally, Subotnik and co-workers developed an extension of FSSH, the augmented-FSSH (A-FSSH), directly from QCLEs to incorporate the decoherence mechanism more rigorously. New dynamical variables are propagated along the nuclear and electronic degrees of freedom to calculate an instantaneous decoherence time. 42,49 Yet propagation with this method is more expensive than FSSH and might not be suitable to study large systems with several hundred molecules.

Trivial crossings and state tracking
The presence of trivial (or unavoided) crossing becomes a substantial limitation in performing charge transport simulations with FSSH. A trivial crossing event occurs when two energy surfaces cross with zero couplings between them, leading to an actual reordering of the state indices (as shown in Fig. 1(B)). Physically, such crossings occur when the adiabats are not interacting, i.e., when the adiabatic wavefunctions are localized in distant regions in space. If the state reordering is not taken care of, the dynamics is continued on the wrong surface, leading to spurious charge transfer that biases any mobility calculation (see Fig. 1(B)). The problem is greatly amplified in systems with many adiabatic states in a narrow energy band (e.g., in OSs), where adiabatic energy surfaces can frequently cross each other.
In practice, the FSSH algorithm is implemented with a finite timestep, meaning that the trivial crossing and index update may be missed. Moreover, the distinction between a trivial crossing and an avoided crossing with a very small energy gap is unclear with a finite timestep. This inherent issue of the FSSH algorithm was often overlooked, as it only arises for systems with many adiabatic states and more strikingly for charge transport. Recently, different solutions emerged in the literature to tackle the missed trivial crossings. Most of them resort to a state tracking algorithm. At each MD timestep, a map is drawn between the indexes of the adiabatic states at time t À Dt and at time t. To build that map, Thiel and co-workers relied on energy criteria and the maximum of overlap between adiabats at time t À Dt and time t, 52 whereas Tretiak and co-workers used the more sophisticated min-cost algorithm. 53,54 Another innovative approach, suggested by Wang and Beljonne, is their flexible surface hopping method (FSH), where the size of the ''active'' region of the OS that transports the charge evolves at each MD timestep. 55 Such an approach not only permits maintaining the relatively small number of adiabatic states (provided that the charge carrier remains localized in space) and diminishing the number of trivial crossings, but also requires new criteria and rules to decide at each MD timestep which part of the OS should be included in the active region. Recently, Wang and co-workers proposed to use the overlap between adiabats at time t À Dt and at time t to classify the surface crossings in different types to determine how to calculate the hopping probability and whether the adiabats' indexes must be updated. 56 They also combined this classification approach with a restriction to hop only to adiabats with a large enough adiabatic population, 57 which strays, however, from the spirit of Tully's original FSSH.
An alternative route would be to improve the calculation of the hopping probability to capture such trivial crossings. A norm-preserving interpolation of the adiabats between time t À Dt and time t can provide a better estimation of the NACV d ad ij . 70 Subotnik and co-workers generalized the norm-preserving interpolation to multiple states crossing using the logarithm of the overlap matrix (eqn (16)). 51 They extended this approach very recently to ensure phase consistency and trivial crossing correction. 71 Wang and Prezhdo proposed a few years ago an alternative straightforward improvement of the probability to hop. 72 They invoked the exact sum rule, to correct the probability to hop to the state the closest in energy, However, there are some issues with eqn (14). First, to derive eqn (14), one needs to integrate eqn (7) from t to t + Dt. If a true trivial crossing is encountered, the adiabatic NACE d ad ij will diverge, the integration will not be permitted and eqn (14) will be invalid. Second, due to the finite timestep, the exact configuration where the surfaces cross is never realized in practice and g ja always remains smaller than 1. Hence, there is no guarantee that the dynamics is continued on the correct surface. Therefore, we suggest that eqn (14) is applied after the trivial crossing problem is accounted for by state tracking/re-indexing of states. In fact, eqn (14) itself cannot correct for a true trivial crossing.
In this work, we opt for a combination of the mapping approach and the self-consistent correction (eqn (15)) for surface hopping. We build the map M between the adiabatic states j at time t and adiabatic states i at time t À Dt with a maximum overlap criterion. First, we calculate the overlap O ij , For each state j = l, we identify state i l with the maximum overlap, |O i l l | = max i |O il |. If |O i l l | 4 1 À e (where e is a constant set to 0.1), we map state l at time t with state i l at time t À Dt, (l) = i l . After that step, all the remaining states j = k at time t that could not be mapped to states at t À Dt (since |O ik | o 1 À e for all unmapped states i) are arranged by index (i.e., by increasing adiabatic energy) and mapped onto one another. As the function map M is a bijection between states at t and states at t À Dt, the reverse map M À1 (which associates states at t À Dt with states at t) is easily found. We can track the index of the active state at t, knowing its value at t À Dt, a t = M À1 (a tÀDt ). This step permits changing the index of the active state without hopping. We stress that our algorithm maps all the states at t with the states at t À Dt, not only the active state, as required by the calculation of the NACEs (eqn (17)). We also note that our mapping criterion produces a unique map and that the algorithm can be run over the states in any order. After the mapping, we make the phase of the eigenvectors consistent along the trajectory by checking the sign of the overlap matrix element O i,M(i) and by reversing the sign of c i if O i,M(i) o 0 (as we discussed in Section 2.1 and underlined in Fig. 2 with the comment ''check sign''). We finally determine the correct hopping probability (eqn (7)), which requires the adiabatic NACEs (eqn (8)) and in particular the second term Â U y _ U Ã ja . As suggested by Hammes-Schiffer and Tully, 67 we take advantage of the anti-symmetry of this term. After mapping, this term now is, Finally, we apply the self-consistent correction to improve the probability to hop (as in eqn (15)) towards the closest state in energy, i.e., the one likely to be affected by numerical inaccuracies due to finite timesteps. The overall efficiency of our approach to remove the trivial crossings will be discussed in Section 4.2.

Spurious charge transfer
Although the decoherence correction is paramount to maintain the internal consistency of the FSSH algorithm as discussed in Section 2.2, it may create an undesirable charge displacement in the cases of large systems and high densities of states.
In those systems, a surface hop between adiabatic states localized in different regions of space is unlikely but still possible in the FSSH algorithm (due to the stochasticity of the hopping algorithm). After an unlikely hop, the decoherence correction scheme (damping-based or collapsing-based) will move the electronic wavefunction C(t) closer to the new active state c a (t), thus leading to unphysical charge transfer, as illustrated in Fig. 1(C). This problem was recently pointed out by our group 26 and Wang and coworkers. 58 We labelled such events decoherence correction-induced spurious charge transfers (DCICTs). They arise because the different decoherence correction schemes act on the adiabatic (i.e., non-local) representation. Note that, since the mean-squared displacement (MSD) depends on the square of the distance, this will have a strong impact on the diffusion coefficient and on the mobility (see Section 2.7).
Recently, Wang and collaborators proposed to switch off the decoherence correction when the surface population is below a certain threshold and showed that the spurious transfer is indeed alleviated. 58 However, this formulation can reduce the internal consistency of surface hopping as some decoherence events are actually removed. Independently from the latter study, we have developed a three-step strategy to remove the DCICTs as illustrated in Fig. 1(C): (i) at each timestep, an ''active'' region that encloses 99.9% of the electronic density |C(t)| 2 is determined, (ii) the decoherence correction is applied and (iii) any change in the diabatic population D|u l | 2 outside the active region is reset to zero, while the diabatic populations inside the active region are scaled accordingly to preserve the norm. We call this strategy spurious charge transfer correction (SCTC, previously termed SPTC in our previous paper 26 ). In practice, it amounts to a local decoherence correction within the active region, while outside the active region the diabatic populations remain unchanged. All DCICTs are removed, while decoherence correction is still applied at each timestep. Note that the propagation of the wave function according to eqn (4) remains unaffected by the presence of the active region.

Elimination of forbidden hops
As discussed in Section 2.1, a key issue in FSSH simulation is the internal inconsistency between the surface population and the wavefunction population. One of the main reasons for this discrepancy is the lack of an intrinsic decoherence mechanism, as detailed in Section 2.2. Another source of internal inconsistency is attributed to forbidden hops, namely attempted hops triggered by the stochastic FSSH algorithm (probability to hop in eqn (7)) but that do not fulfil the conservation energy requirement (eqn (9)). 73 Although frustrated hops are actually essential to maintain detailed balance as pointed out by Tully, 29,74 they can be a source of internal inconsistency. When such forbidden hops are encountered, NACEs between the active state and those classically forbidden states are finite and the TDSE (eqn (4)) will transfer a small amount of the electronic wavefunction C(t) to those states. In some cases an ad hoc decoherence scheme cannot completely correct the problem especially when the electronic coupling is high (i.e., when a large number of forbidden transitions are present because the surfaces are distant from each other). This problem has been partially overlooked in the literature as most applications of surface hopping are related to fast relaxation processes from excited states to low-energy states. 44,48,50 In this case, almost all transitions are downward (thus allowed) and a simple decoherence scheme is sufficient to reach a satisfactory internal consistency.
At the early stage of FSSH development, Hammes-Schiffer and co-workers already suggested a route to ensure internal consistency even in the presence of a great number of forbidden transitions. 44,73 They proposed to remove the undesired population transfers, by setting the corresponding NACEs/NACVs to classically forbidden states to zero in the TDSE. In this work, we call their approach elimination of forbidden hop (EFH). In this approach the FSSH algorithm is modified in the following way: (i) at each timestep, one first determines which transitions are energetically forbidden between the active state and the excited states using the energy conservation criterion in eqn (9), (ii) for such forbidden transitions, the corresponding adiabatic NACEs/ NACVs are exactly set to zero and (iii) one propagates the modified electronic TDSE (analogue to eqn (4)) that, in the adiabatic basis, now reads, in which d ad,efh ij are the modified adiabatic NACEs forming the new matrix D ad,efh . This matrix is now sparse as certain transitions are forbidden. In our implementation, the electronic propagation (eqn (4)) is carried out in the diabatic basis, which gives better numerical stability. 25,75 For this reason, eqn (18) can be transformed in the diabatic basis: in which d efh kl = [D efh ] kl and D efh is the NACE matrix in the diabatic basis. The latter can be written in terms of the adiabatic NACE matrix D ad,efh as (see eqn (8)), The main difficulty arising in eqn (20) is the presence of the time derivative _ U and the adiabatic NACE between two subsequent nuclear time steps. As discussed in Section 2.3, especially near crossing points, NACEs have sharp localized peaks resembling Dirac delta functions that can be easily missed. However, as we show in the following, _ U can be eliminated from eqn (20). In its final form D efh only contains the smooth NACEs between the active state and the energetically forbidden states that are high in energy. To this end, we write D efh in the following way: where D is the matrix of diabatic NACEs with elements that appear in the unmodified Schrödinger, eqn (4). Substituting eqn (20) in eqn (21) and defining D ad ¼ D ad;efh À D ad , _ U is eliminated from eqn (20), We will discuss in Section 4.1 the effects of this alternative electronic propagation on the internal consistency and the equilibrium properties.

Electronic populations
Our objective is to calculate the electronic mobility from FOB-SH simulations. Electronic-based properties are, however, ambiguous in FSSH, as different definitions can be found for the electronic population. Landry and Subotnik provided a detailed account of the existing definitions in ref. 59 and they highlighted that these populations can produce divergent properties. In the most common approach (''Method 1'' in ref. 59 or the ''surface method'') electronic properties are calculated using the active adiabatic state c a (t). To avoid confusion, we prefer to call this population a projected active state population (PAS). The electronic population on site k is Other authors use the intrinsic FSSH wavefunction C(t) (''Method 2'' or the ''wavefunction population'') to obtain electronic properties and the local population on k is 26,76,77 This population definition relies on the propagated electronic wavefunction C(t) that does not observe detailed balance in the absence of decoherence. We will demonstrate in Section 4.2 that this definition fulfils detailed balance as well as Method 1 when decoherence is included. Finally, Landry and Subotnik also suggested using the mixed quantum-classical density 36 (''Method 3'' or the ''MQC population'') and they obtained the following diabatic population: As mentioned before, the projected active state method has the advantage of giving the correct detailed balance distribution (the electronic state distribution follows an approximately Boltzmann population in FSSH). However, this method is also more sensitive to trivial crossings, as any trivial crossing missed will instantaneously modify c a (t). In contrast, C(t) (in ''Method 2'') would not be directly impacted by a missed trivial crossing, although in the long term there will be a bias in the dynamics.

Mean-square displacement
In the literature, the charge mobility m is often obtained from the diffusion coefficient D by means of the Einstein equation: where k B is the Boltzmann constant, T the temperature and q the charge of the carrier. The diffusion coefficient is defined as the slope of the mean square displacement (MSD) at long times, While for a classical particle the MSD is well defined, MSD = hx(t) À x(0) 2 i, with x the position of the classical particle and the brackets referring to a statistical average, various equations for the MSD of a quantum particle can be found in the literature. In our previous papers, 26,60 we propose first to determine the expectation value of the ''position'' of the quantum particle, % x(t) = hC(t)|x|C(t)i, and then to use it within the classical definition to obtain: in which P k,n is the diabatic population of site k and for trajectory n (P k,n could be P PAS k,n , P wf k,n or P MQC k,n ), and x k,n is the distance between the center of mass of molecule k and the molecule initially charged. Finally, the sum runs over trajectories and stands for the statistical average. We denote the resultant MSD in eqn (28) as MSD coc , where coc stands for the center of charge. On the other hand, Elstner and collaborators 76,77 have chosen to define the MSD as the expectation value of the squared displacements (X À % x(0)) 2 (with % x(0) = hC(0)|x|C(0)i) and to average the expectation value of such an operator in the following manner: Other authors 55,58,78 considered the time-evolution of the spread or variance (var) of the wavefunction, The three definitions are of course related, MSD = MSD coc + MSD var . The MSD (eqn (30)) is preferable because it accounts for both the motion of the center of charge and the spreading of the wavefunction. Therefore it can be used for both extremes, small polaron hopping and pure wavefunction spreading, and all intermediate cases. We will analyse the two different contributions MSD coc and MSD var in Section 4.2.
3 Numerical details

Chains of ethylene-like molecules
To investigate the detailed balance of FOB-SH, we used onedimensional chains of two (dimer), three (trimer) and five (pentamer) ethylene-like molecules (ELMs), as shown in Fig. 3(A). The name ''ethylene-like'' stresses that only the nuclear geometries correspond to a real ethylene molecule, while the CT parameters (namely the reorganization energy l and the AOM scaling value C) are chosen freely to explore a large range of physical behaviours. Within a chain, the ELMs are spaced by 4 Å and a weak center of mass restraint potential (force constant = 11 kcal mol À1 Å À2 ) is applied to keep the chain straight. Because the NVE ensemble is used in all our FOB-SH simulations to avoid any dynamical bias introduced by the thermostat, the chain is embedded in a bath of neon atoms that mimics the role of a thermostat and reduces the fluctuations of the temperature (especially for the short chains with few degrees of freedom). The simulation boxes are cubic with size a = 60 Å and contain one chain of ELMs and 123 neon atoms (124 for the dimer). Periodic boundary conditions are applied in all directions of the simulation box. We insist, however, that the electronic propagation occurring within the chain of ELMs is not periodic: when the charge reaches the edges of the chain, it is scattered backward and it does not continue at the other end of the chain.
In the present model, hole transfer is mediated by a set of (orthogonalized) HOMOs of the ethylene molecules, f k , k = 1, M, that are used to construct the electronic Hamiltonian H.
Diagonalization of the Hamiltonian H gives the M adiabatic electronic states. For a detailed explanation of how the orbitals f i (R(t)) are reconstructed along the trajectory we refer to ref. 60. The diagonal elements H kk are calculated with a force field energy function whose parameters for neutral and positively charged ELMs are chosen as in our previous work. 25,60 For charged ELMs, the equilibrium distance of the CQC bond is displaced (1.387 Å) with respect to the one in the neutral state (1.324 Å) corresponding to the reorganization energy for hole transfer between two ELMs of l = 200 meV. Such reorganization energies are typical for organic semiconductors and an order of magnitude smaller than those, e.g., for redox processes in aqueous solution [79][80][81] or oxide materials. 82 Intra-molecular interactions for neutral ELMs are taken from the Generalized Amber Force Field (GAFF). 83 The intermolecular interactions among the ELMs and between ELMs and Ne atoms are modelled by Lennard-Jones terms with parameters taken again from the GAFF database for neutral and charged ELMs and from ref. 84 for Ne and applying the Lorentz-Berthelot mixing rules. Electrostatic interactions in the form of fixed point charges do not significantly alter the energetics of this system because only one ELM carries a net charge and the other ELMs and Ne are charge neutral. Hence, for convenience, electrostatic interactions were switched off in all simulations.
The initial configurations are built with the investigated chain in its energy-minimized geometry and the neon atoms positioned on a regular grid. The system is equilibrated with a 1 ns NVT run at 298 K using a Nosé-Hoover thermostat 85,86 and using a force field energy function where the first molecule of the chain is charged. From the last configuration of the NVT run, 100 ps Born-Oppenheimer molecular dynamics (BOMD) trajectories are started for each adiabatic electronic state (two, three and five states for the dimer, trimer and pentamer, respectively). This is done for each of the six AOM scaling values C that determine the strength of electronic coupling according to H kl = C % We extract Boltzmann-weighted configurations (nuclear coordinates and velocities) from such BOMD runs as starting configurations for the FOB-SH runs, to ensure the correct distribution of excited state populations at the start of the run (at t = 0) and a well-sampled phase space. The electronic wavefunction is initialized in the corresponding adiabatic state i (C(0) = c i ) to ensure perfect internal consistency at t = 0. For each set of parameters (chain length, C value, decoherence correction), we generated 1000 independent FOB-SH trajectories starting from the initial configuration evenly sampled from the corresponding BOMD trajectories. Each trajectory is run for 10 ps in the NVE ensemble. The nuclear dynamics is propagated with the velocity-Verlet algorithm with forces calculated according to eqn (6) and with a MD timestep Dt = 0.1 fs. The wavefunction of the excess charge carrier C(t) was propagated by integrating eqn (4) using the Runge-Kutta algorithm of 4th order and an electronic timestep dt = Dt/5 = 0.02 fs. An interpolation scheme is used to calculate the Hamiltonian matrix elements at each electronic timestep. 60 Error bars were determined by block averaging over the 1000 trajectories with a block size of 200 independent runs.

Chains of embedded anthracenes
To investigate the mobility and IPR of a real system, we have modelled an electron hole transfer in a chain of electronically active anthracene molecules, embedded in a larger crystal comprised of electronically inactive anthracene molecules, see Fig. 3(B). We have compared 4 chains of different lengths: 12, 24, 36 and 48 molecules. The simulation boxes are monoclinic, with angles a = 90.01, b = 124.71 and g = 90.01 and with dimensions a = 8.562 Å and c = 11.184 Å. We have adapted the box length in the b direction to ensure that the distance between the chain and the edges of the box is above 8 Å. Periodic boundary conditions are applied in all directions of the crystal, but, similarly to the chains of ELMs, the electronic dynamics occurs along non-periodic chains. To check the convergence of mobility with respect to the length of the chain, different numbers of molecules of chains are investigated. Table 1 presents the total number of molecules and the length b for the different chain lengths.
As for the ELM model system described in the previous section, we assume that the electron hole transfer is mediated by the HOMOs of the anthracene molecules that form the basis functions for the excess charge expansion (eqn (1)). The M diagonal elements H kk are, again, estimated using M classical force field energy functions. In the kth energy functions, anthracene molecule k is positively charged, while all the others are neutral. Intra-molecular interactions for the neutral anthracene molecule are taken from the Generalized Amber Force Field (GAFF). 83 These intramolecular parameters are used also for the charged anthracene, except for the carbon-carbon bond length which was chosen instead to reproduce the reorganization energy l. The reorganization energy is determined using four DFT calculations on neutral and charge anthracene molecules in both neutral and charged geometries as: where E C/N (R N/C ) is the energy of the charged/neutral molecule in the optimized neutral/charged state and E C/N (R C/N ) is the energy of the charged/neutral molecule in the optimized charged/neutral minimum. The geometries of charged and neutral molecules were optimized with the B3LYP functional 87 and 6-311g(d) basis set. The intermolecular interactions between anthracene molecules are also taken from the GAFF database. As in each classical force field all but one molecule are neutral and the anthracene has zero dipole moment, we did not include electrostatic interactions. The off-diagonal elements of the electronic Hamiltonian H kl are calculated using the AOM. 61 First, the HOMO of anthracene (which is non-degenerate) is projected onto an atomic Slater basis consisting of one atomic p orbital per carbon atom. The calculation of the HOMO and its projection are done using CPMD software 88 using the PBE exchange-correlation functional. 89 Core electrons are described by Goedecker-Teter-Hutter (GTH) pseudo-potentials, 90 and the valence electrons are expanded in plane waves with a reciprocal space plane wave cutoff of 90 Ry. The dimers are centered in a simulation box with dimensions of 12 Â 40 Â 40 Å 3 . After that, the electronic coupling H DFT kl is calculated using the FODFT method 62 for four different dimers extracted from the crystal structure, while the HOMO-HOMO overlap % S kl = hj k |j l i is calculated using the AOM for the same four dimers. The FODFT couplings are scaled by a constant 1.348 as recommended in ref. 62. A linear regression is applied between H DFT kl and % S kl to determine the AOM scaling value C = 3.09 eV. Each FOB-SH simulation involves 1000 independent trajectories initialized from 100 different initial conditions (10 trajectories repeated with a different random seed for each initial condition). Starting from the crystal structure, the system is equilibrated for 500 ps in the NVT ensemble using a Nosé-Hoover thermostat. 85,86 Then a MD run of length 500 ps is carried out in the NVE ensemble from which 100 configurations are chosen at equidistant intervals. These configurations are used as the initial configurations for subsequent FOB-SH runs. The initial wavefunction is fully localized on the first molecule of the chain, C(t = 0) = f 1 (0), and the initial active state is randomly drawn from all adiabatic states with a probability hc i (0)|f 1 (0)i 2 . Each trajectory is then run for 2 ps in the NVE ensemble. We opt for the NVE ensemble to avoid any artificial thermostat that may bias the calculation of the electronic mobility. The large number of degrees of freedom due to the ''inactive'' part (inactive for electronic propagation, as depicted in Fig. 3) of the anthracene crystal plays the role of a thermostat and ensures small temperature fluctuations. The nuclear dynamics

Energy conservation, detailed balance and internal consistency
Before investigating charge transport properties, we focus here on the influence of various decoherence correction schemes and the treatment of classically forbidden hops on key equilibrium properties such as energy conservation, detailed balance and internal consistency. In this section, we exclude analysis of trivial crossings or DCICTs, as these issues arise only when considering large systems and transport properties. 4.1.1 Effect of decoherence correction. The role of decoherence correction in detailed balance has been only partially considered in the literature. 31 To investigate to what extent decoherence correction influences the thermal population of each state, we have carried out FOB-SH simulations using one dimer of ELMs in a bath of neon atoms for several commonly used decoherence correction algorithms: instantaneous decoherence after each attempted hop (IDA), 50 energy-based decoherence correction (EDC, eqn (11)), 47 pure dephasing decoherence correction (PDDC, eqn (12)), force-based decoherence correction 41,43,69 using a damping algorithm (FDC, eqn (13)) and a stochastic algorithm (SC-FDC), and finally the absence of a correction scheme (NO DC). Simulations are initialized as described in Section 3.1. Fig. 4(A) shows the energy drift averaged over 1000 FOB-SH NVE runs as a function of electronic coupling. The general trend is similar to the results previously obtained by us in ref. 25, with a monotonic decrease of the energy drift from 10 À5 Ha per ps per QM atom to 10 À7 Ha per ps per QM atom. We explain this behaviour by observing that, with increasing coupling, the number of successful hops decreases, while the potential energy surface softens. The notable fact is that the energy drift is independent of the decoherence correction scheme. This can be expected as the decoherence only affects the electronic wavefunction, not directly the nuclear degrees of freedom whose total energy is conserved along the simulation.
A similar conclusion holds for the detailed balance. In Fig. 4(B), we show the electronic population of the excited state, averaged over the 1000 trajectories and over time, against the time average electronic coupling. The exact result obtained from the BOMD simulations as described in Section 3.1 is also indicated. Since the work of Tully and collaborators, 29,74 the ability of the ''vanilla'' FSSH (i.e., without decoherence correction) to reach detailed balance is well-known. We recently reinforced the point that the NACV-oriented adjustment of velocities after a hop is paramount for this agreement to hold. 25 Remarkably, we find here that the bias introduced by the decoherence correction in the electronic dynamics is almost negligible in terms of equilibrium distribution. This can be readily explained in the cases of EDC, PDDC and FDC, for which the decoherence time is small (i.e., fast decoherence) far from the crossing region and it is large (i.e., slow decoherence) within the crossing region. For this reason, such corrections have only a minor effect in the proximity of an avoided crossing, which is where the probability for hops sharply increases and the thermal equilibration between the electronic states occurs. Thus, decoherence only affects the dynamics away from the crossing region, where, in any case, the surfaces are quite well separated in energy and the number of hops is small. Therefore, dampingbased schemes maintain the correct flux between states and do not ruin the detailed balance.
It is important to notice that such an argument does not apply to instantaneous decoherence algorithms. These algorithms require the nuclei to be in the crossing region in order to trigger the decoherence event (i.e. there must be either an attempted or a successful hop in order to collapse the wavefunction) and they do not depend on any decoherence time. This explains why for the latter algorithms we can observe larger deviations for both energy drift and excited state population, even though the bias is still small due to the small number of collapsing events with respect to the total number of steps in the dynamics.
We conclude that all the decoherence schemes investigated here can reach approximately the detailed balance, meaning that the bias introduced in the electronic dynamics does not affect the flux between adiabatic states.
While the different decoherence algorithms give virtually identical results for energy drift and detailed balance, they give  (13)), stochastic force-based (SC-FDC), energybased (EDC, eqn (11)), and pure dephasing decoherence times (PDDC, eqn (12)), instant collapse (IDA) and no DC. Exact populations in (B) are obtained as described in Section 3.1. The internal consistency in (C) is measured in terms of the root-mean-square error (RMSE, eqn (34)), and divided by the excited state population P ex 1 in (D).
very different results for internal consistency. We measure the latter by calculating the time-averaged root mean square error between the surface population and the quantum amplitude of the excited state i, P surf i (t) (eqn (10)) and h|c i (t)| 2 i trj , respectively, where hÁ Á Ái trj refers to an average over trajectories. In case of perfect internal consistency, RMSE i = 0 for all i. Fig. 4(C) shows the RMSEs obtained for the usual range of coupling values and Fig. 4(D) shows the RMSEs normalised with respect to the exact excited state population, RMSE i /P surf i . We can observe that RMSEs follow the same trend for all the decoherence methods (an increase from low couplings to medium-sized coupling values (maximum around 50 meV) and a slight decrease thereafter). The damping methods show very good internal consistency, with FDC and DCCP giving the best performance for all coupling strengths, and significantly improving over wavefunction collapse and no decoherence correction ( Fig. 4(C)). Hence, the particular choice of damping time seems rather unimportant for good average internal consistency, i.e. in the long time limit (eqn (34)).
However, Fig. 4(D) reveals that the internal consistency, normalized with respect to the excited state population P ex i , deteriorates with increasing coupling strength. The quantum populations of excited states are generally overestimated in this regime. In our previous paper, 25 we showed that for couplings V 4 k B T/2, adiabatic NACEs still transfer the electronic population from the ground state to the excited state, while attempted hops become increasingly energy-forbidden. Therefore, the wavefunction population in the excited state is overestimated compared to the surface population. While for no DC and collapse the error is substantial to the extent that there is no longer any consistency between quantum and surface amplitudes in the high coupling regime, the damping methods significantly improve on this situation, albeit not perfectly. For a large coupling value of 100 meV, the excited state surface population is about 10 À4 (Fig. 4(B)), while the quantum populations are about 10 À3 , giving RMSE 1 /P ex 1 E 10. While this deviation may not be relevant in many practical situations, it is desirable to investigate further possible improvements to internal consistency such as elimination of energy forbidden hops.

EFH propagation.
To improve the internal consistency at medium and large couplings, we resorted to the EFH scheme suggested by Hammes-Schiffer and collaborators 44,73 and described in detail in Section 2.5 for a diabatic electronic propagation. In EFH the electronic propagation is modified, avoiding electronic population transfer to excited states energetically unreachable by the nuclei (i.e., these states would fail the energy conservation requirement (eqn (9)) with a NACV-oriented adjustment). To test the EFH algorithm, we have carried out FOB-SH simulations using an identical set-up as described before (a dimer of ELMs in a bath of neon atoms), but now the electronic dynamics is propagated using eqn (18). Fig. 5(A) and (B) show, respectively, the excited state populations and RMSE 1 /P ex 1 as functions of electronic couplings for the EFH propagation with decoherence, standard FSSH propagation with a decoherence scheme and standard FSSH propagation without any decoherence scheme. For both EFH and standard FSSH propagation, we use force-based decoherence correction (FDC) as we have shown that all damping-based approaches give identical results for detailed balance and internal consistency. Though EFH propagation biases electronic dynamics and, indirectly, hopping probability and detailed balance, Fig. 5(A) shows that EFH dynamics gives excited state populations in agreement with the exact results. In fact, at small couplings, when the number of frustrated hops is small, EFH yields the same results as the standard FSSH propagation with decoherence. By contrast, when the electronic coupling is high and the effect of EFH becomes active (removal of amplitude transfer to states that are energetically not accessible), EFH propagation shows a larger uncertainty and a larger error bar. This effect will be considered further in the following. Fig. 5(B) shows again that EFH coincides with the standard propagation at small couplings where the internal consistency is good, but performs by far better at larger couplings with respect to standard FSSH propagation with and without decoherence (i.e., RMSE i /P ex i is about two orders of magnitude smaller in the former case).

Trimer and pentamer chains.
To assess the validity of our previous conclusions, we investigate detailed balance and internal consistency for two larger systems: a trimer and a pentamer of ELMs with three and five states, respectively, embedded in a bath of neon atoms. The initialization of these systems is described in Section 3.1. For each system, we have carried out FOB-SH simulations using two different set-ups: standard propagation with FDC to account for decoherence and EFH propagation with FDC. Fig. 6(A) and (C) show the populations P surf i (t) for the different excited states i and Fig. 6(B) and (D) show RMSE i /P ex i against electronic coupling. The standard FSSH propagation produces excited state populations close to the exact ones for all excited states in both systems. This confirms that FSSH with a decoherence scheme can reach detailed balance  (13)). Internal consistency is measured by the relative root-mean-square error defined in eqn (34). The exact result for the excited state population is obtained as described in Section 3.1. Error bars are obtained by block-averaging over five independent blocks of 200 trajectories each. even for larger systems. On the other hand, excited state populations using EFH differ quantitatively from the exact results. We suggest the following explanation. In a situation where the adiabatic NACE is not-negligible between the active state and another state much higher in energy, there can be population transfer to the excited state due to the standard FSSH propagation. If a hop is attempted, it will be rejected due to the large energetic gap between those states. However, the population transfer may increase the probability of hopping a few time steps later, when hypothetically the energy gap might be smaller and therefore this second hop successful. By contrast, using EFH propagation no population transfer is allowed in the first case, leading to a decrease in the hopping probability also in the second case where the energy gap between surfaces would have been small enough to allow the second attempted hop. On average, EFH diminishes the flux to the excited states and thus underestimates their populations. Even though detailed balance in EFH is not as good as in the standard propagation scheme, it massively improves internal consistency for the excited states of both systems in the high coupling region, as can be seen in Fig. 6(B) and (D).
In conclusion, EFH greatly improves the internal consistency; nevertheless, it biases the hopping probability and produces a worse agreement with the exact equilibrium population when compared with the standard FSSH propagation. It is worth noting that in large organic semiconductors the density of states in a given band is quite high, most of the hops are allowed and the most important source of internal inconsistency is the wavefunction branching rather than the presence of frustrated hops. For these reasons we will not consider further this correction to investigate dynamical properties and charge transport.

Charge mobility
We now focus on building and describing the best FOB-SH set-up to calculate the mobility and delocalization of a charge carrier in organic semiconductors. We make use of the system described in Section 3.2, a hole diffusing in a chain of anthracene molecules embedded in a larger crystal. Several parameters and set-ups need to be assessed: the role of the state tracking algorithm and/or the spurious transfer correction, the size of the system, the MD timestep, the decoherence correction scheme, the population definition, the MSD definition and finally the length and the number of trajectories.
4.2.1 MSD definition. We first clarify how to determine the mobility. As we have indicated in Section 2.7, different definitions of the MSD have been used in the literature. We ran a FOB-SH simulation for a chain of 48 anthracene molecules for 2 ps using our reference set-up: state tracking and SCTC are switched on, the MD timestep is 0.1 fs, the decoherence correction scheme is pure-dephasing (PDDC) and the diabatic population P wf is used to calculate the MSD. Fig. 7(A) shows the three MSD definitions against time. MSD coc and MSD quickly increase simultaneously until 200 fs, then they diverge slightly and end up in a linear regime with the same slope, clearly visible after 1 ps. In contrast, MSD var increases more slowly and reaches a plateau regime after 1 ps. In fact, as shown by different authors using simulation tools 26,55,58,91 and by experiments, 92-94 the electronic wavefunction is delocalized over few molecules and the charge spreading remains approximately constant in time.
As described in Section 2.7, the mobility is related via the diffusion coefficient (eqn (26)) to the slope of the MSD at long times (eqn (27)). In Fig. 7(A), the best linear fits are indicated by black dashed lines for all three MSD definitions. We conclude that to determine the mobility, both MSD coc and MSD will give the same value for the diffusion constant, whereas MSD var will give a zero value for this coefficient and so for the mobility.  Besides the mobility, it is also interesting to measure the delocalization of the wavefunction. Rather than using the wavefunction spreading, MSD var , we prefer to follow ref. 55 and to calculate the inverse participation ratio (IPR): In particular, IPR n (t) measures the number of sites over which the wavefunction is delocalized at time t for trajectory n. This quantity can be averaged over time and over trajectories to obtain the equilibrium converged value. In Fig. 7(B) we show the evolution of the IPR against time. At t = 0, the wavefunction is initially localized on the first molecule of the chain and IPR(t = 0) = 1. The IPR increases rapidly during the first few hundreds of femtoseconds before reaching a plateau at around 1 ps. The initial transient increase corresponds to the wavefunction spreading until the equilibrium polaron width is reached. In fact, this behaviour of the IPR mirrors exactly the time evolution of MSD var . We note that in this system the transient behaviour lasts about 1 ps and that trajectories of at least 2 ps are necessary to calculate the mobility. In the following, we always apply a linear fit to the MSD between 1 ps and 2 ps to extract the slope, the diffusion coefficient and the mobility and we calculate the time average of the IPR also between 1 and 2 ps. 4.2.2 Trivial crossings and spurious charge transfer. We now investigate the necessity of the state tracking algorithm and decoherence-induced spurious charge transfer correction described, respectively, in Sections 2.2 and 2.3 in order to obtain convergence for the MSD and the mobility as functions of chain length and timestep. Fig. 8 shows the MSDs for different chain lengths and the wavefunction populations of sites k (P wf k , eqn (24)) for a representative FOB-SH trajectory for three different set-ups: (A) state tracking and SCTC both active, (B) state tracking and no SCTC and (C) no state tracking and no SCTC. If both state tracking and SCTC are switched on, the MSD converges for increasing chain length: up to 500 fs all MSDs are identical for the different chain lengths, as the charge initially explores just a few molecules around its initial position. After 500 fs, the MSD for a chain of 12 molecules is slightly below the MSD for the longer chains, which are all practically identical, as the diffusive charge starts to feel the boundary and to bounce back at the end of the chain. The smooth evolution of the polaronic charge carrier is illustrated for one representative trajectory in Fig. 8(D), where no spurious charge transfer event is present. The spreading of the charge carrier is around 2, in agreement with the IPR in Fig. 7(B). When the SCTC is switched off (panel B), MSDs as function of chain lengths start to diverge after a few femtoseconds, indicating that frequent decoherenceinduced spurious transfer events bias the charge dynamics. It is worth noticing that spurious charge transfers induce a much larger displacement of the charge as in a few time steps the charge can completely change its localization. A spurious transfer is shown in panel (E) for a representative FOB-SH trajectory without SCTC, where the charge carrier ''jumps'' from molecule 7 to molecule 13 in a few femtoseconds. We also note that the order of magnitudes spanned by these MSDs (500-100 Angstroms) and the presence of a linear regime may be deceptive, but the divergence with system size underlines the unphysical aspect of the charge displacement. In Fig. 8(C), we show the MSDs for different chain lengths when both SCTC and state tracking are switched off. In a few femtoseconds, the MSDs reach a plateau that depends on the size of the system (i.e., the larger the system the larger the plateau value). Such dynamics for the charge corresponds to an unphysically fast diffusion in which the numerous missed trivial crossings yield an almost random motion of the charge along the chain. This is well exemplified in Fig. 8(F) for a FOB-SH trajectory without SCTC and state reordering. Missing index updates cause numerous jumps of the charge carrier at long distance (tens of molecules). Only the use of a state-tracking algorithm to detect the trivial crossings and the SCTC to eliminate the decoherence-induced spurious charge transfers leads to a physical MSD independent of system size.

Number of trajectories.
We now consider the convergence of transport properties (mobility m and IPR) as a function of the number of FOB-SH trajectories. We ran FOB-SH simulations using different numbers of trajectories (10,50,100,150,200, 500, 1000) for a chain of 48 anthracene molecules. Fig. 9(A) and (B) show the mobility and the IPR calculated using these different FOB-SH simulations against the number of trajectories. We find that the mobility converges at about For all set-ups, decoherence keeps the charge localized over about 2 molecules (consistently with Fig. 7(B)). Long-range spurious transfer events are highlighted with red arrows in (E) and (F); note that the charge transport in (C) is completely biased by unphysical jumps of the charge. The MD timestep is 0.1 fs and the decoherence correction is damping with pure dephasing decoherence time. Error bars represent standard deviations over five independent blocks of 200 trajectories each. 100 trajectories and the IPR after as few as 10 trajectories. This implies that the IPR distribution is more homogeneous among the trajectories than the mobility one. Based on the convergence of mobility with 100 trajectories, we calculated error bars for mobility and the IPR in the FOB-SH simulation by blockaveraging over five independent blocks of 200 trajectories each.

Timestep.
To confirm that our set-up (decoherence correction, state-tracking and SCTC) also permits the convergence with timestep, we ran a FOB-SH trajectory using four different timesteps (0.025 fs, 0.05 fs, 0.1 fs and 0.5 fs) for a chain of 48 anthracene molecules. We report the mobility and IPR values obtained using these different timesteps in Fig. 9(C) and (D) respectively. Fig. 9(C) and (D) show that the three smallest timesteps (0.025 fs, 0.05 fs, 0.1 fs) give the same mobility results, whereas the mobility is four times higher at 0.5 fs, indicating that the state tracking algorithm fails to detect all trivial crossings at this large timestep. In contrast, the different IPRs presented in Fig. 9(D) are identical for all MD timesteps. The IPR is a static quantity and it is thus insensitive to the failure in the detection of the trivial crossings.
4.2.5 Decoherence and diabatic populations. We now discuss the choice of diabatic population definition and the effect of the decoherence correction. Using the converged set-up described in the paragraph above (state tracking, SCTC and 0.1 fs for the MD timestep), we ran FOB-SH simulations for different chain lengths (12, 24, 36 and 48) either without decoherence or using the pure-dephasing decoherence schemes. For each FOB-SH run, we calculated the mobility and IPR for the three different diabatic population definitions (P wf , P PAS , P MQC , see Section 2.6). The results are reported in Fig. 10 against the number of molecules forming the chain. In the original FSSH implementation, 39 no decoherence correction and P PAS are used to calculate different properties; this set-up corresponds to the green lines in Fig. 10(A) and (B). We immediately see that the mobility does not converge with increasing system size. Without decoherence, in fact, Tully and collaborators 29 already showed that in the long time limit the electronic wavefunction C(t) delocalizes equally over all the available adiabatic states. In that case, the adiabatic population appearing in the denominator of the hopping probability (eqn (7)) is the same for all states, so hops can be attempted between states localized in completely different positions. The larger the chain, the farther apart the charge can jump after such unphysical (but allowed by FSSH) hops. Thus the mobility for all three different diabatic population definitions increases with the number of molecules. Conversely, the IPR for P PAS is independent of system size, showing that the delocalization of the eigenfunction of the active state (c a ) is size independent as well. In contrast, the IPR for P wf increases with chain length, mirroring the delocalization of the electronic wavefunction on the adiabatic states. We note that the IPR for P MQC is very close to the one for P PAS . These results prove that, without decoherence,  (12)). In (C and D), error bars are obtained by block-averaging over five independent blocks of 200 trajectories each. Fig. 10 Importance of decoherence correction for the convergence of charge mobility with respect to system size. No decoherence correction is applied in (A) and (B) and the pure dephasing decoherence correction (PDDC, eqn (12)) scheme is applied in (C) and (D). Results are shown for different choices for the diabatic populations used to calculate the charge mobility m and the inverse participation ratio (IPR): wavefunction (P wf , eqn (24)), active state (P PAS , eqn (23)) and mixed quantum-classical populations (P MQC , eqn (25)). The data were obtained from FOB-SH simulations of hole transport along an embedded chain of anthracene molecules with a MD timestep of 0.1 fs. Error bars are obtained by blockaveraging over five independent blocks of 200 trajectories each. the mobility cannot converge with different system sizes whatever the diabatic population definition.
The mobility and IPR obtained with a decoherence correction (pure-dephasing correction) are shown in Fig. 10(C) and (D) respectively. In contrast to the results obtained without decoherence, the mobility is well converged with respect to chain length. Adding a decoherence correction permits localizing (in adiabatic and diabatic space) the electronic wavefunction C(t), to eliminate the undesired hops present without decoherence and to converge with system size. The decoherence correction also ensures the internal consistency of the method, explaining why the three population definitions behave similarly. The IPR results are similar to the ones for the mobility: convergence for the different system sizes and similar values for all three population definitions. Although, in general, we recommend using the wavefunction population (P wf , eqn (24)) as it is generally less affected by potentially undetected trivial crossings. Based on these results, we conclude that a decoherence correction is mandatory for the calculation of mobility and the IPR.

Conclusion
In the present work, we have employed the fragment-orbital based surface hopping (FOB-SH) approach that is a powerful tool to perform atomistic non-adiabatic dynamics in large realistic molecular systems. We have explored and discussed several possible improvements applicable to any surface hopping code when calculating equilibrium and dynamical properties, i.e., decoherence correction with various decoherence times, spurious charge transfer correction (SCTC), electronic propagation with the elimination of classically forbidden transitions (EFH), trivial crossing correction and state tracking. The correct way to retrieve important observables from FSSH simulations such as electronic populations, mean-square displacement and electronic mobility has also been discussed.
Using an organic semiconductor model formed by chains of ethylene-like molecules, we have first looked at equilibrium properties (i.e., energy conservation, detailed balance and internal consistency) over three orders of magnitude of electronic coupling and different system sizes. We have shown that good energy conservation and detailed balance is obtained regardless of the decoherence time and algorithm used. In fact, generally speaking, the decoherence biases the dynamics only away from the crossing region and it does not significantly modify the flux between adiabatic states. On the other hand, when comparing the effects of different decoherence corrections in restoring the consistency between surface and wavefunction populations, we have shown that the damping-based algorithms with fast decoherence times produce far better results than instantaneous collapsing events and maximize internal consistency across several orders of coupling strengths. However, for small systems, some degrees of internal inconsistency can still be seen when electronic coupling is large due to the presence of a significant number of forbidden transitions and the inability of the ad hoc damping procedure to correct for such cases. Therefore, to further reduce the internal inconsistency in this coupling region, we have explored an alternative electronic propagation (called the elimination of classically forbidden hops -EFH) in which electronic population transfer between adiabatic states is removed in the case of classically forbidden transitions for the nuclei. Although the algorithm massively improves internal consistency at high couplings, the agreement with the detailed balance deteriorates due to the bias introduced in the electronic dynamics.
Then, focussing on charge transport in a real organic crystal (i.e., anthracene), we have studied two fundamental properties related to the actual efficiency of organic semiconductors: the electronic mobility and the inverse participation ratio (the latter measures the size of the charge carrier). We have found that charge carriers propagate through OSs as polarons via diffusive jumps, somewhat in analogy with the diffusion of gas molecules in a complex environment, 95,96 though with the sizes and shapes of the polarons strongly fluctuating in time. We have also found that a state-tracking algorithm is mandatory in the case of a large number of states to detect the trivial crossings and to map the adiabatic states between two different MD timesteps, thus improving the electronic and nuclear dynamics and avoiding spurious long-range charge transfers. Without a state-tracking procedure, the meansquare displacement does not reach a diffusive linear regime, prohibiting mobility calculation. In addition, to ensure the convergence of the electronic mobility with the size of the system and the number of excited states, we have shown that a combination of the decoherence correction scheme and decoherence-induced spurious charge transfer correction is required.
Besides those paramount improvements to the surface hopping algorithm, we have also compared different definitions used in the literature for the mean-square displacement and for the electronic population definition. We have shown that the two commonly used definitions for the mean-square displacement (MSD and MSD coc ) give the same diffusion coefficient and the same mobility, whereas the third one (MSD var ), which is related to the spreading of the wavefunction, rather than to the diffusion of the charge carrier, yields always a zero slope as the polaron reaches a finite equilibrium size and does not grow indefinitely. Regarding the choice of electronic population to use in FSSH, we have compared the three definitions suggested in the literature (P wf , P PAS and P MQC ) and we have shown that these definitions coincide when a decoherence scheme is active.
In conclusion, we have established a well-founded set-up to run fewest switches surface hopping simulation of charge transport that converges electronic mobilities for different timesteps and different system sizes and that achieves detailed balance and good internal consistency.

Conflicts of interest
There are no conflicts to declare.