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