level 8.0 英文手册

CP-663 April 2007 LEVEL 8.0 A Computer Program for Solving the Radial Schro¨dinger Equation for Bound and Quasibound Levels Robert J. Le Roy Guelph-Waterloo Centre for Graduate Work in Chemistry University of Waterloo, Waterloo, Ontario N2L 3G1, Canada Electronic mail: leroy@UWaterloo.ca University of Waterloo Chemical Physics Research Report LEVEL 8.0 A Computer Program for Solving the Radial Schro¨dinger Equation for Bound and Quasibound Levels Robert J. Le Roy Guelph-Waterloo Centre for Graduate Work in Chemistry University of Waterloo, Waterloo, Ontario N2L 3G1, Canada Electronic mail: leroy@UWaterloo.ca This manual describes program LEVEL, which can solve the radial or one-dimensional Schro¨dinger equation for bound and/or quasibound levels of any smooth single or double-minimum potential, and calculate inertial rotation and centrifugal distortion constants, expectation values, and/or Franck-Condon factors and other off-diagonal matrix elements, either for levels of a single potential or between levels of two different potentials. c©Robert J. Le Roy, 2007 1 Introduction Determining the number, energies and properties of vibration-rotational levels of a given one-dimensional or effective radial potential, and calculating matrix elements and transition intensities coupling levels of a single potential or levels of two separate potentials, are ubiquitous problems in chemical physics. The present report describes a robust and flexible computer program for performing such calculations. The original version of this program was based of the famous Franck-Condon intensity program of R.N. Zare [1, 2, 3], but the present version is considerably modified, and incorporates several unique features. In particular: (i) it will automatically locate and calculate the widths of quasibound (or orbiting resonance, or tunneling predissociation) levels; (ii) it can calculate diatomic molecule inertial rotation and centrifugal distortion constants for levels of a given potential; (iii) it can readily locate levels with dominant wave function amplitude over either well of an asymmetric double minimum potential; (iv) it can automatically locate and calculate expectation values for all vibration-rotation levels of any well-behaved single-minimum, “shelf state” or double minimum potential; (v) as an extension of (iv), it can automatically generate (for example) Franck-Condon factors and the radiative lifetimes for all possible discrete transitions allowed by specified (in the input data file) rotational selection rules between the levels of two different potentials, or among the levels of a single potential. While the present version generates the specified matrix elements and calculates Einstein A coefficients using the Ho¨nl-London factors for the case of singlet–singlet electronic transitions, it may be generalized to treat other cases. In the following, Section 2 presents the basic equation being solved, describes how the program func- tions, and outlines some of its options. Section 3 then states the input/output conventions, indicates the units assumed for the physical parameters of interest, and presents a shell to facilitate running the program on a UNIX or Linux system. The program’s operation is controlled by the contents of a data file which is read (on channel–5) during execution. The structure of this data file and the significance of the various read-in parameters are described in Section 4. Section 5 then describes the most significant differences between the current (7.8) and earlier [4] versions of this program. Finally, the Appendices outlines the structure of the program and the roles of its various subroutines, and present listings of illustrative sample data input files and of the resulting output. The current version of the (extensively commented) source code for LEVEL and a ‘pdf’ file of this manual may be obtained from the www site http://leroy.uwaterloo.ca/programs/ . While, there are no charges associated with distribution or use of this program, its use should be acknowledged in publications through reference to this report [5]. Users are also requested not to distribute the program themselves, but to refer other prospective users to the above web site or to the author. The version described herein includes corrections and enhancements incorporated up to 7 April 2007. Individuals currently utilizing older versions of this code [4] will likely find it desirable to obtain the current version since it has some corrections and additional functionality. I would also appreciate having users inform me of any apparent errors or instabilities in the code, or of additional features which might appear desirable for future versions. 1 2 Outline of Program Operation and Options 2.1 Solving the Radial Schro¨dinger Equation The core of the program is concerned with determination of the discrete eigenvalues and eigenfunctions of the radial or (effective) one-dimensional Schro¨dinger equation − ~ 2 2µ d2Ψv,J(r) dr2 + VJ(r) Ψv,J(r) = Ev,J Ψv,J(r) (1) in which µ is the effective or reduced mass of the system, J the rotational quantum number, and the effective one-dimensional potential VJ(r) is a sum of the rotationless (electronic) potential V (r) plus a centrifugal term. For the normal problem of a diatomic molecule rotating in three dimensions, this centrifugal potential has the form [J(J + 1) − Ω2] ~2/2µr2 , where Ω = OMEGA is the projection of the electronic angular momentum onto the internuclear axis. However, for the special case of a diatom rotating in two dimensions, a case invoked by setting the read-in parameter OMEGA > 99 (see the discussion of data input statement #5 in § 4), this term becomes [J2−1/4] ~2/2µr2 . The program also defines the reduced mass appearing in Eq. (1) as Watson’s “charge-modified reduced mass” [6], µ= µW = (MAMB)/(MA +MB −meQ) , where MA and MB are the atomic masses of the two atoms, me is the electron mass, and Q= CHARGE (see input Read statement #1) is the ± integer net charge on the molecule (ion). The core of the calculation is the solution of Eq. (1) to determine the eigenvalues Ev,J and eigen- functions Ψv,J(r) of the potential VJ(r). This is done in subroutine SCHRQ, which is based on the famous Cooley-Cashion-Zare routines SCHR [1, 2, 3, 7, 8], but incorporates special features, such as the ability to automatically locate and calculate the widths of “quasibound” or tunneling-predissociation lev- els [9, 10, 11, 12]. These are metastable levels which lie above the dissociation limit, but whose dissociation is inhibited by a potential energy barrier. The accuracy of the eigenvalues and eigenfunctions obtained is largely determined by the size of the (fixed) radial mesh RH (Read #4 of the input data file) used in the numerical integration of Eq. (1). For potentials that are not too steep or too sharply curved, adequate accuracy is usually obtained using an RH value which yields a minimum of 15 to 30 mesh points between adjacent wavefunction nodes in the classically allowed region. An appropriate mesh size may be estimated using the particle-in-a-box expression RH = pi/ ( NPN× [µ×max{E − V (r)}/16.857 629 20]1/2 ) (2) in which NPN is the selected minimum number of mesh points per wavefunction node (say 25), max{E − V (r)} is the maximum of the local kinetic energy (in cm−1) for the levels under consideration (in general it is . the potential well depth), and the numerical factor is identified below in Section 3. A value of NPN which is too small yields results which are unreliable, while too large a value may require excessive computational effort or cause array dimensions to be exceeded. Thus, while Eq. (2) is a useful guide, a careful user should always try different RH values in order to ensure that the results calculated achieve the accuracy desired for their particular application. The numerical integration of Eq. (1) is performed on the range from RMIN to RMAX (see Read statement #4) using the Numerov algorithm [7, 13]. To initiate this integration, it is necessary to specify initial values of the wave function at two adjacent mesh points at each end of the range. For truly bound states, the wave function at the outer end of the range RMAX is initialized at an arbitrary value (say, unity), while its value at the adjacent inner mesh point is determined using the first-order semiclassical or WKB wavefunction [14] Ψv,J(r) ∝ [VJ(r)− Ev,J ]−1/4 exp ( − √ 2µ/~2 ∫ r [VJ(r ′)− Ev,J ]1/2 dr′ ) (3) At short range, most realistic intermolecular potentials grow very steeply, causing the wavefunctions to die off extremely rapidly with decreasing r . As a result, the wave function at the inner end of the range of 2 integration is normally initialized by placing a node at the lower bound of this range, the read-in distance RMIN. This is effected by setting Ψv,J(r = RMIN) = 0 and giving Ψv,J(r = RMIN + RH) an arbitrary (non- zero) value. This is the normal case for a diatomic molecule problem. Note that one should normally set RMIN>0 , as the centrifugal contribution to the potential becomes singular at r= 0 . A special treatment of the inner boundary condition may be implemented if one is searching for eigen- functions of a symmetric potential whose midpoint is located at RMIN. For asymmetric levels which would have a node at RMIN, the normal treatment described above will suffice. However, another approach must be implemented for symmetric levels of such a potential whose wave functions would have zero slope at RMIN. This option is built into subroutine SCHRQ, and is invoked by setting the control parameter INNOD1≤ 0 . However, since this is an unusual case, varying this parameter is not one of the normal options of the current version of the main program, and a user who wishes to deal with this case may chose to either add parameter INNOD1 to Read #17 (line #367 of the MAIN program), or to recompile the code with the value of INNOD1 and INNOD2 defined on lines #365 & 366 of the code reset to 0. If desired, a hard wall outer boundary condition may be imposed by setting the input integer IV(i) (see Read #18), which would otherwise represent the vibrational quantum number, equal to a large negative number which specifies the radial mesh point where this wall would be placed. In particular, setting IV(i) < −10 causes a hard wall (wave function node) to be placed at mesh point number |IV(i)| for level–i. In the Cooley procedure for finding the eigenvalues of Eq. (1) [7, 8], for any given trial energy the numerical integration proceeds inward from RMAX and outward from RMIN until the two solution segments meet at a chosen matching point rx . The discontinuity in their slopes at rx is then used to estimate the energy correction required to converge on the eigenvalue closest to the given trial energy [15], and this process is repeated until the energy improvement is smaller than the chosen convergence criterion (parameter EPS of Read #4). This procedure usually converges very rapidly, and for a single-minimum potential it is insensitive to the choice of the matching point rx , as long as it lies in the classically-allowed region where the wavefunction amplitude is relatively large. However, to ensure high accuracy of calculated expectation values or matrix elements, EPS should usually be set 2 orders of magnitude smaller than the the actual eigenvalue precision required. For an asymmetric double-well potential, wavefunctions usually have amplitudes of very different mag- nitude over the two wells, and the eigenvalue correction algorithm [15] used by SCHRQ tends to become unstable if the matching distance rx lies in the well where the wavefunction has very small amplitude. As a result, it is usually necessary to require rx to lie in the well where its amplitude is the largest. In the current version of the program, this choice is set by the internal control parameter INNER ( = SINNER), which in turn is set to the appropriate value by the automatic vibrational level-finder subroutine ALF. As a result, calculations involving vibrational levels of a double well or “shelf-state” potential are (usually) performed just as routinely (for the user) as those for a normal single-well potential. In general, the outward and inward numerical integration must start at distances RMIN and RMAX (input via Read #4), respectively, which lie sufficiently far into the classically-forbidden regions (where VJ(r) > Ev,J ) that the wavefunction has decayed by several orders of magnitude relative to its amplitude in the classically-allowed region. The present version of the code prints warning messages if this decay is not by a factor of at least 10−9; if such warnings are printed, a smaller RMIN or larger RMAX value will be needed to ensure the desired accuracy for such cases. On the other hand, if RMIN or RMAX lie sufficiently far into the classically-forbidden regions that [VJ(r)− E] becomes extremely large, the integration algorithm can become numerically unstable for the given mesh size. For realistic diatomic molecule potential curves, this situation is only likely to occur near RMIN. If it does, a warning message is printed and the beginning of the integration range is automatically shifted outward until the problem disappears. However, use of a slightly larger value of RMIN will cause these warning messages to disappear and (marginally) reduce 3 the computational effort. For most diatomic molecules, a reasonable value of RMIN is ca. 0.6 − 0.8 times the smallest inner turning point encountered in the calculation, but for hydrides or other species of low reduced mass, even smaller values may be necessary. The program internally defines the upper bound on the range of numerical integration RMAX as the smaller of the value read in (Read #4) and the largest distance consistent with the specified mesh and the internally-defined (see § 3) potential and distance array dimension NDIMR. As with RMIN, the choice of RMAX is not critical, as long as (for truly bound states) the wave function has decayed to an amplitude much smaller than that in the classically allowed region, and the same amplitude decay test of 10−9 is used for it. However, due to the anharmonicity of typical molecular potential curves, the requisite values of RMAX are much larger for highly excited vibrational levels than for those lying near the potential minimum. In order to reduce computational effort, an integration range upper bound rend(v, J) is therefore determined for each level using the semiclassical result of Eq. (3), which shows that the wavefunction dies off exponentially in the classically forbidden region with an exponent of − √ 2µ/~2 ∫ rend(v,J) r2(v,J) [VJ(r)− Ev,J ]1/2 dr (4) where the turning point r2(v, J) marks the outer end of the classically accessible region at this energy Ev,J . For each level it considers, SCHRQ first locates r2(v, J), and then determines a value of rend(v, J) which is sufficiently large to ensure that this starting amplitude is smaller than that in the classically-allowed region by a factor of at least 10−9. In calculations for levels spanning a wide range of energies, the program’s use of this procedure can reduce the overall computation time by a factor of two or more. 2.2 Locating Quasibound Levels and Determining Their Widths Quasibound or orbiting resonance levels are metastable eigenstates of Eq. (1) which lie at energies above the potential asymptote VJ(r =∞) , but below a maximum in the outer part of the potential. As discussed in Refs. [9] & [10], the most efficient way of treating such levels is to determine their energies by applying an Airy function boundary condition at the third (outermost) turning point, and to calculate their widths using a uniform semiclassical method. The former appears to be the most accurate and efficient bound-state type method of locating quasibound or tunneling predissociation levels proposed to date [9, 10, 11, 16, 17]. It is virtually exact for narrow (long-lived) states, while for the very broadest levels lying marginally below barrier maxima, its differences with any competitive methods is at most a small fraction (. 0.2) of the level width (FWHM). More accurate predictions for such short-lived states would require a detailed simulation of the actual process by which they are observed, since different methods of observing a given quasibound level may yield apparent energies differing by a small fraction of the level width. For example, if such levels are being observed spectroscopically, the peaks in the actual bound→ continuum spectrum would be calculated using a photodissociation simulation code [18, 19]. In the present program, the width calculation is based on Eq. (4.5) of Connor and Smith [11]; a more transparent description of this procedure may be found in §II.B of Ref. [12]. This is a uniform semiclassical procedure in which the predissociation rate may be thought of as being the product of the probability of tunneling past the barrier at the specified energy, times the vibrational frequency (inverse of the vibrational period) for the system trapped in the well behind the barrier. The actual calculation requires the evaluation of an integral of the type seen in Eq. (4) across the interval between the two outermost classical turning points (i.e., with the upper bound rend(v, J) replaced by the outermost turning point r3(v, J)), and of an analogous integral across the classically allowed interval between the two inner turning points:√ 2µ/~2 ∫ r2(v,J) r1(v,J) [Ev,J − VJ(r)]−1/2 dr (5) 4 together with (the energy derivative of) a phase correction factor which takes account of the proximity to the barrier maximum [10, 11]. This procedure yields widths which are expected to be very reliable, particularly for narrow (long-lived) levels, but may have uncertainties of up to ca. 10% or more for the very broadest levels [10, 11]. To obtain more accurate results for such levels would again require one to perform a direct simulations of the process by which they are observed. On a practical note, if the end of the numerical integration range RMAX is smaller than the outermost turning point r3(v, J) of the metastable level of interest, the program attempts to generate a reasonable estimate of the width by completing the quadrature analytically while approximating the potential on the remainder of the interval by a centrifugal-type term C2/r 2 attached to the potential function at RMAX. If this approximation is invoked, warning messages are written to the main channel–6 output file (e.g., see the output for Case 3 in Appendix G). 2.3 Calculating Diatomic Molecule Centrifugal Distortion Constants The rotational sublevels of a given vibrational level of a molecule are conventionally represented by the power series Ev,J = G(v) + Bv[J(J + 1)] − Dv[J(J + 1)]2 + Hv[J(J + 1)]3 + ... (6) If desired, the program will calculate values of the inertial rotation constant Bv = (~2/2µ)〈v, J |1/r2|v, J〉 and of the first six centrifugal distortion constants associated with this expansion (Dv, Hv, Lv, Mv, Nv and Ov). These constants have their normal significance for rotationless (J = 0) vibrational levels [20], and are simply related to derivatives of the energy with respect to [J(J + 1)] when calculated for vibration- rotation levels with J > 0 . Calculation of these constants is invoked by setting input parameter LCDC > 0 (see Read #17), and are performed using a subroutine based on Tellinghuisen’s reformulation [21] of the exact quantum mechanical method of Hutson [22], which has been extended to higher order to allow the calculation of Nv and Ov. To ensure stable, fully convered calculations, it is often necessary to make the eigenvalue convergence parameter EPS of Read #4 quite small (e.g.,