File: rdirac.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 (145 lines) | stat: -rw-r--r-- 3,676 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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

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

!BOP
! !ROUTINE: rdirac
! !INTERFACE:
subroutine rdirac(sol,n,l,k,nr,r,vr,eval,g0,f0)
! !INPUT/OUTPUT PARAMETERS:
!   sol  : speed of light in atomic units (in,real)
!   n    : principal quantum number (in,integer)
!   l    : quantum number l (in,integer)
!   k    : quantum number k (l or l+1) (in,integer)
!   nr   : number of radial mesh points (in,integer)
!   r    : radial mesh (in,real(nr))
!   vr   : potential on radial mesh (in,real(nr))
!   eval : eigenvalue without rest-mass energy (inout,real)
!   g0   : major component of the radial wavefunction (out,real(nr))
!   f0   : minor component of the radial wavefunction (out,real(nr))
! !DESCRIPTION:
!   Finds the solution to the radial Dirac equation for a given potential $v(r)$
!   and quantum numbers $n$, $k$ and $l$. The method involves integrating the
!   equation using the predictor-corrector method and adjusting $E$ until the
!   number of nodes in the wavefunction equals $n-l-1$. The calling routine must
!   provide an initial estimate for the eigenvalue. Note that the arrays
!   {\tt g0} and {\tt f0} represent the radial functions multiplied by $r$.
!
! !REVISION HISTORY:
!   Created September 2002 (JKD)
!EOP
!BOC
implicit none
! arguments
real(8), intent(in) :: sol
integer, intent(in) :: n,l,k,nr
real(8), intent(in) :: r(nr),vr(nr)
real(8), intent(inout) :: eval
real(8), intent(out) :: g0(nr),f0(nr)
! local variables
integer, parameter :: maxit=2000
integer kpa,it,ir,irm
integer nn,nnd,nndp
! energy convergence tolerance
real(8), parameter :: eps=1.d-12
real(8) t1,de
! automatic arrays
real(8) g1(nr),f1(nr),fr(nr)
! external functions
real(8) splint
external splint
if (k.le.0) then
  write(*,*)
  write(*,'("Error(rdirac): k <= 0 : ",I8)') k
  write(*,*)
  stop
end if
if (k.gt.n) then
  write(*,*)
  write(*,'("Error(rdirac): incompatible n and k : ",2I8)') n,k
  write(*,*)
  stop
end if
if ((k.eq.n).and.(l.ne.k-1)) then
  write(*,*)
  write(*,'("Error(rdirac): incompatible n, k and l : ",3I8)') n,k,l
  write(*,*)
  stop
end if
if (k.eq.l) then
  kpa=k
else if (k.eq.l+1) then
  kpa=-k
else
  write(*,*)
  write(*,'("Error(rdirac): incompatible l and k : ",2I8)') l,k
  write(*,*)
  stop
end if
if (nr.lt.4) then
  write(*,*)
  write(*,'("Error(rdirac): nr < 4 : ",I8)') nr
  write(*,*)
  stop
end if
de=1.d0
nndp=0
do it=1,maxit
! integrate the Dirac equation
  call rdiracint(sol,kpa,eval,nr,r,vr,nn,g0,g1,f0,f1)
! check the number of nodes
  nnd=nn-(n-l-1)
  if (nnd.gt.0) then
    eval=eval-de
  else
    eval=eval+de
  end if
  if (it.gt.1) then
    if ((nnd.ne.0).or.(nndp.ne.0)) then
      if (nnd*nndp.le.0) then
        de=de*0.5d0
      else
        de=de*1.1d0
      end if
    end if
  end if
  nndp=nnd
  if (de.lt.eps*(abs(eval)+1.d0)) goto 20
end do
write(*,*)
write(*,'("Warning(rdirac): maximum iterations exceeded")')
20 continue
! find effective infinity and set wavefunction to zero after that point
! major component
irm=nr
do ir=2,nr
  if ((g0(ir-1)*g0(ir).lt.0.d0).or.(g1(ir-1)*g1(ir).lt.0.d0)) irm=ir
end do
g0(irm:nr)=0.d0
! minor component
irm=nr
do ir=2,nr
  if ((f0(ir-1)*f0(ir).lt.0.d0).or.(f1(ir-1)*f1(ir).lt.0.d0)) irm=ir
end do
f0(irm:nr)=0.d0
! normalise
do ir=1,nr
  fr(ir)=g0(ir)**2+f0(ir)**2
end do
t1=splint(nr,r,fr)
t1=sqrt(abs(t1))
if (t1.gt.0.d0) then
  t1=1.d0/t1
else
  write(*,*)
  write(*,'("Error(rdirac): zero wavefunction")')
  write(*,*)
  stop
end if
call dscal(nr,t1,g0,1)
call dscal(nr,t1,f0,1)
return
end subroutine
!EOC