Investigation of geometric phase effects in photodissociation dynamics at a conical intersection

Chemical Physics xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Chemical Physics journal homepage: www.elsevier.com/locate/chemphys

Investigation of geometric phase effects in photodissociation dynamics at a conical intersection Foudhil Bouakline Max Born Institute, Max Born Strasse 2a, 12489 Berlin, Germany

a r t i c l e

i n f o

Article history: Available online xxxx Keywords: Photodissociation Conical intersection Geometric phase Topological unwinding of wavefunctions

a b s t r a c t We investigate the effect of the geometric phase (GP) on photodissociation dynamics at a two-dimensional symmetry-allowed conical intersection (CI). To disentangle the pure effect of the GP from other effects due to non-adiabatic couplings between the two coupled potential energy surfaces, we perform two different calculations, one adopting the diabatic representation which implicitly includes the GP, and another one using the adiabatic picture where GP effects are excluded. To interpret the impact of the GP on nuclear dynamics, we use a recent topological approach (Althorpe et al., 2008 [45]) to completely unwind the nuclear wavefunction from around the CI. This unwinding allows us to extract from the nuclear wavepacket two contributions of reaction paths that wind in different senses around the CI. The solely effect of the GP is to change the sign of the relative phase between their corresponding wavefunctions, and hence to convert any constructive (destructive) interference of the two components, in the asymptotic dissociative limit, into a destructive (constructive) one. This results in a change of the product-state vibrational distribution from only-even (-odd) quanta progression to only-odd (-even) quanta progression. Although our calculations are based on a reduced-dimensionality model Hamiltonian, our observations and conclusions should apply to realistic polyatomic molecules, and could be useful to interpret product-state vibrational distributions of photodissociation experiments. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction Conical intersections (CI) of electronic potential energy surfaces [1] are nowadays widely recognised to be ubiquitous in polyatomic molecules [2–4], and to play a major role in their spectroscopy, photodynamics and also chemical reactivity [4–7]. At a CI, the Born–Oppenheimer approximation (stipulating that electronic and nuclear motions are decoupled) breaks down [8–10]. Molecular systems which exhibit such topologies can easily hop between the two coupled electronic states through the funnel of the CI, giving rise to a rich variety of physico-chemical processes [11–18], such as electron transfer, isomerisation, photoinduced unimolecular decay, and radiationless relaxation of electronic excited states. Another more subtle quantum effect resulting from the presence of conical intersections is the geometric phase (GP) [19,20], also known as the Longuet–Higgins phase, and most commonly as the Berry phase. It is simply the sign change acquired by the electronic wavefunction when the nuclei complete an odd number of loops around the CI. Because the total molecular wavefunction

E-mail address: [email protected]

must be single-valued, the GP introduces a corresponding sign change in the boundary condition of the nuclear wavefunction [21–23], whenever the latter encircles the CI. This has important consequences on molecular spectra and also on molecular collisional processes. For instance, it is well-known that the GP shifts the spectrum of a bound molecular system by altering the pattern of nodes in the nuclear wavefunction [24–31]. Also, substantial geometric phase effects in molecular reactive scattering and photodissociation processes have been reported [32–41]. In order to understand the effect of the geometric phase on non-stationary nuclear wavefunctions, one may adopt a semiclassical picture in which a nuclear wavepacket, approaching the neighbourhood of a conical intersection, splits into two components encircling the CI on opposite sides, and eventually interfere in the asymptotic reactive or dissociative limits. For instance, Köppel and coworkers [42,43] demonstrated numerically that the GP leads to destructive self-interference of the two parts of the wavepacket as it exits the CI region. To further unravel the origin of this destructive interference, and to build a quantum analogue of the above semiclassical picture, Althorpe and coworkers [34,44,45] have recently developed a topological approach to completely unwind the nuclear wavefunction encircling the CI. This method uses simple topological arguments, similar to those

http://dx.doi.org/10.1016/j.chemphys.2014.02.010 0301-0104/Ó 2014 Elsevier B.V. All rights reserved.

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

2

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

previously used in Feynman path integral treatments of the Aharonov–Bohm effect [46–49], to separate the nuclear wavefunction into two contributions, each of which contains all the Feynman paths that wind different numbers of times, and in different senses, around the CI. It is shown that the separation of the wavefunction into even- and odd-looping parts reveals the true effect of the geometric phase on the nuclear wavefunction: It solely changes the sign of the relative phase of these two components. Thus, the recombination and interference of these two parts govern the extent to which dynamical observables are affected by the GP. This unwinding technique applies to the exact nuclear wavefunction, i.e., the paths that it extracts are true Feynman paths, not Newtonian trajectories. To extract the even- and odd-looping components, the only numerical work required is to add and subtract nuclear wavefunctions computed with and without GP boundary conditions. In addition, it has been recently shown that this topological interpretation applies, not only to situations where nuclear dynamics is restricted to the adiabatic electronic ground state, but also to the case where nuclei evolve on both coupled electronic states [45], and could access the conical intersection seam. In this contribution, we aim to exploit this unwinding technique to investigate the effect of the Berry phase in photodissociation dynamics at a symmetry-allowed conical intersection. We use a two-dimensional model Hamiltonian of two coupled electronic states involving one reaction coordinate and one coupling mode. To explore the effect of the geometric phase on the photo-induced nuclear dynamics on the two coupled surfaces, and to disentangle its pure effect from other effects due to non-adiabatic couplings, we perform two sets of calculations, (i) one adopting the diabatic picture which implicitly includes the GP, and (ii) another one using the adiabatic representation where the GP is omitted. Nuclear wavefunctions computed in both representations are used to extract two components that wind in different senses around the CI. These two different reaction paths are then exploited to interpret the observed effects of the GP on the dynamics. The paper is organised as follows. In Section 2, we summarise the basic ingredients of the diabatic and adiabatic representations, used to describe nuclear dynamics at a symmetry-allowed conical intersection, whose main features are given in Section 3. The model Hamiltonian used in our calculations is described in Section 4. In Section 5, we present our results, and discuss the influence of the GP on dynamical observables such as electronic populations, and most importantly its effect on the topology of the wavepacket in the dissociative limit, and on the subsequent product-state vibrational distributions. Finally, in Section 6, we use the topological unwinding technique to interpret the observed effects of the GP on the dynamics. Section 7 concludes the paper.

We want to describe nuclear motion on two coupled electronic potential energy surfaces. To do so, two representations are invoked, the adiabatic and diabatic ones. Within the adiabatic picture, the nuclear Hamiltonian has the following form [50]



Tn þ V

0

0

Tn þ V þ





þ

K11 K21

K12 K22



ð1Þ

where T n represents the nuclear kinetic energy operator, and V  and V þ are the potential energies of the ground and excited adiabatic electronic states, respectively. These states are coupled through the non-adiabatic matrix elements fKij g, which are given in mass-scaled coordinates by

Kij ¼ 

h2 ð2Fij  $ þ Gij Þ 2M

Hd ¼ S Had Sþ

ð3Þ

where the transformation matrix is given by

 S¼

 cosðHÞ sinðHÞ :  sinðHÞ cosðHÞ

ð4Þ

The adiabatic–diabatic mixing angle H, which is a function of the nuclear coordinates {Q}, must be chosen to remove the off-diagonal non-adiabatic kinetic couplings. Although exact diabatic electronic states do not strictly exist [51], several diabatisation schemes [52] have proven to be very efficient to remove the divergent part of the off-diagonal couplings, and to ensure that the residual derivative couplings are vanishingly small. Using these diabatic construction procedures, kinetic derivative couplings are removed and transformed into smooth potential energy couplings, giving rise to the following form of the diabatic nuclear Hamiltonian

Hd ¼ T n



1 0 0 1



 þ

U1

U 12

U 21

U2

 :

ð5Þ

The adiabatic and diabatic pictures are related through the mixing angle

HðQ Þ ¼

1 2U 12 ðQ Þ arctan : 2 U 2 ðQ Þ  U 1 ðQ Þ

ð6Þ

This transformation angle allows to avoid the numerical calculation of the non-adiabatic kinetic couplings, and facilitates their evaluation analytically. After some simple algebra [50,52,53], we can show that

F ¼ $HðQ Þ G ¼ ð$  FÞ þ F  F:

ð7Þ

Although the adiabatic and diabatic representations are linked by a unitary transformation, and hence, should give the same nuclear measured observables, a subtle difference emerges in the presence of a conical intersection. The geometric phase is implicitly taken into account only in the diabatic picture. In the adiabatic representation, to describe nuclear dynamics correctly, the GP has to be included explicitly on both coupled electronic states. 3. Symmetry-allowed conical intersections

2. Quantum nuclear dynamics in the adiabatic and diabatic representations

Had ¼

where Fij ¼ hUi j$jUj i is the (vectorial) first derivative coupling matrix element in the adiabatic electronic basis fUi g, and Gij ¼ hUi jr2 jUj i is the (scalar) second derivative coupling matrix element. Numerical wavepacket propagation in the adiabatic picture is potentially difficult because the off-diagonal non-adiabatic couplings become singular at the CI seam. To overcome these numerical problems, one can convert to a (quasi-) diabatic representation of the wavefunction via a unitary transformation [50]

ð2Þ

Conical intersections are widespread in real polyatomic molecules [4]. Their occurrence requires two different conditions to be fulfilled, the diagonal elements of the diabatic potential matrix given in Eq. (5) must be equal, and the off-diagonal elements must vanish. This gives rise to a conical intersection seam of a dimension ðN  2Þ in the nuclear coordinate space, where N is the number of the internal nuclear degrees of freedom. In general, there is no symmetry element that would determine the location of conical intersections in the nuclear coordinate space. However, there are some specific cases where symmetry plays an important role in the characterisation of CIs. Here, we consider the case of the socalled symmetry-allowed conical intersections [5,8]. They arise from the interaction of two nondegenerate electronic states with different spatial symmetries. For these states to interact, a nuclear distortion of the molecular system along a nonsymmetric nuclear mode is required.

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

In this section, we present a summary of the main features of symmetry-allowed conical intersections described in more details in Ref. [8]. Let us consider two nondegenerate electronic states U1 and U2 having two different spatial symmetries labelled by A1 and A2 , respectively. The product of the latter (A1  A2 ) is denoted by B, and the totally symmetric representation of the point group of the molecular system by A. Normal mode analysis allows to determine the symmetries of the different nuclear coordinates (modes) relevant to the dynamics around the CI. It is clear from symmetry requirements that the two electronic states can only couple through the nontotally symmetric coordinates, called coupling modes [8], which transform as B. The totally symmetric coordinates transforming as A, called tuning modes [8], contribute to the regulation of the energy gap between the two interacting electronic states. Building upon this symmetry analysis, one can easily derive the analytical forms of the diabatic potentials and their coupling. Let us denote the tuning modes by Q 1 , the coupling modes by Q 2 , the remaining other coordinates by Q 3 , and consider one mode of each for the simplicity of the analysis. Knowing that the nuclear Hamiltonian and the kinetic energy operators are totally symmetric, it is obvious that the diabatic potentials, U 1 ¼ hU1 jHd  Tn jU1 i and U 2 ¼ hU2 jHd  Tn jU2 i, are totally symmetric functions of the nuclear coordinates: U 1 ¼ U 1 ðQ 1 ; Q 22 ; Q 23 Þ and U 2 ¼ U 2 ðQ 1 ; Q 22 ; Q 23 Þ. On the other hand, the off-diagonal diabatic potential matrix element U 12 ¼ hU1 jHd  Tn jU2 i transforms as the nontotally symmetric representation B, and has to be an odd function of the coupling mode Q 2 . We can thus write it as U 12 ¼ Q 2 kðQ 1 ; Q 22 ; Q 23 Þ. The diabatic potentials and their coupling can be expanded in Taylor series in the near vicinity of a reference point such as the ground state equilibrium geometry. One can also use the geometry of the conical intersection to perform this development. To the first order of this expansion around the reference point, these quantities read

U i ¼ V i þ ki Q 1 ;

i ¼ 1; 2

ð8Þ

U 12 ¼ U 21 ¼ lQ 2 :

V i are the values of the diabatic potentials U i at the reference point, and ki and l are the slopes at the reference point of the diabatic potentials and their coupling, respectively. From the Taylor development, we can see that the totally symmetric coordinate (the tuning mode) enters the diagonal diabatic potentials linearly, whereas the nontotally symmetric coordinate (the coupling mode) appears linearly in the off-diagonal coupling. A closer look at Eq. (8) shows that the adiabatic potential energy surfaces, which are the eigenvalues of the diabatic potential matrix, form a conical intersection in the fQ 1 ; Q 2 g plane at the position ðQ 1 ¼ ½V 2  V 1 =½k1  k2 ; Q 2 ¼ 0Þ. Notice that nuclear modes of type Q 3 are irrelevant in the neighbourhood of the conical intersection. One of the subtle consequences of CIs is the geometric phase [19,20], which is the sign change acquired by the adiabatic electronic states when the nuclei complete an odd number of loops around the CI. This can be easily shown for the case of a symmetry-allowed conical intersection. After a transformation of the system coordinates fQ 1 ; Q 2 g to the polar coordinates system fQ ; ag centred around the CI, diagonalisation of the diabatic potential matrix of Eq. (8) yields the following adiabatic electronic eigenstates [19]



jU i jUþ i



 ¼

cosða=2Þ

sinða=2Þ

 sinða=2Þ cosða=2Þ



jU1 i jU2 i

 ð9Þ

A comparison between Eqs. (4) and (9) shows that, close to the CI region, the pseudorotation angle a and the adiabatic-diabatic mixing angle H are related through a ¼ 2H. Further inspection of Eq. (9) indicates that the geometrical operation a ! a þ 2p, describing

3

a closed loop around the CI, changes the sign of the adiabatic electronic states jU i. Since the overall molecular wavefunction is single-valued, a compensatory sign change in the nuclear wavefunction has to be introduced whenever the nuclei encircles the CI an odd number of times. Hence, a correct description of nuclear dynamics at a conical intersection within the adiabatic representation requires the explicit inclusion of the geometric phase, either by using nuclear bases and coordinate systems satisfying the GP boundary condition, such as hyperspherical coordinates [26,38], or by using complex electronic wavefunctions which cancel the sign change introduced by the geometric phase, but introduce an additional vector potential in the nuclear Hamiltonian [21,38]. We should emphasise that the Berry phase only appears when using the adiabatic picture. Adoption of the diabatic representation ensures correct boundary conditions because the GP is exactly included through the adiabatic-diabatic mixing angle. 4. Model Hamiltonian and computational details In this section, we describe a simple example of a photodissociation process in the presence of a symmetry-allowed conical intersection between two electronic states. For simplicity of the analysis, we consider a two-dimensional model including one totally symmetric mode and one coupling mode. The diabatic Hamiltonian used in our study is taken from the work of Domcke and coworkers [54,55], featuring the photoinduced dynamics of pyrrole at the 1 A2  S0 conical intersection using a reduced dimensionality model, where only two nuclear modes are taken into account. The mode Q 1 which plays the role of the hydrogen abstraction reaction coordinate, and one normal mode of the pyrrole ring Q 2 , which is the coupling coordinate. In this model, the kinetic energy operator takes the following form

Tn ¼ 

h2 @ 2 1 @2  x 2l @Q 21 2 @Q 22

ð10Þ

where l is the reduced mass corresponding to the motion of the hydrogen atom relative to the ring part of pyrrole. Q 1 is defined as the distance between the H atom and the centre of mass of the ring, and Q 2 is a dimensionless normal mode of the ring with the corresponding frequency x. Both the ground state and excited state diabatic potentials are taken to be harmonic along the coupling coordinate Q 2 . Along the reaction coordinate Q 1 , the ground state is described with a Morse oscillator potential, and the excited state is purely repulsive with a low potential barrier, resulting from its interaction with another excited state at small Q 1 distances. The coupling between the two potentials depends on both coordinates, and is an odd function of the coupling mode Q 2 , i.e., U 12 ¼ Q 2 kðQ 1 Þ. The details of the diabatic potentials can be found in the appendix. Fig. 1 illustrates a one-dimensional cut along the reaction coordinate of the two diabatic potentials, as well as the photoinduced dynamical process under study. An initial wavepacket is created by promoting a vibrational eigenstate of the electronic ground state into the electronic excited state, through a vertical electronic transition. The created wavepacket then moves towards the conical intersection where it separates into two pieces, one staying on the excited state, and another one hopping to the ground state under the influence of non-adiabatic couplings. Most of the wavepacket will eventually dissociate into the exit channels of the two potentials. The photoinduced nuclear dynamics of this model system is treated in the time-dependent picture, where the Schrödinger equation for the nuclear motion on the two coupled potential energy surfaces is solved using the split-operator method [56]. The Hamiltonian and the wavepacket are represented in a basis of distributing approximating functional (DAF) grid functions [57], in

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

Fig. 1. One-dimensional potential energy profiles, along the dissociation coordinate, of the diabatic ground state and the diabatic excited state. The symbol () indicates the location of the conical intersection.

which any local operator is diagonal. The action of the kinetic energy operator on the wavepacket is performed with fast fourier transform (FFT) routines. In the diabatic representation, wavepacket propagation with the split-operator technique necessitates the diagonalisation of the potential energy matrix (Eq. (5)), which we evaluate analytically [58]. On the other hand, in the adiabatic representation, a numerical diagonalisation of the non-adiabatic coupling matrix (Eq. (1)) is required. To avoid unphysical reflections of the wavepacket at the edge of the Q 1 grid, a complex absorbing potential [59] is added to the Hamiltonian in the asymptotic dissociation region. To circumvent the singularity of the non-adiabatic kinetic couplings at the conical intersection, we choose grid functions such that no corresponding grid point coincides with the CI position. A similar technique is often used in atom–diatom scattering [60], where the centrifugal potential in Jacobi coordinates becomes singular when the atom position coincides with the diatom centre of mass. Convergence of the calculations with respect to different parameters, such as the number of grid points along both coordinates, the propagation time, and the time step, has been checked. In order to gauge the effect of the Berry phase on nuclear dynamics, and disentangle its pure effect from other non-adiabatic effects, we compute, in both the diabatic and adiabatic representations, time-dependent adiabatic electronic populations, adiabatic nuclear probability densities, and dissociation probabilities. The latter are defined as the time-accumulated nuclear flux through a dividing surface located in the dissociation asymptotic region at Q 1 ¼ Q flux :

PD ðsÞ ¼

h

Z s

l

0

    @ WðQ 1 ; Q 2 ; tÞ Im WðQ 1 ; Q 2 ; tÞ dt @Q 1 Q 1 ¼Q flux

ð11Þ

where jWðQ 1 ; Q 2 ; tÞi is the time-dependent nuclear wavefunction computed either in the diabatic or adiabatic representation.

vibrational eigenstate ðv Q 1 ¼ 1; v Q 2 ¼ 0Þ is chosen to be the initial wavepacket. Fig. 2 illustrates the time evolution of the electronic population of the adiabatic states, computed in both the diabatic and adiabatic pictures. Both calculations indicate a rapid relaxation of the excited electronic state in the early stage of the dynamics. An extensive population transfer from the excited state to the ground state occurs within only 30 fs, followed by a slow and monotonic decay, and then reaching a plateau at about 150 fs. This fast population transfer is due to the presence of the conical intersection. As the wavepacket approaches the neighbourhood of the CI, it starts to experience the effects of the non-adiabatic couplings between the two coupled surfaces, inducing the nuclear wavefunction to hop to the electronic ground state. At latter stages of the dynamics, most of the wavepacket will dissociate into the lower dissociation limit, and only a small part exits to the upper dissociation channel. In our calculations, both the diabatic and adiabatic representations show qualitatively the same relaxation dynamics of the excited state. However, the geometric phase, only included in the diabatic picture, seems to enhance population transfer as it is shown in Fig. 2. Though, this enhancement could be only specific to our model system, and future work on polyatomic molecules is needed to determine the extent to which the GP hinders or enhances relaxation through conical intersections. To get more detailed insight into the effect of the geometric phase on nuclear dynamics, we follow the time evolution of the nuclear wavepacket on both surfaces, and closely inspect its topology. Figs. 3 and 4 show snapshots of the projections of the nuclear probability densities of the adiabatic excited state and ground state, respectively. These probability densities are illustrated at two different times, t = 10 and 40 fs, and computed in both the diabatic (panels (a) and (b)) and adiabatic (panels (c) and (d)) representations. The dashed horizontal and vertical lines in the figures are shown to locate the position of the CI, and also for further symmetry analysis. At 10 fs, the initial wavepacket has already reached the conical intersection, and a piece of it is transferred to the electronic ground state. At 40 fs, a large part of the wavepacket has escaped the neighbourhood of the conical intersection moving towards the two dissociation exit channels. We observe that wavepackets, computed in the diabatic and adiabatic pictures, do not have the same topology in the vicinity of the CI. The wavefunctions obtained in the adiabatic picture are more delocalised along the coupling mode. This is not a surprise knowing the shapes of the potentials

1

0.8

electronic population

4

0.6

0.4

0.2

5. Results and discussion We first investigate the effect of the geometric phase on the time-dependent populations of the adiabatic electronic states. The initial wavepacket is taken to be an eigenstate of the adiabatic electronic ground state, vertically promoted to the excited electronic state. Throughout this work, the electronic ground state

0

0

50

100

150

200

time (fs) Fig. 2. Population of the lower and upper adiabatic electronic states, computed in the diabatic representation (solid and dashed lines), and in the adiabatic representation (chained and dotted lines).

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

5

Fig. 3. Snapshots of the nuclear probability density in the adiabatic excited state taken at (t = 10 and 40 fs), and computed in the diabatic representation (a and b), and in the adiabatic representation (c and d). The intersection of the dashed vertical and horizontal lines indicates the position of the conical intersection . The solid vertical lines in (b and d) show the Q 1 position at which cuts along Q 2 are made and displayed in Fig. 5.

Fig. 4. Same as in Fig. 3 but for the adiabatic ground state.

and their couplings in this region of space. The diabatic potentials vary smoothly in the whole coordinate space, whereas the corresponding adiabatic potentials (with their double cone shape close to the CI) abruptly change in this region. Moreover, the diabatic matrix coupling U 12 is a local operator of coordinates, while the non-adiabatic coupling K12 is nonlocal. These arguments explain the enhancement of the wavepacket spreading in the adiabatic representation. However, this is not the main important conclusion we draw from Figs. 3 and 4. The most striking observation is about the change in the symmetry of the wavepacket, with respect to the

(Q 2 ¼ 0) plane, when it exits the CI vicinity towards the dissociation asymptotic limits. Let us first consider the results of the diabatic representation, given in Figs. 3b and 4b. The wavepacket transferred to the adiabatic ground state keeps the same symmetry of the initial wavepacket. It is totally symmetric with respect to the Q 2 coordinate. On the other hand, the piece of the wavepacket evolving towards the upper dissociation limit is anti-symmetric, it has a nodal line at Q 2 ¼ 0. Inspection of the adiabatic picture results, given by Figs. 3d and 4d show the opposite trend. The component of the wavepacket evolving on the adiabatic ground state acquires a node at Q 2 ¼ 0, and the other part moving on the

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

6

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

probability density

adiabatic excited state is totally symmetric. To highlight this change of the nodal structure of the wavepacket, we show in Fig. 5 a one-dimensional cut, along the coupling coordinate, of the nuclear probability densities in the asymptotic region. We observe that the bimodal structure of the excited state nuclear probability density, obtained in the diabatic representation (Fig. 5a), is turned into a bell-shaped distribution in the adiabatic picture. The opposite trend happens for the transferred wavepacket into the adiabatic ground state as it is shown in Fig. 5b. We should emphasise that this characteristic nodal structure of the wavepacket in the asymptotic region is observed at all propagation times, and occurs whenever the system passes through the CI towards the dissociation region. The nodal structure of the wavepacket, along the coupling coordiante Q 2 , can be easily interpreted using simple symmetry arguments in the diabatic picture. To ensure dissociation to the lower exit channel, the wavepacket must remain on the same diabatic electronic excited state, which has a well-defined electronic character and symmetry. As the initial wavepacket is an even function of the coupling coordinate Q 2 , it must stay so when it passes through the conical intersection to maintain the same overall vibronic symmetry. On the other hand, dissociation into the upper channel results from the hopping of the initial wavepacket from the diabatic excited state to the diabatic ground state. This surface hopping is induced by the off-diagonal potential matrix element U 12 , which is an odd function of the Q 2 coordinate. Thus, the transferred wavepacket must gain a node along the Q 2 ¼ 0 plane, where the coupling U 12 is zero. This interpretation is based only on symmetry requirements, and could thus be applied to the adiabatic picture. However, the latter shows an opposite wavefunction node pattern to the one observed in the diabatic picture. We will show in the next section that this specific nodal structure of the wavepacket is a clear signature of the geometric phase effect, which is missing in the adiabatic representation. As the wavefunction itself is not an observable, one can argue that the observed nodal structure of the wavepacket is only an artefact of the formalism, and could not be experimentally verified. We show that probe of the dissociation products vibrational distribution, directly reveals the nodal structure of the wavepacket after its passage through the CI. Figs. 6 and 7 show product-state vibrational populations in the upper and lower dissociation channels,

0.0006

(a)

diabatic adiabatic

0.0004

0.0002

0

-4

-2

0

2

4

probability density

Q2 (au)

(b)

diabatic adiabatic

0.0009 0.0006 0.0003 0

-4

-2

0

2

4

Q2 (au) Fig. 5. Cuts of the asymptotic nuclear probability density taken along the coupling coordinate (Q 2 ) at (t = 40 fs), (a) in the upper adiabatic state (at Q 1 ¼ 7:6 au), and (b) in the lower adiabatic state (at Q 1 ¼ 12:4 au). The plots show the results obtained in both the diabatic and adiabatic representations.

respectively. Dissociation probabilities are normalised to the value of the lowest populated vibrational state. We clearly see that the diabatic picture calculations show an odd vibrational quanta progression in the upper dissociation limit. This is a clear proof of a node along the coupling coordinate of the corresponding wavepacket. On the other hand, products collected in the lower dissociation channel show an even quanta progression, which is a clear signature that the vibrational symmetry of the initial wavepacket is preserved throughout the dissociation process. However, the adiabatic picture calculations show the opposite trend, an even quanta progression for products in the upper dissociation limit, and an odd quanta progression in the lower dissociation limit. What is shown in Figs. 6 and 7 is a direct ‘experimental’ mapping of the nodal structure of the wavepacket observed in Figs. 3–5. This selective formation of the products in a specific vibrational quanta progression of the coupling mode, only even or only odd, has been experimentally observed in UV photoinduced dissociation of the O-H bond of phenol [63,64]. The authors interpreted this observation as a dynamical effect of the geometric phase. We should also mention that the change in the nodal structure of the wavepacket passing through a conical intersection has been numerically observed in several works [42,61,62]. However, most of these studies use either symmetry arguments or semiclassical interpretations to explain the effect of the GP on nuclear photodynamics at a conical intersection. Our aim in the next section is to go beyond, and use topological arguments to completely unwind the nuclear wavefunction at the CI, and thus explain the observed discrepancies between the results obtained in the adiabatic and diabatic representations. 6. Topological unwinding of the nuclear wavefunction at a conical intersection In this section, we use simple topological arguments to explain the change in the nodal structure of the wavepacket as it passes through the conical intersection. The time evolution of a quantum system is completely defined by its time-dependent evolution operator or kernel. In Feynman path integral theory, the latter is given by

Kðx; x0 jtÞ ¼

Z

DxðtÞeiSðx;x0 jtÞ=h

ð12Þ

where DxðtÞ defines a sum over all Feynman paths connecting the points x0 and x in the time interval t, and S is the classical action evaluated along each of these paths. In the language of time-dependent wavefunctions, the integral in Eq. (12) has to be performed in the configuration space over all possible paths connecting every point of the initial wavepacket to every point of the final wavepacket. In general, all these Feynman paths are coupled, meaning that if we start with a chosen path, we can go to any other path connecting the same initial and final points. This is achieved by applying tiny distortions to the first path, passing through its immediate neighbours to reach the other final path. In other words, any Feynman path can be continuously deformed to another path connecting the same x0 and x points in the time interval t. This is not exactly the case for molecular systems with conical intersections. The latter geometries is an example of what is called, in the language of topology, a multiply connected space. This means that there is an obstacle or an inaccessible region of space, in the vicinity of the CI, preventing the continuous deformation of Feynman paths. A given path is coupled to only a subset of paths belonging to the same topological homotopy class. The latter contains all Feynman paths connecting the same initial and final points within the same time interval t, but making the same number of loops (and in the same sense) around the conical intersection. Thus, each homotopy class is

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

7

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

1 dissociation probability

diabatic 0.8 0.6 0.4 0.2 0 0

1

2

3

4

5

6

7

8

9

10

1 dissociation probability

adiabatic 0.8 0.6 0.4 0.2 0 0

1

2

3

4

5 v

6

7

8

9

10

Fig. 6. Product-state vibrational distribution in the upper dissociation channel, computed in the diabatic and adiabatic representations. Vibrational populations are normalised to the first populated vibrational state.

1 dissociation probability

diabatic 0.8 0.6 0.4 0.2 0 0

1

2

3

4

5

6

7

8

9

10

1 dissociation probability

adiabatic 0.8 0.6 0.4 0.2 0 0

1

2

3

4

5 v

6

7

8

9

10

Fig. 7. Same as in Fig. 6 but for the lower dissociation channel.

characterised by its winding number n. All Feynman paths can be further grouped into two distinct classes, one corresponding to all even n, and another to all odd n. It has been demonstrated [44] that the solely effect of the geometric phase on the evolution operator is to change the relative sign between the even- and odd-looping kernels K e ðx; x0 jtÞ and K o ðx; x0 jtÞ

K N ðx; x0 jtÞ ¼ K e ðx; x0 jtÞ þ K o ðx; x0 jtÞ K G ðx; x0 jtÞ ¼ K e ðx; x0 jtÞ  K o ðx; x0 jtÞ:

ð13Þ

K N ðx; x0 jtÞ and K G ðx; x0 jtÞ define the kernels computed with normal and GP boundary conditions, respectively. Mapping these kernels to their corresponding nuclear wavefunctions, we get

1

WN ¼ pffiffiffi ½We þ Wo  2 1 p ffiffiffi ½We  Wo  WG ¼ 2

ð14Þ

where the factor p1ffiffi2 is introduced for normalisation. We should mention that this unwinding technique has been first developed for the case where nuclear dynamics is restricted to the adiabatic electronic ground state PES [34], but it has been recently shown [45] that it can be safely applied also to the case where the wavepacket evolves on two coupled electronic states, and has access to the CI seam.

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

8

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

Fig. 8. Schematic diagram showing the effect of the geometric phase on the nodal structure of the wavepacket along the coupling mode.

0.002

(a)

even odd

0.001

0

-4

-2

0

2

constructive interference into a destructive one. This gives rise to a nodal line in the wavefunction along the coupling coordinate as it is illustrated in Fig. 8. On the other hand, if the interference between the even- and odd-looping paths is destructive when using normal boundary conditions, the GP will convert it into a constructive interference. To put the unwinding technique into practice, we apply it to our model Hamiltonian in order to explain the nodal structure of the wavepacket observed in Figs. 3–5. Nuclear probability densities are obtained from the coherent sums of the even- and odd-looping nuclear wavefunctions

1 ½jWe j2 þ jWo j2  þ jWe jjWo j cosð/e  /o Þ 2 1 jWN j2 ¼ ½jWe j2 þ jWo j2   jWe jjWo j cosð/e  /o Þ 2

jWN j2 ¼

0.004

probability density

probability density

Photodissociation at a conical intersection provides an ideal process to check the applicability of this topological approach in the case of two coupled electronic states. As the initial wavepacket approaches the neighbourhood of the CI, it bifurcates into two pieces, each of which going around the CI either in the clockwise sense, or in the counterclockwise sense, and then both parts proceed to the dissociation asymptotic channels. These two components correspond to the even- and odd-looping Feynman paths introduced above. Although the two separate pieces of the wavepacket can hardly make loops around the CI because they both quickly escape this region towards the dissociation channels, the asymptotic recombination of the two bifurcated parts ensures a complete encirclement of the CI as it is shown in Fig. 8. The solely effect of the geometric phase is to change the sign of the overlap integral between the two components, and thus converting their

(c)

even odd

0.002

0

4

-4

-2

Q2 (au)

2

4

2

4

1

(d)

(b) 0.5

cos (φe - φo)

cos (φe - φo)

0

Q2 (au)

1

0 -0.5 -1

ð15Þ

-4

-2

0

Q2 (au)

2

4

0.5 0 -0.5 -1

-4

-2

0

Q2 (au)

Fig. 9. Unwinding of the nuclear wavefunctions displayed in Fig. 5, into even- and odd-looping components. The plots show their nuclear probability densities, and their relative phases. The graphs (a,b) correspond to the upper adiabatic state, and (c,d) to the lower adiabatic state.

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

9

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

where /e and /o are the phases of the even- and odd-looping components of the wavefunctions, respectively. From the last formula, it is clear that, for noticeable GP effects to be observed, a strong overlap between the even- and odd-looping parts is required. Eq. (14) shows that the extraction of the clockwise (counterclockwise) component of the wavefunction is simply achieved by coherently adding (subtracting) the adiabatic wavefunctions computed in the diabatic and adiabatic picture, respectively. Fig. 9 illustrates the application of the unwinding technique to the results shown in Fig. 5. For the upper adiabatic state, the corresponding probability densities of the even- and odd-looping components, and their relative phase are shown in panels (a) and (b), respectively. Panels (c) and (d) show the same, but for the nuclear wavepacket evolving on the lower adiabatic state. Our first observation is that the probability densities of the even- and odd-looping parts of the wavepacket are mirror images of each other with respect to the coupling mode symmetry plane Q 2 ¼ 0. Thus, at this line, there are as many even-looping Feynman trajectories as odd-looping ones, and on both potential energy surfaces, as shown in panels (b) and (d) of Fig. 9. This is a simple consequence of the symmetry of the potentials and the initial wavepacket with respect to Q 2 ¼ 0. The only difference between the two components resides in their phases. On the excited electronic state, Fig. 9b shows that both components are in-phase at the coupling coordinate symmetry plane, giving rise to a complete constructive interference when using normal boundary conditions. The geometric phase, by simply changing the sign of the overlap term (Eq. (15)), turns this interference into a complete destructive one. This is why we observe a nodal line along Q 2 ¼ 0 in the nuclear wavepacket evolving on the upper adiabatic state, when adopting the diabatic picture. On the other hand, the evenand odd-looping components of the nuclear wavefunction on the lower adiabatic state are out-of-phase at Q 2 ¼ 0, inducing a complete destructive interference along the symmetry plane, in the adiabatic representation. The GP changes again the sign of the overlap between the two components, leading to a complete constructive interference.

7. Conclusion In this contribution, we have investigated the effect of the Berry phase in photodissociation dynamics at a symmetry-allowed conical intersection using a simple two-dimensional model Hamiltonian. The main result of this work is the observation of a change in the nodal structure of the nuclear wavefunction as it exits the conical intersection region. The vibrational symmetry of the wavepacket along the coupling mode is preserved when dissociation proceeds to the lower dissociation channel, but it changes when the wavefunction passes through the CI and dissociates into the upper dissociation limit. This effect has important consequences on the vibrational spectra of the dissociation products. The latter are formed with a specific quanta progression in the coupling mode vibration. Depending on the symmetry of the initial wavepacket and the dissociation channel, products with only-odd or only-even quanta in the coupling mode vibration are observed. This effect can be easily verified experimentally, either by probing the vibrational spectra of the photodissociation products, as in recent experiments on the photodissociation of phenol [63,64], or by using pump–probe experiments to analyse the topology of the wavepacket just at the exit of the CI region. We have showed that this main observation is a clear signature of the geometric phase effect, which ensures preservation of the overall vibronic symmetry of the wavefunction. To support our interpretation of these findings, we have used simple topological arguments to completely unwind the nuclear wavepacket around the conical intersection. Two reaction paths, corresponding to the

Table 1 Values of the parameters for the diabatic model Hamiltonian. U2

U1

U 12

D2 ¼ 4:979 eV a1 ¼ 1:293 au r1 ¼ 1:882 au E1 ¼ 4:805 eV b ¼ 2:644 au b2 ¼ 1:325 au r 2 ¼ 2:216 au E2 ¼ 3:956 eV v 12 ¼ 1:248 eV x2 ¼ 0:1096 eV x ¼ 0:0825 eV

D1 ¼ 4:979 eV a0 ¼ 1:137 au r 0 ¼ 1:927 au x2 ¼ 0:1096 eV x ¼ 0:0825 eV

k12 ¼ 0:237 eV a12 ¼ 1:369 au r 12 ¼ 3:679 au

even- and odd-looping Feynman paths, are extracted simply by adding and subtracting nuclear wavefunctions computed in the diabatic and adiabatic representations. The change in the nodal structure of the nuclear wavepacket along the coupling mode symmetry line is brought about by two factors: (1) Because of the symmetry of the initial wavepacket and the potentials with respect to the coupling coordinate, The odd- and even-looping components of the wavefunction have the same amplitude at this symmetry line, where they are always either in-phase or out-of-phase. (2) The solely effect of the GP is to change the relative phase of the two components from being in-phase to out-of-phase, and vice versa. Finally, we should emphasise that these important findings are not restricted to the reduced dimensionality model Hamiltonian used here, but also apply to polyatomic molecules. For the latter, a change in the nodal structure of the wavepacket will not occur only along one coupling coordinate, but in the subspace of all the coupling modes contributing to the dynamics in the vicinity of the CI. In addition, the observed self-destructive (or self-constructive) interference of the wavepacket is not only limited to symmetryallowed CIs, but should also occur for other types of conical intersections. Despite using a simple model Hamiltonian, we believe that this work will stimulate further theoretical and experimental investigations of the observed effects in the photodissociation of realistic polyatomic molecular systems. Acknowledgements The author would like to thank Stuart Althorpe for introducing him to the field of non-adiabatic chemistry, and Olga Smirnova for fruitful discussions. Appendix A In this section, we describe the two-dimensional Hamiltonian (due to Vallet et al. [54]) used in this work. The diabatic potential energy surfaces are modelled by harmonic functions of the coupling coordinate Q 2 and analytical functions of the NH coordinate r (i.e., the distance from the H atom to the N atom of the pyrrole ring). The latter is related to the reaction coordinate Q 1 used in dynamical calculations by (r ¼ Q 1  q0 ), where q0 is the distance between the centre of mass of the ring and the N atom at the equilibrium geometry. The diabatic potentials are given by

1 U 1 ðr; Q 2 Þ ¼ D1 ½1  expða0 ðr  r0 ÞÞ2 þ x1 Q 22 2  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 U 2 ðr; Q 2 Þ ¼ ðv 1 ðrÞ þ v 2 ðrÞ  ððv 1 ðrÞ  v 2 ðrÞÞ2 þ 4v 212 Þ þ x2 Q 22 2 2

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010

ð16Þ

10

F. Bouakline / Chemical Physics xxx (2014) xxx–xxx

where 2

v 1 ðrÞ ¼ D2 ½1  expða1 ðr  r1 ÞÞ v 2 ðrÞ ¼ b expðb2 ðr  r2 ÞÞ þ E2 :

þ E1

ð17Þ

The coupling between the two diabatic potentials is linear in the coupling coordinate

U 12 ðr; Q 2 Þ ¼

   1 r  r 12 Q 2: k12 1  tanh a12 2

ð18Þ

The numerical values of the parameters of this diabatic model are given in Table 1. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

H.C. Longuet-Higgins, Proc. R. Soc. Lond. A 344 (1975) 147. A.G. Atchity, S.S. Xantheas, K. Ruedenberg, J. Chem. Phys. 95 (1991) 1862. D.R. Yarkony, Acc. Chem. Res. 31 (1998) 511. W. Domcke, D.R. Yarkony, H. Köppel (Eds.), Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, World Scientific, River Edge, NJ, 2004. G.A. Worth, L.S. Cederbaum, Annu. Rev. Phys. Chem. 55 (2004) 127. S. Mahapatra, Int. Rev. Phys. Chem. 23 (2004) 483. W. Domcke, D.R. Yarkony, Annu. Rev. Phys. Chem. 63 (2012) 325. H. Köppel, W. Domcke, L.S. Cederbaum, Adv. Chem. Phys. 57 (1984) 59. M. Baer, Beyond Born-Oppenheimer: Electronic Non-Adiabatic Terms and Conical Intersections, Wiley, Hoboken, NJ, 2006. H. Nakamura, Nonadibatic Transition: Concepts, Basic Theories and Applications, World Scientific, Singapore, 2012. D.A. Micha, I. Burghardt (Eds.), Quantum Dynamics of Complex Molecular Systems, Springer, Heidelberg, 2006. S. Matsika, P. Krause, Annu. Rev. Phys. Chem. 62 (2011) 621. M. Barbatti, A.J.A. Aquino, J.J. Szymczak, D. Nachtigallova, P. Hobza, H. Lischka, Proc. Natl. Acad. Sci. USA 107 (2010) 21453. S. Mahapatra, Acc. Chem. Res. 42 (2009) 1004. B.G. Levine, T.J. Martínez, Annu. Rev. Phys. Chem. 58 (2007) 613. I. Burghardt, L.S. Cederbaum, J.T. Hynes, Comput. Phys. Commun. 169 (2005) 95. T.A. Barckholtz, T.A. Miller, Int. Rev. Phys. Chem. 17 (1998) 435. F. Bernardi, M. Olivucci, M.A. Robb, Chem. Soc. Rev. 25 (1996) 321. G. Herzberg, H.C. Longuet-Higgins, Discuss. Faraday Soc. 35 (1963) 77. M.V. Berry, Proc. R. Soc. London Ser. A 392 (1984) 45. C.A. Mead, D.G. Truhlar, J. Chem. Phys. 70 (1979) 2284. C.A. Mead, Rev. Mod. Phys. 64 (1992) 51. M.S. Child, Adv. Chem. Phys. 124 (2002) 1. C.A. Mead, Chem. Phys. 49 (1980) 23. F.S. Ham, Phys. Rev. Lett. 58 (1987) 725. B. Lepetit, Z. Peng, A. Kuppermann, Chem. Phys. Lett. 166 (1990) 572. B.K. Kendrick, Phys. Rev. Lett. 79 (1997) 2431. A.J.C. Varandas, H.G. Yu, J. Chem. Soc. Faraday Trans. 93 (1997) 819.

[29] S. Mahapatra, H. Köppel, Phys. Rev. Lett. 81 (1998) 3116. [30] H. von Busch, V. Dev, H.-A. Eckel, S. Kasahara, J. Wang, W. Demtröder, P. Sebald, W. Meyer, Phys. Rev. Lett. 81 (1998) 4584. [31] D. Babikov, B.K. Kendrick, P. Zhang, K. Morokuma, J. Chem. Phys. 122 (2005) 044315. [32] B.K. Kendrick, J. Chem. Phys. 112 (2000) 5679. [33] B.K. Kendrick, J. Phys. Chem. A 107 (2003) 6739. [34] J.C. Juanes-Marcos, S.C. Althorpe, E. Wrede, Science 309 (2005) 1227. [35] J.C. Juanes-Marcos, S.C. Althorpe, J. Chem. Phys. 122 (2005) 204324. [36] B. Lepetit, R. Abrol, A. Kuppermann, Phys. Rev. A 76 (2007) 04072. [37] F. Bouakline, S.C. Althorpe, D. Peláez-Ruiz, J. Chem. Phys. 128 (2008) 124322. [38] F. Bouakline, B. Lepetit, S.C. Althorpe, A. Kuppermann, in: H. Köppel, D.R. Yarkony, W. Domcke (Eds.), The Jahn–Teller Effect, Fundamentals and Applications in Physics and Chemistry, Springer, Heidelberg, 2009, pp. 201– 237. [39] F. Bouakline, S.C. Althorpe, P. Larregaray, L. Bonnet, Mol. Phys. 108 (2010) 969. [40] J. Jankunas, M. Sneha, R.N. Zare, F. Bouakline, S.C. Althorpe, J. Chem. Phys. 139 (2013) 144316. [41] J. Jankunas, M. Sneha, R.N. Zare, F. Bouakline, S.C. Althorpe, Z. Phys. Chem. 227 (2013) 1281. [42] J. Schön, H. Köppel, Chem. Phys. Lett. 231 (1994) 55. [43] J. Schön, H. Köppel, J. Chem. Phys. 103 (1995) 9292. [44] S.C. Althorpe, J. Chem. Phys. 124 (2006) 084105. [45] S.C. Althorpe, T. Stecher, F. Bouakline, J. Chem. Phys. 129 (2008) 214117. [46] L.S. Schulman, Phys. Rev. 176 (1968) 1558. [47] L.S. Schulman, J. Math. Phys. 12 (1971) 304. [48] M.G.G. Laidlaw, C.M. De Witt, Phys. Rev. D 3 (1971) 1375. [49] G. Morandi, E. Menossi, Eur. J. Phys. 5 (1984) 49. [50] F.T. Smith, Phys. Rev. 179 (1969) 111. [51] C.A. Mead, D.G. Truhlar, J. Chem. Phys. 77 (1982) 6090. [52] H. Köppel, in: W. Domcke, D.R. Yarkony, H. Köppel (Eds.), Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, World Scientific, River Edge, NJ, 2004, pp. 175–204. and references therein. [53] B.K. Kendrick, C.A. Mead, D.G. Truhlar, Chem. Phys. 277 (2002) 31. [54] V. Vallet, Z. Lan, S. Mahapatra, A.L. Sobolewski, W. Domcke, J. Chem. Phys. 123 (2005) 144307. [55] Z. Lan, A. Dupays, V. Vallet, S. Mahapatra, W. Domcke, J. Photochem. Photobiol. A 190 (2007) 177. [56] M.D. Feit, J.A. Fleck, J. Chem. Phys. 78 (1983) 301. [57] D.K. Hoffman, T.L. Marchioro II, M. Arnold, Y. Huang, W. Zhu, D.J. Kouri, J. Math. Chem. 20 (1996) 117. and references therein. [58] P. Schwendner, F. Seyl, R. Schinke, Chem. Phys. 217 (1997) 233. [59] D.E. Manolopoulos, J. Chem. Phys. 117 (2002) 9552. [60] S.C. Althorpe, J. Chem. Phys. 114 (2001) 1601. [61] M. Abe, Y. Ohtsuki, Y. Fujimura, Z. Lan, W. Domcke, J. Chem. Phys. 124 (2006) 224316. [62] T. Rajagopala Rao, S. Mahapatra, J. Chem. Phys. 134 (2011) 204307. [63] M.N.R. Ashfold, A.L. Devine, R.N. Dixon, G.A. King, M.G.D. Nix, T.A.A. Oliver, Proc. Natl. Acad. Sci. USA 105 (2008) 12701. [64] M.G.D. Nix, A.L. Devine, R.N. Dixon, M.N.R. Ashfold, Chem. Phys. Lett. 463 (2008) 305.

Please cite this article in press as: F. Bouakline, Chem. Phys. (2014), http://dx.doi.org/10.1016/j.chemphys.2014.02.010