DECIDE: A Deterministic Mixed Quantum-Classical Dynamics Approach

: Mixed quantum-classical dynamics provides an efﬁcient way of simulating the dynamics of quantum subsystems coupled to many-body environments. Many processes, including proton-transfer reactions, electron-transfer reactions, and vibrational energy transport, for example, take place in such open systems. The most accurate algorithms for performing mixed quantum-classical simulations require very large ensembles of trajectories to obtain converged expectation values, which is computationally prohibitive for quantum subsystems containing even a few degrees of freedom. The recently developed “Deterministic evolution of coordinates with initial decoupled equations” (DECIDE) method has demonstrated high accuracy and low computational cost for a host of model systems; however, these applications relied on representing the equations of motion in subsystem and adiabatic energy bases. While these representations are convenient for certain systems, the position representation is convenient for many other systems, including systems undergoing proton-and electron-transfer reactions. Thus, in this review, we provide a step-by-step derivation of the DECIDE approach and demonstrate how to cast the DECIDE equations in a quantum harmonic oscillator position basis for a simple one-dimensional (1D) hydrogen bond model. After integrating the DECIDE equations of motion on this basis, we show that the total energy of the system is conserved for this model and calculate various quantities of interest. Limitations of casting the equations in an incomplete basis are also discussed.


Introduction
Most systems of interest in chemistry and biology contain very large numbers of degrees of freedom (DOF).An exact simulation of the dynamics of such a system requires a fully quantum mechanical treatment of the entire system, which is computationally prohibitive due to the exponential scaling with the number of DOF.One approach to overcoming this issue is to treat a few DOF of interest quantum mechanically (i.e., subsystem) and the remainder classically (i.e., environment or bath) [1][2][3][4].Within this approach, the quantum subsystem is described in terms of a Hilbert space and the classical environment in terms of a phase space of positions and momenta.Situations in which the quantum subsystem is described by a non-Hermitian Hamiltonian have also been considered [5].Previously, these mixed quantum-classical techniques have been applied to a wide range of phenomena, e.g., proton transfer reactions [6,7], electron transfer reactions [8,9], proton-coupled electron transfer reactions [10][11][12][13][14][15][16], exciton transport in photosynthetic complexes [17][18][19], heat transport in molecular junctions [20][21][22][23][24], metamolecules [25], and magnetic molecules [26,27].
The quantum-classical Liouville equation (QCLE) has served as a starting point for simulating the dynamics of mixed quantum-classical systems [2,28,29].Over the years, a host of techniques have been developed based on approximate solutions of the QCLE [3,30].The most accurate of these techniques are the stochastic surface-hopping solutions [6,[31][32][33][34][35][36], but they require extremely large ensembles of trajectories to obtain converged results and, as a result, will be computationally prohibitive for many systems.On the other hand, the Poisson Bracket Mapping Equation (PBME) approach provides a highly computationally efficient approximate solution of the QCLE, but its applicability is mainly restricted to systems with weak subsystem-bath couplings [8,[37][38][39].Recently, the "Deterministic evolution of coordinates with initial decoupled equations" (DECIDE) method was developed, which offers a favourable balance between computational economy and accuracy [40].DE-CIDE has two main advantages compared to the stochastic QCLE-based methods: (i) It is a deterministic method that requires the integration of L 2 (L 2 − 2 + 2N) coupled differential equations (where L is the number of basis functions used to represent the subsystem and N is the number of environmental DOF).Typically, only a few thousand trajectories are required to obtain converged results, in contrast to several orders of magnitude more trajectories in the case of the stochastic QCLE-based methods; (ii) There is no need to diagonalize the Hamiltonian matrix on-the-fly.DECIDE has demonstrated great promise in solving the QCLE with high accuracy and low computational cost for a number of model systems, including the spin-boson model [40,41], Fenna-Matthews-Olson model [40,41], a three-state photo-induced electron transfer model [40], nonequilibrium spin-boson model [21,22,42], and a quantum battery model [43].
To date, the DECIDE equations have been represented in complete spin, subsystem, and adiabatic energy bases and successfully applied to simulate the dynamics of Hamiltonians expressed in these bases.However, the feasibility of solving the DECIDE equations in an incomplete position basis has yet to be explored.In this review, we explore this by first considering a one-dimensional, quadratic model of a hydrogen bond coupled to a solvent coordinate.
The review is organized as follows.In Section 2, we first discuss the Wigner-Weyl transform, which has been widely used in mixed quantum-classical dynamics, as it provides an invertible mapping between Hilbert space operators and phase space functions.In Section 3, we derive the DECIDE equations of motion and briefly discuss their underlying approximations and energy conservation.In Section 4, we provide a practical demonstration of how to implement the DECIDE approach using an incomplete position basis by considering a simple one-dimensional hydrogen bond model.We also briefly discuss issues that may arise when representing the DECIDE equations of motion for an arbitrary system in an incomplete basis.Finally, we summarize the main points of this review and discuss a future direction in Section 5.

Wigner-Weyl Transform
Before deriving the DECIDE equations of motion, we present an overview of the Wigner-Weyl transform and its properties.The Wigner-Weyl transform, which was first introduced by Wigner in 1932 [44], is a starting point for deriving mixed quantum-classical dynamics methods.Wigner's original goal was to find correction terms that would bridge quantum and classical mechanics [44,45].
The transform can be understood in terms of probability functions.For a wave function expressed in position space, ψ(x), the probability density is |ψ(x)| 2 .The wave function can also be expressed in momentum space as and its associated probability density is |φ(p)| 2 (h = h/2π and h is Planck's constant).Each probability density above depends on either x or p.The Wigner-Weyl transform provides a way to represent the probability distribution in terms of both position (x) and momentum (p).By definition, the Wigner-Weyl transform of an operator Â is where x and p are vectors containing the positions and momenta of all DOF in the system and y is an arbitrary integration variable.Using the above definition, one can obtain several properties of this transform.First, the identity operator 1 is given by Second, for an operator Â that is only a function of the position operator x, the Wigner-Weyl transform is Similarly, for an operator Â that is only a function of p, its Wigner-Weyl transform is A(p).
Third, the trace of the product of two operators Â and B is The Wigner function W(x, p), which corresponds to the probability density of x and p, is defined in terms of the Wigner transform of the density operator ρ = |ψ ψ| as For example, the Wigner functions for the ground and first-excited states of the quantum harmonic oscillator are shown in Figure 1.
The expressions for the Wigner functions are taken from Ref. [45].
In the Wigner-Weyl representation, the expectation value of an operator Â is given by where W(x, p) is normalized, i.e., The Wigner-Weyl transform of the product of two operators is [46] ( In the above equation, * denotes the Moyal star product [47,48] and Λ is the Poisson bracket operator, where an arrow denotes the direction of the operation.(The derivation of Equation ( 9) may be found in Ref. [46].)The exponential in Equation ( 9) may be expanded in the following series,

DECIDE Mixed Quantum-Classical Dynamics
In this section, we derive the DECIDE equations of motion [40].Let us start by considering a fully quantum system with a time-independent Hamiltonian Ĥ = ĤS ( x) + ĤB ( X) + ĤC ( x, X), where ĤS is the subsystem Hamiltonian, ĤB is the bath Hamiltonian, ĤC is the subsystem-bath coupling Hamiltonian, x denotes a set of generalized coordinates that provides a complete description of the state of the subsystem, and X = ( R, P) denotes the set of position and momentum operators of the bath.In the Heisenberg picture, the dynamics of the subsystem and bath coordinates are given by the quantum Liouville equation, namely where e i Lt Â = e i Ĥt/h Âe −i Ĥt/h .In the above equations, the time arguments outside the brackets indicate that one should first evaluate the commutator, then apply the time dependence.Taking a partial Wigner transform over the bath DOF (defined in terms of the bath DOF at t = 0) of the above equations and retaining only zero-order terms in h in the resulting Moyal product expansion yields, Assuming that the subsystem and bath coordinates are initially decoupled, evaluating the partial Wigner transform of the commutators in the above equations, and retaining only first-order terms in h in the resulting Moyal product expansion leads to the DECIDE equations of motion, where the antisymmetrized Poisson bracket is given by { ĤW , ÂW } a = 1 2 { ĤW , ÂW } − The assumption that the subsystem and bath are initially decoupled leads to ([ ) It should be noted that the partially Wigner-transformed Hamiltonian retains an operator character because the transform was not performed on the subsystem DOF.In addition, to account for noncommutativity in the above equations, one has to replace any coupling terms in H C with their Weyl-ordered symmetric forms, e.g., xX would be replaced by 1 2 ( xX + X x).The errors caused by the approximations above become significant when the subsystem dynamics is highly non-Markovian, namely when the dependence on the initial bath coordinates is important [40].This can be the case for very low bath temperatures, very strong subsystembath coupling, and very slow baths, i.e., when memory effects are very pronounced [22,40].
The general form of the partially Wigner-transformed Hamiltonian considered in the next section is given by where r and p are vectors containing the subsystem position and momentum operators, respectively, and m and M are vectors containing the masses of the subsystem and bath DOF, respectively.For this form of Hamiltonian, the DECIDE equations of motion for the subsystem and bath DOF are Finally, one must select a convenient basis for the quantum subsystem to solve these equations.
To show that the total energy of the system is conserved over the course of the dynamics, we start by taking the time derivative of a matrix element (in any basis) of the Hamiltonian, i.e., where H αα = α| ĤW |α .Substituting the expressions for the time derivatives of the subsystem and bath DOFs in Equation ( 16) into Equation ( 17), we find that thereby proving that the energy of the total system is conserved.

Hydrogen Bond Model
We now demonstrate how to implement the DECIDE approach using an incomplete position basis by considering a simple, quadratic, 1D hydrogen bond model.

Model
The Weyl-ordered, partially Wigner-transformed Hamiltonian of the hydrogen bond model is given by where x/ p is the position/momentum operator of the proton, X/P is the position/momentum of the solvent coordinate, and A 0 and k 1 are constants.Plots of the protonic potential for two values of X are shown in Figure 2. A convenient basis for solving this proton transfer problem is a set of quantum harmonic oscillator wave functions, i.e., where H n is the nth Hermite polynomial, b is an arbitrary constant, and n = 0, 1, 2, . ... Note that if X = 0 in Equation ( 19), then the exact solution of the time-independent Schrödinger equation for this model is φ n (x) with b = 2mA 0 /h 2 1/4 .In our simulations, we take b = 8mA 0 /h 2 1/4 .The eigenfunctions of ĤW may be expanded in this basis as where c i n is the coefficient of the nth basis function.For this expansion, we use 12 basis functions (n = 0, . . ., 11).Substituting the eigenfunction expansion into the time-independent Schrödinger equation, we obtain an eigenvalue problem of the form Hc = cE, where H is the Hamiltonian matrix expressed in the chosen basis, c is a matrix containing the c i n 's of the ith eigenfunction, and E is a diagonal matrix containing the eigenvalues E i (X).The Hamiltonian matrix elements are calculated by numerical integration, and the eigenvalue problem is solved (to obtain the c i n 's and eigenvalues) using built-in functions in MATLAB.Figure 3 shows the ground-and first-excited state adiabatic potential energy surfaces [E 0 (X) and E 1 (X), respectively] calculated using this method.
It turns out that the time-independent Schrödinger equation for the Hamiltonian in Equation ( 19) can be solved analytically (and exactly) by first rewriting the potential energy The resulting adiabatic energies are where n labels the energy level.Comparing our variational results to the exact ones, we find that good agreement is obtained when X is close to 0, but differences emerge when |X| > 0.25 (see Figure 3).The convergence of the variational result can be improved by increasing the number of basis functions; however, the level of accuracy achieved with 12 basis functions is sufficient for our purposes.The following parameter values were used to obtain these surfaces: , and m = 1 amu.

Dynamics Simulation Details
Using Equation ( 16) and the Hamiltonian in Equation ( 19), one may obtain the DECIDE equations of motion for the proton transfer model in the quantum harmonic oscillator basis, viz., where, for example, x αα = α| x|α and |α is the α th harmonic oscillator state.Given that the basis set is composed of 12 functions, the total number of coupled differential equations is 12 2 × 4. We assume the initial state of the composite system to be factorized, i.e., ρ(0) = ρS (0)ρ B,W (0), where ρS (0) is the initial density operator of the subsystem (viz., the proton) and ρ B,W (0) is the initial Wigner distribution of the bath (viz., the solvent coordinate).In our simulations, the system is initialized in the adiabatic ground state |ψ 0 of the Hamiltonian in Equation (19), such that ρS (t) = |ψ 0 ψ 0 |.In the harmonic oscillator basis, the initial values of the matrix elements of the subsystem coordinates are given by the following analytical expressions: The initial values of the matrix elements of the bath coordinates are given by X αα (0) = X(0)δ α,α , P αα (0) = P(0)δ α,α , where X(0)/P(0) is the initial value of the bath position/momentum.In practice, X(0) and P(0) may be sampled from Wigner distributions corresponding to thermal equilibrium states.However, in this study, we choose X(0) = −0.2Å and P(0) = 50 a.u., since we only generate a single trajectory for demonstration purposes.Starting from these initial conditions, the DECIDE equations of motion in Equation ( 22) are integrated using the Runge Kutta fourth-order method [49] with a time step of 1 a.u. to simulate the dynamics of the system.In general, a particular integrator is chosen to achieve an optimal balance between efficiency and accuracy for a given system.The expectation value of a property Â(t) may be calculated according to However, since we only consider single trajectories for this demonstration, we take ρ B,W (0) = 1.In particular, we are interested in the conservation of the total energy along a trajectory, which is calculated according to

Results
Figure 4 shows the values of X, adiabatic energy, bath kinetic energy, and total energy along a representative trajectory.For a positive initial momentum P(0), we see that X increases up to about 4000 time steps and then starts to decrease.The adiabatic energy increases while the bath kinetic energy decreases over the course of this trajectory.Finally, as seen in Figure 4d, the total energy of the system is well-conserved along this trajectory.

Limitations of Position Representation
Although the DECIDE method is capable of accurately solving the 1D quadratic hydrogen bond model, problems in energy conservation may arise when integrating the DECIDE equations of motion for more complex models.A detailed analysis has revealed (results not shown) that a breakdown in energy conservation could occur due to a combination of basis set incompleteness (inherent to the position representation) and the presence of higher-order (i.e., cubic, quartic, etc.) terms in the Hamiltonian.In the DECIDE equations of motion, matrix elements of second or higher powers of the coordinates [e.g., (x 2 ) αα ] are not evolved directly.Since only matrix elements of the first powers of the coordinates are evolved directly, the matrix elements of higher-order terms are evaluated by first inserting complete sets of states between each linear term in the product (e.g., (x 2 ) αα = ∑ β x αβ x βα ).Therefore, if the position basis is incomplete, such terms can be evaluated inaccurately, which leads to an accumulation of numerical errors and, in turn, a breakdown in energy conservation.In this connection, even in the 1D quadratic hydrogen bond model, there will be small numerical errors when calculating energies due to the presence of x 2 and p 2 terms in the Hamiltonian [see Equation (19)].

Conclusions
Mixed quantum-classical dynamics is a useful approach for simulating the dynamics of quantum processes occurring in chemical and biological systems with large numbers of DOF.Among the various mixed quantum-classical dynamics methods, the DECIDE method has exhibited a highly favourable balance between computational cost and accuracy for treating model Hamiltonians expressed in complete energy bases.
In this review, we explained the DECIDE method and demonstrated it on a 1D quadratic model of a proton in a hydrogen bond.This required casting the DECIDE equations of motion in an incomplete quantum harmonic oscillator position basis.We showed that with a sufficiently large basis, it is possible to generate trajectories that conserve the total energy of the system well.However, for higher-order models expressed in an incomplete position basis, energy conservation may break down and thus, it would be necessary to re-express the Hamiltonian in an energy basis.This approach will be explored in a future study.

Figure 2 .
Figure 2. Plots of the protonic potential along the proton coordinate x.The blue and red curves represent the potential for two different values of the solvent coordinate X, which modulates the position and depth of the protonic potential.The dashed line corresponds to the ground state energy for a particular value the solvent coordinate.

Figure 3 .
Figure 3. Ground-and first-excited state adiabatic potential energy surfaces of the 1D quadratic proton transfer model.V denotes a variational result, and A denotes an analytical result.The following parameter values were used to obtain these surfaces: A 0 = 3367.6cm −1 Å −2 , k 1 = 1.7 × 10 4 cm −1 Å −2 , and m = 1 amu.