3
$\begingroup$

I have code that attempts to implement a solution to the Schrödinger equation where there is a central potential (more or less im thinking of hydrogen), in 1-D using the numerov method to construct the wave function. It does not take care of the states (N=1,..) as my goal is to understand the link between the theory and the software. Here are some results for some energies, whether or not they are eigenvalues, im not sure:

[6.054999999999915, 6.214999999999912, 6.534999999999905]

and a grapha look at the square of the wave function. The algorithm im trying to implement is noted here in section 6 PDF, but my potential is different. My problem is, the central potential has a 1/r at the origin, how would I deal with that? (I did a sidestep in the code by adding an arbitrary 0.00001 to eliminate 0).

There is a more general code in another question seen here, but I do not follow the shift from equation to 3d matrix.

Here is my code:

import numpy as np
import scipy
import matplotlib.pyplot as plt 

h=np.float64(1)
m=60
c=m/2
estep = 1.0/(100)
delta = 0.0002
x=np.linspace(-1000,1000,num=m)
x=(x - x[m/2]+0.00001 )
v =  1.0/(x**2)
sig=None
E = []

def get_left(wf,k):
    for i in range(2,int(c)):
        wf[i] = 2*(( 1- 1/12* h**2* k[i-2]**2)*wf[i-2] - (1 + 1/12 *h**2* k[i-1]**2)*wf[i-1])/(1 + 1/12* h**2* k[i]**2)
    wf[0:c] = wf[0:c]/(np.sqrt(np.sum(wf[0:c]**2)))
    return wf

def get_right(wf,k):
    mi = m-1
    for i in range(2,int(c)):
        wf[mi-i] = 2*(( 1- 1/12* h**2* k[mi-i+2]**2)*wf[mi-i+2] - (1 + 1/12 *h**2* k[mi-i+1]**2)*wf[mi-i+1])/(1 + 1/12* h**2*k[mi-i]**2)
    wf[c:mi] = wf[c:mi]/(np.sqrt(np.sum(wf[c:mi]**2)))
    return wf

def root_search(epsilonu, epsilonl):
    en = None
    it =0
    while it < 10000:
        phil =  np.zeros(m)
        phir =  np.zeros(m)
        phil[0] =0
        phil[1] = 0.001
        phir[m-1]=0
        phir[m-2]=0.001
        k = (-v + (epsilonu+epsilonl)/2)
        phil = get_left(phil,k)
        phir = get_right(phir,k)
        diffl = np.diff(phil)
        diffr = np.diff(phir)
        erri = diffl[m/2]- diffr[m/2]
        if erri < delta:
            plt.plot(x, (phil+phir)**2 )
            plt.savefig( str('sho'+'.png' ))
            return (epsilonu+epsilonl)/2
        if erri < 0 : 
            epsilonl = (epsilonu+epsilonl)/2
        if erri > 0:
            epsilonu =epsilonl 
            epsilonl = (epsilonu+epsilonl)/2
        it+=1
    return None

def start():
    epsilonl=0.01
    eps = epsilonl
    epsilonu=None
    iteration=0
    sig = None
    while iteration < 100000:
        phil =  np.zeros(m)
        phir =  np.zeros(m)
        phil[0] =0
        phil[1] = 0.001
        phir[m-1]=0
        phir[m-2]=0.001
        k = (-v + eps)
        phil = get_left(phil,k)
        phir = get_right(phir,k)
        diffl = np.diff(phil)
        diffr = np.diff(phir)
        err = diffl[m/2]- diffr[m/2]
        ssig = 1 if err>0 else -1
        if sig:
            if sig != ssig:
                epsilonu= eps 
                E.append(root_search(epsilonu, epsilonl))
                print(E)
        sig = ssig
        epsilonl = eps
        eps = eps+estep
        iteration+=1

if __name__ == "__main__":
    start()

The theory: Im attempting to solve the 1-D Schroedinger equation, $$ \left( -\frac{h^2}{2m}\frac{d^2}{dx^2} + V\left(x\right)\right) \phi\left(x\right) = \epsilon \phi\left( x\right) $$ with V as a central potential (coulombic force, inverse square relation), the final equation looks like (with proper choice of units to give q^2=1, h-1..., energy in hartree) $$ \frac{d^2}{dx^2}\phi\left(x\right) + 2 \left(\epsilon -\frac{1}{x^2} \right) \phi\left(x\right) =0 $$

$\endgroup$
  • 1
    $\begingroup$ If you're thinking about the hydrogen atom, then there shouldn't be a singularity at the origin, because you should do a separation of variables on the 3d PDE, instead of plugging in $V=1/r$ into the 1d PDE, and then there is no singularity at $r=0$ (also, $r\geq0$ then). Otherwise, how to handle a singular potential is more of a mathematics and physics question, rather than numerics: the wavefunction has to vanish so that $V\psi \neq \infty$ is consistent with the equation, and the probability current has to be conserved — that's really a mathematics question. $\endgroup$ – Kirill Dec 31 '15 at 8:44
  • $\begingroup$ See, for example physics.stackexchange.com/q/69889/3696 $\endgroup$ – Kirill Dec 31 '15 at 8:48
  • $\begingroup$ I don't think it's only a mathematics question -- properly handling functions with singularities is also a computational challenge. If possible, you want to rewrite the equation to remove the singularity (as you suggested) or otherwise take an analytic limit, I guess. $\endgroup$ – Jannis Teunissen Dec 31 '15 at 15:27
  • $\begingroup$ I don't know what you mean by: but I do not follow the shift from equation to 3d matrix. If you mean how to translate the differential equation into a system of equations, it was done using Finite Differences. $\endgroup$ – nicoguaro Jan 2 '16 at 21:08
2
$\begingroup$

As mentioned by @Kirill, the differential equation for the hydrogen is slightly different. After the separation of variables, you end up with

$$\left[ - \frac{\hbar^2}{2\mu} \left({1 \over r^2}{\partial \over \partial r}\left(r^2 {\partial R(r)\over \partial r}\right) - {l(l+1)R(r)\over r^2} \right) + V(r)R(r) \right]= E R(r),$$

You can rewrite the equation as

$$\left[ -\frac{1}{2} \frac{d^2}{dy^2} + \frac{1}{2} \frac{l(l+1)}{y^2} - \frac{1}{y}\right] u_l = W u_l$$

The next code (based on my previous post) computes the eigenvalues/eigenvectors for this equation. See that I added a small value to the denominators, so they don't vanish at zero.

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

L = 50*np.pi
N = 1001
xa = 0
xb = L
x = np.linspace(xa, xb, N)
dx = x[1] - x[0]

l = 1
k = 0
T = -0.5*diags([-2., 1., 1.], [0, -1, 1], shape=(N, N))/dx**2
U_vec = 0.5*l*(l + 1)/(x**2 + 1e-6) - 1/(np.abs(x) + 1e-6)
U = diags([U_vec], [0])

H = T + U

vals, vecs = eigsh(H, which='SA')

print(np.round(vals, 6))
print(np.round([-1/(2*n**2) for n in range(k + l + 1, k + l + 7)], 6))

for k in range(5):
    vec = vecs[:, k]
    mag = np.sqrt(np.dot(vecs[:, k],vecs[:, k]))
    vec = vec/mag
    plt.plot(x, vec, label=r"$n=%i$"% (k+1))

plt.xlabel(r"$x$")
plt.ylabel(r"$\psi(x)$")
plt.xlim(xa, xb)
plt.legend()
plt.savefig("eigenvectors.png", dpi=600) 
plt.show()

That gives as solution:

[-0.125106 -0.0556   -0.031272 -0.020012 -0.013896 -0.010209]
[-0.125    -0.055556 -0.03125  -0.02     -0.013889 -0.010204]

enter image description here

For your repulsive potential (different sign than the Laplacian), I used the same code

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

L = 50*np.pi
N = 1001
xa = 0
xb = L
x = np.linspace(xa, xb, N)
dx = x[1] - x[0]

T = -0.5*diags([-2., 1., 1.], [0, -1, 1], shape=(N, N))/dx**2
U_vec = 1/(x**2 + 1e-6)
U = diags([U_vec], [0])

H = T + U

vals, vecs = eigsh(H, which='SA')

print(np.round(vals, 6))

for k in range(5):
    vec = vecs[:, k]
    mag = np.sqrt(np.dot(vecs[:, k],vecs[:, k]))
    vec = vec/mag
    plt.plot(x, vec, label=r"$n=%i$"% (k+1))

plt.xlabel(r"$x$")
plt.ylabel(r"$\psi(x)$")
plt.xlim(xa, xb)
plt.legend()
plt.savefig("eigenvecs.png", dpi=600)  
plt.show()

with answers

[ 0.000408  0.001207  0.002405  0.004001  0.005997  0.008392]

enter image description here

Although, I am not sure what that means. Because it is not letting the electron to be in the origin pushing it away, and the domain should be really unbounded. It is similar to having a square potential with a round edge at the left. For the more realistic case, you can use nonuniform finite differences and increase them exponentially or absorbing boundary conditions.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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