File: potnucl.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (64 lines) | stat: -rw-r--r-- 1,637 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

! 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
integer ir
real(8) rn,t1,t2
! external functions
real(8) radnucl
external radnucl
if (zn.eq.0.d0) then
  vn(:)=0.d0
  return
end if
if (ptnucl) then
! nucleus is taken to be a point particle
  vn(:)=zn/r(:)
else
! approximate nuclear radius
  rn=radnucl(zn)
  t1=zn/(2.d0*rn**3)
  t2=3.d0*rn**2
  do ir=1,nr
    if (r(ir).lt.rn) then
      vn(ir)=t1*(t2-r(ir)**2)
    else
      vn(ir)=zn/r(ir)
    end if
  end do
end if
return
end subroutine
!EOC