1!
2! Copyright (C) 2001-2008 Quantum ESPRESSO 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 drho
11  !-----------------------------------------------------------------------
12  !
13  !    Here we compute, for each mode the change of the charge density
14  !    due to the displacement, at fixed wavefunctions. These terms
15  !    are saved on disk. The orthogonality part is included in the
16  !    computed change.
17  !
18  !
19  !
20  USE kinds,      ONLY : DP
21  USE gvecs,         ONLY : doublegrid
22  USE fft_base,   ONLY : dfftp, dffts
23  USE lsda_mod,   ONLY : nspin
24  USE cell_base,  ONLY : omega
25  USE ions_base,  ONLY : nat
26  USE buffers,    ONLY : save_buffer
27  USE noncollin_module, ONLY : noncolin, npol, nspin_lsda, nspin_mag
28  USE uspp_param, ONLY : upf, nhm
29  USE uspp,       ONLY : okvan, nkb
30  USE wvfct,      ONLY : nbnd
31  USE paw_variables,    ONLY : okpaw
32  USE control_ph, ONLY : ldisp, all_done, rec_code_read
33
34  USE lrus,       ONLY : becp1
35  USE qpoint,     ONLY : nksq
36  USE control_lr, ONLY : lgamma
37
38  USE dynmat,     ONLY : dyn00
39  USE modes,      ONLY : npertx, npert, nirr
40  USE phus,       ONLY : becsumort, alphap
41  USE units_ph,   ONLY : lrdrhous, iudrhous
42
43  USE mp_pools,   ONLY : inter_pool_comm
44  USE mp_bands,   ONLY : intra_bgrp_comm
45  USE mp,         ONLY : mp_sum
46  USE becmod,     ONLY : bec_type, allocate_bec_type, deallocate_bec_type
47  USE fft_interfaces, ONLY : fft_interpolate
48
49  implicit none
50
51  integer :: mode, is, ir, irr, iper, npe, nrstot, nu_i, nu_j, ik, &
52             ipol
53  ! counter on modes
54  ! counter on atoms and polarizations
55  ! counter on atoms
56  ! counter on spin
57  ! counter on perturbations
58  ! the number of points
59  ! counter on modes
60  ! counter on k-point
61  ! counter on coordinates
62
63  real(DP), allocatable :: wgg (:,:,:)
64  ! the weight of each point
65
66
67  complex(DP) :: wdyn (3 * nat, 3 * nat)
68  type (bec_type), pointer :: becq(:), alpq(:,:)
69  complex(DP), allocatable :: dvlocin (:), drhous (:,:,:),&
70       drhoust (:,:,:), dbecsum(:,:,:,:), dbecsum_nc(:,:,:,:,:)
71  ! auxiliary to store bec at k+q
72  ! auxiliary to store alphap at
73  ! the change of the local potential
74  ! the change of the charge density
75  ! the change of the charge density
76  ! the derivative
77
78!
79!  The PAW case requires dbecsumort so we recalculate this starting part
80!  This will be changed soon
81!
82  if (all_done) return
83  if ((rec_code_read >=-20 .and..not.okpaw)) return
84
85  dyn00(:,:) = (0.d0,0.d0)
86  if (.not.okvan) return
87  call start_clock ('drho')
88  !
89  !    first compute the terms needed for the change of the charge density
90  !    due to the displacement of the augmentation charge
91  !
92  call compute_becsum_ph()
93  !
94  call compute_alphasum()
95  !
96  !    then compute the weights
97  !
98  allocate (wgg (nbnd ,nbnd , nksq))
99  if (lgamma) then
100     becq => becp1
101     alpq => alphap
102  else
103     allocate (becq ( nksq))
104     allocate (alpq ( 3, nksq))
105     do ik =1,nksq
106        call allocate_bec_type (  nkb, nbnd, becq(ik))
107        DO ipol=1,3
108           CALL allocate_bec_type (  nkb, nbnd, alpq(ipol,ik))
109        ENDDO
110     end do
111  endif
112  call compute_weight (wgg)
113  !
114  !    becq and alpq are sufficient to compute the part of C^3 (See Eq. 37
115  !    which does not contain the local potential
116  !
117  IF (.not.lgamma) call compute_becalp (becq, alpq)
118  call compute_nldyn (dyn00, wgg, becq, alpq)
119  !
120  !   now we compute the change of the charge density due to the change of
121  !   the orthogonality constraint
122  !
123  allocate (drhous ( dfftp%nnr, nspin_mag , 3 * nat))
124  allocate (dbecsum( nhm * (nhm + 1) /2, nat, nspin_mag, 3 * nat))
125  dbecsum=(0.d0,0.d0)
126  IF (noncolin) THEN
127     allocate (dbecsum_nc( nhm, nhm, nat, nspin, 3 * nat))
128     dbecsum_nc=(0.d0,0.d0)
129     call compute_drhous_nc (drhous, dbecsum_nc, wgg, becq, alpq)
130  ELSE
131     call compute_drhous (drhous, dbecsum, wgg, becq, alpq)
132  ENDIF
133
134  if (.not.lgamma) then
135     do ik=1,nksq
136        call deallocate_bec_type(becq(ik))
137        DO ipol=1,3
138           call deallocate_bec_type(alpq(ipol,ik))
139        ENDDO
140     end do
141     deallocate (becq)
142     deallocate (alpq)
143  endif
144  deallocate (wgg)
145  !
146  !  The part of C^3 (Eq. 37) which contain the local potential can be
147  !  evaluated with an integral of this change of potential and drhous
148  !
149  allocate (dvlocin(dffts%nnr))
150
151  wdyn (:,:) = (0.d0, 0.d0)
152  nrstot = dffts%nr1 * dffts%nr2 * dffts%nr3
153  do nu_i = 1, 3 * nat
154     call compute_dvloc (nu_i, dvlocin)
155     do nu_j = 1, 3 * nat
156        do is = 1, nspin_lsda
157        ! FIXME: use zgemm instead of dot_product
158           wdyn (nu_j, nu_i) = wdyn (nu_j, nu_i) + &
159                dot_product (drhous(1:dffts%nnr,is,nu_j), dvlocin) * &
160                omega / DBLE (nrstot)
161        enddo
162     enddo
163
164  enddo
165  !
166  ! collect contributions from all pools (sum over k-points)
167  !
168  call mp_sum ( dyn00, inter_pool_comm )
169  call mp_sum ( wdyn, inter_pool_comm )
170  !
171  ! collect contributions from nodes of a pool (sum over G & R space)
172  !
173  call mp_sum ( wdyn, intra_bgrp_comm )
174
175  call zaxpy (3 * nat * 3 * nat, (1.d0, 0.d0), wdyn, 1, dyn00, 1)
176  !
177  !     force this term to be hermitean
178  !
179  do nu_i = 1, 3 * nat
180     do nu_j = 1, nu_i
181        dyn00(nu_i,nu_j) = 0.5d0*( dyn00(nu_i,nu_j) + CONJG(dyn00(nu_j,nu_i)))
182        dyn00(nu_j,nu_i) = CONJG(dyn00(nu_i,nu_j))
183     enddo
184  enddo
185  !      call tra_write_matrix('drho dyn00',dyn00,u,nat)
186  !
187  !    add the augmentation term to the charge density and save it
188  !
189  allocate (drhoust(dfftp%nnr, nspin_mag , npertx))
190  drhoust=(0.d0,0.d0)
191  !
192  !  The calculation of dbecsum is distributed across processors (see addusdbec)
193  !  Sum over processors the contributions coming from each slice of bands
194  !
195  IF (noncolin) THEN
196     call mp_sum ( dbecsum_nc, intra_bgrp_comm )
197  ELSE
198     call mp_sum ( dbecsum, intra_bgrp_comm )
199  END IF
200
201  IF (noncolin.and.okvan) CALL set_dbecsum_nc(dbecsum_nc, dbecsum, 3*nat)
202
203  mode = 0
204  if (okpaw) becsumort=(0.0_DP,0.0_DP)
205  do irr = 1, nirr
206     npe = npert (irr)
207     if (doublegrid) then
208        do is = 1, nspin_mag
209           do iper = 1, npe
210              call fft_interpolate (dffts, drhous(:,is,mode+iper), dfftp, drhoust(:,is,iper))
211           enddo
212        enddo
213     else
214        call zcopy (dfftp%nnr*nspin_mag*npe, drhous(1,1,mode+1), 1, drhoust, 1)
215     endif
216
217     call dscal (2*dfftp%nnr*nspin_mag*npe, 0.5d0, drhoust, 1)
218
219     call addusddens (drhoust, dbecsum(1,1,1,mode+1), mode, npe, 1)
220     do iper = 1, npe
221        nu_i = mode+iper
222        call save_buffer (drhoust (1, 1, iper), lrdrhous, iudrhous, nu_i)
223     enddo
224     mode = mode+npe
225  enddo
226   !
227   !  Collect the sum over k points in different pools.
228   !
229   IF (okpaw) call mp_sum ( becsumort, inter_pool_comm )
230
231  deallocate (drhoust)
232  deallocate (dvlocin)
233  deallocate (dbecsum)
234  if (noncolin) deallocate(dbecsum_nc)
235  deallocate (drhous)
236
237  call stop_clock ('drho')
238  return
239end subroutine drho
240