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! 9!-------------------------------------------------------------------------- 10subroutine compute_chi(lam,ikk_in,phi_in,chi_out,xc,e,lbes4) 11 !-------------------------------------------------------------------------- 12 ! 13 ! This routine computes the chi functions: 14 ! |chi> = (\epsilon -T -V_{loc}) |psi> 15 ! 16 use io_global, only : stdout, ionode 17 use kinds, only : DP 18 use radial_grids, only: ndmx, series 19 use ld1inc, only: grid, vpsloc, rho0, verbosity 20 21 implicit none 22 integer :: & 23 ikk_in,& ! the point after which the chi should be zero 24 lam ! the angular momentum 25 logical :: & 26 lbes4 27 28 real(DP) :: & 29 e, & ! input: the energy 30 xc(8), & ! input: the coefficients of the bessel function 31 phi_in(ndmx), & ! input: pseudo wavefunction 32 chi_out(ndmx) ! output: the chi function 33 ! 34 real(DP) :: & 35 j1(ndmx), aux(ndmx), gi(ndmx),& 36 b(4),c(4), arow(ndmx),brow(ndmx),crow(ndmx),drow(ndmx), & 37 b0e, g0, g1, g2, & 38 r_clean, & 39 ddx12, & 40 x4l6, & 41 x6l12, & 42 int_0_inf_dr, & 43 integral 44 45 integer :: & 46 n, nstart, nst 47 ! 48 ! RRKJ: first expand in a taylor series the phis function 49 ! Since we know that the phis functions are a sum of Bessel 50 ! functions with coefficients xc, we compute analytically 51 ! the asymptotic expansion 52 ! 53 ! 54 ddx12=grid%dx*grid%dx/12.0_dp 55 x4l6=4*lam+6 56 x6l12=6*lam+12 57 58 do n=1,6 59 j1(n)=phi_in(n)/grid%r(n)**(lam+1) 60 enddo 61 call seriesbes(j1,grid%r,grid%r2,6,c) 62 ! 63 if (lam == 0) then 64 if(lbes4.or.rho0.eq.0.0_dp)then 65 c(1)=xc(1)+xc(2)+xc(3) 66 c(2)=0.0_dp 67 c(3)=-xc(1)*(xc(4)**2/6.0_dp) & 68 -xc(2)*(xc(5)**2/6.0_dp) & 69 -xc(3)*(xc(6)**2/6.0_dp) 70 c(4)=0.0_dp 71 else 72 c(1)=xc(1)+xc(2)+xc(3)+xc(4) 73 c(2)=0.0_dp 74 c(3)=-xc(1)*(xc(5)**2/6.0_dp) & 75 -xc(2)*(xc(6)**2/6.0_dp) & 76 -xc(3)*(xc(7)**2/6.0_dp) & 77 -xc(4)*(xc(8)**2/6.0_dp) 78 c(4)=0.0_dp 79 endif 80 elseif (lam == 3) then 81 c(1)=xc(1)*(48.0_dp*xc(4)**3/5040.0_dp)+ & 82 xc(2)*(48.0_dp*xc(5)**3/5040.0_dp)+ & 83 xc(3)*(48.0_dp*xc(6)**3/5040.0_dp) 84 c(2)=0.0_dp 85 c(3)=-xc(1)*(192.0_dp*xc(4)**5/362880.0_dp) & 86 -xc(2)*(192.0_dp*xc(5)**5/362880.0_dp) & 87 -xc(3)*(192.0_dp*xc(6)**5/362880.0_dp) 88 c(4)=0.0_dp 89 elseif (lam == 2) then 90 c(1)=xc(1)*(xc(4)**2/15.0_dp)+ & 91 xc(2)*(xc(5)**2/15.0_dp)+ & 92 xc(3)*(xc(6)**2/15.0_dp) 93 c(2)=0.0_dp 94 c(3)=-xc(1)*(xc(4)**4/210.0_dp) & 95 -xc(2)*(xc(5)**4/210.0_dp) & 96 -xc(3)*(xc(6)**4/210.0_dp) 97 c(4)=0.0_dp 98 elseif (lam == 1) then 99 c(1)=xc(1)*(xc(4)/3.0_dp)+ & 100 xc(2)*(xc(5)/3.0_dp)+ & 101 xc(3)*(xc(6)/3.0_dp) 102 c(2)=0.0_dp 103 c(3)=-xc(1)*(xc(4)**3/30.0_dp) & 104 -xc(2)*(xc(5)**3/30.0_dp) & 105 -xc(3)*(xc(6)**3/30.0_dp) 106 c(4)=0.0_dp 107 else 108 call errore('compute_chi','lam not programmed',1) 109 endif 110 ! 111 ! and the potential 112 ! 113 do n=1,4 114 j1(n)=vpsloc(n) 115 enddo 116 if (abs(j1(1)-j1(4))>1.d-12) then 117 call series(j1,grid%r,grid%r2,b) 118 else 119 b=0.0_dp 120 b(1)=j1(1) 121 endif 122 ! 123 ! and compute the taylor expansion of the chis 124 ! 125 b0e=(b(1)-e) 126 127 g0=x4l6*c(3)-b0e*c(1) 128 g1=x6l12*c(4)-c(1)*b(2) 129 g2=-(b0e*c(3)+b(3)*c(1)) 130 nstart=5 131 do n=1,nstart-1 132 chi_out(n)= (g0+grid%r(n)*(g1+g2*grid%r(n)))*grid%r(n)**(lam+3)/grid%sqr(n) 133 enddo 134 do n=1,grid%mesh 135 aux(n)= (g0+grid%r(n)*(g1+g2*grid%r(n))) 136 enddo 137 ! 138 ! set up the equation 139 ! 140 do n=1,grid%mesh 141 gi(n)=phi_in(n)/grid%sqr(n) 142 enddo 143 do n=1,grid%mesh 144 j1(n)=grid%r2(n)*(vpsloc(n)-e)+(lam+0.5_dp)**2 145 j1(n)=1.0_dp-ddx12*j1(n) 146 enddo 147 148 do n=nstart,grid%mesh-3 149 drow(n)= gi(n+1)*j1(n+1) & 150 + gi(n)*(-12.0_dp+10.0_dp*j1(n))+ & 151 gi(n-1)*j1(n-1) 152 153 brow(n)=10.0_dp*ddx12 154 crow(n)=ddx12 155 arow(n)=ddx12 156 enddo 157 drow(nstart)=drow(nstart)-ddx12*chi_out(nstart-1) 158 chi_out(grid%mesh-2)=0.0_dp 159 chi_out(grid%mesh-1)=0.0_dp 160 chi_out(grid%mesh)=0.0_dp 161 ! 162 ! and solve it 163 ! 164 call tridiag(arow(nstart),brow(nstart),crow(nstart), & 165 drow(nstart),chi_out(nstart),grid%mesh-3-nstart) 166 ! 167 ! put the correct normalization and r dependence 168 ! 169 do n=1,grid%mesh 170 chi_out(n)=chi_out(n)*grid%sqr(n)/grid%r2(n) 171 ! if(lam.eq.0) 172 ! + write(stdout,'(5(e20.13,1x))') 173 ! + r(n),chi_out(n),chi_out(n)/r(n)**(lam+1), 174 ! + aux(n),aux(n)*r(n)**(lam+1) 175 enddo 176 ! 177 ! smooth close to the origin with asymptotic expansion 178 ! 179 do n=nstart,grid%mesh 180 if (abs(chi_out(n)/grid%r(n)**(lam+1)-aux(n)) & 181 .lt.1.e-3_dp*abs(aux(n)) ) goto 100 182 chi_out(n)=aux(n)*grid%r(n)**(lam+1) 183 enddo 184 185100 if (n.eq.grid%mesh+1.or.grid%r(min(n,grid%mesh)).gt.0.05_dp)then 186 write(stdout,*) lam,n,grid%mesh,grid%r(min(n,grid%mesh)) 187 call errore('compute_chi','n is too large',1) 188 endif 189! 190! When the input wavefunction is a diverging scattering state 191! this routine might become numerically unstable at large r. 192! Here chi_out should be 0.0 but it is not due to this instability. 193! We now clean chi_out when phi_in > 20.0. 194! Clean also after 7.0 a.u.. 195! 196 r_clean=100.d0 197 do n=1,grid%mesh 198 if (abs(phi_in(n))>20.0_DP) then 199 r_clean=grid%r(n) 200 exit 201 endif 202 enddo 203 r_clean=min(r_clean,7.0_DP) 204 if (r_clean<grid%r(ikk_in)) & 205 call errore ('compute_chi ','phi_in too large before r_c', 1) 206 207 do n=grid%mesh,1,-1 208 if (grid%r(n).lt.r_clean) goto 200 209 chi_out(n)=0.0_dp 210 enddo 211200 continue 212 ! check that the chi are zero beyond ikk 213 nst=0 214 gi=0.0_dp 215 do n=ikk_in+1,grid%mesh 216 gi(n)=chi_out(n)**2 217 enddo 218 do n=min(ikk_in+20,grid%mesh),grid%mesh 219 chi_out(n)=0.0_dp 220 enddo 221 integral=int_0_inf_dr(gi,grid,grid%mesh,nst) 222 if (integral > 2.e-6_dp) then 223 write(stdout, '(5x,'' l='',i4, '' integral='',f15.9, & 224 & '' r(ikk) '',f15.9)') lam, integral, grid%r(ikk_in) 225 IF (verbosity=='high') THEN 226 do n=ikk_in,grid%mesh 227 write(stdout,*) grid%r(n),gi(n) 228 enddo 229 ENDIF 230 call errore ('compute_chi ','chi too large beyond r_c', 1) 231 endif 232 return 233end subroutine compute_chi 234 235 236subroutine tridiag(a,b,c,r,u,n) 237 ! 238 ! See Numerical Recipes. 239 ! 240 use kinds, only : DP 241 implicit none 242 243 integer :: n 244 real(DP) :: a(n),b(n),c(n),r(n),u(n) 245 real(DP) :: gam(n), bet 246 247 integer j 248 249 if (abs(b(1)).lt.1.e-10_DP) & 250 call errore('tridiag','b(1) is too small',1) 251 252 bet=b(1) 253 u(1)=r(1)/bet 254 do j=2,n 255 gam(j)=c(j-1)/bet 256 bet=b(j)-a(j)*gam(j) 257 if (abs(bet) < 1.e-10_DP) & 258 call errore('tridiag','bet is too small',1) 259 u(j)=(r(j)-a(j)*u(j-1))/bet 260 enddo 261 do j=n-1,1,-1 262 u(j)=u(j)-gam(j+1)*u(j+1) 263 enddo 264 return 265end subroutine tridiag 266