TABLE OF CONTENTS

ABINIT/m_ephwg [ Modules ]

[ Top ] [ Modules ]

NAME

 m_ephwg

FUNCTION

  Tools and objects to compute the weights used for the BZ integration of EPH quantities.
  More specifically the integration of quantities such as the imaginary part of the self-energy
  involving delta functions. Different approaches are available.

COPYRIGHT

  Copyright (C) 2008-2024 ABINIT group (MG, HM)
  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

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 module m_ephwg
25 
26  use defs_basis
27  use m_abicore
28  use m_errors
29  use m_xmpi
30  use m_copy
31  use m_dtset
32  use m_htetra
33  use m_nctk
34 #ifdef HAVE_NETCDF
35  use netcdf
36 #endif
37  use m_crystal
38  use m_ifc
39  use m_lgroup
40  use m_ebands
41  use m_eph_double_grid
42  use m_krank
43 
44  use defs_datatypes,    only : ebands_t
45  use m_time,            only : cwtime, cwtime_report
46  use m_symtk,           only : matr3inv
47  use m_numeric_tools,   only : arth, inrange, wrap2_pmhalf
48  use m_special_funcs,   only : gaussian
49  use m_fstrings,        only : strcat, ltoa, itoa, ftoa, ktoa, sjoin
50  use m_simtet,          only : sim0onei, SIM0TWOI
51  use m_kpts,            only : kpts_timrev_from_kptopt, kpts_ibz_from_kptrlatt, kpts_map
52  use m_occ,             only : occ_fd, occ_be
53 
54  implicit none
55 
56  private

m_ephwg/ephwg_double_grid_setup_kpoint [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_setup_kpoint

FUNCTION

  Set internal tables and object required to compute integration weights for a given k-point
  using the double grid routines to map the different k-points.
  This version should be more efficient than its counterpart without the double grid.

INPUTS

  kpoint(3): k-point in reduced coordinates.
  prtvol: Verbosity level
  comm: MPI communicator

OUTPUT

SOURCE

457 subroutine ephwg_double_grid_setup_kpoint(self, eph_doublegrid, kpoint, prtvol, comm)
458 
459 !Arguments ------------------------------------
460 !scalars
461  class(ephwg_t),target,intent(inout) :: self
462  type(eph_double_grid_t),intent(inout) :: eph_doublegrid
463  integer,intent(in) :: prtvol, comm
464 !arrays
465  real(dp),intent(in) :: kpoint(3)
466 
467 !Local variables-------------------------------
468 !scalars
469  integer,parameter :: timrev0 = 0
470  integer :: ierr,ii,ik_idx
471  character(len=80) :: errorstring
472  !character(len=500) :: msg
473  type(crystal_t),pointer :: cryst
474 !arrays
475  integer,allocatable :: lgkibz2bz(:) !indkk(:,:),
476  integer,allocatable :: bz2lgkibz(:), bz2lgkibzkq(:) !, bz2bz(:), mapping(:,:)
477  !real(dp) :: kpt(3), wrap_kpt(3), shift
478 !----------------------------------------------------------------------
479 
480  cryst => self%cryst
481 
482  ! Get little group of the (external) kpoint.
483  call self%lgk%free()
484  self%lgk = lgroup_new(self%cryst, kpoint, self%timrev, self%nbz, self%bz, self%nibz, self%ibz, comm)
485  if (prtvol > 0) call self%lgk%print()
486  self%nq_k = self%lgk%nibz
487 
488  ! get dg%bz --> self%lgrp%ibz
489  ABI_REMALLOC(eph_doublegrid%bz2lgkibz, (eph_doublegrid%dense_nbz))
490 
491  ! Old version using all crystal symmetries
492  !call eph_doublegrid%bz2ibz(self%lgk%ibz, self%lgk%nibz,&
493  !                            cryst%symrel, cryst%nsym, &
494  !                            eph_doublegrid%bz2lgkibz, has_timrev=1)
495  call eph_doublegrid%bz2ibz(self%lgk%ibz, self%lgk%nibz,&
496                             self%lgk%symrec_lg, self%lgk%nsym_lg, &
497                             eph_doublegrid%bz2lgkibz, timrev0, use_symrec=.true.)
498 
499  ! self%lgrp%ibz --> dg%bz
500  ABI_ICALLOC(lgkibz2bz, (self%lgk%nibz))
501  do ii=1,self%nbz
502    ik_idx = eph_doublegrid%bz2lgkibz(ii)
503    lgkibz2bz(ik_idx) = ii
504  enddo
505 
506  ! get self%lgrp%ibz --> dg%bz --> self%ibz
507  ABI_REMALLOC(self%lgk2ibz, (self%nq_k))
508  do ii=1,self%nq_k
509    ik_idx = lgkibz2bz(ii)
510    self%lgk2ibz(ii) = eph_doublegrid%bz2ibz_dense(ik_idx)
511  enddo
512 
513  ! calculate k+q
514  do ii=1,self%nq_k
515    self%lgk%ibz(:, ii) = self%lgk%ibz(:, ii) + kpoint
516  end do
517 
518  ! get dg%bz --> lgk%ibz (k+q)
519  ABI_MALLOC(bz2lgkibzkq, (eph_doublegrid%dense_nbz))
520  ! Old version using all crystal symmetries
521  !call eph_doublegrid%bz2ibz(self%lgk%ibz, self%lgk%nibz,&
522  !                            cryst%symrel, cryst%nsym, &
523  !                            bz2lgkibzkq, has_timrev=1)
524  call eph_doublegrid%bz2ibz(self%lgk%ibz, self%lgk%nibz,&
525                             self%lgk%symrec_lg, self%lgk%nsym_lg, &
526                             bz2lgkibzkq, timrev0, use_symrec=.true.)
527 
528  ! self%lgrp%ibz (k+q) --> dg%bz
529  do ii=1,self%nbz
530    ik_idx = bz2lgkibzkq(ii)
531    lgkibz2bz(ik_idx) = ii
532  enddo
533  ABI_FREE(bz2lgkibzkq)
534 
535  ! get self%lgrp%ibz (k+q) --> dg%bz --> self%ibz
536  ABI_REMALLOC(self%kq2ibz, (self%nq_k))
537  do ii=1,self%nq_k
538    ik_idx = lgkibz2bz(ii)
539    self%kq2ibz(ii) = eph_doublegrid%bz2ibz_dense(ik_idx)
540  enddo
541  ABI_FREE(lgkibz2bz)
542 
543  ! revert change
544  do ii=1,self%nq_k
545    self%lgk%ibz(:, ii) = self%lgk%ibz(:, ii) - kpoint
546  end do
547 
548  ! get self%bz --> dg%bz --> self%lgrp%ibz
549  ABI_MALLOC(bz2lgkibz, (self%nbz))
550 
551  do ii=1,self%nbz
552     ! get self%bz --> dg%bz
553     ik_idx = eph_doublegrid%get_index(self%bz(:,ii),2)
554     ! dg%bz --> self%lgrp%ibz
555     bz2lgkibz(ii) = eph_doublegrid%bz2lgkibz(ik_idx)
556  end do
557 
558  ! Build tetrahedron object using IBZ(k) as the effective IBZ
559  ! This means that input data for tetra routines must be provided in lgk%kibz_q
560  call self%tetra_k%free()
561  call htetra_init(self%tetra_k, bz2lgkibz, cryst%gprimd, self%klatt, self%bz, self%nbz, &
562                   self%lgk%ibz, self%nq_k, ierr, errorstring, comm)
563  if (ierr /= 0) then
564    ABI_ERROR(errorstring)
565  end if
566  ABI_FREE(bz2lgkibz)
567 
568 end subroutine ephwg_double_grid_setup_kpoint

m_ephwg/ephwg_free [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_free

FUNCTION

  Deallocate memory

INPUTS

OUTPUT

SOURCE

975 subroutine ephwg_free(self)
976 
977 !Arguments ------------------------------------
978  class(ephwg_t),intent(inout) :: self
979 
980 !----------------------------------------------------------------------
981 
982  ! integer
983  ABI_SFREE(self%kq2ibz)
984 
985  ! Real
986  ABI_SFREE(self%ibz)
987  ABI_SFREE(self%bz)
988  ABI_SFREE(self%lgk2ibz)
989  ABI_SFREE(self%phfrq_ibz)
990  ABI_SFREE(self%eigkbs_ibz)
991 
992  ! types
993  call self%tetra_k%free()
994  call self%lgk%free()
995 
996  ! nullify pointers
997  self%cryst => null()
998 
999 end subroutine ephwg_free

m_ephwg/ephwg_from_ebands [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_from_ebands

FUNCTION

  Convenience constructor to initialize the object from an ebands_t object

m_ephwg/ephwg_get_deltas [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_get_deltas

FUNCTION

 Compute weights for $ \delta(\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
 for a given (band, spin) and phonon mode nu.

INPUTS

 band=band index (global index i.e. unshifted)
 spin=Spin index
 nu=Phonon branch.
 nene=number of energies for DOS
 eminmax=min and  energy in delta (linear mesh)
 bcorr=1 to include Blochl correction else 0.
 comm=MPI communicator
 [broad]=Gaussian broadening

OUTPUT

  deltaw_pm(nene, nq_k, 2)  (plus, minus) including the weights for BZ integration.

SOURCE

637 subroutine ephwg_get_deltas(self, band, spin, nu, nene, eminmax, bcorr, deltaw_pm, comm, &
638                             broad)  ! optional
639 
640 !Arguments ------------------------------------
641 !scalars
642  integer,intent(in) :: band, spin, nu, nene, bcorr, comm
643  class(ephwg_t),intent(in) :: self
644  real(dp),optional,intent(in) :: broad
645 !arrays
646  real(dp),intent(in) :: eminmax(2)
647  real(dp),intent(out) :: deltaw_pm(nene, self%nq_k, 2)
648 
649 !Local variables-------------------------------
650 !scalars
651  integer :: iq,iq_ibz,ikpq_ibz,ib,ie
652  real(dp),parameter :: max_occ1 = one
653  real(dp) :: omega_step
654 !arrays
655  real(dp) :: wme0(nene)
656  real(dp),allocatable :: thetaw(:,:), pme_k(:,:)
657 
658 !----------------------------------------------------------------------
659 
660  ib = band - self%bstart + 1
661 
662  ABI_MALLOC(thetaw, (nene, self%nq_k))
663  ABI_MALLOC(pme_k, (self%nq_k, 2))
664 
665  ! Fill array for e_{k+q, b} +- w_{q,nu)
666  do iq=1,self%nq_k
667    iq_ibz = self%lgk2ibz(iq)   ! IBZ_k --> IBZ
668    ikpq_ibz = self%kq2ibz(iq)  ! k + q --> IBZ
669    pme_k(iq, 1) = self%eigkbs_ibz(ikpq_ibz, ib, spin) - self%phfrq_ibz(iq_ibz, nu)
670    pme_k(iq, 2) = self%eigkbs_ibz(ikpq_ibz, ib, spin) + self%phfrq_ibz(iq_ibz, nu)
671  end do
672 
673  if (present(broad)) then
674    omega_step = (eminmax(2) - eminmax(1)) / (nene - 1)
675    ! Use thetaw as workspace array
676    thetaw(:, 1) = arth(eminmax(1), omega_step, nene)
677    do iq=1,self%nq_k
678      do ie=1,2
679        wme0 = thetaw(:, 1) - pme_k(iq, ie)
680        deltaw_pm(:, iq, ie) = gaussian(wme0, broad)
681      end do
682    end do
683 
684    ! Multiply by weights
685    do ie=1,nene
686      deltaw_pm(ie, :, 1) = deltaw_pm(ie, :, 1) * self%lgk%weights
687      deltaw_pm(ie, :, 2) = deltaw_pm(ie, :, 2) * self%lgk%weights
688    end do
689 
690  else
691    ! TODO Add routine to compute only delta
692    call self%tetra_k%blochl_weights(pme_k(:,1), eminmax(1), eminmax(2), max_occ1, nene, self%nq_k, &
693      bcorr, thetaw, deltaw_pm(:,:,1), comm)
694    call self%tetra_k%blochl_weights(pme_k(:,2), eminmax(1), eminmax(2), max_occ1, nene, self%nq_k, &
695      bcorr, thetaw, deltaw_pm(:,:,2), comm)
696  end if
697 
698  ABI_FREE(thetaw)
699  ABI_FREE(pme_k)
700 
701 end subroutine ephwg_get_deltas

m_ephwg/ephwg_get_deltas_qibzk [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_get_deltas_qibzk

FUNCTION

 Compute weights for $ \delta(\omega - \omega_{q\nu} $ for given nu, using q-opints in the IBZ(k)

INPUTS

 nu=Phonon branch.
 nene=number of energies for DOS
 eminmax=min and  energy in delta (linear mesh)
 bcorr=1 to include Blochl correction else 0.
 comm=MPI communicator
 [with_qweights]= .False. if q-point weights should not be included in dt_weights

OUTPUT

  dt_weights(nene, nq_k, 2)  weights for BZ integration (delta and theta function)

SOURCE

808 subroutine ephwg_get_deltas_qibzk(self, nu, nene, eminmax, bcorr, dt_weights, comm, with_qweights)
809 
810 !Arguments ------------------------------------
811 !scalars
812  integer,intent(in) :: nu, nene, bcorr, comm
813  class(ephwg_t),intent(in) :: self
814  logical,optional,intent(in) :: with_qweights
815 !arrays
816  real(dp),intent(in) :: eminmax(2)
817  real(dp),intent(out) :: dt_weights(nene, self%nq_k, 2)
818 
819 !Local variables-------------------------------
820 !scalars
821  integer :: iq, iq_ibz, ie, ii
822  real(dp),parameter :: max_occ1 = one
823 !arrays
824  real(dp),allocatable :: eigen_in(:)
825 
826 !----------------------------------------------------------------------
827 
828  ABI_MALLOC(eigen_in, (self%nq_k))
829 
830  ! Fill eigen_in
831  do iq=1,self%nq_k
832    iq_ibz = self%lgk2ibz(iq)  ! IBZ_k --> IBZ
833    eigen_in(iq) = self%phfrq_ibz(iq_ibz, nu)
834  end do
835 
836  call self%tetra_k%blochl_weights(eigen_in, eminmax(1), eminmax(2), max_occ1, nene, self%nq_k, &
837    bcorr, dt_weights(:,:,2), dt_weights(:,:,1), comm)
838 
839  if (present(with_qweights)) then
840   if (.not. with_qweights) then
841     do ii=1,2
842       do ie=1,nene
843         dt_weights(ie, :, ii) = dt_weights(ie, :, ii) * self%lgk%weights
844       end do
845     end do
846    end if
847  end if
848 
849  ABI_FREE(eigen_in)
850 
851 end subroutine ephwg_get_deltas_qibzk

m_ephwg/ephwg_get_deltas_wvals [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_get_deltas_wvals

FUNCTION

 Compute weights for $ \delta(\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
 for a given (band, spin) and phonon mode nu.

INPUTS

 band=band index (global index i.e. unshifted)
 spin=Spin index
 nu=Phonon branch.
 nene=number of energies for DOS
 eminmax=min and  energy in delta (linear mesh)
 bcorr=1 to include Blochl correction else 0.
 comm=MPI communicator
 [broad]=Gaussian broadening

OUTPUT

  deltaw_pm(nene, nq_k, 2)  (plus, minus) including the weights for BZ integration.

SOURCE

730 subroutine ephwg_get_deltas_wvals(self, band, spin, nu, neig, eig, bcorr, deltaw_pm, comm, &
731                                   broad)  ! optional
732 
733 !Arguments ------------------------------------
734 !scalars
735  integer,intent(in) :: band, spin, nu, neig, bcorr, comm
736  real(dp),intent(in) :: eig(neig)
737  class(ephwg_t),intent(in) :: self
738  real(dp),optional,intent(in) :: broad
739 !arrays
740  real(dp),intent(out) :: deltaw_pm(neig,self%nq_k, 2)
741 
742 !Local variables-------------------------------
743 !scalars
744  real(dp),parameter :: max_occ1 = one
745  integer :: iq, iq_ibz, ikpq_ibz, ib, nprocs, my_rank
746  real(dp) :: wme0(neig)
747 !arrays
748  real(dp),allocatable :: pme_k(:,:)
749 
750 !----------------------------------------------------------------------
751 
752  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
753  ib = band - self%bstart + 1
754  deltaw_pm = zero
755 
756  ABI_MALLOC(pme_k, (self%nq_k, 2))
757 
758  ! Fill array for e_{k+q, b} +- w_{q,nu)
759  do iq=1,self%nq_k
760    iq_ibz = self%lgk2ibz(iq)   ! IBZ_k --> IBZ
761    ikpq_ibz = self%kq2ibz(iq)  ! k + q --> IBZ
762    pme_k(iq, 1) = self%eigkbs_ibz(ikpq_ibz, ib, spin) - self%phfrq_ibz(iq_ibz, nu)
763    pme_k(iq, 2) = self%eigkbs_ibz(ikpq_ibz, ib, spin) + self%phfrq_ibz(iq_ibz, nu)
764  end do
765 
766  ! Compute the tetrahedron or gaussian weights
767  if (present(broad)) then
768    do iq_ibz=1,self%nq_k
769      if (mod(iq_ibz, nprocs) /= my_rank) cycle ! MPI parallelism
770      wme0 = eig - pme_k(iq_ibz, 1)
771      deltaw_pm(:,iq_ibz,1) = gaussian(wme0, broad) * self%lgk%weights(iq_ibz)
772      wme0 = eig - pme_k(iq_ibz, 2)
773      deltaw_pm(:,iq_ibz,2) = gaussian(wme0, broad) * self%lgk%weights(iq_ibz)
774    end do
775  else
776    call self%tetra_k%wvals_weights_delta(pme_k(:, 1), neig, eig, max_occ1, self%nq_k, bcorr, deltaw_pm(:,:,1), comm)
777    call self%tetra_k%wvals_weights_delta(pme_k(:, 2), neig, eig, max_occ1, self%nq_k, bcorr, deltaw_pm(:,:,2), comm)
778  end if
779 
780  ABI_FREE(pme_k)
781 
782 end subroutine ephwg_get_deltas_wvals

m_ephwg/ephwg_get_zinv_weights [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_get_zinv_weights

FUNCTION

 Compute weights for a given (kpoint, qpoint, spin) for all phonon modes.

INPUTS

 nz: Number of frequencies
 nbcalc=Number of bands in self-energy matrix elements.
 zvals(nw): z-values
 iband_sum = band index in self-energy sum. (global index i.e. unshifted)
 spin=Spin index
 nu=Phonon branch index
 zinv_opt:
   1 for S. Kaprzyk routines,
   2 for Lambin-Vigneron.
 comm=MPI communicator
 [use_bzsum]= By default the weights are multiplied by the Nstar(q) / Nq where
   Nstar(q) is the number of points in the star of the q-point (using the symmetries of the little group of k)
   If use_bzsum is set to True, the Nstar(q) coefficient is removed so that the caller can
   integrate over the BZ without using symmetries.
  [erange(2)]: if present, weights are computed with an approximated asyntotic expression if
   real(z) is outside of this interval and with tetra if inside.

OUTPUT

  cweights(nz, 2, nbcalc, %nq_k)  (plus, minus)
  include weights for BZ integration.

SOURCE

887 subroutine ephwg_get_zinv_weights(self, nz, nbcalc, zvals, iband_sum, spin, nu, zinv_opt, cweights, comm, use_bzsum, erange)
888 
889 !Arguments ------------------------------------
890 !scalars
891  integer,intent(in) :: iband_sum, spin, nu, nz, nbcalc, zinv_opt, comm
892  class(ephwg_t),intent(in) :: self
893  logical, optional, intent(in) :: use_bzsum
894 !arrays
895  complex(dpc),intent(in) :: zvals(nz, nbcalc)
896  complex(dpc),intent(out) :: cweights(nz, 2, nbcalc, self%nq_k)
897  real(dp),optional,intent(in) :: erange(2)
898 
899 !Local variables-------------------------------
900 !scalars
901  integer,parameter :: master = 0
902  integer :: iq_ibz, ikpq_ibz, ib, ii, iq, nprocs, my_rank !, ierr
903  real(dp),parameter :: max_occ1 = one
904  !real(dp) :: emin, emax
905  logical :: use_bzsum_
906 !arrays
907  real(dp) :: my_erange(2)
908  real(dp),allocatable :: pme_k(:,:)
909  complex(dp),allocatable :: cweights_tmp(:,:)
910 !----------------------------------------------------------------------
911 
912  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
913 
914  use_bzsum_ = .False.; if (present(use_bzsum)) use_bzsum_ = use_bzsum
915  my_erange = [-huge(one), huge(one)]; if (present(erange)) my_erange = erange
916 
917  ! Allocate array for e_{k+q, b} +- w_{q,nu)
918  ABI_MALLOC(pme_k, (self%nq_k, 2))
919 
920  ib = iband_sum - self%bstart + 1
921  do iq=1,self%nq_k
922    iq_ibz = self%lgk2ibz(iq)   ! IBZ_k --> IBZ
923    ikpq_ibz = self%kq2ibz(iq)  ! k + q --> IBZ
924    pme_k(iq, 1) = self%eigkbs_ibz(ikpq_ibz, ib, spin) - self%phfrq_ibz(iq_ibz, nu)
925    pme_k(iq, 2) = self%eigkbs_ibz(ikpq_ibz, ib, spin) + self%phfrq_ibz(iq_ibz, nu)
926  end do
927 
928  ! As this part is quite demanding, especially when nz is large, use input z-mesh
929  ! when we are inside the window in which the denominator can blow up (+- some tolerance)
930  ! Outside the window, downsample the mesh use to compute the weights and spline the results.
931  !if (zinv_opt == 2) then
932  !  emin = minval(self%eigkbs_ibz(:, ib, spin))
933  !  emax = maxval(self%eigkbs_ibz(:, ib, spin))
934  !  !my_erange = [emin - half * abs(emin), emax + half * abs(emax)]
935  !  my_erange = [emin - tol2 * abs(emin), emax + tol2 * abs(emax)]
936  !end if
937 
938  cweights = zero
939  ABI_MALLOC(cweights_tmp, (nz, self%nq_k))
940 
941  do ib=1,nbcalc
942    do ii=1,2
943      call self%tetra_k%weights_wvals_zinv(pme_k(:, ii), nz, zvals(:, ib), max_occ1, self%nq_k, zinv_opt, &
944                                           cweights_tmp, comm, erange=my_erange)
945      do iq=1,self%nq_k
946        cweights(:, ii, ib, iq) = cweights_tmp(:, iq)
947      end do
948    end do
949  end do
950 
951  ABI_FREE(cweights_tmp)
952 
953  ! Rescale weights so that the caller can sum over the full BZ.
954  !if (use_bzsum_) cweights = cweights / ( self%lgk%weights(iqlk) * self%nbz )
955 
956  ABI_FREE(pme_k)
957  !call xmpi_sum(cweights, comm, ierr)
958 
959 end subroutine ephwg_get_zinv_weights

m_ephwg/ephwg_new [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_new

FUNCTION

  Initialize the object from the electronic eigenvalues given in the IBZ.

INPUTS

  cryst<cryst_t>=Crystalline structure.
  ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
  bstart=Index of the first band to be included.
  nbcount=Number of bands included
  kptopt=Option for the k-point generation.
  kptrlatt(3,3)=k-point lattice specification
  nshiftk= number of shift vectors.
  shiftk(3,nshiftk)=shift vectors for k point generation
  nkibz=Number of points in the IBZ
  kibz(3,nkibz)=Reduced coordinates of the k-points in the IBZ.
  nsppol=Number of independent spin polarizations.
  eig_ibz(nbcount, nkibz, nsppol) = Electron eigenvalues for nbcount states

OUTPUT

SOURCE

205 type(ephwg_t) function ephwg_new( &
206    cryst, ifc, bstart, nbcount, kptopt, kptrlatt, nshiftk, shiftk, nkibz, kibz, nsppol, eig_ibz, comm) result(new)
207 
208 !Arguments ------------------------------------
209 !scalars
210  integer,intent(in) :: kptopt, nshiftk, nkibz, bstart, nbcount, nsppol, comm
211  type(crystal_t),target,intent(in) :: cryst
212  type(ifc_type),intent(in) :: ifc
213 !arrays
214  integer,intent(in) :: kptrlatt(3,3)
215  real(dp),intent(in) :: shiftk(3, nshiftk), kibz(3, nkibz)
216  real(dp),intent(in) :: eig_ibz(nbcount, nkibz, nsppol)
217 
218 !Local variables-------------------------------
219 !scalars
220  integer :: nprocs, my_rank, ik, ierr, out_nkibz
221  real(dp) :: cpu, wall, gflops
222 !arrays
223  real(dp) :: rlatt(3,3)
224  integer :: out_kptrlatt(3,3)
225  real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom), phfrq(3*cryst%natom)
226  real(dp),allocatable :: out_kibz(:,:), out_wtk(:)
227 
228 !----------------------------------------------------------------------
229 
230  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
231 
232  new%natom3 = ifc%natom * 3
233  new%nsppol = nsppol
234  new%nbcount = nbcount
235  new%bstart = bstart
236  new%kptopt = kptopt
237  new%timrev = kpts_timrev_from_kptopt(new%kptopt)
238  new%nibz = nkibz
239  new%cryst => cryst
240  call alloc_copy(kibz, new%ibz)
241 
242  call cwtime(cpu, wall, gflops, "start")
243 
244  ! Get full BZ (new%nbz, new%bz) and new kptrlatt for tetra.
245  call kpts_ibz_from_kptrlatt(cryst, kptrlatt, kptopt, nshiftk, shiftk, out_nkibz, out_kibz, out_wtk, new%nbz, new%bz, &
246                              new_kptrlatt=out_kptrlatt)
247  call cwtime_report(" ephwg_new: kpts_ibz_from_kptrlatt", cpu, wall, gflops)
248 
249  new%kptrlatt = out_kptrlatt
250  rlatt = out_kptrlatt; call matr3inv(rlatt, new%klatt)
251 
252  ABI_CHECK(size(out_kibz, dim=2) == new%nibz, "mismatch in nkibz!")
253  ABI_FREE(out_kibz)
254  ABI_FREE(out_wtk)
255 
256  ! Copy eigenvalues in IBZ. Change shape for better performance in other routines.
257  ABI_MALLOC(new%eigkbs_ibz, (new%nibz, new%nbcount, new%nsppol))
258  do ik=1,new%nibz
259    new%eigkbs_ibz(ik, :, :) = eig_ibz(:, ik, :)
260  end do
261 
262  ! Fourier interpolate phonon frequencies on the same mesh.
263  ABI_CALLOC(new%phfrq_ibz, (new%nibz, new%natom3))
264 
265  do ik=1,new%nibz
266    if (mod(ik, nprocs) /= my_rank) cycle ! mpi-parallelism
267    call ifc%fourq(cryst, new%ibz(:, ik), phfrq, displ_cart)
268    new%phfrq_ibz(ik, :) = phfrq
269  end do
270 
271  ! Collect results on each rank
272  call xmpi_sum(new%phfrq_ibz, comm, ierr)
273 
274  !new%max_phfrq = maxval(phfrq_ibz)
275 
276  call cwtime_report(" ephwg_new: ifc_fourq", cpu, wall, gflops)
277 
278 end function ephwg_new

m_ephwg/ephwg_report_stats [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_report_stats

FUNCTION

INPUTS

OUTPUT

SOURCE

585 subroutine ephwg_report_stats(self)
586 
587 !Arguments ------------------------------------
588 !scalars
589  class(ephwg_t),intent(in) :: self
590 
591 !Variables
592  real(dp) :: mem_tot
593 !----------------------------------------------------------------------
594 
595  ! IBZ qpoints
596  mem_tot = 3 * self%nibz * dp
597  ! BZ qpoints
598  mem_tot = mem_tot + 3 * self%nbz * dp
599  ! lgk2ibz and kq2ibz
600  mem_tot = mem_tot + self%nq_k * 2 * 4
601  ! phonon frequencies
602  mem_tot = mem_tot + self%nibz * self%natom3 * dp
603  ! eigenvalues
604  mem_tot = mem_tot + self%nibz * self%nbcount * self%nsppol * dp
605 
606  write(std_out,"(a,f8.1,a)") " Memory allocated for ephwg weights:", mem_tot * b2Mb, " [Mb] <<< MEM"
607 
608 end subroutine ephwg_report_stats

m_ephwg/ephwg_setup_kpoint [ Functions ]

[ Top ] [ m_ephwg ] [ Functions ]

NAME

 ephwg_setup_kpoint

FUNCTION

  Set internal tables and object required to compute integration weights for a given k-point.

INPUTS

  kpoint(3): k-point in reduced coordinates.
  prtvol: Verbosity level
  comm: MPI communicator

OUTPUT

SOURCE

337 subroutine ephwg_setup_kpoint(self, kpoint, prtvol, comm, skip_mapping)
338 
339 !Arguments ------------------------------------
340 !scalars
341  class(ephwg_t),target,intent(inout) :: self
342  integer,intent(in) :: prtvol, comm
343  logical,optional,intent(in) :: skip_mapping
344 !arrays
345  real(dp),intent(in) :: kpoint(3)
346 
347 !Local variables-------------------------------
348 !scalars
349  integer :: ierr,ii
350  logical :: do_mapping
351  real(dp) :: cpu, wall, gflops
352  character(len=80) :: errorstring
353  !character(len=500) :: msg
354  type(crystal_t),pointer :: cryst
355  type(krank_t) :: krank
356 !arrays
357  integer,allocatable :: indkk(:,:)
358 
359 !----------------------------------------------------------------------
360 
361  do_mapping = .true.; if (present(skip_mapping)) do_mapping = .not. skip_mapping
362  cryst => self%cryst
363  call cwtime(cpu, wall, gflops, "start")
364 
365  ! Get little group of the (external) kpoint.
366  call self%lgk%free()
367  self%lgk = lgroup_new(self%cryst, kpoint, self%timrev, self%nbz, self%bz, self%nibz, self%ibz, comm)
368 
369  if (prtvol > 0) call self%lgk%print()
370  self%nq_k = self%lgk%nibz
371 
372  call cwtime_report(" lgroup_new", cpu, wall, gflops)
373 
374  if (do_mapping) then
375    ! TODO: Use symrec conventions although this means that we cannot reuse these tables
376    ! to symmetrize wavefunctions and potentials that require S-1 i.e. the symrel convention.
377 
378    ! Get mapping IBZ_k --> initial IBZ (self%lgk%ibz --> self%ibz)
379    ABI_MALLOC(indkk, (6, self%nq_k))
380 
381    krank = krank_from_kptrlatt(self%nibz, self%ibz, self%kptrlatt, compute_invrank=.False.)
382 
383    if (kpts_map("symrel", self%timrev, cryst, krank, self%nq_k, self%lgk%ibz, indkk) /= 0) then
384      ABI_ERROR("At least one of the points in IBZ(k) could not be generated from a symmetrical one.")
385    end if
386 
387    call krank%free()
388 
389    ABI_SFREE(self%lgk2ibz)
390    call alloc_copy(indkk(1, :), self%lgk2ibz)
391    ABI_FREE(indkk)
392    call cwtime_report(" listkk1", cpu, wall, gflops)
393 
394    ! Get mapping (k + q) --> initial IBZ.
395    do ii=1,self%nq_k
396      self%lgk%ibz(:, ii) = self%lgk%ibz(:, ii) + kpoint
397    end do
398    ABI_MALLOC(indkk, (6, self%nq_k))
399 
400    krank = krank_from_kptrlatt(self%nibz, self%ibz, self%kptrlatt, compute_invrank=.False.)
401 
402    if (kpts_map("symrel", self%timrev, cryst, krank, self%nq_k, self%lgk%ibz, indkk) /= 0) then
403      ABI_ERROR("At least one of the points in IBZ(k) + q could not be generated from a symmetrical one.")
404    end if
405    call krank%free()
406 
407    call cwtime_report(" listkk2", cpu, wall, gflops)
408 
409    ABI_SFREE(self%kq2ibz)
410    call alloc_copy(indkk(1, :), self%kq2ibz)
411    ABI_FREE(indkk)
412 
413    ! Revert changes
414    do ii=1,self%nq_k
415      self%lgk%ibz(:, ii) = self%lgk%ibz(:, ii) - kpoint
416    end do
417  end if
418 
419  ! Get mapping BZ --> IBZ_k (self%bz --> self%lgrp%ibz) required for tetrahedron method
420  ABI_MALLOC(indkk, (self%nbz, 1))
421  indkk(:, 1) = self%lgk%bz2ibz_smap(1, :)
422 
423  ! Build tetrahedron object using IBZ(k) as the effective IBZ
424  ! This means that input data for tetra routines must be provided in lgk%kibz_q
425  call self%tetra_k%free()
426  call htetra_init(self%tetra_k, indkk(:, 1), cryst%gprimd, self%klatt, self%bz, self%nbz, &
427                   self%lgk%ibz, self%nq_k, ierr, errorstring, comm)
428  !call tetra_write(self%tetra_k, self%lgk%nibz, self%lgk%ibz, strcat("tetrak_", ktoa(kpoint)))
429  ABI_CHECK(ierr == 0, errorstring)
430 
431  if (xmpi_comm_rank(comm) == 0) call self%tetra_k%print(std_out)
432  ABI_FREE(indkk)
433 
434  call cwtime_report(" init_tetra", cpu, wall, gflops)
435 
436 end subroutine ephwg_setup_kpoint

m_ephwg/ephwg_t [ Types ]

[ Top ] [ m_ephwg ] [ Types ]

NAME

 ephwg_t

FUNCTION

  Stores electron eigevalues and phonon frequencies in the IBZ (assume same mesh for e and ph).
  Provides tools to compute (e_{k+q} - w{q}) in the IBZ(k)
  and integrate the delta functions for phonon emission/absorption with the tetrahedron method.

SOURCE

 72 type, public :: ephwg_t
 73 
 74   integer :: natom3
 75   ! 3 * natom
 76 
 77   integer :: nsppol
 78   ! Number of independent spin polarizations.
 79 
 80   integer :: nbcount
 81   ! Number of bands treated.
 82 
 83   integer :: bstart
 84   ! The fist band (global index) starts at bstart.
 85   ! Used to select bands around the Fermi level.
 86 
 87   integer :: kptopt
 88   ! Option for k-point generation.
 89 
 90   integer :: timrev
 91   ! 1 if the use of time-reversal is allowed; 0 otherwise
 92 
 93   integer :: nibz, nbz
 94   ! Number of q-points in IBZ and full BZ.
 95 
 96   integer :: nq_k
 97   ! Number of points in IBZ(k) i.e. the irreducible wedge
 98   ! defined by the operations of the little group of k.
 99 
100   !real(dp) :: max_phfrq
101   ! Max Phonon frequency, computed from phfrq_ibz
102 
103   integer :: kptrlatt(3,3)
104    ! Value of kptrlatt after inkpts. So one shift
105 
106   integer, allocatable :: kq2ibz(:)
107   ! kq2ibz(nq_k)
108   ! Mapping (k + q) --> initial IBZ array
109 
110   real(dp),allocatable :: ibz(:,:)
111   ! ibz(3, nibz)
112   ! The initial IBZ.
113 
114   real(dp),allocatable :: bz(:,:)
115   ! bz(3, nbz)
116   ! points in full BZ.
117 
118   real(dp) :: klatt(3, 3)
119   ! Reciprocal of lattice vectors for full kpoint grid. Used by init_tetra
120 
121   integer,allocatable :: lgk2ibz(:)
122   ! lgk2ibz(nq_k)
123   ! Mapping Little-group IBZ_k --> initial IBZ
124   ! TODO: This should be generalized to have the symmetry indices as well so
125   ! that we can use it in sigmaph but then we have to implement similar algo for double grid.
126 
127   real(dp),allocatable :: phfrq_ibz(:,:)
128   ! (nibz, natom3)
129   ! Phonon frequencies in the IBZ
130 
131   real(dp),allocatable :: eigkbs_ibz(:, :, :)
132   ! (nibz, nbcount, nsppol)
133   ! Electron eigenvalues in the IBZ for nbcount states
134   ! (not necessarly equal to global nband, see also bstart and bcount)
135 
136   type(crystal_t), pointer :: cryst => null()
137   ! Pointer to input structure (does not own memory)
138 
139   type(lgroup_t) :: lgk
140   ! Little group of the k-point
141 
142   type(htetra_t) :: tetra_k
143   ! Used to evaluate delta(w - e_{k+q} +/- phw_q) with tetrahedron method.
144 
145  contains
146 
147      procedure :: setup_kpoint => ephwg_setup_kpoint
148      ! Prepare tetrahedron method for given external k-point.
149 
150      procedure :: double_grid_setup_kpoint => ephwg_double_grid_setup_kpoint
151      ! Prepare tetrahedron method for given external k-point using double grid routines.
152 
153      procedure :: report_stats => ephwg_report_stats
154      ! Report how much memory is being used by this object
155 
156      procedure :: get_deltas => ephwg_get_deltas
157      ! Compute weights for $ \int \delta(\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
158 
159      procedure :: get_deltas_wvals => ephwg_get_deltas_wvals
160      ! Compute weights for $ \int \delta(\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
161 
162      procedure :: get_deltas_qibzk => ephwg_get_deltas_qibzk
163      ! Compute weights for $ \delta(\omega - \omega_{q\nu}} $ using IBZ(k) ordering
164 
165      procedure :: get_zinv_weights => ephwg_get_zinv_weights
166      ! Compute weights for $ \int 1 / (\omega - \ee_{k+q, b} \pm \omega_{q\nu} $
167 
168      procedure :: free => ephwg_free
169      ! Free memory
170  end type ephwg_t
171 
172  public :: ephwg_new             ! Basic Constructor
173  public :: ephwg_from_ebands     ! Build object from ebands_t
174 
175 contains