Maximilian A. C.
Saller
a,
Aaron
Kelly
b and
Jeremy O.
Richardson
a
aLaboratory of Physical Chemistry, ETH Zurich, Switzerland. E-mail: maximilian.saller@phys.chem.ethz.ch; jeremy.richardson@phys.chem.ethz.ch
bDepartment of Chemistry, Dalhousie University, Halifax, Canada. E-mail: aaron.kelly@dal.ca
First published on 24th September 2019
The mapping approach addresses the mismatch between the continuous nuclear phase space and discrete electronic states by creating an extended, fully continuous phase space using a set of harmonic oscillators to encode the populations and coherences of the electronic states. Existing quasiclassical dynamics methods based on mapping, such as the linearised semiclassical initial value representation (LSC-IVR) and Poisson bracket mapping equation (PBME) approaches, have been shown to fail in predicting the correct relaxation of electronic-state populations following an initial excitation. Here we generalise our recently published modification to the standard quasiclassical approximation for simulating quantum correlation functions. We show that the electronic-state population operator in any system can be exactly rewritten as a sum of a traceless operator and the identity operator. We show that by treating the latter at a quantum level instead of using the mapping approach, the accuracy of traditional quasiclassical dynamics methods can be drastically improved, without changes to their underlying equations of motion. We demonstrate this approach for the seven-state Frenkel-exciton model of the Fenna–Matthews–Olson light harvesting complex, showing that our modification significantly improves the accuracy of traditional mapping approaches when compared to numerically exact quantum results.
The development of simulation methods for nonadiabatic effects has thus continued to be the focus of considerable research efforts. Methods relying on an explicit expansion and propagation of the wavefunction, often on a grid, have yielded highly accurate results.4–6 However, many models inspired by condensed-phase systems still prove too computationally expensive to treat with these methods, due to their unfavourable exponential scaling with system size. Despite recent efforts to overcome this scaling hurdle,7,8 many systems from the fields of chemistry and biology, especially those in a condensed-phase environment, are simply too large to treat using a wavefunction-based approach. Mixed quantum-classical methods,9–29 though inherently more approximate, are often the only choice when seeking to simulate nonadiabatic dynamics in the condensed phase. As many of these scale linearly with system size, they can readily be applied to large and complex realistic systems, yielding highly valuable insights at reasonable computational costs.
The representation typically chosen for a nonadiabatic process consists of a continuous nuclear phase space and a set of discrete electronic states. The resulting Hamiltonian is given by
(1) |
(2) |
The mismatch between the continuous nuclear phase space and the discrete space of electronic states constitutes a recurring challenge in mixed quantum-classical dynamics. The mapping approach solves this problem by projecting the electronic degrees of freedom into a space of singly-excited harmonic oscillators (SEOs).12,30,31 In the space of the SEO wavefunctions, the representation of state |n〉 is given by
(3) |
(4) |
In addition to the relative simplicity of the mapping approach, the extended phase space grows linearly with the number of electronic states. Furthermore, given the favourable scaling of classical trajectories with respect to the nuclear DoFs, a number of mixed quantum-classical dynamics approaches, aimed specifically at large, realistic systems in the condensed phase, have been developed based on this formalism.10–12,14,24,26 Note that in this work we will use the term quasiclassical to refer to mixed quantum-classical approaches which employ a single set of mapping variables per electronic state as well as a single set of positions and momenta for each nuclear degree of freedom.
Quasiclassical methods yield accurate results for most observables at short times. In the long time limit however, they are well known to degrade in accuracy, especially for the relaxation to thermal equilibrium following an initial electronic excitation. Attempts to address this shortcoming with the use of master equations have shown considerable promise.32,33 Other approaches to improve quasiclassical dynamics have led to the development of related dynamics approaches. For instance, the symmetrical quasiclassical windowing method uses a windowing function to “bin” the electronic populations, ensuring that they have integer values at the beginning and end of each trajectory. This approach has been applied to the benchmark we study below, achieving accuracy comparable to that reported here.34–38 A number of methods which depart from the equations of motion underlying quasiclassical dynamics, but remain close to its overall motivation, have also shown considerable promise in treating multi-state systems. Prominent examples, which have been very successfully applied to the benchmark studied here, include the forward-backward trajectory solution19,20 (FBTS) and the partially linearised density matrix (PLDM) method.16,17
In a recent publication, we have however shown that a simple modification, with a similar motivation as that underlying the use of master equations, can drastically improve the performance of quasiclassical methods, without changing the equations of motion.39 We split the population operator into two parts, one of which is the identity.18 We can then use our understanding of the exact behaviour of this operator to drastically improve traditional quasiclassical methods. The resulting approach has the benefit of retaining all the advantages of these methods, as the underlying equations of motion are unchanged.
Here we extend this approach, which was originally presented for only two electronic states, to an arbitrarily large electronic space. We use the fact that the electronic population, like any Hermitian operator, can be expanded exactly into the identity and a purely traceless component.18 Given that the behaviour of its quantum operator is well understood, we treat the identity exactly, resulting in a simpler phase-space representation of the population operator, involving only the traceless part. Population dynamics calculated using these modified operators are of drastically higher quality than those obtained from the traditional quasiclassical definition.
We apply this general formulation to the challenging benchmark model for the Fenna–Matthews–Olson (FMO) light harvesting complex.40–44 Our results are significantly more accurate than those obtained using the standard operator definitions and in excellent agreement with numerically exact quantum dynamics methods.
(5) |
(6) |
When considering the population operator, there are two representations one can choose to Wigner transform, corresponding to either the left- or right-hand side of eqn (5). Using the left-hand side is equivalent to including a projection on the SEO subspace, which yields an expression in terms of the harmonic oscillator wavefunctions as in eqn (3).18 The resulting phase-space representations, identified as Awn and ASEOn respectively, are
(7a) |
(7b) |
(8) |
Note that crucially, Awn ≠ ASEOn × ϕ. Each of these phase-space representations is derived via the Wigner transform of a formally exact mapping form of the |n〉〈n| operator. Therefore, a clear choice of which to use when calculating observables is not obvious a priori. A more detailed discussion of the possible combinations of phase-space representations and electronic initial conditions can be found in our recent work.39
An observable commonly computed using quasiclassical methods is the population of a given electronic state, |n〉, given that the system was initially in a pure state, |m〉. In quantum mechanics this is defined by
(9) |
There is a unique expansion of the population operator |n〉〈n|, such that:
(10) |
(11) |
Note that in a two-level system, this operator is the Pauli spin matrix, i.e. 1 = z, such that |1〉〈1|=( + z)/2, which was used in our previous work.39 Substituting this definition for the quantum population operator into eqn (9) and expanding yields
(12) |
Following the standard quasiclassical procedure, in order to calculate the value of the population operator given in eqn (12), we Wigner transform the operators in these two constituent correlation functions. The phase-space representation of the traceless operator n is
(13) |
If we had performed the Wigner transform on the projected operator, the phase-space representation would simply be Qn(X,P)ϕ(X,P). Note that neither expression contains constant terms which play the role of ZPE-parameters.
It would be possible to arrive at a phase-space representation of the identity operator via similar Wigner transforms. However, we suggested in our previous work39 that we can instead use our understanding of its behaviour in quantum mechanics, which is to leave its operand unchanged. We therefore simply avoid directly computing the identity altogether.
Starting from the exact expression for in eqn (12), we thus arrive at our final quasiclassical expression for the population of electronic state |n〉, assuming the system was initially in state |m〉,
(14) |
The constituent correlation functions CQn and CQmQn are given by
CQn(t) = 〈ϕa(X,P)Qn(X(t),P(t))〉 | (15a) |
CQmQn(t) = 〈ϕa(X,P)Qm(X,P)Qn(X(t),P(t))〉, | (15b) |
Note that we can include the factors of ϕ(X,P) at time zero, because this function is constant over the course of any trajectory evolving under the Hamiltonian given in eqn (4). Also the two constituent correlation functions can be calculated for all values of m and n in a single simulation. Just as in traditional quasiclassical methods,11,14 the values of these constituent correlation functions, and therefore are exact in the limit of t = 0.
(16a) |
(16b) |
We note that the differences between these two methods do not actually stem from the derivations, but are mere convention. It would in principle be possible to derive a PBME method using two projections. However for convenience, we will use eqn (16) as the definition of PBME and LSC-IVR throughout this work. Finally, it is important to note that at least one of the operators has to be Wigner transformed in its projected form in order to ensure that the dynamics are initialised to the physical subspace.
While both LSC-IVR and PBME, as well as other mixed quantum-classical methods, have been applied to challenging systems with considerable success, their failure to accurately reproduce population dynamics in the long time limit has been well documented.18,39,46,47 As mentioned above, a number of modifications to quasiclassical methods which aim to address this issue have been proposed.34–39,48
In practice, both eqn (16) and (15) are evaluated by averaging over an ensemble of trajectories, propagated with Hamilton’s equations of motion defined by Initial conditions for each trajectory are sampled from ρwb(x,p) for the nuclei and ϕa(X,P) for the mapping variables.
The Frenkel-exciton model for the energy transfer in FMO is a challenging benchmark for quantum dynamics methods. In comparison to the Spin–Boson systems studied in our previous work,39 the FMO Hamiltonian presents a different kind of challenge to quasiclassical dynamics methods: the electronic subsystem is comprised of more electronic states and the system-bath coupling is different. We note that the key challenge resulting from a larger electronic state space is the possibility of reaction chains involving more than two states. As a result this benchmark and the FMO system in general has been extensively studied using a considerable number of approaches.16,20,36–38,48,50–60 In addition, though computationally challenging, numerically exact results are available, e.g. from hierarchical equations of motion (HEOM).42–44,61,62
In the seven-site model (S = 7), the full FMO Hamiltonian is given by
ĤFMO = Ĥs + Ĥsb + Ĥb, | (17) |
(18) |
(19) |
(20) |
(21) |
The frequencies, ωj(n), and coupling coefficients, cj(n), which are identical for each bath, are drawn from a spectral density of the Debye form, given by
(22) |
We define the initial nuclear density matrix b = e−βĤb/Zb, where the partition function Zb is defined such that the trace over bath modes only is Trb[b] = 1. Nuclear positions and momenta were sampled from the thermal Wigner distribution of the uncorrelated bath, given, for any bath, by
(23) |
Fig. 1 Constituent correlation functions of the FMO population, with T = 77 K, and τc = 50 fs and an initial excitation of site |1〉. The solid and dotted lines correspond to the electronic initial conditions having been sampled either from ϕ(X,P) or ϕ2(X,P). In eqn (15) this corresponds to a = 1 and a = 2 respectively. |
Considering the overall expression for the population given in eqn (14), the magnitudes of the constituent correlation functions are as one might expect. Notably the negative values observed for both CQn(t) and CQmQn(t) are not unphysical, as exact quantum mechanics would yield similar magnitudes for both correlation functions. We note that there is a noticeable difference between the correlation functions obtained from the two different initial distributions of the mapping variables we investigated. In order to assess their comparative accuracy however they must be combined, using eqn (14), into a population and compared to exact results.
Fig. 2 FMO site populations, with T = 77 K, and τc = 50 fs. Initial excitations of the |1〉 and |6〉 site are shown in the upper and lower panels respectively. Initial conditions for both constituent correlation functions sampled from ϕ(X,P). Results using our alternative population operator are shown as dashed lines, traditional PBME are dash-dotted while solid lines are the numerical HEOM benchmark.44 |
Our alternative definition of the population operators results in dynamics strikingly close in accuracy to the HEOM benchmark. We reproduce not only the correct ordering of states, even in the long time limit, but also capture all features present in the benchmark. We note that while we do observe negative populations for their magnitude is almost negligible. In addition, our values for are within the same margin of error of the exact HEOM result as every other population. This is especially encouraging when comparing the accuracy of our approach to that achieved with traditional PBME. The latter, while capturing the population dynamics at short times rather well, completely fails to reproduce the long time behaviour. Notably the distribution and ordering of states beyond the short time limit degrades drastically with this approach. Considering that the low-temperature parameter regime of the FMO Hamiltonian poses a considerable challenge to quasiclassical methods, the accuracy of our results is highly encouraging.
Fig. 3 shows populations for the same parameter regime of the FMO Hamiltonian as Fig. 2 however with electronic initial conditions now having been sampled from ϕ2(X,P). This corresponds to adding the dotted lines of Fig. 1 as per eqn (14). Also shown are standard LSC-IVR results, which again are a natural source of comparison as they use the same electronic initial conditions.
Fig. 3 FMO site populations, with T = 77 K, and τc = 50 fs. Initial conditions for both constituent correlation functions sampled from ϕ2(X,P). As in Fig. 2, results using our approach are shown as dashed lines and the numerical HEOM benchmark44 as solid lines. The dash-dotted results are the traditional LSC-IVR approach. |
Comparing to Fig. 2, LSC-IVR is seen to perform better than PBME for this particular parameter regime. Interestingly, the results achieved with our approach, using initial conditions sampled from ϕ2(X,P), do not give a considerable improvement over LSC-IVR. We thus conclude that the most consistently accurate method is that shown in Fig. 2 using our improved population operators and sampling from ϕ(X,P). This is the approach which we shall continue to employ for the remaining calculations.
Fig. 4 and 5 show the population dynamics of the FMO Hamiltonian, calculated with our alternative definition of the population operators and, as in Fig. 2, with electronic initial conditions sampled from ϕ(X,P), for two additional parameter regimes at T = 300 K. Traditional PBME is again shown for comparison along with numerically exact HEOM benchmark results.42–44,61,62 We note that the dynamics in the parameter regimes shown in these two figures are, owing to the higher temperature, less likely to be affected by nuclear quantum effects. Nevertheless it is clear that in the long-time limit, PBME diverges significantly from the HEOM benchmark and yields an incorrect distribution of states.
Fig. 4 FMO site populations, with T = 300 K, and a slower bath, τc = 166 fs. Initial conditions for both constituent correlation functions sampled from ϕ(X,P). Results are presented as in Fig. 2 and compared to numerical HEOM benchmark results.44 |
Fig. 5 FMO site populations, with T = 300 K, and a faster bath, τc = 166 fs. Initial conditions for both constituent correlation functions sampled from ϕ(X,P). Numerical HEOM benchmark results are taken from ref. 42. |
Using our alternative definition of the population operator again drastically improves the traditional quasiclassical result in both cases. Our approach in fact yields dynamics which now approach quantitative accuracy with respect to the exact HEOM benchmark. We furthermore note that the issue of small negative populations observed for our approach observed in Fig. 2 has now disappeared. This is not surprising, given that low temperature systems are well known to constitute a more considerable challenge for quasiclassical approaches. We recognise that quasiclassical approaches are well known not to capture nuclear effects, owing to the fact that the trajectories underlying them are driven by classical equations of motion. We note however that our alternative definition of the population operator is in fact not limited to quasiclassical methods, but may be applicable to other approaches based on mapping, which can capture nuclear quantum effects, such as nonadiabatic ring polymer molecular dynamics.21–23,27–29
Fig. 6 Long time FMO site populations, with T = 77 K, and τc = 50 fs. Initial conditions for both constituent correlation functions sampled from ϕ(X,P). Numerical HEOM benchmark results are taken from ref. 44. |
Comparing our results to the long-time limit of the HEOM benchmark,44 it is clear that our method does not rigorously preserve detailed balance. In addition, one of the populations is unphysically predicted to be slightly negative. We note however that as in Fig. 2 and 3, our negative result is within the same margin of error of the exact result as every other state population. It is well known that the long-time limits of the populations obtained from traditional quasiclassical approaches such as PBME and LSC-IVR can be much worse, predicting more strongly negative populations, and some greater than 1.50
We are therefore encouraged that the equilibrium distribution predicted using our approach is qualitatively similar to the HEOM benchmark. We note furthermore that our approach predicts identical equilibrium distributions, whether site |1〉 or site |6〉 is initially excited. We have not computed the long-time dynamics of the other two parameter regimes we studied above. We do however expect that, owing to their higher temperature and the lower impact nuclear quantum effects are therefore likely to have on them, our approach would perform even better than in the T = 77 K system. The fact that we do not observe negative populations in Fig. 4 and 5 further supports this hypothesis.
Overall we consider these population results, along with the others shown above, to be highly encouraging. They clearly demonstrate that our definition of the population operator can drastically improve the accuracy of traditional quasiclassical approaches11,14 at both intermediate and long times.
We have applied our approach to the challenging seven-state Frenkel-exciton model of the FMO light harvesting complex.40–44 In addition to having been studied extensively with traditional quasiclassical methods,50–52 the fact that numerically exact quantum results are available42–44,61,62 makes this system an ideal benchmark for our modification of the traditional quasiclassical population operators.
Overall we find that using our alternative definition of the electronic population operator drastically improves on the results obtainable with other quasiclassical methods. In addition, our results actually approach the exact quantum benchmark in accuracy for the three parameter regimes we study. Finally, we find that rather encouragingly, our method reproduces the long-time distribution of the electronic states with far higher accuracy than existing quasiclassical approaches.
We recognise that there have been other efforts to fix the well documented shortcomings of traditional quasiclassical dynamics methods. For instance, the symmetrical quasiclassical windowing approach uses “binning” to convert the population into an integer, using a windowing function applied symmetrically at the beginning and end of each classical trajectory.34–38 This approach has been applied to the FMO Hamiltonian we study here with considerable success, yielding results comparable in accuracy to those presented here.36–38
In recent work a post-processing method for the traditional LSC-IVR quasiclassical method has been proposed.48 The dynamics resulting from the traditional approach are shifted by a function which imposes the long-time Boltzmann distribution of the FMO subsystem Hamiltonian. We note that while the results obtained from this long-time correction do constitute an improvement over the traditional approach, they fail to address any inaccuracies at short to medium times. In addition, this approach relies on either having prior knowledge of the correct long-time distribution of states or approximating it.
Other mixed quantum-classical dynamics methods have also been applied to the FMO Hamiltonian with considerable success. In recent work both the forward-backward trajectory solution19,20 and the partially linearised density matrix16,17 approaches have yielded accurate results for the FMO systems studied here. We note however that these methods use a different set of equations of motion to propagate the classical trajectories from the quasiclassical approaches used here. This does however not disqualify them from also benefiting from our alternative definition of the population operator, as the latter is independent of the equations of motion.
While we have shown that our alternative definition of the electronic population operator can drastically improve on the dynamics obtained from traditional quasiclassical methods, it cannot address all their shortcomings. Notably, it cannot capture nuclear quantum effects, which are fundamentally inaccessible to an approach relying on purely classical trajectories to calculate operators. Due to its simplicity and generality however, our approach could be combined with methods which can capture some nuclear quantum effects such as tunnelling. The nonadiabatic ring polymer molecular dynamics21–23,27–29 method would seem to be a logical candidate for benefiting from our definition of the population operator, given that it is also based on the mapping formalism. We anticipate that the low temperature regime of the FMO Hamiltonian studied here may particularly benefit from such a combination, as nuclear quantum effects are likely to have a greater impact in this system.
Overall we have shown that using our alternative definition of the population operator results in a considerable improvement over traditional quasiclassical approaches. Our results in fact approach the accuracy of the numerically exact benchmark for the systems we have studied. This is highly encouraging, given that the traditional methods we compare to have previously been considered inadequate for the simulation of long-time nonadiabatic dynamics. Our modification does not involve changing the equations of motion underlying quasiclassical methods and in fact also scales identically. We therefore hope that this work will further the development of dynamics methods based on the quasiclassical approach, and the progress of mixed quantum-classical dynamics as a whole.
This journal is © The Royal Society of Chemistry 2020 |