1! Copyright (C) 2007 Quantum ESPRESSO group
2! This file is distributed under the terms of the
3! GNU General Public License. See the file `License'
4! in the root directory of the present distribution,
5! or http://www.gnu.org/copyleft/gpl.txt .
6!
7!
8SUBROUTINE pseudo_q (qfunc, qfuncl)
9USE kinds, ONLY : DP
10USE io_global, ONLY : stdout
11USE ld1_parameters, ONLY : nwfsx
12USE ld1inc, ONLY : rcut, lls,  grid, ndmx, lmx2, nbeta, ikk, ecutrho, &
13            rmatch_augfun, rmatch_augfun_nc
14IMPLICIT NONE
15!
16REAL(DP), INTENT(IN) :: qfunc(ndmx,nwfsx,nwfsx)
17REAL(DP), INTENT(OUT) :: qfuncl(ndmx,nwfsx,nwfsx,0:lmx2)
18REAL(DP),  EXTERNAL :: int_0_inf_dr
19!
20! variables for aug. functions generation
21!
22INTEGER  :: irc, ns, ns1, l1, l2, l3, lll, mesh, n, ik
23INTEGER  :: l1_e, l2_e
24REAL(DP) :: aux(ndmx)
25REAL(DP) :: augmom, ecutrhoq, rmatch
26
27ecutrho=0.0_DP
28mesh = grid%mesh
29qfuncl=0.0_DP
30do ns=1,nbeta
31   l1 = lls(ns)
32   do ns1=ns,nbeta
33      l2 = lls(ns1)
34      !
35      !  Find the matching point
36      !
37      ik=0
38      IF (rmatch_augfun_nc) THEN
39         rmatch=min(rcut(ns),rcut(ns1))
40      ELSE
41         rmatch=rmatch_augfun
42      ENDIF
43      do n=1,mesh
44         if (grid%r(n)>rmatch) then
45            ik=n
46            exit
47         endif
48      enddo
49      if (ik==0.or.ik>mesh-20) call errore('pseudo_q','wrong rmatch_augfun',1)
50      !
51      ! Do the pseudization
52      !
53      do l3 = abs(l1-l2), l1+l2, 2
54         CALL compute_q_3bess(l3,l1+l2,ik,qfunc(1,ns,ns1), &
55                                          qfuncl(1,ns,ns1,l3),ecutrhoq)
56         IF (ecutrhoq>ecutrho) then
57            ecutrho=ecutrhoq
58            l1_e=l1
59            l2_e=l2
60         ENDIF
61         qfuncl(1:mesh,ns1,ns,l3)=qfuncl(1:mesh,ns,ns1,l3)
62      end do
63   end do
64end do
65!
66!  Check that multipoles have not changed
67!
68irc = maxval(ikk(1:nbeta))+8
69augmom=0.0_DP
70DO ns=1,nbeta
71   l1=lls(ns)
72   DO ns1=ns,nbeta
73      l2=lls(ns1)
74      DO l3 = abs(l1-l2), l1+l2, 2
75         aux(1:irc) = (qfuncl(1:irc,ns,ns1,l3)-qfunc(1:irc,ns,ns1)) &
76                                               * grid%r(1:irc)**l3
77         lll = l1 + l2 + 2 + l3
78         augmom=int_0_inf_dr(aux(1:irc),grid,irc,lll)
79
80         IF (abs(augmom)>1.d-5) WRITE (stdout,'(5x,a,2i3,a,2i3,a,i3,f15.7)') &
81              " Problem with multipole",ns,l1,":",ns1,l2, " l3=",l3, augmom
82      END DO
83   END DO
84END DO
85WRITE(stdout,'(/,5x, "Q pseudized with Bessel functions")')
86WRITE(stdout,'(5x,"Expected ecutrho= ",f12.4," due to l1=",i3,"   l2=",i3)') &
87                            ecutrho, l1_e, l2_e
88RETURN
89END SUBROUTINE pseudo_q
90