5.73 Lecture 9: Numerov-Cooley Method: 1-D Schrodinger Equation

5.73 Lecture #9 9-1 Numerov-Cooley Method : 1-D Schr. Eq. Last time: Rydberg, Klein, Rees Method and Long-Range Model G(v), B(v) rotation-vibration constants VJ(x) potential energy curve x = R – Re E v ,J ,ψv ,J , all conceivable experiments including wavepackets, wp (x , t ) wp(x,t) = ∑ ai ψ i e − Eit ! Treated in Lectures #19 and #20 i * Initial preparation of wp: ai = ∫ ψ i ⎡⎣ wp(x,0)⎤⎦ dx determined by VJ(x) free evolution of wp Completeness: Any arbitrary wavefunction may be expanded as a linear combination of eigenstates of H. Method for RKR derivation: A(E,J) = area of V(x) below E: used WKB Quantization Condition obtained x±(E,J) – a series of pairs of turning points on a grid of E. Today: What do we do when we have a potential energy curve, VJ(x), (especially when V(x) is not suited for WKB)? No models Solve Schr. Eq. numerically! 15 digit reproducibility cheap This is the final tool we will develop for use in the Schrödinger representation. To summarize the classes of 1–D problem we are able to solve: * piecewise constant potentials (matrix approach for joining at ψ(x) boundaries) * Airy functions (linear potential and joining JWKB across turning point) * JWKB (quantization condition and semi-classical wavefunctions) * numerical integration (today) See Merzbacher Quantum Mechanics pp. 92-103 5.73 Lecture #9 9-2 Numerical Integration of the 1-D Schrödinger Equation WKB OK here NOT HERE widely used incredibly accurate no restrictions on V(x) or on E–V(x) [e.g. nonclassical region, near turning points, double minimum potential, kinks in V(x).] For most 1-D problems, where all one cares about is a set of {Ei, ψi}, where ψi is defined on a grid of points xi, one uses Numerov-Cooley leroy.uwaterloo.ca/programs.html LEVEL: A computer program for solving the radial Schrödinger equation for bound and quasibound levels R.J. Le Roy, Journal of Quantitative Spectroscopy and Radiative Transfer 186, 167-178 (2016). RKR1 – A Computer Program implementing the first-order RKR method for determining diatomic molecule potential energy functions R.J. Le Roy, Journal of Quantitative Spectroscopy and Radiative Transfer 186, 158-166 (2016). J. A.Tellinghuisen, "Potential Fitting RKR Method: Semiclassical vs. Quantal Comparisons". Journal of Molecular Spectroscopy. 330, 20-27 (2016). Supplements for Lecture #9: 1. Cooley, Math. Comput. 15, 363 (1961). 2. Press et. al., Numerical Recipes, Chapters 16 and 17 3. Classic unpublished paper by Zare and Cashion with listing of Fortran program [now better to see LeRoy web site] 4. Tests of Numerov-Cooley vs. other methods by Tellinghuisen 5. LeRoy JQSRT paper Basic Idea: grid method * solve differential equation by starting at some xi and propagating trial solution from one grid point to the next * apply ψ(x) = 0 boundary conditions at x = 0 and ∞ by two different tricks and then force agreement between �left(x) and �right(x) at some intermediate point by adjusting E. 5.73 Lecture #9 9-3 Euler’s Method h • x0 want ψ(x) at a series of grid points x0, x1, …xn = x0+nh ψi = ψ(xi) call these Need a “generating function,” f(xn, ψn) h • x1 h • x2 • prescription for going n → n + 1 must depend on both xn and ψn. xn samples potential, ψn samples previous value of ψ. ψ n +1 = ψ n + hf ( xn ,ψ n ) ↑$ increment in x xn+1 – xn = h [NOT Planck’s constant] (ψn is a number, not the entire wavefunction.) For the Euler method, the generating function is simply: f ( xn , ψ n ) = dψ ψ − ψ n ψ n +1 − ψ n ≈ n +1 = dx xn xn +1 − x n h The value of this derivative actually comes from the differential equation thatψmust satisfy, not from prior knowledge of ψ(x) (which we do not yet have!) d 2ψ 2µ For the Schrödinger Eqn. = – (E − U(x))ψ !2 dx 2 All constants absorbed in V(x). V(x) must be in units of Å–2. −p(x)2 note that V(x) = −k(x)2 = !2 V (x ) ≡ C [U (x ) − E ] U (x ) is the potential. d 2ψ = V (x )ψ(x ) C = 10 ( 8π cµ h ) h here is not increment of distance, Planck dx 2 h and x are in Å. E and U (x ) are in = 0.0593203146µ dψ ψ − ψi gram ⎛ ⎞ ≡ i +1 , C ⎟ cm units (E /hc ) ⎜⎝ µ in amu = ⎠ mole dx xi h not Planck −16 2 A 12 -1 A d 2ψ dx 2 xi ⎧ ⎡ ψ − ψ i ⎤ ⎡ ψ i − ψ i −1 ⎤ ⎫ = ⎨ ⎢ i +1 −⎢ ⎬ h ⎥ ⎥ h h ⎦ ⎣ ⎦⎭ ⎩⎣ µA = = h −2 [ ψ i +1 − 2ψ i + ψ i −1 ] (Planck’s constant h is in italics, grid spacing h is not italic.) m1m 2 m1 + m 2 5.73 Lecture #9 9-4 Schr. Eq. tells us the rule for propagating ψ. Employing Euler’s method (h is not Planck’s constant): h −2[ ψ i+1 − 2ψ i + ψ i −1 ] = Viψ i rearrange to solve for ψi+1 2 ψ i+1 − 2ψ i + ψ i−1 = h Vi ψ i ψ i+1 = 2ψ i − ψ i−1 + h2 Vi ψ i a recursion relationship. Need both ψi and ψi–1 to get ψi+1. In order to get things started we need two values of ψ starting at either edge of the region where ψ is defined and ψ starts out very small. [See Press et. al. handout (Numerical Recipes, Chapters 16 and 17) for discussion of nth-order Runge-Kutta method. The generator is chosen more cleverly than in the Euler method so that stepping errors are minimized by taking more derivatives at intermediate points in the xi, xi+1 interval.] Cooley specifies yi+1 = 2 yi − yi−1 + h 2Vi ψ i ( ) yi = ⎡⎣1− h 2 12 Vi ⎤⎦ ψ i (and vice versa) * use ψι to get yi * use ψi and yi (and yi–1) to get yi+1 * use yi+1 to get ψi+1 The result is that the error in ψ i+1 is on the order of h6 ψ iVi — smaller error if h is smaller 240 (Cooley’s method is much better than Euler’s, which converges as a smaller power of h) 5.73 Lecture #9 9-5 So what do we do? 0 –De 0 R–(E) Re R+(E) R turning points [ e.g., VMORSE (R) = De 1 − e ] −β (R−R e ) 2 − De [zero of V at dissociation] V(∞ ) = 0, V(R e ) = −De or, more generally than Morse, use VJ,RKR. x = R – Re (displacement from equilibrium) at R = 0 x = –Re ψ(–Re) = 0 ⎫ R=∞ ψ(∞) = 0 ⎬boundary conditions ⎭ The two boundary conditions are handled differently, because we want to define a finite # of equally spaced grid points (not actually necessary — see Press: he uses a variable grid spacing, which is needed to sample an infinite range of x with a finite number of grid points) experience — doesn't really matter for ψ1 * at R = 0 ψ0 ≡ 0 (required) –20 (arbitrarily chosen small ψ1 = 10 number, which will be corrected later upon normalization) Use this to start the integration outward. When we discover later that we made a wrong choice for ψ1, this is corrected merely by dividing all ψi i ≥ 1 by an i-independent correction factor. 5.73 Lecture #9 9-6 At large R (in the classically forbidden region), choose ψn at the last grid point, xn, to be small and use WKB only once to compute the next to last grid point. We do this because we have no reason to extend the calculation to x → ∞. ψn = 10–30 (the final grid point) ψ WKB = p −1/2 (based on experience) ⎡ 1 x ⎤ exp ⎢ − ∫ p dx ⎥ ⎣ ! R+ ( E ) ⎦ integrating inward from outer turning point. p2 2m V (x ) ≡ − 2 ( E − U (x )) = 2 ! ! p n = ⎡⎣ 2 m ( E − U n ) ⎤⎦ U n = U (xn ) 1/2 ≡ !Vn1/2 WKB gives us the ratio ψ n−1 = ψn ⎡ 1 xn−1 ⎤ exp ⎢ − ∫ pn−1dx ⎥ ⎣ ! R+ ( E ) ⎦ −1/2 ⎡ 1 xn ⎤ pn exp ⎢ − ∫ pn dx ⎥ ⎣ ! R+ ( E ) ⎦ pn−1 −1/2 pn–1 ≈ pn so pre-exponential terms cancel. Evaluate integral in exp as a sum. ∫ x+ x− f (x )dx = i ( x+ ) ∑ hf (x ) i i ( x− ) 5.73 Lecture #9 9-7 fi+2 fi+3 fi+1 fi x h and e a +b +c h h h h h e a +b +c 1 = e e e , and a +b +c +d = d . e e a b c But we chose ψn = 10–30 thus ψ n −1 = 10 −30 exp ⎡⎣Vn1/2 h ⎤⎦ So we have ψn and ψn–1 and we are off to the races. ψ 0 → ψ 1 → ψ 2 →→→ ψ m ψ m ←←← ψ n −1 ← ψ n m meet in the middle at ψm. 5.73 Lecture #9 9-8 Stop the inward propagation of ψ when a point is reached where, for the first time, |ψm| > |ψm+1|. This locates the first maximum of ψ inside the outer turning point. Since |ψi| is exponentially increasing from 10–30 at i=n until it reaches its first maximum inside the classically allowed region, this outer lobe of ψ is also the most important feature of ψ (because most of the probability resides in it). outermost lobe of ψ$ F maximum of |ψ| at ψ ( xm ) V(x) asymptotic approach to 0 xm Use the outermost lobe because this is the global maximum of ψ(x). This choice minimizes the problem of precision being limited by the finite number of significant figures in the computer. Set value of ψm = 1.0 by renormalizing both functions replace each ψ i by ψ i ψ m * ψ from n, n − 1, …m ! (from the right) * ψ from i = 0, 1, …m for all i down to m. ! replace each ψ i by ψ i ψ m (from the left) for all i up to m. The renormalized ψ's are denoted by ψ′. ψi ψ ′i = ψm ψ ′m = 1 ψ′ must be continuous, especially at the joining grid point, m. 5.73 Lecture #9 9-9 This ensures that ψ(x) is continuous everywhere and that it satisfies grid form of Schr. Eq. everywhere except at i = m. ( ) 0 = –ψ i+1 + 2ψ i − ψ i−1 + h 2Vi ψ i In order to satisfy Schr. Eq. for i = m, it is necessary to adjust E. The above equation can be viewed as a nonlinear requirement on E. At the crucial grid point i = m, define an error function, F(E). ( ) ( ) E E E E F E = –ψ m+1 + 2ψ m − ψ m−1 + h 2VmE ψ m where, by varying E, we want to search for the zeroes of F(E). Assume that F(E) can be expanded about E1 (E1 is the initial, randomly chosen value of E.) F ( E ) = F ( E1 ) + dF dE (E − E ) + 1 E1 discard higher terms (linearization approximation): Newton-Raphson and solve for the value of E where F(E) = 0. Call this E2 0 = F (E 1 ) + dF (E 2 − E1 ) dE E 1 F(E1 ) E2 = − + E1 (dF dE )E 1 Correction to E1 This gives an estimate of where the zero of F(E) nearest E1 is located. 5.73 Lecture #9 9 - 10 Usual approach: compute dF dE = E1 F ( E1 + δ ) − F ( E1 ) δ Once the derivative is known, use it to compute the trial correction to E1 (assuming linearity). Newton-Raphson method for solving nonlinear equation E2 = E1 + Δ Δ≡– F ( E1 ) ( dF dE )E Are we done? Not quite done if F(E) is not sufficiently linear near Ei. Keep iterating! 1 Iterate until the correction, ∆, to E is smaller than a pre-set convergence criterion, ε. This seems like a lot of work, but the computer does not complain. Now we have one not-yet-normalized eigenfunction of H and the corresponding eigenvalue, E. Normalize ψ E by dividing by n ∑ ∫ ψ * ψ dx = box normalized: ( ) ψ E xi = i=0 ψi 1/2 ∫ ψ * ψ dx 2 ψi h 1/2 = NE Recall that normalization integral is evaluated by summation over grid points. ψreal for bound 1-D function ⎡ ⎤ 2 ⎢∑ ψ j h ⎥ ⎢⎣ j ⎥⎦ This procedure has been used and tested by many workers. The latest and extensively documented version, “Level 16” is obtainable at Robert LeRoy’s web site: http://dx.doi.org/10.1016/j.jqsrt.2016.05.028 LEVEL: A computer program for solving the radial Schrödinger equation for bound and quasi-bound levels, Robert J. LeRoy. I will assign some problems based on Numerov-Cooley method for integrating the 1-D Schr. Eq.