potnucl.f90

Go to the documentation of this file.

! Copyright (C) 2008 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: potnucl
! !INTERFACE:
subroutine potnucl(ptnucl,nr,r,zn,vn)
! !INPUT/OUTPUT PARAMETERS:
!   ptnucl : .true. if the nucleus is a point charge (in,logical)
!   nr     : number of radial mesh points (in,integer)
!   r      : radial mesh (in,real(nr))
!   zn     : nuclear charge (in,real)
!   vn     : potential on radial mesh (out,real(nr))
! !DESCRIPTION:
!   Computes the nuclear Coulomb potential on a radial mesh. The nuclear radius
!   $R$ is estimated from the nuclear charge $Z$ and the potential is given by
!   $$ V(r)=\begin{cases}
!    Z(3R^2-r^2)/2R^3 & r<R \\
!    Z/r & r\ge R\end{cases} $$
!   assuming that the nucleus is a uniformly charged sphere. If {\tt ptnucl} is
!   {\tt .true.} then the nucleus is treated as a point particle.
!
! !REVISION HISTORY:
!   Created January 2009 (JKD)
!EOP
!BOC
implicit none
! arguments
logical, intent(in) :: ptnucl
integer, intent(in) :: nr
real(8), intent(in) :: r(nr),zn
real(8), intent(out) :: vn(nr)
! local variables
real(8) rn,t1,t2
! external functions
real(8), external :: radnucl
if (zn == 0.d0) then
  vn(1:nr)=0.d0
  return
end if
if (ptnucl) then
! nucleus is taken to be a point particle
  vn(1:nr)=zn/r(1:nr)
else
! approximate nuclear radius
  rn=radnucl(zn)
  t1=zn/(2.d0*rn**3)
  t2=3.d0*rn**2
  where (r(:) < rn)
    vn(:)=t1*(t2-r(:)**2)
  elsewhere
    vn(:)=zn/r(:)
  end where
end if
end subroutine
!EOC