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