|
1. |
Mixed quantum classical steps: a DVR hopping method |
|
PhysChemComm,
Volume 3,
Issue 7,
2000,
Page 29-35
Adolfo Bastida,
Preview
|
|
摘要:
Mixed quantum classical steps: a DVR hopping method Adolfo Bastida,a José Zúñiga,a Alberto Requena,a Nadine Halberstadtb and J. Alberto Beswick*b a Departamento de Química Física, Universidad de Murcia, 30100 Murcia, Spain b IRSAMC, Université Paul Sabatier, 31062 Toulouse, France. E-mail: beswick@irsamc1.ups-tlse.fr Received 6th April 2000, Accepted 11th May 2000, Published 22nd May 2000 A new method that allows to mix a quantum and a classical description in a system where the quantum dynamics involves both bound and continuum states is presented. The description of the quantum degrees of freedom is based on a discrete variable representation (DVR), and the coupling with the classical degrees of freedom follows the surface hopping method using the fewest switches algorithm developed by Tully10 (J.C. Tully, J. Chem. Phys., 1990, 93, 1061). The method is applied to the simple model of a quantum particle (wave packet) colliding with a "classical" particle adsorbed on a solid surface, which has been studied by other methods. where the system can evolve on several quite different quantum states, a situation that frequently occurs in photodissociation for instance, this method does not describe properly the correlation between classical and quantum motions. In surface hopping methods, on the other hand, classical– quantum correlations are included by allowing a given trajectory to bifurcate properly into different branches, each governed by a particular quantum state and weighted by the amplitude of the state.The strengths and weaknesses of these two methods have been discussed in detail,4,9–12 including an exhaustive comparison between them by Tully.13 One difficulty with implementing surface hopping methods is that the quantum wave function must be expanded in an appropriate basis. The coefficients in this basis are examined at specified intervals of time in order to determine if a state-to-state transition should occur. This procedure becomes intractable if the basis includes a continuum of states, as is, for example, the case in the scattering of a wave packet from a non-rigid surface. This situation can also arise in photodissociation or collisional dissociation. Sholl and Tully14 have recently proposed a method to overcome some of these difficulties.Their method is an extension of the "molecular dynamics with quantum transitions" (MDQT) method of Tully10,15 to include a continuum of quantum states. The quantum wave function is treated as a truncated expansion of orthonormal basis functions together with a "lumped" state that represents all the continuum states. The authors have tested their generalized surface hopping (GSH) method on a simple model of a wave packet scattering from a non-rigid surface. They have shown that although considerably more computer demanding, the GSH method corrects for some of the qualitative failings of the mean-field dynamics, like the fact that some portions of the wave packet remain permanently trapped in resonances at low scattering energy. They also encountered a problem with classically forbidden hops.The problem comes from the fact that the continuum is represented by only one state, with energy equal to the average energy of the continuum states populated. When a 1 Introduction Accurate quantum mechanical simulations of the dynamics of molecular systems are usually not feasible for systems involving a large number of degrees of freedom. The multiconfiguration time-dependent Hartree (MCTDH)1 approach for instance, which has been applied successfully to the short-time dynamics of multidimensional quantum wave packets up to 24 nuclear degrees of freedom, becomes prohibitively computationally expensive for larger systems and long-time dynamics.Quantum approximate treatment, such as Miller’s initial value representation (IVR),2 has been receiving increasing interest lately3 and might become an alternative in the near future. On the other hand, conventional molecular dynamics can easily handle multidimensional problems, but they have severe limitations. The first one is due to the Born– Oppenheimer separation of electronic and atomic degrees of freedom, reducing the dynamics to atomic motion on a single adiabatic potential energy surface (PES). The second limitation stems from the use of classical mechanics to describe atomic motion. Mixed quantum–classical methods4–8 have been developed to address both of these problems. They are targeted at systems for which a few crucial degrees of freedom have to be treated quantum mechanically, but the rest, usually a very large number, can be treated classically.In particular, two approaches attempt to treat the interactions between the quantum and classical systems in a self-consistent way, in the sense that both the influence of the classical motion on the dynamics of the quantum degrees of freedom and the back reaction of the quantum system on the classical one are taken into account. The mean-field method, in which the classical degrees of freedom evolve on a PES averaged over the wave packet of the quantum degrees of freedom, has several interesting properties: the method is invariant to the choice of quantum representation (adiabatic or diabatic); it properly conserves total (quantum plus classical) energy; and it can provide accurate quantum transition probabilities when phase interference effects are unimportant and the energy difference of the relevant quantum states is negligible compared to the kinetic energy in the classical degrees of freedom.However, in cases DOI: 10.1039/b002769n PhysChemComm, 2000, 7transition from one of the bound states to the continuum should occur, it often happens that the energy of the continuum state is larger than the total semiclassical energy, and the hop is classically forbidden. To avoid unphysical results, the authors allowed these hops, thereby accepting that the total semiclassical energy be nonconserved. In this paper we present an alternative scheme, based on the surface hopping approach, that allows us to treat continuum states in quite a natural way.In this scheme, the quantum wave function is expanded in a discrete variable representation16 (DVR) rather than in the usual finite (spectral) basis representation (FBR). This allows treating of continuum states on a grid of points, as is usual in wave packet dynamics, using well-known absorption techniques to avoid spurious reflections at the limit of the grid. The interpretation of the dynamics of the quantum degree of freedom in terms of stochastic "hops" from one grid point to another, governed by the time evolution of the wave packet, is at the heart of the method, and is presented in section 2. The method itself, which we call "mixed quantum classical steps" (MQCS), is described in section 3.In section 4, we present its application to the same model system studied by Sholl and Tully,14 which represents the scattering of a wave packet from a non-rigid solid surface, and compare our results to the fully quantum mechanical solutions and to the generalized surface hopping method of Sholl and Tully.14 We conclude in section 5 with a discussion of our technique. (2) 2 Presentation of the method: quantum model Let us examine the propagation of a one-dimensional wave packet in a DVR basis set: (1) j x is a DVR basis function with where the wave packet Y(x, t) propagates with time along the x coordinate, and eigenvalue xj of the position operator. The values of the time-dependent coefficients can be obtained by solving the time-dependent Schrödinger equation using any conventional method (split operator, etc.).We now turn to an interpretation of the propagation of the wave packet that will allow the introduction of the mixed quantum-classical steps method. Initially, the wave packet is a coherent superposition of DVR states. The time evolution of the coefficients of the wave packet in these DVR states is given by where we have used the fact that the non-diagonal matrix elements of the Hamiltonian reduce to the matrix elements of the kinetic energy operator T in the DVR basis and the potential energy matrix is diagonal. Eqn. (2) determines the transition probabilities between DVR states. Let us introduce the density matrix elements (3) jj(t) represents the probability for the system to such that r be in state j at time t.The time evolution of rjj(t) is deduced from eqn. (2) (4) where bjk is given by (5) where stands for the imaginary part. In eqn. (4), the bjk coefficients can be interpreted as the population flux from state j to k, with the property that bkj(t) = –bjk(t). Hence if the system is in state j at time t, the probability to jump to state k at time t + dt is given by (6) We use this as follows. An initial state j is selected at random, following a probability distribution given by j jj t x x,0 r = ( )0 º ( )2. For each initial state we run a sequence of events by testing at each time step the probability for transition to any other state k.Suppose that the sequence of events has led to state j at time t. A random number g is generated between 0 and 1. If it is smaller than bj1, then at time t + dt the state will be 1; if bj1 < g < bj1 + bj2 then the state will be state 2, etc.; and if åk� jbjk < g, then the state at time t + dt will still be j. At the end of the propagation (t = tf), the recombination of the final states of all these chains of events gives the wave packet at time tf. This approach is similar to the way the quantum PES for the classical particle dynamics is selected in Tully's fewest switches method10 for mixed quantum–classical dynamics. However the origin of the coupling is different. In Tully's case the basis set is adiabatic and the coupling is the kinetic energy operator.In our case the basis set is diabatic and in the FBR representation the coupling would be provided by the potential energy operator. Because we are using the DVR representation the non-diagonal matrix elements of H in eqn. (5) are again due to the kinetic energy operator. This simulation can only provide a probabilistic presentation of the wave packet propagation, since it requires the knowledge of all the time-dependent coefficients cj(t), that is, of the accurate wavepacket Y(x, t), in order to evaluate the transition probabilities. However, it will be useful when introducing mixed quantum–classical dynamics, since the DVR representation allows us to treat dynamical problems involving a continuum in a very natural way.Let us now suppose that the one-dimensional wave packet defined in eqn. (1) is governed by a potential containing both bound and continuum states. The DVR grid should, in principle, extend to infinity and has to be truncated. Theexpansion of the wave packet on the DVR basis can be written as (7) where N is the number of DVR basis functions, i.e. the number of DVR grid points, and Y¥(x, t) represents the part of the wave packet leaking out in the asymptotic region (x > xN), where the potential has died out. This is the same as the separation of the wave packet in two parts, one in the interaction region on which propagation is performed and the other in the product region that can be stored for later analysis, first introduced by Heather and Metiu17 and extended to several dimensions by Pernot and Lester.18 We now follow the treatment of Sholl and Tully,14 but in the DVR representation.This gives a different meaning to the part of the wave packet which is not included in the linear combination of discrete states. In Sholl and Tully’s treatment, YC(x, t) stands for the continuum part of the wave function, which is then represented by a sole state, whereas Y¥(x, t) here represents the part of the wave packet that is no longer in the interaction region and should no longer influence the dynamics. The time dependence of the cj coefficients is obtained from the time-dependent Schrödinger equation for j £ N (8) From the definition of the density matrix coefficients [eqn.(3)], the time-dependence of the population in DVR state j, j îN, is given by (9) where (10) as before, and (11) In practice, it is rather difficult to evaluate bj¥(t) from eqn. (11). It is easier to use (12) with d rj j /dt computed by finite differences using eqn. (3). The probability for the system to jump from state j to state k is (13) and the probability to jump from state j to the asymptotic region is given by for j £ N (14) For completeness, the evolution of the asymptotic state population can also be specified. Since the total wave function is normalized, even if Y¥ is not orthogonal to the j x , we can define From this we derive (15) (16) where we have used eqn.(9). Since bk j = –bjk, the first term on the right hand side is zero and d r¥¥/dt can be written as (17) where we have defined b ¥j(t) = –bj ¥(t). Hence the probability to jump from the asymptotic region back to DVR state j is given by (18) The evaluation of P¥® j is very time consuming since it is necessary to calculate the b¥ j coefficients given by eqn. (12). In practice, Y¥(x, t) corresponds to the part of the wave packet escaping the grid, and the limit of the grid is set at a distance where the interaction potential is negligible. Hence neglecting the probability to jump from the asymptotic region back to the interaction region, P¥®j 0, is a very good approximation, as we could check in various examples.(19) 3 Presentation of the method: mixed quantum–classical steps In this section, we use the DVR propagation for the quantum coordinates in a problem where a number of degrees of freedom R a, a = 1, ... Nc, are treated classically and one or a few degrees of freedom (here we take one for the sake of simplicity) r are treated quantum mechanically. The total Hamiltonian for the system is written as where {R a} represents the Nc classical coordinates R a, T({R a}) the kinetic energy operator associated with these coordinates, and Hr(r; {R a(t)}) the rest of the Hamiltonian, including R-dependent potentials. The semicolon is a reminder that the {R a} are only parameters in the reduced Hamiltonian, and (t) denotes that they are a function oftime, hence Hr is a time-dependent Hamiltonian. In this paper, R a are cartesian coordinates.The quantum wave function evolves according to (20) where no explicit R dependence is introduced in the wave function. The wave function is then expanded in the DVR basis as in eqn. (7) (21) where j r , j = 1, ..., N, is a suitable DVR basis set and Y¥(r, t) represents the asymptotic parts of the wave packet that are leaking out of the grid. We then follow the procedure of the surface hopping method. Classical dynamics is described by an ensemble of trajectories, each trajectory being governed by the forces associated with a single basis function, with instantaneous switches ("hops") from one to another. The resulting classical equations of motion are (22) where focc is one of 1 r , ..., N r .If focc = Y¥, then the trajectory is stopped and the system is considered as dissociated. If focc is one of the DVR basis functions, the trajectories obeying eqn (22) conserve the total semiclassical energy (23) We now describe how the identity of focc changes with time. We follow the same procedure as in the preceding section. This is very close to the procedure of the generalized surface hopping method of Sholl and Tully,14 except that we are working in the DVR basis set. The probability for the system to be in the DVR state j at time t is given by (24) where bjk(t) is given by eqn. (10). The bj¥ coefficient appearing in eqn. (24) is calculated the same way as in eqn.(12). As before, the probability for the system to jump from state j to state k is (25) and the probability to jump from state j to the asymptotic region is for j £ N (26) We use this as in the preceding section, to evaluate the probability for the quantum degrees of freedom to change state at each time step [see the presentation following eqn. (6)]. At time t = 0, an initial state j is selected at random, following a probability distribution given by rj j(t = 0), and initial conditions for the classical coordinates can be sampled wh the usual techniques. Note that no special care is taken to conserve energy (or center of mass position) during the hops. It is not possible to follow Tully’s procedure10 based on Pechukas’ force,19 since we are not working in the adiabatic basis set.Neither is it possible to follow the procedure that we have used before20–25 for propagation in a diabatic basis set, where the change in momentum necessary to conserve energy is distributed among the coordinates as (27) for a transition occurring from state focc to state focc�. This is because the matrix element of Hr in eqn. (27) is only due to the kinetic energy operator associated with r, and for cartesian coordinates it does not depend on R a. However, in the same conditions, the one-dimensional wave packet propagation in the preceding section reproduced quite well the exact wave packet propagation. Note also that although the generalized surface hopping method does, in principle, conserve the semiclassical energy, Sholl and Tully14 had to abandon energy conservation for about 50% of the hops in their application to the simple model presented in the following section in order to obtain meaningful results.We have found similar behaviour in our calculations but a more thorough study of the amount of energy non-conservation is left for future work. 4 Application to a model case In order to demonstrate that the MQCS method is operational, and to compare its results to other existing methods, as well as to accurate results, we have applied it to the same simple model used by Sholl and Tully14 to test their generalized surface hopping model. It consists in a highly simplified representation of a light molecule colliding with a heavier molecule that is bound to a static, bulk surface. The model is restricted to one dimension, perpendicular to the surface.The coordinate r of the light particle with mass m is treated quantum mechanically, and the coordinate R of the heavy particle with mass M is treated classically. The Hamiltonian, defining this system, is written as (28)where Tr and TR are the kinetic energy operators. The interaction of the bound (heavy) particle with the surface is assumed to be harmonic the interaction of the light particle with the surface is described by a Morse potential (29) (30) and the light–heavy particle interaction is taken as exponentially repulsive (31) The heavy particle was assumed to be initially in the ground harmonic oscillator state and the light particle defined as a Gaussian wave packet, (32) The DVR points are equally spaced and identical to the ones used by Sholl and Tully.14 The coefficients of the wave packet on the DVR basis functions are determined by the transformation method of Aiani and Hutchinson.26 During the interaction of the light particle with the surface being taken as a Morse potential, the light particle can have bound states as well as continuum states.Due to its repulsive interaction with the heavy particle, the number and energy of these bound states vary with time. They represent resonances which can strongly affect the scattering of the light particle. Hence the discussion focuses on the time-dependent scattering probability of the model, which is defined as14 (33) where rs corresponds to the beginning of the product region of the grid.We used the same parameters in our simulations as Sholl and Tully.14 In Fig. 1 to 4, we present the time-dependent scattering probabilities obtained with the MQCS method, compared with the exact results, the mean field results, and the results of the generalized surface hopping method, the last three reproduced from Sholl and Tully, for four different values of the incident energy: E0 = 10, 20, 40 and 80 kJ mol–1. As noted by Sholl and Tully,14 both the mean field and the generalized surface hopping methods give the same behaviour at short time: the rapid increase in Ps occurs Fig. 1 Time-dependent scattering probability of the model of Sholl and Tully14 for E0 = 10 kJ mol–1, computed using exact quantum dynamics (dotted curve), mean field dynamics (dashed curve), generalized trajectory surface hopping (solid curve) and the MQCS method of this work (dotted–dashed curve).Fig. 2 Time-dependent scattering probability of the model of Sholl and Tully14 for E0 = 20 kJ mol–1, computed using exact quantum dynamics (dotted curve), mean field dynamics (dashed curve), generalized trajectory surface hopping (solid curve) and the MQCS method of this work (dotted–dashed curve). sooner than in the fully quantum simulations. This corresponds to the fact that the energy gained by the heavy particle during the collision is consistently less when using mixed quantum–classical dynamics: at low incident energies the heavy particle actually loses energy in the scattering process.This is a consequence of the classical treatment of the corresponding coordinate: it can use the energy initially given as zero-point energy, which is of course impossible in the fully quantum treatment. Since the three methods suffer from this same problem, the short time disagreements between the mixed quantum–classical and fully quantum calculations exist for all of them.Fig. 3 Time-dependent scattering probability of the model of Sholl and Tully14 for E0 = 40 kJ mol–1, computed using exact quantum dynamics (dotted curve), mean field dynamics (dashed curve), generalized trajectory surface hopping (solid curve) and the MQCS method of this work (dotted–dashed curve).Fig. 4 Time-dependent scattering probability of the model of Sholl and Tully14 for E0 = 80 kJ mol–1, computed using exact quantum dynamics (dotted curve), mean field dynamics (dashed curve), generalized trajectory surface hopping (solid curve) and the MQCS method of this work (dotted–dashed curve). Sholl and Tully14 have also discussed the asymptotic behaviour of the mean field dynamics, in which the scattering probability does not tend to one. In this method, the response of the classical degree of freedom to the quantum wave function is dominated by the majority quantum channels, in this case the continuum states. In this particular example, as the bulk of the wave packet is scattered, the kinetic energy of the classical particle is reduced almost to zero.Once the classical particle has no kinetic energy, any remaining bound portion of the wave packet is permanently trapped. The effect is strongest at the lower incident energies where the bound states play an important role, and disappears at E0 = 80 eV, where population of the bound states is negligible. In contrast, both our method and the generalized surface hopping method behave well at all incident energies. For the generalized surface hopping method, this is due to the fact that surface hopping allows the dynamics observed, when the system is in the majority channel or in the minority channel, to be quite different. Despite the fact that our method does not impose conservation of the semiclassical energy at the hops, it behaves equally well as the generalized surface hopping method of Sholl and Tully14 in the test case discussed here.The main advantage of our method is that it allows to treat continua of the quantum coordinate in quite a natural way. The asymptotic state Y¥ is only used to ensure normalization of the whole wave packet, and in order to avoid spurious reflection at the end of the grid, as in quantum wave packet propagation. It should be noted that in order to obtain physically meaningful results for this model, Sholl and Tully had to abandon the conservation of the total semiclassical energy. Because the continuum states are treated in a mean field way in their approach, many hops to the continuum are energetically unfeasible and a significant fraction of the trajectories remained permanently trapped even though the quantum amplitudes of the corresponding bound states decayed to negligible values at long times.Hence they allowed hops that do not conserve the total semiclassical energy. Non–conservation of the semiclassical energy at the hops does not violate energy conservation however, because the semiclassic energy involves only the classical particle and the occupied quantum state. 5 Conclusion We have presented a method that allows to treat mixed quantum–classical dynamics of systems involving a continuum of quantum states. This method could be characterized as the implementation of the surface hopping method in the discrete variable representation, rather than in the usual finite basis representation.The main advantage of the method is that it treats continua of the quantum coordinate in a very natural way, as in the propagation of a wave packet on a grid. We have demonstrated one possible use of the mixed quantum–classical steps method on the same simple model of a quantum wave packet scattering from a surface as Sholl and Tully.14 This allowed us to compare the results of our method with the ones of the generalized surface hopping mehod, and the mean field method, in addition to the exact results. As the generalized surface hopping method, the MQCS method corrects for some of the qualitative failings of the mean field method.It also does not conserve the semiclassical energy at the hops. In order to ascertain the validity of this approach, more tests have to be performed. It can already be pointed out that from the computational point of view, the method is quite performable, since the only matrix elements that need to be calculated are the ones of the kinetic energy operator in the quantum coordinate, which can be calculated beforehand and tabulated for use in the propagation, the potential matrix elements being diagonal and simply evaluated at the DVR points. This is much easier than calculating the average of the Hamiltonian over the total wave packet (mean-field method), or than calculating the adiabatic states at each propagation step like in the generalized surface hopping method.6 Acknowledgements We would like to express our thanks to Christoph Meier for very useful discussions. We are particularly grateful to Sholl and Tully for providing us with the detailed results of their calculations, and for many enlightening comments. A collaboration grant from the Spanish and French governments, Acción Integrada Picasso no.98155, is gratefully acknowledged. This work was partially supported by the Spanish Ministerio de Educación y Cultura under project PB 97-1041 and by the Fundación Séneca del Centro de Coordinación de la Investigación de la region de Murcia under project 01733/CV/98. References 1 H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73.2 W. H. Miller, J. Chem. Phys., 1970, 53, 3578. 3 M. F. Herman and E. Kluk, Chem. Phys.,1984, 91, 27; E. J. Heller, J. Chem. Phys.,1991, 94, 2723; X. Sun and W. H. Miller, J. Chem. Phys.,1998, 108, 8870. 4 J. C. Tully, in Modern methods for multidimensional dynamics computations in chemistry, ed. D. L. Thompson, World Scientific, Singapore, 1998, p. 34. 5 J. C. Tully, Int. J. Quantum Chem., 1991, 25, 299. 6 M. F. Herman, Annu. Rev. Phys. Chem., 1994, 45, 83. 7 G. D. Billing, Int. Rev. Phys. Chem., 1994, 13, 309. 8 M. Ben-Nun and and R. D. Levine, Chem. Phys., 1995, 201, 163; T. J. Martinez, M. Ben-Nun and G. Ashkenazi, J. Chem. Phys., 1996, 104, 2847. 9 U. Müller and G. Stock, J. Chem. Phys., 1997, 107, 6230. 10 J. C. Tully, J. Chem. Phys., 1990, 93, 1061. 11 C. C. Martens and J.-Y. Fang, J. Chem. Phys., 1997, 106, 4918. 12 G. D. Billing, J. Chem. Phys., 1997, 107, 4286. 13 J. C. Tully, Faraday Discuss., 1998, 110, 1. 14 D. S. Sholl and J. C. Tully, J. Chem. Phys., 1998, 109, 7702. 15 S. Hammes-Schiffer and J. C. Tully, J. Chem. Phys., 1994, 101, 4657. 16 Z. Bacic and J. C. Light, J. Chem. Phys., 1986, 85, 4594. 17 R. Heather and H. Metiu, J. Chem. Phys., 1987, 86, 5009. 18 P. Pernot and W. A. Lester, Jr., Int. J. Quantum Chem., 1991, 40, 577. 19 P. Pechukas, Phys. Rev., 1969, 181, 174. 20 A. Bastida, J. Zúñiga, A. Requena, I. Sola, N. Halberstadt and J. A. Beswick, Chem. Phys. Lett., 1997, 280, 185. 21 S. Fernandez Alberti, N. Halberstadt, J. A. Beswick and J. Echave, J. Chem. Phys., 1998, 109, 2844. 22 A. Bastida, J. Zúñiga, A. Requena, N. Halberstadt and J. A. Beswick, J. Chem. Phys., 1998, 109, 6320. 23 A. Bastida, J. Zúñiga, A. Requena, N. Halberstadt and J. A. Beswick, Chem. Phys., 1999, 240, 229. 24 S. Fernandez Alberti, N. Halberstadt, J. A. Beswick, A. Bastida, J. Zúñiga and A. Requena, J. Chem. Phys., 1999, 111, 239. 25 A. Bastida, B. Miguel, J. Zúñiga, A. Requena , N. Halberstadt and K. C. Janda, J. Chem. Phys., 1999, 111, 4577. 26 K. E. Aiani and J. S. Hutchinson, Comput. Phys. Commun., 1992, 69, 46. PhysChemComm © The Royal Society of Chemistry 2000
ISSN:1460-2733
DOI:10.1039/b002769n
出版商:RSC
年代:2000
数据来源: RSC
|
|