src/rdirac.f90

   1
   2 ! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
   3 ! This file is distributed under the terms of the GNU Lesser General Public
   4 ! License. See the file COPYING for license details.
   5
   6 !BOP
   7 ! !ROUTINE: rdirac
   8 ! !INTERFACE:
   9 subroutine rdirac(n,l,k,np,nr,r,vr,eval,g0,f0)
  10 ! !INPUT/OUTPUT PARAMETERS:
  11 !   n    : principal quantum number (in,integer)
  12 !   l    : quantum number l (in,integer)
  13 !   k    : quantum number k (l or l+1) (in,integer)
  14 !   np   : order of predictor-corrector polynomial (in,integer)
  15 !   nr   : number of radial mesh points (in,integer)
  16 !   r    : radial mesh (in,real(nr))
  17 !   vr   : potential on radial mesh (in,real(nr))
  18 !   eval : eigenvalue without rest-mass energy (inout,real)
  19 !   g0   : major component of the radial wavefunction (out,real(nr))
  20 !   f0   : minor component of the radial wavefunction (out,real(nr))
  21 ! !DESCRIPTION:
  22 !   Finds the solution to the radial Dirac equation for a given potential $v(r)$
  23 !   and quantum numbers $n$, $k$ and $l$. The method involves integrating the
  24 !   equation using the predictor-corrector method and adjusting $E$ until the
  25 !   number of nodes in the wavefunction equals $n-l-1$. The calling routine must
  26 !   provide an initial estimate for the eigenvalue. Note that the arrays
  27 !   {\tt g0} and {\tt f0} represent the radial functions multiplied by $r$.
  28 !
  29 ! !REVISION HISTORY:
  30 !   Created September 2002 (JKD)
  31 !EOP
  32 !BOC
  33 implicit none
  34 ! arguments
  35 integer, intent(in) :: n
  36 integer, intent(in) :: l
  37 integer, intent(in) :: k
  38 integer, intent(in) :: np
  39 integer, intent(in) :: nr
  40 real(8), intent(in) :: r(nr)
  41 real(8), intent(in) :: vr(nr)
  42 real(8), intent(inout) :: eval
  43 real(8), intent(out) :: g0(nr)
  44 real(8), intent(out) :: f0(nr)
  45 ! local variables
  46 integer, parameter :: maxit=2000
  47 integer kpa,it,nn,ir,irm,nnd,nndp
  48 ! energy convergence tolerance
  49 real(8), parameter :: eps=1.d-10
  50 real(8) t1,de
  51 ! automatic arrays
  52 real(8) g1(nr),f1(nr),fr(nr),gr(nr),cf(3,nr)
  53 if (k.le.0) then
  54   write(*,*)
  55   write(*,'("Error(rdirac): k <= 0 : ",I8)') k
  56   write(*,*)
  57   stop
  58 end if
  59 if (k.gt.n) then
  60   write(*,*)
  61   write(*,'("Error(rdirac): incompatible n and k : ",2I8)') n,k
  62   write(*,*)
  63   stop
  64 end if
  65 if ((k.eq.n).and.(l.ne.k-1)) then
  66   write(*,*)
  67   write(*,'("Error(rdirac): incompatible n, k and l : ",3I8)') n,k,l
  68   write(*,*)
  69   stop
  70 end if
  71 if (k.eq.l) then
  72   kpa=k
  73 else if (k.eq.l+1) then
  74   kpa=-k
  75 else
  76   write(*,*)
  77   write(*,'("Error(rdirac): incompatible l and k : ",2I8)') l,k
  78   write(*,*)
  79   stop
  80 end if
  81 de=1.d0
  82 nndp=0
  83 do it=1,maxit
  84 ! integrate the Dirac equation
  85   call rdiracdme(0,kpa,eval,np,nr,r,vr,nn,g0,g1,f0,f1)
  86 ! check the number of nodes
  87   nnd=nn-(n-l-1)
  88   if (nnd.gt.0) then
  89     eval=eval-de
  90   else
  91     eval=eval+de
  92   end if
  93   if (it.gt.1) then
  94     if ((nnd.ne.0).or.(nndp.ne.0)) then
  95       if (nnd*nndp.le.0) then
  96         de=de*0.5d0
  97       else
  98         de=de*1.1d0
  99       end if
 100     end if
 101   end if
 102   nndp=nnd
 103   if (de.lt.eps*(abs(eval)+1.d0)) goto 20
 104 end do
 105 write(*,*)
 106 write(*,'("Error(rdirac): maximum iterations exceeded")')
 107 write(*,*)
 108 stop
 109 20 continue
 110 ! find effective infinity and set wavefunction to zero after that point
 111 ! major component
 112 irm=nr
 113 do ir=2,nr
 114   if ((g0(ir-1)*g0(ir).lt.0.d0).or.(g1(ir-1)*g1(ir).lt.0.d0)) irm=ir
 115 end do
 116 g0(irm:nr)=0.d0
 117 ! minor component
 118 irm=nr
 119 do ir=2,nr
 120   if ((f0(ir-1)*f0(ir).lt.0.d0).or.(f1(ir-1)*f1(ir).lt.0.d0)) irm=ir
 121 end do
 122 f0(irm:nr)=0.d0
 123 ! normalise
 124 do ir=1,nr
 125   fr(ir)=g0(ir)**2+f0(ir)**2
 126 end do
 127 call fderiv(-1,nr,r,fr,gr,cf)
 128 t1=sqrt(abs(gr(nr)))
 129 if (t1.gt.0.d0) then
 130   t1=1.d0/t1
 131 else
 132   write(*,*)
 133   write(*,'("Error(rdirac): zero wavefunction")')
 134   write(*,*)
 135   stop
 136 end if
 137 g0(:)=t1*g0(:)
 138 f0(:)=t1*f0(:)
 139 return
 140 end subroutine
 141 !EOC