### HYDROGEN-POTENTIAL

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

!Program name: HYDROGEN-POTENTIAL

!Objective: plot the effective potential energy as a function of radial distance

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

!The Coulomb potential energy equation is:

! V(r) = -KZq2/r = -(Kq2)(Z/r)

! K(Coulomb constant)= 9*109 Nm2C-2

! q = - 1.602 x 10-19C

! Kq2=2.3 x 10-28

! The unit is Joule. We have to change into Rydberg.

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

! The V(centrifugal) equation is:

! Vcent=6.1 x 10-39 l(l+1)/r2 (see text)

! The unit is Joule. We have to change into Rydberg

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

! The unit of the radial distance, r, is meter.

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

Program hydrogen

Implicit none

integer :: mesh, n, l, i, zeta

Real*8 rmin, dr, rmax

double precision :: oldr, e, kq2, h2m

double precision, allocatable :: r(:), r2(:), vpot(:), vc(:), vcent(:), veff(:)

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

zeta = 1

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

! Initial data

rmax= 100. ! Bohr

rmin= 3.E-4 ! Bohr

dr= 0.005

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

mesh = 1000

Allocate (r(mesh))

Allocate (r2(mesh))

Allocate (vc(mesh))

Allocate (vcent(mesh))

Allocate (vpot(mesh))

Allocate (veff(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

Do i=0,mesh

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

r(i)= (exp(oldr))*5.23E-11 !r(i) is changed to SI unit

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

End do

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

! Initialize the potential energy grid

kq2=2.3E-28

h2m=6.1E-39

If (l==0) then

Do i=0,mesh

vc(i)= -(kq2*zeta)/r(i)

veff(i)=vc(i) !veff is in Joule unit

!Let us change vpot unit to Rydberg

vpot(i)=veff(i)*0.458E+18 ! vpot is in Ry

End do

Else if (l > 0) then

Do i=0,mesh

vc(i)= -(kq2*zeta)/r(i)

vcent(i)=((h2m)*(l*(l+1))/r2(i))

veff(i)= vc(i)+vl(i) !veff is in Joule unit

!Let us change vpot unit to Rydberg

vpot(i)=veff(i)*0.458E+18 ! vpot is in Ry

End do

End if

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

!Printing results r(i) and vpot(i) in the 'fileout' (optional)

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

do i=0,mesh

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

End do

write (7,'(/)')

!Close(7)

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

Deallocate (r)

Deallocate (r2)

Deallocate (vc)

Deallocate (vcent)

deallocate (vpot)

Deallocate (veff)

End program