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 : !!***
|
