Article pubs.acs.org/JPCA
Quantum Nuclear Extension of Electron Nuclear Dynamics on Folded Effective-Potential Surfaces Benjamin Hall,†,‡ Erik Deumens,† Yngve Ö hrn,† and John R. Sabin*,†,§ †
Quantum Theory Project, Departments of Physics and Chemistry, University of Florida, Gainesville, Florida 32611, United States Tampa Preparatory School, Tampa, Florida 33606, United States § Department of Physics, Chemistry, and Pharmacy, University of Southern Denmark, Odense 5230, Denmark ‡
ABSTRACT: A perennial problem in quantum scattering calculations is accurate theoretical treatment of low energy collisions. We propose a method of extracting a folded, nonadiabatic, effective potential energy surface from electron nuclear dynamics (END) trajectories; we then perform nuclear wave packet dynamics on that surface and calculate differential cross sections for two-center, one (active) electron systems.
1. INTRODUCTION Quantum scattering calculations involving heavy projectiles have only rarely been completed using the fully coupled Schrödinger equation, in either its time-dependent (TDSE) or time-independent (TISE) form, as most approaches approximate some part of the problem. Common approximations include those based on the Born−Oppenheimer (BO) separation of slowly varying (nuclear) and quickly varying (electronic) degrees of freedom1−5 and, more recently, the approach of Deumens et al.6 to calculate the fully coupled dynamics of electrons with zero width, i.e., classical, nuclei, known as electron nuclear dynamics (END). It must be stressed at this point that all the approximation schemes, within their region of validity, must approach the same values. At sufficiently high projectile velocities (corresponding to energy between about 10 eV/amu and the ionization threshold), substantial agreement is found between all of the common methods and implementations. At low projectile velocity (from a few tens of eV/amu down to thermal energies), this broad agreement is lost. The methods diverge from each other, often even qualitatively. This disagreement is exacerbated by the absence of reliable experimental data, especially for differential electron-transfer cross sections. Again, this is due to a breakdown in the fundamental assumptions of the approximation: the BO separation can only be reliably truncated to low order for high collision energies.7 Similarly, the restriction to classical nuclei restricts the accuracy of END calculations at low projectile energies. A dominating approximation that has its roots in early work is the use of adiabatic or diabatic potential surfaces. The construction of potential energy surfaces was possible by fixing © 2014 American Chemical Society
the nuclei and only considering the electronic eigenfunctions at each discrete configuration. Each eigenstate corresponds to a different, independent potential energy surface. Formally, this is not an approximation if the complete set of eigenstates is used; in practice, only a truncated set can be used. Rearrangement collisions require a change of state and thus nonadiabatic behavior; this is introduced via off-diagonal Hamiltonian coupling terms that mediate transfers between surfaces. This is the starting point for BO based methods. The use of adiabatic potential energy surfaces requires a coordinate system in which the various types of motion can be efficiently decoupled, for example, internal coordinates (combined with angular momentum eigenstates) for each collision partner, with the interpartner distance being the sole continuous coordinate. The internuclear distance is the reaction coordinate for the scattering process and all other degrees of freedom are expanded in discrete states (e.g., partial waves for the rotational coordinates). The expansion in partial waves requires spherical central scattering potentials. Using the internuclear distance as the reaction coordinate leads to difficulties for electron transfer problems--without momentum in the wave function, the physical localization of the electrons at asymptotic separations cannot be maintained. Use of electron translation factors (ETFs) and switching functions can alleviate the difficulties, but at the cost of significantly complicating the coupled equations for the nuclear wave functions. Special Issue: Franco Gianturco Festschrift Received: January 16, 2014 Revised: March 20, 2014 Published: March 24, 2014 6385
dx.doi.org/10.1021/jp500532d | J. Phys. Chem. A 2014, 118, 6385−6394
The Journal of Physical Chemistry A
Article
have different values at each spatial point. This folding is analogous to the separate adiabatic surfaces, but the folded surfaces are independent except at the fold, instead of being separate solutions of the electronic Hamiltonian that are coupled everywhere by off-diagonal matrix elements of the molecular Hamiltonian. Once the potential energy surfaces are constructed, we use standard wave packet propagation tools to consider an evolving nuclear wave packet on the folded effective potential surfaces and the extraction of differential cross sections. This work is a first attempt at implementing these ideas in a usable form. We show examples of both the effective potential energy surfaces and the differential and total cross sections for a simple model: a two-center, one-effective-electron system. Thus, in this report, we propose a scheme for removing the limitation of END to classical nuclei and complete the methods and calculations described in ref 8. The proposal includes the nonclassical physics seen in even the simplest tunneling calculations and joins it to electronic behavior described by END. An alternate viewpoint for the quantum scattering problem is proposed, and an illustrative implementation of this framework is presented. The current work combines the intuitive conceptual ideas of the reaction coordinate description with the simplicity of the Cartesian calculations. The numerical implementation is suggestive of the accuracy possible with the present method.
The reaction coordinate description carries most of the intuitive simplicity of the classic one-dimensional tunneling models with a system-dependent coordinate system, which requires approximations to the underlying equations of motion and a priori knowledge of the final outcomes. This greatly restricts the range of applicability of the theory; extending it to new systems requires significant trial and error and the complexity hinders comparability between different calculations. In the mid-1990s, Deumens, Ö hrn, and co-workers6 proposed an explicitly time-dependent non-Born−Oppenheimer theory that includes full interaction among nuclei and electrons. The theory, electron nuclear dynamics (END), incorporates electronic momentum and nonadiabatic transitions in space-fixed Cartesian coordinates. No special systemdependent coordinates are needed, coupling elements are not computed separately (the equivalents of the BO based coupling elements are automatically generated from the equations of motion of END), and the theory extends without substantial modification to any nonrelativistic collision system (except ionization, which none of the common methods handle well).6 The original electron nuclear dynamics does have limited accuracy at low collision energies: To maintain the full coupling between nuclei and electrons, the nuclei are restricted to behaving as classical particles (formally, the zero-width limit of Gaussian wave packets). The trajectories are computed dynamically from the instantaneous forces exerted by the electrons and the other nuclei. This allows unrestricted fragmentation (without specification of the possible fragmentation channels) but limits the accuracy of the angular distributions of the scattered particles, as no interference is possible between adjacent classical trajectories. Here, we propose a scheme which makes accurate calculations at low energy possible. Although potential energy surfaces are generally taken to be adiabatic (or diabatic at most), it is possible to define an effective potential energy surface. This is the potential energy actually experienced at each point in space, including the instantaneous nonadiabatic effects; it is a linear combination of the possible electronic energies plus the nuclear repulsion. The fully coupled quantum-nuclear version of this effective potential can, in principle, be constructed; in practice, it is prohibitively difficult and amounts to solving the entire scattering problem at once. Instead, just as adiabatic potential surfaces are created by taking fixed nuclear geometries and solving the electronic eigenvalue problem at each geometry, it is possible to construct the effective potential energy surface from a set of dynamic (time-dependent) trajectories of classical nuclei responding to instantaneous forces. END provides these trajectories and, though the electronic energy is not an input for END, it is calculated along each trajectory. The classical trajectories then serve as sampling points on the effective potential energy surface. It must be stressed that this surface is not adiabatic. In fact, it cannot be easily decomposed into purely adiabatic terms--it naturally includes the excitations. Because the change in the END electronic state is smooth, the surfaces generated thereby will also be smooth. The extraction and characterization of these surfaces is the primary focus of the present work; we show that they intrinsically incorporate all of the information of the reaction-coordinate description by creating a complicated topology in space-fixed relative coordinates. We talk of these surfaces as having been “folded” so that at various times in the history of each trajectory the system energy may
2. METHODOLOGY 2.1. Electron Nuclear Dynamics (END). Electron nuclear dynamics (END) is an application of the time-dependent variational principle (TDVP), which additionally includes explicit coupling between electrons and nuclei at all times. No potential-energy surfaces are used and neither the Born− Oppenheimer nor the adiabatic approximation is enforced. The theory underlying END and its applications to physical systems has been explained in detail elsewhere.6,9,10 As a result, we present an only an overview of the characteristics of the method. Analysis of other details of the scheme, such as surface hopping, will be addressed in further publications. As described in the above references, END is an alternative formulation of the semiclassical scattering problem that explicitly couples a wave function description of the electrons to classical nuclear trajectories.6 The trajectories are dynamically determined on the basis of the internuclear forces and the instantaneous electronic configuration. The trajectories are not based on precomputed potential energy surfaces. Also, without the Born approximation or adiabatic surfaces, all instantaneously available electronic states are coupled without explicit Hamiltonian coupling elements.6 Without these coupling elements, the calculations can be done in Cartesian coordinates and electron-translation factors1 are included intrinsically without switching functions or other corrections. Electron nuclear dynamics uses classical nuclei. Formally, the nuclei are modeled as the zero-width limit of a Gaussian wave packet and thus all nuclear coordinates and conjugate momenta are replaced by their expectation values. As the electronic wave function depends parametrically on the nuclear coordinates and momenta (and the nuclear trajectories are coupled to the electronic state), relaxing the zero-width approximation, though performing well for energies greater than 10 eV/amu, becomes inadequate at lower energies. At low energy ( 0.5 hartree, the reflection (dotted) is unity, and the transmission (solid) is zero.
potentials near the fold must control the multisurface scattering as well. Thus one can think in terms of reaction coordinates but calculate in terms of Cartesian coordinates. Thus, it is not necessary to compute the reaction coordinate and rewrite the Schrödinger equation in those coordinates; one can work in familiar center-of-mass Cartesian coordinates. As is discussed in previous publications,6,9,10 the evolving wave function is at all times a single determinant, but it is not a single determinantal solution to any stationary-state equations such as the SCF equations for the Hartree−Fock problem at the given nuclear geometry. Rather the determinantal wave functions, through the time-dependent coefficients z(t), are propagated independently of, and in interaction with, the nuclear geometry evolution. The electronic wave function is therefore neither multiconfigurational nor multireference. To develop this correspondence further, we examine both the one-dimensional tunneling case and the one-dimensional folded potential scattering case for a head on collision, where the projectile return path is the same as the incoming trajectory, but with a different energy. The solution of the onedimensional scattering problem can be found in many introductory quantum mechanics texts and need not be presented here. We will sketch the major steps using the method proposed here, and discuss the solution for a specific, illustrative case. This example shows all the characteristics of the method proposed. Consider a simple, one-dimensional, two particle collision where the interparticle potential is ⎡0 for x ≤ 0 ⎢ for 0 < x ≤ a V (x) = ⎢ V0 ⎢ ⎣ 0 < V1 < V0 for x > a
Figure 2. Transmission and reflection from a square barrier. Transmission (solid) and reflection (dotted) off the square barrier shown in Figure 1 for a monoenergetic incident plane wave. The asymptotic potential difference is 0.5 hartree, causing a sharp threshold for transmission.
The transmission then rapidly increases to near unity. Thus the only data needed to create the transmission and reflection coefficients are the width and height (relative both to reactant and to product regions) of the potential surface. For a nonsquare barrier, this reduces to the maximum height of the barrier (above the entry channel) and the relative slopes of the potential near the peak. Physically, the gradient of the potential (in 2+ dimensions) defines the forces on the piece of the wave packet at that location. The effective scattering potential (including internuclear repulsion) in center-of-mass Cartesian coordinates can be mapped onto the reaction-coordinate scattering problem. In Cartesian coordinates the projectile has the opportunity to scatter back over its initial trajectory, forming multiple uncoupled surfaces. The surfaces may be degenerate, but the transitions between these surfaces only occur near a specific area in physical space, namely the fold between the surfaces. This fold is the point (in one dimension) or, in general, a surface of dimensionality one less than that of the system at which the scattering event can be described and is similar in concept to a Riemann sheet in comkplex analysis, where the motion continues smoothly on the surface. In Cartesian space and one dimension, the fold will always occur at the smallest internuclear distance attained; in two dimensions it tends to follow the boundary of the energetically forbidden region quite tightly. At this fold, ensuring that the wave function is continuous and smooth at all points guarantees the proper transfer of wave function density from one surface to another. At all other points, away from the fold, the wave function can evolve freely. Folds with high local gradients, and thus forces, shunt wave function density away from the fold, either preventing it from transferring or keeping it on the new surface. Folds with low local gradients allow the wave function to oscillate many times before leaving the area, resulting in significant transitions. The foregoing discussion is in space-fixed Cartesian coordinates, leading to the complex topologically, but computationally simple potential surfaces. A characteristic of space-fixed coordinates is the existence of forbidden regions,
(3)
There are thus three regions in space: an asymptotic region with zero interaction potential energy, another asymptotic region (on the other side) with potential energy V1, and a barrier of height V0 and width a, as illustrated in Figure 1.
Figure 1. One-dimensional scattering potentials: realistic exponential potential (dotted) and square barrier approximation (solid).
The same figure also shows a simple physical potential to which this square barrier can be thought of as an approximation. In terms of the projectile energy E, there are three regions of interest as well: 0 < E < V0 where reflection should dominate, V0 < E < V1 where we expect significant amounts of both reflection and transmission, and E > V1 where 6388
dx.doi.org/10.1021/jp500532d | J. Phys. Chem. A 2014, 118, 6385−6394
The Journal of Physical Chemistry A
Article
regions in space where the wave function is energetically forbidden from penetrating. Figure 3 shows a selection of arbitrary trajectories approaching a scattering center (solid circle).
Figure 5. Potential surface along a single trajectory in Figure 4. The potential cusp is located along the dashed curve. V0 represents asymptotic separation below which no transfer of density can occur.
Figure 3. Cartesian representation of the scattering process. Trajectories (solid lines) approach the boundary of the forbidden region (dashed curve) and continue (dotted lines), possibly on a different fold of the potential surface.
representation, the entire process can be described as multidimensional tunneling through and/or scattering over a barrier. For energies less than the asymptotic separation (E < V0 in Figure 5), no transmission (transfer to the “upper” surface) can occur. For E > V0, both transmission and reflection (elastic scattering in the Cartesian picture) occur, with angledependent ratios that depend on both the incoming energy and the relative gradients of the potential surfaces near the cusp, or forbidden region, just as in standard time-independent wave function tunneling. From a dynamic perspective, the incoming wave packet approaches the barrier and splits, part of it either continuing over the barrier or tunneling through it, the rest reflecting back toward ξ = −∞. The part of the density that traverses the barrier becomes the product and the rest remains reactant-like. This alternate description is conceptually simple and straightforwardly generalizes to N > 2 dimensions but is computationally difficult. Many common methods require the tedious transformation to the appropriate reaction coordinates, and the choice of coordinates alters the form of the coupled equations to be solved. Often ill-suited coordinates must be chosen for computational purposes, which then lead to problems such as spurious coupling elements at long-range or origin-dependence in the results. The Cartesian approach is complicated conceptually, but quite simple computationally, as no coordinate transformations must be performed, nor do the equations have to be modified. We therefore propose to think in the reaction coordinate but calculate in the Cartesian frame. This is possible because the two are related by a one-to-one transformation, implying that the results must be the same in either frame. One particularly simple “unfolded” coordinate system for two or more dimensions is motivated by a technique used to control the onset of time-dependence in perturbation calculations. In essence we transform from Cartesian coordinates (x,z) at one extreme to impact parameter and time coordinates (b,t) at the other. Here, Z(t) is a timedependent vector of expansion coefficients of the dynamic orbital in terms of the reference orbitals and completely characterizes the time evolving determinantal wave function. We then define a new coordinate system (x′(t), z′(t))by the transformation
The dashed parabolic curve shows the border of the forbidden region, and thus the transition point between surfaces. Trajectories that encounter this boundary are deflected backward, possibly at different energies (dotted lines). An alternate, but completely equivalent, description is from the point of view of the projectile. This is a transformation of the coordinate system into one where the same trajectories shown in Figure 3 become straight lines in the new (nonCartesian) coordinate system. This is shown in Figures 4 and 5. Figure 4 shows the same setup as Figure 3, transformed into the new coordinate system (i.e., from (xrel, zrel) to (ς, ξ). The potential along each trajectory is now described by forms similar to those of Figure 5, where the cusp in Figure 5 corresponds to the intersection of the trajectory (solid line) with the (transformed) border of the forbidden region. In this
Figure 4. Reaction coordinate representation of the scattering (tunneling) process. Trajectories (solid lines) pass through the barrier (peak represented by dashed curve) and are either reflected or transmitted. 6389
dx.doi.org/10.1021/jp500532d | J. Phys. Chem. A 2014, 118, 6385−6394
The Journal of Physical Chemistry A
Article
Figure 6. END-based dynamic potential surface for He2+ + H under a linear transformation from (b, v0t) coordinates (α = 0) to (x, y) coordinates (α = 1) showing selected trajectories and potential energy surfaces.
on trajectory backscatters while the others scatter at various angles. At α = 1 (the normal Cartesian coordinates, Figure 6d), the trajectories bend around the well to conserve angular momentum, leading to a forbidden region, the white area in the middle of the grid. 2.3. Wave Packet Evolution. We use the second-order split-operator method15 for time evolution, with a Fourier transform representation of the kinetic energy operator. In future implementations, more advanced iterators, such as the Lanczos16 method, can be included for increased accuracy. Absorbing imaginary potentials17 are placed at the edge of the main computational region to absorb spurious reflections and periodicities. Wave function density in the asymptotic region is transferred to a separate asymptotic grid and allowed to evolve freely until the total integrated population of the primary grid is less than a threshold, usually 1%. Because the individual sheets of the folded surface are not coupled except at the juncture, the wave packet can only move between surfaces through the fold lines (fold surfaces for higher dimensional problems). Population transfer between surfaces is enforced by ensuring continuity and first-order smoothness across the join. Because this causes duplication of the wave packet due to Fourier aliasing (see section 3.2 for discussion), absorbing potentials are placed inside the forbidden regions to eliminate the spurious extra packet.
x′(t ) = (1 − α)b + αx(t ) z′(t ) = (1 − α)v0t + αz(t )
(4)
where α is a continuous parameter with values between zero and unity that describes how unfolded the coordinates are and x represents all coordinates perpendicular to the incoming beam direction k̂ (t = 0), b are the impact parameters, and v0 = (2E/m)1/2 is the initial speed of the projectile’s classical trajectories. The significance of the parameter α is to continuously convert a single, unfolded surface with no forbidden regions (α = 0) into a strictly-Cartesian, multiplevalued surface (α = 1). At intermediate values of α, the surface has characteristics of both, allowing the visualization of exactly how the surfaces transform into each other. A particularly striking example of the value of this simple transformation is shown in Figure 6. Here we see the surface generated by colliding He2+ with ground-state atomic hydrogen at an energy of 2.5 eV/amu. The collision energy is such that there is no electron transfer (all state-to-state probabilities are smaller than 10−5). Shown are four values of α: α = 0.0, 0.1, 0.32, 1.0, chosen to illustrate the evolution of the unfolding process. Plotted in black on top of the surfaces are five representative trajectories with impact parameters b = 0, ± 1.0, ± 2.5. When α = 0 (Figure 6a), the head-on trajectory goes down into the attractive region, up over the repulsive core, and out the other side. The trajectories at higher impact parameters enter the well later, leading to the asymmetry in z′. At higher values of α (Figures 6b,c), the head6390
dx.doi.org/10.1021/jp500532d | J. Phys. Chem. A 2014, 118, 6385−6394
The Journal of Physical Chemistry A
Article
3. RESULTS AND DISCUSSION The wave packet evolution scheme described above as well as the extraction of the potentials from the END trajectory data (described in more detail in ref 8) has been implemented, and some test results are reported below. This is a reference implementation and can be much improved; the results are suggestive of the method rather than being perfectly precise. The remainder of this section describes some test cases. 3.1. Square Barrier. A simple test case for which analytic solutions are available is the one-dimensional asymmetric square barrier, where the potential energy is defined by eq 3, where V1 = 1.0 hartree, V0 = 0.5 hartree, and a = 2.0 au. The analytic form of transmission function T(E) for an incident monoenergetic point particle, represented by a plane wave is shown in Figure 2. Note the sharp onset of transmission when E = V0 and the decaying oscillations for E → ∞. The first is a consequence of conservation of energy: no propagating states exist for E < V0 where x > a. The maxima of the oscillations occur when [2m(E − V0)]1/2 = nπ, for integer n; the resulting perfect fit of the overbarrier segment of the wave function results in perfect transmission. A monoenergetic incident free particle, represented by a plane wave, is an approximation to a physical particle, especially for heavier particles such as atoms or molecules. To examine the effects of an incident particle with uncertain momentum, as well as to compare better with the results from results calculated using the scheme described here, we calculated the transmission coefficient as a function of the average energy of a Gaussian wave packet, T (E ̅ )
Figure 7. Transmission of Gaussian wave packet through asymmetric square barrier (eq 3): ENDwave folded surface (+ points), ENDwave single-surface tunneling (× points), and exact (solid line).
3.2. Aliasing. A difficulty inherent in any finite-space discrete Fourier transform is aliasing. Aliasing is the name given to the spurious duplication of any signal that occurs at wave numbers (for position-space transformations) higher than the Nyquist frequency18 (half the sampling rate) for the grid. The signal thus duplicated appears in the negative wavenumber portion of the transformed grid and interferes with any signal properly occurring there. Of course, signal from the negative wavenumber side can be aliased into the high positive wavenumber region as well. It must be emphasized that aliasing is an inherent characteristic of discrete Fourier transforms. Completely avoiding aliasing requires that the momentum-space signal (the wave function in this case) is band limited to less than the Nyquist frequency in both positive and negative directions. For most of the evolution of the wave packet, this requirement is met easily: the packet is sharply peaked around a mean value. Two points in the evolution pose problems, however. First is the transition between surfaces at the fold. Because the transfer between surfaces occurs over a narrow region of the grid, this reduces the local sampling rate below the physical limits of the transformed wave packet. As a result, a copy of the wave packet is created with the opposite direction of momentum from the physical packet. This must be absorbed by the optical potential in the forbidden region, as described above. A second, more serious concern is the transfer to the asymptotic grid for analysis. Again, the limitation to only a small section of the primary grid induces aliasing on the transferred packet. In one dimension, this merely shifts the peak of the wave packet toward zero momentum, but with a dense enough grid the problem is minimal. In two or more dimensions the aliasing causes more significant problems, as will be shown in subsection 3.3. Various techniques exist18 to limit the damage caused by this aliasing. In this implementation, only the “brute force” method (adding more grid points in the affected regions) has been implemented. It should be noted that there have been several attempts in this regard in scattering theory, for example the work by Kouri et al.19 3.3. Two-Dimensional Analytic Potential. We present the results for two-dimensional Yukawa potentials; the folded structure is modeled with two separate intersecting potential surfaces.
k
∫ dk A k1 exp(−α(k − k ̅)2 )Tpw(k)
(5)
where k = (2mE)1/2 is the momentum of the incident plane wave component and k 1 = [2m(E − V 0 )] 1/2 is the corresponding plane wave component in the exit channel. The effect of integrating over the wave packet energies is to incorporate plane waves with both higher and lower energies than the peak. This means that for all but the narrowest packets the onset of transmission is no longer sharp; even at nearly zero mean energy some of the packet has enough energy to tunnel through the barrier. Because there is no longer a single wavenumber, the perfect wavelength matching that resulted in the oscillatory peaks is washed out by the other components, resulting in a slow monotonic rise to complete transmission as the energy increases above the barrier. These effects can be seen in Figure 7. For comparison, Figure 7 also shows the results of two calculations using the proposed scheme. The first, shown as dashed lines and points, is a single-surface wave packet scattering calculation with the same barrier as above; the computational region is 20 au wide and consisting of 512 grid points. This serves as a check that the propagation, without folded surfaces, is working properly. The second, shown as disconnected points, is a folded-surface representation of the same problem. The potential curve has been hinged about the point x = 1 au and rotated back on itself to form a 10 au long folded surface with a transfer point at x = 1 au. Note that these points lie very close to the previously calculated ones, both theoretical and computed, showing that the folded surface dynamics accurately models the physical ones. The accuracy is within the expected numerical accuracy, given the number of points and time step chosen. 6391
dx.doi.org/10.1021/jp500532d | J. Phys. Chem. A 2014, 118, 6385−6394
The Journal of Physical Chemistry A Vi =
gi exp( −kiR )
Article
corresponds to the direct transfer, but the low projectile mass means that the wave packet has spread greatly by the time it reaches the asymptotic region. Interference and aliasing provide the circular banding. Note that only the first quadrant of the whole map is shown: the second quadrant (negative px and positive pz) is the mirror image of this quadrant, and all negative pz points can only be reached by aliasing. The first quadrant corresponds to the angular region 0 ≤ ϑ ≤ 2π. Figure 10 shows the elastic and transfer doubly differential cross sections at the mean packet energy of 408 eV with a
+ dVi
(6) R The parameters for the two potentials are g1 = 2.0 hartree, g2 = 1.0 hartree, k1 = k2 = 2.0 au−1, dV1 = 0 hartree, and dV2 = 0.1 hartree. This produces two purely repulsive potentials that rise as 1/R near R = 0; the two overlap at R = 1.0 au and have an asymptotic difference of 0.1 hartree. Figures 8 and 9 show color maps of the elastic and transfer, respectively, momentum distribution of the asymptotic out-
Figure 8. Color map of the elastic momentum distribution for transfer between two Yukawa potentials. Note aliasing (banding) due to insufficient momentum resolution. The forward peak is the portions that were just reflected from the barrier; the peak near zero momentum has been transferred back from the upper surface.
Figure 10. Elastic (solid, red) and transfer (dashed, green) doubly differential cross section in au2 for two Yukawa potentials17 as a function of angle at the mean packet energy E = 15 hartree (408 eV).
reduced mass of 1. As the barrier is only 1 hartree in height, the transferred cross section (dashed) is between five and many orders of magnitude higher than the elastic. It does not fall off strongly with angle; part of this is due to the aliasing problem discussed above, but part of it is that the potential has no forbidden region for R = (x2 + z2)1/2 > 1 where the potentials overlap. This leads to almost spherical transfer. The elastic cross section falls off rapidly with angle away from the central peak into the side lobes of the aliasing. The structures are predominantly aliasing-related. This system shows that improper or inadequate choice of grids leads to scattering that is qualitatively correct, but whose cross sections cannot be accurate due to aliasing effects. In general, the current implementation can be used to do the scattering, but a better implementation of the cross-section-determination routines is required before the cross sections can be fully trusted. 3.4. Test Case: H+ + H(2s). As an example of a relatively complicated folded surface, we consider a collision of a proton with a hydrogen atom in the first excited state.20 END trajectories were calculated from b = 0.0 au to b = 10.0 au using the triple-ζ augmented basis of Dunning21 the basis set reproduces the lowest few atomic electronic states adequately for this test. The initial state of the target hydrogen was 2s, leading to an initially spherical electronic wave function and obviating the need for multiple orientations. This would not have been possible with any of the p-type hydrogen functions. Color maps of the asymptotic momentum distributions (elastic and transfer respectively) are shown in Figures 11 and 12. The elastic scattering (Figure 11) shows the effects of aliasing the same as all the other systems examined but otherwise is confined to a small area of momentum space centered at the incident energy, as should be expected for
Figure 9. Color map of the transfer momentum distribution for transfer between two Yukawa potentials. Note the large spread in momentum due to the light mass (m = 1) and the circular banding patterns (aliasing) due to a narrow transfer window.
going wave function. As can be seen in Figure 8, two structures appear: one centered near pz = 5 au−1 and the other right near pz = 0. The forward peak is the reflected pulse from the potential barrier, only slightly scattered from its initial configuration. The fainter bands on either side are the effects of Fourier aliasing caused by the use of too narrow a window to accumulate the asymptotic wave function. It is presented here as an example of the pathologies that can occur when the initial parameters are set without care. Better use of a greater number of grid points and increasing the width of the “asymptotic” region would be expected to partially alleviate this effect, but nothing can completely remove it. Figure 9 shows the transferred portion of the same wave function. The peak 6392
dx.doi.org/10.1021/jp500532d | J. Phys. Chem. A 2014, 118, 6385−6394
The Journal of Physical Chemistry A
Article
Figure 11. Color map of the momentum distribution of the elastically scattered wave function for H+ + H(2s). Figure 13. Elastic (solid, red) and transfer (dashed, green) differential cross section for H+ + H(2s) from ENDwave for E = 0.16 hartree (4.4 eV).
removed and the potential surfaces are adiabatic (or diabatic in practice) with conical intersections and Hamiltonian coupling elements. A third way, proposed in this work, is to describe the potential surfaces on which the nuclei move in terms of the time-dependent potentials calculated from END trajectories. This naturally leads to a center-of-mass Cartesian coordinate frame where the potential surfaces fold back on themselves near the classical turning points. Nuclear wave packet scattering can then be conducted easily and cheaply; the only requirements for transfer between different folds of the surface are the usual continuity and smoothness conditions on the wave function. An initial implementation of the wave packet scattering ENDwave has been demonstrated to provide qualitatively correct differential cross sections. More refinement of the implementation is required to predict the differential cross sections accurately, but the promise of the theory has been demonstrated.
Figure 12. Color map of the momentum distribution of the transferred (inelastic) wave function for H+ + H(2s).
elastic scattering. The width of the scattered packet is the same as that of the incident wave packet. The wave function on the second (nonincident) surface (Figure 12) shows a much different result. Here, most of the scattered packet is off-axis (the beam that is unscattered is also not transferred at all) and much broadened by the transfer; the center of the band is about 9.0 au−1, whereas on the elastic surface it is centered near 8.5 au−1. This conforms to expectations, as the higher momentum portions of the beam have greater probabilities of transferring. The stippling is the result of aliasing and mutual interference. As discussed in the previous sections, not much can be said about the differential cross sections (Figure 13). At these low energies (approximately 4.4 eV mean packet energy), not much transfer occurs. The elastic cross section falls off slowly but consistently with angle; the transfer is effectively constant (barring noise) out to quite high angles before falling off very rapidly. The effects of aliasing make it hard to tell which if any of the smaller structures are real and which are artifacts.
■ ■ ■
AUTHOR INFORMATION
Notes
The authors declare no competing financial interest.
ACKNOWLEDGMENTS We are grateful to the University of Florida department of Research Computing for providing the computational facilities. REFERENCES
(1) Delos, J. B. Theory of Electronic Transitions in Slow Atomic Collisions. Rev. Mod. Phys. 1981, 53, 287−357. (2) Hutson, J. M. Coupled-Channel Methods for Solving the BoundState Schrodinger-Equation. Comput. Phys. Commun. 1994, 84, 1−18. (3) Maierle, C. S.; Schatz, G. C.; Gordon, M. S.; McCabec, P.; Connorc, J. N. L. Coupled Potential-Energy Surfaces and Quantum Reactive Scattering for the Cl(2P) + HCl -> ClH + Cl(2P) Reaction. J. Chem. Soc., Faraday Trans. 1997, 93, 709−720. (4) Le, A.-T.; Lin, C. D.; Errea, L. F.; et al. Comparison of HyperSpherical versus Common-Reaction-Coordinate Close-Coupling Methods for Ion-Atom Collisions at Low Energies. Phys. Rev. A 2004, 69, 062703. (5) Solov’ev, E. A. The Advanced Adiabatic Aproach and Inelastic Transitions via Hidden Crossings. J. Phys. B: At. Mol. Opt. Phys. 2005, 38, R153−R194.
4. CONCLUSIONS As shown, the same scattering information can be described in multiple complementary ways. One description is that of the reaction coordinate, favored by chemists, where the reactant and product valleys are separated by a barrier of some type. This is conceptually simple, but the needed coordinate transformations are unwieldy and nonunique, leading to computational difficulties. Another description is the Born− Oppenheimer separation, where time dependence has been 6393
dx.doi.org/10.1021/jp500532d | J. Phys. Chem. A 2014, 118, 6385−6394
The Journal of Physical Chemistry A
Article
(6) Deumens, E.; Diz, A.; Longo, R.; Ö hrn, Y. Time Dependent Theoretical Treatments of the Dynamics of Electrons and Nuclei in Molecular Systems. Rev. Mod. Phys. 1994, 66, 917−983. (7) Newton, R. G. Scattering Theory of Waves and Particles, 2nd ed.; Texts and Monographs in Physics; Springer-Verlag: Berlin, 1982. (8) Hall, B.; Deumens, B.; Ö hrn, Y.; Sabin, J. R. Wave Packet Dynamics on Multiply Valued Potential Surfaces: Report on Work in Progress. Int. J. Quantum Chem. 2012, 112, 247−252. (9) Cabrera-Trujillo, R.; Ö hrn, Y.; Deumens, E.; Sabin, J. R. Absolute Differential and Total Cross Sections for Direct and Charge-Transfer Scattering of keV Protons by O2. Phys. Rev. A 2004, 70, 042705. (10) Malinovskaya, S.; Cabrera-Trujillo, R.; Sabin, J. R.; Deumens, E.; Ö hrn, Y. Dynamics of Proton-Acetylene Collisions at 30 eV. J. Chem. Phys. 2002, 117, 1103−1108. (11) Schiff, L. I. Approximation Method for High-Energy Potential Scattering. Phys. Rev. 1956, 103, 443−453. (12) Killian, B. On Electronic Representations in Molecular Reaction Dynamics. Ph.D. thesis, University of Florida, 2005. (13) Born, M.; Oppenheimer, J. R. Quantum Theory of Molecules. Ann. Phys. 1927, 389, 457−484. (14) Thouless, D. Stability Conditions and Nuclear Rotations in Hartree-Fock Theory. Nucl. Phys. 1960, 21, 225−232. (15) Bandrauk, A. D.; Shen, H. Improved Exponential Split Operator Method for Solving the Time-Dependent Schrodinger-Equation. Chem. Phys. Lett. 1991, 176, 428−432. (16) Park, T. J.; Light, J. C. Unitary Quantum Time Evolution by Iterative Lanczos Reduction. J. Chem. Phys. 1986, 85, 5870−5876. (17) Fevens, T.; Jiang, H. Absorbing Boundary Conditions for the Schrodinger Equation. SIAM J. Sci. Comput. 1999, 21, 255−282. (18) Proakis, J.; Manolakis, D. Digital Signal Processing: Principles, Algorithms, and Applications; Macmillan Publishing Co.: New York, 1992; p 401. (19) Iyengar, S. S.; Kouri, D. J.; Hoffman, D. K. Particular and Homogeneous Solutions of Time-independent Wavepacket Schrodinger Equations: Calculations Using a Subset of Eigenstates of Undamped or Damped Hamiltonians. Theor. Chem. Acc. 2000, 104, 471−483. (20) See Pindzola, M. S.; Lee, T. G.; Minami, T.; Schultz, D. R. Excitation and Charge Transfer in p + H(2s) Collisions. Phys. Rev. A 2005, 72, 062703 (2005), and references therein. (21) Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. 1. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023.
6394
dx.doi.org/10.1021/jp500532d | J. Phys. Chem. A 2014, 118, 6385−6394
