Branch data 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:Gygi1986]] 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, have_spin
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 : : ! With MPI-shared epsm1, epsm1%rotate_iqbz contains collective fences on epsm1_qbz_win and the (single)
464 : : ! frequency is rotated by only one proc of shared_comm. Hence all procs must enter the k-loop for the same
465 : : ! q-points, even if they have no occupied band to sum for this (k, spin). Procs without bands then skip the
466 : : ! band loop. Procs only store the GW states of their own spin, hence they load them only if they have bands
467 : : ! to sum for this spin.
468 [ + + ]: 10 : do spin=1,nsppol
469 [ - + ]: 5 : if (epsm1%use_mpi_shared_win) then
470 [ # # # # : 0 : if (ALL(proc_distrb(:,:,spin) == xmpi_undefined_rank)) CYCLE
# # # # ]
471 : : else
472 [ + - + - : 5 : if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE
- + + - ]
473 : : end if
474 [ + - + - : 5 : have_spin = ANY(proc_distrb(:,:,spin) == Wfd%my_rank)
- + ]
475 : :
476 [ - + + - : 20 : ABI_MALLOC(wfr_bdgw,(gwc_nfftot*nspinor,ib1:ib2))
- + + - +
- + - - +
- + ]
477 [ + - ]: 5 : if (have_spin) then
478 [ - + + + : 49 : call wfd%get_many_ur([(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw)
+ - + + ]
479 : : else
480 [ # # # # ]: 0 : wfr_bdgw = czero_gw
481 : : end if
482 : :
483 [ - + ]: 5 : if (Wfd%usepaw==1) then
484 : : ! Load cprj for GW states, note the indexing.
485 : 0 : dimcprj_gw=nspinor*(ib2-ib1+1)
486 [ # # # # : 0 : ABI_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1))
# # # # #
# # # # #
# # # # #
# # # ]
487 : 0 : call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm)
488 : 0 : ibsp=ib1
489 [ # # ]: 0 : if (have_spin) then
490 [ # # ]: 0 : do jb=ib1,ib2
491 : 0 : call wfd%get_cprj(jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
492 : 0 : call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
493 : 0 : call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
494 : 0 : ibsp=ibsp+nspinor
495 : : end do
496 : : end if
497 : : end if
498 : :
499 [ + + ]: 269 : do ik_bz=1,Kmesh%nbz
500 : : ! Parallelization over k-points and spin
501 : : ! For the spin there is another check in the inner loop.
502 [ - + ]: 264 : if (epsm1%use_mpi_shared_win) then
503 [ # # # # : 0 : if (ALL(proc_distrb(:,ik_bz,spin) == xmpi_undefined_rank)) CYCLE
# # ]
504 : : else
505 [ + - - + : 264 : if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE
- + ]
506 : : end if
507 : :
508 : 264 : call timab(443,1,tsec) ! csigme (initq)
509 [ + - + + : 264 : print_time = wfd%my_rank == 0 .and. (ik_bz < LOG_MODK .or. mod(ik_bz, LOG_MODK) == 0)
+ + ]
510 : 69 : if (print_time) call cwtime(cpu_k, wall_k, gflops_k, "start")
511 : :
512 : : ! Find the corresponding irreducible k-point
513 : 264 : call kmesh%get_BZ_item(ik_bz,ksum,ik_ibz,isym_ki,iik,ph_mkt)
514 [ + + ]: 1320 : spinrot_kbz(:)=Cryst%spinrot(:,isym_ki)
515 : :
516 : : ! Identify q and G0 where q+G0=k_GW-k_i
517 [ + + ]: 1056 : kgw_m_ksum=kgw-ksum
518 : 264 : call findqg0(iq_bz,g0,kgw_m_ksum,Qmesh%nbz,Qmesh%bz,Sigp%mG0)
519 : :
520 : : ! Symmetrize the matrix elements.
521 : : ! Sum only q"s in IBZ_k. In this case elements are weighted
522 : : ! according to wtqp and wtqm. wtqm is for time-reversal.
523 : 264 : wtqp=1; wtqm=0
524 [ - + ]: 264 : if (can_symmetrize(spin)) then
525 [ # # ]: 0 : if (Ltg_k%ibzq(iq_bz)/=1) CYCLE
526 : 0 : wtqp=0; wtqm=0
527 [ # # ]: 0 : do isym=1,Ltg_k%nsym_sg
528 : 0 : wtqp=wtqp+Ltg_k%wtksym(1,isym,iq_bz)
529 : 0 : wtqm=wtqm+Ltg_k%wtksym(2,isym,iq_bz)
530 : : end do
531 : : end if
532 : :
533 : : !%write(msg,'(2(a,i4),a,i3)')' csigme : ik_bz ',ik_bz,'/',Kmesh%nbz,' done by processor ',Wfd%my_rank
534 : : !%call wrtout(std_out,msg)
535 : :
536 : : ! Find the corresponding irred q-point.
537 : 264 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
538 : 264 : q_is_gamma = (normv(qbz, Cryst%gmet, "G") < GW_TOLQ0)
539 : :
540 : : ! Tables for the FFT of the oscillators.
541 : : ! a) FFT index of the G-G0.
542 : : ! b) gw_gbound table for the zero-padded FFT performed in rhotwg.
543 [ + - - + : 1056 : ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2))
+ - + - +
- + - - +
- + - + ]
544 : 264 : call Gsph_c%fft_tabs(g0,gwc_mgfft,gwc_ngfft,use_padfft,gw_gbound,igfftcg0)
545 [ - + ]: 264 : if ( ANY(gwc_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
546 [ - + ]: 264 : if (use_padfft==0) then
547 [ # # ]: 0 : ABI_FREE(gw_gbound)
548 [ # # # # : 0 : ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2*use_padfft))
# # # # ]
549 : : end if
550 : :
551 [ - + ]: 264 : if (Psps%usepaw==1) then
552 : : ! Get PAW oscillator matrix elements $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form.
553 [ # # # # : 0 : ABI_MALLOC(Pwij_qg,(Psps%ntypat))
# # # # #
# # # # #
# # ]
554 : 0 : q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
555 : 0 : call pawpwij_init(Pwij_qg, npwc, q0, Gsph_c%gvec, Cryst%rprimd, Psps, Pawtab, Paw_pwff)
556 : : end if
557 : :
558 [ - + ]: 264 : if (epsm1%mqmem==0) then
559 : : ! Read q-slice of epsilon^{-1}|chi0 in epsm1%epsm1(:,:,:,1) (much slower but less memory).
560 : 0 : call epsm1%get_epsm1(Vcp,0,0,iomode,xmpi_comm_self,iqibzA=iq_ibz)
561 : : end if
562 : :
563 : : ! Only omega==0 for SEX or COHSEX
564 : 264 : call epsm1%rotate_iqbz(iq_bz, 1, npwc, Gsph_c, Qmesh, .True.) !, epsm1_qbz)
565 : :
566 : : ! Get Fourier components of the Coulomb interaction in the BZ.
567 : : ! In 3D systems, neglecting umklapp, vc(Sq,sG)=vc(q,G)=4pi/|q+G|
568 : : ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
569 [ + + ]: 13536 : do ig=1,npwc
570 : 13536 : vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q)) = Vcp%vc_sqrt(ig,iq_ibz)
571 : : end do
572 : :
573 : 264 : call timab(443,2,tsec) ! csigme (initq)
574 : :
575 : : ! Sum over bands.
576 [ + + ]: 1768 : do ib=1,Sigp%nbnds
577 : : ! Parallelism over spin
578 : : ! This processor has this k-point but what about spin?
579 [ + + ]: 1504 : if (proc_distrb(ib,ik_bz,spin)/=Wfd%my_rank) CYCLE
580 : :
581 : : ! Skip empty state ib for HF, SEX, and COHSEX.
582 [ - + ]: 1088 : if (qp_occ(ib,ik_ibz,spin)<tol_empty) CYCLE
583 : :
584 : 1088 : theta_mu_minus_e0i=fact_spin*qp_occ(ib,ik_ibz,spin)
585 : :
586 : 1088 : call wfd%get_ur(ib,ik_ibz,spin,ur_sum)
587 : :
588 [ - + ]: 1088 : if (Psps%usepaw==1) then
589 : : ! Load cprj for point ksum, this spin or spinor and *THIS* band.
590 : : ! TODO MG I could avoid doing this but I have to exchange spin and bands ???
591 : : ! For sure there is a better way to do this!
592 : 0 : call wfd%get_cprj(ib,ik_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
593 : 0 : call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
594 : : end if
595 : :
596 [ + + ]: 4032 : do jb=ib1,ib2
597 : : ! Get all <k-q,ib,s|e^{-i(q+G).r}|s,jb,k>, at once.
598 : : call rho_tw_g(nspinor,npwc,gwc_nfftot,ndat1,gwc_ngfft,1,use_padfft,igfftcg0,gw_gbound,&
599 : : ur_sum ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, &
600 : : wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,&
601 : 2944 : nspinor,rhotwg_ki(:,jb))
602 : :
603 [ - + ]: 2944 : if (Psps%usepaw==1) then
604 : : ! Add on-site contribution, projectors are already in BZ !TODO Recheck this!
605 [ # # ]: 0 : i2=jb; if (nspinor==2) i2=(2*jb-1)
606 : 0 : spad=(nspinor-1)
607 : : call paw_rho_tw_g(cryst,Pwij_qg, npwc,nspinor,nspinor,Gsph_c%gvec,&
608 : 0 : Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),rhotwg_ki(:,jb))
609 : : end if
610 : :
611 : : ! Multiply by the square root of the Coulomb term.
612 : : ! In 3-D systems, the factor sqrt(4pi) is included)
613 [ + + ]: 6784 : do ii=1,nspinor
614 : 3840 : spad = (ii-1) * npwc
615 [ + + ]: 159616 : rhotwg_ki(spad+1:spad+npwc,jb) = rhotwg_ki(spad+1:spad+npwc,jb)*vc_sqrt_qbz(1:npwc)
616 : : end do
617 : :
618 : : ! === Treat analytically the case q --> 0 ===
619 : : ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma
620 : : ! while the Colulomb term is integrated out
621 : : ! * In the scalar case we have nonzero contribution only if ib==jb
622 : : ! * For nspinor==2 evalute <ib,up|jb,up> and <ib,dwn|jb,dwn>,
623 : : ! impose orthonormalization since npwwfn might be < npwvec.
624 [ + + ]: 4032 : if (ik_bz == jk_bz) then
625 [ + + ]: 144 : if (nspinor == 1) then
626 : 32 : rhotwg_ki(1, jb) = czero_gw
627 [ + + ]: 32 : if (ib==jb) rhotwg_ki(1, jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp)
628 : : else
629 : 112 : npw_k = Wfd%npwarr(ik_ibz)
630 : 112 : rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero
631 [ + + ]: 112 : if (ib == jb) then
632 [ - + ]: 8 : ABI_CHECK(wfd%get_wave_ptr(ib, ik_ibz, spin, wave_sum, msg) == 0, msg)
633 : 8 : cg_sum => wave_sum%ug
634 [ - + ]: 8 : ABI_CHECK(wfd%get_wave_ptr(jb, jk_ibz, spin, wave_jb, msg) == 0, msg)
635 : 8 : cg_jb => wave_jb%ug
636 : 8 : ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1)
637 : 8 : rhotwg_ki(1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
638 : 8 : ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1)
639 : 8 : rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
640 : : ! PAW is missing
641 : :
642 : : !rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half)
643 : : !rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * sqrt(half)
644 : : end if
645 : : end if
646 : : end if
647 : : end do !jb Got all matrix elements from minbnd up to maxbnd.
648 : :
649 [ + + ]: 4296 : do kb=ib1,ib2
650 : : ! Get the ket \Sigma|\phi_{k,kb}> according to the method.
651 [ + + ]: 155776 : rhotwgp(:) = rhotwg_ki(:,kb)
652 [ + + + + ]: 158720 : sigc_ket = czero_gw
653 : :
654 : : ! SEX part. TODO add check on theta_mu_minus_e0i
655 [ + + ]: 6784 : do ispinor=1,nspinor
656 : 3840 : spadc = (ispinor-1) * npwc
657 [ - + - - : 3840 : call XGEMV('N',npwc,npwc,cone_gw,epsm1%epsm1_qbz(:,:,1),npwc,rhotwgp(1+spadc:),1,czero_gw,sigsex,1)
- - - - -
+ ]
658 : :
659 [ + + ]: 156672 : sigsex(:)= -theta_mu_minus_e0i*sigsex(:)
660 : :
661 [ + + ]: 10624 : do io=1,nomega_tot ! nomega==1 as SEX is energy independent.
662 [ + + ]: 160512 : sigc_ket(spadc+1:spadc+npwc,io) = sigsex(:)
663 : : end do
664 : : end do
665 : :
666 : : ! Loop over the non-zero row elements of this column.
667 : : ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS.
668 : : ! 2) If gwcalctyp>=20:
669 : : ! * Only off-diagonal elements connecting states with same character.
670 : : ! * Only the upper triangle if HF, SEX, or COHSEX.
671 [ + + ]: 7392 : do irow=1,Sigcij_tab(spin)%col(kb)%size1
672 : 2944 : jb = Sigcij_tab(spin)%col(kb)%bidx(irow)
673 [ - + + + ]: 158720 : rhotwg = rhotwg_ki(:,jb)
674 : :
675 : : ! Calculate <\phi_j|\Sigma_c|\phi_k>
676 : : ! Different freqs according to method (AC or Perturbative), see nomega_sigc.
677 [ + + ]: 8576 : do iab=1,Sigp%nsig_ab
678 : 5632 : spadc1=spinor_padc(1,iab); spadc2=spinor_padc(2,iab)
679 [ + + ]: 14208 : do io=1,nomega_sigc
680 : 11264 : sigctmp(io,iab) = XDOTC(npwc,rhotwg(spadc1+1:),1,sigc_ket(spadc2+1:,io),1)
681 : : end do
682 : : end do
683 : :
684 : : ! TODO: save wf1swf2_g to avoid having to recalculate it at each q-point.
685 [ + + ]: 2944 : if (mod10==SIG_COHSEX) then
686 : : ! Evaluate Static COH. TODO add spinor.
687 [ + + ]: 1920 : if (coh_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank) then
688 : : ! COH term is done only once for each k-point.
689 : : ! It does not depend on the index ib summed over.
690 : 368 : coh_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank
691 : :
692 : : #if 1
693 : : call calc_wfwfg(ktabr(:,jk_ibz), jik, spinrot_kgw, & ! TODO why jk_ibz?
694 : 368 : & gwc_nfftot, nspinor, gwc_ngfft, wfr_bdgw(:,jb), wfr_bdgw(:,kb), wf1swf2_g)
695 : : #else
696 : : ABI_CHECK(jik==1,"jik")
697 : : call calc_wfwfg(ktabr(:,jk_bz), jik, spinrot_kgw, &
698 : : gwc_nfftot, nspinor, gwc_ngfft, wfr_bdgw(:,jb), wfr_bdgw(:,kb), wf1swf2_g)
699 : : #endif
700 : :
701 [ - + ]: 368 : if (Psps%usepaw==1) then
702 : 0 : i1=jb; i2=kb
703 [ # # ]: 0 : if (nspinor==2) then
704 : 0 : i1=(2*jb-1); i2=(2*kb-1)
705 : : end if
706 : 0 : spad=(nspinor-1)
707 : : call paw_rho_tw_g(cryst,Pwij_fft,gwc_nfftot,Sigp%nsig_ab,nspinor,&
708 : 0 : gw_gfft,Cprj_kgw(:,i1:i1+spad),Cprj_kgw(:,i2:i2+spad),wf1swf2_g)
709 : : end if
710 : :
711 : : call calc_coh(nspinor,Sigp%nsig_ab,gwc_nfftot,gwc_ngfft,npwc,Gsph_c%gvec,wf1swf2_g,epsm1%epsm1_qbz(:,:,1),&
712 [ - + - - : 368 : vc_sqrt_qbz,Vcp%i_sz,iq_ibz,(jb==kb),sigcohme)
- - - - -
+ ]
713 : :
714 [ + + ]: 736 : do io=1,nomega_sigc ! Should be 1
715 [ + + ]: 1440 : sigctmp(io,:) = sigctmp(io,:)+sigcohme(:)
716 : : end do
717 : :
718 : : end if
719 : : end if ! COHSEX
720 : :
721 : : ! Accumulate and, in case, symmetrize matrix elements of Sigma_c.
722 [ + + ]: 11520 : do iab=1,Sigp%nsig_ab
723 [ + + ]: 5632 : is_idx = spin; if (nspinor==2) is_idx=iab
724 : :
725 : : sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + &
726 [ + + ]: 11264 : (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab))
727 : :
728 [ + + ]: 11264 : sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp* sigctmp(:,iab)
729 [ + + ]: 14208 : sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab))
730 : : ! TODO this should be the contribution coming from the anti-hermitian part.
731 : : end do
732 : : end do !jb used to calculate matrix elements of $\Sigma$
733 : :
734 : : end do !kb to calculate matrix elements of $\Sigma$
735 : : end do !ib
736 : :
737 [ + + ]: 264 : if (print_time) then
738 : 69 : write(msg,'(3(a,i0))')' cohsex: ik_bz: ',ik_bz,'/',Kmesh%nbz,", spin: ",spin
739 [ + + ]: 69 : call cwtime_report(msg, cpu_k, wall_k, gflops_k); if (ik_bz == LOG_MODK) call wrtout(std_out, " ...")
740 : : end if
741 : :
742 : : ! Deallocate k-dependent quantities.
743 [ - + ]: 264 : ABI_FREE(gw_gbound)
744 [ + - ]: 269 : if (Psps%usepaw==1) then
745 : 0 : call pawpwij_free(Pwij_qg)
746 [ # # # # : 0 : ABI_FREE(Pwij_qg)
# # # # ]
747 : : end if
748 : : end do !ik_bz
749 : :
750 [ - + ]: 5 : ABI_FREE(wfr_bdgw)
751 [ + - ]: 10 : if (Wfd%usepaw==1) then
752 : 0 : call pawcprj_free(Cprj_kgw)
753 [ # # # # : 0 : ABI_FREE(Cprj_kgw)
# # # # #
# ]
754 : : end if
755 : : end do !spin
756 : :
757 : : ! If epsm1 is MPI-shared, we have to close the RMA epoch.
758 [ - + ]: 5 : if (epsm1%use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1%epsm1_win, ierr)
759 : 5 : call epsm1%free_epsm1_qbz()
760 : :
761 [ - + ]: 5 : ABI_FREE(igfftcg0)
762 : :
763 : : ! Gather contributions from all the CPUs.
764 : 5 : call xmpi_sum(sigcme_tmp, wfd%comm, ierr)
765 : 5 : call xmpi_sum(sigc, wfd%comm, ierr)
766 : :
767 : : ! Multiply by constants
768 : : ! For 3D systems sqrt(4pi) is included in vc_sqrt_qbz.
769 [ + + + + : 1677 : sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz)
+ + + + ]
770 [ + + + + : 3277 : sigc = sigc /(Cryst%ucvol*Kmesh%nbz)
+ + + + +
+ ]
771 : :
772 : : ! If we have summed over the IBZ_q now we have to average over degenerate states.
773 : : ! Presently only diagonal terms are considered
774 : : ! TODO it does not work if nspinor==2.
775 [ + + ]: 10 : do spin=1,nsppol
776 [ - + ]: 10 : if (can_symmetrize(spin)) then
777 [ # # # # : 0 : ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, sigp%nsig_ab))
# # # # #
# # # # #
# # # # #
# # # #
# ]
778 [ # # # # : 0 : sym_cme=czero
# # # # ]
779 : :
780 : : ! Average over degenerate diagonal elements.
781 : : ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results.
782 : : ! another good reason to use a strict criterion for the tolerance on eigenvalues.
783 [ # # ]: 0 : do ib=ib1,ib2
784 : 0 : ndegs=0
785 [ # # ]: 0 : do jb=ib1,ib2
786 [ # # ]: 0 : if (degtab(ib,jb,spin)==1) then
787 [ # # ]: 0 : if (nspinor == 1) then
788 [ # # # # ]: 0 : sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:, :, jb, jb, spin), dim=1)
789 : : else
790 [ # # ]: 0 : do ii=1,sigp%nsig_ab
791 [ # # # # ]: 0 : sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:, :, jb, jb, ii), dim=1)
792 : : end do
793 : : end if
794 : :
795 : : end if
796 : 0 : ndegs = ndegs + degtab(ib,jb,spin)
797 : : end do
798 [ # # # # ]: 0 : sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs
799 : : end do
800 : :
801 [ # # ]: 0 : if (Sigp%gwcalctyp >= 20) then
802 [ # # # # : 0 : call QP_sym(spin)%symmetrize_mels(ib1,ib2,sigc(:,1,:,:,spin),sym_cme(1,:,:,1))
# # # # #
# # # # #
# # # # #
# # # ]
803 : : end if
804 : :
805 : : ! Copy symmetrized values.
806 [ # # ]: 0 : do ib=ib1,ib2
807 [ # # ]: 0 : do jb=ib1,ib2
808 [ # # ]: 0 : if (nspinor == 1) then
809 [ # # ]: 0 : sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1)
810 : : else
811 [ # # # # ]: 0 : sigcme_tmp(:,ib,jb,:) = sym_cme(:,ib,jb,:)
812 : : end if
813 : : end do
814 : : end do
815 [ # # ]: 0 : ABI_FREE(sym_cme)
816 : : end if
817 : : end do
818 : :
819 : : ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX)
820 [ - + - - ]: 5 : if (Sigp%gwcalctyp>=20 .and. sigp%is_herm() ) then
821 [ # # ]: 0 : ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!")
822 [ # # ]: 0 : do spin=1,nsppol
823 [ # # ]: 0 : do io=1,nomega_sigc
824 : 0 : call hermitianize(sigcme_tmp(io,:,:,spin),"Upper")
825 : : end do
826 : : end do
827 : : end if
828 : :
829 : : ! ===========================
830 : : ! ==== Deallocate memory ====
831 : : ! ===========================
832 [ - + ]: 5 : if (Psps%usepaw==1) then
833 [ # # ]: 0 : ABI_SFREE(gw_gfft)
834 : 0 : call pawcprj_free(Cprj_ksum)
835 [ # # # # : 0 : ABI_FREE(Cprj_ksum)
# # # # #
# ]
836 [ # # ]: 0 : if (allocated(Pwij_fft)) then
837 : 0 : call pawpwij_free(Pwij_fft)
838 [ # # # # : 0 : ABI_FREE(Pwij_fft)
# # # # ]
839 : : end if
840 : : end if
841 : :
842 [ - + ]: 5 : ABI_FREE(ktabr)
843 [ - + ]: 5 : ABI_FREE(ur_sum)
844 [ - + ]: 5 : ABI_FREE(rhotwg_ki)
845 [ - + ]: 5 : ABI_FREE(rhotwg)
846 [ - + ]: 5 : ABI_FREE(rhotwgp)
847 [ - + ]: 5 : ABI_FREE(vc_sqrt_qbz)
848 [ - + ]: 5 : ABI_FREE(sigc_ket)
849 : : !ABI_FREE(epsm1_qbz)
850 [ - + ]: 5 : ABI_FREE(sigctmp)
851 [ - + ]: 5 : ABI_FREE(sigc)
852 [ - + ]: 5 : ABI_FREE(sigsex)
853 [ - + ]: 5 : ABI_FREE(proc_distrb)
854 [ + + ]: 5 : ABI_SFREE(wf1swf2_g)
855 [ + + ]: 5 : ABI_SFREE(coh_distrb)
856 [ - + ]: 5 : ABI_SFREE(degtab)
857 : :
858 : 5 : call timab(495,2,tsec) ! csigme(SigC)
859 : 5 : call timab(491,2,tsec)
860 : 5 : call timab(423,2,tsec) ! cohsex_me
861 : 5 : call cwtime_report(" cohsex_me", cpu_all, wall_all, gflops_all)
862 : :
863 : : DBG_EXIT("COLL")
864 : :
865 [ - + - - : 20 : end subroutine cohsex_me
- - - + -
- - - - +
- + - + -
- - - - -
- + - - -
- - - ]
866 : : !!***
867 : :
868 : : !!****f* ABINIT/calc_coh
869 : : !! NAME
870 : : !! calc_coh
871 : : !!
872 : : !! FUNCTION
873 : : !! Calculates the partial contribution to the COH part of the COHSEX self-energy for a given q-point.
874 : : !!
875 : : !! INPUTS
876 : : !! iqibz=index of the irreducible q-point in the array qibz, point which is
877 : : !! related by a symmetry operation to the point q summed over (see csigme).
878 : : !! This index is also used to treat the integrable coulombian singularity at q=0
879 : : !! ngfft(18)=contain all needed information about 3D FFT for GW wavefuntions,
880 : : !! see ~abinit/doc/variables/vargs.htm#ngfft
881 : : !! nsig_ab=Number of components in the self-energy operator (1 for collinear magnetism)
882 : : !! npwc=number of plane waves in $\tilde epsilon^{-1}$
883 : : !! nspinor=Number of spinorial components.
884 : : !! nfftot=number of points in real space
885 : : !! i_sz=contribution arising from the integrable coulombian singularity at q==0
886 : : !! (see csigme for the method used), note that in case of 3-D systems the factor
887 : : !! 4pi in the coulombian potential is included in the definition of i_sz
888 : : !! gvec(3,npwc)=G vectors in reduced coordinates
889 : : !! vc_sqrt(npwc)= square root of the coulombian matrix elements for this q-point
890 : : !! epsm1q_o(npwc,npwc)= contains $\tilde epsilon^{-1}(q,w=0) - \delta_{G Gp}$ for
891 : : !! the particular q-point considered in the sum
892 : : !! wfg2_jk(nsig_ab*nfftot)= Fourier Transform of $\u_{jb k}^*(r) u_{kb k}$
893 : : !! jb,kb=left and righ band indices definining the left and right states where the
894 : : !! partial contribution to the matrix element of $\Sigma_{COH}$ is evaluated
895 : : !!
896 : : !! OUTPUT
897 : : !! sigcohme=partial contribution to the matrix element of $<jb k \sigma|\Sigma_{COH} | kb k \sigma>$
898 : : !! coming from this single q-point
899 : : !!
900 : : !! SOURCE
901 : :
902 : 368 : subroutine calc_coh(nspinor,nsig_ab,nfftot,ngfft,npwc,gvec,wfg2_jk,epsm1q_o,vc_sqrt,i_sz,iqibz,same_band,sigcohme)
903 : :
904 : : !Arguments ------------------------------------
905 : : !scalars
906 : : integer,intent(in) :: iqibz,nfftot,npwc,nsig_ab,nspinor
907 : : real(dp),intent(in) :: i_sz
908 : : logical,intent(in) :: same_band
909 : : !arrays
910 : : integer,intent(in) :: gvec(3,npwc),ngfft(18)
911 : : complex(gwp),intent(in) :: epsm1q_o(npwc,npwc),vc_sqrt(npwc)
912 : : complex(gwp),intent(in) :: wfg2_jk(nfftot*nsig_ab)
913 : : complex(gwp),intent(out) :: sigcohme(nsig_ab)
914 : :
915 : : !Local variables-------------------------------
916 : : !scalars
917 : : integer,save :: enough=0
918 : : integer :: ig,ig4,ig4x,ig4y,ig4z,igp,igmin,ispinor,spad,outofbox
919 : : !arrays
920 : : integer :: g2mg1(3)
921 : : ! *************************************************************************
922 : :
923 : : DBG_ENTER("COLL")
924 : :
925 : : ! === Partial contribution to the matrix element of Sigma_c ===
926 : : ! * For nspinor==2, the closure relation reads:
927 : : ! $\sum_s \psi_a^*(1)\psi_b(2) = \delta_{ab} \delta(1-2)$
928 : : ! where a,b are the spinor components. As a consequence, Sigma_{COH} is always
929 : : ! diagonal in spin-space and only diagonal matrix elements have to be calculated.
930 : : ! MG TODO wfg2_jk should be calculated on an augmented FFT box to avoid spurious wrapping of G1-G2.
931 : : ! MG: One has to make sure G1-G2 is still in the FFT mesh for each G1 and G2 in chi0 (not always true)
932 : : ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic!
933 : :
934 : : ! * Treat the case q --> 0 adequately.
935 : : ! TODO Better treatment of wings, check cutoff in the coulombian interaction.
936 [ + + ]: 368 : igmin=1; if (iqibz==1) igmin=2
937 : :
938 [ + + ]: 1072 : sigcohme(:)=czero_gw
939 : :
940 [ + + ]: 848 : do ispinor=1,nspinor
941 : 480 : spad=(ispinor-1)*nfftot
942 : 480 : outofbox=0
943 : :
944 [ + + ]: 19552 : do igp=igmin,npwc
945 [ + + ]: 846816 : do ig=igmin,npwc
946 : :
947 [ + + ]: 3309056 : g2mg1 = gvec(:,igp)-gvec(:,ig)
948 [ + + + - : 6618112 : if (ANY(g2mg1(:)>ngfft(1:3)/2) .or. ANY(g2mg1(:)<-(ngfft(1:3)-1)/2)) then
+ + + - -
+ ]
949 : 0 : outofbox = outofbox+1; CYCLE
950 : : end if
951 : :
952 : 827264 : ig4x=MODULO(g2mg1(1),ngfft(1))
953 : 827264 : ig4y=MODULO(g2mg1(2),ngfft(2))
954 : 827264 : ig4z=MODULO(g2mg1(3),ngfft(3))
955 : 827264 : ig4= 1+ig4x+ig4y*ngfft(1)+ig4z*ngfft(1)*ngfft(2)
956 : :
957 : 846336 : sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*epsm1q_o(ig,igp)*vc_sqrt(ig)*vc_sqrt(igp)
958 : : end do !ig
959 : : end do !igp
960 : :
961 [ + + + - ]: 848 : if (iqibz ==1 .and. same_band) sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*epsm1q_o(1,1)*i_sz
962 : : end do !ispinor
963 : :
964 [ - + ]: 368 : if (outofbox/=0) then
965 : 0 : enough=enough+1
966 [ # # ]: 0 : if (enough<=50) then
967 : 0 : ABI_WARNING(sjoin(' Number of G1-G2 pairs outside the G-sphere for Wfns:', itoa(outofbox)))
968 [ # # ]: 0 : if (enough==50) then
969 : 0 : call wrtout(std_out,' ========== Stop writing Warnings ==========')
970 : : end if
971 : : end if
972 : : end if
973 : :
974 : : DBG_EXIT("COLL")
975 : :
976 : 368 : end subroutine calc_coh
977 : : !!***
978 : :
979 : : end module m_cohsex
980 : : !!***
|
