Line data Source code
1 : !!****m* ABINIT/m_sigc
2 : !! NAME
3 : !! m_sigc
4 : !!
5 : !! FUNCTION
6 : !! Compute matrix elements of the correlated part of the e-h self-energy
7 : !!
8 : !! SOURCE
9 :
10 : #if defined HAVE_CONFIG_H
11 : #include "config.h"
12 : #endif
13 :
14 : #include "abi_common.h"
15 :
16 : module m_sigc
17 :
18 : use, intrinsic :: iso_c_binding
19 : use defs_basis
20 : use m_gwdefs
21 : use m_abicore
22 : use m_xmpi
23 : use m_xomp
24 : use m_defs_ptgroups
25 : use m_errors
26 : use m_splines
27 : use m_dtset
28 :
29 : use m_io_tools, only : iomode_from_fname
30 : use defs_datatypes, only : pseudopotential_type
31 : use m_dtfil, only : datafiles_type
32 : use m_array, only : array2_gwpc_t
33 : use m_time, only : timab, cwtime, cwtime_report
34 : use m_hide_blas, only : xdotc, xgemv, xgemm, xherk
35 : use m_hide_lapack, only : xheev
36 : use m_numeric_tools, only : hermitianize, imin_loc, coeffs_gausslegint
37 : use m_fstrings, only : sjoin, itoa
38 : use m_geometry, only : normv
39 : use m_crystal, only : crystal_t
40 : use m_bz_mesh, only : kmesh_t, findqg0, littlegroup_t
41 : use m_gsphere, only : gsphere_t
42 : use m_fft_mesh, only : get_gfft, rotate_fft_mesh, cigfft
43 : use m_vcoul, only : vcoul_t
44 : use m_wfd, only : wfdgw_t, wave_t
45 : use m_oscillators, only : rho_tw_g, calc_wfwfg
46 : use m_screening, only : epsm1_t
47 : use m_ppmodel, only : ppmodel_t
48 : use m_screen, only : em1_symmetrize_op
49 : use m_sigma, only : sigma_t, sigma_distribute_bks
50 : use m_esymm, only : esymm_t
51 : use m_pawang, only : pawang_type
52 : use m_pawtab, only : pawtab_type
53 : use m_pawfgrtab, only : pawfgrtab_type
54 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
55 : use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
56 : use m_paw_sym, only : paw_symcprj
57 : use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t
58 : use m_hide_lapack, only : xheev
59 : use m_occ, only : get_fact_spin_tol_empty
60 : use m_ebands, only : ebands_t
61 : use m_pstat, only : pstat_proc
62 : use m_io_screening, only : read_screening ! hscr_t, get_hscr_qmesh_gsph,
63 :
64 : implicit none
65 :
66 : private
67 : !!***
68 :
69 : public :: calc_sigc_me
70 : !!***
71 :
72 : integer,parameter :: LOG_MODK = 5
73 :
74 : contains
75 : !!***
76 :
77 : !!****f* ABINIT/calc_sigc_me
78 : !! NAME
79 : !! calc_sigc_me
80 : !!
81 : !! FUNCTION
82 : !! Calculate diagonal and off-diagonal matrix elements of the self-energy operator.
83 : !!
84 : !! INPUTS
85 : !! sigmak_ibz=Index of the k-point in the IBZ.
86 : !! minbnd, maxbnd= min and Max band index for GW correction (for this k-point)
87 : !! Dtset <type(dataset_type)>=all input variables in this dataset
88 : !! epsm1 <epsm1_t> (see the definition of this structured datatype)
89 : !! Gsph_c<gsphere_t>= info on G-sphere for Sigma_c
90 : !! Gsph_Max<gsphere_t>= info on biggest G-sphere
91 : !! ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
92 : !! Ltg_k datatype containing information on the little group
93 : !! Kmesh <kmesh_t>
94 : !! gwc_ngfft(18)=Information about 3D FFT for the oscillator strengths used for the correlation part,
95 : !! Vcp <vcoul_t datatype> containing information on the cutoff technique
96 : !! Pawtab(Psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
97 : !! Pawang <type(pawang_type)>=paw angular mesh and related data
98 : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
99 : !! Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
100 : !! Sigp <sigparams_t> (see the definition of this structured datatype)
101 : !! Cryst<crystal_t>=Info on unit cell and symmetries
102 : !! PPm<ppmodel_t>= Datatype gathering information on the Plasmonpole technique (see also ppm_get_qbz).
103 : !! QP_BSt<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
104 : !! eig(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=KS or QP energies for k-points, bands and spin
105 : !! occ(Sigp%nbnds,Kmesh%nibz,Wfd%nsppol)=occupation numbers, for each k point in IBZ, each band and spin
106 : !! Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
107 : !! allQP_sym(Wfd%nkibz,Wfd%nsppol)<esymm_t>=Datatype collecting data on the irreducible representations of the
108 : !! little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
109 : !! Sr=sigma_t (see the definition of this structured datatype)
110 : !! use_aerhor=1 is aepaw_rhor is used, 0 otherwise.
111 : !! aepaw_rhor(rho_nfftot,Wfd%nspden*use_aerhor)=AE PAW density used to generate PPmodel parameters if mqmem==0
112 : !!
113 : !! OUTPUT
114 : !!
115 : !! NOTES
116 : !! 1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) [[cite:Gigy1986]] is included.
117 : !! 2) The calculation of energy derivative is based on finite elements.
118 : !! 3) On the symmetrization of Sigma matrix elements ***/
119 : !! If Sk = k+G0 then M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G} (k,q)
120 : !! If -Sk = k+G0 then M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)
121 : !!
122 : !! Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out
123 : !! that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a
124 : !! similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase
125 : !! has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ.
126 : !!
127 : !! 4) The unitary transformation relating wavefunctions
128 : !! at symmetric k-points should be taken into account during the symmetrization
129 : !! of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
130 : !! it is possible to make an invariant by just including all the degenerate states and
131 : !! averaging the final results over the degenerate subset.
132 : !!
133 : !! SOURCE
134 :
135 455 : subroutine calc_sigc_me(sigmak_ibz,ikcalc,nomega_sigc,minbnd,maxbnd,&
136 : Dtset, dtfil, Cryst,QP_BSt,Sigp,Sr,epsm1,Gsph_Max,Gsph_c,Vcp,Kmesh,Qmesh,Ltg_k,&
137 455 : PPm,Pawtab,Pawang,Paw_pwff,Pawfgrtab,Paw_onsite,Psps,Wfd,Wfdf,allQP_sym,&
138 455 : gwc_ngfft,rho_ngfft,rho_nfftot,rhor,use_aerhor,aepaw_rhor,sigcme_tmp)
139 :
140 : !Arguments ------------------------------------
141 : !scalars
142 : integer,intent(in) :: sigmak_ibz,ikcalc,rho_nfftot,nomega_sigc,minbnd,maxbnd
143 : integer,intent(in) :: use_aerhor
144 : type(crystal_t),intent(in) :: Cryst
145 : type(ebands_t),target,intent(in) :: QP_BSt
146 : type(kmesh_t),intent(in) :: Kmesh,Qmesh
147 : type(vcoul_t),intent(in) :: Vcp
148 : type(dataset_type),intent(in) :: Dtset
149 : type(datafiles_type),intent(in) :: dtfil
150 : type(epsm1_t),intent(inout) :: epsm1
151 : type(gsphere_t),intent(in) :: Gsph_Max,Gsph_c
152 : type(littlegroup_t),intent(in) :: Ltg_k
153 : type(ppmodel_t),intent(inout) :: PPm
154 : type(Pseudopotential_type),intent(in) :: Psps
155 : type(pawang_type),intent(in) :: pawang
156 : type(sigparams_t),target,intent(in) :: Sigp
157 : type(sigma_t),intent(in) :: Sr
158 : type(wfdgw_t),target,intent(inout) :: Wfd,Wfdf
159 : !arrays
160 : integer,intent(in) :: gwc_ngfft(18),rho_ngfft(18)
161 : real(dp),intent(in) :: rhor(rho_nfftot,Wfd%nspden)
162 : real(dp),intent(in) :: aepaw_rhor(rho_nfftot,Wfd%nspden*use_aerhor)
163 : complex(dp),intent(out) :: sigcme_tmp(nomega_sigc,minbnd:maxbnd,minbnd:maxbnd,Wfd%nsppol*Sigp%nsig_ab)
164 : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
165 : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat*Psps%usepaw)
166 : type(esymm_t),target,intent(in) :: allQP_sym(Wfd%nkibz,Wfd%nsppol)
167 : type(pawfgrtab_type),intent(inout) :: Pawfgrtab(Cryst%natom*Psps%usepaw)
168 : type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom*Psps%usepaw)
169 :
170 : !Local variables ------------------------------
171 : !scalars
172 : integer,parameter :: spinor_idxs(2,4)=RESHAPE([1,1,2,2,1,2,2,1], [2,4])
173 : integer,parameter :: tim_fourdp2=2,ndat1=1
174 : integer :: npw_k,iab,ib,ib1,ib2,ierr,ig,ii,iik,itim_q,i1,i2,npls,ib_sum
175 : integer :: ik_bz,ik_ibz,io,iiw,isym_q,iq_bz,iq_ibz,spin,isym,jb,is_idx,iiab,jiab
176 : integer :: band,band1,band2,idle,rank,jik,jk_bz,jk_ibz,kb,nspinor
177 : integer :: nomega_tot,nq_summed,ibsp,dimcprj_gw,npwc, scr_iomode
178 : integer :: spad,spadc1,spadc2,irow,my_nbks,ndegs,wtqm,wtqp,mod10, iwc,ifft
179 : integer :: isym_kgw,isym_ki,gwc_mgfft,use_padfft,gwc_fftalga,gwc_nfftot,nfftf,mgfftf,use_padfftf
180 : integer :: ilwrk, neigmax, ac_epsm1cqwz2_win
181 : integer(kind=XMPI_ADDRESS_KIND) :: ad_count
182 : type(c_ptr) :: void_ptr
183 : real(dp) :: cpu_all, wall_all, gflops_all, cpu_k, wall_k, gflops_k
184 : real(dp) :: e0i,fact_spin,theta_mu_minus_e0i,tol_empty,tol_empty_in,en_high,gw_gsq,w_localmax,w_max
185 : complex(dp) :: ctmp,omegame0i2_ac,omegame0i_ac,ph_mkgwt,ph_mkt
186 : logical :: iscompatibleFFT, q_is_gamma, print_time
187 : character(len=500) :: msg
188 : type(wave_t),pointer :: wave_sum, wave_jb
189 455 : complex(gwp),allocatable :: botsq(:,:),otq(:,:),eig(:,:)
190 : !arrays
191 910 : integer :: g0(3),spinor_padc(2,4),got(Wfd%nproc), units(2)
192 455 : integer,allocatable :: proc_distrb(:,:,:),extrapolar_distrb(:,:,:,:),degtab(:,:,:)
193 455 : integer,allocatable :: igfftcg0(:),gw_gfft(:,:),gw_gbound(:,:),irottb(:,:),ktabr(:,:)
194 455 : integer,allocatable :: igfftfcg0(:),gboundf(:,:),ktabrf(:,:),npoles_missing(:)
195 : real(dp) :: ksum(3),kgw(3),kgw_m_ksum(3),q0(3),tsec(2),qbz(3)
196 : real(dp) :: spinrot_kbz(4),spinrot_kgw(4)
197 1820 : real(dp) :: gl_knots(epsm1%nomega_i),gl_wts(epsm1%nomega_i), omegap(epsm1%nomega_i), omegap2(epsm1%nomega_i), tmp_omegap(epsm1%nomega_i)
198 1365 : real(dp) :: conv_gl_knots(epsm1%nomega_i_conv),conv_gl_wts(epsm1%nomega_i_conv),conv_omegap(epsm1%nomega_i_conv),conv_omegap2(epsm1%nomega_i_conv), tmp_conv_omegap(epsm1%nomega_i_conv)
199 455 : real(dp),contiguous, pointer :: qp_ene(:,:,:),qp_occ(:,:,:)
200 455 : real(dp),allocatable :: omegame0i(:), w_maxval(:)
201 1365 : complex(gwp) :: sigcohme(Sigp%nsig_ab), omegap_cplx(epsm1%nomega_i)
202 455 : complex(gwp),allocatable :: vc_sqrt_qbz(:),rhotwg(:),rhotwgp(:), coh_g(:)
203 455 : complex(gwp),allocatable :: botsq_conjg_transp(:,:)
204 : complex(gwp),pointer, contiguous :: ac_epsm1cqwz2(:,:,:) => null()
205 455 : complex(gwp),allocatable :: epsm1_trcc_qbz(:,:,:), epsm1_tmp(:,:)
206 455 : complex(gwp),allocatable :: sigc_ket(:,:),ket1(:,:),ket2(:,:)
207 455 : complex(gwp),allocatable :: herm_sigc_ket(:,:),aherm_sigc_ket(:,:), rhotwg_ki(:,:)
208 455 : complex(gwp),allocatable :: sigcme2(:,:),sigcme_3(:),sigcme_new(:),sigctmp(:,:)
209 455 : complex(gwp),allocatable :: wfr_bdgw(:,:),ur_ibz(:),wf1swf2_g(:),usr_bz(:)
210 455 : complex(gwp),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:)
211 455 : complex(gwp),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:), otq_transp(:,:), epsm1_ggw(:,:,:)
212 455 : complex(gwp),contiguous, pointer :: cg_jb(:),cg_sum(:)
213 455 : complex(dp),allocatable :: sym_cme(:,:,:,:),sigc(:,:,:,:,:)
214 1365 : logical :: rank_mask(Wfd%nproc),can_symmetrize(Wfd%nsppol)
215 : !logical :: me_calc_poles(Sr%nomega_r+Sr%nomega4sd)
216 455 : type(sigijtab_t),pointer :: Sigcij_tab(:)
217 455 : type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
218 455 : type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
219 455 : type(esymm_t),pointer :: QP_sym(:)
220 910 : integer :: neig(epsm1%nomega_i)
221 455 : real(gwp),allocatable :: epsm1_eig(:)
222 455 : complex(gwp),allocatable :: epsm1_sqrt_rhotw(:,:), rhotw_eqsm1_sqrt(:,:), rhotw_epsm1_rhotw(:,:,:,:), conv_rhotw_epsm1_rhotw(:,:,:,:)
223 910 : complex(dp) :: tmp_rhotw_epsm1_rhotw(epsm1%nomega_i), tmp_conv_rhotw_epsm1_rhotw(epsm1%nomega_i_conv)
224 455 : type(array2_gwpc_t),allocatable :: wc0_qibz(:)
225 455 : complex(gwp),allocatable :: wc0_qbz(:,:) !, vec_coh_nk(:,:), vec_coh_mkq(:,:),
226 : !************************************************************************
227 :
228 : DBG_ENTER("COLL")
229 :
230 1365 : units = [std_out, ab_out]
231 :
232 : ! Initial check
233 455 : ABI_CHECK_IEQ(Sr%nomega_r, Sigp%nomegasr, "")
234 455 : ABI_CHECK_IEQ(Sr%nomega4sd, Sigp%nomegasrd, "")
235 455 : ABI_CHECK_IEQ(Sigp%npwc, Gsph_c%ng, "")
236 455 : ABI_CHECK_IEQ(Sigp%npwvec, Gsph_Max%ng, "")
237 :
238 455 : mod10 = MOD(Sigp%gwcalctyp,10)
239 :
240 455 : call timab(424,1,tsec) ! calc_sigc_me
241 455 : call timab(431,1,tsec) ! calc_sigc_me
242 455 : call timab(432,1,tsec) ! Init
243 455 : call cwtime(cpu_all, wall_all, gflops_all,"start")
244 :
245 455 : qp_ene => QP_BSt%eig; qp_occ => QP_BSt%occ
246 :
247 : ! Extract the symmetries of the bands for this k-point
248 455 : QP_sym => allQP_sym(sigmak_ibz,1:Wfd%nsppol)
249 :
250 : ! Index of the GW point in the BZ array, its image in IBZ and time-reversal
251 455 : jk_bz=Sigp%kptgw2bz(ikcalc)
252 455 : call kmesh%get_BZ_item(jk_bz,kgw,jk_ibz,isym_kgw,jik,ph_mkgwt)
253 : !%call kmesh%get_IBZ_item(jk_ibz,kibz,wtk)
254 :
255 : ! TODO: the new version based of get_uug won't suppporte kptgw vector that are not in
256 : ! the IBZ since one should perform the rotation before entering the band loop
257 : ! In the old version, the rotation was done in rho_tw_g
258 : !ABI_CHECK(jik==1,"jik!=1")
259 : !ABI_CHECK(isym_kgw==1,"isym_kgw!=1")
260 : !ABI_CHECK((ABS(ph_mkgwt - cone) < tol12),"ph_mkgwt!")
261 :
262 2275 : spinrot_kgw=Cryst%spinrot(:,isym_kgw)
263 455 : ib1=minbnd; ib2=maxbnd
264 :
265 455 : write(msg,'(2a,3f8.3,2a,2(i3,a))')ch10,&
266 455 : ' Calculating <nk|Sigma_c(omega)|nk> at k = ',kgw(:),ch10,&
267 910 : ' bands n = from ',ib1,' to ',ib2,ch10
268 455 : call wrtout(std_out, msg)
269 :
270 : ! Today we use the same number of eigenvectors irrespective to iw'. Tomorrow we might optimize this further
271 455 : neigmax = Sigp%npwc
272 455 : if (Dtset%gwaclowrank > 0) neigmax = MIN(Dtset%gwaclowrank, Sigp%npwc)
273 :
274 4602 : ABI_CALLOC(w_maxval,(minbnd:maxbnd))
275 :
276 1694 : if (ANY(gwc_ngfft(1:3) /= Wfd%ngfft(1:3))) then
277 42 : call Wfd%change_ngfft(Cryst,Psps,gwc_ngfft)
278 42 : if (dtset%userie == 456) call wfdf%change_ngfft(Cryst, Psps, gwc_ngfft)
279 : end if
280 1820 : gwc_mgfft = MAXVAL(gwc_ngfft(1:3))
281 455 : gwc_fftalga = gwc_ngfft(7)/100 !; gwc_fftalgc=MOD(gwc_ngfft(7),10)
282 :
283 455 : if (Dtset%pawcross==1) mgfftf = MAXVAL(rho_ngfft(1:3))
284 :
285 916 : can_symmetrize = .FALSE.
286 455 : if (Sigp%symsigma>0) then
287 215 : can_symmetrize = .TRUE.
288 106 : if (Sigp%gwcalctyp >= 20) then
289 0 : do spin=1,Wfd%nsppol
290 0 : can_symmetrize(spin) = .not. QP_sym(spin)%failed()
291 0 : if (.not.can_symmetrize(spin)) then
292 : write(msg,'(a,i0,4a)')&
293 0 : " Symmetrization cannot be performed for spin: ",spin,ch10,&
294 0 : " band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg)
295 0 : ABI_WARNING(msg)
296 : end if
297 : end do
298 : end if
299 : ! Symmetrization is available for nspinor = 2, and the difference from symsigma = 0 lies with degenerate bands.
300 : ! By setting symsigma_de 0.000, the difference can be eliminated.
301 106 : if (Wfd%nspinor == 2) then
302 3 : ABI_WARNING("Symmetrization with nspinor = 2 not implemented")
303 : end if
304 : end if
305 :
306 : ! Print type of calculation.
307 455 : call wrtout(std_out, sigma_type_from_key(mod10))
308 :
309 : ! Set up logical flags for Sigma calculation.
310 455 : if (mod10 == SIG_GW_AC) then
311 28 : ABI_CHECK_IEQ(Sigp%gwcomp, 0, "gwcomp with AC not implemented")
312 28 : if (Sigp%gwcalctyp /= 1) then
313 21 : ABI_CHECK(Sigp%gwcalctyp == 21, "gwcalctyp /= 21 not implemented")
314 21 : write(msg,'(a34,i9)')'Constructing Sigma_c(iw) for k = ',ikcalc
315 21 : call wrtout(units, msg)
316 : end if
317 :
318 : write(msg,'(3a,i0,a,i0)')&
319 28 : ' Using a low-rank formula for AC', ch10, &
320 56 : ' Number of epsm1 eigenvectors retained: ',neigmax,' over: ',Sigp%npwc
321 28 : call wrtout(std_out, msg)
322 : end if
323 :
324 : ! Initialize some values
325 455 : nspinor = Wfd%nspinor; npwc = Sigp%npwc
326 4095 : spinor_padc(:,:)=RESHAPE([0, 0, npwc, npwc, 0, npwc, npwc, 0], [2, 4])
327 :
328 1365 : ABI_MALLOC(npoles_missing, (minbnd:maxbnd))
329 3692 : npoles_missing=0
330 :
331 : ! Set tolerance used to decide if a band is empty
332 : ! and normalization of theta_mu_minus_esum. If nsppol == 2, qp_occ $\in [0,1]$
333 455 : tol_empty_in = 0.01
334 455 : call get_fact_spin_tol_empty(wfd%nsppol, wfd%nspinor, tol_empty_in, fact_spin, tol_empty)
335 :
336 : ! Allocate arrays used to accumulate the matrix elements of \Sigma_c over
337 : ! k-points and bands. Note that for AC requires only the imaginary frequencies
338 : !
339 : ! === Define the G-G0 shifts for the FFT of the oscillators ===
340 : ! * Sigp%mG0 gives the MAX G0 component to account for umklapp.
341 : ! * Note the size MAX(Sigp%npwx,npwc).
342 : !
343 : ! === Precalculate the FFT index of $(R^{-1}(r-\tau))$ ===
344 : ! * S=\transpose R^{-1} and k_BZ = S k_IBZ
345 : ! * irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
346 1820 : gwc_nfftot = PRODUCT(gwc_ngfft(1:3))
347 1820 : ABI_MALLOC(irottb,(gwc_nfftot,Cryst%nsym))
348 455 : call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,gwc_ngfft,irottb,iscompatibleFFT)
349 455 : if (.not.iscompatibleFFT) then
350 0 : ABI_WARNING("FFT mesh is not compatible with symmetries. Results might be affected by large errors!")
351 : end if
352 :
353 1820 : ABI_MALLOC(ktabr, (gwc_nfftot,Kmesh%nbz))
354 25523 : do ik_bz=1,Kmesh%nbz
355 25068 : isym = Kmesh%tabo(ik_bz)
356 121257437 : do ifft=1,gwc_nfftot
357 121256982 : ktabr(ifft,ik_bz) = irottb(ifft,isym)
358 : end do
359 : end do
360 455 : ABI_FREE(irottb)
361 :
362 455 : if (Psps%usepaw==1 .and. Dtset%pawcross==1) then
363 0 : nfftf = PRODUCT(rho_ngfft(1:3))
364 0 : ABI_MALLOC(irottb,(nfftf,Cryst%nsym))
365 0 : call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,rho_ngfft,irottb,iscompatibleFFT)
366 :
367 0 : ABI_MALLOC(ktabrf,(nfftf,Kmesh%nbz))
368 0 : do ik_bz=1,Kmesh%nbz
369 0 : isym=Kmesh%tabo(ik_bz)
370 0 : do ifft=1,nfftf
371 0 : ktabrf(ifft,ik_bz)=irottb(ifft,isym)
372 : end do
373 : end do
374 0 : ABI_FREE(irottb)
375 : end if
376 :
377 455 : Sigcij_tab => Sigp%Sigcij_tab(ikcalc,1:Wfd%nsppol)
378 :
379 1302 : got=0
380 2275 : ABI_MALLOC(proc_distrb,(Wfd%mband,Kmesh%nbz,Wfd%nsppol))
381 455 : call sigma_distribute_bks(Wfd,Kmesh,Ltg_k,Qmesh,Wfd%nsppol,can_symmetrize,kgw,Sigp%mg0, my_nbks,proc_distrb,got,global=.TRUE.)
382 :
383 455 : write(msg,'(a,i0,a)')" Will sum ",my_nbks," (b,k,s) states in Sigma_c."
384 455 : call wrtout(std_out, msg)
385 :
386 455 : if (Sigp%gwcomp == 1) then
387 127 : en_high=MAXVAL(qp_ene(Sigp%nbnds,:,:)) + Sigp%gwencomp
388 23 : write(msg,'(6a,e11.4,a)')ch10,&
389 23 : ' Using the extrapolar approximation to accelerate convergence',ch10,&
390 23 : ' with respect to the number of bands included',ch10,&
391 46 : ' with extrapolar energy: ',en_high*Ha_eV,' [eV]'
392 23 : call wrtout(std_out, msg)
393 69 : ABI_MALLOC(wf1swf2_g, (gwc_nfftot*nspinor))
394 :
395 : ! Setup of MPI table for extrapolar contributions.
396 138 : ABI_MALLOC(extrapolar_distrb,(ib1:ib2,ib1:ib2,Kmesh%nbz,Wfd%nsppol))
397 5304 : extrapolar_distrb = xmpi_undefined_rank
398 :
399 46 : do spin=1,Wfd%nsppol
400 580 : do ik_bz=1,Kmesh%nbz
401 4887 : if (ANY(proc_distrb(:,ik_bz,spin) /= xmpi_undefined_rank) ) then ! This BZ point will be calculated.
402 658 : rank_mask = .FALSE. ! The set of node that will treat (k,s).
403 6473 : do band=1,Wfd%mband
404 6158 : rank = proc_distrb(band,ik_bz,spin)
405 6473 : if (rank /= xmpi_undefined_rank) rank_mask(rank+1)=.TRUE.
406 : end do
407 815 : do band2=ib1,ib2
408 1945 : do irow=1,Sigcij_tab(spin)%col(band2)%size1 ! Looping over the non-zero elements of sigma_ij.
409 1130 : band1 = Sigcij_tab(spin)%col(band2)%bidx(irow)
410 1130 : idle = imin_loc(got,mask=rank_mask)
411 1130 : got(idle) = got(idle)+1
412 1630 : extrapolar_distrb(band1,band2,ik_bz,spin) = idle-1
413 : end do
414 : end do
415 : end if
416 : end do
417 : end do
418 :
419 5304 : write(msg,'(a,i0,a)')" Will treat ",COUNT(extrapolar_distrb==Wfd%my_rank)," extrapolar terms."
420 23 : call wrtout(std_out, msg)
421 : end if
422 :
423 455 : if (Sigp%gwcomp == 2) then
424 : ! Allocate memory to store static screening in the IBZ and workspace array for the value at qq_bz in the BZ.
425 : ! TODO: Store only the iq_ibz needed by this MPI rank to save memory.
426 0 : call wrtout(units, " Activating COH remainder technique")
427 :
428 0 : ABI_MALLOC(wc0_qibz, (qmesh%nibz))
429 0 : ABI_MALLOC(wc0_qbz, (npwc, npwc))
430 0 : ABI_MALLOC(epsm1_ggw, (npwc, npwc, 1))
431 :
432 0 : scr_iomode = iomode_from_fname(dtfil%fnameabi_scr)
433 0 : do iq_ibz=1,qmesh%nibz
434 : call read_screening("inverse_dielectric_function", dtfil%fnameabi_scr, &
435 0 : npwc, 1, 1, epsm1_ggw, scr_iomode, wfd%comm, iqiA=iq_ibz)
436 :
437 : ! Keep static limit of Wc in memory.
438 : !do ii=1,npwc
439 : ! epsm1_ggw(ii,ii,1) = epsm1_ggw(ii,ii,1) - one
440 : !end do
441 0 : ABI_MALLOC(wc0_qibz(iq_ibz)%vals, (npwc, npwc))
442 0 : wc0_qibz(iq_ibz)%vals = epsm1_ggw(:,:,1)
443 : end do
444 :
445 0 : ABI_FREE(epsm1_ggw)
446 : !call hscr%free()
447 0 : ABI_MALLOC(coh_g, (npwc))
448 : end if
449 :
450 1820 : ABI_MALLOC(rhotwg_ki, (npwc*nspinor, minbnd:maxbnd))
451 122463 : rhotwg_ki=czero_gw
452 1365 : ABI_MALLOC(rhotwg, (npwc*nspinor))
453 910 : ABI_MALLOC(rhotwgp, (npwc*nspinor))
454 1365 : ABI_MALLOC(vc_sqrt_qbz, (npwc))
455 :
456 455 : if (epsm1%mqmem == 0) then
457 : ! Use out-of-core solution for epsilon.
458 19 : ABI_COMMENT('Reading q-slices from file. Slower but less memory.')
459 : end if
460 :
461 : ! Additional allocations for PAW.
462 455 : if (Psps%usepaw==1) then
463 84 : ABI_MALLOC(Cprj_ksum,(Cryst%natom,nspinor))
464 9 : call pawcprj_alloc(Cprj_ksum,0,Wfd%nlmn_atm)
465 : !
466 : ! For the extrapolar method we need the onsite terms of the PW in the FT mesh.
467 : ! gw_gfft is the set of plane waves in the FFT Box for the oscillators.
468 9 : if (Sigp%gwcomp==1) then
469 0 : ABI_MALLOC(gw_gfft,(3,gwc_nfftot))
470 0 : q0=zero
471 0 : call get_gfft(gwc_ngfft,q0,Cryst%gmet,gw_gsq,gw_gfft)
472 0 : ABI_MALLOC(Pwij_fft,(Psps%ntypat))
473 0 : call pawpwij_init(Pwij_fft,gwc_nfftot, [zero,zero,zero], gw_gfft,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
474 : end if
475 : end if ! usepaw==1
476 :
477 455 : if (mod10 == SIG_GW_AC) then
478 : ! Calculate Gauss-Legendre quadrature knots and weights for analytic continuation.
479 168 : ABI_MALLOC(rhotw_epsm1_rhotw, (minbnd:maxbnd, minbnd:maxbnd, epsm1%nomega_i, Sigp%nsig_ab))
480 168 : ABI_MALLOC(conv_rhotw_epsm1_rhotw, (minbnd:maxbnd, minbnd:maxbnd, epsm1%nomega_i_conv, Sigp%nsig_ab))
481 :
482 28 : select case (epsm1%hscr%iw_mesh_type)
483 : case ("gauss_legendre")
484 28 : call coeffs_gausslegint(zero, one, gl_knots, gl_wts, epsm1%nomega_i)
485 :
486 : ! To calculate \int_0^\infty domegap f(omegap), we calculate \int_0^1 dz f(1/z-1)/z^2.
487 1728 : omegap(:) = one / gl_knots(:) - one
488 1728 : omegap2(:) = omegap(:) ** 2
489 :
490 28 : ierr = 0
491 1728 : do io=1,epsm1%nomega_i
492 : ! First frequencies are always real
493 1728 : if (ABS(AIMAG(one*epsm1%omega(epsm1%nomega_r+io))-omegap(io)) > 0.0001) then
494 0 : ierr = ierr + 1
495 0 : if (Wfd%my_rank == Wfd%master) then
496 0 : if (io == 1) write(std_out, "(a)")"omega_file, gauss_legendre_omega (ev)"
497 0 : write(std_out,*)io, AIMAG(epsm1%omega(epsm1%nomega_r+io)) * Ha_eV, omegap(io) * Ha_eV
498 : end if
499 : end if
500 : end do
501 :
502 28 : if (ierr /= 0) then
503 0 : write(std_out, *)"epsm1%nomega_r:", epsm1%nomega_r, "epsm1%nomega_i:", epsm1%nomega_i
504 : write(msg,'(3a)')&
505 0 : 'Frequencies in the SCR file are not compatible with the analytic continuation with gauss-legendre mesh.',ch10,&
506 0 : 'Verify the frequencies in the SCR file. '
507 0 : ABI_ERROR(msg)
508 : end if
509 :
510 : case ("minimax")
511 : ! Nothing to do at this level
512 0 : omegap_cplx(:) = epsm1%omega(epsm1%nomega_r+1: epsm1%nomega_r+epsm1%nomega_i)
513 0 : omegap(:) = AIMAG(omegap_cplx(:))
514 : !write(std_out,*) "omegap_cplx", omegap_cplx
515 : !omegap2(:) = zero
516 : !write(std_out,*)"omega_i", Sr%omega_i(:)
517 :
518 : case default
519 28 : ABI_ERROR(sjoin("Invalid iw_mesh_type:" , epsm1%hscr%iw_mesh_type))
520 : end select
521 :
522 28 : if (epsm1%nomega_i_conv > 0) then
523 1 : call coeffs_gausslegint(zero, one, conv_gl_knots, conv_gl_wts, epsm1%nomega_i_conv)
524 31 : conv_omegap(:) = one / conv_gl_knots(:) - one
525 31 : conv_omegap2(:) = conv_omegap(:) ** 2
526 : end if
527 :
528 28 : if (epsm1%use_mpi_shared_win) then
529 : #define _MOK(integer) int(integer, kind=XMPI_OFFSET_KIND)
530 0 : ad_count = _MOK(2*npwc) * _MOK(npwc) * _MOK(epsm1%nomega_i)
531 0 : call epsm1%shared_comm%allocate_shared_master(ad_count, gwp, xmpi_info_null, void_ptr, ac_epsm1cqwz2_win)
532 0 : call c_f_pointer(void_ptr, ac_epsm1cqwz2, shape=[npwc, npwc, epsm1%nomega_i])
533 : else
534 3043880 : ABI_CALLOC(ac_epsm1cqwz2, (npwc, npwc, epsm1%nomega_i))
535 : end if
536 : end if
537 :
538 : ! Calculate total number of frequencies and allocate related arrays.
539 : ! sigcme2 is used to accumulate the diagonal matrix elements over k-points and
540 : ! GW bands, used only in case of ppmodel 3 and 4 (TODO save memory)
541 455 : nomega_tot = Sr%nomega_r + Sr%nomega4sd
542 25658 : ABI_CALLOC(sigcme2, (nomega_tot, ib1:ib2))
543 4798 : ABI_CALLOC(sigcme_3, (nomega_tot))
544 6216 : ABI_CALLOC(sigctmp,(nomega_sigc, Sigp%nsig_ab))
545 455 : if (mod10 /= SIG_GW_AC) then
546 1708 : ABI_MALLOC(sigc_ket, (npwc*nspinor, nomega_sigc))
547 : end if
548 :
549 : #if 0
550 : !TODO gmatteo: these arrays are never used in practice. Should we remove them?
551 : ! Arrays storing the contribution given by the Hermitian/anti-Hermitian part of \Sigma_c
552 : ABI_MALLOC(aherm_sigc_ket, (npwc*nspinor, nomega_sigc))
553 : ABI_MALLOC(herm_sigc_ket, (npwc*nspinor, nomega_sigc))
554 : #endif
555 :
556 247042 : sigcme_tmp = czero
557 :
558 673765 : ABI_CALLOC(sigc, (2,nomega_sigc,ib1:ib2,ib1:ib2,Wfd%nsppol*Sigp%nsig_ab))
559 :
560 455 : if (any(mod10 == [SIG_QPGW_PPM, SIG_QPGW_CD])) then
561 752 : ABI_MALLOC(ket1, (npwc*nspinor, nomega_tot))
562 564 : ABI_MALLOC(ket2, (npwc*nspinor, nomega_tot))
563 : end if
564 1365 : ABI_MALLOC(omegame0i, (nomega_tot))
565 :
566 : ! Here we divide the states where the QP energies are required into degenerate groups
567 : ! Note however that this approach is not based on group theory, and it might lead to
568 : ! spurious results in case of accidental degeneracies.
569 455 : nq_summed=Kmesh%nbz
570 455 : if (Sigp%symsigma > 0) then
571 212 : call Ltg_k%print([std_out], prtvol=Dtset%prtvol)
572 6846 : nq_summed = sum(Ltg_k%ibzq(:))
573 : !
574 : ! Find number of degenerate subspaces and number of bands in each subspace
575 : ! The tolerance is a little bit arbitrary (0.001 eV)
576 : ! It could be reduced, in particular in case of nearly accidental degeneracies.
577 530 : ABI_MALLOC(degtab,(ib1:ib2,ib1:ib2,Wfd%nsppol))
578 5999 : degtab=0
579 215 : do spin=1,Wfd%nsppol
580 884 : do ib=ib1,ib2
581 5893 : do jb=ib1,ib2
582 5784 : if (abs(qp_ene(ib,jk_ibz,spin)-qp_ene(jb,jk_ibz,spin)) < dtset%symsigma_de) degtab(ib,jb,spin)=1
583 : end do
584 : end do
585 : end do
586 : end if !symsigma
587 :
588 455 : write(msg,'(2a,i0,a)')ch10,' calculation status ( ',nq_summed,' to be completed):'
589 455 : call wrtout(std_out, msg)
590 :
591 : ! Here we have a problem in case of CD or AC since epsm1q might be huge.
592 : ! For this reason we use MPI shared memory for epsm1%epsm1_qbz(npwc, npwc, epsm1%nomega) inside shared_comm.
593 : ! TODO if single q (ex molecule) dont allocate epsm1q, avoid waste of memory
594 455 : if (ANY(mod10 == [SIG_GW_AC, SIG_GW_CD, SIG_QPGW_CD])) then
595 127 : if (.not. (mod10==SIG_GW_CD .and. epsm1%mqmem == 0)) then
596 127 : call epsm1%malloc_epsm1_qbz(npwc, epsm1%nomega)
597 : end if
598 : end if
599 :
600 : ! TODO epsm1_trcc_qbz is needed for SIG_GW_CD with symmetries since
601 : ! the Hermitian and the anti-Hermitian part have to be symmetrized in a different way.
602 : ! Clearly this increases significantly the memory requirements.
603 455 : if (mod10 == SIG_QPGW_CD) then
604 80 : ABI_MALLOC_OR_DIE(epsm1_trcc_qbz, (npwc, npwc, epsm1%nomega), ierr)
605 64 : ABI_MALLOC(epsm1_tmp, (npwc, npwc))
606 : end if
607 :
608 1365 : ABI_MALLOC(igfftcg0, (Gsph_Max%ng))
609 1365 : ABI_MALLOC(ur_ibz, (gwc_nfftot*nspinor))
610 910 : ABI_MALLOC(usr_bz, (gwc_nfftot*nspinor))
611 :
612 455 : if (Dtset%pawcross==1) then
613 0 : ABI_MALLOC(igfftfcg0, (Gsph_c%ng))
614 0 : ABI_MALLOC(ur_ae_sum, (nfftf*nspinor))
615 0 : ABI_MALLOC(ur_ae_onsite_sum, (nfftf*nspinor))
616 0 : ABI_MALLOC(ur_ps_onsite_sum, (nfftf*nspinor))
617 : end if
618 455 : call timab(432,2,tsec) ! Init
619 455 : call pstat_proc%print(_PSTAT_ARGS_)
620 :
621 : ! If epsm1 is MPI-shared, we have to close the RMA epoch.
622 455 : if (epsm1%use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOPRECEDE, epsm1%epsm1_win, ierr)
623 :
624 : ! Loop over collinear spins.
625 916 : do spin=1,Wfd%nsppol
626 6600 : if (ALL(proc_distrb(:,:,spin)/=Wfd%my_rank)) CYCLE
627 461 : call timab(433,1,tsec) ! Init spin
628 :
629 : ! Load wavefunctions for GW corrections
630 : ! TODO: Rotate the functions here instead of calling rho_tw_g
631 1844 : ABI_MALLOC(wfr_bdgw, (gwc_nfftot*nspinor, ib1:ib2))
632 :
633 461 : if (dtset%userie == 456) then
634 0 : call wrtout(std_out, "Taking states from Sigma^c_nk from supercell WFK file")
635 0 : call wfdf%get_many_ur([(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw)
636 : else
637 6959 : call wfd%get_many_ur([(jb, jb=ib1,ib2)], jk_ibz, spin, wfr_bdgw)
638 : end if
639 :
640 461 : if (wfd%usepaw == 1) then
641 : ! Load cprj for GW states, note the indexing.
642 9 : dimcprj_gw=nspinor*(ib2-ib1+1)
643 216 : ABI_MALLOC(Cprj_kgw,(Cryst%natom,ib1:ib1+dimcprj_gw-1))
644 9 : call pawcprj_alloc(Cprj_kgw,0,Wfd%nlmn_atm)
645 9 : ibsp=ib1
646 55 : do jb=ib1,ib2
647 46 : call Wfd%get_cprj(jb,jk_ibz,spin,Cryst,Cprj_ksum,sorted=.FALSE.)
648 46 : call paw_symcprj(jk_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
649 46 : call pawcprj_copy(Cprj_ksum,Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
650 55 : ibsp=ibsp+nspinor
651 : end do
652 9 : if (Dtset%pawcross==1) then
653 0 : ABI_MALLOC(ur_ae_bdgw, (nfftf*nspinor,ib1:ib2))
654 0 : ABI_MALLOC(ur_ae_onsite_bdgw, (nfftf*nspinor,ib1:ib2))
655 0 : ABI_MALLOC(ur_ps_onsite_bdgw, (nfftf*nspinor,ib1:ib2))
656 0 : do jb=ib1,ib2
657 : call Wfdf%paw_get_aeur(jb,jk_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
658 0 : ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
659 0 : ur_ae_bdgw(:,jb)=ur_ae_sum
660 0 : ur_ae_onsite_bdgw(:,jb)=ur_ae_onsite_sum
661 0 : ur_ps_onsite_bdgw(:,jb)=ur_ps_onsite_sum
662 : end do
663 : end if
664 : end if ! usepaw
665 :
666 461 : call timab(433,2,tsec) ! Init spin
667 :
668 : ! ==========================================
669 : ! ==== Fat loop over k_i in the full BZ ====
670 : ! ==========================================
671 :
672 25913 : do ik_bz=1,Kmesh%nbz
673 : ! Parallelization over k-points and spin. For the spin there is another check in the inner loop
674 184381 : if (ALL(proc_distrb(:,ik_bz,spin)/=Wfd%my_rank)) CYCLE
675 :
676 19839 : call timab(434,1,tsec) ! initq
677 19839 : print_time = wfd%my_rank == 0 .and. (ik_bz < LOG_MODK .or. mod(ik_bz, LOG_MODK) == 0)
678 3658 : if (print_time) call cwtime(cpu_k, wall_k, gflops_k, "start")
679 :
680 : ! Find the corresponding irreducible k-point
681 19839 : call kmesh%get_BZ_item(ik_bz, ksum, ik_ibz, isym_ki, iik, ph_mkt)
682 99195 : spinrot_kbz(:)=Cryst%spinrot(:,isym_ki)
683 :
684 : ! Identify q and G0 where q + G0 = k_GW - k_i
685 79356 : kgw_m_ksum=kgw-ksum
686 19839 : call findqg0(iq_bz, g0, kgw_m_ksum, Qmesh%nbz, Qmesh%bz, Sigp%mG0)
687 :
688 : ! If symsigma, symmetrize the matrix elements.
689 : ! Sum only q"s in IBZ_k. In this case elements are weighted
690 : ! according to wtqp and wtqm. wtqm is for time-reversal.
691 19839 : wtqp = 1; wtqm = 0
692 19839 : if (can_symmetrize(spin)) then
693 1319 : if (Ltg_k%ibzq(iq_bz)/=1) CYCLE
694 1319 : wtqp = 0; wtqm = 0
695 72287 : do isym=1,Ltg_k%nsym_sg
696 70968 : wtqp = wtqp + Ltg_k%wtksym(1,isym,iq_bz)
697 72287 : wtqm = wtqm + Ltg_k%wtksym(2,isym,iq_bz)
698 : end do
699 : end if
700 :
701 : ! Find the corresponding irred q-point.
702 19839 : call qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym_q, itim_q)
703 19839 : q_is_gamma = normv(qbz, Cryst%gmet, "G") < GW_TOLQ0
704 :
705 : ! Tables for the FFT of the oscillators.
706 : ! a) FFT index of the G-G0.
707 : ! b) gw_gbound table for the zero-padded FFT performed in rhotwg.
708 79356 : ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2))
709 19839 : call Gsph_c%fft_tabs(g0,gwc_mgfft,gwc_ngfft,use_padfft,gw_gbound,igfftcg0)
710 :
711 19839 : if (ANY(gwc_fftalga == [2, 4])) use_padfft=0 ! Pad-FFT is not coded in rho_tw_g
712 19839 : if (use_padfft == 0) then
713 33 : ABI_FREE(gw_gbound)
714 66 : ABI_MALLOC(gw_gbound,(2*gwc_mgfft+8,2*use_padfft))
715 : end if
716 :
717 19839 : if (Dtset%pawcross==1) then
718 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2))
719 0 : call Gsph_c%fft_tabs(g0,mgfftf,rho_ngfft,use_padfftf,gboundf,igfftfcg0)
720 0 : if (any(gwc_fftalga == [2, 4]) ) use_padfftf=0
721 0 : if (use_padfftf == 0) then
722 0 : ABI_FREE(gboundf)
723 0 : ABI_MALLOC(gboundf, (2*mgfftf+8,2*use_padfftf))
724 : end if
725 : end if
726 :
727 : ! Evaluate oscillator matrix elements
728 : ! $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form.
729 19839 : if (Psps%usepaw==1) then
730 774 : ABI_MALLOC(Pwij_qg,(Psps%ntypat))
731 192 : q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
732 192 : call pawpwij_init(Pwij_qg,npwc,q0,Gsph_c%gvec,Cryst%rprimd,Psps,Pawtab,Paw_pwff)
733 : end if
734 :
735 19839 : if (epsm1%mqmem == 0) then
736 : ! Read q-slice of epsilon^{-1}|chi0 in epsm1%epsm1(:,:,:,1) (much slower but less memory).
737 290 : call epsm1%get_epsm1(Vcp,0,0,Dtset%iomode,xmpi_comm_self,iqibzA=iq_ibz)
738 290 : if (sigp%needs_ppm()) then
739 290 : if (Wfd%usepaw==1.and.PPm%userho==1) then
740 : ! Use PAW AE rhor.
741 : call PPm%setup(Cryst,Qmesh,epsm1%npwe,epsm1%nomega,epsm1%omega,epsm1%epsm1,rho_nfftot,Gsph_c%gvec,&
742 0 : rho_ngfft,aepaw_rhor(:,1),iqiA=iq_ibz)
743 : else
744 : call PPm%setup(Cryst,Qmesh,epsm1%npwe,epsm1%nomega,epsm1%omega,epsm1%epsm1,rho_nfftot,Gsph_c%gvec,&
745 290 : rho_ngfft,rhor(:,1),iqiA=iq_ibz)
746 : end if
747 : end if
748 : end if
749 :
750 : ! Symmetrize PPM parameters and epsm1 (q_IBZ --> q_BZ):
751 : ! NOTE:
752 : ! - We are not considering umklapp with G0/=0. In this case,
753 : ! indeed the equation is different since we have to use G-G0.
754 : ! A check, however, is performed in sigma.
755 : ! - If gwcomp==1 and mod10 in [1,2,9], one needs both to set up botsq and epsm1_q
756 19839 : if (sigp%needs_ppm()) then
757 14122 : call PPm%get_qbz(Gsph_c, Qmesh, iq_bz, botsq, otq, eig)
758 : end if
759 :
760 19839 : if (any(mod10 == [SIG_GW_AC, SIG_GW_CD, SIG_QPGW_CD])) then
761 :
762 : ! Numerical integration or model GW with contour deformation or Analytic Continuation
763 : ! TODO In case of AC we should symmetrize only the imaginary frequencies
764 5717 : if (mod10==SIG_GW_CD .and. epsm1%mqmem == 0) then
765 : ! Do in-place symmetrization.
766 0 : call epsm1%rotate_iqbz_inplace(iq_bz, epsm1%nomega, npwc, Gsph_c, Qmesh, remove_exchange=.TRUE.)
767 : else
768 : ! This call sets the value of epsm1%epsm1_qbz(npwc, npwc, epsm1%nomega)
769 5717 : call epsm1%rotate_iqbz(iq_bz, epsm1%nomega, npwc, Gsph_c, Qmesh, remove_exchange=.TRUE.)
770 : end if
771 :
772 5717 : if (mod10 == SIG_GW_AC) then
773 726 : call timab(444,1,tsec)
774 : ! ac_lrk_diag
775 : ! Important to set to zero here for all procs since we're going to use a dirty reduction later
776 : ! The reduction 'xmpi_sum' does not induce a significant performance loss in the tested systems
777 51726 : neig(:) = 0
778 2178 : ABI_MALLOC(epsm1_eig, (Sigp%npwc))
779 726 : if (epsm1%use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOPRECEDE, ac_epsm1cqwz2_win, ierr) ! Start the RMA epoch.
780 :
781 51726 : do iiw=1,epsm1%nomega_i
782 :
783 51000 : if (.not. epsm1%use_mpi_shared_win) then
784 : ! Use the MPI procs in wfd%comm to parallelize over iw'.
785 51000 : if (dtset%gwpara == 2 .and. MODULO(iiw-1, Wfd%nproc) /= Wfd%my_rank) then
786 0 : ac_epsm1cqwz2(:,:,iiw) = zero
787 : CYCLE
788 : end if
789 : else
790 : ! Use the MPI procs in shared_comm to parallelize over iw'.
791 0 : if (dtset%gwpara == 2 .and. epsm1%shared_comm%skip(iiw-1)) CYCLE
792 : end if
793 :
794 51000 : select case (epsm1%hscr%iw_mesh_type)
795 : case ("gauss_legendre")
796 : ! Prepare the integration weights w_i 1/z_i^2 f(1/z_i-1)..
797 : ! The first frequencies are always real, skip them.
798 70071240 : ac_epsm1cqwz2(:,:,iiw) = epsm1%epsm1_qbz(:,:,epsm1%nomega_r+iiw)
799 :
800 : case ("minimax")
801 : ! ac_epsm1cqwz2(:,:,iiw) = epsm1%hscr%omega_wgs(epsm1%nomega_r+iiw) * epsm1%epsm1_qbz(:,:,epsm1%nomega_r+iiw)
802 0 : ac_epsm1cqwz2(:,:,iiw) = epsm1%epsm1_qbz(:,:,epsm1%nomega_r+iiw)
803 :
804 : case default
805 51000 : ABI_ERROR(sjoin("Invalid iw_mesh_type:", epsm1%hscr%iw_mesh_type))
806 : end select
807 :
808 : ! (epsm1-1) has negative eigenvalues, after diago, they will be sorted starting from the most negative.
809 51000 : call xheev('V','L', npwc, ac_epsm1cqwz2(:,:,iiw), epsm1_eig)
810 :
811 : ! Eliminate the spurious positive eigenvalues that may occur in harsh conditions.
812 1554000 : neig(iiw) = MIN(COUNT(epsm1_eig(:) < -1.0e-10_dp), neigmax)
813 : !write(std_out, *) "neig(iiw): ", neig(iiw), ", neigmax: ", neigmax
814 :
815 1532328 : do ilwrk=1,neig(iiw)
816 68717504 : ac_epsm1cqwz2(:,ilwrk,iiw) = ac_epsm1cqwz2(:,ilwrk,iiw) * SQRT( -epsm1_eig(ilwrk) )
817 : end do
818 : end do ! iiw
819 :
820 726 : ABI_FREE(epsm1_eig)
821 :
822 726 : if (Dtset%gwpara == 2) then
823 : ! FIXME: It seems that non all the procs get here if nband is small!
824 : !call wrtout(std_out, "AC xmpi_sum begin")
825 726 : if (.not. epsm1%use_mpi_shared_win) then
826 726 : call xmpi_sum(neig, Wfd%comm, ierr)
827 726 : call xmpi_sum(ac_epsm1cqwz2, Wfd%comm, ierr)
828 : else
829 : ! No need to MPI_SUM ac_epsm1cqwz2_win as we're using MPI shared memory.
830 0 : call xmpi_sum(neig, epsm1%shared_comm%value, ierr)
831 0 : call xmpi_win_fence(XMPI_MODE_NOSUCCEED, ac_epsm1cqwz2_win, ierr) ! Close the RMA epoch.
832 : end if
833 : !call wrtout(std_out, "AC xmpi_sum end")
834 : end if
835 :
836 726 : call timab(444,2,tsec) ! ac_lrk_diag
837 : end if
838 :
839 5717 : if (mod10==SIG_QPGW_CD) then
840 : ! For model GW we need transpose(conjg(epsm1_qbz))
841 4752 : do io=1,epsm1%nomega
842 1041120 : epsm1_tmp(:,:) = GWPC_CONJG(epsm1%epsm1_qbz(:,:,io))
843 1041552 : epsm1_trcc_qbz(:,:,io) = TRANSPOSE(epsm1_tmp)
844 : end do
845 : end if
846 : end if ! gwcalctyp
847 :
848 19839 : if (dtset%gwcomp == 2) then
849 : ! Compute static limit at pp_bz from the symmetrical image in the IBZ
850 0 : call em1_symmetrize_op(iq_bz, npwc, 1, gsph_c, qmesh, wc0_qibz(iq_ibz)%vals, wc0_qbz)
851 : end if
852 :
853 : ! Get Fourier components of the Coulomb interaction in the BZ
854 : ! In 3D systems, neglecting umklapp: vc(Sq,sG) = vc(q,G) = 4pi/|q+G|**2
855 : ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
856 644726 : do ig=1,npwc
857 644726 : vc_sqrt_qbz(Gsph_c%rottb(ig,itim_q,isym_q)) = Vcp%vc_sqrt(ig,iq_ibz)
858 : end do
859 :
860 19839 : call timab(434,2,tsec) ! initq
861 19839 : call timab(445,1,tsec) ! loop
862 :
863 : ! Sum over bands
864 297087 : do ib_sum=1,Sigp%nbnds
865 :
866 : ! MPI Parallelism over spin.
867 277248 : if (proc_distrb(ib_sum,ik_bz,spin) /= wfd%my_rank) CYCLE
868 :
869 212400 : call wfd%get_ur(ib_sum, ik_ibz, spin, ur_ibz)
870 :
871 212400 : if (Psps%usepaw==1) then
872 : ! Load cprj for point ksum, this spin or spinor and *THIS* band.
873 : ! TODO MG I could avoid doing this but I have to exchange spin and bands ???
874 : ! For sure there is a better way to do this!
875 4818 : call Wfd%get_cprj(ib_sum, ik_ibz, spin, Cryst, Cprj_ksum, sorted=.FALSE.)
876 4818 : call paw_symcprj(ik_bz,nspinor,1,Cryst,Kmesh,Pawtab,Pawang,Cprj_ksum)
877 4818 : if (Dtset%pawcross==1) then
878 : call Wfdf%paw_get_aeur(ib_sum,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,&
879 0 : ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
880 : end if
881 : end if
882 :
883 212400 : call timab(436,2,tsec) ! (1)
884 212400 : call timab(437,1,tsec) ! rho_tw_g
885 :
886 : ! =====================================================
887 : ! Get all <k-q,ib_sum,s|e^{-i(q+G).r}|s,jb,k>, at once
888 : ! =====================================================
889 1494118 : do jb=ib1,ib2
890 :
891 : call rho_tw_g(nspinor,npwc,gwc_nfftot,ndat1,gwc_ngfft,1,use_padfft,igfftcg0,gw_gbound,&
892 : ur_ibz ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, &
893 : wfr_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw,&
894 1281718 : nspinor,rhotwg_ki(:,jb))
895 :
896 1281718 : if (Psps%usepaw==1) then
897 : ! Add on-site contribution, projectors are already in BZ !TODO Recheck this!
898 30936 : i2=jb; if (nspinor==2) i2=(2*jb-1)
899 30936 : spad=(nspinor-1)
900 : call paw_rho_tw_g(cryst,Pwij_qg,npwc,nspinor,nspinor,Gsph_c%gvec, &
901 30936 : Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),rhotwg_ki(:,jb))
902 :
903 30936 : if (Dtset%pawcross==1) then ! Add paw cross term
904 : call paw_cross_rho_tw_g(nspinor,npwc,nfftf,rho_ngfft,1,use_padfftf,igfftfcg0,gboundf,&
905 : ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum,iik,ktabrf(:,ik_bz),ph_mkt,spinrot_kbz,&
906 : ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
907 0 : nspinor,rhotwg_ki(:,jb))
908 : end if
909 : end if
910 :
911 : ! Multiply by the square root of the Coulomb term
912 : ! In 3-D systems, the factor sqrt(4pi) is included)
913 2577964 : do ii=1,nspinor
914 1296246 : spad = (ii-1) * npwc
915 55017434 : rhotwg_ki(spad+1:spad+npwc,jb) = rhotwg_ki(spad+1:spad+npwc,jb)*vc_sqrt_qbz(1:npwc)
916 : end do
917 :
918 : ! === Treat analytically the case q --> 0 ===
919 : ! * The oscillator is evaluated at q=O as it is considered constant in the small cube around Gamma
920 : ! while the Colulomb term is integrated out.
921 : ! * In the scalar case we have nonzero contribution only if ib_sum==jb
922 : ! * For nspinor==2 evaluate <ib_sum,up|jb,up> and <ib_sum,dwn|jb,dwn>,
923 : ! impose orthonormalization since npwwfn might be < npwvec.
924 1494118 : if (ik_bz==jk_bz) then
925 50419 : if (nspinor==1) then
926 47907 : rhotwg_ki(1,jb)=czero_gw
927 47907 : if (ib_sum==jb) rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp)
928 :
929 : else
930 2512 : npw_k = Wfd%npwarr(ik_ibz)
931 2512 : rhotwg_ki(1, jb) = zero; rhotwg_ki(npwc+1, jb) = zero
932 2512 : if (ib_sum == jb) then
933 82 : ABI_CHECK(Wfd%get_wave_ptr(ib_sum, ik_ibz, spin, wave_sum, msg) == 0, msg)
934 82 : cg_sum => wave_sum%ug
935 82 : ABI_CHECK(Wfd%get_wave_ptr(jb, jk_ibz, spin, wave_jb, msg) == 0, msg)
936 82 : cg_jb => wave_jb%ug
937 82 : ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1)
938 82 : rhotwg_ki(1,jb)=CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
939 82 : ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1)
940 82 : rhotwg_ki(npwc+1,jb) = CMPLX(SQRT(Vcp%i_sz),0.0_gwp) * real(ctmp)
941 : ! PAW is missing
942 : end if
943 : end if
944 : end if
945 :
946 : end do ! jb Got all matrix elements from ib1 up to ib2.
947 :
948 212400 : theta_mu_minus_e0i = fact_spin*qp_occ(ib_sum,ik_ibz,spin)
949 :
950 : ! Starting point to evaluate the derivative of Sigma and the Spectral function
951 212400 : e0i = qp_ene(ib_sum,ik_ibz,spin)
952 :
953 : ! Frequencies for the spectral function, e0i=qp_ene(ib_sum,ik_ibz,spin)
954 : ! FIXME the interval is not centered on eoi ! WHY?
955 212400 : if (Sr%nomega_r>0) then
956 3435720 : omegame0i(1:Sr%nomega_r)=DBLE(Sr%omega_r(1:Sr%nomega_r))-e0i
957 : !omegame0i = arth(e0i - gwr%wr_step * (sr%nomega_r / 2), gwr%wr_step, sr%nomega_r)
958 : end if
959 :
960 212400 : call timab(437,2,tsec) ! rho_tw_g
961 :
962 212400 : if (mod10 == SIG_GW_AC) then
963 6756 : call timab(443,1,tsec) ! ac_lrk_appl
964 28354712 : rhotw_epsm1_rhotw(:,:,:,:) = czero_gw
965 13512 : do iab=1,Sigp%nsig_ab
966 433716 : do iiw=1,epsm1%nomega_i
967 1707840 : ABI_MALLOC(epsm1_sqrt_rhotw, (neig(iiw), minbnd:maxbnd))
968 1280880 : ABI_MALLOC(rhotw_eqsm1_sqrt, (minbnd:maxbnd, neig(iiw)))
969 : ! epsm1_sqrt_rhotw = SQRT(epsm1) * rho_tw
970 : ! rhotw_eqsm1_sqrt = rho_tw^H * SQRT(epsm1)
971 426960 : iiab = spinor_idxs(1,iab); jiab = spinor_idxs(2,iab)
972 :
973 : call xgemm('C','N',neig(iiw),maxbnd-minbnd+1,npwc,cone_gw,ac_epsm1cqwz2(:,:,iiw),npwc,&
974 426960 : rhotwg_ki((iiab-1)*npwc+1:iiab*npwc,:),npwc,czero_gw,epsm1_sqrt_rhotw,neig(iiw))
975 : call xgemm('C','N',maxbnd-minbnd+1,neig(iiw),npwc,cone_gw,&
976 : rhotwg_ki((jiab-1)*npwc+1:jiab*npwc,:),npwc,ac_epsm1cqwz2(:,:,iiw),npwc,czero_gw,&
977 426960 : rhotw_eqsm1_sqrt,maxbnd-minbnd+1)
978 : call xgemm('N','N',maxbnd-minbnd+1, maxbnd-minbnd+1, neig(iiw), cone_gw, &
979 : rhotw_eqsm1_sqrt, maxbnd-minbnd+1, epsm1_sqrt_rhotw, neig(iiw), &
980 426960 : czero_gw, rhotw_epsm1_rhotw(:,:,iiw,iab), maxbnd-minbnd+1)
981 : ! call xherk('L','C',maxbnd-minbnd+1,neig(iiw),one_gw,epsm1_sqrt_rhotw,neig(iiw),zero_gw,&
982 : ! rhotw_epsm1_rhotw(:,:,iiw,iab), maxbnd-minbnd+1)
983 :
984 : ! Get the upper part of rhotw_epsm1_rhotw that is hermitian by construction
985 3628840 : do jb=minbnd,maxbnd
986 14384080 : do kb=jb+1,maxbnd
987 13957120 : rhotw_epsm1_rhotw(jb,kb,iiw,iab) = CONJG(rhotw_epsm1_rhotw(kb,jb,iiw,iab))
988 : end do
989 : end do
990 426960 : ABI_FREE(epsm1_sqrt_rhotw)
991 433716 : ABI_FREE(rhotw_eqsm1_sqrt)
992 : end do
993 13512 : if (epsm1%nomega_i_conv > 0) then
994 1620 : do jb=minbnd,maxbnd
995 13140 : do kb=minbnd,maxbnd
996 11520 : select case (epsm1%hscr%iw_mesh_type)
997 : case ("gauss_legendre")
998 241920 : tmp_rhotw_epsm1_rhotw = rhotw_epsm1_rhotw(jb,kb,epsm1%nomega_i:1:-1,iab)
999 241920 : tmp_omegap = omegap(epsm1%nomega_i:1:-1)
1000 357120 : tmp_conv_omegap = conv_omegap(epsm1%nomega_i_conv:1:-1)
1001 : call spline_c(epsm1%nomega_i, epsm1%nomega_i_conv, &
1002 : tmp_omegap, tmp_conv_omegap, &
1003 : tmp_conv_rhotw_epsm1_rhotw, &
1004 : tmp_rhotw_epsm1_rhotw, &
1005 11520 : extrapolate=.TRUE.)
1006 : case ("minimax")
1007 0 : tmp_rhotw_epsm1_rhotw = rhotw_epsm1_rhotw(jb,kb,:,iab)
1008 0 : tmp_omegap = omegap(:)
1009 0 : tmp_conv_omegap = conv_omegap(epsm1%nomega_i_conv:1:-1)
1010 : call spline_c(epsm1%nomega_i, epsm1%nomega_i_conv, &
1011 : tmp_omegap, tmp_conv_omegap, &
1012 : tmp_conv_rhotw_epsm1_rhotw, &
1013 : tmp_rhotw_epsm1_rhotw, &
1014 11520 : extrapolate=.TRUE.)
1015 : end select
1016 358560 : conv_rhotw_epsm1_rhotw(jb,kb,:,iab) = tmp_conv_rhotw_epsm1_rhotw(epsm1%nomega_i_conv:1:-1)
1017 : end do
1018 : end do
1019 : end if
1020 : end do ! iab
1021 6756 : call timab(443,2,tsec) ! ac_lrk_appl
1022 : end if
1023 :
1024 1513957 : do kb=ib1,ib2
1025 1281718 : call timab(438,1,tsec) ! (2)
1026 :
1027 : ! Get frequencies $\omega$-\epsilon_in$ to evaluate $d\Sigma/dE$, note the spin
1028 : ! subtract e_KS since we have stored e_KS+ Delta \omega in Sr%omega4sd, not required for AC
1029 8561820 : do io=Sr%nomega_r+1,nomega_tot
1030 8561820 : omegame0i(io)=DBLE(Sr%omega4sd(kb,jk_ibz,io-Sr%nomega_r,spin)) - e0i
1031 : end do
1032 :
1033 : ! Get the ket \Sigma|\phi_{k,kb}> according to the method.
1034 53721188 : rhotwgp(:)=rhotwg_ki(:,kb)
1035 :
1036 : select case (mod10)
1037 : case (SIG_GW_PPM)
1038 : ! GW WITH Plasmon-Pole Model.
1039 : ! Note that ppmodel 3 or 4 work only in case of standard perturbative approach!
1040 : ! Moreover, for ppmodel 3 and 4, spinorial case is not allowed
1041 507066502 : sigc_ket = czero_gw
1042 : call PPm%calc_sigc(nspinor, npwc, nomega_tot, rhotwgp, botsq, otq, &
1043 496082 : omegame0i, Sigp%zcut, theta_mu_minus_e0i, eig, npwc, sigc_ket, sigcme_3)
1044 :
1045 496082 : if (PPm%model==3 .or. PPm%model==4) then
1046 33280 : sigcme2(:,kb)=sigcme2(:,kb) + (wtqp+wtqm)*DBLE(sigcme_3(:)) + (wtqp-wtqm)*j_gw*AIMAG(sigcme_3(:))
1047 : end if
1048 :
1049 : case (SIG_GW_AC)
1050 : ! GW with Analytic continuation. This part is so optimized for AC that there is nothing to do here!
1051 :
1052 : case (SIG_GW_CD)
1053 : ! GW with contour deformation.
1054 : ! Check if pole contributions need to be summed. This avoids unnecessary splint calls and saves time.
1055 : !me_calc_poles = .TRUE.
1056 88622398 : sigc_ket = czero_gw
1057 2449180 : do io=1,nomega_tot
1058 2449180 : if (omegame0i(io)>=zero.AND.(ABS(one-theta_mu_minus_e0i)>zero)) then
1059 : !me_calc_poles(io) = .TRUE.
1060 183948 : if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io))
1061 2058714 : else if (omegame0i(io)<zero.AND.(ABS(theta_mu_minus_e0i)>zero)) then
1062 : !me_calc_poles(io) = .TRUE.
1063 447917 : if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io))
1064 : end if
1065 : end do
1066 :
1067 : ! Check memory saving
1068 206518 : if (epsm1%mqmem == 0) then
1069 : call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,epsm1%nomega,epsm1%nomega_r,epsm1%nomega_i,rhotwgp,&
1070 : epsm1%omega,epsm1%epsm1(:,:,:,1),omegame0i,theta_mu_minus_e0i,sigc_ket,Dtset%ppmfrq,npoles_missing(kb),&
1071 0 : method=Dtset%cd_frqim_method)
1072 : else
1073 : call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,epsm1%nomega,epsm1%nomega_r,epsm1%nomega_i,rhotwgp,&
1074 : epsm1%omega,epsm1%epsm1_qbz,omegame0i,theta_mu_minus_e0i,sigc_ket,Dtset%ppmfrq,npoles_missing(kb),&
1075 206518 : method=Dtset%cd_frqim_method)
1076 : end if
1077 :
1078 : #if 0
1079 : if (wtqm/=0) then
1080 : call calc_sigc_cd(npwc,npwc,nspinor,,nomega_tot,epsm1%nomega,epsm1%nomega_r,epsm1%nomega_i,rhotwgp,&
1081 : epsm1%omega,epsm1_trcc_qbz,omegame0i,theta_mu_minus_e0i,aherm_sigc_ket,Dtset%ppmfrq,npoles_missing(kb),&
1082 : method=Dtset%cd_frqim_method)
1083 :
1084 : herm_sigc_ket = half*(sigc_ket + aherm_sigc_ket)
1085 : aherm_sigc_ket = half*(sigc_ket - aherm_sigc_ket)
1086 : else
1087 : herm_sigc_ket = sigc_ket
1088 : aherm_sigc_ket = czero_gw
1089 : end if
1090 : #endif
1091 :
1092 : case (SIG_QPGW_PPM)
1093 : ! MODEL GW calculation WITH PPm TODO Spinor not tested.
1094 : ! Calculate \Sigma(E_k) |k> to obtain <j|\Sigma(E_k)|k>
1095 1426032 : ABI_MALLOC(sigcme_new, (nomega_tot))
1096 40748304 : sigc_ket = czero_gw; ket1 = czero_gw; ket2 = czero_gw
1097 :
1098 : call PPm%calc_sigc(nspinor, npwc, nomega_tot, rhotwgp, botsq,otq, &
1099 475344 : omegame0i, Sigp%zcut, theta_mu_minus_e0i, eig, npwc, ket1, sigcme_new)
1100 :
1101 475344 : if (Sigp%gwcalctyp==28) then
1102 474984 : if (PPm%model/=1.and.PPm%model/=2) then
1103 : ! This is needed to have npwc=PPm%dm2_botsq=PPm%dm2_otq
1104 : write(msg,'(3a)')&
1105 0 : 'For the time being, gwcalctyp=28 cannot be used with ppmodel=3,4.',ch10,&
1106 0 : 'Use another Plasmon Pole Model when gwcalctyp=28.'
1107 0 : ABI_ERROR(msg)
1108 : end if
1109 1899936 : ABI_MALLOC(botsq_conjg_transp,(PPm%dm2_botsq,npwc))
1110 370015728 : botsq_conjg_transp=TRANSPOSE(botsq) ! Keep these two lines separated, otherwise gfortran messes up
1111 369540744 : botsq_conjg_transp=CONJG(botsq_conjg_transp)
1112 1899936 : ABI_MALLOC(otq_transp,(PPm%dm2_otq,PPm%npwc))
1113 370015728 : otq_transp=TRANSPOSE(otq)
1114 :
1115 : call PPm%calc_sigc(nspinor, npwc, nomega_tot, rhotwgp, botsq_conjg_transp, otq_transp, &
1116 474984 : omegame0i, Sigp%zcut, theta_mu_minus_e0i, eig, npwc, ket2, sigcme_3)
1117 :
1118 474984 : ABI_FREE(botsq_conjg_transp)
1119 474984 : ABI_FREE(otq_transp)
1120 14364208 : sigc_ket= half*(ket1+ket2)
1121 : else
1122 10800 : sigc_ket= ket1
1123 : end if
1124 :
1125 475344 : ABI_FREE(sigcme_new)
1126 :
1127 : case (SIG_QPGW_CD)
1128 : ! MODEL GW with numerical integration.
1129 : ! Check if pole contributions need to be summed. this avoids unnecessary splint calls and saves time
1130 : !me_calc_poles = .TRUE.
1131 2709504 : sigc_ket = czero_gw; ket1 = czero_gw; ket2 = czero_gw
1132 :
1133 110592 : do io=1,nomega_tot
1134 110592 : if (omegame0i(io)>=zero.AND.(ABS(one-theta_mu_minus_e0i)>zero)) then
1135 : !me_calc_poles(io) = .TRUE.
1136 3191 : if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io))
1137 52105 : else if (omegame0i(io)<zero.AND.(ABS(theta_mu_minus_e0i)>zero)) then
1138 : !me_calc_poles(io) = .TRUE.
1139 2820 : if ( w_maxval(kb) < ABS(omegame0i(io)) ) w_maxval(kb) = ABS(omegame0i(io))
1140 : end if
1141 : end do
1142 :
1143 : ! Calculate \Sigma(E_k)|k> to obtain <j|\Sigma(E_k)|k>
1144 : call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,epsm1%nomega,epsm1%nomega_r,epsm1%nomega_i,rhotwgp,&
1145 : epsm1%omega,epsm1%epsm1_qbz,omegame0i,theta_mu_minus_e0i,ket1,Dtset%ppmfrq,npoles_missing(kb),&
1146 55296 : method=Dtset%cd_frqim_method)
1147 :
1148 55296 : if (Sigp%gwcalctyp==29) then
1149 : ! Calculate \Sigma^*(E_k)|k> to obtain <k|\Sigma(E_k)|j>^*
1150 : call calc_sigc_cd(npwc,npwc,nspinor,nomega_tot,epsm1%nomega,epsm1%nomega_r,epsm1%nomega_i,rhotwgp,&
1151 : epsm1%omega,epsm1_trcc_qbz,omegame0i,theta_mu_minus_e0i,ket2,Dtset%ppmfrq,npoles_missing(kb),&
1152 55296 : method=Dtset%cd_frqim_method)
1153 995328 : sigc_ket = half*(ket1+ket2)
1154 : else
1155 0 : sigc_ket = ket1
1156 : end if
1157 :
1158 : case default
1159 1281718 : ABI_ERROR(sjoin("Unsupported value for mod10:", itoa(mod10)))
1160 : end select
1161 :
1162 1281718 : if (Sigp%gwcomp==1) then
1163 : ! TODO spinor not implemented
1164 : call calc_sig_ppm_comp(npwc,nomega_tot,rhotwgp,botsq,otq,DBLE(Sr%egw(kb,jk_ibz,spin)-en_high),&
1165 7620 : Sigp%zcut,theta_mu_minus_e0i,sigc_ket,PPm%model,npwc,PPm%dm2_botsq,PPm%dm2_otq)
1166 : end if
1167 :
1168 1281718 : if (Sigp%gwcomp == 2) then
1169 : !wc0_qbz => epsm1%epsm1_qbz(:,:,1)
1170 0 : coh_g = matmul(wc0_qbz, rhotwg_ki(:,kb))
1171 : end if
1172 :
1173 1281718 : call timab(438,2,tsec) !
1174 1281718 : call timab(439,1,tsec) ! sigma_me
1175 :
1176 : ! Loop over the non-zero row elements of this column.
1177 : ! 1) If gwcalctyp<20 : only diagonal elements since QP==KS.
1178 : ! 2) If gwcalctyp>=20: only off-diagonal elements connecting states with same character.
1179 7984404 : do irow=1,Sigcij_tab(spin)%col(kb)%size1
1180 6702686 : jb = Sigcij_tab(spin)%col(kb)%bidx(irow)
1181 192780578 : rhotwg = rhotwg_ki(:,jb)
1182 :
1183 : ! Calculate <\phi_j|\Sigma_c|\phi_k>
1184 : ! Different freqs according to method (AC or Perturbative), see nomega_sigc.
1185 6702686 : if (mod10==SIG_GW_AC) then
1186 9297064 : sigctmp(:,:) = czero_gw
1187 585404 : do iab=1,Sigp%nsig_ab
1188 9297064 : do io=1,nomega_sigc
1189 8711660 : omegame0i_ac = Sr%omega_i(io) - qp_ene(ib_sum,ik_ibz,spin)
1190 9004362 : if (epsm1%nomega_i_conv > 0) then
1191 28800 : omegame0i2_ac = omegame0i_ac*omegame0i_ac
1192 892800 : do iiw=1,epsm1%nomega_i_conv
1193 : sigctmp(io,iab) = sigctmp(io,iab) + &
1194 : piinv * ((wtqp+wtqm)*DBLE(conv_rhotw_epsm1_rhotw(jb,kb,iiw,iab)) + (wtqp-wtqm)*j_gw*AIMAG(conv_rhotw_epsm1_rhotw(jb,kb,iiw,iab))) * &
1195 892800 : omegame0i_ac / (omegame0i2_ac + conv_omegap2(iiw)) * conv_gl_wts(iiw) / conv_gl_knots(iiw)**2
1196 : end do
1197 : else
1198 8682860 : select case (epsm1%hscr%iw_mesh_type)
1199 : case ("gauss_legendre")
1200 8682860 : omegame0i2_ac = omegame0i_ac*omegame0i_ac
1201 826026460 : do iiw=1,epsm1%nomega_i
1202 : sigctmp(io,iab) = sigctmp(io,iab) + &
1203 : piinv * ((wtqp+wtqm)*DBLE(rhotw_epsm1_rhotw(jb,kb,iiw,iab)) + (wtqp-wtqm)*j_gw*AIMAG(rhotw_epsm1_rhotw(jb,kb,iiw,iab))) * &
1204 826026460 : omegame0i_ac / (omegame0i2_ac + omegap2(iiw)) * gl_wts(iiw) / gl_knots(iiw)**2
1205 : end do
1206 : case ("minimax")
1207 : ! NB: Sigma_c along the imag. axis has a -1/2pi factor.
1208 : ! Here the -1 factor disappears because we have performed an EIGEN decomposition of -(epsm1-1).
1209 0 : do iiw=1,epsm1%nomega_i
1210 : sigctmp(io,iab) = sigctmp(io,iab) + &
1211 : (piinv / two) * ((wtqp+wtqm)*DBLE(rhotw_epsm1_rhotw(jb,kb,iiw,iab)) + (wtqp-wtqm)*j_gw*AIMAG(rhotw_epsm1_rhotw(jb,kb,iiw,iab))) * &
1212 : ((one / (omegame0i_ac + omegap_cplx(iiw))) + (one / (omegame0i_ac - omegap_cplx(iiw)))) * &
1213 0 : epsm1%hscr%omega_wgs(epsm1%nomega_r+iiw)
1214 : end do
1215 : case default
1216 8682860 : ABI_ERROR(sjoin("Invalid iw_mesh_type:", epsm1%hscr%iw_mesh_type))
1217 : end select
1218 : end if
1219 : end do
1220 : end do
1221 :
1222 : else
1223 : ! All other cases (not AC)
1224 12863552 : do iab=1,Sigp%nsig_ab
1225 6453568 : spadc1 = spinor_padc(1, iab); spadc2 = spinor_padc(2, iab)
1226 52467536 : do io=1,nomega_sigc
1227 46057552 : sigctmp(io,iab) = XDOTC(npwc,rhotwg(spadc1+1:),1,sigc_ket(spadc2+1:,io),1)
1228 : end do
1229 : end do
1230 : end if
1231 :
1232 6702686 : if (Sigp%gwcomp==1) then
1233 : ! Evaluate Extrapolar term TODO this does not work with spinor
1234 10860 : if (extrapolar_distrb(jb,kb,ik_bz,spin) == Wfd%my_rank) then
1235 : ! Do it once as it does not depend on the ib_sum index being summed over.
1236 690 : extrapolar_distrb(jb,kb,ik_bz,spin) = xmpi_undefined_rank
1237 : call calc_wfwfg(ktabr(:,jk_ibz),jik, spinrot_kgw, & ! TODO: why jk_ibz?
1238 690 : gwc_nfftot,nspinor,gwc_ngfft,wfr_bdgw(:,jb),wfr_bdgw(:,kb),wf1swf2_g)
1239 :
1240 690 : if (Psps%usepaw==1) then
1241 0 : i1=jb; i2=kb
1242 0 : if (nspinor==2) then
1243 0 : i1=(2*jb-1); i2=(2*kb-1)
1244 : end if
1245 0 : spad=(nspinor-1)
1246 : call paw_rho_tw_g(cryst, Pwij_fft,gwc_nfftot,Sigp%nsig_ab,nspinor, &
1247 0 : gw_gfft,Cprj_kgw(:,i1:i1+spad),Cprj_kgw(:,i2:i2+spad),wf1swf2_g)
1248 :
1249 0 : if (Dtset%pawcross==1) then ! Add paw cross term
1250 : call paw_cross_rho_tw_g(nspinor,npwc,nfftf,rho_ngfft,1,use_padfftf,igfftfcg0,gboundf,&
1251 : ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
1252 : ur_ae_bdgw(:,kb),ur_ae_onsite_bdgw(:,kb),ur_ps_onsite_bdgw(:,kb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
1253 0 : nspinor,wf1swf2_g)
1254 : end if
1255 : end if
1256 :
1257 : ! The static contribution from completeness relation is calculated once.
1258 : call calc_coh_comp(iq_ibz,Vcp%i_sz,(jb==kb),nspinor,Sigp%nsig_ab,DBLE(Sr%egw(kb,jk_ibz,spin)-en_high),&
1259 690 : npwc,Gsph_c%gvec,gwc_ngfft,gwc_nfftot,wf1swf2_g,vc_sqrt_qbz,botsq,otq,sigcohme)
1260 :
1261 4260 : do io=1,nomega_sigc
1262 7830 : sigctmp(io,:) = sigctmp(io,:)+sigcohme(:)
1263 : end do
1264 : end if
1265 : end if ! gwcomp==1
1266 :
1267 6702686 : if (Sigp%gwcomp == 2) then
1268 : !sigcohme(1) = quarter * dot_product(conjg(rhotwg), coh_g)
1269 0 : sigcohme(1) = quarter * dot_product(rhotwg, coh_g)
1270 0 : do io=1,nomega_sigc
1271 0 : sigctmp(io,:) = sigctmp(io,:) + sigcohme(1)
1272 : end do
1273 : end if
1274 :
1275 : ! Accumulate and, in case, symmetrize matrix elements of Sigma_c
1276 14730674 : do iab=1,Sigp%nsig_ab
1277 6746270 : is_idx=spin; if (nspinor==2) is_idx=iab
1278 :
1279 6702686 : select case (mod10)
1280 : case (SIG_GW_AC)
1281 9004362 : sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + sigctmp(:,iab)
1282 9004362 : sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + sigctmp(:,iab)
1283 9004362 : sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + czero
1284 : case default
1285 : sigcme_tmp(:,jb,kb,is_idx)=sigcme_tmp(:,jb,kb,is_idx) + &
1286 46057552 : (wtqp+wtqm)*DBLE(sigctmp(:,iab)) + (wtqp-wtqm)*j_gw*AIMAG(sigctmp(:,iab))
1287 46057552 : sigc(1,:,jb,kb,is_idx)=sigc(1,:,jb,kb,is_idx) + wtqp* sigctmp(:,iab)
1288 52803822 : sigc(2,:,jb,kb,is_idx)=sigc(2,:,jb,kb,is_idx) + wtqm*CONJG(sigctmp(:,iab))
1289 : end select
1290 : ! TODO this should be the contribution coming from the anti-hermitian part.
1291 : end do
1292 : end do ! irow used to calculate matrix elements of $\Sigma$
1293 :
1294 : ! shaltaf (030406): this has to be done in a clean way later.
1295 : ! TODO does not work with spinor.
1296 1281718 : if (mod10==SIG_GW_PPM.and.(PPm%model==3.or.PPm%model==4)) then
1297 33280 : sigcme_tmp(:,kb,kb,spin)= sigcme2(:,kb)
1298 33280 : sigc(1,:,kb,kb,spin)= sigcme2(:,kb)
1299 33280 : sigc(2,:,kb,kb,spin)= czero
1300 : end if
1301 :
1302 1558966 : call timab(439,2,tsec) ! csigme(SigC)
1303 : end do !kb to calculate matrix elements of $\Sigma$
1304 : end do !ib_sum
1305 :
1306 19839 : call timab(445,2,tsec) ! csigme(SigC)
1307 19839 : if (print_time) then
1308 3658 : write(msg,'(3(a,i0))')' Sigma_c: ik_bz: ',ik_bz,'/',Kmesh%nbz,", spin: ",spin
1309 3658 : call cwtime_report(msg, cpu_k, wall_k, gflops_k); if (ik_bz == LOG_MODK) call wrtout(std_out, " ...")
1310 : end if
1311 :
1312 : ! Deallocate k-dependent quantities.
1313 19839 : ABI_FREE(gw_gbound)
1314 19839 : if (Dtset%pawcross==1) then
1315 0 : ABI_FREE(gboundf)
1316 : end if
1317 :
1318 19839 : if (sigp%needs_ppm()) then
1319 14122 : ABI_FREE(botsq)
1320 14122 : ABI_FREE(otq)
1321 14122 : ABI_FREE(eig)
1322 : end if
1323 20300 : if (Psps%usepaw==1) then
1324 192 : call pawpwij_free(Pwij_qg)
1325 390 : ABI_FREE(Pwij_qg)
1326 : end if
1327 :
1328 : end do ! ik_bz
1329 :
1330 461 : ABI_FREE(wfr_bdgw)
1331 916 : if (Wfd%usepaw==1) then
1332 9 : call pawcprj_free(Cprj_kgw)
1333 143 : ABI_FREE(Cprj_kgw)
1334 9 : if (Dtset%pawcross==1) then
1335 0 : ABI_FREE(ur_ae_bdgw)
1336 0 : ABI_FREE(ur_ae_onsite_bdgw)
1337 0 : ABI_FREE(ur_ps_onsite_bdgw)
1338 : end if
1339 : end if
1340 : end do ! spin
1341 :
1342 : ! If epsm1 is MPI-shared, we have to close the RMA epoch.
1343 455 : if (epsm1%use_mpi_shared_win) call xmpi_win_fence(XMPI_MODE_NOSUCCEED, epsm1%epsm1_win, ierr)
1344 455 : call epsm1%free_epsm1_qbz()
1345 :
1346 455 : ABI_FREE(sigcme2)
1347 455 : ABI_FREE(sigcme_3)
1348 455 : ABI_FREE(igfftcg0)
1349 455 : if (Dtset%pawcross==1) then
1350 0 : ABI_FREE(igfftfcg0)
1351 : end if
1352 :
1353 : ! Gather contributions from all the CPUs
1354 455 : call timab(440,1,tsec) ! wfd_barrier
1355 455 : call timab(440,2,tsec) ! wfd_barrier
1356 455 : call timab(441,1,tsec) ! xmpi_sum
1357 :
1358 455 : call xmpi_sum(sigcme_tmp, Wfd%comm, ierr)
1359 455 : call xmpi_sum(sigc, Wfd%comm, ierr)
1360 455 : call timab(441,2,tsec) ! xmpi_sum
1361 :
1362 : ! Multiply by constants. In 3D systems sqrt(4pi) is included in vc_sqrt_qbz.
1363 247042 : sigcme_tmp = sigcme_tmp /(Cryst%ucvol*Kmesh%nbz)
1364 671490 : sigc = sigc /(Cryst%ucvol*Kmesh%nbz)
1365 :
1366 : ! If we have summed over the IBZ_q now we have to average over degenerate states
1367 : ! Presently only diagonal terms are considered
1368 : ! TODO QP-SCGW required a more involved approach, there is a check in sigma
1369 : ! TODO it does not work if nspinor==2.
1370 455 : call timab(442,1,tsec) ! final ops
1371 :
1372 916 : do spin=1,Wfd%nsppol
1373 916 : if (can_symmetrize(spin)) then
1374 109 : if (mod10==SIG_GW_AC) then ! FIXME here there is a problem in case of AC with symmetries
1375 42 : ABI_MALLOC(sym_cme, (Sr%nomega_i, ib1:ib2, ib1:ib2, Sigp%nsig_ab))
1376 : else
1377 612 : ABI_MALLOC(sym_cme, (nomega_tot, ib1:ib2, ib1:ib2, Sigp%nsig_ab))
1378 : end if
1379 69289 : sym_cme=czero
1380 :
1381 : ! Average over degenerate diagonal elements
1382 : ! NOTE: frequencies for \Sigma_c(\omega) should be equal to avoid spurious results.
1383 : ! another good reason to use a strict criterion for the tolerance on eigenvalues.
1384 778 : do ib=ib1,ib2
1385 669 : ndegs=0
1386 5784 : do jb=ib1,ib2
1387 5115 : if (degtab(ib,jb,spin)==1) then
1388 1553 : if (nspinor == 1) then
1389 39755 : sym_cme(:, ib, ib, 1) = sym_cme(:, ib, ib, 1) + SUM(sigc(:,:,jb,jb,spin), DIM=1)
1390 : else
1391 840 : do ii=1,Sigp%nsig_ab
1392 18984 : sym_cme(:, ib, ib, ii) = sym_cme(:, ib, ib, ii) + SUM(sigc(:,:,jb,jb,ii), dim=1)
1393 : end do
1394 : end if
1395 : end if
1396 5784 : ndegs = ndegs + degtab(ib,jb,spin)
1397 : end do
1398 8801 : sym_cme(:,ib,ib,:) = sym_cme(:,ib,ib,:) / ndegs
1399 : end do
1400 :
1401 109 : if (Sigp%gwcalctyp >= 20) then
1402 0 : do iwc=1,nomega_sigc
1403 0 : call QP_sym(spin)%symmetrize_mels(ib1,ib2,sigc(:,iwc,:,:,spin),sym_cme(iwc,:,:,1))
1404 : end do
1405 : end if
1406 :
1407 : ! Copy symmetrized values
1408 778 : do ib=ib1,ib2
1409 5893 : do jb=ib1,ib2
1410 : !if (mod10==SIG_GW_AC.and.average_real) CYCLE ! this is to check another scheme in case of AC
1411 5784 : if (nspinor == 1) then
1412 46833 : sigcme_tmp(:,ib,jb,spin) = sym_cme(:,ib,jb,1)
1413 : else
1414 2680 : do ii=1,Sigp%nsig_ab
1415 21976 : sigcme_tmp(:,ib,jb,ii) = sym_cme(:,ib,jb,ii)
1416 : end do
1417 : end if
1418 : end do
1419 : end do
1420 109 : ABI_FREE(sym_cme)
1421 : end if
1422 : end do
1423 :
1424 : ! Reconstruct the full sigma matrix from the upper triangle (only for HF, SEX and COHSEX)
1425 : !if (Sigp%gwcalctyp>=20 .and. sigma_is_herm(Sigp) ) then
1426 : ! ABI_CHECK(nspinor==1,"cannot hermitianize non-collinear sigma!")
1427 : ! do spin=1,Wfd%nsppol
1428 : ! do io=1,nomega_sigc
1429 : ! call hermitianize(sigcme_tmp(io,:,:,spin),"Upper")
1430 : ! end do
1431 : ! end do
1432 : !end if
1433 :
1434 : ! GW with contour deformation: check on the number of poles not included.
1435 455 : if (ANY(mod10 == [SIG_GW_CD, SIG_QPGW_CD])) then
1436 99 : call xmpi_sum(npoles_missing, Wfd%comm, ierr)
1437 828 : npls = SUM(npoles_missing)
1438 99 : if (npls>0) then
1439 11 : ABI_WARNING(sjoin("Total number of missing poles for contour deformation method:", itoa(npls)))
1440 36 : do band=minbnd,maxbnd
1441 25 : npls = npoles_missing(band)
1442 36 : if (npls > 0) then
1443 24 : write(msg,'(a,2(i0,a))')" For band ",band," there are ",npls," missing poles"
1444 24 : call wrtout(std_out, msg)
1445 : end if
1446 : end do
1447 : end if
1448 : ! Print data on the maximum value needed for the screening along the real axis
1449 927 : w_localmax = MAXVAL(w_maxval)
1450 99 : call xmpi_max(w_localmax,w_max, Wfd%comm, ierr)
1451 99 : write(msg,'(a,f12.5,a)') ' Max omega value used in W(omega): ',w_max*Ha_eV,' [eV]'
1452 99 : call wrtout(std_out, msg)
1453 : end if
1454 455 : call timab(442,2,tsec) ! final ops
1455 :
1456 : ! ===========================
1457 : ! ==== Deallocate memory ====
1458 : ! ===========================
1459 455 : if (Psps%usepaw==1) then
1460 9 : ABI_SFREE(gw_gfft)
1461 9 : call pawcprj_free(Cprj_ksum)
1462 48 : ABI_FREE(Cprj_ksum)
1463 9 : if (allocated(Pwij_fft)) then
1464 0 : call pawpwij_free(Pwij_fft)
1465 0 : ABI_FREE(Pwij_fft)
1466 : end if
1467 9 : if (Dtset%pawcross==1) then
1468 0 : ABI_FREE(ur_ae_sum)
1469 0 : ABI_FREE(ur_ae_onsite_sum)
1470 0 : ABI_FREE(ur_ps_onsite_sum)
1471 0 : ABI_FREE(ktabrf)
1472 : end if
1473 : end if
1474 :
1475 455 : ABI_SFREE(npoles_missing)
1476 455 : ABI_SFREE(ur_ibz)
1477 455 : ABI_SFREE(usr_bz)
1478 455 : ABI_SFREE(ktabr)
1479 455 : ABI_SFREE(rhotwg_ki)
1480 455 : ABI_SFREE(rhotwg)
1481 455 : ABI_SFREE(rhotwgp)
1482 455 : ABI_SFREE(vc_sqrt_qbz)
1483 455 : ABI_SFREE(omegame0i)
1484 455 : ABI_SFREE(sigctmp)
1485 455 : ABI_SFREE(sigc)
1486 455 : ABI_SFREE(w_maxval)
1487 455 : ABI_SFREE(sigc_ket)
1488 455 : ABI_SFREE(ket1)
1489 455 : ABI_SFREE(ket2)
1490 455 : ABI_SFREE(epsm1_trcc_qbz)
1491 455 : ABI_SFREE(epsm1_tmp)
1492 455 : ABI_SFREE(degtab)
1493 455 : ABI_SFREE(rhotw_epsm1_rhotw)
1494 455 : ABI_SFREE(conv_rhotw_epsm1_rhotw)
1495 : ABI_SFREE(aherm_sigc_ket)
1496 : ABI_SFREE(herm_sigc_ket)
1497 455 : ABI_SFREE(wf1swf2_g)
1498 455 : ABI_SFREE(extrapolar_distrb)
1499 455 : ABI_SFREE(proc_distrb)
1500 :
1501 455 : if (mod10 == SIG_GW_AC) then
1502 28 : if (epsm1%use_mpi_shared_win) then
1503 0 : call xmpi_win_free(ac_epsm1cqwz2_win, ierr)
1504 : else
1505 28 : ABI_SFREE_PTR(ac_epsm1cqwz2)
1506 : end if
1507 : end if
1508 :
1509 455 : if (dtset%gwcomp == 2) then
1510 0 : ABI_SFREE(wc0_qbz)
1511 0 : do iq_ibz=1,qmesh%nibz
1512 0 : call wc0_qibz(iq_ibz)%free()
1513 : end do
1514 0 : ABI_FREE(wc0_qibz)
1515 0 : ABI_FREE(coh_g)
1516 : end if
1517 :
1518 455 : call timab(431,2,tsec)
1519 455 : call timab(424,2,tsec) ! calc_sigc_me
1520 :
1521 455 : call cwtime_report(" calc_sigc_me", cpu_all, wall_all, gflops_all)
1522 :
1523 : DBG_EXIT("COLL")
1524 :
1525 2730 : end subroutine calc_sigc_me
1526 : !!***
1527 :
1528 : !----------------------------------------------------------------------
1529 :
1530 : !!****f* ABINIT/calc_coh_comp
1531 : !! NAME
1532 : !! calc_coh_comp
1533 : !!
1534 : !! FUNCTION
1535 : !! Calculates the COH-like contribution to the self-energy when
1536 : !! the extrapolar technique and the closure relation is used to
1537 : !! reduce the number of empty states to be summed over in the Green
1538 : !! function entering the definition of the GW self-energy.
1539 : !!
1540 : !! INPUTS
1541 : !! iqibz=index of the irreducible q-point in the array qibz, point which is
1542 : !! related by a symmetry operation to the point q summed over (see csigme).
1543 : !! This index is also used to treat the integrable coulombian singularity at q=0
1544 : !! ngfft(18)=contain all needed information about 3D FFT for GW wavefuntions,
1545 : !! see ~abinit/doc/variables/vargs.htm#ngfft
1546 : !! nsig_ab=Number of components in the self-energy operator (1 for collinear magnetism)
1547 : !! npwc=number of plane waves in $\tilde epsilon^{-1}$
1548 : !! nspinor=Number of spinorial components.
1549 : !! i_sz=contribution arising from the integrable coulomb singularity at q==0
1550 : !! (see csigme for the method used), note that in case of 3-D systems the factor
1551 : !! 4pi in the Coulomb potential is included in the definition of i_sz
1552 : !! gvec(3,npwc)=G vectors in reduced coordinates
1553 : !! vc_sqrt(npwc)= square root of the coulombian matrix elements for this q-point
1554 : !! botsq = Plasmon-pole parameters
1555 : !! otq = PPm parameters
1556 : !!
1557 : !! OUTPUT
1558 : !! sigcohme=partial contribution to the matrix element of $<jb k|\Sigma_{COH}| kb k>$
1559 : !! coming from this single q-point for completeness trick
1560 : !!
1561 : !! SOURCE
1562 :
1563 690 : subroutine calc_coh_comp(iqibz,i_sz,same_band,nspinor,nsig_ab,ediff,npwc,gvec,&
1564 690 : ngfft,nfftot,wfg2_jk,vc_sqrt,botsq,otq,sigcohme)
1565 :
1566 : !Arguments ------------------------------------
1567 : !scalars
1568 : integer,intent(in) :: iqibz,npwc,nsig_ab,nspinor,nfftot
1569 : real(dp),intent(in) :: i_sz,ediff
1570 : logical,intent(in) :: same_band
1571 : !arrays
1572 : integer,intent(in) :: gvec(3,npwc),ngfft(18)
1573 : complex(gwp),intent(in) :: botsq(npwc,npwc),otq(npwc,npwc)
1574 : complex(gwp),intent(in) :: vc_sqrt(npwc)
1575 : complex(gwp),intent(in) :: wfg2_jk(nfftot*nsig_ab)
1576 : complex(gwp),intent(out) :: sigcohme(nsig_ab)
1577 :
1578 : !Local variables-------------------------------
1579 : !scalars
1580 : integer,save :: enough=0
1581 : integer :: ig,ig4,ig4x,ig4y,ig4z,igp,igmin,ispinor,ngfft1,ngfft2,ngfft3,spad,outofbox
1582 : !arrays
1583 : integer :: g2mg1(3)
1584 : ! *************************************************************************
1585 :
1586 : DBG_ENTER("COLL")
1587 :
1588 : ! === Treat the case q --> 0 adequately ===
1589 : ! TODO Better treatment of wings
1590 690 : igmin=1 ; if (iqibz==1) igmin=2
1591 : !
1592 : ! === Partial contribution to the matrix element of Sigma_c ===
1593 : ! * For nspinor==2, the closure relation reads:
1594 : ! $\sum_s \psi_a^*(1)\psi_b(2) = \delta_{ab} \delta(1-2)$
1595 : ! where a,b are the spinor components. As a consequence, Sigma_{COH} is always
1596 : ! diagonal in spin-space and only diagonal matrix elements have to be calculated.
1597 : ! MG TODO wfg2_jk should be calculated on an augmented FFT box to avoid spurious wrapping of G1-G2.
1598 : !
1599 690 : ngfft1 = ngfft(1); ngfft2 = ngfft(2); ngfft3 = ngfft(3)
1600 1380 : sigcohme(:) = czero_gw
1601 :
1602 1380 : do ispinor=1,nspinor
1603 690 : spad=(ispinor-1)*nfftot
1604 690 : outofbox=0
1605 :
1606 30968 : do igp=igmin,npwc
1607 1649010 : do ig=igmin,npwc
1608 :
1609 6472168 : g2mg1 = gvec(:,igp)-gvec(:,ig)
1610 12944336 : if (ANY(g2mg1(:)>ngfft(1:3)/2) .or. ANY(g2mg1(:)<-(ngfft(1:3)-1)/2)) then
1611 0 : outofbox = outofbox+1; CYCLE
1612 : end if
1613 :
1614 1618042 : ig4x=MODULO(g2mg1(1),ngfft1)
1615 1618042 : ig4y=MODULO(g2mg1(2),ngfft2)
1616 1618042 : ig4z=MODULO(g2mg1(3),ngfft3)
1617 1618042 : ig4= 1+ig4x+ig4y*ngfft1+ig4z*ngfft1*ngfft2
1618 :
1619 : !MG where is neta here, ediff, otq might be close to zero depending on gwecomp
1620 : sigcohme(ispinor) = sigcohme(ispinor) + &
1621 1648320 : & half*wfg2_jk(spad+ig4)*vc_sqrt(ig)*vc_sqrt(igp) * botsq(ig,igp) / ( otq(ig,igp) * ( ediff -otq(ig,igp) ) )
1622 : end do
1623 : end do
1624 :
1625 1380 : if (iqibz==1.and.same_band) then
1626 82 : sigcohme(ispinor) = sigcohme(ispinor) + half*wfg2_jk(spad+ig4)*i_sz*botsq(1,1) / ( otq(1,1) * (ediff -otq(1,1)) )
1627 : end if
1628 : end do !ispinor
1629 :
1630 690 : if (outofbox/=0) then
1631 0 : enough=enough+1
1632 0 : if (enough<=50) then
1633 0 : ABI_WARNING(sjoin('Number of G1-G2 pairs outside the G-sphere for Wfns: ',itoa(outofbox)))
1634 0 : if (enough==50) call wrtout(std_out,' ========== Stop writing Warnings ==========')
1635 : end if
1636 : end if
1637 :
1638 : DBG_EXIT("COLL")
1639 :
1640 690 : end subroutine calc_coh_comp
1641 : !!***
1642 :
1643 : !!****f* ABINIT/calc_sigc_cd
1644 : !! NAME
1645 : !! calc_sigc_cd
1646 : !!
1647 : !! FUNCTION
1648 : !! Calculate contributions to the self-energy operator with the contour deformation method.
1649 : !!
1650 : !! INPUTS
1651 : !! nomega=Total number of frequencies where $\Sigma_c$ matrix elements are evaluated.
1652 : !! nomegae=Number of frequencies where $\epsilon^{-1}$ has been evaluated.
1653 : !! nomegaei=Number of imaginary frequencies for $\epsilon^{-1}$ (non zero).
1654 : !! nomegaer=Number of real frequencies for $\epsilon^{-1}$
1655 : !! npwc=Number of G vectors for the correlation part.
1656 : !! npwx=Number of G vectors in rhotwgp for each spinorial component.
1657 : !! nspinor=Number of spinorial components.
1658 : !! theta_mu_minus_e0i=1 if e0i is occupied, 0 otherwise. Fractional occupancy in case of metals.
1659 : !! omegame0i(nomega)=Contains $\omega-\epsilon_{k-q,b1,\sigma}$
1660 : !! epsm1q(npwc,npwc,nomegae)=Symmetrized inverse dielectric matrix (exchange part is subtracted).
1661 : !! omega(nomegae)=Set of frequencies for $\epsilon^{-1}$.
1662 : !! rhotwgp(npwx*nspinor)=Matrix elements: $<k-q,b1,\sigma|e^{-i(q+G)r} |k,b2,\sigma>*vc_sqrt$
1663 : !!
1664 : !! OUTPUT
1665 : !! ket(npwc,nomega)=Contains \Sigma_c(\omega)|\phi> in reciprocal space.
1666 : !!
1667 : !! SIDE EFFECTS
1668 : !! npoles_missing=Incremented with the number of poles whose contribution has not been taken into account due to
1669 : !! limited frequency mesh used for W.
1670 : !!
1671 : !! SOURCE
1672 :
1673 317110 : subroutine calc_sigc_cd(npwc,npwx,nspinor,nomega,nomegae,nomegaer,nomegaei,rhotwgp,&
1674 317110 : omega,epsm1q,omegame0i,theta_mu_minus_e0i,ket,plasmafreq,npoles_missing,&
1675 0 : calc_poles, method) ! optional
1676 :
1677 : !Arguments ------------------------------------
1678 : !scalars
1679 : integer,intent(in) :: nomega,nomegae,nomegaei,nomegaer,npwc,npwx,nspinor
1680 : integer,intent(inout) :: npoles_missing
1681 : real(dp),intent(in) :: theta_mu_minus_e0i,plasmafreq
1682 : !arrays
1683 : real(dp),intent(in) :: omegame0i(nomega)
1684 : complex(dp),intent(in) :: omega(nomegae)
1685 : complex(gwp),intent(in) :: epsm1q(npwc,npwc,nomegae)
1686 : complex(gwp),intent(in) :: rhotwgp(npwx*nspinor)
1687 : complex(gwp),intent(inout) :: ket(nspinor*npwc,nomega)
1688 : logical, intent(in), optional :: calc_poles(nomega)
1689 : integer, intent(in), optional :: method
1690 :
1691 : !Local variables-------------------------------
1692 : !scalars
1693 : integer, parameter :: FABIEN=1,TRAPEZOID=2,NSPLINE=3
1694 : integer :: ii,ig,io,ios,ispinor,spadc,spadx,my_err,ierr,GK_LEVEL,INTMETHOD
1695 : integer :: i,j
1696 : real(dp) :: rt_imag,rt_real,local_one,local_zero
1697 : real(dp) :: intsign,temp1,temp2,temp3,temp4
1698 : real(dp) :: alpha,inv_alph,beta,alphsq,betasq,inv_beta
1699 : real(dp) :: re_intG,re_intK,im_intG,im_intK,GKttab,tau,ttil
1700 : real(dp) :: ref,imf,r,s,r2,s2
1701 : complex(dp) :: ct,domegaleft,domegaright
1702 : complex(gwp) :: fact
1703 : !arrays
1704 634220 : real(dp) :: omegame0i_tmp(nomega),tmp_x(2),tmp_y(2)
1705 634220 : real(dp) :: left(nomega),right(nomega)
1706 634220 : real(dp) :: tbeta(nomega),tinv_beta(nomega),tbetasq(nomega)
1707 951330 : real(dp) :: atermr(nomega),aterml(nomega),logup(nomega),logdown(nomega)
1708 634220 : real(dp) :: rtmp_r(nomegaer),rtmp_i(nomegaer)
1709 951330 : real(dp) :: ftab(nomegaei+2),ftab2(nomegaei+2),xtab(nomegaei+2),y(3,nomegaei+2)
1710 634220 : real(dp) :: work(nomegaei+2),work2(nomegaei+2),y2(3,nomegaei+2)
1711 634220 : complex(dp) :: omega_imag(nomegaei+1)
1712 634220 : complex(gwp) :: epsrho(npwc,nomegae),epsrho_imag(npwc,nomegaei+1)
1713 634220 : complex(gwp) :: tfone(npwc,nomegaei+1),tftwo(npwc,nomegaei+1)
1714 634220 : complex(gwp) :: weight(nomegaei+1,nomega)
1715 634220 : complex(gwp) :: weight2(nomegaei,nomega)
1716 634220 : logical :: my_calc_poles(nomega)
1717 317110 : real(dp), allocatable :: KronN(:),KronW(:),GaussW(:),fint(:),fint2(:)
1718 : !*************************************************************************
1719 :
1720 2670364 : my_calc_poles=.TRUE.; my_err=0
1721 :
1722 : ! Set integration method for imaginary axis
1723 317110 : INTMETHOD = FABIEN
1724 317110 : if (present(method)) then
1725 : if (method==1) INTMETHOD = FABIEN
1726 317110 : if (method==2) INTMETHOD = TRAPEZOID
1727 317110 : if (method>2) then
1728 2520 : INTMETHOD = NSPLINE
1729 2520 : if (method==3) then
1730 840 : GK_LEVEL = 15
1731 840 : ABI_MALLOC(KronN,(GK_LEVEL))
1732 840 : ABI_MALLOC(KronW,(GK_LEVEL))
1733 840 : ABI_MALLOC(GaussW,(GK_LEVEL-8))
1734 840 : ABI_MALLOC(fint,(GK_LEVEL))
1735 840 : ABI_MALLOC(fint2,(GK_LEVEL))
1736 31920 : KronN(:) = Kron15N(:); KronW(:) = Kron15W(:); GaussW(:) = Gau7W(:)
1737 1680 : else if (method==4) then
1738 840 : GK_LEVEL = 23
1739 840 : ABI_MALLOC(KronN,(GK_LEVEL))
1740 840 : ABI_MALLOC(KronW,(GK_LEVEL))
1741 840 : ABI_MALLOC(GaussW,(GK_LEVEL-12))
1742 840 : ABI_MALLOC(fint,(GK_LEVEL))
1743 840 : ABI_MALLOC(fint2,(GK_LEVEL))
1744 48720 : KronN(:) = Kron23N(:); KronW(:) = Kron23W(:); GaussW(:) = Gau11W(:)
1745 : else if (method>4) then
1746 840 : GK_LEVEL = 31
1747 840 : ABI_MALLOC(KronN,(GK_LEVEL))
1748 840 : ABI_MALLOC(KronW,(GK_LEVEL))
1749 840 : ABI_MALLOC(GaussW,(GK_LEVEL-16))
1750 840 : ABI_MALLOC(fint,(GK_LEVEL))
1751 840 : ABI_MALLOC(fint2,(GK_LEVEL))
1752 65520 : KronN(:) = Kron31N(:); KronW(:) = Kron31W(:); GaussW(:) = Gau15W(:)
1753 : end if
1754 : end if
1755 : end if
1756 :
1757 : ! Avoid divergences in $\omega - \omega_s$.
1758 2670364 : omegame0i_tmp(:)=omegame0i(:)
1759 2670364 : do ios=1,nomega
1760 2670364 : if (ABS(omegame0i_tmp(ios))<tol6) omegame0i_tmp(ios)=sign(tol6,omegame0i_tmp(ios))
1761 : end do
1762 :
1763 642028 : do ispinor=1,nspinor
1764 324918 : spadx=(ispinor-1)*npwx; spadc=(ispinor-1)*npwc
1765 :
1766 : ! Calculate $ \sum_{Gp} (\epsilon^{-1}_{G Gp}(\omega)-\delta_{G Gp}) \rhotwgp(Gp) $
1767 : !$omp parallel do
1768 4275278 : do io=1,nomegae
1769 4275278 : call XGEMV('N',npwc,npwc,cone_gw,epsm1q(:,:,io),npwc,rhotwgp(1+spadx:),1,czero_gw,epsrho(:,io),1)
1770 : end do
1771 :
1772 : ! Integrand along the imaginary axis.
1773 9382600 : epsrho_imag(:,1)=epsrho(:,1)
1774 39518366 : epsrho_imag(:,2:nomegaei+1)=epsrho(:,nomegaer+1:nomegae)
1775 :
1776 : ! Frequency mesh for integral along the imaginary axis.
1777 324918 : omega_imag(1)=omega(1)
1778 1672956 : omega_imag(2:nomegaei+1)=omega(nomegaer+1:nomegae)
1779 :
1780 : ! Original implementation -- saved here for reference during development
1781 : ! === Perform integration along the imaginary axis ===
1782 : !do io=1,nomegaei+1
1783 : ! if (io==1) then
1784 : ! domegaleft = omega_imag(io)
1785 : ! domegaright =(omega_imag(io+1)-omega_imag(io ))*half
1786 : ! else if (io==nomegaei+1) then
1787 : ! domegaleft =(omega_imag(io )-omega_imag(io-1))*half
1788 : ! domegaright =(omega_imag(io )-omega_imag(io-1))*half
1789 : ! else
1790 : ! domegaleft =(omega_imag(io )-omega_imag(io-1))*half
1791 : ! domegaright =(omega_imag(io+1)-omega_imag(io ))*half
1792 : ! end if
1793 : ! do ios=1,nomega
1794 : ! omg2 = -AIMAG(omega_imag(io)+domegaright)/REAL(omegame0i_tmp(ios))
1795 : ! omg1 = -AIMAG(omega_imag(io)-domegaleft )/REAL(omegame0i_tmp(ios))
1796 : ! fact = ATAN(omg2)-ATAN(omg1)
1797 : ! ket(spadc+1:spadc+npwc,ios)=ket(spadc+1:spadc+npwc,ios)+epsrho_imag(:,io)*fact
1798 : ! end do
1799 : !end do !io
1800 :
1801 : !ket(spadc+1:spadc+npwc,:)=ket(spadc+1:spadc+npwc,:)/pi
1802 : ! ---------------- end of original implementation -----------------------
1803 :
1804 321558 : select case (INTMETHOD)
1805 : case (FABIEN)
1806 : ! Hopefully more effective implementation MS 04.08.2011
1807 : ! Perform integration along imaginary axis using BLAS
1808 : ! First calculate first and last weights
1809 2714844 : weight(1,:) = ATAN(-half*AIMAG(omega_imag(2))/REAL(omegame0i_tmp(:)))
1810 321558 : domegaleft = (three*omega_imag(nomegaei+1)-omega_imag(nomegaei))
1811 321558 : domegaright = (omega_imag(nomegaei+1)+omega_imag(nomegaei))
1812 2714844 : right(:) = -AIMAG(omega_imag(nomegaei+1)-omega_imag(nomegaei))*REAL(omegame0i_tmp(:))
1813 : left(:) = quarter*AIMAG(domegaleft)*AIMAG(domegaright) &
1814 2714844 : +REAL(omegame0i_tmp(:))*REAL(omegame0i_tmp(:))
1815 2714844 : do ios=1,nomega
1816 2714844 : weight(nomegaei+1,ios) = ATAN(right(ios)/left(ios))
1817 : end do
1818 : ! Calculate the rest of the weights
1819 1321158 : do io=2,nomegaei
1820 999600 : domegaleft = (omega_imag(io )+omega_imag(io-1))
1821 999600 : domegaright = (omega_imag(io+1)+omega_imag(io ))
1822 8493792 : right(:) = -half*AIMAG(omega_imag(io+1)-omega_imag(io-1))*REAL(omegame0i_tmp(:))
1823 : left(:) = REAL(omegame0i_tmp(:))*REAL(omegame0i_tmp(:)) &
1824 8493792 : +quarter*AIMAG(domegaleft)*AIMAG(domegaright)
1825 8815350 : do ios=1,nomega
1826 8493792 : weight(io,ios) = ATAN(right(ios)/left(ios))
1827 : end do
1828 : end do
1829 :
1830 : ! Use BLAS call to perform matrix-matrix multiplication and accumulation
1831 321558 : fact = CMPLX(piinv,zero)
1832 :
1833 : call xgemm('N','N',npwc,nomega,nomegaei+1,fact,epsrho_imag,npwc,&
1834 8192022 : weight,nomegaei+1,cone_gw,ket(spadc+1:spadc+npwc,:),npwc)
1835 :
1836 : case (TRAPEZOID)
1837 : ! Trapezoidal rule Transform omega coordinates
1838 840 : alpha = plasmafreq
1839 840 : alphsq = alpha*alpha
1840 840 : inv_alph = one/alpha
1841 :
1842 8400 : xtab(1:nomegaei+1) = AIMAG(omega_imag(:))/(AIMAG(omega_imag(:)) + alpha)
1843 840 : xtab(nomegaei+2) = one
1844 :
1845 : ! Efficient trapezoidal rule with BLAS calls
1846 8400 : tbeta(:) = REAL(omegame0i_tmp(:))
1847 8400 : tbetasq(:) = tbeta(:)*tbeta(:)
1848 8400 : tinv_beta(:) = one/tbeta(:)
1849 :
1850 7560 : do io=1,nomegaei
1851 67200 : atermr(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(io+1)-tbetasq(:))
1852 67200 : aterml(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(io )-tbetasq(:))
1853 67200 : right(:) = ATAN((atermr(:)-aterml(:))/(one+atermr(:)*aterml(:)))
1854 : logup(:) = ABS(((alphsq+tbetasq(:))*xtab(io+1)-two*tbetasq(:)) &
1855 67200 : *xtab(io+1)+tbetasq(:))
1856 : logdown(:) = ABS(((alphsq+tbetasq(:))*xtab(io )-two*tbetasq(:)) &
1857 67200 : *xtab(io )+tbetasq(:))
1858 : ! Trapezoid integration weights
1859 : weight(io,:) = CMPLX(-(half*alpha*tbeta(:)*LOG(logup(:)/logdown(:)) + tbetasq(:) &
1860 67200 : *right(:))/(alphsq+tbetasq(:)),zero)
1861 67200 : weight2(io,:) = CMPLX(-right(:),zero)
1862 : ! Linear interpolation coefficients for each section (sum over ig)
1863 : tfone(:,io) = (epsrho_imag(:,io+1)-epsrho_imag(:,io)) &
1864 107520 : /(xtab(io+1)-xtab(io))
1865 108360 : tftwo(:,io) = epsrho_imag(:,io) - tfone(:,io)*xtab(io)
1866 : end do
1867 :
1868 : ! Calculate weights for asymptotic behaviour
1869 8400 : atermr(:) = alpha*tinv_beta(:)
1870 8400 : aterml(:) = inv_alph*tinv_beta(:)*((alphsq+tbetasq(:))*xtab(nomegaei+1)-tbetasq(:))
1871 8400 : logup(:) = alphsq*xtab(nomegaei+1)*xtab(nomegaei+1)
1872 : logdown(:) = ABS(((alphsq+tbetasq(:))*xtab(nomegaei+1)-two*tbetasq(:)) &
1873 8400 : *xtab(nomegaei+1)+tbetasq(:))
1874 8400 : right(:) = ATAN((atermr(:)-aterml(:))/(one+atermr(:)*aterml(:)))
1875 : weight (nomegaei+1,:) = CMPLX(-(half*(alphsq*tinv_beta(:)*LOG(logdown(:)/logup(:)) &
1876 8400 : - tbeta(:)*LOG(xtab(nomegaei+1)*xtab(nomegaei+1))) - alpha*right(:)),zero)
1877 : tfone(:,nomegaei+1) = -(zero-epsrho_imag(:,nomegaei+1)*AIMAG(omega_imag(nomegaei+1))) &
1878 13440 : /(one-xtab(nomegaei+1))
1879 :
1880 : ! Use BLAS call to perform matrix-matrix multiplication and accumulation
1881 840 : fact = CMPLX(piinv,zero)
1882 :
1883 : ! FIXME: Temporary copy of ket
1884 : call xgemm('N','N',npwc,nomega,nomegaei+1,fact,tfone,npwc,&
1885 840 : weight ,nomegaei+1,cone_gw,ket(spadc+1:spadc+npwc,:),npwc)
1886 : call xgemm('N','N',npwc,nomega,nomegaei ,fact,tftwo,npwc,&
1887 840 : weight2,nomegaei ,cone_gw,ket(spadc+1:spadc+npwc,:),npwc)
1888 :
1889 : case (NSPLINE)
1890 : ! Natural spline followed by Gauss-Kronrod
1891 : ! Transform omega coordinates
1892 2520 : alpha = plasmafreq
1893 2520 : alphsq = alpha*alpha
1894 2520 : inv_alph = one/alpha
1895 :
1896 25200 : xtab(1:nomegaei+1) = AIMAG(omega_imag(:))/(AIMAG(omega_imag(:)) + alpha)
1897 2520 : xtab(nomegaei+2) = one
1898 :
1899 : ! Gauss-Kronrod integration of spline fit of f(t)/(1-t) in transformed space
1900 : ! *** OPENMP SECTION *** Added by MS
1901 : !!$OMP PARALLEL DO PRIVATE(ig,ftab,ftab2,s,s2,r,r2,y,y2,work,work2,beta,betasq,inv_beta, &
1902 : !!$OMP intsign,io,ii,i,j,re_intG,re_intK,im_intG,im_intK,temp1,temp2,temp3,temp4, &
1903 : !!$OMP ttil,tau,ref,fint,imf,fint2,GKttab)
1904 365238 : do ig=1,npwc
1905 : ! Spline fit
1906 378000 : ftab (1:nomegaei+1) = REAL(epsrho_imag(ig,1:nomegaei+1))/(one-xtab(1:nomegaei+1))
1907 378000 : ftab2(1:nomegaei+1) = AIMAG(epsrho_imag(ig,1:nomegaei+1))/(one-xtab(1:nomegaei+1))
1908 37800 : ftab (nomegaei+2) = zero; ftab2(nomegaei+2) = zero
1909 : ! Explicit calculation of spline coefficients
1910 37800 : s = zero; s2 = zero
1911 378000 : do i = 1, nomegaei+2-1
1912 340200 : r = ( ftab (i+1) - ftab (i) ) / ( xtab(i+1) - xtab(i) )
1913 340200 : r2 = ( ftab2(i+1) - ftab2(i) ) / ( xtab(i+1) - xtab(i) )
1914 340200 : y (2,i) = r - s; y2(2,i) = r2 - s2
1915 378000 : s = r; s2 = r2
1916 : end do
1917 37800 : s = zero; s2 = zero
1918 37800 : r = zero; r2 = zero
1919 37800 : y(2,1) = zero; y2(2,1) = zero
1920 37800 : y(2,nomegaei+2) = zero; y2(2,nomegaei+2) = zero
1921 340200 : do i = 2, nomegaei+2-1
1922 302400 : y (2,i) = y (2,i) + r * y (2,i-1)
1923 302400 : y2(2,i) = y2(2,i) + r2 * y2(2,i-1)
1924 302400 : work (i) = two * ( xtab(i-1) - xtab(i+1) ) - r * s
1925 302400 : work2(i) = two * ( xtab(i-1) - xtab(i+1) ) - r2 * s2
1926 302400 : s = xtab(i+1) - xtab(i)
1927 302400 : s2 = s
1928 302400 : r = s / work (i)
1929 340200 : r2 = s2 / work2(i)
1930 : end do
1931 340200 : do j = 2, nomegaei+2-1
1932 302400 : i = nomegaei+2+1-j
1933 302400 : y (2,i) = ( ( xtab(i+1) - xtab(i) ) * y (2,i+1) - y (2,i) ) / work (i)
1934 340200 : y2(2,i) = ( ( xtab(i+1) - xtab(i) ) * y2(2,i+1) - y2(2,i) ) / work2(i)
1935 : end do
1936 378000 : do i = 1, nomegaei+2-1
1937 340200 : s = xtab(i+1) - xtab(i); s2 = s;
1938 340200 : r = y(2,i+1) - y(2,i); r2 = y2(2,i+1) - y2(2,i);
1939 340200 : y(3,i) = r / s; y2(3,i) = r2 / s2;
1940 340200 : y(2,i) = three * y(2,i); y2(2,i) = three * y2(2,i);
1941 340200 : y (1,i) = ( ftab (i+1) - ftab (i) ) / s - ( y (2,i) + r ) * s
1942 378000 : y2(1,i) = ( ftab2(i+1) - ftab2(i) ) / s2 - ( y2(2,i) + r2 ) * s2
1943 : end do
1944 : ! End of spline interpolation
1945 380520 : do ios=1,nomega
1946 340200 : beta = REAL(omegame0i_tmp(ios))
1947 340200 : betasq = beta*beta
1948 340200 : inv_beta = one/beta
1949 340200 : intsign = sign(half*piinv,beta)
1950 340200 : beta = ABS(beta)
1951 340200 : io = 1; re_intG = zero; re_intK = zero; im_intG = zero; im_intK = zero
1952 8202600 : do ii=1,GK_LEVEL
1953 : do
1954 10527150 : GKttab = two*alpha*xtab(io+1)/(beta-(beta-alpha)*xtab(io+1))-one
1955 10527150 : if (GKttab > KronN(ii)) EXIT
1956 7824600 : io = io + 1
1957 : end do
1958 7824600 : temp1 = half*(KronN(ii)+one)
1959 7824600 : temp2 = temp1 - half
1960 7824600 : temp3 = temp2*temp2
1961 7824600 : temp4 = half/(temp3 + quarter)
1962 7824600 : ttil = beta*temp1/(alpha-(alpha-beta)*temp1)
1963 7824600 : tau = ttil - xtab(io)
1964 7824600 : ref = ftab (io) + tau*(y (1,io)+tau*(y (2,io)+tau*y (3,io)))
1965 7824600 : fint (ii) = -ref*(one-ttil)*temp4
1966 7824600 : imf = ftab2(io) + tau*(y2(1,io)+tau*(y2(2,io)+tau*y2(3,io)))
1967 7824600 : fint2(ii) = -imf*(one-ttil)*temp4
1968 7824600 : re_intK = KronW(ii)*fint (ii)
1969 7824600 : im_intK = KronW(ii)*fint2(ii)
1970 8164800 : ket(spadc+ig,ios) = ket(spadc+ig,ios)+intsign*CMPLX(re_intK,im_intK)
1971 : end do ! ii
1972 : end do !ios
1973 : end do !ig
1974 : !!$OMP END PARALLEL DO
1975 :
1976 : end select !intmethod
1977 :
1978 324918 : local_one = one
1979 324918 : local_zero = zero
1980 :
1981 : ! ============================================
1982 : ! ==== Add contribution coming from poles ====
1983 : ! ============================================
1984 : ! First see if the contribution has been checked before the routine is entered
1985 324918 : if (present(calc_poles)) then
1986 0 : my_calc_poles = calc_poles
1987 : else
1988 : ! Otherwise check locally if there is a contribution
1989 2748444 : do ios=1,nomega
1990 2423526 : if (omegame0i_tmp(ios)>tol12) then
1991 615145 : if ((local_one-theta_mu_minus_e0i)<tol12) my_calc_poles(ios) = .FALSE.
1992 : end if
1993 2748444 : if (omegame0i_tmp(ios)<-tol12) then
1994 1808381 : if (theta_mu_minus_e0i<tol12) my_calc_poles(ios) = .FALSE.
1995 : end if
1996 : end do !ios
1997 : end if
1998 :
1999 2418133 : if (ANY(my_calc_poles(:))) then ! Make sure we only enter if necessary
2000 : ! *** OPENMP SECTION *** Added by MS
2001 : !!OMP !write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',xomp_get_num_threads()
2002 : !$OMP PARALLEL SHARED(npwc,nomega,nomegaer,theta_mu_minus_e0i,spadc,local_one,local_zero, &
2003 : !$OMP omega,epsrho,omegame0i_tmp,ket,my_calc_poles) &
2004 : !$OMP PRIVATE(ig,ios,rtmp_r,rtmp_i,tmp_x,tmp_y,rt_real,rt_imag,ct,ierr) REDUCTION(+:my_err)
2005 : !!OMP $ write(std_out,'(a,i0)') ' Entering openmp loop. Number of threads: ',xomp_get_num_threads()
2006 : !$OMP DO
2007 2237464 : do ig=1,npwc
2008 : ! Prepare the spline interpolation by filling at once the arrays rtmp_r, rtmp_i
2009 41062147 : call spline(DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),nomegaer,local_zero,local_zero,rtmp_r)
2010 41062147 : call spline(DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),nomegaer,local_zero,local_zero,rtmp_i)
2011 : ! call spline_complex( DBLE(omega(1:nomegaer)), epsrho(ig,1:nomegaer), nomegaer, zero, zero, rtmp )
2012 :
2013 32462203 : do ios=1,nomega
2014 30224739 : if (.NOT.my_calc_poles(ios)) CYCLE
2015 :
2016 : ! Interpolate real and imaginary part of epsrho at |omegame0i_tmp|.
2017 25061273 : tmp_x(1) = ABS(omegame0i_tmp(ios))
2018 500704547 : call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(epsrho(ig,1:nomegaer)),rtmp_r,1,tmp_x,tmp_y,ierr=ierr)
2019 25061273 : if (ig==1.and.ispinor==1) my_err = my_err + ierr
2020 25061273 : rt_real = tmp_y(1)
2021 :
2022 : tmp_x(1) = ABS(omegame0i_tmp(ios))
2023 500704547 : call splint(nomegaer,DBLE(omega(1:nomegaer)),DBLE(AIMAG(epsrho(ig,1:nomegaer))),rtmp_i,1,tmp_x,tmp_y)
2024 25061273 : rt_imag = tmp_y(1)
2025 : !!call splint_complex(nomegaer,DBLE(omega(1:nomegaer)),epsrho(ig,1:nomegaer),rtmp,1,tmp_x,yfit)
2026 :
2027 25061273 : ct=DCMPLX(rt_real,rt_imag)
2028 :
2029 25061273 : if (omegame0i_tmp(ios)>tol12) then
2030 6859194 : ket(spadc+ig,ios)=ket(spadc+ig,ios)+ct*(local_one-theta_mu_minus_e0i)
2031 : end if
2032 27231204 : if (omegame0i_tmp(ios)<-tol12) then
2033 18202079 : ket(spadc+ig,ios)=ket(spadc+ig,ios)-ct*theta_mu_minus_e0i
2034 : end if
2035 :
2036 : end do !ios
2037 : end do !ig
2038 : !$OMP END DO
2039 : !$OMP END PARALLEL
2040 : end if ! ANY(my_calc_poles)
2041 : end do !ispinor
2042 :
2043 317110 : npoles_missing = npoles_missing + my_err
2044 :
2045 317110 : if (INTMETHOD>2) then
2046 2520 : ABI_FREE(KronN)
2047 2520 : ABI_FREE(KronW)
2048 2520 : ABI_FREE(GaussW)
2049 2520 : ABI_FREE(fint)
2050 2520 : ABI_FREE(fint2)
2051 : end if
2052 :
2053 634220 : end subroutine calc_sigc_cd
2054 : !!***
2055 :
2056 : !!****f* ABINIT/calc_sig_ppm_comp
2057 : !!
2058 : !! NAME
2059 : !! calc_sig_ppm_comp
2060 : !!
2061 : !! FUNCTION
2062 : !! Calculate contributions to self-energy operator using a plasmon-pole model.
2063 : !!
2064 : !! INPUTS
2065 : !! nomega=number of frequencies to consider
2066 : !! npwc= number of G vectors in the plasmon pole
2067 : !! npwc1= 1 if ppmodel==3, =npwc if ppmodel== 4, 1 for all the other cases
2068 : !! npwc2= 1 if ppmodel==3, =1 if ppmodel== 4, 1 for all the other cases
2069 : !! npwx=number of G vectors in rhotwgp
2070 : !! ppmodel=plasmon pole model
2071 : !! theta_mu_minus_e0i= $\theta(\mu-\epsilon_{k-q,b1,s}), defines if the state is occupied or not
2072 : !! zcut=small imaginary part to avoid the divergence. (see related input variable)
2073 : !! omegame0i(nomega)=frequencies where evaluate \Sigma_c ($\omega$ - $\epsilon_i$
2074 : !! otq(npwc,npwc2)=plasmon pole parameters for this q-point
2075 : !! botsq(npwc,npwc1)=plasmon pole parameters for this q-point
2076 : !! eig(npwc,npwc)=the eigvectors of the symmetrized inverse dielectric matrix for this q point
2077 : !! (first index for G, second index for bands)
2078 : !! rhotwgp(npwx)=oscillator matrix elements divided by |q+G| i.e
2079 : !! $\frac{\langle b1 k-q s | e^{-i(q+G)r | b2 k s \rangle}{|q+G|}$
2080 : !!
2081 : !! OUTPUT
2082 : !! sigcme(nomega) (to be described), only relevant if ppm3 or ppm4
2083 : !!
2084 : !! ket(npwc,nomega):
2085 : !!
2086 : !! In case of ppmodel==1,2 it contains
2087 : !!
2088 : !! ket(G,omega) = Sum_G2 conjg(rhotw(G)) * Omega(G,G2) * rhotw(G2)
2089 : !! ---------------------------------------------------
2090 : !! 2 omegatw(G,G2) (omega-E_i + omegatw(G,G2)(2f-1))
2091 : !!
2092 : !! NOTES
2093 : !! Taken from old routine
2094 : !!
2095 : !! SOURCE
2096 :
2097 7620 : subroutine calc_sig_ppm_comp(npwc,nomega,rhotwgp,botsq,otq,omegame0i_io,zcut,theta_mu_minus_e0i,ket,ppmodel,npwx,npwc1,npwc2)
2098 :
2099 : !Arguments ------------------------------------
2100 : !scalars
2101 : integer,intent(in) :: nomega,npwc,npwc1,npwc2,npwx,ppmodel
2102 : real(dp),intent(in) :: omegame0i_io,theta_mu_minus_e0i,zcut
2103 : !arrays
2104 : complex(gwp),intent(in) :: botsq(npwc,npwc1),rhotwgp(npwx),otq(npwc,npwc2)
2105 : complex(gwp),intent(inout) :: ket(npwc,nomega)
2106 :
2107 : !Local variables-------------------------------
2108 : !scalars
2109 : integer :: ig,igp,io
2110 : real(dp) :: den,otw,twofm1_zcut
2111 : complex(gwp) :: num,rhotwgdp_igp
2112 : logical :: fully_occupied,totally_empty
2113 : character(len=500) :: msg
2114 : !arrays
2115 7620 : complex(gwp),allocatable :: ket_comp(:)
2116 : !*************************************************************************
2117 :
2118 7620 : if (ppmodel/=1.and.ppmodel/=2) then
2119 0 : write(msg,'(a,i0,a)')' The completeness trick cannot be used when ppmodel is ',ppmodel,' It should be set to 1 or 2. '
2120 0 : ABI_ERROR(msg)
2121 : end if
2122 :
2123 22860 : ABI_MALLOC(ket_comp,(npwc))
2124 458760 : ket_comp(:)=0.d0
2125 :
2126 7620 : fully_occupied=(abs(theta_mu_minus_e0i-1.)<0.001)
2127 7620 : totally_empty=(abs(theta_mu_minus_e0i)<0.001)
2128 :
2129 7620 : if(.not.(totally_empty)) then ! not totally empty
2130 90600 : twofm1_zcut=zcut
2131 90600 : do igp=1,npwc
2132 89040 : rhotwgdp_igp=rhotwgp(igp)
2133 5725920 : do ig=1,npwc
2134 5635320 : otw=DBLE(otq(ig,igp)) ! in principle otw -> otw - ieta
2135 5635320 : num = botsq(ig,igp)*rhotwgdp_igp
2136 :
2137 5635320 : den = omegame0i_io-otw
2138 5724360 : if (den**2>zcut**2) then
2139 5635320 : ket_comp(ig) = ket_comp(ig) - num/(den*otw)*theta_mu_minus_e0i
2140 : end if
2141 : end do !ig
2142 : end do !igp
2143 : end if ! not totally empty
2144 :
2145 7620 : if(.not.(fully_occupied)) then ! not fully occupied
2146 6100 : twofm1_zcut=-zcut
2147 :
2148 368560 : do igp=1,npwc
2149 362460 : rhotwgdp_igp=rhotwgp(igp)
2150 23170660 : do ig=1,npwc
2151 22802100 : otw=DBLE(otq(ig,igp)) ! in principle otw -> otw - ieta
2152 22802100 : num = botsq(ig,igp)*rhotwgdp_igp
2153 :
2154 22802100 : den = omegame0i_io-otw
2155 23164560 : if (den**2>zcut**2) then
2156 22802100 : ket_comp(ig) = ket_comp(ig) - num/(den*otw)*(1.-theta_mu_minus_e0i)
2157 : end if
2158 : end do !ig
2159 : end do !igp
2160 : end if ! not fully occupied
2161 :
2162 70440 : do io=1,nomega
2163 3975180 : ket(:,io)=ket(:,io)+0.5*ket_comp(:)
2164 : end do
2165 :
2166 7620 : ABI_FREE(ket_comp)
2167 :
2168 7620 : end subroutine calc_sig_ppm_comp
2169 : !!***
2170 :
2171 : end module m_sigc
2172 : !!***
|
