Line data Source code
1 : !!****m* ABINIT/m_chi0
2 : !! NAME
3 : !! m_chi0
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (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_chi0
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_xmpi
27 : use m_errors
28 : use m_hide_blas
29 : use m_time
30 : use m_wfd
31 : use m_dtset
32 :
33 : use defs_datatypes, only : pseudopotential_type
34 : use defs_abitypes, only : MPI_type
35 : use m_fstrings, only : ftoa, sjoin, itoa
36 : use m_gwdefs, only : GW_TOL_DOCC, GW_TOL_W0, czero_gw, em1params_t, g0g0w
37 : use m_numeric_tools, only : imin_loc, print_arr
38 : use m_geometry, only : normv, vdotw
39 : use m_crystal, only : crystal_t
40 : use m_fft_mesh, only : rotate_FFT_mesh, get_gfft
41 : use m_occ, only : getnel
42 : use m_ebands, only : ebands_t, pack_eneocc, unpack_eneocc
43 : use m_bz_mesh, only : kmesh_t, littlegroup_t
44 : use m_gsphere, only : gsphere_t
45 : use m_io_tools, only : flush_unit
46 : use m_oscillators, only : rho_tw_g, calc_wfwfg
47 : !use m_ddk, only : ddkop_t
48 : use m_vkbr, only : vkbr_t, vkbr_free, vkbr_init, nc_ihr_comm
49 : use m_chi0tk, only : hilbert_transform, setup_spectral, assemblychi0_sym, assemblychi0sf, symmetrize_afm_chi0, &
50 : approxdelta, completechi0_deltapart, accumulate_chi0sumrule, make_transitions, &
51 : chi0_bbp_mask, accumulate_chi0_q0, accumulate_sfchi0_q0, hilbert_transform_headwings
52 : use m_pawang, only : pawang_type
53 : use m_pawrad, only : pawrad_type
54 : use m_pawtab, only : pawtab_type
55 : use m_paw_ij, only : paw_ij_type
56 : use m_pawfgrtab, only : pawfgrtab_type
57 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
58 : use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
59 : use m_paw_sym, only : paw_symcprj
60 : use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t
61 : use m_paw_hr, only : pawhur_t, pawhur_free, pawhur_init, paw_ihr, paw_cross_ihr_comm
62 : use m_read_plowannier, only : read_plowannier
63 : use m_plowannier, only : plowannier_type
64 : use m_pstat, only : pstat_proc
65 :
66 : implicit none
67 :
68 : private
69 : !!***
70 :
71 : public :: cchi0q0
72 : public :: cchi0
73 : public :: chi0q0_intraband
74 : !!***
75 :
76 : integer,parameter :: LOG_MODK = 5
77 :
78 : contains
79 : !!***
80 :
81 : !!****f* ABINIT/cchi0q0
82 : !! NAME
83 : !! cchi0q0
84 : !!
85 : !! FUNCTION
86 : !! Calculate chi0 in the limit q --> 0
87 : !!
88 : !! INPUTS
89 : !! use_tr=If .TRUE. Wfs_val are allocate and only resonant transitions are evaluated (assumes time reversal symmetry)
90 : !! Dtset <type(dataset_type)>=all input variables in this dataset
91 : !! Ep= datatype gathering differening parameters related to the calculation of the inverse dielectric matrix
92 : !! Gsph_epsG0<gvectors_data_type>: Info on the G-sphere used to describe chi0/espilon (including umklapp)
93 : !! %ng=number of G vectors
94 : !! %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G in the array gvec
95 : !! %phmGt(ng,nsym)=phase factor e^{-iG.\tau} needed to symmetrize oscillator matrix elements and chi0
96 : !! %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
97 : !! %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
98 : !! Ep%inclvkb=flag to include (or not) the grad of Vkb
99 : !! Ltg_q= little group datatype
100 : !! nbvw=number of bands in the arrays wfrv,wfgv
101 : !! Kmesh<kmesh_t> The k-point mesh
102 : !! %kbz(3,nbz)=k-point coordinates, full Brillouin zone
103 : !! %tab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding
104 : !! irreducible point (kIBZ), where kBZ= (IS) kIBZ and I is either the inversion or the identity
105 : !! %tabi(nbzx)= for each point in the BZ defines whether inversion has to be
106 : !! considered in the relation kBZ=(IS) kIBZ (1 => only S; -1 => -S)
107 : !! %tabo(nbzx)= the symmetry operation S that takes kIBZ to each kBZ
108 : !! %tabp(nbzx)= phase factor associated to tnons e^{-i 2 \pi k\cdot R{^-1}t}
109 : !! ktabr(nfftot_gw,Kmesh%nbz) index of R^-(r-t) in the FFT array, where k_BZ = (IS) k_IBZ and S = \transpose R^{-1}
110 : !! Ep%nbnds=number of bands
111 : !! ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths.
112 : !! Ep%nomega=number of frequencies
113 : !! Cryst<crystal_t>= data type gathering info on symmetries and unit cell
114 : !! %natom=number of atoms
115 : !! %nsym=number of symmetry operations
116 : !! %symrec(3,3,nsym)=symmetry operations in reciprocal space
117 : !! %typat(natom)=type of each atom
118 : !! %xred(3,natom)=reduced coordinated of atoms
119 : !! %rprimd(3,3)=dimensional primitive translations in real space (bohr)
120 : !! %timrev=2 if time-reversal symmetry can be used, 1 otherwise
121 : !! Ep%npwe=number of planewaves for sigma exchange (input variable)
122 : !! nfftot_gw=Total number of points in the GW FFT grid
123 : !! Ep%omega(Ep%nomega)=frequencies
124 : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
125 : !! %mpsang=1+maximum angular momentum for nonlocal pseudopotential
126 : !! Pawang<pawang_type> angular mesh discretization and related data:
127 : !! Pawrad(ntypat*usepaw)<Pawrad_type>=paw radial mesh and related data
128 : !! Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels
129 : !! qp_ebands<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities)
130 : !! %mband=MAX number of bands over k-points and spin (==Ep%nbnds)
131 : !! %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band
132 : !! %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes
133 : !! ks_ebands<ebands_t>=KS energies and occupations.
134 : !! %eig(mband,nkpt,nsppol)=KS energies
135 : !! Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
136 : !! Wfd<wfdgw_t>=Object used to access the wavefunctions
137 : !!
138 : !! OUTPUT
139 : !! chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix for wavevector qq,
140 : !! and frequencies defined by Ep%omega
141 : !! chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)= Lower wings
142 : !! chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)= Upper wings
143 : !! chi0_head(3,3,Ep%nomega)=Head of chi0.
144 : !!
145 : !! NOTES
146 : !! *) The terms "head", "wings" and "body" of chi(G,Gp) refer to
147 : !! G=Gp=0, either G or Gp=0, and neither=0 respectively
148 : !!
149 : !! *) Symmetry conventions:
150 : !! 1) symmetry in real space is defined as: R_t f(r) = f(R^-1(r-t))
151 : !! 2) S=\transpose R^-1
152 : !! 3) kbz=S kibz
153 : !!
154 : !! The wavefunctions for the k-point in the BZ are (assuming nondegenerate states):
155 : !!
156 : !! u(G,b, Sk) = u ( S^-1G,b,k)* e^{-i(Sk+G)*t)
157 : !! u(G,b,-Sk) = u*(-S^-1G,b,k)* e^{ i(Sk-G)*t)
158 : !!
159 : !! u(r,b, Sk) = u (R^-1(r-t),b,k) e^{-iSk*t}
160 : !! u(r,b,-Sk) = u*(R^-1(r-t),b,k) e^{ iSK*t}
161 : !!
162 : !! The gradient of Vnl(K,Kp) for the k-point in the BZ should be:
163 : !!
164 : !! gradvnl(SG,SGp,Sk)=S gradvnl(G,Gp,kibz)
165 : !!
166 : !! TODO
167 : !! Check npwepG0 before activating umklapps
168 : !!
169 : !! SOURCE
170 :
171 146 : subroutine cchi0q0(use_tr,Dtset,Cryst,Ep,Psps,Kmesh,qp_ebands,ks_ebands,Gsph_epsG0,&
172 146 : Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,Pawfgrtab,Paw_onsite,ktabr,ktabrf,nbvw,ngfft_gw,&
173 146 : nfftot_gw,ngfftf,nfftf_tot,chi0,chi0_head,chi0_lwing,chi0_uwing,Ltg_q,chi0_sumrule,Wfd,Wfdf,wan) !,mpi_enreg
174 :
175 : !Arguments ------------------------------------
176 : !scalars
177 : integer,intent(in) :: nbvw,nfftot_gw,nfftf_tot
178 : logical,intent(in) :: use_tr
179 : type(ebands_t),target,intent(in) :: qp_ebands,ks_ebands
180 : type(crystal_t),intent(in) :: Cryst
181 : type(Dataset_type),intent(in) :: Dtset
182 : type(littlegroup_t),intent(in) :: Ltg_q
183 : type(em1params_t),intent(in) :: Ep
184 : type(kmesh_t),intent(in) :: Kmesh
185 : type(gsphere_t),intent(in) :: Gsph_epsG0
186 : type(Pseudopotential_type),intent(in) :: Psps
187 : type(Pawang_type),intent(in) :: Pawang
188 : type(wfdgw_t),target,intent(inout) :: Wfd,Wfdf
189 : ! type(MPI_type),intent(in) :: mpi_enreg
190 : !arrays
191 : integer,intent(in) :: ktabr(nfftot_gw,Kmesh%nbz),ktabrf(nfftf_tot*Dtset%pawcross,Kmesh%nbz)
192 : integer,intent(in) :: ngfft_gw(18),ngfftf(18)
193 : real(dp),intent(out) :: chi0_sumrule(Ep%npwe)
194 : complex(gwp),intent(out) :: chi0(Ep%npwe,Ep%npwe,Ep%nomega)
195 : complex(dp),intent(out) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)
196 : complex(dp),intent(out) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)
197 : complex(dp),intent(out) :: chi0_head(3,3,Ep%nomega)
198 : type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw)
199 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
200 : type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw)
201 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
202 : type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
203 : type(plowannier_type),intent(inout) :: wan
204 : type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
205 :
206 : !Local variables ------------------------------
207 : !scalars
208 : integer,parameter :: tim_fourdp=1, enough=10, two_poles=2, one_pole=1, ndat1=1
209 : integer :: bandinf,bandsup,lcor,nspinor,npw_k,istwf_k,mband,nfft,band1c,band2c, mpw
210 : integer :: band1,band2,iat1,iat2,iat,ig,ig1,ig2,itim_k,ik_bz,ik_ibz,io,iqlwl,ispinor1,ispinor2,isym_k,il1,il2
211 : integer :: itypatcor,m1,m2,nkpt_summed,dim_rtwg,use_padfft,gw_fftalga,use_padfftf,mgfftf
212 : integer :: my_nbbp,my_nbbpks,spin,nsppol,iq,nq
213 : integer :: comm,ierr,my_wl,my_wr,iomegal,iomegar,gw_mgfft,dummy
214 : real(dp) :: cpu_all, wall_all, gflops_all, cpu_k, wall_k, gflops_k
215 : real(dp) :: fac,fac1,fac2,fac3,fac4,spin_fact,deltaf_b1b2,weight,factor
216 : real(dp) :: max_rest,min_rest,my_max_rest,my_min_rest, qlen
217 : real(dp) :: en_high,deltaeGW_enhigh_b2,wl,wr,numerator,deltaeGW_b1b2,gw_gsq,memreq
218 : complex(dp) :: deltaeKS_b1b2
219 : logical :: qzero, luwindow, is_metallic, print_time !, use_ddk
220 : character(len=500) :: msg_tmp,msg,allup
221 146 : type(gsphere_t) :: Gsph_FFT
222 : type(wave_t),pointer :: wave1, wave2
223 : !type(ddkop_t) :: ddkop
224 : !arrays
225 146 : integer,contiguous, pointer :: kg_k(:,:)
226 292 : integer :: ucrpa_bands(2), got(Wfd%nproc)
227 292 : integer :: wtk_ltg(Kmesh%nbz)
228 146 : integer,allocatable :: tabr_k(:),tabrf_k(:), igffteps0(:),gspfft_igfft(:),igfftepsG0f(:)
229 146 : integer,allocatable :: gw_gfft(:,:),gw_gbound(:,:),dummy_gbound(:,:),gboundf(:,:), bbp_ks_distrb(:,:,:,:)
230 : real(dp) :: kbz(3),spinrot_kbz(4),q0(3)
231 146 : real(dp),contiguous, pointer :: ks_eig(:,:,:),qp_eig(:,:,:),qp_occ(:,:,:)
232 146 : real(dp),allocatable :: omegasf(:), qdirs(:,:)
233 292 : complex(gwp) :: rhotwx(3,Wfd%nspinor**2)
234 146 : complex(gwp),allocatable :: rhotwg(:)
235 146 : complex(dp),allocatable :: green_w(:),green_enhigh_w(:)
236 146 : complex(dp),allocatable :: sf_lwing(:,:,:),sf_uwing(:,:,:),sf_head(:,:,:)
237 : complex(dp) :: chq(3), wng(3)
238 : complex(dp) :: ph_mkt
239 146 : complex(dp),allocatable :: kkweight(:,:)
240 146 : complex(gwp),allocatable :: ur1_kibz(:),ur2_kibz(:), usr1_k(:),ur2_k(:), wfwfg(:), sf_chi0(:,:,:)
241 146 : complex(gwp),allocatable :: ur_ae1(:),ur_ae_onsite1(:),ur_ps_onsite1(:)
242 146 : complex(gwp),allocatable :: ur_ae2(:),ur_ae_onsite2(:),ur_ps_onsite2(:)
243 146 : complex(gwp),contiguous, pointer :: ug1(:),ug2(:)
244 146 : complex(dp), allocatable :: coeffW_BZ(:,:,:,:,:,:), head_qvals(:)
245 292 : logical :: gradk_not_done(Kmesh%nibz)
246 146 : logical,allocatable :: bbp_mask(:,:)
247 146 : type(pawcprj_type),allocatable :: Cprj1_bz(:,:),Cprj2_bz(:,:), Cprj1_ibz(:,:),Cprj2_ibz(:,:)
248 146 : type(pawpwij_t),allocatable :: Pwij(:),Pwij_fft(:)
249 146 : type(pawhur_t),allocatable :: Hur(:)
250 146 : type(vkbr_t),allocatable :: vkbr(:)
251 : !************************************************************************
252 :
253 : DBG_ENTER("COLL")
254 :
255 146 : call cwtime(cpu_all, wall_all, gflops_all, "start")
256 :
257 : ! Change FFT mesh if needed
258 230 : if (ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3))) call wfd%change_ngfft(Cryst,Psps,ngfft_gw)
259 :
260 584 : gw_mgfft = MAXVAL(ngfft_gw(1:3)); gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10)
261 146 : if (Dtset%pawcross==1) mgfftf = MAXVAL(ngfftf(1:3))
262 :
263 : ! Copy important variables.
264 146 : comm = Wfd%comm; nsppol = Wfd%nsppol; nspinor = Wfd%nspinor; mband = Wfd%mband; nfft = Wfd%nfft
265 146 : ABI_CHECK(Wfd%nfftot == nfftot_gw, "Wrong nfftot_gw")
266 : mpw = maxval(wfd%npwarr)
267 146 : dim_rtwg = 1 !; if (nspinor==2) dim_rtwg=2 ! Can reduce size depending on Ep%nI and Ep%nj
268 :
269 146 : is_metallic = qp_ebands%has_metal_scheme()
270 146 : ucrpa_bands(1)=dtset%ucrpa_bands(1)
271 146 : ucrpa_bands(2)=dtset%ucrpa_bands(2)
272 146 : luwindow=.false.
273 146 : if (abs(dtset%ucrpa_window(1)+1_dp)>tol8.or.(abs(dtset%ucrpa_window(2)+1_dp)>tol8)) luwindow=.true.
274 :
275 : ! For cRPA calculation of U: read forlb.ovlp
276 146 : if(dtset%ucrpa>=1 .AND. dtset%plowan_compute<10) then
277 : call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor,Kmesh,lcor,luwindow,&
278 0 : nspinor,nsppol,pawang,dtset%prtvol,ucrpa_bands)
279 : endif
280 :
281 146 : ks_eig => ks_ebands%eig
282 146 : qp_eig => qp_ebands%eig; qp_occ => qp_ebands%occ
283 :
284 255384 : chi0_lwing = czero; chi0_uwing = czero; chi0_head = czero
285 :
286 146 : if (Psps%usepaw == 0) then
287 141 : if (Ep%inclvkb /= 0) then
288 : ! Include the term <n,k|[Vnl,iqr]|n"k>' for q -> 0.
289 84 : ABI_CHECK(nspinor == 1, "nspinor with inclvkb not coded")
290 : else
291 57 : ABI_WARNING('Neglecting <n,k|[Vnl,iqr]|m,k>')
292 : end if
293 :
294 : else
295 : ! For PAW+DFT+U, precalculate <\phi_i|[Hu,r]|phi_j\>
296 46 : ABI_MALLOC(HUr, (Cryst%natom))
297 5 : if (Dtset%usepawu /= 0) then
298 0 : call pawhur_init(hur,nsppol,Dtset%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij)
299 : end if
300 : end if
301 :
302 : ! TODO: use ddkop instead of commutator so that we can handle SOC terms.
303 : ! Unfortunately, PAW is not supported yet. Also m_ddk should be relocated below 72_response
304 : !use_ddk = .False.
305 : !!use_ddk = .True.
306 : !use_ddk = dtset%userie == 432
307 : !if (use_ddk) then
308 : ! call wrtout(std_out, " Using DDK to compute the commutator matrix elements.")
309 : ! call ddkop%init(dtset, cryst, pawtab, psps, mpi_enreg, mpw, ngfft_gw)
310 : !end if
311 :
312 : ! Initialize the completeness correction.
313 438 : ABI_MALLOC(green_enhigh_w, (Ep%nomega))
314 1056 : green_enhigh_w = czero
315 :
316 146 : if (Ep%gwcomp == 1) then
317 77 : en_high = MAXVAL(qp_eig(Ep%nbnds,:,:))+Ep%gwencomp
318 19 : write(msg,'(a,f8.2,a)')' Using completeness correction with energy ',en_high*Ha_eV,' [eV] '
319 19 : call wrtout(std_out, msg)
320 57 : ABI_MALLOC(wfwfg,(nfft*nspinor**2))
321 :
322 : ! Init the largest G-sphere contained in the FFT box for the wavefunctions.
323 19 : call Gsph_FFT%in_fftbox(Cryst,Wfd%ngfft)
324 38 : call Gsph_FFT%print([std_out], dtset%prtvol)
325 :
326 57 : ABI_MALLOC(gspfft_igfft,(Gsph_FFT%ng))
327 76 : ABI_MALLOC(dummy_gbound,(2*gw_mgfft+8,2))
328 :
329 : ! Mapping between G-sphere and FFT box.
330 19 : call Gsph_FFT%fft_tabs([0, 0, 0],Wfd%mgfft,Wfd%ngfft,dummy,dummy_gbound,gspfft_igfft)
331 19 : ABI_FREE(dummy_gbound)
332 :
333 19 : if (Psps%usepaw==1) then
334 : ! Prepare the onsite contributions on the GW FFT mesh.
335 0 : ABI_MALLOC(gw_gfft,(3,nfft))
336 0 : q0=zero
337 0 : call get_gfft(ngfft_gw,q0,Cryst%gmet,gw_gsq,gw_gfft) ! The set of plane waves in the FFT Box.
338 0 : ABI_MALLOC(Pwij_fft,(Psps%ntypat))
339 0 : call pawpwij_init(Pwij_fft,nfft,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
340 : end if
341 : end if
342 :
343 : ! Setup weight (2 for spin unpolarized systems, 1 for polarized).
344 : ! spin_fact is used to normalize the occupation factors to one.
345 : ! Consider also the AFM case.
346 144 : select case (nsppol)
347 : case (1)
348 144 : weight = two / Kmesh%nbz; spin_fact = half
349 144 : if (Wfd%nspden == 2) then
350 1 : weight = one / Kmesh%nbz; spin_fact = half
351 : end if
352 144 : if (nspinor == 2) then
353 2 : weight = one / Kmesh%nbz; spin_fact = one
354 : end if
355 : case (2)
356 2 : weight = one / Kmesh%nbz; spin_fact = one
357 : case default
358 146 : ABI_BUG(sjoin("Wrong nsppol:", itoa(nsppol)))
359 : end select
360 :
361 : ! k-weights for points in the IBZ_q
362 7018 : wtk_ltg(:) = 1
363 146 : if (Ep%symchi == 1) then
364 6498 : do ik_bz=1,Ltg_q%nbz
365 6366 : wtk_ltg(ik_bz) = 0
366 6366 : if (Ltg_q%ibzq(ik_bz) /= 1) CYCLE ! Only k in IBZ_q
367 110568 : wtk_ltg(ik_bz) = sum(Ltg_q%wtksym(:,:,ik_bz))
368 : end do
369 : end if
370 :
371 146 : write(msg,'(a,i3,a)')' Q-points for long wave-length limit. # ',Ep%nqlwl,ch10
372 292 : do iqlwl=1,Ep%nqlwl
373 146 : write(msg_tmp,'(1x,i5,a,2x,3f12.6,a)') iqlwl,')',Ep%qlwl(:,iqlwl),ch10
374 292 : msg=TRIM(msg)//msg_tmp
375 : end do
376 146 : call wrtout(std_out, msg)
377 :
378 : write(msg,'(a,i2,2a,i2)')&
379 146 : ' Using spectral method for the imaginary part = ',Ep%spmeth,ch10,&
380 292 : ' Using symmetries to sum only over the IBZ_q = ',Ep%symchi
381 146 : call wrtout(std_out, msg)
382 :
383 146 : if (use_tr) then
384 135 : call wrtout(std_out, ' Using faster algorithm based on time reversal symmetry.')
385 : else
386 11 : call wrtout(std_out, ' Using slow algorithm without time reversal symmetry.')
387 : end if
388 :
389 : ! Evaluate oscillator matrix elements btw partial waves. Note q=Gamma
390 146 : if (Psps%usepaw == 1) then
391 24 : ABI_MALLOC(Pwij,(Psps%ntypat))
392 5 : call pawpwij_init(Pwij,Ep%npwepG0, [zero, zero, zero],Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
393 :
394 56 : ABI_MALLOC(Cprj1_bz,(Cryst%natom,nspinor))
395 5 : call pawcprj_alloc(Cprj1_bz,0,Wfd%nlmn_atm)
396 51 : ABI_MALLOC(Cprj2_bz,(Cryst%natom,nspinor))
397 5 : call pawcprj_alloc(Cprj2_bz,0,Wfd%nlmn_atm)
398 :
399 51 : ABI_MALLOC(Cprj1_ibz,(Cryst%natom,nspinor))
400 5 : call pawcprj_alloc(Cprj1_ibz,0,Wfd%nlmn_atm)
401 51 : ABI_MALLOC(Cprj2_ibz,(Cryst%natom,nspinor))
402 5 : call pawcprj_alloc(Cprj2_ibz,0,Wfd%nlmn_atm)
403 5 : if (Dtset%pawcross==1) then
404 0 : ABI_MALLOC(ur_ae1,(nfftf_tot*nspinor))
405 0 : ABI_MALLOC(ur_ae_onsite1,(nfftf_tot*nspinor))
406 0 : ABI_MALLOC(ur_ps_onsite1,(nfftf_tot*nspinor))
407 0 : ABI_MALLOC(ur_ae2,(nfftf_tot*nspinor))
408 0 : ABI_MALLOC(ur_ae_onsite2,(nfftf_tot*nspinor))
409 0 : ABI_MALLOC(ur_ps_onsite2,(nfftf_tot*nspinor))
410 0 : ABI_MALLOC(igfftepsG0f,(Ep%npwepG0))
411 0 : ABI_MALLOC(tabrf_k,(nfftf_tot))
412 : end if
413 : end if
414 :
415 438 : ABI_MALLOC(rhotwg,(Ep%npwe*dim_rtwg))
416 438 : ABI_MALLOC(tabr_k,(nfft))
417 438 : ABI_MALLOC(ur1_kibz,(nfft*nspinor))
418 292 : ABI_MALLOC(ur2_kibz,(nfft*nspinor))
419 292 : ABI_MALLOC(usr1_k,(nfft*nspinor))
420 292 : ABI_MALLOC(ur2_k,(nfft*nspinor))
421 : !
422 : ! Tables for the FFT of the oscillators.
423 : ! a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere).
424 : ! b) gw_gbound table for the zero-padded FFT performed in rhotwg.
425 :
426 438 : ABI_MALLOC(igffteps0,(Gsph_epsG0%ng))
427 584 : ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2))
428 146 : call Gsph_epsG0%fft_tabs([0, 0, 0], gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igffteps0)
429 146 : if (ANY(gw_fftalga == [2, 4])) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
430 : !use_padfft = 0
431 146 : if (use_padfft==0) then
432 0 : ABI_FREE(gw_gbound)
433 0 : ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft))
434 : end if
435 146 : if (Dtset%pawcross==1) then
436 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2))
437 0 : call Gsph_epsG0%fft_tabs((/0,0,0/),mgfftf,ngfftf,use_padfftf,gboundf,igfftepsG0f)
438 0 : if ( ANY(gw_fftalga == (/2,4/)) ) use_padfftf=0
439 0 : if (use_padfftf==0) then
440 0 : ABI_FREE(gboundf)
441 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf))
442 : end if
443 : end if
444 :
445 : ! TODO this table can be calculated for each k-point
446 390 : my_nbbpks=0; allup="All"; got=0
447 876 : ABI_MALLOC(bbp_ks_distrb,(mband,mband,Kmesh%nbz,nsppol))
448 : call wrtout(std_out, sjoin(' Memory needed for bbp_ks_distrb: ', &
449 146 : ftoa(four*mband**2*Kmesh%nbz*nsppol*b2Mb, fmt="f8.1"), ' [Mb] <<< MEM'))
450 :
451 584 : ABI_MALLOC(bbp_mask, (mband, mband))
452 :
453 294 : do spin=1,nsppol
454 7294 : do ik_bz=1,Kmesh%nbz
455 :
456 7000 : if (Ep%symchi==1) then
457 6430 : if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q
458 : end if
459 1354 : ik_ibz=Kmesh%tab(ik_bz)
460 :
461 : call chi0_bbp_mask(ik_ibz, ik_ibz, spin, spin_fact, use_tr, &
462 1354 : ep%gwcomp, ep%spmeth, ep%nbnds, mband, qp_ebands, bbp_mask)
463 :
464 1354 : call wfd%distribute_bbp(ik_ibz,spin,allup,my_nbbp,bbp_ks_distrb(:,:,ik_bz,spin),got=got,bbp_mask=bbp_mask)
465 7148 : my_nbbpks = my_nbbpks + my_nbbp
466 : end do
467 : end do
468 :
469 146 : ABI_FREE(bbp_mask)
470 :
471 146 : write(msg,'(a,i0,a)')" Will sum ",my_nbbpks," (b,b',k,s) states in chi0q0."
472 146 : call wrtout(std_out, msg)
473 :
474 286 : SELECT CASE (Ep%spmeth)
475 : CASE (0)
476 140 : call wrtout(std_out,' Calculating chi0(q=(0,0,0),omega,G,G")')
477 420 : ABI_MALLOC(green_w, (Ep%nomega))
478 :
479 : CASE (1, 2)
480 6 : call wrtout(std_out,' Calculating Im chi0(q=(0,0,0),omega,G,G")')
481 : !
482 : ! === Find max and min resonant transitions for this q, report values for this processor ===
483 : call make_transitions(Wfd,1,Ep%nbnds,nbvw,nsppol,Ep%symchi,Cryst%timrev,GW_TOL_DOCC,&
484 6 : max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,qp_eig,qp_occ, [zero,zero,zero], bbp_ks_distrb)
485 :
486 : ! === Calculate frequency dependent weights for Kramers Kronig transform ===
487 18 : ABI_MALLOC(omegasf, (Ep%nomegasf))
488 24 : ABI_MALLOC(kkweight, (Ep%nomegasf,Ep%nomega))
489 : !my_wl=1; my_wr=Ep%nomegasf
490 : call setup_spectral(Ep%nomega,Ep%omega,Ep%nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,&
491 6 : 0,Ep%zcut,zero,my_wl,my_wr,kkweight)
492 :
493 6 : ABI_CHECK(use_tr, 'Hilbert transform requires time-reversal')
494 :
495 : ! allocate head and wings of the spectral function.
496 18 : ABI_MALLOC(sf_head,(3,3,my_wl:my_wr))
497 30 : ABI_MALLOC(sf_lwing,(Ep%npwe,my_wl:my_wr,3))
498 18 : ABI_MALLOC(sf_uwing,(Ep%npwe,my_wl:my_wr,3))
499 261742 : sf_head=czero; sf_lwing=czero; sf_uwing=czero
500 :
501 6 : memreq = two*gwp*Ep%npwe**2*(my_wr-my_wl+1)*b2Gb
502 6 : write(msg,'(a,f10.4,a)')' memory required per spectral point: ',two*gwp*Ep%npwe**2*b2Mb,' [Mb]'
503 6 : call wrtout(std_out, msg)
504 6 : write(msg,'(a,f10.4,a)')' memory required by sf_chi0q0: ',memreq,' [Gb]'
505 6 : call wrtout(std_out, msg)
506 6 : if (memreq > two) then
507 0 : ABI_WARNING(' Memory required for sf_chi0q0 is larger than 2.0 Gb!')
508 : end if
509 30 : ABI_MALLOC_OR_DIE(sf_chi0,(Ep%npwe,Ep%npwe,my_wl:my_wr), ierr)
510 1525306 : sf_chi0=czero_gw
511 :
512 : CASE DEFAULT
513 146 : ABI_BUG("Wrong spmeth")
514 : END SELECT
515 :
516 146 : nkpt_summed = Kmesh%nbz
517 146 : if (Ep%symchi /= 0) then
518 132 : nkpt_summed = Ltg_q%nibz_ltg
519 264 : call Ltg_q%print([std_out], prtvol=Dtset%prtvol)
520 : end if
521 146 : call wrtout(std_out, sjoin(' Calculation status: ', itoa(nkpt_summed), ' k-points to be completed'))
522 146 : call pstat_proc%print(_PSTAT_ARGS_)
523 :
524 1287 : ABI_MALLOC(vkbr, (Kmesh%nibz))
525 995 : gradk_not_done = .TRUE.
526 :
527 : ! ============================================
528 : ! === Begin big fat loop over transitions ====
529 : ! ============================================
530 3348948 : chi0 = czero_gw; chi0_sumrule = zero
531 :
532 : ! Loop on spin to calculate $\chi_{\up,\up} + \chi_{\down,\down}$
533 294 : do spin=1,nsppol
534 1113 : if (ALL(bbp_ks_distrb(:,:,:,spin) /= Wfd%my_rank)) CYCLE
535 :
536 : ! Loop over k-points in the BZ.
537 7294 : do ik_bz=1,Kmesh%nbz
538 7000 : if (Ep%symchi == 1) then
539 6430 : if (Ltg_q%ibzq(ik_bz) /= 1) CYCLE ! Only IBZ_q
540 : end if
541 :
542 8401 : if (ALL(bbp_ks_distrb(:,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE
543 :
544 1354 : print_time = wfd%my_rank == 0 .and. (ik_bz < LOG_MODK .or. mod(ik_bz, LOG_MODK) == 0)
545 392 : if (print_time) call cwtime(cpu_k, wall_k, gflops_k, "start")
546 :
547 : ! Get ik_ibz, non-symmorphic phase and symmetries from ik_bz.
548 1354 : call kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt)
549 6904175 : tabr_k=ktabr(:,ik_bz) ! Table for rotated FFT points
550 6770 : spinrot_kbz(:)=Cryst%spinrot(:,isym_k)
551 1354 : if (Dtset%pawcross==1) tabrf_k(:) = ktabrf(:,ik_bz)
552 :
553 1354 : istwf_k = Wfd%istwfk(ik_ibz)
554 1354 : npw_k = Wfd%npwarr(ik_ibz)
555 1354 : kg_k => Wfd%Kdata(ik_ibz)%kg_k
556 :
557 1354 : if (psps%usepaw == 0 .and. Ep%inclvkb /= 0 .and. gradk_not_done(ik_ibz)) then
558 : ! Include term <n,k|[Vnl,iqr]|n"k>' for q -> 0.
559 494 : call vkbr_init(vkbr(ik_ibz), Cryst, Psps, Ep%inclvkb, istwf_k, npw_k, Kmesh%ibz(:,ik_ibz), kg_k)
560 494 : gradk_not_done(ik_ibz) = .FALSE.
561 : end if
562 :
563 : !if (use_ddk) then
564 : ! call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, Kmesh%ibz(:,ik_ibz), istwf_k, npw_k, kg_k)
565 : ! ABI_CHECK(istwf_k == 1, "istwfk_k1 not coded")
566 : !end if
567 :
568 : ! Loop over "conduction" states.
569 31904 : do band1=1,Ep%nbnds
570 411866 : if (ALL(bbp_ks_distrb(band1,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE
571 :
572 22810 : ABI_CHECK(wfd%get_wave_ptr(band1, ik_ibz, spin, wave1, msg) == 0, msg)
573 22810 : ug1 => wave1%ug
574 22810 : call wfd%get_ur(band1,ik_ibz,spin,ur1_kibz)
575 :
576 : !if (use_ddk) then
577 : ! ! Compute DH_DK |psi_k,bi>, store results in ddk_ug1
578 : ! ddk_ug1(1,:,1) = real(ug1)
579 : ! ddk_ug1(2,:,1) = aimag(ug1)
580 : ! call ddkop%apply(eig_nk, npw_k, nspinor, ddk_ug1(:,:,1), cwaveprj)
581 : ! ddk_ug1(:,:,:) = ddkop%gh1c
582 : !end if
583 :
584 22810 : if (Psps%usepaw==1) then
585 356 : call wfd%get_cprj(band1,ik_ibz,spin,Cryst,Cprj1_ibz,sorted=.FALSE.)
586 356 : call pawcprj_copy(Cprj1_ibz,Cprj1_bz)
587 356 : call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_bz)
588 356 : if (Dtset%pawcross==1) then
589 0 : call wfdf%paw_get_aeur(band1,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae1,ur_ae_onsite1,ur_ps_onsite1)
590 : end if
591 : end if
592 :
593 : ! Loop over "valence" states.
594 1076294 : do band2=1,Ep%nbnds
595 :
596 : if (luwindow.AND.dtset%ucrpa==1 &
597 : .AND.((ks_ebands%eig(band1,ik_ibz,spin)-ks_ebands%fermie)<=dtset%ucrpa_window(2)) &
598 : .AND.((ks_ebands%eig(band1,ik_ibz,spin)-ks_ebands%fermie)>=dtset%ucrpa_window(1)) &
599 : .AND.((ks_ebands%eig(band2,ik_ibz,spin)-ks_ebands%fermie)<=dtset%ucrpa_window(2)) &
600 1052130 : .AND.((ks_ebands%eig(band2,ik_ibz,spin)-ks_ebands%fermie)>=dtset%ucrpa_window(1))) CYCLE
601 : !----------------- cRPA for U
602 :
603 1052130 : if (bbp_ks_distrb(band1,band2,ik_bz,spin) /= Wfd%my_rank) CYCLE
604 :
605 109578 : deltaeKS_b1b2 = ks_eig(band1, ik_ibz, spin) - ks_eig(band2, ik_ibz, spin)
606 109578 : deltaf_b1b2 = spin_fact * (qp_occ(band1, ik_ibz, spin) - qp_occ(band2, ik_ibz, spin))
607 109578 : deltaeGW_b1b2 = qp_eig(band1, ik_ibz, spin) - qp_eig(band2, ik_ibz, spin)
608 :
609 109578 : if (Ep%gwcomp == 0) then
610 : ! Skip negligible transitions.
611 106126 : if (abs(deltaf_b1b2) < GW_TOL_DOCC) CYCLE
612 : else
613 : ! when the completeness trick is used, we need to also consider transitions with vanishing deltaf
614 : ! Rangel Correction for metals
615 3452 : if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC .and. ( ABS(deltaf_b1b2)< GW_TOL_DOCC .or. band1 < band2)) CYCLE
616 : end if
617 :
618 109578 : ABI_CHECK(wfd%get_wave_ptr(band2, ik_ibz, spin, wave2, msg) == 0, msg)
619 109578 : ug2 => wave2%ug
620 109578 : call wfd%get_ur(band2,ik_ibz,spin,ur2_kibz)
621 :
622 109578 : if (Psps%usepaw==1) then
623 1799 : call wfd%get_cprj(band2,ik_ibz,spin,Cryst,Cprj2_ibz,sorted=.FALSE.)
624 1799 : call pawcprj_copy(Cprj2_ibz,Cprj2_bz)
625 1799 : call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj2_bz)
626 1799 : if (Dtset%pawcross==1) then
627 : call wfdf%paw_get_aeur(band2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
628 0 : ur_ae2,ur_ae_onsite2,ur_ps_onsite2)
629 : end if
630 : end if
631 :
632 103914 : SELECT CASE (Ep%spmeth)
633 : CASE (0)
634 : ! Adler-Wiser expression.
635 : ! Add small imaginary of the Time-Ordered response function but only for non-zero real omega
636 : ! FIXME What about metals?
637 :
638 103914 : if (.not. use_tr) then
639 : ! Adler-Wiser without time-reversal.
640 169492 : do io=1,Ep%nomega
641 169492 : green_w(io) = g0g0w(Ep%omega(io), deltaf_b1b2, deltaeGW_b1b2, Ep%zcut, GW_TOL_W0, one_pole)
642 : end do
643 :
644 : else
645 79352 : if (Ep%gwcomp == 0) then ! cannot be completely skipped in case of completeness correction
646 75900 : if (band1 < band2) CYCLE ! Here we GAIN a factor ~2
647 : end if
648 :
649 357819 : do io=1,Ep%nomega
650 : ! Rangel: In metals, the intra-band transitions term does not contain the antiresonant part
651 : ! if(abs(deltaeGW_b1b2)>GW_TOL_W0) green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,Ep%zcut,GW_TOL_W0)
652 278467 : if (band1 == band2) green_w(io) = g0g0w(Ep%omega(io), deltaf_b1b2, deltaeGW_b1b2, Ep%zcut, GW_TOL_W0, one_pole)
653 278467 : if (band1 /= band2) green_w(io) = g0g0w(Ep%omega(io), deltaf_b1b2, deltaeGW_b1b2, Ep%zcut, GW_TOL_W0, two_poles)
654 :
655 357819 : if (Ep%gwcomp == 1) then
656 : ! Calculate the completeness correction
657 16504 : numerator= -spin_fact * qp_occ(band2,ik_ibz,spin)
658 16504 : deltaeGW_enhigh_b2 = en_high - qp_eig(band2,ik_ibz,spin)
659 : ! Completeness correction is NOT valid for real frequencies
660 16504 : if (REAL(Ep%omega(io)) < GW_TOL_W0) then
661 16504 : green_enhigh_w(io) = g0g0w(Ep%omega(io),numerator,deltaeGW_enhigh_b2,Ep%zcut,GW_TOL_W0,two_poles)
662 : else
663 0 : green_enhigh_w(io) = czero_gw
664 : endif
665 : !
666 : ! Rangel Correction for metals
667 16504 : if (band1 >= band2 .and. abs(deltaf_b1b2) > GW_TOL_DOCC) then
668 14004 : green_w(io)= green_w(io) - green_enhigh_w(io)
669 : else
670 : ! Disregard green_w, since it is already accounted for through the time-reversal
671 2500 : green_w(io)= - green_enhigh_w(io)
672 : end if
673 : end if !gwcomp == 1
674 : end do !io
675 :
676 79352 : if (Ep%gwcomp == 1 .and. band1 == band2) then
677 : ! Add the "delta part", symmetrization is done inside the routine.
678 152 : call calc_wfwfg(tabr_k,itim_k,spinrot_kbz,nfft,nspinor,ngfft_gw,ur2_kibz,ur2_kibz,wfwfg)
679 :
680 152 : if (Psps%usepaw==1) then
681 0 : call paw_rho_tw_g(cryst, Pwij_fft, nfft,dim_rtwg,nspinor,gw_gfft,Cprj2_bz,Cprj2_bz,wfwfg)
682 :
683 : ! Add PAW cross term
684 0 : if (Dtset%pawcross==1) then
685 : call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,&
686 : ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,&
687 : ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,&
688 0 : dim_rtwg,wfwfg)
689 : end if
690 : end if
691 :
692 152 : qzero = .TRUE.
693 : call completechi0_deltapart(ik_bz,qzero,Ep%symchi,Ep%npwe,Gsph_FFT%ng,Ep%nomega,nspinor,&
694 152 : nfft,ngfft_gw,gspfft_igfft,Gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0)
695 : end if
696 : end if ! use_tr
697 :
698 : CASE (1, 2)
699 : ! Spectral method, here time-reversal is always assumed.
700 5664 : if (deltaeGW_b1b2 < 0) CYCLE
701 115242 : call approxdelta(Ep%nomegasf,omegasf,deltaeGW_b1b2,Ep%spsmear,iomegal,iomegar,wl,wr,Ep%spmeth)
702 : END SELECT
703 :
704 : ! FFT of u^*_{b1,k}(r) u_{b2,k}(r) and (q,G=0) limit using small q and k.p perturbation theory
705 : call rho_tw_g(nspinor,Ep%npwe,nfft,ndat1,ngfft_gw,1,use_padfft,igffteps0,gw_gbound,&
706 : ur1_kibz,itim_k,tabr_k,ph_mkt,spinrot_kbz,&
707 : ur2_kibz,itim_k,tabr_k,ph_mkt,spinrot_kbz,&
708 109578 : dim_rtwg,rhotwg)
709 :
710 109578 : if (psps%usepaw == 0) then
711 : ! Matrix elements of i[H,r] for NC pseudopotentials.
712 : rhotwx = nc_ihr_comm(vkbr(ik_ibz), cryst, psps, npw_k, nspinor, istwf_k, Ep%inclvkb, &
713 107779 : Kmesh%ibz(:,ik_ibz), ug1, ug2, kg_k)
714 :
715 : else
716 : ! 1) Add PAW onsite contribution, projectors are already in the BZ.
717 1799 : call paw_rho_tw_g(cryst, Pwij, Ep%npwe,dim_rtwg,nspinor,Gsph_epsG0%gvec,Cprj1_bz,Cprj2_bz,rhotwg)
718 :
719 : ! 2) Matrix elements of i[H,r] for PAW.
720 1799 : rhotwx = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug1,ug2,kg_k,Cprj1_ibz,Cprj2_ibz,HUr)
721 :
722 : ! Add PAW cross term
723 1799 : if (Dtset%pawcross==1) then
724 : call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,&
725 : ur_ae1,ur_ae_onsite1,ur_ps_onsite1,itim_k,tabrf_k,ph_mkt,spinrot_kbz,&
726 : ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k,tabrf_k,ph_mkt,spinrot_kbz,&
727 0 : dim_rtwg,rhotwg)
728 :
729 : ! Add cross-term contribution to the commutator
730 0 : if (Dtset%userib/=111) then
731 : call paw_cross_ihr_comm(rhotwx,nspinor,nfftf_tot,Cryst,Pawfgrtab,Paw_onsite,&
732 0 : ur_ae1,ur_ae2,ur_ae_onsite1,ur_ae_onsite2,Cprj1_ibz,Cprj2_ibz)
733 : end if
734 : end if
735 : end if
736 :
737 : ! Treat a possible degeneracy between v and c.
738 109578 : if (abs(deltaeKS_b1b2) > GW_TOL_W0) then
739 567396 : rhotwx = -rhotwx / deltaeKS_b1b2
740 : else
741 1230 : rhotwx = czero_gw
742 : end if
743 :
744 : !if (use_ddk) then
745 : ! cg2_dp(1,:) = real(ug2)
746 : ! cg2_dp(2,:) = aimag(ug2)
747 : ! do iab=1,gwr%nspinor**2
748 : ! do idir=1,3
749 : ! ! DH_DK operator is Hermitian.
750 : ! spad1 = spinor_pad(1,iab); spad2 = spinor_pad(2,iab)
751 : ! rtmp = cg_zdotc(npw_ki, ddk_ug1(:,spad1+1,idir), cg2_dp(:,spad2+1))
752 : ! new_rhotwx(idir, iab) = rtmp(1) + j_dpc * rtmp(2)
753 : ! end do ! iab
754 : ! end do ! idir
755 :
756 : ! if (abs(deltaeKS_b1b2) > GW_TOLQ0) then
757 : ! new_rhotwx = -new_rhotwx / deltaeKS_b1b2
758 : ! do iab=1,gwr%nspinor**2
759 : ! new_rhotwx(:, iab) = matmul(cryst%rmet, new_rhotwx(:, iab)) / (two_pi ** 2)
760 : ! end do
761 : ! else
762 : ! new_rhotwx = zero
763 : ! end if
764 : !end if
765 :
766 103914 : SELECT CASE (Ep%spmeth)
767 : CASE (0)
768 : ! ---------------- Ucrpa (begin)
769 103914 : if(dtset%ucrpa>=1.and..not.luwindow) then
770 0 : fac=one
771 0 : fac1=zero
772 0 : fac2=zero
773 0 : fac3=zero
774 0 : fac4=one
775 0 : m1=-1
776 0 : m2=-1
777 0 : if(dtset%ucrpa<=2) then
778 0 : call flush_unit(std_out)
779 : if ( band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)&
780 0 : & .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then
781 : ! if(dtset%prtvol>=10)write(6,*)"calculation is in progress",band1,band2,ucrpa_bands(1),ucrpa_bands(2)
782 0 : if (dtset%plowan_compute>=10) then
783 0 : band1c=band1-wan%bandi_wan+1
784 0 : band2c=band2-wan%bandi_wan+1
785 0 : do iat1=1,wan%natom_wan
786 0 : do ispinor1=1,wan%nspinor
787 0 : do il1=1,wan%nbl_atom_wan(iat1)
788 0 : do m1=1,2*(wan%latom_wan(iat1)%lcalc(il1))+1
789 : fac1=fac1 + real(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*&
790 0 : &conjg(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))
791 : fac2=fac2 + real(wan%psichi(ik_bz,band2c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*&
792 0 : &conjg(wan%psichi(ik_bz,band2c,iat1)%atom(il1)%matl(m1,spin,ispinor1))
793 0 : do iat2=1,wan%natom_wan
794 0 : do ispinor2=1,wan%nspinor
795 0 : do il2=1,wan%nbl_atom_wan(iat2)
796 0 : do m2=1,2*(wan%latom_wan(iat2)%lcalc(il2))+1
797 : fac=fac - real(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1)*&
798 : & conjg(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*&
799 : & wan%psichi(ik_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)*&
800 0 : & conjg(wan%psichi(ik_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)))
801 : fac3=fac3+ real(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*&
802 : &conjg(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*&
803 : &wan%psichi(ik_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)*&
804 0 : &conjg(wan%psichi(ik_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2))
805 : enddo !m2
806 : enddo !il2
807 : enddo !ispinor2
808 : enddo !iat2
809 : enddo !m1
810 : enddo !il1
811 : enddo !ispinor1
812 : enddo !iat
813 : else !plowan_compute
814 0 : do iat=1, cryst%nattyp(itypatcor)
815 0 : do ispinor1=1,nspinor
816 0 : do m1=1,2*lcor+1
817 : fac1=fac1+ real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
818 0 : &conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)))
819 : fac2=fac2+ real(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)*&
820 0 : & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)))
821 0 : do ispinor2=1,nspinor
822 0 : do m2=1,2*lcor+1
823 : fac=fac - real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
824 : & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))*&
825 : & coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)*&
826 0 : & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)))
827 : fac3=fac3 + real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
828 : & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* &
829 : & coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)*&
830 0 : & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor2,m2)))
831 : ! if(dtset%prtvol>=10)write(6,*) fac,fac3
832 : enddo !m2
833 : enddo !ispinor2
834 : ! if(dtset%prtvol>=10)write(6,*) fac,fac3,fac1,fac2,fac1*fac2
835 : enddo !m1
836 : enddo !ispinor1
837 : enddo !iat
838 : endif !plowan_compute>=10
839 0 : fac4=fac
840 : ! fac=zero
841 0 : if(dtset%ucrpa==1) fac=zero
842 : ! write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0"
843 : endif
844 0 : else if (dtset%ucrpa==3) then
845 0 : if (band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)) then
846 0 : do iat=1, cryst%nattyp(itypatcor)
847 0 : do ispinor1=1,nspinor
848 0 : do m1=1,2*lcor+1
849 : fac2=fac2-real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
850 0 : & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)))
851 : enddo
852 : enddo
853 : enddo
854 : if(dtset%ucrpa==4) fac2=zero
855 : endif
856 0 : if (band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then
857 0 : do iat=1, cryst%nattyp(itypatcor)
858 0 : do ispinor1=1,nspinor
859 0 : do m1=1,2*lcor+1
860 : fac3=fac3-real(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)*&
861 0 : & conjg(coeffW_BZ(iat,spin,band2,ik_bz,ispinor1,m1)))
862 : enddo
863 : enddo
864 : enddo
865 : if(dtset%ucrpa==4) fac3=zero
866 : endif
867 0 : fac=real(fac2*fac3)
868 : endif
869 : !if(dtset%prtvol>=10) write(6,'(6i4,e15.5,a)') ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0"
870 : !if(abs(fac-one)>0.00001) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac," q==0"
871 : ! if(dtset%prtvol>=10) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ik_bz,band1,band2,m1,m2,fac4," q==0"
872 0 : green_w=green_w*fac
873 : endif
874 : ! ---------------- Ucrpa (end)
875 :
876 : ! Adler-Wiser expression, to be consistent here we use the KS eigenvalues (?)
877 : call accumulate_chi0_q0(is_metallic,ik_bz,isym_k,itim_k,Ep%gwcomp,nspinor,Ep%npwepG0,Ep,&
878 103914 : Cryst,Ltg_q,Gsph_epsG0,chi0,rhotwx,rhotwg,green_w,green_enhigh_w,deltaf_b1b2,chi0_head,chi0_lwing,chi0_uwing)
879 :
880 : CASE (1, 2)
881 : ! Spectral method, to be consistent here we use the KS eigenvalues.
882 : call accumulate_sfchi0_q0(ik_bz,isym_k,itim_k,nspinor,Ep%symchi,Ep%npwepG0,Ep%npwe,Cryst,Ltg_q,&
883 : Gsph_epsG0,deltaf_b1b2,my_wl,iomegal,wl,my_wr,iomegar,wr,rhotwx,rhotwg,Ep%nomegasf,&
884 5664 : sf_chi0,sf_head,sf_lwing,sf_uwing)
885 :
886 : CASE DEFAULT
887 109578 : ABI_BUG("Wrong spmeth")
888 : END SELECT
889 :
890 : ! Accumulating the sum rule on chi0. Eq. (5.284) in G.D. Mahan Many-Particle Physics 3rd edition. [[cite:Mahan2000]]
891 109578 : factor = spin_fact * qp_occ(band2,ik_ibz,spin)
892 :
893 : call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_b1b2,&
894 109578 : Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule)
895 :
896 140128 : if (Ep%gwcomp == 1) then
897 : ! Include also the completeness correction in the sum rule.
898 3452 : factor=-spin_fact*qp_occ(band2,ik_ibz,spin)
899 : call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_enhigh_b2,&
900 3452 : Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule)
901 3452 : if (band1 == Ep%nbnds) then
902 13404 : chi0_sumrule(:) = chi0_sumrule(:) + wtk_ltg(ik_bz)*spin_fact*qp_occ(band2,ik_ibz,spin)*deltaeGW_enhigh_b2
903 : end if
904 : end if
905 :
906 : end do ! band2
907 : end do ! band1
908 :
909 1354 : if (Psps%usepaw == 0 .and. Ep%inclvkb /= 0 .and. Ep%symchi == 1) then
910 450 : call vkbr_free(vkbr(ik_ibz)) ! Not need anymore as we loop only over IBZ.
911 : end if
912 :
913 1502 : if (print_time) then
914 392 : write(msg,'(3(a,i0))')' chi0q0: ik_bz: ',ik_bz,'/',Kmesh%nbz,", spin: ",spin
915 392 : call cwtime_report(msg, cpu_k, wall_k, gflops_k); if (ik_bz == LOG_MODK) call wrtout(std_out, " ...")
916 : end if
917 :
918 : end do !ik_bz
919 : end do !spin
920 :
921 146 : ABI_FREE(igffteps0)
922 :
923 146 : call vkbr_free(vkbr)
924 995 : ABI_FREE(vkbr)
925 : !call ddkop%free()
926 :
927 : ! === After big fat loop over transitions, now MPI ===
928 : ! * Master took care of the contribution in case of (metallic|spin) polarized systems.
929 140 : select case (Ep%spmeth)
930 : case (0)
931 : ! Adler-Wiser expression. Sum contributions from each proc.
932 : ! Looping on frequencies to avoid problems with the size of the MPI packet.
933 946 : do io=1,Ep%nomega
934 946 : call xmpi_sum(chi0(:,:,io),comm,ierr)
935 : end do
936 :
937 : case (1, 2)
938 : ! Spectral method.
939 6 : call hilbert_transform(Ep%npwe,Ep%nomega,Ep%nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,Ep%spmeth)
940 :
941 6 : ABI_SFREE(sf_chi0)
942 :
943 : ! Sum contributions from each proc
944 : ! Looping on frequencies to avoid problems with the size of the MPI packet
945 110 : do io=1,Ep%nomega
946 110 : call xmpi_sum(chi0(:,:,io),comm,ierr)
947 : end do
948 :
949 : call hilbert_transform_headwings(Ep%npwe,Ep%nomega,Ep%nomegasf,&
950 : my_wl,my_wr,kkweight,sf_lwing,sf_uwing,sf_head,chi0_lwing,&
951 6 : chi0_uwing,chi0_head,Ep%spmeth)
952 :
953 : case default
954 146 : ABI_BUG(sjoin("Wrong spmeth:", itoa(ep%spmeth)))
955 : end select
956 :
957 : ! Divide by the volume
958 : !$OMP PARALLEL WORKSHARE
959 3340854 : chi0 = chi0 * weight / Cryst%ucvol
960 : !$OMP END PARALLEL WORKSHARE
961 :
962 : ! Collect sum rule. pi comes from Im[1/(x-ieta)] = pi delta(x)
963 146 : call xmpi_sum(chi0_sumrule, comm, ierr)
964 8240 : chi0_sumrule = chi0_sumrule * pi * weight / Cryst%ucvol
965 :
966 : ! Collect head and wings.
967 146 : call xmpi_sum(chi0_head, comm, ierr)
968 146 : call xmpi_sum(chi0_lwing, comm, ierr)
969 146 : call xmpi_sum(chi0_uwing, comm, ierr)
970 :
971 11976 : chi0_head = chi0_head * weight / cryst%ucvol
972 : ! Tensor in terms of reciprocal lattice vectors.
973 1056 : do io=1,Ep%nomega
974 58386 : chi0_head(:,:,io) = matmul(chi0_head(:,:,io), cryst%gmet) * (two_pi**2)
975 : end do
976 121850 : chi0_lwing = chi0_lwing * weight / cryst%ucvol
977 121850 : chi0_uwing = chi0_uwing * weight / cryst%ucvol
978 :
979 : ! ===============================================
980 : ! ==== Symmetrize chi0 in case of AFM system ====
981 : ! ===============================================
982 : ! Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$.
983 : ! Works only in the case of magnetic group Shubnikov type IV.
984 146 : if (cryst%use_antiferro) then
985 : call symmetrize_afm_chi0(Cryst, Gsph_epsG0, Ltg_q, Ep%npwe, Ep%nomega, chi0=chi0, &
986 1 : chi0_head=chi0_head, chi0_lwing=chi0_lwing, chi0_uwing=chi0_uwing)
987 : end if
988 :
989 : ! ===================================================
990 : ! ==== Construct head and wings from the tensor =====
991 : ! ===================================================
992 1056 : do io=1,Ep%nomega
993 39512 : do ig=2,Ep%npwe
994 154408 : wng = chi0_uwing(ig,io,:)
995 38602 : chi0(1,ig,io) = vdotw(Ep%qlwl(:,1), wng, Cryst%gmet,"G")
996 154408 : wng = chi0_lwing(ig,io,:)
997 39512 : chi0(ig,1,io) = vdotw(Ep%qlwl(:,1), wng, Cryst%gmet,"G")
998 : end do
999 14560 : chq = matmul(chi0_head(:,:,io), Ep%qlwl(:,1))
1000 1056 : chi0(1,1,io) = vdotw(Ep%qlwl(:,1), chq, Cryst%gmet,"G") ! Use user-defined small q
1001 : end do
1002 :
1003 146 : if (wfd%my_rank == 0 .and. dtset%prtvol > 20) then
1004 0 : qlen = tol3
1005 0 : call cryst%get_redcart_qdirs(nq, qdirs, qlen=qlen)
1006 0 : ABI_MALLOC(head_qvals, (nq))
1007 0 : call wrtout([std_out, ab_out], "Head of the irreducible polarizability for q --> 0", pre_newlines=1)
1008 0 : call wrtout([std_out, ab_out], sjoin(" q0_len:", ftoa(qlen), "(Bohr^-1)"))
1009 0 : write(msg, "(*(a14))") "omega_re (eV)", "omega_im (eV)", "[100]", "[010]", "[001]", "x", "y", "z"
1010 0 : call wrtout([std_out, ab_out], msg)
1011 0 : do io=1,Ep%nomega
1012 0 : do iq=1,nq
1013 0 : chq = matmul(chi0_head(:,:,io), qdirs(:,iq))
1014 0 : head_qvals(iq) = vdotw(qdirs(:, iq), chq, cryst%gmet, "G")
1015 : end do
1016 0 : write(msg, "(*(es12.5,2x))") ep%omega(io) * Ha_eV, real(head_qvals(:))
1017 0 : call wrtout([std_out, ab_out], msg)
1018 : ! Write imag part to std_out
1019 0 : write(msg, "(*(es12.5,2x))") ep%omega(io) * Ha_eV, aimag(head_qvals(:))
1020 0 : call wrtout(std_out, msg)
1021 : end do
1022 0 : call wrtout([std_out, ab_out], " ")
1023 0 : ABI_FREE(qdirs)
1024 0 : ABI_FREE(head_qvals)
1025 : end if
1026 : !stop
1027 :
1028 : ! Impose Hermiticity (valid only for zero or purely imaginary frequencies)
1029 : ! MG: what about metals, where we have poles around zero?
1030 : ! FB: because of the intraband term, chi0 is never hermitian in case of metals
1031 146 : if (.not. is_metallic) then
1032 1028 : do io=1,Ep%nomega
1033 1028 : if (ABS(REAL(Ep%omega(io))) < 0.00001) then
1034 32538 : do ig2=1,Ep%npwe
1035 1538823 : do ig1=1,ig2-1
1036 1538205 : chi0(ig2,ig1,io)=GWPC_CONJG(chi0(ig1,ig2,io))
1037 : end do
1038 : end do
1039 : end if
1040 : end do
1041 : end if
1042 :
1043 : ! =====================
1044 : ! ==== Free memory ====
1045 : ! =====================
1046 146 : ABI_FREE(bbp_ks_distrb)
1047 146 : ABI_FREE(rhotwg)
1048 146 : ABI_FREE(tabr_k)
1049 146 : ABI_FREE(ur1_kibz)
1050 146 : ABI_FREE(ur2_kibz)
1051 146 : ABI_FREE(usr1_k)
1052 146 : ABI_FREE(ur2_k)
1053 146 : ABI_FREE(gw_gbound)
1054 :
1055 146 : if (Dtset%pawcross==1) then
1056 0 : ABI_FREE(gboundf)
1057 : end if
1058 :
1059 146 : ABI_SFREE(green_enhigh_w)
1060 146 : ABI_SFREE(gw_gfft)
1061 146 : ABI_SFREE(wfwfg)
1062 146 : ABI_SFREE(kkweight)
1063 146 : ABI_SFREE(omegasf)
1064 146 : ABI_SFREE(green_w)
1065 146 : ABI_SFREE(sf_head)
1066 146 : ABI_SFREE(sf_lwing)
1067 146 : ABI_SFREE(sf_uwing)
1068 146 : ABI_SFREE(gspfft_igfft)
1069 :
1070 146 : call Gsph_FFT%free()
1071 :
1072 146 : if (Psps%usepaw==1) then
1073 : ! deallocation for PAW.
1074 5 : call pawcprj_free(Cprj1_bz)
1075 36 : ABI_FREE(Cprj1_bz)
1076 5 : call pawcprj_free(Cprj2_bz)
1077 36 : ABI_FREE(Cprj2_bz)
1078 5 : call pawcprj_free(Cprj1_ibz)
1079 36 : ABI_FREE(Cprj1_ibz)
1080 5 : call pawcprj_free(Cprj2_ibz)
1081 36 : ABI_FREE(Cprj2_ibz)
1082 5 : call pawpwij_free(Pwij)
1083 14 : ABI_FREE(Pwij)
1084 5 : if (allocated(Pwij_fft)) then
1085 0 : call pawpwij_free(Pwij_fft)
1086 0 : ABI_FREE(Pwij_fft)
1087 : end if
1088 5 : call pawhur_free(Hur)
1089 36 : ABI_FREE(Hur)
1090 5 : if (Dtset%pawcross==1) then
1091 0 : ABI_FREE(ur_ae1)
1092 0 : ABI_FREE(ur_ae_onsite1)
1093 0 : ABI_FREE(ur_ps_onsite1)
1094 0 : ABI_FREE(ur_ae2)
1095 0 : ABI_FREE(ur_ae_onsite2)
1096 0 : ABI_FREE(ur_ps_onsite2)
1097 0 : ABI_FREE(tabrf_k)
1098 0 : ABI_FREE(gboundf)
1099 0 : ABI_FREE(igfftepsG0f)
1100 : end if
1101 : end if
1102 :
1103 146 : if(dtset%ucrpa >= 1 .AND. dtset%plowan_compute < 10) then
1104 0 : ABI_FREE(coeffW_BZ)
1105 : endif
1106 :
1107 146 : call cwtime_report("cchi0q0", cpu_all, wall_all, gflops_all)
1108 :
1109 : DBG_EXIT("COLL")
1110 :
1111 292 : end subroutine cchi0q0
1112 : !!***
1113 :
1114 :
1115 : !!****f* ABINIT/cchi0
1116 : !! NAME
1117 : !! cchi0
1118 : !!
1119 : !! FUNCTION
1120 : !! Main calculation of the independent-particle susceptibility chi0 for qpoint != 0
1121 : !!
1122 : !! INPUTS
1123 : !! use_tr=If .TRUE. valence states are stored in Wfs_val and only resonant transitions are calculated
1124 : !! (time reversal is assumed)
1125 : !! Dtset <type(dataset_type)>=all input variables in this dataset
1126 : !! Cryst<crystal_t>= data type gathering info on symmetries and unit cell
1127 : !! %natom=number of atoms
1128 : !! %nsym=number of symmetries
1129 : !! %xred(3,natom)=reduced coordinated of atoms
1130 : !! %typat(natom)=type of each atom
1131 : !! %rprimd(3,3)=dimensional primitive translations in real space (bohr)
1132 : !! %timrev= 2 if time reversal can be used, 1 otherwise
1133 : !! qpoint(3)=reciprocal space coordinates of the q wavevector
1134 : !! Ep<type(em1params_t_type)>= Parameters related to the calculation of the inverse dielectric matrix.
1135 : !! %nbnds=number of bands summed over
1136 : !! %npwe=number of planewaves for the irreducible polarizability X^0_GGp
1137 : !! %npwvec=maximum number of G vectors
1138 : !! used to define the dimension of some arrays e.g igfft
1139 : !! %nsppol=1 for unpolarized, 2 for spin-polarized
1140 : !! %nomega=total number of frequencies in X^0 (both real and imaginary)
1141 : !! %nomegasf=number of real frequencies used to sample the imaginary part of X^0 (spectral method)
1142 : !! %spmeth=1 if we use the spectral method, 0 for standard Adler-Wiser expression
1143 : !! %spsmear=gaussian broadening used to approximate the delta distribution
1144 : !! %zcut=small imaginary shift to avoid poles in X^0
1145 : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
1146 : !! Kmesh <kmesh_t>= datatype gathering parameters related to the k-point sampling
1147 : !! %nibz=number of k-points in the IBZ
1148 : !! %nbz=number of k-points in the BZ
1149 : !! %bz(3,nbz)=reduced coordinates for k-points in the full Brillouin zone
1150 : !! %ibz(3,nibz)=reduced coordinates for k-points in the irreducible wedge
1151 : !! %tab(nbz)=mapping between a kpt in the BZ (array bz) and the irred point in the array ibz
1152 : !! %tabi(nbz)= -1 if inversion is needed to obtain this particular kpt in the BZ, 1 means identity
1153 : !! %tabo(nbz)= for each point in the BZ, the index of the symmetry operation S in reciprocal
1154 : !! space which rotates k_IBZ onto \pm k_BZ (depending on tabi)
1155 : !! %tabp(nbz)= For each k_BZ, it gives the phase factors associated to non-symmorphic operations, i.e
1156 : !! e^{-i 2 \pi k_IBZ \cdot R{^-1}t} == e{-i 2\pi k_BZ cdot t} where :
1157 : !! \transpose R{-1}=S and (S k_IBZ) = \pm k_BZ (depending on ktabi)
1158 : !! %tabr(nfftot,nbz) For each point r on the real mesh and for each k-point in the BZ, tabr
1159 : !! gives the index of (R^-1 (r-t)) in the FFT array where R=\transpose S^{-1} and k_BZ=S k_IBZ.
1160 : !! t is the fractional translation associated to R
1161 : !! Gsph_epsG0<gsphere_t data type> The G-sphere used to describe chi0/eps. (including umklapp G0 vectors)
1162 : !! %ng=number of G vectors for chi0
1163 : !! %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G
1164 : !! %phmGt(Ep%npwe,nsym)=phase factors e^{-iG \cdot t} needed to symmetrize oscillator matrix elements and epsilon
1165 : !! %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
1166 : !! %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
1167 : !! nbvw=number of bands in the arrays wfrv
1168 : !! ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths (see input variable)
1169 : !! nfftot_gw=Total number of points in the GW FFT grid
1170 : !! Ltg_q<Little group>=Data type gathering information on the little group of the q-points.
1171 : !! Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
1172 : !! Pawang<pawang_type> angular mesh discretization and related data:
1173 : !! qp_ebands<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities)
1174 : !! %mband=MAX number of bands over k-points and spin (==Ep%nbnds)
1175 : !! %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band
1176 : !! %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes
1177 : !! Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
1178 : !! Wfd<wfdgw_t>=Object used to access the wavefunctions
1179 : !!
1180 : !! OUTPUT
1181 : !! chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix at wavevector qpoint and
1182 : !! each frequency defined by Ep%omega and Ep%nomega.
1183 : !!
1184 : !! SOURCE
1185 :
1186 776 : subroutine cchi0(use_tr,Dtset,Cryst,qpoint,Ep,Psps,Kmesh,qp_ebands,Gsph_epsG0,&
1187 776 : Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,nbvw,ngfft_gw,nfftot_gw,ngfftf,nfftf_tot,&
1188 776 : chi0,ktabr,ktabrf,Ltg_q,chi0_sumrule,Wfd,Wfdf,wan)
1189 :
1190 : !Arguments ------------------------------------
1191 : !scalars
1192 : integer,intent(in) :: nbvw,nfftot_gw,nfftf_tot
1193 : logical,intent(in) :: use_tr
1194 : type(ebands_t),target,intent(in) :: qp_ebands
1195 : type(kmesh_t),intent(in) :: Kmesh
1196 : type(crystal_t),intent(in) :: Cryst
1197 : type(Dataset_type),intent(in) :: Dtset
1198 : type(em1params_t),intent(in) :: Ep
1199 : type(gsphere_t),intent(in) :: Gsph_epsG0
1200 : type(littlegroup_t),intent(in) :: Ltg_q
1201 : type(Pawang_type),intent(in) :: Pawang
1202 : type(Pseudopotential_type),intent(in) :: Psps
1203 : type(wfdgw_t),target,intent(inout) :: Wfd,Wfdf
1204 : !arrays
1205 : integer,intent(in) :: ktabr(nfftot_gw,Kmesh%nbz),ktabrf(nfftf_tot*Dtset%pawcross,Kmesh%nbz)
1206 : integer,intent(in) :: ngfft_gw(18),ngfftf(18)
1207 : real(dp),intent(in) :: qpoint(3)
1208 : real(dp),intent(out) :: chi0_sumrule(Ep%npwe)
1209 : complex(gwp),intent(out) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega)
1210 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
1211 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
1212 : type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
1213 : type(plowannier_type),intent(inout) :: wan
1214 : type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom)
1215 :
1216 : !Local variables ------------------------------
1217 : !scalars
1218 : integer,parameter :: tim_fourdp1 = 1, two_poles = 2, one_pole = 1, ndat1 = 1
1219 : integer :: bandinf,bandsup,dim_rtwg,band1,band2,ierr,band1c,band2c
1220 : integer :: ig1,ig2,iat1,iat2,iat,ik_bz,ik_ibz,ikmq_bz,ikmq_ibz
1221 : integer :: io,iomegal,iomegar,ispinor1,ispinor2,isym_k,itypatcor,nfft,il1,il2
1222 : integer :: isym_kmq,itim_k,itim_kmq,m1,m2,my_wl,my_wr,size_chi0
1223 : integer :: nfound,nkpt_summed,nspinor,nsppol,mband
1224 : integer :: comm,gw_mgfft,use_padfft,gw_fftalga,lcor,mgfftf,use_padfftf
1225 : integer :: my_nbbp,my_nbbpks,spin,nbmax,dummy
1226 : real(dp) :: cpu_all, wall_all, gflops_all, cpu_k, wall_k, gflops_k
1227 : real(dp) :: deltaeGW_b1kmq_b2k,deltaeGW_enhigh_b2k,deltaf_b1kmq_b2k
1228 : real(dp) :: e_b1_kmq,en_high,fac,fac2,fac3,f_b1_kmq,factor,max_rest,min_rest,my_max_rest
1229 : real(dp) :: my_min_rest,numerator,spin_fact,weight,wl,wr
1230 : real(dp) :: gw_gsq,memreq
1231 : complex(dp) :: ph_mkmqt,ph_mkt
1232 : complex(gwp) :: local_czero_gw
1233 : logical :: qzero,isirred_k,isirred_kmq,luwindow,is_metallic, print_time
1234 : character(len=500) :: msg,allup
1235 776 : type(gsphere_t) :: Gsph_FFT
1236 : !arrays
1237 : integer :: G0(3),umklp_k(3),umklp_kmq(3), ucrpa_bands(2)
1238 2328 : integer :: wtk_ltg(Kmesh%nbz),got(Wfd%nproc)
1239 776 : integer,allocatable :: tabr_k(:),tabr_kmq(:),tabrf_k(:),tabrf_kmq(:)
1240 776 : integer,allocatable :: igfftepsG0(:),gspfft_igfft(:),igfftepsG0f(:)
1241 776 : integer,allocatable :: gw_gfft(:,:),gw_gbound(:,:),dummy_gbound(:,:),gboundf(:,:)
1242 776 : integer,allocatable :: bbp_ks_distrb(:,:,:,:)
1243 : real(dp) :: kbz(3),kmq_bz(3),spinrot_k(4),spinrot_kmq(4),q0(3),tsec(2)
1244 776 : real(dp),contiguous, pointer :: qp_eig(:,:,:),qp_occ(:,:,:)
1245 776 : real(dp),allocatable :: omegasf(:)
1246 776 : complex(dp),allocatable :: green_enhigh_w(:),green_w(:),kkweight(:,:)
1247 776 : complex(gwp),allocatable :: sf_chi0(:,:,:),rhotwg(:)
1248 776 : complex(gwp),allocatable :: ur1_kmq_ibz(:),ur2_k_ibz(:),wfwfg(:)
1249 776 : complex(gwp),allocatable :: usr1_kmq(:),ur2_k(:)
1250 776 : complex(gwp),allocatable :: ur_ae1(:),ur_ae_onsite1(:),ur_ps_onsite1(:)
1251 776 : complex(gwp),allocatable :: ur_ae2(:),ur_ae_onsite2(:),ur_ps_onsite2(:)
1252 776 : complex(dp), allocatable :: coeffW_BZ(:,:,:,:,:,:)
1253 776 : logical,allocatable :: bbp_mask(:,:)
1254 776 : type(pawcprj_type),allocatable :: Cprj1_kmq(:,:),Cprj2_k(:,:)
1255 776 : type(pawpwij_t),allocatable :: Pwij(:),Pwij_fft(:)
1256 : !************************************************************************
1257 :
1258 : DBG_ENTER("COLL")
1259 :
1260 776 : call timab(331,1,tsec) ! cchi0
1261 776 : call cwtime(cpu_all, wall_all, gflops_all, "start")
1262 :
1263 776 : nsppol = Wfd%nsppol; nspinor = Wfd%nspinor
1264 776 : is_metallic = qp_ebands%has_metal_scheme()
1265 :
1266 776 : ucrpa_bands(1)=dtset%ucrpa_bands(1)
1267 776 : ucrpa_bands(2)=dtset%ucrpa_bands(2)
1268 776 : luwindow=.false.
1269 776 : if(abs(dtset%ucrpa_window(1)+1_dp)>tol8.or.(abs(dtset%ucrpa_window(2)+1_dp)>tol8)) then
1270 0 : luwindow=.true.
1271 : endif
1272 : !write(6,*)"ucrpa_bands",ucrpa_bands; write(6,*)"ucrpa_window",dtset%ucrpa_window; write(6,*)"luwindow",luwindow
1273 :
1274 : ! For cRPA calculation of U: read forlb.ovlp
1275 776 : if(dtset%ucrpa>=1 .AND. dtset%plowan_compute <10) then
1276 : call read_plowannier(Cryst,bandinf,bandsup,coeffW_BZ,itypatcor,Kmesh,lcor,luwindow,&
1277 0 : nspinor,nsppol,pawang,dtset%prtvol,ucrpa_bands)
1278 : endif
1279 : ! End of reading forlb.ovlp
1280 :
1281 3098 : if ( ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst,Psps,ngfft_gw)
1282 3104 : gw_mgfft = MAXVAL(ngfft_gw(1:3))
1283 776 : gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10)
1284 :
1285 776 : if (Dtset%pawcross==1) mgfftf = MAXVAL(ngfftf(1:3))
1286 :
1287 : ! == Copy some values ===
1288 776 : comm = Wfd%comm
1289 776 : mband = Wfd%mband
1290 776 : nfft = Wfd%nfft
1291 776 : ABI_CHECK(Wfd%nfftot==nfftot_gw,"Wrong nfftot_gw")
1292 :
1293 776 : dim_rtwg=1 !; if (nspinor==2) dim_rtwg=2 ! can reduce size depending on Ep%nI and Ep%nj
1294 776 : size_chi0 = Ep%npwe*Ep%nI*Ep%npwe*Ep%nJ*Ep%nomega
1295 :
1296 776 : qp_eig => qp_ebands%eig; qp_occ => qp_ebands%occ
1297 :
1298 : ! Initialize the completeness correction
1299 776 : if (Ep%gwcomp==1) then
1300 170 : en_high=MAXVAL(qp_eig(Ep%nbnds,:,:)) + Ep%gwencomp
1301 20 : write(msg,'(a,f8.2,a)')' Using completeness correction with the energy ',en_high*Ha_eV,' [eV]'
1302 20 : call wrtout(std_out, msg)
1303 :
1304 : ! Allocation of wfwfg and green_enhigh_w moved inside openmp loop
1305 : ! Init the largest G-sphere contained in the FFT box for the wavefunctions.
1306 20 : call Gsph_FFT%in_fftbox(Cryst,Wfd%ngfft)
1307 :
1308 : !call Gsph_FFT%print(unit=std_out,prtvol=10)
1309 :
1310 60 : ABI_MALLOC(gspfft_igfft,(Gsph_FFT%ng))
1311 80 : ABI_MALLOC(dummy_gbound,(2*gw_mgfft+8,2))
1312 :
1313 : ! Mapping between G-sphere and FFT box.
1314 20 : call Gsph_FFT%fft_tabs((/0,0,0/),Wfd%mgfft,Wfd%ngfft,dummy,dummy_gbound,gspfft_igfft)
1315 20 : ABI_FREE(dummy_gbound)
1316 :
1317 20 : if (Psps%usepaw==1) then ! * Prepare the onsite contributions on the GW FFT mesh.
1318 0 : ABI_MALLOC(gw_gfft,(3,nfft))
1319 0 : q0=zero
1320 0 : call get_gfft(ngfft_gw,q0,Cryst%gmet,gw_gsq,gw_gfft) ! Get the set of plane waves in the FFT Box.
1321 0 : ABI_MALLOC(Pwij_fft,(Psps%ntypat))
1322 0 : call pawpwij_init(Pwij_fft,nfft,(/zero,zero,zero/),gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
1323 : end if
1324 : end if
1325 :
1326 : ! Setup weights (2 for spin unpolarized system, 1 for polarized).
1327 : ! spin_fact is used to normalize the occupation factors to one. Consider also the AFM case.
1328 758 : select case (nsppol)
1329 : case (1)
1330 758 : weight = two / Kmesh%nbz; spin_fact = half
1331 758 : if (Wfd%nspden==2) then
1332 9 : weight = one / Kmesh%nbz; spin_fact = half
1333 : end if
1334 758 : if (nspinor==2) then
1335 4 : weight = one / Kmesh%nbz; spin_fact = one
1336 : end if
1337 : case (2)
1338 18 : weight = one / Kmesh%nbz; spin_fact = one
1339 : case default
1340 776 : ABI_BUG("Wrong nsppol")
1341 : end select
1342 :
1343 : ! Weight for points in the IBZ_q.
1344 70814 : wtk_ltg(:) = 1
1345 776 : if (Ep%symchi == 1) then
1346 67401 : do ik_bz=1,Ltg_q%nbz
1347 66696 : wtk_ltg(ik_bz) = 0
1348 66696 : if (Ltg_q%ibzq(ik_bz) /= 1) CYCLE ! Only k-points in the IBZ_q.
1349 4251159 : wtk_ltg(ik_bz)=SUM(Ltg_q%wtksym(:,:,ik_bz))
1350 : end do
1351 : end if
1352 :
1353 : write(msg,'(a,i2,2a,i2)')&
1354 776 : ' Using spectral method for the imaginary part = ',Ep%spmeth,ch10,&
1355 1552 : ' Using symmetries to sum only over the IBZ_q = ',Ep%symchi
1356 776 : call wrtout(std_out, msg)
1357 :
1358 776 : if (use_tr) then
1359 : ! Special care has to be taken in metals and/or spin dependent systems
1360 : ! as Wfs_val might contain unoccupied states.
1361 718 : call wrtout(std_out,' Using faster algorithm based on time reversal symmetry.')
1362 : else
1363 58 : call wrtout(std_out,' Using slow algorithm without time reversal symmetry.')
1364 : end if
1365 :
1366 : ! TODO this table can be calculated for each k-point
1367 2168 : my_nbbpks=0; allup="All"; got=0
1368 4656 : ABI_MALLOC(bbp_ks_distrb,(mband,mband,Kmesh%nbz,nsppol))
1369 : call wrtout(std_out, sjoin(' Memory needed for bbp_ks_distrb: ', &
1370 776 : ftoa(four*mband**2*Kmesh%nbz*nsppol*b2Mb, fmt="f8.1"), ' [Mb] <<< MEM'))
1371 :
1372 3104 : ABI_MALLOC(bbp_mask,(mband, mband))
1373 :
1374 1570 : do spin=1,nsppol
1375 72760 : do ik_bz=1,Kmesh%nbz
1376 71190 : if (Ep%symchi == 1) then
1377 67272 : if (Ltg_q%ibzq(ik_bz) /= 1) CYCLE ! Only IBZ_q
1378 : end if
1379 :
1380 : ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz.
1381 33707 : call kmesh%get_BZ_item(ik_bz, kbz, ik_ibz, isym_k, itim_k)
1382 :
1383 : ! Get index of k-q in the BZ, stop if not found as the weight=one/nkbz is not correct.
1384 33707 : call kmesh%get_BZ_diff(kbz, qpoint, ikmq_bz, g0, nfound)
1385 33707 : ABI_CHECK(nfound == 1, "Check kmesh")
1386 :
1387 : ! Get ikmq_ibz, non-symmorphic phase, ph_mkmqt, and symmetries from ikmq_bz.
1388 33707 : call kmesh%get_BZ_item(ikmq_bz, kmq_bz, ikmq_ibz, isym_kmq, itim_kmq)
1389 :
1390 : call chi0_bbp_mask(ikmq_ibz, ik_ibz, spin, spin_fact, use_tr, &
1391 33707 : ep%gwcomp, ep%spmeth, ep%nbnds, mband, qp_ebands, bbp_mask)
1392 :
1393 33707 : call wfd%distribute_kb_kpbp(ikmq_ibz,ik_ibz,spin,allup,my_nbbp,bbp_ks_distrb(:,:,ik_bz,spin),got=got,bbp_mask=bbp_mask)
1394 71984 : my_nbbpks = my_nbbpks + my_nbbp
1395 : end do
1396 : end do
1397 :
1398 776 : ABI_FREE(bbp_mask)
1399 :
1400 776 : write(msg,'(a,i0,a)')" Will sum ",my_nbbpks," (b,b',k,s) states in chi0."
1401 776 : call wrtout(std_out, msg)
1402 :
1403 776 : if (Psps%usepaw==1) then
1404 97 : ABI_MALLOC(Pwij,(Psps%ntypat))
1405 23 : call pawpwij_init(Pwij,Ep%npwepG0,qpoint,Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
1406 : ! Allocate statements moved to inside openmp loop
1407 : end if
1408 :
1409 1504 : SELECT CASE (Ep%spmeth)
1410 : CASE (0)
1411 728 : call wrtout(std_out,' Calculating chi0(q,omega,G,G")')
1412 : ! Allocation of green_w moved inside openmp loop
1413 :
1414 : CASE (1, 2)
1415 48 : call wrtout(std_out,' Calculating Im chi0(q,omega,G,G")')
1416 :
1417 : ! Find Max and min resonant transitions for this q, report also treated by this proc.
1418 : call make_transitions(Wfd,1,Ep%nbnds,nbvw,nsppol,Ep%symchi,Cryst%timrev,GW_TOL_DOCC,&
1419 48 : max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,qp_eig,qp_occ,qpoint,bbp_ks_distrb)
1420 : !
1421 : ! Calculate frequency dependent weights for Hilbert transform.
1422 144 : ABI_MALLOC(omegasf,(Ep%nomegasf))
1423 192 : ABI_MALLOC(kkweight,(Ep%nomegasf,Ep%nomega))
1424 : !my_wl=1; my_wr=Ep%nomegasf
1425 : call setup_spectral(Ep%nomega,Ep%omega,Ep%nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,&
1426 48 : 0,Ep%zcut,zero,my_wl,my_wr,kkweight)
1427 :
1428 48 : if (.not. use_tr) then
1429 0 : ABI_BUG('spectral method requires time-reversal')
1430 : end if
1431 :
1432 48 : memreq = two*gwp*Ep%npwe**2*(my_wr-my_wl+1)*b2Gb
1433 48 : write(msg,'(a,f10.4,a)')' memory required per spectral point: ',two*gwp*Ep%npwe**2*b2Mb,' [Mb]'
1434 48 : call wrtout(std_out,msg)
1435 48 : write(msg,'(a,f10.4,a)')' memory required by sf_chi0: ',memreq,' [Gb]'
1436 48 : call wrtout(std_out,msg)
1437 48 : if (memreq > two) then
1438 0 : ABI_WARNING(' Memory required for sf_chi0 is larger than 2.0 Gb!')
1439 : end if
1440 240 : ABI_MALLOC_OR_DIE(sf_chi0,(Ep%npwe,Ep%npwe,my_wl:my_wr), ierr)
1441 16608048 : sf_chi0=czero_gw
1442 :
1443 : CASE DEFAULT
1444 776 : ABI_BUG("Wrong spmeth")
1445 : END SELECT
1446 :
1447 776 : nkpt_summed=Kmesh%nbz
1448 776 : if (Ep%symchi == 1) then
1449 705 : nkpt_summed = Ltg_q%nibz_ltg
1450 1410 : call Ltg_q%print([std_out], prtvol=Dtset%prtvol)
1451 : end if
1452 :
1453 776 : write(msg,'(a,i0,a)')' Calculation status: ',nkpt_summed,' k-points to be completed'
1454 776 : call wrtout(std_out, msg)
1455 776 : call pstat_proc%print(_PSTAT_ARGS_)
1456 :
1457 : ! ============================================
1458 : ! === Begin big fat loop over transitions ===
1459 : ! ============================================
1460 14722809 : chi0=czero_gw; chi0_sumrule=zero
1461 :
1462 : ! === Loop on spin to calculate trace $\chi_{up,up}+\chi_{down,down}$ ===
1463 : ! Only $\chi_{up,up} for AFM.
1464 1570 : do spin=1,nsppol
1465 5666 : if (ALL(bbp_ks_distrb(:,:,:,spin) /= Wfd%my_rank)) CYCLE
1466 :
1467 : ! Allocation of arrays that are private to loop
1468 794 : if (Ep%gwcomp==1) then
1469 60 : ABI_MALLOC(wfwfg,(nfft*nspinor**2))
1470 : end if
1471 794 : if (Ep%gwcomp==1) then
1472 60 : ABI_MALLOC(green_enhigh_w,(Ep%nomega))
1473 : end if
1474 794 : if (Ep%spmeth==0) then
1475 2238 : ABI_MALLOC(green_w,(Ep%nomega))
1476 : end if
1477 794 : if (Psps%usepaw==1) then
1478 182 : ABI_MALLOC(Cprj2_k ,(Cryst%natom,nspinor))
1479 23 : call pawcprj_alloc(Cprj2_k, 0,Wfd%nlmn_atm)
1480 159 : ABI_MALLOC(Cprj1_kmq,(Cryst%natom,nspinor))
1481 23 : call pawcprj_alloc(Cprj1_kmq,0,Wfd%nlmn_atm)
1482 23 : if (Dtset%pawcross==1) then
1483 0 : ABI_MALLOC(ur_ae1,(nfftf_tot*nspinor))
1484 0 : ABI_MALLOC(ur_ae_onsite1,(nfftf_tot*nspinor))
1485 0 : ABI_MALLOC(ur_ps_onsite1,(nfftf_tot*nspinor))
1486 0 : ABI_MALLOC(ur_ae2,(nfftf_tot*nspinor))
1487 0 : ABI_MALLOC(ur_ae_onsite2,(nfftf_tot*nspinor))
1488 0 : ABI_MALLOC(ur_ps_onsite2,(nfftf_tot*nspinor))
1489 0 : ABI_MALLOC(igfftepsG0f,(Ep%npwepG0))
1490 0 : ABI_MALLOC(tabrf_k,(nfftf_tot))
1491 0 : ABI_MALLOC(tabrf_kmq,(nfftf_tot))
1492 : end if
1493 : end if
1494 :
1495 2382 : ABI_MALLOC(rhotwg,(Ep%npwepG0*dim_rtwg))
1496 2382 : ABI_MALLOC(tabr_k,(nfft))
1497 1588 : ABI_MALLOC(tabr_kmq,(nfft))
1498 2382 : ABI_MALLOC(ur1_kmq_ibz,(nfft*nspinor))
1499 1588 : ABI_MALLOC(ur2_k_ibz,(nfft*nspinor))
1500 1588 : ABI_MALLOC(usr1_kmq,(nfft*nspinor))
1501 1588 : ABI_MALLOC(ur2_k, (nfft*nspinor))
1502 2382 : ABI_MALLOC(igfftepsG0,(Ep%npwepG0))
1503 :
1504 : ! Loop over k-points in the BZ.
1505 71984 : do ik_bz=1,Kmesh%nbz
1506 :
1507 71190 : if (Ep%symchi==1) then
1508 67272 : if (Ltg_q%ibzq(ik_bz) /= 1) CYCLE ! Only IBZ_q
1509 : end if
1510 :
1511 201094 : if (ALL(bbp_ks_distrb(:,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE
1512 :
1513 33707 : print_time = wfd%my_rank == 0 .and. (ik_bz < LOG_MODK .or. mod(ik_bz, LOG_MODK) == 0)
1514 6122 : if (print_time) call cwtime(cpu_k, wall_k, gflops_k, "start")
1515 :
1516 : ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz.
1517 33707 : call kmesh%get_BZ_item(ik_bz, kbz, ik_ibz, isym_k, itim_k, ph_mkt, umklp_k, isirred_k)
1518 :
1519 33707 : call kmesh%get_BZ_diff(kbz, qpoint, ikmq_bz, G0, nfound)
1520 33707 : if (nfound == 0) then
1521 0 : ABI_ERROR("Cannot find kbz - qpoint in Kmesh")
1522 : end if
1523 :
1524 : ! Get ikmq_ibz, non-symmorphic phase, ph_mkmqt, and symmetries from ikmq_bz.
1525 33707 : call kmesh%get_BZ_item(ikmq_bz, kmq_bz, ikmq_ibz, isym_kmq, itim_kmq, ph_mkmqt, umklp_kmq, isirred_kmq)
1526 :
1527 : !BEGIN DEBUG
1528 : !if (ANY(umklp_k /=0)) then
1529 : ! write(msg,'(a,3i2)')" umklp_k /= 0 ",umklp_k
1530 : ! ABI_ERROR(msg)
1531 : !end if
1532 : !if (ANY( g0 /= -umklp_kmq + umklp_k) ) then
1533 : !if (ANY( g0 /= -umklp_kmq ) ) then
1534 : ! write(msg,'(a,3(1x,3i2))')" g0 /= -umklp_kmq + umklp_k ",g0, umklp_kmq, umklp_k
1535 : ! ABI_ERROR(msg)
1536 : !end if
1537 : !g0 = -umklp_k + umklp_kmq
1538 : !g0 = +umklp_k - umklp_kmq
1539 : !if (ANY (ABS(g0) > Ep%mg0) ) then
1540 : ! write(msg,'(a,6(1x,i0))')" ABS(g0) > Ep%mg0 ",g0,Ep%mg0
1541 : ! ABI_ERROR(msg)
1542 : !end if
1543 : !END DEBUG
1544 :
1545 : ! Copy tables for rotated FFT points
1546 212390575 : tabr_k(:) =ktabr(:,ik_bz)
1547 168535 : spinrot_k(:)=Cryst%spinrot(:,isym_k)
1548 :
1549 212390575 : tabr_kmq(:)=ktabr(:,ikmq_bz)
1550 168535 : spinrot_kmq(:)=Cryst%spinrot(:,isym_kmq)
1551 :
1552 33707 : if (Dtset%pawcross==1) then
1553 0 : tabrf_k(:) =ktabrf(:,ik_bz)
1554 0 : tabrf_kmq(:)=ktabrf(:,ikmq_bz)
1555 : end if
1556 : !
1557 : ! Tables for the FFT of the oscillators.
1558 : ! a) FFT index of G-G0.
1559 : ! b) gw_gbound table for the zero-padded FFT performed in rhotwg.
1560 134828 : ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2))
1561 33707 : call Gsph_epsG0%fft_tabs(g0,gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igfftepsG0)
1562 33707 : if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
1563 : !use_padfft = 0
1564 33707 : if (use_padfft==0) then
1565 107 : ABI_FREE(gw_gbound)
1566 214 : ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft))
1567 : end if
1568 :
1569 33707 : if (Dtset%pawcross==1) then
1570 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2))
1571 0 : call Gsph_epsG0%fft_tabs(g0,mgfftf,ngfftf,use_padfftf,gboundf,igfftepsG0f)
1572 0 : if (ANY(gw_fftalga == [2, 4])) use_padfftf=0
1573 0 : if (use_padfftf==0) then
1574 0 : ABI_FREE(gboundf)
1575 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf))
1576 : end if
1577 : end if
1578 :
1579 33707 : nbmax=Ep%nbnds
1580 610049 : do band1=1,nbmax ! Loop over "conduction" states.
1581 6776356 : if (ALL(bbp_ks_distrb(band1,:,ik_bz,spin) /= Wfd%my_rank)) CYCLE
1582 :
1583 344287 : call wfd%get_ur(band1,ikmq_ibz,spin,ur1_kmq_ibz)
1584 :
1585 344287 : if (Psps%usepaw==1) then
1586 23748 : call wfd%get_cprj(band1,ikmq_ibz,spin,Cryst,Cprj1_kmq,sorted=.FALSE.)
1587 23748 : call paw_symcprj(ikmq_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_kmq)
1588 23748 : if (Dtset%pawcross==1) then
1589 0 : call wfdf%paw_get_aeur(band1,ikmq_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae1,ur_ae_onsite1,ur_ps_onsite1)
1590 : end if
1591 : end if
1592 :
1593 344287 : e_b1_kmq = qp_eig(band1,ikmq_ibz,spin)
1594 344287 : f_b1_kmq = qp_occ(band1,ikmq_ibz,spin)
1595 :
1596 11836304 : do band2=1,nbmax ! Loop over "valence" states.
1597 : if (luwindow.AND.dtset%ucrpa==1 &
1598 : .AND.((qp_ebands%eig(band1,ik_ibz ,spin)-qp_ebands%fermie)<=dtset%ucrpa_window(2)) &
1599 : .AND.((qp_ebands%eig(band1,ik_ibz ,spin)-qp_ebands%fermie)>=dtset%ucrpa_window(1)) &
1600 : .AND.((qp_ebands%eig(band2,ikmq_ibz,spin)-qp_ebands%fermie)<=dtset%ucrpa_window(2)) &
1601 11458310 : .AND.((qp_ebands%eig(band2,ikmq_ibz,spin)-qp_ebands%fermie)>=dtset%ucrpa_window(1))) CYCLE
1602 :
1603 11458310 : if (bbp_ks_distrb(band1,band2,ik_bz,spin) /= Wfd%my_rank) CYCLE
1604 :
1605 1503924 : deltaf_b1kmq_b2k = spin_fact * (f_b1_kmq - qp_occ(band2, ik_ibz, spin))
1606 :
1607 1503924 : if (Ep%gwcomp==0) then ! Skip negligible transitions.
1608 1470484 : if (abs(deltaf_b1kmq_b2k) < GW_TOL_DOCC) CYCLE
1609 :
1610 : else
1611 : ! When the completeness correction is used,
1612 : ! we need to also consider transitions with vanishing deltaf
1613 : !if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC) CYCLE
1614 : !
1615 : ! Rangel This is to compute chi correctly when using the extrapolar method
1616 33440 : if (qp_occ(band2,ik_ibz,spin) < GW_TOL_DOCC .and. (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC .or. band1<band2)) CYCLE
1617 : end if
1618 :
1619 1503924 : deltaeGW_b1kmq_b2k = e_b1_kmq - qp_eig(band2,ik_ibz,spin)
1620 :
1621 1503924 : call wfd%get_ur(band2,ik_ibz,spin,ur2_k_ibz)
1622 :
1623 1503924 : if (Psps%usepaw==1) then
1624 95832 : call wfd%get_cprj(band2,ik_ibz,spin,Cryst,Cprj2_k,sorted=.FALSE.)
1625 95832 : call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj2_k)
1626 95832 : if (Dtset%pawcross==1) then
1627 0 : call wfdf%paw_get_aeur(band2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae2,ur_ae_onsite2,ur_ps_onsite2)
1628 : end if
1629 : end if
1630 :
1631 1251604 : SELECT CASE (Ep%spmeth)
1632 : CASE (0)
1633 : ! Standard Adler-Wiser expression.
1634 : ! Add the small imaginary of the Time-Ordered RF only for non-zero real omega ! FIXME What about metals?
1635 1251604 : if (.not. use_tr) then
1636 : ! Have to sum over all possible resonant and anti-resonant transitions.
1637 1148200 : do io=1,Ep%nomega
1638 1148200 : green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,one_pole)
1639 : end do
1640 :
1641 : else
1642 1081776 : if (Ep%gwcomp==0) then ! cannot be completely skipped in case of completeness correction
1643 1048336 : if (band1 < band2) CYCLE ! Here we GAIN a factor ~2
1644 : end if
1645 :
1646 4688106 : do io=1,Ep%nomega
1647 : !Rangel: In metals, the intra-band transitions term does not contain the antiresonant part
1648 : !green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0)
1649 3606330 : if (band1 == band2) then
1650 10288 : green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,one_pole)
1651 : else
1652 3596042 : green_w(io) = g0g0w(Ep%omega(io),deltaf_b1kmq_b2k,deltaeGW_b1kmq_b2k,Ep%zcut,GW_TOL_W0,two_poles)
1653 : end if
1654 :
1655 4688106 : if (Ep%gwcomp==1) then ! Calculate the completeness correction
1656 86080 : numerator= -spin_fact*qp_occ(band2,ik_ibz,spin)
1657 86080 : deltaeGW_enhigh_b2k = en_high-qp_eig(band2,ik_ibz,spin)
1658 :
1659 86080 : if (REAL(Ep%omega(io))<GW_TOL_W0) then ! Completeness correction is NOT valid for real frequencies
1660 86080 : green_enhigh_w(io) = g0g0w(Ep%omega(io),numerator,deltaeGW_enhigh_b2k,Ep%zcut,GW_TOL_W0,two_poles)
1661 : else
1662 0 : green_enhigh_w(io) = local_czero_gw
1663 : end if
1664 : !
1665 : !Rangel Correction for metals
1666 : !if (deltaf_b1kmq_b2k<0.d0) then
1667 86080 : if (band1>=band2 .and. ABS(deltaf_b1kmq_b2k) > GW_TOL_DOCC ) then
1668 70388 : green_w(io)= green_w(io) - green_enhigh_w(io)
1669 : else ! Disregard green_w, since it is already accounted for through the time-reversal
1670 15692 : green_w(io)= - green_enhigh_w(io)
1671 : end if
1672 : end if !gwcomp==1
1673 : end do !io
1674 :
1675 1081776 : if (Ep%gwcomp==1.and.band1==band2) then
1676 : ! Add the "delta part" of the extrapolar method. TODO does not work for spinor
1677 1632 : call calc_wfwfg(tabr_k,itim_k,spinrot_k,nfft,nspinor,ngfft_gw,ur2_k_ibz,ur2_k_ibz,wfwfg)
1678 :
1679 1632 : if (Psps%usepaw==1) then
1680 0 : call paw_rho_tw_g(cryst,Pwij_fft, nfft,dim_rtwg,nspinor,gw_gfft,Cprj2_k,Cprj2_k,wfwfg)
1681 :
1682 : ! Add PAW cross term
1683 0 : if (Dtset%pawcross==1) then
1684 : call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,&
1685 : ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_kmq,tabrf_kmq,ph_mkmqt,spinrot_kmq,&
1686 0 : ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k ,tabrf_k ,ph_mkt ,spinrot_k,dim_rtwg,wfwfg)
1687 : end if
1688 : end if
1689 :
1690 1632 : qzero=.FALSE.
1691 : call completechi0_deltapart(ik_bz,qzero,Ep%symchi,Ep%npwe,Gsph_FFT%ng,Ep%nomega,nspinor,&
1692 1632 : nfft,ngfft_gw,gspfft_igfft,gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0)
1693 :
1694 : end if
1695 : end if ! use_tr
1696 :
1697 : CASE (1, 2)
1698 : ! Spectral method, WARNING time-reversal here is always assumed!
1699 252320 : if (deltaeGW_b1kmq_b2k<0) CYCLE
1700 1756244 : call approxdelta(Ep%nomegasf,omegasf,deltaeGW_b1kmq_b2k,Ep%spsmear,iomegal,iomegar,wl,wr,Ep%spmeth)
1701 : END SELECT
1702 :
1703 : ! Form rho-twiddle(r)=u^*_{b1,kmq_bz}(r) u_{b2,kbz}(r) and its FFT transform.
1704 : call rho_tw_g(nspinor,Ep%npwepG0,nfft,ndat1,ngfft_gw,1,use_padfft,igfftepsG0,gw_gbound,&
1705 : ur1_kmq_ibz,itim_kmq,tabr_kmq,ph_mkmqt,spinrot_kmq,&
1706 1503924 : ur2_k_ibz, itim_k ,tabr_k ,ph_mkt ,spinrot_k,dim_rtwg,rhotwg)
1707 :
1708 1503924 : if (Psps%usepaw==1) then
1709 : ! Add PAW on-site contribution, projectors are already in the BZ.
1710 95832 : call paw_rho_tw_g(cryst, Pwij, Ep%npwepG0,dim_rtwg,nspinor,Gsph_epsG0%gvec,Cprj1_kmq,Cprj2_k,rhotwg)
1711 :
1712 : ! Add PAW cross term
1713 95832 : if (Dtset%pawcross==1) then
1714 : call paw_cross_rho_tw_g(nspinor,Ep%npwepG0,nfftf_tot,ngfftf,1,use_padfftf,igfftepsG0f,gboundf,&
1715 : ur_ae1,ur_ae_onsite1,ur_ps_onsite1,itim_kmq,tabrf_kmq,ph_mkmqt,spinrot_kmq,&
1716 0 : ur_ae2,ur_ae_onsite2,ur_ps_onsite2,itim_k ,tabrf_k ,ph_mkt ,spinrot_k,dim_rtwg,rhotwg)
1717 : end if
1718 : end if
1719 :
1720 1251604 : SELECT CASE (Ep%spmeth)
1721 :
1722 : CASE (0) ! Adler-Wiser.
1723 : !debug if(dtset%ucrpa==2) then
1724 1251604 : if(dtset%ucrpa>=1.and..not.luwindow) then
1725 0 : fac=one
1726 0 : fac2=one
1727 0 : fac3=one
1728 0 : m1=-1
1729 0 : m2=-1
1730 0 : call flush_unit(std_out)
1731 0 : if(dtset%ucrpa<=2) then
1732 : if ( band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)&
1733 0 : & .AND.band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then
1734 0 : if (dtset%plowan_compute >=10) then
1735 0 : band1c=band1-wan%bandi_wan+1
1736 0 : band2c=band2-wan%bandi_wan+1
1737 0 : do iat1=1, wan%natom_wan
1738 0 : do iat2=1, wan%natom_wan
1739 0 : do ispinor1=1,wan%nspinor
1740 0 : do ispinor2=1,wan%nspinor
1741 0 : do il1=1,wan%nbl_atom_wan(iat1)
1742 0 : do il2=1,wan%nbl_atom_wan(iat2)
1743 0 : do m1=1,2*wan%latom_wan(iat1)%lcalc(il1)+1
1744 0 : do m2=1,2*wan%latom_wan(iat2)%lcalc(il2)+1
1745 : fac=fac - real(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1)*&
1746 : &conjg(wan%psichi(ik_bz,band1c,iat1)%atom(il1)%matl(m1,spin,ispinor1))*&
1747 : &wan%psichi(ikmq_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)*&
1748 0 : &conjg(wan%psichi(ikmq_bz,band2c,iat2)%atom(il2)%matl(m2,spin,ispinor2)))
1749 : enddo !m2
1750 : enddo !m1
1751 : enddo !il2
1752 : enddo !il1
1753 : enddo !ispinor2
1754 : enddo !isspinor1
1755 : enddo !iat2
1756 : enddo !iat1
1757 : else !plowan_compute>=10
1758 0 : do iat=1, cryst%nattyp(itypatcor)
1759 0 : do ispinor1=1,nspinor
1760 0 : do ispinor2=1,nspinor
1761 0 : do m1=1,2*lcor+1
1762 0 : do m2=1,2*lcor+1
1763 : fac=fac - real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
1764 : & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1))* &
1765 : & coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor2,m2)*&
1766 0 : & conjg(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor2,m2)))
1767 : enddo !m2
1768 : enddo !m1
1769 : enddo !ispinor2
1770 : enddo !ispinor1
1771 : enddo !iat
1772 : endif !plowan_compute>=10
1773 0 : if(dtset%ucrpa==1) fac=zero
1774 : endif
1775 : else if (dtset%ucrpa>=3) then
1776 0 : if (band1<=ucrpa_bands(2).AND.band1>=ucrpa_bands(1)) then
1777 0 : do iat=1, cryst%nattyp(itypatcor)
1778 0 : do ispinor1=1,nspinor
1779 0 : do m1=1,2*lcor+1
1780 : fac2=fac2-real(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)*&
1781 0 : & conjg(coeffW_BZ(iat,spin,band1,ik_bz,ispinor1,m1)))
1782 : enddo
1783 : enddo
1784 : enddo
1785 0 : if(dtset%ucrpa==4) fac2=zero
1786 : endif
1787 0 : if (band2<=ucrpa_bands(2).AND.band2>=ucrpa_bands(1)) then
1788 0 : do iat=1, cryst%nattyp(itypatcor)
1789 0 : do ispinor1=1,nspinor
1790 0 : do m1=1,2*lcor+1
1791 : fac3=fac3-real(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor1,m1)*&
1792 0 : & conjg(coeffW_BZ(iat,spin,band2,ikmq_bz,ispinor1,m1)))
1793 : enddo
1794 : enddo
1795 : enddo
1796 0 : if(dtset%ucrpa==4) fac3=zero
1797 : endif
1798 0 : fac=real(fac2*fac3)
1799 : endif
1800 :
1801 : ! if(dtset%prtvol>=10) write(6,'(6i3,e15.5,a)') ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q=/0"
1802 : ! if(dtset%prtvol>=10.and.abs(fac-one)>0.00001) &
1803 : !& write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q/=0"
1804 : ! if(dtset%prtvol>=10) write(6,'(a,6i4,e15.5,a)') "*****FAC*********",ik_bz,ikmq_bz,band1,band2,m1,m2,fac," q/=0"
1805 0 : green_w=green_w*fac
1806 : endif
1807 :
1808 1251604 : call assemblychi0_sym(is_metallic,ik_bz,nspinor,Ep,Ltg_q,green_w,Ep%npwepG0,rhotwg,Gsph_epsG0,chi0)
1809 :
1810 : CASE (1, 2)
1811 : ! Spectral method (not yet adapted for nspinor=2)
1812 : call assemblychi0sf(ik_bz,Ep%symchi,Ltg_q,Ep%npwepG0,Ep%npwe,rhotwg,Gsph_epsG0,&
1813 252320 : & deltaf_b1kmq_b2k,my_wl,iomegal,wl,my_wr,iomegar,wr,Ep%nomegasf,sf_chi0)
1814 :
1815 : CASE DEFAULT
1816 1503924 : ABI_BUG("Wrong spmeth")
1817 : END SELECT
1818 :
1819 : ! Accumulating the sum rule on chi0. Eq. (5.284) in G.D. Mahan Many-Particle Physics 3rd edition. [[cite:Mahan2000]]
1820 : ! TODO Does not work with spinor
1821 1503924 : factor=spin_fact*qp_occ(band2,ik_ibz,spin)
1822 : call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_b1kmq_b2k,&
1823 1503924 : & Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule)
1824 :
1825 : ! Include also the completeness correction in the sum rule
1826 2080266 : if (Ep%gwcomp==1) then
1827 33440 : factor=-spin_fact*qp_occ(band2,ik_ibz,spin)
1828 : call accumulate_chi0sumrule(ik_bz,Ep%symchi,Ep%npwe,factor,deltaeGW_enhigh_b2k,&
1829 33440 : & Ltg_q,Gsph_epsG0,Ep%npwepG0,rhotwg,chi0_sumrule)
1830 33440 : if (band1==Ep%nbnds) then
1831 114240 : chi0_sumrule(:)=chi0_sumrule(:) + wtk_ltg(ik_bz)*spin_fact*qp_occ(band2,ik_ibz,spin)*deltaeGW_enhigh_b2k
1832 : end if
1833 : end if
1834 :
1835 : end do !band2
1836 : end do !band1
1837 :
1838 33707 : if (print_time) then
1839 6122 : write(msg,'(3(a,i0))')' chi0: ik_bz: ',ik_bz,'/',Kmesh%nbz,", spin: ",spin
1840 6122 : call cwtime_report(msg, cpu_k, wall_k, gflops_k); if (ik_bz == LOG_MODK) call wrtout(std_out, " ...")
1841 : end if
1842 :
1843 33707 : ABI_FREE(gw_gbound)
1844 34501 : if (Dtset%pawcross==1) then
1845 0 : ABI_FREE(gboundf)
1846 : end if
1847 : end do !ik_bz
1848 :
1849 : ! Deallocation of arrays private to the spin loop.
1850 794 : ABI_FREE(igfftepsG0)
1851 794 : ABI_FREE(ur1_kmq_ibz)
1852 794 : ABI_FREE(ur2_k_ibz)
1853 794 : ABI_FREE(usr1_kmq)
1854 794 : ABI_FREE(ur2_k)
1855 794 : ABI_FREE(rhotwg)
1856 794 : ABI_FREE(tabr_k)
1857 794 : ABI_FREE(tabr_kmq)
1858 :
1859 794 : ABI_SFREE(green_w)
1860 794 : ABI_SFREE(wfwfg)
1861 794 : ABI_SFREE(green_enhigh_w)
1862 1570 : if (Psps%usepaw==1) then
1863 23 : call pawcprj_free(Cprj2_k)
1864 90 : ABI_FREE(Cprj2_k)
1865 23 : call pawcprj_free(Cprj1_kmq)
1866 90 : ABI_FREE(Cprj1_kmq)
1867 23 : if (Dtset%pawcross==1) then
1868 0 : ABI_FREE(ur_ae1)
1869 0 : ABI_FREE(ur_ae_onsite1)
1870 0 : ABI_FREE(ur_ps_onsite1)
1871 0 : ABI_FREE(ur_ae2)
1872 0 : ABI_FREE(ur_ae_onsite2)
1873 0 : ABI_FREE(ur_ps_onsite2)
1874 0 : ABI_FREE(tabrf_k)
1875 0 : ABI_FREE(tabrf_kmq)
1876 0 : ABI_FREE(gboundf)
1877 0 : ABI_FREE(igfftepsG0f)
1878 : end if
1879 : end if
1880 : end do !spin
1881 :
1882 : ! After big loop over transitions, now MPI
1883 : ! Master took care of the contribution in case of metallic|spin polarized systems.
1884 728 : SELECT CASE (Ep%spmeth)
1885 : CASE (0)
1886 : ! Adler-Wiser
1887 : ! Collective sum of the contributions of each node.
1888 : ! Looping on frequencies to avoid problems with the size of the MPI packet
1889 4846 : do io=1,Ep%nomega
1890 4846 : call xmpi_sum(chi0(:,:,io),comm,ierr)
1891 : end do
1892 :
1893 : CASE (1, 2)
1894 : ! Spectral method.
1895 48 : call hilbert_transform(Ep%npwe,Ep%nomega,Ep%nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,Ep%spmeth)
1896 :
1897 : ! Deallocate here before xmpi_sum
1898 48 : ABI_SFREE(sf_chi0)
1899 :
1900 : ! Collective sum of the contributions.
1901 : ! Looping over frequencies to avoid problems with the size of the MPI packet
1902 727 : do io=1,Ep%nomega
1903 727 : call xmpi_sum(chi0(:,:,io),comm,ierr)
1904 : end do
1905 :
1906 : CASE DEFAULT
1907 776 : ABI_BUG("Wrong spmeth")
1908 : END SELECT
1909 :
1910 : ! Divide by the volume
1911 : !$OMP PARALLEL WORKSHARE
1912 14682915 : chi0 = chi0 * weight / Cryst%ucvol
1913 : !$OMP END PARALLEL WORKSHARE
1914 :
1915 : ! === Collect the sum rule ===
1916 : ! * The pi factor comes from Im[1/(x-ieta)] = pi delta(x)
1917 776 : call xmpi_sum(chi0_sumrule,comm,ierr)
1918 40670 : chi0_sumrule=chi0_sumrule*pi*weight/Cryst%ucvol
1919 : !
1920 : ! *************************************************
1921 : ! **** Now each node has chi0(q,G,Gp,Ep%omega) ****
1922 : ! *************************************************
1923 :
1924 : ! Impose Hermiticity (valid only for zero or purely imaginary frequencies)
1925 : ! MG what about metals, where we have poles around zero?
1926 : ! FB because of the intraband term, chi0 is never hermitian in case of metals
1927 : ! FIXME: as of today, hermitianity is also enforced for metallic systems
1928 : !if (.not. is_metallic) then
1929 5573 : do io=1,Ep%nomega
1930 5573 : if (ABS(REAL(Ep%omega(io))) <0.00001) then
1931 159034 : do ig2=1,Ep%npwe
1932 6418785 : do ig1=1,ig2-1
1933 6415553 : chi0(ig2,ig1,io) = GWPC_CONJG(chi0(ig1,ig2,io))
1934 : end do
1935 : end do
1936 : end if
1937 : end do
1938 : !endif
1939 :
1940 : ! === Symmetrize chi0 in case of AFM system ===
1941 : ! Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$.
1942 : ! Works only in case of magnetic group Shubnikov type IV.
1943 776 : if (Cryst%use_antiferro) then
1944 9 : call symmetrize_afm_chi0(Cryst,Gsph_epsG0,Ltg_q,Ep%npwe,Ep%nomega,chi0=chi0)
1945 : end if
1946 :
1947 : ! =====================
1948 : ! ==== Free memory ====
1949 : ! =====================
1950 776 : ABI_FREE(bbp_ks_distrb)
1951 :
1952 776 : ABI_SFREE(gw_gfft)
1953 776 : ABI_SFREE(kkweight)
1954 776 : ABI_SFREE(omegasf)
1955 776 : ABI_SFREE(gspfft_igfft)
1956 :
1957 776 : call Gsph_FFT%free()
1958 :
1959 : ! deallocation for PAW.
1960 776 : if (Psps%usepaw==1) then
1961 23 : call pawpwij_free(Pwij)
1962 51 : ABI_FREE(Pwij)
1963 23 : if (allocated(Pwij_fft)) then
1964 0 : call pawpwij_free(Pwij_fft)
1965 0 : ABI_FREE(Pwij_fft)
1966 : end if
1967 : end if
1968 :
1969 776 : if(dtset%ucrpa>=1 .AND. dtset%plowan_compute<10) then
1970 0 : ABI_FREE(coeffW_BZ)
1971 : endif
1972 :
1973 776 : call timab(331,2,tsec)
1974 776 : call cwtime_report("cchi0", cpu_all, wall_all, gflops_all)
1975 :
1976 : DBG_EXIT("COLL")
1977 :
1978 1552 : end subroutine cchi0
1979 : !!***
1980 :
1981 :
1982 : !!****f* ABINIT/chi0q0_intraband
1983 : !! NAME
1984 : !! chi0q0_intraband
1985 : !!
1986 : !! FUNCTION
1987 : !! Calculate chi0 in the limit q-->0
1988 : !!
1989 : !! INPUTS
1990 : !! use_tr=If .TRUE. Wfs_val are allocate and only resonant transitions are evaluated (assumes time reversal symmetry)
1991 : !! Ep= datatype gathering differening parameters related to the calculation of the inverse dielectric matrix
1992 : !! Gsph_epsG0<gvectors_data_type>: Info on the G-sphere used to describe chi0/espilon (including umklapp)
1993 : !! %ng=number of G vectors
1994 : !! %rottbm1(ng,2,nsym)=contains the index (IS^{-1}) G in the array gvec
1995 : !! %phmGt(ng,nsym)=phase factor e^{-iG.\tau} needed to symmetrize oscillator matrix elements and chi0
1996 : !! %gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
1997 : !! %gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1)
1998 : !! Ep%nbnds=number of bands
1999 : !! ngfft_gw(18)= array containing all the information for 3D FFT for the oscillator strengths.
2000 : !! Ep%nomega=number of frequencies
2001 : !! Cryst<crystal_t>= data type gathering info on symmetries and unit cell
2002 : !! %natom=number of atoms
2003 : !! %nsym=number of symmetry operations
2004 : !! %symrec(3,3,nsym)=symmetry operations in reciprocal space
2005 : !! %typat(natom)=type of each atom
2006 : !! %xred(3,natom)=reduced coordinated of atoms
2007 : !! %rprimd(3,3)=dimensional primitive translations in real space (bohr)
2008 : !! %timrev=2 if time-reversal symmetry can be used, 1 otherwise
2009 : !! Ep%npwe=number of planewaves for sigma exchange (input variable)
2010 : !! Ep%nsppol=1 for unpolarized, 2 for spin-polarized
2011 : !! Ep%omega(Ep%nomega)=frequencies
2012 : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
2013 : !! %mpsang=1+maximum angular momentum for nonlocal pseudopotential
2014 : !! Pawang<pawang_type> angular mesh discretization and related data:
2015 : !! Pawrad(ntypat*usepaw)<Pawrad_type>=paw radial mesh and related data
2016 : !! Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels
2017 : !! BSt<ebands_t>=Quasiparticle energies and occupations (for the moment real quantities)
2018 : !! %mband=MAX number of bands over k-points and spin (==Ep%nbnds)
2019 : !! %occ(mband,nkpt,nsppol)=QP occupation numbers, for each k point in IBZ, and each band
2020 : !! %eig(mband,nkpt,nsppol)=GW energies, for self-consistency purposes
2021 : !! Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
2022 : !!
2023 : !! OUTPUT
2024 : !! chi0(Ep%npwe,Ep%npwe,Ep%nomega)=independent-particle susceptibility matrix for wavevector qq,
2025 : !! and frequencies defined by Ep%omega
2026 : !!
2027 : !! NOTES
2028 : !! *) The terms "head", "wings" and "body" of chi(G,Gp) refer to
2029 : !! G=Gp=0, either G or Gp=0, and neither=0 respectively
2030 : !!
2031 : !! TODO
2032 : !! Check npwepG0 before Switching on umklapp
2033 : !!
2034 : !! SOURCE
2035 :
2036 0 : subroutine chi0q0_intraband(Wfd,Cryst,Ep,Psps,BSt,Gsph_epsG0,Pawang,Pawrad,Pawtab,Paw_ij,Paw_pwff,use_tr,usepawu,&
2037 0 : & ngfft_gw,chi0,chi0_head,chi0_lwing,chi0_uwing)
2038 :
2039 : !Arguments ------------------------------------
2040 : !scalars
2041 : integer,intent(in) :: usepawu
2042 : logical,intent(in) :: use_tr
2043 : type(ebands_t),intent(in) :: BSt
2044 : type(crystal_t),intent(in) :: Cryst
2045 : type(em1params_t),intent(in) :: Ep
2046 : type(gsphere_t),intent(in) :: Gsph_epsG0
2047 : type(Pseudopotential_type),intent(in) :: Psps
2048 : type(Pawang_type),intent(in) :: Pawang
2049 : type(wfdgw_t),target,intent(inout) :: Wfd
2050 : !arrays
2051 : integer,intent(in) :: ngfft_gw(18)
2052 : complex(gwp),intent(out) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega)
2053 : complex(dp),intent(out) :: chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)
2054 : complex(dp),intent(out) :: chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)
2055 : complex(dp),intent(out) :: chi0_head(3,3,Ep%nomega)
2056 : type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat*Psps%usepaw)
2057 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
2058 : type(Paw_ij_type),intent(in) :: Paw_ij(Cryst%natom*Psps%usepaw)
2059 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
2060 :
2061 : !Local variables ------------------------------
2062 : !scalars
2063 : integer,parameter :: tim_fourdp1=1,two_poles=2,one_pole=1,ndat1=1
2064 : integer,parameter :: unitdos0=0,option1=1,NOMEGA_PRINTED=15
2065 : integer :: nqlwl,nband_k,iomega,istwf_k,npw_k,my_nband,lbidx
2066 : integer :: band,itim_k,ik_bz,ik_ibz,io,isym_k,spin,iqlwl !ig,ig1,ig2,my_nbbp,my_nbbpks
2067 : integer :: nkpt_summed,dim_rtwg,use_padfft,gw_fftalga,ifft
2068 : integer :: kptopt,isym,nsppol,nspinor
2069 : integer :: comm,ierr,gw_mgfft,use_umklp,inclvkb
2070 : real(dp) :: spin_fact,deltaf_b1b2,weight
2071 : real(dp) :: deltaeGW_b1b2,zcut
2072 : real(dp),parameter :: dummy_dosdeltae=HUGE(zero)
2073 : real(dp) :: o_entropy,o_nelect,maxocc
2074 : complex(dp) :: ph_mkt
2075 : logical :: iscompatibleFFT, is_metallic !, print_time
2076 : character(len=500) :: msg,msg_tmp !,allup
2077 0 : type(kmesh_t) :: Kmesh
2078 0 : type(littlegroup_t) :: Ltg_q
2079 0 : type(vkbr_t) :: vkbr
2080 : type(wave_t),pointer :: wave
2081 : !arrays
2082 0 : integer :: my_band_list(Wfd%mband)
2083 0 : integer,contiguous, pointer :: kg_k(:,:)
2084 0 : integer,allocatable :: ktabr(:,:),irottb(:,:)
2085 : !integer :: got(Wfd%nproc)
2086 0 : integer,allocatable :: tabr_k(:),igffteps0(:),gw_gbound(:,:)
2087 : real(dp),parameter :: q0(3)=(/zero,zero,zero/)
2088 : real(dp) :: kpt(3),dedk(3),kbz(3),spinrot_kbz(4)
2089 : !real(dp),contiguous, pointer :: ks_eig(:,:,:),qp_eig(:,:,:),qp_occ(:,:,:)
2090 0 : real(dp) :: shift_ene(BSt%mband,BSt%nkpt,BSt%nsppol)
2091 0 : real(dp) :: delta_occ(BSt%mband,BSt%nkpt,BSt%nsppol)
2092 : !real(dp) :: eigen_vec(BSt%bantot)
2093 0 : real(dp) :: o_doccde(BSt%bantot)
2094 0 : real(dp) :: eigen_pdelta_vec(BSt%bantot),eigen_mdelta_vec(BSt%bantot)
2095 0 : real(dp) :: o_occ_pdelta(BSt%bantot),o_occ_mdelta(BSt%bantot)
2096 0 : real(dp) :: delta_ene(BSt%mband,BSt%nkpt,BSt%nsppol)
2097 0 : real(dp) :: test_docc(BSt%mband,BSt%nkpt,BSt%nsppol)
2098 0 : real(dp),allocatable :: qlwl(:,:)
2099 0 : complex(gwp) :: comm_kbbs(3,Wfd%nspinor**2)
2100 0 : complex(dp),allocatable :: ihr_comm(:,:,:,:,:)
2101 0 : complex(gwp),allocatable :: rhotwg(:)
2102 0 : complex(dp) :: green_w(Ep%nomega)
2103 0 : complex(gwp),allocatable :: ur1(:)
2104 0 : complex(gwp),contiguous, pointer :: ug(:)
2105 0 : logical :: bmask(Wfd%mband)
2106 0 : type(pawcprj_type),allocatable :: Cprj1_bz(:,:),Cprj1_ibz(:,:),Cp_bks(:,:)
2107 0 : type(pawpwij_t),allocatable :: Pwij(:)
2108 0 : type(pawhur_t),allocatable :: Hur(:)
2109 : !************************************************************************
2110 :
2111 : DBG_ENTER("COLL")
2112 :
2113 0 : nsppol = Wfd%nsppol
2114 0 : nspinor = Wfd%nspinor
2115 0 : is_metallic = BSt%has_metal_scheme()
2116 :
2117 0 : gw_mgfft = MAXVAL(ngfft_gw(1:3))
2118 0 : gw_fftalga = ngfft_gw(7)/100 !; gw_fftalgc=MOD(ngfft_gw(7),10)
2119 :
2120 : ! Calculate <k,b1|i[H,r]|k',b2>.
2121 0 : inclvkb=2; if (Wfd%usepaw==1) inclvkb=0
2122 0 : ABI_MALLOC(ihr_comm,(3,nspinor**2,Wfd%mband,Wfd%nkibz,nsppol))
2123 0 : ihr_comm = czero
2124 :
2125 0 : if (Wfd%usepaw==1) then
2126 0 : ABI_MALLOC(Cp_bks,(Cryst%natom,nspinor))
2127 0 : call pawcprj_alloc(Cp_bks,0,Wfd%nlmn_atm)
2128 0 : ABI_MALLOC(HUr,(Cryst%natom))
2129 0 : if (usepawu/=0) then ! For PAW+DFT+U, precalculate <\phi_i|[Hu,r]|phi_j\>.
2130 0 : call pawhur_init(hur,nsppol,Wfd%pawprtvol,Cryst,Pawtab,Pawang,Pawrad,Paw_ij)
2131 : end if
2132 : end if
2133 :
2134 0 : do spin=1,nsppol
2135 0 : do ik_ibz=1,Wfd%nkibz
2136 0 : npw_k = Wfd%npwarr(ik_ibz)
2137 0 : nband_k= Wfd%nband(ik_ibz,spin)
2138 0 : kpt = Wfd%kibz(:,ik_ibz)
2139 0 : kg_k => Wfd%Kdata(ik_ibz)%kg_k
2140 0 : istwf_k = Wfd%istwfk(ik_ibz)
2141 0 : ABI_CHECK(istwf_k==1,"istwf_k/=1 not coded")
2142 :
2143 : ! Distribute bands.
2144 0 : bmask=.FALSE.; bmask(1:nband_k)=.TRUE. ! TODO only bands around EF should be included.
2145 0 : call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,bmask=bmask)
2146 0 : if (my_nband==0) CYCLE
2147 :
2148 0 : if (Wfd%usepaw==0.and.inclvkb/=0) then ! Include term <n,k|[Vnl,iqr]|n"k>' for q->0.
2149 0 : call vkbr_init(vkbr,Cryst,Psps,inclvkb,istwf_k,npw_k,kpt,kg_k)
2150 : end if
2151 :
2152 0 : do lbidx=1,my_nband
2153 0 : band=my_band_list(lbidx)
2154 :
2155 0 : ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg)
2156 0 : ug => wave%ug
2157 :
2158 0 : if (Wfd%usepaw==0) then
2159 : ! Matrix elements of i[H,r] for NC pseudopotentials.
2160 0 : comm_kbbs = nc_ihr_comm(vkbr,cryst,psps,npw_k,nspinor,istwf_k,inclvkb,Kmesh%ibz(:,ik_ibz),ug,ug,kg_k)
2161 : else
2162 : ! Matrix elements of i[H,r] for PAW.
2163 0 : call wfd%get_cprj(band,ik_ibz,spin,Cryst,Cp_bks,sorted=.FALSE.)
2164 0 : comm_kbbs = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug,ug,kg_k,Cp_bks,Cp_bks,HUr)
2165 : end if
2166 :
2167 0 : ihr_comm(:,:,band,ik_ibz,spin) = comm_kbbs
2168 : end do
2169 :
2170 0 : call vkbr_free(vkbr) ! Not need anymore as we loop only over IBZ.
2171 : end do
2172 : end do
2173 : !
2174 : ! Gather the commutator on each node.
2175 0 : call xmpi_sum(ihr_comm,Wfd%comm,ierr)
2176 :
2177 0 : if (Wfd%usepaw==1) then
2178 0 : call pawcprj_free(Cp_bks)
2179 0 : ABI_FREE(Cp_bks)
2180 0 : call pawhur_free(Hur)
2181 0 : ABI_FREE(Hur)
2182 : end if
2183 :
2184 0 : nqlwl=1
2185 0 : ABI_MALLOC(qlwl,(3,nqlwl))
2186 : !qlwl = GW_Q0_DEFAULT(3)
2187 0 : qlwl(:,1) = (/0.00001_dp, 0.00002_dp, 0.00003_dp/)
2188 : !
2189 0 : write(msg,'(a,i3,a)')' Q-points for long wave-length limit in chi0q_intraband. # ',nqlwl,ch10
2190 0 : do iqlwl=1,nqlwl
2191 0 : write(msg_tmp,'(1x,i5,a,2x,3f12.6,a)') iqlwl,')',qlwl(:,iqlwl),ch10
2192 0 : msg=TRIM(msg)//msg_tmp
2193 : end do
2194 0 : call wrtout(std_out, msg)
2195 : !
2196 : ! delta_ene = e_{b,k-q} - e_{b,k} = -q. <b,k| i[H,r] |b,k> + O(q^2).
2197 0 : delta_ene = zero
2198 0 : do spin=1,nsppol
2199 0 : do ik_ibz=1,Wfd%nkibz
2200 0 : do band=1,Wfd%nband(ik_ibz,spin)
2201 0 : dedk = REAL(ihr_comm(:,1,band,ik_ibz,spin))
2202 0 : delta_ene(band,ik_ibz,spin) = -vdotw(qlwl(:,1),dedk,Cryst%gmet,"G")
2203 : end do
2204 : end do
2205 : end do
2206 :
2207 0 : maxocc=two/(nsppol*nspinor)
2208 :
2209 : ! Calculate the occupations at f(e+delta/2).
2210 0 : shift_ene = BSt%eig + half*delta_ene
2211 :
2212 0 : call pack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,BSt%bantot,shift_ene,eigen_pdelta_vec)
2213 :
2214 0 : if (BSt%occopt < 9) then
2215 : call getnel(o_doccde,dummy_dosdeltae,eigen_pdelta_vec,o_entropy,BSt%fermie,BSt%fermie,maxocc,BSt%mband,BSt%nband,&
2216 0 : & o_nelect,BSt%nkpt,BSt%nsppol,o_occ_pdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk,1,BSt%nband(1))
2217 : else
2218 0 : ABI_ERROR('occopt 9 not implemented for GW calculations')
2219 : end if
2220 : !
2221 : ! Calculate the occupations at f(e-delta/2).
2222 0 : shift_ene = BSt%eig - half*delta_ene
2223 :
2224 0 : call pack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,BSt%bantot,shift_ene,eigen_mdelta_vec)
2225 :
2226 0 : if (BSt%occopt < 9) then
2227 : call getnel(o_doccde,dummy_dosdeltae,eigen_mdelta_vec,o_entropy,BSt%fermie,BSt%fermie,maxocc,BSt%mband,BSt%nband,&
2228 0 : & o_nelect,BSt%nkpt,BSt%nsppol,o_occ_mdelta,BSt%occopt,option1,BSt%tphysel,BSt%tsmear,unitdos0,BSt%wtk,1,BSt%nband(1))
2229 0 : write(std_out,*)"nelect2: ",o_nelect
2230 : else
2231 0 : ABI_ERROR("occopt 9 not implemented for GW calculations")
2232 : end if
2233 :
2234 : ! f(e-delta/2) - f(e+delta/2).
2235 0 : o_occ_pdelta = o_occ_mdelta - o_occ_pdelta
2236 :
2237 0 : call unpack_eneocc(BSt%nkpt,BSt%nsppol,BSt%mband,BSt%nband,o_occ_pdelta,delta_occ)
2238 : !
2239 : ! Expand f(e-delta/2) - f(e+delta/2) up to the first order in the small q.
2240 0 : do spin=1,nsppol
2241 0 : do ik_ibz=1,Wfd%nkibz
2242 0 : do band=1,Wfd%nband(ik_ibz,spin)
2243 0 : dedk = REAL(ihr_comm(:,1,band,ik_ibz,spin))
2244 0 : test_docc(band,ik_ibz,spin) = +vdotw(qlwl(:,1),dedk,Cryst%gmet,"G") * BSt%doccde(band,ik_ibz,spin)
2245 0 : write(std_out,'(a,3(i0,1x),1x,3es16.8)')" spin,ik_ibz,band, delta_occ: ",&
2246 0 : & spin,ik_ibz,band,delta_occ(band,ik_ibz,spin),&
2247 0 : & test_docc(band,ik_ibz,spin),delta_occ(band,ik_ibz,spin)-test_docc(band,ik_ibz,spin)
2248 : end do
2249 : end do
2250 : end do
2251 :
2252 : ! ABI_ERROR("DONE")
2253 : ! do spin=1,nsppol
2254 : ! do ik_ibz=1,Wfd%nkibz
2255 : ! nband_k = Wfd%nband(ik_ibz,spin)
2256 : ! do band=1,nband_k
2257 : ! write(std_out,'(a,3i3,2es14.6)')" spin, band, ik_ibz, delta_ene, delta_occ ",&
2258 : !& spin,band,ik_ibz,delta_ene(band,ik_ibz,spin),delta_occ(band,ik_ibz,spin)
2259 : ! end do
2260 : ! end do
2261 : ! end do
2262 :
2263 0 : ABI_FREE(ihr_comm)
2264 0 : ABI_FREE(qlwl)
2265 :
2266 0 : if ( ANY(ngfft_gw(1:3) /= Wfd%ngfft(1:3)) ) call wfd%change_ngfft(Cryst,Psps,ngfft_gw)
2267 :
2268 : ! TODO take into account the case of random k-meshes.
2269 0 : kptopt=3
2270 0 : call Kmesh%init(Cryst,Wfd%nkibz,Wfd%kibz,kptopt)
2271 : !
2272 : !=== Get the FFT index of $ (R^{-1}(r-\tau)) $ ===
2273 : !* S= $\transpose R^{-1}$ and k_BZ = S k_IBZ
2274 : !* irottb is the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Sk.
2275 0 : ABI_MALLOC(irottb,(Wfd%nfftot,Cryst%nsym))
2276 :
2277 0 : call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,irottb,iscompatibleFFT)
2278 0 : ABI_CHECK(iscompatibleFFT,"FFT mesh not compatible with symmetries")
2279 :
2280 0 : ABI_MALLOC(ktabr,(Wfd%nfftot,Kmesh%nbz))
2281 0 : do ik_bz=1,Kmesh%nbz
2282 0 : isym=Kmesh%tabo(ik_bz)
2283 0 : do ifft=1,Wfd%nfftot
2284 0 : ktabr(ifft,ik_bz)=irottb(ifft,isym)
2285 : end do
2286 : end do
2287 0 : ABI_FREE(irottb)
2288 : !
2289 : ! === Setup weight (2 for spin unpolarized systems, 1 for polarized) ===
2290 : ! spin_fact is used to normalize the occupation factors to one.
2291 : ! Consider also the AFM case.
2292 0 : SELECT CASE (nsppol)
2293 : CASE (1)
2294 0 : weight=two/Kmesh%nbz; spin_fact=half
2295 0 : if (Wfd%nspden==2) then
2296 0 : weight=one/Kmesh%nbz; spin_fact=half
2297 : end if
2298 0 : if (nspinor==2) then
2299 0 : weight=one/Kmesh%nbz; spin_fact=one
2300 : end if
2301 : CASE (2)
2302 0 : weight=one/Kmesh%nbz; spin_fact=one
2303 : CASE DEFAULT
2304 0 : ABI_BUG("Wrong nsppol")
2305 : END SELECT
2306 :
2307 0 : use_umklp = 0
2308 0 : call Ltg_q%init(q0, Kmesh%nbz, Kmesh%bz, Cryst, use_umklp, Ep%npwepG0, gvec=Gsph_epsG0%gvec)
2309 :
2310 0 : write(msg,'(a,i2)')' Using symmetries to sum only over the IBZ_q = ',Ep%symchi
2311 0 : call wrtout(std_out, msg)
2312 : !
2313 : ! Evaluate oscillator matrix elements btw partial waves. Note that q=Gamma is used.
2314 0 : if (Psps%usepaw==1) then
2315 0 : ABI_MALLOC(Pwij,(Psps%ntypat))
2316 0 : call pawpwij_init(Pwij,Ep%npwepG0, [zero,zero,zero], Gsph_epsG0%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
2317 :
2318 0 : ABI_MALLOC(Cprj1_bz ,(Cryst%natom,nspinor))
2319 0 : call pawcprj_alloc(Cprj1_bz, 0,Wfd%nlmn_atm)
2320 0 : ABI_MALLOC(Cprj1_ibz,(Cryst%natom,nspinor))
2321 0 : call pawcprj_alloc(Cprj1_ibz,0,Wfd%nlmn_atm)
2322 : end if
2323 :
2324 0 : ABI_MALLOC(rhotwg,(Ep%npwe*nspinor**2))
2325 0 : ABI_MALLOC(tabr_k,(Wfd%nfftot))
2326 0 : ABI_MALLOC(ur1,(Wfd%nfft*nspinor))
2327 : !
2328 : ! Tables for the FFT of the oscillators.
2329 : ! a) FFT index of the G sphere (only vertical transitions, unlike cchi0, no need to shift the sphere).
2330 : ! b) gw_gbound table for the zero-padded FFT performed in rhotwg.
2331 0 : ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2))
2332 0 : ABI_MALLOC(igffteps0,(Gsph_epsG0%ng))
2333 :
2334 0 : call Gsph_epsG0%fft_tabs([0, 0, 0], gw_mgfft,ngfft_gw,use_padfft,gw_gbound,igffteps0)
2335 0 : if ( ANY(gw_fftalga == [2, 4]) ) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
2336 0 : if (use_padfft==0) then
2337 0 : ABI_FREE(gw_gbound)
2338 0 : ABI_MALLOC(gw_gbound,(2*gw_mgfft+8,2*use_padfft))
2339 : end if
2340 :
2341 0 : nkpt_summed=Kmesh%nbz
2342 0 : if (Ep%symchi/=0) then
2343 0 : nkpt_summed=Ltg_q%nibz_ltg
2344 0 : call Ltg_q%print([std_out], prtvol=Wfd%prtvol)
2345 : end if
2346 : !
2347 : ! ============================================
2348 : ! === Begin big fat loop over transitions ====
2349 : ! ============================================
2350 0 : chi0 = czero_gw
2351 0 : chi0_head = czero_gw; chi0_lwing = czero_gw; chi0_uwing = czero_gw
2352 0 : dim_rtwg=1; if (nspinor==2) dim_rtwg=2 !can reduce size depending on Ep%nI and Ep%nj
2353 :
2354 : zcut = Ep%zcut
2355 0 : zcut = 0.1/Ha_eV
2356 0 : write(std_out,*)" using zcut ",zcut*Ha_eV," [eV]"
2357 :
2358 : ! Loop on spin to calculate $ \chi_{\up,\up} + \chi_{\down,\down}
2359 0 : do spin=1,nsppol
2360 : ! Loop over k-points in the BZ.
2361 0 : do ik_bz=1,Kmesh%nbz
2362 0 : if (Ep%symchi==1) then
2363 0 : if (Ltg_q%ibzq(ik_bz)/=1) CYCLE ! Only IBZ_q
2364 : end if
2365 :
2366 : ! Get ik_ibz, non-symmorphic phase and symmetries from ik_bz.
2367 0 : call kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt)
2368 0 : tabr_k=ktabr(:,ik_bz) ! Table for rotated FFT points
2369 0 : spinrot_kbz(:)=Cryst%spinrot(:,isym_k)
2370 0 : nband_k=Wfd%nband(ik_ibz,spin)
2371 :
2372 : ! Distribute bands.
2373 0 : bmask=.FALSE.; bmask(1:nband_k)=.TRUE. ! TODO only bands around EF should be included.
2374 0 : call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,bmask=bmask)
2375 0 : if (my_nband==0) CYCLE
2376 :
2377 0 : write(msg,'(2(a,i0),a,i0,a,i0)')' ik: ',ik_bz,' / ',Kmesh%nbz,' spin: ',spin,' done by rank: ',Wfd%my_rank
2378 0 : call wrtout(std_out, msg)
2379 :
2380 0 : do lbidx=1,my_nband
2381 : ! Loop over bands treated by this node.
2382 0 : band=my_band_list(lbidx)
2383 0 : call wfd%get_ur(band,ik_ibz,spin,ur1)
2384 :
2385 0 : if (Psps%usepaw==1) then
2386 0 : call wfd%get_cprj(band,ik_ibz,spin,Cryst,Cprj1_ibz,sorted=.FALSE.)
2387 0 : call pawcprj_copy(Cprj1_ibz,Cprj1_bz)
2388 0 : call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj1_bz)
2389 : end if
2390 :
2391 0 : deltaf_b1b2 = spin_fact*delta_occ(band,ik_ibz,spin)
2392 0 : deltaeGW_b1b2= delta_ene(band,ik_ibz,spin)
2393 :
2394 : ! Add small imaginary of the Time-Ordered resp function but only for non-zero real omega FIXME What about metals?
2395 0 : if (.not.use_tr) then
2396 0 : do io=1,Ep%nomega
2397 : !green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,-one,one_pole)
2398 0 : green_w(io) = g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,GW_TOL_W0,one_pole)
2399 : end do
2400 : else
2401 0 : do io=1,Ep%nomega ! This expression implements time-reversal even when the input k-mesh breaks it.
2402 : !green_w(io) = half * g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,-one,two_poles)
2403 0 : green_w(io) = half * g0g0w(Ep%omega(io),deltaf_b1b2,deltaeGW_b1b2,zcut,GW_TOL_W0,two_poles)
2404 : end do !io
2405 : end if ! use_tr
2406 :
2407 : ! FFT of u^*_{b1,k}(r) u_{b2,k}(r).
2408 : call rho_tw_g(nspinor,Ep%npwe,Wfd%nfft,ndat1,ngfft_gw,1,use_padfft,igffteps0,gw_gbound,&
2409 : & ur1,itim_k,tabr_k,ph_mkt,spinrot_kbz,&
2410 : & ur1,itim_k,tabr_k,ph_mkt,spinrot_kbz,&
2411 0 : & dim_rtwg,rhotwg)
2412 :
2413 0 : if (Psps%usepaw==1) then
2414 : ! Add PAW onsite contribution, projectors are already in the BZ.
2415 0 : call paw_rho_tw_g(cryst,Pwij,Ep%npwe,dim_rtwg,nspinor,Gsph_epsG0%gvec,Cprj1_bz,Cprj1_bz,rhotwg)
2416 : end if
2417 :
2418 : ! ==== Adler-Wiser expression, to be consistent here we use the KS eigenvalues (?) ====
2419 0 : call assemblychi0_sym(is_metallic,ik_bz,nspinor,Ep,Ltg_q,green_w,Ep%npwepG0,rhotwg,Gsph_epsG0,chi0)
2420 : end do !band
2421 : end do !ik_bz
2422 : end do !spin
2423 :
2424 : ! Collect body, head and wings within comm
2425 0 : comm=Wfd%comm
2426 0 : do io=1,Ep%nomega
2427 0 : call xmpi_sum(chi0(:,:,io),comm,ierr)
2428 : end do
2429 0 : call xmpi_sum(chi0_head,comm,ierr)
2430 0 : call xmpi_sum(chi0_lwing,comm,ierr)
2431 0 : call xmpi_sum(chi0_uwing,comm,ierr)
2432 :
2433 : ! Divide by the volume
2434 0 : chi0 = chi0 * weight/Cryst%ucvol
2435 0 : chi0_head = chi0_head * weight/Cryst%ucvol
2436 0 : do io=1,Ep%nomega ! Tensor in the basis of the reciprocal lattice vectors.
2437 0 : chi0_head(:,:,io) = MATMUL(chi0_head(:,:,io),Cryst%gmet) * (two_pi**2)
2438 : end do
2439 0 : chi0_lwing = chi0_lwing * weight/Cryst%ucvol
2440 0 : chi0_uwing = chi0_uwing * weight/Cryst%ucvol
2441 : !
2442 : ! ===============================================
2443 : ! ==== Symmetrize chi0 in case of AFM system ====
2444 : ! ===============================================
2445 : ! * Reconstruct $chi0{\down,\down}$ from $chi0{\up,\up}$.
2446 : ! * Works only in the case of magnetic group Shubnikov type IV.
2447 0 : if (Cryst%use_antiferro) then
2448 : call symmetrize_afm_chi0(Cryst, Gsph_epsG0, Ltg_q, Ep%npwe, Ep%nomega, chi0=chi0, &
2449 0 : chi0_head=chi0_head, chi0_lwing=chi0_lwing, chi0_uwing=chi0_uwing)
2450 : end if
2451 : !
2452 : ! ==================================================
2453 : ! ==== Construct head and wings from the tensor ====
2454 : ! ==================================================
2455 : !do io=1,Ep%nomega
2456 : ! do ig=2,Ep%npwe
2457 : ! wng = chi0_uwing(ig,io,:)
2458 : ! chi0(1,ig,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G")
2459 : ! wng = chi0_lwing(ig,io,:)
2460 : ! chi0(ig,1,io) = vdotw(Ep%qlwl(:,1),wng,Cryst%gmet,"G")
2461 : ! end do
2462 : ! chq = MATMUL(chi0_head(:,:,io), Ep%qlwl(:,1))
2463 : ! chi0(1,1,io) = vdotw(Ep%qlwl(:,1),chq,Cryst%gmet,"G") ! Use user-defined small q
2464 : !end do
2465 : !call wfd_barrier(Wfd)
2466 :
2467 : ! Impose Hermiticity (valid only for zero or purely imaginary frequencies)
2468 : ! MG what about metals, where we have poles around zero?
2469 : !do io=1,Ep%nomega
2470 : ! if (ABS(REAL(Ep%omega(io)))<0.00001) then
2471 : ! do ig2=1,Ep%npwe
2472 : ! do ig1=1,ig2-1
2473 : ! chi0(ig2,ig1,io)=CONJG(chi0(ig1,ig2,io))
2474 : ! end do
2475 : ! end do
2476 : ! end if
2477 : !end do
2478 :
2479 0 : do iomega=1,MIN(Ep%nomega,NOMEGA_PRINTED)
2480 0 : write(msg,'(1x,a,i4,a,2f9.4,a)')' chi0_intra(G,G'') at the ',iomega,' th omega',Ep%omega(iomega)*Ha_eV,' [eV]'
2481 0 : call wrtout(std_out, msg)
2482 0 : call print_arr([std_out], chi0(:,:,iomega))
2483 : end do
2484 :
2485 : ! =====================
2486 : ! ==== Free memory ====
2487 : ! =====================
2488 0 : ABI_FREE(rhotwg)
2489 0 : ABI_FREE(tabr_k)
2490 0 : ABI_FREE(ur1)
2491 0 : ABI_FREE(gw_gbound)
2492 0 : ABI_FREE(ktabr)
2493 0 : ABI_FREE(igffteps0)
2494 :
2495 : ! deallocation for PAW.
2496 0 : if (Psps%usepaw==1) then
2497 0 : call pawcprj_free(Cprj1_bz)
2498 0 : ABI_FREE(Cprj1_bz)
2499 0 : call pawcprj_free(Cprj1_ibz)
2500 0 : ABI_FREE(Cprj1_ibz)
2501 0 : call pawpwij_free(Pwij)
2502 0 : ABI_FREE(Pwij)
2503 : end if
2504 :
2505 0 : call Ltg_q%free()
2506 0 : call Kmesh%free()
2507 :
2508 : DBG_EXIT("COLL")
2509 :
2510 0 : end subroutine chi0q0_intraband
2511 : !!***
2512 :
2513 910 : end module m_chi0
2514 : !!***
|
