TABLE OF CONTENTS


ABINIT/calc_sigx_me [ Functions ]

[ Top ] [ Functions ]

NAME

 calc_sigx_me

FUNCTION

 Calculate diagonal and off-diagonal matrix elements of the exchange part of the self-energy operator.

INPUTS

 sigmak_ibz=Index of the k-point in the IBZ.
 minbnd, maxbnd= min and Max band index for GW correction (for this k-point)
 Gsph_x<gsphere_t>= Info on the G-sphere used for Sigma_x
    %nsym=number of symmetry operations
    %rottb(ng,timrev,nsym)=index of (IS) G where I is the identity or the inversion
      operation and G is one of the ng vectors in reciprocal space
    %timrev=2 if time-reversal symmetry is used, 1 otherwise
    %gvec(3,ng)=integer coordinates of each plane wave in reciprocal space
 ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
 Ltg_k datatype containing information on the little group
 Kmesh <kmesh_t>
    %nbz=Number of points in the BZ
    %nibz=Number of points in IBZ
    %kibz(3,nibz)=k-point coordinates, irreducible Brillouin zone
    %kbz(3,nbz)=k-point coordinates, full Brillouin zone
    %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
    %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered
    %ktabp(nbz)= phase factor associated to tnons
 gwx_ngfft(18)=Information about 3D FFT for the oscillator strengths, see ~abinit/doc/variables/vargs.htm#ngfft
 gwx_nfftot=number of points of the FFT grid for GW wavefunctions
 Vcp <vcoul_t datatype> containing information on the cutoff technique
    %vc_sqrt(npwx,nqibz)= square-root of the coulombian potential for q-points in the IBZ
 Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
 Pawang <type(pawang_type)>=paw angular mesh and related data
 Psps <type(pseudopotential_type)>=variables related to pseudopotentials
    %usepaw=1 for PAW, 0 for NC pseudopotentials.
 Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
    %ibz=q points where $\tilde\epsilon^{-1}$ has been computed
    %bz(3,nqbz)=coordinates of all q-points in BZ
 Sigp <sigparams_t> (see the definition of this structured datatype)
 Cryst<crystal_t>=Info on unit cell and symmetries
    %natom=number of atoms in unit cell
    %ucvol=unit cell volume
    %nsym=number of symmetry operations
    %typat(natom)=type of each atom
  much slower but it requires less memory
 QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
  eig(Sigp%nbnds,Kmesh%nibz,%nsppol)=KS or QP energies for k-points, bands and spin
  occ(Sigp%nbnds,Kmesh%nibz,nsppol)=occupation numbers, for each k point in IBZ, each band and spin
  Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
  allQP_sym(%nkibz, %nsppol)<esymm_t>=Datatype collecting data on the irreducible representations of the
    little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
 prtvol=Flags governing verbosity level.

OUTPUT

  Sr%x_mat(minbnd:maxbnd,minbnd:maxbnd,%nsppol*Sigp%nsig_ab)=Matrix elements of Sigma_x.

NOTES

  1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) [[cite:Gigy1986]] is included.

  2) On the symmetrization of Sigma matrix elements
     If  Sk = k+G0 then  M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G}   (k,q)
     If -Sk = k+G0 then  M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)

 Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out
 that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a
 similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase
 has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ.

  3) the unitary transformation relating wavefunctions
     at symmetric k-points should be taken into account during the symmetrization
     of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
     it is possible to make an invariant by just including all the degenerate states and
     averaging the final results over the degenerate subset. Here we divide the states
     where the QP energies are required into complexes. Note however that this approach is not
     based on group theory, and it might lead to spurious results in case of accidental degeneracies.

PARENTS

      sigma

CHILDREN

      cwtime,esymm_symmetrize_mels,findqg0,get_bz_item,gsph_fft_tabs
      hermitianize,littlegroup_print,paw_cross_rho_tw_g,paw_rho_tw_g
      paw_symcprj,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawmknhat_psipsi
      pawpwij_free,pawpwij_init,rho_tw_g,rotate_fft_mesh,sigma_distribute_bks
      timab,wfd_change_ngfft,wfd_get_cprj,wfd_get_many_ur,wfd_get_ur
      wfd_paw_get_aeur,wrtout,xmpi_sum

SOURCE

128 subroutine calc_sigx_me(sigmak_ibz,ikcalc,minbnd,maxbnd,Cryst,QP_BSt,Sigp,Sr,Gsph_x,Vcp,Kmesh,Qmesh,&
129 & Ltg_k,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,allQP_sym,gwx_ngfft,ngfftf,&
130 & prtvol,pawcross)
131 
132  use defs_basis
133  use defs_datatypes
134  use m_abicore
135  use m_gwdefs
136  use m_xmpi
137  use m_defs_ptgroups
138  use m_errors
139  use m_time
140 
141  use m_hide_blas,     only : xdotc, xgemv
142  use m_numeric_tools, only : hermitianize
143  use m_geometry,      only : normv
144  use m_crystal,       only : crystal_t
145  use m_fft_mesh,      only : rotate_FFT_mesh, cigfft
146  use m_bz_mesh,       only : kmesh_t, get_BZ_item, findqg0, littlegroup_t, littlegroup_print
147  use m_gsphere,       only : gsphere_t, gsph_fft_tabs
148  use m_vcoul,         only : vcoul_t
149  use m_pawpwij,       only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
150  use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t
151  use m_pawang,        only : pawang_type
152  use m_pawtab,        only : pawtab_type
153  use m_pawfgrtab,     only : pawfgrtab_type
154  use m_pawcprj,       only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
155  use m_paw_nhat,      only : pawmknhat_psipsi
156  use m_paw_sym,       only : paw_symcprj
157  use m_wfd,           only : wfd_t, wfd_get_ur, wfd_get_cprj, wfd_change_ngfft, wfd_paw_get_aeur, wfd_get_many_ur,&
158 &                            wfd_sym_ur
159  use m_sigma,         only : sigma_t, sigma_distribute_bks
160  use m_oscillators,   only : rho_tw_g
161  use m_esymm,         only : esymm_t, esymm_symmetrize_mels, esymm_failed
162  use m_ptgroups,      only : sum_irreps
163 
164 !This section has been created automatically by the script Abilint (TD).
165 !Do not modify the following lines by hand.
166 #undef ABI_FUNC
167 #define ABI_FUNC 'calc_sigx_me'
168 !End of the abilint section
169 
170  implicit none
171 
172 !Arguments ------------------------------------
173 !scalars
174  integer,intent(in) :: sigmak_ibz,ikcalc,prtvol,minbnd,maxbnd,pawcross
175  type(crystal_t),intent(in) :: Cryst
176  type(ebands_t),target,intent(in) :: QP_BSt
177  type(kmesh_t),intent(in) :: Kmesh,Qmesh
178  type(vcoul_t),intent(in) :: Vcp
179  type(gsphere_t),intent(in) :: Gsph_x
180  type(littlegroup_t),intent(in) :: Ltg_k
181  type(Pseudopotential_type),intent(in) :: Psps
182  type(sigparams_t),target,intent(in) :: Sigp
183  type(sigma_t),intent(inout) :: Sr
184  type(pawang_type),intent(in) :: Pawang
185  type(wfd_t),target,intent(inout) :: Wfd,Wfdf
186 !arrays
187  integer,intent(in) :: gwx_ngfft(18),ngfftf(18)
188  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
189  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
190  type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol)
191  type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
192  type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom*Psps%usepaw)
193 
194 !Local variables ------------------------------
195 !scalars
196  integer,parameter :: tim_fourdp=2,ndat1=1
197  integer,parameter :: use_pawnhat=0,ider0=0
198  integer :: gwcalctyp,izero,iab,ib_sum,ib,ib1,ib2,ierr,ig,ig_rot,ii,iik,itim_q,i2
199  integer :: ik_bz,ik_ibz,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx
200  integer :: jik,jk_bz,jk_ibz,kb,nspinor,nsppol,ifft
201  integer :: nq_summed,ibsp,dimcprj_gw,dim_rtwg
202  integer :: spad,spadx1,spadx2,irow,npw_k,ndegs,wtqm,wtqp,my_nbks
203  integer :: isym_kgw,isym_ki,gwx_mgfft,use_padfft,use_padfftf,gwx_fftalga,gwx_fftalgb
204  integer :: gwx_nfftot,nfftf,mgfftf,nhat12_grdim,npwx
205  real(dp) :: cpu_time,wall_time,gflops
206  real(dp) :: fact_sp,theta_mu_minus_esum,tol_empty
207  complex(dpc) :: ctmp,ph_mkgwt,ph_mkt
208  complex(gwpc) :: gwpc_sigxme
209  logical :: iscompatibleFFT,q_is_gamma
210  character(len=500) :: msg
211 !arrays
212  integer :: g0(3),spinor_padx(2,4)
213  integer,allocatable :: igfftxg0(:),igfftfxg0(:)
214  integer,allocatable :: degtab(:,:,:)
215  integer,allocatable :: gwx_gfft(:,:),gwx_gbound(:,:),gboundf(:,:)
216  integer,allocatable ::  ktabr(:,:),irottb(:,:),ktabrf(:,:)
217  integer,allocatable :: proc_distrb(:,:,:)
218  real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),qbz(3),q0(3),tsec(2)
219  real(dp) :: spinrot_kbz(4),spinrot_kgw(4)
220  real(dp),ABI_CONTIGUOUS pointer :: qp_ene(:,:,:),qp_occ(:,:,:)
221  real(dp),allocatable :: nhat12(:,:,:),grnhat12(:,:,:,:)
222  complex(gwpc),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:)
223  complex(gwpc),allocatable :: rhotwg_ki(:,:)
224  complex(gwpc),allocatable :: wfr_bdgw(:,:),ur_ibz(:)
225  complex(gwpc),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:)
226  complex(gwpc),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:)
227  complex(dpc),allocatable :: sigxme_tmp(:,:,:),sym_sigx(:,:,:),sigx(:,:,:,:)
228  complex(gwpc),ABI_CONTIGUOUS pointer :: cg_jb(:),cg_sum(:)
229  logical :: can_symmetrize(Wfd%nsppol)
230  logical,allocatable :: bks_mask(:,:,:)
231  type(esymm_t),pointer :: QP_sym(:)
232  type(sigijtab_t),pointer :: Sigxij_tab(:)
233  type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
234  type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
235 
236 !************************************************************************
237 
238  DBG_ENTER("COLL")
239 
240  call timab(430,1,tsec) ! csigme (SigX)
241  call cwtime(cpu_time,wall_time,gflops,"start")
242 
243  ! Initialize some values.
244  gwcalctyp=Sigp%gwcalctyp
245  nspinor = Wfd%nspinor; nsppol = Wfd%nsppol; npwx = sigp%npwx
246  dim_rtwg = 1; if (nspinor == 2) dim_rtwg = 2
247  spinor_padx = RESHAPE([0, 0, npwx, npwx, 0, npwx, npwx, 0], [2, 4])
248  ABI_CHECK(Sigp%npwx==Gsph_x%ng,"")
249 
250  qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ
251 
252  ! Exctract the symmetries of the bands for this k-point
253  QP_sym => allQP_sym(sigmak_ibz,1:nsppol)
254 
255  ib1=minbnd; ib2=maxbnd
256 
257  ! Index of the GW point in the BZ array, its image in IBZ and time-reversal
258  jk_bz = Sigp%kptgw2bz(ikcalc)
259  call get_BZ_item(Kmesh,jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
260 
261  ! TODO: the new version based of get_uug won't support kptgw vectors that are not in
262  ! the IBZ since one should perform the rotation before entering the band loop
263  ! In the old version, the rotation was done in rho_tw_g
264  !ABI_CHECK(jik==1,"jik!=1")
265  !ABI_CHECK(isym_kgw==1,"isym_kgw!=1")
266  !ABI_CHECK((ABS(ph_mkgwt - cone) < tol12),"ph_mkgwt!")
267  !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
268  spinrot_kgw(:)=Cryst%spinrot(:,isym_kgw)
269  write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,&
270 &  ' Calculating <nk|Sigma_x|nk> at k= ',kgw,ch10,&
271 &  ' bands from ',ib1,' to ',ib2,ch10
272  call wrtout(std_out,msg,'COLL')
273 
274  if (ANY(gwx_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call wfd_change_ngfft(Wfd,Cryst,Psps,gwx_ngfft)
275  gwx_mgfft = MAXVAL(gwx_ngfft(1:3))
276  gwx_fftalga = gwx_ngfft(7)/100; gwx_fftalgb = MOD(gwx_ngfft(7),100)/10
277 
278  if (pawcross==1) mgfftf = MAXVAL(ngfftf(1:3))
279 
280  can_symmetrize = .FALSE.
281  if (Sigp%symsigma>0) then
282    can_symmetrize = .TRUE.
283    if (gwcalctyp >= 20) then
284      do spin=1,Wfd%nsppol
285        can_symmetrize(spin) = .not. esymm_failed(QP_sym(spin))
286        if (.not.can_symmetrize(spin)) then
287          write(msg,'(a,i0,4a)')&
288 &          "Symmetrization cannot be performed for spin: ",spin,ch10,&
289 &          "band classification encountered the following problem: ",ch10,&
290 &          TRIM(QP_sym(spin)%err_msg)
291          MSG_WARNING(msg)
292        end if
293      end do
294    end if
295    if (nspinor == 2) then
296      MSG_WARNING('Symmetrization with nspinor=2 not implemented')
297    end if
298  end if
299 
300  ABI_MALLOC(rhotwg_ki, (npwx*nspinor, minbnd:maxbnd))
301  rhotwg_ki=czero_gw
302  ABI_MALLOC(rhotwg, (npwx*nspinor))
303  ABI_MALLOC(rhotwgp, (npwx*nspinor))
304  ABI_MALLOC(vc_sqrt_qbz, (npwx))
305 
306  ! Normalization of theta_mu_minus_esum
307  ! If nsppol==2, qp_occ $\in [0,1]$
308  SELECT CASE (nsppol)
309  CASE (1)
310    fact_sp=half; tol_empty=0.01   ! below this value the state is assumed empty
311    if (Sigp%nspinor==2) then
312     fact_sp=one; tol_empty=0.005  ! below this value the state is assumed empty
313    end if
314  CASE (2)
315    fact_sp=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic
316  CASE DEFAULT                    ! spin unpolarized system is treated using nsppol==2
317    MSG_BUG('Wrong nsppol')
318  END SELECT
319 
320  ! Table for \Sigmax_ij matrix elements.
321  Sigxij_tab => Sigp%Sigxij_tab(ikcalc,1:nsppol)
322 
323  ! Remove empty states from the list of states that will be distributed.
324  ABI_MALLOC(bks_mask,(Wfd%mband,Kmesh%nbz,nsppol))
325  bks_mask=.FALSE.
326  do spin=1,nsppol
327    do ik_bz=1,Kmesh%nbz
328      ik_ibz = Kmesh%tab(ik_bz)
329      do ib_sum=1,Sigp%nbnds
330        bks_mask(ib_sum,ik_bz,spin) = qp_occ(ib_sum,ik_ibz,spin) >= tol_empty
331      end do
332    end do
333  end do
334 
335  ! Distribute tasks.
336  ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,nsppol))
337  call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,proc_distrb,bks_mask=bks_mask)
338  ABI_FREE(bks_mask)
339  write(msg,'(2(a,i4),a)')" Will sum ",my_nbks ," (b, k, s) occupied states in Sigma_x."
340  call wrtout(std_out,msg)
341 
342  ! The index of G-G0 in the FFT mesh the oscillators
343  ! Sigp%mG0 gives the MAX G0 component to account for umklapp.
344  ! Note the size MAX(npwx, Sigp%npwc).
345  ABI_MALLOC(igfftxg0, (Gsph_x%ng))
346 
347  ! Precalculate the FFT index of $ R^{-1}(r-\tau)$
348  ! S = \transpose R^{-1} and k_BZ = S k_IBZ
349  ! irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
350  gwx_nfftot = PRODUCT(gwx_ngfft(1:3))
351  ABI_MALLOC(irottb,(gwx_nfftot,Cryst%nsym))
352  call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwx_ngfft,irottb,iscompatibleFFT)
353  if (.not. iscompatibleFFT) then
354    MSG_WARNING("FFT mesh is not compatible with symmetries. Results might be affected by large errors!")
355  end if
356 
357  ABI_MALLOC(ktabr,(gwx_nfftot,Kmesh%nbz))
358  do ik_bz=1,Kmesh%nbz
359    isym=Kmesh%tabo(ik_bz)
360    do ifft=1,gwx_nfftot
361      ktabr(ifft,ik_bz)=irottb(ifft,isym)
362    end do
363  end do
364  ABI_FREE(irottb)
365 
366  if (Psps%usepaw==1 .and. pawcross==1) then
367    nfftf = PRODUCT(ngfftf(1:3))
368    ABI_MALLOC(irottb,(nfftf,Cryst%nsym))
369    call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,ngfftf,irottb,iscompatibleFFT)
370 
371    ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz))
372    do ik_bz=1,Kmesh%nbz
373      isym=Kmesh%tabo(ik_bz)
374      do ifft=1,nfftf
375        ktabrf(ifft,ik_bz)=irottb(ifft,isym)
376      end do
377    end do
378    ABI_FREE(irottb)
379  end if
380 
381  ! Additional allocations for PAW.
382  if (Psps%usepaw==1) then
383    ABI_DT_MALLOC(Cprj_ksum,(Cryst%natom,nspinor))
384    call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm)
385 
386    nhat12_grdim=0
387    if (use_pawnhat==1) then
388      ! Compensation charge for \phi_a^*\phi_b
389      call wrtout(std_out,"Using nhat12","COLL")
390      ABI_MALLOC(nhat12  ,(2,gwx_nfftot,nspinor**2))
391      ABI_MALLOC(grnhat12,(2,gwx_nfftot,nspinor**2,3*nhat12_grdim))
392    end if
393  end if
394 
395  ABI_MALLOC(sigxme_tmp, (minbnd:maxbnd, minbnd:maxbnd, Wfd%nsppol*Sigp%nsig_ab))
396  ABI_MALLOC(sigx,(2, ib1:ib2, ib1:ib2, nsppol*Sigp%nsig_ab))
397  sigxme_tmp = czero; sigx = czero
398 
399  nq_summed=Kmesh%nbz
400  if (Sigp%symsigma > 0) then
401    call littlegroup_print(Ltg_k,std_out,prtvol,'COLL')
402    nq_summed=SUM(Ltg_k%ibzq(:))
403 
404    ! Find number of degenerates subspaces and number of bands in each subspace.
405    ! The tolerance is a little bit arbitrary (0.001 eV)
406    ! It could be reduced, in particular in case of nearly accidental degeneracies
407    ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,nsppol))
408    degtab=0
409    do spin=1,nsppol
410      do ib=ib1,ib2
411        do jb=ib1,ib2
412         if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin))<0.001/Ha_ev) degtab(ib,jb,spin)=1
413        end do
414      end do
415    end do
416  end if !symsigma
417 
418  write(msg,'(2a,i0,a)')ch10,' calc_sigx_me: calculation status (', nq_summed, ' to be completed):'
419  call wrtout(std_out,msg,'COLL')
420 
421  ABI_MALLOC(ur_ibz,(gwx_nfftot*nspinor))
422  if (pawcross==1) then
423    ABI_MALLOC(ur_ae_sum,(nfftf*nspinor))
424    ABI_MALLOC(ur_ae_onsite_sum,(nfftf*nspinor))
425    ABI_MALLOC(ur_ps_onsite_sum,(nfftf*nspinor))
426  end if
427 
428  ! =======================================
429  ! ==== Begin loop over k_i in the BZ ====
430  ! =======================================
431  do spin=1,nsppol
432    if (ALL(proc_distrb(:,:,spin) /= Wfd%my_rank)) CYCLE
433 
434    ! Load wavefunctions for GW corrections.
435    ABI_STAT_MALLOC(wfr_bdgw, (gwx_nfftot*nspinor,ib1:ib2), ierr)
436    ABI_CHECK(ierr==0, "Out of memory in wfr_bdgw")
437    call wfd_get_many_ur(Wfd, [(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw)
438 
439    if (Wfd%usepaw==1) then
440      ! Load cprj for GW states, note the indexing.
441      dimcprj_gw=nspinor*(ib2-ib1+1)
442      ABI_DT_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1))
443      call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm)
444      ibsp=ib1
445      do jb=ib1,ib2
446        call wfd_get_cprj(Wfd,jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
447        call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
448        call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
449        ibsp=ibsp+nspinor
450      end do
451      if (pawcross==1) then
452        ABI_MALLOC(ur_ae_bdgw,(nfftf*nspinor,ib1:ib2))
453        ABI_MALLOC(ur_ae_onsite_bdgw,(nfftf*nspinor,ib1:ib2))
454        ABI_MALLOC(ur_ps_onsite_bdgw,(nfftf*nspinor,ib1:ib2))
455        do jb=ib1,ib2
456          call wfd_paw_get_aeur(Wfdf,jb,jk_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
457 &          ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
458          ur_ae_bdgw(:,jb)=ur_ae_sum
459          ur_ae_onsite_bdgw(:,jb)=ur_ae_onsite_sum
460          ur_ps_onsite_bdgw(:,jb)=ur_ps_onsite_sum
461        end do
462      end if
463    end if
464 
465    do ik_bz=1,Kmesh%nbz
466      ! Parallelization over k-points and spin.
467      ! For the spin there is another check in the inner loop
468      if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE
469 
470      ! Find the corresponding irreducible k-point
471      call get_BZ_item(Kmesh,ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt)
472      spinrot_kbz = Cryst%spinrot(:,isym_ki)
473 
474      ! Identify q and G0 where q+G0=k_GW-k_i
475      kgw_m_ksum = kgw - ksum
476      call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0)
477 
478      ! Symmetrize the matrix elements ===
479      ! Sum only q"s in IBZ_k. In this case elements are weighted
480      ! according to wtqp and wtqm. wtqm is for time-reversal.
481      wtqp=1; wtqm=0
482      if (can_symmetrize(spin)) then
483        if (Ltg_k%ibzq(iq_bz)/=1) CYCLE
484        wtqp=0; wtqm=0
485        do isym=1,Ltg_k%nsym_sg
486          wtqp = wtqp + Ltg_k%wtksym(1,isym,iq_bz)
487          wtqm = wtqm + Ltg_k%wtksym(2,isym,iq_bz)
488        end do
489      end if
490 
491      write(msg,'(2(a,i4),a,i3)')' calc_sigx_me: ik_bz ',ik_bz,'/',Kmesh%nbz,' done by mpi-rank: ',Wfd%my_rank
492      call wrtout(std_out,msg,'PERS')
493 
494      ! Find the corresponding irreducible q-point.
495      call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q)
496      q_is_gamma = (normv(qbz,Cryst%gmet,"G") < GW_TOL_W0)
497 
498      ! Tables for the FFT of the oscillators.
499      !  a) FFT index of the G-G0.
500      !  b) gwx_gbound table for the zero-padded FFT performed in rhotwg.
501      ABI_MALLOC(gwx_gbound,(2*gwx_mgfft+8,2))
502      call gsph_fft_tabs(Gsph_x,g0,gwx_mgfft,gwx_ngfft,use_padfft,gwx_gbound,igfftxg0)
503 
504      if (ANY(gwx_fftalga == [2, 4])) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
505      !use_padfft=0
506      if (use_padfft==0) then
507        ABI_FREE(gwx_gbound)
508        ABI_MALLOC(gwx_gbound,(2*gwx_mgfft+8,2*use_padfft))
509      end if
510 
511      if (pawcross==1) then
512        ABI_MALLOC(gboundf,(2*mgfftf+8,2))
513        ABI_MALLOC(igfftfxg0,(Gsph_x%ng))
514        call gsph_fft_tabs(Gsph_x,g0,mgfftf,ngfftf,use_padfftf,gboundf,igfftfxg0)
515        if ( ANY(gwx_fftalga == [2, 4]) ) use_padfftf=0
516        if (use_padfftf==0) then
517          ABI_FREE(gboundf)
518          ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf))
519        end if
520      end if
521 
522      if (Psps%usepaw==1 .and. use_pawnhat==0) then
523        ! Evaluate oscillator matrix elements
524        ! $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form
525        q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
526        ABI_DT_MALLOC(Pwij_qg,(Psps%ntypat))
527        call pawpwij_init(Pwij_qg, npwx, q0, Gsph_x%gvec, Cryst%rprimd, Psps, Pawtab, Paw_pwff)
528      end if
529 
530      ! Get Fourier components of the Coulomb interaction in the BZ
531      ! In 3D systems, neglecting umklapp,  vc(Sq,sG)=vc(q,G)=4pi/|q+G|
532      ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
533      do ig=1,npwx
534        ig_rot = Gsph_x%rottb(ig,itim_q,isym_q)
535        vc_sqrt_qbz(ig_rot)=Vcp%vc_sqrt_resid(ig,iq_ibz)
536      end do
537 
538      ! Sum over (occupied) bands.
539      do ib_sum=1,Sigp%nbnds
540        ! Parallelism over bands
541        ! This processor has this k-point but what about spin?
542        if (proc_distrb(ib_sum,ik_bz,spin)/=Wfd%my_rank) CYCLE
543 
544        ! Skip empty states.
545        if (qp_occ(ib_sum,ik_ibz,spin)<tol_empty) CYCLE
546 
547        call wfd_get_ur(Wfd,ib_sum,ik_ibz,spin,ur_ibz)
548 
549        if (Psps%usepaw==1) then
550          ! Load cprj for point ksum, this spin or spinor and *THIS* band.
551          ! TODO MG I could avoid doing this but I have to exchange spin and bands ???
552          ! For sure there is a better way to do this!
553          call wfd_get_cprj(Wfd,ib_sum,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
554          call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
555          if (pawcross==1) then
556            call wfd_paw_get_aeur(Wfdf,ib_sum,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
557 &              ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
558          end if
559        end if
560 
561        ! Get all <k-q,ib_sum,s|e^{-i(q+G).r}|s,jb,k>
562        do jb=ib1,ib2
563 
564          if (Psps%usepaw==1 .and. use_pawnhat==1) then
565            MSG_ERROR("use_pawnhat is disabled")
566            i2=jb; if (nspinor==2) i2=(2*jb-1)
567            spad=(nspinor-1)
568 
569            izero=0
570            call pawmknhat_psipsi(Cprj_ksum,Cprj_kgw(:,i2:i2+spad),ider0,izero,Cryst%natom,&
571 &            Cryst%natom,gwx_nfftot,gwx_ngfft,nhat12_grdim,nspinor,Cryst%ntypat,Pawang,Pawfgrtab,&
572 &            grnhat12,nhat12,pawtab)
573 
574          else
575            call rho_tw_g(nspinor,npwx,gwx_nfftot,ndat1,gwx_ngfft,1,use_padfft,igfftxg0,gwx_gbound, &
576 &            ur_ibz        ,iik,ktabr(:,ik_bz),ph_mkt  ,spinrot_kbz, &
577 &            wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw, &
578 &            nspinor,rhotwg_ki(:,jb))
579 
580            if (Psps%usepaw==1.and.use_pawnhat==0) then
581              ! Add on-site contribution, projectors are already in BZ.
582              i2=jb; if (nspinor==2) i2=(2*jb-1)
583              spad=(nspinor-1)
584              call paw_rho_tw_g(npwx,nspinor,nspinor,Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xred,Gsph_x%gvec,&
585 &              Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),Pwij_qg,rhotwg_ki(:,jb))
586            end if
587            if (Psps%usepaw==1.and.pawcross==1) then ! Add paw cross term
588              call paw_cross_rho_tw_g(nspinor,npwx,nfftf,ngfftf,1,use_padfftf,igfftfxg0,gboundf,&
589 &             ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum,iik,ktabrf(:,ik_bz),ph_mkt,spinrot_kbz,&
590 &             ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
591 &             nspinor,rhotwg_ki(:,jb))
592            end if
593          end if
594 
595          ! Multiply by the square root of the Coulomb term
596          ! In 3-D systems, the factor sqrt(4pi) is included
597          do ii=1,nspinor
598            spad = (ii-1) * npwx
599            rhotwg_ki(spad+1:spad+npwx,jb) = rhotwg_ki(spad+1:spad + npwx,jb) * vc_sqrt_qbz(1:npwx)
600          end do
601 
602          if (ik_bz == jk_bz) then
603            ! Treat analytically the case q --> 0
604            ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma
605            !   while the Colulomb term is integrated out.
606            ! * In the scalar case we have nonzero contribution only if ib_sum==jb
607            ! * For nspinor==2 evalute <ib_sum,up|jb,up> and <ib_sum,dwn|jb,dwn>,
608            !   impose orthonormalization since npwwfn might be < npwvec.
609 
610            if (nspinor==1) then
611              rhotwg_ki(1,jb) = czero_gw
612              !Note the use of i_sz_resid and not i_sz, to account for the possibility to have generalized KS basis set from hybrid
613              if (ib_sum == jb) rhotwg_ki(1,jb) = CMPLX(SQRT(Vcp%i_sz_resid), 0.0_gwp)
614              !rhotwg_ki(1,jb) = czero_gw ! DEBUG
615            else
616              npw_k = Wfd%npwarr(ik_ibz)
617              rhotwg_ki(1, jb) = zero; rhotwg_ki(npwx+1, jb) = zero
618              if (ib_sum == jb) then
619                cg_sum => Wfd%Wave(ib_sum,ik_ibz,spin)%ug
620                cg_jb  => Wfd%Wave(jb,jk_ibz,spin)%ug
621                ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1)
622                rhotwg_ki(1, jb) = CMPLX(SQRT(Vcp%i_sz_resid), 0.0_gwp) * real(ctmp)
623                ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1)
624                rhotwg_ki(npwx+1, jb) = CMPLX(SQRT(Vcp%i_sz_resid), 0.0_gwp) * real(ctmp)
625              end if
626              !rhotwg_ki(1, jb) = zero; rhotwg_ki(npwx+1, jb) = zero
627              ! PAW is missing
628            end if
629          end if
630        end do ! jb Got all matrix elements from minbnd up to maxbnd.
631 
632        theta_mu_minus_esum = fact_sp * qp_occ(ib_sum,ik_ibz,spin)
633 
634        do kb=ib1,ib2
635          ! Compute the ket Sigma_x |phi_{k,kb}>.
636          rhotwgp(:) = rhotwg_ki(:,kb)
637 
638          ! Loop over the non-zero row elements of this column.
639          ! If gwcalctyp <  20: only diagonal elements since QP == KS.
640          ! If gwcalctyp >= 20:
641          !      * Only off-diagonal elements connecting states with same character.
642          !      * Only the upper triangle if HF, SEX, or COHSEX.
643          do irow=1,Sigxij_tab(spin)%col(kb)%size1
644            jb = Sigxij_tab(spin)%col(kb)%bidx(irow)
645            rhotwg = rhotwg_ki(:,jb)
646 
647            ! Calculate bare exchange <phi_j|Sigma_x|phi_k>.
648            ! Do the scalar product only if ib_sum is occupied.
649            if (theta_mu_minus_esum/fact_sp >= tol_empty) then
650              do iab=1,Sigp%nsig_ab
651                spadx1 = spinor_padx(1, iab); spadx2 = spinor_padx(2, iab)
652                gwpc_sigxme = -XDOTC(npwx, rhotwg(spadx1+1:), 1, rhotwgp(spadx2+1:), 1) * theta_mu_minus_esum
653 
654                ! Accumulate and symmetrize Sigma_x
655                ! -wtqm comes from time-reversal (exchange of band indeces)
656                is_idx = spin; if (nspinor == 2) is_idx = iab
657                sigxme_tmp(jb, kb, is_idx) = sigxme_tmp(jb, kb, is_idx) + &
658 &                (wtqp + wtqm)*DBLE(gwpc_sigxme) + (wtqp - wtqm)*j_gw*AIMAG(gwpc_sigxme)
659 
660                sigx(1, jb, kb, is_idx) = sigx(1, jb, kb, is_idx) + wtqp *      gwpc_sigxme
661                sigx(2, jb, kb, is_idx) = sigx(2, jb, kb, is_idx) + wtqm *CONJG(gwpc_sigxme)
662              end do
663            end if
664          end do ! jb used to calculate matrix elements of Sigma_x
665 
666        end do ! kb to calculate matrix elements of Sigma_x
667      end do ! ib_sum
668 
669      ! Deallocate k-dependent quantities.
670      ABI_FREE(gwx_gbound)
671      if (pawcross==1) then
672        ABI_FREE(gboundf)
673      end if
674 
675      if (Psps%usepaw==1.and.use_pawnhat==0) then
676        call pawpwij_free(Pwij_qg)
677        ABI_DT_FREE(Pwij_qg)
678      end if
679    end do ! ik_bz Got all diagonal (off-diagonal) matrix elements.
680 
681    ABI_FREE(wfr_bdgw)
682    if (Wfd%usepaw==1) then
683      call pawcprj_free(Cprj_kgw)
684      ABI_DT_FREE(Cprj_kgw)
685      if (pawcross==1) then
686        ABI_FREE(ur_ae_bdgw)
687        ABI_FREE(ur_ae_onsite_bdgw)
688        ABI_FREE(ur_ps_onsite_bdgw)
689      end if
690    end if
691  end do !spin
692 
693  ABI_FREE(igfftxg0)
694  if (pawcross==1) then
695    ABI_FREE(igfftfxg0)
696  end if
697 
698  ! Gather contributions from all the CPUs.
699  call xmpi_sum(sigxme_tmp, wfd%comm, ierr)
700  call xmpi_sum(sigx, wfd%comm, ierr)
701 
702  ! Multiply by constants. For 3D systems sqrt(4pi) is included in vc_sqrt_qbz.
703  sigxme_tmp = (one/(Cryst%ucvol*Kmesh%nbz)) * sigxme_tmp * Sigp%sigma_mixing
704  sigx       = (one/(Cryst%ucvol*Kmesh%nbz)) * sigx       * Sigp%sigma_mixing
705 
706  !
707  ! If we have summed over the IBZ_q now we have to average over degenerate states.
708  ! NOTE: Presently only diagonal terms are considered
709  ! TODO QP-SCGW required a more involved approach, there is a check in sigma
710  ! TODO it does not work if spinor==2.
711  do spin=1,nsppol
712    if (can_symmetrize(spin)) then
713      ABI_MALLOC(sym_sigx, (ib1:ib2, ib1:ib2, sigp%nsig_ab))
714      sym_sigx = czero
715 
716      ! Average over degenerate diagonal elements.
717      ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results.
718      ! another good reason to use a strict criterion for the tollerance on eigenvalues.
719      do ib=ib1,ib2
720        ndegs=0
721        do jb=ib1,ib2
722          if (degtab(ib,jb,spin)==1) then
723            if (nspinor == 1) then
724              sym_sigx(ib, ib, 1) = sym_sigx(ib, ib, 1) + SUM(sigx(:,jb,jb,spin))
725            else
726              do ii=1,sigp%nsig_ab
727                sym_sigx(ib, ib, ii) = sym_sigx(ib, ib, ii) + SUM(sigx(:,jb,jb,ii))
728              end do
729            end if
730          end if
731          ndegs = ndegs + degtab(ib,jb,spin)
732        end do
733        sym_sigx(ib,ib,:) = sym_sigx(ib,ib,:) / ndegs
734      end do
735 
736      if (gwcalctyp >= 20) then
737        call esymm_symmetrize_mels(QP_sym(spin),ib1,ib2,sigx(:,:,:,spin),sym_sigx(:,:,1))
738      end if
739 
740      ! Copy symmetrized values.
741      do ib=ib1,ib2
742        do jb=ib1,ib2
743          if (nspinor == 1) then
744            sigxme_tmp(ib,jb,spin) = sym_sigx(ib,jb,1)
745          else
746            do ii=1,sigp%nsig_ab
747              sigxme_tmp(ib,jb,ii) = sym_sigx(ib,jb,ii)
748            end do
749          end if
750        end do
751      end do
752 
753      ABI_FREE(sym_sigx)
754    end if
755  end do
756 
757  if (gwcalctyp>=20) then
758    ! Reconstruct the full sigma_x matrix from the upper triangle.
759    if (nspinor == 1) then
760      do spin=1,nsppol
761        call hermitianize(sigxme_tmp(:,:,spin), "Upper")
762      end do
763    else
764      MSG_WARNING("Should hermitianize non-collinear sigma!")
765    end if
766  end if
767 
768  ! Save diagonal elements or ab components of Sigma_x (hermitian)
769  ! TODO It should be hermitian also if nspinor==2
770  do spin=1,nsppol
771    do jb=ib1,ib2
772      do iab=1,Sigp%nsig_ab
773        is_idx = spin; if (Sigp%nsig_ab > 1) is_idx = iab
774        if (is_idx <= 2) then
775          Sr%sigxme(jb,jk_ibz,is_idx) = DBLE(sigxme_tmp(jb,jb,is_idx))
776        else
777          Sr%sigxme(jb,jk_ibz,is_idx) = sigxme_tmp(jb,jb,is_idx)
778        end if
779      end do
780      !if (Sigp%nsig_ab>1) then
781      !  write(std_out,'(i3,4f8.3,a,f8.3)')jb,Sr%sigxme(jb,jk_ibz,:)*Ha_eV,' Tot ',SUM(Sr%sigxme(jb,jk_ibz,:))*Ha_eV
782      !end if
783    end do
784  end do
785 
786  ! Save full exchange matrix in Sr%
787  Sr%x_mat(minbnd:maxbnd,minbnd:maxbnd,jk_ibz,:) = sigxme_tmp(minbnd:maxbnd,minbnd:maxbnd,:)
788  ABI_FREE(sigxme_tmp)
789 
790  ! ===========================
791  ! ==== Deallocate memory ====
792  ! ===========================
793  if (Psps%usepaw==1) then
794    if (allocated(gwx_gfft))  then
795      ABI_FREE(gwx_gfft)
796    end if
797    call pawcprj_free(Cprj_ksum)
798    ABI_DT_FREE(Cprj_ksum)
799    if (allocated(Pwij_fft)) then
800      call pawpwij_free(Pwij_fft)
801      ABI_DT_FREE(Pwij_fft)
802    end if
803    if (use_pawnhat==1) then
804      ABI_FREE(nhat12)
805      ABI_FREE(grnhat12)
806    end if
807    if (pawcross==1) then
808      ABI_FREE(ur_ae_sum)
809      ABI_FREE(ur_ae_onsite_sum)
810      ABI_FREE(ur_ps_onsite_sum)
811      ABI_FREE(ktabrf)
812    end if
813  end if
814 
815  ABI_FREE(ur_ibz)
816  ABI_FREE(rhotwg_ki)
817  ABI_FREE(rhotwg)
818  ABI_FREE(rhotwgp)
819  ABI_FREE(vc_sqrt_qbz)
820  ABI_FREE(ktabr)
821  ABI_FREE(sigx)
822  ABI_FREE(proc_distrb)
823 
824  if (allocated(degtab)) then
825    ABI_FREE(degtab)
826  end if
827 
828  call timab(430,2,tsec) ! csigme (SigX)
829  call cwtime(cpu_time,wall_time,gflops,"stop")
830  write(std_out,'(2(a,f9.1))')" cpu_time = ",cpu_time,", wall_time = ",wall_time
831 
832  DBG_EXIT("COLL")
833 
834 end subroutine calc_sigx_me

ABINIT/m_sigx [ Modules ]

[ Top ] [ Modules ]

NAME

  m_six

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (FB, GMR, VO, LR, RWG, MG, RShaltaf)
  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_sigx
28 
29  implicit none
30 
31  private