Branch data 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 : 1253116 : 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 : 2506232 : integer :: Sm1_gmG0(Ep%npwe)
132 : 1253116 : complex(gwp),allocatable :: rhotwg_sym(:,:)
133 : : ! *************************************************************************
134 : :
135 : : ABI_UNUSED(nspinor)
136 : :
137 : : nthreads = xomp_get_max_threads()
138 : :
139 : 1253116 : 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 [ + + + + : 121055384 : nsymop = count(Ltg_q%wtksym(:,:,ik_bz)==1)
+ + ]
179 [ - + + - : 3901664 : ABI_MALLOC(rhotwg_sym,(Ep%npwe,nsymop))
- + + - +
- + - - +
- + ]
180 : 975416 : isymop = 0
181 : :
182 : : ! Prepare all the rhotwg at once to use BLAS level 3 routines
183 [ + + ]: 41002072 : do isym=1,Ltg_q%nsym_sg
184 [ + + ]: 120907928 : do itim=1,Ltg_q%timrev
185 [ + + ]: 119932512 : 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 [ + - + + : 308718722 : Sm1_gmG0(1:Ep%npwe) = Gsph_epsG0%rottbm1( Ltg_q%igmG0(1:Ep%npwe,itim,isym), itim,isym)
+ + ]
192 : :
193 : 2123310 : isymop = isymop + 1
194 : 2052188 : SELECT CASE (itim)
195 : : CASE (1)
196 [ + + ]: 152026248 : rhotwg_sym(1:Ep%npwe,isymop) = rhotwg(Sm1_gmG0) * Gsph_epsG0%phmGt(1:Ep%npwe,isym)
197 : : CASE (2)
198 [ + + ]: 3394768 : rhotwg_sym(1:Ep%npwe,isymop) = GWPC_CONJG(rhotwg(Sm1_gmG0))*Gsph_epsG0%phmGt(1:Ep%npwe,isym)
199 : : CASE DEFAULT
200 [ + + - ]: 2123310 : 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 [ + + ]: 4354690 : 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 [ + + + + ]: 4354690 : if( ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic ) then
214 : 2420910 : dr=green_w(io)
215 : 2420910 : 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 [ - + ]: 975416 : ABI_FREE(rhotwg_sym)
223 : :
224 : : CASE DEFAULT
225 [ + + - ]: 1253116 : ABI_BUG(sjoin('Wrong symchi:', itoa(Ep%symchi)))
226 : : END SELECT
227 : :
228 : 1253116 : 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 : 104097 : subroutine accumulate_chi0_q0(is_metallic,ik_bz,isym_kbz,itim_kbz,gwcomp,nspinor,npwepG0,Ep,Cryst,Ltg_q,Gsph_epsG0,&
615 : 104097 : 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 : 104097 : integer,contiguous, pointer :: Sm1G(:)
643 : : complex(dp) :: mir_kbz(3)
644 : 104097 : complex(gwp),allocatable :: rhotwg_sym(:,:)
645 : 104097 : complex(gwp), contiguous, pointer :: phmGt(:)
646 : : !************************************************************************
647 : :
648 : : ABI_UNUSED(deltaf_b1b2)
649 : :
650 : 104097 : npwe = ep%npwe; nomega = ep%nomega
651 : :
652 : 104097 : 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 [ + + + + : 7154311 : nsymop = count(Ltg_q%wtksym(:,:,ik_bz) == 1)
+ + ]
718 [ - + + - : 245692 : ABI_MALLOC(rhotwg_sym, (npwe, nsymop))
- + + - +
- + - - +
- + ]
719 : 61423 : isymop = 0
720 : :
721 : : ! Loop over symmetries of the space group and time-reversal.
722 [ + + ]: 2425719 : do isym=1,Ltg_q%nsym_sg
723 [ + + ]: 7099015 : do itim=1,Ltg_q%timrev
724 : :
725 [ + + ]: 7037592 : 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 : 323181 : phmGt => Gsph_epsG0%phmGt (1:npwe, isym) ! In the 2 lines below note the slicing (1:npwe)
728 : 323181 : Sm1G => Gsph_epsG0%rottbm1(1:npwe, itim, isym)
729 : :
730 : 323181 : isymop = isymop + 1
731 : 310989 : select case (itim)
732 : : case (1)
733 [ + + ]: 25547236 : rhotwg_sym(1:npwe, isymop) = rhotwg(Sm1G(1:npwe)) * phmGt(1:npwe)
734 : : case (2)
735 [ + + ]: 584640 : rhotwg_sym(1:npwe,isymop) = conjg(rhotwg(Sm1G(1:npwe))) * phmGt(1:npwe)
736 : : case default
737 [ + + - ]: 323181 : 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 [ + + ]: 323181 : if (nspinor == 1) then
745 [ + + + + : 8963052 : 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 [ + + + + ]: 359757 : 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 [ + + ]: 1292724 : do idir=1,3
754 [ + + ]: 5144094 : do io=1,nomega
755 [ + + ]: 232872264 : chi0_uwing(:,io,idir) = chi0_uwing(:,io,idir) + green_w(io) * mir_kbz(idir) * conjg(rhotwg_sym(:,isymop))
756 [ + + ]: 232872264 : chi0_lwing(:,io,idir) = chi0_lwing(:,io,idir) + green_w(io) * rhotwg_sym(:,isymop) * conjg(mir_kbz(idir))
757 [ + + ]: 4820913 : 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 [ + + ]: 1606971 : do io=1,nomega
767 [ + + ]: 5458341 : do jdir=1,3
768 [ + + ]: 16689270 : do idir=1,3
769 : 11554110 : chi0_head(idir,jdir,io) = chi0_head(idir,jdir,io) + green_w(io) * mir_kbz(idir) * conjg(mir_kbz(jdir))
770 [ + + ]: 15405480 : 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 [ + + ]: 289184 : 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 [ + + + + ]: 289184 : if (ABS(AIMAG(green_w(io))) < 1.0e-6_dp .and. .not. is_metallic) then
788 : 161041 : dr = green_w(io)
789 : 161041 : 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 [ - + ]: 61423 : 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 [ + + - ]: 104097 : 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 : 931 : 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 [ + + ]: 931 : 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 [ + + ]: 48978 : 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 [ + + ]: 48978 : do ig=igstart,npwe
2044 : 48047 : average = average + chi0sumrule(ig) * real( vc_sqrt(ig)**2 * (epsm1_w0(ig,ig) - 1.0_dp ) )
2045 : 48978 : 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 [ + - ]: 931 : if (abs(norm)>tol8) then
2052 : : write(msg,'(1x,a,i4,a,f10.2,2x,a)')&
2053 : 931 : ' Average fulfillment of the sum rule on Im[epsilon] for q-point ',&
2054 : 1862 : iq,' :',average/norm/(0.5_dp*omegaplasma**2*pi)*100.0_dp,'[%]'
2055 [ + + ]: 2793 : call wrtout([std_out, ab_out], msg)
2056 : : end if
2057 : :
2058 : 931 : 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 : 1652089 : subroutine accumulate_chi0sumrule(ik_bz,symchi,npwe,factor,delta_ene,&
2094 : 1652089 : & 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 : 1652089 : integer,allocatable :: Sm1_gmG0(:)
2112 : 1652089 : integer, contiguous, pointer :: gmG0(:)
2113 : 1652089 : 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 : 1978223 : 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 [ - + + - : 3977865 : ABI_MALLOC(rhotwg_sym,(npwe))
+ - + - -
+ - + ]
2127 [ - + + - : 3977865 : ABI_MALLOC(Sm1_gmG0,(npwe))
+ - + - -
+ - + ]
2128 : :
2129 [ + + ]: 3973641 : do itim=1,Ltg_q%timrev
2130 [ + + ]: 114844121 : do isym=1,Ltg_q%nsym_sg
2131 [ + + ]: 113518166 : 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 : 3166083 : gmG0 => Ltg_q%igmG0(1:npwe,itim,isym)
2135 [ + + ]: 212244252 : Sm1_gmG0(1:npwe)=Gsph_epsG0%rottbm1(gmG0(1:npwe),itim,isym)
2136 [ + + ]: 212244252 : rhotwg_sym(1:npwe)=rhotwg(Sm1_gmG0)
2137 : :
2138 [ + + ]: 212244252 : chi0sumrule(:)=chi0sumrule(:) + factor*delta_ene*ABS(rhotwg_sym(1:npwe))**2
2139 : : end if
2140 : : end do !isym
2141 : : end do !itim
2142 : :
2143 [ - + ]: 1325955 : ABI_FREE(rhotwg_sym)
2144 [ - + ]: 1325955 : ABI_FREE(Sm1_gmG0)
2145 : :
2146 : : CASE DEFAULT
2147 : 0 : ABI_BUG(sjoin('Wrong value for symchi:', itoa(symchi)))
2148 : : END SELECT
2149 : :
2150 [ + + - ]: 1652089 : 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 : 35122 : subroutine chi0_bbp_mask(ikmq_ibz, ik_ibz, spin, spin_fact, use_tr, &
2327 : 35122 : 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 : 35122 : real(dp), contiguous, pointer :: qp_eig(:,:,:),qp_occ(:,:,:)
2344 : : !************************************************************************
2345 : :
2346 : 35122 : qp_eig => ebands%eig; qp_occ => ebands%occ
2347 [ + + + + ]: 19700948 : bbp_mask = .FALSE.
2348 : :
2349 : 69791 : select case (gwcomp)
2350 : : case (0)
2351 : : ! Loop over "conduction" states.
2352 [ + + ]: 632897 : do ib1=1,chi_nband
2353 : 598228 : e_b1_kmq = qp_eig(ib1, ikmq_ibz, spin)
2354 : 598228 : f_b1_kmq = qp_occ(ib1, ikmq_ibz, spin)
2355 : :
2356 : : ! Loop over "valence" states.
2357 [ + + ]: 19492861 : do ib2=1,chi_nband
2358 : 18859964 : deltaf_b1kmq_b2k = spin_fact * (f_b1_kmq - qp_occ(ib2,ik_ibz,spin))
2359 : 18859964 : deltaeGW_b1kmq_b2k = e_b1_kmq - qp_eig(ib2,ik_ibz,spin)
2360 : :
2361 : 598228 : select case (spmeth)
2362 : : case (0)
2363 : : ! Standard Adler-Wiser expression.
2364 [ + + ]: 16664616 : if (ABS(deltaf_b1kmq_b2k) >= GW_TOL_DOCC) then
2365 : 3349604 : bbp_mask(ib1, ib2) = .TRUE.
2366 [ + + + + ]: 3349604 : 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 [ + + - ]: 18859964 : 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 [ + + - ]: 35122 : ABI_ERROR(sjoin("Wrong value of gwcomp:", itoa(gwcomp)))
2411 : : end select
2412 : :
2413 : 35122 : end subroutine chi0_bbp_mask
2414 : : !!***
2415 : :
2416 [ - - + + : 432671 : END MODULE m_chi0tk
- - - - -
- + + ]
2417 : : !!***
|
