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