Nonadiabatic molecular dynamics simulation: An approach based on quantum measurement picture

Mixed-quantum-classical molecular dynamics simulation implies an effective measurement on the electronic states owing to continuously tracking the atomic forces.Based on this insight, we propose a quantum trajectory mean-field approach for nonadiabatic molecular dynamics simulations. The new protocol provides a natural interface between the separate quantum and classical treatments, without invoking artificial surface hopping algorithm. Moreover, it also bridges two widely adopted nonadiabatic dynamics methods, the Ehrenfest mean-field theory and the trajectory surface-hopping method. Excellent agreement with the exact results is illustrated with representative model systems, including the challenging ones for traditional methods.


I. INTRODUCTION
A full quantum mechanical evaluation for molecular dynamics (MD) would quickly become intractable with the increase of atomic degrees of freedom. As alternatives in practice, some mixed-quantum-classical (MQC) MD approaches were developed and proved to be very powerful [1,2]. A typical class of such studies is nonadiabatic MD. Nonadiabatic effects are of crucial importance in the proximity of conical intersection, see Fig. 1(A), where the energy separation between different potential energy surfaces (PESs) becomes comparable with the nonadiabatic coupling. MQC treatment of nonadiabatic MD has a long history. The widely applied schemes include the so-called "Ehrenfest" or "time-dependent-Hartree mean-field" approach [3][4][5], the "trajectory surface-hopping" methods [6][7][8][9][10][11], and their mixed scheme [12][13][14][15]. The former views that the electronic wave function is in general a linear combination of Born-Oppenheimer adiabatic states, and the atomic effective potential (and force) is calculated by averaging the electronic Hamiltonian over such wave function. The trajectory surface-hopping scheme is in a different extreme. It believes that the trajectories should split into branches, i.e., each trajectory should be on one state or another, not somewhere in between. In this type of theories the trajectories distribution is achieved by allowing hops between PESs according to some probability distribution.
Among the trajectory-based surface hopping methods, the most popular one is Tully's fewest switches surface hopping (FSSH) approach [11], together with its variations [2]. In this approach, nonadiabatic dynamics is treated by allowing hops from one PES to another, with * Electronic address: lixinqi@bnu.edu.cn

Electronic States
Atomic Forces the hopping probability determined by the weight change of the respective electronic states that are in a quantum superposition. From observation that the classical atomic motion must decohere the electronic subsystem (from quantum superposition), considerable efforts were pushed towards accounting for the associated decoherence effect [16][17][18][19][20][21][22]. One may notice that the FSSH treatment has an obvious flaw from basic physical point of view. By an analogy with quantum measurement or the popular Schrödinger's Cat paradox, as illustrated in Fig. 1(B), the FSSH scheme simply indicates that, while having found the Cat definitely "alive" or "dead", one is still treating the radioactive decay in "quantum superposition". In addition, it was also noticed that the FSSH scheme involves a man-made hopping algorithm to generate the stochastic surface-switching events [2,18,19]. In this work, by explicitly identifying the role of the atomic motion as a quantum measurement to the electronic subsystem, we propose a novel quantum trajectory mean field (QTMF) approach to nonadiabatic MD simulations. The protocol is developed from an insight that the involved quantum weak measurement actually serves as an interface between the quantum and classical parts of the MQC strategy. Our scheme also naturally unifies the Ehrenfest-type mean field theory and the trajectory surface-hopping method. As illustrated by excellent agreement with the exact results on the three representative models discussed by Tully [11], the QTMF approach can eliminate all the unsatisfactory features of the FSSH method.

Let us start with the electronic Hamiltonian
The potential v(r; R), which includes the nuclear potential, depends on electronic coordinates operator r ≡ {r j } and also atomic configuration, R ≡ {R α (t)}, that is assumed a set of classical dynamics variables. Expand the electronic wavefunction with an orthogonal basis set functions, Ψ(r, t; R) = j c j (t)φ j (r; R). In particular one often exploits the Born-Oppenheimer (BO) adiabatic wavefunctions. These are the instantaneous eigenstates of H el (r; R), satisfying H el (r; R)φ j (r; R) = ε j (R)φ j (r; R), the standard output of quantum chemistry computation. Each BO energy serves as the BO potential energy surface (PES), ε j (R) ≡ V j (R), for nuclear motion.
Without loss of generality, we proceed with the BO representation hereafter. The Schrödinger equation for the coherent electronic evolution readṡ with the nonadiabatic coupling characterized by Treating atomic motion with classical trajectories on individual PESs, {V j (R)}, the highly celebrated FSSH method [11] goes by a Monte-Carlo algorithm as follows. It uses Eq. (2) for the hopping probability from a given V j (R) to another. That is p j (t) = |c j (t)| 2 − |c j (t + ∆t)| 2 /|c j (t)| 2 , the normalized population change in the BO electronic state |φ j (R) . However, this algorithm is of problematic basis. It completely neglects the influence of classical trajectories back onto the electronic state evolution. Atomic motion that experiences a series of single PESs should collapse the electronic state from a quantum superposition, given by Eq. (2), onto the corresponding single BO basis state, due to the entanglement-type correlation between the two subsystems. In other words, atomic motion is continuously extracting information on the electronic state, via the correlation between the PES and BO basis state. For instance, the force experienced by atomic motion plays essentially the same role as the meter's output in quantum measurement process. Based on this insight, we propose to apply the well-established quantum trajectory equation, in replacement of Eq. (2), to account for the backaction effect of the atomic "meter" on the electronic subsystem [23]: In this equation, ρ denotes the reduced density matrix of electronic state, with diagonal elements for BOstate population probabilities, and off-diagonal ones for their coherence. The first term in Eq. (4) describes the same coherent dynamics of Eq. (2), corresponding to the Ehrenfest mean-field approach. The second term accounts for the decoherence effect owing to ensemble average of the nuclear degrees of freedom, with an overall rate Γ. The third term, significantly, reflects the backaction effect of the atomic motion in each single trajectory realization, with a rate γ F for force-mediated information gain and γ ′ for information gain by other means, e.g., the nuclear coordinate (atomic configuration) and velocity (atomic kinetic energy). Here we omitted the PES indices of the rates for brevity and for a general description. In the following Sec. II (A) and (B) we will explain how these rates can be implemented in MQC-MD simulation. Before that, we describe in more detail the decoherence and measurement backaction terms.
| indicates a dephasing between states φ j (R) and φ k (R). The last backaction term, explicitly, is described in terms of an superoperator as . Involved in Eq. (4) for this back-action effect are also the quantum jump (from the Copenhagen postulate) related stochastic noises, {ξ jk (t)}, which satisfy the ensemble average property of . From the quantum trajectory theory [23], the last term in Eq. (4) has a role of collapsing the electronic state from a quantum superposition onto a single BO basis state. Therefore, now, the issue on "devising" hopping algorithms that are not contained in Eq. (2) does no longer exist anymore.

A. Information Gain Rates
In the MQC-MD approach, the nuclear part is treated classically. As a consequence, just like Tully pointed out in his pioneering work [11], the classical force experienced by the atomic motion in the no-transition adiabatic area should come from a single PES. This indicates that, from a measurement perspective, the classical force plays a role of measurement output. Below we analyze this forcemediated information gain rate (γ F ).
We know that the emergence of classicality from a closed quantum system is a fundamental puzzle in quantum mechanics. In essence, this transition is accompanied by quantum jumps. This implies that the classical force has certain stochastic fluctuations. Since the atomic motion is much slower than its electronic counterpart, it would be reasonable to use a coarse-graining force to evolve the Newton equation. Let us denote the coarse-grained fluctuating component byF j (t), which is an average over a characteristic time "τ c " around t as follows:F F j is the BO force associated with the j th PES at R(t).
Notice also that ξ j dt = dW j , the Wiener increment, has a magnitude order and dimension of √ dt [23]. As a result, the coarse-grainedF j (t) is no longer δ-function correlated, but has a correlation function of where E[· · · ] denotes an ensemble average over the stochastic realization ξ j , which satisfies E[ξ j (t)ξ j (t ′ )] = δ(t−t ′ ). Accordingly, the zero-frequency spectrum of the force-force correlation function reads Now we return to the original (stochastic) force of the j th PES, F j (t) =F j (t) + S j ξ j (t). The total average force F j , given by averaging F j over time interval (t, t + τ m ), is a stochastic variable satisfying a Gaussian distribution P with the variance given by D j = S j /τ m . Following the theory for realistic quantum measurements [24], the state distinguishable criterion allows us to extract the measurement time (t m ) which is minimally required to identify the state being |φ j or |φ k . Obviously, t m is given by τ m from Eq. (8) in the case of equality. Then, the information gain rate γ F,jk in Eq. (4) coincides with 1/t m , taking a compact form as, As inferred from the coarse-grained force, the characteristic time τ c physically scales atomic motion that is typically of picosecond. With respect to the femtosecond timescale of the electronic part, in practice we adopt 1/τ c ≃ 10 −3 in (reduced) units of the electronic energies.
Favorably, the information gain rate given by Eq. (9) is of configuration (R) dependence, but it does not need any microscopic information of the nuclear (quantum) wavepackets. It allows thus for a convenient implementation even for simulation on complex molecular systems. Except for the force-mediated information gain discussed above, there exit also other channels of information gain, which are formally accounted for by γ ′ in Eq. (4). The channels include the nuclear coordinate R and velocityṘ in the MQC-MD simulation. For instance, if we performed a microscopic full quantum treatment, the nuclear wavepacket would have distinct spatial extension along different PES. With this "knowledge" in priori, one can infer certain information for the electronic state from the classical "output" R. Another information channel is the nuclear kinetic energy (associated withṘ). For an initial state with specific energy, the distinct kinetic energies on different PESs in MQC-MD simulation exposure also some information of the electronic state. Unfortunately, to our knowledge, the information gain rates through these channels have not yet well developed so far. However, fortunately, as we will elaborate further in next subsection, the specific form of these rates are not important for us to get the correct results.
B. More Remarks on Eqs. (4) and (9) Generally speaking, a desirable MQC-MD approach should satisfy two requirements. One is that the equation for the electronic subsystem should satisfy the ensemble average property, corresponding to averaging the nuclear degrees of freedom from the exact dynamics of the full electronic-plus-nuclear system. Another is that the equation should allow for performing reasonable classical MD simulation (with correct "force" and "kinetic energy") on the nuclear subsystem. Our protocol is to combine Eq. (4) with a classical MD simulation to fulfill these two requirements. That is, the second term of Eq. (4) satisfies the first requirement, and the last term satisfies the second one.
In Eq. (4), we distinguished the information-gain rates γ F and γ ′ from the overall decoherence rate Γ. Formally, we may express Γ = γ F + γ ′ +γ. As discussed above in Sec. II (A), γ F and γ ′ are, respectively, the information gain rates mediated by force and other channels (e.g., R andṘ). Therefore,γ represents the decoherence rate not unraveled (the rate of information loss). For practical applications, we propose the following strategies to implement these rates: 1. In the nonadiabatic crossing region, use the rate γ F given by Eq. (9) to approximate the total information gain and decoherence rates. This approximation is from an insight that, in the conical intersection area, the total decoherence rate should be quite weak, otherwise the result will be strongly distorted from the correct one (the Ehrenfest meanfield approach is a support to this viewpoint). We believe that our coarse-graining argument for obtaining Eq. (9) gives a reasonable order of magnitude for this rate.
2. Apart from the crossing area, add a nonzero γ ′ into Eq. (4), by using a simple phenomenological parameter with similar/stronger magnitude order of 1/τ c , or by certain more sophisticated manner [19,22,25]. We may remark that this different implementation of γ ′ is anticipated to result in slight difference only in the "narrow" region between the nonadiabatic crossing and the no-transition adiabatic areas. It will affect neither the molecular dynamics in the broad adiabatic region, nor the ensemble statistical properties. For the overall decoherence rate Γ, one can set either Γ = γ F + γ ′ , or Γ = γ F + γ ′ +γ by including a more informationloss rate. However,γ will have no effect, since γ F and γ ′ will collapse the system in the adiabatic area onto a single PES in each trajectory realization, implying a mixed state after ensemble average. The role ofγ is simply to facilitate the formation of an ensemble-averaged mixed state.
Finally, we mention a special case that may remind our attention. For parallel PESs, the "force output" reveals no information of the electronic state, thus giving a vanished measurement rate. This is in consistence with the result of Eq. (9). In this case, the rate γ ′ from other informational channels will collapse the system onto a single PES in the adiabatic area. Whether or not collapsing the system onto a single PES in this case will have no effect on the force, but it does affect the kinetic energy (nuclear velocity) that should be of importance in the MD simulation for real molecular systems.

C. Issue of Energy Conservation
In the MQC-MD approach, the atomic motion defines a time-dependent electronic Hamiltonian, which does not conserve the electronic energy. In turn, the electronic energy defines a potential to guide the classical motion of atoms. The sum of the kinetic and potential energies, , is conserved, simply as the situation in the Ehrenfest mean-field approach.
However, Eq. (4) is stochastic. This would lead to a stochastic potential energy V (R). The non-analytic V (R) makes the force not perfectly defined in mathematics, causing thus some errors in determining the nuclear velocity. This would violate slightly the total energy conservation. Noticeably, in the present QTMF approach, this violation is quite weak (particularly if a coarse graining procedure is involved), unlike the drastic violation in the FSSH scheme where an energy calibration must be performed after each hopping event. For practice of the QTMF approach, we propose the following scheme to address this issue: 1. Define the whole simulation region with the criterion V (R) ≤ E 0 , where E 0 is the initial energy of the whole system. Of course, this renders also that V j (R) ≤ E 0 once the system is fully collapsed onto the j th -PES.
2. If V (R 1 ) = E 0 occurs at R 1 , reset the system to its proximity point R 2 , given by V M (R 2 ) = E 0 .
Here, V M is the renormalized Ehrenfest mean-field potential energy and determined as follows: at R 2 , keep the electronic wavefunction unchanged as the one at R 1 ; subtract the lowest PES component and re-normalize the wavefunction; then use the renormalized wavefunction to calculate the Ehrenfest V M (R 2 ).
3. Restart the MD evolution from the determined proximity point R 2 , with the original superposition of BO PESs at R 1 but a newly assigned atomic kinetic energy of E 0 − V (R 2 ) and the momentum direction opposite to that at R 1 .
4. After passing through the nonadiabatic crossing area, check the total energy of the collapsed state (onto a single PES) and make it be E 0 via proper modification to the kinetic energy.

III. ILLUSTRATIVE EXAMPLES
In this section we present our QTMF results versus the exact and FSSH counterparts, on the well-known set of Tully test systems [11], each of them being a onedimensional two-surface model, with an atomic mass of M = 2000 a.u. (all parameters below are in atomic unit). The scheme for exact quantum dynamics simulation was clearly described in Ref. [11]. In the present work, we simply extract the results from Ref. [11] for comparison. In our simulation, we assume an incident energy E 0 = k 2 /(2M ) to initiate the system evolution. And, as discussed earlier, we adopt 1/τ c = 10 −3 . In the nonadiabatic coupling area, we approximate the entire decoherence and information rates by γ F through Eq. (9). In the adiabatic (no-transition) area, we add γ ′ = 10 −2 to account for the backaction effect of other informational channels, and set Γ = γ F + γ ′ . As explained in Sec. II (A) and (B), the choice of γ ′ in the adiabatic area can be rather arbitrary, having almost no influence on the results. For each model, we run 2000 stochastic trajectories and accordingly determine the population probabilities of the final "products". Also, each trajectory begins with the classical particle (atom) on the lower energy surface at x = −10, with an incident momentum to the right, and ends at |x| = 15. Model-I: Single-Avoided Crossing -The diabatic potential matrix elements for this model are Set A = 0.01, B = 1.6, C = 0.005, and D = 1. The adiabatic potential surfaces of this model are plotted in Fig. 2(a), while the results are shown in Fig. 2 Desirably, both our QTMF and the FSSH schemes work very well for this model, being almost in an overall agreement with the exact results. We only make two remarks on the extremely quantum regime. (i) The steep stepfunction behavior at k ∼ 5 is an indicator for the failure of almost all semiclassical MD methods, i.e., failing to predict tunneling through the barrier on the lower surface at very low momentum. Physically, in our case this is caused by setting the semiclassical rule of energyconservation when propagating the state. Giving up this restriction at k ∼ 5, we can actually recover the exact result. (ii) Another interesting quantum regime is k ∼ 8 (7.7 < k < 8.9), which is above the threshold for transmission. Satisfactorily, both QTMF and FSSH captured the essential physics here, e.g., predicting the small amount of particle reflections. This is somehow a challenging test for any semiclassical approaches. Model-II: Dual-Avoided Crossing -This is a more challenging model and contains two avoided crossings. The key feature of this model is the Stueckelberg oscillations, owing to quantum interference between the successive nonadiabatic quantum transitions. The diabatic potentials for this model are given by where A = 0.10, B = 0.28, C = 0.015, D = 0.06, and E = 0.05. The adiabatic potentials of this model and results comparison are shown in Fig. 3. At high incident energies, both the FSSH and QTMF can give correct results in good agreement with the exact ones, all exhibiting the expected Stueckelberg oscillations. However, at low energies, the FSSH method fails to predict both the transmission and reflection probabilities on the lower surface, see Fig. 3 (b) and (d) in the low energy regime. In particular, the FSSH algorithm overestimates the amount of reflection by an order of magnitude (roughly a factor of 10). This overestimation is owing to the artificial hopping algorithm, which results in too many upward transitions. In sharp contrast, our QTMF approach can physically rule out this drawback. Model-III: Extended Coupling -This is the most challenging model for classical mechanics based approach to address, since it involves an extended region of strong nonadiabatic coupling. The diabatic potentials are The parameters are A = 6 × 10 −4 , B = 0.1 and C = 0.9. The adiabatic potentials and comparative results are shown in Fig. 3. We see that, unfortunately, the FSSH algorithm completely fails for the reflection probabilities, to either the upper or lower surface. This failure clearly indicates that the FSSH algorithm breaks down in strong quantum transition region. Again, in sharp contrast, our QTMF approach gives satisfactory results even for this very demanding model.

IV. SUMMARY
To summarize, we have proposed a quantum trajectory mean field (QTMF) approach to the powerful mixedquantum-classical molecular dynamics simulation with surface hopping. Simulations on three nontrivial models are quantitatively satisfactory. While Eq. (9) offers a compact position-dependent measurement rate on atomic motion timescale (τ c ), the present study reveals also an important observation: results are rather insensitive to the choice of decoherence rate, as long as it is weak (∼ 1/τ c ) in the nonadiabatic crossing area. Unraveling any decoherence rate in the no-transition adiabatic area can stochastically collapse the system onto a single potential surface, and gives about the same satisfactory statistical results.
In this context we would like to remark that quantum superposition is rooted in the exact quantum dynamics simulation, but involving not at all the concept of classical atomic trajectory. In Tully's fewest switches algorithm, while the evolution of electronic state is not consistent with the successive complete surface hopping, it keeps the essential feature of quantum superposition. It is merely this reason, in our opinion, that makes the most celebrated FSSH approach be often comparable to the exact results from full quantum dynamics simulation.
Compared with the FSSH approach, the QTMF scheme is founded on a more physical and simpler treatment. For the electronic part, the replacement of the Schrödinger equation with a master equation will increase only negligible amount of computational complexity, since the involved BO states are very few (for instance, only two in most real molecular simulations). On the other hand, the QTMF scheme avoids the hopping algorithm and simplifies the procedures of calibrating the total energy. This will speed the simulation. As a conservative estimate, the time cost of the QTMF scheme is comparable to or less than the FSSH approach (and its many variations). With this computational efficiency together with the satisfactory accuracy, and most importantly its physical foundation, the proposed QTMF scheme is anticipated to be a powerful tool in real MD simulations.