We propose an approximate method for evaluating the importance of non-Born–Oppenheimer effects on the quantum dynamics of nuclei. The method uses a generalization of the dephasing representation (DR) of quantum fidelity to several diabatic potential energy surfaces and its computational cost is the cost of dynamics of a classical phase space distribution. It can be implemented easily into any molecular dynamics program and also can utilize on-the-fly ab initio electronic structure information. We test the methodology on three model problems introduced by Tully and on the photodissociation of NaI. The results show that for dynamics close to the diabatic limit, the decay of fidelity due to nondiabatic effects is described accurately by the DR. In this regime, unlike the mixed quantum-classical methods such as surface hopping or Ehrenfest dynamics, the DR can capture more subtle quantum effects than the population transfer between potential energy surfaces. Hence we propose using the DR to estimate the dynamical importance of diabatic, spin-orbit, or other couplings between potential energy surfaces. The acquired information can help reduce the complexity of a studied system without affecting the accuracy of the quantum simulation.

Introduction. The nonadiabatic effects play an important role in many chemical phenomena and often must be taken into account in accurate calculations of molecular properties such as spectra or reaction rates.1 Many nonadiabatic quantum (QM) simulations are performed in the diabatic (or quasidiabatic) basis, which offers several computational advantages, especially in the vicinity of conical intersections of adiabatic potential energy surfaces (PESs). Usually, more PESs must be included to describe a system of interest accurately. One may ask: Which diabatic surfaces are important? How much is the result affected by neglecting the less important surfaces? The method we propose below quantifies the importance of couplings between PESs and thus can help in answering these questions.

A direct way to determine whether the coupling of diabatic PESs affects the result of a simulation is to compute the desired quantity by running QM dynamics for both the uncoupled and coupled systems and to compare the results. Here, we propose a more general approach, which can give the information about the importance of couplings for all observables at the same time. It consists in comparing the wave functions by computing the QM fidelity, defined as2 

FQM(t)=|f(t)|2=|ψ0(t)|ψp(t)|2.
(1)

In the general setting, ψ0(t) and ψp(t) are wave functions evolved in the unperturbed and perturbed systems, respectively. In our case, the “unperturbed” represents the uncoupled system, “perturbed” the coupled system, and ψ denotes the full molecular wave function, i.e., it includes both nuclear and electronic degrees of freedom. Values of FQM(t) close to unity imply that the perturbation is not important: The uncoupled Hamiltonian can be used in quantitative simulations. Values significantly below unity imply that the perturbation is important: The affected PESs and couplings should be included in the simulation.

Instead of computing fidelity directly from the definition (1), one can use the dephasing representation (DR),3–5 which is an efficient semiclassical approximation of fidelity.6,7 Recently, the DR was used successfully to evaluate the accuracy of QM dynamics on a single but approximate PES.8,9 Below, we generalize this method to several surfaces and use it to evaluate how the dynamics is affected by couplings between the diabatic PESs.

Unlike the computational cost of a direct QM calculation, the cost of the DR does not grow exponentially with the number D of degrees of freedom.5,8 Hence, the DR can be used for many-dimensional systems inaccessible to current methods of QM dynamics. An advantage of the DR in comparison with mixed QM-classical methods for nonadiabatic dynamics, such as the Ehrenfest dynamics, various surface hopping methods,10,11 or methods in which the classical limit is obtained by the linearization of the QM propagator,12 is that the DR, being a semiclassical method, takes the nuclear coherence effects into account at least approximately. Other semiclassical approaches to nonadiabatic dynamics exist, including, among others, the multiple-spawning methods13 or methods14 based on the Miller–Meyer–Stock–Thoss classical electron model. The advantage of the DR is that unlike the majority of semiclassical approaches, the DR does not require the Hessian of the potential energy, which is the most expensive element of first-principles semiclassical dynamics methods (see, e.g., Ref. 15).

Clearly, the above-mentioned advantages do not come for free. First, unlike most other approaches noted above, the DR is not a general dynamical method. It can only describe properties expressible in terms of fidelity. Second, the DR in the diabatic basis is expected to work accurately mainly when the coupling is relatively weak or the nuclear velocities are large, i.e., when the dynamics is close to the diabatic limit. Nevertheless, in many cases, the DR works very well even for strong perturbations that lead to completely different classical phase space distributions for uncoupled and coupled Hamiltonians.

Applications for which the DR method is particularly well suited include the chemical reactions which proceed near the diabatic limit. The photodissociation of bromopropionyl chloride or bromoacetyl chloride16 serve as examples. Another possible application would be to quantify the relative importance of Hamiltonian coupling terms of various origins, e.g., of the spin-orbit versus diabatic couplings. For instance, in the Cl(P2)+H2 reaction, the spin-orbit coupling terms clearly dominate over diabatic couplings, which can therefore be neglected in a simulation.17 

Theory. The starting point is the Hamiltonian describing nuclear motion in a molecule, expressed in the diabatic basis. To compute the decay of fidelity due to the coupling between n PESs, the Hamiltonian n×n matrix is split into the uncoupled (i.e., diagonal) part Ĥdiag and the coupling (i.e., off-diagonal) part ΔV̂,

Ĥ=Ĥdiag+ΔV̂,
(2)
Ĥdiag=T̂+V̂diagandΔV̂=V̂offdiag.
(3)

In Eq. (3), T̂ is the diagonal nuclear kinetic energy matrix and V̂diag contains the diabatic PESs, which are uncoupled in Ĥdiag and coupled in Ĥ by the elements of V̂offdiag. (Bold face denotes n×n matrices; hat ̂ denotes operators.)

To derive the DR, one starts from the expression for QM fidelity amplitude, applicable to both pure and mixed states5 and generalized to the multi-PES setting

fQM(t)=Tr(eiĤdiagt/ρ̂e+iĤt/),
(4)

where ρ̂ is the density operator of the initial state. Generalizing the derivation from Ref. 18, fidelity amplitude can be written exactly as

fQM(t)=TrdxρW(x)(e+iĤt/eiĤdiagt/)W(x),
(5)

where ρW is the Wigner transform of the initial state ρ̂

[ρW]ij(x)=hDdξqξ2|ρ̂|q+ξ2exp(iξp),

and x denotes the point (q,p) in the 2×D-dimensional phase space. Approximating the Wigner transform of the product of the time-evolution operators, one arrives at the DR expression

fDR(t)=Tr[dx0ρW(x0)eiΔS(x0,t)/],
(6)
ΔS(x0,t)=0tdτΔVW[xτ(x0)],
(7)

where ΔS(x0,t) is the action at time t due to Wigner representation ΔVW of ΔV̂ along the trajectory xτ of Hdiag. Note that ΔVW, ΔS, and ρW are matrix quantities. If ΔV̂ contains only diabatic coupling elements, then ΔV̂ΔV(q̂) and ΔVW(x)=ΔV(q). For initial Gaussian wave packets, the Wigner function equals the classical phase space density which is strictly a probability distribution.

Since Trdx0ρW(x0)=1, it follows from Eq. (6) that fDR can be computed as a Monte Carlo average exp(iΔS(x0,t)/)ρW(x0), with initial conditions sampled from the Wigner distribution ρW(x0) of the initial state. However, as ΔS in Eq. (7) is much smaller than the action in other semiclassical methods, the DR alleviates the notorious “sign problem.”

If the initial state lies on a single surface n, i.e., if [ρW]nj=0 and [ρW]jj=0 for jn, then no other elements of V=Vdiag+Voffdiag than Vnn and Vnj,jn enter the calculation. This means that no information about other diagonal elements Vjj is used in the DR calculation. Thus, except for special cases, FDR is expected to approximate FQM accurately only when the detailed structure of the remaining PESs does not significantly affect the dynamics on Vnn.

Results and Discussion: Tully's model problems. First, the method is tested on three one-dimensional model potentials proposed by Tully11 to cover the most important characteristics of nonadiabatic transitions. Diabatic and adiabatic PESs as well as the coupling terms V12 for Tully’s problems A and C are shown in Figs. 1 and 2; further details can be found in Ref. 11.

FIG. 1.

Fidelity in Tully’s problem A. (a) The diabatic and adiabatic PESs and the diabatic coupling. [(b) and (c)] QM fidelity FQM, the QM survival probability P1,QM on the initial PES, and FDR as functions of time for two different values of initial kinetic energy T0.

FIG. 1.

Fidelity in Tully’s problem A. (a) The diabatic and adiabatic PESs and the diabatic coupling. [(b) and (c)] QM fidelity FQM, the QM survival probability P1,QM on the initial PES, and FDR as functions of time for two different values of initial kinetic energy T0.

Close modal
FIG. 2.

Fidelity in Tully’s problem C.

FIG. 2.

Fidelity in Tully’s problem C.

Close modal

In all cases, the initial wave packet was a Gaussian with approximately 10% dispersion in momentum and located on the lower (problems A and B) or upper (problem C) energy diabatic PES in the asymptotic region without diabatic couplings. The equations of motion were integrated until the wave packet left the interaction region.

The simple avoided crossing model (problem A in the original paper11) represents the most often encountered situation. As can be seen in Fig. 1, the survival probability P1,QM=Trρ̂11(t) on the initial PES is nearly equal to the QM fidelity, suggesting that the fidelity decay is caused almost exclusively by the transition to the second diabatic PES. Close to the diabatic limit, FDR accurately reproduces FQM. In practice, the error of FDR comes from the intrinsic error of the approximation and from the statistical error due to the finite number of trajectories. In this case, the intrinsic relative error is 0.5% when fidelity decays by 10% to 0.9, and even less for the wave packets with higher initial momentum and hence fidelity. The relative statistical error remains below 1% even when only a single trajectory is used, since most of the fidelity decay is due to the transitions to the second surface. For slower wave packets, fidelity decreases, and the agreement between FQM and FDR stays qualitative, with FDR always decaying faster than FQM. Finally, for a fixed initial momentum, decreasing the coupling V12 results in a slower decay of fidelity but the relative error of 1FDR, i.e. (FQMFDR)/(1FQM) is approximately independent of the coupling. Conversely, if V12 is increased substantially, fidelity initially decays quickly to zero, which is captured well by the DR. Subsequently, fidelity may rise again, which is usually not well reflected by the DR.

In the dual avoided crossing model (problem B, not shown), the DR works best again for high energy wave packets. For those, the detailed structure of the second PES does not significantly affect the motion on the first PES. At low energies, where a significant transfer of probability density to the diabatic surface V22 and back occurs, the DR fails to reproduce the QM fidelity even qualitatively. Nevertheless, in this case, both FDR and FQM initially decay almost to zero (although they might rise again later), correctly reflecting that the coupling is important and should not be neglected.

A very interesting situation occurs in the extended coupling model (problem C), where the two diabatic PESs V11 and V22 are almost equal (they differ just by a small constant shift) but the adiabatic surfaces are well separated due to the coupling. At very low energies of the wave packet [Fig. 2(b)], fidelity F goes to zero despite that the survival probability P1 on the upper surface remains close to unity. Therefore, in contrast to previously discussed cases, it is not possible to estimate the extent of the nondiabaticity of QM dynamics on the basis of electronic transitions only. Unlike survival probability, fidelity describes the nondiabaticity effects on the coherent nuclear dynamics correctly. It is remarkable that the DR describes F so accurately despite that the DR dynamics on the diabatic surface ignores the reflection of a large part of the QM wave packet from the upper adiabatic surface. The decrease of fidelity is correctly captured by FDR for all energies of the wave packet. At higher energies [Fig. 2(c)], fidelity converges to zero after several oscillations, whereas the survival probability converges to 1/2. At very high energies (not shown), two wave packets moving on the two PESs interfere for a long time; FDR and FQM agree and oscillate between zero and unity. Surface hopping and Ehrenfest dynamics might predict P1 quite well, but since V11V22, surface hopping would incorrectly predict that FP1. While Ehrenfest dynamics would predict a decay of F below P, it is unlikely that it would capture F quantitatively since the QM wave packet splits into faster and slower components, whereas Ehrenfest dynamics uses a single mean-field surface. One of the reasons why the DR performs so well in the extended coupling model is the similarity of the diabatic PESs in the coupling region. If the surfaces are different, e.g., when V22=cx, the fidelity decay at low energies is still well approximated. However, at higher energies, FDR oscillations slowly dephase from FQM, since the DR neglects effects of V22.

Photodissociation of NaI. We also applied the methodology to the photodissociation of NaI using a two-surface model of Engel and Metiu19 (see Fig. 3). In the original experiment of Mokhtari et al.,20 the molecule was excited by the light in the 310–390 nm range, which led to only a weakly nonadiabatic motion of the wave packet on the excited surface. In this regime (corresponding toT00.07a.u. at the plateau of the PES), far from the diabatic limit, FDR does not quantitatively reproduce FQM. However, when the momentum is approximately ten times higher and the dissociation of NaI is almost diabatic, FDR reproduces FQM rather well [Fig. 3(b) and 3(c)].

FIG. 3.

Fidelity in the photodissociation of NaI.

FIG. 3.

Fidelity in the photodissociation of NaI.

Close modal

Computational details. QM dynamics calculations used the first-order split operator method. Classical trajectories were computed using the first-order symplectic Euler algorithm. The time step varied from 0.02 to 1 a.u., depending on the system and initial momentum.

To guarantee convergence, 16 384 classical paths were used to produce the DR plots. The statistical error of fDR for a two-surface system due to a finite number N of paths is given by

σstat2=1N(cos2(φ)cos(φ)2),
(8)

where φ=ΔS12(x0,t)/ is the phase accumulated along a trajectory and cos(φ)=fDR. As demonstrated on Tully’s model A, where a single trajectory was sufficient, a much lower number of trajectories than 16 384 is needed to obtain an accurate estimate of fidelity in cases where fidelity stays close to unity. Equation (8) implies that for a given value of σstat and fDR, the number N of trajectories needed is approximately independent of dimensionality.

Conclusions. Presented results demonstrate the utility of the DR in analyzing the molecular QM dynamics involving multiple PESs. On one hand, in the nearly diabatic regime, FDR accurately approximates FQM. On the other hand, in systems far from the diabatic limit, FDR decays quickly to zero and thus detects the importance of nondiabatic couplings although it may not reproduce FQM accurately. Hence, the method can be used to establish the level of nondiabaticity of QM dynamics, without the need for a QM dynamics simulation. In fact, we propose the condition F1 (instead of the standard requirement of high survival probability) as the rigorous definition of the diabatic limit (see Fig. 2). Nevertheless, for single avoided crossings, e.g., fidelity could be used to estimate the survival probability and hence the branching ratios (see Figs. 1 and 3).

The DR calculation can be performed easily for all systems accessible to classical molecular dynamics and for which coupling elements Vij,ij are available. However, it remains to be verified how the method will perform in higher-dimensional systems and in systems with more than two important surfaces. The first issue was partially addressed in Ref. 8 where the DR was applied to the two-dimensional photodissociation of CO2, confirming that neither the accuracy nor the number of classical trajectories needed are significantly affected by increased dimensionality.

At the moment, we are exploring a related and complementary problem to that of the present paper, namely, whether the DR in the adiabatic basis can estimate the level of nonadiabaticity of QM dynamics near the adiabatic limit.

Acknowledgments. This research was supported by the Swiss NSF (Grant No. 200021_124936/1) and by the EPFL. We thank Ivano Tavernelli for helpful discussions.

1.
L. J.
Butler
,
Annu. Rev. Phys. Chem.
49
,
125
(
1998
);
[PubMed]
G. A.
Worth
and
L. S.
Cederbaum
,
Annu. Rev. Phys. Chem.
55
,
127
(
2004
);
[PubMed]
M. S.
Child
and
M. A.
Robb
, Eds.,
Non-Adiabatic Effects in Chemical Dynamics
,
Faraday Discussions
Vol.
127
(
Royal Society of Chemistry
,
Cambridge, UK
,
2004
).
2.
A.
Peres
,
Phys. Rev. A
30
,
1610
(
1984
).
3.
J.
Vaníček
,
Phys. Rev. E
70
,
055201
(
2004
).
4.
J.
Vaníček
, e-print arXiv:quant-ph/0410205v1.
5.
J.
Vaníček
,
Phys. Rev. E
73
,
046204
(
2006
).
6.
N. R.
Cerruti
and
S.
Tomsovic
,
Phys. Rev. Lett.
88
,
054103
(
2002
).
7.
J.
Vaníček
and
E. J.
Heller
,
Phys. Rev. E
68
,
056208
(
2003
).
8.
B.
Li
,
C.
Mollica
, and
J.
Vaníček
,
J. Chem. Phys.
131
,
041101
(
2009
).
9.
T.
Zimmermann
,
J.
Ruppen
,
B.
Li
, and
J.
Vaníček
, “
Efficient evaluation of the accuracy of molecular quantum dynamics on an approximate anaytical or interpolated ab initio potential energy surface
,”
Int. J. Quantum Chem.
(unpublished).
10.
J. C.
Tully
and
R. K.
Preston
,
J. Chem. Phys.
55
,
562
(
1971
);
S.
Nielsen
,
R.
Kapral
, and
G.
Ciccotti
,
J. Stat. Phys.
101
,
225
(
2000
);
E. J.
Heller
,
B.
Segev
, and
A. V.
Sergeev
,
J. Phys. Chem. B
106
,
8471
(
2002
).
11.
J. C.
Tully
,
J. Chem. Phys.
93
,
1061
(
1990
).
12.
E. R.
Dunkel
,
S.
Bonella
, and
D. F.
Coker
,
J. Chem. Phys.
129
,
114106
(
2008
).
13.
T. J.
Martínez
,
M.
Ben-Nun
, and
R. D.
Levine
,
J. Phys. Chem.
100
,
7884
(
1996
);
S.
Yang
,
J. D.
Coe
,
B.
Kaduk
, and
T. J.
Martínez
,
J. Chem. Phys.
130
,
134113
(
2009
).
[PubMed]
14.
W. H.
Miller
,
J. Phys. Chem. A
113
,
1405
(
2009
).
15.
M.
Ceotto
,
S.
Atahan
,
G. F.
Tantardini
, and
A.
Aspuru-Guzik
,
J. Chem. Phys.
130
,
234113
(
2009
).
16.
G. C. G.
Waschewsky
,
P. W.
Kash
,
T. L.
Myers
,
D. C.
Kitchen
, and
L. J.
Butler
,
J. Chem. Soc., Faraday Trans.
90
,
1581
(
1994
);
R.
Valero
and
D. G.
Truhlar
,
J. Chem. Phys.
125
,
194305
(
2006
).
[PubMed]
17.
M. H.
Alexander
,
G.
Capecchi
, and
H. J.
Werner
,
Faraday Discuss.
127
,
59
(
2004
).
18.
J.
Vaníček
,
C.
Mollica
,
T.
Prosen
, and
W.
Strunz
(to be published).
19.
M. B.
Faist
and
R. D.
Levine
,
J. Chem. Phys.
64
,
2953
(
1976
);
N.
van Veen
,
M. D.
Vries
,
J.
Sokol
,
T.
Baller
, and
A.
de Vries
,
Chem. Phys.
56
,
81
(
1981
);
V.
Engel
and
H.
Metiu
,
J. Chem. Phys.
90
,
6116
(
1989
).
20.
A.
Mokhtari
,
P.
Cong
,
J. L.
Herek
, and
A. H.
Zewail
,
Nature (London)
348
,
225
(
1990
).