TABLE OF CONTENTS


ABINIT/m_chi0tk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_chi0tk

FUNCTION

  This module provides tools for the computation of the irreducible polarizability.

COPYRIGHT

 Copyright (C) 1999-2018 ABINIT group (MG, FB)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 MODULE m_chi0tk
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32  use m_xmpi
33  use m_xomp
34  use m_sort
35  use m_wfd
36 
37  use defs_datatypes, only : ebands_t
38  use m_gwdefs,   only : GW_TOL_DOCC, czero_gw, cone_gw, em1params_t, j_gw
39  use m_fstrings, only : sjoin, itoa
40  use m_hide_blas,     only : xgerc, xgemm
41  use m_crystal,  only : crystal_t
42  use m_gsphere,  only : gsphere_t, gsph_gmg_idx, gsph_gmg_fftidx
43  use m_bz_mesh,  only : littlegroup_t, kmesh_t, has_BZ_item
44 
45  implicit none
46 
47  private
48 
49  public :: assemblychi0_sym
50  public :: symmetrize_afm_chi0
51  public :: accumulate_chi0_q0
52  public :: accumulate_sfchi0_q0
53  public :: assemblychi0sf
54  public :: approxdelta
55  public :: setup_spectral
56  public :: hilbert_transform
57  public :: hilbert_transform_headwings
58  public :: completechi0_deltapart
59  public :: output_chi0sumrule
60  public :: accumulate_chi0sumrule
61  public :: make_transitions
62  public :: chi0_bbp_mask

m_chi0tk/accumulate_chi0_q0 [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 accumulate_chi0_q0

FUNCTION

 Update the independent particle susceptibility at q==0 for the contribution
 of one pair of occupied-unoccupied band, for each frequency.
 This routine takes advantage of the symmetries of the little group of the external q-point
 to symmetrize the contribution arising from the input k-point located in the IBZ_q.
 It computes:

   $ \chi_0(G1,G2,io) = \chi_0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))*green_w(io) $

 where S is a symmetry in reciprocal space.
 The matrix elements of the gradient operator and [V_{nl},r] are symmetrized as well.

INPUTS

  ik_bz=Index of the k-point whose contribution has to be added to chi0.
  isym_kbz=Index of the symmetry such that k_bz = IS k_ibz
  itim_kbz=2 if time-reversal has to be used to obtain k_bz, 1 otherwise.
  npwepG0=Maximum number of G vectors
  rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states.
  rhotwx(3,nspinor**2)=Matrix element of the operator $-i[H,r]/(e1-e2) = -i r$ in reciprocal lattice units.
  green_w(nomega)=Frequency dependent part of the Green function.
  Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point.
    %timrev=2 it time-reversal is used, 1 otherwise.
    %nsym_sg=Number of space group symmetries.
    %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point.
  Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors
    %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg.
    %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion.
    %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations.
  Cryst<crystal_t>=Structure defining the unit cell and its symmetries
    %nsym=Number of symmetries.
    %symrec(3,3,nsym)=Symmetry operation in reciprocal space (reduced coordinates)
  Ep<em1params_t>=Parameters of the chi0 calculation.
     %npwe=number of plane waves in chi0.
     %symchi=1 if symmetrization has to be performed.
     %nomega=number of frequencies in chi0.

OUTPUT

  (see side effects)

SIDE EFFECTS

  chi0(npwe,npwe,nomega)= Updated independent-particle susceptibility matrix in reciprocal space at q==0.
  chi0_head(3,3,Ep%nomega)=Head.
  chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)=Lower wing.
  chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)=Upper wing.

NOTES

  1) Symmetrization of the oscilator matrix elements.
    If  Sq = q then  M_G( Sk,q)= e^{-i(q+G).\tau} M_{ S^-1G}  (k,q)
    If -Sq = q then  M_G(-Sk,q)= e^{-i(q+G).\tau} M_{-S^-1G}^*(k,q)

    In the case of umklapps:
    If  Sq = q+G0 then  M_G( Sk,q)= e^{-i(q+G).\tau} M_{ S^-1(G-G0}   (k,q)
    If -Sq = q+G0 then  M_G(-Sk,q)= e^{-i(q+G).\tau} M_{-S^-1(G-G0)}^*(k,q)

  In the equation below there is no need to take into account the phases due to q.t
  as they cancel each other in the scalar product ==> only phmGt(G,isym)=e^{-iG.\tau} is needed.

  2) Symmetrization of the matrix elements of the position operator.

    <Sk,b|\vec r| Sk,b'> = <k b| R\vec r + \tau|k b'>

     where S is one of the symrec operation, R and \tau is the corresponding
     operation in real space. The term involving the fractional translation is zero provided that b /= b'.

PARENTS

      cchi0q0

CHILDREN

SOURCE

674 subroutine accumulate_chi0_q0(ik_bz,isym_kbz,itim_kbz,gwcomp,nspinor,npwepG0,Ep,Cryst,Ltg_q,Gsph_epsG0,&
675 & chi0,rhotwx,rhotwg,green_w,green_enhigh_w,deltaf_b1b2,chi0_head,chi0_lwing,chi0_uwing)
676 
677 
678 !This section has been created automatically by the script Abilint (TD).
679 !Do not modify the following lines by hand.
680 #undef ABI_FUNC
681 #define ABI_FUNC 'accumulate_chi0_q0'
682 !End of the abilint section
683 
684  implicit none
685 
686 !Arguments ------------------------------------
687 !scalars
688  integer,intent(in) :: ik_bz,isym_kbz,itim_kbz,npwepG0,nspinor,gwcomp
689  real(dp),intent(in) :: deltaf_b1b2
690  type(littlegroup_t),intent(in) :: Ltg_q
691  type(gsphere_t),target,intent(in) :: Gsph_epsG0
692  type(crystal_t),intent(in) :: Cryst
693  type(em1params_t),intent(in) :: Ep
694 !arrays
695  complex(gwpc),intent(in) :: rhotwg(npwepG0)
696  complex(gwpc),intent(in) :: rhotwx(3,nspinor**2)
697  complex(gwpc),intent(inout) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega)
698  complex(dpc),intent(in) :: green_w(Ep%nomega),green_enhigh_w(Ep%nomega)
699  complex(dpc),intent(inout) :: chi0_head(3,3,Ep%nomega)
700  complex(dpc),intent(inout) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)
701  complex(dpc),intent(inout) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)
702 
703 !Local variables-------------------------------
704 !scalars
705  integer :: itim,io,isym,idir,jdir
706  complex(gwpc) :: dd
707  !character(len=500) :: msg
708 !arrays
709  integer,ABI_CONTIGUOUS pointer :: Sm1G(:)
710  complex(dpc) :: mir_kbz(3)
711  complex(gwpc),allocatable :: rhotwg_sym(:)
712  complex(gwpc), ABI_CONTIGUOUS pointer :: phmGt(:)
713 
714 !************************************************************************
715 
716  ABI_UNUSED(deltaf_b1b2)
717 
718  SELECT CASE (Ep%symchi)
719  CASE (0)
720    ! Do not use symmetries.
721    ! Symmetrize rhotwg in the full BZ and accumulate over the full BZ i.e.
722    !   chi0(G1,G2,io) = chi0(G1,G2,io) + (rhotwg(G1)*CONJG(rhotwg(G2)))*green_w(io)
723    !
724    ! The non-analytic term is symmetrized for this k-point in the BZ according to:
725    !    rhotwg(1)= S^-1q * rhotwx_ibz
726    !    rhotwg(1)=-S^-1q * CONJG(rhotwx_ibz) if time-reversal is used.
727 
728    ! Multiply elements G1,G2 of rhotwg_sym by green_w(io) and accumulate in chi0(G1,G2,io)
729 !$omp parallel do private(dd)
730    do io=1,Ep%nomega
731      dd=green_w(io)
732      call XGERC(Ep%npwe,Ep%npwe,dd,rhotwg,1,rhotwg,1,chi0(:,:,io),Ep%npwe)
733    end do
734 
735    ! === Accumulate heads and wings for each small q ===
736    ! FIXME extrapolar method should be checked!!
737    ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
738    if (nspinor == 1) then
739      mir_kbz = (3-2*itim_kbz) * MATMUL(Cryst%symrec(:,:,isym_kbz),rhotwx(:,1))
740    else
741      mir_kbz = (3-2*itim_kbz) * MATMUL(Cryst%symrec(:,:,isym_kbz), sum(rhotwx(:,1:2), dim=2))
742    end if
743    if (itim_kbz==2) mir_kbz=CONJG(mir_kbz)
744 
745    ! here we might take advantage of Hermiticity along Im axis in RPA (see mkG0w)
746    do idir=1,3
747      do io=1,Ep%nomega
748        chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir)*CONJG(rhotwg)
749        chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg*CONJG(mir_kbz(idir))
750        ! Add contribution due to extrapolar technique.
751        !if (gwcomp==1.and.ABS(deltaf_b1b2) >= GW_TOL_DOCC) then
752        if (gwcomp==1) then
753          chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_enhigh_w(io) * mir_kbz(idir)*CONJG(rhotwg)
754          chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_enhigh_w(io) * rhotwg*CONJG(mir_kbz(idir))
755        end if
756      end do
757    end do
758 
759    ! Accumulate the head.
760    do io=1,Ep%nomega
761      do jdir=1,3
762        do idir=1,3
763          chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_w(io) * mir_kbz(idir)*CONJG(mir_kbz(jdir))
764          ! Add contribution due to extrapolar technique.
765          !if (gwcomp==1.and.ABS(deltaf_b1b2) >= GW_TOL_DOCC) then
766          if (gwcomp==1) then
767            chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_enhigh_w(io) * mir_kbz(idir)*CONJG(mir_kbz(jdir))
768          end if
769        end do
770      end do
771    end do
772 
773  CASE (1)
774    ! Use symmetries to reconstruct the integrand.
775    ABI_MALLOC(rhotwg_sym, (Ep%npwe))
776 
777    ! Loop over symmetries of the space group and time-reversal.
778    do isym=1,Ltg_q%nsym_sg
779      do itim=1,Ltg_q%timrev
780 
781        if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
782          ! This operation belongs to the little group and has to be considered to reconstruct the BZ.
783          phmGt => Gsph_epsG0%phmGt  (1:Ep%npwe,isym) ! In the 2 lines below note the slicing (1:npwe)
784          Sm1G  => Gsph_epsG0%rottbm1(1:Ep%npwe,itim,isym)
785 
786          SELECT CASE (itim)
787          CASE (1)
788            rhotwg_sym(1:Ep%npwe)=rhotwg(Sm1G(1:Ep%npwe))*phmGt(1:Ep%npwe)
789          CASE (2)
790            rhotwg_sym(1:Ep%npwe)=CONJG(rhotwg(Sm1G(1:Ep%npwe)))*phmGt(1:Ep%npwe)
791          CASE DEFAULT
792            MSG_BUG(sjoin('Wrong value of itim:', itoa(itim)))
793          END SELECT
794 
795          ! Multiply elements G1,G2 of rhotwg_sym by green_w(io) and accumulate in chi0(G,Gp,io)
796 !$omp parallel do private(dd)
797          do io=1,Ep%nomega
798            dd=green_w(io)
799            call XGERC(Ep%npwe,Ep%npwe,dd,rhotwg_sym,1,rhotwg_sym,1,chi0(:,:,io),Ep%npwe)
800          end do
801 
802          ! === Accumulate heads and wings for each small q ===
803          ! FIXME extrapolar method should be checked!!
804 
805          ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
806          if (nspinor == 1) then
807            mir_kbz =(3-2*itim) * MATMUL(Cryst%symrec(:,:,isym),rhotwx(:,1))
808          else
809            mir_kbz = (3-2*itim) * MATMUL(Cryst%symrec(:,:,isym), sum(rhotwx(:,1:2), dim=2))
810          end if
811 
812          if (itim==2) mir_kbz=CONJG(mir_kbz)
813 
814          ! here we might take advantage of Hermiticity along Im axis in RPA (see mkG0w)
815          do idir=1,3
816            do io=1,Ep%nomega
817              chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir)*CONJG(rhotwg_sym)
818              chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg_sym*CONJG(mir_kbz(idir))
819              ! Add contribution due to extrapolar technique.
820              !if (gwcomp==1.and.ABS(deltaf_b1b2)>=GW_TOL_DOCC) then
821              if (gwcomp==1) then
822                chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_enhigh_w(io) * mir_kbz(idir)*CONJG(rhotwg_sym)
823                chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_enhigh_w(io) * rhotwg_sym*CONJG(mir_kbz(idir))
824              end if
825            end do
826          end do
827 
828          ! Accumulate the head.
829          do io=1,Ep%nomega
830            do jdir=1,3
831              do idir=1,3
832                 chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) +  green_w(io) * mir_kbz(idir)*CONJG(mir_kbz(jdir))
833                 ! Add contribution due to extrapolar technique.
834                 !if (gwcomp==1.and.ABS(deltaf_b1b2) >= GW_TOL_DOCC) then
835                 if (gwcomp==1) then
836                   chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_enhigh_w(io)*mir_kbz(idir)*CONJG(mir_kbz(jdir))
837                 end if
838              end do
839            end do
840          end do
841 
842        end if !wtksym
843      end do !itim
844    end do !isym
845 
846    ABI_FREE(rhotwg_sym)
847 
848  CASE DEFAULT
849    MSG_BUG(sjoin('Wrong value of symchi ',itoa(Ep%symchi)))
850  END SELECT
851 
852 end subroutine accumulate_chi0_q0

m_chi0tk/accumulate_chi0sumrule [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 accumulate_chi0sumrule

FUNCTION

  Accumulate the contribution to the sum rule for Im chi0
  arising from a single transition. Eventually symmetrize
  it using the symmetry operations of the little group of q.

INPUTS

  ik_bz=Index of the k-point in the full BZ whose contribution has to be added and symmetrized.
  symchi=1 if we are summing over IBZ_q and symmetrization has to be performed.
  npwe=Number of G vectors in chi0.
  npwepG0=Number of G vectors in the "enlarged" sphere to treat umklapp.
  factor=factor entering the expression.
  delta_ene=Transition energy.
  Ltg_q=<littlegroup_t>= Structure gathering information on the little group of the external q.
  Gsph_epsG0=<gsphere_t>=Info on the G-sphere for chi0.
  rhotwg(npwepG0)=Fouriet transform of u_{b1 k-q} u_{b2 k} in the "enlarged" sphere.

OUTPUT

  See SIDES EFFECTS

 SIDES EFFECTS
  chi0sumrule(npwe)= In input the sum rule calculated so far,
  In output the contribution of this transition is accounted for, and, eventually, symmetrized.
  using the symmetry operations of the little group of the external q.

PARENTS

      cchi0,cchi0q0

CHILDREN

SOURCE

2132 subroutine accumulate_chi0sumrule(ik_bz,symchi,npwe,factor,delta_ene,&
2133 & Ltg_q,Gsph_epsG0,npwepG0,rhotwg,chi0sumrule)
2134 
2135 
2136 !This section has been created automatically by the script Abilint (TD).
2137 !Do not modify the following lines by hand.
2138 #undef ABI_FUNC
2139 #define ABI_FUNC 'accumulate_chi0sumrule'
2140 !End of the abilint section
2141 
2142  implicit none
2143 
2144 !Arguments ------------------------------------
2145 !scalars
2146  integer,intent(in) :: ik_bz,npwe,npwepG0,symchi
2147  real(dp),intent(in) :: delta_ene,factor
2148  type(gsphere_t),intent(in) :: Gsph_epsG0
2149  type(littlegroup_t),target,intent(in) :: Ltg_q
2150 !arrays
2151  real(dp),intent(inout) :: chi0sumrule(npwe)
2152  complex(gwpc),intent(in) :: rhotwg(npwepG0)
2153 
2154 !Local variables-------------------------------
2155 !scalars
2156  integer :: isym,itim
2157  !character(len=500) :: msg
2158 !arrays
2159  integer,allocatable :: Sm1_gmG0(:)
2160  integer, ABI_CONTIGUOUS pointer :: gmG0(:)
2161  complex(gwpc),allocatable :: rhotwg_sym(:)
2162 
2163 !************************************************************************
2164 
2165  ! Accumulating the sum rule on chi0.
2166  ! Eq.(5.284) in G. D. Mahan Many-Particle Physics 3rd edition [[cite:Mahan2000]]
2167 
2168  SELECT CASE (symchi)
2169  CASE (0)
2170    ! Do not use symmetries, sum is performed in the full BZ.
2171    chi0sumrule(:)=chi0sumrule(:) + factor*delta_ene*ABS(rhotwg(1:npwe))**2
2172 
2173  CASE (1)
2174    ! Symmetrize the contribution in the full BZ.
2175    ABI_ALLOCATE(rhotwg_sym,(npwe))
2176    ABI_ALLOCATE(Sm1_gmG0,(npwe))
2177 
2178    do itim=1,Ltg_q%timrev
2179      do isym=1,Ltg_q%nsym_sg
2180        if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
2181         ! This operation belongs to the little group and has to be used to reconstruct the BZ ===
2182         ! In the following 2 lines mind the slicing (1:npwe)
2183         gmG0  => Ltg_q%igmG0(1:npwe,itim,isym)
2184         Sm1_gmG0(1:npwe)=Gsph_epsG0%rottbm1(gmG0(1:npwe),itim,isym)
2185         rhotwg_sym(1:npwe)=rhotwg(Sm1_gmG0)
2186 
2187         chi0sumrule(:)=chi0sumrule(:) + factor*delta_ene*ABS(rhotwg_sym(1:npwe))**2
2188        end if
2189      end do !isym
2190    end do !itim
2191 
2192    ABI_DEALLOCATE(rhotwg_sym)
2193    ABI_DEALLOCATE(Sm1_gmG0)
2194 
2195  CASE DEFAULT
2196    MSG_BUG(sjoin('Wrong value for symchi:', itoa(symchi)))
2197  END SELECT
2198 
2199 end subroutine accumulate_chi0sumrule

m_chi0tk/accumulate_sfchi0_q0 [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 accumulate_sfchi0_q0

FUNCTION

 Update the spectral function of the independent particle susceptibility at q==0 for the contribution
 of one pair of occupied-unoccupied band, for each frequency.
 If symchi==1, the symmetries belonging to the little group of the external point q are used
 to reconstruct the contributions in the full Brillouin zone. In this case, the equation implented is:

  $ chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))* \delta(\omega -trans) $

 where S is a symmetry belonging to the little group of q.
 The subroutine also performs the symmetrization of the matrix elements of the
 gradient operator and of the commutator [V_{nl},r] with the position operator.

INPUTS

  ikbz=Index in the BZ of the k-point whose contribution to chi0 has to be added,
   if we use symmetries, the contribution to chi0 by this k-point has to be symmetrized.
  isym_kbz=Index of the symmetry such as k_bz = IS k_ibz
  itim_kbz=2 if time-reversal has to be used to obtain k_bz, 1 otherwise.
  my_wl,my_wr=min and Max frequency index treated by this processor.
  npwe=Number of plane waves used to describe chi0.
  npwepG0=Maximum number of G vectors to account for umklapps.
  nomega=Number of frequencies in the imaginary part.
  rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states.
  rhotwx(3,nspinor**2)=Matrix elements of the gradient and of the commutator of the non-local operator with
    the position operator. The second term is present only if inclvkb=1,2.
  Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors
    %ng=number of G vectors in the enlarged sphere. It MUST be equal to the size of rhotwg
    %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion
    %phmGt(ng,nsym)=phase factors associated to non-symmorphic operations
  Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point.
    %timrev=2 it time-reversal is used, 1 otherwise
    %nsym_sg=Number of space group symmetries
    %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point
    %flag_umklp(timrev,nsym)= flag for umklapp processes
     if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise
 Cryst<crystal_t>=Info on unit cell and it symmetries
    %nsym=Number of symmetry operations.
    %symrec(3,3,nsym)=Symmetry operations in reciprocal space (reduced coordinates).

OUTPUT

  (see side effects)

SIDE EFFECTS

  sf_chi0(npwe,npwe,my_wl:my_wr)=Updated spectral function at q==0.
  sf_lwing(npwe,my_wl:my_wr,3)=Updated lower wing of the spectral function.
  sf_uwing(npwe,mw_wl:my_wr,3)=Updated upper wing of the spectral function.
  sf_head(3,3,my_wl:my_wr)=Updated head of the spectral function.

PARENTS

      cchi0q0

CHILDREN

SOURCE

 915 subroutine accumulate_sfchi0_q0(ikbz,isym_kbz,itim_kbz,nspinor,symchi,npwepG0,npwe,Cryst,Ltg_q,Gsph_epsG0,&
 916 & factocc,my_wl,iomegal,wl,my_wr,iomegar,wr,rhotwx,rhotwg,nomegasf,sf_chi0,sf_head,sf_lwing,sf_uwing)
 917 
 918 
 919 !This section has been created automatically by the script Abilint (TD).
 920 !Do not modify the following lines by hand.
 921 #undef ABI_FUNC
 922 #define ABI_FUNC 'accumulate_sfchi0_q0'
 923 !End of the abilint section
 924 
 925  implicit none
 926 
 927 !Arguments ------------------------------------
 928 !scalars
 929  integer,intent(in) :: ikbz,my_wl,my_wr,nomegasf,npwe,npwepG0,nspinor
 930  integer,intent(in) :: isym_kbz,itim_kbz,symchi,iomegal,iomegar
 931  real(dp),intent(in) :: factocc,wl,wr
 932  type(littlegroup_t),intent(in) :: Ltg_q
 933  type(gsphere_t),target,intent(in) :: Gsph_epsG0
 934  type(crystal_t),intent(in) :: Cryst
 935 !arrays
 936  complex(gwpc),intent(in) :: rhotwg(npwepG0)
 937  complex(gwpc),intent(in) :: rhotwx(3,nspinor**2)
 938  complex(gwpc),intent(inout) :: sf_chi0(npwe,npwe,my_wl:my_wr)
 939  complex(dpc),intent(inout) :: sf_head(3,3,my_wl:my_wr)
 940  complex(dpc),intent(inout) :: sf_lwing(npwe,my_wl:my_wr,3)
 941  complex(dpc),intent(inout) :: sf_uwing(npwe,my_wl:my_wr,3)
 942 
 943 !Local variables-------------------------------
 944 !scalars
 945  integer :: itim,isym,idir,jdir
 946  complex(gwpc) :: num
 947  character(len=500) :: msg
 948 !arrays
 949  integer, ABI_CONTIGUOUS pointer :: Sm1G(:)
 950  complex(dpc) :: mir_kbz(3)
 951  complex(gwpc), ABI_CONTIGUOUS pointer :: phmGt(:)
 952  complex(gwpc),allocatable :: rhotwg_sym(:)
 953 
 954 !************************************************************************
 955 
 956  if (iomegal<my_wl .or. iomegar>my_wr) then
 957    write(msg,'(3a,2(a,i0,a,i0,a))')ch10,&
 958 &    'Indices out of boundary ',ch10,&
 959 &    '  my_wl = ',my_wl,' iomegal = ',iomegal,ch10,&
 960 &    '  my_wr = ',my_wr,' iomegar = ',iomegar,ch10
 961    MSG_BUG(msg)
 962  end if
 963 
 964  SELECT CASE (symchi)
 965  CASE (0)
 966    !
 967    ! Calculation without symmetries
 968    ! rhotwg(1)= R^-1q*rhotwx_ibz
 969    ! rhotwg(1)=-R^-1q*conjg(rhotwx_ibz) for inversion
 970    if (wl<huge(0.0_dp)*1.d-11) then
 971      !this is awful but it is still a first coding
 972      ! Num is single precision needed for cgerc check factocc
 973      num=-wl*factocc
 974      call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,sf_chi0(:,:,iomegal),npwe)
 975    end if
 976    ! Last point, must accumulate left point but not the right one
 977    if (iomegar/=nomegasf+1 .and. wr<huge(0.0_dp)*1.d-11) then
 978      num=-wr*factocc
 979      call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,sf_chi0(:,:,iomegar),npwe)
 980    end if
 981 
 982    ! ================================
 983    ! ==== Update heads and wings ====
 984    ! ================================
 985 
 986    ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
 987    if (nspinor == 1) then
 988       mir_kbz =(3-2*itim_kbz) * MATMUL(Cryst%symrec(:,:,isym_kbz),rhotwx(:,1))
 989     else
 990       mir_kbz = (3-2*itim_kbz) * MATMUL(Cryst%symrec(:,:,isym_kbz), sum(rhotwx(:,1:2), dim=2))
 991     end if
 992    if (itim_kbz==2) mir_kbz=CONJG(mir_kbz)
 993 
 994    do jdir=1,3
 995      if (wl<huge(0.0_dp)*1.d-11) then
 996        ! this is awful but it is still a first coding
 997        ! Num is single precision needed for cgerc check factocc
 998        num=-wl*factocc
 999        sf_uwing(:,iomegal,jdir) = sf_uwing(:,iomegal,jdir) + num * mir_kbz(jdir) * CONJG(rhotwg(1:npwepG0))
1000        sf_lwing(:,iomegal,jdir) = sf_lwing(:,iomegal,jdir) + num * rhotwg(1:npwepG0) * CONJG(mir_kbz(jdir))
1001        do idir=1,3
1002          sf_head(idir,jdir,iomegal) = sf_head(idir,jdir,iomegal) + num * mir_kbz(idir) * CONJG(mir_kbz(jdir))
1003        end do
1004      end if
1005 
1006      ! Last point, must accumulate left point but not the right one
1007      if (iomegar/=nomegasf+1 .and. wr<huge(0.0_dp)*1.d-11) then
1008        num=-wr*factocc
1009        sf_uwing(:,iomegar,jdir) = sf_uwing(:,iomegar,jdir) + num * mir_kbz(jdir) * CONJG(rhotwg(1:npwepG0))
1010        sf_lwing(:,iomegar,jdir) = sf_lwing(:,iomegar,jdir) + num * rhotwg(1:npwepG0) * CONJG(mir_kbz(jdir))
1011        do idir=1,3
1012          sf_head(idir,jdir,iomegar) = sf_head(idir,jdir,iomegar) + num * mir_kbz(idir) * CONJG(mir_kbz(jdir))
1013        end do
1014      end if
1015    end do ! jdir
1016 
1017  CASE (1)
1018    ! === Notes on the symmetrization of oscillator matrix elements ===
1019    ! If  Sq = q then  M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1G}  (k,q)
1020    ! If -Sq = q then  M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1G}^*(k,q)
1021    !
1022    ! In case of an umklapp process
1023    ! If  Sq = q+G_o then  M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1(G-G_o}   (k,q)
1024    ! If -Sq = q+G_o then  M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1(G-G-o)}^*(k,q)
1025    !
1026    ! rhotwg(1)= R^-1q*rhotwx_ibz
1027    ! rhotwg(1)=-R^-1q*conjg(rhotwx_ibz) for inversion
1028 
1029    ABI_MALLOC(rhotwg_sym,(npwe))
1030 
1031    ! Loop over symmetries of the space group and time-reversal
1032    do isym=1,Ltg_q%nsym_sg
1033      do itim=1,Ltg_q%timrev
1034 
1035        if (Ltg_q%wtksym(itim,isym,ikbz)==1) then
1036          ! This operation belongs to the little group and has to be considered to reconstruct the BZ
1037          phmGt => Gsph_epsG0%phmGt(1:npwe,isym) ! In these 2 lines mind the slicing (1:npwe)
1038          Sm1G  => Gsph_epsG0%rottbm1(1:npwe,itim,isym)
1039 
1040          SELECT CASE (itim)
1041          CASE (1)
1042            rhotwg_sym(1:npwe)=rhotwg(Sm1G(1:npwe))*phmGt(1:npwe)
1043          CASE (2)
1044            rhotwg_sym(1:npwe)=CONJG(rhotwg(Sm1G(1:npwe)))*phmGt(1:npwe)
1045          CASE DEFAULT
1046            MSG_BUG(sjoin('Wrong value of itim:', itoa(itim)))
1047          END SELECT
1048 
1049          ! Multiply elements G,Gp of rhotwg_sym*num and accumulate in sf_chi0(G,Gp,io)
1050          if (wl<huge(0.0_dp)*1.d-11) then
1051            num=-wl*factocc
1052            call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,sf_chi0(:,:,iomegal),npwe)
1053          end if
1054 
1055          ! Last point, must accumulate left point but not the right one
1056          if (iomegar/=nomegasf+1 .and. wr<huge(0.0_dp)*1.d-11) then
1057            num=-wr*factocc
1058            call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,sf_chi0(:,:,iomegar),npwe)
1059          end if
1060 
1061          ! Accumulate heads and wings.
1062          ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
1063          if (nspinor == 1) then
1064            mir_kbz =(3-2*itim) * MATMUL(Cryst%symrec(:,:,isym),rhotwx(:,1))
1065          else
1066            mir_kbz = (3-2*itim) * MATMUL(Cryst%symrec(:,:,isym), sum(rhotwx(:,1:2), dim=2))
1067          end if
1068          if (itim==2) mir_kbz=CONJG(mir_kbz)
1069 
1070          do jdir=1,3
1071            if (wl<huge(0.0_dp)*1.d-11) then
1072              ! this is awful but it is still a first coding
1073              ! Num is single precision needed for cgerc check factocc
1074              num=-wl*factocc
1075              sf_uwing(:,iomegal,jdir) = sf_uwing(:,iomegal,jdir) + num * mir_kbz(jdir) * CONJG(rhotwg_sym(1:npwe))
1076              sf_lwing(:,iomegal,jdir) = sf_lwing(:,iomegal,jdir) + num * rhotwg_sym(1:npwe) * CONJG(mir_kbz(jdir))
1077              do idir=1,3
1078                sf_head(idir,jdir,iomegal) = sf_head(idir,jdir,iomegal) + num * mir_kbz(idir) * CONJG(mir_kbz(jdir))
1079              end do
1080            end if
1081 
1082            ! Last point, must accumulate left point but not the right one
1083            if (iomegar/=nomegasf+1 .and. wr<huge(0.0_dp)*1.d-11) then
1084              num=-wr*factocc
1085              sf_uwing(:,iomegar,jdir) = sf_uwing(:,iomegar,jdir) + num * mir_kbz(jdir) * CONJG(rhotwg_sym(1:npwe))
1086              sf_lwing(:,iomegar,jdir) = sf_lwing(:,iomegar,jdir) + num * rhotwg_sym(1:npwe) * CONJG(mir_kbz(jdir))
1087              do idir=1,3
1088                sf_head(idir,jdir,iomegar) = sf_head(idir,jdir,iomegar) + num * mir_kbz(idir) * CONJG(mir_kbz(jdir))
1089              end do
1090            end if
1091          end do ! jdir
1092 
1093        end if !wtksym
1094      end do !inv
1095    end do !isym
1096    ABI_FREE(rhotwg_sym)
1097 
1098  CASE DEFAULT
1099    MSG_BUG(sjoin('Wrong value of symchi:', itoa(symchi)))
1100  END SELECT
1101 
1102 end subroutine accumulate_sfchi0_q0

m_chi0tk/approxdelta [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

  approxdelta

FUNCTION

  Approximate the Dirac function using two methods:
  method 1) a triangular funtion centered at the value egwdiff_re, Eq 17 of PRB 74, 035101 (2006) [[cite:Shishkin2006]]
  method 2) a gaussian of witdth ep%spsmear expandended in Taylor series
  (at the moment only the 0-th moments)

  Subroutine needed to implement the calculation
  of the polarizability using the spectral representation as proposed in:
  PRB 74, 035101 (2006) [[cite:Shishkin2006]]
  and PRB 61, 7172 (2000) [[cite:Miyake2000]]

INPUTS

  nomegasf=number of frequencies in the grid for Im \chi_0
  omegasf(0:nomega+1)= frequencies (real)
  egwdiff_re = transition energy where the delta function is centered

  method= 1: a triangular shaped function used to approximated the delta
          2: gaussian approximation with standard deviation (smear)
 smear= used only in case of method==2, defines the width of the gaussian

OUTPUT

  wl = weight associated to omegal (last omega wich is smaller than egwdiff_re
  wr = weight associate to omegar  (first omega larger than egwdff_re
  iomegal= index in the array omegasf of the last frequency < egwdiff
  iomegar= index in the array omegasf of the first frequency > egwdiff

PARENTS

      cchi0,cchi0q0

CHILDREN

SOURCE

1364 subroutine approxdelta(nomegasf,omegasf,egwdiff_re,smear,iomegal,iomegar,wl,wr,spmeth)
1365 
1366 
1367 !This section has been created automatically by the script Abilint (TD).
1368 !Do not modify the following lines by hand.
1369 #undef ABI_FUNC
1370 #define ABI_FUNC 'approxdelta'
1371 !End of the abilint section
1372 
1373  implicit none
1374 
1375 !Arguments ------------------------------------
1376 !scalars
1377  integer,intent(in) :: nomegasf,spmeth
1378  integer,intent(out) :: iomegal,iomegar
1379  real(dp),intent(in) :: egwdiff_re,smear
1380  real(dp),intent(out) :: wl,wr
1381 !arrays
1382  real(dp),intent(in) :: omegasf(nomegasf)
1383 
1384 !Local variables-------------------------------
1385  integer :: io,iomega
1386  real(dp) :: omegal,omegar,deltal,deltar
1387  !character(len=500) :: msg
1388 ! *************************************************************************
1389 
1390  iomega=-999
1391  do io=nomegasf,1,-1
1392    if (omegasf(io)<egwdiff_re) then
1393     iomega=io; EXIT
1394    end if
1395  end do
1396 
1397  iomegal=iomega   ; omegal=omegasf(iomegal)
1398  iomegar=iomegal+1; omegar=omegasf(iomegar)
1399 
1400  SELECT CASE (spmeth)
1401  CASE (1)
1402    ! Weights for triangular shaped function
1403    wr=  (egwdiff_re-omegal)/(omegar-omegal)
1404    wl= -(egwdiff_re-omegar)/(omegar-omegal)
1405 
1406  CASE (2)
1407    ! Weights for gaussian method (0-th moment)
1408    deltal=(egwdiff_re-omegal)/smear
1409    deltar=(omegar-egwdiff_re)/smear
1410    if (deltar>=deltal) then
1411      wl=EXP(-deltal*deltal)
1412      ! this value is used to avoid double counting and speed-up
1413      wr=huge(one)*1.d-10
1414    else
1415      wl=huge(one)*1.d-10
1416      wr=exp(-deltal*deltal)
1417    end if
1418 
1419  CASE DEFAULT
1420    MSG_BUG(sjoin('Wrong value for spmeth:', itoa(spmeth)))
1421  END SELECT
1422 
1423 end subroutine approxdelta

m_chi0tk/assemblychi0_sym [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 assemblychi0_sym

FUNCTION

 Update the independent particle susceptibility for the contribution
 of one pair of occupied-unoccupied band, for each frequency.
 If symchi=1 the expression is symmetrized taking into account the symmetries
 of the little group associated to the external q-point.
 Compute chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S \hat S (rhotwg(G1)*rhotwg*(G2))*green_w(io)
 where S are the symmetries of the little group associated to the external q-point.

INPUTS

  nspinor=Number of spinorial components.
  ik_bz=Index of the k-point in the BZ array whose contribution has to be symmetrized and added to cchi0
  npwepG0=Maximum number of G vectors taking into account possible umklapp G0, ie enlarged sphere G-G0
  rhotwg(npwe)=Oscillator matrix elements for this k-point and the transition that has to be summed
  green_w(nomega)=frequency dependent part coming from the green function
  Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors
    %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg
    %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion
    %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations
  Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point.
    %timrev=2 it time-reversal is used, 1 otherwise
    %nsym_sg=Number of space group symmetries
    %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point
    %flag_umklp(timrev,nsym)= flag for umklapp processes
      if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise
    %igmG0(npwepG0,timrev,nsym) index of G-G0 in the array gvec
  Ep<em1params_t>=Parameters related to the calculation of chi0/epsilon^-1
    %symchi
    %nomega=number of frequencies
    %npwe=number of plane waves for epsilon (input variable)

OUTPUT

  (see side effects)

SIDE EFFECTS

  chi0(npwe,npwe,nomega)=independent-particle susceptibility matrix in reciprocal space

PARENTS

      cchi0,cchi0q0_intraband

CHILDREN

SOURCE

118 subroutine assemblychi0_sym(ik_bz,nspinor,Ep,Ltg_q,green_w,npwepG0,rhotwg,Gsph_epsG0,chi0)
119 
120 
121 !This section has been created automatically by the script Abilint (TD).
122 !Do not modify the following lines by hand.
123 #undef ABI_FUNC
124 #define ABI_FUNC 'assemblychi0_sym'
125 !End of the abilint section
126 
127  implicit none
128 
129 !Arguments ------------------------------------
130 !scalars
131  integer,intent(in) :: ik_bz,npwepG0,nspinor
132  type(gsphere_t),intent(in) :: Gsph_epsG0
133  type(littlegroup_t),intent(in) :: Ltg_q
134  type(em1params_t),intent(in) :: Ep
135 !arrays
136  complex(gwpc),intent(in) :: rhotwg(npwepG0)
137  complex(dpc),intent(in) :: green_w(Ep%nomega)
138  complex(gwpc),intent(inout) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega)
139 
140 !Local variables-------------------------------
141 !scalars
142  integer :: itim,io,isym,ig1,ig2,nthreads
143  complex(gwpc) :: dd
144  !character(len=500) :: msg
145 !arrays
146  integer :: Sm1_gmG0(Ep%npwe)
147  complex(gwpc) :: rhotwg_sym(Ep%npwe)
148  !complex(gwpc),allocatable :: rhotwg_I(:),rhotwg_J(:)
149 
150 ! *************************************************************************
151 
152  ABI_UNUSED(nspinor)
153 
154  nthreads = xomp_get_max_threads()
155 
156  SELECT CASE (Ep%symchi)
157  CASE (0)
158    ! Do not use symmetries
159    if (nthreads==1 .or. Ep%nomega>=nthreads) then
160      ! MKL_10.3 xgerc is not threaded. BTW: single precision is faster (sometimes factor ~2).
161 !$omp parallel do private(dd)
162      do io=1,Ep%nomega
163        dd = green_w(io)
164        call XGERC(Ep%npwe,Ep%npwe,dd,rhotwg,1,rhotwg,1,chi0(:,:,io),Ep%npwe)
165      end do
166    else
167 !$omp parallel private(dd)
168      do io=1,Ep%nomega
169        dd = green_w(io)
170 !$omp do
171        do ig2=1,Ep%npwe
172          do ig1=1,Ep%npwe
173            chi0(ig1,ig2,io) = chi0(ig1,ig2,io) + dd * rhotwg(ig1) * GWPC_CONJG(rhotwg(ig2))
174          end do
175        end do
176 !$omp end do NOWAIT
177      end do
178 !$omp end parallel
179    end if
180 
181  CASE (1)
182    ! Use symmetries to reconstruct the integrand in the BZ.
183    !
184    ! Notes on the symmetrization of the oscillator matrix elements
185    !  If  Sq = q then  M_G^( Sk,q)= e^{-i(q+G).t} M_{ S^-1G}  (k,q)
186    !  If -Sq = q then  M_G^(-Sk,q)= e^{-i(q+G).t} M_{-S^-1G}^*(k,q)
187    !
188    ! In case of an umklapp process
189    !  If  Sq = q+G0 then  M_G( Sk,q)= e^{-i(q+G).t} M_{ S^-1(G-G0}   (k,q)
190    !  If -Sq = q+G0 then  M_G(-Sk,q)= e^{-i(q+G).t} M_{-S^-1(G-G0)}^*(k,q)
191    !
192    ! Ltg_q%igmG0(ig,itim,isym) contains the index of G-G0 where ISq=q+G0
193    ! Note that there is no need to take into account the phases due to q,
194    ! They cancel in the scalar product ==> phmGt(G,isym)=e^{-iG\cdot t}
195    !
196    ! Mind the slicing of %rottbm1(npwepG0,timrev,nsym) and %phmGt(npwepG0,nsym) as
197    ! these arrays, usually, do not conform to rho_twg_sym(npw) !
198    !
199    ! Loop over symmetries of the space group and time-reversal.
200    do isym=1,Ltg_q%nsym_sg
201      do itim=1,Ltg_q%timrev
202        if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
203          ! This operation belongs to the little group and has to be used to reconstruct the BZ ===
204          ! * In the following 3 lines mind the slicing (1:npwe)
205          ! TODO this is a hot-spot, should add a test on the umklapp
206          !
207          !gmG0 => Ltg_q%igmG0(1:Ep%npwe,itim,isym)
208          Sm1_gmG0(1:Ep%npwe) = Gsph_epsG0%rottbm1( Ltg_q%igmG0(1:Ep%npwe,itim,isym), itim,isym)
209 
210          SELECT CASE (itim)
211          CASE (1)
212            rhotwg_sym(1:Ep%npwe) = rhotwg(Sm1_gmG0) * Gsph_epsG0%phmGt(1:Ep%npwe,isym)
213          CASE (2)
214            rhotwg_sym(1:Ep%npwe) = GWPC_CONJG(rhotwg(Sm1_gmG0))*Gsph_epsG0%phmGt(1:Ep%npwe,isym)
215          CASE DEFAULT
216            MSG_BUG(sjoin('Wrong itim:', itoa(itim)))
217          END SELECT
218 
219          ! Multiply rhotwg_sym by green_w(io) and accumulate in chi0(G,Gp,io)
220          if (nthreads==1 .or. Ep%nomega>=nthreads) then
221            ! MKL_10.3 xgerc is not threaded. BTW: single precision is faster (sometimes factor ~2).
222 !$omp parallel do private(dd)
223            do io=1,Ep%nomega
224              dd=green_w(io)
225              call XGERC(Ep%npwe,Ep%npwe,dd,rhotwg_sym,1,rhotwg_sym,1,chi0(:,:,io),Ep%npwe)
226            end do
227 
228          else
229            !write(std_out,*)"in new with nthreads = ",nthreads
230 !$omp parallel private(dd)
231            do io=1,Ep%nomega
232              dd=green_w(io)
233 !$omp do
234              do ig2=1,Ep%npwe
235                do ig1=1,Ep%npwe
236                  chi0(ig1,ig2,io) = chi0(ig1,ig2,io) + dd * rhotwg_sym(ig1) * GWPC_CONJG(rhotwg_sym(ig2))
237                end do
238              end do
239 !$omp end do NOWAIT
240            end do
241 !$omp end parallel
242          end if
243        end if
244 
245      end do
246    end do
247 
248  CASE DEFAULT
249    MSG_BUG(sjoin('Wrong symchi:', itoa(Ep%symchi)))
250  END SELECT
251 
252 end subroutine assemblychi0_sym

m_chi0tk/assemblychi0sf [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 assemblychi0sf

FUNCTION

 Update the spectral function of the irreducible polarizability for the contribution
 of one pair of occupied-unoccupied states, for each frequency.
 If symchi==1, the symmetries of the little group of the external q-point are used
 to symmetrize the contribution in the full Brillouin zone. In this case, the routine computes:

   $ chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))*\delta(w - trans) $

 where S are the symmetries of the little group of the external q-point.

INPUTS

  ik_bz=Index of the k-point in the BZ whose contribution has to be added to the spectral function of chi0
    If symchi=1, the contribution is symmetrized.
  my_wl,my_wr=min and Max frequency index treated by this processor.
  npwe=Number of plane waves used to describe chi0.
  npwepG0=Maximum number of G vectors taking into account umklapp vectors.
  nomegasf=Number of frequencies for the spectral function.
  nspinor=Number of spinorial components.
  symchi=1 if symmetries are used, 0 otherwise
  rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states.
  timrev=if 2, time reversal has to be used to obtain k_bz; 1 otherwise.
  Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors
    %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg
    %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion
    %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations
  Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point.
    %timrev=2 it time-reversal is used, 1 otherwise
    %nsym_sg=Number of space group symmetries
    %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point
    %flag_umklp(timrev,nsym)= flag for umklapp processes
      if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise
    %igmG0(npwepG0,timrev,nsym) index of G-G0 in the array gvec
  factocc=occupation factor=f_occ*(ockp-occk) (see cchi0.F90)
  wl,wr=Weights used to approximate the delta function.

OUTPUT

  (see side effects)

SIDE EFFECTS

  chi0sf(npwe,npwe,my_wl:my_wr)= updated spectral function.

NOTES

  Umklapp processes are not yet implemented

PARENTS

      cchi0

CHILDREN

SOURCE

1162 subroutine assemblychi0sf(ik_bz,symchi,Ltg_q,npwepG0,npwe,rhotwg,Gsph_epsG0,&
1163 & factocc,my_wl,iomegal,wl,my_wr,iomegar,wr,nomegasf,chi0sf)
1164 
1165 
1166 !This section has been created automatically by the script Abilint (TD).
1167 !Do not modify the following lines by hand.
1168 #undef ABI_FUNC
1169 #define ABI_FUNC 'assemblychi0sf'
1170 !End of the abilint section
1171 
1172  implicit none
1173 
1174 !Arguments ------------------------------------
1175 !scalars
1176  integer,intent(in) :: ik_bz,iomegal,iomegar,my_wl,my_wr,nomegasf,npwe,npwepG0
1177  integer,intent(in) :: symchi
1178  real(dp),intent(in) :: factocc,wl,wr
1179  type(gsphere_t),intent(in) :: Gsph_epsG0
1180  type(littlegroup_t),intent(in) :: Ltg_q
1181 !arrays
1182  complex(gwpc),intent(in) :: rhotwg(npwepG0)
1183  complex(gwpc),intent(inout) :: chi0sf(npwe,npwe,my_wl:my_wr)
1184 
1185 !Local variables-------------------------------
1186 !scalars
1187  integer :: isym,itim,ig1,ig2
1188  complex(gwpc) :: num
1189  character(len=500) :: msg
1190 !arrays
1191  integer :: Sm1_gmG0(npwe)
1192  complex(gwpc) :: rhotwg_sym(npwe)
1193 
1194 
1195 ! *************************************************************************
1196 
1197  if (iomegal < my_wl .or. iomegar > my_wr) then
1198    write(msg,'(3a,2(a,i0,a,i0,a))')ch10,&
1199 &    ' Indices out of boundary ',ch10,&
1200 &    '  my_wl = ',my_wl,' iomegal = ',iomegal,ch10,&
1201 &    '  my_wr = ',my_wr,' iomegar = ',iomegar,ch10
1202    MSG_BUG(msg)
1203  end if
1204 
1205  SELECT CASE (symchi)
1206  CASE (0)
1207     ! Do not use symmetries.
1208 
1209 ! MG: This is the best I can do for this part.
1210 !$omp PARALLEL private(num)
1211 !$omp SECTIONS
1212 !$omp SECTION
1213    if (wl<huge(0.0_dp)*1.d-11) then !FIXME this is awful
1214      num=-wl*factocc
1215      call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,chi0sf(:,:,iomegal),npwe)
1216    end if
1217 
1218    ! Last point, must accumulate left point but not the right one
1219 !$omp SECTION
1220    if (iomegar/=nomegasf+1 .and. wr<huge(0.0_dp)*1.d-11) then
1221      num=-wr*factocc
1222      call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,chi0sf(:,:,iomegar),npwe)
1223    end if
1224 !$omp end SECTIONS
1225 !$omp end PARALLEL
1226 
1227  CASE (1)
1228    ! Use symmetries to reconstruct oscillator matrix elements
1229    ! Notes on the symmetrization of the oscillator maxtri elements:
1230    !
1231    ! If  Sq=q then  M_G^( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1G}  (k,q)
1232    ! If -Sq=q then  M_G^(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1G}^*(k,q)
1233    !
1234    ! In case of an umklapp process
1235    ! If  Sq=q+G_o then  M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1(G-G_o}   (k,q)
1236    ! If -Sq=q+G_o then  M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1(G-G_o)}^*(k,q)
1237    !
1238    ! Ltg_q%igmG0(ig,itim,isym) contains the index of G-G0 where ISq=q+G0
1239    ! Note that there is no need to take into account the phases due to q,
1240    ! They cancel in the scalar product ==> phmGt(G,isym)=e^{-iG\cdot t}
1241    !
1242    ! Mind the slicing of %rottbm1(npwepG0,timrev,nsym) and %phgt(npwepG0,nsym) as
1243    ! these arrays, usually, do not conform to rho_twg_sym(npw) !
1244    !
1245    !ABI_MALLOC(rhotwg_sym,(npwe))
1246    !
1247    ! Loop over symmetries of the space group and time-reversal
1248    do isym=1,Ltg_q%nsym_sg
1249      do itim=1,Ltg_q%timrev
1250 
1251        if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
1252          ! This operation belongs to the little group and has to be used to reconstruct BZ.
1253          ! TODO this is a hot-spot, should add a test on the umklapp
1254          !
1255          ! In these 3 lines mind the slicing (1:npwe)
1256          Sm1_gmG0(1:npwe)=Gsph_epsG0%rottbm1( Ltg_q%igmG0(1:npwe,itim,isym), itim,isym)
1257 
1258          SELECT CASE (itim)
1259          CASE (1)
1260            rhotwg_sym(1:npwe)=rhotwg(Sm1_gmG0(1:npwe)) * Gsph_epsG0%phmGt(1:npwe,isym)
1261          CASE (2)
1262            rhotwg_sym(1:npwe)=CONJG(rhotwg(Sm1_gmG0(1:npwe))) * Gsph_epsG0%phmGt(1:npwe,isym)
1263          CASE DEFAULT
1264            MSG_BUG(sjoin('Wrong value for itim:', itoa(itim)))
1265          END SELECT
1266 
1267 #if 0
1268 !! MG: This is the best I can do, at present.
1269 !$omp PARALLEL private(num)
1270 !$omp SECTIONS
1271 
1272 !$omp SECTION
1273          if (wl<huge(0.0_dp)*1.d-11) then
1274            num=-wl*factocc
1275            call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe)
1276          end if
1277 !$omp SECTION
1278          !
1279          ! Last point, must accumulate left point but not the right one
1280          if (iomegar/=nomegasf+1 .and. wr<huge(0.0_dp)*1.d-11) then
1281            num=-wr*factocc
1282            call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegar),npwe)
1283          end if
1284 !$omp end SECTIONS
1285 !$omp end PARALLEL
1286 #else
1287 
1288          if (wl<huge(0.0_dp)*1.d-11) then
1289            !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe)
1290            num=-wl*factocc
1291 !$omp parallel do
1292            do ig2=1,npwe
1293              do ig1=1,npwe
1294                chi0sf(ig1,ig2,iomegal) = chi0sf(ig1,ig2,iomegal) + num * rhotwg_sym(ig1) * CONJG(rhotwg_sym(ig2))
1295              end do
1296            end do
1297          end if
1298 
1299          ! Last point, must accumulate left point but not the right one
1300          if (iomegar/=nomegasf+1 .and. wr<huge(0.0_dp)*1.d-11) then
1301            !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegar),npwe)
1302            num=-wr*factocc
1303 !$omp parallel do
1304            do ig2=1,npwe
1305              do ig1=1,npwe
1306                !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe)
1307                chi0sf(ig1,ig2,iomegar) = chi0sf(ig1,ig2,iomegar) + num * rhotwg_sym(ig1) * CONJG(rhotwg_sym(ig2))
1308              end do
1309            end do
1310          end if
1311 #endif
1312        end if !wtksym
1313 
1314      end do !inv
1315    end do !isym
1316    !ABI_FREE(rhotwg_sym)
1317 
1318  CASE DEFAULT
1319    MSG_BUG(sjoin('Wrong value for symchi:', itoa(symchi)))
1320  END SELECT
1321 
1322 end subroutine assemblychi0sf

m_chi0tk/calc_kkweight [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

  calc_kkweight

FUNCTION

  Calculate frequency dependent weights needed to perform the Hilbert transform

  Subroutine needed to implement the calculation
  of the polarizability using the spectral representation as proposed in:
  PRB 74, 035101 (2006) [[cite:Shishkin2006]]
  and PRB 61, 7172 (2000) [[cite:Miyake2000]]

INPUTS

 nsp=number of frequencies where the imaginary part of the polarizability is evaluated
 ne=number of frequencies for the polarizability (same as in epsilon^-1)
 omegasp(nsp)=real frequencies for the imaginary part of the polarizability
 omegae(ne)= imaginary frequencies for the polarizability
 delta=small imaginary part used to avoid poles, input variables

OUTPUT

 kkweight(nsp,ne)=frequency dependent weights Eq A1 PRB 74, 035101 (2006) [[cite:Shishkin2006]]

PARENTS

      m_chi0

CHILDREN

SOURCE

1458 subroutine calc_kkweight(ne,omegae,nsp,omegasp,delta,omegamax,kkw)
1459 
1460 
1461 !This section has been created automatically by the script Abilint (TD).
1462 !Do not modify the following lines by hand.
1463 #undef ABI_FUNC
1464 #define ABI_FUNC 'calc_kkweight'
1465 !End of the abilint section
1466 
1467  implicit none
1468 
1469 !Arguments ------------------------------------
1470 !scalars
1471  integer,intent(in) :: ne,nsp
1472  real(dp),intent(in) :: delta,omegamax
1473 !arrays
1474  real(dp),intent(in) :: omegasp(nsp)
1475  complex(dpc),intent(in) :: omegae(ne)
1476  complex(dpc),intent(out) :: kkw(nsp,ne)
1477 
1478 !Local variables-------------------------------
1479 !scalars
1480  integer :: isp,je
1481  real(dp) :: eta,xx1,xx2,den1,den2
1482  complex(dpc) :: c1,c2,wt
1483 !************************************************************************
1484 
1485  DBG_ENTER("COLL")
1486 
1487  kkw(:,:)=czero
1488 
1489  do je=1,ne
1490    eta=delta
1491    wt=omegae(je)
1492    ! Not include shift at omega==0, what about metallic systems?
1493    if (abs(real(omegae(je)))<tol6 .and. abs(aimag(wt))<tol6) eta=tol12
1494    !  Not include shift along the imaginary axis
1495    if (abs(aimag(wt))>tol6) eta=zero
1496    do isp=1,nsp
1497      if (isp==1) then
1498        ! Skip negative point, should check that this would not lead to spurious effects
1499        c1=czero
1500        den1=one
1501      else
1502        xx1=omegasp(isp-1)
1503        xx2=omegasp(isp)
1504        den1= xx2-xx1
1505        c1= -(wt-xx1+j_dpc*eta)*log( (wt-xx2+j_dpc*eta)/(wt-xx1+j_dpc*eta) )&
1506 &          +(wt+xx1-j_dpc*eta)*log( (wt+xx2-j_dpc*eta)/(wt+xx1-j_dpc*eta) )
1507        c1= c1/den1
1508      end if
1509      xx1=omegasp(isp)
1510      if (isp==nsp) then
1511        ! Skip last point should check that this would not lead to spurious effects
1512        xx2=omegamax
1513      else
1514        xx2=omegasp(isp+1)
1515      end if
1516      den2=xx2-xx1
1517      c2=  (wt-xx2+j_dpc*eta)*log( (wt-xx2+j_dpc*eta)/(wt-xx1+j_dpc*eta) )&
1518 &        -(wt+xx2-j_dpc*eta)*log( (wt+xx2-j_dpc*eta)/(wt+xx1-j_dpc*eta) )
1519      c2= c2/den2
1520      kkw(isp,je)=  c1/den1 + c2/den2
1521    end do
1522  end do
1523 
1524  DBG_EXIT("COLL")
1525 
1526 end subroutine calc_kkweight

m_chi0tk/chi0_bbp_mask [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

  chi0_bbp_mask

FUNCTION

INPUTS

OUTPUT

PARENTS

      cchi0,cchi0q0

CHILDREN

SOURCE

2394 subroutine chi0_bbp_mask(Ep,use_tr,QP_BSt,mband,ikmq_ibz,ik_ibz,spin,spin_fact,bbp_mask)
2395 
2396 
2397 !This section has been created automatically by the script Abilint (TD).
2398 !Do not modify the following lines by hand.
2399 #undef ABI_FUNC
2400 #define ABI_FUNC 'chi0_bbp_mask'
2401 !End of the abilint section
2402 
2403  implicit none
2404 
2405 !Arguments ------------------------------------
2406 !scalars
2407  integer,intent(in) :: spin,ik_ibz,ikmq_ibz,mband
2408  real(dp),intent(in) :: spin_fact
2409  logical,intent(in) :: use_tr
2410  type(em1params_t),intent(in) :: Ep
2411  type(ebands_t),target,intent(in) :: QP_BSt
2412 !arrays
2413  logical,intent(out) :: bbp_mask(mband,mband)
2414 
2415 !Local variables-------------------------------
2416 !scalars
2417  integer :: ib1,ib2
2418  real(dp) :: deltaeGW_b1kmq_b2k,deltaf_b1kmq_b2k,e_b1_kmq,f_b1_kmq
2419 !arrays
2420  real(dp), ABI_CONTIGUOUS pointer :: qp_energy(:,:,:),qp_occ(:,:,:)
2421 
2422 !************************************************************************
2423 
2424  qp_energy => QP_BSt%eig; qp_occ => QP_BSt%occ
2425  bbp_mask=.FALSE.
2426 
2427  !use_tr = (Ep%awtr==1)
2428  SELECT CASE (Ep%gwcomp)
2429  CASE (0)
2430    do ib1=1,Ep%nbnds
2431      ! Loop over "conduction" states.
2432      e_b1_kmq = qp_energy(ib1,ikmq_ibz,spin)
2433      f_b1_kmq =    qp_occ(ib1,ikmq_ibz,spin)
2434 
2435      do ib2=1,Ep%nbnds ! Loop over "valence" states.
2436        deltaf_b1kmq_b2k   = spin_fact*(f_b1_kmq-qp_occ(ib2,ik_ibz,spin))
2437        deltaeGW_b1kmq_b2k = e_b1_kmq-qp_energy(ib2,ik_ibz,spin)
2438 
2439        SELECT CASE (Ep%spmeth)
2440        CASE (0)
2441          ! Standard Adler-Wiser expression.
2442          if (ABS(deltaf_b1kmq_b2k) >= GW_TOL_DOCC) then
2443            bbp_mask(ib1,ib2)=.TRUE.
2444            if (use_tr .and. ib1<ib2) bbp_mask(ib1,ib2)=.FALSE. ! GAIN a factor ~2 thanks to time-reversal.
2445          end if
2446 
2447        CASE (1,2)
2448          ! Spectral method, WARNING time-reversal here is always assumed!
2449          if (ABS(deltaf_b1kmq_b2k) >= GW_TOL_DOCC) then
2450            bbp_mask(ib1,ib2)=.TRUE.
2451            if (deltaeGW_b1kmq_b2k<zero) bbp_mask(ib1,ib2)=.FALSE. ! Only positive frequencies are needed for the Hilbert transform.
2452            !$if (use_tr .and. ib1<ib2) bbp_mask(ib1,ib2)=.FALSE. ! GAIN a factor ~2 thanks to time-reversal.
2453          end if
2454 
2455        CASE DEFAULT
2456          MSG_ERROR(sjoin(" Wrong value for spmeth:", itoa(Ep%spmeth)))
2457        END SELECT
2458        !write(std_out,*) "bbp_mask(ib1,ib2)",bbp_mask(ib1,ib2)
2459      end do !ib2
2460    end do !ib1
2461 
2462  CASE (1)
2463    ! Extrapolar technique
2464    ABI_CHECK(Ep%spmeth==0,"Hilbert transform and extrapolar method are not compatible")
2465 
2466    ! Loop over "conduction" states.
2467    do ib1=1,Ep%nbnds
2468      e_b1_kmq=qp_energy(ib1,ikmq_ibz,spin)
2469      f_b1_kmq=   qp_occ(ib1,ikmq_ibz,spin)
2470 
2471      ! Loop over "valence" states.
2472      do ib2=1,Ep%nbnds
2473        deltaf_b1kmq_b2k  =spin_fact*(f_b1_kmq-qp_occ(ib2,ik_ibz,spin))
2474        deltaeGW_b1kmq_b2k=e_b1_kmq-qp_energy(ib2,ik_ibz,spin)
2475 
2476        ! When the completeness correction is used,
2477        ! we need to also consider transitions with vanishing deltaf
2478        ! Rangel: This is to compute chi in metals correctly with the extrapolar method.
2479        bbp_mask(ib1,ib2)=.TRUE.
2480        !if (qp_occ(ib2,ik_ibz,is) < GW_TOL_DOCC) CYCLE
2481        if (qp_occ(ib2,ik_ibz,spin) < GW_TOL_DOCC .and. (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC .or. ib1<ib2)) then
2482          bbp_mask(ib1,ib2)=.FALSE.
2483        end if
2484 
2485      end do
2486    end do
2487 
2488   CASE DEFAULT
2489     MSG_ERROR(sjoin("Wrong value of gwcomp:", itoa(Ep%gwcomp)))
2490   END SELECT
2491 
2492 end subroutine chi0_bbp_mask

m_chi0tk/completechi0_deltapart [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 completechi0_deltapart

FUNCTION

  Apply the delta part of the completeness correction to chi0

INPUTS

  ik_bz=Index of the k-point in the full BZ whose contribution has to be added and symmetrized.
  qzero=.TRUE. is long wave-length limit.
  symchi=1 if we are summing over IBZ_q and symmetrization has to be performed.
  npwe=Number of G vectors in chi0.
  npwvec=MAX number of G.
  nomega=Number of frequencies.
  nspinor=Number of spinorial components.
  nfftot=Total Number of points in the FFT
  ngfft(18)=Info on the FFT.
  igfft0(npwvec)=Index of each G in the FFT array.
  Gsph_FFT=<gsphere_t>=Info on the largest G-sphere contained in the FFT box used for wavefunctions.
  Ltg_q=<littlegroup_t>= Structure gathering information on the little group of the external q.
  green_enhigh_w=Approximated frequency dependent part of the Green function entering equation (TODO put reference)
  wfwfg=Fourier components of u_{kb1}.u_{kb2}

OUTPUT

  See SIDES EFFECTS

 SIDES EFFECTS
  chi0(npwe,npwe,nomega)= In input chi0 calculated so far,
  In output the "delta part" of the completeness correction is added.

PARENTS

      cchi0,cchi0q0

CHILDREN

SOURCE

1912 subroutine completechi0_deltapart(ik_bz,qzero,symchi,npwe,npwvec,nomega,nspinor,&
1913 & nfftot,ngfft,igfft0,Gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0)
1914 
1915 
1916 !This section has been created automatically by the script Abilint (TD).
1917 !Do not modify the following lines by hand.
1918 #undef ABI_FUNC
1919 #define ABI_FUNC 'completechi0_deltapart'
1920 !End of the abilint section
1921 
1922  implicit none
1923 
1924 !Arguments ------------------------------------
1925 !scalars
1926  integer,intent(in) :: ik_bz,nfftot,nomega,npwe,npwvec,nspinor,symchi
1927  logical,intent(in) :: qzero
1928  type(gsphere_t),intent(in) :: Gsph_FFT
1929  type(littlegroup_t),intent(in) :: Ltg_q
1930 !arrays
1931  integer,intent(in) :: igfft0(npwvec),ngfft(18)
1932  complex(dpc),intent(in) :: green_enhigh_w(nomega)
1933  complex(gwpc),intent(in) :: wfwfg(nfftot*nspinor**2)
1934  complex(gwpc),intent(inout) :: chi0(npwe,npwe,nomega)
1935 
1936 !Local variables ------------------------------
1937 !scalars
1938  integer,save :: enough=0
1939  integer :: iSm1_g1mg2,iSm1_g1mg2_fft,ig,gmg_sph,gmg_fft
1940  integer :: igp,igstart,isym,itim,outofbox_wfn
1941  complex(gwpc) :: phmGt
1942  !character(len=500) :: msg
1943 
1944 !************************************************************************
1945 
1946  igstart=1; if (qzero) igstart=2
1947  outofbox_wfn=0
1948 
1949  SELECT CASE (symchi)
1950 
1951  CASE (0) ! Do not use symmetries.
1952    ! MG: One has to make sure G1-G2 is still in the FFT mesh for each G1 and G2 in chi0 (not always true)
1953    ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic!
1954    do igp=igstart,npwe
1955      do ig=igstart,npwe
1956        gmg_fft = gsph_gmg_fftidx(Gsph_FFT,ig,igp,ngfft)
1957        if (gmg_fft==0) then
1958          outofbox_wfn=outofbox_wfn+1; CYCLE
1959        end if
1960        chi0(ig,igp,:) = chi0(ig,igp,:) + wfwfg(gmg_fft)*green_enhigh_w(:)
1961      end do
1962    end do
1963 
1964  CASE (1)
1965    ! Symmetrize the integrand in the full BZ.
1966    ! * <Sk b|e^{-i(G1-G2}.r}|b Sk> = e^{-i(G1-G2).\tau} <k b|e^{-i(S^{-1}(G1-G2).r)|b k>
1967    ! * green_enhigh_w in invariant under symmetry
1968    ! * We symmetrize using the operations of the little group of q since this routine
1969    !   is called inside a sum over IBZ_q, it would be possible to symmetrize
1970    !   this term by just summing over the IBZ and rotating the matrix elements.
1971    ! * Time-reversal does not lead to a complex conjugated since bra and ket are the same.
1972    !
1973    do igp=igstart,npwe
1974      do ig=igstart,npwe
1975 
1976       ! Get the index of G1-G2.
1977       gmg_sph = gsph_gmg_idx(Gsph_FFT,ig,igp)
1978       if (gmg_sph==0) then
1979         outofbox_wfn=outofbox_wfn+1; CYCLE
1980       end if
1981 
1982       do itim=1,Ltg_q%timrev
1983         do isym=1,Ltg_q%nsym_sg
1984           if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
1985             ! * This operation belongs to the little group and has to be used to reconstruct the BZ.
1986             ! * Time-reversal in not used to rotate (G1-G2) see comment above.
1987             phmGt          = Gsph_FFT%phmGt  (gmg_sph,isym)
1988             iSm1_g1mg2     = Gsph_FFT%rottbm1(gmg_sph,1,isym)
1989             iSm1_g1mg2_fft = igfft0(iSm1_g1mg2)
1990 
1991             chi0(ig,igp,:) = chi0(ig,igp,:) + phmGt*wfwfg(iSm1_g1mg2_fft)*green_enhigh_w(:)
1992           end if
1993         end do !isym
1994       end do !itim
1995 
1996      end do !igp
1997    end do !ig
1998 
1999  CASE DEFAULT
2000    MSG_BUG("Wrong value of symchi")
2001  END SELECT
2002 
2003  if (outofbox_wfn/=0) then
2004    enough=enough+1
2005    if (enough<=50) then
2006      MSG_WARNING(sjoin(' Number of G1-G2 pairs outside the G-sphere for Wfns: ', itoa(outofbox_wfn)))
2007      if (enough==50) then
2008        call wrtout(std_out,' ========== Stop writing Warnings ==========','COLL')
2009      end if
2010    end if
2011  end if
2012 
2013 end subroutine completechi0_deltapart

m_chi0tk/hilbert_transform [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

  hilbert_transform

FUNCTION

  Compute the hilbert transform.

INPUTS

 nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$
 nomega=number of frequencies in $\chi0(q,\omega)$.
 max_rest,min_res=max and min resonant transition energy (for this q-point)
 my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor

OUTPUT

PARENTS

      cchi0,cchi0q0

CHILDREN

SOURCE

1719 subroutine hilbert_transform(npwe,nomega,nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,spmeth)
1720 
1721 
1722 !This section has been created automatically by the script Abilint (TD).
1723 !Do not modify the following lines by hand.
1724 #undef ABI_FUNC
1725 #define ABI_FUNC 'hilbert_transform'
1726 !End of the abilint section
1727 
1728  implicit none
1729 
1730 !Arguments ------------------------------------
1731 !scalars
1732  integer,intent(in) :: spmeth,nomega,nomegasf,my_wl,my_wr,npwe
1733 !arrays
1734  complex(dpc),intent(in) :: kkweight(nomegasf,nomega)
1735  complex(gwpc),intent(inout) :: sf_chi0(npwe,npwe,my_wl:my_wr)
1736  complex(gwpc),intent(inout) :: chi0(npwe,npwe,nomega)
1737 
1738 !Local variables-------------------------------
1739 !scalars
1740  integer :: ig2,my_nwp
1741  character(len=500) :: msg
1742 !arrays
1743  complex(gwpc),allocatable :: A_g1wp(:,:),H_int(:,:),my_kkweight(:,:)
1744 
1745 !************************************************************************
1746 
1747 #ifdef HAVE_OPENMP
1748  write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform (with OpenMP) using method ',spmeth,' It might take some time...'
1749 #else
1750  write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform using method ',spmeth,' It might take some time...'
1751 #endif
1752  call wrtout(std_out, msg, do_flush=.True.)
1753 
1754  my_nwp = my_wr - my_wl +1
1755 
1756 !$omp parallel private(my_kkweight, A_g1wp, H_int, ig2)
1757  ABI_MALLOC(my_kkweight, (my_wl:my_wr,nomega))
1758  my_kkweight = kkweight(my_wl:my_wr,:)
1759 
1760  ABI_MALLOC(A_g1wp, (npwe, my_nwp))
1761  ABI_MALLOC(H_int, (npwe, nomega))
1762 
1763 !$omp do
1764  do ig2=1,npwe
1765    A_g1wp = sf_chi0(:,ig2,:)
1766 
1767    ! Compute H_int = MATMUL(A_g1wp,my_kkweight)
1768    call XGEMM('N','N',npwe,nomega,my_nwp,cone_gw,A_g1wp,npwe,my_kkweight,my_nwp,czero_gw,H_int,npwe)
1769    chi0(:,ig2,:) = H_int
1770  end do
1771 
1772  ABI_FREE(my_kkweight)
1773  ABI_FREE(A_g1wp)
1774  ABI_FREE(H_int)
1775 !$omp end parallel
1776 
1777 end subroutine hilbert_transform

m_chi0tk/hilbert_transform_headwings [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

  hilbert_transform_headwings

FUNCTION

  Compute the hilbert transform the heads and wings of the polarizability.

INPUTS

 nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$
 nomega=number of frequencies in $\chi0(q,\omega)$.
 max_rest,min_res=max and min resonant transition energy (for this q-point)
 my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor

OUTPUT

PARENTS

      cchi0q0

CHILDREN

SOURCE

1804 subroutine hilbert_transform_headwings(npwe,nomega,nomegasf,my_wl,my_wr,kkweight, &
1805 & sf_lwing,sf_uwing,sf_head,chi0_lwing,chi0_uwing,chi0_head,spmeth)
1806 
1807 
1808 !This section has been created automatically by the script Abilint (TD).
1809 !Do not modify the following lines by hand.
1810 #undef ABI_FUNC
1811 #define ABI_FUNC 'hilbert_transform_headwings'
1812 !End of the abilint section
1813 
1814  implicit none
1815 
1816 !Arguments ------------------------------------
1817 !scalars
1818  integer,intent(in) :: spmeth,nomega,nomegasf,my_wl,my_wr,npwe
1819 !arrays
1820  complex(dpc),intent(in) :: kkweight(nomegasf,nomega)
1821  complex(dpc),intent(inout) :: sf_lwing(npwe,my_wl:my_wr,3)
1822  complex(dpc),intent(inout) :: sf_uwing(npwe,my_wl:my_wr,3)
1823  complex(dpc),intent(inout) :: sf_head(3,3,my_wl:my_wr)
1824  complex(dpc),intent(inout) :: chi0_lwing(npwe,nomega,3)
1825  complex(dpc),intent(inout) :: chi0_uwing(npwe,nomega,3)
1826  complex(dpc),intent(inout) :: chi0_head(3,3,nomega)
1827 
1828 !Local variables-------------------------------
1829 !scalars
1830  integer :: ig1,idir,io,iw
1831  complex(dpc) :: kkw
1832  character(len=500) :: msg
1833 !************************************************************************
1834 
1835 #ifdef HAVE_OPENMP
1836  write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform (with OpenMP) using method ',spmeth,' It might take some time...'
1837 #else
1838  write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform using method ',spmeth,' It might take some time...'
1839 #endif
1840  call wrtout(std_out,msg,'COLL',do_flush=.True.)
1841 
1842  ! Hilbert transform of the head.
1843  do io=1,nomega
1844    chi0_head(1,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,1,my_wl:my_wr))
1845    chi0_head(2,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,1,my_wl:my_wr))
1846    chi0_head(3,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,1,my_wl:my_wr))
1847    chi0_head(1,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,2,my_wl:my_wr))
1848    chi0_head(2,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,2,my_wl:my_wr))
1849    chi0_head(3,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,2,my_wl:my_wr))
1850    chi0_head(1,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,3,my_wl:my_wr))
1851    chi0_head(2,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,3,my_wl:my_wr))
1852    chi0_head(3,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,3,my_wl:my_wr))
1853  end do
1854 
1855  ! Hilbert transform for wings.
1856  ! Partial contributions to chi0 will be summed afterwards.
1857 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(kkw)
1858  do idir=1,3
1859    do io=1,nomega
1860      do iw=my_wl,my_wr
1861        kkw = kkweight(iw,io)
1862        do ig1=1,npwe
1863          chi0_lwing(ig1,io,idir) = chi0_lwing(ig1,io,idir) + kkw*sf_lwing(ig1,iw,idir)
1864          chi0_uwing(ig1,io,idir) = chi0_uwing(ig1,io,idir) + kkw*sf_uwing(ig1,iw,idir)
1865        end do
1866      end do
1867    end do
1868  end do  !idir
1869 
1870 end subroutine hilbert_transform_headwings

m_chi0tk/make_transitions [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 make_transitions

FUNCTION

  Calculate transition energies entering the espression for the irreducible polarizability.

INPUTS

  nsspol=1 for spin unpolarized, 2 for spin polarized calculations
  nbnds=total number of bands
  kmesh<kmesh_t>=datatype gathering info on the k-mesh:
   | %nbz=number of k-points in the full BZ
   | %nibz=number of k-points in the IBZ
   | %tab(nkbz)=table giving for each k-point in the BZ, the corresponding irreducible point in the IBZ array
   | %bz(3,nkbz)=reduced coordinated of k-points
  TOL_DELTA_OCC=tolerance on the difference of the occupation numbers
  gw_energy(nbnds,kmesh%nkibz,nsppol)=quasi-particle energies energies
  occ(nbnds,kmesh%nkibz,nsppol)=occupation numbers
  chi0alg=integer defining the method used to calculate chi0
   0 ==> calculate chi0 using the Adler-Wiser expression
   1 ==> use spectral method
  timrev=if 2, time-reversal symmetry is considered; 1 otherwise

OUTPUT

 my_max_rest,my_min_rest=Maximum and minimum resonant (posite) transition energy.
 max_rest,min_rest=Maximun and minimum resonant (posite) transition energy treated by this node.

PARENTS

      cchi0,cchi0q0

CHILDREN

SOURCE

2236 subroutine make_transitions(Wfd,chi0alg,nbnds,nbvw,nsppol,symchi,timrev,TOL_DELTA_OCC,&
2237 & max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,gw_energy,occ,qpoint,bbp_ks_distrb)
2238 
2239 
2240 !This section has been created automatically by the script Abilint (TD).
2241 !Do not modify the following lines by hand.
2242 #undef ABI_FUNC
2243 #define ABI_FUNC 'make_transitions'
2244 !End of the abilint section
2245 
2246  implicit none
2247 
2248 !Arguments ------------------------------------
2249 !scalars
2250  integer,intent(in) :: chi0alg,nbnds,nbvw,nsppol,symchi,timrev
2251  real(dp),intent(in) :: TOL_DELTA_OCC
2252  real(dp),intent(out) :: max_rest,min_rest
2253  real(dp),intent(out) :: my_max_rest,my_min_rest
2254  type(kmesh_t),intent(in) :: Kmesh
2255  type(littlegroup_t),intent(in) :: Ltg_q
2256  type(wfd_t),intent(in) :: Wfd
2257 !arrays
2258  real(dp),intent(in) :: gw_energy(nbnds,Kmesh%nibz,nsppol)
2259  real(dp),intent(in) :: occ(nbnds,Kmesh%nibz,nsppol),qpoint(3)
2260  integer,intent(in) :: bbp_ks_distrb(Wfd%mband,Wfd%mband,Kmesh%nbz,Wfd%nsppol)
2261 
2262 !Local variables-------------------------------
2263 !scalars
2264  integer :: ib1,ib2,ii,ik_bz,ik_ibz,ikmq_bz,ikmq_ibz,is,nt,ntrans,my_ntrans,iloop
2265  real(dp) :: delta_ene,delta_occ,spin_fact
2266  character(len=500) :: msg
2267 !arrays
2268  integer :: G0(3)
2269  real(dp) :: kmq(3)
2270 
2271 !************************************************************************
2272 
2273  DBG_ENTER("COLL")
2274 
2275  if (chi0alg<0 .or. chi0alg>=2) then
2276    MSG_BUG(sjoin('chi0alg:', itoa(chi0alg),' not allowed'))
2277  end if
2278  if (timrev/=1 .and. timrev/=2) then
2279    MSG_BUG(sjoin('timrev:', itoa(timrev),' not allowed'))
2280  end if
2281 
2282  ABI_UNUSED(nbvw)
2283  !
2284  ! In the first loop calculate total number of transitions for this q-point
2285  ! as well min and max transition without taking into account distribution of bands.
2286  ! In the second iteration calculate min and Max transition for this processor.
2287  !
2288  spin_fact=half; if (nsppol==2) spin_fact=one
2289  my_max_rest=smallest_real; my_min_rest=greatest_real
2290     max_rest=smallest_real;    min_rest=greatest_real
2291 
2292  do iloop=1,2
2293    nt=0
2294    do ik_bz=1,Kmesh%nbz
2295      ik_ibz=Kmesh%tab(ik_bz)
2296      kmq(:)=Kmesh%bz(:,ik_bz)-qpoint(:)
2297 
2298      if (symchi==1) then
2299        if (Ltg_q%ibzq(ik_bz)/=1) cycle ! This point does not belong to the IBZ defined by the little group
2300      end if
2301 
2302      ! Find kp=k-q-G0 and also G0 where kp is in the first BZ
2303      if (.not.has_BZ_item(Kmesh,kmq,ikmq_bz,g0)) then ! Stop as the weight 1.0/nkbz is wrong.
2304        write(msg,'(4a,2(2a,3f12.6),2a)')ch10,&
2305 &        ' make_transitions : ERROR - ',ch10,&
2306 &        ' kp  = k-q-G0 not found in the BZ mesh',ch10,&
2307 &        ' k   = ',(Kmesh%bz(ii,ik_bz),ii=1,3),ch10,&
2308 &        ' k-q = ',(kmq(ii),ii=1,3),ch10,&
2309 &        ' weight in cchi0/cchi0q is wrong '
2310        MSG_ERROR(msg)
2311      end if
2312 
2313      ikmq_ibz=Kmesh%tab(ikmq_bz)
2314      do is=1,nsppol
2315        do ib1=1,nbnds
2316          do ib2=1,nbnds
2317 
2318            if (iloop==2) then
2319              if (bbp_ks_distrb(ib1,ib2,ik_bz,is)/=Wfd%my_rank) cycle
2320            end if
2321 
2322            if (timrev==2 .and. ib1<ib2) cycle ! Thanks to time-reversal we gain a factor ~2.
2323 
2324            delta_occ=spin_fact*(occ(ib1,ikmq_ibz,is)-occ(ib2,ik_ibz,is))
2325            delta_ene=gw_energy(ib1,ikmq_ibz,is)-gw_energy(ib2,ik_ibz,is)
2326 
2327            if (chi0alg==0)  then
2328              ! Adler-Wiser expression. Skip only if factor due to occupation number is smaller than TOL_DELTA_OCC
2329              if (abs(delta_occ) < abs(TOL_DELTA_OCC)) cycle
2330            else if (chi0alg==1) then
2331              ! Spectral method with time-reversal, only resonant transitions
2332              ! This has to changed to include spectral method without time-reversal
2333              if (delta_ene < -abs(TOL_DELTA_OCC) .or. abs(delta_occ) < abs(TOL_DELTA_OCC)) cycle
2334            end if
2335 
2336            ! We have a new transition
2337            nt=nt+1
2338 
2339            if (iloop==1) then
2340              max_rest=MAX(max_rest,zero,delta_ene)
2341              if (delta_ene>=-tol6) min_rest=MIN(min_rest,delta_ene)
2342            end if
2343            if (iloop==2) then
2344              my_max_rest=MAX(my_max_rest,zero,delta_ene)
2345              if (delta_ene>=-tol6) my_min_rest=MIN(my_min_rest,delta_ene)
2346            end if
2347 
2348          end do
2349        end do
2350      end do
2351    end do
2352    if (iloop==1) ntrans=nt
2353    if (iloop==2) my_ntrans=nt
2354  end do !iloop
2355 
2356  write(msg,'(2a,i9,2a,f8.3,3a,f8.3,a)')ch10,&
2357 &  ' Total number of transitions = ',ntrans,ch10,&
2358 &  ' min resonant     = ',min_rest*Ha_eV,' [eV] ',ch10,&
2359 &  ' Max resonant     = ',max_rest*Ha_eV,' [eV] '
2360  call wrtout(std_out,msg,'COLL')
2361 
2362  if (Wfd%nproc/=1) then
2363    write(msg,'(2a,i9,2a,f8.3,3a,f8.3,a)')ch10,&
2364 &    ' Total number of transitions for this processor= ',my_ntrans,ch10,&
2365 &    ' min resonant     = ',my_min_rest*Ha_eV,' [eV] ',ch10,&
2366 &    ' Max resonant     = ',my_max_rest*Ha_eV,' [eV] '
2367    call wrtout(std_out,msg,'PERS')
2368  end if
2369 
2370  DBG_EXIT("COLL")
2371 
2372 end subroutine make_transitions

m_chi0tk/mkrhotwg_sigma [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 mkrhotwg_sigma

FUNCTION

  Helper function used to calculate selected linear combination
  of the oscillator matrix elements in the case of non-collinear magnetism.

INPUTS

  ii=Index selecting the particolar combination of spin components.
  npw=Number of plane-waves in the oscillators.
  nspinor=Number of spinorial components.
  rhotwg(npw*nspinor**2)=OScillator matrix elements.

OUTPUT

  rhotwg_I(npw)=Required linear combination of the oscillator matrix elements.

PARENTS

CHILDREN

SOURCE

280 subroutine mkrhotwg_sigma(ii,nspinor,npw,rhotwg,rhotwg_I)
281 
282 
283 !This section has been created automatically by the script Abilint (TD).
284 !Do not modify the following lines by hand.
285 #undef ABI_FUNC
286 #define ABI_FUNC 'mkrhotwg_sigma'
287 !End of the abilint section
288 
289  implicit none
290 
291 !Arguments ------------------------------------
292 !scalars
293  integer,intent(in) :: ii,npw,nspinor
294 !arrays
295  complex(gwpc),intent(in) :: rhotwg(npw*nspinor**2)
296  complex(gwpc),intent(out) :: rhotwg_I(npw)
297 
298 ! *************************************************************************
299 
300  SELECT CASE (ii)
301  CASE (1)
302    ! $ M_0 = M_{\up,\up} + M_{\down,\down} $
303    rhotwg_I(:) = rhotwg(1:npw) + rhotwg(npw+1:2*npw)
304  CASE (2)
305    ! $ M_z = M_{\up,\up} - M_{\down,\down} $
306    rhotwg_I(:) = rhotwg(1:npw) - rhotwg(npw+1:2*npw)
307  CASE (3)
308    ! $ M_x = M_{\up,\down} + M_{\down,\up} $
309    rhotwg_I(:) = ( rhotwg(2*npw+1:3*npw) + rhotwg(3*npw+1:4*npw) )
310  CASE (4)
311    ! $ M_y = i * (M_{\up,\down} -M_{\down,\up}) $
312    rhotwg_I(:) = (rhotwg(2*npw+1:3*npw) - rhotwg(3*npw+1:4*npw) )*j_gw
313  CASE DEFAULT
314    MSG_BUG(sjoin('Wrong ii value:', itoa(ii)))
315  END SELECT
316 
317 end subroutine mkrhotwg_sigma

m_chi0tk/output_chi0sumrule [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 output_chi0sumrule

FUNCTION

  Calculate and output the value of the sum rule for
  the non-interacting polarizability chi0

INPUTS

OUTPUT

  (for writing routines, no output)
  otherwise, should be described

PARENTS

      screening

CHILDREN

SOURCE

2039 subroutine output_chi0sumrule(qeq0,iq,npwe,omegaplasma,chi0sumrule,epsm1_w0,vc_sqrt)
2040 
2041 
2042 !This section has been created automatically by the script Abilint (TD).
2043 !Do not modify the following lines by hand.
2044 #undef ABI_FUNC
2045 #define ABI_FUNC 'output_chi0sumrule'
2046 !End of the abilint section
2047 
2048  implicit none
2049 
2050 !Arguments ------------------------------------
2051 !scalars
2052  integer,intent(in) :: iq,npwe
2053  real(dp),intent(in) :: omegaplasma
2054  logical,intent(in) :: qeq0
2055 !arrays
2056  real(dp),intent(inout) :: chi0sumrule(npwe)
2057  complex(gwpc),intent(in) :: epsm1_w0(npwe,npwe),vc_sqrt(npwe)
2058 
2059 !Local variables ------------------------------
2060 !scalars
2061  integer :: ig,igstart
2062  real(dp) :: average,norm
2063  character(len=500) :: msg
2064 
2065 !************************************************************************
2066 
2067  igstart=1; if (qeq0) igstart=2
2068  !
2069  ! The sumrule reads:
2070  ! $ \int d\omega \omega v * Im[ \chi_0(\omega) ] = \pi/2 * w_p^2 $.
2071  chi0sumrule(igstart:npwe) = chi0sumrule(igstart:npwe) * vc_sqrt(igstart:npwe)**2
2072  !
2073  ! Calculate a weighted average of the fulfilment of the sumrule on epsilon
2074  ! The weight is given according to the significance of each q+G in the
2075  ! subsequent GW calculation: It is proportional to v * (epsm1 -1 )
2076  average = zero; norm = zero
2077  do ig=igstart,npwe
2078    average = average + chi0sumrule(ig) * real( vc_sqrt(ig)**2 * (epsm1_w0(ig,ig) - 1.0_dp ) )
2079    norm    = norm    +                   real( vc_sqrt(ig)**2 * (epsm1_w0(ig,ig) - 1.0_dp ) )
2080    !average = average + chi0sumrule(ig) * real(  (epsm1_w0(ig,ig) - 1.0_dp ) )
2081    !norm    = norm    +                   real(  (epsm1_w0(ig,ig) - 1.0_dp ) )
2082    !write(203,'(i4,8(2x,e12.6))') ig,1.0_dp/vc_sqrt(ig),chi0sumrule(ig)/ (0.5d0*omegaplasma**2*pi)
2083  end do
2084 
2085  if (abs(norm)>tol8) then
2086    write(msg,'(1x,a,i4,a,f10.2,2x,a)')&
2087     ' Average fulfillment of the sum rule on Im[epsilon] for q-point ',&
2088     iq,' :',average/norm/(0.5_dp*omegaplasma**2*pi)*100.0_dp,'[%]'
2089    call wrtout(std_out,msg,'COLL'); call wrtout(ab_out, msg,'COLL')
2090  end if
2091 
2092 end subroutine output_chi0sumrule

m_chi0tk/setup_spectral [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

  setup_spectral

FUNCTION

 Calculation of \chi_o based on the spectral method as proposed in PRB 74, 035101 (2006) [[cite:Shishkin2006]]
 and PRB 61, 7172 (2000) [[cite:Miyake2000]].
 Setup of the real frequency mesh for $\Im\chi_o$ and of the frequency-dependent weights for
 Hilbert transform. Note that CPU time does not depend dramatically on nomegasf unlike memory.
 spmeth defines the approximant for the delta function:
  ==1 : use Triangular approximant (Kresse method)
  ==2 : use Gaussian method, requiring smearing (Miyake method)

INPUTS

 nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$
 nomega=number of frequencies in $\chi0(q,\omega)$.
 max_rest,min_res=max and min resonant transition energy (for this q-point)
 my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor
 method=integer flag defining the type of frequency mesh used for $\Im chi0$
  | 0 for a linear mesh
  | 1 for a mesh densified around omegaplasma
 omegaplasma=frequency around which the mesh is densifies (usually Drude plasma frequency)
  used only in case of method==1
 zcut=small imaginary shift to avoid pole in chi0

OUTPUT

  kkweight(nomegasf,nomega)=Frequency dependent weight for Hilber transform.
  omegasf(nomegasf+1)=frequencies for imaginary part.

PARENTS

      cchi0,cchi0q0

CHILDREN

SOURCE

1567 subroutine setup_spectral(nomega,omega,nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,&
1568 &  method,zcut,omegaplasma,my_wl,my_wr,kkweight)
1569 
1570 
1571 !This section has been created automatically by the script Abilint (TD).
1572 !Do not modify the following lines by hand.
1573 #undef ABI_FUNC
1574 #define ABI_FUNC 'setup_spectral'
1575 !End of the abilint section
1576 
1577  implicit none
1578 
1579 !Arguments ------------------------------------
1580 !scalars
1581  integer,intent(in) :: method,nomega,nomegasf
1582  integer,intent(out) :: my_wl,my_wr
1583  real(dp),intent(in) :: max_rest,min_rest,omegaplasma,zcut
1584  real(dp),intent(in) :: my_max_rest,my_min_rest
1585 !arrays
1586  real(dp),intent(out) :: omegasf(nomegasf)
1587  complex(dpc),intent(in) :: omega(nomega)
1588  complex(dpc),intent(out) :: kkweight(nomegasf,nomega)
1589 
1590 !Local variables-------------------------------
1591 !scalars
1592  integer :: io,ii
1593  real(dp) :: nu_min,nu_max,nu1,nu2,dd,domegasf,wp,deltat
1594  character(len=500) :: msg
1595 !arrays
1596  integer,allocatable :: insort(:)
1597 !************************************************************************
1598 
1599  ! The mesh must enclose the entire range of transitions.
1600  dd=(max_rest-min_rest)/(nomegasf-1)
1601  domegasf=(max_rest-min_rest+2*dd)/(nomegasf-1)
1602 
1603  write(msg,'(4a,f8.3,3a,i5,2a,f8.5,a)')ch10,&
1604 &  ' === Info on the real frequency mesh for spectral method === ',ch10,&
1605 &  '  maximum frequency = ',max_rest*Ha_eV,' [eV]',ch10,&
1606 &  '  nomegasf = ',nomegasf,ch10,&
1607 &  '  domegasf = ',domegasf*Ha_eV,' [eV]'
1608  call wrtout(std_out,msg,'COLL')
1609 
1610  if (min_rest<tol6) then
1611    MSG_WARNING("System seems to be metallic")
1612  end if
1613 
1614  ! ======================================================
1615  ! === Setup of the w-mesh for the spectral function ====
1616  ! ======================================================
1617  SELECT CASE (method)
1618  CASE (0)
1619    ! Linear mesh.
1620    call wrtout(std_out,' Using linear mesh for Im chi0','COLL')
1621    do io=1,nomegasf
1622      omegasf(io)=(io-1)*domegasf+min_rest-dd
1623    end do
1624 
1625  CASE (1)
1626    ! Non-homogeneous mesh densified around omega_plasma, do not improve results ===
1627    ! WARNING_ this part has to be checked since I modified omegasf
1628    write(msg,'(a,f7.4,a)')' Using mesh densified around ',omegaplasma*Ha_eV,' [eV] '
1629    call wrtout(std_out,msg,'COLL')
1630    wp=omegaplasma ; deltat=max_rest-min_rest
1631    nu_min=zero
1632    if (deltat<wp ) then
1633      nu_max = wp/sqrt2 *   ATAN(sqrt2*deltat*wp/(-deltat**2+wp**2))
1634    else
1635      nu_max = wp/sqrt2 * ( ATAN(sqrt2*deltat*wp/(-deltat**2+wp**2)) + pi)
1636    end if
1637    domegasf=(nu_max-nu_min)/(nomegasf+1)
1638    !write(std_out,*)  -(wp/sqrt2) * atan(sqrt2*deltat*wp/(deltat**2-wp**2))
1639    omegasf(1)=zero ; omegasf(nomegasf+1)=deltat
1640    ii=0
1641    do io=2,nomegasf
1642      nu1=domegasf*(io-1) ; nu2=TAN(-sqrt2*nu1/wp)
1643      if (nu2<0) then
1644        omegasf(io) = wp * (one - SQRT(1+2*nu2**2))/(sqrt2*nu2)
1645      else
1646        omegasf(io) = wp * (one + SQRT(1+2*nu2**2))/(sqrt2*nu2)
1647      end if
1648      if (omegasf(io)> deltat ) then
1649        omegasf(io)= deltat-0.1*ii
1650        ii=ii+1
1651      end if
1652      ! write(102,'(i4,2x,3(f9.4,2x))')io,nu1,nu2,ep%omegasf(io)*Ha_eV
1653    end do
1654 
1655    ! Reorder frequencies in ascending order
1656    ABI_MALLOC(insort,(nomegasf+1))
1657    insort(:)=(/ (io,io=1,nomegasf+1) /)
1658    call sort_dp(nomegasf+1,omegasf,insort,tol14)
1659    ABI_FREE(insort)
1660 
1661  CASE DEFAULT
1662    MSG_BUG(sjoin('Wrong value for method:', itoa(method)))
1663  END SELECT
1664  !write(std_out,*)omegasf(1)*Ha_eV,omegasf(nomegasf)*Ha_eV
1665 
1666  ! Find min and max index in omegasf treated by this processor.
1667  my_wr=-999
1668  do io=1,nomegasf
1669    if (omegasf(io)>my_max_rest) then
1670      my_wr=io; EXIT
1671    end if
1672  end do
1673  if (my_wr==nomegasf+2) my_wr=nomegasf+1
1674  my_wl=-999
1675  do io=nomegasf,1,-1
1676    if (omegasf(io)< my_min_rest) then ! Check metals
1677      my_wl=io; EXIT
1678    end if
1679  end do
1680 
1681  write(msg,'(a,2(1x,i0))')' my_wl and my_wr:',my_wl,my_wr
1682  call wrtout(std_out,msg,'PERS')
1683 
1684  if (my_wl==-999 .or. my_wr==-999) then
1685    write(msg,'(a,2i6)')' wrong value in my_wl and/or my_wr ',my_wl,my_wr
1686    MSG_ERROR(msg)
1687  end if
1688 
1689  ! Calculate weights for Hilbert transform.
1690  call calc_kkweight(nomega,omega,nomegasf,omegasf,zcut,max_rest,kkweight)
1691 
1692 end subroutine setup_spectral

m_chi0tk/symmetrize_afm_chi0tk [ Functions ]

[ Top ] [ m_chi0tk ] [ Functions ]

NAME

 symmetrize_afm_chi0

FUNCTION

  Reconstruct the (down, down) component of the irreducible polarizability
  starting from the (up,up) element in case of systems with AFM symmetries
  (i.e nspden==2 and nsppol=1). Return the trace (up,up)+(down,down) of the
  matrix as required by GW calculations.

INPUTS

  Cryst<crystal_t>= Information on symmetries and unit cell.
  Gsph<gsphere_t>= The G-sphere used to descrive chi0.
  npwe=Number of G-vectors in chi0.
  nomega=number of frequencies.
  Ltg_q<littlegroup_t>=Structure with useful table describing the little group of the q-point.

SIDE EFFECTS

 chi0(npwe,npwe,nomega)= In input the up-up component, in output the trace of chi0.
   The value of matrix elements that should be zero due to AFM symmetry properties are
   forced to be zero (see NOTES below).
 [chi0_lwing(npwe,nomega,3)] = Lower wings, symmetrized in output.
 [chi0_uwing(npwe,nomega,3)] = Upper wings, symmetrized in output.
 [chi0_head(3,3,nomega)    ] = Head of chi0, symmetrized  in output.

NOTES

 In the case of magnetic group Shubnikov type III:
   For each set of paired FM-AFM symmetries, the down-down component of
   a generic response function in reciprocal space can be obtained according to:

    chi^{down,down}_{G1,G2}(q) = chi^{up up}_{G1,G2}(q) e^{iS(G1-G2).(tnonsFM - tnonsAFM)}

   where S is the rotational part common to the FM-AFM pair, tnonsFM and tnonsAFM
   are the fractional translations associated to the ferromagnetic and antiferromagnetic symmetry, respectively.
   Note that, if for a given G1-G2 pair, the phase e^{iS(G1-G2).(tnonsFM - tnonsAFM) depends
   on the FM-AFM symmetry pair, then the corresponding matrix element of chi0 must be zero.
   Actually this is manually enforced in the code because this property might not be
   perfectly satisfied due to round-off errors.

 In the case of magnetic group Shubnikov type III:
   Only the AFM symmetries that preserve the external q-point (with or without time-reversal)
   are used to get the (down, down) component using the fact that:

    chi^{down,down}_{G1,G2}(Sq) = chi^{up up}_{S^{-1}G1,S^{-1}G2}(q) e^{i(G2-G1).tnons_S }

   Actually we perform an average over subset of the little group of q with AFM character in
   order to reduce as much as possible errors due to round off errors. In brief we evaluate:

    1/N_{Ltq} \sum_{S\in Ltg AFM} chi^{up up}_{S^{-1}G1,S^{-1}G2}(q) e^{i(G2-G1).tnons_S }

   where N_{Ltg} is the number of AFM operation in the little group (time reversal included)

TODO

  It is possible to symmetrize chi0 without any the extra allocation for afm_mat.
  More CPU demanding but safer in case of a large chi0 matrix. One might loop over G1 and G2 shells ...

PARENTS

      cchi0,cchi0q0,cchi0q0_intraband

CHILDREN

SOURCE

385 subroutine symmetrize_afm_chi0(Cryst,Gsph,Ltg_q,npwe,nomega,chi0,chi0_head,chi0_lwing,chi0_uwing)
386 
387 
388 !This section has been created automatically by the script Abilint (TD).
389 !Do not modify the following lines by hand.
390 #undef ABI_FUNC
391 #define ABI_FUNC 'symmetrize_afm_chi0'
392 !End of the abilint section
393 
394  implicit none
395 
396 !Arguments ------------------------------------
397 !scalars
398  integer,intent(in) :: npwe,nomega
399  type(gsphere_t),intent(in) :: Gsph
400  type(crystal_t),intent(in) :: Cryst
401  type(littlegroup_t),intent(in) :: Ltg_q
402 !arrays
403  complex(gwpc),intent(inout) :: chi0(npwe,npwe,nomega)
404  complex(dpc),optional,intent(inout) :: chi0_lwing(npwe,nomega,3)
405  complex(dpc),optional,intent(inout) :: chi0_uwing(npwe,nomega,3)
406  complex(dpc),optional,intent(inout) :: chi0_head(3,3,nomega)
407 
408 !Local variables ------------------------------
409 !scalars
410  integer :: io,ig1,ig2,isymf,isyma,isym,ipair,k0g,kg,npairs,nonzero
411  integer :: iSmg1,iSmg2,itim,shubnikov,ntest
412  complex(gwpc) :: phase,phase_old,sumchi,ctmp
413  logical :: found
414  !character(len=500) :: msg
415 !arrays
416  integer :: rotfm(3,3),rotafm(3,3),pairs2sym(2,Cryst%nsym/2)
417  real(dp) :: tfm(3),tafm(3)
418  complex(gwpc),allocatable :: afm_mat(:),chi0_afm(:,:)
419 
420 !************************************************************************
421 
422  ABI_CHECK(ANY(Cryst%symafm==-1),'Not magnetic space group')
423  !
424  ! ==== Find shubnikov type ====
425  ! TODO This info should be stored in Cryst%
426  shubnikov=4; npairs=0
427 
428  do isymf=1,Cryst%nsym
429    if (Cryst%symafm(isymf)==-1) CYCLE
430    rotfm = Cryst%symrec(:,:,isymf)
431    tfm = Cryst%tnons(:,isymf)
432    found = .FALSE.
433 
434    do isyma=1,Cryst%nsym
435      if (Cryst%symafm(isyma)==1) CYCLE
436      rotafm = Cryst%symrec(:,:,isyma)
437 
438      if (ALL(rotfm==rotafm)) then
439        found=.TRUE.
440        tafm = Cryst%tnons(:,isyma)
441        npairs=npairs+1
442        ABI_CHECK(npairs<=Cryst%nsym/2,'Wrong AFM group')
443        pairs2sym(1,npairs)=isymf
444        pairs2sym(2,npairs)=isyma
445      end if
446    end do !isyma
447 
448    if (.not.found) then
449     shubnikov=3; EXIT !isymf
450    end if
451  end do !isymf
452 
453  select case (shubnikov)
454 
455  case (4)
456    call wrtout(std_out,' Found Magnetic group Shubnikov type IV','COLL')
457    ABI_CHECK(npairs==Cryst%nsym/2,'Wrong AFM space group')
458 
459    ABI_MALLOC(afm_mat,(npwe*(npwe+1)/2))
460 
461 ! jmb
462    phase_old=zero
463 
464    do ig2=1,npwe
465      k0g=ig2*(ig2-1)/2
466      do ig1=1,ig2
467        kg=k0g+ig1
468        nonzero=1
469 
470        do ipair=1,Cryst%nsym/2
471          isymf = pairs2sym(1,ipair)
472          isyma = pairs2sym(2,ipair)
473          phase = ( Gsph%phmSGt(ig1,isymf)*CONJG(Gsph%phmSGt(ig1,isyma)) ) * &
474 &                ( Gsph%phmSGt(ig2,isymf)*CONJG(Gsph%phmSGt(ig2,isyma)) )
475          if (ipair>1 .and. (ABS(phase_old-phase) > tol6)) then
476            nonzero=0; EXIT
477          end if
478          phase_old=phase
479        end do !ipair
480 
481        afm_mat(kg)=nonzero*(cone_gw + phase)
482      end do !ig1
483    end do !ig2
484    !
485    ! =======================================================================
486    ! ==== Symmetrize chi0 constructing chi0^{\up\up} + chi0{\down\down} ====
487    ! =======================================================================
488    !
489    !  head^{\down\down} = head^{\up\up}
490    if (PRESENT(chi0_head)) chi0_head = two * chi0_head
491    !
492    ! w^{\down\down}_{0 G'} =  e^{-iSG'.(tFM-tAFM)} w^{\up\up}_{0 G'}.
493    ! w^{\down\down}_{G 0 } =  e^{+iSG .(tFM-tAFM)} w^{\up\up}_{G 0 }.
494    if (PRESENT(chi0_uwing)) then
495      do io=1,nomega
496        do ig2=1,npwe
497          k0g=ig2*(ig2-1)/2
498          kg=k0g+1
499          chi0_uwing(ig2,io,:)=afm_mat(kg)*chi0_uwing(ig2,io,:)
500        end do
501      end do
502    end if
503 
504    if (PRESENT(chi0_lwing)) then
505      do io=1,nomega
506        do ig1=1,npwe
507          k0g=ig1*(ig1-1)/2
508          kg=k0g+1
509          chi0_lwing(ig1,io,:)=CONJG(afm_mat(kg))*chi0_lwing(ig1,io,:)
510        end do
511      end do
512    end if
513 
514    do io=1,nomega
515      ! Take care of diagonal.
516      do ig1=1,npwe
517        chi0(ig1,ig1,io)=two*chi0(ig1,ig1,io)
518      end do
519 
520      ! Upper and lower triangle are treated differently:
521      ! We took advantage of the fact the afm_mat is hermitian to reduce memory.
522      do ig2=2,npwe
523        k0g=ig2*(ig2-1)/2
524        do ig1=1,ig2-1
525          kg=k0g+ig1
526          chi0(ig1,ig2,io)=afm_mat(kg)*chi0(ig1,ig2,io)
527        end do
528      end do
529 
530      do ig1=2,npwe
531        k0g=ig1*(ig1-1)/2
532        do ig2=1,ig1-1
533          kg=k0g+ig2
534          chi0(ig1,ig2,io)=CONJG(afm_mat(kg))*chi0(ig1,ig2,io)
535        end do
536      end do
537 
538    end do !io
539 
540    ABI_FREE(afm_mat)
541 
542  case (3)
543    call wrtout(std_out,' Found Magnetic group Shubnikov type III',"COLL")
544    MSG_ERROR('Shubnikov type III not implemented')
545 
546    ntest=0
547    do itim=1,ltg_q%timrev
548      do isym=1,ltg_q%nsym_sg
549        ! use only afm sym preserving q with and without time-reversal
550        if ( cryst%symafm(isym)==-1 .and. ltg_q%preserve(itim,isym)==1 ) ntest=ntest+1
551      end do
552    end do
553 
554    if (ntest==0) then
555        MSG_WARNING("no symmetry can be used!")
556    end if
557    !RETURN
558    ABI_MALLOC(chi0_afm,(npwe,npwe))
559 
560    do io=1,nomega
561 
562      do ig2=1,npwe
563        do ig1=1,npwe
564          sumchi=czero_gw
565 
566          do itim=1,ltg_q%timrev
567            do isym=1,ltg_q%nsym_sg
568              ! use only afm sym preserving q with and without time-reversal
569              if ( cryst%symafm(isym)==-1 .and. ltg_q%preserve(itim,isym)==1 ) then
570                phase =  Gsph%phmGt(ig1,isym)*CONJG(Gsph%phmGt(ig2,isym))
571                iSmg1=Gsph%rottbm1(ig1,itim,isym)
572                iSmg2=Gsph%rottbm1(ig2,itim,isym)
573                ctmp=chi0(iSmg1,iSmg2,io)*phase !; if (itim==2) ctmp=CONJG(ctmp) !check this
574                sumchi=sumchi+ctmp !chi0(iSmg1,iSmg2,io)*phase
575              end if
576            end do ! isym
577          end do !itim
578 
579          chi0_afm(ig1,ig2)=sumchi/Ltg_q%nsym_ltg  !has to be changed in case of time-reversal
580        end do !ig1
581      end do !ig2
582 
583      ! We want chi_{up,up} +chi_{dwn,dwn}.
584      chi0(:,:,io)=chi0(:,:,io)+chi0_afm(:,:)
585    end do !iomega
586 
587    ABI_FREE(chi0_afm)
588 
589  case default
590    MSG_BUG(sjoin('Wrong value for shubnikov= ', itoa(shubnikov)))
591  end select
592 
593 end subroutine symmetrize_afm_chi0