anchor.diff:  !   Created September 2004 (JKD)
   !EOP
   !BOC
  -implicit none
  -! arguments
  -character(*), intent(in) :: findlinentype
  -integer, intent(in) :: l
  -integer, intent(in) :: k
  -integer, intent(in) :: np
  -integer, intent(in) :: nr
  -real(8), intent(in) :: r(nr)
  -real(8), intent(in) :: vr(nr)
  -real(8), intent(in) :: de0
  -real(8), intent(inout) :: e
  -! local variables
  -! maximum number of steps
  -integer, parameter :: maxstp=1000
  -integer::ie, nn
  -! energy search tolerance
  -real(8), parameter :: eps=1.d-5
  -real(8)::de, et, eb, t, tp
  -! automatic arrays
  -real(8)::p0(nr), p1(nr), q0(nr), q1(nr)
  -
  -
  -character*4, parameter ::        EMAIN='CONT'
  -
  -real(8)   DDNOLD, DNOLD, DU, DUDN, DUPOLD, DUTST, DUUP
  -real(8)   E1, E2, EIDN, EIUP, U, UDN, UPOLD, UTST
  -real(8)   UUP
  -
  -select case(trim(findlinentype))
  -case('simple')
  -
  -tp=0.d0
  -! find the top of the band
  -de=abs(de0)
  -et=e
  -do ie=1, maxstp
  -  et=et+de
  -  call rschroddme(0, l, k, et, np, nr, r, vr, nn, p0, p1, q0, q1)
  -  t=p0(nr)
  -  if (ie.gt.1) then
  -    if (t*tp.le.0.d0) then
  -      if (abs(de).lt.eps) goto 10
  -      de=-0.5d0*de
  -    end if
  -  end if
  -  tp=t
  -end do
  -goto 30
  -10 continue
  -! find the bottom of the band
  -de=-abs(de0)
  -eb=et+5.d0*abs(de0)
  -do ie=1, maxstp
  -  eb=eb+de
  -  call rschroddme(0, l, k, eb, np, nr, r, vr, nn, p0, p1, q0, q1)
  -  t=p1(nr)
  -  if (ie.gt.1) then
  -    if (t*tp.le.0.d0) then
  -      if (abs(de).lt.eps) goto 20
  -      de=-0.5d0*de
  -    end if
  -  end if
  -  tp=t
  -end do
  -goto 30
  -20 continue
  -! set the band energy to the mid-point
  -e=(et+eb)/2.d0
  -30 continue
  -
  -case('advanced')
  - de=de0
  -      call rschroddme(0,l,k,E,np,nr,r,vr,nn,p0,p1,q0,q1)
  -      U=p0(nr)
  -      DU=p1(nr)
  -      EIUP = E
  -      EIDN = E
  -      E1 = -100.0D+0
  -      E2 = -100.0D+0
  -      UPOLD = U
  -      DUPOLD = DU
  -      DNOLD = U
  -      DDNOLD = DU
  -   101 CONTINUE
  -      EIUP = EIUP + DE
  -      call rschroddme(0,l,k,EIUP,np,nr,r,vr,nn,p0,p1,q0,q1)
  -      UUP=p0(nr)
  -      DUUP=p1(nr)
  -      UTST = UPOLD*UUP
  -      DUTST = DUPOLD*DUUP
  -      UPOLD = UUP
  -      DUPOLD = DUUP
  -      IF (UTST .LT. 0.0D+0) THEN
  -         E2 = EIUP
  -         IF (E1 .GT. -30.0D+0) GOTO 301
  -      ENDIF
  -      IF (DUTST .LT. 0.0D+0) THEN
  -         E1 = EIUP
  -         IF (E2 .GT. -30.0D+0) GOTO 301
  -      ENDIF
  -      IF ((E1 .LT. -30.0D+0) .AND. (EIUP .LT. 2.5D+0)) THEN
  -   201    CONTINUE
  -         EIDN = EIDN - DE
  -         call rschroddme(0,l,k,EIDN,np,nr,r,vr,nn,p0,p1,q0,q1)
  -         UDN=p0(nr)
  -         DUDN=p1(nr)
  -         UTST = DNOLD*UDN
  -         DUTST = DDNOLD*DUDN
  -         DNOLD = UDN
  -         DDNOLD = DUDN
  -         IF (UTST .LT. 0.0D+0) THEN
  -            E2 = EIDN
  -            IF (E1 .GT. -30.0D+0) GOTO 301
  -         ENDIF
  -         IF (DUTST .LT. 0.) THEN
  -            E1 = EIDN
  -            IF (E2 .GT. -30.0D+0) GOTO 301
  -         ENDIF
  -         IF (E2 .LT. -30.0D+0) THEN
  -            GOTO 101
  -         ELSEIF (EIDN .GT. (E - 3.0D+0)) THEN
  -            GOTO 201
  -         ENDIF
  -      ELSEIF (EIUP .LT. 0.5D+0) THEN
  -         GOTO 101
  -      ENDIF
  -   301 CONTINUE
  -      IF ((E1 .LT. -30.0D+0) .AND. (E2 .LT. -30.0D+0)) THEN
  -         IF (EMAIN .EQ. 'STOP') THEN
  -            GOTO 900
  -         ELSE
  -            E = 1.0D+0
  -            write(*,*) 'l,ei,e1,e2', l, E, E1, E2
  -         ENDIF
  -      ELSEIF (E2 .LT. -30.0D+0) THEN
  -         IF (EMAIN .EQ. 'STOP') THEN
  -            GOTO 900
  -         ELSE
  -            E = MAX(E1,E)
  -            write(*,*) 'l,ei,e1,e2', l, E, E1, E2
  -         ENDIF
  -      ELSEIF (E1 .LT. -30.0D+0) THEN
  -         GOTO 900
  -      ELSE
  -         E = (E1+E2)*0.5D+0
  -         write(*,*) 'l,ei,e1,e2', l, E, E1, E2
  -      ENDIF
  -      RETURN
  -  900 continue
  -      write(*,'("no energy limits found for L=",i2)') l
  -      write(*,'("E-bottom ",g18.10,3x,"E-top ",g18.10)') e1,e2
  -
  -case default
  -  write(*,*)
  -  write(*,'("Error(findband): No such method for search of linearization energies: ",a)') &
  -    trim(findlinentype)
  -  write(*,*)
  -  stop
  -end select
  -
  -end subroutine
  +  implicit none
  +  ! arguments
  +  character(*), intent(in) :: findlinentype
  +  integer, intent(in) :: l
  +  integer, intent(in) :: k
  +  integer, intent(in) :: np
  +  integer, intent(in) :: nr
  +  real(8), intent(in) :: r(nr)
  +  real(8), intent(in) :: vr(nr)
  +  real(8), intent(in) :: de0
  +  real(8), intent(inout) :: e
  +  ! local variables
  +  ! maximum number of steps
  +  integer, parameter :: maxstp=1000
  +  character*4, parameter :: emain='continue'
  +  integer::ie, nn
  +  ! energy search tolerance
  +  real(8), parameter :: eps=1.d-5
  +  real(8)::de, et, eb, t, tp
  +  real(8)::u, uup, udn, upold, ddnold, dnold, du, dudn, dupold, duup, utst, dutst
  +  real(8)::e1, e2, eidn, eiup
  +  ! automatic arrays
  +  real(8)::p0(nr), p1(nr), q0(nr), q1(nr)
  +  select case(trim(findlinentype))
  +  case('simple')
  +     tp=0.d0
  +     ! find the top of the band
  +     de=abs(de0)
  +     et=e
  +     do ie=1, maxstp
  +        et=et+de
  +        call rschroddme(0, l, k, et, np, nr, r, vr, nn, p0, p1, q0, q1)
  +        t=p0(nr)
  +        if (ie.gt.1) then
  +           if (t*tp.le.0.d0) then
  +              if (abs(de).lt.eps) goto 10
  +              de=-0.5d0*de
  +           end if
  +        end if
  +        tp=t
  +     end do
  +     goto 30
  +10   continue
  +     ! find the bottom of the band
  +     de=-abs(de0)
  +     eb=et+5.d0*abs(de0)
  +     do ie=1, maxstp
  +        eb=eb+de
  +        call rschroddme(0, l, k, eb, np, nr, r, vr, nn, p0, p1, q0, q1)
  +        t=p1(nr)
  +        if (ie.gt.1) then
  +           if (t*tp.le.0.d0) then
  +              if (abs(de).lt.eps) goto 20
  +              de=-0.5d0*de
  +           end if
  +        end if
  +        tp=t
  +     end do
  +     goto 30
  +20   continue
  +     ! set the band energy to the mid-point
  +     e=(et+eb)/2.d0
  +30   continue
  +  case('advanced')
  +     de=de0
  +     call rschroddme(0,l,k,e,np,nr,r,vr,nn,p0,p1,q0,q1)
  +     u=p0(nr)
  +     du=p1(nr)
  +     eiup=e
  +     eidn=e
  +     e1=-1000.0d0
  +     e2=-1000.0d0
  +     upold=u
  +     dupold=du
  +     dnold=u
  +     ddnold=du
  +101  continue
  +     eiup=eiup+de
  +     call rschroddme(0,l,k,eiup,np,nr,r,vr,nn,p0,p1,q0,q1)
  +     uup=p0(nr)
  +     duup=p1(nr)
  +     utst=upold*uup
  +     dutst=dupold*duup
  +     upold=uup
  +     dupold=duup
  +     if (utst.lt.0.d0) then
  +        e2=eiup
  +        if (e1.gt.-30.d0) goto 301
  +     end if
  +     if (dutst.lt.0.d0) then
  +        e1=eiup
  +        if (e2.gt.-30.d0) goto 301
  +     end if
  +     if ((e1.lt.-30.d0).and.(eiup.lt.2.5d0)) then
  +201     continue
  +        eidn=eidn-de
  +        call rschroddme(0,l,k,eidn,np,nr,r,vr,nn,p0,p1,q0,q1)
  +        udn=p0(nr)
  +        dudn=p1(nr)
  +        utst=dnold*udn
  +        dutst=ddnold*dudn
  +        dnold=udn
  +        ddnold=dudn
  +        if (utst.lt.0.d0) then
  +           e2=eidn
  +           if (e1 .gt. -30.d0) goto 301
  +        end if
  +        if (dutst .lt. 0.d0) then
  +           e1=eidn
  +           if (e2 .gt. -30.d0) goto 301
  +        end if
  +        if (e2 .lt. -30.d0) then
  +           goto 101
  +        else if (eidn .gt. (e-3.d0)) then
  +           goto 201
  +        end if
  +     else if (eiup .lt. 0.5d0) then
  +        goto 101
  +     end if
  +301  continue
  +     if ((e1 .lt. -30.d0) .and. (e2 .lt. -30.d0)) then
  +        if (emain .eq. 'stop') then
  +           goto 900
  +        else
  +           e=1.d0
  +           write(*,*) 'l,ei,e1,e2', l, e, e1, e2
  +        end if
  +     else if (e2 .lt. -30.d0) then
  +        if (emain .eq. 'stop') then
  +           goto 900
  +        else
  +           e=max(e1,e)
  +           write(*,*) 'l,ei,e1,e2', l, e, e1, e2
  +        end if
  +     else if (e1 .lt. -30.d0) then
  +        goto 900
  +     else
  +        e=(e1+e2)*0.5d0
  +        write(*,*) 'l,ei,e1,e2', l, e, e1, e2
  +     end if
  +     return
  +900  continue
  +     write(*,*)
  +     write(*,'("Error(findband): no energy limits found for l=",i2)') l
  +     write(*,'("E-bottom ",g18.10,4x,"E-top ",g18.10)') e1,e2
  +     stop
  +  case default
  +     write(*,*)
  +     write(*,'("Error(findband): No such method for search of linearization energies: ",a)') &
  +          trim(findlinentype)
  +     write(*,*)
  +     stop
  +  end select
  +end subroutine findband
   !EOC
  
positive.code: 
  
  
  ! 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
  
negative.code: !BOP
  ! !ROUTINE: findband
  ! !INTERFACE:
  
  subroutine findband(findlinentype,l, k, np, nr, r, vr, de0, e)
  ! !INPUT/OUTPUT PARAMETERS:
  