Issues with eigenvalues and spatial distribution of 3D hydrogen atom

[Ionizing]

Hi,
I just found this wonderful repo and it helped me a lot in understanding the spatial distribution of each level of the hydrogen atom.

But when I ran the script 3D_hydrogen_atom.py, I found two issues with this program:

1. The eigenvalues of it returned seemed to be far from the known values:

The ground state should be -13.6eV, but the script returned -8.26297717eV. Also, I found the grid is 303030 with the actual length of each edge are also 30A30A30A. Should I make a finner grid to solve for a better result?

2. The axis unit is much larger than the known radial distribution. According to the radial distribution of hydrogen 1s state. There is a peak below the Bohr radius (0.53A). But the extent of of the result is much larger

However, if I change the `extent` parameter to 5A, the issue still remains:

Though I came with such issues, I must say thank you for giving a feasible numeric approach to solve the 3D single particle schrodinger equation and giving fancy results to help students understanding the QM.

Looking forward to your reply!

[rafael-fuente]

The problem is the spatial discretization, as you said 30x30x30 is very small, especially considering that the Coulomb potential has a singularity. The problem is fixed using a smoother discretization, but If you try a discretization smoother than (50x50x50) you will soon find that a personal computer cannot handle it.

But I already fixed the problem!

It happens that LOBPCG algorithm can deal with very large grids (200x200x200) and improve the accuracy of the solver considerably. I'm working on implementing it. I will update the repository with the implementation this weekend.

Using a 200x200x200 grid these are energies:

[-13.52, -3.384, -1.502, -0.826, -0.531] eV

while the analytical ones are:

[13.605 , 3.401 , 1.512, 0.850, 0.544] eV

This is the best you will be able to get with a 3D simulation. ( the computation last about 100 seconds on a good pc)

But the hydrogen atom has a clearly azimuthal symmetry and the problem can be reduced to 2D.
I'm currently working on an azimuthal solver too. This is what I get using a 700x700 grid:

[-13.58, -3.394, -1.504, -0.848, -0.543] eV

[rafael-fuente]

Yes, singularity is a big problem. We need a large grid to process accurately it. You may find useful this visualization of the radial distribution of hydrogen ground state, in which I compared the results of a simulation with extent = 16, and N = 200 with the analytical expression:

As you can see it compares quite well with the analytical value. However, there is a little disagreement near the singularity. This doesn't happen for the other orbitals, because its wavefunction distribution is further from the singularity.

Here also a 3d density plot:

I want to ask that if the spatial distribution of eigenstates is zoomed in the visual window because I saw that n=2,3 states can spread to 5 angstroms or more.

It's completely fine, they grow up pretty fast. Exactly, the mean growth of the radius of the hydrogen orbitals grows proportionally to the square of the quantum number n. Here an exact formula that illustrates this fact:

And if I choose extent = 8, will it have an impact on the spatial distribution (e.g. the eigenstates are squeezed)?

Yes. You wouldn't be simulating a hydrogen atom, but a hydrogen atom on a box. The extent must be always larger than the size of the wave function.

I just re-implemented the solving part of 3D_hydrogen_atom in julia, see here. A vast performance improvement was observed (It took only 32 secs to calculate a 50*50*50 system using eigs instead of lobpcg on a laptop). I am not sure if there are still some bugs in it, but the eigen values seemed to be reasonable.

That's pretty nice! Definitely, a port of the solver to Julia may worth it.

[Ionizing]

Okay, it seems all of my problems are already solved. Thank you very much!

[Ionizing]

Could I ask one more question? How to construct the kinetic matrix in Hamiltonian when Nx != Ny != Nz? I read the code and found:

def get_kinetic_matrix(self, H):  # Single particle kinetic matrix

    I = eye(H.N)
    T_ =  diags([-2., 1., 1.], [0,-1, 1] , shape=(H.N, H.N))*-k/(self.m*H.dx**2)
    if H.spatial_ndim ==1:
        T = T_

    elif H.spatial_ndim ==2:
        T =  kron(T_,I) + kron(I,T_)

    elif H.spatial_ndim ==3:
        T =  kron(T_, kron(I, I)) + kron(I, kron(T_, I)) + kron(I, kron(I, T_))

    return T

In my view, I should constuct not only the kinetic matrices for each axis: Tx_, Ty_, Tz_, but also the Identity matrices Ix, Iy, Iz. Then the total kinetic matrix should be T = kron(Tx_, kron(Iy, Iz)) + kron(Ix, kron(Ty_, Iz)) + kron(Ix, kron(Iy, Tz_)). Is there anything wrong with my derivation? Are there any papers or books to help understand the math?

[rafael-fuente]

Could I ask one more question?

Of course! :)

In my view, I should constuct not only the kinetic matrices for each axis: Tx_, Ty_, Tz_, but also the Identity matrices Ix, Iy, Iz. Then the total kinetic matrix should be T = kron(Tx_, kron(Iy, Iz)) + kron(Ix, kron(Ty_, Iz)) + kron(Ix, kron(Iy, Tz_)). Is there anything wrong with my derivation? Are there any papers or books to help understand the math?

We didn't implement it, but different spatial discretization for each axis may be useful for some systems. Your guess is completely right.
Here is some introduction to the math:

Discretizing the Hamiltonian:

• Discretization of 1D continuous Hamiltonian

Using 2D Grids for higher dimensional discrete operators:

• Kronecker sum of discrete Laplacians

[Ionizing]

Thanks for the reply! The supplementary links are very helpful!