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