Improved population operators for multi-state nonadiabatic dynamics with the mixed quantum-classical mapping approach

TL;DR

This paper introduces a modification to the mapping approach in nonadiabatic dynamics that improves the accuracy of population relaxation predictions by treating the identity component of the population operator quantum mechanically.

Contribution

The authors generalize a recent modification to quasiclassical methods, enhancing the accuracy of population dynamics in multi-state systems without altering existing equations.

Findings

Significant improvement in population relaxation predictions for a seven-state model

Enhanced accuracy of traditional mapping methods compared to exact quantum results

Applicable to complex light-harvesting systems like Fenna-Matthews-Olson

Abstract

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.

Equations

$$\hat{H}=\sum _{j=1}^F\frac{{\hat{p}}_j^2}{2m_j}+U(\hat{\mathsf{x}})+\hat{V}(\hat{\mathsf{x}}),$$

$$\hat{V}(\hat{\mathsf{x}})=\sum _{n,m}^S{V}_{nm}(\hat{\mathsf{x}})|n⟩⟨m|,$$

$$⟨\mathsf{X}⟩n=\frac{\sqrt{2}}{{\pi }^{S/4}}{X}_n\exp \left[-\frac{1}{2}\sum _{m=1}^S{X}_m^2\right],$$

$$\mathcal{H}=\sum _{j=1}^F\frac{p_j^2}{2m_j}+U(\mathsf{x})+\frac{1}{2}\sum _{n,m}^S\left(X_n{X}_m+P_n{P}_m-{\delta }_{nm}\right)V_{nm}(\mathsf{x}).$$

$$\Ket{n}\!\Bra{n}\equiv\hat{A}_{n}\mapsto\frac{1}{2}\bigg{(}\hat{X}_{n}^{2}+\hat{P}_{n}^{2}-1\bigg{)}\,,$$

$$O^w(\mathsf{x},\mathsf{p},\mathsf{X},\mathsf{P})=$$

$$\displaystyle\Braket{\mathsf{x}-\frac{\mathsf{y}}{2},\mathsf{X}-\vspace{1em}\frac{\mathsf{Y}}{2}}{\hat{O}}{\mathsf{x}+\frac{\mathsf{y}}{2},\mathsf{X}+\frac{\mathsf{Y}}{2}}\operatorname{\mathrm{d}}\mathsf{y}\,\operatorname{\mathrm{d}}\mathsf{Y}.$$

$$A_n^w(\mathsf{X},\mathsf{P})$$

$$A_n^{\mathrm{SEO}}(\mathsf{X},\mathsf{P})$$

$$\varphi (\mathsf{X},\mathsf{P})=2^{S+2}\exp \left[-\sum _{m=1}^S(X_m^2+P_m^2)\right].$$

$$P_{n\leftarrow m}(t)=\operatorname{Tr}\bigg{[}\hat{\rho}_{\text{b}}\Ket{m}\!\Bra{m}\mathrm{e}^{\operatorname{\mathrm{i}}\hat{H}t}\Ket{n}\!\Bra{n}\mathrm{e}^{-\operatorname{\mathrm{i}}\hat{H}t}\bigg{]}\,,$$

$$\Ket{n}\!\Bra{n}=\frac{1}{S}\bigg{(}\hat{\mathbb{I}}+\hat{Q}_{n}\bigg{)}\,,$$

$${\hat{Q}}_n=(S-1)|n⟩⟨n|-\sum _{m\ne n}^S|m⟩⟨m|.$$

$$\displaystyle P_{n\leftarrow m}(t)=\frac{1}{S^{2}}\bigg{(}S$$

$$\displaystyle+\operatorname{Tr}\left[\hat{\rho}_{\mathrm{b}}\hat{Q}_{m}\mathrm{e}^{\operatorname{\mathrm{i}}\hat{H}t}\hat{Q}_{n}\mathrm{e}^{-\operatorname{\mathrm{i}}\hat{H}t}\right]\bigg{)}\,,$$

$$Q_{n}(\mathsf{X},\mathsf{P})=\frac{1}{2}\bigg{[}(S-1)(X_{n}^{2}+P_{n}^{2})-\sum\limits_{m\neq n}^{S}(X_{m}^{2}+P_{m}^{2})\bigg{]}\,.$$

$$P_{n\leftarrow m}(t)\approx\frac{1}{S^{2}}\bigg{(}S+C_{\mathbb{I}Q_{n}}(t)+C_{Q_{m}Q_{n}}(t)\bigg{)}\,.$$

$$C_{\mathbb{I}{Q}_n}(t)$$

$$C_{Q_m{Q}_n}(t)$$

$$P_{n\leftarrow m}^{\text{PBME}}(t)$$

$$P_{n\leftarrow m}^{\text{LSCIVR}}(t)$$

$${\hat{H}}_{\mathrm{FMO}}={\hat{H}}_{\mathrm{s}}+{\hat{H}}_{\mathrm{sb}}+{\hat{H}}_{\mathrm{b}},$$

$${\hat{H}}_{\mathrm{s}}=\sum _{n=1}^S{\epsilon }_n|n⟩⟨n|+\sum _{n\ne m}^S{\Delta }_{nm}|n⟩⟨m|,$$

$$\!\!\!\!\!\!\!\!\!\!\!\!\hat{H}_{\mathrm{s}}\!=\!\left(\begin{array}[]{rrrrrrr}12410&-87.7&5.5&-5.9&6.7&-13.7&-9.9\\ -87.7&12530&30.8&8.2&0.7&11.8&4.3\\ 5.5&30.8&12210&-53.5&-2.2&-9.6&6.0\\ -5.9&8.2&-53.5&12320&-70.7&-17.0&-63.3\\ 6.7&0.7&-2.2&-70.7&12480&81.1&-1.3\\ -13.7&11.8&-9.6&-17.0&81.1&12630&39.7\\ -9.9&4.3&6.0&-63.3&-1.3&39.7&12440\end{array}\right)$$

$${\hat{H}}_{\mathrm{sb}}=-\sum _{n=1}^S|n⟩⟨n|\sum _{j=1}^B{c}_j^{(n)}{x}_j^{(n)},$$

$${\hat{H}}_{\mathrm{b}}=\frac{1}{2}\sum _{n=1}^S|n⟩⟨n|\sum _{j=1}^B\left(\frac{{\left(p_j^{(n)}\right)}^2}{m_j}+m_j{\left({\omega }_j^{(n)}{x}_j^{(n)}\right)}^2\right),$$

$$J(\omega )=2\lambda \frac{\omega {\omega }_{\text{c}}}{{\omega }^2+{\omega }_{\text{c}}^2},$$

$${\rho }_{\text{b}}^w(\mathsf{x},\mathsf{p})=\prod _{j=1}^B$$

$$\times \exp \left[-\tanh \left(\frac{1}{2}\beta {\omega }_j\right)\left(\frac{p_j^2}{{\omega }_j}+{\omega }_j{x}_j^2\right)\right].$$

Full text

I Introduction

Simulating nonadiabatic effects in quantum dynamics continues to pose a considerable challenge in theoretical chemistry and physics, especially in the condensed phase. Arising when the energies of two or more electronic states approach each other, resulting in the breakdown of the Born–Oppenheimer approximation, these effects have been found to have a profound impact on a wide range of systems spanning physics, chemistry and biology. Marcus (1993); Chandler (1998); Tully (2012)

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. Meyer et al. (1990, 2009); Richings et al. (2015) 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, Ben-Nun and Martínez (1998); Saller and Habershon (2017) 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, Tully (1990); Sun and Miller (1997); Sun et al. (1998); Stock and Thoss (2005); Bonella and Coker (2005); Kim et al. (2008); Ananth and Miller III (2010); Huo and Coker (2011, 2012); Kelly et al. (2012); Hsieh and Kapral (2012, 2013); Ananth (2013); Richardson and Thoss (2013); Richardson et al. (2017); Agostini et al. (2014); Makri (2015); Kapral (2015); Duke and Ananth (2015); Hele and Ananth (2016); Chowdhury and Huo (2017) 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

[TABLE]

where $\hat{p}_{j}$ and $m_{j}$ are the momentum operator and mass of nuclear degree of freedom (DoF) $j$ respectively, $\hat{\mathsf{x}}$ is a vector of length $F$ consisting of the position operators for each nuclear DoF. $U(\hat{\mathsf{x}})$ is the state-independent potential and the state-dependent potential is given by

[TABLE]

where $S$ is the number of electronic states. The diagonal elements of $\hat{V}(\hat{\mathsf{x}})$ are the diabatic potential energy surfaces, while its off-diagonal elements are the couplings between the electronic states. Everything that follows does not rely on a particular choice of potential, i.e. we are not limited to simple harmonic models. Note that we will use reduced units throughout, such that $\hbar=1$.

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 problem by projecting the electronic degrees of freedom into a space of singly-excited harmonic oscillators (SEOs).Meyer and Miller (1979); Stock and Thoss (1997, 2005) In the space of the SEO wavefunctions, the representation of state $\ket{n}$ is given by

[TABLE]

where $\mathsf{X}$ and its conjugate $\mathsf{P}$ are vectors of length $S$, corresponding to the position and momenta of the SEOs. The mapping variables $\set{\mathsf{X},\mathsf{P}}$ extend the nuclear phase space, $\set{\mathsf{x},\mathsf{p}}$. The resulting space, now completely continuous, can be used to propagate classical trajectories evolving under the mapping Hamiltonian, $\mathcal{H}$, given by

[TABLE]

In addition to the relative simplicity of the mapping approach, the extended phase space grows linearly with the number of electronic states. Given furthermore 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. Sun and Miller (1997); Sun et al. (1998); Stock and Thoss (2005); Kim et al. (2008); Agostini et al. (2014); Kapral (2015) 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 master equations have shown considerable promise.Kelly et al. (2015, 2016) 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, insuring 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. Cotton and Miller (2013, 2013); Miller and Cotton (2016); Cotton and Miller (2016, 2019) 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 solution Hsieh and Kapral (2012, 2013) (FBTS) and the partially linearised density matrix (PLDM) method. Huo and Coker (2011, 2012)

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. Saller et al. (2019) We split the population operator into two parts, one of which is the identity.Kelly et al. (2012) 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. Kelly et al. (2012) 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. Fenna and Matthews (1975); Adolphs and Renger (2006); Ishizaki and Fleming (2009); Sarovar et al. (2010) Our results are significantly more accurate than those obtained using the standard operator definitions and in excellent agreement with numerically exact quantum dynamics methods.

II Theory

Here we extend our previous workSaller et al. (2019) by presenting a general formalism, which can be applied to any system of multiple electronic states.

II.1 Quasiclassical population operators

In the mapping formalism, the operator $\Ket{n}\!\Bra{n}$, which measures the population of electronic state $\ket{n}$ can be written as

[TABLE]

where $X_{n}$ and $P_{n}$ are the mapping variables associated with state $\ket{n}$. In the quasiclassical approximation, the Wigner transform is used to define a phase-space representation for the operators of interest. The Wigner transform of a general operator, $\hat{O}$, is given by

[TABLE]

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 Eq. 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 Eq. 3.Kelly et al. (2012) The resulting phase-space representations, identified as $A_{n}^{\operatorname{\mathrm{w}}}$ and $A_{n}^{\operatorname{\mathrm{SEO}}}$ respectively, are

[TABLE]

where

[TABLE]

Note that crucially, $A_{n}^{\operatorname{\mathrm{w}}}\neq A_{n}^{\operatorname{\mathrm{SEO}}}\times\phi$. Each of these phase-space representations is derived via the Wigner transform of a formally exact mapping form of the $\Ket{n}\!\Bra{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. Saller et al. (2019)

An observable commonly computed using quasiclassical methods is the population of a given electronic state, $\ket{n}$, given that the system was initially in a pure state, $\ket{m}$. In quantum mechanics this is defined by

[TABLE]

where $\hat{\rho}_{\text{b}}$ is a density matrix which defines the initial state of the nuclei, normalised such that the trace over nuclear DoFs only is $\operatorname{Tr}_{\mathrm{b}}[\hat{\rho}_{\mathrm{b}}]=1$.

II.2 Traceless projection operators

There are two differences between the phase-space representations given in Eq. 7a and Eq. 7b: the factor of $\phi(\mathsf{X},\mathsf{P})$, which is only present in $A_{n}^{\operatorname{\mathrm{SEO}}}$, and the differing constant terms, which are related to zero-point energy (ZPE) of the mapping DoFs.Müller and Stock (1999) The origin of the latter is that both the projected and unprojected forms of $\Ket{n}\!\Bra{n}$ have a non-zero trace. We propose a form of the quantum population operator in which the trace is shifted to the identity operator, which in turn is treated exactly using quantum mechanics. Saller et al. (2019) The result is a phase-space representation of the quantum population operator which is traceless.

There is a unique expansion of the population operator $\Ket{n}\!\Bra{n}$, such that:

[TABLE]

where $\hat{\mathbb{I}}=\sum_{m=1}^{S}\Ket{m}\!\Bra{m}$ is the identity operator, and $\hat{Q}_{n}$ is, by design, traceless,

[TABLE]

Note that in a two-level system, this operator is the Pauli spin matrix, i.e. $\hat{Q}_{1}=\hat{\sigma}_{z}$, such that $\Ket{1}\!\Bra{1}=(\hat{\mathbb{I}}+\hat{\sigma}_{\mathrm{z}})/2$, which was used in our previous work.Saller et al. (2019) Substituting this definition for the quantum population operator into Eq. 9 and expanding yields

[TABLE]

where we have used $\operatorname{Tr}[\hat{\rho}_{\mathrm{b}}\hat{Q}_{m}]=0$ and $\operatorname{Tr}[\hat{\rho}_{\mathrm{b}}\hat{\mathbb{I}}]=S$. The final two terms in this expression are quantum correlation functions which can be approximated by well-known quasiclassical dynamics methods.

Following the standard quasiclassical procedure, in order to calculate the value of the population operator given in Eq. II.2, we Wigner transform the operators in these two constituent correlation functions. The phase-space representation of the traceless operator $\hat{Q}_{n}$ is

[TABLE]

If we had performed the Wigner transform on the projected operator, the phase-space representation would simply be $Q_{n}(X,P)\phi(X,P)$. Note that either expression contains no 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. We however suggested in our previous work,Saller et al. (2019) 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 $\hat{P}_{n\leftarrow m}(t)$ in Eq. II.2, we thus arrive at our final quasiclassical expression for the population of electronic state $\ket{n}$, assuming the system was initially in state $\ket{m}$,

[TABLE]

The constituent correlation functions $C_{\mathbb{I}Q_{n}}$ and $C_{Q_{m}Q_{n}}$ are given by

[TABLE]

where $\braket{\cdots}=\frac{1}{(2\pi)^{F+S}}\iiiint\rho_{b}^{\operatorname{\mathrm{w}}}(\mathsf{x},\mathsf{p})\cdots\operatorname{\mathrm{d}}\mathsf{x}\operatorname{\mathrm{d}}\mathsf{p}\operatorname{\mathrm{d}}\mathsf{X}\operatorname{\mathrm{d}}\mathsf{P}$ and $\rho_{b}^{\operatorname{\mathrm{w}}}(\mathsf{x},\mathsf{p})$ is the Wigner transformed density matrix of the nuclear DOFs. In practice the values of these correlation functions are averaged over an ensemble of trajectories, with initial conditions for the mapping variables being drawn from either $\phi(\mathsf{X},\mathsf{P})$ or $\phi^{2}(\mathsf{X},\mathsf{P})$, depending on whether the projected forms of one of both of $\hat{Q}_{m}$ and $\hat{Q}_{n}$ were Wigner transformed. This corresponds to $a=1$ and $a=2$ respectively.

Note that we can include the factors of $\phi(\mathsf{X},\mathsf{P})$ at time zero, because this function is constant over the course of any trajectory evolving under the Hamiltonian $\mathcal{H}$, given in Eq. 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,Sun et al. (1998); Kim et al. (2008) the values of these constituent correlation functions, and therefore $P_{n\leftarrow m}(t)$, are exact in the limit of $t=0$.

II.3 Traditional quasiclassical dynamics methods

The traditional quasiclassical approach does not involve treating the identity quantum mechanically as we have done above. There are two standard approaches which differ in whether both population operators are projected onto the subspace, or just one. These methods were derived in different ways Sun et al. (1998); Kim et al. (2008) and are called the Poisson bracket mapping equationKim et al. (2008); Kelly et al. (2012) (PBME) and the linearized semiclassical initial value representation Sun and Miller (1997); Sun et al. (1998) (LSC-IVR) methods. LSC-IVR commonly involves projecting both operators prior to Wigner transforming them, i.e. using $\Ket{m}\!\Bra{m}\mapsto A_{m}^{\operatorname{\mathrm{SEO}}}(\mathsf{X},\mathsf{P})$ and $\Ket{n}\!\Bra{n}\mapsto A_{n}^{\operatorname{\mathrm{SEO}}}(\mathsf{X},\mathsf{P})$. The Wigner transform of each operator yields, as per Eq. 7b, a factor of $\phi(\mathsf{X},\mathsf{P})$. Initial conditions for the mapping variables are therefore sampled from $\phi^{2}(\mathsf{X},\mathsf{P})$. In PBME on the other hand, traditionally only the operator for the initial population is Wigner transformed in its projected from. The operators are therefore $\Ket{m}\!\Bra{m}\mapsto A_{m}^{\operatorname{\mathrm{SEO}}}(\mathsf{X},\mathsf{P})$ and $\Ket{n}\!\Bra{n}\mapsto A_{n}^{\operatorname{\mathrm{w}}}(\mathsf{X},\mathsf{P})$. Only the transform of $\Ket{m}\!\Bra{m}$ yields a factor of $\phi(\mathsf{X},\mathsf{P})$. Consequently, electronic initial conditions are sampled from $\phi(\mathsf{X},\mathsf{P})$. Using these definitions, the electronic population can be calculated from

[TABLE]

We note that the differences between these two methods does not actually stem from the derivations, but is mere convention. It would in principle be possible to derive a PBME method using two projections. However for convenience, we will use Eqs. 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. Müller and Stock (1998); Bonella and Coker (2003); Kelly et al. (2012); Saller et al. (2019) As mentioned above, a number of modifications to quasiclassical methods which aim to address this issue have been proposed. Cotton and Miller (2013, 2013); Miller and Cotton (2016); Cotton and Miller (2016, 2019); Teh and Cheng (2017); Saller et al. (2019)

In practice, both Eqs. 16 and 15 are evaluated by averaging over an ensemble of trajectories, propagated with Hamilton’s equations of motion defined by $\mathcal{H}$. Initial conditions for each trajectory are sampled from $\rho_{\mathrm{b}}^{\operatorname{\mathrm{w}}}(\mathsf{x},\mathsf{p})$ for the nuclei and $\phi^{a}(\mathsf{X},\mathsf{P})$ for the mapping variables.

III Results and Discussion

III.1 The Fenna-Matthews-Olson Hamiltonian

The Fenna-Matthews-Olson complex is a pigment protein biomolecule found in green sulfur bacteria adapted for low-light environments. It consists of three identical trimers, each containing seven bacteriochlorophyll (BChl) pigments supported by a protein backbone. In photosynthesis the task of FMO is to transport the excitation gained from absorbing sunlight to the reaction centre where it is converted into electrochemical energy. Fenna and Matthews (1975); Li et al. (1997); Adolphs and Renger (2006); Ishizaki and Fleming (2009)

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,Saller et al. (2019) 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. Tao and Miller (2010); Huo and Coker (2011); Kelly and Rhee (2011); Kim et al. (2012); Hsieh and Kapral (2013); Miller and Cotton (2016); Cotton and Miller (2016, 2019); Teh and Cheng (2017); Scholes et al. (2007); Adolphs et al. (2007); Plenio and Huelga (2008); Caruso et al. (2010); Fassioli and Olaya-Castro (2010); Kolli et al. (2011); Olaya-Castro and Scholes (2011); Chin et al. (2013) In addition, though computationally challenging, numerically exact results are available, e.g. from hierarchical equations of motion (HEOM). Ishizaki and Fleming (2009); Sarovar et al. (2010); Ishizaki et al. (2010); Zhu et al. (2011)

In the seven-site model ($S=7$), the full FMO Hamiltonian is given by

[TABLE]

where $\hat{H}_{\mathrm{s}}$ is the electronic sub-system Hamiltonian, given by

[TABLE]

where $\varepsilon_{n}$ is the energy of BChl site $\ket{n}$ and $\Delta_{nm}$ is the electronic coupling between sites $\ket{n}$ and $\ket{m}$. The values of site energies and couplings used in the matrix representation of $\hat{H}_{\mathrm{s}}$ are given by

[TABLE]

all energies being in units of cm*-1*. The protein environment around every BChl site is modelled by a bath of harmonic oscillators. The system-bath Hamiltonian, $\hat{H}_{\mathrm{sb}}$, defines the coupling between the electronic sub-system and these baths. It is given by

[TABLE]

where $c_{j}^{(n)}$ is the vibronic coupling coefficient between site $\ket{n}$ and bath mode $j$. $B$ is the number of modes per bath, such that $B=S\times F$. The position coordinate of bath mode $j$ of the $n$th bath is $x_{j}^{(n)}$. Finally, the Hamiltonian for the baths, $\hat{H}_{\mathrm{b}}$, is given by

[TABLE]

where $p_{j}^{(n)}$ and $\omega^{(n)}_{j}$ are the momentum coordinate and frequency of bath mode $j$ associated with site $\ket{n}$. The choice of masses does not affect results, so one can effectively set $m_{j}=1$. Note that, following previous work, Ishizaki et al. (2010); Zhu et al. (2011); Kim et al. (2012); Cotton and Miller (2019) each BChl site is coupled to an identical bath, which in turn is uncoupled from all other baths. The coupling between sites is thus contained purely in $\hat{H}_{\mathrm{s}}$.

The frequencies, $\omega^{(n)}_{j}$, and coupling coefficients, $c_{j}^{(n)}$, which are therefore identical for each bath, are drawn from a spectral density of the Debye form, given by

[TABLE]

where $\omega_{\text{c}}$ is the characteristic frequency of the bath, related inversely to the phonon relaxation time, $\omega_{\text{c}}^{-1}=\tau_{\text{c}}$, and we use $\lambda=35\,\mathrm{cm}^{-1}$ throughout, following previous work. Ishizaki et al. (2010); Zhu et al. (2011); Kim et al. (2012); Cotton and Miller (2019) We discretize this function using a scheme known to reproduce exact reorganisation energies.Wang et al. (2001); Craig et al. (2007)

We define the initial nuclear density matrix $\rho_{\text{b}}=\mathrm{e}^{-\beta\hat{H}_{\text{b}}}/Z_{\text{b}}$, where the partition function $Z_{\text{b}}$ is defined such that the trace over bath modes only is $\operatorname{Tr}_{\text{b}}[\rho_{\text{b}}]=1$. Nuclear positions and momenta were sampled from the thermal Wigner distribution of the uncorrelated bath, given, for any bath, by

[TABLE]

III.2 Simulation parameters

In order to test our alternative definition of the quasiclassical population in Eqs. 15, we investigated three parameter regimes of the FMO Hamiltonian, which have been studied extensively, including using the numerically exact HEOM approach. Ishizaki and Fleming (2009); Sarovar et al. (2010); Ishizaki et al. (2010); Zhu et al. (2011) All our simulations used a timestep of $\delta t=1\text{fs}$, which was found to be numerically converged. The results presented here are averaged over an ensemble of $10^{6}$ trajectories in order to demonstrate the converged performance of our approach. We found however that using as few as $10^{3}$ trajectories was enough to qualitatively capture all significant features of the population dynamics and already exhibits the clear improvement over a fully converged traditional quasiclassical result. We note that in all our simulations, we used the traceless form of the $\hat{V}(\mathsf{x})$ matrix to propagate our trajectories, absorbing the remainder into $U(\mathsf{x})$.

III.3 Constituent Correlation Functions

Figure 1 shows the constituent correlation functions, $C_{\mathbb{I}Q_{n}}(t)$ and $C_{Q_{m}Q_{n}}(t)$, calculated with electronic initial conditions having been sampled from both $\phi(\mathsf{X},\mathsf{P})$ and $\phi^{2}(\mathsf{X},\mathsf{P})$.

Considering the overall expression for the population given in Eq. 14, the magnitudes of the constituent correlation functions are as one might expect. Notably the negative values observed for both $C_{\mathbb{I}Q_{n}}(t)$ and $C_{Q_{m}Q_{n}}(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 Eq. 14, into a population and compared to exact results.

III.4 Population dynamics

Figure 2 shows the population dynamics resulting from combining the constituent correlation functions using initial conditions sampled from $\phi(\mathsf{X},\mathsf{P})$, i.e. the solid lines in Figure 1. In addition to results obtained for an initial excitation of the $\ket{1}$ site, populations starting in site $\ket{6}$ are also shown. In both cases the traditional PBME populations, calculated as per Eq. 16a, are shown for comparison, along with numerically exact HEOM results. Ishizaki and Fleming (2009); Sarovar et al. (2010); Ishizaki et al. (2010); Zhu et al. (2011) The comparison to PBME is a natural one here, as, in practice, the electronic initial conditions of PBME are also sampled from $\phi(\mathsf{X},\mathsf{P})$, i.e. $a=1$. It is worth noting that this parameter regime of the FMO Hamiltonian is the most challenging for quasiclassical methods, due to the significant impact of quantum effects at low temperature.

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 $P_{6\leftarrow 1}(t)$, their magnitude is almost negligible. In addition, our values for $P_{6\leftarrow 1}(t)$ 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.

Figure 3 shows populations for the same parameter regime of the FMO Hamiltonian as Figure 2 however with electronic initial conditions now having been sampled from $\phi^{2}(\mathsf{X},\mathsf{P})$. This corresponds to adding the dotted lines of Figure 1 as per Eq. 14. Also shown are standard LSC-IVR results, which again are a natural source of comparison as they use the same electronic initial conditions.

Comparing these results those shown in Figure 2, where the mapping variables were sampled from $\phi(\mathsf{X},\mathsf{P})$, we observe a slight decrease in accuracy. Nevertheless, our approach retains all qualitative features of the dynamics and encouragingly yields the correct ordering of states throughout. LSC-IVR performs significantly better than PBME for this particular parameter regime, however still fails to yield accurate results beyond the short-time limit. We note in particular the incorrect ordering of states with respect to the HEOM benchmark. Our alternative definition of the population operator, although less accurate when sampling from $\phi^{2}(\mathsf{X},\mathsf{P})$, therefore still improves considerably on the LSC-IVR result, which is highly encouraging.

Figures 4 and 5 show the population dynamics of the FMO Hamiltonian, calculated with our alternative definition of the population operators and, as in Figure 2, with electronic initial conditions sampled from $\phi(\mathsf{X},\mathsf{P})$, for two additional parameter regimes at $T=300$ K. Traditional PBME is again shown for comparison along with numerically exact HEOM benchmark results. Ishizaki and Fleming (2009); Sarovar et al. (2010); Ishizaki et al. (2010); Zhu et al. (2011) 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.

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 Figure 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. Richardson and Thoss (2013); Ananth (2013); Duke and Ananth (2015); Hele and Ananth (2016); Richardson et al. (2017); Chowdhury and Huo (2017)

III.5 Populations in the Long Time Limit

One of the well known failings of quasiclassical dynamics methods is that they do not preserve detailed balance of the populations and are therefore inaccurate at long times. In order to investigate whether our alternative definition of the population operator can improve on this, we have carried out a longer simulation of the parameter regime shown in Figures 2 and 3. We used the same timestep, $\delta t=1\,\text{fs}$, as in the simulations above and averaged over the same number of trajectories ($10^{6}$). Figure 6 shows the FMO site populations, following an initial excitation of either the $\ket{1}$ or the $\ket{6}$ site, up to 10 ps, calculated with our definition of the population operator. Electronic initial conditions for both were sampled from $\phi(\mathsf{X},\mathsf{P})$.

It is clear that our method does not rigorously preserve detailed balance as one of the populations is unphysically predicted to be slightly negative. We note however that as in Figures 2 and 3, our negative result is within the same margin of error of the exact result as every other state population. It is however 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. Tao and Miller (2010)

On the figure, arrows indicate an approximation to the distribution obtained from the diagonal elements of the matrix exponential $\exp[-\beta\hat{H}_{\mathrm{s}}]/Z_{\mathrm{s}}$ where $Z_{\mathrm{s}}=\operatorname{Tr}[\exp(-\beta\hat{H}_{\mathrm{s}})]$. Note however that this is an approximation which neglects coupling to the bath modes. We are nevertheless encouraged that the equilibrium distribution predicted using our approach yields relatively similar results to this approximation. We note furthermore that our approach predicts identical equilibrium distributions, whether site $\ket{1}$ or site $\ket{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 thus the less impact nuclear quantum effects are likely to have in them, our approach would perform even better than in the $T=77$ K system. The fact that we do not observe negative populations in Figures 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 approachesSun et al. (1998); Kim et al. (2008) at both intermediate and long times.

IV Conclusion

We have outlined an extension of our previous work,Saller et al. (2019) presenting an alternative definition of the electronic population operator for any system of multiple electronic states. We rely on the fact that any Hermitian operator can be split into two terms, one of which is the identity, the other being traceless.Kelly et al. (2012) We then use our understanding of the exact behaviour of the quantum identity instead of a quasiclassical treatment. The combination of this splitting and exact treatment of the identity results in a new form of the electronic population operator. Our approach retains the excellent scaling with respect to system size of traditional quasiclassical methods as well as their underlying equations of motion. Notably, as the constituent correlation functions into which our new operator is split can be calculated for all states in a single simulation, our approach is no more computationally expensive than the traditional methods it seeks to improve.

We have applied our approach to the challenging seven-state Frenkel-Exciton model of the FMO light harvesting complex. Fenna and Matthews (1975); Adolphs and Renger (2006); Ishizaki and Fleming (2009); Sarovar et al. (2010) In addition to having been studied extensively with traditional quasiclassical methods, Tao and Miller (2010); Kelly and Rhee (2011); Kim et al. (2012) the fact that numerically exact quantum results are available Ishizaki and Fleming (2009); Sarovar et al. (2010); Ishizaki et al. (2010); Zhu et al. (2011) 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 continuous mapping variables into integers, using a windowing function applied symmetrically at the beginning and end of each classical trajectory. Cotton and Miller (2013, 2013); Miller and Cotton (2016); Cotton and Miller (2016, 2019) This approach has been applied to the FMO Hamiltonian we study here with considerable success, yielding results comparable in accuracy to those presented here.Miller and Cotton (2016); Cotton and Miller (2016, 2019)

In recent work a post-processing method for the traditional LSC-IVR quasiclassical method has been proposed.Teh and Cheng (2017) 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 FMO Hamiltonian with considerable success. In recent work both the forward-backward trajectory solutionHsieh and Kapral (2012, 2013) and the partially linearised density matrix Huo and Coker (2011, 2012) 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 dynamics Richardson and Thoss (2013); Ananth (2013); Duke and Ananth (2015); Hele and Ananth (2016); Richardson et al. (2017); Chowdhury and Huo (2017) 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 of the progress of mixed quantum-classical dynamics as a whole.

Acknowledgements

M. A. C. S. would like to acknowledge financial support via the ETH postdoctoral fellowship. The authors also acknowledge the support from the Swiss National Science Foundation through the NCCR MUST (Molecular Ultrafast Science and Technology) Network.

Notes and references

• Marcus (1993)

R. A. Marcus, Rev. Mod. Phys., 1993, 65, 599.

• Chandler (1998)

D. Chandler, Classical and Quantum Dynamics in Condensed Phase Simulations, World Scientific, Singapore, 1998, ch. 2, pp. 25–49.

• Tully (2012)

J. C. Tully, J. Chem. Phys., 2012, 137, 22A301.

• Meyer et al. (1990)

H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78.

• Meyer et al. (2009)

Multidimensional Quantum Dynamics: MCTDH Theory and Applications, ed. H.-D. Meyer, F. Gatti and G. W. Worth, Wiley-VCH, Weinheim, 2009.

• Richings et al. (2015)

G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt and B. Lasorne, Int. Rev. Phys. Chem., 2015, 34, 269–308.

• Ben-Nun and Martínez (1998)

M. Ben-Nun and T. J. Martínez, J. Chem. Phys., 1998, 108, 7244.

• Saller and Habershon (2017)

M. A. C. Saller and S. Habershon, J. Chem. Theory Comput., 2017, 13, 3085–3096.

• Tully (1990)

J. C. Tully, J. Chem. Phys., 1990, 93, 1061–1071.

• Sun and Miller (1997)

X. Sun and W. H. Miller, J. Chem. Phys., 1997, 106, 6346–6353.

• Sun et al. (1998)

X. Sun, H. Wang and W. H. Miller, J. Chem. Phys., 1998, 109, 7064–7074.

• Stock and Thoss (2005)

G. Stock and M. Thoss, Adv. Chem. Phys., 2005, 131, 243–376.

• Bonella and Coker (2005)

S. Bonella and D. F. Coker, J. Chem. Phys., 2005, 122, 194102.

• Kim et al. (2008)

H. Kim, A. Nassimi and R. Kapral, J. Chem. Phys., 2008, 129, 084102.

• Ananth and Miller III (2010)

N. Ananth and T. F. Miller III, J. Chem. Phys., 2010, 133, 234103.

• Huo and Coker (2011)

P. Huo and D. F. Coker, J. Chem. Phys., 2011, 135, 201101.

• Huo and Coker (2012)

P. Huo and D. F. Coker, Mol. Phys., 2012, 110, 1035–1052.

• Kelly et al. (2012)

A. Kelly, R. van Zon, J. Schofield and R. Kapral, J. Chem. Phys., 2012, 136, 084101.

• Hsieh and Kapral (2012)

C.-Y. Hsieh and R. Kapral, J. Chem. Phys., 2012, 137, 22A507.

• Hsieh and Kapral (2013)

C.-Y. Hsieh and R. Kapral, J. Chem. Phys., 2013, 138, 134110.

• Ananth (2013)

N. Ananth, J. Chem. Phys., 2013, 139, 124102.

• Richardson and Thoss (2013)

J. O. Richardson and M. Thoss, J. Chem. Phys., 2013, 139, 031102.

• Richardson et al. (2017)

J. O. Richardson, P. Meyer, M.-O. Pleinert and M. Thoss, Chem. Phys., 2017, 482, 124–134.

• Agostini et al. (2014)

F. Agostini, A. Abedi and E. K. U. Gross, J. Chem. Phys., 2014, 141, 214101.

• Makri (2015)

N. Makri, Int. J. Quantum Chem., 2015, 115, 1209–1214.

• Kapral (2015)

R. Kapral, J. Phys.-Condens. Mat., 2015, 27, 073201.

• Duke and Ananth (2015)

J. R. Duke and N. Ananth, J. Phys. Chem. Lett., 2015, 6, 4219–4223.

• Hele and Ananth (2016)

T. J. H. Hele and N. Ananth, Faraday Discuss., 2016, 195, 269–289.

• Chowdhury and Huo (2017)

S. N. Chowdhury and P. Huo, J. Chem. Phys., 2017, 147, 214109.

• Meyer and Miller (1979)

H.-D. Meyer and W. H. Miller, J. Chem. Phys., 1979, 70, 3214–3223.

• Stock and Thoss (1997)

G. Stock and M. Thoss, Phys. Rev. Lett., 1997, 78, 578–581.

• Kelly et al. (2015)

A. Kelly, N. Brackbill and T. E. Markland, J. Chem. Phys., 2015, 142, 094110.

• Kelly et al. (2016)

A. Kelly, A. Montoya-Castillo, L. Wang and T. E. Markland, J. Chem. Phys., 2016, 144, 184105.

• Cotton and Miller (2013)

S. J. Cotton and W. H. Miller, J. Phys. Chem. A, 2013, 117, 7190–7194.

• Cotton and Miller (2013)

S. J. Cotton and W. H. Miller, J. Chem. Phys., 2013, 139, 234112.

• Miller and Cotton (2016)

W. H. Miller and S. J. Cotton, Faraday Discuss., 2016, 195, 9–30.

• Cotton and Miller (2016)

S. J. Cotton and W. H. Miller, J. Chem. Theory Comput., 2016, 12, 983–991.

• Cotton and Miller (2019)

S. J. Cotton and W. H. Miller, J. Chem. Phys., 2019, 150, 104101.

• Saller et al. (2019)

M. A. C. Saller, A. Kelly and J. O. Richardson, J. Chem. Phys., 2019, 150, 071101.

• Fenna and Matthews (1975)

R. E. Fenna and B. W. Matthews, Nature, 1975, 258, 573–577.

• Adolphs and Renger (2006)

J. Adolphs and T. Renger, Biophys. J., 2006, 91, 2778–2797.

• Ishizaki and Fleming (2009)

A. Ishizaki and G. R. Fleming, Proc. Nat. Acad. Sci., 2009, 106, 17255–17260.

• Sarovar et al. (2010)

M. Sarovar, A. Ishizaki, G. R. Fleming and K. B. Whaley, Nat. Phys., 2010, 6, 462 EP.

• Müller and Stock (1999)

U. Müller and G. Stock, J. Chem. Phys., 1999, 111, 77.

• Müller and Stock (1998)

U. Müller and G. Stock, J. Chem. Phys., 1998, 108, 7516–7526.

• Bonella and Coker (2003)

S. Bonella and D. F. Coker, J. Chem. Phys., 2003, 118, 4370–4385.

• Teh and Cheng (2017)

H.-H. Teh and Y.-C. Cheng, J. Chem. Phys., 2017, 146, 144105.

• Li et al. (1997)

Y.-F. Li, W. Zhou, R. E. Blankenship and J. P. Allen, *J. Mol. Bio. *, 1997, 271, 456 – 471.

• Tao and Miller (2010)

G. Tao and W. H. Miller, J. Phys. Chem. Lett., 2010, 1, 891–894.

• Kelly and Rhee (2011)

A. Kelly and Y. M. Rhee, J. Phys. Chem. Lett., 2011, 2, 808–812.

• Kim et al. (2012)

H. W. Kim, A. Kelly, J. W. Park and Y. M. Rhee, J. Am. Chem. Soc., 2012, 134, 11640–11651.

• Scholes et al. (2007)

G. D. Scholes, C. Curutchet, B. Mennucci, R. Cammi and J. Tomasi, J. Phys. Chem. B, 2007, 111, 6978–6982.

• Adolphs et al. (2007)

J. Adolphs, F. Müh, M. E.-A. Madjet and T. Renger, Photosyn. Res., 2007, 95, 197.

• Plenio and Huelga (2008)

M. B. Plenio and S. F. Huelga, New J. Phys., 2008, 10, 113019.

• Caruso et al. (2010)

F. Caruso, A. W. Chin, A. Datta, S. F. Huelga and M. B. Plenio, Phys. Rev. A, 2010, 81, 062346.

• Fassioli and Olaya-Castro (2010)

F. Fassioli and A. Olaya-Castro, New J. Phys., 2010, 12, 085006.

• Kolli et al. (2011)

A. Kolli, A. Nazir and A. Olaya-Castro, J. Chem. Phys., 2011, 135, 154112.

• Olaya-Castro and Scholes (2011)

A. Olaya-Castro and G. D. Scholes, Int. Rev. Phys. Chem., 2011, 30, 49–77.

• Chin et al. (2013)

A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga and M. B. Plenio, Nat. Phys., 2013, 9, 113 EP.

• Ishizaki et al. (2010)

A. Ishizaki, T. R. Calhoun, G. S. Schlau-Cohen and G. R. Fleming, Phys. Chem. Chem. Phys., 2010, 12, 7319–7337.

• Zhu et al. (2011)

J. Zhu, S. Kais, P. Rebentrost and A. Aspuru-Guzik, J. Phys. Chem. B, 2011, 115, 1531–1537.

• Wang et al. (2001)

H. Wang, M. Thoss and W. H. Miller, J. Chem. Phys., 2001, 115, 2979.

• Craig et al. (2007)

I. R. Craig, M. Thoss and H. Wang, J. Chem. Phys., 2007, 127, 144503.