HYDROGEN-RADIAL

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

!Program name: HYDROGEN-RADIAL

!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'

Read (*,*) n, l

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

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

! Initialize constant radial grid

 

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'

 

open(7,file='radial.dat',status='replace')

  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