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