A numerical solution of the time-independent Schroedinger equation for hydrogenic atoms with reduced mass m@me.

The program Hatom1.txt is written in QuickBasic.  It demonstrates a method for solving the radial equation for a hydrogenic atom.  This is a boundary value problem.  We solve for the eigenvalues E of the Hamiltonian for which the wavefunction and its derivative are continuous and for which the wavefunction takes on given values at the boundaries of the integration volume.  The program plots the radial functions Rkl(r) corresponding to the calculated eigenvalues.  The nuclear charge Z and the angular momentum quantum number l can easily be varied.

Explanation:

The wavefunction y(r) is of the form We solve the radial equation for ukl(r).

We integrate this equation using the Numerov method.  (See Program 3, QM1.)   For each l we find all eigenvalues Ekl<0.5Z2eV.  We denote the successive energy levels by k=1, 2, 3, ... .  We expect that Ekl depends only on n=k+l.

We assume that ukl(r) is zero at r=0 and at r=rmax=6(l+1)Z2.

As r goes to zero the effective potential varies rapidly, while for large r it is close to zero and varies slowly.  Integration in regions where k varies rapidly requires small steps, while larger step sizes can be used when k varies slowly.  Using a small step size everywhere increases CPU time and numerical inaccuracies, while using a large step size everywhere will not produce an acceptable result.  We divide the region between r=0 and r=rmax into two intervals to demonstrate how variable step sizes can be used.  In the first interval the step size is 100 times smaller than in the second interval.  Ideally we would like to change the step size continuously.  We will demonstrate how this can be done in  another program.

We start integrating at r=dr.  (At r=0  Veff blows up resulting in a "divide by zero" error.)  We choose a starting energy lower than the expected ground state energy.  We integrate outward to rmax, counting the number of nodes in the wavefunction as a function of energy, as we increase the energy in small steps DE.  When the number of nodes changes from one energy to the next the wavefunction is close to zero at r=rmax and we are close to an eigenvalue.  We decrease DE once and then take as the eigenvalue the average value of the energy just before and just after the increase in the number of nodes.  This interpolation can easily be refined, for example by using a secant search.  (See Program 3, QM1.)

The program Hatom1.txt plots the radial functions Rkl(r) from r=0 to r=50 Å corresponding to the calculated eigenvalues.

Link: