TalentDFT

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

• 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.

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 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 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$?