Line data Source code
1 : !!****m* ABINIT/m_chi0tk
2 : !! NAME
3 : !! m_chi0tk
4 : !!
5 : !! FUNCTION
6 : !! This module provides tools for the computation of the irreducible polarizability.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (MG, FB)
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_chi0tk
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xmpi
28 : use m_xomp
29 : use m_sort
30 :
31 : use m_ebands, only : ebands_t
32 : use m_gwdefs, only : GW_TOL_DOCC, czero_gw, cone_gw, one_gw, em1params_t, j_gw
33 : use m_fstrings, only : sjoin, itoa
34 : use m_hide_blas,only : xgerc, xgemm, xherk, xher
35 : use m_crystal, only : crystal_t
36 : use m_gsphere, only : gsphere_t
37 : use m_bz_mesh, only : littlegroup_t, kmesh_t
38 : use m_wfd, only : wfdgw_t
39 :
40 : implicit none
41 :
42 : private
43 :
44 : public :: assemblychi0_sym
45 : public :: symmetrize_afm_chi0
46 : public :: accumulate_chi0_q0
47 : public :: accumulate_head_wings_imagw
48 : public :: accumulate_sfchi0_q0
49 : public :: assemblychi0sf
50 : public :: approxdelta
51 : public :: setup_spectral
52 : public :: hilbert_transform
53 : public :: hilbert_transform_headwings
54 : public :: completechi0_deltapart
55 : public :: output_chi0sumrule
56 : public :: accumulate_chi0sumrule
57 : public :: make_transitions
58 : public :: chi0_bbp_mask
59 : !!***
60 :
61 : CONTAINS !=======================================================================================================
62 : !!***
63 :
64 : !----------------------------------------------------------------------
65 :
66 : !!****f* m_chi0tk/assemblychi0_sym
67 : !! NAME
68 : !! assemblychi0_sym
69 : !!
70 : !! FUNCTION
71 : !! Update the independent particle susceptibility for the contribution
72 : !! of one pair of occupied-unoccupied band, for each frequency.
73 : !! If symchi=1 the expression is symmetrized taking into account the symmetries
74 : !! of the little group associated to the external q-point.
75 : !! Compute chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S \hat S (rhotwg(G1)*rhotwg*(G2))*green_w(io)
76 : !! where S are the symmetries of the little group associated to the external q-point.
77 : !!
78 : !! INPUTS
79 : !! nspinor=Number of spinorial components.
80 : !! ik_bz=Index of the k-point in the BZ array whose contribution has to be symmetrized and added to cchi0
81 : !! npwepG0=Maximum number of G vectors taking into account possible umklapp G0, ie enlarged sphere G-G0
82 : !! rhotwg(npwe)=Oscillator matrix elements for this k-point and the transition that has to be summed
83 : !! green_w(nomega)=frequency dependent part coming from the green function
84 : !! Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors
85 : !! %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg
86 : !! %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion
87 : !! %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations
88 : !! Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point.
89 : !! %timrev=2 it time-reversal is used, 1 otherwise
90 : !! %nsym_sg=Number of space group symmetries
91 : !! %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point
92 : !! %flag_umklp(timrev,nsym)= flag for umklapp processes
93 : !! if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise
94 : !! %igmG0(npwepG0,timrev,nsym) index of G-G0 in the array gvec
95 : !! Ep<em1params_t>=Parameters related to the calculation of chi0/epsilon^-1
96 : !! %symchi
97 : !! %nomega=number of frequencies
98 : !! %npwe=number of plane waves for epsilon (input variable)
99 : !!
100 : !! OUTPUT
101 : !! (see side effects)
102 : !!
103 : !! SIDE EFFECTS
104 : !! chi0(npwe,npwe,nomega)=independent-particle susceptibility matrix in reciprocal space
105 : !!
106 : !! SOURCE
107 :
108 :
109 1251604 : subroutine assemblychi0_sym(is_metallic,ik_bz,nspinor,Ep,Ltg_q,green_w,npwepG0,rhotwg,Gsph_epsG0,chi0)
110 :
111 : !Arguments ------------------------------------
112 : !scalars
113 : logical,intent(in) :: is_metallic
114 : integer,intent(in) :: ik_bz,npwepG0,nspinor
115 : type(gsphere_t),intent(in) :: Gsph_epsG0
116 : type(littlegroup_t),intent(in) :: Ltg_q
117 : type(em1params_t),intent(in) :: Ep
118 : !arrays
119 : complex(gwp),intent(in) :: rhotwg(npwepG0)
120 : complex(dp),intent(in) :: green_w(Ep%nomega)
121 : complex(gwp),intent(inout) :: chi0(Ep%npwe*Ep%nI,Ep%npwe*Ep%nJ,Ep%nomega)
122 :
123 : !Local variables-------------------------------
124 : !scalars
125 : integer :: itim,io,isym,nthreads
126 : integer :: isymop,nsymop
127 : real(gwp) :: dr
128 : complex(gwp) :: dd
129 : !character(len=500) :: msg
130 : !arrays
131 2503208 : integer :: Sm1_gmG0(Ep%npwe)
132 1251604 : complex(gwp),allocatable :: rhotwg_sym(:,:)
133 : ! *************************************************************************
134 :
135 : ABI_UNUSED(nspinor)
136 :
137 : nthreads = xomp_get_max_threads()
138 :
139 1251604 : SELECT CASE (Ep%symchi)
140 : CASE (0)
141 : ! Do not use symmetries
142 :
143 : ! note that single precision is faster (sometimes factor ~2).
144 : ! Rely on MKL threads for OPENMP parallelization
145 :
146 1486152 : do io=1,Ep%nomega
147 : ! Check if green_w(io) is real (=> pure imaginary omega)
148 : ! and that it is not a metal
149 : ! then the corresponding chi0(io) is hermitian
150 1486152 : if( ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic ) then
151 375360 : dr=green_w(io)
152 375360 : call xher('U',Ep%npwe,dr,rhotwg,1,chi0(:,:,io),Ep%npwe)
153 : else
154 833092 : dd=green_w(io)
155 833092 : call xgerc(Ep%npwe,Ep%npwe,dd,rhotwg,1,rhotwg,1,chi0(:,:,io),Ep%npwe)
156 : endif
157 : end do
158 :
159 : CASE (1)
160 : ! Use symmetries to reconstruct the integrand in the BZ.
161 : !
162 : ! Notes on the symmetrization of the oscillator matrix elements
163 : ! If Sq = q then M_G^( Sk,q)= e^{-i(q+G).t} M_{ S^-1G} (k,q)
164 : ! If -Sq = q then M_G^(-Sk,q)= e^{-i(q+G).t} M_{-S^-1G}^*(k,q)
165 : !
166 : ! In case of an umklapp process
167 : ! If Sq = q+G0 then M_G( Sk,q)= e^{-i(q+G).t} M_{ S^-1(G-G0} (k,q)
168 : ! If -Sq = q+G0 then M_G(-Sk,q)= e^{-i(q+G).t} M_{-S^-1(G-G0)}^*(k,q)
169 : !
170 : ! Ltg_q%igmG0(ig,itim,isym) contains the index of G-G0 where ISq=q+G0
171 : ! Note that there is no need to take into account the phases due to q,
172 : ! They cancel in the scalar product ==> phmGt(G,isym)=e^{-iG\cdot t}
173 : !
174 : ! Mind the slicing of %rottbm1(npwepG0,timrev,nsym) and %phmGt(npwepG0,nsym) as
175 : ! these arrays, usually, do not conform to rho_twg_sym(npw) !
176 : !
177 : ! Loop over symmetries of the space group and time-reversal.
178 120945008 : nsymop = count(Ltg_q%wtksym(:,:,ik_bz)==1)
179 3895616 : ABI_MALLOC(rhotwg_sym,(Ep%npwe,nsymop))
180 973904 : isymop = 0
181 :
182 : ! Prepare all the rhotwg at once to use BLAS level 3 routines
183 40964272 : do isym=1,Ltg_q%nsym_sg
184 120797552 : do itim=1,Ltg_q%timrev
185 119823648 : if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
186 : ! This operation belongs to the little group and has to be used to reconstruct the BZ ===
187 : ! * In the following 3 lines mind the slicing (1:npwe)
188 : ! TODO this is a hot-spot, should add a test on the umklapp
189 : !
190 : !gmG0 => Ltg_q%igmG0(1:Ep%npwe,itim,isym)
191 308485202 : Sm1_gmG0(1:Ep%npwe) = Gsph_epsG0%rottbm1( Ltg_q%igmG0(1:Ep%npwe,itim,isym), itim,isym)
192 :
193 2118750 : isymop = isymop + 1
194 2049254 : SELECT CASE (itim)
195 : CASE (1)
196 151951656 : rhotwg_sym(1:Ep%npwe,isymop) = rhotwg(Sm1_gmG0) * Gsph_epsG0%phmGt(1:Ep%npwe,isym)
197 : CASE (2)
198 3350320 : rhotwg_sym(1:Ep%npwe,isymop) = GWPC_CONJG(rhotwg(Sm1_gmG0))*Gsph_epsG0%phmGt(1:Ep%npwe,isym)
199 : CASE DEFAULT
200 2118750 : ABI_BUG(sjoin('Wrong itim:', itoa(itim)))
201 : END SELECT
202 : end if
203 : end do
204 : end do
205 :
206 : ! Multiply rhotwg_sym by green_w(io) and accumulate in chi0(G,Gp,io)
207 : ! note that single precision is faster (sometimes factor ~2).
208 : ! Rely on MKL threads for OPENMP parallelization
209 4350154 : do io=1,Ep%nomega
210 : ! Check if green_w(io) is real (=> pure imaginary omega)
211 : ! and that it is not a metal
212 : ! then the corresponding chi0(io) is hermitian
213 4350154 : if( ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic ) then
214 2417886 : dr=green_w(io)
215 2417886 : call xherk('U','N',Ep%npwe,nsymop,dr,rhotwg_sym,Ep%npwe,one_gw,chi0(:,:,io),Ep%npwe)
216 : else
217 958364 : dd=green_w(io)
218 958364 : call xgemm('N','C',Ep%npwe,Ep%npwe,nsymop,dd,rhotwg_sym,Ep%npwe,rhotwg_sym,Ep%npwe,cone_gw,chi0(:,:,io),Ep%npwe)
219 : endif
220 : end do
221 :
222 973904 : ABI_FREE(rhotwg_sym)
223 :
224 : CASE DEFAULT
225 1251604 : ABI_BUG(sjoin('Wrong symchi:', itoa(Ep%symchi)))
226 : END SELECT
227 :
228 1251604 : end subroutine assemblychi0_sym
229 : !!***
230 :
231 : !----------------------------------------------------------------------
232 :
233 : !!****f* m_chi0tk/mkrhotwg_sigma
234 : !! NAME
235 : !! mkrhotwg_sigma
236 : !!
237 : !! FUNCTION
238 : !! Helper function used to calculate selected linear combination
239 : !! of the oscillator matrix elements in the case of non-collinear magnetism.
240 : !!
241 : !! INPUTS
242 : !! ii=Index selecting the particolar combination of spin components.
243 : !! npw=Number of plane-waves in the oscillators.
244 : !! nspinor=Number of spinorial components.
245 : !! rhotwg(npw*nspinor**2)=OScillator matrix elements.
246 : !!
247 : !! OUTPUT
248 : !! rhotwg_I(npw)=Required linear combination of the oscillator matrix elements.
249 : !!
250 : !! SOURCE
251 :
252 : subroutine mkrhotwg_sigma(ii,nspinor,npw,rhotwg,rhotwg_I)
253 :
254 : !Arguments ------------------------------------
255 : !scalars
256 : integer,intent(in) :: ii,npw,nspinor
257 : !arrays
258 : complex(gwp),intent(in) :: rhotwg(npw*nspinor**2)
259 : complex(gwp),intent(out) :: rhotwg_I(npw)
260 : ! *************************************************************************
261 :
262 : SELECT CASE (ii)
263 : CASE (1)
264 : ! $ M_0 = M_{\up,\up} + M_{\down,\down} $
265 : rhotwg_I(:) = rhotwg(1:npw) + rhotwg(npw+1:2*npw)
266 : CASE (2)
267 : ! $ M_z = M_{\up,\up} - M_{\down,\down} $
268 : rhotwg_I(:) = rhotwg(1:npw) - rhotwg(npw+1:2*npw)
269 : CASE (3)
270 : ! $ M_x = M_{\up,\down} + M_{\down,\up} $
271 : rhotwg_I(:) = ( rhotwg(2*npw+1:3*npw) + rhotwg(3*npw+1:4*npw) )
272 : CASE (4)
273 : ! $ M_y = i * (M_{\up,\down} -M_{\down,\up}) $
274 : rhotwg_I(:) = (rhotwg(2*npw+1:3*npw) - rhotwg(3*npw+1:4*npw) )*j_gw
275 : CASE DEFAULT
276 : ABI_BUG(sjoin('Wrong ii value:', itoa(ii)))
277 : END SELECT
278 :
279 : end subroutine mkrhotwg_sigma
280 : !!***
281 :
282 : !----------------------------------------------------------------------
283 :
284 : !!****f* m_chi0tk/symmetrize_afm_chi0tk
285 : !! NAME
286 : !! symmetrize_afm_chi0
287 : !!
288 : !! FUNCTION
289 : !! Reconstruct the (down, down) component of the irreducible polarizability
290 : !! starting from the (up,up) element in case of systems with AFM symmetries
291 : !! (i.e nspden==2 and nsppol=1). Return the trace (up,up)+(down,down) of the
292 : !! matrix as required by GW calculations.
293 : !!
294 : !! INPUTS
295 : !! Cryst<crystal_t>= Information on symmetries and unit cell.
296 : !! Gsph<gsphere_t>= The G-sphere used to descrive chi0.
297 : !! npwe=Number of G-vectors in chi0.
298 : !! nomega=number of frequencies.
299 : !! Ltg_q<littlegroup_t>=Structure with useful table describing the little group of the q-point.
300 : !!
301 : !! SIDE EFFECTS
302 : !! chi0(npwe,npwe,nomega)= In input the up-up component, in output the trace of chi0.
303 : !! The value of matrix elements that should be zero due to AFM symmetry properties are
304 : !! forced to be zero (see NOTES below).
305 : !! [chi0_lwing(npwe,nomega,3)] = Lower wings, symmetrized in output.
306 : !! [chi0_uwing(npwe,nomega,3)] = Upper wings, symmetrized in output.
307 : !! [chi0_head(3,3,nomega) ] = Head of chi0, symmetrized in output.
308 : !!
309 : !! NOTES
310 : !! In the case of magnetic group Shubnikov type III:
311 : !! For each set of paired FM-AFM symmetries, the down-down component of
312 : !! a generic response function in reciprocal space can be obtained according to:
313 : !!
314 : !! chi^{down,down}_{G1,G2}(q) = chi^{up up}_{G1,G2}(q) e^{iS(G1-G2).(tnonsFM - tnonsAFM)}
315 : !!
316 : !! where S is the rotational part common to the FM-AFM pair, tnonsFM and tnonsAFM
317 : !! are the fractional translations associated to the ferromagnetic and antiferromagnetic symmetry, respectively.
318 : !! Note that, if for a given G1-G2 pair, the phase e^{iS(G1-G2).(tnonsFM - tnonsAFM) depends
319 : !! on the FM-AFM symmetry pair, then the corresponding matrix element of chi0 must be zero.
320 : !! Actually this is manually enforced in the code because this property might not be
321 : !! perfectly satisfied due to round-off errors.
322 : !!
323 : !! In the case of magnetic group Shubnikov type III:
324 : !! Only the AFM symmetries that preserve the external q-point (with or without time-reversal)
325 : !! are used to get the (down, down) component using the fact that:
326 : !!
327 : !! chi^{down,down}_{G1,G2}(Sq) = chi^{up up}_{S^{-1}G1,S^{-1}G2}(q) e^{i(G2-G1).tnons_S }
328 : !!
329 : !! Actually we perform an average over subset of the little group of q with AFM character in
330 : !! order to reduce as much as possible errors due to round off errors. In brief we evaluate:
331 : !!
332 : !! 1/N_{Ltq} \sum_{S\in Ltg AFM} chi^{up up}_{S^{-1}G1,S^{-1}G2}(q) e^{i(G2-G1).tnons_S }
333 : !!
334 : !! where N_{Ltg} is the number of AFM operation in the little group (time reversal included)
335 : !!
336 : !! TODO
337 : !! It is possible to symmetrize chi0 without any the extra allocation for afm_mat.
338 : !! More CPU demanding but safer in case of a large chi0 matrix. One might loop over G1 and G2 shells ...
339 : !!
340 : !! SOURCE
341 :
342 13 : subroutine symmetrize_afm_chi0(Cryst,Gsph,Ltg_q,npwe,nomega,chi0,chi0_head,chi0_lwing,chi0_uwing)
343 :
344 : !Arguments ------------------------------------
345 : !scalars
346 : integer,intent(in) :: npwe,nomega
347 : type(gsphere_t),intent(in) :: Gsph
348 : type(crystal_t),intent(in) :: Cryst
349 : type(littlegroup_t),intent(in) :: Ltg_q
350 : !arrays
351 : complex(gwp),optional,intent(inout) :: chi0(npwe,npwe,nomega)
352 : complex(dp),optional,intent(inout) :: chi0_lwing(npwe,nomega,3)
353 : complex(dp),optional,intent(inout) :: chi0_uwing(npwe,nomega,3)
354 : complex(dp),optional,intent(inout) :: chi0_head(3,3,nomega)
355 :
356 : !Local variables ------------------------------
357 : !scalars
358 : integer :: io,ig1,ig2,isymf,isyma,isym,ipair,k0g,kg,npairs,nonzero
359 : integer :: iSmg1,iSmg2,itim,shubnikov,ntest
360 : complex(gwp) :: phase,phase_old,sumchi,ctmp
361 : logical :: found
362 : !character(len=500) :: msg
363 : !arrays
364 20 : integer :: rotfm(3,3),rotafm(3,3),pairs2sym(2,Cryst%nsym/2)
365 : real(dp) :: tfm(3),tafm(3)
366 10 : complex(gwp),allocatable :: afm_mat(:),chi0_afm(:,:)
367 : !************************************************************************
368 :
369 20 : ABI_CHECK(ANY(Cryst%symafm==-1),'Not magnetic space group')
370 : !
371 : ! ==== Find shubnikov type ====
372 : ! TODO This info should be stored in Cryst%
373 10 : shubnikov=4; npairs=0
374 :
375 970 : do isymf=1,Cryst%nsym
376 960 : if (Cryst%symafm(isymf)==-1) CYCLE
377 6240 : rotfm = Cryst%symrec(:,:,isymf)
378 : tfm = Cryst%tnons(:,isymf)
379 : found = .FALSE.
380 :
381 46560 : do isyma=1,Cryst%nsym
382 46080 : if (Cryst%symafm(isyma)==1) CYCLE
383 299520 : rotafm = Cryst%symrec(:,:,isyma)
384 :
385 55680 : if (ALL(rotfm==rotafm)) then
386 480 : found=.TRUE.
387 : tafm = Cryst%tnons(:,isyma)
388 480 : npairs=npairs+1
389 480 : ABI_CHECK(npairs<=Cryst%nsym/2,'Wrong AFM group')
390 480 : pairs2sym(1,npairs)=isymf
391 480 : pairs2sym(2,npairs)=isyma
392 : end if
393 : end do !isyma
394 :
395 490 : if (.not.found) then
396 0 : shubnikov=3; EXIT !isymf
397 : end if
398 : end do !isymf
399 :
400 10 : select case (shubnikov)
401 :
402 : case (4)
403 10 : call wrtout(std_out,' Found Magnetic group Shubnikov type IV')
404 10 : ABI_CHECK(npairs==Cryst%nsym/2,'Wrong AFM space group')
405 :
406 30 : ABI_MALLOC(afm_mat,(npwe*(npwe+1)/2))
407 :
408 : ! jmb
409 820 : phase_old=zero
410 :
411 820 : do ig2=1,npwe
412 810 : k0g=ig2*(ig2-1)/2
413 34030 : do ig1=1,ig2
414 33210 : kg=k0g+ig1
415 33210 : nonzero=1
416 :
417 1627290 : do ipair=1,Cryst%nsym/2
418 1594080 : isymf = pairs2sym(1,ipair)
419 1594080 : isyma = pairs2sym(2,ipair)
420 : phase = ( Gsph%phmSGt(ig1,isymf)*conjg(Gsph%phmSGt(ig1,isyma)) ) * &
421 1594080 : & ( Gsph%phmSGt(ig2,isymf)*conjg(Gsph%phmSGt(ig2,isyma)) )
422 1594080 : if (ipair>1 .and. (ABS(phase_old-phase) > tol6)) then
423 : nonzero=0; EXIT
424 : end if
425 1627290 : phase_old=phase
426 : end do !ipair
427 :
428 34020 : afm_mat(kg)=nonzero*(cone_gw + phase)
429 : end do !ig1
430 : end do !ig2
431 : !
432 : ! =======================================================================
433 : ! ==== Symmetrize chi0 constructing chi0^{\up\up} + chi0{\down\down} ====
434 : ! =======================================================================
435 : !
436 : ! head^{\down\down} = head^{\up\up}
437 36 : if (PRESENT(chi0_head)) chi0_head = two * chi0_head
438 : !
439 : ! w^{\down\down}_{0 G'} = e^{-iSG'.(tFM-tAFM)} w^{\up\up}_{0 G'}.
440 : ! w^{\down\down}_{G 0 } = e^{+iSG .(tFM-tAFM)} w^{\up\up}_{G 0 }.
441 10 : if (PRESENT(chi0_uwing)) then
442 3 : do io=1,nomega
443 165 : do ig2=1,npwe
444 162 : k0g=ig2*(ig2-1)/2
445 162 : kg=k0g+1
446 650 : chi0_uwing(ig2,io,:)=afm_mat(kg)*chi0_uwing(ig2,io,:)
447 : end do
448 : end do
449 : end if
450 :
451 10 : if (PRESENT(chi0_lwing)) then
452 3 : do io=1,nomega
453 165 : do ig1=1,npwe
454 162 : k0g=ig1*(ig1-1)/2
455 162 : kg=k0g+1
456 650 : chi0_lwing(ig1,io,:)=conjg(afm_mat(kg))*chi0_lwing(ig1,io,:)
457 : end do
458 : end do
459 : end if
460 :
461 10 : if (PRESENT(chi0)) then
462 30 : do io=1,nomega
463 : ! Take care of diagonal.
464 1640 : do ig1=1,npwe
465 1640 : chi0(ig1,ig1,io)=two*chi0(ig1,ig1,io)
466 : end do
467 :
468 : ! Upper and lower triangle are treated differently:
469 : ! We took advantage of the fact the afm_mat is hermitian to reduce memory.
470 1620 : do ig2=2,npwe
471 1600 : k0g=ig2*(ig2-1)/2
472 66420 : do ig1=1,ig2-1
473 64800 : kg=k0g+ig1
474 66400 : chi0(ig1,ig2,io)=afm_mat(kg)*chi0(ig1,ig2,io)
475 : end do
476 : end do
477 :
478 1630 : do ig1=2,npwe
479 1600 : k0g=ig1*(ig1-1)/2
480 66420 : do ig2=1,ig1-1
481 64800 : kg=k0g+ig2
482 66400 : chi0(ig1,ig2,io)=conjg(afm_mat(kg))*chi0(ig1,ig2,io)
483 : end do
484 : end do
485 : end do !io
486 : end if
487 :
488 10 : ABI_FREE(afm_mat)
489 :
490 : case (3)
491 0 : call wrtout(std_out,' Found Magnetic group Shubnikov type III')
492 0 : ABI_ERROR('Shubnikov type III not implemented')
493 :
494 0 : ntest=0
495 0 : do itim=1,ltg_q%timrev
496 0 : do isym=1,ltg_q%nsym_sg
497 : ! use only afm sym preserving q with and without time-reversal
498 0 : if ( cryst%symafm(isym)==-1 .and. ltg_q%preserve(itim,isym)==1 ) ntest=ntest+1
499 : end do
500 : end do
501 :
502 0 : if (ntest==0) then
503 0 : ABI_WARNING("no symmetry can be used!")
504 : end if
505 : !RETURN
506 0 : ABI_MALLOC(chi0_afm,(npwe,npwe))
507 :
508 0 : do io=1,nomega
509 :
510 0 : do ig2=1,npwe
511 0 : do ig1=1,npwe
512 : sumchi=czero_gw
513 :
514 0 : do itim=1,ltg_q%timrev
515 0 : do isym=1,ltg_q%nsym_sg
516 : ! use only afm sym preserving q with and without time-reversal
517 0 : if ( cryst%symafm(isym)==-1 .and. ltg_q%preserve(itim,isym)==1 ) then
518 0 : phase = Gsph%phmGt(ig1,isym)*conjg(Gsph%phmGt(ig2,isym))
519 0 : iSmg1=Gsph%rottbm1(ig1,itim,isym)
520 0 : iSmg2=Gsph%rottbm1(ig2,itim,isym)
521 0 : ctmp=chi0(iSmg1,iSmg2,io)*phase !; if (itim==2) ctmp=conjg(ctmp) !check this
522 0 : sumchi=sumchi+ctmp !chi0(iSmg1,iSmg2,io)*phase
523 : end if
524 : end do ! isym
525 : end do !itim
526 :
527 0 : chi0_afm(ig1,ig2)=sumchi/Ltg_q%nsym_ltg !has to be changed in case of time-reversal
528 : end do !ig1
529 : end do !ig2
530 :
531 : ! We want chi_{up,up} +chi_{dwn,dwn}.
532 0 : chi0(:,:,io)=chi0(:,:,io)+chi0_afm(:,:)
533 : end do !iomega
534 :
535 0 : ABI_FREE(chi0_afm)
536 :
537 : case default
538 10 : ABI_BUG(sjoin('Wrong value for shubnikov= ', itoa(shubnikov)))
539 : end select
540 :
541 13 : end subroutine symmetrize_afm_chi0
542 : !!***
543 :
544 : !----------------------------------------------------------------------
545 :
546 : !!****f* m_chi0tk/accumulate_chi0_q0
547 : !! NAME
548 : !! accumulate_chi0_q0
549 : !!
550 : !! FUNCTION
551 : !! Update the independent particle susceptibility at q==0 for the contribution
552 : !! of one pair of occupied-unoccupied band, for each frequency.
553 : !! This routine takes advantage of the symmetries of the little group of the external q-point
554 : !! to symmetrize the contribution arising from the input k-point located in the IBZ_q.
555 : !! It computes:
556 : !!
557 : !! $ \chi_0(G1,G2,io) = \chi_0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))*green_w(io) $
558 : !!
559 : !! where S is a symmetry in reciprocal space.
560 : !! The matrix elements of the gradient operator and [V_{nl},r] are symmetrized as well.
561 : !!
562 : !! INPUTS
563 : !! ik_bz=Index of the k-point whose contribution has to be added to chi0.
564 : !! isym_kbz=Index of the symmetry such that k_bz = IS k_ibz
565 : !! itim_kbz=2 if time-reversal has to be used to obtain k_bz, 1 otherwise.
566 : !! npwepG0=Maximum number of G vectors
567 : !! rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states.
568 : !! rhotwx(3,nspinor**2)=Matrix element of the operator $-i[H,r]/(e1-e2) = -i r$ in reciprocal lattice units.
569 : !! green_w(nomega)=Frequency dependent part of the Green function.
570 : !! Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point.
571 : !! %timrev=2 it time-reversal is used, 1 otherwise.
572 : !! %nsym_sg=Number of space group symmetries.
573 : !! %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point.
574 : !! Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors
575 : !! %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg.
576 : !! %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion.
577 : !! %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations.
578 : !! Cryst<crystal_t>=Structure defining the unit cell and its symmetries
579 : !! %nsym=Number of symmetries.
580 : !! %symrec(3,3,nsym)=Symmetry operation in reciprocal space (reduced coordinates)
581 : !! Ep<em1params_t>=Parameters of the chi0 calculation.
582 : !! %npwe=number of plane waves in chi0.
583 : !! %symchi=1 if symmetrization has to be performed.
584 : !! %nomega=number of frequencies in chi0.
585 : !!
586 : !! SIDE EFFECTS
587 : !! chi0(npwe,npwe,nomega)= Updated independent-particle susceptibility matrix in reciprocal space at q==0.
588 : !! chi0_head(3,3,Ep%nomega)=Head.
589 : !! chi0_lwing(Ep%npwe*Ep%nI,Ep%nomega,3)=Lower wing.
590 : !! chi0_uwing(Ep%npwe*Ep%nJ,Ep%nomega,3)=Upper wing.
591 : !!
592 : !! NOTES
593 : !!
594 : !! 1) Symmetrization of the oscilator matrix elements.
595 : !! If Sq = q then M_G( Sk,q)= e^{-i(q+G).\tau} M_{ S^-1G} (k,q)
596 : !! If -Sq = q then M_G(-Sk,q)= e^{-i(q+G).\tau} M_{-S^-1G}^*(k,q)
597 : !!
598 : !! In the case of umklapps:
599 : !! If Sq = q+G0 then M_G( Sk,q)= e^{-i(q+G).\tau} M_{ S^-1(G-G0} (k,q)
600 : !! If -Sq = q+G0 then M_G(-Sk,q)= e^{-i(q+G).\tau} M_{-S^-1(G-G0)}^*(k,q)
601 : !!
602 : !! In the equation below there is no need to take into account the phases due to q.t
603 : !! as they cancel each other in the scalar product ==> only phmGt(G,isym)=e^{-iG.\tau} is needed.
604 : !!
605 : !! 2) Symmetrization of the matrix elements of the position operator.
606 : !!
607 : !! <Sk,b|\vec r| Sk,b'> = <k b| R\vec r + \tau|k b'>
608 : !!
609 : !! where S is one of the symrec operation, R and \tau is the corresponding
610 : !! operation in real space. The term involving the fractional translation is zero provided that b /= b'.
611 : !!
612 : !! SOURCE
613 :
614 103914 : subroutine accumulate_chi0_q0(is_metallic,ik_bz,isym_kbz,itim_kbz,gwcomp,nspinor,npwepG0,Ep,Cryst,Ltg_q,Gsph_epsG0,&
615 103914 : chi0,rhotwx,rhotwg,green_w,green_enhigh_w,deltaf_b1b2,chi0_head,chi0_lwing,chi0_uwing)
616 :
617 : !Arguments ------------------------------------
618 : !scalars
619 : logical,intent(in) :: is_metallic
620 : integer,intent(in) :: ik_bz,isym_kbz,itim_kbz,npwepG0,nspinor,gwcomp
621 : real(dp),intent(in) :: deltaf_b1b2
622 : type(littlegroup_t),intent(in) :: Ltg_q
623 : type(gsphere_t),target,intent(in) :: Gsph_epsG0
624 : type(crystal_t),intent(in) :: Cryst
625 : type(em1params_t),intent(in) :: Ep
626 : !arrays
627 : complex(gwp),intent(in) :: rhotwg(npwepG0)
628 : complex(gwp),intent(in) :: rhotwx(3, nspinor**2)
629 : complex(gwp),intent(inout) :: chi0(Ep%npwe*Ep%nI, Ep%npwe*Ep%nJ, Ep%nomega)
630 : complex(dp),intent(in) :: green_w(Ep%nomega), green_enhigh_w(Ep%nomega)
631 : complex(dp),intent(inout) :: chi0_head(3, 3, Ep%nomega)
632 : complex(dp),intent(inout) :: chi0_lwing(Ep%npwe*Ep%nI, Ep%nomega, 3)
633 : complex(dp),intent(inout) :: chi0_uwing(Ep%npwe*Ep%nJ, Ep%nomega, 3)
634 :
635 : !Local variables-------------------------------
636 : !scalars
637 : integer :: itim,io,isym,idir,jdir,isymop,nsymop,npwe,nomega
638 : real(gwp) :: dr
639 : complex(gwp) :: dd
640 : !character(len=500) :: msg
641 : !arrays
642 103914 : integer,contiguous, pointer :: Sm1G(:)
643 : complex(dp) :: mir_kbz(3)
644 103914 : complex(gwp),allocatable :: rhotwg_sym(:,:)
645 103914 : complex(gwp), contiguous, pointer :: phmGt(:)
646 : !************************************************************************
647 :
648 : ABI_UNUSED(deltaf_b1b2)
649 :
650 103914 : npwe = ep%npwe; nomega = ep%nomega
651 :
652 103914 : select case (Ep%symchi)
653 : case (0)
654 : ! Do not use symmetries.
655 : ! Symmetrize rhotwg in the full BZ and accumulate over the full BZ i.e.
656 : ! chi0(G1,G2,io) = chi0(G1,G2,io) + (rhotwg(G1)*conjg(rhotwg(G2)))*green_w(io)
657 : !
658 : ! The non-analytic term is symmetrized for this k-point in the BZ according to:
659 : ! rhotwg(1) = S^-1q * rhotwx_ibz
660 : ! rhotwg(1) = -S^-1q * conjg(rhotwx_ibz) if time-reversal is used.
661 :
662 : ! Multiply elements G1,G2 of rhotwg by green_w(io) and accumulate in chi0(G1,G2,io)
663 :
664 : !$OMP PARALLEL DO PRIVATE(dr, dd) IF (nomega > 2)
665 238676 : do io=1,nomega
666 : ! Check if green_w(io) is real (=> pure imaginary omega)
667 : ! and that it is not a metal
668 : ! then the corresponding chi0(io) is hermitian
669 238676 : if (ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic) then
670 71488 : dr = green_w(io)
671 71488 : call xher('U', npwe, dr, rhotwg, 1, chi0(:,:,io), npwe)
672 : else
673 124514 : dd = green_w(io)
674 124514 : call xgerc(npwe, npwe, dd, rhotwg, 1, rhotwg, 1, chi0(:,:,io), npwe)
675 : endif
676 : end do
677 :
678 : ! === Accumulate heads and wings for each small q ===
679 : ! FIXME extrapolar method should be checked!!
680 : ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
681 42674 : if (nspinor == 1) then
682 1194872 : mir_kbz = (3-2*itim_kbz) * matmul(Cryst%symrec(:,:,isym_kbz), rhotwx(:,1))
683 : else
684 0 : mir_kbz = (3-2*itim_kbz) * matmul(Cryst%symrec(:,:,isym_kbz), sum(rhotwx(:,1:2), dim=2))
685 : end if
686 60818 : if (itim_kbz == 2) mir_kbz = conjg(mir_kbz)
687 :
688 : ! here we might take advantage of Hermiticity along Im axis in RPA (see mkG0w)
689 170696 : do idir=1,3
690 758702 : do io=1,nomega
691 30677160 : chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir)*conjg(rhotwg)
692 30677160 : chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg*conjg(mir_kbz(idir))
693 716028 : if (gwcomp == 1) then
694 : ! Add contribution due to extrapolar technique.
695 5564160 : chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_enhigh_w(io) * mir_kbz(idir)*conjg(rhotwg)
696 5564160 : chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_enhigh_w(io) * rhotwg*conjg(mir_kbz(idir))
697 : end if
698 : end do
699 : end do
700 :
701 : ! Accumulate the head.
702 238676 : do io=1,nomega
703 826682 : do jdir=1,3
704 2548026 : do idir=1,3
705 1764018 : chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_w(io) * mir_kbz(idir)*conjg(mir_kbz(jdir))
706 2352024 : if (gwcomp == 1) then
707 : ! Add contribution due to extrapolar technique.
708 120960 : chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_enhigh_w(io) * mir_kbz(idir)*conjg(mir_kbz(jdir))
709 : end if
710 : end do
711 : end do
712 : end do
713 :
714 : case (1)
715 : ! Use symmetries to reconstruct the integrand.
716 :
717 7140952 : nsymop = count(Ltg_q%wtksym(:,:,ik_bz) == 1)
718 244960 : ABI_MALLOC(rhotwg_sym, (npwe, nsymop))
719 61240 : isymop = 0
720 :
721 : ! Loop over symmetries of the space group and time-reversal.
722 2421144 : do isym=1,Ltg_q%nsym_sg
723 7085656 : do itim=1,Ltg_q%timrev
724 :
725 7024416 : if (Ltg_q%wtksym(itim, isym, ik_bz) == 1) then
726 : ! This operation belongs to the little group and has to be considered to reconstruct the BZ.
727 322053 : phmGt => Gsph_epsG0%phmGt (1:npwe, isym) ! In the 2 lines below note the slicing (1:npwe)
728 322053 : Sm1G => Gsph_epsG0%rottbm1(1:npwe, itim, isym)
729 :
730 322053 : isymop = isymop + 1
731 310245 : select case (itim)
732 : case (1)
733 25530724 : rhotwg_sym(1:npwe, isymop) = rhotwg(Sm1G(1:npwe)) * phmGt(1:npwe)
734 : case (2)
735 573888 : rhotwg_sym(1:npwe,isymop) = conjg(rhotwg(Sm1G(1:npwe))) * phmGt(1:npwe)
736 : case default
737 322053 : ABI_BUG(sjoin('Wrong value of itim:', itoa(itim)))
738 : end select
739 :
740 : ! === Accumulate heads and wings for each small q ===
741 : ! FIXME extrapolar method should be checked!!
742 :
743 : ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
744 322053 : if (nspinor == 1) then
745 8931468 : mir_kbz = (3-2*itim) * matmul(Cryst%symrec(:,:,isym), rhotwx(:,1))
746 : else
747 113664 : mir_kbz = (3-2*itim) * matmul(Cryst%symrec(:,:,isym), sum(rhotwx(:,1:2), dim=2))
748 : end if
749 :
750 357477 : if (itim == 2) mir_kbz = conjg(mir_kbz)
751 :
752 : ! here we might take advantage of Hermiticity along Im axis in RPA (see mkG0w)
753 1288212 : do idir=1,3
754 5132814 : do io=1,nomega
755 232708680 : chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir) * conjg(rhotwg_sym(:,isymop))
756 232708680 : chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg_sym(:,isymop) * conjg(mir_kbz(idir))
757 4810761 : if (gwcomp == 1) then
758 : ! Add contribution due to extrapolar technique.
759 4095936 : chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_enhigh_w(io) * mir_kbz(idir) * conjg(rhotwg_sym(:,isymop))
760 4095936 : chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_enhigh_w(io) * rhotwg_sym(:,isymop) * conjg(mir_kbz(idir))
761 : end if
762 : end do
763 : end do
764 :
765 : ! Accumulate the head.
766 1603587 : do io=1,nomega
767 5448189 : do jdir=1,3
768 16659942 : do idir=1,3
769 11533806 : chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_w(io) * mir_kbz(idir) * conjg(mir_kbz(jdir))
770 15378408 : if (gwcomp == 1) then
771 : ! Add contribution due to extrapolar technique.
772 194256 : chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_enhigh_w(io)*mir_kbz(idir) * conjg(mir_kbz(jdir))
773 : end if
774 : end do
775 : end do
776 : end do
777 :
778 : end if !wtksym
779 : end do !itim
780 : end do !isym
781 :
782 : ! Multiply rhotwg_sym by green_w(io) and accumulate in chi0(G,Gp,io)
783 : !$OMP PARALLEL DO PRIVATE(dr, dd) IF (nomega > 2)
784 288635 : do io=1,nomega
785 : ! Check if green_w(io) is real (=> pure imaginary omega) and that it is not a metal
786 : ! then the corresponding chi0(io) is hermitian
787 288635 : if (ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic) then
788 160675 : dr = green_w(io)
789 160675 : call xherk('U', 'N', npwe, nsymop, dr, rhotwg_sym, npwe, one_gw, chi0(:,:,io), npwe)
790 : else
791 66720 : dd = green_w(io)
792 66720 : call xgemm('N', 'C', npwe, npwe, nsymop, dd, rhotwg_sym, npwe, rhotwg_sym, npwe, cone_gw, chi0(:,:,io), npwe)
793 : endif
794 : end do
795 :
796 61240 : ABI_FREE(rhotwg_sym)
797 :
798 : case default
799 0 : ABI_BUG(sjoin('Wrong value of symchi:', itoa(Ep%symchi)))
800 : end select
801 :
802 103914 : end subroutine accumulate_chi0_q0
803 : !!***
804 :
805 : !----------------------------------------------------------------------
806 :
807 : !!****f* m_chi0tk/accumulate_head_wings_imagw
808 : !! NAME
809 : !! accumulate_head_wings_imagw
810 : !!
811 : !! FUNCTION
812 :
813 0 : subroutine accumulate_head_wings_imagw( &
814 : npwe, nomega, nI, nJ, symchi, &
815 : is_metallic, ik_bz, isym_kbz, itim_kbz, nspinor, Cryst, Ltg_q, Gsph_epsG0,&
816 0 : rhotwx, rhotwg, green_w, chi0_head, chi0_lwing, chi0_uwing)
817 :
818 : !Arguments ------------------------------------
819 : !scalars
820 : integer,intent(in) :: npwe, nomega, nI, nJ, symchi
821 : logical,intent(in) :: is_metallic
822 : integer,intent(in) :: ik_bz,isym_kbz,itim_kbz,nspinor
823 : type(littlegroup_t),intent(in) :: Ltg_q
824 : type(gsphere_t),target,intent(in) :: Gsph_epsG0
825 : type(crystal_t),intent(in) :: Cryst
826 : !arrays
827 : complex(gwp),intent(in) :: rhotwg(npwe)
828 : complex(gwp),intent(in) :: rhotwx(3, nspinor**2)
829 : complex(dp),intent(in) :: green_w(nomega)
830 : complex(dp),intent(inout) :: chi0_head(3, 3, nomega)
831 : complex(dp),intent(inout) :: chi0_lwing(npwe*nI, nomega, 3)
832 : complex(dp),intent(inout) :: chi0_uwing(npwe*nJ, nomega, 3)
833 :
834 : !Local variables-------------------------------
835 : !scalars
836 : integer :: itim,io,isym,idir,jdir,isymop,nsymop
837 : !character(len=500) :: msg
838 : !arrays
839 0 : integer,contiguous, pointer :: Sm1G(:)
840 : complex(dp) :: mir_kbz(3)
841 0 : complex(gwp),allocatable :: rhotwg_sym(:,:)
842 0 : complex(gwp), contiguous, pointer :: phmGt(:)
843 : !************************************************************************
844 :
845 : ABI_UNUSED([is_metallic])
846 :
847 0 : select case (symchi)
848 : case (0)
849 : ! Do not use symmetries.
850 : ! Symmetrize rhotwg in the full BZ and accumulate over the full BZ i.e.
851 : ! chi0(G1,G2,io) = chi0(G1,G2,io) + (rhotwg(G1)*conjg(rhotwg(G2)))*green_w(io)
852 : !
853 : ! The non-analytic term is symmetrized for this k-point in the BZ according to:
854 : ! rhotwg(1) = S^-1q * rhotwx_ibz
855 : ! rhotwg(1) = -S^-1q * conjg(rhotwx_ibz) if time-reversal is used.
856 :
857 : ! Multiply elements G1,G2 of rhotwg by green_w(io) and accumulate in chi0(G1,G2,io)
858 :
859 : ! === Accumulate heads and wings for each small q ===
860 : ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
861 0 : if (nspinor == 1) then
862 0 : mir_kbz = (3-2*itim_kbz) * matmul(Cryst%symrec(:,:,isym_kbz), rhotwx(:,1))
863 : else
864 0 : mir_kbz = (3-2*itim_kbz) * matmul(Cryst%symrec(:,:,isym_kbz), sum(rhotwx(:,1:2), dim=2))
865 : end if
866 0 : if (itim_kbz == 2) mir_kbz = conjg(mir_kbz)
867 :
868 : ! here we might take advantage of Hermiticity along Im axis in RPA (see mkG0w)
869 0 : do idir=1,3
870 0 : do io=1,nomega
871 0 : chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir)*conjg(rhotwg)
872 0 : chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg*conjg(mir_kbz(idir))
873 : end do
874 : end do
875 :
876 : ! Accumulate the head.
877 0 : do io=1,nomega
878 0 : do jdir=1,3
879 0 : do idir=1,3
880 0 : chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_w(io) * mir_kbz(idir)*conjg(mir_kbz(jdir))
881 : end do
882 : end do
883 : end do
884 :
885 : case (1)
886 : ! Use symmetries to reconstruct the integrand.
887 :
888 0 : nsymop = count(Ltg_q%wtksym(:,:,ik_bz) == 1)
889 0 : ABI_MALLOC(rhotwg_sym, (npwe, nsymop))
890 0 : isymop = 0
891 :
892 : ! Loop over symmetries of the space group and time-reversal.
893 0 : do isym=1,Ltg_q%nsym_sg
894 0 : do itim=1,Ltg_q%timrev
895 :
896 0 : if (Ltg_q%wtksym(itim, isym, ik_bz) == 1) then
897 : ! This operation belongs to the little group and has to be considered to reconstruct the BZ.
898 0 : phmGt => Gsph_epsG0%phmGt (1:npwe, isym) ! In the 2 lines below note the slicing (1:npwe)
899 0 : Sm1G => Gsph_epsG0%rottbm1(1:npwe, itim, isym)
900 :
901 0 : isymop = isymop + 1
902 0 : select case (itim)
903 : case (1)
904 0 : rhotwg_sym(1:npwe, isymop) = rhotwg(Sm1G(1:npwe)) * phmGt(1:npwe)
905 : case (2)
906 0 : rhotwg_sym(1:npwe,isymop) = conjg(rhotwg(Sm1G(1:npwe))) * phmGt(1:npwe)
907 : case default
908 0 : ABI_BUG(sjoin('Wrong value of itim:', itoa(itim)))
909 : end select
910 :
911 : ! === Accumulate heads and wings for each small q ===
912 : ! FIXME extrapolar method should be checked!!
913 :
914 : ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
915 0 : if (nspinor == 1) then
916 0 : mir_kbz = (3-2*itim) * matmul(Cryst%symrec(:,:,isym), rhotwx(:,1))
917 : else
918 0 : mir_kbz = (3-2*itim) * matmul(Cryst%symrec(:,:,isym), sum(rhotwx(:,1:2), dim=2))
919 : end if
920 :
921 0 : if (itim == 2) mir_kbz = conjg(mir_kbz)
922 :
923 : ! here we might take advantage of Hermiticity along Im axis in RPA (see mkG0w)
924 0 : do idir=1,3
925 0 : do io=1,nomega
926 0 : chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir) * conjg(rhotwg_sym(:,isymop))
927 0 : chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg_sym(:,isymop) * conjg(mir_kbz(idir))
928 : end do
929 : end do
930 :
931 : ! Accumulate the head.
932 0 : do io=1,nomega
933 0 : do jdir=1,3
934 0 : do idir=1,3
935 0 : chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_w(io) * mir_kbz(idir) * conjg(mir_kbz(jdir))
936 : end do
937 : end do
938 : end do
939 :
940 : end if !wtksym
941 : end do !itim
942 : end do !isym
943 :
944 0 : ABI_FREE(rhotwg_sym)
945 :
946 : case default
947 0 : ABI_BUG(sjoin('Wrong value of symchi:', itoa(symchi)))
948 : end select
949 :
950 0 : end subroutine accumulate_head_wings_imagw
951 : !!***
952 :
953 : !----------------------------------------------------------------------
954 :
955 : !!****f* m_chi0tk/accumulate_sfchi0_q0
956 : !! NAME
957 : !! accumulate_sfchi0_q0
958 : !!
959 : !! FUNCTION
960 : !! Update the spectral function of the independent particle susceptibility at q==0 for the contribution
961 : !! of one pair of occupied-unoccupied band, for each frequency.
962 : !! If symchi==1, the symmetries belonging to the little group of the external point q are used
963 : !! to reconstruct the contributions in the full Brillouin zone. In this case, the equation implented is:
964 : !!
965 : !! $ chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))* \delta(\omega -trans) $
966 : !!
967 : !! where S is a symmetry belonging to the little group of q.
968 : !! The subroutine also performs the symmetrization of the matrix elements of the
969 : !! gradient operator and of the commutator [V_{nl},r] with the position operator.
970 : !!
971 : !! INPUTS
972 : !! ikbz=Index in the BZ of the k-point whose contribution to chi0 has to be added,
973 : !! if we use symmetries, the contribution to chi0 by this k-point has to be symmetrized.
974 : !! isym_kbz=Index of the symmetry such as k_bz = IS k_ibz
975 : !! itim_kbz=2 if time-reversal has to be used to obtain k_bz, 1 otherwise.
976 : !! my_wl,my_wr=min and Max frequency index treated by this processor.
977 : !! npwe=Number of plane waves used to describe chi0.
978 : !! npwepG0=Maximum number of G vectors to account for umklapps.
979 : !! nomega=Number of frequencies in the imaginary part.
980 : !! rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states.
981 : !! rhotwx(3,nspinor**2)=Matrix elements of the gradient and of the commutator of the non-local operator with
982 : !! the position operator. The second term is present only if inclvkb=1,2.
983 : !! Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors
984 : !! %ng=number of G vectors in the enlarged sphere. It MUST be equal to the size of rhotwg
985 : !! %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion
986 : !! %phmGt(ng,nsym)=phase factors associated to non-symmorphic operations
987 : !! Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point.
988 : !! %timrev=2 it time-reversal is used, 1 otherwise
989 : !! %nsym_sg=Number of space group symmetries
990 : !! %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point
991 : !! %flag_umklp(timrev,nsym)= flag for umklapp processes
992 : !! if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise
993 : !! Cryst<crystal_t>=Info on unit cell and it symmetries
994 : !! %nsym=Number of symmetry operations.
995 : !! %symrec(3,3,nsym)=Symmetry operations in reciprocal space (reduced coordinates).
996 : !!
997 : !! OUTPUT
998 : !! (see side effects)
999 : !!
1000 : !! SIDE EFFECTS
1001 : !! sf_chi0(npwe,npwe,my_wl:my_wr)=Updated spectral function at q==0.
1002 : !! sf_lwing(npwe,my_wl:my_wr,3)=Updated lower wing of the spectral function.
1003 : !! sf_uwing(npwe,mw_wl:my_wr,3)=Updated upper wing of the spectral function.
1004 : !! sf_head(3,3,my_wl:my_wr)=Updated head of the spectral function.
1005 : !!
1006 : !! SOURCE
1007 :
1008 5664 : subroutine accumulate_sfchi0_q0(ikbz,isym_kbz,itim_kbz,nspinor,symchi,npwepG0,npwe,Cryst,Ltg_q,Gsph_epsG0,&
1009 5664 : & factocc,my_wl,iomegal,wl,my_wr,iomegar,wr,rhotwx,rhotwg,nomegasf,sf_chi0,sf_head,sf_lwing,sf_uwing)
1010 :
1011 : !Arguments ------------------------------------
1012 : !scalars
1013 : integer,intent(in) :: ikbz,my_wl,my_wr,nomegasf,npwe,npwepG0,nspinor
1014 : integer,intent(in) :: isym_kbz,itim_kbz,symchi,iomegal,iomegar
1015 : real(dp),intent(in) :: factocc,wl,wr
1016 : type(littlegroup_t),intent(in) :: Ltg_q
1017 : type(gsphere_t),target,intent(in) :: Gsph_epsG0
1018 : type(crystal_t),intent(in) :: Cryst
1019 : !arrays
1020 : complex(gwp),intent(in) :: rhotwg(npwepG0)
1021 : complex(gwp),intent(in) :: rhotwx(3,nspinor**2)
1022 : complex(gwp),intent(inout) :: sf_chi0(npwe,npwe,my_wl:my_wr)
1023 : complex(dp),intent(inout) :: sf_head(3,3,my_wl:my_wr)
1024 : complex(dp),intent(inout) :: sf_lwing(npwe,my_wl:my_wr,3)
1025 : complex(dp),intent(inout) :: sf_uwing(npwe,my_wl:my_wr,3)
1026 :
1027 : !Local variables-------------------------------
1028 : !scalars
1029 : integer :: itim,isym,idir,jdir
1030 : complex(gwp) :: num
1031 : character(len=500) :: msg
1032 : !arrays
1033 5664 : integer, contiguous, pointer :: Sm1G(:)
1034 : complex(dp) :: mir_kbz(3)
1035 5664 : complex(gwp), contiguous, pointer :: phmGt(:)
1036 5664 : complex(gwp),allocatable :: rhotwg_sym(:)
1037 : !************************************************************************
1038 :
1039 5664 : if (iomegal<my_wl .or. iomegar>my_wr) then
1040 0 : write(msg,'(3a,2(a,i0,a,i0,a))')ch10,&
1041 0 : 'Indices out of boundary ',ch10,&
1042 0 : ' my_wl = ',my_wl,' iomegal = ',iomegal,ch10,&
1043 0 : ' my_wr = ',my_wr,' iomegar = ',iomegar,ch10
1044 0 : ABI_BUG(msg)
1045 : end if
1046 :
1047 5664 : SELECT CASE (symchi)
1048 : CASE (0)
1049 : !
1050 : ! Calculation without symmetries
1051 : ! rhotwg(1)= R^-1q*rhotwx_ibz
1052 : ! rhotwg(1)=-R^-1q*conjg(rhotwx_ibz) for inversion
1053 0 : if (wl<huge(zero)*1.d-11) then
1054 : !this is awful but it is still a first coding
1055 : ! Num is single precision needed for cgerc check factocc
1056 0 : num=-wl*factocc
1057 0 : call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,sf_chi0(:,:,iomegal),npwe)
1058 : end if
1059 : ! Last point, must accumulate left point but not the right one
1060 0 : if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then
1061 0 : num=-wr*factocc
1062 0 : call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,sf_chi0(:,:,iomegar),npwe)
1063 : end if
1064 :
1065 : ! ================================
1066 : ! ==== Update heads and wings ====
1067 : ! ================================
1068 :
1069 : ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
1070 0 : if (nspinor == 1) then
1071 0 : mir_kbz =(3-2*itim_kbz) * MATMUL(Cryst%symrec(:,:,isym_kbz),rhotwx(:,1))
1072 : else
1073 0 : mir_kbz = (3-2*itim_kbz) * MATMUL(Cryst%symrec(:,:,isym_kbz), sum(rhotwx(:,1:2), dim=2))
1074 : end if
1075 0 : if (itim_kbz==2) mir_kbz=conjg(mir_kbz)
1076 :
1077 0 : do jdir=1,3
1078 0 : if (wl<huge(zero)*1.d-11) then
1079 : ! this is awful but it is still a first coding
1080 : ! Num is single precision needed for cgerc check factocc
1081 0 : num=-wl*factocc
1082 0 : sf_uwing(:,iomegal,jdir) = sf_uwing(:,iomegal,jdir) + num * mir_kbz(jdir) * conjg(rhotwg(1:npwepG0))
1083 0 : sf_lwing(:,iomegal,jdir) = sf_lwing(:,iomegal,jdir) + num * rhotwg(1:npwepG0) * conjg(mir_kbz(jdir))
1084 0 : do idir=1,3
1085 0 : sf_head(idir,jdir,iomegal) = sf_head(idir,jdir,iomegal) + num * mir_kbz(idir) * conjg(mir_kbz(jdir))
1086 : end do
1087 : end if
1088 :
1089 : ! Last point, must accumulate left point but not the right one
1090 0 : if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then
1091 0 : num=-wr*factocc
1092 0 : sf_uwing(:,iomegar,jdir) = sf_uwing(:,iomegar,jdir) + num * mir_kbz(jdir) * conjg(rhotwg(1:npwepG0))
1093 0 : sf_lwing(:,iomegar,jdir) = sf_lwing(:,iomegar,jdir) + num * rhotwg(1:npwepG0) * conjg(mir_kbz(jdir))
1094 0 : do idir=1,3
1095 0 : sf_head(idir,jdir,iomegar) = sf_head(idir,jdir,iomegar) + num * mir_kbz(idir) * conjg(mir_kbz(jdir))
1096 : end do
1097 : end if
1098 : end do ! jdir
1099 :
1100 : CASE (1)
1101 : ! === Notes on the symmetrization of oscillator matrix elements ===
1102 : ! If Sq = q then M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1G} (k,q)
1103 : ! If -Sq = q then M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1G}^*(k,q)
1104 : !
1105 : ! In case of an umklapp process
1106 : ! If Sq = q+G_o then M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1(G-G_o} (k,q)
1107 : ! If -Sq = q+G_o then M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1(G-G-o)}^*(k,q)
1108 : !
1109 : ! rhotwg(1)= R^-1q*rhotwx_ibz
1110 : ! rhotwg(1)=-R^-1q*conjg(rhotwx_ibz) for inversion
1111 :
1112 16992 : ABI_MALLOC(rhotwg_sym,(npwe))
1113 :
1114 : ! Loop over symmetries of the space group and time-reversal
1115 277536 : do isym=1,Ltg_q%nsym_sg
1116 821280 : do itim=1,Ltg_q%timrev
1117 :
1118 815616 : if (Ltg_q%wtksym(itim,isym,ikbz)==1) then
1119 : ! This operation belongs to the little group and has to be considered to reconstruct the BZ
1120 52992 : phmGt => Gsph_epsG0%phmGt(1:npwe,isym) ! In these 2 lines mind the slicing (1:npwe)
1121 52992 : Sm1G => Gsph_epsG0%rottbm1(1:npwe,itim,isym)
1122 :
1123 52992 : SELECT CASE (itim)
1124 : CASE (1)
1125 1947264 : rhotwg_sym(1:npwe)=rhotwg(Sm1G(1:npwe))*phmGt(1:npwe)
1126 : CASE (2)
1127 0 : rhotwg_sym(1:npwe)=conjg(rhotwg(Sm1G(1:npwe)))*phmGt(1:npwe)
1128 : CASE DEFAULT
1129 52992 : ABI_BUG(sjoin('Wrong value of itim:', itoa(itim)))
1130 : END SELECT
1131 :
1132 : ! Multiply elements G,Gp of rhotwg_sym*num and accumulate in sf_chi0(G,Gp,io)
1133 52992 : if (wl<huge(zero)*1.d-11) then
1134 47907 : num=-wl*factocc
1135 47907 : call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,sf_chi0(:,:,iomegal),npwe)
1136 : end if
1137 :
1138 : ! Last point, must accumulate left point but not the right one
1139 52992 : if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then
1140 48573 : num=-wr*factocc
1141 48573 : call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,sf_chi0(:,:,iomegar),npwe)
1142 : end if
1143 :
1144 : ! Accumulate heads and wings.
1145 : ! Symmetrize <r> in full BZ: <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'}
1146 52992 : if (nspinor == 1) then
1147 1440768 : mir_kbz =(3-2*itim) * MATMUL(Cryst%symrec(:,:,isym),rhotwx(:,1))
1148 : else
1149 56832 : mir_kbz = (3-2*itim) * MATMUL(Cryst%symrec(:,:,isym), sum(rhotwx(:,1:2), dim=2))
1150 : end if
1151 52992 : if (itim==2) mir_kbz=conjg(mir_kbz)
1152 :
1153 211968 : do jdir=1,3
1154 158976 : if (wl<huge(zero)*1.d-11) then
1155 : ! this is awful but it is still a first coding
1156 : ! Num is single precision needed for cgerc check factocc
1157 143721 : num=-wl*factocc
1158 5597712 : sf_uwing(:,iomegal,jdir) = sf_uwing(:,iomegal,jdir) + num * mir_kbz(jdir) * conjg(rhotwg_sym(1:npwe))
1159 5597712 : sf_lwing(:,iomegal,jdir) = sf_lwing(:,iomegal,jdir) + num * rhotwg_sym(1:npwe) * conjg(mir_kbz(jdir))
1160 574884 : do idir=1,3
1161 574884 : sf_head(idir,jdir,iomegal) = sf_head(idir,jdir,iomegal) + num * mir_kbz(idir) * conjg(mir_kbz(jdir))
1162 : end do
1163 : end if
1164 :
1165 : ! Last point, must accumulate left point but not the right one
1166 211968 : if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then
1167 145719 : num=-wr*factocc
1168 5629680 : sf_uwing(:,iomegar,jdir) = sf_uwing(:,iomegar,jdir) + num * mir_kbz(jdir) * conjg(rhotwg_sym(1:npwe))
1169 5629680 : sf_lwing(:,iomegar,jdir) = sf_lwing(:,iomegar,jdir) + num * rhotwg_sym(1:npwe) * conjg(mir_kbz(jdir))
1170 582876 : do idir=1,3
1171 582876 : sf_head(idir,jdir,iomegar) = sf_head(idir,jdir,iomegar) + num * mir_kbz(idir) * conjg(mir_kbz(jdir))
1172 : end do
1173 : end if
1174 : end do ! jdir
1175 :
1176 : end if !wtksym
1177 : end do !inv
1178 : end do !isym
1179 5664 : ABI_FREE(rhotwg_sym)
1180 :
1181 : CASE DEFAULT
1182 5664 : ABI_BUG(sjoin('Wrong value of symchi:', itoa(symchi)))
1183 : END SELECT
1184 :
1185 5664 : end subroutine accumulate_sfchi0_q0
1186 : !!***
1187 :
1188 : !----------------------------------------------------------------------
1189 :
1190 : !!****f* m_chi0tk/assemblychi0sf
1191 : !! NAME
1192 : !! assemblychi0sf
1193 : !!
1194 : !! FUNCTION
1195 : !! Update the spectral function of the irreducible polarizability for the contribution
1196 : !! of one pair of occupied-unoccupied states, for each frequency.
1197 : !! If symchi==1, the symmetries of the little group of the external q-point are used
1198 : !! to symmetrize the contribution in the full Brillouin zone. In this case, the routine computes:
1199 : !!
1200 : !! $ chi0(G1,G2,io)=chi0(G1,G2,io)+\sum_S (rhotwg(G1)*rhotwg^\dagger(G2))*\delta(w - trans) $
1201 : !!
1202 : !! where S are the symmetries of the little group of the external q-point.
1203 : !!
1204 : !! INPUTS
1205 : !! ik_bz=Index of the k-point in the BZ whose contribution has to be added to the spectral function of chi0
1206 : !! If symchi=1, the contribution is symmetrized.
1207 : !! my_wl,my_wr=min and Max frequency index treated by this processor.
1208 : !! npwe=Number of plane waves used to describe chi0.
1209 : !! npwepG0=Maximum number of G vectors taking into account umklapp vectors.
1210 : !! nomegasf=Number of frequencies for the spectral function.
1211 : !! nspinor=Number of spinorial components.
1212 : !! symchi=1 if symmetries are used, 0 otherwise
1213 : !! rhotwg(npwepG0)=Oscillator matrix elements corresponding to an occupied-unoccupied pair of states.
1214 : !! timrev=if 2, time reversal has to be used to obtain k_bz; 1 otherwise.
1215 : !! Gsph_epsG0<gsphere_t> Information on the "enlarged" G-sphere used for chi0, it contains umklapp G0 vectors
1216 : !! %ng=number of G vectors in the enlarged sphere, actually MUST be equal to the size of rhotwg
1217 : !! %rottbm1(ng,2,nsym)=index of (IR)^{-1} G where I is the identity or the inversion
1218 : !! %phmGt(ng,nsym)=phase factors associated to non-simmorphic operations
1219 : !! Ltg_q<littlegroup_t_type>=Info on the little group associated to the external q-point.
1220 : !! %timrev=2 it time-reversal is used, 1 otherwise
1221 : !! %nsym_sg=Number of space group symmetries
1222 : !! %wtksym(2,nsym,nkbz)=1 if the symmetry (with or without time-reversal) must be considered for this k-point
1223 : !! %flag_umklp(timrev,nsym)= flag for umklapp processes
1224 : !! if 1 that the particular operation (IS) requires a G_o to preserve Q, 0 otherwise
1225 : !! %igmG0(npwepG0,timrev,nsym) index of G-G0 in the array gvec
1226 : !! factocc=occupation factor=f_occ*(ockp-occk) (see cchi0.F90)
1227 : !! wl,wr=Weights used to approximate the delta function.
1228 : !!
1229 : !! OUTPUT
1230 : !! (see side effects)
1231 : !!
1232 : !! SIDE EFFECTS
1233 : !! chi0sf(npwe,npwe,my_wl:my_wr)= updated spectral function.
1234 : !!
1235 : !! NOTES
1236 : !! Umklapp processes are not yet implemented
1237 : !!
1238 : !! SOURCE
1239 :
1240 252320 : subroutine assemblychi0sf(ik_bz,symchi,Ltg_q,npwepG0,npwe,rhotwg,Gsph_epsG0,&
1241 252320 : & factocc,my_wl,iomegal,wl,my_wr,iomegar,wr,nomegasf,chi0sf)
1242 :
1243 : !Arguments ------------------------------------
1244 : !scalars
1245 : integer,intent(in) :: ik_bz,iomegal,iomegar,my_wl,my_wr,nomegasf,npwe,npwepG0
1246 : integer,intent(in) :: symchi
1247 : real(dp),intent(in) :: factocc,wl,wr
1248 : type(gsphere_t),intent(in) :: Gsph_epsG0
1249 : type(littlegroup_t),intent(in) :: Ltg_q
1250 : !arrays
1251 : complex(gwp),intent(in) :: rhotwg(npwepG0)
1252 : complex(gwp),intent(inout) :: chi0sf(npwe,npwe,my_wl:my_wr)
1253 :
1254 : !Local variables-------------------------------
1255 : !scalars
1256 : integer :: isym,itim,ig1,ig2
1257 : complex(gwp) :: num
1258 : character(len=500) :: msg
1259 : !arrays
1260 504640 : integer :: Sm1_gmG0(npwe)
1261 504640 : complex(gwp) :: rhotwg_sym(npwe)
1262 : ! *************************************************************************
1263 :
1264 252320 : if (iomegal < my_wl .or. iomegar > my_wr) then
1265 0 : write(msg,'(3a,2(a,i0,a,i0,a))')ch10,&
1266 0 : & ' Indices out of boundary ',ch10,&
1267 0 : & ' my_wl = ',my_wl,' iomegal = ',iomegal,ch10,&
1268 0 : & ' my_wr = ',my_wr,' iomegar = ',iomegar,ch10
1269 0 : ABI_BUG(msg)
1270 : end if
1271 :
1272 252320 : SELECT CASE (symchi)
1273 : CASE (0)
1274 : ! Do not use symmetries.
1275 :
1276 : ! MG: This is the best I can do for this part.
1277 : !$omp PARALLEL private(num)
1278 : !$omp SECTIONS
1279 : !$omp SECTION
1280 0 : if (wl<huge(zero)*1.d-11) then !FIXME this is awful
1281 0 : num=-wl*factocc
1282 0 : call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,chi0sf(:,:,iomegal),npwe)
1283 : end if
1284 :
1285 : ! Last point, must accumulate left point but not the right one
1286 : !$omp SECTION
1287 0 : if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then
1288 0 : num=-wr*factocc
1289 0 : call XGERC(npwe,npwe,num,rhotwg,1,rhotwg,1,chi0sf(:,:,iomegar),npwe)
1290 : end if
1291 : !$omp end SECTIONS
1292 : !$omp end PARALLEL
1293 :
1294 : CASE (1)
1295 : ! Use symmetries to reconstruct oscillator matrix elements
1296 : ! Notes on the symmetrization of the oscillator maxtri elements:
1297 : !
1298 : ! If Sq=q then M_G^( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1G} (k,q)
1299 : ! If -Sq=q then M_G^(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1G}^*(k,q)
1300 : !
1301 : ! In case of an umklapp process
1302 : ! If Sq=q+G_o then M_G( Sk,q)= e^{-i(q+G)\cdot t} M_{ S^-1(G-G_o} (k,q)
1303 : ! If -Sq=q+G_o then M_G(-Sk,q)= e^{-i(q+G)\cdot t} M_{-S^-1(G-G_o)}^*(k,q)
1304 : !
1305 : ! Ltg_q%igmG0(ig,itim,isym) contains the index of G-G0 where ISq=q+G0
1306 : ! Note that there is no need to take into account the phases due to q,
1307 : ! They cancel in the scalar product ==> phmGt(G,isym)=e^{-iG\cdot t}
1308 : !
1309 : ! Mind the slicing of %rottbm1(npwepG0,timrev,nsym) and %phgt(npwepG0,nsym) as
1310 : ! these arrays, usually, do not conform to rho_twg_sym(npw) !
1311 : !
1312 : !ABI_MALLOC(rhotwg_sym,(npwe))
1313 : !
1314 : ! Loop over symmetries of the space group and time-reversal
1315 12363680 : do isym=1,Ltg_q%nsym_sg
1316 36586400 : do itim=1,Ltg_q%timrev
1317 :
1318 36334080 : if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
1319 : ! This operation belongs to the little group and has to be used to reconstruct BZ.
1320 : ! TODO this is a hot-spot, should add a test on the umklapp
1321 : !
1322 : ! In these 3 lines mind the slicing (1:npwe)
1323 46058752 : Sm1_gmG0(1:npwe)=Gsph_epsG0%rottbm1( Ltg_q%igmG0(1:npwe,itim,isym), itim,isym)
1324 :
1325 583168 : SELECT CASE (itim)
1326 : CASE (1)
1327 23320960 : rhotwg_sym(1:npwe)=rhotwg(Sm1_gmG0(1:npwe)) * Gsph_epsG0%phmGt(1:npwe,isym)
1328 : CASE (2)
1329 0 : rhotwg_sym(1:npwe)=conjg(rhotwg(Sm1_gmG0(1:npwe))) * Gsph_epsG0%phmGt(1:npwe,isym)
1330 : CASE DEFAULT
1331 583168 : ABI_BUG(sjoin('Wrong value for itim:', itoa(itim)))
1332 : END SELECT
1333 :
1334 : #if 0
1335 : !! MG: This is the best I can do, at present.
1336 : !$omp PARALLEL private(num)
1337 : !$omp SECTIONS
1338 :
1339 : !$omp SECTION
1340 : if (wl<huge(zero)*1.d-11) then
1341 : num=-wl*factocc
1342 : call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe)
1343 : end if
1344 : !$omp SECTION
1345 : !
1346 : ! Last point, must accumulate left point but not the right one
1347 : if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then
1348 : num=-wr*factocc
1349 : call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegar),npwe)
1350 : end if
1351 : !$omp end SECTIONS
1352 : !$omp end PARALLEL
1353 : #else
1354 :
1355 583168 : if (wl<huge(zero)*1.d-11) then
1356 : !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe)
1357 512019 : num=-wl*factocc
1358 : !$omp parallel do
1359 22182576 : do ig2=1,npwe
1360 1108179555 : do ig1=1,npwe
1361 1107667536 : chi0sf(ig1,ig2,iomegal) = chi0sf(ig1,ig2,iomegal) + num * rhotwg_sym(ig1) * conjg(rhotwg_sym(ig2))
1362 : end do
1363 : end do
1364 : end if
1365 :
1366 : ! Last point, must accumulate left point but not the right one
1367 583168 : if (iomegar/=nomegasf+1 .and. wr<huge(zero)*1.d-11) then
1368 : !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegar),npwe)
1369 511757 : num=-wr*factocc
1370 : !$omp parallel do
1371 22178384 : do ig2=1,npwe
1372 1108116413 : do ig1=1,npwe
1373 : !call XGERC(npwe,npwe,num,rhotwg_sym,1,rhotwg_sym,1,chi0sf(:,:,iomegal),npwe)
1374 1107604656 : chi0sf(ig1,ig2,iomegar) = chi0sf(ig1,ig2,iomegar) + num * rhotwg_sym(ig1) * conjg(rhotwg_sym(ig2))
1375 : end do
1376 : end do
1377 : end if
1378 : #endif
1379 : end if !wtksym
1380 :
1381 : end do !inv
1382 : end do !isym
1383 : !ABI_FREE(rhotwg_sym)
1384 :
1385 : CASE DEFAULT
1386 252320 : ABI_BUG(sjoin('Wrong value for symchi:', itoa(symchi)))
1387 : END SELECT
1388 :
1389 252320 : end subroutine assemblychi0sf
1390 : !!***
1391 :
1392 : !----------------------------------------------------------------------
1393 :
1394 : !!****f* m_chi0tk/approxdelta
1395 : !! NAME
1396 : !! approxdelta
1397 : !!
1398 : !! FUNCTION
1399 : !! Approximate the Dirac function using two methods:
1400 : !! method 1) a triangular funtion centered at the value egwdiff_re, Eq 17 of PRB 74, 035101 (2006) [[cite:Shishkin2006]]
1401 : !! method 2) a gaussian of witdth ep%spsmear expandended in Taylor series
1402 : !! (at the moment only the 0-th moments)
1403 : !!
1404 : !! Subroutine needed to implement the calculation
1405 : !! of the polarizability using the spectral representation as proposed in:
1406 : !! PRB 74, 035101 (2006) [[cite:Shishkin2006]]
1407 : !! and PRB 61, 7172 (2000) [[cite:Miyake2000]]
1408 : !!
1409 : !! INPUTS
1410 : !! nomegasf=number of frequencies in the grid for Im \chi_0
1411 : !! omegasf(0:nomega+1)= frequencies (real)
1412 : !! egwdiff_re = transition energy where the delta function is centered
1413 : !!
1414 : !! method= 1: a triangular shaped function used to approximated the delta
1415 : !! 2: gaussian approximation with standard deviation (smear)
1416 : !! smear= used only in case of method==2, defines the width of the gaussian
1417 : !!
1418 : !! OUTPUT
1419 : !! wl = weight associated to omegal (last omega wich is smaller than egwdiff_re
1420 : !! wr = weight associate to omegar (first omega larger than egwdff_re
1421 : !! iomegal= index in the array omegasf of the last frequency < egwdiff
1422 : !! iomegar= index in the array omegasf of the first frequency > egwdiff
1423 : !!
1424 : !! SOURCE
1425 :
1426 257984 : subroutine approxdelta(nomegasf,omegasf,egwdiff_re,smear,iomegal,iomegar,wl,wr,spmeth)
1427 :
1428 : !Arguments ------------------------------------
1429 : !scalars
1430 : integer,intent(in) :: nomegasf,spmeth
1431 : integer,intent(out) :: iomegal,iomegar
1432 : real(dp),intent(in) :: egwdiff_re,smear
1433 : real(dp),intent(out) :: wl,wr
1434 : !arrays
1435 : real(dp),intent(in) :: omegasf(nomegasf)
1436 :
1437 : !Local variables-------------------------------
1438 : integer :: io,iomega
1439 : real(dp) :: omegal,omegar,deltal,deltar
1440 : !character(len=500) :: msg
1441 : ! *************************************************************************
1442 :
1443 257984 : iomega=-999
1444 26714312 : do io=nomegasf,1,-1
1445 26714312 : if (omegasf(io)<egwdiff_re) then
1446 : iomega=io; EXIT
1447 : end if
1448 : end do
1449 :
1450 257984 : iomegal=iomega ; omegal=omegasf(iomegal)
1451 257984 : iomegar=iomegal+1; omegar=omegasf(iomegar)
1452 :
1453 452168 : SELECT CASE (spmeth)
1454 : CASE (1)
1455 : ! Weights for triangular shaped function
1456 194184 : wr= (egwdiff_re-omegal)/(omegar-omegal)
1457 194184 : wl= -(egwdiff_re-omegar)/(omegar-omegal)
1458 :
1459 : CASE (2)
1460 : ! Weights for gaussian method (0-th moment)
1461 63800 : deltal=(egwdiff_re-omegal)/smear
1462 63800 : deltar=(omegar-egwdiff_re)/smear
1463 63800 : if (deltar>=deltal) then
1464 31949 : wl=EXP(-deltal*deltal)
1465 : ! this value is used to avoid double counting and speed-up
1466 31949 : wr=huge(one)*1.d-10
1467 : else
1468 31851 : wl=huge(one)*1.d-10
1469 31851 : wr=exp(-deltal*deltal)
1470 : end if
1471 :
1472 : CASE DEFAULT
1473 257984 : ABI_BUG(sjoin('Wrong value for spmeth:', itoa(spmeth)))
1474 : END SELECT
1475 :
1476 257984 : end subroutine approxdelta
1477 : !!***
1478 :
1479 : !----------------------------------------------------------------------
1480 :
1481 : !!****f* m_chi0tk/calc_kkweight
1482 : !! NAME
1483 : !! calc_kkweight
1484 : !!
1485 : !! FUNCTION
1486 : !! Calculate frequency dependent weights needed to perform the Hilbert transform
1487 : !!
1488 : !! Subroutine needed to implement the calculation
1489 : !! of the polarizability using the spectral representation as proposed in:
1490 : !! PRB 74, 035101 (2006) [[cite:Shishkin2006]]
1491 : !! and PRB 61, 7172 (2000) [[cite:Miyake2000]]
1492 : !!
1493 : !! INPUTS
1494 : !! nsp=number of frequencies where the imaginary part of the polarizability is evaluated
1495 : !! ne=number of frequencies for the polarizability (same as in epsilon^-1)
1496 : !! omegasp(nsp)=real frequencies for the imaginary part of the polarizability
1497 : !! omegae(ne)= imaginary frequencies for the polarizability
1498 : !! delta=small imaginary part used to avoid poles, input variables
1499 : !!
1500 : !! OUTPUT
1501 : !! kkweight(nsp,ne)=frequency dependent weights Eq A1 PRB 74, 035101 (2006) [[cite:Shishkin2006]]
1502 : !!
1503 : !! SOURCE
1504 : !!
1505 :
1506 54 : subroutine calc_kkweight(ne,omegae,nsp,omegasp,delta,omegamax,kkw)
1507 :
1508 : !Arguments ------------------------------------
1509 : !scalars
1510 : integer,intent(in) :: ne,nsp
1511 : real(dp),intent(in) :: delta,omegamax
1512 : !arrays
1513 : real(dp),intent(in) :: omegasp(nsp)
1514 : complex(dp),intent(in) :: omegae(ne)
1515 : complex(dp),intent(out) :: kkw(nsp,ne)
1516 :
1517 : !Local variables-------------------------------
1518 : !scalars
1519 : integer :: isp,je
1520 : real(dp) :: eta,xx1,xx2,den1,den2
1521 : complex(dp) :: c1,c2,wt
1522 : !************************************************************************
1523 :
1524 : DBG_ENTER("COLL")
1525 :
1526 190187 : kkw(:,:)=czero
1527 :
1528 837 : do je=1,ne
1529 783 : eta=delta
1530 783 : wt=omegae(je)
1531 : ! Not include shift at omega==0, what about metallic systems?
1532 783 : if (abs(real(omegae(je)))<tol6 .and. abs(aimag(wt))<tol6) eta=tol12
1533 : ! Not include shift along the imaginary axis
1534 783 : if (abs(aimag(wt))>tol6) eta=zero
1535 190187 : do isp=1,nsp
1536 189350 : if (isp==1) then
1537 : ! Skip negative point, should check that this would not lead to spurious effects
1538 : c1=czero
1539 : den1=one
1540 : else
1541 188567 : xx1=omegasp(isp-1)
1542 188567 : xx2=omegasp(isp)
1543 188567 : den1= xx2-xx1
1544 : c1= -(wt-xx1+j_dpc*eta)*log( (wt-xx2+j_dpc*eta)/(wt-xx1+j_dpc*eta) )&
1545 188567 : & +(wt+xx1-j_dpc*eta)*log( (wt+xx2-j_dpc*eta)/(wt+xx1-j_dpc*eta) )
1546 188567 : c1= c1/den1
1547 : end if
1548 189350 : xx1=omegasp(isp)
1549 189350 : if (isp==nsp) then
1550 : ! Skip last point should check that this would not lead to spurious effects
1551 783 : xx2=omegamax
1552 : else
1553 188567 : xx2=omegasp(isp+1)
1554 : end if
1555 189350 : den2=xx2-xx1
1556 : c2= (wt-xx2+j_dpc*eta)*log( (wt-xx2+j_dpc*eta)/(wt-xx1+j_dpc*eta) )&
1557 189350 : & -(wt+xx2-j_dpc*eta)*log( (wt+xx2-j_dpc*eta)/(wt+xx1-j_dpc*eta) )
1558 189350 : c2= c2/den2
1559 190133 : kkw(isp,je)= c1/den1 + c2/den2
1560 : end do
1561 : end do
1562 :
1563 : DBG_EXIT("COLL")
1564 :
1565 54 : end subroutine calc_kkweight
1566 : !!***
1567 :
1568 : !----------------------------------------------------------------------
1569 :
1570 : !!****f* m_chi0tk/setup_spectral
1571 : !! NAME
1572 : !! setup_spectral
1573 : !!
1574 : !! FUNCTION
1575 : !! Calculation of \chi_o based on the spectral method as proposed in PRB 74, 035101 (2006) [[cite:Shishkin2006]]
1576 : !! and PRB 61, 7172 (2000) [[cite:Miyake2000]].
1577 : !! Setup of the real frequency mesh for $\Im\chi_o$ and of the frequency-dependent weights for
1578 : !! Hilbert transform. Note that CPU time does not depend dramatically on nomegasf unlike memory.
1579 : !! spmeth defines the approximant for the delta function:
1580 : !! ==1 : use Triangular approximant (Kresse method)
1581 : !! ==2 : use Gaussian method, requiring smearing (Miyake method)
1582 : !!
1583 : !! INPUTS
1584 : !! nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$
1585 : !! nomega=number of frequencies in $\chi0(q,\omega)$.
1586 : !! max_rest,min_res=max and min resonant transition energy (for this q-point)
1587 : !! my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor
1588 : !! method=integer flag defining the type of frequency mesh used for $\Im chi0$
1589 : !! | 0 for a linear mesh
1590 : !! | 1 for a mesh densified around omegaplasma
1591 : !! omegaplasma=frequency around which the mesh is densifies (usually Drude plasma frequency)
1592 : !! used only in case of method==1
1593 : !! zcut=small imaginary shift to avoid pole in chi0
1594 : !!
1595 : !! OUTPUT
1596 : !! kkweight(nomegasf,nomega)=Frequency dependent weight for Hilber transform.
1597 : !! omegasf(nomegasf+1)=frequencies for imaginary part.
1598 : !!
1599 : !! SOURCE
1600 :
1601 54 : subroutine setup_spectral(nomega,omega,nomegasf,omegasf,max_rest,min_rest,my_max_rest,my_min_rest,&
1602 54 : & method,zcut,omegaplasma,my_wl,my_wr,kkweight)
1603 :
1604 : !Arguments ------------------------------------
1605 : !scalars
1606 : integer,intent(in) :: method,nomega,nomegasf
1607 : integer,intent(out) :: my_wl,my_wr
1608 : real(dp),intent(in) :: max_rest,min_rest,omegaplasma,zcut
1609 : real(dp),intent(in) :: my_max_rest,my_min_rest
1610 : !arrays
1611 : real(dp),intent(out) :: omegasf(nomegasf)
1612 : complex(dp),intent(in) :: omega(nomega)
1613 : complex(dp),intent(out) :: kkweight(nomegasf,nomega)
1614 :
1615 : !Local variables-------------------------------
1616 : !scalars
1617 : integer :: io,ii
1618 : real(dp) :: nu_min,nu_max,nu1,nu2,dd,domegasf,wp,deltat
1619 : character(len=500) :: msg
1620 : !arrays
1621 54 : integer,allocatable :: insort(:)
1622 : !************************************************************************
1623 :
1624 : ! The mesh must enclose the entire range of transitions.
1625 54 : dd=(max_rest-min_rest)/(nomegasf-1)
1626 54 : domegasf=(max_rest-min_rest+2*dd)/(nomegasf-1)
1627 :
1628 54 : write(msg,'(4a,f8.3,3a,i5,2a,f8.5,a)')ch10,&
1629 54 : ' === Info on the real frequency mesh for spectral method === ',ch10,&
1630 54 : ' maximum frequency = ',max_rest*Ha_eV,' [eV]',ch10,&
1631 54 : ' nomegasf = ',nomegasf,ch10,&
1632 108 : ' domegasf = ',domegasf*Ha_eV,' [eV]'
1633 54 : call wrtout(std_out,msg)
1634 :
1635 54 : if (min_rest<tol6) then
1636 0 : ABI_WARNING("System seems to be metallic")
1637 : end if
1638 :
1639 : ! ======================================================
1640 : ! === Setup of the w-mesh for the spectral function ====
1641 : ! ======================================================
1642 108 : SELECT CASE (method)
1643 : CASE (0)
1644 : ! Linear mesh.
1645 54 : call wrtout(std_out, ' Using linear mesh for Im chi0')
1646 10354 : do io=1,nomegasf
1647 10354 : omegasf(io)=(io-1)*domegasf+min_rest-dd
1648 : end do
1649 :
1650 : CASE (1)
1651 : ! Non-homogeneous mesh densified around omega_plasma, do not improve results ===
1652 : ! WARNING_ this part has to be checked since I modified omegasf
1653 0 : write(msg,'(a,f7.4,a)')' Using mesh densified around ',omegaplasma*Ha_eV,' [eV] '
1654 0 : call wrtout(std_out, msg)
1655 0 : wp=omegaplasma ; deltat=max_rest-min_rest
1656 0 : nu_min=zero
1657 0 : if (deltat<wp ) then
1658 0 : nu_max = wp/sqrt2 * ATAN(sqrt2*deltat*wp/(-deltat**2+wp**2))
1659 : else
1660 0 : nu_max = wp/sqrt2 * ( ATAN(sqrt2*deltat*wp/(-deltat**2+wp**2)) + pi)
1661 : end if
1662 0 : domegasf=(nu_max-nu_min)/(nomegasf+1)
1663 : !write(std_out,*) -(wp/sqrt2) * atan(sqrt2*deltat*wp/(deltat**2-wp**2))
1664 0 : omegasf(1)=zero ; omegasf(nomegasf+1)=deltat
1665 0 : ii=0
1666 0 : do io=2,nomegasf
1667 0 : nu1=domegasf*(io-1) ; nu2=TAN(-sqrt2*nu1/wp)
1668 0 : if (nu2<0) then
1669 0 : omegasf(io) = wp * (one - SQRT(1+2*nu2**2))/(sqrt2*nu2)
1670 : else
1671 0 : omegasf(io) = wp * (one + SQRT(1+2*nu2**2))/(sqrt2*nu2)
1672 : end if
1673 0 : if (omegasf(io)> deltat ) then
1674 0 : omegasf(io)= deltat-0.1*ii
1675 0 : ii=ii+1
1676 : end if
1677 : ! write(102,'(i4,2x,3(f9.4,2x))')io,nu1,nu2,ep%omegasf(io)*Ha_eV
1678 : end do
1679 :
1680 : ! Reorder frequencies in ascending order
1681 0 : ABI_MALLOC(insort,(nomegasf+1))
1682 0 : insort(:)=(/ (io,io=1,nomegasf+1) /)
1683 0 : call sort_dp(nomegasf+1,omegasf,insort,tol14)
1684 0 : ABI_FREE(insort)
1685 :
1686 : CASE DEFAULT
1687 54 : ABI_BUG(sjoin('Wrong value for method:', itoa(method)))
1688 : END SELECT
1689 : !write(std_out,*)omegasf(1)*Ha_eV,omegasf(nomegasf)*Ha_eV
1690 :
1691 : ! Find min and max index in omegasf treated by this processor.
1692 54 : my_wr=-999
1693 10300 : do io=1,nomegasf
1694 10300 : if (omegasf(io)>my_max_rest) then
1695 54 : my_wr=io; EXIT
1696 : end if
1697 : end do
1698 54 : if (my_wr==nomegasf+2) my_wr=nomegasf+1
1699 54 : my_wl=-999
1700 10300 : do io=nomegasf,1,-1
1701 10300 : if (omegasf(io)< my_min_rest) then ! Check metals
1702 54 : my_wl=io; EXIT
1703 : end if
1704 : end do
1705 :
1706 54 : write(msg,'(a,2(1x,i0))')' my_wl and my_wr:',my_wl,my_wr
1707 54 : call wrtout(std_out, msg)
1708 :
1709 54 : if (my_wl==-999 .or. my_wr==-999) then
1710 0 : write(msg,'(a,2i6)')' wrong value in my_wl and/or my_wr ',my_wl,my_wr
1711 0 : ABI_ERROR(msg)
1712 : end if
1713 :
1714 : ! Calculate weights for Hilbert transform.
1715 54 : call calc_kkweight(nomega,omega,nomegasf,omegasf,zcut,max_rest,kkweight)
1716 :
1717 54 : end subroutine setup_spectral
1718 : !!***
1719 :
1720 : !----------------------------------------------------------------------
1721 :
1722 : !!****f* m_chi0tk/hilbert_transform
1723 : !! NAME
1724 : !! hilbert_transform
1725 : !!
1726 : !! FUNCTION
1727 : !! Compute the hilbert transform.
1728 : !!
1729 : !! INPUTS
1730 : !! nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$
1731 : !! nomega=number of frequencies in $\chi0(q,\omega)$.
1732 : !! max_rest,min_res=max and min resonant transition energy (for this q-point)
1733 : !! my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor
1734 : !!
1735 : !! OUTPUT
1736 : !!
1737 : !! SOURCE
1738 :
1739 54 : subroutine hilbert_transform(npwe,nomega,nomegasf,my_wl,my_wr,kkweight,sf_chi0,chi0,spmeth)
1740 :
1741 : !Arguments ------------------------------------
1742 : !scalars
1743 : integer,intent(in) :: spmeth,nomega,nomegasf,my_wl,my_wr,npwe
1744 : !arrays
1745 : complex(dp),intent(in) :: kkweight(nomegasf,nomega)
1746 : complex(gwp),intent(inout) :: sf_chi0(npwe,npwe,my_wl:my_wr)
1747 : complex(gwp),intent(inout) :: chi0(npwe,npwe,nomega)
1748 :
1749 : !Local variables-------------------------------
1750 : !scalars
1751 : integer :: ig2,my_nwp
1752 : character(len=500) :: msg
1753 : !arrays
1754 54 : complex(gwp),allocatable :: A_g1wp(:,:),H_int(:,:),my_kkweight(:,:)
1755 : !************************************************************************
1756 :
1757 : #ifdef HAVE_OPENMP
1758 : write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform (with OpenMP) using method ',spmeth,' It might take some time...'
1759 : #else
1760 54 : write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform using method ',spmeth,' It might take some time...'
1761 : #endif
1762 54 : call wrtout(std_out, msg, do_flush=.True.)
1763 :
1764 54 : my_nwp = my_wr - my_wl +1
1765 :
1766 : !$omp parallel private(my_kkweight, A_g1wp, H_int, ig2)
1767 216 : ABI_MALLOC(my_kkweight, (my_wl:my_wr,nomega))
1768 190241 : my_kkweight = kkweight(my_wl:my_wr,:)
1769 :
1770 216 : ABI_MALLOC(A_g1wp, (npwe, my_nwp))
1771 216 : ABI_MALLOC(H_int, (npwe, nomega))
1772 :
1773 : !$omp do
1774 1796 : do ig2=1,npwe
1775 18126484 : A_g1wp = sf_chi0(:,ig2,:)
1776 :
1777 : ! Compute H_int = MATMUL(A_g1wp,my_kkweight)
1778 1742 : call XGEMM('N','N',npwe,nomega,my_nwp,cone_gw,A_g1wp,npwe,my_kkweight,my_nwp,czero_gw,H_int,npwe)
1779 1429412 : chi0(:,ig2,:) = H_int
1780 : end do
1781 :
1782 54 : ABI_FREE(my_kkweight)
1783 54 : ABI_FREE(A_g1wp)
1784 54 : ABI_FREE(H_int)
1785 : !$omp end parallel
1786 :
1787 54 : end subroutine hilbert_transform
1788 : !!***
1789 :
1790 : !----------------------------------------------------------------------
1791 :
1792 : !!****f* m_chi0tk/hilbert_transform_headwings
1793 : !! NAME
1794 : !! hilbert_transform_headwings
1795 : !!
1796 : !! FUNCTION
1797 : !! Compute the hilbert transform the heads and wings of the polarizability.
1798 : !!
1799 : !! INPUTS
1800 : !! nomegasf=number of points for the imaginary part of $\chi0(q,\omega)$
1801 : !! nomega=number of frequencies in $\chi0(q,\omega)$.
1802 : !! max_rest,min_res=max and min resonant transition energy (for this q-point)
1803 : !! my_max_rest,my_min_rest=max and min resonant transition energy treated by this processor
1804 : !!
1805 : !! OUTPUT
1806 : !!
1807 : !! SOURCE
1808 :
1809 6 : subroutine hilbert_transform_headwings(npwe,nomega,nomegasf,my_wl,my_wr,kkweight, &
1810 6 : & sf_lwing,sf_uwing,sf_head,chi0_lwing,chi0_uwing,chi0_head,spmeth)
1811 :
1812 : !Arguments ------------------------------------
1813 : !scalars
1814 : integer,intent(in) :: spmeth,nomega,nomegasf,my_wl,my_wr,npwe
1815 : !arrays
1816 : complex(dp),intent(in) :: kkweight(nomegasf,nomega)
1817 : complex(dp),intent(inout) :: sf_lwing(npwe,my_wl:my_wr,3)
1818 : complex(dp),intent(inout) :: sf_uwing(npwe,my_wl:my_wr,3)
1819 : complex(dp),intent(inout) :: sf_head(3,3,my_wl:my_wr)
1820 : complex(dp),intent(inout) :: chi0_lwing(npwe,nomega,3)
1821 : complex(dp),intent(inout) :: chi0_uwing(npwe,nomega,3)
1822 : complex(dp),intent(inout) :: chi0_head(3,3,nomega)
1823 :
1824 : !Local variables-------------------------------
1825 : !scalars
1826 : integer :: ig1,idir,io,iw
1827 : complex(dp) :: kkw
1828 : character(len=500) :: msg
1829 : !************************************************************************
1830 :
1831 : #ifdef HAVE_OPENMP
1832 : write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform (with OpenMP) using method ',spmeth,' It might take some time...'
1833 : #else
1834 6 : write(msg,'(2a,i3,a)')ch10,' Performing Hilbert transform using method ',spmeth,' It might take some time...'
1835 : #endif
1836 6 : call wrtout(std_out,msg, do_flush=.True.)
1837 :
1838 : ! Hilbert transform of the head.
1839 110 : do io=1,nomega
1840 25704 : chi0_head(1,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,1,my_wl:my_wr))
1841 25704 : chi0_head(2,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,1,my_wl:my_wr))
1842 25704 : chi0_head(3,1,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,1,my_wl:my_wr))
1843 25704 : chi0_head(1,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,2,my_wl:my_wr))
1844 25704 : chi0_head(2,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,2,my_wl:my_wr))
1845 25704 : chi0_head(3,2,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,2,my_wl:my_wr))
1846 25704 : chi0_head(1,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(1,3,my_wl:my_wr))
1847 25704 : chi0_head(2,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(2,3,my_wl:my_wr))
1848 25710 : chi0_head(3,3,io) = SUM(kkweight(my_wl:my_wr,io)*sf_head(3,3,my_wl:my_wr))
1849 : end do
1850 :
1851 : ! Hilbert transform for wings.
1852 : ! Partial contributions to chi0 will be summed afterwards.
1853 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(kkw)
1854 24 : do idir=1,3
1855 336 : do io=1,nomega
1856 77130 : do iw=my_wl,my_wr
1857 76800 : kkw = kkweight(iw,io)
1858 2402112 : do ig1=1,npwe
1859 2325000 : chi0_lwing(ig1,io,idir) = chi0_lwing(ig1,io,idir) + kkw*sf_lwing(ig1,iw,idir)
1860 2401800 : chi0_uwing(ig1,io,idir) = chi0_uwing(ig1,io,idir) + kkw*sf_uwing(ig1,iw,idir)
1861 : end do
1862 : end do
1863 : end do
1864 : end do !idir
1865 :
1866 6 : end subroutine hilbert_transform_headwings
1867 : !!***
1868 :
1869 : !----------------------------------------------------------------------
1870 :
1871 : !!****f* m_chi0tk/completechi0_deltapart
1872 : !! NAME
1873 : !! completechi0_deltapart
1874 : !!
1875 : !! FUNCTION
1876 : !! Apply the delta part of the completeness correction to chi0
1877 : !!
1878 : !! INPUTS
1879 : !! ik_bz=Index of the k-point in the full BZ whose contribution has to be added and symmetrized.
1880 : !! qzero=.TRUE. is long wave-length limit.
1881 : !! symchi=1 if we are summing over IBZ_q and symmetrization has to be performed.
1882 : !! npwe=Number of G vectors in chi0.
1883 : !! npwvec=MAX number of G.
1884 : !! nomega=Number of frequencies.
1885 : !! nspinor=Number of spinorial components.
1886 : !! nfftot=Total Number of points in the FFT
1887 : !! ngfft(18)=Info on the FFT.
1888 : !! igfft0(npwvec)=Index of each G in the FFT array.
1889 : !! Gsph_FFT=<gsphere_t>=Info on the largest G-sphere contained in the FFT box used for wavefunctions.
1890 : !! Ltg_q=<littlegroup_t>= Structure gathering information on the little group of the external q.
1891 : !! green_enhigh_w=Approximated frequency dependent part of the Green function entering equation (TODO put reference)
1892 : !! wfwfg=Fourier components of u_{kb1}.u_{kb2}
1893 : !!
1894 : !! OUTPUT
1895 : !! See SIDES EFFECTS
1896 : !!
1897 : !! SIDES EFFECTS
1898 : !! chi0(npwe,npwe,nomega)= In input chi0 calculated so far,
1899 : !! In output the "delta part" of the completeness correction is added.
1900 : !!
1901 : !! SOURCE
1902 :
1903 1784 : subroutine completechi0_deltapart(ik_bz,qzero,symchi,npwe,npwvec,nomega,nspinor,&
1904 1784 : & nfftot,ngfft,igfft0,Gsph_FFT,Ltg_q,green_enhigh_w,wfwfg,chi0)
1905 :
1906 : !Arguments ------------------------------------
1907 : !scalars
1908 : integer,intent(in) :: ik_bz,nfftot,nomega,npwe,npwvec,nspinor,symchi
1909 : logical,intent(in) :: qzero
1910 : type(gsphere_t),intent(in) :: Gsph_FFT
1911 : type(littlegroup_t),intent(in) :: Ltg_q
1912 : !arrays
1913 : integer,intent(in) :: igfft0(npwvec),ngfft(18)
1914 : complex(dp),intent(in) :: green_enhigh_w(nomega)
1915 : complex(gwp),intent(in) :: wfwfg(nfftot*nspinor**2)
1916 : complex(gwp),intent(inout) :: chi0(npwe,npwe,nomega)
1917 :
1918 : !Local variables ------------------------------
1919 : !scalars
1920 : integer,save :: enough=0
1921 : integer :: iSm1_g1mg2,iSm1_g1mg2_fft,ig,gmg_sph,gmg_fft
1922 : integer :: igp,igstart,isym,itim,outofbox_wfn
1923 : complex(gwp) :: phmGt
1924 : !character(len=500) :: msg
1925 : !************************************************************************
1926 :
1927 1784 : igstart=1; if (qzero) igstart=2
1928 1784 : outofbox_wfn=0
1929 :
1930 1976 : SELECT CASE (symchi)
1931 :
1932 : CASE (0) ! Do not use symmetries.
1933 : ! MG: One has to make sure G1-G2 is still in the FFT mesh for each G1 and G2 in chi0 (not always true)
1934 : ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic!
1935 26432 : do igp=igstart,npwe
1936 3612608 : do ig=igstart,npwe
1937 3586176 : gmg_fft = Gsph_FFT%gmg_fftidx(ig,igp,ngfft)
1938 3586176 : if (gmg_fft==0) then
1939 0 : outofbox_wfn=outofbox_wfn+1; CYCLE
1940 : end if
1941 28715648 : chi0(ig,igp,:) = chi0(ig,igp,:) + wfwfg(gmg_fft)*green_enhigh_w(:)
1942 : end do
1943 : end do
1944 :
1945 : CASE (1)
1946 : ! Symmetrize the integrand in the full BZ.
1947 : ! * <Sk b|e^{-i(G1-G2}.r}|b Sk> = e^{-i(G1-G2).\tau} <k b|e^{-i(S^{-1}(G1-G2).r)|b k>
1948 : ! * green_enhigh_w in invariant under symmetry
1949 : ! * We symmetrize using the operations of the little group of q since this routine
1950 : ! is called inside a sum over IBZ_q, it would be possible to symmetrize
1951 : ! this term by just summing over the IBZ and rotating the matrix elements.
1952 : ! * Time-reversal does not lead to a complex conjugated since bra and ket are the same.
1953 : !
1954 101060 : do igp=igstart,npwe
1955 6523900 : do ig=igstart,npwe
1956 :
1957 : ! Get the index of G1-G2.
1958 6422840 : gmg_sph = Gsph_FFT%gmg_idx(ig,igp)
1959 6422840 : if (gmg_sph==0) then
1960 0 : outofbox_wfn=outofbox_wfn+1; CYCLE
1961 : end if
1962 :
1963 19367988 : do itim=1,Ltg_q%timrev
1964 327741864 : do isym=1,Ltg_q%nsym_sg
1965 321319024 : if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
1966 : ! * This operation belongs to the little group and has to be used to reconstruct the BZ.
1967 : ! * Time-reversal in not used to rotate (G1-G2) see comment above.
1968 17254456 : phmGt = Gsph_FFT%phmGt (gmg_sph,isym)
1969 17254456 : iSm1_g1mg2 = Gsph_FFT%rottbm1(gmg_sph,1,isym)
1970 17254456 : iSm1_g1mg2_fft = igfft0(iSm1_g1mg2)
1971 :
1972 51763368 : chi0(ig,igp,:) = chi0(ig,igp,:) + phmGt*wfwfg(iSm1_g1mg2_fft)*green_enhigh_w(:)
1973 : end if
1974 : end do !isym
1975 : end do !itim
1976 :
1977 : end do !igp
1978 : end do !ig
1979 :
1980 : CASE DEFAULT
1981 1784 : ABI_BUG("Wrong value of symchi")
1982 : END SELECT
1983 :
1984 1784 : if (outofbox_wfn/=0) then
1985 0 : enough=enough+1
1986 0 : if (enough<=50) then
1987 0 : ABI_WARNING(sjoin(' Number of G1-G2 pairs outside the G-sphere for Wfns: ', itoa(outofbox_wfn)))
1988 0 : if (enough==50) then
1989 0 : call wrtout(std_out,' ========== Stop writing Warnings ==========')
1990 : end if
1991 : end if
1992 : end if
1993 :
1994 1784 : end subroutine completechi0_deltapart
1995 : !!***
1996 :
1997 : !----------------------------------------------------------------------
1998 :
1999 : !!****f* m_chi0tk/output_chi0sumrule
2000 : !! NAME
2001 : !! output_chi0sumrule
2002 : !!
2003 : !! FUNCTION
2004 : !! Calculate and output the value of the sum rule for
2005 : !! the non-interacting polarizability chi0
2006 : !!
2007 : !! INPUTS
2008 : !!
2009 : !! OUTPUT
2010 : !! (for writing routines, no output)
2011 : !! otherwise, should be described
2012 : !!
2013 : !! SOURCE
2014 :
2015 922 : subroutine output_chi0sumrule(qeq0,iq,npwe,omegaplasma,chi0sumrule,epsm1_w0,vc_sqrt)
2016 :
2017 : !Arguments ------------------------------------
2018 : !scalars
2019 : integer,intent(in) :: iq,npwe
2020 : real(dp),intent(in) :: omegaplasma
2021 : logical,intent(in) :: qeq0
2022 : !arrays
2023 : real(dp),intent(inout) :: chi0sumrule(npwe)
2024 : complex(gwp),intent(in) :: epsm1_w0(npwe,npwe),vc_sqrt(npwe)
2025 :
2026 : !Local variables ------------------------------
2027 : !scalars
2028 : integer :: ig,igstart
2029 : real(dp) :: average,norm
2030 : character(len=500) :: msg
2031 : !************************************************************************
2032 :
2033 922 : igstart=1; if (qeq0) igstart=2
2034 : !
2035 : ! The sumrule reads:
2036 : ! $ \int d\omega \omega v * Im[ \chi_0(\omega) ] = \pi/2 * w_p^2 $.
2037 48764 : chi0sumrule(igstart:npwe) = chi0sumrule(igstart:npwe) * vc_sqrt(igstart:npwe)**2
2038 : !
2039 : ! Calculate a weighted average of the fulfilment of the sumrule on epsilon
2040 : ! The weight is given according to the significance of each q+G in the
2041 : ! subsequent GW calculation: It is proportional to v * (epsm1 -1 )
2042 : average = zero; norm = zero
2043 48764 : do ig=igstart,npwe
2044 47842 : average = average + chi0sumrule(ig) * real( vc_sqrt(ig)**2 * (epsm1_w0(ig,ig) - 1.0_dp ) )
2045 48764 : norm = norm + real( vc_sqrt(ig)**2 * (epsm1_w0(ig,ig) - 1.0_dp ) )
2046 : !average = average + chi0sumrule(ig) * real( (epsm1_w0(ig,ig) - 1.0_dp ) )
2047 : !norm = norm + real( (epsm1_w0(ig,ig) - 1.0_dp ) )
2048 : !write(203,'(i4,8(2x,e12.6))') ig,1.0_dp/vc_sqrt(ig),chi0sumrule(ig)/ (0.5d0*omegaplasma**2*pi)
2049 : end do
2050 :
2051 922 : if (abs(norm)>tol8) then
2052 : write(msg,'(1x,a,i4,a,f10.2,2x,a)')&
2053 922 : ' Average fulfillment of the sum rule on Im[epsilon] for q-point ',&
2054 1844 : iq,' :',average/norm/(0.5_dp*omegaplasma**2*pi)*100.0_dp,'[%]'
2055 2766 : call wrtout([std_out, ab_out], msg)
2056 : end if
2057 :
2058 922 : end subroutine output_chi0sumrule
2059 : !!***
2060 :
2061 : !----------------------------------------------------------------------
2062 :
2063 : !!****f* m_chi0tk/accumulate_chi0sumrule
2064 : !! NAME
2065 : !! accumulate_chi0sumrule
2066 : !!
2067 : !! FUNCTION
2068 : !! Accumulate the contribution to the sum rule for Im chi0
2069 : !! arising from a single transition. Eventually symmetrize
2070 : !! it using the symmetry operations of the little group of q.
2071 : !!
2072 : !! INPUTS
2073 : !! ik_bz=Index of the k-point in the full BZ whose contribution has to be added and symmetrized.
2074 : !! symchi=1 if we are summing over IBZ_q and symmetrization has to be performed.
2075 : !! npwe=Number of G vectors in chi0.
2076 : !! npwepG0=Number of G vectors in the "enlarged" sphere to treat umklapp.
2077 : !! factor=factor entering the expression.
2078 : !! delta_ene=Transition energy.
2079 : !! Ltg_q=<littlegroup_t>= Structure gathering information on the little group of the external q.
2080 : !! Gsph_epsG0=<gsphere_t>=Info on the G-sphere for chi0.
2081 : !! rhotwg(npwepG0)=Fouriet transform of u_{b1 k-q} u_{b2 k} in the "enlarged" sphere.
2082 : !!
2083 : !! OUTPUT
2084 : !! See SIDES EFFECTS
2085 : !!
2086 : !! SIDES EFFECTS
2087 : !! chi0sumrule(npwe)= In input the sum rule calculated so far,
2088 : !! In output the contribution of this transition is accounted for, and, eventually, symmetrized.
2089 : !! using the symmetry operations of the little group of the external q.
2090 : !!
2091 : !! SOURCE
2092 :
2093 1650394 : subroutine accumulate_chi0sumrule(ik_bz,symchi,npwe,factor,delta_ene,&
2094 1650394 : & Ltg_q,Gsph_epsG0,npwepG0,rhotwg,chi0sumrule)
2095 :
2096 : !Arguments ------------------------------------
2097 : !scalars
2098 : integer,intent(in) :: ik_bz,npwe,npwepG0,symchi
2099 : real(dp),intent(in) :: delta_ene,factor
2100 : type(gsphere_t),intent(in) :: Gsph_epsG0
2101 : type(littlegroup_t),target,intent(in) :: Ltg_q
2102 : !arrays
2103 : real(dp),intent(inout) :: chi0sumrule(npwe)
2104 : complex(gwp),intent(in) :: rhotwg(npwepG0)
2105 :
2106 : !Local variables-------------------------------
2107 : !scalars
2108 : integer :: isym,itim
2109 : !character(len=500) :: msg
2110 : !arrays
2111 1650394 : integer,allocatable :: Sm1_gmG0(:)
2112 1650394 : integer, contiguous, pointer :: gmG0(:)
2113 1650394 : complex(gwp),allocatable :: rhotwg_sym(:)
2114 : !************************************************************************
2115 :
2116 : ! Accumulating the sum rule on chi0.
2117 : ! Eq.(5.284) in G. D. Mahan Many-Particle Physics 3rd edition [[cite:Mahan2000]]
2118 :
2119 1976528 : SELECT CASE (symchi)
2120 : CASE (0)
2121 : ! Do not use symmetries, sum is performed in the full BZ.
2122 15288808 : chi0sumrule(:)=chi0sumrule(:) + factor*delta_ene*ABS(rhotwg(1:npwe))**2
2123 :
2124 : CASE (1)
2125 : ! Symmetrize the contribution in the full BZ.
2126 3972780 : ABI_MALLOC(rhotwg_sym,(npwe))
2127 3972780 : ABI_MALLOC(Sm1_gmG0,(npwe))
2128 :
2129 3968556 : do itim=1,Ltg_q%timrev
2130 114757676 : do isym=1,Ltg_q%nsym_sg
2131 113433416 : if (Ltg_q%wtksym(itim,isym,ik_bz)==1) then
2132 : ! This operation belongs to the little group and has to be used to reconstruct the BZ ===
2133 : ! In the following 2 lines mind the slicing (1:npwe)
2134 3160395 : gmG0 => Ltg_q%igmG0(1:npwe,itim,isym)
2135 212097948 : Sm1_gmG0(1:npwe)=Gsph_epsG0%rottbm1(gmG0(1:npwe),itim,isym)
2136 212097948 : rhotwg_sym(1:npwe)=rhotwg(Sm1_gmG0)
2137 :
2138 212097948 : chi0sumrule(:)=chi0sumrule(:) + factor*delta_ene*ABS(rhotwg_sym(1:npwe))**2
2139 : end if
2140 : end do !isym
2141 : end do !itim
2142 :
2143 1324260 : ABI_FREE(rhotwg_sym)
2144 1324260 : ABI_FREE(Sm1_gmG0)
2145 :
2146 : CASE DEFAULT
2147 0 : ABI_BUG(sjoin('Wrong value for symchi:', itoa(symchi)))
2148 : END SELECT
2149 :
2150 1650394 : end subroutine accumulate_chi0sumrule
2151 : !!***
2152 :
2153 : !!****f* m_chi0tk/make_transitions
2154 : !! NAME
2155 : !! make_transitions
2156 : !!
2157 : !! FUNCTION
2158 : !! Calculate transition energies entering the espression for the irreducible polarizability.
2159 : !!
2160 : !! INPUTS
2161 : !! nsspol=1 for spin unpolarized, 2 for spin polarized calculations
2162 : !! nbnds=total number of bands
2163 : !! kmesh<kmesh_t>=datatype gathering info on the k-mesh:
2164 : !! | %nbz=number of k-points in the full BZ
2165 : !! | %nibz=number of k-points in the IBZ
2166 : !! | %tab(nkbz)=table giving for each k-point in the BZ, the corresponding irreducible point in the IBZ array
2167 : !! | %bz(3,nkbz)=reduced coordinated of k-points
2168 : !! TOL_DELTA_OCC=tolerance on the difference of the occupation numbers
2169 : !! gw_energy(nbnds,kmesh%nkibz,nsppol)=quasi-particle energies energies
2170 : !! occ(nbnds,kmesh%nkibz,nsppol)=occupation numbers
2171 : !! chi0alg=integer defining the method used to calculate chi0
2172 : !! 0 ==> calculate chi0 using the Adler-Wiser expression
2173 : !! 1 ==> use spectral method
2174 : !! timrev=if 2, time-reversal symmetry is considered; 1 otherwise
2175 : !!
2176 : !! OUTPUT
2177 : !! my_max_rest,my_min_rest=Maximum and minimum resonant (posite) transition energy.
2178 : !! max_rest,min_rest=Maximun and minimum resonant (posite) transition energy treated by this node.
2179 : !!
2180 : !! SOURCE
2181 :
2182 54 : subroutine make_transitions(Wfd,chi0alg,nbnds,nbvw,nsppol,symchi,timrev,TOL_DELTA_OCC,&
2183 54 : max_rest,min_rest,my_max_rest,my_min_rest,Kmesh,Ltg_q,gw_energy,occ,qpoint,bbp_ks_distrb)
2184 :
2185 : !Arguments ------------------------------------
2186 : !scalars
2187 : integer,intent(in) :: chi0alg,nbnds,nbvw,nsppol,symchi,timrev
2188 : real(dp),intent(in) :: TOL_DELTA_OCC
2189 : real(dp),intent(out) :: max_rest,min_rest, my_max_rest,my_min_rest
2190 : type(kmesh_t),intent(in) :: Kmesh
2191 : type(littlegroup_t),intent(in) :: Ltg_q
2192 : type(wfdgw_t),intent(in) :: Wfd
2193 : !arrays
2194 : real(dp),intent(in) :: gw_energy(nbnds,Kmesh%nibz,nsppol)
2195 : real(dp),intent(in) :: occ(nbnds,Kmesh%nibz,nsppol),qpoint(3)
2196 : integer,intent(in) :: bbp_ks_distrb(Wfd%mband,Wfd%mband,Kmesh%nbz,Wfd%nsppol)
2197 :
2198 : !Local variables-------------------------------
2199 : !scalars
2200 : integer :: ib1,ib2,ii,ik_bz,ik_ibz,ikmq_bz,ikmq_ibz,is,nt,ntrans,my_ntrans,iloop
2201 : real(dp) :: delta_ene,delta_occ,spin_fact
2202 : character(len=500) :: msg
2203 : !arrays
2204 : integer :: G0(3)
2205 : real(dp) :: kmq(3)
2206 : !************************************************************************
2207 :
2208 : DBG_ENTER("COLL")
2209 :
2210 54 : if (chi0alg < 0 .or. chi0alg >= 2) then
2211 0 : ABI_BUG(sjoin('chi0alg:', itoa(chi0alg),' not allowed'))
2212 : end if
2213 54 : if (timrev /= 1 .and. timrev /= 2) then
2214 0 : ABI_BUG(sjoin('timrev:', itoa(timrev),' not allowed'))
2215 : end if
2216 :
2217 : ABI_UNUSED(nbvw)
2218 : !
2219 : ! In the first loop, we calculate total number of transitions for this q-point
2220 : ! as well the min and max transition without taking into account distribution of bands.
2221 : ! In the second iteration, we calculate the min and Max transition treated by this processor.
2222 : !
2223 54 : spin_fact =half; if (nsppol == 2) spin_fact = one
2224 :
2225 54 : my_max_rest=smallest_real; my_min_rest=greatest_real
2226 54 : max_rest=smallest_real; min_rest=greatest_real
2227 :
2228 162 : do iloop=1,2
2229 : nt=0
2230 16076 : do ik_bz=1,Kmesh%nbz
2231 15968 : ik_ibz=Kmesh%tab(ik_bz)
2232 63872 : kmq(:)=Kmesh%bz(:,ik_bz)-qpoint(:)
2233 :
2234 15968 : if (symchi == 1) then
2235 15968 : if (Ltg_q%ibzq(ik_bz) /= 1) cycle ! This point does not belong to the IBZ defined by the little group
2236 : end if
2237 :
2238 : ! Find kp=k-q-G0 and also G0 where kp is in the first BZ
2239 6604 : if (.not. kmesh%has_BZ_item(kmq,ikmq_bz,g0)) then
2240 : ! Stop as the weight 1.0/nkbz is wrong.
2241 0 : write(msg,'(4a,2(2a,3f12.6),2a)')ch10,&
2242 0 : ' make_transitions : ERROR - ',ch10,&
2243 0 : ' kp = k-q-G0 not found in the BZ mesh',ch10,&
2244 0 : ' k = ',(Kmesh%bz(ii,ik_bz),ii=1,3),ch10,&
2245 0 : ' k-q = ',(kmq(ii),ii=1,3),ch10,&
2246 0 : ' weight in cchi0/cchi0q is wrong '
2247 0 : ABI_ERROR(msg)
2248 : end if
2249 :
2250 6604 : ikmq_ibz = Kmesh%tab(ikmq_bz)
2251 13316 : do is=1,nsppol
2252 177540 : do ib1=1,nbnds
2253 4552268 : do ib2=1,nbnds
2254 :
2255 4390696 : if (iloop == 2) then
2256 2195348 : if (bbp_ks_distrb(ib1,ib2,ik_bz,is)/=Wfd%my_rank) cycle
2257 : end if
2258 :
2259 2453332 : if (timrev == 2 .and. ib1 < ib2) cycle ! Thanks to time-reversal we gain a factor ~2.
2260 :
2261 1394400 : delta_occ = spin_fact * (occ(ib1,ikmq_ibz,is) - occ(ib2,ik_ibz,is))
2262 1394400 : delta_ene = gw_energy(ib1,ikmq_ibz,is) - gw_energy(ib2,ik_ibz,is)
2263 :
2264 1394400 : if (chi0alg == 0) then
2265 : ! Adler-Wiser expression. Skip only if factor due to occupation number is smaller than TOL_DELTA_OCC
2266 0 : if (abs(delta_occ) < abs(TOL_DELTA_OCC)) cycle
2267 1394400 : else if (chi0alg==1) then
2268 : ! Spectral method with time-reversal, only resonant transitions
2269 : ! This has to be changed to include spectral method without time-reversal
2270 1394400 : if (delta_ene < -abs(TOL_DELTA_OCC) .or. abs(delta_occ) < abs(TOL_DELTA_OCC)) cycle
2271 : end if
2272 :
2273 : ! We have a new transition
2274 515968 : nt=nt+1
2275 :
2276 515968 : if (iloop==1) then
2277 257984 : max_rest=MAX(max_rest,zero,delta_ene)
2278 257984 : if (delta_ene>=-tol6) min_rest=MIN(min_rest,delta_ene)
2279 : end if
2280 670936 : if (iloop==2) then
2281 257984 : my_max_rest=MAX(my_max_rest,zero,delta_ene)
2282 257984 : if (delta_ene>=-tol6) my_min_rest=MIN(my_min_rest,delta_ene)
2283 : end if
2284 :
2285 : end do
2286 : end do
2287 : end do
2288 : end do
2289 108 : if (iloop==1) ntrans=nt
2290 108 : if (iloop==2) my_ntrans=nt
2291 : end do !iloop
2292 :
2293 54 : write(msg,'(2a,i9,2a,f8.3,3a,f8.3,a)')ch10,&
2294 54 : ' Total number of transitions = ',ntrans,ch10,&
2295 54 : ' min resonant = ',min_rest*Ha_eV,' [eV] ',ch10,&
2296 108 : ' Max resonant = ',max_rest*Ha_eV,' [eV] '
2297 54 : call wrtout(std_out, msg)
2298 :
2299 54 : if (Wfd%nproc/=1) then
2300 0 : write(msg,'(2a,i9,2a,f8.3,3a,f8.3,a)')ch10,&
2301 0 : ' Total number of transitions for this processor= ',my_ntrans,ch10,&
2302 0 : ' min resonant = ',my_min_rest*Ha_eV,' [eV] ',ch10,&
2303 0 : ' Max resonant = ',my_max_rest*Ha_eV,' [eV] '
2304 0 : call wrtout(std_out, msg)
2305 : end if
2306 :
2307 : DBG_EXIT("COLL")
2308 :
2309 54 : end subroutine make_transitions
2310 : !!***
2311 :
2312 : !----------------------------------------------------------------------
2313 :
2314 : !!****f* m_chi0tk/chi0_bbp_mask
2315 : !! NAME
2316 : !! chi0_bbp_mask
2317 : !!
2318 : !! FUNCTION
2319 : !!
2320 : !! INPUTS
2321 : !!
2322 : !! OUTPUT
2323 : !!
2324 : !! SOURCE
2325 :
2326 35061 : subroutine chi0_bbp_mask(ikmq_ibz, ik_ibz, spin, spin_fact, use_tr, &
2327 35061 : gwcomp, spmeth, chi_nband, mband, ebands, bbp_mask)
2328 :
2329 : !Arguments ------------------------------------
2330 : !scalars
2331 : integer,intent(in) :: spin,ik_ibz,ikmq_ibz,mband,gwcomp, spmeth, chi_nband
2332 : real(dp),intent(in) :: spin_fact
2333 : logical,intent(in) :: use_tr
2334 : type(ebands_t),target,intent(in) :: ebands
2335 : !arrays
2336 : logical,intent(out) :: bbp_mask(mband, mband)
2337 :
2338 : !Local variables-------------------------------
2339 : !scalars
2340 : integer :: ib1, ib2
2341 : real(dp) :: deltaeGW_b1kmq_b2k,deltaf_b1kmq_b2k,e_b1_kmq,f_b1_kmq
2342 : !arrays
2343 35061 : real(dp), contiguous, pointer :: qp_eig(:,:,:),qp_occ(:,:,:)
2344 : !************************************************************************
2345 :
2346 35061 : qp_eig => ebands%eig; qp_occ => ebands%occ
2347 19693077 : bbp_mask = .FALSE.
2348 :
2349 69669 : select case (gwcomp)
2350 : case (0)
2351 : ! Loop over "conduction" states.
2352 632182 : do ib1=1,chi_nband
2353 597574 : e_b1_kmq = qp_eig(ib1, ikmq_ibz, spin)
2354 597574 : f_b1_kmq = qp_occ(ib1, ikmq_ibz, spin)
2355 :
2356 : ! Loop over "valence" states.
2357 19484990 : do ib2=1,chi_nband
2358 18852808 : deltaf_b1kmq_b2k = spin_fact * (f_b1_kmq - qp_occ(ib2,ik_ibz,spin))
2359 18852808 : deltaeGW_b1kmq_b2k = e_b1_kmq - qp_eig(ib2,ik_ibz,spin)
2360 :
2361 597574 : select case (spmeth)
2362 : case (0)
2363 : ! Standard Adler-Wiser expression.
2364 16657460 : if (ABS(deltaf_b1kmq_b2k) >= GW_TOL_DOCC) then
2365 3346214 : bbp_mask(ib1, ib2) = .TRUE.
2366 3346214 : if (use_tr .and. ib1 < ib2) bbp_mask(ib1,ib2) = .FALSE. ! GAIN a factor ~2 thanks to time-reversal.
2367 : end if
2368 :
2369 : case (1,2)
2370 : ! Spectral method, WARNING time-reversal here is always assumed!
2371 2195348 : if (ABS(deltaf_b1kmq_b2k) >= GW_TOL_DOCC) then
2372 515968 : bbp_mask(ib1,ib2)=.TRUE.
2373 515968 : if (deltaeGW_b1kmq_b2k<zero) bbp_mask(ib1,ib2)=.FALSE. ! Only positive frequencies are needed for the Hilbert transform.
2374 : !$if (use_tr .and. ib1<ib2) bbp_mask(ib1,ib2)=.FALSE. ! GAIN a factor ~2 thanks to time-reversal.
2375 : end if
2376 :
2377 : case default
2378 18852808 : ABI_ERROR(sjoin("Wrong value for spmeth:", itoa(spmeth)))
2379 : end select
2380 : !write(std_out,*) "bbp_mask(ib1,ib2)",bbp_mask(ib1,ib2)
2381 : end do !ib2
2382 : end do !ib1
2383 :
2384 : case (1)
2385 : ! Extrapolar technique
2386 453 : ABI_CHECK(spmeth == 0, "Hilbert transform and extrapolar method are not compatible")
2387 :
2388 : ! Loop over "conduction" states.
2389 9771 : do ib1=1,chi_nband
2390 9318 : e_b1_kmq=qp_eig(ib1,ikmq_ibz,spin)
2391 9318 : f_b1_kmq= qp_occ(ib1,ikmq_ibz,spin)
2392 :
2393 : ! Loop over "valence" states.
2394 208087 : do ib2=1,chi_nband
2395 198316 : deltaf_b1kmq_b2k = spin_fact*(f_b1_kmq-qp_occ(ib2,ik_ibz,spin))
2396 198316 : deltaeGW_b1kmq_b2k= e_b1_kmq-qp_eig(ib2,ik_ibz,spin)
2397 :
2398 : ! When the completeness correction is used,
2399 : ! we need to also consider transitions with vanishing deltaf
2400 : ! Rangel: This is to compute chi in metals correctly with the extrapolar method.
2401 198316 : bbp_mask(ib1,ib2)=.TRUE.
2402 : !if (qp_occ(ib2,ik_ibz,is) < GW_TOL_DOCC) CYCLE
2403 207634 : if (qp_occ(ib2,ik_ibz,spin) < GW_TOL_DOCC .and. (ABS(deltaf_b1kmq_b2k) < GW_TOL_DOCC .or. ib1<ib2)) then
2404 161328 : bbp_mask(ib1,ib2)=.FALSE.
2405 : end if
2406 : end do
2407 : end do
2408 :
2409 : case default
2410 35061 : ABI_ERROR(sjoin("Wrong value of gwcomp:", itoa(gwcomp)))
2411 : end select
2412 :
2413 35061 : end subroutine chi0_bbp_mask
2414 : !!***
2415 :
2416 431543 : END MODULE m_chi0tk
2417 : !!***
|
