A symmetrical quasi-classical windowing model for the molecular dynamics treatment of non-adiabatic processes involving many electronic states

Electronically non-adiabatic processes—i.e., those involving molecular motion influenced by multiple potential energy surfaces (PESs) corresponding to multiple electronic states and the transitions between them—are of critical importance in many areas of molecular science. Examples include electronic excitation energy transfer in light-harvesting biological complexes, charge-migration in semiconductors, the electrical dynamics of energy storage devices, the dynamics of excited electronic states in atmospheric photochemistry, and so forth.

Nevertheless, though a variety of methods exist for treating these systems (of varying levels of complexity), there remain significant open questions as to the preferred approach in the context of the classical mechanical molecular dynamics (MD) framework routinely applied in computational chemistry for large molecular systems. The challenge for theory and simulation presented by such systems is that they would appear to be intrinsically highly non-classical: as stated, they involve the dynamics of electronic quantum states, the transitions between them, and the nuclear evolution on multiple PESs corresponding to those states. Quantum mechanics (QM), of course, provides the correct description, but the cost of a QM dynamics simulation scales exponentially with the number of degrees of freedom (DOF), which limits feasibility to relatively small system sizes, by chemistry standards.

Our work has thus focused on the development and application of a completely general, atomistic, classical trajectory-based method for treating electronically non-adiabatic molecular processes, which is thereby seamlessly compatible with classical MD simulation. The approach comprises two basic components: (i) Meyer and Miller’s (MM) representation of the electronic DOF in non-adiabatic processes as a collection of classical harmonic oscillators and (ii) a symmetrical quasiclassical (SQC) windowing model for quantizing these “electronic oscillator” DOF.

The approach is simple—it is essentially classical MD for the full vibronic dynamics in electronically non-adiabatic processes—but it has been shown to give remarkably quantum-like results in a fairly comprehensive set of application contexts^1–5 including gas and condensed phase systems, with electronic coupling ranging from extremely strong to moderate to weak to very weak, and in these disparate regime classes, accurately recovering electronic “quantum” coherence and decoherence effects, while in principle providing a reasonable treatment of detailed-balance⁶ (though not in all cases⁷). The approach is, furthermore, equally applicable in both diabatic or adiabatic⁸ representations, depending on what is convenient for the application at hand (though, in an adiabatic picture, a quasi/crude-diabatization can give improved numerical efficiency^9,10), and recent work by Geva and co-workers has also demonstrated¹¹ how classical MM oscillator DOF plus SQC quantization may additionally be used to represent nuclear quantum states, thereby providing a dynamically consistent means for incorporating nuclear quantum effects into the standard SQC/MM approach (though efficiency would generally dictate this be done selectively for a few key nuclear DOF).

Other recent applications of the SQC/MM approach include the work of Tao¹² treating single fission, the work of Provazza and Coker¹³ demonstrating the calculation of vibronic spectra (which employed a procedure¹⁴ for calculating the full electronic density matrix within the SQC framework), and our own work, in collaboration with Burghardt and co-workers, to treat a 20 monomer segment of oligothiophene, a model polymeric semiconductor. This application⁵ to oligothiophene—building on the prior work of Burghardt et al.¹⁵ to construct and treat QM’ly an ab initio-parametrized Frenkel-exciton Hamiltonian for this system—illustrated the capability of the SQC/MM approach to accurately simulate, classically, an elementary step in the excitation energy transfer processes in these systems, including concerted motion of the electronic and nuclear DOF and polaronic exciton formation and stabilization, all of which have physical significance to the potential use of these and similar materials in electronic devices.

This paper provides a next step in the development of the SQC/MM methodology in the form of a significantly improved SQC windowing model. The basic idea of the SQC component of the theory is to provide a semiclassical-like quantization of the electronic DOF without actually invoking any of the complicated (and potentially expensive) machinery of a true semiclassical calculation—i.e., the SQC/MM approach, just like any other classical MD approach, relies solely on ordinary classical mechanics.

In our work explaining how to use the SQC framework to calculate the full electronic density matrix, we showed¹⁴ how the SQC windowing model can be viewed as similar to a classical Wigner model where, instead of traditional Wigner functions, one employs “direct” Wigner functions, meaning that they are calculated by doing the Wigner transform of the electronic density matrix directly in action-angle variables. This gives δ-functions centered at the quantum values of the electronic actions (integer for the diagonal elements of the density matrix and half-integer for the off-diagonals), which are then broadened to give the SQC windowing functions. The question is, of course, what they should be broadened to, and this is where the modelistic aspect of the approach enters because experience shows that using narrow δ-function-like window functions gives poor results (even if one were not concerned with the difficulty of the root-search problem implied by the δ-function boundary condition). However, while the selection of the SQC windowing function might be viewed as an ambiguity, we do not want to give the impression that any window will do or that the choice is not crucial. Quite the contrary. The fact is that SQC/MM calculations are very sensitive to the form of the windowing function and our experience is that there are not many choices that work well at all, which leads into the topic of this paper.

In a previous work,¹⁶ a new SQC model based on triangle-shaped window functions was developed in order to extend the SQC approach to the regime of weak-coupling between the electronic states and, in so doing, the new SQC model was in fact shown to be universally superior to the original in all coupling regimes. This new SQC triangle model was therefore viewed as superseding the original square/histogram SQC windowing model with respect to the MD treatment of 2 electronic states. However, as illustrated below, the generalization of this model to more than 2 electronic states/DOF was not obvious: none was presented in the original paper and none has since been proposed.

The purpose of this paper is therefore to provide such a generalization, to explain its development and advantageous characteristics, and to illustrate its application to several standard benchmark problems. The organization is as follows: Section II A provides a review of the relevant theory, including the details of the SQC triangle window functions developed in Ref. 16 for 2 electronic states. Section II B then describes the issues encountered in extending the triangle model to more than 2 states, with the sought generalization given in Sec. II C and an optimal algorithm (and implementation) given in Sec. II D. The remainder of the paper is devoted to validating this algorithm through various statistical tests in Sec. III and example applications in Sec. IV. Section V provides conclusions.

The details of the SQC/MM approach have been reviewed in many recent papers. Comprehensive discussions are given, for instance, in Refs. 3–5 and 8. The focus of this paper is solely on the definition of the symmetrical quasi-classical (SQC) windowing functions which represent the electronic states and define the SQC quantization condition.

Briefly, in the MM model, the electronic DOF corresponding to a set of electronic states are mapped to a set of classical harmonic oscillators, which are coupled to each other (and to the nuclear DOF), and whose excitations determine which electronic state is occupied. Specifically, for a set of F diabatic electronic states, the MM Hamiltonian is given by

where {x_k, p_k} are the coordinates and momenta of the F “electronic oscillators,” (R, P) are the coordinates and momenta of the nuclear DOF, {H_kk′(R)} is an F × F diabatic electronic Hamiltonian matrix which depends on the nuclear coordinates, and γ is a zero point energy (ZPE) parameter. The time evolution of the classical oscillators in Eq. (1) thus describes the electronic configuration in the MM model and, in particular, the classical actions associated with each oscillator

(i.e., the classical analogues of the vibrational quantum numbers) represent the occupation of each state [each H_kk is multiplied by n_k in Eq. (1)], and they are specifically what are quantized via the SQC windowing protocol. Taking the time-derivative of Eq. (2), substituting Hamilton’s equations, and summing show that ∑_kn_k is conserved, i.e., each trajectory moves on a “polyad” of constant total action.

In the original SQC work treating electronically non-adiabatic processes,¹ the SQC windowing function for electronic state k of F electronic states was given by

where

i.e., as a product of F simple histogram windows requiring action n_k to be within ±γ of unity (the integer quantum value for state k) and the other actions {n_k′≠k} to be within ±γ of zero.

For instance, for F = 2 electronic states, the window functions defining electronic states 1 and 2 are just square boxes centered at (n₁, n₂) = (1, 0) and (n₁, n₂) = (0, 1), respectively, as shown schematically in Fig. 1. The use of simple histograms [Eq. (3)] is analogous to what had been done previously¹⁷ to symmetrically quantize the vibrational mode in the reactive scattering of H + H₂, and it was an obvious choice to try.

What was not so obvious was that the amount of ZPE built into the MM Hamiltonian, as specified by the parameter γ in Eq. (1), needed to be matched to the width of the SQC windowing functions—i.e., the same parameter γ appears in Eqs. (1) and (3b) (the actions being required to be within ±γ of the center of the window). It was also not obvious that this γ-parameter would have what appears to be an optimal universal value of about $1 3$ (though we later realized, we were not the first to have found better results with the MM Hamiltonian by reducing γ from its quantum value of $1 2$⁠, see Refs. 18 and 19).

The only real advantage of using the window functions defined by Eq. (3) is that they are the simplest one can envision (and that work). The triangular-shaped windows described in Ref. 16 and shown in Fig. 2 are a much better choice (for F = 2 states), one that is adapted to the dynamics generated by the MM Hamiltonian, which is essentially exact from a statistical standpoint, and which deals with Eq. (3)’s manifest inability to treat weak-coupling [because, as shown in Fig. 1, there is a gap between the windows of Eq. (3) in the action space].

In addition to the foregoing practical benefits (touching at a point in the action space in order to probe the weak-coupling limit, being statistically exact for a predetermined nuclear trajectory), they also intrinsically provide the optimal value of γ simply because $1 3$ is the centroid of a unit-width right triangle (i.e., the “correct” value of $γ= 1 3$ is essentially a geometric parameter arising from the shape of the windowing functions).

As stated, the purpose of this paper is to generalize the triangle-shaped windowing functions, essentially Fig. 2, to more than 2 electronic DOF/states. One possible generalization analogous to Eq. (3) is a set of F-dimensional “orthogonal simplex” windows given by

where

and

Equation (4) differs from Eq. (3) (for the histogram windows) in that the upper limits of the individual histogram factors have been removed in Eq. (4b), and instead the factor b(n) of Eq. (4c) appears as an upper bound to the sum of the actions in Eq. (4). For F = 2, this reduces to Eq. (15) from Ref. 16 which generates Fig. 2, and for F = 3, Eq. (4) generates the pleasing looking array of 3 three-dimensional scatter plots shown in Fig. 3.

The b(n) factor in Eq. (4) makes these window functions non-separable in the individual {n_k} and defines a line/plane of maximum total action and total ZPE in the ensemble (summed over the F electronic DOF). Note that the trajectory in the ensemble with no ZPE is the one initialized at exactly the bottom right-angled corner of whatever window function corresponds to the initial state (in Figs. 1, 2, or 3) and its time-evolution maps out the no-ZPE polyad which is furthest from b(n) and this is, in fact, the sole Ehrenfest trajectory. Also note that when binning the final actions with the W_k(n) of Eq. (4), there is no need to explicitly evaluate b(n) [Eq. (4c)] because ∑_kn_k is a constant of the motion. Nor is it required to check the values of the other action variables {n_k′≠k} [with Eq. (4b)] as they cannot fall below −γ dynamically. In sum, Eq. (4) is one possible generalization of the F = 2 triangle windows and to apply it one samples initial actions from Eq. (4) and windows/bins the final actions simply by assigning each trajectory to a state k when

of which there will be only one state k for a given trajectory.

For any F > 2, however, there is a problem with the value of γ implied by Eq. (4) which is a simple consequence of the geometry in higher dimensions: the centroid of an F-dimensional unit orthogonal simplex is (F + 1)⁻¹. Thus, for example, the implied γ value of the 3-dimensional windows shown in Fig. 3 is not $1 3$ but $1 4$⁠, as indicated in Fig. 3 ²⁰ and, in general, with Eq. (4), as the number of states F gets large, the γ-parameter goes to zero. Not only does this strike one as being a poor model, but SQC/MM calculations have been done using the original square histogram windows for F > 2 and, at least up to F = 7, it has been verified that $γ≈ 1 3$ is still effective. See, for example, the SQC/MM calculations for the Fenna-Matthews-Olson (FMO) complex in Ref. 3 or those for the oligothiophene 20-mer in Ref. 5.

The other potential issue evident from Fig. 3 is that there is no point in the action space where the windows for all three states touch. It is hard to conclude definitively that this is a fatal defect but consider, in view of Fig. 3, the F = 3 scenario of electron donor and acceptor states only coupled together indirectly through a third bridge state. In the weak coupling limit, where trajectories move only slightly in the action space, there would be some trajectories that move from state 1 (Fig. 3, red window) to state 2 (green window), and some trajectories that move (again, only slightly) from state 2 (green window) to state 3 (blue window), but they are not the same trajectories. For a single trajectory to traverse all 3 windows by moving only slightly, the windows for all 3 states would at least need to touch at a central point in the action space. This reasoning just generalizes the reasoning of Ref. 16 for having the F = 2 triangle windows touch at a point. Again, it is hard to claim definitively that the effectiveness of a SQC windowing model absolutely depends on having all the windowing functions meet at a central point in F-dimensions, but nevertheless, this characteristic, if present, does give a more intuitively satisfying picture of the partitioning of the classical action space into regions which characterize the quantum states.

In discovering a high-dimensional SQC windowing scheme which possesses the foregoing characteristics, it has been helpful to invoke a more stringent criterion than simply the requirement that the generalization reduces to the triangle windows for F = 2 states. Instead, it is useful to require that the generalization preserves the triangle window statistics between all pairs of states, regardless of the total dimensionality F of the electronic space. Obviously, this is not the case for the triangle/simplex scheme of Eq. (4) because it does not even keep the same value of γ independent of F: lim_F→∞γ = 0.

It turns out that a windowing scheme of arbitrary dimension F which possesses these properties—meeting at a central point, having $γ= 1 3$ (for all F), and, moreover, preserving the triangle window’s statistical treatment between all pairs of states—is found, surprisingly, not in the shape of a unit right/orthogonal simplex, but in the shape of a unit right pyramid. Such a generalization is given by

where

and

Figure 4 illustrates these windowing functions for F = 3. As shown in Fig. 4, the central point at which these generalized triangle windows meet (to handle weak-coupling) is

(where I is the F × F identify matrix) and, as claimed, Eq. (6) preserves the value of $γ= 1 3$ for any dimension F, though, to achieve this, it is noteworthy that there is a density weight factor present in Eq. (6b).

To understand the role of the density weight factor, it is useful to verify the value of γ. It is easiest to first calculate the normalization factor, determined by evaluating the density weighted volume, which (as shown) is independent of F,

which, letting $n k → n k +1−γ$⁠,

Next, the task is to prove that $⟨ n k ′ ⟩ W k = δ k k ′$⁠, i.e., the expectation value of action n_k over the window for state k is 1, that the other action variables n_k′≠k average (over W_k) to 0, and that this is so if and only if $γ= 1 3$⁠,

which, letting $n k → n k +1−γ$⁠,

and is 1, if and only if $γ= 1 3$⁠. Likewise, for the other actions {n_k′≠k}

which, letting $n k → n k +1−γ$⁠,

and is 0, if and only if $γ= 1 3$⁠. This analysis is, of course, equivalent to the statement that the centroid of the generalized windowing functions specified by Eq. (6) [with the density weight factor in Eq. (6b)] is $1 3$⁠, and indeed proving this is even slightly simpler if one locates the windowing function at the origin and then performs the foregoing integrations. In any case, the integrations reveal that the role of the density weight factor in Eq. (6b) is to balance the statistics of the expanded number of DOF (i.e., by cancelling the integrations over these other DOF). Of course, for the case of F = 2, the density weight factor vanishes and Eq. (6) does reduce to the original triangle windowing scheme.

This leaves the question of the more stringent criteria described above, for which having $γ= 1 3$ is a necessary but not sufficient condition. The simplest way to understand that the criteria is satisfied is to see it graphically. To preserve the triangle statistics between all pairs of states of F total states is to require that trajectories be initialized so that the projection of their F action variables onto any 2-dimensional sub-plane n_i–n_j gives precisely a 2-dimensional triangle of uniform density (i.e., what is shown in Fig. 2 for the original F = 2 triangle windows). Figure 5 shows such projections for the three sub-planes in the F = 3 action space: panel (a) contains the original 3-dimensional scatter plot for the 3 states, and panels (b)–(d) give the same plot (and same data), rotated to show the 3 possible n_i–n_j sub-planes. For all 3 orthogonal views, a pair of triangles of exactly uniform density is revealed. Further statistical validation is given in Sec. III.

To make it absolutely straightforward to take advantage of this new generalized SQC triangle windowing protocol, provided in this section is a practical implementation given in C++ for concreteness. Algorithm 1 (below) is for Monte Carlo sampling the initial actions. Algorithm 2 (below) is for binning the actions at each desired final time.

With regards to Algorithm 2’s binning procedure, first, note that there is no density weighting, and this is because the correct (triangle) statistics are given by only weighting the initial conditions.²¹ (Also, a lack of density weighting in the binning procedure is likely to converge with fewer trajectories.) The “correctness” of the sampling and binning procedure in terms of statistical accuracy will be confirmed in Sec. III. The physical significance of the initial sampling is that it sets the ensemble’s distribution of ZPE, which is conserved dynamically because ∑_kn_k is a constant of the motion; experience has shown this to be the critical aspect.

Second, note Algorithm 2’s binning procedure for assigning trajectories to the various states: A trajectory is assigned to state k if

This windowing/binning condition is totally analogous to that given in Eq. (5) for the triangle/simplex windowing scheme except that, here, the presence of the second inequality is necessary to guarantee that each trajectory only contributes to a single state k. It would seem that this is the most efficient “in-or-out” division of the dynamically allowed action space into quantum states that can be envisioned (and still be accurate/correct). Once again, a strict application of Eq. (6) which imposes additional upper bounds on the actions has been tested, and it appears to be just as accurate as using Eq. (8), but it is significantly less efficient for large F. In short, use Eq. (8) as implemented in Algorithm 2.

The windowing procedure set forth in Sec. II is simple to implement and perhaps even intuitive but, because it is not obvious how one could derive it, ab initio, the task of this section is to provide some statistical validation. For this purpose, it is useful to return again to the “universal” scattering model from the original Meyer-Miller paper, which was also used to justify the choice of triangle windows for the F = 2 case in Ref. 16.

The universal model problem is based on nothing more than the assumption that the initial and final electronic amplitudes in a non-adiabatic process may be given by a single unitary transformation

independent of the nuclear DOF. Thus, it provides a sort of statistical validation of a computational method with respect to a baseline process in which the nuclei are considered fixed or, equivalently, of a sufficiently short-time/instantaneous process (where the nuclei do not have time to move) or, more generally, of a process where all the nuclei are assumed to follow the same predetermined trajectory.

Recapping what was done in Ref. 16 for 2 states: a 2 × 2 scattering/S-matrix in Eq. (9) may be parametrized in terms of two angles, α and β, and a third parameter p₁₂

The parameter p₁₂ obviously sets the transition amplitude P_2←1 affected by application of the matrix S to the initial electronic amplitudes c(0), and so by systematically varying parameter p₁₂, all non-adiabatic transition regimes may be “universally” evaluated. In the context of testing the statistics of the SQC protocol, this amounts to the following:

• sampling initial actions n(0) using Algorithm 1;

• calculating c(0) from n(0) using the relation

(with angles q(0) selected randomly);

• (iii)

  selecting a p₁₂ with α and β given arbitrary values (since the averaging over q removes the dependence on α and β);

• (iv)

  applying Eqs. (9) and (10) for the selected p₁₂ to give the final amplitudes c(∞);

• (v)

  calculating the final actions n(∞) using Eq. (2); and

• (vi)

  binning n(∞) using Algorithm 2.

The result of this procedure is given in Fig. 6 and, as explained in Ref. 16, it is seen that the triangle windows give the exact result (the diagonal line P_2←1 = p₁₂) for all possible S-matrices and transitions, whereas the original square histograms only approximately fit the exact result and fall off completely in the weak-coupling limit.

The goal, then, is to perform a similar procedure for the generalized triangle windows, for which F = 3 is the simplest case they can be evaluated. For this purpose, rather than employing an expression for a completely general 3 × 3 unitary matrix, it is more straightforward to consider 3-state transitions as compositions of transitions between pairs of states and, because only transition probabilities will be evaluated, to ignore any phase-dependence, i.e., to employ the real matrices

and

and to sequentially apply them to generate 3-state transitions over ranges of values for the parameters p₁₂, p₂₃, and p₁₃.

It is, however, still difficult to visualize (and scrutinize) results dependent on 3 independent parameters, so specifically considered will be the SQC calculation of P_3←1 for the two possible distinct (non-trivial) ways of applying any two of the forgoing S-matrices as follows:

• a donor-bridge-acceptor type transition

for which the exact result is

and

• (ii)

  an initial “depletion” of state 1 via transition to state 2, followed by a transition from state 1 to state 3,

for which the exact result is

SQC results using the generalized triangle windows for the 3-state donor-bridge-acceptor model of Eq. (13a) are plotted in Fig. 7(a) for the full range of possible values of p₁₂ and p₂₃. The results exhibit a linear dependence on p₂₃ with a slope equal to p₁₂ in accord with the exact result of Eq. (13b). The results are not perfectly linear, but they are indeed very close.

Likewise, SQC results for the case of the 3-state initial “depletion” model of Eq. (14a) are plotted in Fig. 8(a), for this case in terms of p₁₂ and p₁₃. Here, from the exact result of Eq. (14b), what is expected in Fig. 8(a) are linear functions of p₁₃ with slopes equal to 1 − p₁₂ and that is what is seen. Again, the functional dependence is not perfectly linear, but it is nevertheless very close. Note that the plot-line in Fig. 8(a) for which p₁₂ = 0 is exactly the 2-state universal problem in the presence of an uncoupled third state. The agreement is perfect which is, of course, a consequence of the sub-plane projection analysis given in Sec. II.

In addition, it is reassuring to verify that the results for the 3-state problems shown in Figs. 7(a) and 8(a) are not affected by the presence of additional uncoupled electronic states/DOF. Thus, while panels (a) of Figs. 7 and 8 show results for 3-state transitions in an electronic space of dimension F = 3 (i.e., 3 action variables, n₁, n₂, and n₃), panels (b) and (c) go further showing analogous results in an expanded electronic space of dimension F = 10 and F = 100, respectively (i.e., having action variables n₁, n₂, … n₁₀ and n₁, n₂, … n₁₀₀). Comparing panels (a)–(c) of Figs. 7 and 8, the results are seen to be identical up to sampling variation, explicitly illustrating the desired invariance to electronic space dimension F, even for up to 100 electronic states and action variables. Of course, in view of the sub-plane projection analysis/arguments made in Sec. II and, particularly, in view of Fig. 5, this comes as no surprise. It should also be noted in all these results that the weak-coupling limit is handled properly.²² 

In summary, the results of Figs. 7 and 8, taken together, illustrate the robust statistical map provided by the new windowing protocol.

The most important validation of the new SQC windowing protocol is its application to realistic benchmark models with significant coupling between the electronic and nuclear DOF. This is now done for the following example applications, many of which are Frenkel/site-exciton Hamiltonians parameterized to model light-harvesting complexes. Some of these models were treated previously in Ref. 3 with the original square histogram SQC windows and, thus, as an initial matter, it is reassuring to verify similar results here.

Section III demonstrated that the new windowing model extends the advantages of the triangle windows—a robust and accurate statistical description of the electronic DOF and the seamless handling of the weak-coupling limit—to the treatment of higher numbers of electronic states. For the following applications, however, weak-coupling is not an issue. What is an issue is the efficiency of the windowing model as the number of states F becomes large and it turns out that vastly improved efficiency for large F is the most important benefit of the new model.

Accordingly, the accuracy and efficiency of the new generalized SQC triangle model is now evaluated for non-trivial model problems having significant coupling between the electronic and nuclear DOF and progressively increasing numbers of electronic states. Beginning with the case of F = 3, Fig. 9 shows the treatment of the 3-state spin-boson problem from Ref. 23 with SQC results plotted versus QM path integral results, and the agreement is reasonable. Note that an SQC treatment of this same 3-state spin-boson problem was given in Ref. 4 but, as explained in Sec. II, that treatment was not generalizable to more than 3 states.

Next, the new scheme is validated for 4-state, 5-state, and 7-state electronic models. For the purpose of progressively increasing dimensionality, a simple approach is to take a well-studied higher-dimensional problem and to remove some of the less important states. This has been done for the well-known 7-state model of the Fenna-Matthews-Olson (FMO) complex from Ref. 24, which was reasonably handled by the original square histogram SQC model in Ref. 3.

At a temperature of 77 K, when the FMO system is initialized with the excitation at site/pigment 1 of the chromophore complex, only 4 of the 7 site-specific electronic states play a significant role in the dynamics. Deleting the other 3 states from the model gives a reasonable 4-site/state benchmark (see the  Appendix for the exact electronic site-energy/coupling matrix). Figure 10 shows SQC results using the new windowing protocol versus benchmark hierarchical equations of motion (HEOM) results for both 4- and 7-site versions of the problem. The results are interesting.

For the standard 7-site version, panel (a) of Fig. 10, the result is essentially the same as that calculated in Ref. 3 for the square histogram windows: the agreement with the exact HEOM result is reasonable (and the similarity to Ref. 3 is reassuring), but there is clearly too much decay from the initial state 1, with about half of the un-retained population transferred to pigment 3, and the rest, presumably, ending up distributed amongst pigments 5–7 (which otherwise play little role). Figure 10(b) shows results for the analogous 4-site model with pigments 5–7 deleted, and for this reduced dimensional model, the SQC calculations give results in very close agreement with the exact HEOM calculations. This is an important validation of the windowing scheme, but it is also interesting because it suggests that it is the SQC/MM treatment of the “unimportant” states that is causing the deviation from the exact results. It is also surprising that the exact HEOM results actually show less decay from initial state 1 in the presence of states 5–7, rather than the dynamics just exhibiting invariance to their presence.

In any event, one plausible explanation for the SQC/MM model’s discrepancy when pigments 5–7 are present is that, at the low simulation temperature of 77 K, nuclear quantum effects are playing a role. In particular, nuclear ZPE effects may be “locking out” certain vibrational DOF of pigments 5–7, whereas all vibrational modes are open for energy exchange in fully classical MM dynamics. The very good agreement between SQC and HEOM in the 4-state version certainly supports this view, but an additional confirmation comes from running the same simulations at 300 K. Figure 10(c) shows such results and the deviation between the SQC and HEOM treatments is significantly reduced. The SQC/MM approach is designed with classical MD-level efficiency in mind and, although electronic coherence effects and the like are captured, it is not intended to deal with nuclear quantum effects. In any event, the task here is validate the new windowing scheme and, nuclear quantum effects aside, panels (a)–(c) of Fig. 10 do collectively indicate that the statistical description of the electronic DOF is robust and accurate.

Further confirmation is provided in Fig. 11 which shows SQC results for the population dynamics of the FMO complex when site/pigment 6 is initially excited. In this case, the dynamics primarily involves 5 pigments, so 2 site-local electronic states may be removed to provide a reasonable 5-state benchmark problem (see the  Appendix for the specific electronic site-energy/coupling matrix). Analogously to Fig. 10, panel (a) of Fig. 11 shows the full 7-site/state problem at 77 K, panel (b) the 5-site version of the model at 77 K, and in panel (c), results for the full 7-site model at 300 K. In all cases, the SQC/MM approach replicates quite well the exact HEOM results. This may be unsurprising since the original square histogram windowing model also provided excellent results for the full 7-site model initialized in state 6 at 77 K. Nevertheless, it is important to demonstrate this for the new windowing scheme and for the particular value of F = 5—in addition to F = 3, 4, and 7 above—and particularly so in view of the density weight factor in Eq. (6b) having an explicit exponential-in-F dependence. In sum, the new SQC windowing protocol gives reasonable results for relevant systems of 3, 4, 5, and 7 electronic states/DOF.

Other than for systems where treating weak-coupling is crucial, the most important benefit of the generalized SQC triangle windows is in handling problems involving many electronic states, where its superior statistical efficiency becomes paramount. Table I gives trajectory statistics for the various model problems treated in this paper, specifically, the approximate average fraction of trajectories which are windowed at any point in time towards the end of the given SQC/MM simulation²⁵ using both the original square histogram windows and the generalized triangle windows. Table I shows that, in all cases, the triangle windows give superior numbers, but that for problems having many relevant electronic states—i.e., double-digit numbers of states—so few trajectories are windowed with the original square histogram SQC model that the calculations effectively become infeasible. The generalized SQC triangle windowing scheme, on the other hand, is able to effectively draw statistics from sufficient numbers of trajectories to enable successful calculations for every example tried, even for the model system having nearly 100 electronic states.

Effective treatment of high numbers of electronic states is first demonstrated with respect to the FMO trimer of 3 × 8 = 24 states (as indicated in Table I). The 7-site model considered above (of a single FMO monomer) is actually missing a later appreciated 8th pigment. The matrix of site-energies and intra-monomer electronic couplings associated with this 8-state model is given in Ref. 26. In conjunction with the 8 × 8 inter-monomer coupling matrix given in Ref. 27, the full 24-state Frenkel/site-excitonic model of the FMO trimer may be assembled as

as described in Ref. 28. Reference 28 also provides the effectively exact, course-grained quantum path-integral results which are taken here as the benchmark. The monomer excitonic site-energy/coupling matrix H_mono and inter-monomer couplings H_inter are reproduced in the  Appendix.

Furthermore, as a means of assessing the feasibility of very high-dimensional calculations with the new approach, one can go a step further and construct fictitious FMO multi-mers having many times more electronic states. Thus, as indicated in Table I, hexamer and dodecamer FMO models have been constructed with 8 × 6 = 48 and 8 × 12 = 96 electronic states, respectively. For example, the fictitious FMO hexamer model treated below is given by

Figure 12(a) displays the SQC treatment of the FMO trimer model with the exact QM path integral results alongside in Fig. 12(b). In all calculations here involving FMO multi-mers, the simulation temperature is 300 K and pigment 1 of monomer 2 is given the initial electronic excitation. The plots show the resulting short time dynamics (out to 500 fs, top plots) involving, almost exclusively, intra-monomer excitation energy transfer amongst the pigments of monomer 2. Inter-monomer energy transfer occurs on the time scale of picoseconds, and this is shown in the middle and lower plots of Fig. 12 with different vertical-axis magnifications and color-coding of the pigments.

Overall, the SQC results of Fig. 12(a) compare very well with the exact QM path-integral results of Fig. 12(b), including the SQC calculation¹⁴ of the largest off-diagonal element of the electronic density matrix, ρ₁₂, which decays in a few hundred femtoseconds.²⁹ There are some discrepancies, most notably with the 5× magnified vertical-axis of the lower plots, one discerns a little less intra-monomer energy transfer to pigments 3 and 8 of monomer 2 (versus the exact result) and a little more inter-monomer energy transfer to monomers 1 and 3. Nevertheless, relative to unit probabilities, these differences look to be roughly 2%, which, given the crudeness of the model Hamiltonians to begin with, seems quite reasonable. Moreover, as has been emphasized before, because the SQC/MM approach is totally classical and provides an atomistic MD simulation of the motion of the nuclei (here, the bath modes associated with each pigment/site), there are no restrictions placed on the form of the Hamiltonian which may be treated dynamically with the MM/SQC approach: long non-Markovian bath memory times are not an issue; site-specific spectral densities are not an issue; realistic, highly peaked spectral densities are, if anything, simply more efficient to discretize; and, again, because the SQC/MM method is basically classical MD, one has wide choices available as to the analytic form of the atom-to-atom interaction potentials.

It is also important to take note of the considerable efficiency of the SQC/MM approach versus, in this case, a QM path-integral calculation. Quantum mechanically, in the path integral view, there is an exponential increase in the number of viable paths whenever there is a non-vanishing probability of switching electronic states. The SQC/MM approach is efficient because instead of dealing with this exponential path expansion, the dynamics just evolves on the occupation-weighted average of the PES associated with the different electronic states. This tremendous simplification is born out in the single-processor run-times reported in Table II. Table II shows that what takes days to compute quantum mechanically (even with the substantially improved course-grained approach of Ref. 28) is reduced to mere minutes using the SQC/MM approach. Moreover, as stated, these are single-processor run-times. The classical SQC/MM approach is absolutely trivial to parallelize over multiple processors (without shared memory) and so system size can generally be extended linearly with processor number, which is not generally the case with approaches that cannot leverage the locality built into classical mechanics. Of course, the question is whether classical mechanics and, specifically, classical mean-field Meyer-Miller (MM) vibronic dynamics is a good enough approximation. Here, for a problem and regime of obvious interest to molecular bioscience, it does indeed appear (from Fig. 12) that the SQC/MM approach quite reasonably emulates the true electronically non-adiabatic/excitonic quantum dynamics.

As stated above, one would also like to assess the feasibility of handling systems having even higher numbers of relevant electronic states, e.g., say, 50 or 100, and one simple procedure, as described above, is to construct a fictitious FMO hexamer of 48 states and a dodecamer model of 96 states. Figure 13 shows the SQC/MM results, which look very similar to those for the FMO trimer model, the primary difference being the stronger quenching of monomer 2’s initial electronic excitation (note the difference in vertical axis scales in the lower panels), which is easily explained by the fact that there are many more pigment channels amongst which the electronic excitation may be distributed. These results are thus just what one would expect and are, perhaps, by themselves, not too significant. What is significant are the trajectory statistics listed in Table I and their interpretation: While windowed trajectory fractions of 12% and 7% for 48 and 96 states, respectively, are not large percentages, they are at least a factor or 10³ greater than when using the square histogram windows and not so small as to render the calculations impracticable. Thus, the results in Fig. 13 look reasonable, and they were generated with the same sized trajectory ensemble as all the other results presented in this section (i.e., 10 000, see below).

Table I shows, however, that even with the new windowing scheme, the windowing fraction does tend towards zero as the number of states gets larger and larger. It is therefore natural to wonder whether there is some limit where even the new scheme becomes intractable.

A definitive answer to that question will have to wait, but there are a few encouraging points worth noting: First, for a large system with many electronic states coupled to each other and to many dispersive nuclear DOF—such as our fictitious FMO dodacamer—at long times, when the electronic DOF have equilibrated, the electronic energy will tend to be distributed amongst the many sites, which is manifested by an ensemble of trajectories that tend to have their total electronic action split amongst their different electronic modes, thus tending to be clustered in the center of the electronic action space. According to Eq. (8), the central region of the action space where n_k < 1 − γ, ∀k does not fall within any window function (the outer edges of this region are actually visible in Fig. 4 for the 3 state case), and hence, one anticipates potentially poor trajectory statistics for high-dimensional problems at long time. What is encouraging, however, is that because all the windows touch at a central point, there would seem to always be some statistical probing of the population dynamics, even in this limit where the electronic energy is roughly evenly distributed amongst the different states. The second encouraging point is that, in the long-time limit, the interesting dynamics has likely already occurred, and therefore, it probably does not take great statistics to assess the outcome, which is precisely what one sees in Fig. 13, with the electronic population distributed across a sea of states/pigments and being essentially constant between 15 and 20 ps.

One final note regarding efficiency: practitioners often wonder how many classical trajectories are necessary to achieve reasonable results. Every result presented in Sec. IV was calculated from an ensemble of 10 000 trajectories (in some cases, applying a Bezier smoothing with the plotting software to beautify the raw results). SQC calculations can look a bit “ragged” because the binning criterion is a binary operation, but experience shows the statistics to be robust and, generally, even results that appear noisy are not unreasonable. In view of these points, the ensemble size was kept fixed at 10 000 regardless of the problem—whether treating 3 electronic states (Fig. 9) or 96 electronic states (Fig. 13)—which is enough to be confident in the results, but not so many as to smooth away imperfections. As a practical matter, these would likely be the relevant considerations when treating realistic problems where trajectory calculations may be expensive.

Finally, all of the results reported in this paper make use of an unrelated “γ-adjustment” procedure which is the exclusive subject of a contemporaneous paper. For most of the model problems reported here, the procedure has only a small effect (if any). Nevertheless, because this additional protocol can, it some cases, give markedly improved results, it is our standard practice to employ it and, thus, it is used consistently throughout this paper. The goal, of course, is to have a standard, universal, classical algorithm applicable in all regimes without any problem-specific tunings or modifications.

This paper has provided a generalization of the 2-state triangle SQC windowing scheme of Ref. 16 to higher numbers of electronic states and illustrated its application to several benchmark problems. Surprisingly, the preferred generalization takes the form of density-weighted right-angle pyramidal-shaped windows, rather than a high-dimensional simplex geometry as one might initially suppose.

By construction, this new generalized triangle SQC windowing scheme carries forward the benefits of the 2-state triangle windows—seamless treatment of the limit of weak coupling between the electronic states and improved accuracy in all coupling regimes—to the treatment of systems having more than 2 relevant electronic states. More importantly, however, the new scheme provides vastly improved trajectory statistics (e.g., by factors of 10³) when treating large numbers of electronic states, thereby bringing complex/high-dimensional electronically non-adiabatic and excitation energy transfer problems within practical reach of the SQC/MM approach. Furthermore, though its development has manifested some subtleties, the new method/algorithm is actually extremely simple to implement, and a full implementation, requiring just a few lines of C++ code, is given in Sec. II D.

With these benefits and without apparent drawbacks, we recommend this new windowing model as the preferred SQC scheme in all circumstances. Of course, it is also true that although the new algorithm is simple, accurate, efficient, and powerful, and developed to satisfy a comprehensive and pragmatic set of criteria (and it is not easy to envision other models that would do so), it is nevertheless conceivable that even better procedures are still possible, though we believe in a spirit similar to what has been developed here.

The authors thank Addison Schile, Ph.D. student in the Pitzer Theoretical Center at UC Berkeley, for performing the HEOM calculations for the for the 4- and 5-state reduced-dimensional FMO models used as benchmarks in Figs. 10 and 11.

The authors also thank Dr. Martin Richter, formerly of the Max Born Institute, Berlin, and now at the Friedrich Schiller University, Jena, for providing the course-grained QM path integral calculations of Ref. 28, used as a benchmark in Fig. 12.

This work was supported by the National Science Foundation under Grant No. CHE-1464647 and by the Director, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

In addition, this research utilized computation resources provided by the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The standard 7-site FMO excitonic energy and coupling matrix from Ref. 24 is

where the site energies {ϵ_k} and non-adiabatic couplings {J_k,k′} are in cm⁻¹.

Since states 5–7 are relatively unimportant when the initial excitation is on site/pigment 1, a reasonable 4-state benchmark problem is given by cancelling the rows and columns of the matrix corresponding to states/pigments 5, 6, and 7

Likewise, when the initial excitation is on site/pigment 6, the dynamics does not significantly involve excitations on pigments 2 and 3. Thus, cancelling rows and columns corresponding to states 2 and 3 gives a reasonable 5 state benchmark

For the FMO trimer, following Ref. 28, the 8-state monomer site/coupling matrix is taken from Ref. 26 

and the inter-monomer couplings from Ref. 27