1!
2! Copyright (C) 2001-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!
8SUBROUTINE compute_qdipol( dpqq )
9  !! This routine computes the term dpqq, i.e. the dipole moment of the
10  !! augmentation charge. The output is given on cartesian coordinates.
11  !
12  USE kinds,       ONLY: DP
13  USE constants,   ONLY: fpi
14  USE atom,        ONLY: rgrid
15  USE ions_base,   ONLY: ntyp => nsp
16  USE uspp,        ONLY: nhtol, nhtolm, indv, nlx, ap
17  USE uspp_param,  ONLY: upf, nbetam, nh, nhm
18  !
19  IMPLICIT NONE
20  !
21  REAL(DP) :: dpqq(nhm, nhm, 3, ntyp)
22  !! Dipole moment of augmentation charge
23  REAL(DP), ALLOCATABLE :: qrad2(:,:,:), qtot(:,:,:), aux(:)
24  REAL(DP) :: fact
25  INTEGER :: nt, l, ir, nb, mb, ijv, ilast, ipol, ih, ivl, jh, jvl, lp, ndm
26  !
27  CALL start_clock( 'cmpt_qdipol' )
28  !
29  ndm = MAXVAL( upf(1:ntyp)%kkbeta )
30  ALLOCATE( qrad2( nbetam , nbetam, ntyp) )
31  ALLOCATE( aux( ndm)                     )
32  ALLOCATE( qtot( ndm, nbetam, nbetam)    )
33
34  qrad2(:,:,:)=0.d0
35  dpqq=0.d0
36
37  DO nt = 1, ntyp
38     IF ( upf(nt)%tvanp ) THEN
39        l=1
40        !
41        !   Only l=1 terms enter in the dipole of Q
42        !
43        DO nb = 1, upf(nt)%nbeta
44           DO mb = nb, upf(nt)%nbeta
45              ijv = mb * (mb-1) /2 + nb
46              IF ( ( l >= ABS(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .AND. &
47                   ( l <=     upf(nt)%lll(nb) + upf(nt)%lll(mb)  ) .AND. &
48                   (MOD(l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) THEN
49                 qtot(1:upf(nt)%kkbeta,nb,mb) = upf(nt)%qfuncl(1:upf(nt)%kkbeta,ijv,l)
50              ENDIF
51           ENDDO
52        ENDDO
53        do nb=1, upf(nt)%nbeta
54           !
55           !    the Q are symmetric with respect to indices
56           !
57           do mb=nb, upf(nt)%nbeta
58              if ( ( l >= ABS(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .AND. &
59                   ( l <=     upf(nt)%lll(nb) + upf(nt)%lll(mb)  ) .AND. &
60                   (MOD(l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) THEN
61                 do ir = 1, upf(nt)%kkbeta
62                    aux(ir)=rgrid(nt)%r(ir)*qtot(ir, nb, mb)
63                 ENDDO
64                 call simpson ( upf(nt)%kkbeta, aux, rgrid(nt)%rab, &
65                                qrad2(nb,mb,nt) )
66              ENDIF
67           ENDDO
68        ENDDO
69     ENDIF
70     ! ntyp
71  ENDDO
72
73  DO ipol = 1, 3
74     fact = -SQRT(fpi/3.d0)
75     IF (ipol == 1) lp=3
76     IF (ipol == 2) lp=4
77     IF (ipol == 3) THEN
78        lp=2
79        fact=-fact
80     ENDIF
81     DO nt = 1, ntyp
82        IF ( upf(nt)%tvanp ) THEN
83           DO ih = 1, nh(nt)
84              ivl = nhtolm(ih, nt)
85              mb = indv(ih, nt)
86              DO jh = ih, nh(nt)
87                 jvl = nhtolm(jh, nt)
88                 nb = indv(jh,nt)
89                 IF (ivl > nlx) CALL errore( 'compute_qdipol',' ivl > nlx', ivl )
90                 IF (jvl > nlx) CALL errore( 'compute_qdipol',' jvl > nlx', jvl )
91                 IF (nb > nbetam) &
92                    CALL errore( 'compute_qdipol',' nb out of bounds', nb )
93                 IF (mb > nbetam) &
94                    CALL errore( 'compute_qdipol',' mb out of bounds', mb )
95                 IF (mb > nb) CALL errore( 'compute_qdipol',' mb > nb', 1 )
96                 dpqq(ih,jh,ipol,nt) = fact * ap(lp,ivl,jvl) * qrad2(mb,nb,nt)
97                 dpqq(jh,ih,ipol,nt) = dpqq(ih,jh,ipol,nt)
98                 ! WRITE( stdout,'(3i5,2f15.9)') ih,jh,ipol,dpqq(ih,jh,ipol,nt)
99              ENDDO
100           ENDDO
101        ENDIF
102     ENDDO
103  ENDDO
104  !
105  DEALLOCATE( qtot  )
106  DEALLOCATE( aux   )
107  DEALLOCATE( qrad2 )
108  !
109  CALL stop_clock( 'cmpt_qdipol' )
110  !
111  RETURN
112  !
113END SUBROUTINE compute_qdipol
114