1!
2! Copyright (C) 2004 PWSCF group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!--------------------------------------------------------------------------
9subroutine find_qi(logderae,xc,ik,lam,ncn,flag,iok)
10  !--------------------------------------------------------------------------
11  !
12  !      This routine finds three values of q such that the
13  !      functions f_l have a logarithmic derivative equal to
14  !      logderae at the point ik
15  !
16  !      f_l=j_l(r) * r**flag
17  !
18  use kinds, only:dp
19  use ld1inc, only: grid
20  implicit none
21
22  integer,parameter ::  ncmax=10   ! maximum allowed nc
23
24  integer ::      &
25       ik,    & ! input: the point corresponding to rcut
26       ncn,   & ! input: the number of qi to compute
27       flag,  & ! input: the type of function
28       iok,   & ! output: if 0 the calculation in this routine is ok
29       lam      ! input: the angular momentum
30
31  real(DP) :: &
32       xc(ncn),& ! output: the values of qi
33       logderae  ! input: the logarithmic derivative
34
35
36  real(DP) ::   &
37       j1(ncmax),& ! the bessel function in three points
38       qmax,qmin,& ! the limits of the q search
39       logdermax,logdermin,& ! the maximum and minimum logder
40       logder, & ! the actual logder
41       jlmin, jlmax, & ! the value of jl in qmin and qmax
42       compute_log, &! function for log derivative
43       dq, dq_0 ! the step to braket the q
44
45  integer ::    &
46       nc,  &    ! counter on the q found
47       icount, &  ! too many iterations
48       icount1, & ! too many iterations
49       imax,&   ! maximum number of iteration to braket
50       iq      ! counter on iteration
51
52  iok=0
53  if (ncn.gt.ncmax) &
54       call errore('find_qi','ncn is too large',1)
55
56  if (flag.eq.0.and.lam.ne.0) &
57       call errore('find_qi','lam too large for this iflag',1)
58
59  if (lam.gt.6) &
60       call errore('find_qi','l not programmed',1)
61  !
62  !    fix deltaq and the maximum step number
63  !
64  dq_0=0.05_dp
65  imax=600
66  !
67  !    prepare for the first iteration. For too small q the function could
68  !    have noise.
69  !
70  qmax=0.5_dp
71  call sph_bes(7,grid%r(ik-3),qmax,lam,j1)
72  j1(1:7) = j1(1:7)*grid%r(ik-3:ik+3)**flag
73  logdermax=compute_log(j1,grid%r(ik),grid%dx)-logderae
74  jlmax = j1(4)
75
76  icount=0
77  do nc=1,ncn
78     !
79     !    bracket the zero
80     !
81     icount1=0
82200  dq=dq_0
83     qmin=qmax
84     logdermin=logdermax
85     jlmin = jlmax
86     do iq=1,imax
87        qmax=qmin+dq
88        call sph_bes(7,grid%r(ik-3),qmax,lam,j1)
89        j1(1:7) = j1(1:7)*grid%r(ik-3:ik+3)**flag
90        logdermax=compute_log(j1,grid%r(ik),grid%dx)-logderae
91        jlmax = j1(4)
92        !
93        !    the zero has been bracketed?
94        !
95        if ( jlmin * jlmax .gt. 0.0_dp ) then ! far from an asintote
96           !
97           if (logdermax*logdermin.lt.0.0_dp) goto 100
98           !
99           qmin = qmax
100           logdermin=logdermax
101           jlmin = jlmax
102        else
103           if (logdermax*logdermin.lt.0.0_dp) then
104              qmin = qmax
105              logdermin=logdermax
106              jlmin = jlmax
107           else
108              dq = 0.5_dp * dq
109           end if
110        end if
111     enddo
112     call infomsg ('find_qi','qmax not found ')
113     iok=1
114     return
115100  continue
116     !
117     !      start bisection loop
118     !
119     xc(nc)=qmin+(qmax-qmin)/2.0_dp
120     call sph_bes(7,grid%r(ik-3),xc(nc),lam,j1)
121     j1(1:7) = j1(1:7)*grid%r(ik-3:ik+3)**flag
122     logder=compute_log(j1,grid%r(ik),grid%dx)-logderae
123     if (logder*logdermin.lt.0.0_dp) then
124        qmax=xc(nc)
125        logdermax=logder
126     else
127        qmin=xc(nc)
128        logdermin=logder
129     endif
130     !
131     !    avoid the asintotes
132     !
133     if (abs(logdermin-logdermax).gt.1.e3_dp) then
134        qmax=xc(nc)
135        logdermax=logder
136        icount1=icount1+1
137        if (icount1<20) then
138           goto 200
139        else
140           call errore('find_q','problem finding q',1)
141        endif
142     endif
143     !
144     !    check for convergence
145     !
146     icount=icount+1
147     if (icount>1000) call errore('find_q','too many iterations',1)
148     if (abs(logdermax-logdermin).gt.1.e-8_dp) goto 100
149  enddo
150
151  return
152end subroutine find_qi
153
154function compute_log(j1,rj,dx)
155  use kinds, only : DP
156  implicit none
157  real(DP) ::   &
158       compute_log, &
159       deriv_7pts,  &
160       dx,          &
161       j1(7),       &
162       rj
163
164  !      compute_log=(j1(2)-j1(1))*2.0_dp/( (rj(2)-rj(1))*(j1(2)+j1(1)) )
165  compute_log=deriv_7pts(j1,4,rj,dx)/j1(4)
166
167  return
168end function compute_log
169