===== Numerical Exercises for Tuesday July 15 =====


In today's exercises, we complete the tasks from yesterday, which unfortunately had ambiguous and confusing instructions!

==== The easy stuff ====

  * Install Git, and try out some of the commands covered in Morten or Nicolas's lecture slides ({{:computing.pdf|}} {{:talentdftguides.pdf|}}) for your codes in the following problems.

  * In your favorite programming language, make a program to construct a real symmetric $NxN$ matrix. Diagonalize it using the appropriate LAPACK or GSL routine, and write out some number of the lowest eigenvalues. (Suggestion: You might find it useful to use Mathematica (available on the ECT* computers) to diagonalize a small matrix that you can benchmark against.) This will help you test that you've linked to the GSL or LAPACK library and things are properly installed.


=== Harmonic Oscillator wf's and matrix elements to solve the Hydrogen atom ===

Here, our goal is to do some simple exercises to help us build useful tools such as subroutines to calculate HO wf's and matrix elements $\langle n'l|V|nl\rangle$ that are useful building blocks for more complicated calculations (e.g., Hartree-Fock calculations of many-body systems). Our ultimate test will be to see if we can get the ground state energy of the hydrogen atom (-.5 in atomic units) when we solve the Schr\"odinger equation by matrix diagonalization in the HO basis.

  * In the code {{:coulomboscrelme.f90.zip|here}}, find the subroutine that computes the generalized laguerre polynomials $L_n^{l+1/2}(x)$, "laguerre_general". NOTE: this is not a standalone code. Rather, the purpose is to see the simple recursion algorithm for calculating the Laguerre polynomials, and by extension $R_{nl}$, which you can then implement according to your own taste.

Recall that the HO wavefunctions are given by 
\begin{equation}
R_{nl}(r) = \frac{A_{nl}}{b^{3/2}}\xi^le^{-\xi^2/2}L_{n}^{l+1/2}(\xi^2)\,,
\end{equation}

where the oscillator length $b\equiv \sqrt{\hbar/(m\omega)}$, and the dimensionless variable $\xi\equiv r/b$. The normalization constant is given by

\begin{equation}
A_{nl}=\sqrt{\frac{2^{n+l+2}n!}{\pi^{1/2}(2n+2l+1)!!}}
\end{equation}

  * Using the uploaded code for guidance, make your own subroutine to calculate  $L_{n}^{\alpha}(x)$. Then, using the definition of the HO wf's, make a subroutine/function to calculate $R_{nl}(r)$. 

  * Why should you not calculate the factorial and double factorial functions appearing in $A_{nl}$ naively as they are written? Hint: What happens for large values of $n$ and $l$? How might you avoid this problem? Hint: What is $\log{ABC\cdots}$? 

  * As a check on your code, check the orthogonality of the $R_{nl}$,

\begin{equation}
\int_{0}^{\infty}R_{nl}(r)R_{n'l}(r)r^2dr= \delta_{nn'}
\end{equation}

  * As discussed in the HO notes, the "best" type of Gaussian quadrature over HO wf's is Gauss-Laguerre. However, as you see in the code {{:coulomboscrelme.f90.zip|here}}, you can use Gauss-Legendre as well. You need to play around to make sure your upper limit $r_{max}$ of your discretization range is sufficiently large. (Note: there are plenty of freely downloadable Gauss-Legendre subroutines. Please feel free to use one instead of writing it from scratch.) In this way, the orthogonality check takes the form
 
\begin{equation}
\int_{0}^{\infty}R_{nl}(r)R_{n'l}(r)r^2dr\approx \sum_{i=1}^{N}w_i r_i^2 R_{nl}(r_i)R_{n'l}(r_i)\,,
\end{equation}

where the sum is over the N mesh points and $w_{i}$ are the quadrature weights.

  * Now, write a subroutine to evaluate the matrix elements $\langle nl|V|n'l\rangle$ for some arbitrary central two-body potential $V(r)$. 

  * To test your $\langle nl|V|n'l\rangle$ routine, apply it to the coulomb potential in the hydrogen atom problem. Use atomic units ($\hbar=e=1/4\pi\epsilon_0=m_e=1$).  Now, build the Hamiltonian matrix $\langle nl|H|n'l\rangle$ (you may take $l=0$) for the hydrogen atom, keeping states $n,n'< N_{max}$. Use the analytic expressions for the K.E. matrix elements in {{:ho_spherical.pdf|}}. Now diagonalize the matrix to get the ground state energy. 

  * For a given $N_{max}$, repeat the calculation for many different values of oscillator length and plot the ground state versus $b$. Repeat for larger values of $N_{max}$ and put the $E$ versus $b$ curves on the same plot. Are they behaving according to expectation? (Hint: the diagonalization in a finite basis is variational.) Should results depend on $b$ as $N_{max}\rightarrow\infty$?