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:Gigy1986]] 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 660 : subroutine calc_sigx_me(sigmak_ibz, ikcalc, bmin, bmax, cryst, qp_ebands, dtset, Sigp, Sr, Gsph_x, Vcp, Kmesh, Qmesh, &
126 660 : ltg_k, Pawtab, Pawang, Paw_pwff, Pawfgrtab, Paw_onsite, psps, wfd, Wfdf, &
127 660 : 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 660 : integer,allocatable :: igfftxg0(:), igfftfxg0(:), x_gbound(:,:), gboundf(:,:)
173 660 : 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 660 : real(dp),contiguous, pointer :: qp_ene(:,:,:), qp_occ(:,:,:)
176 660 : real(dp),allocatable :: nhat12(:,:,:),grnhat12(:,:,:,:)
177 660 : complex(gwp),allocatable :: vc_sqrt_qbz(:), rhotwg(:), rhotwgp(:), rhotwg_ki(:,:), ur_bdgw(:,:), ur_ibz(:)
178 660 : complex(dp),allocatable :: sigxcme_tmp(:,:), sigxme_tmp(:,:,:), sigx(:,:,:,:)
179 660 : complex(gwp),allocatable :: ur_ae_sum(:),ur_ae_onsite_sum(:),ur_ps_onsite_sum(:)
180 660 : complex(gwp),allocatable :: ur_ae_bdgw(:,:),ur_ae_onsite_bdgw(:,:),ur_ps_onsite_bdgw(:,:)
181 660 : complex(gwp),contiguous, pointer :: cg_jb(:),cg_sum(:)
182 1320 : logical :: can_symmetrize(wfd%nsppol)
183 660 : logical,allocatable :: bks_mask(:,:,:)
184 660 : type(esymm_t),pointer :: QP_sym(:)
185 660 : type(sigijtab_t),pointer :: Sigxij_tab(:)
186 660 : type(pawcprj_type),allocatable :: Cprj_kgw(:,:),Cprj_ksum(:,:)
187 660 : type(pawpwij_t),allocatable :: Pwij_qg(:),Pwij_fft(:)
188 : !************************************************************************
189 :
190 : DBG_ENTER("COLL")
191 :
192 660 : call timab(430,1,tsec) ! csigme (SigX)
193 660 : call cwtime(cpu_all, wall_all, gflops_all, "start")
194 :
195 : ! Initialize some values.
196 660 : gwcalctyp = Sigp%gwcalctyp; nspinor = wfd%nspinor; nsppol = wfd%nsppol; npwx = sigp%npwx
197 660 : dim_rtwg = 1; if (nspinor == 2) dim_rtwg = 2
198 5940 : nsig_ab = sigp%nsig_ab; spinor_padx = reshape([0, 0, npwx, npwx, 0, npwx, npwx, 0], [2, 4])
199 660 : ABI_CHECK(Sigp%npwx == Gsph_x%ng, "Sigp%npwx != Gsph_x%ng")
200 :
201 660 : qp_ene => qp_ebands%eig; qp_occ => qp_ebands%occ
202 :
203 : ! Exctract the symmetries of the bands for this k-point
204 660 : 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 660 : jk_bz = Sigp%kptgw2bz(ikcalc)
208 660 : call kmesh%get_BZ_item(jk_bz, kgw, jk_ibz, isym_kgw, jik, ph_mkgwt)
209 3300 : spinrot_kgw(:) = cryst%spinrot(:,isym_kgw)
210 :
211 660 : write(msg,'(6a)') ch10, &
212 2640 : ' Calculating <nk|Sigma_x|nk> at k: ',trim(ktoa(kgw)), ", for bands: ", trim(ltoa([bmin, bmax])),ch10
213 660 : call wrtout(std_out, msg)
214 :
215 2211 : if (any(x_ngfft(1:3) /= wfd%ngfft(1:3)) ) then
216 143 : call wfd%change_ngfft(cryst, psps, x_ngfft)
217 143 : if (dtset%userie == 456) call wfdf%change_ngfft(Cryst, Psps, x_ngfft)
218 : end if
219 4620 : x_nfft = product(x_ngfft(1:3)); x_mgfft = maxval(x_ngfft(1:3)); x_fftalga = x_ngfft(7) / 100
220 :
221 660 : if (pawcross==1) mgfftf = MAXVAL(ngfftf(1:3))
222 :
223 : ! Define whether we can use symmetries to sum over the IBZ_kgw
224 1328 : can_symmetrize = .FALSE.
225 660 : if (Sigp%symsigma > 0) then
226 239 : can_symmetrize = .TRUE.
227 118 : 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 118 : 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 660 : call get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
247 :
248 : ! Table for \Sigmax_ij matrix elements.
249 660 : Sigxij_tab => Sigp%Sigxij_tab(ikcalc, 1:nsppol)
250 :
251 : ! Remove empty states from the list of states that will be distributed.
252 3300 : ABI_MALLOC(bks_mask, (wfd%mband, Kmesh%nbz, nsppol))
253 557588 : bks_mask = .FALSE.
254 :
255 1328 : do spin=1,nsppol
256 35376 : do ik_bz=1,Kmesh%nbz
257 34048 : ik_ibz = Kmesh%tab(ik_bz)
258 556928 : do band_sum=1,Sigp%nbnds
259 556260 : 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 2640 : ABI_MALLOC(proc_distrb, (wfd%mband, Kmesh%nbz, nsppol))
266 660 : call sigma_distribute_bks(wfd,Kmesh,ltg_k,Qmesh,nsppol,can_symmetrize,kgw,Sigp%mg0,my_nbks,proc_distrb,bks_mask=bks_mask)
267 660 : ABI_FREE(bks_mask)
268 660 : 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 1980 : 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 2640 : ABI_MALLOC(irottb, (x_nfft, cryst%nsym))
279 660 : call rotate_FFT_mesh(cryst%nsym, cryst%symrel, cryst%tnons, x_ngfft, irottb, iscompatibleFFT)
280 660 : 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 2640 : ABI_MALLOC(ktabr, (x_nfft, Kmesh%nbz))
285 34322 : do ik_bz=1,Kmesh%nbz
286 33662 : isym = Kmesh%tabo(ik_bz)
287 175031268 : do ifft=1,x_nfft
288 175030608 : ktabr(ifft,ik_bz) = irottb(ifft,isym)
289 : end do
290 : end do
291 660 : ABI_FREE(irottb)
292 :
293 660 : 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 660 : if (psps%usepaw == 1) then
310 84 : ABI_MALLOC(Cprj_ksum, (cryst%natom, nspinor))
311 9 : call pawcprj_alloc(Cprj_ksum, 0, wfd%nlmn_atm)
312 :
313 9 : 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 660 : nq_summed = Kmesh%nbz
323 660 : if (Sigp%symsigma > 0) then
324 236 : call ltg_k%print([std_out], prtvol=prtvol)
325 7698 : nq_summed = sum(ltg_k%ibzq(:))
326 : end if ! symsigma
327 :
328 660 : write(msg,'(2a,i0,a)')ch10,' calc_sigx_me: calculation status (', nq_summed, ' to be completed):'
329 660 : call wrtout(std_out, msg)
330 :
331 1980 : ABI_MALLOC(ur_ibz, (x_nfft * nspinor))
332 2640 : ABI_MALLOC(rhotwg_ki, (npwx * nspinor, bmin:bmax))
333 1980 : ABI_MALLOC(rhotwg, (npwx * nspinor))
334 1320 : ABI_MALLOC(rhotwgp, (npwx * nspinor))
335 1980 : ABI_MALLOC(vc_sqrt_qbz, (npwx))
336 :
337 72061 : ABI_CALLOC(sigxme_tmp, (bmin:bmax, bmin:bmax, nsppol * nsig_ab))
338 9392 : ABI_CALLOC(sigxcme_tmp, (bmin:bmax, nsppol * nsig_ab))
339 196079 : ABI_CALLOC(sigx, (2, bmin:bmax, bmin:bmax, nsppol * nsig_ab))
340 :
341 660 : 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 660 : call pstat_proc%print(_PSTAT_ARGS_)
348 :
349 1328 : do spin=1,nsppol
350 6936 : if (ALL(proc_distrb(:,:,spin) /= wfd%my_rank)) CYCLE ! Spin parallelism.
351 :
352 : ! ===============================================
353 : ! Load wavefunctions for Sigma_x matrix elements
354 : ! ===============================================
355 2632 : ABI_MALLOC_OR_DIE(ur_bdgw, (x_nfft * nspinor, bmin:bmax), ierr)
356 :
357 658 : if (dtset%userie == 456) 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 11716 : call wfd%get_many_ur([(jb, jb=bmin, bmax)], jk_ibz, spin, ur_bdgw)
362 : end if
363 :
364 658 : if (wfd%usepaw == 1) then
365 : ! Load cprj for GW states, note the indexing.
366 9 : dimcprj_gw = nspinor * (bmax - bmin + 1)
367 216 : ABI_MALLOC(Cprj_kgw, (cryst%natom, bmin:bmin+dimcprj_gw-1))
368 9 : call pawcprj_alloc(Cprj_kgw, 0, wfd%nlmn_atm)
369 9 : ibsp = bmin
370 55 : do jb=bmin,bmax
371 46 : call wfd%get_cprj(jb, jk_ibz, spin, cryst, Cprj_ksum, sorted=.FALSE.)
372 46 : call paw_symcprj(jk_bz, nspinor, 1, cryst, Kmesh, Pawtab, Pawang, Cprj_ksum)
373 46 : call pawcprj_copy(Cprj_ksum, Cprj_kgw(:,ibsp:ibsp+(nspinor-1)))
374 55 : ibsp = ibsp + nspinor
375 : end do
376 9 : 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 34696 : do ik_bz=1,Kmesh%nbz
394 :
395 : ! Parallelization over k-points and spin.
396 199205 : if (ALL(proc_distrb(:,ik_bz,spin) /= wfd%my_rank)) CYCLE
397 :
398 : ! Find the symmetrical image of ksum in the IBZ
399 27674 : call kmesh%get_BZ_item(ik_bz, ksum, ik_ibz, isym_ki, iik, ph_mkt)
400 138370 : spinrot_kbz = cryst%spinrot(:,isym_ki)
401 :
402 : ! Identify q and G0 where q + G0 = k_GW - ksum
403 110696 : kgw_m_ksum = kgw - ksum
404 27674 : 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 27674 : wtqp = 1; wtqm = 0
409 27674 : if (can_symmetrize(spin)) then
410 1408 : if (ltg_k%ibzq(iq_bz) /= 1) cycle
411 1408 : wtqp = 0; wtqm = 0
412 76648 : do isym=1,ltg_k%nsym_sg
413 75240 : wtqp = wtqp + ltg_k%wtksym(1, isym, iq_bz)
414 76648 : wtqm = wtqm + ltg_k%wtksym(2, isym, iq_bz)
415 : end do
416 : end if
417 :
418 27674 : print_time = wfd%my_rank == 0 .and. (ik_bz < LOG_MODK .or. mod(ik_bz, LOG_MODK) == 0)
419 5910 : 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 27674 : call qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym_q, itim_q)
424 27674 : 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 110696 : ABI_MALLOC(x_gbound, (2*x_mgfft+8, 2))
430 27674 : call Gsph_x%fft_tabs(g0, x_mgfft, x_ngfft, use_padfft, x_gbound, igfftxg0)
431 :
432 27674 : if (any(x_fftalga == [2, 4])) use_padfft = 0 ! Padded-FFT is not coded in rho_tw_g
433 : !use_padfft = 0
434 27674 : 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 27674 : 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 27674 : 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 192 : q0 = qbz !;if (q_is_gamma) q0 = (/0.00001_dp,0.00001_dp,0.00001_dp/) ! GW_Q0_DEFAULT
454 774 : ABI_MALLOC(Pwij_qg, (psps%ntypat))
455 192 : 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 2387720 : do ig=1,npwx
462 2360046 : ig_rot = Gsph_x%rottb(ig, itim_q, isym_q)
463 2387720 : vc_sqrt_qbz(ig_rot) = Vcp%vc_sqrt_resid(ig, iq_ibz)
464 : end do
465 :
466 : ! ==========================
467 : ! Sum over (occupied) bands
468 : ! ==========================
469 395540 : do band_sum=1,Sigp%nbnds
470 :
471 : ! Parallelism over bands.
472 367866 : if (proc_distrb(band_sum, ik_bz, spin) /= wfd%my_rank) CYCLE
473 :
474 : ! Skip empty states. MRM: allow negative occ numbers.
475 86383 : if (abs(qp_occ(band_sum, ik_ibz, spin)) < tol_empty) CYCLE
476 :
477 86383 : call wfd%get_ur(band_sum, ik_ibz, spin, ur_ibz)
478 :
479 86383 : 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 888 : call wfd%get_cprj(band_sum, ik_ibz, spin, cryst, Cprj_ksum, sorted=.FALSE.)
484 888 : call paw_symcprj(ik_bz, nspinor, 1, cryst, Kmesh, Pawtab, Pawang, Cprj_ksum)
485 888 : 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 821387 : 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 735004 : nspinor,rhotwg_ki(:,jb))
509 :
510 735004 : if (psps%usepaw == 1 .and. use_pawnhat0 == 0) then
511 : ! Add on-site contribution, projectors are already in BZ.
512 5184 : i2=jb; if (nspinor==2) i2=(2*jb-1)
513 5184 : spad = nspinor - 1
514 : call paw_rho_tw_g(cryst,Pwij_qg,npwx,nspinor,nspinor,Gsph_x%gvec,&
515 5184 : Cprj_ksum(:,:),Cprj_kgw(:,i2:i2+spad),rhotwg_ki(:,jb))
516 : end if
517 735004 : 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 1499336 : do ii=1,nspinor
528 764332 : spad = (ii-1) * npwx
529 98001140 : rhotwg_ki(spad+1:spad+npwx, jb) = rhotwg_ki(spad+1:spad + npwx, jb) * vc_sqrt_qbz(1:npwx)
530 : end do
531 :
532 821387 : 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 20017 : if (nspinor == 1) then
544 18865 : rhotwg_ki(1, jb) = czero_gw
545 18865 : 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 86383 : theta_mu_minus_esum = fact_spin * qp_occ(band_sum, ik_ibz, spin)
568 86383 : theta_mu_minus_esum2 = sqrt(abs(fact_spin * qp_occ(band_sum, ik_ibz, spin))) ! MBB Nat. orb. funct. approx. sqrt(occ)
569 :
570 114057 : if (abs(theta_mu_minus_esum / fact_spin) >= tol_empty) then ! MRM: allow negative occ numbers
571 821387 : do kb=bmin,bmax
572 : ! Copy the ket Sigma_x |phi_{k,kb}>.
573 97236808 : 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 4394930 : do irow=1,Sigxij_tab(spin)%col(kb)%size1
582 3573543 : jb = Sigxij_tab(spin)%col(kb)%bidx(irow)
583 341731546 : 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 7970074 : do iab=1,nsig_ab
588 3661527 : spadx1 = spinor_padx(1, iab); spadx2 = spinor_padx(2, iab)
589 3661527 : xdot_tmp = -XDOTC(npwx, rhotwg(spadx1+1:), 1, rhotwgp(spadx2+1:), 1)
590 3661527 : gwpc_sigxme = xdot_tmp * theta_mu_minus_esum
591 3661527 : 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 3661527 : is_idx = spin; if (nspinor == 2) is_idx = iab
596 : sigxme_tmp(jb, kb, is_idx) = sigxme_tmp(jb, kb, is_idx) + &
597 3661527 : (wtqp + wtqm) * DBLE(gwpc_sigxme) + (wtqp - wtqm) * j_gw * AIMAG(gwpc_sigxme)
598 3661527 : if (jb == kb) then
599 : sigxcme_tmp(jb, is_idx) = sigxcme_tmp(jb, is_idx) + &
600 822988 : (wtqp + wtqm) * DBLE(gwpc_sigxme2) + (wtqp - wtqm) *j_gw * AIMAG(gwpc_sigxme2)
601 : end if
602 :
603 3661527 : sigx(1, jb, kb, is_idx) = sigx(1, jb, kb, is_idx) + wtqp * gwpc_sigxme
604 7235070 : 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 27674 : if (print_time) then
613 5910 : write(msg,'(3(a,i0))')' sigx: ik_bz: ',ik_bz,'/',Kmesh%nbz,", spin: ",spin
614 5910 : 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 27674 : ABI_FREE(x_gbound)
619 27674 : if (pawcross==1) then
620 0 : ABI_FREE(gboundf)
621 : end if
622 :
623 28332 : if (psps%usepaw==1 .and. use_pawnhat0 == 0) then
624 192 : call pawpwij_free(Pwij_qg)
625 390 : ABI_FREE(Pwij_qg)
626 : end if
627 : end do ! ik_bz Got all diagonal (off-diagonal) matrix elements.
628 :
629 658 : ABI_FREE(ur_bdgw)
630 1318 : if (wfd%usepaw == 1) then
631 9 : call pawcprj_free(Cprj_kgw)
632 143 : ABI_FREE(Cprj_kgw)
633 9 : 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 660 : ABI_FREE(igfftxg0)
642 660 : if (pawcross==1) then
643 0 : ABI_FREE(igfftfxg0)
644 : end if
645 :
646 : ! Gather contributions from all the CPUs.
647 660 : call xmpi_sum(sigxme_tmp, wfd%comm, ierr)
648 660 : call xmpi_sum(sigxcme_tmp, wfd%comm, ierr)
649 660 : call xmpi_sum(sigx, wfd%comm, ierr)
650 :
651 : ! Multiply by constants. For 3D systems sqrt(4pi) is included in vc_sqrt_qbz.
652 69421 : sigxme_tmp = (one / (cryst%ucvol * Kmesh%nbz)) * sigxme_tmp * Sigp%sigma_mixing
653 7412 : sigxcme_tmp = (one / (cryst%ucvol * Kmesh%nbz)) * sigxcme_tmp * Sigp%sigma_mixing
654 193439 : 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 1328 : do spin=1,nsppol
662 668 : if (.not. can_symmetrize(spin)) cycle
663 1328 : 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 660 : 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 1328 : do spin=1,nsppol
680 6947 : do jb=bmin,bmax
681 12338 : do iab=1,nsig_ab
682 6051 : is_idx = spin; if (nsig_ab > 1) is_idx = iab
683 11670 : if (is_idx <= 2) then
684 5763 : Sr%sigxme(jb,jk_ibz,is_idx) = DBLE( sigxme_tmp(jb,jb,is_idx))
685 5763 : 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 69421 : Sr%x_mat(bmin:bmax, bmin:bmax, jk_ibz, :) = sigxme_tmp(bmin:bmax, bmin:bmax,:)
699 660 : ABI_FREE(sigxme_tmp)
700 660 : ABI_FREE(sigxcme_tmp)
701 :
702 : ! ===========================
703 : ! ==== Deallocate memory ====
704 : ! ===========================
705 660 : if (psps%usepaw == 1) then
706 9 : call pawcprj_free(Cprj_ksum)
707 48 : ABI_FREE(Cprj_ksum)
708 9 : 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 9 : 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 660 : ABI_FREE(ur_ibz)
725 660 : ABI_FREE(rhotwg_ki)
726 660 : ABI_FREE(rhotwg)
727 660 : ABI_FREE(rhotwgp)
728 660 : ABI_FREE(vc_sqrt_qbz)
729 660 : ABI_FREE(ktabr)
730 660 : ABI_FREE(sigx)
731 660 : ABI_FREE(proc_distrb)
732 :
733 660 : call timab(430,2,tsec) ! csigme (SigX)
734 660 : call cwtime_report(" calc_sigx_me:", cpu_all, wall_all, gflops_all)
735 :
736 : DBG_EXIT("COLL")
737 :
738 3300 : 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 121 : 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 121 : integer,allocatable :: degtab(:,:)
759 121 : 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 7253 : ABI_ICALLOC(degtab, (bmin:bmax, bmin:bmax))
767 886 : do ib=bmin,bmax
768 6769 : do jb=bmin,bmax
769 6648 : 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 9111 : ABI_CALLOC(sym_sigx, (bmin:bmax, bmin:bmax, nsig_ab))
774 :
775 : ! Average over degenerate diagonal elements.
776 886 : do ib=bmin,bmax
777 : ndegs=0
778 6648 : do jb=bmin,bmax
779 5883 : if (degtab(ib,jb)==1) then
780 1757 : if (nspinor == 1) then
781 4767 : 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 6648 : ndegs = ndegs + degtab(ib,jb)
789 : end do
790 1771 : 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 886 : do ib=bmin,bmax
797 6769 : do jb=bmin,bmax
798 6648 : if (nspinor == 1) then
799 5347 : 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 121 : ABI_FREE(sym_sigx)
809 121 : ABI_FREE(degtab)
810 :
811 121 : end subroutine sigx_symmetrize
812 : !!***
813 :
814 : end module m_sigx
815 : !!***
|
