1!
2! Copyright (C) 2001 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!
8!
9!----------------------------------------------------------------------
10subroutine addusddenseq (drhoscf, dbecsum)
11  !----------------------------------------------------------------------
12  !
13  !  This routine adds to the change of the charge and magnetization
14  !  densities due to an electric field perturbation
15  !  the part due to the US augmentation.
16  !  It assumes that the array dbecsum has already accumulated the
17  !  change of the becsum term.
18  !  The expression implemented is given in Eq. B32 of PRB 64, 235118
19  !  (2001) with b=c=0.
20  !
21
22  USE kinds, only : DP
23  USE ions_base, ONLY : nat, ityp, ntyp => nsp
24  USE cell_base, ONLY : tpiba
25  use fft_base,  only: dfftp
26  use fft_interfaces, only: invfft
27  USE gvect, ONLY : g, gg, ngm, eigts1, eigts2, eigts3, mill
28  USE uspp, ONLY: okvan
29  USE uspp_param, ONLY: upf, lmaxq, nh, nhm
30  USE noncollin_module, ONLY : nspin_mag
31
32  USE qpoint, ONLY : xq, eigqts
33  implicit none
34  !
35  !   the dummy variables
36  !
37
38  ! input: if zero does not compute drho
39  ! input: the number of perturbations
40
41  complex(DP) :: drhoscf(dfftp%nnr,nspin_mag,1), &
42                 dbecsum(nhm*(nhm+1)/2,nat,nspin_mag,1)
43
44  ! inp/out: change of the charge density
45  ! input: sum over kv of bec
46  !
47  !     here the local variables
48  !
49
50  integer :: ig, na, nt, ih, jh, ijh, is
51
52  ! counters
53
54  real(DP), allocatable  :: qmod(:), qpg(:,:), ylmk0(:,:)
55  ! the modulus of q+G
56  ! the spherical harmonics
57
58  complex(DP) :: zsum
59  complex(DP), allocatable ::  sk (:), qg (:), qgm (:), aux (:,:)
60  ! the structure factor
61  ! work space
62
63  if (.not.okvan) return
64  call start_clock ('addusddenseq')
65  allocate (aux(  ngm, nspin_mag))
66  allocate (sk (  ngm))
67  allocate (qg (  dfftp%nnr))
68  allocate (ylmk0(ngm , lmaxq * lmaxq))
69  allocate (qgm  (ngm))
70  allocate (qmod (ngm))
71  allocate (qpg(3, ngm))
72  !
73  !  And then we compute the additional charge in reciprocal space
74  !
75  call setqmod (ngm, xq, g, qmod, qpg)
76  call ylmr2 (lmaxq * lmaxq, ngm, qpg, qmod, ylmk0)
77  do ig = 1, ngm
78     qmod (ig) = sqrt (qmod (ig) ) * tpiba
79  enddo
80
81  aux (:,:) = (0.d0, 0.d0)
82  do nt = 1, ntyp
83     if (upf(nt)%tvanp ) then
84        ijh = 0
85        do ih = 1, nh (nt)
86           do jh = ih, nh (nt)
87              call qvan2 (ngm, ih, jh, nt, qmod, qgm, ylmk0)
88              ijh = ijh + 1
89              do na = 1, nat
90                 if (ityp (na) == nt) then
91                    !
92                    ! calculate the structure factor
93                    !
94                    do ig = 1, ngm
95                       sk(ig)=eigts1(mill(1,ig),na)*eigts2(mill(2,ig),na) &
96                             *eigts3(mill(3,ig),na)*eigqts(na)*qgm(ig)
97                    enddo
98                    !
99                    !  And qgmq and becp and dbecq
100                    !
101                    do is=1,nspin_mag
102                       zsum = dbecsum (ijh, na, is, 1)
103                       call zaxpy(ngm,zsum,sk,1,aux(1,is),1)
104                    enddo
105                 endif
106              enddo
107           enddo
108        enddo
109     endif
110  enddo
111  !
112  !     convert aux to real space
113  !
114  do is=1,nspin_mag
115     qg (:) = (0.d0, 0.d0)
116     qg (dfftp%nl (:) ) = aux (:, is)
117     CALL invfft ('Rho', qg, dfftp)
118     drhoscf(:,is,1) = drhoscf(:,is,1) + 2.d0*qg(:)
119  enddo
120
121  deallocate (qpg)
122  deallocate (qmod)
123  deallocate (qgm)
124  deallocate (ylmk0)
125  deallocate (qg)
126  deallocate (sk)
127  deallocate (aux)
128
129  call stop_clock ('addusddenseq')
130  return
131end subroutine addusddenseq
132