Hydrogen Atom

Hydrogen Atom#

Author:
Date:
Time spent on this assignment:

Throughout this exercise, we are going to work in atomic units where \(\hbar=m=e=1/(4\pi\epsilon_0)=1\). This gives us the Bohr radius which is also equal to 1 in these units (and \(0.529 \times 10^{-10}\) m in SI units). All distances in this unit will be measured in Bohr radii.

Exercise 1. Plotting the Hydrogen Atom Orbitals#

In this exercise, we are going to plot the orbitals of the Hydrogen atom starting from the relevant solutions that you found analytically. In future exercises, we will see how to get these results from scratch without using the analytic solutions.

Recall that the Hydrogen atom wave-function is

\[ \Psi_{nlm}(r,\theta,\phi) = R_{nl}(r) Y_{lm}(\theta,\phi)\]

a. Spherical Harmonics#

Our first task is to look at and plot the spherical harmonics \(Y_{lm}(\theta,\phi)\).

To get the spherical Harmonics in python we can call

sp.sph_harm(m,l,theta,phi)

where the \(m\) and \(l\) are indexing the eigenstates of the angular momentum operator. (Notice that this is backwards from how you typical think about the quantum numbers for the hydrogen atom). \(\theta\) and \(\phi\) are in spherical coordinates.

It’s going to turn out that we will want the spherical harmonics on many different values of \(\theta\) and \(\phi\) at once. To do that, we are going to have to build some grids.

The first grid, we will build is a two-dimensional grid of \(\theta\) and \(\phi\) values using

phi_1d=np.linspace(0,np.pi,100)
theta_1d=np.linspace(0,2*np.pi,100)
theta_2d,phi_2d=np.meshgrid(theta_1d,phi_1d)

Using these 2d grids, we can now compute the values of \(Y_{lm}\) for every \((\theta,\phi)\) pair at once by calling the spherical harmonic function on these 2d grids. You should get out a \(100 \times 100\) matrix (let’s call it Ylm_2d).

Now we want to figure out how to plot these spherical harmonics. The spherical harmonics are a function of two variables (\(\theta\) and \(\phi\)). So one thing we could do is let the r coordinate represent the value of the Ylm(\(\theta\),\(\phi\)) and then plot \((r,\theta,\phi) = (Y_{lm}(\theta,\phi),\theta,\phi)\).

To do this, write a function x,y,z=SphericalToCartesian(r,theta,phi) which takes spherical coordinates to cartesian coordinates - i.e x=r * np.sin(theta)* np.cos(phi) - and then send your grid as well as the Ylm_2d result (i.e. our \(r\) values) to get the (x,y,z) coordinates.

We can now plot

fig = plt.figure()
ax = fig.add_subplot( 111 , projection='3d')
ax.plot_surface(x, y, z, linewidth = 0.5, edgecolors = 'k')
plt.axis("equal")

Plot spherical harmonics for

\((\ell,m) = (1,0)\)
\((\ell,m) = (2,0)\)
\((\ell,m) = (3,2)\)

b. Radial Function#

Now that we understand how to get the spherical Harmonics, our next goal is to get the radial wave-function of the Hydrogen atom working. Recall that the radial wave-function is

\[ R_{nl}(r) = \frac{u_{nl}(r)}{r} = N \exp(-\rho / 2) \rho^l L^{n-l-1}_{2l+1}(\rho) \]

where

\[\rho = \frac{2r}{na_0}\]

and

\[N = \sqrt{\left(\frac{2 }{n a_0}\right)^3 \frac{(n - l - 1)!}{2n (n + l)!}}\]

and \(L\) is a special function which can be generated by doing L=sp.genlaguerre(n-l-1,2*l+1). Note that this just returns the python function. Later you will have to call L(rho).

Write a function radial_function(n, l, r) which takes \(n\), \(l\), and \(r\) and returns \(R_{nl}(r)\). Recall we are in units where \(a_0=1\).

Let’s start by evaluating this on a one-dimensional grid of positions r=np.linspace(0.001,45,1000).

Plot on the same plot, \(u_{nl}(r)\) (notice this is different from \(R_{nl}\) by a factor of \(r\)) for

\((n,\ell) = (1,0)\)
\((n,\ell) = (2,0)\)
\((n,\ell) = (3,0)\)

How do the number of nodes (times when it hits zero) correspond to the value of \(n\)?

Now also plot

\((n,\ell) = (3,2)\)
\((n,\ell) = (3,1)\)
\((n,\ell) = (3,0)\)

Notice how the larger values of \(l\) are pushed further away from the origin. This is coming from the additional force due to the fake potential.

It should also be that the different \(u_{nl}(r)\) are orthogonal to each other. Check this by verifying that \(u_{21} u_{20}^\dagger=0\) You can do this by doing u21 @ u20.T.conjugate()

c. Hydrogen Orbitals#

We can now combine the previous two sections to write down the hydrogen orbitals. They are

\[\Psi_{nlm}(x,y,z) = R_{nl}(x,y,z)Y_{lm}(\theta,\phi)\]

We are going to want to evaluate these orbitals on the entire space. Therefore, we will start by putting together a three-dimensional grid of the (x,y,z) positions.

def HydrogenGrid(g=10,numGridPoints=50):
    p=np.linspace(-g,g,numGridPoints)
    h=p[1]-p[0]
    x,y,z=np.meshgrid(p,p,p)
    return x,y,z,h

You’ll adjust the grid length g depending on which orbital you’re trying to plot.

Write a function Cartesian2Spherical(x,y,z) which returns \(r\), \(\theta\), and \(\phi\). Recall that

\[\theta = \tan^{-1}(y/x) +\pi \]

(use np.arctan(y,x))

and

\[\phi = \cos^{-1}(z/r)\]

You should now have access to \(r,\theta,\phi\) on your grid from which you can call both your radial function and spherical harmonics. Write a function HydrogenAtom(n, l, m,x,y,z) which computes the hydrogen orbital on this grid and returns psi.

To plot the wave-function use the following

def PlotMe(psi,x,y,z):
    isosurface = go.Isosurface(
        x=x.flatten(),
        y=y.flatten(),
        z=z.flatten(),
        value=(np.real(psi)**2).flatten(),  # Use the absolute value for intensity
        surface_count=20, # Number of isosurfaces to display
        opacity=0.5,  # Set the opacity of the isosurface
        surface_fill=1.0,
        caps=dict(x_show=False, y_show=False,z_show=False)
    )
    fig = go.Figure(data=[isosurface])
    fig.show()

Plot \((n,\ell,m) = (1,0,0); (2,1,0); (3,2,2); (3,2,0)\) where the grid-length I used to get a good figure was

g=dict()
g[(1,0,0)]=4
g[(2,1,0)]=7.7
g[(3,2,2)]=18
g[(3,2,0)]=20

d. Getting the Energy of the Hydrogen Atom#

Now we are going to compute the energy of the Hydrogen atom. Recall that the Hamiltonian for the Hydrogen atom is

\[ H = -\frac{1}{2} \nabla^2 + \frac{1}{r} \]

where \(-\nabla^2 = \partial^2/\partial x^2 + \partial^2/\partial y^2 + \partial^2/\partial z^2\) and \(r\) is the distance between the nucleus and the electron.

Let’s figure out how to compute each of these terms.

This second term can be computed by taking the integral

\[ \int \Psi(x,y,z) \left(-\frac{1}{r}\right) \Psi^*(x,y,z) dx dx dz \]

If we evaluate our wave-function psi on the grid (which is how we got our orbitals that we visualize above), then we simply multiply these three terms, sum them up (np.sum) and then multiply by our grid spacing cubed.

For the \((n,l,m)=(1,0,0)\) you should get a potential Energy of -0.9897509703021152 using the following grid x,y,z,delta_x=HydrogenGrid(20,200)

Implement this function and make sure you get the correct energy.

Now we need to get the next term

\[ -\frac{1}{2} \int \Psi^*(x,y,z)\frac{\partial^2}{\partial x^2} \Psi(x,y,z) dx dx dz \]

We can compute this by finite differences. The standard formula for the second derivative is

\[ \frac{\partial^2}{\partial x^2} = \frac{f(x+dx)+f(x-dx)-2f(x)}{dx^2} \]

Write a function Laplacian(n,l,m) which calls HydrogenWaveFunction three times (per direction) returning the laplacian; for example for the x-direction it calls x+dx, x-dx, and x. For dx use 1e-3.

Compute the kinetic energy of the \((n,l,m)=(1,0,0)\) state. You should get 0.4897835961971023 total with each of the three kinetic terms being the same (because of the symmetry of this state)

Now sum up the kinetic and potential energies and multiply by 27.2 to get it into eV. Verify that you get the right answer for the (1,0,0) term. Also compute the energies for the (2,0,0), (2,1,0), (2,1,1) and (3,2,1) terms. Check that each of them are close to the known value.

Exercise 2. Solving the Hydrogen Atom Radial Equation Numerically#

Recall that the Hydrogen atom Hamiltonian is

\[H = -\nabla^2 + \frac{1}{r} \]

We solve the hydrogen atom using separation of variables eventually finding that

\[\Psi(R,\theta,\phi) = R_{nl}(r) Y_{lm}(\theta,\phi)\]

where \(R_{nl}(r) = u(r)/r\) and \(u(r)\) obeys the raidial eigenvalue equation

\[ \left[-\frac{\hbar^2}{2m}\left(\frac{d^2}{dr^2} - \frac{l(l+1)}{r^2} \right) - \frac{e^2}{4\pi\epsilon_0 r} \right] u(r) = E u(r)\]

which we can write as

\[H_{eff} u(r) = E u(r)\]

This problem looks a lot like the particle in the box. We can think of \(r\) here just like the \(x\) coordinate. Then the Hamiltonian has three terms:

the standard kinetic energy term
a potential term coming from the angular momentum
a potential term coming from the coulomb potential.

To solve this analytically required a significant amount of work. We will see that solving it numerically is relatively straightforward. For each of the terms in \(H_{eff}\) we need to write out a \(N \times N\) Hamiltonian matrix (where \(N\) is the number of points in our grid which discretizes space).

a. The Kinetic Term#

We will start by working with a discretized grid which is \(N=200\) points defined by

N=200
r = np.linspace(40, 0.0, N, endpoint=False)

Later to actually get accurate answers, we are going to want to increase \(N\) which will actually require us to make some changes (moving to sparse matrices) but let’s start with the smaller dense matrices.

We are building \(N \times N\) matrices where each row (respectively column) corresponds to a values of \(r\) on the grid - i.e. the “third” row corresponds to r[2] and the “fourth” column corresponds to the value of r[3].

Let’s start with the term

\[\frac{d^2}{dr^2} = \frac{2 f(r) - f(r+dr ) - f(r-dr)}{dr^2} \]

As a matrix this means:

the elements that connect the r’th column to the r’th row should get a \(2/dr^2\). This is just the diagonal of the matrix
the elements that connect the (r+1)’st column to the r’th row should get a \(-1/dr^2\). This is just one row above the diagonal
the elements that connect the (r-1)’st column to the r’th row should get a \(-1/dr^2\). This is just one row below the diagonal.

Write a function Laplacian(r) which takes the grid r and returns a matrix corresponding to the second derivative. You will find the np.diag function useful (even the off-diagonal term can be built using it - look at the documentation!). You should essentially be able to build this matrix in three such function calls.

Working in units where \(\hbar =1\) and \(m=1\), go ahead and compute the full kinetic term. To diagonalize a matrix, you can call e,v=np.linalg.eigh(L). Verify that the lowest four energies are [0.00305358 0.01221356 0.02747771 0.04884231].

#answer here

b. The potential terms#

The two potential terms are a function of \(r\) and ther matrices for these two terms therefore only connect the r’th row to the r’th column. This means both additional terms only exist on the diagonal. Build a function V(r,l) that takes the r grid and returns the matrix corresponding to the potential

\[V(r,l) = \frac{l(l+1)}{r^2} + \frac{1}{r}\]

Again use np.diag

#answer here

c. Putting it together#

Now we simply want to add the full Hamiltonian matrix together and diagonalize it using e,v=np.linalg.eigh(H)

Both print and plot the lowest 6 energies - plt.plot(e[0:6]*27.2). Also on your plot make some straight lines - plt.axhline(the_energy) where the energies for the n’th eignenstate should be - i.e. \(-13.6/n^2\). How close are you getting?

Also plot the lowest three eigenstates and see if the radial functions match the ones from above. To normalize the plot correctly you have to divide your eigenvectors by \(\sqrt{h}\).

#answer here

d. Finer Grids#

We are not getting an accuracy as high as we might desire. What we need is a better grid: finer and going out to larger \(r\). Unfortunately, if there are more grid points, the matrix gets larger and using dense matrices and solving for all the eigenstates becomes intractable. Instead we need to work with sparse matrices.

To do this, you need to do the following:

Change V and Laplacian to return sparse matrices. You can do this by just having the return be return scipy.sparse.csr_matrix(oldMatrix)
When you diagonalize, take e,v=scipy.sparse.linalg.eigsh(H,k=6,which='SA') This will give you sparse versions of the lowest six eigenstates.

Use now a grid of \(N=4000\) from \(L=-60\) to \(L=60\)

Go ahead and plot the same curves as above and verify you get more reasonable results.

#answer here

Exercise 3. Hydrogen Atom on a Grid#

So far we’ve plotted the orbitals for the Hydrogen atom by

using the formula for the orbitals that one derives analytically
using the radial equation you derived analytically (along with the spherical harmonics)

In both cases, we (somebody?) had to do a lot of analytic work to get something that we then worked with computationally. Now we just want to envision a situation where we have a Hamiltonian and want to directly work with and diagonalize this on a grid.

a. The Hamiltonian#

Our goal here is to be able to solve for the eigenstates and eigenvalues of the Hydrogen atom.

The hydrogen atom has a single electron and a Hamiltonian of

\[H = -\nabla^2 + V(R) \]

where

\[V(R) = \frac{1}{R}\]

where \(R\) is the distance between the proton and electron of the Hydrogen atom.

Let’s learn some things to put this together:

First recall that \(\nabla^2 = \frac{d^2}{dx^2}+\frac{d^2}{dy^2}+ \frac{d^2}{dz^2}\). We need to generate a matrix for this laplacian.

There’s something a little subtle here. \(d^2/dx^2\) looks like the same term that we saw for the particle in the box but it now needs to exist over three-dimensions. If \(L\) was the one-dimensional derivative, we instead now write this as \(L \otimes I \otimes I\) where \(I\) is the identity of the same size matrix.

Let’s think how large a matrix we are making here. We have a matrix that is number of grid points \(\times\) grid points. If we use 80 grid points in every direction, we get a matrix of size \(80^3 \times 80^3\) or \(512,000 \times 512,000\). This is much too large to be densely represented on a computer. To get around this, we are going to have to use sparse matrices. The trick here is to

first take your one-dimensional Laplacian (call it L) and turn it into a sparse matrix (i.e. scipy.sparse.lil_matrix(L)) and
then when you produce your three-dimensional laplacian make sure that you do a kroneker product - i.e. scipy.sparse.kron(scipy.sparse.kron(L,I),I) These two steps will give you \(\frac{d^2}{dx^2}\). Then you need to also get \(\frac{d^2}{dy^2}\) and \(\frac{d^2}{dz^2}\). Go ahead and generate the matrix then for \(\nabla^2\) (make sure in your laplacian you rememember to divide by \(1/h^2\) where \(h\) is the grid spacing).

Let’s check that our Laplacian is giving us the correct thing. We should be able to diagonalize it and get the energies for a particle in a three-dimensional box.

You can’t actually afford to get all the eigenvalues of such a large matrix. Instead, we are going to only get the first couple. To do this, you can use e,v=scipy.sparse.linalg.eigsh(-0.5*L) and show that the first few energies are correct for a grid

HydrogenGrid(12,numGridPoints=80)

After you have your laplacian in place, you want to generate the Coulomb potential piece of your Hamiltonian (\(1/R\)) on a grid.

The potential \(V(R)\) is a diagonal term of the Hamiltonian (i.e. it’s on the diagonal of the matrix). First compute the potential \(V(R)\) on your three-dimensional grid (call it Vext). To put it on the diagonal, you need to reshape your three-dimensional grid onto a one-dimensional vector which you can then make the diagonal of your matrix

  m=np.shape(Vext)[0]
  Vext=np.reshape(Vext,m*m*m)
  potential=scipy.sparse.diags(Vext)

Write a function ExternalPotential which generates this matrix.

You now have your laplacian term and your potential term. You can add these two terms together to get your full Hamiltonian.

#answer here

b. Eigenvalues#

Diagonalize your full Hamiltonian using scipy.sparse.linalg.eigsh on your \(12 \times 12 \times 12\) grid.

We would like to see that you get the expected eigenenergies and eigenstates.

Recall the standard hydrogen orbitals from chemistry each with an energy of \(-13.6/n^2\) eV

The lowest 6 energies for the Hydrogen atom should be

1s
2s
2px, 2py, 2pz
3s

Let’s start with looking at the energies of the 1s, 2s, and 3s states which we expect it to be -13.6 eV, -13.6/4 eV and -13.6/9 eV respectively. Obviously because we are using a grid, we aren’t going to get the exact answer. But as we increase the grid size we should get closer and closer to the true answer. Starting with 20 grid points per direction and increasing the grid in units of 20 up to 100 grid points per direction, graph the energy of these three eigenstates as a function of number of grid points per direction. Also draw horizontal lines at the expected value of the eigenstates. You should see the energies approach the target energies at large grid-size. Hint: Your g=20 energy for 1s should be -10.71473301. (There are some eigenstates where you need to increase the size of the grid to capture the full wave-function)

Note:. You need to be careful about units. To convert into eV you have to multiply your energy by 27.21.

c. Eigenvectors#

Now we can also look at the eigenstates. For this section, let’s use a grid with 80 points.

Let’s start by looking at the s-states. First, we will make some two-dimensional slices of the three-dimensional wave-functions we’ve found. To get your 2d visualization, you first need to reshape your wave-function into a 3d grid - if you have a \(80 \times 80 \times 80\) grid this can be done for the \(i\) th eigenstate by doing A=psi[:,i].reshape(80,80,80). Now you can go ahead and plot it by plt.matshow(A[40,:,:]) where 40 is the slice that you want. You probably want to use plt.colorbar() so you can see the magnitude and sign of things.

If you look at the 0 (1s),1 (2s) or 5 (3s) eigenstate you will find that it should have a clear “S” character - symmetric around the origin.

You can also graph the probability density as a function of the radius. You could average but because it is symmetric the easiest thing to do is just plot a one dimensional line starting at the origin. To match up with earlier results, actually we’d like to plot not the radial function \(|R_{nl}(r)|^2\) but \(|u_{nl}(r)|^2\) so multiply by \(r\) - i.e. plt.plot(r,r*A[40,40,40:]). Notice we are starting at 40: for that last term because that’s where zero starts. For the different shells, count how many nodes you see (i.e. places where the probability density goes to zero). To normalize you have to divide by \(1/(4\pi\sqrt{h})\)

You can also look at the 2’nd eigenstate (and third and fourth) eigenstate on a 2d plot. This should show clearly the p-character of your state.

c. Three-dimensional plots#

Now we’d like to make some three-dimensional plots.

Make some three-dimensional plots. To do this, you need to put your wave-function back on a grid - i.e. psi[:,i].reshape(80,80,80). Then use your PlotMe function from above. Compare these orbitals to the ones you got by directly plotting the analytic formula.

#answer here