Branch data Line data Source code
1 : : !!****m* ABINIT/m_screen
2 : : !! NAME
3 : : !! m_screen
4 : : !!
5 : : !! FUNCTION
6 : : !! Screening object used in the BSE/GWPT code.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2014-2026 ABINIT group (MG)
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 : : !! TODO
15 : : !! Integrate these objects with quartic GW.
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : module m_screen
26 : :
27 : : use defs_basis
28 : : use m_xmpi
29 : : use m_hide_blas
30 : : use m_errors
31 : : use m_splines
32 : : use m_abicore
33 : : use m_kxc
34 : : use m_screening
35 : : use m_nctk
36 : : use m_sort
37 : : use m_yaml
38 : :
39 : : use m_gwdefs, only : GW_TOLQ0, czero_gw
40 : : use m_fstrings, only : firstchar, endswith, strcat, itoa, sjoin
41 : : use m_numeric_tools, only : print_arr
42 : : use m_geometry, only : normv
43 : : use m_crystal, only : crystal_t
44 : : use m_bz_mesh, only : kmesh_t
45 : : use m_gsphere, only : gsphere_t
46 : : use m_vcoul, only : vcoul_t
47 : : use m_io_screening, only : read_screening, hscr_t, ncname_from_id, em1_ncname
48 : : use m_ppmodel, only : ppmodel_t, PPM_NONE, PPM_NOTAB
49 : : use m_pstat, only : pstat_proc
50 : :
51 : : implicit none
52 : :
53 : : private
54 : :
55 : : public :: em1_symmetrize_ip
56 : : public :: em1_symmetrize_op
57 : :
58 : : ! Flags defining the content of the %mat buffer in the fgg_t type.
59 : : integer,public,parameter :: MAT_NOTYPE = 0
60 : : integer,public,parameter :: MAT_CHI0 = 1
61 : : integer,public,parameter :: MAT_CHI = 2
62 : : integer,public,parameter :: MAT_EPSILON = 3
63 : : integer,public,parameter :: MAT_INV_EPSILON = 4
64 : : integer,public,parameter :: MAT_INV_EPSILON_M1 = 5
65 : : integer,public,parameter :: MAT_W = 6
66 : : integer,public,parameter :: MAT_W_M1 = 7
67 : : !
68 : : ! Family vertex.
69 : : integer,public,parameter :: VTX_FAMILY_NONE = 0 ! No vertex correction.
70 : : integer,public,parameter :: VTX_FAMILY_TDDFT = 1 ! TDDFT-based vertex.
71 : : integer,public,parameter :: VTX_FAMILY_ADA = 2 ! ADA vertex.
72 : : !
73 : : ! Test charge or test particle.
74 : : integer,public,parameter :: VTX_TEST_CHARGE = 0
75 : : integer,public,parameter :: VTX_TEST_PARTICLE = 1
76 : : !
77 : : ! Named constants for the frequency mesh.
78 : : integer,private,parameter :: WMESH_LINEAR = 1
79 : : integer,private,parameter :: WMESH_GAUSS_LEG = 2
80 : : integer,private,parameter :: WMESH_TAN_GRID = 3
81 : : !
82 : : ! Method used for the frequency integration.
83 : : integer,public,parameter :: WINT_NONE = 0
84 : : integer,public,parameter :: WINT_PPMODEL = 1
85 : : integer,public,parameter :: WINT_CONTOUR = 2
86 : : integer,public,parameter :: WINT_AC = 3
87 : : !
88 : : ! Parameters used for the model dielectric function.
89 : : integer,public,parameter :: MDL_NONE = 0
90 : : integer,public,parameter :: MDL_BECHSTEDT = 1
91 : : !
92 : : ! Flags giving the status of the local buffers defined in fgg_t.
93 : : integer,private,parameter :: MAT_NODATA = 0
94 : : integer,private,parameter :: MAT_ALLOCATED = 1
95 : : integer,private,parameter :: MAT_STORED = 2
96 : : !
97 : : ! Flags giving the status of the local buffers defined in fgg_t.
98 : : integer,private,parameter :: FGG_QBZ_ISPOINTER =1 ! Fgg_qbz is used to store the address in memory.
99 : : integer,private,parameter :: FGG_QBZ_ISALLOCATED=2 ! Fgg_qbz is used as an allocable array.
100 : : !!***
101 : :
102 : : !----------------------------------------------------------------------
103 : :
104 : : !!****t* m_screen/screen_info_t
105 : : !! NAME
106 : : !! screen_info_t
107 : : !!
108 : : !! FUNCTION
109 : : !! Container storing the parameters used to initialize a screen_t datatype or to
110 : : !! calculate a new SCR file from the SUSC file containing the independent-particle
111 : : !! polarizability.
112 : : !!
113 : : !! NOTES
114 : : !! The list of parameters passed to screening_init is copied in W%Info.
115 : : !! At present there is no need to provide a copy method since the structure does
116 : : !! not contain pointers but such a method must be defined and used if
117 : : !! dynamic entities are added to the datatype.
118 : : !!
119 : : !! SOURCE
120 : :
121 : : type,public :: screen_info_t
122 : :
123 : : integer :: mat_type = MAT_NOTYPE
124 : : ! Matrix identifier. See MAT_* flags.
125 : :
126 : : integer :: vtx_family = VTX_FAMILY_NONE
127 : : ! Vertex correction family.
128 : :
129 : : integer :: invalid_freq = 0
130 : : ! Sets the procedure to follow when a ppm frequency is invalid (negative or imaginary),
131 : : ! see input variable gw_invalid_freq
132 : :
133 : : integer :: ixc = 0
134 : : ! XC functional used for the TDDFT-based vertex.
135 : :
136 : : integer :: use_ada = 0
137 : : ! >0 if ADA vertex is used.
138 : :
139 : : integer :: use_mdf = MDL_NONE
140 : : ! >0 if model dielectric function is used.
141 : :
142 : : integer :: use_ppm = PPM_NONE
143 : : ! >0 if ppmodel is used.
144 : :
145 : : integer :: vtx_test = VTX_TEST_CHARGE
146 : : ! test charge or test particle.
147 : :
148 : : integer :: wint_method = WINT_NONE
149 : : ! Defines the frequency integration technique. See WIN_ flags.
150 : : ! NOTE that this flag can be changed at run time. For example
151 : : ! one can switch from the CD to the PPm if the ppmodel parameters are in memory
152 : :
153 : : real(dp) :: ada_kappa = 2.1_dp
154 : : ! Inverse smearing length used for ADA.
155 : :
156 : : real(dp) :: eps_inf = 12.0_dp
157 : : ! Dielectric constant used for the model dielectric function.
158 : :
159 : : real(dp) :: drude_plsmf = zero
160 : : ! Drude plasma frequency used for PPmodel 1.
161 : :
162 : : contains
163 : :
164 : : procedure :: print => screen_info_print
165 : :
166 : : end type screen_info_t
167 : : !!***
168 : :
169 : : !----------------------------------------------------------------------
170 : :
171 : : !!****t* m_screen/fgg_t
172 : : !! NAME
173 : : !! fgg_t
174 : : !!
175 : : !! FUNCTION
176 : : !! Object used to store F(G,G')(w) for a given q-point.
177 : : !!
178 : : !! SOURCE
179 : :
180 : : type,public :: fgg_t
181 : :
182 : : integer :: nomega
183 : : ! Number of frequencies.
184 : :
185 : : integer :: npw
186 : : ! Number of G vectors.
187 : :
188 : : integer :: nqlwl
189 : : ! Number of points for the treatment of the long wave-length limit.
190 : :
191 : : integer :: has_mat = MAT_NODATA
192 : : ! Flag giving the status of mat.
193 : :
194 : : complex(gwp),allocatable :: mat(:,:,:)
195 : : ! (npw, npw, nomega)
196 : : ! The component of the two-point function $F_{G,G',w}$ for a given q.
197 : :
198 : : !complex(dp),allocatable :: head(:,:,:)
199 : : ! head(3,3,nomega)
200 : :
201 : : !complex(dp),allocatable :: lwing(:,:,:)
202 : : ! lwing(3,npwe,nomega)
203 : : ! Lower wings
204 : :
205 : : !complex(dp),allocatable :: uwing(:,:,:)
206 : : ! uwing(3,npwe,nomega)
207 : : ! Upper wings.
208 : :
209 : : contains
210 : :
211 : : procedure :: init => fgg_init ! Creation method.
212 : : !procedure :: free => fgg_free
213 : : end type fgg_t
214 : :
215 : : public :: fgg_free ! Free memory.
216 : : !!***
217 : :
218 : : interface fgg_free
219 : : module procedure fgg_free_0D
220 : : module procedure fgg_free_1D
221 : : end interface fgg_free
222 : :
223 : : !----------------------------------------------------------------------
224 : :
225 : : !!****t* m_screen/screen_t
226 : : !! NAME
227 : : !! screen_t
228 : : !!
229 : : !! FUNCTION
230 : : !! Object used to store the screening matrix in reciprocal space.
231 : : !!
232 : : !! SOURCE
233 : :
234 : : type,public :: screen_t
235 : :
236 : : ! scalars
237 : : integer :: iomode ! Flag defining the IO mode.
238 : : integer :: debug_level=0 ! Internal Flag defining the debug level.
239 : : integer :: mqmem ! =0 for out-of-core solution, =nqibz if entire matrix is stored in memory.
240 : : integer :: nI,nJ ! Number of components (rows,columns) in chi|eps^-1. (1,1) if collinear.
241 : : integer :: nqibz ! Number of q-points in the IBZ used.
242 : : integer :: nqlwl ! Number of points used for the treatment of the long wave-length limit.
243 : : integer :: nomega ! Total Number of frequencies used.
244 : : integer :: nomega_i ! Number of purely imaginary frequencies used.
245 : : integer :: nomega_r ! Number of real frequencies used.
246 : : integer :: npw ! Number of G vectors.
247 : : integer :: prtvol ! Verbosity level.
248 : : integer :: has_ppmodel ! 1 if PPmodel tables are stored.
249 : : integer :: has_fgg ! 1 if Fgg tables are stored.
250 : : integer :: nfftf_tot
251 : : integer :: nspden
252 : :
253 : : ! arrays
254 : : integer :: ngfftf(18) ! Info on the FFT mesh used for ae_rhor (used for the model dielectric function)
255 : :
256 : : real(dp),allocatable :: ae_rhor(:,:)
257 : : ! ae_rhor(nfft,nspden)
258 : : ! Density in real space used to construct the TDDFT kernel or the model dielectric function.
259 : : ! NOTE that ae_rhor is given on the dense mesh as it contains the PAW onsite contribution.
260 : :
261 : : character(len=fnlen) :: fname = ABI_NOFILE ! Name of the file used for the out-of-core solution.
262 : :
263 : : real(dp),allocatable :: qibz(:,:)
264 : : ! (3,nqibz)
265 : : ! q-points in reduced coordinates
266 : :
267 : : real(dp),allocatable :: qlwl(:,:)
268 : : ! (3,nqlwl)
269 : : ! q-points used for the long wave-length limit treatment.
270 : :
271 : : complex(dp),allocatable :: omega(:)
272 : : ! (nomega)
273 : : ! List of frequencies. Real frequencies are packed first.
274 : :
275 : : integer,allocatable :: gvec(:,:)
276 : : ! (3,npw)
277 : : ! G-vectors used to describe the two-point function (r.l.u.).
278 : :
279 : : logical,allocatable :: keep_qibz(:)
280 : : ! (nqibz)
281 : : ! Storage strategy: keep or not keep Em1(q) in memory.
282 : :
283 : : type(fgg_t),pointer :: Fgg(:) => null()
284 : : ! (nqibz)
285 : : ! F_{G,G'}(q,w) for q in the IBZ.
286 : :
287 : : integer :: fgg_qbz_stat = FGG_QBZ_ISPOINTER
288 : : ! Status of Fgg_qbz
289 : :
290 : : integer :: fgg_qbz_idx = 0
291 : : ! The index of the q-point in BZ pointed by Fgg_qbz. Used for debugging purpose.
292 : :
293 : : type(fgg_t),pointer :: Fgg_qbz => null()
294 : : ! Buffer used for storing F_GG' at the point q_bz in the BZ
295 : : ! If q_bz is in the IBZ, Fgg_qbz *points* to Fgg(iq_ibz)
296 : : ! If q_bz is not in the IBZ, Fgg_qbz is *allocated* and used to store the symmetrized matrix.
297 : :
298 : : type(ppmodel_t) :: ppm
299 : : ! Structure storing the plasmon-pole parameters.
300 : :
301 : : type(screen_info_t) :: info
302 : : ! Parameters used to construct the screening.
303 : :
304 : : contains
305 : :
306 : : procedure :: nullify => screen_nullify
307 : : ! Nullify all pointers before use.
308 : :
309 : : procedure :: init => screen_init
310 : : ! Creation method.
311 : :
312 : : procedure :: print => screen_print
313 : : ! Print info on object
314 : :
315 : : procedure :: free => screen_free
316 : : ! Free dynamic memory
317 : :
318 : : procedure :: rotate_iqbz => screen_rotate_iqbz
319 : : ! Prepare the object for applying W_qbz.
320 : :
321 : : procedure :: w0gemv => screen_w0gemv
322 : : ! Matrix vector multiplication \sum_{G'} F_{G,G') |u(G')>.
323 : :
324 : : procedure :: calc_sigc => screen_calc_sigc
325 : : ! Compute the frequency convolution.
326 : :
327 : : procedure :: ihave_fgg => screen_ihave_fgg
328 : : ! Inquire the processor whether it has a particular F_{GG')(q) and with which status.
329 : :
330 : : end type screen_t
331 : : !!***
332 : :
333 : : contains
334 : : !----------------------------------------------------------------------
335 : :
336 : : !!****f* m_screen/screen_info_print
337 : : !! NAME
338 : : !! screen_info_print
339 : : !!
340 : : !! FUNCTION
341 : : !! Printout object
342 : : !!
343 : : !! INPUTS
344 : : !! units=Unit numbers for output
345 : : !! header=String to be printed as header for additional info.
346 : : !!
347 : : !! SOURCE
348 : :
349 : 58 : subroutine screen_info_print(W_info, units, header)
350 : :
351 : : !Arguments ------------------------------------
352 : : !scalars
353 : : class(screen_info_t),intent(in) :: W_info
354 : : integer,intent(in) :: units(:)
355 : : character(len=*),optional,intent(in) :: header
356 : :
357 : : !Local variables-------------------------------
358 : : character(len=500) :: msg
359 : : ! *********************************************************************
360 : :
361 : 29 : msg=' ==== Info on the screen_info_t% object ==== '
362 [ + - + - : 29 : if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
+ - ]
363 : 29 : call wrtout(units, msg)
364 : :
365 : : !integer
366 : 29 : write(msg,'(a,i3)')" mat_type ",W_info%mat_type
367 : 29 : call wrtout(units, msg)
368 : 29 : write(msg,'(a,i3)')" vtx_family ",W_info%vtx_family
369 : 29 : call wrtout(units, msg)
370 : 29 : write(msg,'(a,i3)')" invalid_freq",W_info%invalid_freq
371 : 29 : call wrtout(units, msg)
372 : 29 : write(msg,'(a,i3)')" ixc ",W_info%ixc
373 : 29 : call wrtout(units, msg)
374 : 29 : write(msg,'(a,i3)')" use_ada ",W_info%use_ada
375 : 29 : call wrtout(units, msg)
376 : 29 : write(msg,'(a,i3)')" use_mdf ",W_info%use_mdf
377 : 29 : call wrtout(units, msg)
378 : 29 : write(msg,'(a,i3)')" use_ppm ",W_info%use_ppm
379 : 29 : call wrtout(units, msg)
380 : 29 : write(msg,'(a,i3)')" vtx_test ",W_info%vtx_test
381 : 29 : call wrtout(units, msg)
382 : 29 : write(msg,'(a,i3)')" wint_method ",W_info%wint_method
383 : 29 : call wrtout(units, msg)
384 : :
385 : : !real
386 : 29 : write(msg,'(a,f8.3)')" ada_kappa ",W_info%ada_kappa
387 : 29 : call wrtout(units, msg)
388 : 29 : write(msg,'(a,f8.3)')" eps_inf ",W_info%eps_inf
389 : 29 : call wrtout(units, msg)
390 : 29 : write(msg,'(a,f8.3)')" drude_plsmf ",W_info%drude_plsmf
391 : 29 : call wrtout(units, msg)
392 : :
393 [ - + ]: 29 : end subroutine screen_info_print
394 : : !!***
395 : :
396 : : !----------------------------------------------------------------------
397 : :
398 : : !!****f* m_screen/fgg_free_0D
399 : : !! NAME
400 : : !! fgg_free_0D
401 : : !!
402 : : !! FUNCTION
403 : : !! Free dynamic memory.
404 : : !!
405 : : !! SOURCE
406 : :
407 : 2155 : subroutine fgg_free_0D(Fgg)
408 : :
409 : : !Arguments ------------------------------------
410 : : type(fgg_t),intent(inout) :: Fgg
411 : : ! *************************************************************************
412 : :
413 [ + + ]: 2155 : ABI_SFREE(Fgg%mat)
414 : 2155 : Fgg%has_mat = MAT_NODATA
415 : :
416 : 2155 : end subroutine fgg_free_0D
417 : : !!***
418 : :
419 : : !----------------------------------------------------------------------
420 : :
421 : : !!****f* m_screen/fgg_free_1D
422 : : !! NAME
423 : : !! fgg_free_1D
424 : : !!
425 : : !! FUNCTION
426 : : !! Deallocate all the memory.
427 : : !!
428 : : !! INPUT
429 : : !! [keep_qibz(:)]=Optional logical mask used to select the q-points that are deallocated.
430 : : !!
431 : : !! SOURCE
432 : :
433 [ - + - + : 29 : subroutine fgg_free_1D(Fgg, keep_qibz)
- - - - ]
434 : :
435 : : !Arguments ------------------------------------
436 : : !scalars
437 : : type(fgg_t),intent(inout) :: Fgg(:)
438 : : logical,optional,intent(in) :: keep_qibz(:)
439 : :
440 : : !Local variables ------------------------------
441 : : !scalars
442 : : integer :: iq_ibz
443 : : logical :: keep_it
444 : : ! *************************************************************************
445 : :
446 [ + - + + ]: 234 : do iq_ibz=LBOUND(Fgg,DIM=1),UBOUND(Fgg,DIM=1)
447 [ - + ]: 176 : keep_it = .FALSE.; if (PRESENT(keep_qibz)) keep_it = keep_qibz(iq_ibz)
448 [ - - ]: 205 : if (.not. keep_it) call fgg_free_0D(Fgg(iq_ibz))
449 : : end do
450 : :
451 : 29 : end subroutine fgg_free_1D
452 : : !!***
453 : :
454 : : !----------------------------------------------------------------------
455 : :
456 : : !!****f* m_screen/fgg_init
457 : : !! NAME
458 : : !! fgg_init
459 : : !!
460 : : !! FUNCTION
461 : : !! Initialize the structure allocating the memory and initializing the internal variables.
462 : : !!
463 : : !! INPUT
464 : : !! npw
465 : : !! nqlwl
466 : : !!
467 : : !! SOURCE
468 : :
469 : 2149 : subroutine fgg_init(Fgg, npw, nomega, nqlwl)
470 : :
471 : : !Arguments ------------------------------------
472 : : !scalars
473 : : class(fgg_t),intent(inout) :: Fgg
474 : : integer,intent(in) :: npw, nqlwl, nomega
475 : :
476 : : !Local variables ------------------------------
477 : : integer :: ierr
478 : : ! *************************************************************************
479 : :
480 : 2149 : Fgg%nomega = nomega; Fgg%npw = npw; Fgg%nqlwl = nqlwl
481 : :
482 [ + - + - ]: 2149 : if (npw > 0 .and. nomega > 0) then
483 [ - + + - : 10745 : ABI_MALLOC_OR_DIE(Fgg%mat, (npw, npw, nomega), ierr)
- + + - -
+ + - + -
+ - + - -
+ + - -
+ ]
484 : 2149 : Fgg%has_mat = MAT_ALLOCATED
485 : : end if
486 : :
487 : 2149 : end subroutine fgg_init
488 : : !!***
489 : :
490 : : !----------------------------------------------------------------------
491 : :
492 : : !!****f* m_screen/screen_fgg_qbz_set
493 : : !! NAME
494 : : !! screen_fgg_qbz_set
495 : : !!
496 : : !! FUNCTION
497 : : !! Helper function used to perform the setup W%Fgg_qbz setting also the internal
498 : : !! flag that defines its status.
499 : : !!
500 : : !! INPUTS
501 : : !! iq_bz=Index of the q-point in the BZ.
502 : : !! nqlwl=Number of wings wanted.
503 : : !! how= "Pointer" is a true pointer is wanted.
504 : : !! "Allocated" if memory has to be allocated.
505 : : !!
506 : : !! NOTES
507 : : !! iq_bz and nqlwl are not used if how="Pointer".
508 : : !!
509 : : !! SOURCE
510 : :
511 : 22740 : subroutine screen_fgg_qbz_set(screen, iq_bz, nqlwl, how)
512 : :
513 : : !Arguments ------------------------------------
514 : : !scalars
515 : : class(screen_t),intent(inout) :: screen
516 : : integer,intent(in) :: iq_bz,nqlwl
517 : : character(len=*),intent(in) :: how
518 : :
519 : : !Local variables ------------------------------
520 : : !scalars
521 : : !character(len=500) :: msg
522 : : !************************************************************************
523 : :
524 : : ! Save the index of the q-point in the BZ.
525 : 22740 : screen%fgg_qbz_idx = iq_bz
526 : :
527 [ + + + + ]: 45480 : if (firstchar(how, (/"P"/)) ) then
528 : : ! We want a pointer.
529 : 4607 : select case (screen%fgg_qbz_stat)
530 : : case (FGG_QBZ_ISALLOCATED)
531 : 1978 : call fgg_free_0D(screen%Fgg_qbz)
532 [ - + - + ]: 1978 : ABI_FREE(screen%Fgg_qbz)
533 : : nullify(screen%Fgg_qbz)
534 : 1978 : screen%fgg_qbz_stat = FGG_QBZ_ISPOINTER
535 : :
536 : : case (FGG_QBZ_ISPOINTER)
537 : : ! Set it to null().
538 : 651 : nullify(screen%Fgg_qbz)
539 : :
540 : : case default
541 [ + + - ]: 2629 : ABI_ERROR(sjoin("Wrong status:", itoa(screen%fgg_qbz_stat)))
542 : : end select
543 : :
544 [ + + + - ]: 40222 : else if (firstchar(how, (/"A"/)) ) then
545 : : ! We want an allocatable array.
546 : :
547 : 22090 : select case (screen%fgg_qbz_stat)
548 : : case (FGG_QBZ_ISPOINTER)
549 : : ! Allocate memory
550 : : nullify(screen%Fgg_qbz)
551 [ - + ]: 1979 : ABI_MALLOC(screen%Fgg_qbz,)
552 : :
553 : 1979 : call screen%Fgg_qbz%init(screen%npw, screen%nomega, nqlwl)
554 : 1979 : screen%fgg_qbz_stat = FGG_QBZ_ISALLOCATED
555 : :
556 : : case (FGG_QBZ_ISALLOCATED)
557 : 18132 : screen%Fgg_qbz%has_mat = MAT_ALLOCATED ! STORED --> ALLOCATED
558 : :
559 : : case default
560 [ + + - ]: 20111 : ABI_ERROR(sjoin("Wrong status:", itoa(screen%fgg_qbz_stat)))
561 : : end select
562 : :
563 : : else
564 : 0 : ABI_BUG(sjoin("Wrong how:", how))
565 : : end if
566 : :
567 : 22740 : end subroutine screen_fgg_qbz_set
568 : : !!***
569 : :
570 : : !----------------------------------------------------------------------
571 : :
572 : : !!****f* m_screen/screen_ihave_fgg
573 : : !! NAME
574 : : !! screen_ihave_fgg
575 : : !!
576 : : !! FUNCTION
577 : : !! Inquire the processor whether it has a particular F_{GG')(q) and with which status.
578 : : !!
579 : : !! INPUTS
580 : : !! iq_ibz=k-point index
581 : : !! [how]=string defining which status is checked. By default the function returns
582 : : !! .TRUE. if the wave is either MAT_ALLOCATED or MAT_STORED.
583 : : !! Possible mutually exclusive values: "Allocated", "Stored".
584 : : !! Only the first character is checked (no case-sensitive)
585 : : !!
586 : : !! NOTES
587 : : !! A zero index can be used to inquire the status of the full set of q-points.
588 : : !!
589 : : !! SOURCE
590 : :
591 : 23268 : logical pure function screen_ihave_fgg(screen, iq_ibz, how)
592 : :
593 : : !Arguments ------------------------------------
594 : : !scalars
595 : : class(screen_t),intent(in) :: screen
596 : : integer,intent(in) :: iq_ibz
597 : : character(len=*),optional,intent(in) :: how
598 : :
599 : : !Local variables ------------------------------
600 : : integer :: ii, check(2)
601 : : !************************************************************************
602 : :
603 : 23268 : check = [MAT_ALLOCATED, MAT_STORED]
604 [ + - ]: 23268 : if (PRESENT(how)) then
605 [ + + + + ]: 69804 : if (firstchar(how, (/"A","a"/))) check = [MAT_ALLOCATED, MAT_ALLOCATED]
606 [ + + + + ]: 69804 : if (firstchar(how, (/"S","s"/))) check = [MAT_STORED, MAT_STORED]
607 : : end if
608 : :
609 [ + - ]: 23268 : if (iq_ibz > 0) then
610 : : screen_ihave_fgg = (screen%Fgg(iq_ibz)%has_mat == check(1) .or.&
611 [ + + - + ]: 23268 : screen%Fgg(iq_ibz)%has_mat == check(2) )
612 : : else
613 : : ! check the status of the full set of q-tables.
614 : 0 : screen_ihave_fgg=.TRUE.
615 [ # # ]: 0 : do ii=1,screen%nqibz
616 : : screen_ihave_fgg = screen_ihave_fgg .and. &
617 : : (screen%Fgg(ii)%has_mat == check(1) .or.&
618 [ # # # # : 0 : screen%Fgg(ii)%has_mat == check(2) )
# # ]
619 : : end do
620 : : end if
621 : :
622 : 23268 : end function screen_ihave_fgg
623 : : !!***
624 : :
625 : : !----------------------------------------------------------------------
626 : :
627 : : !!****f* m_screen/screen_nullify
628 : : !! NAME
629 : : !! screen_nullify
630 : : !!
631 : : !! FUNCTION
632 : : !! Initialize the pointers to null()
633 : : !!
634 : : !! SOURCE
635 : :
636 : 58 : subroutine screen_nullify(screen)
637 : :
638 : : !Arguments ------------------------------------
639 : : class(screen_t),intent(inout) :: screen
640 : : ! *************************************************************************
641 : :
642 : 58 : nullify(screen%Fgg_qbz); screen%fgg_qbz_stat=FGG_QBZ_ISPOINTER ! Needed since the initial status is undefined.
643 : 58 : nullify(screen%Fgg)
644 : :
645 : 58 : end subroutine screen_nullify
646 : : !!***
647 : :
648 : : !----------------------------------------------------------------------
649 : :
650 : : !!****f* m_screen/screen_free
651 : : !! NAME
652 : : !! screen_free
653 : : !!
654 : : !! FUNCTION
655 : : !! Free the memory allocated in the datatype.
656 : : !!
657 : : !! SOURCE
658 : :
659 : 29 : subroutine screen_free(screen)
660 : :
661 : : !Arguments ------------------------------------
662 : : class(screen_t),intent(inout) :: screen
663 : : ! *************************************************************************
664 : :
665 : : ! integer
666 [ + - ]: 29 : ABI_SFREE(screen%gvec)
667 : :
668 : : !real
669 [ + - ]: 29 : ABI_SFREE(screen%ae_rhor)
670 [ + - ]: 29 : ABI_SFREE(screen%qibz)
671 [ + - ]: 29 : ABI_SFREE(screen%qlwl)
672 : :
673 : : !complex
674 [ + - ]: 29 : ABI_SFREE(screen%omega)
675 : :
676 : : ! logical
677 [ + - ]: 29 : ABI_SFREE(screen%keep_qibz)
678 : :
679 : : ! types
680 : : ! Here be careful with dangling pointers.
681 : : ! First Fgg_qbz that might point to one of the %Fgg then %Fgg.
682 : 30 : select case (screen%fgg_qbz_stat)
683 : :
684 : : case (FGG_QBZ_ISALLOCATED)
685 : 1 : call fgg_free_0D(screen%Fgg_qbz)
686 [ - + - + ]: 1 : ABI_FREE(screen%Fgg_qbz)
687 : : nullify(screen%Fgg_qbz)
688 : 1 : screen%fgg_qbz_stat=FGG_QBZ_ISPOINTER
689 : :
690 : : case (FGG_QBZ_ISPOINTER)
691 : 28 : nullify(screen%Fgg_qbz)
692 : 28 : screen%fgg_qbz_stat=FGG_QBZ_ISPOINTER
693 : :
694 : : case default
695 [ + + - ]: 29 : continue
696 : : end select
697 : :
698 : : ! Free the Fgg matrices.
699 [ + - ]: 29 : if (associated(screen%Fgg)) then
700 : 29 : call fgg_free(screen%Fgg)
701 [ - + ]: 29 : ABI_FREE(screen%Fgg)
702 : : end if
703 : :
704 : : ! Free the plasmon pole tables.
705 : 29 : call screen%ppm%free()
706 : :
707 : 29 : end subroutine screen_free
708 : : !!***
709 : :
710 : : !----------------------------------------------------------------------
711 : :
712 : : !!****f* m_screen/screen_print
713 : : !! NAME
714 : : !! screen_print
715 : : !!
716 : : !! FUNCTION
717 : : !! Print info on the object.
718 : : !!
719 : : !! SOURCE
720 : :
721 : 0 : subroutine screen_print(screen, units, header)
722 : :
723 : : !Arguments ------------------------------------
724 : : class(screen_t),intent(in) :: screen
725 : : integer,intent(in) :: units(:)
726 : : character(len=*),optional,intent(in) :: header
727 : :
728 : : !Local variables-------------------------------
729 : : character(len=500) :: msg
730 : : type(yamldoc_t) :: ydoc
731 : : ! *************************************************************************
732 : :
733 [ # # # # ]: 0 : msg = ' ==== Info on the screen_t object ==== '; if (present(header)) msg=' ==== '//trim(adjustl(header))//' ==== '
734 : 0 : call wrtout(units, msg)
735 : :
736 : 0 : ydoc = yamldoc_open('screen_params') !, width=11, real_fmt='(3f8.3)')
737 : 0 : call ydoc%add_int("nomega_r", screen%nomega_r)
738 : 0 : call ydoc%add_int("nomega_i", screen%nomega_i)
739 : : !call ydoc%add_real("drude_plsmf", ppm%drude_plsmf)
740 : : !call ydoc%add_int1d("has_qibz", ppm%has_qibz)
741 : :
742 : 0 : call ydoc%write_units_and_free(units)
743 : :
744 [ # # # # ]: 0 : end subroutine screen_print
745 : : !!***
746 : :
747 : : !----------------------------------------------------------------------
748 : :
749 : : !!****f* m_screen/screen_init
750 : : !! NAME
751 : : !! screen_init
752 : : !!
753 : : !! FUNCTION
754 : : !! Initialize basic dimensions and other important arrays
755 : : !! starting from a file containing either epsilon^{-1} (_SCR) or chi0 (_SUSC).
756 : : !!
757 : : !! INPUTS
758 : : !! W_Info<screen_info_t>=The list of parameters used to construct the screen function.
759 : : !! Cryst<crystal_t>=Info on the unit cell.
760 : : !! Qmesh<kmesh_t>=Info on the Q-mesh.
761 : : !! Gsph<gsphere_t>=Info on the plane-wave basis set used for the two-point function.
762 : : !! Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction.
763 : : !! ifname=The name of the external file used to read the matrix.
764 : : !! id_required=Identifier used to specify the type of two-point function that is wanted.
765 : : !! iomode=Option defining the file format of the external file.
766 : : !! mqmem=0 for out-of-core solution, /=0 if entire matrix has to be stored in memory.
767 : : !! npw_asked=Number of G-vector to be used in the calculation, if <=0 use Max allowed number.
768 : : !! ngfftf(18)=Info on the (fine) mesh used for the density.
769 : : !! nfftf_tot=Total number of point in the FFT mesh for ae_rhor
770 : : !! nsppol=Number of independent spin polarizations.
771 : : !! nspden=Number of spin density components in ae_rhor
772 : : !! ae_rhor(nfftf_tot,nspden)
773 : : !! prtvol=Verbosity level.
774 : : !! comm=MPI communicator.
775 : : !!
776 : : !! SOURCE
777 : :
778 : 29 : subroutine screen_init(screen, W_Info, Cryst, Qmesh, Gsph, Vcp, ifname, mqmem, npw_asked, &
779 : 29 : iomode, ngfftf, nfftf_tot, nsppol, nspden, ae_rhor, prtvol, comm)
780 : :
781 : : !Arguments ------------------------------------
782 : : !scalars
783 : : class(screen_t),intent(out) :: screen
784 : : integer,intent(in) :: mqmem,iomode,npw_asked,comm,prtvol,nsppol, nfftf_tot,nspden
785 : : character(len=fnlen),intent(in) :: ifname
786 : : type(crystal_t),intent(in) :: Cryst
787 : : type(gsphere_t),intent(in) :: Gsph
788 : : type(vcoul_t),intent(in) :: Vcp
789 : : type(kmesh_t),intent(in) :: Qmesh
790 : : type(screen_info_t),intent(in) :: W_Info
791 : : !arrays
792 : : integer,intent(in) :: ngfftf(18)
793 : : real(dp),intent(in) :: ae_rhor(nfftf_tot,nspden)
794 : :
795 : : !Local variables-------------------------------
796 : : !scalars
797 : : integer,parameter :: master=0
798 : : integer :: option_test,approx_type,ixc_required,id_required !,nkxc
799 : : integer :: fform,my_rank,mat_type_read, nqibz,nomega,iq_ibz,npw,nqlwl
800 : : integer :: nI,nJ,iq_bz,mdf_type,ppmodel,ierr, iw,qsort,ii
801 : : real(dp) :: eps_inf,drude_plsmf
802 : : logical :: free_Fgg,found,from_file,is_qeq0,remove_dgg !only_one_kpt,
803 : : character(len=500) :: msg
804 : : character(len=fnlen) :: sus_fname,scr_fname
805 : : character(len=nctk_slen) :: varname
806 [ + + ]: 116 : type(hscr_t) :: Hscr
807 : : !arrays
808 : 58 : integer :: units(2), g0(3), iperm(Qmesh%nibz)
809 [ + - ]: 29 : real(dp) :: wt_list(Qmesh%nibz)
810 : : !complex(gwp),contiguous, pointer :: em1_ggw(:,:,:)
811 : : ! *********************************************************************
812 : :
813 : : DBG_ENTER("COLL")
814 : : ABI_UNUSED(nsppol)
815 : :
816 : 29 : my_rank = xmpi_comm_rank(comm); units = [std_out, ab_out]
817 : 29 : call screen%nullify()
818 : :
819 : : ! Initialize basic parameters
820 : 29 : screen%info = w_info
821 [ + + ]: 58 : call screen%info%print([std_out], header="W info")
822 : :
823 : 29 : id_required = W_Info%mat_type
824 : 29 : approx_type = W_Info%vtx_family
825 : 29 : option_test = W_Info%vtx_test
826 : 29 : ixc_required = W_Info%ixc
827 : 29 : varname = ncname_from_id(id_required)
828 : :
829 [ - + ]: 29 : if (all(id_required /= [MAT_INV_EPSILON])) then
830 : 0 : ABI_ERROR(sjoin("id_required:", itoa(id_required), " not available"))
831 : : end if
832 : :
833 : : ! This part must be rationalized.
834 : 29 : remove_dgg = (id_required == MAT_W_M1)
835 : :
836 [ + + ]: 29 : if (screen%info%use_mdf == MDL_NONE) screen%fname = ifname
837 : 29 : screen%nI = 1; screen%nJ = 1
838 : :
839 : : ! The q-point sampling is initialized from qmesh.
840 : 29 : screen%nqibz = Qmesh%nibz
841 [ - + + - : 87 : ABI_MALLOC(screen%qibz, (3, screen%nqibz))
+ - + - -
+ - + -
+ ]
842 [ + - + - : 762 : screen%qibz= Qmesh%ibz
- + - - -
- - - + +
+ + ]
843 : :
844 [ + + ]: 29 : screen%mqmem = mqmem; if (screen%mqmem /= 0) screen%mqmem = screen%nqibz !; screen%mqmem = 0
845 : :
846 [ - + + - : 87 : ABI_MALLOC(screen%keep_qibz, (screen%nqibz))
+ - + - -
+ - + -
+ ]
847 [ + + + + : 211 : screen%keep_qibz = .TRUE.; if (screen%mqmem == 0) screen%keep_qibz = .False.
+ + ]
848 : :
849 [ + + - + ]: 29 : if (screen%mqmem /= 0 .and. screen%mqmem < screen%nqibz) then
850 : : ! Keep in memory the most representative q-points.
851 [ # # ]: 0 : screen%keep_qibz = .FALSE.
852 [ # # # # : 0 : wt_list = Qmesh%wt; iperm = (/(ii,ii=1,Qmesh%nibz)/)
# # # # ]
853 : 0 : call sort_dp(Qmesh%nibz, wt_list, iperm, tol12)
854 [ # # ]: 0 : do qsort=Qmesh%nibz,Qmesh%nibz-mqmem+1,1
855 : 0 : iq_ibz = iperm(qsort)
856 : 0 : screen%keep_qibz(iq_ibz) = .TRUE.
857 : : end do
858 : : end if
859 : :
860 : 29 : screen%fgg_qbz_idx = 0
861 : 29 : screen%iomode = iomode
862 : 29 : screen%prtvol = prtvol
863 [ - + ]: 29 : screen%has_ppmodel = 0; if (screen%info%use_ppm /= PPM_NONE) screen%has_ppmodel = 1
864 : :
865 : : ! Copy the AE density for the model dielectric function or for the vertex corrections.
866 : 29 : screen%nspden = nspden
867 [ + + ]: 551 : screen%ngfftf = ngfftf
868 : 29 : screen%nfftf_tot = nfftf_tot
869 : :
870 [ - + + - : 116 : ABI_MALLOC(screen%ae_rhor, (nfftf_tot, nspden))
- + + - +
- + - - +
- + - + ]
871 [ - + + + : 286992 : screen%ae_rhor = ae_rhor
+ + ]
872 : :
873 : 29 : free_Fgg = .FALSE.
874 : :
875 [ - + ]: 29 : screen%has_fgg = 0; if (ANY(screen%info%wint_method == [WINT_CONTOUR, WINT_AC])) screen%has_fgg = 1
876 : :
877 [ - + - - ]: 29 : if (screen%has_fgg > 0 .and. screen%has_ppmodel > 0) then
878 : 0 : ABI_WARNING("Both PPmodel tables and F_(GG')(q,w) are stored in memory")
879 : : end if
880 : :
881 : : ! Default values used if external file is not read.
882 : 29 : nqlwl = 0; nomega = 1
883 : 29 : screen%npw = npw_asked
884 : :
885 : : ! Model dielectric function does not require any external file.
886 : 29 : from_file = (screen%info%use_mdf == MDL_NONE)
887 : :
888 [ + + ]: 29 : if (from_file) then
889 : : ! Open file and check its content.
890 [ + - ]: 11 : if (endswith(screen%fname, ".nc")) screen%iomode = IO_MODE_ETSF
891 : 11 : call hscr%from_file(screen%fname, fform, comm)
892 : : ! Echo of the header
893 [ + - - + : 11 : if (my_rank == master .and. screen%prtvol > 0) call hscr%print([std_out], 0)
- - ]
894 : :
895 : 11 : mat_type_read = Hscr%id
896 : 11 : nqlwl = Hscr%nqlwl
897 : 11 : nomega = Hscr%nomega
898 : : end if
899 : :
900 : 29 : screen%nqlwl = nqlwl
901 : 29 : screen%nomega = nomega
902 : :
903 [ + + + - : 87 : ABI_MALLOC(screen%qlwl, (3, screen%nqlwl))
+ - + + -
+ - + -
+ ]
904 [ - + + - : 87 : ABI_MALLOC(screen%omega, (screen%nomega))
+ - + - -
+ - + -
+ ]
905 : :
906 [ + + ]: 29 : if (from_file) then
907 [ + - + - : 66 : screen%qlwl = Hscr%qlwl
- + - - -
- - - + +
+ + ]
908 [ + - - + : 33 : screen%omega = Hscr%omega
- - - - +
+ ]
909 : :
910 : : ! G-vectors.
911 : 11 : screen%npw = Hscr%npwe
912 [ + - ]: 11 : if (npw_asked > 0) then
913 [ - + ]: 11 : if (npw_asked > Hscr%npwe) then
914 : : write(msg,'(a,i0,2a,i0)') &
915 : 0 : 'The number of G-vectors saved on file is less than the value required: ',npw_asked,ch10,&
916 : 0 : 'Calculation will proceed with the Max available npw: ',Hscr%npwe
917 : 0 : ABI_WARNING(msg)
918 : : else
919 : 11 : screen%npw = npw_asked ! Redefine the no. of G"s for W.
920 : : write(msg,'(a,i0,2a,i0)')&
921 : 11 : 'The number of G-vectors saved on file is larger than the value required: ',npw_asked,ch10,&
922 : 22 : 'Calculation will proceed with npw: ',screen%npw
923 : 11 : ABI_COMMENT(msg)
924 : : end if
925 : : end if
926 : :
927 : : ! Here consistency check on G-vectors and q-points.
928 [ + + + + : 2303 : if (ANY(Hscr%gvec(:,1:screen%npw) /= Gsph%gvec(:,1:screen%npw))) then
+ - - + ]
929 : : !write(std_out) W%gvec, Gsph%gvec
930 : 0 : ABI_ERROR("Hscr%gvec /= Gsph%gvec(1:W%npw)")
931 : : end if
932 [ - + ]: 11 : ABI_CHECK(Hscr%nqibz == Qmesh%nibz, "Mismatch in the number of q-points in the IBZ")
933 : 11 : ierr = 0
934 [ + + ]: 78 : do iq_ibz=1,Hscr%nqibz
935 [ + + + - : 279 : if (ANY(ABS(Qmesh%ibz(:,iq_ibz) - Hscr%qibz(:,iq_ibz)) > tol6) ) then
- + ]
936 : 0 : ierr = ierr + 1
937 : 0 : write(std_out,'(i0,2(3f7.3,1x))')iq_ibz, Qmesh%ibz(:,iq_ibz), Hscr%qibz(:,iq_ibz)
938 : : end if
939 : : end do
940 [ - + ]: 11 : ABI_CHECK(ierr == 0, "Wrong ordering in q-point list, Aborting now")
941 : : end if
942 : :
943 [ - + + - : 87 : ABI_MALLOC(screen%gvec, (3, screen%npw))
+ - + - -
+ - + -
+ ]
944 [ + - - + : 6246 : screen%gvec = Gsph%gvec(:,1:screen%npw)
- - + + +
+ ]
945 : :
946 : : ! Frequency mesh.
947 : 29 : screen%nomega_r = 1; screen%nomega_i = 0
948 [ - + ]: 29 : if (screen%nomega == 2 ) then
949 : 0 : screen%nomega_r = 1; screen%nomega_i = 1
950 : : else
951 : : ! Real frequencies are packed in the first locations.
952 : : screen%nomega_r = 1
953 [ + + ]: 58 : do iw=1,screen%nomega
954 [ - + ]: 58 : if (DBLE(screen%omega(iw))>0.001*Ha_eV) screen%nomega_r=iw
955 : : end do
956 : 29 : screen%nomega_i = screen%nomega - screen%nomega_r
957 : : end if
958 : :
959 : : ! ------------------------------ Initialization completed --------------------------------
960 : : !
961 : : ! Just to keep the code below more readable.
962 : 29 : npw = screen%npw
963 : 29 : nqibz = screen%nqibz
964 : 29 : nomega = screen%nomega
965 : 29 : nI = screen%ni
966 : 29 : nJ = screen%nj
967 [ - + + - : 263 : ABI_MALLOC(screen%Fgg, (nqibz))
+ - + - -
+ - + +
+ ]
968 : :
969 [ + + ]: 29 : if (from_file) then
970 : :
971 : : ! Read ab-initio em1 from file.
972 : 11 : select case (mat_type_read)
973 : : case (MAT_INV_EPSILON)
974 : 11 : call wrtout(std_out, strcat(" Em1 will be initialized from SCR file: ", screen%fname))
975 : :
976 : : case (MAT_CHI0)
977 : : ! Should Write new SCR file.
978 : 0 : ABI_ERROR("Not coded yet")
979 : 0 : sus_fname = screen%fname; scr_fname="TESTING_SUS2SCR"
980 : :
981 : 0 : screen%fname = scr_fname ! Change the name of the file associated to W.
982 : :
983 : : case default
984 : 0 : write(msg,'(a,i0)')" Unsupported conversion from mat_type ",mat_type_read
985 [ + - - ]: 11 : ABI_ERROR(msg)
986 : : end select
987 : :
988 : : ! Begin reading.
989 [ + + ]: 78 : do iq_ibz=1,nqibz
990 [ - + ]: 67 : if (.not. screen%keep_qibz(iq_ibz)) then
991 : : !call wrtout(std_out, strcat("Skipping iq_ibz: ",itoa(iq_ibz)))
992 : : CYCLE
993 : : end if
994 : :
995 : 67 : nqlwl = 0; is_qeq0 = (normv(screen%qibz(:,iq_ibz),Cryst%gmet,'G') < GW_TOLQ0)
996 [ + + ]: 67 : if (is_qeq0) nqlwl=screen%nqlwl
997 : :
998 : : ! Allocate F_{GG'}(w)
999 : 67 : call screen%Fgg(iq_ibz)%init(npw, nomega, nqlwl)
1000 : :
1001 : : ! Read data from file (use MPI-IO if possible)
1002 [ - + ]: 67 : if (screen%iomode /= IO_MODE_ETSF .and. xmpi_mpiio == 1) then
1003 : : !call wrtout(std_out, "read_screening with MPI_IO")
1004 : 0 : call read_screening(varname, screen%fname, npw, 1, nomega, screen%Fgg(iq_ibz)%mat, IO_MODE_MPI, comm, iqiA=iq_ibz)
1005 : : else
1006 : 67 : call read_screening(varname, screen%fname, npw, 1, nomega, screen%Fgg(iq_ibz)%mat, screen%iomode, comm, iqiA=iq_ibz)
1007 : : end if
1008 : :
1009 : : ! W contains Em1 and is ready to use.
1010 : 78 : screen%Fgg(iq_ibz)%has_mat = MAT_STORED
1011 : : end do
1012 : :
1013 : : else
1014 : :
1015 : : ! Model dielectric function. Only epsm-1 is supported here.
1016 : 18 : call wrtout(std_out," Calculating model dielectric function... ")
1017 [ - + ]: 18 : ABI_CHECK(screen%nomega == 1, "Cannot use nomega > 1 in model dielectric function")
1018 : :
1019 [ + + ]: 127 : do iq_ibz=1,nqibz
1020 [ + + ]: 109 : if (.not.screen%keep_qibz(iq_ibz)) CYCLE
1021 : :
1022 : : ! The wings are not used here.
1023 : 103 : nqlwl=0; is_qeq0= (normv(screen%qibz(:,iq_ibz),Cryst%gmet,'G')<GW_TOLQ0)
1024 : :
1025 : : ! Calculate the model. Note that mdielf awaits an index in the BZ.
1026 : 103 : found = qmesh%has_bz_item(Qmesh%ibz(:,iq_ibz),iq_bz,g0)
1027 [ + + + - : 412 : if (.not.found .or. any(g0 /= 0)) then
- + ]
1028 : 0 : ABI_ERROR("Problem in retrieving ibz points")
1029 : : end if
1030 : :
1031 : : ! Allocate F_{GG'}(w).
1032 : 103 : call screen%Fgg(iq_ibz)%init(npw, nomega, nqlwl)
1033 : :
1034 : 103 : eps_inf = screen%info%eps_inf
1035 : 103 : mdf_type = screen%info%use_mdf
1036 : : !em1_ggw => screen%Fgg(iq_ibz)%mat
1037 : :
1038 : : ! Construct W TODO check the new implementation.
1039 : : call screen_mdielf(iq_bz,npw,nomega,mdf_type,eps_inf,Cryst,Qmesh,Vcp,Gsph,&
1040 : 103 : nspden,nfftf_tot,ngfftf,ae_rhor,"EM1",screen%Fgg(iq_ibz)%mat,comm)
1041 : :
1042 : 103 : screen%Fgg(iq_ibz)%has_mat = MAT_STORED
1043 : :
1044 [ - + ]: 121 : if (screen%prtvol > 0) then
1045 [ # # ]: 0 : do iw=1,nomega
1046 : 0 : write(msg,'(a,i3,a,i4,a)')' Model symmetrical e^{-1} (q=',iq_ibz,', omega=',iw,', G,G'')'
1047 : 0 : call wrtout(std_out,msg)
1048 [ # # ]: 0 : call print_arr([std_out], screen%Fgg(iq_ibz)%mat(:,:,iw))
1049 : : end do
1050 : : end if
1051 : : end do ! iq_ibz
1052 : : end if
1053 : :
1054 : : ! Init plasmon-pole parameters from em1.
1055 [ - + ]: 29 : if (screen%has_ppmodel > 0) then
1056 : 0 : call wrtout(std_out, " Calling ppm_init ...")
1057 : 0 : ppmodel = screen%info%use_ppm; drude_plsmf = screen%info%drude_plsmf
1058 : 0 : call screen%ppm%init(screen%mqmem, screen%nqibz, screen%npw, ppmodel, drude_plsmf, screen%info%invalid_freq)
1059 : : !call screen%ppm%print(units)
1060 : :
1061 [ # # ]: 0 : do iq_ibz=1,nqibz
1062 [ # # ]: 0 : if (screen%ihave_fgg(iq_ibz, how="Stored")) then
1063 : : !call wrtout(std_out, sjoin(" Calling ppm%new_setup for iq_ibz:", itoa(iq_ibz)))
1064 : : call screen%ppm%new_setup(iq_ibz, Cryst, Qmesh, npw, nomega, screen%omega, &
1065 : 0 : screen%Fgg(iq_ibz)%mat, nfftf_tot, Gsph%gvec, ngfftf, screen%ae_rhor(:,1))
1066 : : end if
1067 : : end do
1068 : :
1069 : : end if
1070 : : !stop
1071 : :
1072 : : ! Deallocate Fgg if the matrices are not needed anymore.
1073 : : if (free_Fgg) then
1074 : : call screen_fgg_qbz_set(screen, 0, 0, "Pointer") ! Avoid dangling pointer.
1075 : : call fgg_free(screen%Fgg, keep_qibz=screen%keep_qibz)
1076 : : end if
1077 : :
1078 [ + + ]: 29 : if (from_file) call Hscr%free()
1079 : :
1080 : 29 : call pstat_proc%print(_PSTAT_ARGS_)
1081 : :
1082 : : DBG_EXIT("COLL")
1083 : :
1084 [ - + - + : 29 : end subroutine screen_init
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - - -
- - - - -
- - ]
1085 : : !!***
1086 : :
1087 : : !----------------------------------------------------------------------
1088 : :
1089 : : !!****f* m_screen/screen_rotate_iqbz
1090 : : !! NAME
1091 : : !! screen_rotate_iqbz
1092 : : !!
1093 : : !! FUNCTION
1094 : : !! Modify the status of the object so that the symmetrized component F(q_bz)_GG' is calculated
1095 : : !! (if needed) and is made available in the internal buffer. This routine must be called before
1096 : : !! performing any operation that involves the symmetrized component of the two-point function.
1097 : : !!
1098 : : !! INPUTS
1099 : : !! iq_bz=Index of the q-point in the BZ where F(q_bz)_GG' is wanted.
1100 : : !! Cryst=Crystal structure.
1101 : : !! Gsph=The G-sphere
1102 : : !! Qmesh=Structure defining the q-mesh used for sampling the BZ.
1103 : : !!
1104 : : !! SIDE EFFECTS
1105 : : !! screen%ppm
1106 : : !! screen%Fgg_qbz
1107 : : !!
1108 : : !! NOTES
1109 : : !! In the present implementation, we are not considering a possible umklapp vector G0 in the
1110 : : !! expression Sq = q+G0. Treating this case would require some changes in the G-sphere
1111 : : !! since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required
1112 : : !! to reconstruct the BZ.
1113 : : !!
1114 : : !! SOURCE
1115 : :
1116 : 22740 : subroutine screen_rotate_iqbz(screen, iq_bz, Cryst, Gsph, Qmesh, Vcp)
1117 : :
1118 : : !Arguments ------------------------------------
1119 : : !scalars
1120 : : class(screen_t),intent(inout) :: screen
1121 : : integer,intent(in) :: iq_bz
1122 : : type(crystal_t),intent(in) :: Cryst
1123 : : type(gsphere_t),intent(in) :: Gsph
1124 : : type(kmesh_t),intent(in) :: Qmesh
1125 : : type(vcoul_t),intent(in) :: Vcp
1126 : :
1127 : : !Local variables-------------------------------
1128 : : !scalars
1129 : : integer,parameter :: nqlwl0=0
1130 : : integer :: iq_ibz,isym_q,itim_q,npw,nomega,nqibz,mdf_type
1131 : : real(dp) :: eps_inf
1132 : : logical :: q_isirred
1133 : : !character(len=500) :: msg
1134 : : !arrays
1135 : : real(dp) :: qbz(3)
1136 : : ! *********************************************************************
1137 : :
1138 : : DBG_ENTER("COLL")
1139 : :
1140 : 22740 : npw = screen%npw; nqibz = screen%nqibz; nomega = screen%nomega
1141 : 22740 : call qmesh%get_bz_item(iq_bz, qbz, iq_ibz, isym_q, itim_q, isirred=q_isirred)
1142 : :
1143 : : ! ========================================================
1144 : : ! ==== Branching for in-core or out-of-core solutions ====
1145 : : ! ========================================================
1146 [ + + ]: 22740 : if (screen%ihave_fgg(iq_ibz, how="Stored")) then
1147 : :
1148 [ + + ]: 22212 : if (q_isirred) then
1149 : : ! Symmetrization is not needed. Target the data in memory.
1150 : 2629 : call screen_fgg_qbz_set(screen, iq_bz, nqlwl0, "Pointer")
1151 : 2629 : screen%Fgg_qbz => screen%Fgg(iq_ibz)
1152 : : else
1153 : : ! Allocate space. ! TODO Wings are not symmetrized but oh well
1154 : 19583 : call screen_fgg_qbz_set(screen, iq_bz, nqlwl0, "Allocate") ! Dimensions should not be changed.
1155 : :
1156 : : ! Out-of-place symmetrization.
1157 : : !em1_qibz => screen%Fgg(iq_ibz)%mat; em1_qbz => screen%Fgg_qbz%mat
1158 : 19583 : call em1_symmetrize_op(iq_bz, npw, nomega, Gsph, Qmesh, screen%Fgg(iq_ibz)%mat, screen%Fgg_qbz%mat)
1159 : : end if
1160 : :
1161 [ - + ]: 22212 : if (screen%has_ppmodel > 0) then
1162 : : ! Symmetrize the ppmodel tables: em1_qibz => W%Fgg(iq_ibz)%mat
1163 : : call screen%ppm%rotate_iqbz(iq_bz, Cryst, Qmesh, Gsph, npw, nomega, screen%omega, screen%Fgg(iq_ibz)%mat, &
1164 : 0 : screen%nfftf_tot, screen%ngfftf, screen%ae_rhor(:,1))
1165 : :
1166 : : !call screen%ppm_get_qbz(Gsph, Qmesh, iq_bz, botsq_qbz, otq_qbz, eig_qbz)
1167 : : end if
1168 : :
1169 [ - + ]: 528 : else if (screen%ihave_fgg(iq_ibz, how="Allocated")) then
1170 : 0 : ABI_ERROR("Fgg_iqibz is allocated but not initialized!")
1171 : :
1172 : : else
1173 : : ! Out of core branch
1174 : :
1175 [ + - ]: 528 : if (screen%fgg_qbz_idx /= iq_bz) then
1176 : : ! Must compute em1_qbz here. em1_qbz => W%Fgg_qbz%mat
1177 : : ! Allocate the BZ buffer.
1178 : 528 : call screen_fgg_qbz_set(screen, iq_bz, nqlwl0, "Allocate")
1179 : :
1180 [ + - ]: 528 : if (screen%info%use_mdf /= MDL_NONE) then
1181 : : ! Compute the model-dielectric function at qbz on-the fly and in sequential
1182 : : !call wrtout(std_out,"Will compute MDF on the fly")
1183 : : call screen_mdielf(iq_bz,npw,nomega,screen%info%use_mdf,screen%info%eps_inf,Cryst,Qmesh,Vcp,Gsph,&
1184 : : screen%nspden,screen%nfftf_tot,screen%ngfftf,screen%ae_rhor,"EM1", &
1185 : 528 : screen%Fgg_qbz%mat,xmpi_comm_self)
1186 : :
1187 : : else
1188 : : ! Read W(q_ibz) and symmetrize it (do this only if we don't have the correct q_bz in memory).
1189 : 0 : call wrtout(std_out,sjoin("Out of core with file: ",screen%fname))
1190 : : call read_screening(em1_ncname, screen%fname, npw, 1, nomega, screen%Fgg_qbz%mat, &
1191 : 0 : screen%iomode, xmpi_comm_self, iqiA=iq_ibz)
1192 : :
1193 : : ! In-place symmetrization to get the q-point in the BZ.
1194 [ # # ]: 0 : if (.not. q_isirred) then
1195 : 0 : call em1_symmetrize_ip(iq_bz, npw, nomega, Gsph, Qmesh, screen%Fgg_qbz%mat)
1196 : : end if
1197 : : end if
1198 : :
1199 : 528 : screen%Fgg_qbz%has_mat = MAT_STORED
1200 : : end if
1201 : :
1202 [ - + ]: 528 : ABI_CHECK(screen%Fgg_qbz%has_mat == MAT_STORED, "Wrong has_mat")
1203 : :
1204 : : ! Ppmodel calculations with ppm tables in memory.
1205 : : ! TODO treat the case in which IBZ tables are stored in memory.
1206 [ - + ]: 528 : if (screen%has_ppmodel > 0) then
1207 : 0 : ABI_ERROR("Not implemented error")
1208 : : ! Symmetrize the ppmodel using em1_qibz.
1209 : : call screen%ppm%rotate_iqbz(iq_bz, Cryst, Qmesh, Gsph, npw, nomega, screen%omega, &
1210 : 0 : screen%Fgg_qbz%mat, screen%nfftf_tot, screen%ngfftf, screen%ae_rhor(:,1))
1211 : : end if
1212 : : end if
1213 : :
1214 : : ! Calculate model dielectric function for this q-point in the BZ.
1215 : 22740 : eps_inf = screen%info%eps_inf; mdf_type = screen%info%use_mdf
1216 : :
1217 : : ! Model dielectric function. Only epsm-1 is supported here.
1218 : : !call wrtout(std_out," Calculating model dielectric function... ")
1219 : : !ABI_CHECK(W%nomega==1,"Cannot use nomega>1 in model dielectric function")
1220 : :
1221 : : ! screen%Fgg_qbz%mat
1222 : :
1223 : : !% call screen_mdielf(iq_bz,npw,nomega,mdf_type,eps_inf,Cryst,Qmesh,Vcp,Gsph,&
1224 : : !% & screen%nspden,screen%nfftf_tot,screen%ngfftf,screen%ae_rhor,"EM1",em1_qbz,xmpi_comm_self)
1225 : :
1226 : : ! Store the index of the q-point in the BZ for checking purpose.
1227 : 22740 : screen%fgg_qbz_idx = iq_bz
1228 : :
1229 : : DBG_EXIT("COLL")
1230 : :
1231 : 22740 : end subroutine screen_rotate_iqbz
1232 : : !!***
1233 : :
1234 : : !----------------------------------------------------------------------
1235 : :
1236 : : !!****f* m_screen/screen_w0gemv
1237 : : !! NAME
1238 : : !! screen_w0gemv
1239 : : !!
1240 : : !! FUNCTION
1241 : : !! Perform the matrix multiplication W x vector in reciprocal space.
1242 : : !!
1243 : : !! INPUTS
1244 : : !! in_npw=Number of G vectors in in_ket
1245 : : !! nspinor=Number of spinorial components.
1246 : : !! in_ket(in_npw)= |\phi> in reciprocal space.
1247 : : !! trans= On entry, TRANS specifies the operation to be performed as follows:
1248 : : !! TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
1249 : : !! TRANS = 'T' or 't' y := alpha*A**T*x + beta*y.
1250 : : !! TRANS = 'C' or 'c' y := alpha*A**H*x + beta*y.
1251 : : !!
1252 : : !! OUTPUT
1253 : : !! out_ket(in_npw)= W |\phi\> in reciprocal space.
1254 : : !! ZGEMV performs one of the matrix-vector operations
1255 : : !! *
1256 : : !! * y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y, or
1257 : : !! *
1258 : : !! * y := alpha*A**H*x + beta*y,
1259 : : !! *
1260 : : !! * where alpha and beta are scalars, x and y are vectors and A is an m by n matrix.
1261 : : !!
1262 : : !! SOURCE
1263 : :
1264 : 217920 : subroutine screen_w0gemv(screen, trans, in_npw, nspinor, only_diago, alpha, beta, in_ket, out_ket)
1265 : :
1266 : : !Arguments ------------------------------------
1267 : : !scalars
1268 : : class(screen_t),intent(in) :: screen
1269 : : integer,intent(in) :: in_npw,nspinor
1270 : : complex(gwp),intent(in) :: alpha,beta
1271 : : logical,intent(in) :: only_diago
1272 : : character(len=*),intent(in) :: trans
1273 : : !arrays
1274 : : complex(gwp),intent(in) :: in_ket(in_npw*nspinor)
1275 : : complex(gwp),intent(out) :: out_ket(in_npw*nspinor)
1276 : :
1277 : : !Local variables-------------------------------
1278 : : !scalars
1279 : : integer :: ig,lda
1280 : : !arrays
1281 : 217920 : complex(gwp),contiguous, pointer :: em1_qbz(:,:)
1282 : : ! *************************************************************************
1283 : :
1284 : 217920 : lda = screen%npw; em1_qbz => screen%Fgg_qbz%mat(:,:,1)
1285 : :
1286 : 217920 : if (.not.only_diago) then
1287 : 190880 : call xgemv(trans,in_npw,in_npw,alpha,em1_qbz,lda,in_ket,1,beta,out_ket,1)
1288 : :
1289 : : else
1290 [ - + ]: 27040 : if (beta /= czero_gw) then
1291 [ # # # # ]: 0 : if (firstchar(trans, (/"C"/))) then
1292 [ # # ]: 0 : do ig=1,in_npw
1293 : 0 : out_ket(ig) = alpha * CONJG(em1_qbz(ig,ig)) * in_ket(ig) + beta * out_ket(ig)
1294 : : end do
1295 [ # # # # ]: 0 : else if (firstchar(trans, (/"N","T"/))) then
1296 [ # # ]: 0 : do ig=1,in_npw
1297 : 0 : out_ket(ig) = alpha * em1_qbz(ig,ig) * in_ket(ig) + beta * out_ket(ig)
1298 : : end do
1299 : : else
1300 : 0 : ABI_ERROR(sjoin("Wrong trans:", trans))
1301 : : end if
1302 : :
1303 : : else
1304 : : ! beta == 0
1305 [ + + + - ]: 54080 : if (firstchar(trans, (/"C"/)) ) then
1306 [ + + ]: 757120 : do ig=1,in_npw
1307 : 757120 : out_ket(ig) = alpha * CONJG(em1_qbz(ig,ig)) * in_ket(ig)
1308 : : end do
1309 [ # # # # ]: 0 : else if (firstchar(trans, (/"N","T"/))) then
1310 [ # # ]: 0 : do ig=1,in_npw
1311 : 0 : out_ket(ig) = alpha * em1_qbz(ig,ig) * in_ket(ig)
1312 : : end do
1313 : : else
1314 : 0 : ABI_ERROR(sjoin("Wrong trans:", trans))
1315 : : end if
1316 : : end if
1317 : : end if
1318 : :
1319 [ + + ]: 217920 : end subroutine screen_w0gemv
1320 : : !!***
1321 : :
1322 : : !----------------------------------------------------------------------
1323 : :
1324 : : !!****f* m_screen/screen_calc_sigc
1325 : : !! NAME
1326 : : !! screen_calc_sigc
1327 : : !!
1328 : : !! FUNCTION
1329 : : !!
1330 : : !! INPUTS
1331 : : !! npw_c=Number of G vectors in in_ket
1332 : : !! nspinor=Number of spinorial components.
1333 : : !! in_ket(npw_c)= |\phi> in reciprocal space.
1334 : : !! trans= On entry, TRANS specifies the operation to be performed as follows:
1335 : : !! TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
1336 : : !! TRANS = 'T' or 't' y := alpha*A**T*x + beta*y.
1337 : : !! TRANS = 'T' or 't' y := alpha*A**T*x + beta*y.
1338 : : !!
1339 : : !! OUTPUT
1340 : : !!
1341 : : !! SOURCE
1342 : :
1343 : 0 : subroutine screen_calc_sigc(screen, trans, nomega, omegame0i, theta_mu_minus_e0i, zcut, &
1344 : 0 : nspinor, npw_x, npw_c, rhotwgp, out_ket, sigcme)
1345 : :
1346 : : !Arguments ------------------------------------
1347 : : !scalars
1348 : : class(screen_t),intent(in) :: screen
1349 : : character(len=*),intent(in) :: trans
1350 : : integer,intent(in) :: nomega, nspinor, npw_x, npw_c
1351 : : real(dp),intent(in) :: theta_mu_minus_e0i, zcut
1352 : : !arrays
1353 : : real(dp),intent(in) :: omegame0i(nomega)
1354 : : complex(gwp),intent(in) :: rhotwgp(npw_x*nspinor)
1355 : : complex(gwp),intent(inout) :: out_ket(npw_c*nspinor, nomega)
1356 : : complex(gwp),intent(out) :: sigcme(nomega)
1357 : :
1358 : : !Local variables-------------------------------
1359 : 0 : complex(gwp),allocatable :: botsq_conjg_transp(:,:), otq_transp(:,:)
1360 : : ! *************************************************************************
1361 : :
1362 [ # # # # ]: 0 : out_ket = czero_gw
1363 : :
1364 : 0 : select case (screen%info%wint_method)
1365 : :
1366 : : case (WINT_PPMODEL)
1367 : : !call wrtout(std_out, "screen in WINT_PPMODEL")
1368 [ # # ]: 0 : ABI_CHECK_IGE(screen%has_ppmodel, 0, "has_ppmodel should be > 0")
1369 : :
1370 : 0 : select case (trans)
1371 : : case ("N")
1372 : : associate (botsq => screen%ppm%bigomegatwsq_qbz_vals, &
1373 : : otq => screen%ppm%omegatw_qbz_vals, &
1374 : 0 : eig => screen%ppm%eigpot_qbz_vals)
1375 : :
1376 : : call screen%ppm%calc_sigc(nspinor, npw_c, nomega, rhotwgp, botsq, otq, &
1377 : 0 : omegame0i, zcut, theta_mu_minus_e0i, eig, npw_x, out_ket, sigcme)
1378 : :
1379 : : end associate
1380 : :
1381 : : case ("T")
1382 : : associate (ppm => screen%ppm, &
1383 : : botsq => screen%ppm%bigomegatwsq_qbz_vals, &
1384 : : otq => screen%ppm%omegatw_qbz_vals, &
1385 : 0 : eig => screen%ppm%eigpot_qbz_vals)
1386 : :
1387 [ # # # # : 0 : ABI_MALLOC(botsq_conjg_transp,(PPm%dm2_botsq,npw_c))
# # # # #
# # # #
# ]
1388 [ # # # # : 0 : botsq_conjg_transp=TRANSPOSE(botsq) ! Keep these two lines separated, otherwise gfortran messes up
# # # # #
# # # ]
1389 : : !botsq_conjg_transp=CONJG(botsq_conjg_transp)
1390 [ # # # # : 0 : ABI_MALLOC(otq_transp,(ppm%dm2_otq, ppm%npwc))
# # # # #
# # # # #
# # ]
1391 [ # # # # : 0 : otq_transp=TRANSPOSE(otq)
# # # # #
# # # ]
1392 : :
1393 : : call screen%ppm%calc_sigc(nspinor, npw_c, nomega, rhotwgp, botsq_conjg_transp, otq_transp, &
1394 : 0 : omegame0i, zcut, theta_mu_minus_e0i, eig, npw_x, out_ket, sigcme)
1395 : :
1396 [ # # ]: 0 : ABI_FREE(botsq_conjg_transp)
1397 [ # # # # ]: 0 : ABI_FREE(otq_transp)
1398 : : end associate
1399 : :
1400 : : case default
1401 [ # # # ]: 0 : ABI_ERROR(sjoin("Invalid trans:", trans))
1402 : : end select
1403 : :
1404 : : case default
1405 [ # # ]: 0 : ABI_ERROR(sjoin("Unsupported wint_method:", itoa(screen%info%wint_method)))
1406 : : end select
1407 : :
1408 : 0 : end subroutine screen_calc_sigc
1409 : : !!***
1410 : :
1411 : : !----------------------------------------------------------------------
1412 : :
1413 : :
1414 : : !!****f* m_screen/em1_symmetrize_ip
1415 : : !! NAME
1416 : : !! em1_symmetrize_ip
1417 : : !!
1418 : : !! FUNCTION
1419 : : !! Symmetrizes the two-point function in G-space. Symmetrization is done
1420 : : !! inplace through an auxiliary work array of dimension (npw_c,npw_c)
1421 : : !!
1422 : : !! INPUTS
1423 : : !! nomega=All frequencies from 1 up to nomega are symmetrized.
1424 : : !! npw_c=Number of G vectors in the symmetrized matrix.
1425 : : !! Gsph<gsphere_t>=data related to the G-sphere
1426 : : !! Qmesh<kmesh_t>=Structure defining the q-mesh used for Er.
1427 : : !! iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.
1428 : : !!
1429 : : !! SIDE EFFECTS
1430 : : !! epsm1(npw_c,npw_c,nomega)
1431 : : !! input: filled with the matrix at the q-point that has to be symmetrized.
1432 : : !! output: symmetrised matrix.
1433 : : !!
1434 : : !! NOTES
1435 : : !! In the present implementation we are not considering a possible umklapp vector G0 in the
1436 : : !! expression Sq = q+G0. Treating this case would require some changes in the G-sphere
1437 : : !! since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required
1438 : : !! to reconstruct the BZ.
1439 : : !!
1440 : : !! * Remember the symmetry properties of E
1441 : : !! If q_bz=Sq_ibz+G0:
1442 : : !!
1443 : : !! $ E_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau} E_{G1,G2)}(q)
1444 : : !!
1445 : : !! The invariance under exchange of the real space position E(1,2) = E(2,1) leads to:
1446 : : !! $ E_{-G2,-G1}(-q) = E_{G1,G2)
1447 : : !!
1448 : : !! SOURCE
1449 : :
1450 : 0 : subroutine em1_symmetrize_ip(iq_bz, npw_c, nomega, Gsph, Qmesh, epsm1)
1451 : :
1452 : : !Arguments ------------------------------------
1453 : : !scalars
1454 : : integer,intent(in) :: iq_bz,nomega,npw_c
1455 : : type(gsphere_t),intent(in) :: Gsph
1456 : : type(kmesh_t),intent(in) :: Qmesh
1457 : : !arrays
1458 : : complex(gwp),intent(inout) :: epsm1(npw_c,npw_c,nomega)
1459 : :
1460 : : !Local variables-------------------------------
1461 : : !scalars
1462 : : integer :: iw,g1,g2,isg1,isg2,iq_ibz,itim_q,isym_q,ierr
1463 : : logical :: q_isirred
1464 : : complex(gwp) :: phmsg1t,phmsg2t_star
1465 : : !character(len=500) :: msg
1466 : : !arrays
1467 : : real(dp) :: qbz(3)
1468 : 0 : complex(gwp),allocatable :: work(:,:)
1469 : : ! *********************************************************************
1470 : :
1471 : : ! Get iq_ibz, and symmetries from iq_ibz.
1472 : 0 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred)
1473 : :
1474 [ # # ]: 0 : if (q_isirred) RETURN ! Nothing to do
1475 : :
1476 : : !write(msg,'(a,f8.2,a)')" out of memory in work , requiring ",npw_c**2*gwp*b2Mb," Mb"
1477 [ # # # # : 0 : ABI_MALLOC_OR_DIE(work, (npw_c,npw_c), ierr)
# # # # #
# # # # #
# # # # ]
1478 : :
1479 : : !$OMP PARALLEL DO PRIVATE(isg2,isg1,phmsg1t,phmsg2t_star,work) IF (nomega > 1)
1480 [ # # ]: 0 : do iw=1,nomega
1481 [ # # ]: 0 : do g2=1,npw_c
1482 : 0 : isg2 = Gsph%rottb(g2,itim_q,isym_q)
1483 : 0 : phmsg2t_star = CONJG(Gsph%phmSGt(g2,isym_q))
1484 [ # # ]: 0 : do g1=1,npw_c
1485 : 0 : isg1 = Gsph%rottb(g1,itim_q,isym_q)
1486 : 0 : phmsg1t = Gsph%phmSGt(g1,isym_q)
1487 : 0 : work(isg1,isg2) = epsm1(g1,g2,iw) * phmsg1t * phmsg2t_star
1488 : : end do
1489 : : end do
1490 [ # # # # ]: 0 : epsm1(:,:,iw) = work
1491 : : end do
1492 : :
1493 [ # # ]: 0 : ABI_FREE(work)
1494 : :
1495 : : ! Account for time-reversal
1496 [ # # ]: 0 : if (itim_q==2) then
1497 : : !$OMP PARALLEL DO IF (nomega > 1)
1498 [ # # ]: 0 : do iw=1,nomega
1499 : 0 : call sqmat_itranspose(npw_c, epsm1(:,:,iw))
1500 : : end do
1501 : : end if
1502 : :
1503 : : end subroutine em1_symmetrize_ip
1504 : : !!***
1505 : :
1506 : : !----------------------------------------------------------------------
1507 : :
1508 : : !!****f* m_screen/em1_symmetrize_op
1509 : : !! NAME
1510 : : !! em1_symmetrize_op
1511 : : !!
1512 : : !! FUNCTION
1513 : : !! Symmetrizes the two-point function in G-space. Symmetrization is done outofplace.
1514 : : !!
1515 : : !! INPUTS
1516 : : !! nomega=All frequencies from 1 up to nomega are symmetrized.
1517 : : !! npw_c=Number of G vectors in the symmetrized matrix.
1518 : : !! Gsph<gsphere_t>=data related to the G-sphere
1519 : : !! Qmesh<kmesh_t>=Structure defining the q-mesh used for Er.
1520 : : !! iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.
1521 : : !! in_epsm1(npw_c,npw_c,nomega)
1522 : : !!
1523 : : !! OUTPUT
1524 : : !! out_epsm1(npw_c,npw_c,nomega)
1525 : : !!
1526 : : !! NOTES
1527 : : !! In the present implementation we are not considering a possible umklapp vector G0 in the
1528 : : !! expression Sq = q+G0. Treating this case would require some changes in the G-sphere
1529 : : !! since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required
1530 : : !! to reconstruct the BZ.
1531 : : !!
1532 : : !! * Remember the symmetry properties of E
1533 : : !! If q_bz=Sq_ibz+G0:
1534 : : !!
1535 : : !! $ E_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau} E_{G1,G2)}(q)
1536 : : !!
1537 : : !! The invariance under exchange of the real space position E(1,2) = E(2,1) leads to:
1538 : : !! $ E_{-G2,-G1}(-q) = E_{G1,G2)
1539 : : !!
1540 : : !! SOURCE
1541 : :
1542 : 19583 : subroutine em1_symmetrize_op(iq_bz, npw_c, nomega, Gsph, Qmesh, in_epsm1, out_epsm1)
1543 : :
1544 : : !Arguments ------------------------------------
1545 : : !scalars
1546 : : integer,intent(in) :: iq_bz,nomega,npw_c
1547 : : type(gsphere_t),target,intent(in) :: Gsph
1548 : : type(kmesh_t),intent(in) :: Qmesh
1549 : : !arrays
1550 : : complex(gwp),intent(in) :: in_epsm1(npw_c,npw_c,nomega)
1551 : : complex(gwp),intent(out) :: out_epsm1(npw_c,npw_c,nomega)
1552 : :
1553 : : !Local variables-------------------------------
1554 : : !scalars
1555 : : integer :: iw,g1,g2,isg1,isg2,iq_ibz,itim_q,isym_q
1556 : : logical :: q_isirred
1557 : : complex(gwp) :: phmsg1t,phmsg2t_star
1558 : : !arrays
1559 : : real(dp) :: qbz(3)
1560 : : ! *********************************************************************
1561 : :
1562 : : ! Get iq_ibz, and symmetries from iq_ibz.
1563 : 19583 : call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred)
1564 : :
1565 [ - + ]: 19583 : if (q_isirred) then
1566 [ # # # # : 0 : out_epsm1 = in_epsm1; return
# # ]
1567 : : end if
1568 : :
1569 : : ! grottb is a 1-1 mapping.
1570 : : !$OMP PARALLEL DO PRIVATE(isg1,isg2,phmsg1t,phmsg2t_star) COLLAPSE(2) IF (nomega > 1)
1571 [ + + ]: 39166 : do iw=1,nomega
1572 [ + + ]: 921363 : do g2=1,npw_c
1573 : 882197 : isg2 = Gsph%rottb(g2,itim_q,isym_q)
1574 : 882197 : phmsg2t_star = CONJG(Gsph%phmSGt(g2,isym_q))
1575 [ + + ]: 44918267 : do g1=1,npw_c
1576 : 44016487 : isg1 = Gsph%rottb(g1,itim_q,isym_q)
1577 : 44016487 : phmsg1t = Gsph%phmSGt(g1,isym_q)
1578 : 44898684 : out_epsm1(isg1,isg2,iw) = in_epsm1(g1,g2,iw) * phmsg1t * phmsg2t_star
1579 : : end do
1580 : : end do
1581 : : end do
1582 : :
1583 : : ! Account for time-reversal
1584 [ - + ]: 19583 : if (itim_q==2) then
1585 : : !$OMP PARALLEL DO IF (nomega > 1)
1586 [ # # ]: 0 : do iw=1,nomega
1587 : 0 : call sqmat_itranspose(npw_c, out_epsm1(:,:,iw))
1588 : : end do
1589 : : end if
1590 : :
1591 : : end subroutine em1_symmetrize_op
1592 : : !!***
1593 : :
1594 : : !----------------------------------------------------------------------
1595 : :
1596 [ - + - - : 87 : end module m_screen
+ + - + +
- - + - +
+ - - + -
+ - + - +
- + - + -
+ - + - +
- + - - -
- - + - +
- - - - -
+ - + - -
- - - + ]
1597 : : !!***
|
