Branch data Line data Source code
1 : : !!****m* ABINIT/m_sigx
2 : : !! NAME
3 : : !! m_sigx
4 : : !!
5 : : !! FUNCTION
6 : : !! Calculate diagonal and off-diagonal matrix elements of the exchange part of the self-energy operator.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1999-2026 ABINIT group (FB, GMR, VO, LR, RWG, MG, RShaltaf)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_sigx
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_gwdefs
27 : : use m_xmpi
28 : : use m_defs_ptgroups
29 : : use m_errors
30 : :
31 : : use defs_datatypes, only : pseudopotential_type
32 : : use m_dtset, only : dataset_type
33 : : use m_time, only : timab, cwtime, cwtime_report
34 : : use m_fstrings, only : itoa, sjoin, ktoa, ltoa
35 : : use m_hide_blas, only : xdotc, xgemv
36 : : use m_numeric_tools, only : hermitianize
37 : : use m_geometry, only : normv
38 : : use m_crystal, only : crystal_t
39 : : use m_fft_mesh, only : rotate_FFT_mesh, cigfft
40 : : use m_bz_mesh, only : kmesh_t, findqg0, littlegroup_t
41 : : use m_gsphere, only : gsphere_t
42 : : use m_vcoul, only : vcoul_t
43 : : use m_pawpwij, only : pawpwff_t, pawpwij_t, pawpwij_init, pawpwij_free, paw_rho_tw_g, paw_cross_rho_tw_g
44 : : use m_paw_pwaves_lmn,only : paw_pwaves_lmn_t
45 : : use m_pawang, only : pawang_type
46 : : use m_pawtab, only : pawtab_type
47 : : use m_pawfgrtab, only : pawfgrtab_type
48 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, paw_overlap
49 : : use m_paw_nhat, only : pawmknhat_psipsi
50 : : use m_paw_sym, only : paw_symcprj
51 : : use m_wfd, only : wfdgw_t, wave_t
52 : : use m_sigma, only : sigma_t, sigma_distribute_bks
53 : : use m_oscillators, only : rho_tw_g
54 : : use m_esymm, only : esymm_t
55 : : use m_occ, only : get_fact_spin_tol_empty
56 : : use m_ebands, only : ebands_t
57 : : use m_pstat, only : pstat_proc
58 : :
59 : : implicit none
60 : :
61 : : private
62 : : !!***
63 : :
64 : : public :: calc_sigx_me
65 : : public :: sigx_symmetrize ! Symmetrize Sig_x matrix elements
66 : : !!***
67 : :
68 : : integer,parameter :: LOG_MODK = 5
69 : :
70 : : contains
71 : : !!***
72 : :
73 : : !!****f* ABINIT/calc_sigx_me
74 : : !! NAME
75 : : !! calc_sigx_me
76 : : !!
77 : : !! FUNCTION
78 : : !! Calculate diagonal and off-diagonal matrix elements of the exchange part of the self-energy operator.
79 : : !!
80 : : !! INPUTS
81 : : !! sigmak_ibz=Index of the k-point in the IBZ.
82 : : !! bmin, bmax= min and Max band index for GW correction (for this k-point)
83 : : !! Gsph_x<gsphere_t>= Info on the G-sphere used for Sigma_x
84 : : !! ikcalc=index in the array Sigp%kptgw2bz of the k-point where GW corrections are calculated
85 : : !! ltg_k datatype containing information on the little group
86 : : !! Kmesh <kmesh_t>
87 : : !! x_ngfft(18)=Information about 3D FFT for the oscillator strengths, see ~abinit/doc/variables/vargs.htm#ngfft
88 : : !! Vcp <vcoul_t datatype> containing information on the cutoff technique
89 : : !! Pawtab(psps%ntypat) <type(pawtab_type)>=paw tabulated starting data
90 : : !! Pawang <type(pawang_type)>=paw angular mesh and related data
91 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
92 : : !! Qmesh <kmesh_t> : datatype gathering information of the q-mesh used
93 : : !! Sigp <sigparams_t> (see the definition of this structured datatype)
94 : : !! cryst<crystal_t>=Info on unit cell and symmetries
95 : : !! qp_ebands<ebands_t>=Datatype gathering info on the QP energies (KS if one shot)
96 : : !! Paw_pwff<pawpwff_t>=Form factor used to calculate the onsite mat. elements of a plane wave.
97 : : !! allQP_sym(%nkibz, %nsppol)<esymm_t>=Datatype collecting data on the irreducible representations of the
98 : : !! little group of kcalc in the KS representation as well as the symmetry of the bdgw_k states.
99 : : !! prtvol=Flags governing verbosity level.
100 : : !!
101 : : !! OUTPUT
102 : : !! Sr%x_mat(bmin:bmax,bmin:bmax,%nsppol*Sigp%nsig_ab)=Matrix elements of Sigma_x.
103 : : !!
104 : : !! NOTES
105 : : !! 1) The treatment of the divergence of Gygi+Baldereschi (PRB 1986) [[cite:Gygi1986]] is included.
106 : : !!
107 : : !! 2) On the symmetrization of Sigma matrix elements
108 : : !! If Sk = k+G0 then M_G(k, Sq)= e^{-i( Sq+G).t} M_{ S^-1(G} (k,q)
109 : : !! If -Sk = k+G0 then M_G(k,-Sq)= e^{-i(-Sq+G).t} M_{-S^-1(G)}^*(k,q)
110 : : !!
111 : : !! Notice the absence of G0 in the expression. Moreover, when we sum over the little group, it turns out
112 : : !! that there is a cancellation of the phase factor associated to the non-symmorphic operations due to a
113 : : !! similar term coming from the symmetrization of \epsilon^{-1}. Mind however that the nonsymmorphic phase
114 : : !! has to be considered when epsilon^-1 is reconstructed starting from the q-points in the IBZ.
115 : : !!
116 : : !! 3) the unitary transformation relating wavefunctions
117 : : !! at symmetric k-points should be taken into account during the symmetrization
118 : : !! of the oscillator matrix elements. In case of G_oW_o and GW_o calculations, however,
119 : : !! it is possible to make an invariant by just including all the degenerate states and
120 : : !! averaging the final results over the degenerate subset. Here we divide the states
121 : : !! where the QP energies are required into complexes. Note however that this approach is not
122 : : !! based on group theory, and it might lead to spurious results in case of accidental degeneracies.
123 : : !!
124 : :
125 : 665 : subroutine calc_sigx_me(sigmak_ibz, ikcalc, bmin, bmax, cryst, qp_ebands, dtset, Sigp, Sr, Gsph_x, Vcp, Kmesh, Qmesh, &
126 : 665 : ltg_k, Pawtab, Pawang, Paw_pwff, Pawfgrtab, Paw_onsite, psps, wfd, Wfdf, &
127 : 665 : allQP_sym, x_ngfft, ngfftf, prtvol, pawcross, tol_empty_in)
128 : :
129 : : !Arguments ------------------------------------
130 : : !scalars
131 : : integer,intent(in) :: sigmak_ibz,ikcalc,prtvol,bmin,bmax,pawcross
132 : : type(crystal_t),intent(in) :: cryst
133 : : type(ebands_t),target,intent(in) :: qp_ebands
134 : : type(kmesh_t),intent(in) :: Kmesh,Qmesh
135 : : type(vcoul_t),intent(in) :: Vcp
136 : : type(gsphere_t),intent(in) :: Gsph_x
137 : : type(littlegroup_t),intent(in) :: ltg_k
138 : : type(Pseudopotential_type),intent(in) :: psps
139 : : type(dataset_type),intent(in) :: dtset
140 : : type(sigparams_t),target,intent(in) :: Sigp
141 : : type(sigma_t),intent(inout) :: Sr
142 : : type(pawang_type),intent(in) :: Pawang
143 : : type(wfdgw_t),target,intent(inout) :: wfd,Wfdf
144 : : real(dp),intent(in) :: tol_empty_in
145 : : !arrays
146 : : integer,intent(in) :: x_ngfft(18),ngfftf(18)
147 : : type(Pawtab_type),intent(in) :: Pawtab(psps%ntypat)
148 : : type(pawpwff_t),intent(in) :: Paw_pwff(psps%ntypat*psps%usepaw)
149 : : type(esymm_t),target,intent(in) :: allQP_sym(wfd%nkibz, wfd%nsppol)
150 : : type(pawfgrtab_type),intent(inout) :: Pawfgrtab(cryst%natom*psps%usepaw)
151 : : type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(cryst%natom*psps%usepaw)
152 : :
153 : : !Local variables ------------------------------
154 : : !scalars
155 : : integer,parameter :: ndat1 = 1, use_pawnhat0 = 0, ider0 = 0
156 : : integer :: gwcalctyp,izero,iab,band_sum,ierr,ig,ig_rot,ii,iik,itim_q,i2
157 : : integer :: ik_bz, ik_ibz, isym_q, iq_bz, iq_ibz, spin, isym, jb, is_idx
158 : : integer :: jik,jk_bz,jk_ibz,kb,nspinor,nsppol,ifft
159 : : integer :: nq_summed,ibsp,dimcprj_gw,dim_rtwg, isym_kgw, isym_ki
160 : : integer :: spad, spadx1, spadx2, irow, npw_k, wtqm, wtqp
161 : : integer :: npwx, x_nfft, x_mgfft, x_fftalga, nsig_ab
162 : : integer :: nfftf, mgfftf, nhat12_grdim, my_nbks, use_padfft, use_padfftf
163 : : real(dp) :: cpu_all, wall_all, gflops_all, cpu_k, wall_k, gflops_k
164 : : real(dp) :: fact_spin, theta_mu_minus_esum, theta_mu_minus_esum2, tol_empty
165 : : complex(dp) :: ctmp,ph_mkgwt,ph_mkt
166 : : complex(gwp) :: gwpc_sigxme,gwpc_sigxme2,xdot_tmp
167 : : logical :: iscompatibleFFT, q_is_gamma, print_time
168 : : character(len=5000) :: msg
169 : : type(wave_t),pointer :: wave_sum, wave_jb
170 : : !arrays
171 : : integer :: g0(3), spinor_padx(2,4)
172 : 665 : integer,allocatable :: igfftxg0(:), igfftfxg0(:), x_gbound(:,:), gboundf(:,:)
173 : 665 : integer,allocatable :: ktabr(:,:),irottb(:,:),ktabrf(:,:), proc_distrb(:,:,:)
174 : : real(dp) :: ksum(3), kgw(3), kgw_m_ksum(3), qbz(3), q0(3), spinrot_kbz(4), spinrot_kgw(4), tsec(2)
175 : 665 : real(dp),contiguous, pointer :: qp_ene(:,:,:), qp_occ(:,:,:)
176 : 665 : real(dp),allocatable :: nhat12(:,:,:),grnhat12(:,:,:,:)
177 : 665 : complex(gwp),allocatable :: vc_sqrt_qbz(:), rhotwg(:), rhotwgp(:), rhotwg_ki(:,:), ur_bdgw(:,:), ur_ibz(:)
178 : 665 : complex(dp),allocatable :: sigxcme_tmp(:,:), sigxme_tmp(:,:,:), sigx(:,:,:,:)
179 : 665 : complex(gwp),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:)
180 : 665 : complex(gwp),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:)
181 : 665 : complex(gwp),contiguous, pointer :: cg_jb(:),cg_sum(:)
182 : 1330 : logical :: can_symmetrize(wfd%nsppol)
183 : 665 : logical,allocatable :: bks_mask(:,:,:)
184 : 665 : type(esymm_t),pointer :: QP_sym(:)
185 : 665 : type(sigijtab_t),pointer :: Sigxij_tab(:)
186 : 665 : type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
187 : 665 : type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
188 : : !************************************************************************
189 : :
190 : : DBG_ENTER("COLL")
191 : :
192 : 665 : call timab(430,1,tsec) ! csigme (SigX)
193 : 665 : call cwtime(cpu_all, wall_all, gflops_all, "start")
194 : :
195 : : ! Initialize some values.
196 : 665 : gwcalctyp = Sigp%gwcalctyp; nspinor = wfd%nspinor; nsppol = wfd%nsppol; npwx = sigp%npwx
197 : 665 : dim_rtwg = 1; if (nspinor == 2) dim_rtwg = 2
198 [ + + ]: 5985 : nsig_ab = sigp%nsig_ab; spinor_padx = reshape([0, 0, npwx, npwx, 0, npwx, npwx, 0], [2, 4])
199 [ - + ]: 665 : ABI_CHECK(Sigp%npwx == Gsph_x%ng, "Sigp%npwx != Gsph_x%ng")
200 : :
201 : 665 : qp_ene => qp_ebands%eig; qp_occ => qp_ebands%occ
202 : :
203 : : ! Exctract the symmetries of the bands for this k-point
204 : 665 : QP_sym => allQP_sym(sigmak_ibz,1:nsppol)
205 : :
206 : : ! Index of sigma_k k-point in the BZ array, its image in the IBZ and symmetries
207 : 665 : jk_bz = Sigp%kptgw2bz(ikcalc)
208 : 665 : call kmesh%get_BZ_item(jk_bz, kgw, jk_ibz, isym_kgw, jik, ph_mkgwt)
209 [ + + ]: 3325 : spinrot_kgw(:) = cryst%spinrot(:,isym_kgw)
210 : :
211 : 665 : write(msg,'(6a)') ch10, &
212 [ + - + + : 2660 : ' Calculating <nk|Sigma_x|nk> at k: ',trim(ktoa(kgw)), ", for bands: ", trim(ltoa([bmin, bmax])),ch10
+ - ]
213 : 665 : call wrtout(std_out, msg)
214 : :
215 [ + + + + : 2216 : if (any(x_ngfft(1:3) /= wfd%ngfft(1:3)) ) then
+ + ]
216 : 148 : call wfd%change_ngfft(cryst, psps, x_ngfft)
217 [ - + ]: 148 : if (dtset%userie == 456 .or. dtset%userie == 458) call wfdf%change_ngfft(Cryst, Psps, x_ngfft)
218 : : end if
219 [ + + + + ]: 4655 : x_nfft = product(x_ngfft(1:3)); x_mgfft = maxval(x_ngfft(1:3)); x_fftalga = x_ngfft(7) / 100
220 : :
221 [ - + - - ]: 665 : if (pawcross==1) mgfftf = MAXVAL(ngfftf(1:3))
222 : :
223 : : ! Define whether we can use symmetries to sum over the IBZ_kgw
224 [ + + ]: 1338 : can_symmetrize = .FALSE.
225 [ + + ]: 665 : if (Sigp%symsigma > 0) then
226 [ + + ]: 243 : can_symmetrize = .TRUE.
227 [ - + ]: 120 : if (gwcalctyp >= 20) then
228 [ # # ]: 0 : do spin=1,nsppol
229 : 0 : can_symmetrize(spin) = .not. QP_sym(spin)%failed()
230 [ # # ]: 0 : if (.not.can_symmetrize(spin)) then
231 : : write(msg,'(a,i0,4a)')&
232 : 0 : "Symmetrization cannot be performed for spin: ",spin,ch10,&
233 [ # # ]: 0 : "band classification encountered the following problem: ",ch10,TRIM(QP_sym(spin)%err_msg)
234 : 0 : ABI_WARNING(msg)
235 : : end if
236 : : end do
237 : : end if
238 [ + + ]: 120 : if (nspinor == 2) then
239 : 3 : ABI_WARNING('Symmetrization with nspinor=2 not implemented')
240 : : end if
241 : : end if
242 : :
243 : : ! MRM allow lower occ numbers
244 : : ! Set tolerance used to decide if a band is empty
245 : : ! and normalization of theta_mu_minus_esum. If nsppol == 2, qp_occ $\in [0,1]$
246 : 665 : call get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
247 : :
248 : : ! Table for \Sigmax_ij matrix elements.
249 : 665 : Sigxij_tab => Sigp%Sigxij_tab(ikcalc, 1:nsppol)
250 : :
251 : : ! Remove empty states from the list of states that will be distributed.
252 [ - + + - : 3325 : ABI_MALLOC(bks_mask, (wfd%mband, Kmesh%nbz, nsppol))
- + + - -
+ + - + -
+ - - + -
+ ]
253 [ + + + + : 558102 : bks_mask = .FALSE.
+ + ]
254 : :
255 [ + + ]: 1338 : do spin=1,nsppol
256 [ + + ]: 35426 : do ik_bz=1,Kmesh%nbz
257 : 34088 : ik_ibz = Kmesh%tab(ik_bz)
258 [ + + ]: 557437 : do band_sum=1,Sigp%nbnds
259 : 556764 : bks_mask(band_sum, ik_bz, spin) = (abs(qp_occ(band_sum, ik_ibz, spin)) >= tol_empty) ! MRM allow negative occ
260 : : end do
261 : : end do
262 : : end do
263 : :
264 : : ! Distribute tasks.
265 [ - + + - : 2660 : ABI_MALLOC(proc_distrb, (wfd%mband, Kmesh%nbz, nsppol))
- + + - -
+ + - - +
- + ]
266 : 665 : call sigma_distribute_bks(wfd,Kmesh,ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,proc_distrb,bks_mask=bks_mask)
267 [ - + ]: 665 : ABI_FREE(bks_mask)
268 : 665 : call wrtout(std_out, sjoin(" Will sum ", itoa(my_nbks) ," (b, k, s) occupied states in Sigma_x."))
269 : :
270 : : ! The index of G-G0 in the FFT mesh for the oscillators
271 : : ! Sigp%mG0 gives the MAX G0 component to account for umklapp.
272 [ - + + - : 1995 : ABI_MALLOC(igfftxg0, (Gsph_x%ng))
+ - + - -
+ - + ]
273 : :
274 : : ! Precompute the FFT index of $ R^{-1}(r-\tau)$
275 : : ! S = \transpose R^{-1} and k_BZ = S k_IBZ
276 : : ! irottb is the FFT index of $R^{-1} (r-\tau)$ used to symmetrize u_Sk.
277 : :
278 [ - + + - : 2660 : ABI_MALLOC(irottb, (x_nfft, cryst%nsym))
- + + - +
- + - - +
- + ]
279 : 665 : call rotate_FFT_mesh(cryst%nsym, cryst%symrel, cryst%tnons, x_ngfft, irottb, iscompatibleFFT)
280 [ - + ]: 665 : if (.not. iscompatibleFFT) then
281 : 0 : ABI_WARNING("FFT mesh is not compatible with symmetries. Results might be affected by large errors!")
282 : : end if
283 : :
284 [ - + + - : 2660 : ABI_MALLOC(ktabr, (x_nfft, Kmesh%nbz))
- + + - +
- + - - +
- + ]
285 [ + + ]: 34367 : do ik_bz=1,Kmesh%nbz
286 : 33702 : isym = Kmesh%tabo(ik_bz)
287 [ + + ]: 175150113 : do ifft=1,x_nfft
288 : 175149448 : ktabr(ifft,ik_bz) = irottb(ifft,isym)
289 : : end do
290 : : end do
291 [ - + ]: 665 : ABI_FREE(irottb)
292 : :
293 [ + + - + ]: 665 : if (psps%usepaw == 1 .and. pawcross == 1) then
294 [ # # ]: 0 : nfftf = PRODUCT(ngfftf(1:3))
295 [ # # # # : 0 : ABI_MALLOC(irottb, (nfftf, cryst%nsym))
# # # # #
# # # # #
# # ]
296 : 0 : call rotate_FFT_mesh(cryst%nsym, cryst%symrel, cryst%tnons, ngfftf, irottb, iscompatibleFFT)
297 : :
298 [ # # # # : 0 : ABI_MALLOC(ktabrf,(nfftf, Kmesh%nbz))
# # # # #
# # # # #
# # ]
299 [ # # ]: 0 : do ik_bz=1,Kmesh%nbz
300 : 0 : isym=Kmesh%tabo(ik_bz)
301 [ # # ]: 0 : do ifft=1,nfftf
302 : 0 : ktabrf(ifft,ik_bz)=irottb(ifft,isym)
303 : : end do
304 : : end do
305 [ # # ]: 0 : ABI_FREE(irottb)
306 : : end if
307 : :
308 : : ! Additional allocations for PAW.
309 [ + + ]: 665 : if (psps%usepaw == 1) then
310 [ - + + - : 98 : ABI_MALLOC(Cprj_ksum, (cryst%natom, nspinor))
- + + - +
- + - - +
- + + + +
+ ]
311 : 11 : call pawcprj_alloc(Cprj_ksum, 0, wfd%nlmn_atm)
312 : :
313 : 11 : nhat12_grdim = 0
314 : : if (use_pawnhat0 == 1) then
315 : : ! Compensation charge for \phi_a^*\phi_b
316 : : call wrtout(std_out, "Using nhat12")
317 : : ABI_MALLOC(nhat12 ,(2, x_nfft, nspinor**2))
318 : : ABI_MALLOC(grnhat12,(2, x_nfft, nspinor**2, 3*nhat12_grdim))
319 : : end if
320 : : end if
321 : :
322 : 665 : nq_summed = Kmesh%nbz
323 [ + + ]: 665 : if (Sigp%symsigma > 0) then
324 [ + + ]: 240 : call ltg_k%print([std_out], prtvol=prtvol)
325 [ + + ]: 7716 : nq_summed = sum(ltg_k%ibzq(:))
326 : : end if ! symsigma
327 : :
328 : 665 : write(msg,'(2a,i0,a)')ch10,' calc_sigx_me: calculation status (', nq_summed, ' to be completed):'
329 : 665 : call wrtout(std_out, msg)
330 : :
331 [ - + + - : 1995 : ABI_MALLOC(ur_ibz, (x_nfft * nspinor))
+ - + - -
+ - + ]
332 [ - + + - : 2660 : ABI_MALLOC(rhotwg_ki, (npwx * nspinor, bmin:bmax))
- + + - +
- + - - +
- + ]
333 [ - + + - : 1995 : ABI_MALLOC(rhotwg, (npwx * nspinor))
+ - + - -
+ - + ]
334 [ - + + - : 1330 : ABI_MALLOC(rhotwgp, (npwx * nspinor))
- + - + ]
335 [ - + + - : 1995 : ABI_MALLOC(vc_sqrt_qbz, (npwx))
+ - + - -
+ - + ]
336 : :
337 [ - + + - : 72149 : ABI_CALLOC(sigxme_tmp, (bmin:bmax, bmin:bmax, nsppol * nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + ]
338 [ - + + - : 9431 : ABI_CALLOC(sigxcme_tmp, (bmin:bmax, nsppol * nsig_ab))
- + + - +
- + - - +
- + + + +
+ ]
339 [ - + + - : 196255 : ABI_CALLOC(sigx, (2, bmin:bmax, bmin:bmax, nsppol * nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + + + ]
340 : :
341 [ - + ]: 665 : if (pawcross==1) then
342 [ # # # # : 0 : ABI_MALLOC(ur_ae_sum,(nfftf*nspinor))
# # # # #
# # # ]
343 [ # # # # : 0 : ABI_MALLOC(ur_ae_onsite_sum,(nfftf*nspinor))
# # # # ]
344 [ # # # # : 0 : ABI_MALLOC(ur_ps_onsite_sum,(nfftf*nspinor))
# # # # ]
345 : : end if
346 : :
347 : 665 : call pstat_proc%print(_PSTAT_ARGS_)
348 : :
349 [ + + ]: 1338 : do spin=1,nsppol
350 [ + + + + : 6941 : if (ALL(proc_distrb(:,:,spin) /= wfd%my_rank)) CYCLE ! Spin parallelism.
+ + + + ]
351 : :
352 : : ! ===============================================
353 : : ! Load wavefunctions for Sigma_x matrix elements
354 : : ! ===============================================
355 [ - + + - : 2652 : ABI_MALLOC_OR_DIE(ur_bdgw, (x_nfft * nspinor, bmin:bmax), ierr)
- + + - +
- + - + -
- + - + ]
356 : :
357 [ - + ]: 663 : if (dtset%userie == 456 .or. dtset%userie == 458) then
358 : 0 : call wrtout(std_out, "Taking states from Sigma^x_nk from supercell WFK file")
359 [ # # # # : 0 : call wfdf%get_many_ur([(jb, jb=bmin, bmax)], jk_ibz, spin, ur_bdgw)
# # # # ]
360 : : else
361 [ - + + + : 11749 : call wfd%get_many_ur([(jb, jb=bmin, bmax)], jk_ibz, spin, ur_bdgw)
+ - + + ]
362 : : end if
363 : :
364 [ + + ]: 663 : if (wfd%usepaw == 1) then
365 : : ! Load cprj for GW states, note the indexing.
366 : 11 : dimcprj_gw = nspinor * (bmax - bmin + 1)
367 [ - + + - : 248 : ABI_MALLOC(Cprj_kgw, (cryst%natom, bmin:bmin+dimcprj_gw-1))
- + + - +
- + - - +
- + - + +
+ + + ]
368 : 11 : call pawcprj_alloc(Cprj_kgw, 0, wfd%nlmn_atm)
369 : 11 : ibsp = bmin
370 [ + + ]: 65 : do jb=bmin,bmax
371 : 54 : call wfd%get_cprj(jb, jk_ibz, spin, cryst, Cprj_ksum, sorted=.FALSE.)
372 : 54 : call paw_symcprj(jk_bz, nspinor, 1, cryst, Kmesh, Pawtab, Pawang, Cprj_ksum)
373 : 54 : call pawcprj_copy(Cprj_ksum, Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
374 : 65 : ibsp = ibsp + nspinor
375 : : end do
376 [ - + ]: 11 : if (pawcross ==1) then
377 [ # # # # : 0 : ABI_MALLOC(ur_ae_bdgw,(nfftf*nspinor,bmin:bmax))
# # # # #
# # # # #
# # # # ]
378 [ # # # # : 0 : ABI_MALLOC(ur_ae_onsite_bdgw,(nfftf*nspinor,bmin:bmax))
# # # # #
# # # #
# ]
379 [ # # # # : 0 : ABI_MALLOC(ur_ps_onsite_bdgw,(nfftf*nspinor,bmin:bmax))
# # # # #
# # # #
# ]
380 [ # # ]: 0 : do jb=bmin,bmax
381 : : call wfdf%paw_get_aeur(jb,jk_ibz,spin,cryst,Paw_onsite,psps,Pawtab,Pawfgrtab,&
382 : 0 : ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
383 [ # # ]: 0 : ur_ae_bdgw(:,jb)=ur_ae_sum
384 [ # # ]: 0 : ur_ae_onsite_bdgw(:,jb)=ur_ae_onsite_sum
385 [ # # ]: 0 : ur_ps_onsite_bdgw(:,jb)=ur_ps_onsite_sum
386 : : end do
387 : : end if
388 : : end if
389 : :
390 : : ! ==============================
391 : : ! ==== Sum over k in the BZ ====
392 : : ! ==============================
393 [ + + ]: 34741 : do ik_bz=1,Kmesh%nbz
394 : :
395 : : ! Parallelization over k-points and spin.
396 [ + + + + : 199385 : if (ALL(proc_distrb(:,ik_bz,spin) /= wfd%my_rank)) CYCLE
+ + ]
397 : :
398 : : ! Find the symmetrical image of ksum in the IBZ
399 : 27704 : call kmesh%get_BZ_item(ik_bz, ksum, ik_ibz, isym_ki, iik, ph_mkt)
400 [ + + ]: 138520 : spinrot_kbz = cryst%spinrot(:,isym_ki)
401 : :
402 : : ! Identify q and G0 where q + G0 = k_GW - ksum
403 [ + + ]: 110816 : kgw_m_ksum = kgw - ksum
404 : 27704 : call findqg0(iq_bz, g0, kgw_m_ksum, Qmesh%nbz, Qmesh%bz, Sigp%mG0)
405 : :
406 : : ! If symmetries are exploited only q-points in the IBZ_k are computed.
407 : : ! In this case elements are weighted according to wtqp and wtqm. wtqm is for time-reversal.
408 : 27704 : wtqp = 1; wtqm = 0
409 [ + + ]: 27704 : if (can_symmetrize(spin)) then
410 [ + - ]: 1414 : if (ltg_k%ibzq(iq_bz) /= 1) cycle
411 : 1414 : wtqp = 0; wtqm = 0
412 [ + + ]: 76798 : do isym=1,ltg_k%nsym_sg
413 : 75384 : wtqp = wtqp + ltg_k%wtksym(1, isym, iq_bz)
414 : 76798 : wtqm = wtqm + ltg_k%wtksym(2, isym, iq_bz)
415 : : end do
416 : : end if
417 : :
418 [ + + + + : 27704 : print_time = wfd%my_rank == 0 .and. (ik_bz < LOG_MODK .or. mod(ik_bz, LOG_MODK) == 0)
+ + ]
419 : 5929 : if (print_time) call cwtime(cpu_k, wall_k, gflops_k, "start")
420 : :
421 : : ! Find the corresponding irreducible q-point.
422 : : ! NB: non-zero umklapp G_o is not allowed. There's a check in setup_sigma
423 : 27704 : call qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym_q, itim_q)
424 : 27704 : q_is_gamma = normv(qbz, cryst%gmet, "G") < GW_TOLQ0
425 : :
426 : : ! Tables for the FFT of the oscillators.
427 : : ! a) FFT index of G-G0.
428 : : ! b) x_gbound table for the zero-padded FFT performed in rhotwg.
429 [ + - - + : 110816 : ABI_MALLOC(x_gbound, (2*x_mgfft+8, 2))
+ - + - +
- + - - +
- + - + ]
430 : 27704 : call Gsph_x%fft_tabs(g0, x_mgfft, x_ngfft, use_padfft, x_gbound, igfftxg0)
431 : :
432 [ - + ]: 27704 : if (any(x_fftalga == [2, 4])) use_padfft = 0 ! Padded-FFT is not coded in rho_tw_g
433 : : !use_padfft = 0
434 [ - + ]: 27704 : if (use_padfft == 0) then
435 [ # # ]: 0 : ABI_FREE(x_gbound)
436 [ # # # # : 0 : ABI_MALLOC(x_gbound, (2*x_mgfft+8, 2*use_padfft))
# # # # ]
437 : : end if
438 : :
439 [ - + ]: 27704 : if (pawcross==1) then
440 [ # # # # : 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2))
# # # # #
# # # # #
# # # # ]
441 [ # # # # : 0 : ABI_MALLOC(igfftfxg0,(Gsph_x%ng))
# # # # #
# # # #
# ]
442 : 0 : call Gsph_x%fft_tabs(g0,mgfftf,ngfftf,use_padfftf,gboundf,igfftfxg0)
443 [ # # ]: 0 : if ( ANY(x_fftalga == [2, 4]) ) use_padfftf=0
444 [ # # ]: 0 : if (use_padfftf==0) then
445 [ # # ]: 0 : ABI_FREE(gboundf)
446 [ # # # # : 0 : ABI_MALLOC(gboundf,(2*mgfftf+8,2*use_padfftf))
# # # # ]
447 : : end if
448 : : end if
449 : :
450 [ + + ]: 27704 : if (psps%usepaw==1 .and. use_pawnhat0 == 0) then
451 : : ! Evaluate oscillator matrix elements
452 : : ! $ <phj/r|e^{-i(q+G)}|phi/r> - <tphj/r|e^{-i(q+G)}|tphi/r> $ in packed form
453 : 198 : q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
454 [ - + + - : 804 : ABI_MALLOC(Pwij_qg, (psps%ntypat))
+ - + - -
+ - + - +
+ + ]
455 : 198 : call pawpwij_init(Pwij_qg, npwx, q0, Gsph_x%gvec, cryst%rprimd, psps, Pawtab, Paw_pwff)
456 : : end if
457 : :
458 : : ! Get Fourier components of the Coulomb interaction in the BZ
459 : : ! In 3D systems, neglecting umklapp, vc(Sq,sG)=vc(q,G)=4pi/|q+G|
460 : : ! The same relation holds for 0-D systems, but not in 1-D or 2D systems. It depends on S.
461 [ + + ]: 2390256 : do ig=1,npwx
462 : 2362552 : ig_rot = Gsph_x%rottb(ig, itim_q, isym_q)
463 : 2390256 : vc_sqrt_qbz(ig_rot) = Vcp%vc_sqrt_resid(ig, iq_ibz)
464 : : end do
465 : :
466 : : ! ==========================
467 : : ! Sum over (occupied) bands
468 : : ! ==========================
469 [ + + ]: 395894 : do band_sum=1,Sigp%nbnds
470 : :
471 : : ! Parallelism over bands.
472 [ + + ]: 368190 : if (proc_distrb(band_sum, ik_bz, spin) /= wfd%my_rank) CYCLE
473 : :
474 : : ! Skip empty states. MRM: allow negative occ numbers.
475 [ + - ]: 86533 : if (abs(qp_occ(band_sum, ik_ibz, spin)) < tol_empty) CYCLE
476 : :
477 : 86533 : call wfd%get_ur(band_sum, ik_ibz, spin, ur_ibz)
478 : :
479 [ + + ]: 86533 : if (psps%usepaw == 1) then
480 : : ! Load cprj for point ksum, this spin or spinor and *THIS* band.
481 : : ! TODO MG I could avoid doing this but I have to exchange spin and bands ???
482 : : ! For sure there is a better way to do this!
483 : 942 : call wfd%get_cprj(band_sum, ik_ibz, spin, cryst, Cprj_ksum, sorted=.FALSE.)
484 : 942 : call paw_symcprj(ik_bz, nspinor, 1, cryst, Kmesh, Pawtab, Pawang, Cprj_ksum)
485 [ - + ]: 942 : if (pawcross==1) then
486 : : call wfdf%paw_get_aeur(band_sum,ik_ibz,spin,cryst,Paw_onsite,psps,Pawtab,Pawfgrtab,&
487 : 0 : ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum)
488 : : end if
489 : : end if
490 : :
491 : : ! Get all <k-q,band_sum,s|e^{-i(q+G).r}|s,jb,k>
492 [ + + ]: 821945 : do jb=bmin,bmax
493 : :
494 : : if (Psps%usepaw==1 .and. use_pawnhat0 == 1) then
495 : : ABI_ERROR("use_pawnhat is disabled")
496 : : i2=jb; if (nspinor==2) i2=(2*jb-1)
497 : : spad = nspinor - 1
498 : :
499 : : izero=0
500 : : call pawmknhat_psipsi(Cprj_ksum,Cprj_kgw(:,i2:i2+spad),ider0,izero,cryst%natom,&
501 : : cryst%natom,x_nfft,x_ngfft,nhat12_grdim,nspinor,cryst%ntypat,1,1,Pawang,Pawfgrtab,&
502 : : grnhat12,nhat12,pawtab)
503 : :
504 : : else
505 : : call rho_tw_g(nspinor,npwx,x_nfft,ndat1,x_ngfft,1,use_padfft,igfftxg0,x_gbound, &
506 : : ur_ibz ,iik,ktabr(:,ik_bz),ph_mkt ,spinrot_kbz, &
507 : : ur_bdgw(:,jb),jik,ktabr(:,jk_bz),ph_mkgwt,spinrot_kgw, &
508 : 735412 : nspinor,rhotwg_ki(:,jb))
509 : :
510 [ + + ]: 735412 : if (psps%usepaw == 1 .and. use_pawnhat0 == 0) then
511 : : ! Add on-site contribution, projectors are already in BZ.
512 [ - + ]: 5400 : i2=jb; if (nspinor==2) i2=(2*jb-1)
513 : 5400 : spad = nspinor - 1
514 : : call paw_rho_tw_g(cryst,Pwij_qg,npwx,nspinor,nspinor,Gsph_x%gvec,&
515 : 5400 : Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),rhotwg_ki(:,jb))
516 : : end if
517 [ + + - + ]: 735412 : if (psps%usepaw==1.and.pawcross==1) then ! Add paw cross term
518 : : call paw_cross_rho_tw_g(nspinor,npwx,nfftf,ngfftf,1,use_padfftf,igfftfxg0,gboundf,&
519 : : ur_ae_sum,ur_ae_onsite_sum,ur_ps_onsite_sum,iik,ktabrf(:,ik_bz),ph_mkt,spinrot_kbz,&
520 : : ur_ae_bdgw(:,jb),ur_ae_onsite_bdgw(:,jb),ur_ps_onsite_bdgw(:,jb),jik,ktabrf(:,jk_bz),ph_mkgwt,spinrot_kgw,&
521 : 0 : nspinor,rhotwg_ki(:,jb))
522 : : end if
523 : : end if
524 : :
525 : : ! Multiply by the square root of the Coulomb term
526 : : ! In 3-D systems, the factor sqrt(4pi) is included
527 [ + + ]: 1500152 : do ii=1,nspinor
528 : 764740 : spad = (ii-1) * npwx
529 [ + + ]: 98065516 : rhotwg_ki(spad+1:spad+npwx, jb) = rhotwg_ki(spad+1:spad + npwx, jb) * vc_sqrt_qbz(1:npwx)
530 : : end do
531 : :
532 [ + + ]: 821945 : if (ik_bz == jk_bz) then
533 : : ! Treat analytically the case q --> 0:
534 : : !
535 : : ! * The oscillator is evaluated at q = 0 as it is considered constant in the small cube around Gamma
536 : : ! while the Colulomb term is integrated out.
537 : : ! * If nspinor == 1, we have nonzero contribution only if band_sum == jb
538 : : ! * If nspinor == 2, we evaluate <band_sum,up|jb,up> and <band_sum,dwn|jb,dwn>,
539 : : ! and impose orthonormalization since npwwfn might be < npwvec.
540 : : ! * Note the use of i_sz_resid and not i_sz, to account for the possibility
541 : : ! to have generalized KS basis set from hybrid.
542 : :
543 [ + + ]: 20113 : if (nspinor == 1) then
544 : 18961 : rhotwg_ki(1, jb) = czero_gw
545 [ + + ]: 18961 : if (band_sum == jb) rhotwg_ki(1,jb) = cmplx(sqrt(Vcp%i_sz_resid), 0.0_gwp)
546 : : !rhotwg_ki(1,jb) = czero_gw ! DEBUG
547 : :
548 : : else
549 : 1152 : npw_k = wfd%npwarr(ik_ibz)
550 : 1152 : rhotwg_ki(1, jb) = zero; rhotwg_ki(npwx+1, jb) = zero
551 [ + + ]: 1152 : if (band_sum == jb) then
552 [ - + ]: 88 : ABI_CHECK(wfd%get_wave_ptr(band_sum, ik_ibz, spin, wave_sum, msg) == 0, msg)
553 : 88 : cg_sum => wave_sum%ug
554 [ - + ]: 88 : ABI_CHECK(wfd%get_wave_ptr(jb, jk_ibz, spin, wave_jb, msg) == 0, msg)
555 : 88 : cg_jb => wave_jb%ug
556 : 88 : ctmp = xdotc(npw_k, cg_sum(1:), 1, cg_jb(1:), 1)
557 : 88 : rhotwg_ki(1, jb) = cmplx(sqrt(Vcp%i_sz_resid), 0.0_gwp) * real(ctmp)
558 : 88 : ctmp = xdotc(npw_k, cg_sum(npw_k+1:), 1, cg_jb(npw_k+1:), 1)
559 : 88 : rhotwg_ki(npwx+1, jb) = cmplx(sqrt(Vcp%i_sz_resid), 0.0_gwp) * real(ctmp)
560 : : end if
561 : : !rhotwg_ki(1, jb) = zero; rhotwg_ki(npwx+1, jb) = zero
562 : : ! PAW is missing
563 : : end if
564 : : end if
565 : : end do ! jb Got all matrix elements from bmin up to bmax.
566 : :
567 : 86533 : theta_mu_minus_esum = fact_spin * qp_occ(band_sum, ik_ibz, spin)
568 : 86533 : theta_mu_minus_esum2 = sqrt(abs(fact_spin * qp_occ(band_sum, ik_ibz, spin))) ! MBB Nat. orb. funct. approx. sqrt(occ)
569 : :
570 [ + - ]: 114237 : if (abs(theta_mu_minus_esum / fact_spin) >= tol_empty) then ! MRM: allow negative occ numbers
571 [ + + ]: 821945 : do kb=bmin,bmax
572 : : ! Copy the ket Sigma_x |phi_{k,kb}>.
573 [ + + ]: 97300776 : rhotwgp(:) = rhotwg_ki(:, kb)
574 : :
575 : : ! Loop over the non-zero row elements of this column.
576 : : ! If gwcalctyp < 20: only diagonal elements since QP == KS.
577 : : ! If gwcalctyp >= 20:
578 : : ! * Only off-diagonal elements connecting states with same character.
579 : : ! * Only the upper triangle if HF, SEX, or COHSEX.
580 : :
581 [ + + ]: 4395896 : do irow=1,Sigxij_tab(spin)%col(kb)%size1
582 : 3573951 : jb = Sigxij_tab(spin)%col(kb)%bidx(irow)
583 [ + + ]: 341795514 : rhotwg(:) = rhotwg_ki(:,jb)
584 : :
585 : : ! Calculate bare exchange <phi_jb|Sigma_x|phi_kb>.
586 : : ! Do the scalar product only if band_sum is occupied.
587 [ + + ]: 7971298 : do iab=1,nsig_ab
588 : 3661935 : spadx1 = spinor_padx(1, iab); spadx2 = spinor_padx(2, iab)
589 : 3661935 : xdot_tmp = -XDOTC(npwx, rhotwg(spadx1+1:), 1, rhotwgp(spadx2+1:), 1)
590 : 3661935 : gwpc_sigxme = xdot_tmp * theta_mu_minus_esum
591 : 3661935 : gwpc_sigxme2 = xdot_tmp * theta_mu_minus_esum2
592 : :
593 : : ! Accumulate and symmetrize Sigma_x matrix elements.
594 : : ! -wtqm comes from time-reversal (exchange of band indices)
595 [ + + ]: 3661935 : is_idx = spin; if (nspinor == 2) is_idx = iab
596 : : sigxme_tmp(jb, kb, is_idx) = sigxme_tmp(jb, kb, is_idx) + &
597 : 3661935 : (wtqp + wtqm) * DBLE(gwpc_sigxme) + (wtqp - wtqm) * j_gw * AIMAG(gwpc_sigxme)
598 [ + + ]: 3661935 : if (jb == kb) then
599 : : sigxcme_tmp(jb, is_idx) = sigxcme_tmp(jb, is_idx) + &
600 : 823396 : (wtqp + wtqm) * DBLE(gwpc_sigxme2) + (wtqp - wtqm) *j_gw * AIMAG(gwpc_sigxme2)
601 : : end if
602 : :
603 : 3661935 : sigx(1, jb, kb, is_idx) = sigx(1, jb, kb, is_idx) + wtqp * gwpc_sigxme
604 : 7235886 : sigx(2, jb, kb, is_idx) = sigx(2, jb, kb, is_idx) + wtqm *CONJG(gwpc_sigxme)
605 : : end do
606 : : end do ! jb
607 : : end do ! kb
608 : : end if
609 : :
610 : : end do ! band_sum
611 : :
612 [ + + ]: 27704 : if (print_time) then
613 : 5929 : write(msg,'(3(a,i0))')' sigx: ik_bz: ',ik_bz,'/',Kmesh%nbz,", spin: ",spin
614 [ + + ]: 5929 : call cwtime_report(msg, cpu_k, wall_k, gflops_k); if (ik_bz == LOG_MODK) call wrtout(std_out, " ...")
615 : : end if
616 : :
617 : : ! Deallocate k-dependent quantities.
618 [ - + ]: 27704 : ABI_FREE(x_gbound)
619 [ - + ]: 27704 : if (pawcross==1) then
620 [ # # ]: 0 : ABI_FREE(gboundf)
621 : : end if
622 : :
623 [ + + ]: 28367 : if (psps%usepaw==1 .and. use_pawnhat0 == 0) then
624 : 198 : call pawpwij_free(Pwij_qg)
625 [ + - + + : 408 : ABI_FREE(Pwij_qg)
- + - + ]
626 : : end if
627 : : end do ! ik_bz Got all diagonal (off-diagonal) matrix elements.
628 : :
629 [ - + ]: 663 : ABI_FREE(ur_bdgw)
630 [ + + ]: 1328 : if (wfd%usepaw == 1) then
631 : 11 : call pawcprj_free(Cprj_kgw)
632 [ + - + + : 161 : ABI_FREE(Cprj_kgw)
- + - + -
+ ]
633 [ - + ]: 11 : if (pawcross==1) then
634 [ # # ]: 0 : ABI_FREE(ur_ae_bdgw)
635 [ # # ]: 0 : ABI_FREE(ur_ae_onsite_bdgw)
636 [ # # ]: 0 : ABI_FREE(ur_ps_onsite_bdgw)
637 : : end if
638 : : end if
639 : : end do !spin
640 : :
641 [ - + ]: 665 : ABI_FREE(igfftxg0)
642 [ - + ]: 665 : if (pawcross==1) then
643 [ # # ]: 0 : ABI_FREE(igfftfxg0)
644 : : end if
645 : :
646 : : ! Gather contributions from all the CPUs.
647 : 665 : call xmpi_sum(sigxme_tmp, wfd%comm, ierr)
648 : 665 : call xmpi_sum(sigxcme_tmp, wfd%comm, ierr)
649 : 665 : call xmpi_sum(sigx, wfd%comm, ierr)
650 : :
651 : : ! Multiply by constants. For 3D systems sqrt(4pi) is included in vc_sqrt_qbz.
652 [ + + + + : 69489 : sigxme_tmp = (one / (cryst%ucvol * Kmesh%nbz)) * sigxme_tmp * Sigp%sigma_mixing
+ + ]
653 [ + + + + ]: 7436 : sigxcme_tmp = (one / (cryst%ucvol * Kmesh%nbz)) * sigxcme_tmp * Sigp%sigma_mixing
654 [ + + + + : 193595 : sigx = (one / (cryst%ucvol * Kmesh%nbz)) * sigx * Sigp%sigma_mixing
+ + + + ]
655 : :
656 : : ! If we have summed over the IBZ_q, we have to average over degenerate states.
657 : : ! NOTE: Presently only diagonal terms are considered
658 : : ! TODO QP-SCGW required a more involved approach, there is a check in sigma
659 : : ! TODO it does not work if spinor == 2.
660 : :
661 [ + + ]: 1338 : do spin=1,nsppol
662 [ + + ]: 673 : if (.not. can_symmetrize(spin)) cycle
663 : 1338 : call sigx_symmetrize(jk_ibz, spin, bmin, bmax, nsppol, nspinor, nsig_ab, dtset%symsigma_de, qp_ene, sigx, sigxme_tmp)
664 : : end do
665 : :
666 [ + + ]: 665 : if (gwcalctyp >= 20) then
667 : : ! Reconstruct the full sigma_x matrix from the upper triangle.
668 [ + - ]: 400 : if (nspinor == 1) then
669 [ + + ]: 802 : do spin=1,nsppol
670 : 802 : call hermitianize(sigxme_tmp(:,:,spin), "Upper")
671 : : end do
672 : : else
673 : 0 : ABI_WARNING("Should hermitianize non-collinear sigma!")
674 : : end if
675 : : end if
676 : :
677 : : ! Save diagonal elements or ab components of Sigma_x (Hermitian)
678 : : ! TODO It should be hermitian also if nspinor == 2
679 [ + + ]: 1338 : do spin=1,nsppol
680 [ + + ]: 6971 : do jb=bmin,bmax
681 [ + + ]: 12371 : do iab=1,nsig_ab
682 [ + + ]: 6065 : is_idx = spin; if (nsig_ab > 1) is_idx = iab
683 [ + + ]: 11698 : if (is_idx <= 2) then
684 : 5777 : Sr%sigxme(jb,jk_ibz,is_idx) = DBLE( sigxme_tmp(jb,jb,is_idx))
685 : 5777 : Sr%sigxcnofme(jb,jk_ibz,is_idx) = DBLE(sigxcme_tmp(jb,is_idx))
686 : : else
687 : 288 : Sr%sigxme(jb,jk_ibz,is_idx) = sigxme_tmp(jb,jb,is_idx)
688 : 288 : Sr%sigxcnofme(jb,jk_ibz,is_idx) = sigxcme_tmp(jb,is_idx)
689 : : end if
690 : : end do
691 : : !if (nsig_ab > 1) then
692 : : ! write(std_out,'(i3,4f8.3,a,f8.3)')jb,Sr%sigxme(jb,jk_ibz,:)*Ha_eV,' Tot ',SUM(Sr%sigxme(jb,jk_ibz,:))*Ha_eV
693 : : !end if
694 : : end do
695 : : end do
696 : :
697 : : ! Save full exchange matrix in Sr%
698 [ + + + + : 69489 : Sr%x_mat(bmin:bmax, bmin:bmax, jk_ibz, :) = sigxme_tmp(bmin:bmax, bmin:bmax,:)
+ + ]
699 [ - + ]: 665 : ABI_FREE(sigxme_tmp)
700 [ - + ]: 665 : ABI_FREE(sigxcme_tmp)
701 : :
702 : : ! ===========================
703 : : ! ==== Deallocate memory ====
704 : : ! ===========================
705 [ + + ]: 665 : if (psps%usepaw == 1) then
706 : 11 : call pawcprj_free(Cprj_ksum)
707 [ + - + + : 54 : ABI_FREE(Cprj_ksum)
- + - + -
+ ]
708 [ - + ]: 11 : if (allocated(Pwij_fft)) then
709 : 0 : call pawpwij_free(Pwij_fft)
710 [ # # # # : 0 : ABI_FREE(Pwij_fft)
# # # # ]
711 : : end if
712 : : if (use_pawnhat0 == 1) then
713 : : ABI_FREE(nhat12)
714 : : ABI_FREE(grnhat12)
715 : : end if
716 [ + - ]: 11 : if (pawcross == 1) then
717 [ # # ]: 0 : ABI_FREE(ur_ae_sum)
718 [ # # ]: 0 : ABI_FREE(ur_ae_onsite_sum)
719 [ # # ]: 0 : ABI_FREE(ur_ps_onsite_sum)
720 [ # # ]: 0 : ABI_FREE(ktabrf)
721 : : end if
722 : : end if
723 : :
724 [ - + ]: 665 : ABI_FREE(ur_ibz)
725 [ - + ]: 665 : ABI_FREE(rhotwg_ki)
726 [ - + ]: 665 : ABI_FREE(rhotwg)
727 [ - + ]: 665 : ABI_FREE(rhotwgp)
728 [ - + ]: 665 : ABI_FREE(vc_sqrt_qbz)
729 [ - + ]: 665 : ABI_FREE(ktabr)
730 [ - + ]: 665 : ABI_FREE(sigx)
731 [ - + ]: 665 : ABI_FREE(proc_distrb)
732 : :
733 : 665 : call timab(430,2,tsec) ! csigme (SigX)
734 : 665 : call cwtime_report(" calc_sigx_me:", cpu_all, wall_all, gflops_all)
735 : :
736 : : DBG_EXIT("COLL")
737 : :
738 [ - + - + : 3325 : end subroutine calc_sigx_me
- + - + -
+ - + - +
- + - - -
- - + - -
- - - + -
+ - + - +
- - - - -
- - + - -
- - - - ]
739 : : !!***
740 : :
741 : : !!****f* ABINIT/sigx_symmetrize
742 : : !! NAME
743 : : !! sigx_symmetrize
744 : : !!
745 : : !! FUNCTION
746 : : !! Symmetrize Sig_x matrix elements
747 : : !!
748 : :
749 [ - + ]: 123 : subroutine sigx_symmetrize(jk_ibz, spin, bmin, bmax, nsppol, nspinor, nsig_ab, symsigma_de, qp_ene, sigx, sigxme_tmp)
750 : :
751 : : integer,intent(in) :: jk_ibz, spin, bmin, bmax, nsppol, nspinor, nsig_ab
752 : : real(dp),intent(in) :: symsigma_de, qp_ene(:,:,:)
753 : : complex(dp),intent(in) :: sigx(2, bmin:bmax, bmin:bmax, nsppol * nsig_ab)
754 : : complex(dp),intent(inout) :: sigxme_tmp(bmin:bmax, bmin:bmax, nsppol * nsig_ab)
755 : :
756 : : !Local variables ------------------------------
757 : : integer :: ib, jb, ndegs, ii
758 : 123 : integer,allocatable :: degtab(:,:)
759 : 123 : complex(dp),allocatable :: sym_sigx(:,:,:)
760 : : !************************************************************************
761 : :
762 : : ! Find number of degenerates subspaces and number of bands in each subspace.
763 : : ! The tolerance is a little bit arbitrary (0.001 eV)
764 : : ! It could be reduced, in particular in case of nearly accidental degeneracies
765 : :
766 [ - + + - : 7303 : ABI_ICALLOC(degtab, (bmin:bmax, bmin:bmax))
- + + - +
- + - - +
+ - + + +
+ ]
767 [ + + ]: 896 : do ib=bmin,bmax
768 [ + + ]: 6811 : do jb=bmin,bmax
769 [ + + ]: 6688 : if (abs(qp_ene(ib, jk_ibz, spin) - qp_ene(jb, jk_ibz, spin)) < symsigma_de) degtab(ib, jb) = 1
770 : : end do
771 : : end do
772 : :
773 [ - + + - : 9163 : ABI_CALLOC(sym_sigx, (bmin:bmax, bmin:bmax, nsig_ab))
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + ]
774 : :
775 : : ! Average over degenerate diagonal elements.
776 [ + + ]: 896 : do ib=bmin,bmax
777 : : ndegs=0
778 [ + + ]: 6688 : do jb=bmin,bmax
779 [ + + ]: 5915 : if (degtab(ib,jb)==1) then
780 [ + + ]: 1777 : if (nspinor == 1) then
781 [ + + ]: 4827 : sym_sigx(ib, ib, 1) = sym_sigx(ib, ib, 1) + sum(sigx(:,jb,jb,spin))
782 : : else
783 [ + + ]: 840 : do ii=1,nsig_ab
784 [ + + ]: 2184 : sym_sigx(ib, ib, ii) = sym_sigx(ib, ib, ii) + sum(sigx(:,jb,jb,ii))
785 : : end do
786 : : end if
787 : : end if
788 : 6688 : ndegs = ndegs + degtab(ib,jb)
789 : : end do
790 [ + + ]: 1789 : sym_sigx(ib,ib,:) = sym_sigx(ib,ib,:) / ndegs
791 : : end do
792 : :
793 : : !if (gwcalctyp >= 20) call QP_sym(spin)%symmetrize_mels(,bmin,bmax,sigx(:,:,:,spin),sym_sigx(:,:,1))
794 : :
795 : : ! Copy symmetrized values.
796 [ + + ]: 896 : do ib=bmin,bmax
797 [ + + ]: 6811 : do jb=bmin,bmax
798 [ + + ]: 6688 : if (nspinor == 1) then
799 : 5379 : sigxme_tmp(ib,jb,spin) = sym_sigx(ib,jb,1)
800 : : else
801 [ + + ]: 2680 : do ii=1,nsig_ab
802 : 2680 : sigxme_tmp(ib,jb,ii) = sym_sigx(ib,jb,ii)
803 : : end do
804 : : end if
805 : : end do
806 : : end do
807 : :
808 [ - + ]: 123 : ABI_FREE(sym_sigx)
809 [ - + ]: 123 : ABI_FREE(degtab)
810 : :
811 : 123 : end subroutine sigx_symmetrize
812 : : !!***
813 : :
814 : : end module m_sigx
815 : : !!***
|
