TABLE OF CONTENTS
- ABINIT/m_ephtk
- m_ephtk/ephtk_gam_atm2qnu
- m_ephtk/ephtk_gkknu_from_atm
- m_ephtk/ephtk_mkqtabs
- m_ephtk/ephtk_set_pertables
- m_ephtk/ephtk_set_phmodes_skip
- m_ephtk/ephtk_update_ebands
ABINIT/m_ephtk [ Modules ]
NAME
m_ephtk
FUNCTION
Helper functions common to e-ph calculations.
COPYRIGHT
Copyright (C) 2008-2022 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_ephtk 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 use m_dtset 29 use m_ebands 30 use m_crystal 31 use m_krank 32 use m_xmpi 33 34 use m_fstrings, only : itoa, sjoin, ltoa, ftoa, ktoa 35 use m_bz_mesh, only : isamek 36 use defs_datatypes, only : ebands_t 37 38 implicit none 39 40 private 41 42 public :: ephtk_set_phmodes_skip ! Setup a mask to skip accumulating the contribution of certain phonon modes. 43 public :: ephtk_set_pertables ! Set tables for parallelism over perturbations from my_npert and comm 44 public :: ephtk_mkqtabs ! Build tables with correspondence between q-points as needed by complete_gamma. 45 public :: ephtk_gam_atm2qnu ! Compute phonon linewidths from gamma matrix in reduced coordinates. 46 public :: ephtk_gkknu_from_atm ! Transform the gkk matrix elements from (atom, red_direction) basis to phonon-mode basis. 47 public :: ephtk_update_ebands ! Update ebands according to dtset%occopt, tsmear, mbpt_sciss, eph_fermie, eph_extrael
m_ephtk/ephtk_gam_atm2qnu [ Functions ]
[ Top ] [ m_ephtk ] [ Functions ]
NAME
ephtk_gam_atm2qnu
FUNCTION
This routine takes the gamma matrices in the atomic representation and multiplies them by the displ_red matrices. Based on gam_mult_displ
INPUTS
natom3 = number of phonon branches (3*natom) displ_red = phonon mode displacement vectors in reduced coordinates. gam_bare = bare gamma matrices before multiplication
OUTPUT
gam_now = output gamma matrices multiplied by displacement matrices
SOURCE
286 subroutine ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, gam_qnu) 287 288 !Arguments ------------------------------- 289 integer, intent(in) :: natom3 290 real(dp), intent(in) :: displ_red(2,natom3,natom3) 291 real(dp), intent(in) :: gam_atm(2,natom3,natom3) 292 real(dp), intent(out) :: gam_qnu(natom3) 293 294 !Local variables ------------------------- 295 integer,save :: enough = 0 296 integer :: nu 297 character(len=500) :: msg 298 real(dp) :: zgemm_tmp_mat(2,natom3,natom3), gam_now(2,natom3,natom3) 299 300 ! ********************************************************************* 301 302 call zgemm('c','n',natom3, natom3, natom3, cone, displ_red, natom3, gam_atm, natom3, czero, zgemm_tmp_mat, natom3) 303 304 gam_now = zero 305 call zgemm('n','n',natom3,natom3,natom3,cone,zgemm_tmp_mat,natom3,displ_red,natom3,czero,gam_now,natom3) 306 307 ! Extract gamma(q,nu) 308 do nu=1,natom3 309 gam_qnu(nu) = gam_now(1, nu, nu) 310 if (abs(gam_now(2, nu, nu)) > tol8) then 311 enough = enough + 1 312 if (enough <= 30) then 313 write (msg,'(a,i0,a,es16.8)')' non-zero imaginary part for branch: ',nu,', img: ',gam_now(2, nu, nu) 314 ABI_WARNING(msg) 315 end if 316 end if 317 end do 318 319 end subroutine ephtk_gam_atm2qnu
m_ephtk/ephtk_gkknu_from_atm [ Functions ]
[ Top ] [ m_ephtk ] [ Functions ]
NAME
ephtk_gkknu_from_atm
FUNCTION
Transform the gkk matrix elements from (atom, red_direction) basis to phonon-mode basis.
INPUTS
nb1,nb2=Number of bands in gkq_atm matrix. nk=Number of k-points (usually 1) natom=Number of atoms. gkq_atm(2,nb1,nb2,3*natom)=EPH matrix elements in the atomic basis. phfrq(3*natom)=Phonon frequencies in Ha displ_red(2,3*natom,3*natom)=Phonon displacement in reduced coordinates.
OUTPUT
gkq_nu(2,nb1,nb2,3*natom)=EPH matrix elements in the phonon-mode basis.
SOURCE
344 subroutine ephtk_gkknu_from_atm(nb1, nb2, nk, natom, gkq_atm, phfrq, displ_red, gkq_nu) 345 346 !Arguments ------------------------------------ 347 !scalars 348 integer,intent(in) :: nb1, nb2, nk, natom 349 !arrays 350 real(dp),intent(in) :: phfrq(3*natom),displ_red(2,3*natom,3*natom) 351 real(dp),intent(in) :: gkq_atm(2,nb1,nb2,nk,3*natom) 352 real(dp),intent(out) :: gkq_nu(2,nb1,nb2,nk,3*natom) 353 354 !Local variables------------------------- 355 !scalars 356 integer :: nu,ipc 357 358 ! ************************************************************************* 359 360 gkq_nu = zero 361 362 ! Loop over phonon branches. 363 do nu=1,3*natom 364 ! Ignore negative or too small frequencies 365 if (phfrq(nu) < EPHTK_WTOL) cycle 366 367 ! Transform the gkk from (atom, reduced direction) basis to phonon mode representation 368 do ipc=1,3*natom 369 gkq_nu(1,:,:,:,nu) = gkq_nu(1,:,:,:,nu) & 370 + gkq_atm(1,:,:,:,ipc) * displ_red(1,ipc,nu) & 371 - gkq_atm(2,:,:,:,ipc) * displ_red(2,ipc,nu) 372 gkq_nu(2,:,:,:,nu) = gkq_nu(2,:,:,:,nu) & 373 + gkq_atm(1,:,:,:,ipc) * displ_red(2,ipc,nu) & 374 + gkq_atm(2,:,:,:,ipc) * displ_red(1,ipc,nu) 375 end do 376 377 gkq_nu(:,:,:,:,nu) = gkq_nu(:,:,:,:,nu) / sqrt(two * phfrq(nu)) 378 end do 379 380 end subroutine ephtk_gkknu_from_atm
m_ephtk/ephtk_mkqtabs [ Functions ]
[ Top ] [ m_ephtk ] [ Functions ]
NAME
ephtk_mkqtabs
FUNCTION
Build tables with correspondence between q-points in the IBZ/BZ as needed by complete_gamma. INPUT cryst<crystal_t>=Crystal structure. nqibz, qibz = Points in the IBZ nqbz, qbz = Points in the BZ
OUTPUT
qirredtofull(nqibz) = mapping irred to full qpoints qpttoqpt(2, cryst%nsym, nqbz)) = qpoint index mapping under symops.
SOURCE
203 subroutine ephtk_mkqtabs(cryst, nqibz, qibz, nqbz, qbz, qirredtofull, qpttoqpt) 204 205 !Arguments ------------------------------------ 206 type(crystal_t),intent(in) :: cryst 207 integer,intent(in) :: nqibz, nqbz 208 !arrays 209 real(dp),intent(in) :: qibz(3, nqibz), qbz(3, nqbz) 210 integer,allocatable :: qirredtofull(:),qpttoqpt(:,:,:) 211 212 !Local variables ------------------------------ 213 !scalars 214 integer :: iq_bz, iq_ibz, isq_bz, isym 215 type(krank_t) :: qrank 216 !arrays 217 integer :: g0(3) 218 real(dp) :: qirr(3), tmp_qpt(3) 219 220 ! ************************************************************************* 221 222 qrank = krank_new(nqbz, qbz) 223 224 ! Compute index of IBZ q-point in the BZ array 225 ABI_CALLOC(qirredtofull, (nqibz)) 226 227 do iq_ibz=1,nqibz 228 qirr = qibz(:,iq_ibz) 229 iq_bz = qrank%get_index(qirr) 230 if (iq_bz /= -1) then 231 ABI_CHECK(isamek(qirr, qbz(:,iq_bz), g0), "isamek") 232 qirredtofull(iq_ibz) = iq_bz 233 else 234 ABI_ERROR(sjoin("Full BZ does not contain IBZ q-point:", ktoa(qirr))) 235 end if 236 end do 237 238 ! Build qpttoqpt table. See also mkqptequiv 239 ABI_MALLOC(qpttoqpt, (2, cryst%nsym, nqbz)) 240 qpttoqpt = -1 241 do iq_bz=1,nqbz 242 do isym=1,cryst%nsym 243 tmp_qpt = matmul(cryst%symrec(:,:,isym), qbz(:,iq_bz)) 244 245 isq_bz = qrank%get_index(tmp_qpt) 246 if (isq_bz == -1) then 247 ABI_ERROR("Looks like no kpoint equiv to q by symmetry without time reversal!") 248 end if 249 qpttoqpt(1,isym,isq_bz) = iq_bz 250 251 ! q --> -q 252 tmp_qpt = -tmp_qpt 253 isq_bz = qrank%get_index(tmp_qpt) 254 if (isq_bz == -1) then 255 ABI_ERROR("Looks like no kpoint equiv to q by symmetry with time reversal!") 256 end if 257 qpttoqpt(2,isym,isq_bz) = iq_bz 258 end do 259 end do 260 261 call qrank%free() 262 263 end subroutine ephtk_mkqtabs
m_ephtk/ephtk_set_pertables [ Functions ]
[ Top ] [ m_ephtk ] [ Functions ]
NAME
ephtk_set_pertables
FUNCTION
Build tables for parallelism over perturbations from my_npert and comm INPUT natom: Number of atoms my_npert: Number of atomic perturbations or phonon modes treated by this MPI rank. comm: MPI communicator for parallelism over atomic perturbations.
OUTPUT
integer,allocatable :: my_pinfo(:,:) my_pinfo(3, my_npert) my_pinfo(1, ip) gives the `idir` index of the ip-th perturbation. my_pinfo(2, ip) gives the `ipert` index of the ip-th perturbation. my_pinfo(3, ip) gives `pertcase`=idir + (ipert-1)*3 integer,allocatable :: pert_table(:,:) pert_table(2, natom3) pert_table(1, npert): rank of the processor treating this atomic perturbation. pert_table(2, npert): imyp index in my_pinfo table, -1 if this rank is not treating ipert.
SOURCE
142 subroutine ephtk_set_pertables(natom, my_npert, pert_table, my_pinfo, comm) 143 144 !Arguments ------------------------------------ 145 integer,intent(in) :: natom, my_npert, comm 146 !arrays 147 integer,allocatable :: pert_table(:,:), my_pinfo(:,:) 148 149 !Local variables ------------------------------ 150 !scalars 151 integer :: iatom, idir, pertcase, bstart, bstop, ii, ip, natom3, my_rank, nproc 152 !arrays 153 integer :: all_pinfo(3, natom*3) 154 155 ! ************************************************************************* 156 157 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 158 159 ! Build table with list of perturbations treated by this CPU. 160 natom3 = natom * 3 161 ABI_MALLOC(my_pinfo, (3, my_npert)) 162 ABI_MALLOC(pert_table, (2, natom3)) 163 164 do iatom=1,natom 165 do idir=1,3 166 pertcase = idir + (iatom-1) * 3 167 all_pinfo(:, pertcase) = [idir, iatom, pertcase] 168 pert_table(1, pertcase) = (pertcase - 1) / (natom3 / nproc) 169 end do 170 end do 171 bstart = (natom3 / nproc) * my_rank + 1 172 bstop = bstart + my_npert - 1 173 my_pinfo = all_pinfo(:, bstart:bstop) 174 175 pert_table(2, :) = -1 176 do ii=1,my_npert 177 ip = my_pinfo(3, ii) 178 pert_table(2, ip) = ii 179 end do 180 !write(std_out,*)"my_npert", my_npert, "nproc", nproc; write(std_out,*)"my_pinfo", my_pinfo 181 182 end subroutine ephtk_set_pertables
m_ephtk/ephtk_set_phmodes_skip [ Functions ]
[ Top ] [ m_ephtk ] [ Functions ]
NAME
ephtk_set_phmodes_skip
FUNCTION
Setup a mask to skip accumulating the contribution of certain phonon modes. INPUT eph_phrange=Abinit input variable.
OUTPUT
phmodes_skip(natom3) For each mode: 1 to skip the contribution given by this phonon branch else 0
SOURCE
73 subroutine ephtk_set_phmodes_skip(natom, eph_phrange, phmodes_skip) 74 75 !Arguments ------------------------------------ 76 integer,intent(in) :: natom 77 !arrays 78 integer,intent(in) :: eph_phrange(2) 79 integer,allocatable,intent(out) :: phmodes_skip(:) 80 81 !Local variables ------------------------------ 82 !scalars 83 integer :: natom3 84 85 ! ************************************************************************* 86 87 ! Setup a mask to skip accumulating the contribution of certain phonon modes. 88 ! By default do not skip, if set skip all but specified 89 natom3 = natom * 3 90 ABI_MALLOC(phmodes_skip, (natom3)) 91 phmodes_skip = 0 92 93 if (all(eph_phrange /= 0)) then 94 if (minval(abs(eph_phrange)) < 1 .or. & 95 maxval(abs(eph_phrange)) > natom3 .or. & 96 abs(eph_phrange(2)) < abs(eph_phrange(1))) then 97 ABI_ERROR('Invalid range for eph_phrange. Should be between [1, 3*natom] and eph_modes(2) > eph_modes(1)') 98 end if 99 if (all(eph_phrange > 0)) then 100 call wrtout(std_out, sjoin(" Including phonon modes between [", & 101 itoa(eph_phrange(1)), ',', itoa(eph_phrange(2)), "]")) 102 phmodes_skip = 1 103 phmodes_skip(eph_phrange(1):eph_phrange(2)) = 0 104 else if (all(eph_phrange < 0)) then 105 call wrtout(std_out, sjoin(" Excluding phonon modes between [", & 106 itoa(abs(eph_phrange(1))), ',', itoa(abs(eph_phrange(2))), "]")) 107 phmodes_skip = 0 108 phmodes_skip(abs(eph_phrange(1)):abs(eph_phrange(2))) = 1 109 else 110 ABI_ERROR(sjoin("Invalid eph_phrange: ", itoa(eph_phrange(1)), ',', itoa(eph_phrange(2)))) 111 end if 112 end if 113 114 end subroutine ephtk_set_phmodes_skip
m_ephtk/ephtk_update_ebands [ Functions ]
[ Top ] [ m_ephtk ] [ Functions ]
NAME
ephtk_update_ebands
FUNCTION
Update ebands according to dtset%occopt, tsmear, mbpt_sciss, eph_fermie, eph_extrael
INPUTS
dtset<dataset_type>=All input variables for this dataset.
SOURCE
397 subroutine ephtk_update_ebands(dtset, ebands, header) 398 399 !Arguments ------------------------------------ 400 !scalars 401 type(dataset_type),intent(in) :: dtset 402 type(ebands_t),intent(inout) :: ebands 403 character(len=*),intent(in) :: header 404 405 !Local variables------------------------- 406 !scalars 407 real(dp),parameter :: nholes = zero 408 character(len=500) :: msg 409 integer :: unts(2) 410 411 ! ************************************************************************* 412 413 unts = [std_out, ab_out] 414 415 if (dtset%occopt /= ebands%occopt .or. abs(dtset%tsmear - ebands%tsmear) > tol12) then 416 !if (.True.) then 417 write(msg,"(2a,2(a,i0,a,f14.6,a))")& 418 " Changing occupation scheme as input occopt and tsmear differ from those read from WFK file.",ch10,& 419 " From WFK file: occopt = ",ebands%occopt,", tsmear = ",ebands%tsmear,ch10,& 420 " From input: occopt = ",dtset%occopt,", tsmear = ",dtset%tsmear,ch10 421 call wrtout(unts, msg) 422 call ebands_set_scheme(ebands, dtset%occopt, dtset%tsmear, dtset%spinmagntarget, dtset%prtvol) 423 424 if (abs(dtset%mbpt_sciss) > tol6) then 425 ! Apply the scissor operator 426 call wrtout(unts, sjoin(" Applying scissors operator to the conduction states with value: ", & 427 ftoa(dtset%mbpt_sciss * Ha_eV, fmt="(f6.2)"), " (eV)")) 428 call ebands_apply_scissors(ebands, dtset%mbpt_sciss) 429 end if 430 end if 431 432 ! Default value of eph_fermie is zero hence no tolerance is used! 433 if (dtset%eph_fermie /= zero) then 434 ABI_CHECK(dtset%eph_extrael == zero, "eph_fermie and eph_extrael are mutually exclusive") 435 call wrtout(unts, sjoin(" Fermi level set by the user at:", ftoa(dtset%eph_fermie * Ha_eV, fmt="(f6.2)"), " (eV)")) 436 call ebands_set_fermie(ebands, dtset%eph_fermie, msg) 437 call wrtout(unts, msg) 438 439 else if (abs(dtset%eph_extrael) > zero) then 440 call wrtout(unts, sjoin(" Adding eph_extrael:", ftoa(dtset%eph_extrael), "to input nelect:", ftoa(ebands%nelect))) 441 call ebands_set_scheme(ebands, dtset%occopt, dtset%tsmear, dtset%spinmagntarget, dtset%prtvol, update_occ=.False.) 442 call ebands_set_extrael(ebands, dtset%eph_extrael, nholes, dtset%spinmagntarget, msg) 443 call wrtout(unts, msg) 444 end if 445 446 ! Recompute occupations. This is needed if WFK files have been produced in a NSCF run 447 ! since occ are set to zero, and fermie is taken from the previous density. 448 if (dtset%kptopt > 0) then 449 call ebands_update_occ(ebands, dtset%spinmagntarget, prtvol=dtset%prtvol) 450 call ebands_print(ebands, header=header, prtvol=dtset%prtvol) 451 end if 452 453 end subroutine ephtk_update_ebands