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