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