Symmetrical Windowing for Quantum States in Quasi-Classical

Article pubs.acs.org/JPCA

Symmetrical Windowing for Quantum States in Quasi-Classical Trajectory Simulations Stephen J. Cotton and William H. Miller* Department of Chemistry and Kenneth S. Pitzer Center for Theoretical Chemistry, University of California, Berkeley, California 94720, United States, and Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, United States ABSTRACT: A microscopically reversible approach toward computing reaction probabilities via classical trajectory simulation has been developed that bins trajectories symmetrically on the basis of their initial and final classical actions. The symmetrical quasi-classical (SQC) approach involves defining a classical action window function centered at integer quantum values of the action, choosing a width parameter that is less than unit quantum width, and applying the window function to both initial reactant and final product vibrational states. Calculations were performed using flat histogram windows and Gaussian windows over a range of width parameters. Use of the Wigner distribution function was also investigated as a possible choice. It was demonstrated for collinear H + H2 reactive scattering on the BKMP2 potential energy surface that reaction probabilities computed via the SQC methodology using a Gaussian window function of 1/2 unit width produces good agreement with quantum mechanical results over the 0.4−0.6 eV energy range relevant to the ground vibrational state to the ground vibrational state reactive transition.

I. INTRODUCTION AND HISTORY To keep the presentation below as simple as possible, so as to focus on the essential idea, all discussion will refer explicitly to a collinear A + BC(N1) → A + BC(N2) or AB(N2) + C inelastic or reactive scattering process, where N is the vibrational quantum number of the diatomic molecule. It should be clear how the procedure may be generalized. In the standard “quasi-classical” prescription for using classical trajectory calculations to describe state-to-state inelastic and reactive scattering processes,1 one begins trajectories with the vibrational energy of the diatom BC(N1) being equal to the quantum mechanical energy level of the initial vibrational state N1, and then “assigns” the trajectory to the final quantum state N2 whose energy level is closest to that of the actual final classical vibrational energy. It occurred to one of us many years ago2 that this unsymmetrical treatmenti.e., that the initial vibrational energy of the molecule is equal to the quantum mechanical value corresponding to the initial state, whereas this is not true for the final vibrational energy of the moleculemight be improved by having a symmetrical prescription (i.e., that obeys microscopic reversibility) for associating classical trajectories with initial and final quantum states. To accomplish this it is useful to begin with the “primitive” semiclassical expression3 for the state-to-state transition probability, for which both initial and final vibrational states have their correct (quantum) vibrational energy, ⎛ ∂n2(n1 ,q1) PN2 ← N1 = ⎜⎜2π ∂q1 ⎝

⎞ ⎟ ⎟ ⎠

where (n, q) are the classical action-angle variables for the vibrational degree of freedom; (n1, q1) are their initial conditions, which (together with initial conditions for the translational degree of freedom, which are not indicated explicitly) determine the final value of the action variable, n2(n1, q1). [Note: as usual in earlier work, the classical action is actually n + 1/2, so that (with ℏ = 1 throughout) n is the classical counterpart to the vibrational quantum number; i.e., the vibrational energy levels are the energies for which n = 0, 1, 2, ....] The appropriate initial conditions for the N1 → N2 transition probability in eq 1 are n1 = N1, with q1 determined by

n2(N1 ,q1) = N2

(2)

The Jacobian, ∂n2 (N1,q1)/∂q1, in eq 1 is evaluated at the root of eq 2 (with a sum over multiple roots). It is not apparent from eq 1 that this gives a symmetric transition probability matrix, but it is well-known from semiclassical theory to be so; we refer to earlier work3 for more details of this. The standard (asymmetrical) quasi-classical model can be obtained by averaging eq 1 over a quantum number width of the final action variable about its quantum value, Special Issue: Joel M. Bowman Festschrift

−1

© 2013 American Chemical Society

Received: January 30, 2013 Revised: February 20, 2013 Published: February 24, 2013

(1) 7190

dx.doi.org/10.1021/jp401078u | J. Phys. Chem. A 2013, 117, 7190−7194

The Journal of Physical Chemistry A

Article

Figure 1. Ground-state to ground-state reaction probability, P0←0, versus total energy (eV) for the collinear H + H2 reaction over the threshold energy region: exact QM (+); traditional (asymmetric) QC (×); symmetric QC (∗).

PNQC 2 ← N1

=

N2 + 1/2

∫N −1/2 2

⎛ ⎞−1 ∂ n 2 ⎟ dn2 ⎜⎜2π ∂q1 ⎟⎠ ⎝

average over the initial Cartesian vibrational coordinate and momentum rather than their action-angle variables.] Another reason for being interested in this symmetrical “double averaging” approach was the thought that it might smooth out sharp classical reaction thresholds somewhat and perhaps mimic, at least qualitatively, the threshold smoothing that quantum tunneling provides. At the Faraday Discussion where this suggestion was first made, results for the prototypical H + H2(0) → H2(0) + H collinear reaction were presented in the General Discussion5 following the published papers. Figure 1 shows an updated version of these results recalculated for the BKMP2 H + H2 potential energy surface (PES)6 over the threshold energy region. The standard quasi-classical result (× points) shows a sharp reaction threshold, which is smoothed out in the correct quantum result (+ points) due to tunneling. The symmetrical, double-averaged quasi-classical result (∗ points) is indeed seen to smooth out the sharp classical threshold, but much too much, and the approach was thus not pursued. [Note that over this threshold energy region, P0←0(E) is equal to N(E), the cumulative reaction probability, because the only open reactive channel is between ground vibrational states. Thus, the QM curve displayed in Figure 1 was generated via a direct calculation of N(E) using a discrete variable representation (DVR) of the QM Green’s function with absorbing boundary conditions (ABC).7−10]

(3a)

and by changing the variable of integration from n2 to q1, this gives = PNQC 2 ← N1

1 2π

∫0

2π

⎛1 ⎞ dq1 h⎜ − |n2(N1 ,q1) − N2|⎟ ⎝2 ⎠

(3b)

where h(x) = 0 (x < 0) or 1 (x > 0) is the usual Heaviside function. (The reader familiar with recent developments in semiclassical theory will recognize this change of integration variables to be the same “trick” that leads to the initial value representation4 of SC theory.) Thus, the initial action variable is fixed at its quantum value N1, and the conjugate angle variable q1 is averaged over. (It is clear that one can compute transition probabilities to all final states N2 from the same ensemble of trajectories by simply accumulating the trajectories that land in the various histogram “boxes”.) It is rather obvious, though, that one can obtain a symmetrical result by averaging eq 1 also over the initial action variable by a quantum number width about the quantum value N1, just as done for the final action; this gives2 PNSQC 2 ← N1

=

N1+ 1/2

∫N −1/2 1

dn1

N2 + 1/2

∫N −1/2 2

−1 ⎛ ∂n2 ⎞⎟ ⎜ dn2 ⎜2π ∂q1 ⎟⎠ ⎝

∞ 2π ⎛1 ⎞ 1 dq1 dn1 h⎜ − |n1 − N1|⎟ · = ⎝ ⎠ −1/2 2π 0 2 ⎛1 ⎞ h⎜ − |n2(n1 ,q1) − N2|⎟ ⎝2 ⎠

∫

II. SYMMETRICAL WINDOWING IN GENERAL

∫

In more recent years, however, Bonnet and Rayez11 have shown that considerably improved quasi-classical results for product state distributions for inelastic and reactive scattering can be obtained by using “Gaussian binning” rather than the traditional histogram bins. Czako and Bowman12 have also made good use of such methods. These authors typically find the best results when the bin sizes are significantly smaller than the unit quantum number width used in eqs 3 and 4. The symmetrical version of this approach is simply to use a window functionGaussian or otherwisefor both the initial and the final vibrational state. The purpose of this paper is to

(4)

which is seen to be an average over the initial phase space of the vibrational degree of freedom, with the integrand being the product of normalized window f unctions for the initial and final vibrational states. [Note that 2π (equivalent to 2πℏ here because ℏ = 1) is the usual “volume of phase space” for one degree of freedom, and also that because phase space averages are invariant to a canonical transformation, one could do the phase space 7191



Article

Figure 2. (A) Same as Figure 1, but with results for P0←0(E) = N(E) (equal at these energies) calculated via the symmetric quasi-classical (SQC) approach using histogram window functions of various widths: exact QM (+); SQC with Δn = 0.4 (×), Δn = 0.5 (∗), and Δn = 0.6 (□). (B) Same as Figure 2A, but showing N(E) over an expanded energy (eV) range (where P0←0(E) ≠ N(E)): exact QM (+); SQC with Δn = 0.4 (×), Δn = 0.5 (∗), and Δn = 0.6 (□).

show how this works for the prototypical H + H2(0) → H2(0) + H reaction for various possible choices of window functions.

(again calculated using the BKMP2 PES6), and one sees that a value of Δn = 1/2 gives quite reasonable results over the threshold energy region shown in Figure 2A: the maximum reaction probability in Figure 2A is in reasonable agreement with the quantum value, as is the general shape of the reaction probability curve over the threshold region. The reaction probabilities shown in Figure 2A for various Δn are shown over a wider energy range in Figure 2B. More precisely, Figure 2B plots the cumulative reaction probability N(E), which, of course, is just equal to the sum of all state-to-state reaction probabilities Pp←r(E), which are nonzero over the higher energies of this expanded range. Here, one sees that the results given by the symmetric binning procedure using various width parameters exhibit a bit of an oscillation about the plateau region of the QM N(E) curve before climbing in the high-energy region where additional reactive channels open up. Next, following the lead of Bonnet, Bowman, and their respective co-workers, we employ in eq 6 a normalized Gaussian window function defined by

III. EXAMPLES OF SYMMETRICAL WINDOWING APPLIED TO H + H2 REACTIVE SCATTERING We first consider a square histogram window function as above, but with various widths Δn about the quantum value; the normalized histogram window function for state N is thus WN (n) =

⎞ 1 ⎛⎜ Δn h − |n − N |⎟ ⎝ ⎠ 2 Δn

(5)

and the N1 → N2 transition probability is given by 1 2π

∫0

2π

∞

dq1

∫−1/2 dn1 WN (n1)·WN (n2(n1,q1)) 1

2

(6)

The 40-year old result recalculated in Figure 1 corresponds to Δn = 1, the maximum possible value, which was seen to be too large. Figure 2A,B above show results for smaller values of Δn 7192



Article

Figure 3. (A) N(E) calculated via the SQC method over the expanded energy (eV) range (where P0←0(E) ≠ N(E)) using histogram (+) and Gaussian (×) window functions with width parameter Δn = 0.5. (B) Analogue of Figure 2B: N(E) calculated via the symmetric quasi-classical (SQC) approach using various width parameters Δn, but here using Gaussian window functions instead of histogram windows: exact QM (+); SQC with Δn = 0.4 (×), Δn = 0.5 (∗), and Δn = 0.6 (□). (C) Results for P0←0(E) = N(E) calculated via the symmetric quasi-classical (SQC) approach using Gaussian window functions having various width parameters Δn: exact QM (+); SQC with Δn = 0.4 (×), Δn = 0.5 (∗), and Δn = 0.6 (□). (D) Results for P0←0(E) = N(E): exact QM (+); traditional (asymmetric) QC (×); SQC with histogram window of width Δn = 1 (∗); SQC with Gaussian window of width Δn = 0.5 (□).

Figure 4. Results for calculation of P0←0(E) using the SQC method with the Wigner distribution function as the classical action window (×) compared with the exact QM result (+).

7193


The Journal of Physical Chemistry A W NG(n) =

Article

2 1 e−[(n − N )/(Δn /2)] π (Δn/2)

to “quantize” the most relevant and quantum-like degree(s) of freedom in the process via binning of initial and final action variables, e.g., the H atom vibration in a H atom transfer reaction in a large molecular system.

(7)

The result for Δn = 1/2 is shown in Figure 3A and compared to the result obtained by using the symmetric histogram binning procedure, also with Δn = 1/2. Although the use of Gaussian windows does seem to give somewhat smoother results than the use of histogram windows, the choice of the optimal width parameter Δn [see, e.g., the variation in Figure 2B (histogram windows) and Figure 3B (Gaussian windows)] is seen to be much more important than the choice of the functional form of the window itself for the symmetric quasi-classical method to give the best approximation to the true QM results. Nevertheless, despite the more significant effect of varying width parameter, the best overall SQC approximation found here was obtained by choosing a Gaussian window having an optimal width of Δn = 1/2. Over the reaction threshold energy region shown in Figure 3C, results for P0←0(E) computed using the SQC method with this window function show strikingly good agreement with the true QM result: in terms of the onset of the tunneling region, the location and value of maximum reaction probability, as well as the overall curve shape P0←0(E). Figure 3D overlays this result for Δn = 1/2 with the analogue of the 40-year old result displayed in Figure 1. Finally, we also considered using the Wigner distribution for the initial and final action window functions. The Wigner distribution (in Cartesian variables) for a given quantum state is ∞

WN (p ,x) =

⎛

⎞

⎛

■

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Science Foundation Grant No. CHE-1148645 and by the Director, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, U.S. Department of Energy, under Contract No. DE-AC02-05CH11231.

■

⎟

⎜

⎞ ⎟

(8)

where ϕN(x) is the wave function corresponding to quantum state N. For a 1D-harmonic oscillator, expressed in terms of action-angle variables, eq 8 gives W0(n) = 2e−2[n + (1/2)]

REFERENCES

(1) See, for example: Karplus, M.; Porter, R. N.; Sharma, R. D. Exchange Reactions with Activation Energy. I. Simple Barrier Potential for (H, H2). J. Chem. Phys. 1965, 43, 3259−3287. (2) Miller, W. H.; Raczkowski, A. W. Partial Averaging in Classical SMatrix Theory. Vibrational Excitation of H2 by He. Faraday Discuss. Chem. Soc. 1973, 55, 45−50. See specifically p 47. (3) Miller, W. H. Classical S Matrix: Numerical Application to Inelastic Collisions. J. Chem. Phys. 1970, 53, 3578−3587. (4) Miller, W. H. The Semiclassical Initial Value Representation: A Potentially Practical Way for Adding Quantum Effects to Classical Molecular Dynamics Simulations. J. Phys. Chem. A 2001, 105, 2942− 2955. (5) Miller, W. H.; Raczkowski, A. W. Partial Averaging in Classical SMatrix Theory. Vibrational Excitation of H2 by He. Faraday Discuss. Chem. Soc. 1973, 55, 71−72. (6) Boothroyd, A. I.; Keogh, W. J.; Martin, P. G.; Peterson, M. R. A Refined H3 Potential Energy Surface. J. Chem. Phys. 1996, 104, 7139− 7152. (7) Details of the DVR-ABC Green’s function method are found in refs 8, 9, and 10. In particular, see eq 2.9(d) of ref 8 for the expression used to directly calculate N(E). Equation 3.2 of ref 10 was used as the analytical form of the reactant and product adsorbing potentials. (8) Seideman, T.; Miller, W. H. Quantum Mechanical Reaction Probabilities Via a Discrete Variable Representation-Absorbing Boundary Condition Green’s Function. J. Chem. Phys. 1992, 97, 2499−2514. (9) Colbert, D. T.; Miller, W. H. A Novel Discrete Variable Representation for Quantum Mechanical Reactive Scattering via the S-Matrix Kohn Method. J. Chem. Phys. 1992, 96, 1982−1991. (10) Manthe, U.; Miller, W. H. The Cumulative Reaction Probability as Eigenvalue Problem. J. Chem. Phys. 1993, 99, 3411−3419. (11) Bonnet, L.; Rayez, J. C. Quasiclassical Trajectory Method for Molecular Scattering Processes: Necessity of a Weighted Binning Approach. Chem. Phys. Lett. 1997, 277, 183−190. (12) Czako, G.; Bowman, J. M. Quasiclassical Trajectory Calculations of Correlated Product Distributions for the F+CHD3(v1 = 0,1) Reactions Using an Ab Initio Potential Energy Surface. J. Chem. Phys. 2009, 131, 244302.1−18.

∫−∞ dΔx eipΔxϕN ⎝x + Δ2x ⎠ ϕN ⎝x − Δ2x ⎠ ⎜

AUTHOR INFORMATION

Corresponding Author

(9)

for the ground vibrational state. The result of using this Wigner window function in eq 6 is shown in Figure 4. Because eq 9 gives a very broad distribution for the action variable, and one weighted heavily at small valuesrather than having the action localized about the integer quantum value (as is the case above for Gaussian or histogram window functions with Δn = 1/2) the results unsurprisingly give a rather poor approximation to the correct quantum transition probabilities.

IV. SUMMARY A microscopically reversible approach toward computing reaction probabilities via classical trajectory simulation has been developed that bins trajectories symmetrically on the basis of their initial and final classical actions. The symmetrical quasiclassical (SQC) approach involves defining a classical action window function centered at integer quantum values of the action, choosing a width parameter that is less than unit quantum width, and applying the window function to both initial reactant and final product vibrational states. It is demonstrated, for collinear H + H2 reactive scattering, that reaction probabilities computed via the SQC methodology using a Gaussian window function of 1/2 unit width produces good agreement with true quantum mechanical results over the energy region relevant to the ground-vibrational state to ground-vibrational state reactive transition. Finally, we note that when applying these approaches to a large (polyatomic) molecular system, it is likely that one will only need 7194


Symmetrical Windowing for Quantum States in Quasi-Classical

Recommend Documents