Line data Source code
1 : !!****m* ABINIT/m_cohsex
2 : !! NAME
3 : !! m_cohsex
4 : !!
5 : !! FUNCTION
6 : !! Calculate diagonal and off-diagonal matrix elements of the SEX or COHSEX self-energy operator.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (FB, GMR, VO, LR, RWG, MG, RShaltaf)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_cohsex
23 :
24 : use defs_basis
25 : use m_defs_ptgroups
26 : use m_gwdefs
27 : use m_xmpi
28 : use m_errors
29 : use m_abicore
30 :
31 : use defs_datatypes, only : pseudopotential_type
32 : use m_dtset, only : dataset_type
33 : use m_time, only : timab, cwtime, cwtime_report
34 : use m_fstrings, only : sjoin, itoa
35 : use m_hide_blas, only : xdotc, xgemv
36 : use m_numeric_tools, only : hermitianize, imin_loc
37 : use m_geometry, only : normv
38 : use m_crystal, only : crystal_t
39 : use m_ebands, only : ebands_t
40 : use m_bz_mesh, only : kmesh_t, findqg0, littlegroup_t
41 : use m_gsphere, only : gsphere_t
42 : use m_fft_mesh, only : get_gfft, rotate_fft_mesh, cigfft
43 : use m_vcoul, only : vcoul_t
44 : use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g
45 : use m_wfd, only : wfdgw_t, wave_t
46 : use m_oscillators, only : rho_tw_g, calc_wfwfg
47 : use m_screening, only : epsm1_t
48 : use m_esymm, only : esymm_t
49 : use m_sigma, only : sigma_t, sigma_distribute_bks
50 : use m_pawang, only : pawang_type
51 : use m_pawtab, only : pawtab_type
52 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
53 : use m_paw_sym, only : paw_symcprj
54 :
55 : implicit none
56 :
57 : private
58 : !!***
59 :
60 : public :: cohsex_me
61 : !!***
62 :
63 : integer,parameter :: LOG_MODK = 5
64 :
65 : contains
66 : !!***
67 :
68 : !!****f* ABINIT/cohsex_me
69 : !! NAME
70 : !! cohsex_me
71 : !!
72 : !! FUNCTION
73 : !! Calculate diagonal and off-diagonal matrix elements of the SEX or COHSEX self-energy operator.
74 : !!
75 : !! COPYRIGHT
76 : !! Copyright (C) 1999-2026 ABINIT group (FB, GMR, VO, LR, RWG, MG, RShaltaf)
77 : !! This file is distributed under the terms of the
78 : !! GNU General Public License, see ~abinit/COPYING
79 : !! or http://www.gnu.org/copyleft/gpl.txt .
80 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
81 : !!
82 : !! INPUTS
83 : !! sigmak_ibz=Index of the k-point in the IBZ.
84 : !! minbnd, maxbnd= min and Max band index for GW correction (for this k-point)
85 : !! iomode=Option defining the file format of the SCR file (Fortran, NETCDF)
86 : !! epsm1 <epsm1_t> (see the definition of this structured datatype)
87 : !! %mqmem=if 0 use out-of-core method in which a single q-slice of espilon is read inside the loop over k
88 : !! %nomega_i=Number of imaginary frequencies.
89 : !! %nomega_r=Number of real frequencies.
90 : !! %nomega=Total number of frequencies.
91 : !! Gsph_c<gsphere_t>= info on the G-sphere for Sigma_x
92 : !! %nsym=number of symmetry operations
93 : !! %rottb(ng,timrev,nsym)=index of (IS) G where I is the identity or the inversion
94 : !! operation and G is one of the ng vectors in reciprocal space
95 : !! %timrev=2 if time-reversal symmetry is used, 1 otherwise
96 : !! %gvec(3,Sigp%npwc)=integer coordinates of each plane wave in reciprocal space
97 : !! ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
98 : !! Ltg_k datatype containing information on the little group
99 : !! Kmesh <kmesh_t>
100 : !! %nbz=Number of points in the BZ
101 : !! %nibz=Number of points in IBZ
102 : !! %kibz(3,nibz)=k-point coordinates, irreducible Brillouin zone
103 : !! %kbz(3,nbz)=k-point coordinates, full Brillouin zone
104 : !! %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
105 : !! %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered
106 : !! %ktabp(nbz)= phase factor associated to tnons
107 : !! gwc_ngfft(18)=Information about 3D FFT for the oscillator strengths used for the correlation part,
108 : !! Vcp <vcoul_t datatype> containing information on the cutoff technique
109 : !! %vc_sqrt(npwc,nqibz)= square-root of the coulombian potential for q-points in the IBZ
110 : !! Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
111 : !! Pawang <type(pawang_type)>=paw angular mesh and related data
112 : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
113 : !! %usepaw=1 for PAW, 0 for NC pseudopotentials.
114 : !! Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
115 : !! %ibz=q points where $\tilde\epsilon^{-1}$ has been computed
116 : !! %bz(3,nqbz)=coordinates of all q-points in BZ
117 : !! Sigp <sigparams_t> (see the definition of this structured datatype)
118 : !! Cryst<crystal_t>=Info on unit cell and symmetries
119 : !! %natom=number of atoms in unit cell
120 : !! %ucvol=unit cell volume
121 : !! %nsym=number of symmetry operations
122 : !! %typat(natom)=type of each atom
123 : !! much slower but it requires less memory
124 : !! QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
125 : !! eig(Sigp%nbnds,Kmesh%nibz,%nsppol)=KS or QP energies for k-points, bands and spin
126 : !! occ(Sigp%nbnds,Kmesh%nibz,%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
127 : !! Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
128 : !! allQP_sym(%nkibz,%nsppol)<esymm_t>=Datatype collecting data on the irreducible representaions of the
129 : !! little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
130 : !! Sr=sigma_t (see the definition of this structured datatype)
131 : !!
132 : !! OUTPUT
133 : !!
134 : !! NOTES
135 : !! 1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) [[cite:Gigy1986]] is included.
136 : !! 2) The calculation of energy derivative is based on finite elements.
137 : !! 3) On the symmetrization of Sigma matrix elements ***/
138 : !! If Sk = k+G0 then M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G} (k,q)
139 : !! If -Sk = k+G0 then M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)
140 : !!
141 : !! Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out
142 : !! that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a
143 : !! similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase
144 : !! has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ.
145 : !!
146 : !! 4) The unitary transformation relating wavefunctions
147 : !! at symmetric k-points should be taken into account during the symmetrization
148 : !! of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
149 : !! it is possible to make an invariant by just including all the degenerate states and
150 : !! averaging the final results over the degenerate subset.
151 : !!
152 : !! SOURCE
153 :
154 5 : subroutine cohsex_me(sigmak_ibz,ikcalc,nomega_sigc,minbnd,maxbnd,dtset, Cryst,QP_BSt,Sigp,Sr,epsm1,Gsph_c,Vcp,&
155 5 : & Kmesh,Qmesh,Ltg_k,Pawtab,Pawang,Paw_pwff,Psps,Wfd,allQP_sym,gwc_ngfft,iomode,prtvol,sigcme_tmp)
156 :
157 : !Arguments ------------------------------------
158 : !scalars
159 : integer,intent(in) :: sigmak_ibz,ikcalc,prtvol,iomode,nomega_sigc,minbnd,maxbnd
160 : type(dataset_type),intent(in) :: dtset
161 : type(crystal_t),intent(in) :: Cryst
162 : type(ebands_t),target,intent(in) :: QP_BSt
163 : type(kmesh_t),intent(in) :: Kmesh,Qmesh
164 : type(vcoul_t),intent(in) :: Vcp
165 : type(epsm1_t),intent(inout) :: epsm1
166 : type(gsphere_t),intent(in) :: Gsph_c
167 : type(littlegroup_t),intent(in) :: Ltg_k
168 : type(Pseudopotential_type),intent(in) :: Psps
169 : type(pawang_type),intent(in) :: pawang
170 : type(sigparams_t),target,intent(in) :: Sigp
171 : type(sigma_t),intent(in) :: Sr
172 : type(wfdgw_t),target,intent(inout) :: Wfd
173 : !arrays
174 : integer,intent(in) :: gwc_ngfft(18)
175 : complex(dp),intent(out) :: sigcme_tmp(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Wfd%nsppol*Sigp%nsig_ab)
176 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
177 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
178 : type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol)
179 :
180 : !Local variables ------------------------------
181 : !scalars
182 : integer,parameter :: ndat1=1
183 : integer :: iab,ib,ib1,ib2,ierr,ig,ii,iik,itim_q,i1,i2,npwc
184 : integer :: ik_bz,ik_ibz,io,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx
185 : integer :: band,band1,band2,idle,rank
186 : integer :: jik,jk_bz,jk_ibz,kb,nspinor,nsppol
187 : integer :: nomega_tot,nq_summed,ispinor,ibsp,dimcprj_gw
188 : integer :: spad,spadc,spadc1,spadc2,irow,my_nbks
189 : integer :: ndegs,wtqm,wtqp,mod10
190 : integer :: isym_kgw,isym_ki,gwc_mgfft,use_padfft,gwc_fftalga,gwc_nfftot,ifft,npw_k
191 : real(dp) :: fact_spin,theta_mu_minus_e0i,tol_empty,gw_gsq
192 : real(dp) :: cpu_all, wall_all, gflops_all, cpu_k, wall_k, gflops_k
193 : complex(dp) :: ctmp,ph_mkgwt,ph_mkt
194 : logical :: iscompatibleFFT, q_is_gamma, print_time
195 : character(len=500) :: msg
196 : type(wave_t),pointer :: wave_sum, wave_jb
197 : !arrays
198 10 : integer :: g0(3),spinor_padc(2,4),nbv_ks(Kmesh%nibz,Wfd%nsppol)
199 5 : integer,allocatable :: proc_distrb(:,:,:),coh_distrb(:,:,:,:),degtab(:,:,:)
200 5 : integer,allocatable :: igfftcg0(:),gw_gfft(:,:),gw_gbound(:,:),irottb(:,:),ktabr(:,:)
201 10 : integer :: got(Wfd%nproc)
202 : real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),q0(3),tsec(2),qbz(3),spinrot_kbz(4),spinrot_kgw(4)
203 5 : real(dp),contiguous, pointer :: qp_ene(:,:,:),qp_occ(:,:,:)
204 10 : complex(gwp) :: sigcohme(Sigp%nsig_ab)
205 5 : complex(gwp),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:),sigsex(:)
206 5 : complex(gwp),allocatable :: sigc_ket(:,:) ! epsm1_qbz(:,:,:),
207 5 : complex(gwp),allocatable :: rhotwg_ki(:,:), sigctmp(:,:)
208 5 : complex(gwp),allocatable :: wfr_bdgw(:,:),ur_sum(:),wf1swf2_g(:)
209 5 : complex(gwp),contiguous, pointer :: cg_jb(:),cg_sum(:)
210 5 : complex(dp),allocatable :: sym_cme(:,:,:,:),sigc(:,:,:,:,:)
211 15 : logical :: rank_mask(Wfd%nproc),can_symmetrize(Wfd%nsppol)
212 5 : logical,allocatable :: bks_mask(:,:,:)
213 5 : type(sigijtab_t),pointer :: Sigcij_tab(:)
214 5 : type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
215 5 : type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
216 5 : type(esymm_t),pointer :: QP_sym(:)
217 : !************************************************************************
218 :
219 : DBG_ENTER("COLL")
220 :
221 5 : call timab(423,1,tsec) ! cohsex_me
222 5 : call cwtime(cpu_all, wall_all, gflops_all,"start")
223 :
224 : ! Initial check
225 5 : ABI_CHECK(Sr%nomega_r == Sigp%nomegasr,"")
226 5 : ABI_CHECK(Sr%nomega4sd == Sigp%nomegasrd,"")
227 : !ABI_CHECK(Sigp%npwc==Gsph_c%ng,"")
228 :
229 : ! Initialize some values
230 5 : nspinor = Wfd%nspinor; nsppol = Wfd%nsppol
231 5 : npwc = sigp%npwc
232 45 : spinor_padc = RESHAPE([0, 0, npwc, npwc, 0, npwc, npwc,0], [2, 4])
233 :
234 5 : qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ
235 :
236 : ! Extract the symmetries of the bands for this k-point
237 5 : QP_sym => allQP_sym(sigmak_ibz,1:nsppol)
238 :
239 : ! Index of the GW point in the BZ array, its image in IBZ and time-reversal ===
240 5 : jk_bz=Sigp%kptgw2bz(ikcalc)
241 5 : call kmesh%get_BZ_item(jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
242 : !%call get_IBZ_item(Kmesh,jk_ibz,kibz,wtk)
243 25 : spinrot_kgw=Cryst%spinrot(:,isym_kgw)
244 5 : ib1 = minbnd; ib2 = maxbnd
245 :
246 5 : write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,&
247 5 : ' Calculating <nk|Sigma_c(omega)|nk> at k = ',kgw(:),ch10,&
248 10 : ' bands n = from ',ib1,' to ',ib2,ch10
249 5 : call wrtout(std_out,msg)
250 :
251 20 : if (ANY(gwc_ngfft(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst,Psps,gwc_ngfft)
252 20 : gwc_mgfft = MAXVAL(gwc_ngfft(1:3))
253 5 : gwc_fftalga = gwc_ngfft(7)/100 !; gwc_fftalgc=MOD(gwc_ngfft(7),10)
254 :
255 10 : can_symmetrize = .FALSE.
256 5 : if (Sigp%symsigma>0) then
257 0 : can_symmetrize = .TRUE.
258 0 : if (Sigp%gwcalctyp >= 20) then
259 0 : do spin=1,Wfd%nsppol
260 0 : can_symmetrize(spin) = .not. QP_sym(spin)%failed()
261 0 : if (.not.can_symmetrize(spin)) then
262 0 : write(msg,'(a,i0,4a)')" Symmetrization cannot be performed for spin: ",spin,ch10,&
263 0 : " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg)
264 0 : ABI_WARNING(msg)
265 : end if
266 : end do
267 : end if
268 0 : if (nspinor == 2) ABI_WARNING('Symmetrization with nspinor=2 not implemented')
269 : end if
270 :
271 5 : mod10=MOD(Sigp%gwcalctyp, 10)
272 :
273 5 : call timab(491,1,tsec) ! csigme(tot) Overall clock. TODO check this
274 5 : call timab(495,1,tsec) ! csigme (SigC)
275 :
276 : ! Normalization of theta_mu_minus_e0i
277 : ! If nsppol==2, qp_occ $\in [0,1]$
278 : !call get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
279 5 : SELECT CASE (nsppol)
280 : CASE (1)
281 5 : fact_spin=half; tol_empty=0.01 ! below this value the state is assumed empty
282 5 : if (nspinor==2) then
283 1 : fact_spin=one; tol_empty=0.005 ! below this value the state is assumed empty
284 : end if
285 : CASE (2)
286 0 : fact_spin=one; tol_empty=0.005 ! to be consistent and obtain similar results if a metallic
287 : CASE DEFAULT ! spin unpolarized system is treated using nsppol==2
288 5 : ABI_BUG('Wrong nsppol')
289 : END SELECT
290 :
291 5 : call timab(442,1,tsec) ! csigme(init0)
292 :
293 : ! Precalculate the FFT index of $(R^{-1}(r-\tau))$ ===
294 : ! S=\transpose R^{-1} and k_BZ = S k_IBZ
295 : ! irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
296 20 : gwc_nfftot = PRODUCT(gwc_ngfft(1:3))
297 20 : ABI_MALLOC(irottb,(gwc_nfftot,Cryst%nsym))
298 5 : call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwc_ngfft,irottb,iscompatibleFFT)
299 5 : if (.not.iscompatibleFFT) then
300 0 : ABI_WARNING("FFT mesh is not compatible with symmetries. Results might be affected by large errors!")
301 : end if
302 :
303 20 : ABI_MALLOC(ktabr,(gwc_nfftot, Kmesh%nbz))
304 269 : do ik_bz=1,Kmesh%nbz
305 264 : isym=Kmesh%tabo(ik_bz)
306 2112269 : do ifft=1,gwc_nfftot
307 2112264 : ktabr(ifft,ik_bz)=irottb(ifft,isym)
308 : end do
309 : end do
310 5 : ABI_FREE(irottb)
311 :
312 : ! The number of occupied states for each point in the IBZ and spin.
313 10 : do spin=1,nsppol
314 45 : do ik_ibz=1,Kmesh%nibz
315 284 : nbv_ks(ik_ibz,spin) = COUNT(qp_occ(:,ik_ibz,spin)>=tol_empty)
316 : end do
317 : end do
318 :
319 : ! (b,k,s) mask for MPI distribution of the sum over occupied states in the BZ.
320 25 : ABI_MALLOC(bks_mask,(Wfd%mband,Kmesh%nbz,nsppol))
321 1778 : bks_mask=.FALSE.
322 10 : do spin=1,nsppol
323 274 : do ik_bz=1,Kmesh%nbz
324 264 : ik_ibz = Kmesh%tab(ik_bz)
325 1357 : bks_mask(1:nbv_ks(ik_ibz,spin),ik_bz,spin) = .TRUE.
326 : end do
327 : end do
328 :
329 : ! Distribute the individual terms of the sum over the BZ taking into account symmetries and MPI memory distribution.
330 : ! got is used to optimize the distribution if more than one node can calculate the same (b,k,s) element.
331 10 : got=0
332 20 : ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,nsppol))
333 : call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,&
334 5 : proc_distrb,got,bks_mask,global=.TRUE.)
335 :
336 5 : ABI_FREE(bks_mask)
337 :
338 5 : write(msg,'(a,i0,a)')" Will sum ",my_nbks," (b,k,s) occupied states in (COHSEX|SEX)."
339 5 : call wrtout(std_out,msg)
340 :
341 5 : Sigcij_tab => Sigp%Sigcij_tab(ikcalc,1:nsppol)
342 :
343 5 : if (mod10==SIG_COHSEX) then
344 : ! Distribute the COHSEX terms, taking into account the symmetries of the Sigma_ij matrix.
345 18 : ABI_MALLOC(coh_distrb,(ib1:ib2,ib1:ib2,Kmesh%nbz,nsppol))
346 :
347 2590 : coh_distrb = xmpi_undefined_rank
348 6 : do spin=1,nsppol
349 142 : do ik_bz=1,Kmesh%nbz
350 139 : if (ANY(proc_distrb(:,ik_bz,spin) /= xmpi_undefined_rank) ) then ! This BZ point will be calculated.
351 272 : rank_mask = .FALSE. ! To select only those nodes that will treat (k,s).
352 1000 : do band=1,Wfd%mband
353 864 : rank = proc_distrb(band,ik_bz,spin)
354 1000 : if (rank /= xmpi_undefined_rank) rank_mask(rank+1)=.TRUE.
355 : end do
356 504 : do band2=ib1,ib2
357 872 : do irow=1,Sigcij_tab(spin)%col(band2)%size1 ! Looping over the upper triangle of sigma_ij with non-zero elements.
358 368 : band1 = Sigcij_tab(spin)%col(band2)%bidx(irow)
359 368 : idle = imin_loc(got,mask=rank_mask)
360 368 : got(idle) = got(idle)+1
361 736 : coh_distrb(band1,band2,ik_bz,spin) = idle-1
362 : end do
363 : end do
364 : end if
365 : end do
366 : end do
367 :
368 2590 : write(msg,'(a,i0,a)')" will treat ",COUNT(coh_distrb==Wfd%my_rank)," COH terms."
369 3 : call wrtout(std_out,msg)
370 : end if
371 :
372 20 : ABI_MALLOC(rhotwg_ki, (npwc * nspinor, minbnd:maxbnd))
373 1191 : rhotwg_ki=czero_gw
374 15 : ABI_MALLOC(rhotwg, (npwc * nspinor))
375 10 : ABI_MALLOC(rhotwgp ,(npwc * nspinor))
376 15 : ABI_MALLOC(vc_sqrt_qbz, (npwc))
377 :
378 : ! Additional allocations for PAW
379 5 : if (Psps%usepaw==1) then
380 0 : ABI_MALLOC(Cprj_ksum,(Cryst%natom,nspinor))
381 0 : call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm)
382 :
383 : ! For COHSEX we need the onsite terms of the PW on the FFT mesh.
384 : ! gw_gfft is the set of plane waves in the FFT Box for the oscillators.
385 0 : if (mod10==SIG_COHSEX) then
386 0 : ABI_MALLOC(gw_gfft,(3,gwc_nfftot))
387 0 : q0=zero
388 0 : call get_gfft(gwc_ngfft,q0,Cryst%gmet,gw_gsq,gw_gfft)
389 0 : ABI_MALLOC(Pwij_fft,(Psps%ntypat))
390 0 : call pawpwij_init(Pwij_fft,gwc_nfftot,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
391 : end if
392 : end if ! usepaw==1
393 :
394 : ! === Calculate total number of frequencies and allocate related arrays ===
395 : ! sigcme2 is used to accumulate the diagonal matrix elements over k-points and
396 : ! GW bands, used only in case of ppmodel 3 and 4 (TODO save memory)
397 5 : nomega_tot=Sr%nomega_r+Sr%nomega4sd
398 :
399 20 : ABI_MALLOC(sigctmp, (nomega_sigc,Sigp%nsig_ab))
400 21 : sigctmp = czero_gw
401 20 : ABI_MALLOC(sigc_ket, (npwc*nspinor, nomega_sigc))
402 :
403 5 : if (mod10==SIG_COHSEX) then
404 9 : ABI_MALLOC(wf1swf2_g,(gwc_nfftot*nspinor))
405 : end if
406 :
407 : ! Arrays storing the contribution given by the Hermitian/anti-Hermitian part of \Sigma_c
408 : !allocate(aherm_sigc_ket(npwc*nspinor,nomega_sigc))
409 : !allocate( herm_sigc_ket(npwc*nspinor,nomega_sigc))
410 15 : ABI_MALLOC(sigsex,(npwc))
411 1677 : sigcme_tmp=czero
412 :
413 30 : ABI_MALLOC(sigc,(2,nomega_sigc,ib1:ib2,ib1:ib2,nsppol*Sigp%nsig_ab))
414 3277 : sigc=czero
415 :
416 : ! Here we divide the states where the QP energies are required into complexes. Note however that this approach is not
417 : ! based on group theory, and it might lead to spurious results in case of accidental degeneracies.
418 5 : nq_summed=Kmesh%nbz
419 5 : if (Sigp%symsigma > 0) then
420 0 : call Ltg_k%print([std_out], prtvol=prtvol)
421 0 : nq_summed=SUM(Ltg_k%ibzq(:))
422 : !
423 : ! Find number of degenerate states and number of bands in each subspace
424 : ! The tolerance is a little bit arbitrary (0.001 eV)
425 : ! It could be reduced, in particular in case of nearly accidental degeneracies
426 0 : ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,nsppol))
427 0 : degtab=0
428 0 : do spin=1,nsppol
429 0 : do ib=ib1,ib2
430 0 : do jb=ib1,ib2
431 0 : if (ABS(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin)) < dtset%symsigma_de) then
432 0 : degtab(ib,jb,spin)=1
433 : end if
434 : end do
435 : end do
436 : end do
437 : end if !symsigma
438 :
439 5 : write(msg,'(2a,i6,a)')ch10,' calculation status ( ',nq_summed,' to be completed):'
440 5 : call wrtout(std_out,msg)
441 :
442 : ! TODO if single q (ex molecule) dont allocate epsm1q, avoid waste of memory
443 : !ABI_MALLOC_OR_DIE(epsm1_qbz, (npwc, npwc, 1), ierr)
444 5 : call epsm1%malloc_epsm1_qbz(npwc, 1)
445 :
446 15 : ABI_MALLOC(igfftcg0,(Gsph_c%ng))
447 :
448 : ! Out-of-core solution for epsilon.
449 5 : if (epsm1%mqmem==0) then
450 0 : ABI_COMMENT('Reading q-slices from file. Slower but less memory.')
451 : end if
452 :
453 : ! If epsm1 is MPI-shared, we have to start the RMA epoch. Note that epsm1%epsm1 is read-only.
454 5 : if (epsm1%use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1%epsm1_win, ierr)
455 :
456 5 : call timab(442,2,tsec)
457 :
458 : ! ==========================================
459 : ! ==== Fat loop over k_i in the full BZ ====
460 : ! ==========================================
461 15 : ABI_MALLOC(ur_sum,(gwc_nfftot*nspinor))
462 :
463 10 : do spin=1,nsppol
464 5 : if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE
465 :
466 20 : ABI_MALLOC(wfr_bdgw,(gwc_nfftot*nspinor,ib1:ib2))
467 49 : call wfd%get_many_ur([(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw)
468 :
469 5 : if (Wfd%usepaw==1) then
470 : ! Load cprj for GW states, note the indexing.
471 0 : dimcprj_gw=nspinor*(ib2-ib1+1)
472 0 : ABI_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1))
473 0 : call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm)
474 0 : ibsp=ib1
475 0 : do jb=ib1,ib2
476 0 : call wfd%get_cprj(jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
477 0 : call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
478 0 : call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
479 0 : ibsp=ibsp+nspinor
480 : end do
481 : end if
482 :
483 269 : do ik_bz=1,Kmesh%nbz
484 : ! Parallelization over k-points and spin
485 : ! For the spin there is another check in the inner loop.
486 264 : if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE
487 :
488 264 : call timab(443,1,tsec) ! csigme (initq)
489 264 : print_time = wfd%my_rank == 0 .and. (ik_bz < LOG_MODK .or. mod(ik_bz, LOG_MODK) == 0)
490 69 : if (print_time) call cwtime(cpu_k, wall_k, gflops_k, "start")
491 :
492 : ! Find the corresponding irreducible k-point
493 264 : call kmesh%get_BZ_item(ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt)
494 1320 : spinrot_kbz(:)=Cryst%spinrot(:,isym_ki)
495 :
496 : ! Identify q and G0 where q+G0=k_GW-k_i
497 1056 : kgw_m_ksum=kgw-ksum
498 264 : call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0)
499 :
500 : ! Symmetrize the matrix elements.
501 : ! Sum only q"s in IBZ_k. In this case elements are weighted
502 : ! according to wtqp and wtqm. wtqm is for time-reversal.
503 264 : wtqp=1; wtqm=0
504 264 : if (can_symmetrize(spin)) then
505 0 : if (Ltg_k%ibzq(iq_bz)/=1) CYCLE
506 0 : wtqp=0; wtqm=0
507 0 : do isym=1,Ltg_k%nsym_sg
508 0 : wtqp=wtqp+Ltg_k%wtksym(1,isym,iq_bz)
509 0 : wtqm=wtqm+Ltg_k%wtksym(2,isym,iq_bz)
510 : end do
511 : end if
512 :
513 : !%write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
514 : !%call wrtout(std_out,msg)
515 :
516 : ! Find the corresponding irred q-point.
517 264 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
518 264 : q_is_gamma = (normv(qbz, Cryst%gmet, "G") < GW_TOLQ0)
519 :
520 : ! Tables for the FFT of the oscillators.
521 : ! a) FFT index of the G-G0.
522 : ! b) gw_gbound table for the zero-padded FFT performed in rhotwg.
523 1056 : ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2))
524 264 : call Gsph_c%fft_tabs(g0,gwc_mgfft,gwc_ngfft,use_padfft,gw_gbound,igfftcg0)
525 264 : if ( ANY(gwc_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
526 264 : if (use_padfft==0) then
527 0 : ABI_FREE(gw_gbound)
528 0 : ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2*use_padfft))
529 : end if
530 :
531 264 : if (Psps%usepaw==1) then
532 : ! Get PAW oscillator matrix elements $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form.
533 0 : ABI_MALLOC(Pwij_qg,(Psps%ntypat))
534 0 : q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
535 0 : call pawpwij_init(Pwij_qg, npwc, q0, Gsph_c%gvec, Cryst%rprimd, Psps, Pawtab, Paw_pwff)
536 : end if
537 :
538 264 : if (epsm1%mqmem==0) then
539 : ! Read q-slice of epsilon^{-1}|chi0 in epsm1%epsm1(:,:,:,1) (much slower but less memory).
540 0 : call epsm1%get_epsm1(Vcp,0,0,iomode,xmpi_comm_self,iqibzA=iq_ibz)
541 : end if
542 :
543 : ! Only omega==0 for SEX or COHSEX
544 264 : call epsm1%rotate_iqbz(iq_bz, 1, npwc, Gsph_c, Qmesh, .True.) !, epsm1_qbz)
545 :
546 : ! Get Fourier components of the Coulomb interaction in the BZ.
547 : ! In 3D systems, neglecting umklapp, vc(Sq,sG)=vc(q,G)=4pi/|q+G|
548 : ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
549 13536 : do ig=1,npwc
550 13536 : vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q)) = Vcp%vc_sqrt(ig,iq_ibz)
551 : end do
552 :
553 264 : call timab(443,2,tsec) ! csigme (initq)
554 :
555 : ! Sum over bands.
556 1768 : do ib=1,Sigp%nbnds
557 : ! Parallelism over spin
558 : ! This processor has this k-point but what about spin?
559 1504 : if (proc_distrb(ib,ik_bz,spin)/=Wfd%my_rank) CYCLE
560 :
561 : ! Skip empty state ib for HF, SEX, and COHSEX.
562 1088 : if (qp_occ(ib,ik_ibz,spin)<tol_empty) CYCLE
563 :
564 1088 : theta_mu_minus_e0i=fact_spin*qp_occ(ib,ik_ibz,spin)
565 :
566 1088 : call wfd%get_ur(ib,ik_ibz,spin,ur_sum)
567 :
568 1088 : if (Psps%usepaw==1) then
569 : ! Load cprj for point ksum, this spin or spinor and *THIS* band.
570 : ! TODO MG I could avoid doing this but I have to exchange spin and bands ???
571 : ! For sure there is a better way to do this!
572 0 : call wfd%get_cprj(ib,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
573 0 : call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
574 : end if
575 :
576 4032 : do jb=ib1,ib2
577 : ! Get all <k-q,ib,s|e^{-i(q+G).r}|s,jb,k>, at once.
578 : call rho_tw_g(nspinor,npwc,gwc_nfftot,ndat1,gwc_ngfft,1,use_padfft,igfftcg0,gw_gbound,&
579 : ur_sum ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, &
580 : wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,&
581 2944 : nspinor,rhotwg_ki(:,jb))
582 :
583 2944 : if (Psps%usepaw==1) then
584 : ! Add on-site contribution, projectors are already in BZ !TODO Recheck this!
585 0 : i2=jb; if (nspinor==2) i2=(2*jb-1)
586 0 : spad=(nspinor-1)
587 : call paw_rho_tw_g(cryst,Pwij_qg, npwc,nspinor,nspinor,Gsph_c%gvec,&
588 0 : Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),rhotwg_ki(:,jb))
589 : end if
590 :
591 : ! Multiply by the square root of the Coulomb term.
592 : ! In 3-D systems, the factor sqrt(4pi) is included)
593 6784 : do ii=1,nspinor
594 3840 : spad = (ii-1) * npwc
595 159616 : rhotwg_ki(spad+1:spad+npwc,jb) = rhotwg_ki(spad+1:spad+npwc,jb)*vc_sqrt_qbz(1:npwc)
596 : end do
597 :
598 : ! === Treat analytically the case q --> 0 ===
599 : ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma
600 : ! while the Colulomb term is integrated out
601 : ! * In the scalar case we have nonzero contribution only if ib==jb
602 : ! * For nspinor==2 evalute <ib,up|jb,up> and <ib,dwn|jb,dwn>,
603 : ! impose orthonormalization since npwwfn might be < npwvec.
604 4032 : if (ik_bz == jk_bz) then
605 144 : if (nspinor == 1) then
606 32 : rhotwg_ki(1, jb) = czero_gw
607 32 : if (ib==jb) rhotwg_ki(1, jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp)
608 : else
609 112 : npw_k = Wfd%npwarr(ik_ibz)
610 112 : rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero
611 112 : if (ib == jb) then
612 8 : ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, spin, wave_sum, msg) == 0, msg)
613 8 : cg_sum => wave_sum%ug
614 8 : ABI_CHECK(wfd%get_wave_ptr(jb, jk_ibz, spin, wave_jb, msg) == 0, msg)
615 8 : cg_jb => wave_jb%ug
616 8 : ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1)
617 8 : rhotwg_ki(1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
618 8 : ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1)
619 8 : rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
620 : ! PAW is missing
621 :
622 : !rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half)
623 : !rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half)
624 : end if
625 : end if
626 : end if
627 : end do !jb Got all matrix elements from minbnd up to maxbnd.
628 :
629 4296 : do kb=ib1,ib2
630 : ! Get the ket \Sigma|\phi_{k,kb}> according to the method.
631 155776 : rhotwgp(:) = rhotwg_ki(:,kb)
632 158720 : sigc_ket = czero_gw
633 :
634 : ! SEX part. TODO add check on theta_mu_minus_e0i
635 6784 : do ispinor=1,nspinor
636 3840 : spadc = (ispinor-1) * npwc
637 3840 : call XGEMV('N',npwc,npwc,cone_gw,epsm1%epsm1_qbz(:,:,1),npwc,rhotwgp(1+spadc:),1,czero_gw,sigsex,1)
638 :
639 156672 : sigsex(:)= -theta_mu_minus_e0i*sigsex(:)
640 :
641 10624 : do io=1,nomega_tot ! nomega==1 as SEX is energy independent.
642 160512 : sigc_ket(spadc+1:spadc+npwc,io) = sigsex(:)
643 : end do
644 : end do
645 :
646 : ! Loop over the non-zero row elements of this column.
647 : ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS.
648 : ! 2) If gwcalctyp>=20:
649 : ! * Only off-diagonal elements connecting states with same character.
650 : ! * Only the upper triangle if HF, SEX, or COHSEX.
651 7392 : do irow=1,Sigcij_tab(spin)%col(kb)%size1
652 2944 : jb = Sigcij_tab(spin)%col(kb)%bidx(irow)
653 158720 : rhotwg = rhotwg_ki(:,jb)
654 :
655 : ! Calculate <\phi_j|\Sigma_c|\phi_k>
656 : ! Different freqs according to method (AC or Perturbative), see nomega_sigc.
657 8576 : do iab=1,Sigp%nsig_ab
658 5632 : spadc1=spinor_padc(1,iab); spadc2=spinor_padc(2,iab)
659 14208 : do io=1,nomega_sigc
660 11264 : sigctmp(io,iab) = XDOTC(npwc,rhotwg(spadc1+1:),1,sigc_ket(spadc2+1:,io),1)
661 : end do
662 : end do
663 :
664 : ! TODO: save wf1swf2_g to avoid having to recalculate it at each q-point.
665 2944 : if (mod10==SIG_COHSEX) then
666 : ! Evaluate Static COH. TODO add spinor.
667 1920 : if (coh_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank) then
668 : ! COH term is done only once for each k-point.
669 : ! It does not depend on the index ib summed over.
670 368 : coh_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank
671 :
672 : #if 1
673 : call calc_wfwfg(ktabr(:,jk_ibz), jik, spinrot_kgw, & ! TODO why jk_ibz?
674 368 : & gwc_nfftot, nspinor, gwc_ngfft, wfr_bdgw(:,jb), wfr_bdgw(:,kb), wf1swf2_g)
675 : #else
676 : ABI_CHECK(jik==1,"jik")
677 : call calc_wfwfg(ktabr(:,jk_bz), jik, spinrot_kgw, &
678 : gwc_nfftot, nspinor, gwc_ngfft, wfr_bdgw(:,jb), wfr_bdgw(:,kb), wf1swf2_g)
679 : #endif
680 :
681 368 : if (Psps%usepaw==1) then
682 0 : i1=jb; i2=kb
683 0 : if (nspinor==2) then
684 0 : i1=(2*jb-1); i2=(2*kb-1)
685 : end if
686 0 : spad=(nspinor-1)
687 : call paw_rho_tw_g(cryst,Pwij_fft,gwc_nfftot,Sigp%nsig_ab,nspinor,&
688 0 : gw_gfft,Cprj_kgw(:,i1:i1+spad),Cprj_kgw(:,i2:i2+spad),wf1swf2_g)
689 : end if
690 :
691 : call calc_coh(nspinor,Sigp%nsig_ab,gwc_nfftot,gwc_ngfft,npwc,Gsph_c%gvec,wf1swf2_g,epsm1%epsm1_qbz(:,:,1),&
692 368 : vc_sqrt_qbz,Vcp%i_sz,iq_ibz,(jb==kb),sigcohme)
693 :
694 736 : do io=1,nomega_sigc ! Should be 1
695 1440 : sigctmp(io,:) = sigctmp(io,:)+sigcohme(:)
696 : end do
697 :
698 : end if
699 : end if ! COHSEX
700 :
701 : ! Accumulate and, in case, symmetrize matrix elements of Sigma_c.
702 11520 : do iab=1,Sigp%nsig_ab
703 5632 : is_idx = spin; if (nspinor==2) is_idx=iab
704 :
705 : sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + &
706 11264 : (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab))
707 :
708 11264 : sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp* sigctmp(:,iab)
709 14208 : sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab))
710 : ! TODO this should be the contribution coming from the anti-hermitian part.
711 : end do
712 : end do !jb used to calculate matrix elements of $\Sigma$
713 :
714 : end do !kb to calculate matrix elements of $\Sigma$
715 : end do !ib
716 :
717 264 : if (print_time) then
718 69 : write(msg,'(3(a,i0))')' cohsex: ik_bz: ',ik_bz,'/',Kmesh%nbz,", spin: ",spin
719 69 : call cwtime_report(msg, cpu_k, wall_k, gflops_k); if (ik_bz == LOG_MODK) call wrtout(std_out, " ...")
720 : end if
721 :
722 : ! Deallocate k-dependent quantities.
723 264 : ABI_FREE(gw_gbound)
724 269 : if (Psps%usepaw==1) then
725 0 : call pawpwij_free(Pwij_qg)
726 0 : ABI_FREE(Pwij_qg)
727 : end if
728 : end do !ik_bz
729 :
730 5 : ABI_FREE(wfr_bdgw)
731 10 : if (Wfd%usepaw==1) then
732 0 : call pawcprj_free(Cprj_kgw)
733 0 : ABI_FREE(Cprj_kgw)
734 : end if
735 : end do !spin
736 :
737 : ! If epsm1 is MPI-shared, we have to close the RMA epoch.
738 5 : if (epsm1%use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1%epsm1_win, ierr)
739 5 : call epsm1%free_epsm1_qbz()
740 :
741 5 : ABI_FREE(igfftcg0)
742 :
743 : ! Gather contributions from all the CPUs.
744 5 : call xmpi_sum(sigcme_tmp, wfd%comm, ierr)
745 5 : call xmpi_sum(sigc, wfd%comm, ierr)
746 :
747 : ! Multiply by constants
748 : ! For 3D systems sqrt(4pi) is included in vc_sqrt_qbz.
749 1677 : sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz)
750 3277 : sigc = sigc /(Cryst%ucvol*Kmesh%nbz)
751 :
752 : ! If we have summed over the IBZ_q now we have to average over degenerate states.
753 : ! Presently only diagonal terms are considered
754 : ! TODO it does not work if nspinor==2.
755 10 : do spin=1,nsppol
756 10 : if (can_symmetrize(spin)) then
757 0 : ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, sigp%nsig_ab))
758 0 : sym_cme=czero
759 :
760 : ! Average over degenerate diagonal elements.
761 : ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results.
762 : ! another good reason to use a strict criterion for the tolerance on eigenvalues.
763 0 : do ib=ib1,ib2
764 0 : ndegs=0
765 0 : do jb=ib1,ib2
766 0 : if (degtab(ib,jb,spin)==1) then
767 0 : if (nspinor == 1) then
768 0 : sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:, :, jb, jb, spin), dim=1)
769 : else
770 0 : do ii=1,sigp%nsig_ab
771 0 : sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:, :, jb, jb, ii), dim=1)
772 : end do
773 : end if
774 :
775 : end if
776 0 : ndegs = ndegs + degtab(ib,jb,spin)
777 : end do
778 0 : sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs
779 : end do
780 :
781 0 : if (Sigp%gwcalctyp >= 20) then
782 0 : call QP_sym(spin)%symmetrize_mels(ib1,ib2,sigc(:,1,:,:,spin),sym_cme(1,:,:,1))
783 : end if
784 :
785 : ! Copy symmetrized values.
786 0 : do ib=ib1,ib2
787 0 : do jb=ib1,ib2
788 0 : if (nspinor == 1) then
789 0 : sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1)
790 : else
791 0 : sigcme_tmp(:,ib,jb,:) = sym_cme(:,ib,jb,:)
792 : end if
793 : end do
794 : end do
795 0 : ABI_FREE(sym_cme)
796 : end if
797 : end do
798 :
799 : ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX)
800 5 : if (Sigp%gwcalctyp>=20 .and. sigp%is_herm() ) then
801 0 : ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!")
802 0 : do spin=1,nsppol
803 0 : do io=1,nomega_sigc
804 0 : call hermitianize(sigcme_tmp(io,:,:,spin),"Upper")
805 : end do
806 : end do
807 : end if
808 :
809 : ! ===========================
810 : ! ==== Deallocate memory ====
811 : ! ===========================
812 5 : if (Psps%usepaw==1) then
813 0 : ABI_SFREE(gw_gfft)
814 0 : call pawcprj_free(Cprj_ksum)
815 0 : ABI_FREE(Cprj_ksum)
816 0 : if (allocated(Pwij_fft)) then
817 0 : call pawpwij_free(Pwij_fft)
818 0 : ABI_FREE(Pwij_fft)
819 : end if
820 : end if
821 :
822 5 : ABI_FREE(ktabr)
823 5 : ABI_FREE(ur_sum)
824 5 : ABI_FREE(rhotwg_ki)
825 5 : ABI_FREE(rhotwg)
826 5 : ABI_FREE(rhotwgp)
827 5 : ABI_FREE(vc_sqrt_qbz)
828 5 : ABI_FREE(sigc_ket)
829 : !ABI_FREE(epsm1_qbz)
830 5 : ABI_FREE(sigctmp)
831 5 : ABI_FREE(sigc)
832 5 : ABI_FREE(sigsex)
833 5 : ABI_FREE(proc_distrb)
834 5 : ABI_SFREE(wf1swf2_g)
835 5 : ABI_SFREE(coh_distrb)
836 5 : ABI_SFREE(degtab)
837 :
838 5 : call timab(495,2,tsec) ! csigme(SigC)
839 5 : call timab(491,2,tsec)
840 5 : call timab(423,2,tsec) ! cohsex_me
841 5 : call cwtime_report(" cohsex_me", cpu_all, wall_all, gflops_all)
842 :
843 : DBG_EXIT("COLL")
844 :
845 20 : end subroutine cohsex_me
846 : !!***
847 :
848 : !!****f* ABINIT/calc_coh
849 : !! NAME
850 : !! calc_coh
851 : !!
852 : !! FUNCTION
853 : !! Calculates the partial contribution to the COH part of the COHSEX self-energy for a given q-point.
854 : !!
855 : !! INPUTS
856 : !! iqibz=index of the irreducible q-point in the array qibz, point which is
857 : !! related by a symmetry operation to the point q summed over (see csigme).
858 : !! This index is also used to treat the integrable coulombian singularity at q=0
859 : !! ngfft(18)=contain all needed information about 3D FFT for GW wavefuntions,
860 : !! see ~abinit/doc/variables/vargs.htm#ngfft
861 : !! nsig_ab=Number of components in the self-energy operator (1 for collinear magnetism)
862 : !! npwc=number of plane waves in $\tilde epsilon^{-1}$
863 : !! nspinor=Number of spinorial components.
864 : !! nfftot=number of points in real space
865 : !! i_sz=contribution arising from the integrable coulombian singularity at q==0
866 : !! (see csigme for the method used), note that in case of 3-D systems the factor
867 : !! 4pi in the coulombian potential is included in the definition of i_sz
868 : !! gvec(3,npwc)=G vectors in reduced coordinates
869 : !! vc_sqrt(npwc)= square root of the coulombian matrix elements for this q-point
870 : !! epsm1q_o(npwc,npwc)= contains $\tilde epsilon^{-1}(q,w=0) - \delta_{G Gp}$ for
871 : !! the particular q-point considered in the sum
872 : !! wfg2_jk(nsig_ab*nfftot)= Fourier Transform of $\u_{jb k}^*(r) u_{kb k}$
873 : !! jb,kb=left and righ band indices definining the left and right states where the
874 : !! partial contribution to the matrix element of $\Sigma_{COH}$ is evaluated
875 : !!
876 : !! OUTPUT
877 : !! sigcohme=partial contribution to the matrix element of $<jb k \sigma|\Sigma_{COH} | kb k \sigma>$
878 : !! coming from this single q-point
879 : !!
880 : !! SOURCE
881 :
882 368 : subroutine calc_coh(nspinor,nsig_ab,nfftot,ngfft,npwc,gvec,wfg2_jk,epsm1q_o,vc_sqrt,i_sz,iqibz,same_band,sigcohme)
883 :
884 : !Arguments ------------------------------------
885 : !scalars
886 : integer,intent(in) :: iqibz,nfftot,npwc,nsig_ab,nspinor
887 : real(dp),intent(in) :: i_sz
888 : logical,intent(in) :: same_band
889 : !arrays
890 : integer,intent(in) :: gvec(3,npwc),ngfft(18)
891 : complex(gwp),intent(in) :: epsm1q_o(npwc,npwc),vc_sqrt(npwc)
892 : complex(gwp),intent(in) :: wfg2_jk(nfftot*nsig_ab)
893 : complex(gwp),intent(out) :: sigcohme(nsig_ab)
894 :
895 : !Local variables-------------------------------
896 : !scalars
897 : integer,save :: enough=0
898 : integer :: ig,ig4,ig4x,ig4y,ig4z,igp,igmin,ispinor,spad,outofbox
899 : !arrays
900 : integer :: g2mg1(3)
901 : ! *************************************************************************
902 :
903 : DBG_ENTER("COLL")
904 :
905 : ! === Partial contribution to the matrix element of Sigma_c ===
906 : ! * For nspinor==2, the closure relation reads:
907 : ! $\sum_s \psi_a^*(1)\psi_b(2) = \delta_{ab} \delta(1-2)$
908 : ! where a,b are the spinor components. As a consequence, Sigma_{COH} is always
909 : ! diagonal in spin-space and only diagonal matrix elements have to be calculated.
910 : ! MG TODO wfg2_jk should be calculated on an augmented FFT box to avoid spurious wrapping of G1-G2.
911 : ! MG: One has to make sure G1-G2 is still in the FFT mesh for each G1 and G2 in chi0 (not always true)
912 : ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic!
913 :
914 : ! * Treat the case q --> 0 adequately.
915 : ! TODO Better treatment of wings, check cutoff in the coulombian interaction.
916 368 : igmin=1; if (iqibz==1) igmin=2
917 :
918 1072 : sigcohme(:)=czero_gw
919 :
920 848 : do ispinor=1,nspinor
921 480 : spad=(ispinor-1)*nfftot
922 480 : outofbox=0
923 :
924 19552 : do igp=igmin,npwc
925 846816 : do ig=igmin,npwc
926 :
927 3309056 : g2mg1 = gvec(:,igp)-gvec(:,ig)
928 6618112 : if (ANY(g2mg1(:)>ngfft(1:3)/2) .or. ANY(g2mg1(:)<-(ngfft(1:3)-1)/2)) then
929 0 : outofbox = outofbox+1; CYCLE
930 : end if
931 :
932 827264 : ig4x=MODULO(g2mg1(1),ngfft(1))
933 827264 : ig4y=MODULO(g2mg1(2),ngfft(2))
934 827264 : ig4z=MODULO(g2mg1(3),ngfft(3))
935 827264 : ig4= 1+ig4x+ig4y*ngfft(1)+ig4z*ngfft(1)*ngfft(2)
936 :
937 846336 : sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*epsm1q_o(ig,igp)*vc_sqrt(ig)*vc_sqrt(igp)
938 : end do !ig
939 : end do !igp
940 :
941 848 : if (iqibz ==1 .and. same_band) sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*epsm1q_o(1,1)*i_sz
942 : end do !ispinor
943 :
944 368 : if (outofbox/=0) then
945 0 : enough=enough+1
946 0 : if (enough<=50) then
947 0 : ABI_WARNING(sjoin(' Number of G1-G2 pairs outside the G-sphere for Wfns:', itoa(outofbox)))
948 0 : if (enough==50) then
949 0 : call wrtout(std_out,' ========== Stop writing Warnings ==========')
950 : end if
951 : end if
952 : end if
953 :
954 : DBG_EXIT("COLL")
955 :
956 368 : end subroutine calc_coh
957 : !!***
958 :
959 : end module m_cohsex
960 : !!***
|
