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) = -KZq²/r = -(Kq²)(Z/r)

! K(Coulomb constant)= 9*10⁹ Nm²C^-2

! q = - 1.602 x 10^-19C

! Kq²=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)/r² (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