1
$\begingroup$

I'm trying to solve the Schrödinger equation for the hydrogen atom in the following form numerically:

$$\left[-\frac{\hbar^2}{2m}\frac{d^2}{dr^2}+V(r)+\frac{\hbar^2l(l+1)}{2mr^2}\right]R(r)=ER(r).$$

With the Coulomb-Potential $$V(r) = -\frac{e^2}{4\pi\epsilon_0r^2}.$$ By introducing the dimensionless variables $x = r/a_0$ and $\epsilon = E/E_0$, $E_0=1\text{Ry}$, the Schrödinger equation can be written as follows:

$$\left[-\frac{d^2}{dx^2}-\frac{2}{x}+\frac{l(l+1)}{x^2}\right]R(x)=\epsilon R(x)$$

I used the method of Leandro M. described in this post to solve the Schrödinger equation by discretizing the derivative of $R$ and transforming it into a Matrix equation.

I implemented this solution in Python and I think it worked, the only problem I have is that I don't know how to get from the Eigenvalues and -vectors of the resulting Matrix to the Energies and the corresponding Wavefunctions of the Hydrogen Atom for $n=1,2,3,...$ and $l = 0$.

Edit: To clarify, in the last line of my code i do get some values for the Eigenenergies, they are however in the range of a few 10000 instead of 1 for $n=1$ for example. The problem is that I don't know how to get from the values that I calculated to the correct values.

Here is my code:

import numpy as np
from numpy.linalg import eig
from matplotlib import pyplot as plt

d = 0.001
# set values for r
N = 3000
rmax = 10
r = np.linspace(1e-20, rmax, N)

# create first matrix of schrodinger eq corresponding to the derivative with the following shape:
# (-2  1                   )
# ( 1 -2  1                )
# (    1 -2  1             )  * (-1/(d**2))
# (         ...            )
# (                1  -2  1)
# (                    1 -2)


m = np.empty((N, N))

x = -1 / d ** 2

m[0, 0] = -2 * x
m[0, 1] = 1 * x
m[N - 1, N - 1] = -2 * x
m[N - 1, N - 2] = 1 * x

for i in range(1, N - 1):
    m[i, i - 1] = 1 * x
    m[i, i] = -2 * x
    m[i, i + 1] = 1 * x

for i in range(2, 10):
    m[0, i] = 0

# create matrix corresponding to the potential with the following shape:
# (V1                 )
# (   V2              )
# (      V3           )
# (       ...         )
# (            VN-1   )
# (                 VN)

vdiag = []
l = 0

for i in range(0, N - 1):
    vdiag.append(-2 / r[i] + l * (l + 1) / (r[i] ** 2))

v = np.empty((N, N))
np.fill_diagonal(v, vdiag)

# add matrices to get eigenvalue equation: H*R(x) = E*R(x)
H = np.empty((N, N))

for i in range(0, len(v[0] - 1)):
    for j in range(0, len(v[1] - 1)):
        H[i, j] = m[i, j] + v[i, j]

# setting boundary conditions R_1 = R_N = 0
H[:, 0] = H[:, N - 1] = 0

# determine eigenvalues and eigenvectors
energies, wavefcts = eig(H)
Improve this question
$\endgroup$
8
  • 2
    $\begingroup$ Well, actually eigenvalues are indeed energies and eigenvectors are the wavefunctions or orbitals. $\endgroup$
    – Mithridates the Great
    Nov 18, 2021 at 18:52
  • $\begingroup$ Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. $\endgroup$
    – Community Bot
    Nov 18, 2021 at 18:52
  • $\begingroup$ My Problem is that i don't really know how to get the Eigenenergy for the Quantum Numbers n=1, l=0 for example (which should be approximately 1), since the vector with the eigenvalues only contains values that are way too high. $\endgroup$
    – mysterion123
    Nov 18, 2021 at 19:01
  • $\begingroup$ Your lowest eigenvalue for l=0 should be 0.5. If you don't get close, there's an error in your implementation. Maybe you shouldnt start your r-grid from 1e-20 bur rather some reasonably larger number. Or, you let it start for r=0, but then impose the Dirichlet condition R(r=0)=0. $\endgroup$
    – davidhigh
    Nov 18, 2021 at 20:05
  • $\begingroup$ And wait -- your Rmax, your N and your d are dependent on esch other. You can't simply choose all three separately as you do in your code, but rather calculate the third with the first two. $\endgroup$
    – davidhigh
    Nov 18, 2021 at 23:36

1 Answer 1

Reset to default
2
$\begingroup$

I finally found the solution to my problem, the Eigenenergies where inside the vector of the Eigenvalues, but not in the right order. With the following code I got the right result, if someone is interested:

import numpy as np
from numpy.linalg import eig
from matplotlib import pyplot as plt


# set values for x
xmax = 500
N = 3000
d = xmax/N
x = np.linspace(1e-20, xmax, N)

# create first matrix of schrodinger eq corresponding to the derivative with the following shape:
# (-2  1                   )
# ( 1 -2  1                )
# (    1 -2  1             )  * (-1/(d**2))
# (         ...            )
# (                1  -2  1)
# (                    1 -2)


m = np.zeros((N, N))

g = -1 / d ** 2

m[0, 0] = -2 * g
m[0, 1] = 1 * g
m[N - 1, N - 1] = -2 * g
m[N - 1, N - 2] = 1 * g

for i in range(1, N - 1):
    m[i, i - 1] = 1 * g
    m[i, i] = -2 * g
    m[i, i + 1] = 1 * g

for i in range(2, 10):
    m[0, i] = 0

# create matrix corresponding to the potential with the following shape:
# (V1                 )
# (   V2              )
# (      V3           )
# (       ...         )
# (            VN-1   )
# (                 VN)

vdiag = []
l = 0

for i in range(0, N - 1):
    vdiag.append(-2 / x[i] + l * (l + 1) / (x[i] ** 2))

v = np.zeros((N, N))
np.fill_diagonal(v, vdiag)

# add matrices to get eigenvalue equation: H*R(x) = E*R(x)
H = np.zeros((N, N))

H = m + v



# determine eigenvalues and eigenvectors
energies, wavefcts = eig(H)

energies[0] = energies[N-1] = 0

E = np.sort(energies)
Improve this answer
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.