TABLE OF CONTENTS
ABINIT/psp5nl [ Functions ]
NAME
psp5nl
FUNCTION
Make Kleinman-Bylander form factors f_l(q) for each l from 0 to lmax; Vloc is assumed local potential.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, FrD, GZ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
al=grid spacing in exponent for radial grid
lmax=maximum ang momentum for which nonlocal form factor is desired.
Usually lmax=1, sometimes = 0 (e.g. for oxygen); lmax <= 2 allowed.
mmax=number of radial grid points for atomic grid
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mqgrid=number of grid points for q grid
qgrid(mqgrid)=values at which form factors are returned
rad(mmax)=radial grid values
vloc(mmax)=local pseudopotential on radial grid
vpspll(mmax,3)=nonlocal pseudopotentials for each l on radial grid
wfll(mmax,3)=reference state wavefunctions on radial grid
mmax and mqgrid
OUTPUT
ekb(mpsang)=Kleinman-Bylander energy, {{\\ \begin{equation} \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]} {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]} \end{equation} }} for each l ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum
NOTES
u_l(r) is reference state wavefunction (input as wf);
j_l(q) is a spherical Bessel function;
dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
f_l(q) = $ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
where dvms = $\int_0^{rmax} [(u_l(r) dV_l(r))^2 dr]$ is the mean
square value of the nonlocal correction for angular momentum l.
Xavier Gonze s E_KB = $ dvms/\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]$.
This is the eigenvalue of the Kleinman-Bylander operator and sets
the energy scale of the nonlocal psp corrections.
PARENTS
psp5in,psp6in
CHILDREN
ctrap,spline
SOURCE
60 #if defined HAVE_CONFIG_H 61 #include "config.h" 62 #endif 63 64 #include "abi_common.h" 65 66 subroutine psp5nl(al,ekb,ffspl,lmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll) 67 68 use defs_basis 69 use m_profiling_abi 70 use m_splines 71 use m_errors 72 73 use m_numeric_tools, only : ctrap 74 75 !This section has been created automatically by the script Abilint (TD). 76 !Do not modify the following lines by hand. 77 #undef ABI_FUNC 78 #define ABI_FUNC 'psp5nl' 79 !End of the abilint section 80 81 implicit none 82 83 !Arguments ------------------------------------ 84 !scalars 85 integer,intent(in) :: lmax,mmax,mpsang,mqgrid 86 real(dp),intent(in) :: al 87 !arrays 88 real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax),vpspll(mmax,mpsang) 89 real(dp),intent(in) :: wfll(mmax,mpsang) 90 real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang) 91 92 !Local variables------------------------------- 93 !DEBUG 94 !real(dp) :: norm,wf 95 !ENDDEBUG 96 !scalars 97 integer,parameter :: dpsang=5 98 integer :: iq,ir,lp1 99 real(dp) :: arg,bessel,dvwf,qr,result,yp1,ypn,ztor1 100 character(len=500) :: message 101 !arrays 102 real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang) 103 real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:) 104 105 !************************************************************************* 106 107 !l=0,1,2 and 3 spherical Bessel functions 108 !The accuracy of the bes1, bes2, bes3 functions for small arguments 109 !may be insufficient. In the present version 110 !of the routines, some care is taken with the value of the argument. 111 !If smaller than 1.d-3, a two terms 112 !Taylor series expansion is prefered. 113 ! bes0(arg)=sin(arg)/arg 114 ! bes1(arg)=(sin(arg)-arg*cos(arg))/arg**2 115 ! bes2(arg)=( (3.0d0-arg**2)*sin(arg)-& 116 !& 3.0d0*arg*cos(arg) ) /arg**3 117 118 ! bes3(arg)=(15.d0*sin(arg)-15.d0*arg*cos(arg) & 119 !& -6.d0*arg**2*sin(arg)+arg**3*cos(arg) )/arg**4 120 121 !Zero out Kleinman-Bylander energies ekb 122 ekb(:)=0.0d0 123 124 ABI_ALLOCATE(work1,(mmax)) 125 ABI_ALLOCATE(work2,(mmax)) 126 ABI_ALLOCATE(work3,(mmax)) 127 ABI_ALLOCATE(work4,(mmax)) 128 129 !Allow for no nonlocal correction (lmax=-1) 130 if (lmax/=-1) then 131 132 ! Check that lmax is within allowed range 133 if (lmax<0.or.lmax>3) then 134 write(message, '(a,i12,a,a,a,a,a,a,a)' )& 135 & 'lmax=',lmax,' is not an allowed value.',ch10,& 136 & 'Allowed values are -1 for no nonlocal correction or else',ch10,& 137 & '0, 1,2 or 3 for maximum l nonlocal correction.',ch10,& 138 & 'Action : check the input atomic psp data file for lmax.' 139 MSG_ERROR(message) 140 end if 141 142 ! Compute normalizing integrals eta=<dV> and mean square 143 ! nonlocal psp correction dvms=<dV^2> 144 ! "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction 145 do lp1=1,lmax+1 146 147 ! integral from 0 to r1 148 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1) 149 ztor1=(wfll(1,lp1)*dvwf)*rad(1)/dble(2*(lp1-1)+3) 150 ! integrand for r1 to r(mmax) (incl extra factor of r) 151 do ir=1,mmax 152 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 153 work1(ir)=rad(ir)*(wfll(ir,lp1)*dvwf) 154 end do 155 ! do integral by corrected trapezoidal integration 156 call ctrap(mmax,work1,al,result) 157 eta(lp1)=ztor1+result 158 159 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1) 160 ztor1=dvwf**2*rad(1)/dble(2*(lp1-1)+3) 161 do ir=1,mmax 162 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 163 work1(ir)=rad(ir)*(dvwf**2) 164 end do 165 call ctrap(mmax,work1,al,result) 166 dvms(lp1)=ztor1+result 167 168 ! DEBUG 169 ! Compute the norm of wfll 170 ! wf=wfll(1,lp1) 171 ! ztor1=wf**2*rad(1)/dble(2*(lp1-1)+3) 172 ! do ir=1,mmax 173 ! wf=wfll(ir,lp1) 174 ! work1(ir)=rad(ir)*(wf**2) 175 ! end do 176 ! call ctrap(mmax,work1,al,result) 177 ! norm=ztor1+result 178 ! write(std_out,*)' lp1, norm',lp1,norm 179 ! ENDDEBUG 180 181 ! If dvms is not 0 for any given angular momentum l, 182 ! compute Xavier Gonze s definition of the Kleinman-Bylander 183 ! energy E_KB = dvms/eta. In this case also renormalize 184 ! the projection operator to u_KB(r)=$u_l(r)*dV(r)/\sqrt{dvms}$. 185 ! This means dvwf gets multiplied by the normalization factor 186 ! "renorm"=$1/\sqrt{dvms}$ as seen below. 187 if (dvms(lp1)/=0.0d0) then 188 ekb(lp1)=dvms(lp1)/eta(lp1) 189 renorm(lp1)=1.0d0/sqrt(dvms(lp1)) 190 ! ckb is Kleinman-Bylander "cosine" (Xavier Gonze) 191 ckb(lp1)=eta(lp1)/sqrt(dvms(lp1)) 192 else 193 ekb(lp1)=0.0d0 194 end if 195 196 end do 197 198 ! l=0 form factor if ekb(1) not 0 (lmax always at least 0) 199 if (ekb(1)/=0.0d0) then 200 201 ! do q=0 separately 202 lp1=1 203 ! 0 to r1 integral 204 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 205 ztor1=(rad(1)*dvwf)*rad(1)/3.0d0 206 ! integrand 207 do ir=1,mmax 208 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 209 work1(ir)=rad(ir)*(rad(ir)*dvwf) 210 end do 211 call ctrap(mmax,work1,al,result) 212 ffspl(1,1,1)=ztor1+result 213 214 ! do rest of q points 215 do iq=2,mqgrid 216 arg=two_pi*qgrid(iq) 217 ! 0 to r1 integral 218 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 219 ztor1=(bes0_psp5(arg*rad(1))*rad(1)*dvwf)*rad(1)/3.0d0 220 ! integrand 221 do ir=1,mmax 222 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 223 work1(ir)=rad(ir)*(rad(ir)*bes0_psp5(arg*rad(ir))*dvwf) 224 end do 225 call ctrap(mmax,work1,al,result) 226 ffspl(iq,1,1)=ztor1+result 227 end do 228 229 ! Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid) 230 ! yp1=0 for l=0 231 yp1=0.0d0 232 ! ypn=$ \int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$ 233 arg=two_pi*qgrid(mqgrid) 234 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 235 qr=arg*rad(1) 236 if(qr<1.d-3)then 237 bessel=(10.d0-qr*qr)*qr/30.0d0 238 else 239 bessel=bes1_psp5(qr) 240 end if 241 ! ztor1=(-bes1(arg*rad(1))*two_pi*rad(1)*r(1)*dvwf)*rad(1)/5.0d0 242 ztor1=(-bessel*two_pi*rad(1)*rad(1)*dvwf)*rad(1)/5.0d0 243 do ir=1,mmax 244 qr=arg*rad(ir) 245 if(qr<1.d-3)then 246 bessel=(10.d0-qr*qr)*qr/30.0d0 247 else 248 bessel=bes1_psp5(qr) 249 end if 250 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 251 ! work(ir)=rad(ir)*(-bes1(arg*rad(ir))*two_pi*rad(ir)*rad(ir)*dvwf) 252 work1(ir)=rad(ir)*(-bessel*two_pi*rad(ir)*rad(ir)*dvwf) 253 end do 254 call ctrap(mmax,work1,al,result) 255 ypn=ztor1+result 256 257 ! Fit spline to get second derivatives by spline fit 258 call spline(qgrid,ffspl(1,1,1),mqgrid,yp1,ypn,ffspl(1,2,1)) 259 260 else 261 ! or else put nonlocal correction at l=0 to 0 262 ffspl(:,:,1)=0.0d0 263 end if 264 265 ! Finished if lmax=0 (highest nonlocal correction) 266 ! Do l=1 form factor if ekb(2) not 0 and lmax>=1 267 if (lmax>0)then 268 if(ekb(2)/=0.0d0) then 269 270 lp1=2 271 ! do q=0 separately: f_1(q=0) vanishes ! 272 ffspl(1,1,2)=0.0d0 273 274 ! do rest of q points 275 do iq=2,mqgrid 276 arg=two_pi*qgrid(iq) 277 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 278 qr=arg*rad(1) 279 if(qr<1.d-3)then 280 bessel=(10.d0-qr*qr)*qr/30.0d0 281 else 282 bessel=bes1_psp5(qr) 283 end if 284 ! ztor1=(bes1(arg*rad(1))*rad(1)*dvwf)*rad(1)/5.0d0 285 ztor1=(bessel*rad(1)*dvwf)*rad(1)/5.0d0 286 287 do ir=1,mmax 288 qr=arg*rad(ir) 289 if(qr<1.d-3)then 290 bessel=(10.d0-qr*qr)*qr/30.0d0 291 else 292 bessel=bes1_psp5(qr) 293 end if 294 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 295 work2(ir)=rad(ir)*(rad(ir)*bessel*dvwf) 296 end do 297 298 call ctrap(mmax,work2,al,result) 299 ffspl(iq,1,2)=ztor1+result 300 end do 301 302 ! Compute yp1,ypn for l=1 303 ! yp1=$\displaystyle \int [2\pi r^2 wf(r) dV(r)]/3$ 304 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 305 ztor1=((two_pi*rad(1)**2)*dvwf)*rad(1)/(3.0d0*5.0d0) 306 do ir=1,mmax 307 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 308 work2(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf/3.0d0) 309 end do 310 call ctrap(mmax,work2,al,result) 311 yp1=ztor1+result 312 ! ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$ 313 ! where x=2 Pi qgrid(mqgrid) r 314 arg=two_pi*qgrid(mqgrid) 315 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 316 qr=arg*rad(1) 317 if(qr<1.d-3)then 318 bessel=(10.d0-3.0d0*qr*qr)/30.0d0 319 else 320 bessel=bes0_psp5(qr)-2.d0*bes1_psp5(qr)/qr 321 end if 322 ! ztor1=( (two_pi*rad(1)**2)*dvwf* (bes0(arg*rad(1))- 323 ! 2.0d0*bes1(arg*rad(1))/(arg*rad(1))) ) * rad(1)/5.0d0 324 ztor1=( (two_pi*rad(1)**2)*dvwf*bessel)* rad(1)/5.0d0 325 326 do ir=1,mmax 327 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 328 qr=arg*rad(ir) 329 if(qr<1.d-3)then 330 bessel=(10.d0-3.0d0*qr*qr)/30.0d0 331 else 332 bessel=bes0_psp5(qr)-2.d0*bes1_psp5(qr)/qr 333 end if 334 ! work(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf* 335 ! (bes0(arg*rad(ir))-2.d0*bes1(arg*rad(ir))/(arg*rad(ir))) ) 336 work2(ir)=rad(ir)*(two_pi*rad(ir)**2)*dvwf*bessel 337 end do 338 call ctrap(mmax,work2,al,result) 339 ypn=ztor1+result 340 341 ! Fit spline for l=1 Kleinman-Bylander form factor 342 call spline(qgrid,ffspl(1,1,2),mqgrid,yp1,ypn,ffspl(1,2,2)) 343 344 else 345 ! or else put form factor to 0 for l=1 346 ffspl(:,:,2)=0.0d0 347 end if 348 ! Endif condition of lmax>0 349 end if 350 351 ! Finished if lmax=1 (highest nonlocal correction) 352 ! Do l=2 nonlocal form factor if eta(3) not 0 and lmax>=2 353 if (lmax>1)then 354 if(ekb(3)/=0.0d0) then 355 356 lp1=3 357 ! do q=0 separately; f_2(q=0) vanishes 358 ffspl(1,1,3)=0.0d0 359 360 ! do rest of q points 361 do iq=2,mqgrid 362 arg=two_pi*qgrid(iq) 363 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 364 qr=arg*rad(1) 365 if(qr<1.d-3)then 366 bessel=qr*qr/15.0d0-qr**4/210.0d0 367 else 368 bessel=bes2_psp5(qr) 369 end if 370 ! ztor1=(bes2(arg*rad(1))*rad(1)*dvwf)*rad(1)/7.0d0 371 ztor1=(bessel*rad(1)*dvwf)*rad(1)/7.0d0 372 do ir=1,mmax 373 qr=arg*rad(ir) 374 if(qr<1.d-3)then 375 bessel=qr*qr/15.0d0-qr**4/210.0d0 376 else 377 bessel=bes2_psp5(qr) 378 end if 379 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 380 ! work(ir)=rad(ir)*(r(ir)*bes2(arg*rad(ir))*dvwf) 381 work3(ir)=rad(ir)*(rad(ir)*bessel*dvwf) 382 end do 383 call ctrap(mmax,work3,al,result) 384 ffspl(iq,1,3)=ztor1+result 385 end do 386 387 ! Compute yp1,ypn for l=2 388 ! yp1=0 for l=2 389 yp1=0.0d0 390 ! ypn=$\int [2 \pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$ 391 ! where x=2 Pi qgrid(mqgrid) r 392 arg=two_pi*qgrid(mqgrid) 393 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 394 qr=arg*rad(1) 395 if(qr<1.d-3)then 396 bessel=qr*2.0d0/15.0d0-qr**3*4.0d0/210.0d0 397 else 398 bessel=bes1_psp5(qr)-3.0d0*bes2_psp5(qr)/qr 399 end if 400 ! ztor1=( (two_pi*rad(1)**2)*dvwf* (bes1(arg*rad(1))- 401 ! 3.0d0*bes2(arg*rad(1))/(arg*rad(1))) ) * rad(1)/7.0d0 402 ztor1=( (two_pi*rad(1)**2)*dvwf* bessel ) * rad(1)/7.0d0 403 do ir=1,mmax 404 qr=arg*rad(ir) 405 if(qr<1.d-3)then 406 bessel=qr*2.0d0/15.0d0-qr**3*4.0d0/210.0d0 407 else 408 bessel=bes1_psp5(qr)-3.0d0*bes2_psp5(qr)/qr 409 end if 410 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 411 ! work3(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf* 412 ! (bes1(arg*rad(ir))-3.d0*bes2(arg*rad(ir))/(arg*rad(ir))) ) 413 work3(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*bessel) 414 end do 415 call ctrap(mmax,work3,al,result) 416 ypn=ztor1+result 417 418 ! Fit spline for l=2 Kleinman-Bylander form factor 419 call spline(qgrid,ffspl(1,1,3),mqgrid,yp1,ypn,ffspl(1,2,3)) 420 421 else 422 ! or else put form factor to 0 for l=1 423 ffspl(:,:,3)=0.0d0 424 end if 425 ! Endif condition of lmax>1 426 end if 427 428 ! Finished if lmax=2 (highest nonlocal correction) 429 ! Do l=3 nonlocal form factor if eta(4) not 0 and lmax>=3 430 if (lmax>2)then 431 if(ekb(4)/=0.0d0) then 432 433 lp1=4 434 ! do q=0 separately; f_3(q=0) vanishes 435 ffspl(1,1,4)=0.0d0 436 437 ! do rest of q points 438 do iq=2,mqgrid 439 arg=two_pi*qgrid(iq) 440 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 441 qr=arg*rad(1) 442 if(qr<1.d-3)then 443 bessel=qr*qr*qr/105.0d0-qr**5/1890.0d0+qr**7/83160.0d0 444 else 445 bessel=bes3_psp5(qr) 446 end if 447 ! ztor1=(bes3(arg*rad(1))*rad(1)*dvwf)*rad(1)/9.0d0 448 ztor1=(bessel*rad(1)*dvwf)*rad(1)/9.0d0 449 do ir=1,mmax 450 qr=arg*rad(ir) 451 if(qr<1.d-3)then 452 bessel=qr*qr*qr/105.0d0-qr**5/1890.0d0+qr**7/83160.0d0 453 else 454 bessel=bes3_psp5(qr) 455 end if 456 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 457 ! work(ir)=rad(ir)*(rad(ir)*bes3(arg*rad(ir))*dvwf) 458 work4(ir)=rad(ir)*(rad(ir)*bessel*dvwf) 459 end do 460 call ctrap(mmax,work4,al,result) 461 ffspl(iq,1,4)=ztor1+result 462 end do 463 464 ! Compute yp1,ypn for l=3 465 ! yp1=0 for l=3 466 yp1=0.0d0 467 ! ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$ 468 ! where x=2 Pi qgrid(mqgrid) r 469 arg=two_pi*qgrid(mqgrid) 470 dvwf=(vpspll(1,lp1)-vloc(1))*wfll(1,lp1)*renorm(lp1) 471 qr=arg*rad(1) 472 if(qr<1.d-3)then 473 bessel=3.d0*qr**2/105.0d0-5.d0*qr**4/1890.0d0+7.d0*qr**6/83160.0d0 474 else 475 bessel=bes2_psp5(qr)-4.0d0*bes3_psp5(qr)/qr 476 end if 477 ! ztor1=( (two_pi*rad(1)**2)*dvwf* (bes2(arg*rad(1))- 478 ! 3.0d0*bes3(arg*rad(1))/(arg*rad(1))) ) * rad(1)/9.0d0 479 ztor1=( (two_pi*rad(1)**2)*dvwf* bessel ) * rad(1)/9.0d0 480 do ir=1,mmax 481 qr=arg*rad(ir) 482 if(qr<1.d-3)then 483 bessel=3.d0*qr**2/105.0d0-5.d0*qr**4/1890.0d0+7.d0*qr**6/83160.0d0 484 else 485 bessel=bes2_psp5(qr)-4.0d0*bes3_psp5(qr)/qr 486 end if 487 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 488 ! work4(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf* 489 ! (bes2(arg*rad(ir))-4.d0*bes3(arg*rad(ir))/(arg*rad(ir))) ) 490 work4(ir)=rad(ir)*((two_pi*rad(ir)**2)*dvwf*bessel) 491 end do 492 call ctrap(mmax,work4,al,result) 493 ypn=ztor1+result 494 495 ! Fit spline for l=3 Kleinman-Bylander form factor 496 call spline(qgrid,ffspl(1,1,4),mqgrid,yp1,ypn,ffspl(1,2,4)) 497 498 else 499 ! or else put form factor to 0 for l=3 500 ffspl(:,:,4)=0.0d0 501 end if 502 ! Endif condition of lmax>2 503 end if 504 505 ! Endif condition lmax/=-1 506 end if 507 508 !DEBUG 509 !write(std_out,*) 'EKB=',(ekb(iq),iq=1,3) 510 !write(std_out,*) 'COSKB=',(ckb(iq),iq=1,3) 511 !ENDDEBUG 512 513 ABI_DEALLOCATE(work1) 514 ABI_DEALLOCATE(work2) 515 ABI_DEALLOCATE(work3) 516 ABI_DEALLOCATE(work4) 517 518 contains 519 520 function bes0_psp5(arg) 521 522 523 !This section has been created automatically by the script Abilint (TD). 524 !Do not modify the following lines by hand. 525 #undef ABI_FUNC 526 #define ABI_FUNC 'bes0_psp5' 527 !End of the abilint section 528 529 real(dp) :: bes0_psp5,arg 530 bes0_psp5=sin(arg)/arg 531 end function bes0_psp5 532 533 function bes1_psp5(arg) 534 535 536 !This section has been created automatically by the script Abilint (TD). 537 !Do not modify the following lines by hand. 538 #undef ABI_FUNC 539 #define ABI_FUNC 'bes1_psp5' 540 !End of the abilint section 541 542 real(dp) :: bes1_psp5,arg 543 bes1_psp5=(sin(arg)-arg*cos(arg))/arg**2 544 end function bes1_psp5 545 546 function bes2_psp5(arg) 547 548 549 !This section has been created automatically by the script Abilint (TD). 550 !Do not modify the following lines by hand. 551 #undef ABI_FUNC 552 #define ABI_FUNC 'bes2_psp5' 553 !End of the abilint section 554 555 real(dp) :: bes2_psp5,arg 556 bes2_psp5=( (3.0d0-arg**2)*sin(arg)-& 557 & 3.0d0*arg*cos(arg) ) /arg**3 558 end function bes2_psp5 559 560 function bes3_psp5(arg) 561 562 563 !This section has been created automatically by the script Abilint (TD). 564 !Do not modify the following lines by hand. 565 #undef ABI_FUNC 566 #define ABI_FUNC 'bes3_psp5' 567 !End of the abilint section 568 569 real(dp) :: bes3_psp5, arg 570 bes3_psp5=(15.d0*sin(arg)-15.d0*arg*cos(arg) & 571 & -6.d0*arg**2*sin(arg)+arg**3*cos(arg) )/arg**4 572 end function bes3_psp5 573 574 end subroutine psp5nl