1! 2! Copyright (C) 2004 Vanderbilt's group at Rutgers University, NJ 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 calc_btq( ql, qr_k, idbes ) 10 !---------------------------------------------------------------------- 11 !! Calculates the Bessel-transform (or its derivative if idbes=1) 12 !! of the augmented qrad charges at a given ql point. 13 !! Rydberg atomic units are used. 14 ! 15 USE kinds, ONLY: DP 16 USE atom, ONLY: rgrid 17 USE ions_base, ONLY: ntyp => nsp 18 USE cell_base, ONLY: omega 19 USE constants, ONLY: fpi 20 USE uspp_param, ONLY: upf, nbetam, lmaxq 21 ! 22 IMPLICIT NONE 23 ! 24 REAL(DP) :: ql, qr_k(nbetam,nbetam,lmaxq,ntyp) 25 INTEGER :: idbes 26 ! 27 INTEGER :: i, np, l, ilmin, ilmax, iv, jv, ijv, ilast 28 REAL(DP) :: qrk 29 REAL(DP), ALLOCATABLE :: jl(:), aux(:) 30 ! 31 DO np = 1, ntyp 32 ! 33 IF ( upf(np)%tvanp ) THEN 34 ! 35 ALLOCATE( jl(upf(np)%kkbeta), aux(upf(np)%kkbeta) ) 36 DO iv = 1, upf(np)%nbeta 37 DO jv = iv, upf(np)%nbeta 38 ijv = jv * (jv-1) / 2 + iv 39 ilmin = ABS( upf(np)%lll(iv) - upf(np)%lll(jv) ) 40 ilmax = upf(np)%lll(iv) + upf(np)%lll(jv) 41 ! only need to calculate for l=lmin,lmin+2 ...lmax-2,lmax 42 DO l = ilmin,ilmax,2 43 aux(:) = 0.0_DP 44 aux(1:upf(np)%kkbeta) = upf(np)%qfuncl(1:upf(np)%kkbeta,ijv,l) 45 IF (idbes == 1) THEN 46 ! 47 CALL sph_dbes( upf(np)%kkbeta, rgrid(np)%r, ql, l, jl ) 48 ! 49 ELSE 50 ! 51 CALL sph_bes( upf(np)%kkbeta, rgrid(np)%r, ql, l, jl ) 52 ! 53 ENDIF 54 ! 55 ! jl is the Bessel function (or its derivative) calculated at ql 56 ! now integrate qfunc*jl*r^2 = Bessel transform of qfunc 57 ! 58 DO i = 1, upf(np)%kkbeta 59 aux(i) = jl(i)*aux(i) 60 ENDDO 61 ! if (tlog(np)) then 62 CALL simpson( upf(np)%kkbeta, aux, rgrid(np)%rab, qrk ) 63 ! 64 qr_k(iv,jv,l+1,np) = qrk * fpi / omega 65 qr_k(jv,iv,l+1,np) = qr_k(iv,jv,l+1,np) 66 ! 67 ENDDO 68 ENDDO 69 ENDDO 70 DEALLOCATE( aux, jl ) 71 ENDIF 72 ENDDO 73 ! 74 RETURN 75 ! 76END SUBROUTINE calc_btq 77