Quantum Molecular Dynamics of the Topological Properties of the

ARTICLE pubs.acs.org/JPCA

Quantum Molecular Dynamics of the Topological Properties of the Electron Density: Charge Transfer in H3+ and LiF Rodrigo Chavez-Calvillo* and Jesus Hernandez-Trujillo* Departamento de Física y Química Teorica, Facultad de Química, UNAM, Mexico, D.F. 04510 Mexico ABSTRACT: This paper presents a method to analyze the time evolution of electron density descriptors defined by the quantum theory of atoms in molecules. The wave packet nuclear dynamics was followed solving the time-dependent Schr€odinger equation. The time evolution of the nuclear wave packets was combined with the electronic wave functions to follow the time dependence of the average values of topological electron density descriptors. The method was applied to the reactive collision of H+ + H2 under different initial conditions and the photodissociation of LiF for either diabatic or adiabatic processes, with emphasis on the information provided by the time evolution of the atomic charges. These examples illustrate how this approach allows for a detailed analysis of the electronic structure in the time domain.

’ INTRODUCTION The time evolution of chemical processes at the molecular level is a relevant field of study in current chemistry. Available sources of information such as femtosecond pumpprobe laser experiments allow for the direct observation of chemical reactions as they proceed in time. In addition, the detailed knowledge obtained with these methods allows the control of elementary reactions by selectively breaking and making of chemical bonds.14 The theoretical methods have also reached enough development for the study of the reaction dynamics in chemistry, biology, and material science. As a consequence, the interplay between theory and experiment allows for a detailed microscopic characterization of chemical processes and also to propose reactivity models. The properties of the electron density play a major role in the understanding of molecular structure and chemical reactivity. In this respect, the quantum theory of atoms in molecules5,6 (QTAIM) provides us with a number of tools for quantifying useful chemical notions such as charge transfer and chemical bonding. This formalism has been applied to a wide number of systems using experimental or theoretical electron densities.7 In particular, electron density descriptors obtained from stationary wave functions have been used to understand the mechanism of chemical transformation along the reaction coordinate.8,9 The information concerning the time evolution of the electron density is not as abundant in the literature. The study of the coupled electronicnuclear dynamics including the time dependence of the electron density has been carried out for simple model systems.10 There are also reports of molecular dynamics simulations, using the method of density matrix propagation along classical trajectories in combination with a QTAIM analysis, that account for the oscillation of the electron density for molecules fluctuating around the equilibrium r 2011 American Chemical Society

geometries.11 Closer to the approach followed in this work, Grønager et al.12 analyzed the photodissociaton of the NaI molecule with the electron transfer taking place in the process being followed graphically by the electron density contributions of two valence orbitals. The theoretical analysis of how the electron density and related properties unfold in real time and the wealth of information this analysis could provide with emphasis on the QTAIM is a field that has not been extensively exploited in the literature and is the subject of this contribution. The small molecules H3+ and LiF were selected as test cases to carry out this study. These are well-known molecules that allow for the possibility of studying transition states, tunneling, change of molecular structure, or interaction with electromagnetic radiation as a means to illustrate the use of electron density descriptors for the characterization of chemical transformations using the QTAIM in real time.

’ TIME EVOLUTION OF THE ELECTRON DENSITY The state of a molecular system can be expressed in terms of a wave function Ψ(X,x,t) that depends on the nuclear coordinates, collectively denoted by X, on the position r and spin ω electronic coordinates, x = {r,ω}, and on the time, t. The wave function can be represented as a weighted superposition of electronic states ΨðX, x, tÞ ¼

n

∑i χi ðX, tÞψi ðx; XÞ

ð1Þ

Special Issue: Richard F. W. Bader Festschrift Received: May 29, 2011 Revised: September 5, 2011 Published: September 26, 2011 13036

dx.doi.org/10.1021/jp205003y | J. Phys. Chem. A 2011, 115, 13036–13044

The Journal of Physical Chemistry A

ARTICLE

In this expression, the summation runs over the n coupled states relevant for the process under study, and χi(X,t) is a timedependent nuclear wave packet whose square gives the probability of the system being in the i-th electronic state ψi(x;X) at X; note also that the electronic wave functions depend only parametrically on the nuclear coordinates. The time evolution of the system can be described in terms of Ψ(X,x,t) by solving the time-dependent Shr€odinger equation. In addition, the matrix representation of the electronic Hamiltonian can be diagonalized to yield a set {ψai } of adiabatic states whose eigenvalues {Eai } provide a potential for the nuclear dynamics. Using this representation, the corresponding wavepackets, {χai }, form a basis for a nondiagonal nuclear Hamiltonian with kinetic coupling terms which are difficult to calculate; these nonadiabatic coupling terms are usually small whenever the separation of the energy levels is sufficiently large3 but become relevant in the case of near degeneracy. A way to bypass this difficulty is to carry out a unitary transformation, U, of the electronic states to the so-called diabatic basis, {ψdi }, in which the coupling terms only contain multiplicative potential energy contributions easier to evaluate. In this work, the time evolution of selected properties of the electron density and of the topological atoms defined in the QTAIM are discussed. For nonstationary wave functions, the time evolution of the average value of scalar fields such as the electron density F(r) and its laplacian at r ∈ R 3 are obtained from eq 1 Z n ij ð2Þ ÆFðrÞæðtÞ ¼ dX χi ðX, tÞχj ðX, tÞFX ðrÞ

∑ij

Z Æ∇2 FðrÞæðtÞ ¼

n

dX

∑ij χi ðX, tÞχjðX, tÞ∇2FijX ðrÞ

ð3Þ

where the transition density between the electronic states i and j at the nuclear configuration X, FijX(r), is given by Z ij FX ðrÞ ¼ N dτ0 ψi ðx; XÞψj ðx; XÞ ð4Þ R In this expression, N is the number of electrons, and dτ0 denotes the integration over all the electronic spin coordinates and over the spatial coordinates of all of the electrons but one. As an additional example, the atomic charge of an atom Ω in a molecule can be conveniently denoted as qiX(Ω) in the sense that its value depends on both the nuclear configuration and the electronic state. The evaluation of this property involves the integration of |ψi(x;X)|2 over the basin of the atom Ω in the following way (in atomic units) Z Z i qX ðΩÞ ¼ ZΩ N dr dτ0 jψi ðx; XÞj2 ð5Þ Ω

ZΩ being the atomic number. For nonstationary wave functions, the time evolution of the average value of an atomic charge is obtained using eq 1 and eq 5 Z Z ÆqðΩÞæðtÞ ¼ ZΩ dX dτ0 jΨj2

’ COMPUTATIONAL METHODS The wave packet propagation procedure used was the Chebychev method with Fourier mapping for kinetic energy operators.1416 In this approach, the space of nuclear configurations is discretized on a uniform grid, and the wave packet amplitudes are sampled over this grid through time; the dynamic sampling is achieved by advancing an initial condition for the wave packet by finite time steps with the evolution operator exp(iH^δt), which is presented as a Chebychev expansion in the Hamiltonian governing the system motions. Time evolution is thus modeled by successive mappings of the Hamiltonian operator over the wave packet amplitudes on the grid. For the one-dimensional problems analyzed in this work, the grid has the form {Xj = jΔ, j = 1, ..., M}, where Δ is the step size on the nuclear degree of freedom and M is the number of nodes. The discretized potentials for the nuclear motion were obtained at each point of the spatial grid from electronic structure calculations of the corresponding electronic wave functions and their energies. Whereas the potential energy contribution is directly computed on the spatial grid, the kinetic energy is more conveniently calculated in the following steps:14 (1) transformation of the amplitude of the discretized wave packets to momentum space by means of a discrete Fourier transform;17 (2) the application of the multiplicative kinetic energy operator on the corresponding grid in momentum space; and (3) application of an inverse Fourier transform to return to position space. This procedure has the advantages of high accuracy and invariance of the commutation relations of quantum mechanics.14 The eigenfunctions of the Hamiltonian operator, which are used to set the initial condition in the LiF photodissociation dynamics below, are obtained by direct diagonalization of the Hamiltonian matrix in the discrete position representation induced by the grid.18 A computer program in C language was generated to carry out the wave packet time evolution involving arbitrary potentials in one or two degrees of freedom with the possibility of including dynamic potential coupling. The reactive scattering of the H3+ cation and the photodissociation of LiF were chosen as test examples for which the necessary potentials were generated from electronic structure calculations using the Molpro19 program. The performance of the code was tested and compared against reported results on model systems such as the wave packet dynamics over a potential step, a harmonic well, and the Morse potential.14,18 Practical versions of eq 2, eq 3, and eq 6 can be obtained when the discretization procedure is used. In the case of one degree of freedom and a regular spatial grid {Xj} with M nodes and step size Δ, the time evolution of the atomic charge is expressed as follows

Ω

Z ¼ ZΩ

R density, qijX(Ω) = ΩdrFijX(r), in the manner discussed in ref 13. In this notation, qiiX qiX. Notice that the time dependence of the average values given by eq 2, eq 3, and eq 6 is contained in the wave packets. The most important contributions to these values will come from the stationary properties at the nuclear configurations for which the wave packets have large amplitudes. In addition, these formulas can be used for either adiabatic or diabatic representations of the stationary electronic wave functions.

n

dX

∑ij

ij χi ðX, tÞχj ðX, tÞqX ðΩÞ

hqðΩÞiðtÞ ¼

ð6Þ

Note that this equation involves not only the atomic charges defined in eq 5 but also the basin integrations of the transition

M

n

∑k Δ ∑ij χi ðXk, tÞχjðXk, tÞqijX ðΩÞ k

ð7Þ

The graphical representation of the molecular graphs was carried out with Opendx,20 and Aimstudio21 was used for the calculation 13037

dx.doi.org/10.1021/jp205003y |J. Phys. Chem. A 2011, 115, 13036–13044

The Journal of Physical Chemistry A Scheme 1. Atomic Labels and Jacobi Coordinates for the H3+ Cation

Figure 1. (a) Ground state potential energy curve and (b) charges of the atoms in H3+. The horizontal solid lines in (a) denote the total energy of the wave packets for the reactive dynamics D1 and D2 that use two different initial conditions. The atomic labels are defined in Scheme 1.

of the stationary integrated properties. The atomic integration of the transition densities for the state averaged natural orbitals of LiF was carried out with a modified version of the AIMPAC22 program.

’ RESULTS AND DISCUSSION Reactive Scattering of H3+. H3+ is the simplest polyatomic

molecule. It has been the subject of numerous experimental and theoretical studies because of its interestelar abundance, its formation in hydrogen plasmas, and its use as a model for the neutral H3 molecule in high energetic states.23 The dissociation of this molecule can yield either H2 + H+, H2+ + H, or H+ + 2H. The global potential energy surface and analytical representations encompassing all of these processes have been reported because of the interest in describing the molecular dynamics and the rovibrational spectra of this highly anharmonic molecule.2426 One of the difficulties on accounting for its long-range behavior arises because of the presence of an avoided crossing of electronic states involving the two first dissociation channels when the internuclear distance between two of the protons is ca. 1.3 Å and the third proton is separated by at least 2.6 Å.27 From the point of view of the electron density analysis, this system provides the possibility of following in real time some interesting features of chemical reactivity such as the breaking and making of chemical bonds. In addition, the reactive collision

ARTICLE

of H+ + H2 is one of the simplest examples of charge transfer. The nuclear dynamics was studied restricted to the C2v channel which involves two equivalent D3h equilibrium structures connected by the linear arrangement of the transition state and the entrance/exit channels of H+ + H2 along the adiabatic ground electronic state. This portion of the potential energy surface does not involve the region of the avoided crossing of states which is relevant in the case of the linear approach. The electronic structure calculations were carried out at the CISD level (which for this two-electron system is equivalent to full CI) using the aug-cc-pVTZ28 basis set. It is convenient to represent the H3+ molecule in terms of the (r,R,θ) Jacobi coordinates, i.e., the HbHc separation, the distance from Ha to the HbHc center of mass, and the angle formed by the corresponding lines, respectively, as shown in Scheme 1. Hence, for the C2v channel, the intrinsic reaction coordinate shown in Figure 1(a) is defined in terms of the R coordinate with θ = 90 and allowing r to freely adjust for each R value. This means that even though the process involves two Jacobi coordinates the nuclear dynamics is described by a one-dimensional nuclear potential. At the equilibrium equilateral geometry, the Jacobi coordinates are Req = (0.758 and req = 0.875 Å, the latter being the HH internuclear distance. These values agree with other theoretical reports (see, for example, ref 25). In the case of the linear transition state, the corresponding values are RTS = 0 and rTS = 1.632 Å. In addition, the dissociation energy and the barrier to surpass the transition state are 4.591 and 1.761 eV, respectively. The stationary state atomic charges in H3+ were obtained using eq 5 and are shown in Figure 1(b) as a function of R. They are the atomic charges obtained at fixed nuclear configurations along the potential energy surface. At very large absolute values of R, the charges of the Ha and Hb atoms are +1 e and 0, respectively, in agreement with a proton far removed from the H2 molecule. Considering the negative values of R as the entrance region for the formation of H3+, as the proton approaches the H2 molecule, its atomic charge q(Ha) decreases because it gains electron density to the point that at the equilibrium geometry the three atoms bear identical atomic charges of +0.38 e each. Going beyond the energy minimum, an abrupt charge transfer continues toward Ha as the molecule becomes linear, q(Ha) = 0.27 e when R = 0; Figure 1(b) also reveals some hydridic character for this atom in the vicinity of the transition state. Figure 2 shows the evolution of the molecular structure of H3+ by the molecular graphs defined in terms of the gradient vector field of the electron density5 for some relevant molecular configurations along the C2v reaction coordinate. There are two structural catastrophes on each branch of the potential. The first of them appears at ca. R = 0.779 Å for the transition from a structurally unstable configuration to a cyclic arrangement by means of a bifurcation mechanism. Interestingly, the resulting structure involves the appearance of a non-nuclear attractor in the centroid of the triangle of H atoms with a charge of ca. 0.14 e at the equilibrium geometry, as shown in Figure 2. Further increase of R very soon reaches the point of the second catastrophe, ca. R = 0.757 Å, where the non-nuclear attractor vanishes by a bifurcation mechanism again to yield an open angular structure until the linear transition state is finally reached. Hence, the non-nuclear attractor persists only in a rather narrow interval of R values of width 0.02 Å enclosing the equilibrium geometry. This is a consequence of the steep electron transfer 13038



ARTICLE

Figure 2. Evolution of the molecular structure of H3+ along the potential energy curve. The light gray and red spheres denote the H attractors and the bond critical points, respectively; the non-nuclear attractor is displayed as a blue sphere. The atomic charges of the H atoms are also displayed; note the presence of a non-nuclear atractor at the equilibrium configuration. The approximate positions (in Å) of the structural catastrophes are shown in (b) and (d).

from the H2 molecule to the approaching proton at those R values. The next step is to analyze the time evolution of electron density descriptors as dictated by the nuclear dynamics for the reactive dispersion of a Gaussian wave packet by the potential shown in Figure 1(a). The effect of the initial conditions on the development of the dynamics was evaluated for two cases, denoted by the horizontal lines D1 and D2 in Figure 1(a), for the entrance from negative R values. Initially the wave packets were defined as χ(R) = Ng exp[(R R0)2/2σ20], Ng being the normalization constant, with zero kinetic energy and width σ0 = 0.42 Å centered at the classical turning points R0 = 2.65 and 1.06 Å for D1 and D2, respectively. They were propagated during 248 fs with a time step of 1 fs. In each case, the spatial grid was defined for the R interval from 5.92 to +5.92 Å, with M = 300 grid points. Whereas D1 is a high-energy case near the energy continuum, D2 allows us to analyze the tunneling through the transition state in the passage from the negative to the positive branch of the potential. Figure 3 shows some instances of the time evolution of the wave packets during these nuclear dynamics. For example, when 11 fs has passed on the D1 dynamics, Figure 3(a) shows that the wave packet already exhibits the dispersion provoked by the presence of the minimum in the potential energy well, and for t = 15 fs there is a maximum near the transition state. At t = 21 fs, the maximum of the wave packet is already in the R > 0 portion of the potential. In addition, the longest time displayed, 148 fs, is an instance showing the large nuclear delocalization over the reaction coordinate that can take place. The nuclear dynamics D2 evolves in time in a different manner. Figure 3(b) shows that even though the wave packet remains mostly localized around

Figure 3. Amplitude of selected wave packets dispersed by the potential for the C2v reactive collision H+ + H2 during the nuclear dynamics D1 in (a) and D2 in (b). The arrows on the abscissa show the positions of the energy minima.

the equilibrium geometry for the three times selected tunneling past the transition state can also take place. As an example, for t = 23 fs the tunneling probability is ca. 18%, obtained by the integration of |χ|2 for R g 0.5 Å. Figure 4(a) shows the time evolution of ÆRæ, the average value of R for the wave packets. The high-energy nuclear dynamics shows that the system behaves as a damped oscillator that oscillates about ÆRæ = 0 during the time evolution of the system. This does not mean that the system reaches and stays near the transition state when enough time has passed but that the wave 13039



Figure 4. Time evolution of (a) average wave packet position and (b) average atomic charge of the atom Ha during the reactive colision H+ + H2 for the dynamics D1 and D2. The horizontal lines in (a) denote the positions of the energy minimum and the transition state and in (b) the atomic charge of Ha for the equilibrium geometry.

packet eventually becomes symmetrically delocalized over the entire R axis, as in the example shown for t = 148 fs in Figure 3(a). This assertion is supported by the standard deviation for the wave packet, σ = (ÆR2æ ÆRæ2)1/2, which has a value of ca. 2 Å for t > 200 fs, thus encompassing molecular geometries beyond the equilibrium configuration. In contrast, in the case of the dynamics D2, the time evolution of ÆRæ shows that the wave packet virtually remains trapped on the left-hand well of the potential describing an anharmonic behavior all of the time; accordingly, the width σ of the wave packet oscillates with time around a much smaller value of ≈0.48 Å compared to the high-energy dynamics D1, reflecting the fact that less population is allowed to overpass the transition state region and reach the rightmost potential well. It is also interesting to observe in Figure 4(a) that the wave packet centroid coincides with the energy minimum for the first time much earlier in the case of D2 because of the initial conditions imposed. The time evolution of the average value of the atomic charge of Ha in the H3+ cation is shown in Figure 4(b) for the two nuclear dynamics. These values were obtained using eq 7 for the one state 2 case (n = 1). For example, Æq(Ha)æ(t) = ∑M j Δ |χ(Rj, t)| qRj (Ha) is a summation of the atomic charges of Ha weighted by the probabilities of the molecule being at the corresponding nuclear configurations {Rj} at the time t, multiplied by the step size. The average value of q(Ha) during the reactive collision of H+ + H2 in the high-energy case D1 shows an oscillating behavior about ca. +0.65 e. This value is bigger than the stationary atomic charge for the equilibrium geometry because the wave packet accumulates at longer distances, beyond the energy minima, in agreement with the σ value discussed above. As a consequence, the Ha atom retains a highly cationic character most of the time in the highenergy case. It is therefore possible to suppose that at long times the molecular graph that prevails is that of Figure 2(a) in the case D1. The situation is different when the system evolves in time on a state bound by the well of the energy minimum, as in the

ARTICLE

dynamics D2. In this case, Æq(Ha)æ oscillates around ca. +0.2 e, i.e., between the stationary values for the equilibrium geometry and the transition state, in the same manner as ÆRæ mostly oscillates between 0 and Re, Figure 4(a), because most of the time the system oscillates in the region of the energy minimum. In this case, the amplitude of the oscillations of the average values of both R and q(Ha) and the narrow stability interval for the molecular graph of the equilibrium geometry show the dynamical nature of the chemical bonding of the H3+ cation: the atomic interaction lines are formed and broken all the time, and in average, the bent molecular arrangement of Figure 2(e) dominates the dynamics. Photodissociation of LiF. The femtosecond photodissociation dynamics of alkali halides has been the subject of numerous experimental and theoretical studies.29 There are many aspects of these reactions that make them important for the understanding of the dynamical nature of chemical reactivity. In this respect, the NaI molecule has played a prominent role in the efforts to characterize the nature of excited state dissociative channels and the potential energy surfaces involved. For example, the adiabatic trapping of the molecule in the first excited electronic state, the electronic energy transfer on the diabatic dissociation, and the selective experimental control of the dynamics are important aspects that have been reported. In this work, the LiF molecule was selected as a model compound feasible of an accurate theoretical treatment at the quantum mechanical level with emphasis on the time evolution of electron density descriptors during the photodissociation dynamics of the molecule. A remarkable feature of the potential energy curves of LiF is the presence of an avoided crossing of the two lowest-lying electronic states of 1Σ+ symmetry that leads to an important coupling that has to be taken into account. To achieve a proper description of the system in the region of the curve crossing, the state-averaged multireference CI method including single and double excitations (MRCISD)30 was used with the Sadlej PVTZ orbital basis set31 giving equal weights to both states. The associated CASSCF calculations involve an active space defined by six electrons and the highest three occupied and the lowest three unnocupied natural orbitals. With this procedure, the adiabatic electronic wave functions ψai and energies Eai were obtained for the ground and first excited (i = 1,2) electronic states. The nuclear dynamics is more conveniently solved, however, by switching to the diabatic representation for the molecular Hamiltonian.4,12 The unitary transformation 0 1 cos αðRÞ sin αðRÞ A U ¼@ ð8Þ sin αðRÞ cos αðRÞ that defines the change from the adiabiatic to the diabatic basis ψd1 ¼ ψa1 U11 þ ψa2 U21 , ψd2 ¼ ψa1 U12 þ ψa2 U22

ð9Þ

is such that the corresponding dipole matrices are related by μa = UμdU†, and μd is diagonal.32 The parameter α(R) is obtained by diagonalization of μa for each internuclear separation, R. The coupling term is given by Ed12 ¼ U11 Ea1 U12 þ U21 Ea2 U22 Ea1

Ea2

ð10Þ ψa1

ψa2,

where and are the eigenvalues of and respectively, obtained from the multireference quantum chemical calculations. The diabatic potential energy surfaces obtained with this procedure are displayed in Figure 5(a) for both states. The ground state 13040



ARTICLE

dissociation energy and equilibrium separation are 5.76 eV and 1.60 Å, respectively, in agreement with the experimental values,33,34 5.96 eV and 1.56 Å, and with other theoretical results.32,35 At the equilibrium internuclear separation, Req, the ground state equilibrium geometry is dominated by the Li+F dipolar species and the first excited state by the so-called covalent Li + F configuration. At the crossing point of the diabatic states, Rc = 6.85 Å, the value of the parameter α is 41.088, denoting a nearly equitative mixture of the adiabatic states. It can also be seen in the inset of Figure 5(a) that the adiabatic curves virtually overlap the diabatic ones at short distances, but when R is close to Rc the adiabatic curves do not cross; this figure also shows that for R > Rc the ground electronic state corresponds to the covalent configuration; dissociation in this case involves an electron transfer, Li+ + F f Li + F. The properties of the adiabatic ground state electron density have been studied in detail before for the formation of LiF going beyond the equilibrium geometry to extremely short internuclear separations.8,36 The analysis in this work includes the first excited electronic state, the use of diabatic wave functions, and the time evolution of the system during the photodissociation. Figure 5(b) shows the behavior of the stationary state charge of the Li atom along the potential energy surfaces of both diabatic electronic states. These values were obtained replacing the wave functions ψd1 and ψd2 from eq 9 in eq 5 for each electronic state. In the case of the Li atom, the charges for the diabatic states are expressed in terms of the corresponding adiabatic charges, the basin integrations of the transition density, and the transformation matrix elements defined in eq 8, as follows qid ðLiÞ ¼ ZLi

2

∑jk Uij†q jk, aðLiÞUki ,

i ¼ 1, 2

ð11Þ

They agree with the dominance of the ionic configuration for the ground state and the covalent for the excited state at short distances with a reversed role beyond Rc. The corresponding adiabatic charges are also represented by the dashed lines shown in Figure 5(b). In this case, beyond the avoided crossing of states, there is an abrupt change in the atomic charge because the covalent configuration becomes the ground state and the ionic changes to the excited state in a small interval around Rc. In principle, the experimental conditions can favor either the diabatic or adiabatic dissociation mechanisms. In this work, the photodissociation process favors the diabatic channel because of the small difference between the electron affinity of F and the ionization potential of Li (which is responsible for the positive values of the ionic curve for R > Rc). This behavior should be contrasted with that of NaI where the excited state at long distances has a steeper energy increase of the positive values and the system can oscillate in a trapping well.37,38 For the LiF lightinduced dissociation, the full Hamiltonian governing the evolution of the wave packets is written ! Ed11 þ T Ed12 þ μ12 εðtÞ ð12Þ H ¼ Ed22 þ T Ed12 þ μ12 εðtÞ using the dipole approximation for the lightmatter interaction and T being the kinetic energy operator. In this expression, Ed12 is given by eq 10, and Ed11 and Ed22 are obtained using similar equations. For this work, we use a Gaussian electromagnetic pulse of the form ε(i) = exp[t2/2τ2]cos(ω12 t) of width τ = 10 fs tuned to the Bohr frequency ω12 between the electronic states

Figure 5. Diabatic (a) potential energy curves and (b) Li atomic charges for the two lowest 1Σ+ electronic states of LiF. The inset in (a) shows both the diabatic and adiabatic curves near the crossing of states. The diabatic Li atomic charges were obtained in terms of the adiabatic ones (shown as dashed lines) using eq 11.

at the equilibrium geometry. The field amplitude is chosen such that |μ12E0| = 0.06 hartree at this same geometry. The initial condition for the molecular state, expressed in terms of the nuclear wave packets, was taken as ! ! χ1 j1 ðRÞ ¼ ð13Þ χ2 0 where j1(R) is the ground vibrational state supported by the lowest diabatic potential energy curve, Ed11(R). This eigenfunction was obtained by the Fourier Grid Hamiltonian method.18 The vertical transition to the repulsive excited state is modeled by explicit propagation of this condition with the Hamiltonian defined by eq 12, and time-dependent perturbation theory is bypassed. The maximum excitation probability observed within the pulse duration was 0.57, measured by the integration of |χ2|2 over the R interval. The one-dimensional diabatic potentials of Figure 5(a) were used with M = 300 grid points, for R from 0.53 to 11.11 Å. Wavepacket propagation was performed during 148 fs with a time step of 1 fs. In addition, it is also illustrative to compare the time evolution of the photodissociation when an adiabatic process is favored. For this reason, the nuclear dynamics was also analyzed modeling the process with the adiabatic potentials and the same initial conditions without involving any coupling terms. This is a manner to mimic the behavior of another system with a larger transition probability to the ionic channel. Figure 6(a) shows several wave packets dispersed by the diabatic potentials during the photodissociation of LiF at selected times. For t = 6 fs, very early after the electronic transition, the wave packets of both the ground and excited states are localized and nearly superimposed on each other. At longer times, the wave packet of the covalent state has important amplitudes only in a narrow interval along the R window. This behavior has to be 13041



ARTICLE

Figure 6. Amplitude of the wavepackets dispersed by the ground and excited state (a) diabatic and (b) adiabatic potentials for the photodissociation of LiF. The dashed and solid lines correspond to the ionic and covalent electronic states, respectively.

contrasted with that for H3+ in Figure 3 for which the wave packets become highly dispersed. It can also be seen that even though the wave packet for the ground state remains centered near Req there is a small but not negliglible probability for the molecule in this configuration to dissociate because of the coupling of the electronic states. Note also that the excited state wave packet reaches the region of the crossing of states at ca. 94 fs; this time scale is similar to that found experimentally for the photodissociation of other alkali halides. Figure 6(b) shows the corresponding wave packets for the nuclear dynamics governed by the adiabatic potentials. In this case, they are still highly localized, and the excited state evolves in a similar manner as before; however, the ground state molecule remains trapped around the equilibrium geometry. The time evolution of ÆRæ of LiF is shown in Figure 7(a). The increasing behavior for the diabatic excited values correctly describes the dissociation processes and is accompanied by a moderate increase for the ionic configuration which has a local maximum near the time when the curve crossing is reached. This behavior has to be contrasted with that of the adiabatic process: as the excited state dissociates, the ground state oscillates about the equilibrium geometry. Figure 7(b) shows Æq(Li)æ for the photodissociation processes. These values are obtained replacing the stationary atomic charges (diabatic or adiabatic) in eq 7 for the two-state case (n = 2) for each internuclear separation on the spatial grid. The time evolution for the diabatic case (solid line) is discussed first. At the beginning of the process, Æq(Li)æ = +0.93 e, virtually the same as for the ground state equilibrium geometry, +0.94 e. This is a consequence of the peaked form of the wave packet which at t = 0 fs corresponds to the ionic ground state at Req. During the first 10 fs, the time when the light pulse promotes most of the population transfer to the excited state, there is an abrupt decrease of the average charge to +0.48 e, i.e., between those values for the ionic and covalent configurations (see Figure 5(b)). For some time, this average value remains nearly constant to later increase when the system reaches the avoided crossing of states (for t between 90 and 100 fs) because of a population transfer back to the ionic channel accounted for by the potential coupling terms. When

Figure 7. (a) Time evolution of (a) the average wave packet position for the ground and excited state diabatic and adiabatic potentials and (b) average atomic charges for the photodissociation dynamics of LiF. The solid and dashed lines denote the diabatic and adiabatic processes, respectively. The inset in (a) shows the oscillations about the equilibrium position of the ground state wave packet trapped on the adiabatic potential.

the system goes beyond Rc the average atomic charges become stable again. A remarkable feature of this photodissociation process is the fact that the biggest wave packet amplitudes remain highly centered around two characteristic nuclear configurations, namely Req and ÆRæ, the latter value increasing at constant rate with time for the covalent state; the corresponding packet χ2 broadens only slightly up to a width of ≈0.75 Å in a fashion that has been previously reported.2 This behavior has a direct connection with the observed evolution of the average value of the atomic charge (see eq 7), which now can be understood as being dictated mainly by the contributions from the stationary 13042


The Journal of Physical Chemistry A charges at the relevant nuclear configurations. This should be contrasted with the dynamics for the reactive collision of H3+ reported above, where the spread of the wave packets prevents identifying a dominant contribution to the average atomic charge values from any specific nuclear arrangement. The adiabatic process displays a behavior similar to the diabatic one before the excited wave packet arrives at Rc, the only difference being the larger charge values, Æq(Li)æ ≈ +0.57 e, in this case. However, after the curve crossing, Æq(Li)æ increases to +0.96 e. At these long times both χa1 and χa2 have their largest amplitudes at regions of ionic character, the former about the equilibrium geometry and the latter past the crossing point, showing that the adiabatic process is clearly one of collisional ionization.

’ CONCLUSIONS Molecular reaction dynamics is very much concerned with the description of elementary reaction mechanisms and the making and breaking of chemical bonds in real time. QTAIM on the other hand is a rigorous physical theory that provides the tools for the study of the molecular electronic structure and its change. The combination of the two can be bridged by the method of wave packet dynamics in the manner proposed in this paper. The construction of dynamic reactivity models based on the time evolution of the electron density descriptors can be an important complement to the richness of information that available experimental methods such as femtochemistry can provide. The examples analyzed in this work illustrate some possible applications of this approach. For example, the dynamical nature of the chemical bond in the H3+ cation, accounted for by the molecular graphs, can involve different molecular structures if two processes take place on the same potential energy surface but under different conditions. In the case of the photodissociation of LiF, the temporal evolution of the atomic charges is consistent with the diabatic or adiabatic nature of the process in question. These and other dynamical features of chemical reactivity can be modeled quantitatively in the time domain in terms of the topological properties of the electron density. ’ AUTHOR INFORMATION Corresponding Author

*E-mail: crystalinfl[email protected]; [email protected].

’ ACKNOWLEDGMENT The authors thank DGAPA-UNAM for financial support (project PAPIIT IN105609). The Supercomputing Department of DGTIC-UNAM is acknowledged for providing computational resources. R. Ch-C thanks CONACYT-Mexico for financial support (scholarship 235230). The authors also thank Prof. Richard F. W. Bader for allowing us to modify his AIMPAC program to calculate the atomic integrals of the transition density. ’ REFERENCES (1) (a) Carley, R.; Heesel, E.; Fielding, H. Chem. Soc. Rev. 2005, 34, 949–969. (b) Henriksen, N. E. Chem. Soc. Rev. 2002, 31, 37–42. (c) Kleyn, A. W. Chem. Soc. Rev. 2003, 32, 87–95. (2) Zewail, A. H. J. Phys. Chem. A 2000, 104, 5660–5694. (3) Levine, R. D. MolecularReaction Dynamics; Cambridge University Press: Cambridge U. K., 2009.

ARTICLE

(4) Tannor, D. J. Introduction to Quantum Mechanics. A time Dependent Perspective; University Science Books: Sausalito, CA, 2007. (5) Bader, R. F. W. Atoms in Molecules. A Quantum Theory; International Series of Monographs; Oxford University Press: New York, 1994. (6) Popelier, P. A., Ed. Atoms in Molecules: An Introduction; Prentice Hall: Harlow, England, 2000. (7) Matta, C. F., Boyd, R. J., Eds. The Quantum Theory of Atoms in Molecules. From Solid State to DNA and Drug Design; Wiley-VCH: Weinheim, 2007. (8) Hernandez-Trujillo, J.; Bader, R. F. W. J. Phys. Chem. A 2000, 104, 1779–1794. (9) (a) Gatti, C.; Barzaghi, M.; Bonati, L.; Pitea, D. In Quantum Chemistry: Basic Aspects, Actual Trends; Carbo, R., Ed.; Studies in Physical and Theoretical Chemistry; Elsevier: Amsterdam, 1989; Vol. 62, pp 401427. (b) Reid, D. L.; Hernandez-Trujillo, J.; Warkentin, J. J. Phys. Chem. A 2000, 104, 3398–3405. (c) Salinas-Olvera, J. P.; Gomez, R. M.; Cortes-Guzman, F. J. Phys. Chem. A 2008, 112, 2906–2912. (d) García-Revilla, M.; Hernandez-Trujillo, J. Phys. Chem. Chem. Phys. 2009, 11, 8425–8432. (10) (a) Erdmann, M.; Marquetand, P.; Engel, V. J. Chem. Phys. 2003, 119, 672–679. (b) Erdmann, M.; Engel, V. J. Chem. Phys. 2004, 120, 158–164. (c) Jhala, C.; Lein, M. Phys. Rev. A 2010, 81, 063421. (11) (a) Chevreau, H. Chem. Phys. Lett. 2004, 400, 59–61. (b) Chevreau, H. J. Chem. Phys. 2005, 122, 244316. (c) Soler, P.; Berges, J.; Fuster, F.; Chevreau, H. Chem. Phys. Lett. 2005, 411, 117–120. (12) Grønager, M.; Henriksen, N. E. J. Chem. Phys. 1998, 109, 4335–4341. (13) Bader, R. F. W.; Bayles, D.; Heard, G. L. J. Chem. Phys. 2000, 112, 10095–10105. (14) Kosloff, D.; Kosloff, R. J. Comput. Phys. 1983, 52, 35–53. (15) Kosloff, R. J. Phys. Chem. 1988, 92, 2087–2100. (16) Tal-Ezer, H.; Kosloff, R. J. Chem. Phys. 1984, 81, 3967–3971. (17) Cooley, J. W.; Tukey, J. W. Math. Comput. 1965, 19, 297–301. (18) Marston, C. C.; Balintkurti, G. G. J. Chem. Phys. 1989, 91, 3571–3576. (19) Werner, H.-J. et al. MOLPRO, version 2010.1; 2010; http:// www.molpro.net. Accessed on May 25, 2011. (20) The Open Source Software Project Based on IBM Visualization Data Explorer, http://www.opendx.org. Accessed on May 25, 2011. (21) Keith, T. A. AIMALL, version 10.07.01; 2010; http://aim. tkgristmill.com. Accessed on May 25, 2011. (22) Biegler-k€onig, F. W.; Bader, R. F. W.; Tang, T.-H. J. Comput. Chem. 1982, 3, 317–328. (23) (a) Burt, J. A.; Dunn, J. L.; McEwan, M. J.; Sutton, M.; Roche, A. E.; Schiff, H. I. J. Chem. Phys. 1970, 52, 6062–6075. (b) Oka, T. Phys. Rev. Lett. 1980, 45, 531–534. (c) Larsson, M. Int. J. Astrobiol. 2008, 7, 237–241. (d) McKenna, J.; Sayler, A. M.; Gaire, B.; Johnson, N. G.; Carnes, K. D.; Esry, B. D.; Ben-Itzhak, I. Phys. Rev. Lett. 2009, 103, 103004. (24) Carney, G. D.; Porter, R. N. J. Chem. Phys. 1974, 60, 4251–4264. (25) Meyer, W.; Botschwina, P.; Burton, P. J. Chem. Phys. 1986, 84, 891–900. (26) Rohse, R.; Kutzelnigg, W.; Jaquet, R.; Klopper, W. J. Chem. Phys. 1994, 101, 2231–2243. (27) (a) Bauschlicher, C. W.; O’Neil, S. V.; Preston, R. K.; Schaefer, H. F.; Bender, C. F. J. Chem. Phys. 1973, 59, 1286–1292. (b) Prosmiti, R.; Polyansky, O. L.; Tennyson, J. Chem. Phys. Lett. 1997, 273, 107–114. (c) Ichihara, A.; Yokoyama, K. J. Chem. Phys. 1995, 103, 2109–2112. (28) Thom, H.; Dunning, J. J. Chem. Phys. 1989, 90, 1007–1023. (29) (a) Rose, T. S.; Rosker, M. J.; Zewail, A. H. J. Chem. Phys. 1989, 91, 7415–7436. (b) Jouvet, C.; Martrenchard, S.; Solgadi, D.; DedonderLardeux, C.; Mons, M.; Gregoire, G.; Dimicoli, I.; Piuzzi, F.; Visticot, J. P.; Mestdagh, J. M.; D’Oliveira, P.; Meynadier, P.; Perdrix, M. J. Phys. Chem. A 1997, 101, 2555–2560. (c) Charron, E.; Suzor-Weiner, A. J. Chem. Phys. 1998, 108, 3922–3931. (d) Hartmann, M.; Pittner, J.; van Dam, H.; Bona^cic-Koutecky, V. Eur. Phys. J. D 1999, 9, 393–397. (e) Abu-samha, M.; Madsen, L. B. Phys. Rev. A 2010, 82, 043413. 13043



ARTICLE

(30) (a) Knowles, P. J.; Werner, H. J. Chem. Phys. Lett. 1988, 145, 514–522. (b) Knowles, P. J.; Werner, H. J. Theor. Chim. Acta 1992, 84, 95–103. (c) Werner, H. J.; Knowles, P. J. J. Chem. Phys. 1988, 89, 5803–5814. (31) Sadlej, A. J. Collect. Czech. Chem. Commun. 1988, 53, 1995–2016. (32) Werner, H.-J.; Meyer, W. J. Chem. Phys. 1981, 74, 5802–5807. (33) Brewer, L.; Brackett, E. Chem. Rev. 1961, 61, 425–432. (34) Pearson, E. F.; Gordy, W. Phys. Rev. 1969, 177, 52–58. (35) Bauschlicher, C. W.; Langhoff, S. R. J. Chem. Phys. 1988, 89, 4246–4254. (36) Hernandez-Trujillo, J.; Cortes-Guzman, F.; Fang, D. C.; Bader, R. F. W. Faraday Discuss. 2007, 135, 79–95. (37) Engel, V.; Metiu, H. J. Chem. Phys. 1989, 90, 6116–6128. (38) Rose, T. S.; Rosker, M. J.; Zewail, A. H. J. Chem. Phys. 1988, 88, 6672–6673.

13044


Quantum Molecular Dynamics of the Topological Properties of the

Recommend Documents