Connecting the Dunham Expansion to the Dissociation Limit for

Subscriber access provided by University of Sunderland

Article

Connecting the Dunham Expansion to the Dissociation Limit for Interatomic Potentials: Application to Lennard-Jones m-n Potentials Randall B. Shirts J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.8b08095 • Publication Date (Web): 08 Oct 2018 Downloaded from http://pubs.acs.org on October 14, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Connecting the Dunham Expansion to the Dissociation Limit for Interatomic Potentials: Application to Lennard-Jones m-n Potentials Randall B. Shirts* Department of Chemistry and Biochemistry Brigham Young University Provo, Utah 84602 [email protected] Abstract Dunham generated the expansion for energy levels of a rotating, vibrating diatomic molecule from an expansion of the potential about the equilibrium position. For partition functions, however, the energy levels are needed all the way to dissociation. Analytic Morse oscillator energies are not very useful because the exponential decay of the Morse potential is much too short-ranged for any physical system. The longer-range Lennard-Jones 12-6 potential could be used, but quantum energies have not previously been conveniently fit. I show how Dunham coefficients begin a set of asymptotic functions for any interaction potential, one function arising from each successive term in the WKB expansion. I apply this to the family of Lennard-Jones m-n (LJ m-n) potentials with R-m repulsive term and R-n attractive term (m > n) and demonstrate how m can be used as a parameter to adjust either the equilibrium distance or harmonic frequency. I present an empirical parameterization of LJ m-n vibrotor energies starting with Dunham coefficients generated from four terms in the WKB expansion. This information is combined with data from numerically solved energies and asymptotic limits to fit the functions all the way to dissociation. One can also treat exp-6 and similar model potentials with different repulsive parts using the same method because the expansion form is controlled by the longrange part of the potential.

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 29

2 I. Introduction Molecular vibration-rotation state energies as a function of quantum number are required for the calculation of gas-phase partition functions needed for statistical calculation of thermodynamic quantities including equilibrium constants between reactive species. 1 Such information is essential for modeling studies of thermochemistry, air pollution, and climate change, and model potentials can play an important role when high quality potential energy surfaces are not available. Molecular vibrations were first approximated by independent harmonic normal modes; 2 however, all chemical bonds are anharmonic and have finite dissociation energy. The Morse oscillator is the usual anharmonic model potential for stretching vibrations because the harmonic frequency, bond length, and dissociation energy are all convenient parameters of the model. In addition, Morse oscillator quantum energies and wave functions are known analytically. 3,4 However, the exponential decay of the Morse potential is much too short-ranged for any physical system. The weak dispersion potential between nonpolar molecules is known to have a long-range decay proportional to R-6, and the dipole-dipole potential between polar molecules decays as R-3 at large distance. 5 The Lennard-Jones 12-6 6 potential was formulated to have the correct long-range behavior for nonpolar species and yet be computationally efficient because the repulsive term is proportional to the square of the attractive term. Even though the LennardJones 12-6 potential has a long history, quantum energies have not previously been conveniently fit. Other inverse power potentials with longer range could also be convenient model potentials, but the quantum energies of such potentials have not been conveniently fit either. I seek a way to accurately fit energy levels of long-range intermolecular potentials with an analytic form.


Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3 Dunham 7 showed how to generate the quantum energies of a rotating-vibrating diatomic molecule by expanding the interatomic potential energy function about the equilibrium geometry using an expansion in the variables 𝑣𝑣 + 12 and 𝐽𝐽(𝐽𝐽 + 1) : 𝑙𝑙

1 ∞ 𝑗𝑗 𝐸𝐸𝑙𝑙,𝑗𝑗 = ∑∞ 𝑗𝑗=0 ∑𝑙𝑙=0 𝑌𝑌𝑙𝑙,𝑗𝑗 �𝑣𝑣 + 2� [𝐽𝐽(𝐽𝐽 + 1)]

(1)

where Yl,j are the Dunham coefficients, which are algebraic functions of the potential expansion parameters and are expansions in the variable Be/νe, where Be is the rotational constant and νe is the harmonic vibrational frequency. Dunham used the first two terms in the WKB asymptotic expansion. It is now plausible to incorporate four or more terms in the WKB series. Electronic energy levels must be calculated using quantum calculations. However, semiclassical methods such as the WKB expansion, Born-Sommerfeld quantization, or the RKR inversion, have proven useful for rotational and vibrational energy levels of molecules because nuclei are considerably heavier than electrons and move much more slowly. Their motion is mostly classical. Therefore, the quantum energies associated with rotation and vibration are much closer together and quite accurately quantized using the classical action variables. My early career was spent exploring these methods, and this paper returns to semiclassical methods for describing diatomic vib-rotors and the stretching vibrations of weakly bonded systems. Another purpose of this article is to demonstrate the usefulness of Lennard-Jones m-n potentials (LJ m-n), where m = 12 and n = 6 is only one member of this larger family. In particular, the LJ m-3 potentials are of interest for understanding the weakly bound states of dipole-dipole complexes, including hydrogen-bonded species such as the water dimer. Other



Page 4 of 29

4 5

model potentials such as the exp-6 and MSV (Morse-Spline-van der Waals)

8,9

potentials can

be treated similarly. The LJ m-n potential is a two-parameter potential energy function given by 𝑉𝑉(𝑅𝑅) =

𝐷𝐷𝑒𝑒

𝑚𝑚−𝑛𝑛

𝑅𝑅

𝑚𝑚

𝑅𝑅

𝑛𝑛

�𝑛𝑛 � 𝑒𝑒� − 𝑚𝑚 � 𝑒𝑒� � , 𝑅𝑅

𝑅𝑅

(2)

where Re is the equilibrium distance at which the potential has a minimum value of –De. The potential approaches zero for large R values, so De is the classical dissociation energy. To be physical meaningful, m > n and n > 0, but neither m nor n need be an integer. The harmonic frequency for low energy oscillation is controlled by V″(Re) = mnDe/Re2. For simplicity, I choose reduced time units so the harmonic frequency at the bottom of the well is ω = (V″(Re)/µ)1/2 = 1, where µ is the reduced mass parameter. I also choose reduced energy units so ω = 1 and reduced length units so V″(Re) = 1 (Re = (mnDe)1/2), where D is hereafter measured in multiples of ω (D = De/ω). The form of LJ m-n potentials is shown in Fig. 1. Figure 1(a) shows plots of LJ 2n-n potentials in reduced units for n = 2.5, 3, 4, 5, 6, and 7 for D = 10. Note that they all have the desired form with a repulsive wall at small R values, a minimum with identical dissociation energy and harmonic force constant, and all decay to zero at large distance. They differ, however, in their equilibrium distance, Re, and this is required to equalize the harmonic force constant because, for given n and m, the potential depends on only two parameters, Re and De. The Morse potential, on the other hand, has three parameters, and Re can be independently specified with a given well depth and harmonic force constant. Figure 1(b) shows plots of LJ m6 potentials with several values of m. It is difficult to see that higher m values have steeper repulsive walls, but the increase of Re with m is obvious. On the other hand, Figure 1(c) shows


Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5 LJ 8-6 and LJ 12-4 potentials with D = 10. These potentials have the same Re, and it is easier to see that the LJ 12-4 potential has longer range for large R and a steeper repulsive wall at small R, even though the harmonic force constant is identical for the two potentials. Thus, it is useful to consider the LJ m-n family of potentials to be a three-parameter family with n chosen to have the correct long-range behavior and m adjusted to give a correct value for Re (or alternatively m can be used to adjust ωe with prescribed Re). As long as only bound quantum states are of interest, the steepness of the repulsive wall dependent on the m parameter is of little concern except in high-energy collisions that sample the potential at smaller values of R.

Fig. 1. (a) Plots of LJ 2n-n potentials in reduced units for n = 2.5 (dot-dash), 3 (solid), 4 (dot), 5 (dash), 6 (gray solid), and 7 (gray dash) for D = 10. (b) Plots of LJ m-6 potentials in reduced units for m = 6.001 (dot-dash), 8 (solid), 10 (dot)), 12 (dash), 14 (gray solid), 16 (gray dot), 18 (gray dash) with D = 10. (c) Plots of LJ 8-6 (solid) and LJ 12-4 (dash) potentials, which have the same Re. Again D =10 in reduced units.



Page 6 of 29

6 In this paper, I give numerical fits for the quantum energies eigenvalues for pure vibration of LJ m-n potentials and for an exp-6 potential adjusted to have the same harmonic frequency as a LJ 12-6 potential. In addition, I present analytic formulas for vib-rotor Dunham coefficients of LJ m-n systems up through terms with v + j < 7. These are all the values accessible using four terms in the WKB series (ninth-order in ℏ). In these parameterizations, D

is a parameter, so the fits should be appropriate for any value of D, not just for those used in the fit. II. Methods 1. Vibrational Energies The leading term in the fit of quantum energy values of LJ m-n (and similar) potentials is given by the prediction of Bohr-Sommerfeld quantization as discussed in the supporting information. This is equivalent to using just one term in the WKB expansion. This procedure gives an approximation for the energies as a function of the vibrational quantum number that is proportional to D and that I represent as 𝐸𝐸𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 (𝑣𝑣) = −𝐷𝐷(1 − 𝑧𝑧)𝑝𝑝 𝑓𝑓1,0 (𝑧𝑧) ,

(3)

where p, the exponent controlling how fast eigenvalues approach zero, was determined by LeRoy et al. 10,11 to have the value p = 2n/(n – 2). Here, f1,0 (z) is a smooth function of z independent of D, and 𝑧𝑧 =

ℏ𝜔𝜔(𝑣𝑣+12) 𝐷𝐷𝑒𝑒 𝐺𝐺𝑚𝑚𝑚𝑚

,

(4)

where v is the vibrational quantum number and −1/2

𝐺𝐺𝑚𝑚𝑚𝑚 = (2𝜋𝜋)

−3/2

(𝑚𝑚 − 𝑛𝑛)

𝑛𝑛

𝑚𝑚−2𝑛𝑛+2 2(𝑚𝑚−𝑛𝑛)

𝑚𝑚

2𝑚𝑚−𝑛𝑛−2 2(𝑚𝑚−𝑛𝑛)

Γ�

Gmn has a simpler form for m = 2n given by


𝑛𝑛−2

Γ�2(𝑚𝑚−𝑛𝑛)�

3𝑚𝑚−2𝑛𝑛−2 � 2(𝑚𝑚−𝑛𝑛)

.

(5)

Page 7 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝐺𝐺𝑛𝑛 =

1 1− 2 𝑛𝑛

√𝜋𝜋

7 1 Γ�12−𝑛𝑛 � 1 Γ�2−𝑛𝑛�

.

(6)

The derivation of Gmn given in Eq. (5) is provided in the supporting information. Here Gmn determines the maximum number of bound energy states, and the continuous variable z varies from 0 at the bottom of the well to 1 at dissociation, but only for integer values of v do z values represent vibrational quantum states. Energy eigenvalues, then, are equally spaced in the variable z with the ground state located at half that spacing, and the last bound state is near z = 1. The maximum value of v is given by rounding down to an integer the expression 𝑣𝑣𝑚𝑚𝑚𝑚𝑚𝑚 =

𝐺𝐺𝑚𝑚𝑚𝑚 𝐷𝐷𝑒𝑒 ℏ𝜔𝜔

− 12 ,

(7)

Thus, when D is given in units of ω, Gmn represents a multiplier telling how many more states are bound compared to the states a harmonic oscillator has below energy D. For example when m = 2n, in the limit n → ∞, Gn → 2, which is same behavior as for the exponentially decaying Morse oscillator. In the limit n → 2, Gn → ∞, meaning that there are an infinite number of bound states for potentials that fall off as slow or slower than R–2. Other values of physical interest are G6 = 2.86296… for LJ 12-6, and G3 = 5.52223… for LJ 6-3. A quantity like Gmn giving the number of bound states also exists for potentials not of the LJ m-n form, such as exp6, and may be determined numerically as described in the supporting information. Figure 2 shows a plot of Gn-1 and p-1 as a function of n-1. In this plot, p-1 is linear. Both functions go to zero at n–1 = 1/2 (n = 2) showing that the R–2 has an infinite number of bound states and that the eigenvalues consequently decay too slowly for Gn to be finite. Both functions approach 1/2 as n– 1

→ 0 (Gn → 2 for n → ∞). Again, this is the result for the Morse oscillator. Figure 2 shows the

close relationship of Gn and p. Note, however, that the value of p depends only on n, the power of the attractive term. In fact, the value of Gmn depends only weakly on m. Figure 3 shows how



Page 8 of 29

8 Gmn varies with m for n = 3, 4, 5, 6, and 8. Gm3 is greater than other plotted values of Gmn, which means the attractive R-3 potential, because of its longer range, binds significantly more states than other attractive potentials with the same value of D. Each curve in Fig. 3 has a minimum at a value of m given by m = n + 2 + ∆, where ∆ increases slowly from zero to 0.53 as m increases from 2 to 10. Decreasing m from the minimum, Gmn increases to the limiting value lim 𝐺𝐺𝑛𝑛+𝜀𝜀,𝑛𝑛 = � 𝜀𝜀→0

1 1

− 2𝑛𝑛 3/2 𝑒𝑒 2 𝑛𝑛

�

𝑛𝑛+2

√2𝜋𝜋

.

(8)

The plots of Gmn for higher values of n in Fig. 3 are flatter and approach 2.

Fig. 2. Plot of Gn–1 (solid) and p–1 (dash) versus n–1 showing the close relationship between Gn and p.

Fig. 3. Plot of Gmn versus m for n = 3 (solid), 4 (dot), 5 (dash), 6 (gray solid) and 8 (gray dotted).


Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


9 To achieve greater accuracy than Eq. (3), one may include more terms in the WKB asymptotic expansion. An expression including four KWB terms is given by Kesarwani and Varshni. 12 Each term in the WKB series allows one to add an additional term to Eq. (3), arriving at the final expression given below as Eq. (10). The WKB integrals needed for higher order terms are contour integrals in the complex R plane with integrands containing various derivatives of V(R) and powers of the quantity [E − V(R)]k/2, where k is an odd integer.12 Dunham’s method7 involves expanding the potential energy function about the equilibrium geometry in the following form: 𝑉𝑉(𝑅𝑅) = 𝑎𝑎0 𝜉𝜉 2 (1 + 𝑎𝑎1 𝜉𝜉 + 𝑎𝑎2 𝜉𝜉 2 + 𝑎𝑎3 𝜉𝜉 3 + ⋯ )

(9)

where 𝜉𝜉 = (𝑅𝑅 − 𝑅𝑅𝑒𝑒 )/𝑅𝑅𝑒𝑒 and the coefficients ai are algebraic functions of the potential

parameters. Dunham then evaluates the WKB contour integrals by the method of residues. The contours of the integral circle the zeros of E − V(R) in the complex R plane in the counterclockwise direction. Dunham transformed the integrals to the complex V plane, in which case the contour circles the zeros of E − V(R) twice. Dunham expanded the integrands in powers of V1/2 and expressed the integrals in terms of the coefficients (residues) of V-1 in the expansion. This method allows evaluation of the integrals as an expansion in powers of E. This expansion can be inverted to obtain E as an expansion in powers of (𝜐𝜐 + 12) involving algebraic expressions in terms of the ai coefficients, which gives the Dunham coefficients for J = 0 directly or can be

algebraically converted to give the expansion coefficients for f1,0 (z) in Eq. (3) and for a few higher-order corrections included in Eq. (10) below. Dunham used two terms in the WKB expansion, and extension to include the four terms given by Kesarwani and Varshni12 is straightforward though algebraically demanding.



Page 10 of 29

10 For E approaching dissociation, however, the expansion coefficients are not easily extrapolated, and evaluating the WKB contour integrals numerically near dissociation is computationally unstable. Numerical evaluation of the integrals is complicated by singularities at the turning points. Kirschner and LeRoy 13 show how to avoid these singularities using energy derivatives of nonsingular integrals. However, highly accurate numerical differentiation is also unstable. I found it easier just to solve the Schrödinger equation numerically and fit the computed eigenvalues to the functional form below. I used the Numerov method as described in the supporting information to numerically determine several thousand eigenvalues and fit them using the form given below in Eq. (10). I fit eigenvalues for n = 2.5 with m = 5, n = 3 with m = {3.5,4,5,6,7,8,9,10,11,12}, n = 4 with m = {4.5,5,6,7,8,9,10,11,12} n = 5 with m = {5.5,6,7,8,9,10,11,12}, n = 6 with m = {6.5,7,8,9,10,11,12,13,14,15}, n = 7 with m = {7.5,8,9,10,11,12,13,14,15}, n = 8 with m = {8.5,9,10,11,12,13,14,15,16,17,18}, and n = 10 with m = {10.5,11,12,13,14,15,16,17,18,19,20}, and also for an exp-6 potential fit to a LJ 12-6. Some results are shown below, but mostly in the supporting information. Including corrections found empirically from numerically determined eigenvalues, vibrational energy levels of potentials of the LJ m-n form were fit with the following expansion:

𝐶𝐶4,0 𝐷𝐷5

𝐸𝐸(𝑣𝑣) = −𝐷𝐷(1 − 𝑧𝑧)𝑝𝑝 𝑓𝑓1,0 (𝑧𝑧) +

(1 − 𝑧𝑧)𝑝𝑝−6 𝑓𝑓4,0 (𝑧𝑧) ⋯ ,

𝐶𝐶2,0 𝐷𝐷

(1 − 𝑧𝑧)𝑝𝑝−2 𝑓𝑓2,0 (𝑧𝑧) +

𝐶𝐶3,0 𝐷𝐷3

(1 − 𝑧𝑧)𝑝𝑝−4 𝑓𝑓3,0 (𝑧𝑧) +

(10)

where, in each successive term, the power of (1 – z) decreases by 2, and the power of D in the denominator increases by 2. In the chosen reduced units, D-1 = 2mnBe , where Be is in units of hve; thus an expansion in D-1 is equivalent to one in Be/ve as expressed by Dunham. As the first term in Eq. (10) comes from the Born-Sommerfeld quantization of the classical action at half-


Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


11 integer values, each additional term in Eq. (10) arises from the subsequent term in the WKB expansion (with earlier terms convolved). The first term in the WKB expansion is equivalent to Born-Sommerfeld quantization. The second subscript zero in Eq. (10) indicates J = 0, and nonzero rotations will be included later. When no comma nor second subscript is included in the Ci,J values, fi,J(z) functions or their polynomial coefficients, J = 0 is assumed. The initial parameters in this expansion can be determined by Dunham’s7 method. Using the four terms in the expansion given by Kasarwani and Varshni12 allows the exact determination of Cj cofficients through j = 4 and expansion coefficients up through fij where j = 9 – 2i. For example, terms up through f17 and f41 can be determined exactly. Further terms extrapolating to dissociation can be determined by iteratively fitting the numerically determined eigenvalues. First, using calculated values for C2, C3, C4, calculated expansions for f2, f3, and f4, and neglecting higher order terms not explicitly displayed in Eq. (10), one numerically solves for E(v), v = 0, 1, 2, . . .vmax for several values of D and plots the following versus z. 𝑓𝑓1𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 (𝑧𝑧) =

𝐶𝐶 𝐶𝐶 𝐶𝐶 −𝐸𝐸(𝑣𝑣)+ 𝐷𝐷2 (1−𝑧𝑧)𝑝𝑝−2 𝑓𝑓2 (𝑧𝑧)+ 33 (1−𝑧𝑧)𝑝𝑝−4 𝑓𝑓3 (𝑧𝑧)+ 45 (1−𝑧𝑧)𝑝𝑝−6 𝑓𝑓4 (𝑧𝑧) 𝐷𝐷

𝐷𝐷(1−𝑧𝑧)𝑝𝑝

𝐷𝐷

.

(10)

The largest errors in f1 when doing this are proportional to D-2 and (1 – z)-2. Typically, f1(z) is well determined for most of the range (0,1) except near z = 1, and one can extrapolate f1(1) and f1′(1) from the behavior of f1approx for the largest value of D used. The appearance of such a plot is shown in Fig. 4 with m = 12, n = 6 and D = 4, 16, and 64 with f2 and f3 functions set to zero. This figure shows how the f1(z) values computed for each quantum state fall along the same curve independent of D except for the neglect of f2 corrections, which are proportional to D-2 and (1 – z)-2. The curve for D = 4 has 11 equally spaced points (v = 0 through 10) as predicted by Eq. (4) with G6 = 2.86296. . . The curves for D = 16 and D = 64 likewise have 46 and 183 equally



Page 12 of 29

12 spaced points respectively. It easy to extrapolate from Fig. 4 that f1(1) is close to 1.10 with zero derivative. After subsequent determination of and correction for f2, f3, and f4, all points will fall on the same curve to numerical accuracy. The numerical value for f1(1) can be determined by the method of Paulsson et al., 14 and the empirically fit value agrees with the analytic value to several decimal points. Agreement between empirically fit values of fk(1) for k = 2, 3, 4 and values determined using the methods of Ref. 14 becomes poorer as k increases, and the source of this disagreement may be either the approximation used in Ref. 14 or uncertainties in the numerical data. Analytic formulas for fk(1) for k = 1, 2, 3, 4 using Ref. 14 are contained in the supporting information.

Fig. 4. Plot showing initial determination of f1(z) for the LJ 12-6 potential using Eq. (11) with C2 = 0 and C3 = 0, etc. Values for f1approx for D = 4 (gray triangles), 16 (squares), and 64 (diamonds) are plotted. Note all points fall on the same curve for f1(z) on the left side. Deviations on the right are smaller for higher values of D, and subsequent determination of f2(z) and f3(z) correct this deviation, and points eventually fall on the same f1(z) line for the full range of z.

After initial determination of f1approx, one solves Eq. (10) for f2(z) and determines f2approx using the just-determined f1approx in Eq. (10). After determining f2(z), one then determines f3approx, etc. After initial determination of the fiapprox functions, one then returns and refines the approximation for f1, then f2, etc. Two iterations are usually sufficient.


Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


13 The form in Eq. (10) can only be asymptotic for large values of D because the powers of (1 – z) become negative after a few terms. For example, for n = 3, p = 6, and this means that the term containing f5 blows up as z → 1 and E → 0 near dissociation. The term containing f5, however, is too small to be determined without higher accuracy than obtained so far. On the other hand, for n = 6, p = 3, and the term containing f3 is singular as z → 1. However, these corrections are so small for larger values of D having many bound states that the singularity is not important, and for smaller values of D for which the third term is larger, fewer bound eigenvalues mean few states near the singularity at z = 1. If the form in Eq. (10) is assumed, one only need determine the scale factors, Cj, and the functions fj(z) to accurately fit energy eigenvalues for arbitrary values of D. To do this, I have adopted for convenience the convention that fj(0) = 1, and thus the magnitude and sign of each term is determined by the Cj coefficients, whose magnitudes are typically much less than one, indicating that the expansion converges quickly. The Cj scale factors up through j = 4, and the first few polynomial expansion coefficients fij were determined using Dunham’s methods. These are algebraic functions of the ai coefficients, and thus m and n. For example, C2 = (7 a12 – 12 a2)/(32a0) = [2m2 + m(12 − 5n) + 2n2 + 12n + 36]/(576nm), which simplifies to C2 = (n + 1)/32n2 for m = 2n. The correction term proportional to C2 may be considered a quantum correction for tunneling at the outer turning point, which tunneling is more prominent for lower values of n because the potential has longer range. The potential range also explains why the value of the parameter p is independent of m, as it describes energy levels near dissociation, and wave functions have the largest magnitude near the outer turning point where the repulsive potential term is negligible. Using this normalization, the form of fi,J(z) can easily be taken to have polynomial form.



2

3

𝑓𝑓𝑖𝑖,𝐽𝐽 = 1 + 𝑓𝑓𝑖𝑖1,𝐽𝐽 𝑧𝑧 + 𝑓𝑓𝑖𝑖2,𝐽𝐽 𝑧𝑧 + 𝑓𝑓𝑖𝑖3,𝐽𝐽 𝑧𝑧 + ⋯ .

Page 14 of 29

14 (11)

The fij,J coefficients are fitting parameters; however, f11,0 = Gmn – p is required to maintain the condition that the harmonic frequency at the bottom of the well is unity in reduced units. Note that f11,0 is quite small because of the approximate equality of Gmn and p (see Fig. 2). The values of f12,0 and f13,0, together with the value of p determine first and second anharmonicities respectively, etc. In most cases, f12,0 and f13,0 are even smaller than f11,0, indicating that, in addition to describing the high energy approach to E = 0, most of the vibrational anharmonicity is described by the term −D(1 − z)p. However, the functions fj(z) can be interpreted to be functions which interpolate from behavior at the bottom of the well (z = 0) to behavior at dissociation (z = 1). For n = 3 and m > 8, I found that C3 becomes very small (C3 is zero at approximately m = 9.8232, thus making f32>>1 and f33>>1. In this case, one could drop C3 and solve for f3 directly with a variable f30 placed in Eq. (11), if desired. Longer range LJ m-n potentials with n = 3 and n = 4, which have many states near dissociation, the functions f1 and to some extent f2 are sufficiently flat as z → 1 that it is easier to expand around z = 1, giving the following form 𝑓𝑓𝑖𝑖 = 𝑓𝑓̃𝑖𝑖0 + (1 − 𝑧𝑧)𝑘𝑘 �𝑓𝑓̃𝑖𝑖𝑖𝑖 + (1 − 𝑧𝑧)𝑓𝑓̃𝑖𝑖,𝑘𝑘+1 + (1 − 𝑧𝑧)2 𝑓𝑓̃𝑖𝑖,𝑘𝑘+2 + ⋯� ,

(12)

where k – 1 is the number of zero derivatives at z = 1. When this parameterization is used, algebraic constraints are applied so that

and

𝑓𝑓𝑖𝑖 (0) = 1 = 𝑓𝑓̃𝑖𝑖0 + 𝑓𝑓̃𝑖𝑖𝑖𝑖 + 𝑓𝑓̃𝑖𝑖,𝑘𝑘+1 + 𝑓𝑓̃𝑖𝑖,𝑘𝑘+2 + ⋯ ,

(13)

𝑓𝑓1′ (0) = 𝐺𝐺𝑚𝑚𝑚𝑚 – 𝑝𝑝 = −𝑘𝑘𝑓𝑓̃1𝑘𝑘 − (𝑘𝑘 + 1)𝑓𝑓̃1,𝑘𝑘+1 − (𝑘𝑘 + 2)𝑓𝑓̃1,𝑘𝑘+2 ⋯ ,

(14)


Page 15 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


15 etc. to fix coefficients known from the Dunham expansion. A plot of fit f1 functions for LJ 2n-n potentials is shown in Fig. 5(a). Note that in all cases f1(z) is monotonically increasing, which means that anharmonicity lowers the energy of higher excited states more than predicted by the simple power law given by −D(1 − z)p. Also, the maximum value of f1(z) decreases with increasing n, indicating that corrections to the power law are greater for a longer range potential. However, for n < 5, the function is best fit by Eq. (12) (where k = 9 for n = 2.5, k = 5 for n = 3, and k = 3 for n = 4) (Theoretical considerations suggest that k = p 1). In this case where n is small, p is large and many states are near dissociation, and the same value of f1 corrects for increased binding over a larger range of vibrational quantum number even though it is a small energy range. The function f1 is easily calculated using Eq. (11) to an accuracy of 10-11. Thus, tabulation and interpolation can be more accurate than a polynomial fit.



Page 16 of 29

16 Fig. 5. Plots of fi(z) functions for LJ 2n-n potentials for n = 2.5 (dot-dash), 3 (solid), 4 (dot), 5 (dash) , 6 (gray solid), 8 (gray dot), and 10 gray dash). (a) f1(z) (b) f2(z) (c) f3(z) (d) f4(z) (approximate).

A plot of f2 functions of LJ 2n-n potentials is shown in Fig. 5(b). Again, f2(z) values decrease with increasing n. In fact, only functions for n = 2.5 and n = 3 differ significantly from unity. Determined f2(z) values are only accurate to an accuracy of about 10-6 due to 10-15 level round-off errors in the eigenvalues. Plots of f3 functions for LJ 2n-n potentials are shown in Fig. 5(c). Again, only functions for n = 2.5 (which limits to 9.6) and n = 3 (which limits to 3.0) differ significantly from unity. Determined f3(z) values are only accurate to about 0.001. Plots of f4 functions, shown in Fig. 5(d) are accurate only to about 50%. The value of f4(1) in Fig. 5(d) is about -12 for n = 2.5, but this may not be accurate due to uncertainties in estimating a small correction. Plots of f1, f2, f3, and f4 functions for other LJ m-n potentials are shown in Figures S2–S8 (supporting information). In the plots, it is of interest to note that the curve with the highest value of fi(1) is, with few exceptions, the one with m nearest the minimum of Gmn (usually m = n + 2). Also of note, when m = n + ½, the f1(z) curves typically bend down to lower values on the right, perhaps anticipating a negative singularity as m approaches n. For these reasons, I only recommend using m ≥ n + 1. This range of m guarantees that the fi(z) functions are well-behaved and with an accurate fit. Analytic formulas for C3 and C4 and tables of computed C3 (Table S1) and C4 (Table S2) values are included in the supporting information. Note also that all tabulated values of C3 are less than 10-4, indicating that the series in Eq. (10) converges quickly, especially for high values of D.


Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


17 -15

The 10

level of accuracy is sufficient to determine corrections only though f4.

Higher order terms appear to contribute for smaller values of D, but I have not been able to determine f5 to any degree of accuracy. In fact, the uncertainty in values of f4 is significant. Only f41 is determined analytically, and the function for z > 0.5 is, at best, a guess. I used the iterative method described earlier for fitting eigenvalues with the following comments: 1. The Numerov method used to calculate the eigenvalues of a given potential as described in the supporting information requires initial estimates of the E, Rmin and Rmax values for each eigenvalue integration. These values were then refined iteratively. I used D values of 2, 4, 8, 16, 32, 64, and 128 as explained in comment #3 below. Additional comments on application of the Numerov method are found in the supporting information. 2. I determined values of f1(z) using Eq. (11) and polynomial fits for f1(z) using data for D = 128, where errors from inaccurate values of f2(z), f3(z) and f4(z) near z = 1 are smallest. 3. I determined values for f2(z) using interpolated values of f1(z) from D = 128 because interpolation was more accurate than my polynomial fits. For D = 128, computed points are closer together for interpolation, and errors due to inaccurate f2 and f3 values are smaller. This meant that D = 64 and D = 32 data were used to determine f2(z). By choosing powers of 2 for D values, interpolated z values were halfway between z values for computed f1(z) values for all lower values of D. I used Lagrangian 4-point interpolation. 4. I determined f3(z) values using interpolated f1(z) data from D = 128 and fit f2(z) data from D = 64 and D = 32. The f3(z) fits used data from D ≤ 16 where f3(z) corrections were larger and errors due to inaccurate values of f1(z) and f2(z) are smaller. In determining f3(z), fitting



Page 18 of 29

18 2

errors in f2(z) are proportional to D , so it is useful to see the calculated f3approx(z) points for D = 16, 8, 4, and 2 to look for convergence to the correct value. 5. I determined f4(z) values using interpolated f1(z) data from D = 128, fit f2(z) data from D = 64 and D = 32, and fit f3(z) data from D = 4 and D = 8. The f4(z) fits used data from D = 2, but most data is insufficiently defined for accurate determination.

2. Using the Results Three different levels of approximation may be useful in applications. The first level of approximation uses the fact that, for most values of m and n, f1(z) and f2(z) are approximately one. Thus, if not too high accuracy is needed, the energy levels of the Lennard-Jones m-n potential can be approximated by 𝐸𝐸𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎1 (𝑛𝑛) = −𝐷𝐷(1 − 𝑧𝑧)𝑝𝑝 +

𝐶𝐶2 𝐷𝐷

(1 − 𝑧𝑧)𝑝𝑝−2 ,

(15)

where C2 is given by the analytic form given above Eq. (10). The maximum error in energy eigenvalues computed using this approximation is tabulated in Table S3, and is approximately D × 2.8 ×10-2 for m = n + 2 and decreases to below D × 1.0 ×10-2 for m >> n. The error is proportional to D because the error is dominated by approximating f1(z) by 1, and the term containing f1(z) is proportional to D. The error is maximum for m = n + 2 because f1(1) is maximum for m = n + 2, as seen in Figs. S1-S7. Errors in the approximated energies for m and n values not included in the data (for example, non-integer values of n) may be somewhat larger. If higher accuracy is needed, it may be sufficient to approximate the functions f1, f2 and f3 by a quartic polynomial that reproduces the values of fi1, fi2 and fi(1) and fi′(1). However, because the terms containing f2(z) and f3(z) are small, it is sufficient to approximate f1(z) by a quartic polynomial, giving


Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝐸𝐸𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎2 (𝑛𝑛) = −𝐷𝐷(1 −

(4) 𝑧𝑧)𝑝𝑝 𝑓𝑓1 (𝑧𝑧)

+

𝐶𝐶2 𝐷𝐷

19 (1 − 𝑧𝑧)

𝑝𝑝−2

,

(16)

Except for m = n + ½, where fitting is poor, maximum fitting errors using this approximation were D × 5 ×10-3 for m = 12 and n = 10, in which case a quartic polynomial cannot adequately follow a late rise in f1(z) near z = 1 (see Fig. S7(a)). More typical maximum errors using this level of approximation are D × 1 × 10-3 or less. The maximum fitting errors found for individual values of m and n are tabulated in Table S3. To determine the polynomial coefficients of f1(4)(z) using the supporting information, f10 = 1, f11 = Gmn – p, f13 = − f1′(1) + 4 (f1(1) − 1 ) −3 f11 − 2 f12, and f14 = f1(1) − 1 − f11 − f12 − f13. If maximum accuracy is desired, then one can use the tabulated f1(z) values and leastsquared polynomial approximations for the f2, f3 and f4 functions, all given in the supporting information. In most cases, one can achieve a relative error of 10-8 × D by using these parameterizations. Exceptions include the last bound states for small values of D. In this case, the terms containing f3(z) and f4(z) are singular near z = 1, and the asymptotic series diverges, especially so for higher value of n. For n = 3 and n = 4 with small values of D where the divergences of the terms containing f3(z) and f4(z) are less severe, the term containing f5(z) could give better accuracy, as Ck are larger and the series does not converge as quickly as for higher n values. However, except in rare cases where a state very near z = 1 for a particular value of D, the approximate energy values using fit functions of only the first four WKB terms are accurate enough for almost any purpose. 3. Exp-6 and other model potentials Quantum energies for potentials that have a repulsive term which is not of the LennardJones m-n form, but which have an attractive term dominated by R−n at large distances, may be



Page 20 of 29

20 fit with the same form given by Eq. (10). Potentials that have additional attractive terms of shorter range than a R-n term may also be treated. However, the value of G, the scaling factor for the number of bound states required in Eq. (4), is no longer given by Eq. (5). In this case, the value of G may be determined numerically as described in the supporting information. One or both of the methods described should work just as well if the attractive potential is a sum of two or more terms rather than a single power as in MSV or ESMSV potentials.8,9 The exp-n potential can be taken to be of the following form:5 𝑉𝑉𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 =

𝐷𝐷

𝛾𝛾−𝑛𝑛

�𝑛𝑛𝑒𝑒

𝑅𝑅

−𝛾𝛾(𝑅𝑅 −1) 𝑒𝑒

𝑅𝑅

𝑛𝑛

− 𝛾𝛾 � 𝑒𝑒� � . 𝑅𝑅

(17)

This model potential was developed because the repulsive term in the LJ m-n potential was considered too hard. However, the exp-n repulsive exponential term is finite at small distances, so it is always overwhelmed by the singularity in the attractive term as R → 0. Thus, the exp-n potential is nonphysical by having a maximum and becoming negative at small distances. To be useful, the maximum must be at sufficiently small R values and with Vmax sufficiently high that tunneling to small R can be ignored. To agree with a LJ m-n potential with V(Re) = −D, V′(Re) = 0 and V″(Re) = 1 at Re = (mnD)1/2, one must choose γ as the + root of 𝛾𝛾 2 − 𝛾𝛾(𝑛𝑛 + 𝑚𝑚 + 1) − 𝑛𝑛𝑛𝑛 = 0 .

(18)

To be physically reasonable, a good rule of thumb is m > n + 2 and n > 2.4. For example, with n = 6 and m = 12, then γ is approximately 13.772. For this value of γ, the value of G was determined to be 2.78807303, a value measurably lower than for the LJ 12-6 potential and in keeping with a softer repulsive wall (see Fig. 3), as expected. Values for the eigenvalue fitting parameters of this potential are given in Table S4 for this value of γ up through f2(z), where I used a polynomial of degree 8 in z for f1(z) and a polynomial of degree 5 for f2(z). Neither f3(z)


Page 21 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


21 nor f4(z) was determined because I determined the eigenvalues for only a limited number of D values. However, this example is provided to illustrate the utility, in principle, of using Eq. (10) to fit eigenvalues, without an exhaustive survey of its applications.

4. Rotational Terms The rotational energies and centrifugal terms are determined by adding the centrifugal potential to the Lennard-Jones m-n potential: 𝑉𝑉𝑒𝑒𝑒𝑒𝑒𝑒 (𝑅𝑅) =

𝐷𝐷

𝑚𝑚−𝑛𝑛

𝑅𝑅

𝑚𝑚

𝑅𝑅

𝑛𝑛

�𝑛𝑛 � 𝑒𝑒� − 𝑚𝑚 � 𝑒𝑒� � + 𝑅𝑅

𝑅𝑅

ℏ2 𝐽𝐽(𝐽𝐽+1) 2𝜇𝜇𝑅𝑅2

.

(19)

Using Dunham’s method,7 one determines a perturbed equilibrium position in powers of J(J+1) and then calculates perturbed ak coefficients in powers of J(J+1). I determined Dunham coefficients to the maximum order possible using four WKB terms for Yv,J expressions for v+J < 7. Analytic expressions for these are given in the supporting information. Then I converted the Dunham coefficients to a form suggested by Eq. (10) and given as Eq. (20) below. So far, I have only determined the rotational energies only for m = 2n for n = {2.5, 3, 4, 5, 6, 8}, and additionally for n = 3 with m = {5,7,8}, and for n = 6, with m = 10,11,13. For each one, I determined 60–65 energy eigenvalues using the Numerov method as a function of C = 2J(J+1)/2µ, and I fit the eigenvalues to obtain the expansion coefficients up through terms proportional to C7. The values of these coefficients at the bottom of the well can be determined by expanding the potential about the minimum of Veff(R) using Dunham’s method, but the vdependence limiting to z = 1 was determined by empirical fitting.

and

𝐸𝐸(𝐶𝐶, 𝑣𝑣) = 𝐸𝐸0 (𝑣𝑣) + 𝑏𝑏1,𝑣𝑣 𝐽𝐽(𝐽𝐽 + 1) − 𝑏𝑏2,𝑣𝑣 𝐽𝐽2 (𝐽𝐽 + 1)2 − 𝑏𝑏3,𝑣𝑣 𝐽𝐽3 (𝐽𝐽 + 1)3 − ⋯ ,


(20)


(21)

𝑏𝑏𝐽𝐽,𝑣𝑣 =

𝐶𝐶1,𝐽𝐽

𝐷𝐷2𝐽𝐽−1

[𝑓𝑓1,𝐽𝐽 (𝑧𝑧)(1 − 𝑧𝑧)𝑝𝑝−2𝐽𝐽 +

𝐶𝐶2,𝐽𝐽 𝑓𝑓2,𝐽𝐽 (𝑧𝑧)(1−𝑧𝑧)𝑝𝑝−2𝐽𝐽−2 𝐷𝐷2

+

Page 22 of 29

𝐶𝐶3,𝐽𝐽 𝑓𝑓3,𝐽𝐽 (𝑧𝑧)(1−𝑧𝑧)𝑝𝑝−2𝐽𝐽−4 𝐷𝐷4

22 + ⋯] ,

and similarly for terms higher order in J(J + 1). The powers of (1 - z) multiplying f1,J(z) were determined by LeRoy and Barwell, 15 and the powers of (1 - z) multiplying f2,J(z), etc. follow the experience of Eq. (10). Here, b2,v is the centrifugal distortion term (usually denoted Dv); b3,v takes the place of the Hv term, etc. As before, the functions, fi,J, are taken to be polynomials in z whose constant term is 1 by factoring out the Ck,J. Note again that these expressions are only asymptotic because p is finite. The first few terms in the expansions were determined from algebraically determined Yn,J expressions, and the values of fk,J(1) were estimated from numerically determined eigenvalues. Higher-order expansion coefficients not determined algebraically were fit to computed numerical eigenvalues. The estimated uncertainty in the numerically fit coefficients increases as k and J increase. Analytic formulas and fit data for Ck,J, and fk,J(z), etc. are included in the supporting information. The procedure for calculating rotational terms in the energy could be done for any LJ m-n potential, and the fitting parameters would be expressed as algebraic functions of m and n. For other interatomic potentials, the fitting parameters would be algebraic functions of the potential parameters. III. Conclusion I argued in the introduction that the m value in the LJ m-n potential can adjust the equilibrium distance of the potential minimum. However, the LJ m-n potential family cannot be used universally. For example, if you try to use a LJ m-n potential to fit the H2 diatomic molecule with experimental spectroscopic constants ωe = 4395 cm-1, D0 =4.476 eV, and Re = 74.16 pm, 16 you calculate a harmonic length unit of (h/νµ)1/2 = 12.3 pm and in reduced units D =


Page 23 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


23 1/2

8.71 and Re = 6.01 pm . If you now require n = 6 in Re = (mnD) , it means that m = 0.69, which is less than n and disallowed. This same situation appears for all strongly bound small molecules like HCl, HF, O2, N2, and even for more weakly bound molecules like Li2 and Na2. In fact, an exponential Morse potential may be appropriate in these because covalent binding in molecular orbital theory can be viewed as being described by overlap and exchange integrals that have exponential dependence on R. On the other hand, a LJ m-n potential can be used profitably in more weakly bound and larger molecules like rare gas dimers, 17,18,19,20 where n = 6 and the computed value of m is between 12 and 13.2. Another area for profitable application of the LJ m-n potential is for complexes of polar species. One of the vibrational modes of such a complex correlates to separating the centers of mass. In this case, n = 3, and the weak bond and large equilibrium distance of the energy minimum easily puts the value of m greater than 5. I have successfully applied this model to one of the vibrational motions of water dimer and other hydrogen bonded species. 21 Fitting of the fi(z) functions becomes precarious for m too close to n because fi(z) functions develop singularities at z = 1. I recommend using m ≥ n + 1 except when absolutely necessary. This choice also avoids the parameter region when determining Dcrit(v) becomes difficult in the numerical determination of G described in the supporting information. This recommendation, however, limits the range over which Re may be adjusted. If one were to use a harmonic approximation or even a Morse approximation for the quantum states of a floppy stretching vibration, it could seriously undercount the populations of excited states. For example in the water dimer, using a LJ m-3 model for a dipole-dipole potential includes more that five times as many bound states below dissociation than a harmonic



Page 24 of 29

24 potential has below that same energy (see Fig. 3), and compared to a Morse potential with the same dissociation energy, a LJ m-3 potential has more than 2.5 times as many bound states. Since the hydrogen bond in water dimer is relatively weak, these extra states could be significantly populated at some temperatures important to atmospheric physics or chemistry. Consequently, using long-range model potential energies where appropriate could be important in partition function calculations. This work is similar to previous work tabulating Dunham coefficients to higher order than Dunham’s7 original paper using either higher order WKB results 22 or quantum perturbation theory. 23,24,25 However, the polynomial potential expansion about the equilibrium geometry will not converge for a potential that becomes constant at large separations, and the expressions for the higher order terms become increasingly complicated. It seems reasonable to break off the WKB series at four terms where they are still tractable and yet sufficiently accurate for most purposes. Using the interpolating functions, fi(z), to connect to the correct dissociation limit is a useful and intuitive idea that deserves consideration. It might be possible to achieve this using numerical quadrature as suggested by Barnwell et al. 26 However, I chose empirical fitting to make this connection, which appears to be successful. These results should increase the ability of researchers to use simple long-range model potentials to accurately take into account the quantum states of weakly bound systems.

Acknowledgment This work was supported by NSF grant AGS 1238947. Supporting Information 1. Notes on Born-Sommerfeld quantization and the Gmn factor


Page 25 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


25 2. Notes on the Numerov calculation of quantum eigenvalues 3. Notes on numerical determination of G factor which determines the maximum number of bound states, Figure S1. G factor fitting parameters. 4. Table S1. Table of C3 values 5. Table S2. Table of C4 values 6. Table S3. Quartic fit parameters and maximum errors using the approximations of Eq. (15) and Eq. (16). 7. Table S4. Exp-6 fitting parameters 8. Figures S2-S8. Plots of fi(z) functions for Lennard-Jones m-n potentials. 9. Formulas for ai coefficients for Lennard-Jones m-n potentials. 10. Analytic expressions for the Ak contour integration expansion variables in terms of the ak potential expansion variables. 11. Analytic expressions for the Bk contour integration expansion variables in terms of the ak potential expansion variables. 12. Analytic expressions for the Ck contour integration expansion variables in terms of the ak potential expansion variables. 13. Analytic expressions for the Dk contour integration expansion variables in terms of the ak potential expansion variables. 14. Analytic expressions for the Fk contour integration expansion variables in terms of the ak potential expansion variables. 15. Analytic expressions for the Ek contour integration expansion variables in terms of the ak potential expansion variables. 16. Analytic expressions for the Hk contour integration expansion variables in terms of the ak potential expansion variables. 17. Analytic expressions for the Gk contour integration expansion variables in terms of the ak potential expansion variables. 18. Analytic expressions for Y(v,j) Dunham coefficients calculable with four WKB terms in terms of the ak potential expansion variables. 19. Analytic expressions for the Ck,J coefficients, scaling constants in the expansion given in Eq. (20-21) 20. Analytic expressions for the fij,J coefficients calculable with four WKB terms. 21. Analytic expressions for fi,0(1) values. 22. Tabulated values of f1(z) functions 23. Tabulated expansion coefficients for fitting functions used in Eq. 20-21. References 1

McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976, Ch. 9.

2

Herzberg, G. Infrared and Raman Spectra of Polyatomic Molecules; Van Nostrand Reinhold: New York, 1945; Ch. 2.



Page 26 of 29

26 3

Morse, P. M. Diatomic Molecules According to the Wave Mechanics. II. Vibrational Levels. Phys. Rev. 1929, 34, 57–64.

4

M. L. Sage, Morse Oscillator Transition Probabilities for Molecular Bond Modes. Chem. Phys. 1978, 35, 375–380.

5

Cooksy, A. Quantum Chemistry and Molecular Interactions; Pearson: Boston, 2014, Ch. 10.

6

Lennard-Jones, J. E. On the Determination of Molecular Fields. II. From the Equation of State of a Gas. Proc. R. Soc. Lond. A 1924, 106, 463–477.

7

Dunham, J. L. The Energy Levels of a Rotating Vibrator. J. Chem. Phys. 1932, 41, 721–731.

8

Siska, P. E.; Parson, J. M.; Schafer, T. P.; Lee, Y. T. Intermolecular Potentials from Crossed Beam Differential Elastic Scattering Measurements, III. He + He and Ne + Ne. J. Chem. Phys. 1971, 55, 5762–5770.

9

Parson, J. M.; Siska, P. E.; Lee, Y. T. Intermolecular Potentials from Crossed Beam Differential Elastic Scattering Measurements, IV. Ar + Ar. J. Chem. Phys. 1972, 56, 1511– 1516.

10

LeRoy, R. J. ; Bernstein, R. B. Dissociation Energies of Diatomic Molecules from Vibrational Spacings of Higher Levels: Application to the Halogens. Chem. Phys. Letters 1970, 5, 42–44.

11

LeRoy, R. J.; Bernstein, R. B. Dissociation Energies and Long-range Potential of Diatomic Molecules from Vibrational Spacings of Higher Levels. J. Chem. Phys. 1970, 52, 3869–3879;

12

Kesarwani, R. N.; Varshni, Y. P. Four-term WKBJ Approximation and its Application to Varshni V and Lennard-Jones Potentials. Can. J. Phys. 1980, 58, 363–369.


Page 27 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


27 13

Kirschner, S. M.; LeRoy, R. J. On the Application, Breakdown, and Near-dissociation Behavior of the Higher-order JWKB Quantization Condition, J. Chem. Phys. 1978, 68, 3139– 3148.

14

R. Paulsson, R.; F. Karlsson, F.; LeRoy, R. J. Reliability of High-order Phase-integral Eigenvalues for Single and Double Minimum Potentials. J. Chem. Phys. 1983, 79, 4346– 4354.

15

LeRoy, R. J.; Barwell, M. G. Ground State D2 Dissociation Energy from the Near-dissociation Behavior of Rotational Energy Spacings, Can. J. Phys. 1975, 53, 1983–1990.

16

Herzberg, G. Infrared and Raman Spectra of Polyatomic Molecules; Van Nostrand Reinhold: New York, 1945; pp. 531–532.

17

Aziz, R. A.; Chen, H. H. An Accurate Intermolecular Potential for Argon. J. Chem. Phys. 1977, 67, 5719–5726.

18

Aziz, R. A.; Nain, V. P. S.; Carley, J. S.; Taylor, W. L.; McConville, G. T. An Accurate Intermolecular Potential for Helium. J. Chem. Phys. 1979, 70, 4330–4342.

19

Aziz, R. A.; Slaman, M. J. The Ne-Ne Interatomic Potential Revisited. Chem. Phys. 1989, 130, 187–194.

20

Dham, A. K.; Meath, W. J.; Allnatt, A. R.; Aziz, R. A.; Slaman, M. J. XC and HFD-B Potential Energy Curves for Xe-Xe and Related Physical Properties. Chem. Phys. 1990, 142, 173–189.

21

Shirts, R. B.; Kumbhani, S. R.; Burrell, E.; Hanson, J. C. An Improved Model to Calculate Equilibrium Constants for Formation of Peroxy Radical–Water Complexes. Theor. Chem. Accounts. 2018, 137, 96.



Page 28 of 29

28 22

Kesarwani, R. N.; Varshni, Y. P. Five Term WKBJ Approximation, J. Math. Phys. 1980, 21, 90–92.

23

Bouanich, J. P. Higher-order Contributions to the Dunham Coefficients: Application to the ground state of CO, J. Quant. Spectrosc. Radiat. Transfer, 1978, 19, 381–386.

24

Ogilvie, J. F.; Bouanich, J. P. Further Dunham Energy Coefficients of Diatomic Molecules, J. Quant. Spectrosc. Radiat. Transfer, 1982, 22, 481–482.

25

Rytel, M.; Rytel, T. Dunham Series Coefficients up to 20th Order, J. Mol. Spectroscopy, 1997, 185, 417–419.

26

Barnwell, M. G.; LeRoy, R. J.; Pajunen, P.; Child M. S. Evaluation of High-order JWKB Phase Integrals, J. Chem. Phys. 1979, 71, 2618–2622.


Page 29 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic Manuscript jp-2018-080957 R. B. Shirts

1.12

f1 LJ12-6

1.08

D=64

1.04

D=16 D=4

1 0

0.2

0.4

0.6

0.8

1

z


Connecting the Dunham Expansion to the Dissociation Limit for

Recommend Documents