TABLE OF CONTENTS
ABINIT/integrate_gamma_alt [ Functions ]
NAME
integrate_gamma_alt
FUNCTION
This routine integrates the electron phonon coupling matrix over the kpoints on the fermi surface. A dependency on qpoint remains for gamma_qpt
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
elph_ds = elphon datastructure with data and dimensions elph_ds%qpt_full = qpoint coordinates elph_ds%nqptirred = number of irred qpoints elph_ds%qirredtofull = indexing of the GKK qpoints found elph_tr_ds = elphon transport datastructure elph_tr_ds nrpt = number of real space points for FT
OUTPUT
elph_ds = modified elph_ds%gamma_qpt with much larger size elph_tr_ds = modified elph_ds%gamma_qpt with much larger size
PARENTS
CHILDREN
complete_gamma,ftgam,ftgam_init,get_rank_1kpt,k_neighbors,wrtout
SOURCE
39 #if defined HAVE_CONFIG_H 40 #include "config.h" 41 #endif 42 43 #include "abi_common.h" 44 45 subroutine integrate_gamma_alt(elph_ds,elph_tr_ds,Cryst,gprim,kptrlatt,& 46 & natom,nrpt,nsym,qpttoqpt,rpt,wghatm) 47 48 use defs_basis 49 use defs_elphon 50 use m_profiling_abi 51 use m_errors 52 use m_kptrank 53 use m_xmpi 54 55 use m_dynmat, only : ftgam_init, ftgam 56 use m_numeric_tools, only : interpol3d 57 use m_crystal, only : crystal_t 58 59 !This section has been created automatically by the script Abilint (TD). 60 !Do not modify the following lines by hand. 61 #undef ABI_FUNC 62 #define ABI_FUNC 'integrate_gamma_alt' 63 use interfaces_14_hidewrite 64 use interfaces_77_ddb, except_this_one => integrate_gamma_alt 65 !End of the abilint section 66 67 implicit none 68 69 !Arguments ------------------------------------ 70 !scalars 71 integer, intent(in) :: nrpt 72 integer, intent(in) :: natom, nsym 73 type(crystal_t),intent(in) :: Cryst 74 type(elph_type),intent(inout) :: elph_ds 75 type(elph_tr_type), intent(inout) :: elph_tr_ds 76 !arrays 77 integer,intent(in) :: qpttoqpt(2,nsym,elph_ds%nqpt_full) 78 real(dp), intent(in) :: wghatm(natom,natom,nrpt), rpt(3,nrpt) 79 real(dp), intent(in) :: gprim(3,3) 80 integer, intent(in) :: kptrlatt(3,3) 81 82 !Local variables------------------------------- 83 !scalars 84 integer :: ikpt_phon,ib1,ib2,ibeff,iqpt,iqpt_fullbz,isppol 85 integer :: irec, qtor, rtoq, ibranch,nbranch,ngkkband,ierr 86 integer :: symrankkpt 87 integer :: ikptfine,ikptqfine 88 integer :: ineighbor 89 integer :: counter_fine 90 integer :: fib1, fib2 91 integer :: itensor, icomp, jcomp,sz1,sz2 92 real(dp) :: sd1, sd2, tolfs 93 94 !arrays 95 integer :: kpt_phon_indices (8) 96 real(dp) :: etaout, etain 97 real(dp) :: qpt(3), rel_kpt(3) 98 real(dp) :: kpt(3) 99 real(dp) :: gam_now(2,elph_ds%nbranch*elph_ds%nbranch) 100 real(dp) :: elvelock(3), elvelockpq(3) 101 102 character(len=500) :: message 103 character(len=fnlen) :: fname 104 105 real(dp),allocatable :: tmp_gkk(:,:,:,:,:) 106 real(dp),allocatable :: tmp_gkk_1q(:,:,:,:,:) 107 real(dp),allocatable :: tmp_gkr(:,:,:) 108 real(dp),allocatable :: coskr_full(:,:) 109 real(dp),allocatable :: sinkr_full(:,:) 110 real(dp),allocatable :: coskr_fine(:,:) 111 real(dp),allocatable :: sinkr_fine(:,:) 112 113 ! ************************************************************************* 114 115 write (message,'(3a)')ch10,' entering integrate_gamma_alt ',ch10 116 call wrtout(std_out,message,'COLL') 117 118 if (xmpi_comm_size(xmpi_world) > 1) then 119 write (message, '(a)') ' alternative integration scheme does not work in parallel yet' 120 MSG_ERROR(message) 121 end if 122 123 tolfs = 1.e-11 124 125 if (elph_ds%ep_keepbands == 0) then 126 write (message,'(3a)')& 127 & ' Cannot call integrate_gamma_alt without all bands in memory',ch10,& 128 & ' ACTION: set ep_keepbands = 1 in anaddb input file' 129 MSG_ERROR(message) 130 end if 131 132 nbranch = elph_ds%nbranch 133 ngkkband = elph_ds%ngkkband 134 135 !allocations 136 ABI_STAT_ALLOCATE(elph_ds%gamma_qpt,(2,nbranch**2,elph_ds%nsppol,elph_ds%k_fine%nkpt), ierr) 137 ABI_CHECK(ierr==0, 'integrate_gamma_alt trying to allocate array elph_ds%gamma_qpt') 138 elph_ds%gamma_qpt(:,:,:,:) = zero 139 140 !transport case 141 if (elph_tr_ds%ifltransport==1) then 142 ABI_STAT_ALLOCATE(elph_tr_ds%gamma_qpt_trin,(2,9,nbranch**2,elph_ds%nsppol,elph_ds%k_fine%nkpt), ierr) 143 ABI_CHECK(ierr==0, 'trying to allocate array elph_tr_ds%gamma_qpt_trin') 144 elph_tr_ds%gamma_qpt_trin = zero 145 146 ABI_STAT_ALLOCATE(elph_tr_ds%gamma_qpt_trout,(2,9,nbranch**2,elph_ds%nsppol,elph_ds%k_fine%nkpt), ierr) 147 ABI_CHECK(ierr==0, "trying to allocate array elph_tr_ds%gamma_qpt_trout") 148 elph_tr_ds%gamma_qpt_trout = zero 149 end if 150 151 152 !temp variables 153 ABI_STAT_ALLOCATE(tmp_gkk_1q ,(2,ngkkband**2,nbranch**2,8,elph_ds%nsppol), ierr) 154 ABI_CHECK(ierr==0, 'trying to allocate tmp_gkk_1q') 155 156 sz1=elph_ds%ngkkband*elph_ds%ngkkband 157 sz2=elph_ds%nbranch*elph_ds%nbranch 158 ABI_STAT_ALLOCATE(tmp_gkk ,(2,sz1,sz2,elph_ds%nsppol,elph_ds%nqpt_full), ierr) 159 ABI_CHECK(ierr==0, 'trying to allocate tmp_gkk') 160 161 ABI_STAT_ALLOCATE(tmp_gkr ,(2,elph_ds%nbranch*elph_ds%nbranch,nrpt), ierr) 162 ABI_CHECK(ierr==0, 'trying to allocate array tmp_gkr') 163 164 if (elph_ds%gkqwrite == 0) then 165 call wrtout(std_out,' integrate_gamma_alt : getting gamma matrices from memory','COLL') 166 else if (elph_ds%gkqwrite == 1) then 167 fname=trim(elph_ds%elph_base_name) // '_GKKQ' 168 write (message,'(2a)')' integrate_gamma_alt : reading gamma matrices from file ',trim(fname) 169 call wrtout(std_out,message,'COLL') 170 else 171 write (message,'(a,i0)')' Wrong value for gkqwrite = ',elph_ds%gkqwrite 172 MSG_BUG(message) 173 end if 174 175 qtor = 1 ! q --> r 176 rtoq = 0 ! r --> q 177 178 counter_fine = 0 179 180 !FIXME: group fine kpt in batches with the same corners in sparse k-grid, to do only 1 181 !retrieval of gkk for each batch. 182 do ikptfine=1,elph_ds%k_fine%nkpt 183 ! if we have found a good k-point on the FS 184 if (sum(abs(elph_ds%k_fine%wtk(:,ikptfine,:))) < tolfs) cycle 185 counter_fine = counter_fine+1 186 187 call k_neighbors (elph_ds%k_fine%kpt(:,ikptfine), kptrlatt, elph_ds%k_phon%kptrank_t, & 188 & rel_kpt, kpt_phon_indices) 189 190 ! FIXME: this read in of 8 matrices is done several times for fine kpt inside the same phon_kpt cell 191 ! more generally could keep 4 matrices when moving to the neighboring cell, but bookkeeping 192 ! will be a pain. 193 tmp_gkk = zero 194 do iqpt=1,elph_ds%nqptirred 195 do ineighbor = 1, 8 196 ikpt_phon = kpt_phon_indices(ineighbor) 197 if (elph_ds%gkqwrite == 0) then 198 tmp_gkk_1q(:,:,:,ineighbor,:) = elph_ds%gkk_qpt(:,:,:,ikpt_phon,:,iqpt) 199 else if (elph_ds%gkqwrite == 1) then 200 irec = (iqpt-1)*elph_ds%k_phon%nkpt+ikpt_phon 201 read (elph_ds%unitgkq,REC=irec) tmp_gkk_1q(:,:,:,ineighbor,:) 202 end if 203 end do ! neighbor ikpt_phon 204 205 ! interpolate gkk matrix elements wrt k vector 206 ! FIXME: make a batch version of interpol3d for all matrix elements at once 207 do isppol = 1, elph_ds%nsppol 208 do ibranch = 1, elph_ds%nbranch*elph_ds%nbranch 209 do ibeff = 1, elph_ds%ngkkband*elph_ds%ngkkband 210 211 tmp_gkk(1,ibeff,ibranch,isppol,elph_ds%qirredtofull(iqpt)) = & 212 & interpol3d(rel_kpt,2,2,2,tmp_gkk_1q(1,ibeff,ibranch,:,isppol)) 213 214 tmp_gkk(2,ibeff,ibranch,isppol,elph_ds%qirredtofull(iqpt)) = & 215 & interpol3d(rel_kpt,2,2,2,tmp_gkk_1q(2,ibeff,ibranch,:,isppol)) 216 end do 217 end do 218 end do 219 end do ! do iqpt 220 ! now tmpgkk is filled for present fine kpt and all irred qpoints on sparse grid 221 222 ABI_ALLOCATE(coskr_full, (elph_ds%nqpt_full,nrpt)) 223 ABI_ALLOCATE(sinkr_full, (elph_ds%nqpt_full,nrpt)) 224 call ftgam_init(gprim, elph_ds%nqpt_full,nrpt, elph_ds%qpt_full, rpt, coskr_full, sinkr_full) 225 226 ABI_ALLOCATE(coskr_fine, (elph_ds%k_fine%nkpt,nrpt)) 227 ABI_ALLOCATE(sinkr_fine, (elph_ds%k_fine%nkpt,nrpt)) 228 call ftgam_init(gprim, elph_ds%k_fine%nkpt, nrpt, elph_ds%k_fine%kpt, rpt, coskr_fine, sinkr_fine) 229 230 231 ! find spin and band for FS contributing state 232 do isppol=1,elph_ds%nsppol 233 do ib1=1,elph_ds%ngkkband 234 if (abs(elph_ds%k_fine%wtk(ib1,ikptfine,isppol)) < tolfs) cycle 235 sd1 = elph_ds%k_fine%wtk(ib1,ikptfine,isppol) 236 237 if (elph_tr_ds%ifltransport==1) then 238 fib1 = ib1 + elph_ds%minFSband - 1 239 elvelock(:)=elph_tr_ds%el_veloc(ikptfine,fib1,:,isppol) 240 end if 241 242 ! loop over k+q, n' point 243 ! FIXME: check if we need a second spin index here... looks like there is none in abinit matrix elements... 244 ! does this presume no SO operator in the perturbed potential? Otherwise the deltaV is spin diagonal... 245 do ib2=1,elph_ds%ngkkband 246 ibeff = ib2+(ib1-1)*elph_ds%ngkkband 247 248 ! do FT of gkq matrix to gkR 249 if (elph_ds%symgkq ==1) then 250 ! FIXME: this does both sppol by default, which is a waste here inside the loop on isppol 251 call complete_gamma(Cryst,elph_ds%nbranch,elph_ds%nsppol,elph_ds%nqptirred,elph_ds%nqpt_full,& 252 & elph_ds%ep_scalprod,elph_ds%qirredtofull,qpttoqpt,tmp_gkk(:,ibeff,:,:,:)) 253 end if 254 255 ! Now FT to real space too 256 call ftgam(wghatm,tmp_gkk(:,ibeff,:,isppol,:),tmp_gkr(:,:,:),natom,& 257 & elph_ds%nqpt_full,nrpt,qtor, coskr_full, sinkr_full) 258 259 ! run _qpt_ over the dense grid of kpt_fine 260 ! FIXME: could at least use time reversal sym, and probably more, to do fewer loop iterations here 261 do iqpt_fullbz = 1, elph_ds%k_fine%nkpt 262 263 ! find bloody index of k + q 264 kpt(:) = elph_ds%k_fine%kpt(:,ikptfine) + elph_ds%k_fine%kpt(:,iqpt_fullbz) 265 ! NB: this indexing is not the same order as the normal full qpts (elph_ds%k_fine%kpt(:,ikpt_fine) + qpt(:,iqpt_fullbz) 266 ! and it runs over all kpt_fine, not just the sparse input qpts. 267 268 ! which kpt is k+q among the full FS kpts 269 call get_rank_1kpt (kpt,symrankkpt,elph_ds%k_fine%kptrank_t) 270 ikptqfine = elph_ds%k_fine%kptrank_t%invrank(symrankkpt) 271 272 if (elph_tr_ds%ifltransport==1) then 273 fib2 = ib2 + elph_ds%minFSband - 1 274 elvelockpq(:)=elph_tr_ds%el_veloc(ikptqfine,fib2,:,isppol) 275 end if 276 277 278 if (abs(elph_ds%k_fine%wtk(ib2,ikptqfine,isppol)) < tolfs) cycle 279 sd2 = elph_ds%k_fine%wtk(ib2,ikptqfine,isppol) 280 qpt = elph_ds%k_fine%kpt(:,iqpt_fullbz) 281 282 ! interpolate gkR to required gkq 283 call ftgam(wghatm,gam_now,tmp_gkr(:,:,:),natom,1,nrpt,rtoq, & 284 & coskr_fine(iqpt_fullbz,:), sinkr_fine(iqpt_fullbz,:) ) 285 286 ! add to gamma(q) = sum over k and bands 287 elph_ds%gamma_qpt(:,:,isppol,iqpt_fullbz) = elph_ds%gamma_qpt(:,:,isppol,iqpt_fullbz) & 288 & + gam_now*sd1*sd2 289 290 if (elph_tr_ds%ifltransport==1) then 291 do icomp = 1, 3 292 do jcomp = 1, 3 293 itensor = (icomp-1)*3+jcomp 294 ! FIXME: could use symmetry i <-> j 295 296 etain = elvelock(icomp)*elvelockpq(jcomp) 297 etaout = elvelock(icomp)*elvelock(jcomp) 298 299 elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,iqpt_fullbz) = & 300 & elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,iqpt_fullbz) + & 301 & gam_now(:,:)*etain*sd1*sd2 302 303 elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,iqpt_fullbz) = & 304 & elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,iqpt_fullbz) + & 305 & gam_now(:,:)*etaout*sd1*sd2 306 end do ! jcomp 307 end do ! icomp 308 end if 309 310 end do !iqpt_fullbz 311 end do !ib2 312 end do !ib1 313 end do !isppol 314 end do !ikpt_fine 315 316 ABI_DEALLOCATE(tmp_gkk) 317 ABI_DEALLOCATE(tmp_gkr) 318 319 write (message,'(a,E20.10)') 'fraction of fine kpt used : ', dble(counter_fine) / elph_ds%k_fine%nkpt 320 call wrtout(std_out,message,'COLL') 321 322 !need prefactor of 1/nkpt for each integration over 1 kpoint index. 323 !NOT INCLUDED IN elph_ds%k_fine%wtk 324 elph_ds%gamma_qpt = elph_ds%gamma_qpt * elph_ds%occ_factor / elph_ds%k_fine%nkpt 325 elph_tr_ds%gamma_qpt_trin = elph_tr_ds%gamma_qpt_trin * elph_ds%occ_factor / elph_ds%k_fine%nkpt 326 elph_tr_ds%gamma_qpt_trout = elph_tr_ds%gamma_qpt_trout * elph_ds%occ_factor / elph_ds%k_fine%nkpt 327 328 !normalization in transport case 329 if (elph_tr_ds%ifltransport==1) then 330 do isppol = 1, elph_ds%nsppol 331 do icomp = 1, 3 332 do jcomp = 1, 3 333 itensor = (icomp-1)*3+jcomp 334 elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,:) = elph_tr_ds%gamma_qpt_trin(:,itensor,:,isppol,:) / & 335 & sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol)) 336 elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,:) = elph_tr_ds%gamma_qpt_trout(:,itensor,:,isppol,:) / & 337 & sqrt(elph_tr_ds%FSelecveloc_sq(icomp,isppol)*elph_tr_ds%FSelecveloc_sq(jcomp,isppol)) 338 end do 339 end do 340 end do 341 end if 342 343 ABI_DEALLOCATE(coskr_full) 344 ABI_DEALLOCATE(sinkr_full) 345 ABI_DEALLOCATE(coskr_fine) 346 ABI_DEALLOCATE(sinkr_fine) 347 348 write (message,'(a,a)') ' integrate_gamma_alt : gamma matrices have been calculated',& 349 & 'for recip space and irred qpoints ' 350 call wrtout(std_out,message,'COLL') 351 352 end subroutine integrate_gamma_alt