!---------------------------------------------------------------

!Radial equation of the hydrogen-like atom

!Objective: plot the radial wave function for each n,l pair

! The equation is atomic units

!---------------------------------------------------------------

!The Coulomb potential energy equation is:

! V(r) = -2Z/r

! The ceentrifugal potential equation is:

! VL=l(l+1)/r^2

! The energy is:

! E=-2*Z/n^2

!---------------------------------------------------------------

! Use of Numerov method to find the wave function

! Function fn from Numerov equation is

! Fn = 1 + K2*h^2/12 , where K2=-2(Veff - E)

! In the code, h is changed into dr

!---------------------------------------------------------------

Program hydrogen

Implicit none

integer :: mesh, n, l, i, j, zeta, m, maxiter=100

Real*8 rmin, dr, zmesh, y_out_m

double precision :: e, fh12, norm, eps=1.0E-6

double precision, allocatable :: r(:), r2(:), y(:), vc(:), vcent(:), &

veff(:), k2(:), fn(:), f(:)

!----------------------------------------------------

zeta = 1

!------------------------------------------------------

! Initial data

zmesh= zeta !zmesh is zeta as a real number

rmin= 1.E-5

dr= 0.005d0

!----------------------------------------------------------------

! Number of points for the calculation (mesh)=1000 for Z=1

mesh = 1000

Allocate (r(mesh))

Allocate (r2(mesh))

Allocate (y(mesh))

Allocate (vc(mesh))

Allocate (f(mesh))

Allocate (fn(mesh))

Allocate (vcent(mesh))

Allocate (veff(mesh))

Allocate (k2(mesh))

!----------------------------------------------------------------

!User entering data

Write (*,*) 'Enter n and l. n > 0 and l=n-1'

If (n<1) then

Write (*,*) 'n < 1. Unphysical.'

Stop

Else if (n < l+1) then

Write (*,*) 'n < l+1. It has to be n=l+1!'

Stop

Else if (l < 0) then

Write (*,*) 'l < 0. Unphysical.'

Stop

End if

!-----------------------------------------------------------------

! Initialize exponential radial grid only for the potential

! This grid does not work for 1s and 2p states

Do i=0,mesh

oldr = rmin + dr*i !oldr is in Bohr unit

r(i)= (exp(oldr))

r2(i)= (r(i) * r(i))

End do

!-----------------------------------------------------

! Initialize the potential energy grid for the y(i)

Do i=0,mesh

vc(i)= -2*zmesh/r(i)

vcent(i)=(l*(l+1.d0))/(r2(i))

veff(i)= vc(i)+vcent(i)

End do

!----------------------------------------------------------------------

dr= 0.02d0

Do i=0,mesh

r(i) = rmin + dr*i !oldr is in Bohr unit

r2(i)= (r(i) * r(i))

End do

!----------------------------------------------------------------------

!Set the  eigenvalue

e= -zmesh/n**2

!-----------------------------------------------------------------------

!calculate part of Numerov function f, fh12 and K2

fh12=dr*dr/12.0d0

k2= vc-e+vcent

f(0)= fh12*k2(0)

!-----------------------------------------------------------------

!Start Numerov integration

!----------------------------------------------------------------

!set boundary conditions

y(0)=r(0)**(l+1)

y(1)=r(1)**(l+1)

y(mesh)=dr

m=0

!---------------------------------------------------------------

!find the matching point m checking the change of sign of f

do i=1,mesh

f(i)= fh12*k2(i)

if (f(i) == 0.d0) f(i)=1.E-20

if (f(i) /= sign(f(i),f(i-1))) m=i

end do

!--------------------------------------------------------------

! With the values of f, obtain fn

!obtain the y(mesh-1) for inward integration

fn=1.d0-f

y(mesh-1)=y(mesh)*(12.d0-10.d0*fn(mesh))/fn(mesh-1)

!--------------------------------------------------------------

!Start Numerov outward integration

do i=1,m-1

y(i+1)=y(i)*(12.d0-10.d0*fn(i))-(fn(i-1)*y(i-1))/fn(i+1)

if (y(i) == 0.0d0) y(i)=1.d-20

end do

y_out_m = y(m)

!print *, y(m)

!------------------------------------------------------------------

! Start Numerov inward integration

do i = mesh-1,m+1,-1

y(i-1)=y(i)*(12.d0-10.d0*fn(i))-(fn(i+1)*y(i+1))/fn(i-1)

if (y(i-1) > 1.E10) then

do j=mesh,i-1,-1

y(j)=y(j)/y(i-1)

end do

end if

end do

!------------------------------------------------------------------

! rescale function to match at the turning point

y(:m-1) = y(:m-1)/y_out_m

y(m:) = y(m:)/y(m)

!----------------------------------------------------------------

!Normalization process

norm = sqrt(dot_product(y,y))

y=y/norm

!----------------------------------------------------------------

!Print the result of the wave function

do i=0,mesh

print *, r(i),y(i)

End do

!-------------------------------------------------------------

!Printing results r(i) and y(i) in the 'radial'

do i=0,mesh

write (7,'(3e16.8,f12.6,3e16.8,f12.6)') r(i),y(i)

End do

write (7,'(/)')

Close(7)

!--------------------------------------------------------------------------

Deallocate (r)

Deallocate (r2)

Deallocate (vc)

Deallocate (vcent)

Deallocate (veff)

Deallocate (k2)

Deallocate (y)

Deallocate (f)

Deallocate (fn)

End program