A new symmetrical quasi-classical model for electronically non-adiabatic processes: Application to the case of weak non-adiabatic coupling

In recent years we have been pursuing the question of how well some very simple classical molecular dynamics (MD)-based approaches can describe electronically non-adiabatic processes in molecular systems.^1–7 The practical motivation for this is obvious: only classical MD simulations can treat very large molecular systems with essentially arbitrary interaction potentials, and electronically non-adiabatic processes are ubiquitous in many important areas of molecular science. The foundation of the approach is to treat electronic and nuclear degrees of freedom (DOFs) equivalently, so as to describe the dynamics of their interaction consistently, albeit classically, and the first step in this regard is to use the Meyer-Miller⁸ (MM) classical model for the electronic DOFs, to go along with the standard classical treatment of the nuclear DOFs. The resulting classical vibronic (i.e., nuclear + electronic) Hamiltonian (which can be derived in various ways, heuristic or more rigorous⁹) is actually an exact Hamiltonian in the sense that if the classical nuclear and electronic coordinates and momenta were all replaced by quantum mechanical operators in the standard way, and the resulting Hamiltonian operator used in the Schrödinger equation, this would describe the quantum vibronic dynamics exactly. The approximation is that we treat both nuclear and electronic DOFs classically, i.e., by computing classical trajectories in these nuclear and electronic variables.

The second aspect of the approach is how electronic state information is extracted from such classical trajectories, for example, to calculate transition probabilities from one electronic state to another. Semiclassical (SC) theory^10,11 provides a general and well-defined prescription for how to do this, but a rigorous implementation of SC approaches leads to a much more involved calculation than a straight-forward classical MD simulation. For this reason we have focused on very simple quasi-classical (QC) models—which originated¹² in the 1960’s for treating state-to-state vibrational transitions—to see how well they can work for the “electronic harmonic oscillators” of the MM Hamiltonian, especially when implemented with some of the useful “tweaks” that have been added to the QC technique in the intervening years: the “symmetrical” QC (SQC) idea suggested many years ago¹³ but not pursued, Stock’s suggestion^14,15 of including less than the full zero point energy (ZPE) in the electronic oscillators of the MM model, and Bonnet and Rayez’s¹⁶ introduction of window functions of reduced size.

The “standard” version of this SQC/MM approach that we have focused on has been described in detail in our recent papers (for example, Refs. 2 and 6) and is summarized briefly in Sec. II. It has shown a remarkably good agreement with accurate quantum benchmark calculations for a variety of examples involving electronic transitions, from 1-D scattering problems to a variety of models of complex molecular systems (e.g., various spin-boson type models). It has been seen to quantitatively describe both coherence effects in the dynamics, as well as the quenching/“de-coherence” of such effects, all coming out of a straightforward classical MD simulation with no additional inputs nor ad hoc additions to the theory. Though the classical dynamics resulting from the MM Hamiltonian was noted from its origin to be “Ehrenfest dynamics”—i.e., the force on the nuclei at any time is the coherent average of over all electronic states—we have explained how the boundary conditions implicit in the SQC procedure eliminate the Ehrenfest method’s well-known deficiency of emerging from a region of electronic coupling in an intermediate electronic configuration; relatedly, we have shown how the boundary conditions further ensure that detailed balance is described in a reasonable way⁴ (though not necessarily exactly). Most recently, we have shown⁷ that in addition to accurately describing the time-evolution of the electronic state populations—i.e., the diagonal elements of the electronic density matrix—off-diagonal elements of the density matrix may also be extracted from a standard SQC/MM MD simulation simply by defining additional window functions corresponding to these off-diagonal elements.

However, the SQC model is such a simple and primitive way for “quantizing” the electronic DOFs that it has been clear to us from the beginning that it could not always perform as well as we have observed. One situation where we have anticipated it to fail is for the case of very weak coupling between the electronic states. This is easy to understand by referring to Fig. 1 which depicts the square histogram window functions of the standard SQC approach for the case of two electronic states (and thus two electronic oscillators in the MM model). In the SQC model the initial electronic action variables (n₁, n₂) begin in one of the histograms, and the probability of being in the other state at time t is the fraction of such trajectories that have their action variables in the other histogram at that time. If the electronic coupling were zero, the initial action variables would of course retain their initial values for all times and the transition probability would be zero for all times, which is correct. However, if the coupling is very small but not zero, the true transition probability will be very small (but not absolutely zero), yet for no trajectories will the action variables have changed enough to fall within the bounds of the other histogram, thus yielding a zero transition probability. Our supposition, therefore, has been that our “standard” SQC approach would fail in the case of exponentially small transition probabilities. (This has been a well-recognized problem with the QC approximation for vibrationally inelastic transitions in earlier years.)

The purpose of this paper is to describe a new SQC model—a variation on our original—which does in fact describe exponentially small transition probabilities quite well, and as a bonus it describes transitions in the “normal” regime (i.e., probabilities ≳0.1) even better than our previous standard model. Section II begins with a brief overview of the classical Meyer-Miller (MM) vibronic Hamiltonian and a short description of how electronic states are determined within the symmetrical quasi-classical (SQC) framework. Section III then describes a model problem utilized in the original MM paper (a universal 2-channel scattering problem) and treats it via our “standard” SQC/MM approach to illustrate (and quantify) the problem with very weak transitions. A solution to the problem is then given in Section IV, which presents a new windowing scheme seen to treat the weak coupling regime as well as normal coupling regime exceptionally well. Section V presents our results of applying this new SQC/MM approach to more realistic problems, beginning with validation of the new approach in the normal coupling regime by treating a set of challenging spin-boson benchmarks for which the original SQC/MM approach did very well. More realistic examples in the weak coupling regime are then explored including a version of Tully’s¹⁷ 1-D single avoided crossing scattering problem (Tully #1), but modified with exponentially reduced values of the non-adiabatic coupling, and also a treatment of a very weakly coupled 2-state site exciton model similar to those treated in Ref. 6. Our conclusions are summarized in Section VI.

We have recently provided a thorough review of the SQC/MM methodology in conjunction with our work⁶ treating the light-harvesting Fenna-Mathews-Olson (FMO) complex, so only the bare essentials are covered here. The two core ingredients of the methodology are (i) the Meyer-Miller (MM) classical vibronic Hamiltonian⁸ and (ii) a symmetrical quasi-classical (SQC) windowing procedure^1,2 for “quantizing” the electronic DOFs embodied as harmonic oscillators in the MM Hamiltonian.

The MM Hamiltonian was originally given in terms of action-angle (a-a) variables, a pair for each electronic state,

where the first term is the nuclear kinetic energy, and {H_k,k′ (R)} is the diabatic electronic matrix (which depends parametrically on the nuclear coordinates R); there is an equivalent adiabatic version of the Hamiltonian. Calculations (i.e., numerical integration of Hamilton’s equations), however, are usually more conveniently carried out in the Cartesian electronic variables which are related to the a-a variables (by canonical transformation) as follows:

or inversely,

and in which terms the MM Hamiltonian for say F electronic states is

The MM representation thus characterizes the F electronic states as F (coupled) harmonic oscillators, the excitation of which represents the occupation of that electronic state, and electronic transitions emerge in this picture as “vibrational” transitions between these “electronic oscillators.” A classical simulation thus adds F vibrational-like DOFs to the set of (perhaps very many) nuclear DOFs, and is thus typically a very modest addition to a MD simulation with only the nuclear DOFs—provided, of course, that one has the potential energy surfaces (PESs) for the various electronic states and their couplings, which may come from rigorous “quantum chemistry” (possibly computed on-the-fly), less-rigorous density functional theory (DFT), or from a semi-empirical “molecular mechanics” force field as often used in bio-molecular simulations.

The parameter γ in Eqs. (1) and (4) sets the zero point energy (ZPE) of the electronic oscillators in the MM model. Originally, MM took it to have the QM value of $12$⁠, but Stock observed^14,15 that better results could be obtained in a classical simulation by using a value less than $12$⁠. In the SQC/MM model, we have thus also taken γ to be an empirical parameter, though (as described below) we have essentially found the same value to be optimal in every application.

There is also a symmetrized version of the MM Hamiltonian which has proved convenient (and is used in all calculations in this paper). In terms of the Cartesian electronic variables, it is given by

and is obtained^2,6 from Eq. (4) by adding and subtracting the average of the diagonal electronic matrix elements, $H̄(R)≡1F∑kFHkk(R)$⁠. Note that in this symmetrized version, the ZPE γ-parameter of Eq. (4) does not appear, however it still enters the calculation through the initial and final actions using Eq. (3a). Eq. (5) also makes it apparent that the local frequencies of the electronic oscillators (i.e., at a given R) are the differences of the diabatic PES, which is essentially why the vibrations of these electronic oscillators couple effectively to the nuclear DOFs.

To determine the population of individual electronic states (initially and at later times) within a simple classical framework, we have employed an updated version of the quasi-classical (QC) model that was used many years ago for determining state-to-state vibrational transitions from classical trajectory simulations. The underlying idea here is that quantum states correspond to integer values of the classical action variables {n_k} (i.e., Bohr-Sommerfeld or semiclassical quantization), but rather than imposing this quantization condition strictly it is replaced by requiring the action variables to be within a “window” about the relevant integer values of the actions. This is not only easier to implement than requiring the actions to be exactly integer values, but also the smoothing effect this brings about leads to classical results that are in much better agreement with the correct quantum results.¹⁸ (There are a variety of ways to derive this SQC model, as reviewed in Ref. 6.) The symmetrical aspect of the SQC procedure is that the quantization condition (the “windowing”) be imposed equivalently on both the initial and final actions of each computed trajectory. This builds in microscopic reversibility (unlike the standard QC prescription) and, as will be illustrated below, the symmetrical aspect is critical for properly handling the limit of weak non-adiabatic coupling.

In all previous work, we have applied the SQC quantization procedure using histogram window functions, e.g.,

where

which is basically a box of width 2γ centered at the QM value of the classical action, N = 1 or 0 (for occupied or unoccupied, respectively). In the limit γ → 0, these window functions become delta functions, δ(n − N), which is the exact limit of Bohr-Sommerfeld quantization but which, as noted above, is not what is desired. These window functions may thus be thought of as “pre-limit” delta functions and are thus the essential parameters of the SQC model. There is, of course, no unique pre-limit delta function, and we have to date applied Occam’s razor in choosing the simplest model that makes sense, realizing that the window functions for the different states should not overlap and that each window should be centered on the quantum values of the action variables for its corresponding state (e.g., as in Fig. 1).

In addition, we have found that much better results are obtained if the width parameter γ in Eq. (6) has the same value as the zero point energy (ZPE) parameter γ in Eq. (4) (see also Eqs. (2) and (3)), and that this single γ-parameter be chosen to have the value $(3−1)/2≈0.366$⁠. In practice, we have found this particular value of γ to be optimal (or nearly so) for every application of the SQC/MM methodology that we have reported, and this choice also has some theoretical justification, as described in the Appendix of Ref. 2. Nevertheless, we have always been careful to point out that γ is essentially still just an empirical parameter (QM’ly it should be 1/2), and that other values of γ may work better in future applications or in certain regimes yet to be explored (as will be demonstrated below).

The histogram window function of Eq. (6) “quantizes” a single electronic DOF (constraining a single action n_k to lie within an interval $N−γ,N+γ$⁠, N = 0 or 1). To window all the actions associated with a full electronic configuration, one may construct a joint window function by taking the product of 1-D window functions. E.g., for the case of F = 2 electronic states, the complete configuration of actions for state 1 is (n₁, n₂) = (1, 0), and likewise, for state 2 it is (n₁, n₂) = (0, 1); therefore, joint window functions based on Eq. (6) may be written as

Likewise, for an arbitrary F number of electronic states/DOFs, a joint window function for the full configuration of actions corresponding to a state k may be written as

Thus, in the SQC/MM approach, a non-adiabatic transition probability between an initial electronic state i and final state f is calculated by Monte Carlo evaluation of the following phase space integral over the initial values of the coordinates and momenta of all DOFs (nuclear + electronic), which is carried out simultaneously for all energetically feasible final states k,

In evaluating Eq. (9), the window function for the initial state, W_i(n), is used to sample the F initial actions; likewise, ρ(P, R) is used as the sampling function for the G nuclear DOF. The full set of action window functions {W_k(n)}, corresponding to all possible final states, is used to “bin” the time-evolved actions n(t) at each desired final time t over the ensemble of trajectories. This gives a set of raw transition probabilities ${P̃k←i(t)}$⁠, one for each possible final state k, which are renormalized via

to give the final SQC/MM result for the transition probability.

While the symmetrical windowing is done in terms of the electronic actions, trajectories are usually run in Cartesian oscillator variables using the symmetrized Cartesian version of the MM Hamiltonian given in Eq. (5). Thus, after selecting initial actions using W_i(n), the corresponding Cartesian coordinates and momenta are calculated from Eq. (2), where the q_k’s—the classical “angle” variables conjugate to the electronic actions—are chosen randomly from the interval [0, 2π]. At each final time t, Eq. (3a) is used to calculate the set of final actions associated with each trajectory, n(t), which are then “quantized” in Eq. (9) with {W_k}.

The primary goals of this paper are (i) to apply the SQC/MM approach to the problem of very weak electronically non-adiabatic coupling, (ii) to gauge the extent to which our standard SQC/MM model fails in the weak-coupling regime (as we have suspected it will), and (iii) to present a new version of the SQC approach which expands the model to treat this regime.

To begin, we apply our standard SQC model (to date) to a universal 2-channel time-independent scattering model that was considered and used as a benchmark in the original MM paper.⁸ 

A general 2 × 2 unitary matrix can be written in terms of three parameters, e.g., two angles and one scalar, as

which we take to be the 2 × 2 S-matrix for a 2-channel scattering process (i.e., two electronic states with a given time-dependent trajectory for the nuclei), with the final electronic amplitudes related to their initial values by

The parameter p in Eq. (11) thus sets the transition probability between the two electronic states, i.e., $S1,22=p$⁠. Writing the initial and final electronic amplitudes in terms of classical action-angle variables (as done by MM),

gives the final electronic action variables in terms of the initial action and angle variables as

This expression may then be used to calculate the P_2←1 transition probability within the SQC framework of Eq. (9), though here there are no nuclear DOFs (P, R) to consider. The phase space average is thus only over the initial actions and angles of the electronic DOFs (n₀ q₀), with the initial actions sampled using the state 1 window function of Eq. (7). The final actions are “binned” in the usual way (using Eq. (7)), but here are calculated directly from the initial actions and angles using Eq. (14) (i.e., there are no trajectories to integrate in order to determine the final variables as a function of their initial values).

This 2-channel scattering problem is universal in the sense that it encompasses all possible electronic transition probabilities (as specified by the parameter p in Eq. (11)). The limitation is that it assumes a given nuclear trajectory R(t) (i.e., the “classical path approximation”) restricting its scope to cases where the details of the nuclear dynamics do not affect the overall electronic transition probabilities—e.g., for high collisional velocities. We are not concerned with such details here, however, our focus being on using this universal model problem as a tool for assessing the performance of the SQC/MM methodology for all possible transition probabilities simultaneously, and most importantly, for analyzing the regime of exponentially weak non-adiabatic coupling.

Note also that in Eq. (14) the S-matrix phase parameters α and β enter the result only as an additive constant to the initial angle variables q₁ and q₂. However, since q₁ and q₂ are each averaged uniformly over the interval [0, 2π], α and β do not affect the final results and can thus be ignored; i.e., the only parameter of this universal 2-channel model embodied in Eq. (11) is p which parametrizes the transition probability.

Fig. 2 shows the results obtained by applying the standard SQC windowing methodology to this universal model problem. Specifically, plotted in Fig. 2 are the P_2←1 transition probabilities calculated using different choices of the γ-parameter over a range of possible values of S-matrix parameter p. (Parameter p can, of course, vary from 0 to 1, but by symmetry, P_2←1 over $p∈[0,12]$ is equivalent to 1 − P_2←1 over $p∈[1,12]$⁠.)

The left and right panels of Fig. 2 show the same results, but on the right the results are plotted using a log-log scale. The exact result is that P_2←1 = p, as indicated by the diagonal line in both panels. It is seen that for no value of the γ-parameter does the SQC methodology give good agreement with the exact results over the full range of p. The main conclusions from the left panel are (i) that small values of γ (⁠$<14$⁠) perform quite poorly (in addition to being less efficient due to the more restrictive windowing) and (ii) that our preferred value of γ = 0.366 (⁠$≈ (3−1)/2$⁠) does in fact appear, on average, to be about the optimal value.

The log-log scale of the right panel provides a detailed view of the weak-coupling (i.e., small p) limit. As we anticipated based on Fig. 1, if the SQC window functions do not touch (i.e, $γ<12$⁠) the SQC methodology fails to properly treat the weak-coupling limit. Fig. 2 illustrates that even for values of p ≈ 0.1 (i.e., not too small), the SQC results begin to show significant deviation from the exact results, and by p ≈ 0.01 the SQC results have essentially fallen to zero.

The exception is the result corresponding to $γ=12$⁠. Although the $γ=12$ windows do a relatively poor job for values of p > 0.1, as seen in the left panel of Fig. 2, it is quite surprising to see how well they perform for values of p < 0.1, as seen in the right hand panel, even down to p ≈ 10⁻⁵. The reason for this is that the $γ=12$ windows touch at a point in the electronic action space thereby effectively probing the perturbative limit of the non-adiabatic process. The notion that the windows must touch to describe this limit properly also underscores another advantage of the SQC approach versus the traditional QC method: in the latter, because trajectories are always initialized with integer action, they can never be perturbatively close to another (weakly coupled) state.

The challenge is to find a SQC windowing scheme which provides the best of both worlds—i.e., performs as well as using γ = 0.366 in the “normal” coupling limit, but as well as using $γ=12$ in the weak coupling limit.

We have indeed discovered a simple alternative SQC windowing scheme which properly treats the limit of exponentially weak coupling, and as a bonus, it performs even better than our standard methodology in the “normal” coupling regime while being no more complicated to implement.

The solution is to use triangular-shaped window functions as shown in Fig. 3. This may at first seem strange, but the logic leading to this development is straightforward. (Also, keep in mind that these window functions are actually pre-limit delta functions about the integer action values which define the quantum states, and there is no unique pre-limit delta function; this is the modelistic aspect of the SQC approach.) From the analysis in Section III of the data in Fig. 2, it is clear that a set of improved SQC window functions should satisfy several criteria:

• The window functions corresponding to the different states must touch at a point (the key feature of the $γ=12$ square histograms).

• A value of γ in the neighborhood of 0.366 should be retained, as experience has repeatedly shown that this value is approximately optimal (aside from the weak-coupling limit).

• The lower/outer edge of each window function should be matched/tied to the ZPE, i.e., be defined by n > − γ (this is the remedy for the problem of ZPE leakage).

• The centroid of the window functions should be the integer action values which define the quantum states, i.e., (1, 0) and (0, 1) (because the SQC approach is to be viewed as an approximation to semiclassical Bohr-Sommerfeld quantization).

It is criteria (iii), and what is missing from it, which gives the important clue: there does not appear to be any reason why the upper limits of n₁ and n₂ imposed by each window need to be tied to the ZPE/γ-parameter. Hence, a logical tactic is to take the original two square histogram windows (of γ ≈ 0.366) and to extend their inner edges towards each other until they touch (specifically, by relaxing the upper limit of n₂ along the inner n₁ = 1 − γ edge of the state 1 window, and by relaxing the upper limit of n₁ along the inner n₂ = 1 − γ edge of the state 2 window such that both edges extend to and touch at the point (1 − γ, 1 − γ)), thus satisfying criteria (i).

Following this logic further, the corresponding outer limits of each window must be adjusted so that the windows are balanced, and then—in an attempt to preserve the value of γ we know to be effective (criteria (ii))—the upper portion of the window functions is cutoff with a diagonal line which extends through the point (1 − γ, 1 − γ) where the windows touch.

Equations for these triangle window functions (analogous to Eq. (7) for square histograms) are thus given (with normalization) by

where, again,

Note that the h(2 − 2γ − n₁ − n₂) factor in each of Eqs. (15a) and (15b) serves as an upper bound to the actions and makes the triangular windows non-separable in the electronic DOFs. However, one notes that separability is not a requirement for a multi-dimensional pre-limit delta function, and so Eqs. (15a) and (15b) may still be viewed as approximately implementing Bohr-Sommerfeld quantization.

The diagonal line, n₁ + n₂ = 2 − 2γ, used in Eq. (15) is a polyad of constant total n, and because the sum of the n’s (n₁ + n₂) is conserved by the MM Hamiltonian, all trajectories which start on this line move along this line, and trajectories which start away from it nevertheless move parallel to it. In other words, one could say that the triangular shapes of these window functions “fit” the dynamics of the electronic action space better than the square histogram windows.

A more subtle point, and potentially a practical consideration from the standpoint of numerical convergence, is that the extension of the outer edges of these windows (in conjunction with the elimination of trajectories having the highest range of total n) actually patches a gap which was present in the scheme provided by the square histogram windows. Specifically, the ensemble of trajectories—each of which moves diagonally in the action space on a polyad of constant total n—in the most extreme case, has its turning points (in the action space) at the lines given by n₁ = − γ and n₂ = − γ, which are exactly the outer edges of these triangle windows. As a consequence, these new windows provide a more complete solution to the problem of ZPE leakage than the square-histogram windows by virtue of their triangular shape. It is also pleasing to note (although perhaps of little consequence) that the triangle windows, in a symmetrical fashion, exactly divide the dynamically accessible electronic action space into three regions: the regions of the two states and a triangular-shaped intermediate region.

The remaining considerations are criteria (ii): the selection of a value for the γ-parameter close to 0.366; and, relatedly, the satisfaction of criteria (iv): that the centroids of the window functions be the integer action values which define the Bohr-Sommerfeld quantum states. Fig. 3 shows a set of triangle (and square histogram) window functions having $γ=13$ (which is reasonably close to 0.366). It turns out that this is the unique value of γ for which the triangle windows simultaneously satisfy both criteria (i), that the windows touch at a point, and criteria (iv), that the windows are properly centered at the quantum values.

To see this, note that triangle windows for states 1 and 2, as defined in Eq. (15), always touch at the single point (1 − γ, 1 − γ). Thus, rather than altering the width of these windows, varying γ instead has the effect of translating both windows in concert along a 45° line (n₂ = n₁) extending from the origin. For example, Fig. 4 shows the triangle windows (as defined by Eq. (15)) for extreme values of the γ-parameter, γ = 0 and $γ=12$⁠. For γ = 0, the windows are just above the n₁ + n₂ = 1 polyad (line of constant total n) which passes through the quantum values. The $γ=12$ windows border this same polyad from below. For the intermediate value of $γ=13$⁠, the windows are correctly positioned, centered on the quantum values of the actions, as shown in Fig. 3, and this is simply because the centroid of an isosceles right triangle, short sides of unit length, is exactly $13$⁠. Thus, unlike the square histogram scheme where γ was viewed as an empirical parameter whose optimal value was arrived at through numerical “experimentation” (and separately justified on independent theoretical grounds²), for these triangle windows, the “right” value of γ (i.e., $13$⁠) is uniquely determined from the geometric considerations required to satisfy the four criteria listed above.

The triangle window functions have been designed to address the weak-coupling limit (by touching at a point), and to otherwise attempt to preserve the excellent results we have observed in the “normal” coupling regime (inter alia, by approximately retaining our preferred value for the γ-parameter). Fig. 5 shows the results of applying the SQC triangle windowing scheme to the universal model problem described in Section III, along with the prior results taken from Fig. 2 for square histogram windows of γ = 0.366 and of $γ=12$⁠. The log-log scale plot of Fig. 5 (right panel) confirms that the triangle scheme does in fact do an excellent job of reproducing the exact results in the limit of very weak coupling, all the way down to transition probabilities of 10⁻⁵. This is not too surprising: the triangle window scheme was constructed to address this limit, and it is of course true that the $γ=12$ square histogram windows (which also touch at a point) do an excellent job in this regime as well.

What is truly surprising about this new windowing scheme, however, are the results one obtains with it over the full range of “normal” coupling strengths. As shown in Fig. 5 (left panel), the $γ=12$ square histogram windows generally do a relatively poor job for anything other than the weak coupling regime, significantly worse than the γ = 0.366 square histogram results (also shown), which are generally quite good, representing a sort of “best fit” compromise (although the γ = 0.366 square histograms are not perfect either). The triangles, on the other hand, are seen in Fig. 5 (left panel) to reproduce the exact result essentially perfectly over the full range of possible coupling strengths (and this is without any alteration relative to the weak-coupling case, γ already being $13$⁠).

In summary, Fig. 5 shows that the triangle windows elegantly handle the limit of exponentially weak coupling (as they were designed to do) but as a bonus also essentially exactly reproduce the correct P_2←1 transition probabilities over all possible coupling strengths—i.e., for any 2-channel S-matrix within this simple universal model (where the dynamics of the nuclear DOFs do not enter).

While the universal model treated in Secs. III and IV has served its purpose as a useful discovery tool, the important question is how well the new SQC triangle window scheme performs in more realistic applications where the dynamics of the nuclear DOFs are intertwined with the electronic dynamics, and where the previous SQC square histogram model (with γ = 0.366) has been seen to provide reliable performance.

We begin by validating the new SQC triangle windowing model against a set of benchmark spin-boson problems that we have treated in recent papers^2,5 using our standard SQC approach with square histogram windows.

The spin-boson problem models an electronically non-adiabatic process in a condensed phase medium which acts to dissipate/absorb electronic energy and modulate the decoherence in the electronic DOFs. The model consists of a pair of diabatic potential energy surfaces (PESs) which are spatially offset multi-dimensional harmonic oscillators in the nuclear DOFs, energetically biased by 2ϵ, and coupled together by a non-adiabatic coupling constant Δ (which is independent of nuclear coordinates),

The SB problem comes in symmetric (ϵ = 0) and asymmetric versions (ϵ ≠ 0), the latter generally being considered far more challenging for simple (and inexpensive) methodologies to describe accurately. The most difficult aspect of the problem, of course, is accurately describing the time-dependence of the quantum coherence effects in the electronic DOFs (including the decoherence of such effects). The electronic coherences are most pronounced when the temperature T (⁠$=1kBβ$⁠) is low (relative to the coupling Δ), as determined by selecting an appropriate thermal distribution for the initial oscillator coordinates {Q_k} and momenta {P_k},

where $ak=2ωktanhβ ωk2$⁠.

The frequencies {ω_k} and coupling constants {c_k} are chosen according to some spectral density (SD) function,

which characterizes the distribution and coupling strengths of local vibrational modes in a specific condensed phase environment; the coupling constants in Eq. (16) are thus given in terms of the SD by

For simplicity, the SD is often taken to have a continuous functional form which cuts off at high-frequency, and for the benchmark problems considered here, the common exponentially damped Ohmic form is used given by

where ω_c is the SD’s characteristic frequency and α is the bath coupling (or friction) parameter.

The same four cases are considered here as in Refs. 2 and 4—symmetric and asymmetric versions of the SB problem at high and low temperatures. Results using standard square histograms and the triangle windows are shown in Fig. 6, with the parameters given in the figure caption for the four cases. Specifically, the results correspond to the time-dependent population difference between electronic states 1 and 2,

after the system is initialized in electronic state 1 with the appropriate thermal distribution for the nuclear DOFs (as given in Eq. (17)). SQC results are plotted against the exact QM benchmark results from Ref. 19 for the symmetric SB problems and from Ref. 20 for the asymmetric versions. All the SQC calculations employed G = 100 bath modes (see Eqs. (16) and (9)), so of a dimensionality on the order of what would be relevant to an electronic transition in the condensed phase.

These are not weak-coupling problems, and the prior results obtained using square histogram window functions were impressive to begin with. Fig. 6 basically serves to confirm that the new triangle scheme is able to maintain the good results previously seen with the SQC/MM approach.

There is one small point of improvement: at very short time in Fig. 6(a) there is a discernible delay in the initial population decay when calculated with the γ = 0.366 histogram window functions. This is again because these windows do not touch (and therefore a certain amount of time must pass before any trajectories have crossed the gap between the windows). In contrast, the results in Fig. 6(a) calculated using the SQC triangle windows show a smooth quadratic decay right from time zero. While it is only a small effect, it does remedy a slight qualitative discrepancy seen in our earlier results, and for some problems, this may turn out to be more pronounced.

Finally, with respect to classical trajectory simulations of QM processes, a great deal of interest is often expressed in the number of trajectories required to obtain converged results, as well as how few trajectories may be used to obtain reasonable results (even if not perfect). Therefore, provided in Fig. 15 of the  Appendix are results analogous to those shown in Fig. 6 but with only 1000 (1K) trajectories (versus a trajectory number on the order of 10⁴ used to obtain the essentially fully converged definite results shown in Fig. 6). The 1K results are not perfect, but they are indeed still reasonable and would likely be quite acceptable for modeling a “real” molecular system; and though these results have been presented in raw form in Fig. 15 (to show the level of statistical noise), one could always apply a simple smoothing procedure to improve their aesthetics if that is desired.

Additional SB results for cases we had not previously considered are shown in Fig. 7. The SQC results are benchmarked against exact QM results from Ref. 21 (as calculated using a combination of the quantum-classical path integral (QCPI), full path sum, and full “blip”²² sum methodologies, as detailed therein). These two SB examples are different from the prior examples in Fig. 6 in that the coupling to the bath is much stronger: α = 2 for the symmetric problem and α = 4 for the asymmetric version (ϵ = 5) (the other parameters as indicated in the figure caption).

The case presented in Fig. 7(b) is strongly asymmetric (with ϵ = 5), exhibiting monotonic decay down to the final state 1 population of ≈ 25%. It also has the stronger bath coupling (of α = 4). But the decay curve is simple, and both SQC methodologies (triangles and square histograms) reproduce the exact QM results quantitatively.

The symmetric case presented in Fig. 7(a) seems to be the more challenging of the two, the QM decay curve having a bit more structure to it, and the shorter time scale illustrating the quadratic decay from time zero. Here, there is once again apparent a slight qualitative discrepancy exhibited in the square histogram results at short time (as discussed above with respect to Fig. 6(a)). On the other hand, the SQC triangle results perfectly capture the initial quadratic decay, and also track the QM result with discernible improvement in the region of t ≈ 1 AU.

The reasonable conclusion to take away from these SB results—in combination with those from Sec. V A—is that the good performance of the triangle windowing scheme is not limited to the weak-coupling regime. Rather, the new scheme seems to do at least as well as the original square histogram methodology for realistic problems having higher coupling strengths, and in some cases, even slightly better. This is, of course, what one expects based on the analysis of the universal model in Sections III and IV.

With the new scheme suitably validated against previous benchmark problems (and some new ones), the next task is to apply the triangle windowing to more realistic weakly coupled problems where we know treatment with the square histogram windows must fail. A simple example having a single nuclear DOF interacting with two coupled electronic states is the well-known Tully #1 problem¹⁷ which we previously treated² with good results. Here, we revisit Tully 1, and use it to examine the limit of very weak coupling by modifying the non-adiabatic coupling relative to the standard version by many orders of magnitude.

Briefly, Tully 1 is a simple 1-D model of the avoided crossing of two potential energy surfaces (PES) which is given by

where, in the original version, A = 0.01, B = 1.6, C = 0.005, and D = 1 (all expressed in atomic units). Here we consider various exponentially reduced values of the C-parameter down to 10⁻⁵ thereby modifying the off-diagonal non-adiabatic coupling elements, H₁₂(R) and H₂₁(R), of the diabatic Hamiltonian matrix given in Eq. (21).

In its original form, Tully 1 is a very strongly coupled problem, and so we begin again by validating the new SQC triangle windowing scheme against prior square histogram results and a relevant QM benchmark calculation. Originally, we compared our standard SQC square histogram windowing methodology against QM benchmark results from time-domain wave-packet scattering calculations (as was done in Tully’s original paper¹⁷ for testing the fewest switches surface-hopping algorithm). Here, because the aim is to consider very weak transition probabilities, rigorous time-independent energy-domain scattering calculations were performed and then energy-averaged with Gaussian-weighting (with ΔP = 1) to produce results analogous to the original Gaussian wavepacket calculations. An example calculation is shown in Fig. 8 which corresponds to the original version of the Tully 1 problem (i.e., with parameter C = 0.005). The time-independent energy-domain scattering calculations in Fig. 8 exhibit pronounced “Feshbach resonances” which correspond to a weakly bound state on the upper adiabatic surface. These resonances, however, are smoothed away by the Gaussian averaging, as also shown in Fig. 8, which gives us a benchmark result very similar (if not indistinguishable) to what one obtains from an analogous wavepacket calculation with a similar spread in momenta.

Fig. 9 shows SQC/MM results using both the triangle and square histogram windowing schemes versus the exact (Gaussian-weighted) QM benchmark result taken from Fig. 8. The SQC square histogram result is a recomputed version of what we have reported previously² and it agrees very well with the QM result. Results shown as computed with the triangle windowing scheme are also quite good, with some slight differences. The triangle windowing does seem to over-smooth the threshold region to a slightly greater extent than the square histogram windows do (probably because the latter are somewhat more localized in the action space as shown in Fig. 3). Above threshold, however, while the square histogram results are good, the triangle scheme gives results which are in even better agreement with the QM transmission probabilities, over the entire momentum range from 9 to 29 AU. (Note that, for consistency, the same Gaussian-weighted averaging has been applied to the SQC/MM results, though it has little effect as the SQC/MM model already has some intrinsic energy averaging due to the dispersion in electronic energy from the windowing.)

Our real interest here, though, is the weak coupling regime. Fig. 10 shows a calculation analogous to Fig. 9 (same Gaussian-weighted averaging, etc.) but with the coupling constant C in Eq. (21) reduced by a factor of 50 to 10⁻⁴, and the resulting T_2←1 transmission probabilities reduced by more than two orders of magnitude. Unsurprisingly, the standard histogram scheme captures almost none of the non-adiabatic transition probability, while the triangle scheme, though not perfect, does capture a good portion of the non-adiabatic transition probability throughout the main portion of the sharply peaked region (momenta ≈ 7 AU) and slightly beyond (≈ 8 AU) where there is significant structure in the raw QM results; they also describe the threshold region (i.e., 4–7 AU) quite well. More impressive are results generated with the triangle methodology in the region of higher collisional momenta (from 9–29 AU) where the agreement with the QM calculations (throughout this range) is essentially exact. (In this region, as seen in Fig. 8, there are no QM resonance-related phenomena in the nuclear DOFs likely complicating matters.)

Finally, similar to what was done with the universal model, a broad range of non-adiabatic coupling strengths was investigated by varying parameter C (Eq. (21)) —from its original value of 0.005 down to 10⁻⁵—while fixing as constant the collisional momentum (P = 15 AU, i.e., a value in the middle of Figs. 9 and 10 was chosen). Fig. 11 shows the SQC-calculated T_2←1 transmission probabilities plotted in log-log scale versus the exact QM results and what is seen is reminiscent of what was seen with the universal model in Fig. 5: the original/standard SQC square histogram scheme begins to show some deviation from the exact results for non-adiabatic transition probabilities around 0.1, significant deviation around 0.01, and below that, the original methodology is unable to capture any of the non-adiabatic effects. Once again, however, the new triangle scheme is seen to capture²³ these effects, seemingly with near perfect quantitative accuracy.

In recent work,⁶ we used our standard (square histogram) SQC/MM model to treat a suite of 2-state site-exciton (SE) models (as a preamble to treating a far more complicated 7-state SE model of the light-harvesting Fenna-Mathews-Olson (FMO) complex). Here, to give an example of weak electronic coupling in the condensed phase, we revisit one of these 2-site models, but reduce the electronic coupling from Δ = 100 cm⁻¹ to 10 cm⁻¹, and then to 1 cm⁻¹ (i.e., to 10% and then to 1% of its original value).

The SE Hamiltonian²⁴ is very similar to the SB Hamiltonian given above in Eq. (16), the primary differences being that, in the SE model, there are different sites each of which may be electronically excited, and that each site is coupled to its own independent harmonic bath. Considering the set of electronic states where only one site is electronically excited at a time (the notion being that the electronic excitation hops from site to site), the diabatic Hamiltonian matrix is given by

where {ϵ_k} are the site energies, Δ is the non-adiabatic coupling constant, and

In the standard SE model, the bath vibrational modes {ω_ξ} are characterized by a Debye spectral density (SD),

ω_c being the characteristic frequency and λ the “reorganization energy,”²⁵ which is used (like Eq. (20) above) to select the bath frequencies and coupling strengths. Specifically, the coupling constants D ≡ {D_ξ} in Eq. (22) are given by

which is almost the same as Eq. (19) above for the SB problem, the difference being one of form and consistency with Eq. (22). See Refs. 6 and 24 for further details. As in Ref. 6, calculations here employed G = 200 bath modes per site (see Eq. (23)) and thus 400 nuclear DOFs.

In Ref. 6 we treated 8 separate parameter regimes for this problem. Here, to consider weak-coupling, we treat as an example the case which exhibited the strongest electronic coherence (because drastically reducing the coupling strength will reduce the coherence).

Fig. 12 shows the results corresponding to the original problem with non-adiabatic coupling Δ = 100 cm⁻¹, and for exponentially reduced values of the coupling, Δ = 10 cm⁻¹ and Δ = 1 cm⁻¹ (with the other parameters as given in the figure caption). To show treatment of the weak-coupling limit, the figure plots the population transfer to the ground state P_2←1 on a log-scale. SQC/MM results using the original square histogram and triangle SQC windowing schemes are shown versus results calculated via Ishizaki and Fleming’s hierarchical equations of motion (HEOM) methodology, which is considered essentially exact for this site-exciton problem.

The results in Fig. 12 indicate that what we have learned from the universal model in Section III (and specifically, Fig. 5) likely holds for non-trivial electronically non-adiabatic processes involving hundreds of nuclear DOFs. For the original version of the problem (Δ = 100 cm⁻¹)—which is in the “normal” coupling regime (i.e., transition probabilities ≳ 0.1)—both the triangle and square histograms show almost perfect agreement with the benchmark HEOM results (the triangle scheme perhaps doing a slightly better job of gauging the magnitude of the first recurrence in the transition probability).

However, when the non-adiabatic coupling is reduced to 10% of its original value (Δ = 10 cm⁻¹), the HEOM benchmark shows that the true P_2←1 transition probabilities never reach 0.1. In this regime, the square histogram result begins to show significant deviation from the exact answer—just the same behavior as seen in Fig. 5—while the SQC triangle calculation is still able to reproduce the true dynamics.

Moreover, the conclusion generalizes to the weaker coupling strength of Δ = 1 cm⁻¹ (i.e., 1% of its original value): in this very weak-coupling regime, the true (HEOM-calculated) P_2←1 transition probability never exceeds 10⁻³, and here, as expected, the square histogram result is exactly zero (and so on this log-scale plot does not even appear). In contrast, the SQC triangle scheme is still able to retain its original accuracy—again, precisely what was seen in Fig. 5 for the universal model system (and also with Tully 1 in Fig. 11).

Therefore, the qualitative behavior exhibited by the universal (classical path) model treated in Section III is seen to hold for the more realistic benchmark systems treated in this section, and the new SQC triangle window scheme is demonstrated capable of describing these more complex systems quite well in both the normal coupling regime and the weak-coupling limit.

We recently showed⁷ how the SQC approach may be viewed as an approximate implementation of a new Wigner model where the Wigner functions are calculated directly (using the Wigner transform) in action-angle variables (instead of being calculated in Cartesian variables and then converted to action-angle variables, as is traditionally done), and that the practical consequence of this view is that it provides an extremely simple prescription for calculating the full electronic density matrix (i.e., the off-diagonal electronic coherences as well as the electronic populations (diagonals)) within the standard SQC framework, and to do so at essentially no additional computational expense.

As a concrete example, Ref. 7 demonstrated one way this may be done using square histogram window functions and treated the same SE model above (in Fig. 12) with Ishizaki and Fleming’s original electronic coupling of Δ = 100 cm⁻¹ (i.e., not weak coupling). Here, we provide an analogous demonstration (of calculating the full density matrix) for the same SE model except using the new triangle windowing scheme. Fig. 13 depicts a natural extension of the triangle scheme to the calculation of the full density matrix, where the red triangle window centered at $(n1,n2)=(12,12)$ is used to collect the electronic phase information encoded in the angle variables {q_i} as trajectories pass through it. (The trajectory calculation necessarily includes propagation of the angles variables anyway—but at time t they are simply discarded if one is only interested in the populations.)

The results of such a calculation to estimate the full 2 × 2 electronic density matrix for the SE model of Eq. (22) are shown in Fig. 14 plotted against exact HEOM results (as in Ref. 7). One can see that the off-diagonal electronic coherences (right panel) agree quite well with the exact results, though not perfectly. (The diagonal elements (left panel) are the same populations plotted in Fig. 12 (but not on a log scale).) We emphasize that this is a preliminary calculation based on what seems a natural extension of the triangle windowing scheme following the principles of Ref. 7, however, whether further optimizations are feasible will be a topic of future work covering a wider range of examples.

The primary accomplishment of this paper is the new SQC triangle windowing scheme described in Section IV. It was motivated by the goal of being able to treat very weak non-adiabatic transitions, and the critical aspect of being able to do so was revealed by considering the universal 2-channel scattering system (with a given nuclear path) in Section III. Section IV A discusses the intuitive insights that led to the final result. As a bonus, the new triangle windows are seen to describe not only the weak coupling regime of exponentially small transition probabilities extremely well, but also the “normal” regime (of moderate to strong non-adiabatic coupling) even better than our previous “standard” model (square histogram windows with the parameter γ = 0.366). At present we thus consider it to be our new “standard” version of the SQC/MM approach (at least for the case of 2 electronic states, as applications involving higher electronic dimensionality (more states) will be considered in later work).

For none of the benchmark systems treated herein (nor any that we have yet encountered) have we seen this new windowing scheme fail to work reasonably well (and often it is seen to work spectacularly well). Although we feel there is elegance in the simplicity of this new windowing geometry, there will surely be cases where it will fail to provide good results. The spirit of this work is to push forward the limits of what one can accomplish (in the treatment of non-adiabatic processes) within the framework of a straightforward classical MD approach, but in a theoretically sound and intuitive way, and as universally as possible.

The formal origin of the window functions in these SQC approaches is as “pre-limit” delta functions that impose an approximate version of Bohr-Sommerfeld (i.e., semiclassical) quantization (i.e., quantum states defined by integer values of the classical action variables). Though it has been realized for a long time that pre-limit delta functions introduce a smoothing effect that often makes a classical approximation work better than the actual delta functions themselves (and it certainly makes the calculations much simpler), there is nevertheless no unique pre-limit delta function, and this is where the modelistic aspect of the approach enters. One thus looks for models that are as simple and well-defined as possible, with as few empirical parameters as possible, and which cover as wide a range of systems as possible. We feel that this work represents significant progress along these lines, but it undoubtedly is not the end.

The authors thank Professor David Manolopoulos for providing source code for the Tully #1 benchmark QM scattering calculations, and Professor Nancy Makri for providing the benchmark QM path integral results shown in Fig. 7.

The authors also thank Mr. Herman Chan of the Whaley Research Group at Berkeley for providing the HEOM results for the site exciton problem shown in Figs. 12 and 14.

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 purpose of this work is to present the new SQC triangle windowing scheme and to provide an assessment of its performance, particularly in the limit of very weak electronic coupling. Therefore, the focus has been on definitive comparisons to benchmark results and no attempt has been made (in general) to minimize the number of trajectories used in the calculations presented here. Nevertheless (as mentioned in Section V), there is often interest in what can be obtained with relatively few trajectories for a given methodology. Thus, presented in Fig. 15 are results for the same four versions of the spin-boson (SB) problem presented in Fig. 6, but here running only 1000 (1K) trajectories, and one sees that quite reasonable results are obtained, although obviously not perfect. They show more statistical noise, but do track the true results with relatively small deviations at any fixed point in time (despite their oscillations); this would likely be adequate if the formulation of the model for a given molecular system is not that detailed and/or accurate to begin with. It is certainly sufficient for gauging the general character of the results and if, based on the noise, one desires additional confirmation, more trajectories can be added at a later time.

An additional detail regarding Fig. 15 versus Fig. 6 is that G = 50 (instead of 100) oscillators were used to discretize the harmonic bath of Eq. (23). This is clearly a sufficient number (based on Fig. 15), but it was observed that better 1K results were obtained using the more coarsely discretized bath, underscoring that it is the number of nuclear DOFs (each of which is thermalized) that dictates the number of trajectories required for reasonable convergence. Even with the smaller G = 50 discretization, the calculation in Fig. 15 still utilizes only 20 trajectories per nuclear DOF.