TABLE OF CONTENTS


ABINIT/m_screen [ Modules ]

[ Top ] [ Modules ]

NAME

  m_screen

FUNCTION

   Screening object used in the BSE code.

COPYRIGHT

  Copyright (C) 2014-2022 ABINIT group (MG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 MODULE m_screen
23 
24  use defs_basis
25  use m_xmpi
26  use m_hide_blas
27  use m_errors
28  use m_splines
29  use m_abicore
30  use m_kxc
31  use m_screening
32  use m_nctk
33  use m_sort
34 
35  use m_gwdefs,         only : GW_TOLQ0, czero_gw
36  use m_fstrings,       only : firstchar, endswith, strcat, itoa, sjoin
37  use m_numeric_tools,  only : print_arr
38  use m_geometry,       only : normv
39  use m_crystal,        only : crystal_t
40  use m_bz_mesh,        only : kmesh_t
41  use m_gsphere,        only : gsphere_t
42  use m_vcoul,          only : vcoul_t
43  use m_io_screening,   only : read_screening, hscr_t, ncname_from_id, em1_ncname
44  use m_ppmodel,        only : ppmodel_t, ppm_init, ppm_free, ppm_nullify, PPM_NONE, new_setup_ppmodel, ppm_symmetrizer
45 
46  implicit none
47 
48  private
49 
50  ! Flags defining the content of the %mat buffer in the fgg_t type.
51  integer,public,parameter :: MAT_NOTYPE         = 0
52  integer,public,parameter :: MAT_CHI0           = 1
53  integer,public,parameter :: MAT_CHI            = 2
54  integer,public,parameter :: MAT_EPSILON        = 3
55  integer,public,parameter :: MAT_INV_EPSILON    = 4
56  integer,public,parameter :: MAT_INV_EPSILON_M1 = 5
57  integer,public,parameter :: MAT_W              = 6
58  integer,public,parameter :: MAT_W_M1           = 7
59  !
60  ! Family vertex.
61  integer,public,parameter :: VTX_FAMILY_NONE  = 0   ! No vertex correction.
62  integer,public,parameter :: VTX_FAMILY_TDDFT = 1   ! TDDFT-based vertex.
63  integer,public,parameter :: VTX_FAMILY_ADA   = 2   ! ADA vertex.
64  !
65  ! Test charge or test particle.
66  integer,public,parameter :: VTX_TEST_CHARGE   = 0
67  integer,public,parameter :: VTX_TEST_PARTICLE = 1
68  !
69  ! Named constants for the frequency mesh.
70  integer,private,parameter :: WMESH_LINEAR    = 1
71  integer,private,parameter :: WMESH_GAUSS_LEG = 2
72  integer,private,parameter :: WMESH_TAN_GRID  = 3
73  !
74  ! Method used for the frequency integration.
75  integer,public,parameter ::  WINT_NONE    = 0
76  integer,public,parameter ::  WINT_PPMODEL = 1
77  integer,public,parameter ::  WINT_CONTOUR = 2
78  integer,public,parameter ::  WINT_AC      = 3
79  !
80  ! Parameters used for the model dielectric function.
81  integer,public,parameter :: MDL_NONE      = 0
82  integer,public,parameter :: MDL_BECHSTEDT = 1
83  !
84  ! Flags giving the status of the local buffers defined in fgg_t.
85  integer,private,parameter :: MAT_NODATA    = 0
86  integer,private,parameter :: MAT_ALLOCATED = 1
87  integer,private,parameter :: MAT_STORED    = 2
88  !
89  ! Flags giving the status of the local buffers defined in fgg_t.
90  integer,private,parameter :: FGG_QBZ_ISPOINTER  =1 ! Fgg_qbz is used to store the address in memory.
91  integer,private,parameter :: FGG_QBZ_ISALLOCATED=2 ! Fgg_qbz is used as an allocable array.

m_screen/em1_symmetrize_ip [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

  em1_symmetrize_ip

FUNCTION

  Symmetrizes the two-point function in G-space. Symmetrization is done
  inplace thorugh an auxiliary work array of dimension (npwc,npwc)

INPUTS

  nomega=All frequencies from 1 up to nomega are symmetrized.
  npwc=Number of G vectors in the symmetrized matrix.
  Gsph<gsphere_t>=data related to the G-sphere
  Qmesh<kmesh_t>=Structure defining the q-mesh used for Er.
  iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.

SIDE EFFECTS

  epsm1(npwc,npwc,nomega)
   input:  filled with the matrix at the q-point that has to be symmetrized.
   output: symmetrised matrix.

NOTES

  In the present implementation we are not considering a possible umklapp vector G0 in the
  expression Sq = q+G0. Treating this case would require some changes in the G-sphere
  since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required
  to reconstruct the BZ.

  * Remember the symmetry properties of E
    If q_bz=Sq_ibz+G0:

    $ E_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau} E_{G1,G2)}(q)

    The invariance under exchange of the real space position E(1,2) = E(2,1) leads to:
    $ E_{-G2,-G1}(-q) = E_{G1,G2)

SOURCE

1378 subroutine em1_symmetrize_ip(iq_bz,npwc,nomega,Gsph,Qmesh,epsm1)
1379 
1380 !Arguments ------------------------------------
1381 !scalars
1382  integer,intent(in) :: iq_bz,nomega,npwc
1383  type(gsphere_t),intent(in) :: Gsph
1384  type(kmesh_t),intent(in) :: Qmesh
1385 !arrays
1386  complex(gwpc),intent(inout) :: epsm1(npwc,npwc,nomega)
1387 
1388 !Local variables-------------------------------
1389 !scalars
1390  integer :: iw,g1,g2,isg1,isg2,iq_ibz,itim_q,isym_q,ierr
1391  logical :: q_isirred
1392  complex(gwpc) :: phmsg1t,phmsg2t_star
1393  !character(len=500) :: msg
1394 !arrays
1395  real(dp) :: qbz(3)
1396  complex(gwpc),allocatable :: work(:,:)
1397 
1398 ! *********************************************************************
1399 
1400  ! * Get iq_ibz, and symmetries from iq_ibz.
1401  call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred)
1402 
1403  if (q_isirred) RETURN ! Nothing to do
1404 
1405  !write(msg,'(a,f8.2,a)')" out of memory in work , requiring ",npwc**2*gwpc*b2Mb," Mb"
1406  ABI_MALLOC_OR_DIE(work,(npwc,npwc), ierr)
1407 
1408 !$OMP PARALLEL DO PRIVATE(isg2,isg1,phmsg1t,phmsg2t_star,work) IF (nomega > 1)
1409  do iw=1,nomega
1410    do g2=1,npwc
1411      isg2 = Gsph%rottb(g2,itim_q,isym_q)
1412      phmsg2t_star = CONJG(Gsph%phmSGt(g2,isym_q))
1413      do g1=1,npwc
1414        isg1 = Gsph%rottb(g1,itim_q,isym_q)
1415        phmsg1t = Gsph%phmSGt(g1,isym_q)
1416        work(isg1,isg2) = epsm1(g1,g2,iw) * phmsg1t  * phmsg2t_star
1417      end do
1418    end do
1419    epsm1(:,:,iw) = work
1420  end do
1421 
1422  ABI_FREE(work)
1423  !
1424  ! Account for time-reversal ----------------------
1425  if (itim_q==2) then
1426 !$OMP PARALLEL DO IF (nomega > 1)
1427    do iw=1,nomega
1428      call sqmat_itranspose(npwc,epsm1(:,:,iw))
1429    end do
1430  end if
1431 
1432 end subroutine em1_symmetrize_ip

m_screen/em1_symmetrize_op [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

  em1_symmetrize_op

FUNCTION

  Symmetrizes the two-point function in G-space. Symmetrization is done outofplace.

INPUTS

  nomega=All frequencies from 1 up to nomega are symmetrized.
  npwc=Number of G vectors in the symmetrized matrix.
  Gsph<gsphere_t>=data related to the G-sphere
  Qmesh<kmesh_t>=Structure defining the q-mesh used for Er.
  iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.
  in_epsm1(npwc,npwc,nomega)

OUTPUT

  out_epsm1(npwc,npwc,nomega)

NOTES

  In the present implementation we are not considering a possible umklapp vector G0 in the
  expression Sq = q+G0. Treating this case would require some changes in the G-sphere
  since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required
  to reconstruct the BZ.

  * Remember the symmetry properties of E
    If q_bz=Sq_ibz+G0:

    $ E_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau} E_{G1,G2)}(q)

    The invariance under exchange of the real space position E(1,2) = E(2,1) leads to:
    $ E_{-G2,-G1}(-q) = E_{G1,G2)

SOURCE

1471 subroutine em1_symmetrize_op(iq_bz,npwc,nomega,Gsph,Qmesh,in_epsm1,out_epsm1)
1472 
1473 !Arguments ------------------------------------
1474 !scalars
1475  integer,intent(in) :: iq_bz,nomega,npwc
1476  type(gsphere_t),target,intent(in) :: Gsph
1477  type(kmesh_t),intent(in) :: Qmesh
1478 !arrays
1479  complex(gwpc),intent(in) :: in_epsm1(npwc,npwc,nomega)
1480  complex(gwpc),intent(out) :: out_epsm1(npwc,npwc,nomega)
1481 
1482 !Local variables-------------------------------
1483 !scalars
1484  integer :: iw,g1,g2,isg1,isg2,iq_ibz,itim_q,isym_q
1485  logical :: q_isirred
1486  complex(gwpc) :: phmsg1t,phmsg2t_star
1487 !arrays
1488  real(dp) :: qbz(3)
1489 
1490 ! *********************************************************************
1491 
1492  ! Get iq_ibz, and symmetries from iq_ibz.
1493  call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred)
1494 
1495  if (q_isirred) then
1496    out_epsm1 = in_epsm1; return
1497  end if
1498 
1499 ! grottb is a 1-1 mapping.
1500 !$OMP PARALLEL DO PRIVATE(isg1,isg2,phmsg1t,phmsg2t_star) COLLAPSE(2) IF (nomega > 1)
1501  do iw=1,nomega
1502    do g2=1,npwc
1503      isg2=Gsph%rottb(g2,itim_q,isym_q)
1504      phmsg2t_star = CONJG(Gsph%phmSGt(g2,isym_q))
1505      do g1=1,npwc
1506        isg1=Gsph%rottb(g1,itim_q,isym_q)
1507        phmsg1t = Gsph%phmSGt(g1,isym_q)
1508        out_epsm1(isg1,isg2,iw) = in_epsm1(g1,g2,iw) * phmsg1t * phmsg2t_star
1509      end do
1510    end do
1511  end do
1512  !
1513  ! Account for time-reversal ----------------------
1514  if (itim_q==2) then
1515 !$OMP PARALLEL DO IF (nomega > 1)
1516    do iw=1,nomega
1517      call sqmat_itranspose(npwc,out_epsm1(:,:,iw))
1518    end do
1519  end if
1520 
1521 end subroutine em1_symmetrize_op

m_screen/fgg_free_0D [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

 fgg_free_0D

FUNCTION

  Free memory.

SOURCE

401 subroutine fgg_free_0D(Fgg)
402 
403 !Arguments ------------------------------------
404 !scalars
405  type(fgg_t),intent(inout) :: Fgg
406 ! *************************************************************************
407 
408  !@fgg_t
409  !complex
410  if (allocated(Fgg%mat))  then
411    ABI_FREE(Fgg%mat)
412  end if
413  Fgg%has_mat = MAT_NODATA
414 
415 end subroutine fgg_free_0D

m_screen/fgg_free_1D [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

 fgg_free_1D

FUNCTION

 Deallocate all the pointers in the structure that result to be associated.

 INPUT
  [keep_q(:)]=Optional logical mask used to select the q-points that are deallocated.

SOURCE

432 subroutine fgg_free_1D(Fgg,keep_q)
433 
434 !Arguments ------------------------------------
435 !scalars
436  type(fgg_t),intent(inout) :: Fgg(:)
437  logical,optional,intent(in) :: keep_q(:)
438 
439 !Local variables ------------------------------
440 !scalars
441  integer :: iq_ibz
442  logical :: keep_it
443 ! *************************************************************************
444 
445  do iq_ibz=LBOUND(Fgg,DIM=1),UBOUND(Fgg,DIM=1)
446    keep_it = .FALSE.; if (PRESENT(keep_q)) keep_it = keep_q(iq_ibz)
447    if (.not. keep_it) call fgg_free_0D(Fgg(iq_ibz))
448  end do
449 
450 end subroutine fgg_free_1D

m_screen/fgg_init [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

 fgg_init

FUNCTION

 Initialize the structure allocating the memory and initializing the internal variables.

 INPUT
  npw
  nqlwl

SIDE EFFECTS

  Fgg, See function description

SOURCE

471 subroutine fgg_init(Fgg,npw,nomega,nqlwl)
472 
473 !Arguments ------------------------------------
474 !scalars
475  integer,intent(in) :: npw,nqlwl,nomega
476  type(fgg_t),intent(inout) :: Fgg
477 
478 !Local variables ------------------------------
479 !scalars
480  integer :: ierr
481 
482 ! *************************************************************************
483 
484  !@fgg_t
485  Fgg%nomega= nomega; Fgg%npw   = npw
486  Fgg%nqlwl = nqlwl
487 
488  if (npw>0.and.nomega>0) then
489    ABI_MALLOC_OR_DIE(Fgg%mat,(npw,npw,nomega), ierr)
490    Fgg%has_mat = MAT_ALLOCATED
491  end if
492 
493 end subroutine fgg_init

m_screen/fgg_t [ Types ]

[ Top ] [ m_screen ] [ Types ]

NAME

 fgg_t

FUNCTION

  Object used to store F(G,G')(w) for a given q-point.

SOURCE

169  type,public :: fgg_t
170 
171   integer :: nomega
172   ! Number of frequencies.
173 
174   integer :: npw
175   ! Number of G vectors.
176 
177   integer :: nqlwl
178   ! Number of points for the treatment of the long wave-length limit.
179 
180   integer :: has_mat = MAT_NODATA
181   ! Flag giving the status of mat.
182 
183   !arrays
184   complex(gwpc),allocatable :: mat(:,:,:)
185   ! mat(npw,npw,nomega)
186   ! The component of the two-point function $F_{G,G',w}$ for a given q.
187 
188   !complex(dpc),allocatable :: head(:,:,:)
189   ! head(3,3,nomega)
190 
191   !complex(dpc),allocatable :: lwing(:,:,:)
192   ! lwing(3,npwe,nomega)
193   ! Lower wings
194 
195   !complex(dpc),allocatable :: uwing(:,:,:)
196   ! uwing(3,npwe,nomega)
197   ! Upper wings.
198 
199  end type fgg_t
200 
201  public :: fgg_init   ! Creation method.
202  public :: fgg_free   ! Free memory.

m_screen/screen_fgg_qbz_set [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

  screen_fgg_qbz_set

FUNCTION

  Helper function used to perform the setup W%Fgg_qbz setting also the internal
  flag that defines its status.

INPUTS

  W<screen_t>=Object describing the screened interaction.
  iq_bz=Index of the q-point in the BZ.
  nqlwl=Number of wings wanted.
  how= "Pointer" is a true pointer is wanted.
       "Allocated" if memory has to be allocated.

NOTES

  iq_bz and nqlwl are not used if how="Pointer".

SOURCE

518 subroutine screen_fgg_qbz_set(W,iq_bz,nqlwl,how)
519 
520 !Arguments ------------------------------------
521 !scalars
522  integer,intent(in) :: iq_bz,nqlwl
523  character(len=*),intent(in) :: how
524  type(screen_t),intent(inout) :: W
525 
526 !Local variables ------------------------------
527 !scalars
528  !character(len=500) :: msg
529 
530 !************************************************************************
531 
532  W%fgg_qbz_idx = iq_bz ! Save the index of the q-point in the BZ.
533 
534  if (firstchar(how,(/"P"/)) ) then
535    ! We want a pointer.
536    select case (W%fgg_qbz_stat)
537    case (FGG_QBZ_ISALLOCATED)
538      call fgg_free_0D(W%Fgg_qbz)
539      ABI_FREE(W%Fgg_qbz)
540      nullify(W%Fgg_qbz)
541      W%fgg_qbz_stat = FGG_QBZ_ISPOINTER
542 
543    case (FGG_QBZ_ISPOINTER)
544      ! Set it to null().
545      nullify(W%Fgg_qbz)
546 
547    case default
548      ABI_ERROR(sjoin("Wrong status:", itoa(W%fgg_qbz_stat)))
549    end select
550 
551  else if (firstchar(how,(/"A"/)) ) then
552    ! We want an allocatable array.
553 
554    select case (W%fgg_qbz_stat)
555    case (FGG_QBZ_ISPOINTER)
556      ! Allocate memory
557      nullify(W%Fgg_qbz)
558      ABI_MALLOC(W%Fgg_qbz,)
559 
560      call fgg_init(W%Fgg_qbz,W%npw,W%nomega,nqlwl)
561      W%fgg_qbz_stat = FGG_QBZ_ISALLOCATED
562 
563    case (FGG_QBZ_ISALLOCATED)
564      W%Fgg_qbz%has_mat = MAT_ALLOCATED  ! STORED --> ALLOCATED
565 
566    case default
567      ABI_ERROR(sjoin("Wrong status:", itoa(W%fgg_qbz_stat)))
568    end select
569 
570  else
571    ABI_BUG(sjoin("Wrong how:", how))
572  end if
573 
574 end subroutine screen_fgg_qbz_set

m_screen/screen_free [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

 screen_free

FUNCTION

 Free the memory allocate in the datatype.

SIDE EFFECTS

 W<screen_t>=The data structure to be nullified.

OUTPUT

SOURCE

686 subroutine screen_free(W)
687 
688 !Arguments ------------------------------------
689 !scalars
690  type(screen_t),intent(inout) :: W
691 
692 ! *************************************************************************
693 
694  !@screen_t
695  !
696  ! integer
697  ABI_SFREE(W%gvec)
698  !
699  !real
700  ABI_SFREE(W%ae_rhor)
701  ABI_SFREE(W%qibz)
702  ABI_SFREE(W%qlwl)
703  !
704  !complex
705  ABI_SFREE(W%omega)
706 
707  ! logical
708  ABI_SFREE(W%keep_q)
709  !
710  ! types ------------------------------------------
711  !
712  ! Here be careful with dangling pointers.
713  ! First Fgg_qbz that might point to one of the %Fgg then %Fgg.
714  select case (W%fgg_qbz_stat)
715  case (FGG_QBZ_ISALLOCATED)
716    call fgg_free_0D(W%Fgg_qbz)
717    ABI_FREE(W%Fgg_qbz)
718    nullify(W%Fgg_qbz)
719    W%fgg_qbz_stat=FGG_QBZ_ISPOINTER
720 
721  case (FGG_QBZ_ISPOINTER)
722    nullify(W%Fgg_qbz)
723    W%fgg_qbz_stat=FGG_QBZ_ISPOINTER
724 
725  case default
726    continue
727  end select
728  !
729  ! Free the Fgg matrices.
730  if (associated(W%Fgg)) then
731    call fgg_free(W%Fgg)
732    ABI_FREE(W%Fgg)
733  end if
734  !
735  ! Free the plasmon pole tables.
736  call ppm_free(W%PPm)
737 
738 end subroutine screen_free

m_screen/screen_ihave_fgg [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

  screen_ihave_fgg

FUNCTION

  Inquire the processor whether it has a particular F_{GG')(q) and with which status.

INPUTS

  W<screen_t>=Object describing the screened interaction.
  iq_ibz=k-point index
  [how]=string defining which status is checked. By default the function returns
     .TRUE. if the wave is either MAT_ALLOCATED or MAT_STORED.
     Possible mutually exclusive values: "Allocated", "Stored".
     Only the first character is checked (no case-sensitive)

NOTES

   A zero index can be used to inquire the status of the full set of q-points.

SOURCE

599 pure function screen_ihave_fgg(W,iq_ibz,how)
600 
601 !Arguments ------------------------------------
602 !scalars
603  integer,intent(in) :: iq_ibz
604  logical :: screen_ihave_fgg
605  character(len=*),optional,intent(in) :: how
606  type(screen_t),intent(in) :: W
607 
608 !Local variables ------------------------------
609 !scalars
610  integer :: ii
611 !arrays
612  integer :: check(2)
613 
614 !************************************************************************
615 
616  check = (/MAT_ALLOCATED, MAT_STORED/)
617  if (PRESENT(how)) then
618    if (firstchar(how,(/"A","a"/))) check = (/MAT_ALLOCATED, MAT_ALLOCATED/)
619    if (firstchar(how,(/"S","s"/))) check = (/MAT_STORED, MAT_STORED/)
620  end if
621 
622  if (iq_ibz>0) then
623    screen_ihave_fgg = (W%Fgg(iq_ibz)%has_mat==check(1) .or.&
624                        W%Fgg(iq_ibz)%has_mat==check(2) )
625  else
626    ! check the status of the full set of q-tables.
627    screen_ihave_fgg=.TRUE.
628    do ii=1,W%nqibz
629      screen_ihave_fgg = screen_ihave_fgg .and. &
630                      (W%Fgg(ii)%has_mat==check(1) .or.&
631                       W%Fgg(ii)%has_mat==check(2) )
632    end do
633  end if
634 
635 end function screen_ihave_fgg

m_screen/screen_info_print [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

  screen_info_print

FUNCTION

  Printout of screen_info_t.

INPUTS

  W_info<screen_info_t>=The structure.
  [unit]=Unit number for output
  [prtvol]=Verbosity level
  [mode_paral]=Either "COLL" or "PERS"
  [header]=String to be printed as header for additional info.

OUTPUT

  Only printing

SOURCE

335 subroutine screen_info_print(W_info,header,unit,mode_paral,prtvol)
336 
337 !Arguments ------------------------------------
338 !scalars
339  integer,optional,intent(in) :: unit,prtvol
340  character(len=4),optional,intent(in) :: mode_paral
341  character(len=*),optional,intent(in) :: header
342  type(screen_info_t),intent(in) :: W_info
343 
344 !Local variables-------------------------------
345  integer :: my_unt,my_prtvol
346  character(len=4) :: my_mode
347  character(len=500) :: msg
348 ! *********************************************************************
349 
350  !@screen_info_t
351  my_unt   =std_out; if (PRESENT(unit      )) my_unt   =unit
352  my_prtvol=0      ; if (PRESENT(prtvol    )) my_prtvol=prtvol
353  my_mode  ='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral
354 
355  msg=' ==== Info on the screen_info_t% object ==== '
356  if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
357  call wrtout(my_unt,msg,my_mode)
358 
359 !integer
360  write(msg,'(a,i3)')" mat_type    ",W_info%mat_type
361  call wrtout(my_unt,msg,my_mode)
362  write(msg,'(a,i3)')" vtx_family  ",W_info%vtx_family
363  call wrtout(my_unt,msg,my_mode)
364  write(msg,'(a,i3)')" invalid_freq",W_info%invalid_freq
365  call wrtout(my_unt,msg,my_mode)
366  write(msg,'(a,i3)')" ixc         ",W_info%ixc
367  call wrtout(my_unt,msg,my_mode)
368  write(msg,'(a,i3)')" use_ada     ",W_info%use_ada
369  call wrtout(my_unt,msg,my_mode)
370  write(msg,'(a,i3)')" use_mdf     ",W_info%use_mdf
371  call wrtout(my_unt,msg,my_mode)
372  write(msg,'(a,i3)')" use_ppm     ",W_info%use_ppm
373  call wrtout(my_unt,msg,my_mode)
374  write(msg,'(a,i3)')" vtx_test    ",W_info%vtx_test
375  call wrtout(my_unt,msg,my_mode)
376  write(msg,'(a,i3)')" wint_method ",W_info%wint_method
377  call wrtout(my_unt,msg,my_mode)
378 
379 !real
380  write(msg,'(a,f8.3)')" ada_kappa   ",W_info%ada_kappa
381  call wrtout(my_unt,msg,my_mode)
382  write(msg,'(a,f8.3)')" eps_inf     ",W_info%eps_inf
383  call wrtout(my_unt,msg,my_mode)
384  write(msg,'(a,f8.3)')" drude_plsmf ",W_info%drude_plsmf
385  call wrtout(my_unt,msg,my_mode)
386 
387 end subroutine screen_info_print

m_screen/screen_info_t [ Types ]

[ Top ] [ m_screen ] [ Types ]

NAME

 screen_info_t

FUNCTION

  Container storing the parameters used to initialize a screen_t datatype or to
  calculate a new SCR file from the SUSC file containing the independent-particle
  polarizability.

NOTES

  The list of parameters passed to screening_init is copied in W%Info.
  At present there is no need to provide a copy method since the structure does
  not contain pointers but such a method must be defined and used if the
  dynamic entities are added to the datatype.

SOURCE

113 type,public :: screen_info_t
114 
115   integer :: mat_type = MAT_NOTYPE
116   ! Matrix identifier. See MAT_* flags.
117 
118   integer :: vtx_family = VTX_FAMILY_NONE
119   ! Vertex correction family.
120 
121   integer :: invalid_freq=0
122   ! Sets the procedure to follow when a PPM frequency is invalid (negative or imaginary), see input variable gw_invalid_freq
123 
124   integer :: ixc=0
125   ! XC functional used for the TDDFT-based vertex.
126 
127   integer :: use_ada=0
128   ! >0 if ADA vertex is used.
129 
130   integer :: use_mdf = MDL_NONE
131   ! >0 if model dielectric function is used.
132 
133   integer :: use_ppm = PPM_NONE
134   ! >0 if ppmodel is used.
135 
136   integer :: vtx_test = VTX_TEST_CHARGE
137   ! test charge or test particle.
138 
139   integer :: wint_method = WINT_NONE
140   ! Defines the frequency integration technique. See WIN_ flags.
141   ! NOTE that this flag can be changed at run time. For example
142   ! one can switch from the CD to the PPm if the ppmodel parameters are in memory
143 
144   real(dp) :: ada_kappa = 2.1_dp
145   ! Inverse smearing length used for ADA.
146 
147   real(dp) :: eps_inf = 12.0_dp
148   ! Dielectric constant used for the model dielectric function.
149 
150   real(dp) :: drude_plsmf = zero
151   ! Drude plasma frequency used for PPmodel 1.
152 
153 end type screen_info_t
154 
155  public :: screen_info_print

m_screen/screen_init [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

  screen_init

FUNCTION

  Initialize basic dimensions and the important (small) arrays in an Epsilonm1_results data type
  starting from a file containing either epsilon^{-1} (_SCR) or chi0 (_SUSC).

INPUTS

  W_Info<screen_info_t>=The list of parameters used to construct the screen function.
  Cryst<crystal_t>=Info on the unit cell.
  Qmesh<kmesh_t>=Info on the Q-mesh.
  Gsph<gsphere_t>=Info on the plane-wave basis set used for the two-point function.
  Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction.
  ifname=The name of the external file used to read the matrix.
  id_required=Identifier used to specify the type of two-point function that is wanted.
  iomode=Option defining the file format of the external file.
  mqmem=0 for out-of-core solution, /=0 if entire matrix has to be stored in memory.
  npw_asked=Number of G-vector to be used in the calculation, if <=0 use Max allowed number.
  ngfftf(18)=Info on the (fine) mesh used for the density.
  nfftf_tot=Total number of point in the FFT mesh for ae_rhor
  nsppol=Numer of independent spin polarizations.
  nspden=Number of spin density components in ae_rhor
  ae_rhor(nfftf_tot,nspden)
  prtvol=Verbosity level.
  comm=MPI communicator.

OUTPUT

  W<screen_t>=The structure initialized with basic dimensions and arrays.

SOURCE

 775 subroutine screen_init(W,W_Info,Cryst,Qmesh,Gsph,Vcp,ifname,mqmem,npw_asked,&
 776 &  iomode,ngfftf,nfftf_tot,nsppol,nspden,ae_rhor,prtvol,comm)
 777 
 778 !Arguments ------------------------------------
 779 !scalars
 780  integer,intent(in) :: mqmem,iomode,npw_asked,comm,prtvol,nsppol
 781  integer,intent(in) :: nfftf_tot,nspden
 782  character(len=fnlen),intent(in) :: ifname
 783  type(crystal_t),intent(in) :: Cryst
 784  type(gsphere_t),intent(in) :: Gsph
 785  type(vcoul_t),intent(in) :: Vcp
 786  type(kmesh_t),intent(in) :: Qmesh
 787  type(screen_info_t),intent(in) :: W_Info
 788  type(screen_t),intent(out) :: W
 789 !arrays
 790  integer,intent(in) :: ngfftf(18)
 791  real(dp),intent(in) :: ae_rhor(nfftf_tot,nspden)
 792 
 793 !Local variables-------------------------------
 794 !scalars
 795  integer,parameter :: master=0
 796  integer :: option_test,approx_type,ixc_required,id_required !,nkxc
 797  integer :: fform,my_rank,mat_type_read
 798  integer :: nqibz,nomega,iq_ibz,npw,nqlwl
 799  integer :: nI,nJ,iq_bz,mdf_type,ppmodel,ierr
 800  integer :: iw,qsort,ii
 801  real(dp) :: eps_inf,drude_plsmf
 802  logical :: deallocate_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  type(hscr_t) :: Hscr
 807 !arrays
 808  integer :: g0(3),iperm(Qmesh%nibz)
 809  real(dp) :: wt_list(Qmesh%nibz)
 810  !complex(gwpc),ABI_CONTIGUOUS pointer :: em1_ggw(:,:,:)
 811 
 812 ! *********************************************************************
 813 
 814  DBG_ENTER("COLL")
 815 
 816  ABI_UNUSED(nsppol)
 817 
 818  !@screen_t
 819  my_rank = xmpi_comm_rank(comm)
 820 
 821  call screen_nullify(W)
 822  !
 823  ! Initialize basic parameters ----------------------------------------
 824  !
 825  !@screen_info_t
 826  W%Info = W_info ! Copy basic parameters.
 827  call screen_info_print(W%Info,header="W info",unit=std_out)
 828 
 829  id_required  = W_Info%mat_type
 830  approx_type  = W_Info%vtx_family
 831  option_test  = W_Info%vtx_test
 832  ixc_required = W_Info%ixc
 833  varname = ncname_from_id(id_required)
 834 
 835  !if (ALL(id_required /= (/MAT_W_M1, MAT_W/)) ) then
 836  if (all(id_required /= [MAT_INV_EPSILON])) then
 837    ABI_ERROR(sjoin("id_required:",itoa(id_required),"not available"))
 838  end if
 839 
 840  ! This part must be rationalized.
 841  remove_dgg = (id_required==MAT_W_M1)
 842 
 843  if (W%Info%use_mdf == MDL_NONE) W%fname = ifname
 844  W%nI = 1; W%nJ = 1
 845  !
 846  ! The Q-point sampling is inited from Qmesh.
 847  W%nqibz=Qmesh%nibz
 848  ABI_MALLOC(W%qibz,(3,W%nqibz))
 849  W%qibz= Qmesh%ibz
 850 
 851  W%mqmem=mqmem!;  W%mqmem = 0
 852  if (W%mqmem/=0) W%mqmem=W%nqibz
 853 
 854  ABI_MALLOC(W%keep_q,(W%nqibz))
 855  W%keep_q=.TRUE.; if (W%mqmem == 0) W%keep_q = .False.
 856 
 857  if (W%mqmem/=0 .and. W%mqmem<W%nqibz) then
 858    ! Keep in memory the most representative q-points.
 859    W%keep_q=.FALSE.
 860    wt_list = Qmesh%wt; iperm = (/(ii,ii=1,Qmesh%nibz)/)
 861    call sort_dp(Qmesh%nibz,wt_list,iperm,tol12)
 862    do qsort=Qmesh%nibz,Qmesh%nibz-mqmem+1,1
 863      iq_ibz = iperm(qsort)
 864      W%keep_q(iq_ibz) = .TRUE.
 865    end do
 866  end if
 867  W%fgg_qbz_idx = 0
 868 
 869  W%iomode = iomode
 870  W%prtvol = prtvol
 871 
 872  W%has_ppmodel=0; if (W%Info%use_ppm/=PPM_NONE) W%has_ppmodel=1
 873  !
 874  ! Copy the AE density for the model dielectric function or for the vertex corrections.
 875  W%nspden     = nspden
 876  W%ngfftf     = ngfftf
 877  W%nfftf_tot  = nfftf_tot
 878 
 879  ABI_MALLOC(W%ae_rhor,(nfftf_tot,nspden))
 880  W%ae_rhor = ae_rhor
 881 
 882  deallocate_Fgg=.FALSE.
 883 
 884  W%has_fgg=0; if ( ANY(W%Info%wint_method == (/WINT_CONTOUR, WINT_AC/)) ) W%has_fgg=1
 885 
 886  if (W%has_fgg>0 .and. W%has_ppmodel>0) then
 887    ABI_WARNING("Both PPmodel tables and F_(GG')(q,w) are stored in memory")
 888  end if
 889 
 890  ! G-sphere for W and Sigma_c is initialized from the SCR file.
 891  !%% only_one_kpt=.FALSE.
 892  !%% call gsph_init(Gsph,Cryst,Sigp%npwc,gvec=Er%gvec)
 893  !%% deallocate(gvec_p)
 894 
 895  ! Default values used if external file is not read.
 896  nqlwl  = 0
 897  nomega = 1
 898  W%npw  = npw_asked
 899 
 900  ! Model dielectric function does not require any external file.
 901  from_file = (W%Info%use_mdf==MDL_NONE)
 902 
 903  if (from_file) then
 904    ! Open file and check its content.
 905    if (endswith(W%fname, ".nc")) W%iomode = IO_MODE_ETSF
 906 
 907    call hscr%from_file(W%fname, fform, comm)
 908    ! Echo of the header
 909    if (my_rank == master .and. W%prtvol>0) call hscr%print()
 910 
 911    ! Communicate the header and copy basic parameters.
 912    !call hscr_bcast(Hscr,master,my_rank,comm)
 913    mat_type_read = Hscr%id
 914    nqlwl         = Hscr%nqlwl
 915    nomega        = Hscr%nomega
 916  end if
 917 
 918  W%nqlwl  = nqlwl
 919  W%nomega = nomega
 920 
 921  ABI_MALLOC(W%qlwl,(3,W%nqlwl))
 922  ABI_MALLOC(W%omega,(W%nomega))
 923 
 924  if (from_file) then
 925    W%qlwl  = Hscr%qlwl
 926    W%omega = Hscr%omega
 927    !
 928    ! G-vectors.
 929    W%npw=Hscr%npwe
 930    if (npw_asked>0) then
 931      if (npw_asked>Hscr%npwe) then
 932        write(msg,'(a,i8,2a,i8)')&
 933         'The number of G-vectors saved on file is less than the value required = ',npw_asked,ch10,&
 934         'Calculation will proceed with the Max available npw = ',Hscr%npwe
 935        ABI_WARNING(msg)
 936      else
 937        W%npw=npw_asked ! Redefine the no. of G"s for W.
 938        write(msg,'(a,i8,2a,i8)')&
 939         'The Number of G-vectors saved on file is larger than the value required = ',npw_asked,ch10,&
 940         'Calculation will proceed with npw = ',W%npw
 941        ABI_COMMENT(msg)
 942      end if
 943    end if
 944    !
 945    ! consistency check on G-vectors and q-points.
 946    if ( ANY(Hscr%gvec(:,1:W%npw) /= Gsph%gvec(:,1:W%npw)) ) then
 947      !write(std_out) W%gvec, Gsph%gvec
 948      ABI_ERROR("Hscr%gvec /= Gsph%gvec(1:W%npw)")
 949    end if
 950    ABI_CHECK(Hscr%nqibz==Qmesh%nibz,"Mismatch in the number of q-points")
 951    ierr=0
 952    do iq_ibz=1,Hscr%nqibz
 953      if (ANY( ABS(Qmesh%ibz(:,iq_ibz) - Hscr%qibz(:,iq_ibz)) > tol6) ) then
 954        ierr=ierr+1
 955        write(std_out,'(i0,2(3f7.3,1x))')iq_ibz,Qmesh%ibz(:,iq_ibz),Hscr%qibz(:,iq_ibz)
 956      end if
 957    end do
 958    ABI_CHECK(ierr==0,"Wrong ordering in q-point list, aborting now")
 959  end if
 960 
 961  ABI_MALLOC(W%gvec,(3,W%npw))
 962  W%gvec = Gsph%gvec(:,1:W%npw)
 963 
 964  ! Frequency mesh.
 965  W%nomega_r=1; W%nomega_i=0
 966  if (W%nomega==2) then
 967    W%nomega_r=1; W%nomega_i=1
 968  else
 969    ! Real frequencies are packed in the first locations.
 970    W%nomega_r=1
 971    do iw=1,W%nomega
 972      if (DBLE(W%omega(iw))>0.001*Ha_eV) W%nomega_r=iw
 973    end do
 974    W%nomega_i=W%nomega-W%nomega_r
 975  end if
 976  !
 977  ! ------------------------------ Initialization completed --------------------------------
 978  !
 979  ! Just to keep the code readable.
 980  npw    = W%npw
 981  nqibz  = W%nqibz
 982  nomega = W%nomega
 983  nI     = W%ni
 984  nJ     = W%nj
 985  ABI_MALLOC(W%Fgg,(nqibz))
 986 
 987  if (from_file) then
 988 
 989    ! Read the ab-initio em1 from file.
 990    select case (mat_type_read)
 991    case (MAT_INV_EPSILON)
 992      call wrtout(std_out,strcat("Em1 will be initialized from SCR file: ",W%fname),"COLL")
 993      !
 994    case  (MAT_CHI0)  ! Write new SCR file.
 995      ABI_ERROR("Not coded yet")
 996      sus_fname = W%fname; scr_fname="TESTING_SUS2SCR"
 997 
 998      W%fname = scr_fname  ! Change the name of the file associated to W.
 999      !
1000    case default
1001      write(msg,'(a,i0)')" Unsupported conversion from mat_type ",mat_type_read
1002      ABI_ERROR(msg)
1003    end select
1004 
1005    do iq_ibz=1,nqibz
1006      if (.not.W%keep_q(iq_ibz)) then
1007         call wrtout(std_out,strcat("Skipping iq_ibz: ",itoa(iq_ibz)),"COLL")
1008         CYCLE
1009      end if
1010      !
1011      nqlwl=0; is_qeq0=(normv(W%qibz(:,iq_ibz),Cryst%gmet,'G')<GW_TOLQ0)
1012      if (is_qeq0) nqlwl=W%nqlwl
1013      !
1014      ! Allocate F_{GG'}(w)
1015      call fgg_init(W%Fgg(iq_ibz),npw,nomega,nqlwl)
1016      !
1017      ! Read data from file (use MPI-IO if possible)
1018      if (W%iomode /= IO_MODE_ETSF .and. xmpi_mpiio==1) then
1019        call wrtout(std_out, "read_screening with MPI_IO")
1020        call read_screening(varname,W%fname,npw,1,nomega,W%Fgg(iq_ibz)%mat,IO_MODE_MPI,comm,iqiA=iq_ibz)
1021      else
1022        call read_screening(varname,W%fname,npw,1,nomega,W%Fgg(iq_ibz)%mat,W%iomode,comm,iqiA=iq_ibz)
1023      end if
1024      W%Fgg(iq_ibz)%has_mat = MAT_STORED
1025      ! W contains Em1 and is ready to use.
1026    end do
1027 
1028  else
1029    !
1030    ! Model dielectric function. Only epsm-1 is supported here.
1031    call wrtout(std_out," Calculating model dielectric function... ","COLL")
1032    ABI_CHECK(W%nomega==1,"Cannot use nomega>1 in model dielectric function")
1033    !
1034    do iq_ibz=1,nqibz
1035      !
1036      if (.not.W%keep_q(iq_ibz)) CYCLE
1037      !
1038      ! The wings are not used here.
1039      nqlwl=0; is_qeq0= (normv(W%qibz(:,iq_ibz),Cryst%gmet,'G')<GW_TOLQ0)
1040      !
1041      ! Calculate the model. Note that mdielf awaits an index in the BZ.
1042      found = qmesh%has_bz_item(Qmesh%ibz(:,iq_ibz),iq_bz,g0)
1043      if (.not.found.or.ANY(g0/=0)) then
1044        ABI_ERROR("Problem in retrieving ibz point")
1045      end if
1046      !
1047      ! Allocate F_{GG'}(w).
1048      call fgg_init(W%Fgg(iq_ibz),npw,nomega,nqlwl)
1049 
1050      eps_inf  =  W%Info%eps_inf
1051      mdf_type =  W%Info%use_mdf
1052      !em1_ggw  => W%Fgg(iq_ibz)%mat
1053 
1054      ! Construct W TODO check the new implementation.
1055      call screen_mdielf(iq_bz,npw,nomega,mdf_type,eps_inf,Cryst,Qmesh,Vcp,Gsph,&
1056 &      nspden,nfftf_tot,ngfftf,ae_rhor,"EM1",W%Fgg(iq_ibz)%mat,comm)
1057 
1058      W%Fgg(iq_ibz)%has_mat = MAT_STORED
1059 
1060      if (W%prtvol > 0) then
1061        do iw=1,nomega
1062          write(msg,'(a,i3,a,i4,a)')'  Model symmetrical e^{-1} (q=',iq_ibz,', omega=',iw,', G,G'')'
1063          call wrtout(std_out,msg,'COLL')
1064          call print_arr(W%Fgg(iq_ibz)%mat(:,:,iw))
1065        end do
1066      end if
1067      !
1068    end do ! iq_ibz
1069    !
1070  end if
1071  !
1072  ! Init plasmon-pole parameters.
1073  if (W%has_ppmodel>0) then
1074    ABI_WARNING("Calculating PPmodel parameters")
1075    ppmodel = W%Info%use_ppm; drude_plsmf = W%Info%drude_plsmf
1076    call ppm_init(W%PPm,W%mqmem,W%nqibz,W%npw,ppmodel,drude_plsmf,W%Info%invalid_freq)
1077 
1078    do iq_ibz=1,nqibz
1079      if (screen_ihave_fgg(W,iq_ibz,how="Stored")) then
1080        !em1_ggw => W%Fgg(iq_ibz)%mat
1081        call new_setup_ppmodel(W%PPm,iq_ibz,Cryst,Qmesh,npw,nomega,W%omega,W%Fgg(iq_ibz)%mat,&
1082           nfftf_tot,Gsph%gvec,ngfftf,W%ae_rhor(:,1))
1083      end if
1084    end do
1085  end if
1086  !
1087  ! Deallocate Fgg if the matrices are not needed anymore.
1088  if (deallocate_Fgg) then
1089    call screen_fgg_qbz_set(W,0,0,"Pointer") ! Avoid dangling pointer.
1090    call fgg_free(W%Fgg,keep_q=W%keep_q)
1091  end if
1092 
1093  if (from_file) call Hscr%free()
1094 
1095  DBG_EXIT("COLL")
1096 
1097 end subroutine screen_init

m_screen/screen_nullify [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

 screen_nullify

FUNCTION

 Initialize the pointers to null()

SIDE EFFECTS

 W<screen_t>=The data structure to be nullified.

SOURCE

652 subroutine screen_nullify(W)
653 
654 !Arguments ------------------------------------
655 !scalars
656  type(screen_t),intent(inout) :: W
657 ! *************************************************************************
658 
659  !@screen_t
660  !
661  !types
662  nullify(W%Fgg_qbz); W%fgg_qbz_stat=FGG_QBZ_ISPOINTER ! Nedded since the initial status is undefined.
663  nullify(W%Fgg)
664 
665  call ppm_nullify(W%PPm)
666 
667 end subroutine screen_nullify

m_screen/screen_symmetrizer [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

  screen_symmetrizer

FUNCTION

  Modify the status of the object so that the symmetrized component F(q_bz)_GG' is calculated
  (if needed) and is made available in the internal buffer. This routine must be called before
  performing any operation that involves the symmetrized component of the two-point function.

INPUTS

  W<screen_t>=Data structure used to represent the two-point function
  Cryst<crystal_t>=Info on the crystal structure.
  Gsph<gsphere_t>=data related to the G-sphere
  Qmesh<kmesh_t>=Structure defining the q-mesh used for sample the BZ.
  iq_bz=Index of the q-point in the BZ where F(q_bz)_GG' is wanted.

OUTPUT

SIDE EFFECTS

   W%PPm
   W%Fgg_qbz

NOTES

  In the present implementation we are not considering a possible umklapp vector G0 in the
  expression Sq = q+G0. Treating this case would require some changes in the G-sphere
  since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required
  to reconstruct the BZ.

SOURCE

1132 subroutine screen_symmetrizer(W,iq_bz,Cryst,Gsph,Qmesh,Vcp)
1133 
1134 !Arguments ------------------------------------
1135 !scalars
1136  integer,intent(in) :: iq_bz
1137  type(screen_t),intent(inout) :: W
1138  type(crystal_t),intent(in) :: Cryst
1139  type(gsphere_t),intent(in) :: Gsph
1140  type(kmesh_t),intent(in) :: Qmesh
1141  type(vcoul_t),intent(in) :: Vcp
1142 
1143 !Local variables-------------------------------
1144 !scalars
1145  integer,parameter :: nqlwl0=0
1146  integer :: iq_ibz,isym_q,itim_q,npw,nomega,nqibz,mdf_type
1147  real(dp) :: eps_inf
1148  logical :: q_isirred
1149  !character(len=500) :: msg
1150 !arrays
1151  real(dp) :: qbz(3)
1152 
1153 ! *********************************************************************
1154 
1155  DBG_ENTER("COLL")
1156 
1157  npw = W%npw; nqibz = W%nqibz; nomega = W%nomega
1158 
1159  call qmesh%get_bz_item(iq_bz,qbz,iq_ibz,isym_q,itim_q,isirred=q_isirred)
1160  !
1161  ! ========================================================
1162  ! ==== Branching for in-core or out-of-core solutions ====
1163  ! ========================================================
1164  if (screen_ihave_fgg(W,iq_ibz,how="Stored")) then
1165 
1166    if (q_isirred) then ! Symmetrization is not needed. Point the data in memory.
1167      call screen_fgg_qbz_set(W,iq_bz,nqlwl0,"Pointer")
1168      W%Fgg_qbz => W%Fgg(iq_ibz)
1169    else
1170      ! Allocate space. ! TODO Wings are not symmetrized but oh well
1171      call screen_fgg_qbz_set(W,iq_bz,nqlwl0,"Allocate")  ! Dimensions should not be changed.
1172 
1173      ! Out-of-place symmetrization.
1174      !em1_qibz => W%Fgg(iq_ibz)%mat; em1_qbz  => W%Fgg_qbz%mat
1175      call em1_symmetrize_op(iq_bz,npw,nomega,Gsph,Qmesh,W%Fgg(iq_ibz)%mat,W%Fgg_qbz%mat)
1176    end if
1177 
1178    if (W%has_ppmodel>0) then
1179      ! Symmetrize the ppmodel tables: em1_qibz => W%Fgg(iq_ibz)%mat
1180      call ppm_symmetrizer(W%PPm,iq_bz,Cryst,Qmesh,Gsph,npw,nomega,W%omega,W%Fgg(iq_ibz)%mat,&
1181 &      W%nfftf_tot,W%ngfftf,W%ae_rhor(:,1))
1182    end if
1183 
1184  else if (screen_ihave_fgg(W,iq_ibz,how="Allocated")) then
1185    ABI_ERROR("Fgg_iqibz is allocated but not initialized!")
1186 
1187  else
1188    if (W%fgg_qbz_idx /= iq_bz) then
1189      ! Must compute em1_qbz here. em1_qbz => W%Fgg_qbz%mat
1190 
1191      ! Allocate the BZ buffer.
1192      call screen_fgg_qbz_set(W,iq_bz,nqlwl0,"Allocate")
1193 
1194      if (W%Info%use_mdf /= MDL_NONE) then
1195        ! Compute the model-dielectric function at qbz on-the fly and in sequential
1196        !call wrtout(std_out,"Will compute MDF on the fly","COLL")
1197        call screen_mdielf(iq_bz,npw,nomega,W%Info%use_mdf,W%Info%eps_inf,Cryst,Qmesh,Vcp,Gsph,&
1198 &         W%nspden,W%nfftf_tot,W%ngfftf,W%ae_rhor,"EM1",W%Fgg_qbz%mat,xmpi_comm_self)
1199 
1200      else
1201        ! Read W(q_ibz) and symmetrize it (do this only if we don't have the correct q_bz in memory).
1202        call wrtout(std_out,sjoin("Out of core with file: ",W%fname),"COLL")
1203        call read_screening(em1_ncname,W%fname,npw,1,nomega,W%Fgg_qbz%mat,W%iomode,xmpi_comm_self,iqiA=iq_ibz)
1204 
1205        ! In-place symmetrization to get the q-point in the BZ.
1206        if (.not.q_isirred) then
1207          call em1_symmetrize_ip(iq_bz,npw,nomega,Gsph,Qmesh,W%Fgg_qbz%mat)
1208        end if
1209      end if
1210 
1211      W%Fgg_qbz%has_mat = MAT_STORED
1212    end if
1213 
1214    ABI_CHECK(W%Fgg_qbz%has_mat == MAT_STORED, "Wrong has_mat")
1215    !
1216    ! Ppmodel calculations with ppm tables in memory.
1217    ! TODO treat the case in which IBZ tables are stored in memory.
1218    if (W%has_ppmodel>0) then
1219      ABI_ERROR("Not implemented error")
1220      ! Symmetrize the ppmodel using em1_qibz.
1221      !%call ppm_symmetrizer(W%PPm,iq_bz,Cryst,Qmesh,Gsph,npw,nomega,W%omega,em1_qibz,W%nfftf_tot,W%ngfftf,W%ae_rhor(:,1))
1222      call ppm_symmetrizer(W%PPm,iq_bz,Cryst,Qmesh,Gsph,npw,nomega,W%omega,W%Fgg_qbz%mat ,W%nfftf_tot,W%ngfftf,W%ae_rhor(:,1))
1223    end if
1224  end if
1225 
1226  ! Calculate model dielectric function for this q-point in the BZ.
1227  eps_inf = W%Info%eps_inf; mdf_type = W%Info%use_mdf
1228 
1229  ! Model dielectric function. Only epsm-1 is supported here.
1230  !call wrtout(std_out," Calculating model dielectric function... ","COLL")
1231  !ABI_CHECK(W%nomega==1,"Cannot use nomega>1 in model dielectric function")
1232 
1233  ! W%Fgg_qbz%mat
1234 
1235  !%  call screen_mdielf(iq_bz,npw,nomega,mdf_type,eps_inf,Cryst,Qmesh,Vcp,Gsph,&
1236  !% &  W%nspden,W%nfftf_tot,W%ngfftf,W%ae_rhor,"EM1",em1_qbz,xmpi_comm_self)
1237 
1238  ! Store the index of the q-point in the BZ for checking purpose.
1239  W%fgg_qbz_idx = iq_bz
1240 
1241  DBG_EXIT("COLL")
1242 
1243 end subroutine screen_symmetrizer

m_screen/screen_t [ Types ]

[ Top ] [ m_screen ] [ Types ]

NAME

 screen_t

FUNCTION

  Object used to store the screening matrix in reciprocal space.

SOURCE

221  type,public :: screen_t
222 
223 ! scalars
224   integer :: iomode               ! Flag defining the IO mode.
225   integer :: debug_level=0        ! Internal Flag defining the debug level.
226   integer :: mqmem                ! =0 for out-of-core solution, =nqibz if entire matrix is stored in memory.
227   integer :: nI,nJ                ! Number of components (rows,columns) in chi|eps^-1. (1,1) if collinear.
228   integer :: nqibz                ! Number of q-points in the IBZ used.
229   integer :: nqlwl                ! Number of points used for the treatment of the long wave-length limit.
230   integer :: nomega               ! Total Number of frequencies used.
231   integer :: nomega_i             ! Number of purely imaginary frequencies used.
232   integer :: nomega_r             ! Number of real frequencies used.
233   integer :: npw                  ! Number of G vectors.
234   integer :: prtvol               ! Verbosity level.
235   integer :: has_ppmodel          ! 1 if PPmodel tables are stored.
236   integer :: has_fgg              ! 1 if Fgg tables are stored.
237   !%integer :: wmesh_type
238   !%real(dp) :: ecut
239   integer :: nfftf_tot
240   integer :: nspden
241 
242 !arrays
243   integer :: ngfftf(18)          ! Info on the FFT mesh used for ae_rhor (used for the model dielectric function)
244 
245   real(dp),allocatable :: ae_rhor(:,:)
246   ! ae_rhor(nfft,nspden)
247   ! Density in real space used to construct the TDDFT kernel or the model dielectric function.
248   ! NOTE that ae_rhor is given on the dense mesh as it contains the PAW onsite contribution.
249 
250   character(len=fnlen) :: fname=ABI_NOFILE  ! Name of the file used for the out-of-core solution.
251 
252   real(dp),allocatable :: qibz(:,:)
253   ! qibz(3,nqibz)
254   ! q-points in reduced coordinates
255 
256   !type(kmesh_t) :: Qmesh
257   ! Info the q-point sampling.
258 
259   real(dp),allocatable :: qlwl(:,:)
260   ! qlwl(3,nqlwl)
261   ! q-points used for the long wave-length limit treatment.
262 
263   complex(dpc),allocatable :: omega(:)
264   ! omega(nomega)
265   ! List of frequencies. Real frequencies are packed first.
266 
267   integer,allocatable :: gvec(:,:)
268   ! gvec(3,npw)
269   ! G-vectors used to describe the two-point function (r.l.u.).
270 
271   !type(gsphere_t) :: Gsphere
272   ! Info on the G-sphere. Note that the basis set does not depend on q.
273   ! The G-vectors are ordered in shells to facilitate the application of symmetry operations.
274   ! See m_gsphere.F90.
275 
276   logical,allocatable :: keep_q(:)
277    ! keep_q(nqibz)
278    ! Storage strategy: keep or not keep Em1(q) in memory.
279 
280   type(fgg_t),pointer :: Fgg(:)
281   ! Fgg(nqibz)
282   ! F_{G,G'}(q,w) for q in the IBZ.
283 
284   integer :: fgg_qbz_stat=FGG_QBZ_ISPOINTER
285   ! Status of Fgg_qbz
286 
287   integer :: fgg_qbz_idx=0
288   ! The index of the q-point in BZ pointed by Fgg_qbz. Used for debugging purpose.
289 
290   type(fgg_t),pointer :: Fgg_qbz  => null()
291   ! Buffer used for storing F_GG' at the point q_bz in the BZ
292   ! If q_bz is in the IBZ, Fgg_qbz *points* to Fgg(iq_ibz)
293   ! If q_bz is not in the IBZ, Fgg_qbz is *allocated* and used to store the symmetrized matrix.
294 
295   type(ppmodel_t) :: PPm
296   ! Structure storing the plasmon-pole parameters.
297 
298   type(screen_info_t) :: Info
299   ! Parameters used to construct the screening.
300 
301  end type screen_t
302 
303  public :: screen_nullify       ! Nullify all pointers before use.
304  public :: screen_init          ! Creation method.
305  !public :: screen_print         ! Print info
306  public :: screen_free          ! Free all associated pointers.
307  public :: screen_symmetrizer   ! Prepare the object for applying W_qbz.
308  public :: screen_w0gemv        ! Matrix vector multiplication \sum_{G'} F_{G,G') |u(G')>.
309  !public :: screen_times_ket     ! Compute \Sigma_c(\omega)|\phi> in reciprocal space.

m_screen/screen_w0gemv [ Functions ]

[ Top ] [ m_screen ] [ Functions ]

NAME

 screen_w0gemv

FUNCTION

  Perform the matrix multiplication W x vector in reciprocal space.

INPUTS

  W<screen_t>=
  in_npw=Number of G vectors in in_ket
  nspinor=Number of spinorial components.
  in_ket(in_npw)= |\phi> in reciprocal space.
  trans= On entry, TRANS specifies the operation to be performed as follows:
  TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
  TRANS = 'T' or 't'   y := alpha*A**T*x + beta*y.
  TRANS = 'C' or 'c'   y := alpha*A**H*x + beta*y.

OUTPUT

   out_ket(in_npw)= W |\phi\> in reciprocal space.
   ZGEMV  performs one of the matrix-vector operations
   *
   *     y := alpha*A*x + beta*y,   or   y := alpha*A**T*x + beta*y,   or
   *
   *     y := alpha*A**H*x + beta*y,
   *
   *  where alpha and beta are scalars, x and y are vectors and A is an m by n matrix.

SOURCE

1277 subroutine screen_w0gemv(W,trans,in_npw,nspinor,only_diago,alpha,beta,in_ket,out_ket)
1278 
1279 !Arguments ------------------------------------
1280 !scalars
1281  integer,intent(in) :: in_npw,nspinor
1282  complex(gwpc),intent(in) :: alpha,beta
1283  logical,intent(in) :: only_diago
1284  character(len=*),intent(in) ::  trans
1285  type(screen_t),intent(in) :: W
1286 !arrays
1287  complex(gwpc),intent(in) :: in_ket(in_npw*nspinor)
1288  complex(gwpc),intent(out) :: out_ket(in_npw*nspinor)
1289 
1290 !Local variables-------------------------------
1291 !scalars
1292  integer :: ig,lda
1293 !arrays
1294  complex(gwpc),ABI_CONTIGUOUS pointer :: em1_qbz(:,:)
1295 
1296 ! *************************************************************************
1297 
1298  lda = W%npw
1299  em1_qbz => W%Fgg_qbz%mat(:,:,1)
1300 
1301  if (.not.only_diago) then
1302    !out_ket = MATMUL(TRANSPOSE(CONJG(em1_qbz(1:in_npw,1:in_npw)),in_ket)
1303    call xgemv(trans,in_npw,in_npw,alpha,em1_qbz,lda,in_ket,1,beta,out_ket,1)
1304 
1305  else
1306    if (beta /= czero_gw) then
1307      !
1308      if ( firstchar(trans,(/"C"/)) ) then
1309        do ig=1,in_npw
1310          out_ket(ig) = alpha * CONJG(em1_qbz(ig,ig)) * in_ket(ig) + beta * out_ket(ig)
1311        end do
1312      else if ( firstchar(trans,(/"N","T"/)) ) then
1313        do ig=1,in_npw
1314          out_ket(ig) = alpha * em1_qbz(ig,ig) * in_ket(ig) + beta * out_ket(ig)
1315        end do
1316      else
1317        ABI_ERROR(sjoin("Wrong trans:", trans))
1318      end if
1319 
1320    else
1321      ! beta==0
1322      if ( firstchar(trans,(/"C"/)) ) then
1323        do ig=1,in_npw
1324          out_ket(ig) = alpha * CONJG(em1_qbz(ig,ig)) * in_ket(ig)
1325        end do
1326      else if ( firstchar(trans,(/"N","T"/)) ) then
1327        do ig=1,in_npw
1328          out_ket(ig) = alpha * em1_qbz(ig,ig) * in_ket(ig)
1329        end do
1330      else
1331        ABI_ERROR(sjoin("Wrong trans:", trans))
1332      end if
1333 
1334    end if
1335  end if
1336 
1337 end subroutine screen_w0gemv