""" This script solves the hydrogen problem using scipy routine
integrate.odeint to integrate Schroedinger differential equation
It coputes the electron density for heaviear atoms in the
non-interacting limit.
"""
from scipy import *
from scipy import integrate, optimize
from pylab import *
def Schrod_deriv(y, r, l, E, Z):
""" Routine which is called by odeint to solve the
Schroedinger equation of the Hydrogen atom.
We are solving the Eq:
u'' + (E -l(l+1)/r^2 + 2Z/r ) * u = 0
which is rewritten as
y == [u, u']
dy/dr == [u', (l(l+1)/r^2 - 2Z/r -E ) u ]
"""
du2 = y[0]*(l*(l+1)/r**2 - 2*Z/r - E)
return [y[1], du2]
def SolveSchroedinger(E,R,l,Z):
""" Integrates Schroedinger equation given a energy E, and angular momentum
eigenvalue l, and charge Z.
For stability, we integrate Schroedinger equation starting at infinity
and integrating down to zero.
We use initial condition at infinity : u(infinity)=0, and u'(infinity)=small.
In order that solution is of the order of unity, we normalize it to (roughly) the peak
value of the function.
"""
R0 = R[::-1]
y0=[0, -1e-4] # [value, derivative]
y = integrate.odeint(Schrod_deriv, y0, R0, args=(l,E,Z) )
norm = max(y[:len(R)*0.9,0])
u = y[::-1,0]/norm
return u
def Shoot(E, R, l, Z):
""" Given a energy E (momentum eigenvalue l, and charge Z), and initial condition
at infinity (u(infinity)=0, u'(infinity)=small) it finds the value of u(r=0).
For the state to be bound state, we have the following condition: u(r=0)==0.
"""
u = SolveSchroedinger(E,R,l,Z)
return u[0] + (u[1]-u[0])*(0.0-R[0])/(R[1]-R[0])
def FindBoundStates(R,lmax,Z,nmax,MaxSteps):
""" For each momentum eigenvalue l=[0,...lmax], it find first nmax-l bound-states.
It first brackets roots of function Shoot using a logarithmic mesh dense at zero.
It later refinds the roots by optimize.brentq method.
"""
Eb=[]
for l in range(lmax+1):
E0 = -1.2*Z**2
dE = 0.1*Z**2
uo = Shoot(E0, R, l, Z)
nfound=0
for i in range(MaxSteps):
E0 += dE
un = Shoot(E0, R, l, Z)
if uo*un < 0.0:
Eb0 = optimize.brentq(Shoot, E0-dE, E0, args=(R, l, Z),xtol=1e-17)
Eb.append([l,Eb0])
nfound+=1
print 'Found bound state at ', Eb0
if (nfound>=nmax-l): break
uo = un
dE /= 1.091
return Eb
def cmpb(x,y):
"""This subroutine is used for sorting bound states.
It sorts them according to eigenvalue, but if two eigenstates are
degenerate, it sorts according to angular momentum l.
"""
if abs(x[1]-y[1])>1e-4:
return cmp(x[1],y[1])
else:
return cmp(x[0],y[0])
R = logspace(-6,2.2,500) # Logarithmic mesh for our integration: [10^-6,...10^2.2]
Z=1 # Charge of the nucleous
nmax = 4 # Maximum principal state
lmax = 3 # Maximum angular momentum
MaxSteps=1000 # Maximum number of steps in searching for the bound state
Ntot = 36 # nuber of all electrons we put into these hydrogen-like states.
# Computes eigenstates requested by lmax and nmax
bound_states = FindBoundStates(R,lmax,Z,nmax,MaxSteps)
# Sorts eigenestates according to their energy
bound_states.sort(cmpb)
N=0 # Total number of electrons in the charge density
rho = zeros(len(R),dtype=float) # Total charge density
for state in bound_states:
l = state[0] # l of current state
Ei = state[1] # Energy of current state
u = SolveSchroedinger(Ei,R,l,Z) # u(r) for this bound state
u2 = u*u # normalizing bound state
norm = integrate.simps(u2,R) # normalization factor
u2 *= 1/norm
dN = 2*(2*l+1) # degeneracy of this particular bound state
if N+dN