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