TABLE OF CONTENTS


ABINIT/m_ephtk [ Modules ]

[ Top ] [ 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