8212
J. Phys. Chem. 1991, 95, 8212-8215
Complex Log Derivative Method for Nonreactive Coupled-Channel Scattering Calculations Danny L. Yeagert and William H. Miller* Department of Chemistry, University of California, and Chemical Sciences Division, Lawrence Berkeley Laboratory, Berkeley, California 94720 (Received: January 15, 1991)
Downloaded by CENTRAL MICHIGAN UNIV on September 11, 2015 | http://pubs.acs.org Publication Date: October 1, 1991 | doi: 10.1021/j100174a037
A new version of Johnson's log derivative method for integrating the (nonreactive) Schrodinger equation is described. The new feature is that the log derivativeequation is integrated inward (from large r to small r) with initial conditionsthat correspond to asymptotically incoming (or outgoing) radial waves. The result is that the log derivative function Q(r) (which is now a complex function of r) is a very smooth, nonsingular function (compared to the oscillatory, singular behavior of Q(r) with the usual real boundary conditions), meaning that many fewer grid points are required for its numerical integration. The number of grid points necessary to achieve a given level of accuracy is also only a weak function of the scattering energy E or the mass of the particles. Test calculations for a one-dimensional and a multichannel example illustrate this behavior.
I. Introduction
In the process of recent work in our group on reactive scattering processes in molecular collision^,^-^ a new approach has been suggested to us for carrying out calculations for the simpler process of nonreactive (i.e., inelastic) scattering. The purpose of this paper is to describe this new (so far as we are aware) approach. The method we present is based on the log derivative approach of Johnson,6s7 the difference being that we integrate the radial equation inward (from large r to small r ) with initial conditions (at large r) that correspond to purely incoming (or outgoing) radial waves. The log derivative is therefore a complex function of r. The important feature of using these boundary conditions is that the log derivative is an extremely smooth function of r (rather than oscillatory and singular, as in the usual approach), meaning that many (- an order of magnitude) fewer grid points are needed for the numerical integration. The procedure for extracting the S matrix from the calculations is also different from that in the usual approach. Section I1 describes the approach for the case of one-channel (i.e., potential) scattering, which illustrates its essential features, and then section 111 shows how it generalizes to the multichannel case. Section I1 also describes an application to potential scattering that illustrates the dramatic reduction in grid points that is possible, and section 111 describes an application to inelastic scattering. 11. Method: Potential Scattering A. Theory. For s wave potential scattering the Schrodinger
equation is of the standard form
-hZ d2
l-0
--
and the asymptotic form as r m is lim \k(r) -e-ikr + e W
Q'(r) = -Q(r)*
lim u(r)
-
e+
+ W(r)
(2.7a)
where 2m
W )= j - p r )
- E)
(2.7b)
From eq 2.5, one sees that the large r limit of Q(r) is lim Q(r) = -ik
(2.8)
I---
The procedure, therefore, is to solve eq 2.7a for Q(r), with the boundary condition eq 2.8 at r = m. (More will be said about this below.) The function u(r) is then obtained by integrating eq 2.6
W = 4 r 1 )exp[
1'
dr'
QW]
(2.9a)
where rl is a large value of r; thus u ( r l ) = e-ik
(2.9b)
u(r) is one solution of the Schrodinger equation; u(r)* (it's complex conjugate) is another linearly independent solution. The asymptotic form of u(r)* corresponds to a purely outgoing wave
lim u(r)* I---
- dkr
(2.10)
The physically correct solution of the Schrodinger equation, *(r), is a linear combination of u(r) and u(r)*, and in light of the various asymptotic forms (eqs 2.3, 2.5, and 2.10) it is clear that this linear combination is @ ( r ) = -u(r) u(r)*S (2.1 1) (2.3)
The ( 1 X 1) S matrix value is a number of unit modulus S = e2in (2.4) where 9 is the (real) phase shift. Rather than attempting to solve for \k(r) directly, we seek first the (irregular) solution of eq 2.1, u(r),that is asymptotically a purely incoming wave, Le. r-+-
Q(r) = W ) / W (2.6) and as Johnson shows, the Schrcdinger equation (eq 2.1) for u(r) leads to the following equation for Q(r)
(2.1)
where \k(r) is required to be regular at r = 0, i.e. lim \k = 0
?.-
We utilize Johnson's log derivative algorithm6v7to solve for u(r). The log derivative Q(r) is defined as
(2.5)
'Permanent address: Chemistry Department, Texas A & M University, College Station, TX 77843-3255.
+
-
The value of S is obtained by requiring that \k(r) be regular as r 0 (Le., eq 2.2) whereby eq 2.1 1 implies that S = lim u(r)/u(r)* (2.12) f-0
(1) Zhang, J. 2.H.; Chu, S. 1.; Miller, W. H. J . Chem. Phys. 1988, 88, 6233. ( 2 ) Miller, W. H.; Jansen op de Harr, B. M. D. D. J . Chem. Phys. 1987, 86, 6213. (3) Miller, W. H. J . Phys. B: At., Mol. Opt. Phys. 1988, 22, 115. (4) Zhang, J. H. Z.; Miller, W. H. Chem. Phys. Lett. 1987, 140, 329. (5) Zhang, J. H. Z.; Miller, W . H. J . Chem. Phys. 1989, 91, 1528. (6) Johnson, B. R. J . Comput. Phys. 1973, 13, 445. (7) Johnson, B. R. Proceedings of the NRCC Workshop on Algorithms and Computer Codes in Atomic and Molecular Quantum Scattering Theory; Lawrence Berkeley Laboratory Report LBL-9501; 1979; Vol. 1.
0022-3654/91/2095-8212$02.50/0 0 1991 American Chemical Society
The Journal of Physical Chemistry, Vol. 95, NO.21, 1991 8213
Nonreactive Coupled-Channel Scattering Calculations With u(r) given by eq 2.9, eq 2.12 becomes S
+
= lim exp[-2ikrl ro
1'
dr' (Q(r?
=
I
+2il:
lim exp[-2ikrl f-0
i.e., the phase shift
-
- Q(r)*)]
TABLE I: Phase Shift for tbe One-Dimensional Potential V = e-" step size qa qb step size qa ,Ib
9
9
dr' Im Q(r')] (2.13a)
of eq 2.4 is
+
= -krl
-
1'
dr' Im Q(r')
(2.13b)
with rl and r 0. In summary, therefore, one solves eq 2.7a for the log derivative Q(r), with the boundary condition eq 2.8, and then the S matrix is given by eq 2.13. To see the potential power of this approach, we note that an approximate solution of eqs 2.7 with boundary condition (2.8) is Q(r) -ik(r) (2.14a)
-
Downloaded by CENTRAL MICHIGAN UNIV on September 11, 2015 | http://pubs.acs.org Publication Date: October 1, 1991 | doi: 10.1021/j100174a037
where k(r) is the local wave vector (2.1 4b)
Within this approximation one thus has Im Q(r) = -k(r) r > r,
=0 r
-0.67342 -0.673 43 -0.67358 -0.67408
-0.155 14 -0.155 15 -0.15531 -0.15810
-0.2 -0.3 -0.4
-0,67604 -0.34987 -0.67942 -0.68627
'k = Suo-'. b k = 50ao-'.
I
k(r) = [2m(E - V(r))/h2]1/2
4.001 -6.05 -0.01 -0.1
< rTp
where rTp is the classical turning point, and one recognizes that eq 2.13b then gives the WKB approximation for the phase shift. This shows that Q(r) can be expected to be a very smoothfunction of r, and also that the limit r 0 in eq 2.13 is achieved very quickly as soon as r is less than rTp. We note that this behavior of Q(r) is in strong contrast to that found if the usual (real) boundary conditions are used.6 In the usual case, Q(r) is an oscillatory function of r, and in fact has singularities (two per wavelength). We thus expect that the log derivative equation, eq 2.7, can be integrated with many fewer grid points with eq 2.8 as the boundary condition rather than the usual ones. To conclude this section, we describe an efficient algorithm for integrating the log derivative equation, eq 2.7. Using the trapezoid gives rule to integrate eq 2.7 from r, to Ar Q ~ + I - Qn = +T(-Qn+i2 - QnZ + wn+l + wn) (2.15)
-
of H atom, E = 0.185 eV) and k = Souo-' (e.g., m = mass of argon, E = 0.463 eV). The results obtained by use of the complex log derivative method are given in Table I. For the former (smaller) value of k the complex log derivative method (integrated via the trapezoid algorithm, eq 2.16) yields the phase shift corrm to 0.1% with a step size 1Ar1 = 0 . 1 and ~ ~a phase shift correct to better than 1% with a step size lArl = 0 . 3 ~In ~contrast, the renormalized Numerov algorithm (which is essentially equivalent to the log derivative with the usual real boundary conditions) requires lArl = 0 . 0 5 ~to~yield an accuracy in the phase shift correct to 0.1%. For the latter (larger) value of k the complex log derivative algorithm yields a phase shift correct to 0.1%with a step size of 1 4 = 0 . 0 5 and ~ ~ a phase shift correct to better than 2% with a step size of 0.1%. The renormalized Numerov algorithm requires 1Ar1 = 0 . 0 0 5 ~for~ k = 50ao-' for an accuracy of 0.1%. (The number of grid points per wavelength of the wave function, -2~/klAr(, is -25 in both cases for the renormalized Numerov approach.) The above example illustrates another important practical feature of the complex log derivative method: the required number of grid points does not change dramatically with energy or mass-it is essentially the number of points required to evaluate the integral
1; dr [2m(E - V(r))/h2I1l2
(2.19)
by the trapezoid rule. The usual log derivative method with real boundary conditions, on the other hand, requires a given number of grid points per wavelength, so that the number required increases
111. Generalization to Multichannel Scattering A. Theory. Everything in the previous section generalizes to the multichannel case, provided no closed channels are included. The coupled-channel Schrodinger equation is of the form
where Ar = r,+l - r,, and
= w(rk) for k = 1, 2, .... Equation 2.1 5 is a quadratic equation for QHI, which can be solved to give Qk
Q(rk)
Qn+l = 2A/[1
wk
+ (1 + 2ArA)1/2]
where V(r) = (Vlp(r)}is the potential matrix, E = E,& (El is the translational energy in channel 1), and *(r) has the standard boundary conditions
(2.16a) (3.lb)
with Ar
(2.16b)
The (irregular) function matrix u(r) is thus chosen to have the asymptotic form
Equation 2.16 defines a two-point iteration scheme: one begins at rl = large (2.17a)
(3.2)
A = Q n + ~ ( w n + i + Wn-QnZ)
QI = -ik
(2.17b)
and then propagates inward: thus Ar = r,,+l-r, = 4 4 in eq 2.16. (We have tried a number of other integration schemes, e.g. Simpson's rule, but find the one described above to be the most stable.) B. Example. To illustrate the effectiveness of this complex log derivative method, we tested it first on a simple repulsive potential V(r) = Ae-2'
(2.18)
with two values of k = (2mE/h2)lI2,k = 5ao-I (e.g., m = mass
so that *(r) is then given by W(r) = -u(r) Regularity of q ( r ) as r
-
+ u(r)*S
(3.3)
0, Le.
lim * ( r ) = 0 l-0
implies that
S = lim u(r)*-Iu(r) l-0
(3.4)
The log derivative matrix Q(r) is defined by Q(r) = u'(r) u(r)-I
(3.5)
Yeager and Miller
8214 The Journal of Physical Chemistry, Vol. 95, No. 21, 1991
and the equation for Q'(r) is
TABLE 11: Inelastic Transition Probabilities for the Potenthi of Eq
QW = *w2 + W ( r )
(3.6a)
where 2m W ( r ) = Y ( V ( r )- E)
h
The boundary condition for Q ( r ) at r
-
m
(3.6b) is
lim Q ( r ) = -ik
F
(3.7)
m
where ( k ) , = blr ( 2 m E l / h 2 ) 1 / 2Q. ( r ) is a complex symmetric matrix. Equation 3.6 is thus to be integrated inward with eq 3.7 as the initial condition at rl = large. Integrating (3.6a) by the trapezoid rule gives Ar Q ~ +-IQn = - Q? + Wn+I + Wn) (3.8a)
z(Qn+12
'
3.1V
step size -0.01 -0.1 -0.2 -0.3 -0.4
Kohnb step size -0.01 -0.1 -0.2 -0.3 -0.4
Kohnb
Pol
Po2
0.892(-1) 0.893(-1) 0.894(-1) 0.897(-1) 0.901(-1) 0.892(-1)
0.441(-1) 0.441(-1) 0.442(-1) 0.443(-1) 0.444(-1) 0.441(-1)
0.202(-1) 0.202(-1) 0.203(-1) 0.203(-1) 0.204(-1) 0.202(-1)
PO4
p12
0.134(-1) 0.134(-1) 0.134(-1) 0.134(-1) 0.135(-1) 0.134(-1)
0.108 0.108 0.109 0.109 0.109 0.108
PI3
PI4
p13
p24
5 4
0.634(-1) 0.634(-1) 0.635(-1) 0.636(-1) 0.638(-1) 0.634(-1)
0.416(-1) 0.416(-1) 0.416(-1) 0.416(-1) 0.417(-1) 0.416(-1)
0.163 0.164 0.164 0.164 0.164 0.163
0.123 0.123 0.123 0.123 0.123 0.123
0.334 0.334 0.334 0.334 0.334 0.334
"All channels are open with E = 0.25. The on-diagonal E elements are (1.25, 1.0, 0.75, 0.50, 0.25) with corresponding channel momenta (1.58, 1.41, 1.22, 1.00, 0.71). *Reference 8.
Downloaded by CENTRAL MICHIGAN UNIV on September 11, 2015 | http://pubs.acs.org Publication Date: October 1, 1991 | doi: 10.1021/j100174a037
where note that in practice r,, has to be taken only so small for all the eigenvalues of Q , to be essentially real, because in this limit
Wk = W(rk) = Q(rk) This can be rearranged to give Qk
ei&P)Q.*e(b/2)Qn =1
In the iteration scheme the right side of eq 3.8b is known, and Q H l is to be determined. Let U H Ibe the (complex) orthogonal transformation that diagonalizes the (complex symmetric) matrix Qn+l. Un+Ialso diagonalizes Qn+12,so that Un+lis the transformation that diagonalizes the left side of eq 3.8b. Thus Un+l must also be the transformation that diagonalizes the right side of eq 3.8b. The iteration procedure is now clear: with Q, known, one diagonalizes the right side of eq 3.8b to obtain eigenvalues aHI(o and the transformation matrix U H l .If qn+l(oare the eigenvalues of Qn+l,then the diagonalized version of eq 3.8b is Ar (3.8~) qn+I(" + l q n + I ( ' = an+I('
and this factor thus does not contribute to the S matrix for larger values of n (smaller values of r,,). Finally, we note that since the transformation matrix Uk that diagonalizes Qk is known (from the iterative solution of the log derivative equation), the exponentials of the matrices in eq 3.12 are conveniently evaluated as ehrQk = U k e h U k T
(3.14)
e-hrQk' = Uk*e-P"lk'Uk*T
(3.15)
and
B. Example. As an example of the complex log derivative method we consider an attractive five-channel potential first used by Lucchese for scattering calculations using the complex Kohn variational principle.8 For this example the potential is
which can be solved for the eigenvalues qn+l(')
Vqs(r)= Apse-'
+ + 2 A r ~ , , + ~ ( o ) ~ /(3.9) ~]
qn+l(o = 2 ~ , , + ~ ( o / [ l (1
= Un+lqn+lUn+lT
un = u(rn)
then successive use of the first-order Magnus approximation gives un = e(&/2)QnehrQ~l,.,eMe(&/2)Qlul (3.1 la) with uI = e-hv-l/2
(3.11b)
where v = h k / m and k are diagonal. The S matrix, defined by eq 3.4, is thus given by S = lim r.4
0; -2.00
(3.10)
where qn+l is the diagonal matrix of eigenvalues qn+l(o. To summarize the trapezoid integration algorithm for the log derivative matrix Q(r): the initial value (at r l ) is given by eq 3.7; to obtain Q H l from Q,,,one diagonalizes the right side of eq 3.8b, yielding the transformation matrix U,,+l (i.e., the matrix with the eigenvectors as columns) and eigenvalues Q d I is then given by eqs 3.9 and 3.10. The wavefunction u(r) is then obtained by integrating eq 3.5. If
ei&/2Woe(&/2)Qa.. .~(&/~)QI~-W (3.1 V 2) -I/~
where we remember that Ar = -1Arl and r, = rl + (n - 1)Ar. We
(3.16)
with the potential coupling coefficients given by
The matrix Qn+Iis then given by Qn+l
(3.13)
A =
0.10
0.40 -1.75 0.40
0.30 0.20
0.30
0.20
0.40
0.30
-1.50 0.40 0.30
0.40 -1.25 0.40
0.10 0.20 0.30 0.40 -1.00
)
(3.17)
and m = 1.0. The diagonal E matrix elements are chosen as (1.25, 1 .OO,0.75, 0.50, 0.25) and (2.00, 1.75, 1.50, 1.25, 1.00). (These correspond to E = 0.25 and 1.OO respectively in Lucchese's nomenclature.) The channel momentum is given by kl = (2EI)1/2or (1.58, 1.41, 1.22, 1.00, 0.71) and (2.00, 1.87, 1.73, 1.58, 1.41), respectively. The results obtained by use of the complex log derivative method are listed in Tables I1 and 111. With E = 0.25, the maximum error in an inelastic probability matrix element is about 1% with a step size of -0.4. With E = 1.0, the maximum error is 2.3% with a step size of -0.4. As for the potential scattering case in section IIB, we would not expect the step size necessary to achieve this level of accuracy to change significantly for higher energies or masses (Le., larger wavevectors). C. Asymptotically Closed Channels. The above theory must be modified in order to include asymptotically closed channels. The wavefunction matrix is thus written as * ( r ) = -ua(r) + ub(r)S (3.18) (8) Lucchese, R. L. Phys. Rev. A 1989, 40, 6879.
J. Phys. Chem. 1991,95,8215-8219 TABLE Ilk Iaclrstic Transition Probabilities for the Potential of Eq
(3.23a)
3.16”
step size -0.01 -0.1 -0.2 -0.3
-0.4 Kohnb
Pol
Po2
0.755(-1) 0.408(-I) 0.755(-1) 0.408(-1) 0.757(-1) 0.409(-1) 0.760(-1) 0.41 I(-1) 0.764(-1) 0.413(-1) 0.7551-1) 0.4081-1 I
p03
PM
PI 2
0.189(-1) 0.189(-1) 0.190(-1) 0.191(-1) 0.192(-1) 0.189(-1)
0.854(-2) 0.855(-2) 0.859(-2) 0.865(-2) 0.874(-2) 0.854(-2)
0.923(-1) 0.924(-1) 0.926(-1) 0.930(-1) 0.934(-1) 0.923(-1)
~~~
step size -0.01 -0.1 -0.2 -0.3 -0.4
Kohnb
PI3
0.552(-1) 0.553(-1) 0.554(-1)
0.557(-1) 0.560(-1) 0.553(-1)
p14
pz3
p14
4 4
0.284(-1) 0.284(-1) 0.285(-1) 0.287(-1) 0.289(-1) 0.284(-1)
0.126 0.126 0.126 0.127 0.127 0.126
0.795(-1) 0.796(-1) 0.798(-1) 0.802(-1) 0.808(-1) 0.795(-1)
0.182 0.182 0.183 0.183 0.185 0.182
Downloaded by CENTRAL MICHIGAN UNIV on September 11, 2015 | http://pubs.acs.org Publication Date: October 1, 1991 | doi: 10.1021/j100174a037
“All channels are open with E = 1.00. The on-diagonal E elements are (2.00, 1.75, 1.50, I .25, 1 .00)with corresponding channel momenta (2.00, 1.87, 1.73, 1.58, 1.41). bReference 8.
where the function matrix u,(r) is chosen to have the asymptotic form (3.19) The upper left block of u, is the open channekpen channel block and has the asymptotic form given by eq 3.1 lb. The linearly independent function ub has the asymtotic form (3.20) where for both functions the closed-channel wave vector is
Thus, u, and ub are not complex conjugate matrices and this means that they must be individually propagated inward from large r, and the S matrix is then given by
S = lim Ub(f)-ha(r) 4
(3.22)
-
The algorithm for propagation of the log derivative of either u, or ub is identical with eqs 3.8-3.10. The boundary condition for Q ( r ) as r rl (Le. m) is
8215
for ua, and (3.23b) for ub In some test applications, however, we have had difficulty in obtaining satisfactory results when closed channels are included. Integration of the equations for the log derivative matrices Q, and Qb, with the asymptotic initial conditions of eq 2.23, is quite well-behaved, but numerical instabilities arise when constructing the S matrix via eq 3.12. For the present we have not pursued this problem because for current purposes-Le., constructing inelastic distorted waves for the “free” functions u&) in the S matrix Kohn method for reactive scatterings-we are only interested in the open-channel problem. Further work would be useful, however, in developing a stable procedure for including closed channels in this approach. IV. Concluding Remarks The essential advantage of the complex log derivative method described in this paper is that the log derivative function Q(r) is a very smooth, nonsingular function of r, much more so than the log derivative function that results with the usual real boundary conditions. The equation for Q(r) can thus be integrated with fewer grid points than is necessary in the usual case, and the number of grid points required for a certain degree of accuracy is approximately independent of the energy and the mass of the particles. The algorithm generalizes naturally to the multichannel case. Test calculations for a one-channel example and a multichannel problem (open channels only) illustrate these features. Further work is needed to obtain a numerically stable version of the method when closed channels are included in a multichannel calculation. Acknowledgment. D.L.Y. thanks Professor John Zhang and Professor Moshe Shapiro for very stimulating discussions on the subject of this paper. We also thank Professor Robin Lucchese for providing detailed numerical results obtained by use of the Kohn variational principle in Tables I1 and 111. We acknowledge the National Science Foundation (Grant No. CHE-8920690) for generous support. D.L.Y.also thanks the Robert A. Welch Foundation (Grant A-770) for research support).
Anisotropic Integral Cross Sections for Li2(Y =0, Y =20)-Na Collisions H . 4 . Rubahn MPI fur Striimungsforschung, Bunsenstrasse IO. 0-3400 Gottingen, Germany (Received: January 30, 1991) Integral cross sections for Li2(u=0,v=20;J)-Na elastic, inelastic, and reactive scattering at E = 300 and 500 meV have been calculated via quasiclassical trajectories on a semiempirical valence bond energy surface. The cross sections show a pronounced dependence on alignment of the diatomic bond axis with respect to the initial relative velocity vector. The effect increases with increasing vibrational quantum number and is also observed in the exit channels of the collisions. Experimental data on the alignment dependence of the total Li2(o=O)-Na integral cross section show good agreement with the theoretical results. 1. Introduction Anisotropies in the collisional potentials between atoms and molecules in the vibrational ground state have investigated thoroughly in the past via the measurement of rotationally inel&c differential cross sections.’ Those studies lead mainly to a refined (1) Faubel, M. Ado. Ar. Mol. Phys. 1983, 19, 345.
0022-3654/91/2095-8215$02.50/0
understanding of the short-range anisotropic part of the interaction potential. Recently, it was shown in a laser moleculear beam study that it is possible to obtain this kind of information also for collisions involving highly vibrationally excited (U = 31) molecules.2 (2) Ziegler, G.;Kumar, S. V. K.; Rubahn, H.-G.; Kuhn, A.; Sun, B.; Bergmann, K . J . Chem. Phys. 1991, 94,4252.
0 1991 American Chemical Society
