TABLE OF CONTENTS


ABINIT/m_screening [ Modules ]

[ Top ] [ Modules ]

NAME

  m_screening

FUNCTION

  This module contains the definition of the object used to deal
  with the inverse dielectric matrix as well as related methods.

COPYRIGHT

 Copyright (C) 2008-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

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_screening
24 
25  use defs_basis
26  use m_abicore
27  use m_hide_blas
28  use m_linalg_interfaces
29  use m_xmpi
30  use m_errors
31  use m_copy
32  use m_splines
33  use m_lebedev
34  use m_spectra
35  use m_nctk
36  use m_distribfft
37  use netcdf
38 
39  use defs_abitypes,     only : MPI_type
40  use m_gwdefs,          only : GW_TOLQ0, czero_gw, GW_Q0_DEFAULT
41  use m_fstrings,        only : toupper, endswith, sjoin, itoa
42  use m_io_tools,        only : open_file
43  use m_numeric_tools,   only : print_arr, hermitianize
44  use m_special_funcs,   only : k_fermi, k_thfermi
45  use m_geometry,        only : normv, vdotw, metric
46  use m_hide_lapack,     only : xginv
47  use m_crystal,         only : crystal_t
48  use m_bz_mesh,         only : kmesh_t, box_len
49  use m_fft_mesh,        only : g2ifft
50  use m_fftcore,         only : kgindex
51  use m_fft,             only : fourdp
52  use m_gsphere,         only : gsphere_t
53  use m_vcoul,           only : vcoul_t
54  use m_io_screening,    only : hscr_io, read_screening, write_screening, &
55                                HSCR_LATEST_HEADFORM, hscr_t, ncname_from_id, em1_ncname
56  use m_paw_sphharm,     only : ylmc
57  use m_mpinfo,          only : destroy_mpi_enreg, initmpi_seq
58 
59  implicit none
60 
61  private

m_screening/atddft_hyb_symepsm1 [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 atddft_hyb_symepsm1

FUNCTION

  Calculate $\tilde\epsilon^{-1}$ using ALDA within TDDFT

  Based on atddft_symepsm1, modified for the LR+ALDA hybrid kernel

  Output the electron energy loss function and the macroscopic dielectric function with and
  without local field effects (only if non-zero real frequencies are available)

INPUTS

  iqibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
  npwe=Number of G-vectors in chi0.
  nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
  dim_wing=Dimension of the wings (0 or 3 if q-->0)
  kxcg_mat_sr=Short-range fxc kernel used in the TE epsilon^-1
  option_test=Only for TDDFT:
   == 0 for TESTPARTICLE ==
   == 1 for TESTELECTRON ==
  Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
   %nqibz=Number of q-points.
   %qibz(3,nqibz)=q-points in the IBZ.
  comm=MPI communicator.
  chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
  chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
  chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)

OUTPUT

SIDE EFFECTS

  chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output
   the symmetrized inverse dielectric matrix.

SOURCE

2258 subroutine atddft_hyb_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,kxcg_mat_sr,option_test,my_nqlwl,dim_wing,omega,&
2259 & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
2260 
2261 !Arguments ------------------------------------
2262 !scalars
2263  integer,intent(in) :: iqibz,nI,nJ,npwe,dim_wing,my_nqlwl
2264  integer,intent(in) :: option_test,comm
2265  type(vcoul_t),target,intent(in) :: Vcp
2266 !arrays
2267  complex(gwpc),intent(in) :: kxcg_mat(npwe,npwe), kxcg_mat_sr(npwe,npwe)
2268  complex(dpc),intent(in) :: omega
2269  complex(dpc),intent(inout) :: chi0_lwing(npwe*nI,dim_wing)
2270  complex(dpc),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing)
2271  complex(dpc),intent(inout) :: chi0_head(dim_wing,dim_wing)
2272  complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2273  real(dp),intent(out) :: eelf(my_nqlwl)
2274  complex(dpc),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2275 
2276 !Local variables-------------------------------
2277 !scalars
2278  integer,parameter :: master=0,prtvol=0
2279  integer :: ig1,ig2,my_rank,nprocs,ierr
2280  real(dp) :: ucvol, Zr
2281  logical :: is_qeq0
2282  character(len=500) :: msg
2283 !arrays
2284  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2285  complex(gwpc),allocatable :: chitmp(:,:)
2286  complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:)
2287 
2288 ! *************************************************************************
2289 
2290  ABI_UNUSED(chi0_head(1,1))
2291  ABI_UNUSED(chi0_lwing(1,1))
2292  ABI_UNUSED(chi0_uwing(1,1))
2293 
2294  if (nI/=1.or.nJ/=1) then
2295    ABI_ERROR("nI or nJ=/1 not yet implemented")
2296  end if
2297 
2298  ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
2299  ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded")
2300 
2301  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2302 
2303  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2304 
2305  is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0)
2306 
2307  if (iqibz==1) then
2308    !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl)  ! Use Coulomb term for q-->0
2309    vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! TODO add treatment of non-Analytic behavior
2310  else
2311    vc_sqrt => Vcp%vc_sqrt(:,iqibz)
2312  end if
2313 
2314  write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwpc*b2Mb," Mb"
2315  ABI_MALLOC_OR_DIE(chitmp,(npwe,npwe), ierr)
2316  !
2317  ! * Calculate chi0*fxc.
2318  chitmp = MATMUL(chi0,kxcg_mat)
2319  ! * First calculate the NLF contribution
2320  do ig1=1,npwe
2321    do ig2=1,npwe
2322      chitmp(ig1,ig2)=-chitmp(ig1,ig2)
2323    end do
2324    chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2325  end do
2326 
2327  call xginv(chitmp,npwe,comm=comm)
2328 
2329  chitmp = MATMUL(chitmp,chi0)
2330  !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2331  chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1)
2332  chitmp(1,1)=chitmp(1,1)+one
2333 
2334  epsm_nlf(1)=chitmp(1,1)
2335 
2336  chitmp = MATMUL(chi0,kxcg_mat)
2337  ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp.
2338  do ig1=1,npwe
2339    do ig2=1,npwe
2340      chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2
2341    end do
2342    chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2343  end do
2344 
2345  ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0.
2346  call xginv(chitmp,npwe,comm=comm)
2347  chitmp=MATMUL(chitmp,chi0)
2348 
2349  ! * Save result, now chi0 contains chi.
2350  chi0=chitmp
2351 
2352  select case (option_test)
2353  case (0)
2354    ! Symmetrized TESTPARTICLE epsilon^-1
2355    call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi')
2356    do ig1=1,npwe
2357      chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)
2358      chi0(ig1,ig1)=one+chi0(ig1,ig1)
2359    end do
2360 
2361  case (1)
2362    ! Symmetrized TESTELECTRON epsilon^-1
2363    call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + Vc*chi + Zr*Kxc_sr*chi')
2364    chitmp=MATMUL(kxcg_mat_sr,chi0)
2365    Zr = 0.78
2366 
2367    ! Perform hermitianization, only valid along the imaginary axis.
2368    if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2369 
2370    do ig1=1,npwe
2371      chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:)*Zr
2372      chi0(ig1,ig1)=one+chi0(ig1,ig1)
2373    end do
2374 
2375  case default
2376    ABI_BUG(sjoin('Wrong option_test:',itoa(option_test)))
2377  end select
2378 
2379  ABI_FREE(chitmp)
2380  !
2381  ! === chi0 now contains symmetrical epsm1 ===
2382  ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w).
2383  epsm_lf(1) =  one/chi0(1,1)
2384  eelf   (1) = -AIMAG(chi0(1,1))
2385 
2386  if (prtvol > 0) then
2387    write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]'
2388    call wrtout(std_out,msg)
2389    call print_arr(chi0,unit=std_out)
2390  end if
2391 
2392 end subroutine atddft_hyb_symepsm1

m_screening/atddft_symepsm1 [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 atddft_symepsm1

FUNCTION

  Calculate $\tilde\epsilon^{-1}$ using ALDA within TDDFT

  2) Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
     calculating these quantities for different small q-directions specified by the user
     (Not yet operative)

  Output the electron energy loss function and the macroscopic dielectric function with and
  without local field effects (only if non-zero real frequencies are available)

INPUTS

  iqibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
  npwe=Number of G-vectors in chi0.
  nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
  dim_wing=Dimension of the wings (0 or 3 if q-->0)
  option_test=Only for TDDFT:
   == 0 for TESTPARTICLE ==
   == 1 for TESTELECTRON ==
  Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
   %nqibz=Number of q-points.
   %qibz(3,nqibz)=q-points in the IBZ.
  comm=MPI communicator.
  chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
  chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
  chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)

OUTPUT

SIDE EFFECTS

  chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output
   the symmetrized inverse dielectric matrix.

SOURCE

2083 subroutine atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,option_test,my_nqlwl,dim_wing,omega,&
2084 & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
2085 
2086 !Arguments ------------------------------------
2087 !scalars
2088  integer,intent(in) :: iqibz,nI,nJ,npwe,dim_wing,my_nqlwl
2089  integer,intent(in) :: option_test,comm
2090  type(vcoul_t),target,intent(in) :: Vcp
2091 !arrays
2092  complex(gwpc),intent(in) :: kxcg_mat(npwe,npwe)
2093  complex(dpc),intent(in) :: omega
2094  complex(dpc),intent(inout) :: chi0_lwing(npwe*nI,dim_wing)
2095  complex(dpc),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing)
2096  complex(dpc),intent(inout) :: chi0_head(dim_wing,dim_wing)
2097  complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2098  real(dp),intent(out) :: eelf(my_nqlwl)
2099  complex(dpc),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2100 
2101 !Local variables-------------------------------
2102 !scalars
2103  integer,parameter :: master=0,prtvol=0
2104  integer :: ig1,ig2,my_rank,nprocs,ierr
2105  real(dp) :: ucvol
2106  logical :: is_qeq0
2107  character(len=500) :: msg
2108 !arrays
2109  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2110  complex(gwpc),allocatable :: chitmp(:,:)
2111  complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:)
2112 
2113 ! *************************************************************************
2114 
2115  ABI_UNUSED(chi0_head(1,1))
2116  ABI_UNUSED(chi0_lwing(1,1))
2117  ABI_UNUSED(chi0_uwing(1,1))
2118 
2119  if (nI/=1.or.nJ/=1) then
2120    ABI_ERROR("nI or nJ=/1 not yet implemented")
2121  end if
2122 
2123  ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
2124  ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded")
2125 
2126  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2127 
2128  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2129 
2130  is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0)
2131 
2132  if (iqibz==1) then
2133    !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl)  ! Use Coulomb term for q-->0
2134    vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! TODO add treatment of non-Analytic behavior
2135  else
2136    vc_sqrt => Vcp%vc_sqrt(:,iqibz)
2137  end if
2138 
2139  write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwpc*b2Mb," Mb"
2140  ABI_MALLOC_OR_DIE(chitmp, (npwe,npwe), ierr)
2141  !
2142  ! Calculate chi0*fxc.
2143  chitmp = MATMUL(chi0,kxcg_mat)
2144  ! First calculate the NLF contribution
2145  do ig1=1,npwe
2146    do ig2=1,npwe
2147      chitmp(ig1,ig2)=-chitmp(ig1,ig2)
2148    end do
2149    chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2150  end do
2151 
2152  call xginv(chitmp,npwe,comm=comm)
2153 
2154  chitmp = MATMUL(chitmp,chi0)
2155  !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2156  chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1)
2157  chitmp(1,1)=chitmp(1,1)+one
2158 
2159  epsm_nlf(1)=chitmp(1,1)
2160 
2161  chitmp = MATMUL(chi0,kxcg_mat)
2162  ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp.
2163  do ig1=1,npwe
2164    do ig2=1,npwe
2165      chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2
2166    end do
2167    chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2168  end do
2169 
2170  ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0.
2171  call xginv(chitmp,npwe,comm=comm)
2172  chitmp=MATMUL(chitmp,chi0)
2173 
2174  ! * Save result, now chi0 contains chi.
2175  chi0=chitmp
2176 
2177  select case (option_test)
2178  case (0)
2179    ! Symmetrized TESTPARTICLE epsilon^-1
2180    call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi')
2181    do ig1=1,npwe
2182      chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)
2183      chi0(ig1,ig1)=one+chi0(ig1,ig1)
2184    end do
2185 
2186  case (1)
2187    ! Symmetrized TESTELECTRON epsilon^-1
2188    call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + (Vc + fxc)*chi')
2189    chitmp=MATMUL(kxcg_mat,chi0)
2190 
2191    ! Perform hermitianization, only valid along the imaginary axis.
2192    if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2193 
2194    do ig1=1,npwe
2195      chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:)
2196      chi0(ig1,ig1)=one+chi0(ig1,ig1)
2197    end do
2198 
2199  case default
2200    ABI_BUG(sjoin('Wrong option_test:',itoa(option_test)))
2201  end select
2202 
2203  ABI_FREE(chitmp)
2204  !
2205  ! === chi0 now contains symmetrical epsm1 ===
2206  ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w).
2207  epsm_lf(1) =  one/chi0(1,1)
2208  eelf   (1) = -AIMAG(chi0(1,1))
2209 
2210  if (prtvol > 0) then
2211    write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]'
2212    call wrtout(std_out,msg)
2213    call print_arr(chi0,unit=std_out)
2214  end if
2215 
2216 end subroutine atddft_symepsm1

m_screening/chi_free [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 chi_free

FUNCTION

INPUTS

OUTPUT

SOURCE

3069 subroutine chi_free(chi)
3070 
3071 !Arguments ------------------------------------
3072 !scalars
3073  type(chi_t),intent(inout) :: chi
3074 
3075 ! *************************************************************************
3076 
3077  if (allocated(chi%mat)) then
3078    ABI_FREE(chi%mat)
3079  end if
3080  if (allocated(chi%head)) then
3081    ABI_FREE(chi%head)
3082  end if
3083  if (allocated(chi%lwing)) then
3084    ABI_FREE(chi%lwing)
3085  end if
3086  if (allocated(chi%uwing)) then
3087    ABI_FREE(chi%uwing)
3088  end if
3089 
3090 end subroutine chi_free

m_screening/chi_new [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 chi_new

FUNCTION

INPUTS

OUTPUT

SOURCE

3037 type(chi_t) function chi_new(npwe, nomega) result(chi)
3038 
3039 !Arguments ------------------------------------
3040  integer,intent(in) :: npwe,nomega
3041 
3042 ! *************************************************************************
3043 
3044  chi%nomega = nomega; chi%npwe = npwe
3045 
3046  ABI_MALLOC(chi%mat, (npwe,npwe,nomega))
3047 
3048  ABI_MALLOC(chi%head, (3,3,nomega))
3049  ABI_MALLOC(chi%lwing, (npwe, nomega,3))
3050  ABI_MALLOC(chi%uwing, (npwe, nomega,3))
3051 
3052 end function chi_new

m_screening/chi_t [ Types ]

[ Top ] [ m_screening ] [ Types ]

NAME

  chi_t

FUNCTION

  This object contains the head and the wings of the polarizability
  These quantities are used to treat the q-->0 limit

SOURCE

179  type,public :: chi_t
180 
181    integer :: npwe
182    ! Number of G vectors.
183 
184    integer :: nomega
185    ! Number of frequencies
186 
187    complex(gwpc),allocatable :: mat(:,:,:)
188    ! mat(npwe, npwe, nomega)
189 
190    complex(dpc),allocatable :: head(:,:,:)
191    ! head(3,3,nomega)
192 
193    complex(dpc),allocatable :: lwing(:,:,:)
194    ! lwing(3,npwe,nomega)
195    ! Lower wings
196 
197    complex(dpc),allocatable :: uwing(:,:,:)
198    ! uwing(3,npwe,nomega)
199    ! Upper wings.
200 
201  end type chi_t
202 
203  public :: chi_new    ! Create new object (allocate memory)
204  public :: chi_free   ! Free memory.

m_screening/decompose_epsm1 [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 decompose_epsm1

FUNCTION

 Decompose the complex symmetrized dielectric

INPUTS

OUTPUT

SOURCE

1179 subroutine decompose_epsm1(Er,iqibz,eigs)
1180 
1181 !Arguments ------------------------------------
1182 !scalars
1183  integer,intent(in) :: iqibz
1184  type(Epsilonm1_results),intent(in) :: Er
1185 !arrays
1186  complex(dpc),intent(out) :: eigs(Er%npwe,Er%nomega)
1187 
1188 !Local variables-------------------------------
1189 !scalars
1190  integer :: info,lwork,iw,negw,ig1,ig2,idx,sdim,npwe,ierr
1191  character(len=500) :: msg
1192 !arrays
1193  real(dp),allocatable :: ww(:),rwork(:)
1194  complex(dpc),allocatable :: work(:),Adpp(:),eigvec(:,:),Afull(:,:),vs(:,:),wwc(:)
1195  logical,allocatable :: bwork(:)
1196  logical :: sortcplx !BUG in abilint
1197 
1198 ! *********************************************************************
1199 
1200  ABI_CHECK(Er%mqmem/=0,'mqmem==0 not implemented')
1201 
1202  npwe = Er%npwe
1203 
1204  do iw=1,Er%nomega
1205 
1206    if (ABS(REAL(Er%omega(iw)))>0.00001) then
1207      ! Eigenvalues for a generic complex matrix
1208      lwork=4*2*npwe
1209      ABI_MALLOC(wwc,(npwe))
1210      ABI_MALLOC(work,(lwork))
1211      ABI_MALLOC(rwork,(npwe))
1212      ABI_MALLOC(bwork,(npwe))
1213      ABI_MALLOC(vs,(npwe,npwe))
1214      ABI_MALLOC(Afull,(npwe,npwe))
1215 
1216      Afull=Er%epsm1(:,:,iw,iqibz)
1217 
1218      !for the moment no sort, maybe here I should sort using the real part?
1219      call ZGEES('V','N',sortcplx,npwe,Afull,npwe,sdim,wwc,vs,npwe,work,lwork,rwork,bwork,info)
1220      if (info/=0) then
1221        ABI_ERROR(sjoin("ZGEES returned info:",itoa(info)))
1222      end if
1223 
1224      eigs(:,iw)=wwc(:)
1225 
1226      ABI_FREE(wwc)
1227      ABI_FREE(work)
1228      ABI_FREE(rwork)
1229      ABI_FREE(bwork)
1230      ABI_FREE(vs)
1231      ABI_FREE(Afull)
1232 
1233    else
1234      ! Hermitian version.
1235      lwork=2*npwe-1
1236      ABI_MALLOC(ww,(npwe))
1237      ABI_MALLOC(work,(lwork))
1238      ABI_MALLOC(rwork,(3*npwe-2))
1239      ABI_MALLOC(eigvec,(npwe,npwe))
1240      ABI_MALLOC_OR_DIE(Adpp,(npwe*(npwe+1)/2), ierr)
1241 
1242      idx=0 ! Pack the matrix
1243      do ig2=1,npwe
1244        do ig1=1,ig2
1245          idx=idx+1
1246          Adpp(idx)=Er%epsm1(ig1,ig2,iw,iqibz)
1247        end do
1248      end do
1249 
1250      ! For the moment we require also the eigenvectors.
1251      call ZHPEV('V','U',npwe,Adpp,ww,eigvec,npwe,work,rwork,info)
1252      if (info/=0) then
1253        ABI_ERROR(sjoin('ZHPEV returned info=', itoa(info)))
1254      end if
1255 
1256      negw=(COUNT((REAL(ww)<tol6)))
1257      if (negw/=0) then
1258        write(msg,'(a,i5,a,i3,a,f8.4)')&
1259         'Found negative eigenvalues. No. ',negw,' at iqibz= ',iqibz,' minval= ',MINVAL(REAL(ww))
1260        ABI_WARNING(msg)
1261      end if
1262 
1263      eigs(:,iw)=ww(:)
1264 
1265      ABI_FREE(ww)
1266      ABI_FREE(work)
1267      ABI_FREE(rwork)
1268      ABI_FREE(eigvec)
1269      ABI_FREE(Adpp)
1270    end if
1271  end do !iw
1272 
1273 ! contains
1274 ! function sortcplx(carg) result(res)
1275 !  implicit none
1276 !  complex(dpc),intent(in) :: carg
1277 !  logical :: res
1278 !  res=.TRUE.
1279 ! end function sortcplx
1280 
1281 end subroutine decompose_epsm1

m_screening/em1results_free [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 em1results_free

FUNCTION

 Deallocate all the pointers in Er that result to be associated.
 Perform also a cleaning of the Header.

INPUTS

OUTPUT

SOURCE

274 subroutine em1results_free(Er)
275 
276 !Arguments ------------------------------------
277 !scalars
278  type(Epsilonm1_results),intent(inout) :: Er
279 ! *************************************************************************
280 
281  !@Epsilonm1_results
282  !integer
283  ABI_SFREE(Er%gvec)
284 
285  !real
286  ABI_SFREE(Er%qibz)
287  ABI_SFREE(Er%qlwl)
288 
289  !complex
290  ABI_SFREE(Er%epsm1)
291  ABI_SFREE(Er%omega)
292 
293  !datatypes
294  call Er%Hscr%free()
295 
296 end subroutine em1results_free

m_screening/em1results_print [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  em1results_print

FUNCTION

 Print the basic dimensions and the most important
 quantities reported in the Epsilonm1_results data type.

INPUTS

  Er<Epsilonm1_results>=The data type.
  unit[optional]=the unit number for output.
  prtvol[optional]=verbosity level.
  mode_paral[optional]=either COLL or PERS.

OUTPUT

  Only printing.

SOURCE

320 subroutine em1results_print(Er,unit,prtvol,mode_paral)
321 
322 !Arguments ------------------------------------
323  integer,optional,intent(in) :: unit,prtvol
324  character(len=4),optional,intent(in) :: mode_paral
325  type(Epsilonm1_results),intent(in) :: Er
326 
327 !Local variables-------------------------------
328  integer :: iw,iqibz,iqlwl,unt,my_prtvol
329  character(len=50) :: rfname,rforder,rfapprox,rftest,kxcname
330  character(len=500) :: msg
331  character(len=4) :: mode
332 ! *************************************************************************
333 
334  unt = std_out; if (present(unit)) unt = unit
335  my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol
336  mode   ='COLL'; if (present(mode_paral)) mode = mode_paral
337 
338  ! === chi0 or \epsilon^{-1} ? ===
339  SELECT CASE (Er%ID)
340  CASE (0)
341    rfname = 'Undefined'
342  CASE (1)
343    rfname = 'Irreducible Polarizability'
344  CASE (2)
345    rfname = 'Polarizability'
346  CASE (3)
347    rfname = 'Symmetrical Dielectric Matrix'
348  CASE (4)
349    rfname = 'Symmetrical Inverse Dielectric Matrix'
350  CASE DEFAULT
351    ABI_BUG(sjoin('Wrong Er%ID:',itoa(Er%ID)))
352  END SELECT
353 
354  ! For chi, \espilon or \epsilon^{-1}, define the approximation.
355  rfapprox='None'
356  if (Er%ID>=2.or.Er%ID<=4) then
357    if (Er%ikxc==0) then
358      rfapprox='RPA'
359    else if (Er%ikxc>0) then
360      rfapprox='Static TDDFT'
361    else
362      rfapprox='TDDFT'
363    end if
364  end if
365 
366  ! === If TDDFT and \epsilon^{-1}, define the type ===
367  rftest='None'
368 ! if (Er%ID==0) then
369 !  if (Er%test_type==0) then
370 !   rftest='TEST-PARTICLE'
371 !  else if (Er%test_type==1) then
372 !   rftest='TEST-ELECTRON'
373 !  else
374 !   write(msg,'(4a,i3)')ch10,&
375 !&   ' em1results_print : BUG - ',ch10,&
376 !&   ' Wrong value of Er%test_type = ',Er%test_type
377 !   ABI_ERROR(msg)
378 !  end if
379 ! end if
380 
381  ! === Define time-ordering ===
382  rforder='Undefined'
383  if (Er%Tordering==1) then
384    rforder='Time-Ordered'
385  else if (Er%Tordering==2) then
386    rforder='Advanced'
387  else if (Er%Tordering==3) then
388    rforder='Retarded'
389  else
390    ABI_BUG(sjoin('Wrong er%tordering= ',itoa(Er%Tordering)))
391  end if
392 
393  kxcname='None'
394  if (Er%ikxc/=0) then
395    !TODO Add function to retrieve kxc name
396    ABI_ERROR('Add function to retrieve kxc name')
397    kxcname='XXXXX'
398  end if
399 
400  write(msg,'(6a,5(3a))')ch10,&
401 &  ' ==== Info on the Response Function ==== ',ch10,&
402 &  '  Associated File ................  ',TRIM(Er%fname),ch10,&
403 &  '  Response Function Type .......... ',TRIM(rfname),ch10,&
404 &  '  Type of Approximation ........... ',TRIM(rfapprox),ch10,&
405 &  '  XC kernel used .................. ',TRIM(kxcname),ch10,&
406 &  '  Type of probing particle ........ ',TRIM(rftest),ch10,&
407 &  '  Time-Ordering ................... ',TRIM(rforder),ch10
408  call wrtout(unt,msg,mode)
409  write(msg,'(a,2i4,a,3(a,i4,a),a,3i4,2a,i4,a)')&
410    '  Number of components ............ ',Er%nI,Er%nJ,ch10,&
411    '  Number of q-points in the IBZ ... ',Er%nqibz,ch10,&
412    '  Number of q-points for q-->0 .... ',Er%nqlwl,ch10,&
413    '  Number of G-vectors ............. ',Er%npwe,ch10,&
414    '  Number of frequencies ........... ',Er%nomega,Er%nomega_r,Er%nomega_i,ch10,&
415    '  Value of mqmem .................. ',Er%mqmem,ch10
416  call wrtout(unt,msg,mode)
417 
418  if (Er%nqlwl/=0) then
419    write(msg,'(a,i3)')' q-points for long wavelength limit: ',Er%nqlwl
420    call wrtout(unt,msg,mode)
421    do iqlwl=1,Er%nqlwl
422      write(msg,'(1x,i5,a,3es16.8)')iqlwl,') ',Er%qlwl(:,iqlwl)
423      call wrtout(unt,msg,mode)
424    end do
425  end if
426 
427  if (my_prtvol>0) then
428    ! Print out head and wings in the long-wavelength limit.
429    ! TODO add additional stuff.
430    write(msg,'(a,i4)')' Calculated Frequencies: ',Er%nomega
431    call wrtout(unt,msg,mode)
432    do iw=1,Er%nomega
433      write(msg,'(i4,es14.6)')iw,Er%omega(iw)*Ha_eV
434      call wrtout(unt,msg,mode)
435    end do
436 
437    write(msg,'(a,i4)')' Calculated q-points: ',Er%nqibz
438    call wrtout(unt,msg,mode)
439    do iqibz=1,Er%nqibz
440      write(msg,'(1x,i4,a,3es16.8)')iqibz,') ',Er%qibz(:,iqibz)
441      call wrtout(unt,msg,mode)
442    end do
443  end if ! my_prtvol>0
444 
445 end subroutine em1results_print

m_screening/epsilonm1_results [ Types ]

[ Top ] [ m_screening ] [ Types ]

NAME

 epsilonm1_results

FUNCTION

 For the GW part of ABINIT, the epsilonm1_results structured datatype
 gather the results of screening: the inverse dielectric matrix, and the omega matrices.

SOURCE

 76  type,public :: Epsilonm1_results
 77 
 78   integer :: id
 79   ! Matrix identifier: O if not yet defined, 1 for chi0,
 80   ! 2 for chi, 3 for epsilon, 4 for espilon^{-1}, 5 for W.
 81 
 82   integer :: ikxc
 83   ! Kxc kernel used, 0 for None (RPA), >0 for static TDDFT (=ixc), <0 for TDDFT
 84 
 85   integer :: fform
 86   ! File format: 1002 for SCR|SUSC files.
 87 
 88   integer :: mqmem
 89   ! =0 for out-of-core solution, =nqibz if entire matrix is stored in memory.
 90 
 91   integer :: nI,nJ
 92   ! Number of components (rows,columns) in chi|eps^-1. (1,1) if collinear.
 93 
 94   integer :: nqibz
 95   ! Number of q-points in the IBZ used.
 96 
 97   integer :: nqlwl
 98   ! Number of point used for the treatment of the long wave-length limit.
 99 
100   integer :: nomega
101   ! Total number of frequencies.
102 
103   integer :: nomega_i
104   ! Number of purely imaginary frequencies used.
105 
106   integer :: nomega_r
107   ! Number of real frequencies used.
108 
109   integer :: npwe
110   ! Number of G vectors.
111 
112   integer :: test_type
113   ! 0 for None, 1 for TEST-PARTICLE, 2 for TEST-ELECTRON (only for TDDFT)
114 
115   integer :: tordering
116   ! 0 if not defined, 1 for Time-Ordered, 2 for Advanced, 3 for Retarded.
117 
118   character(len=fnlen) :: fname
119   ! Name of the file from which epsm1 is read.
120 
121   integer,allocatable :: gvec(:,:)
122   ! gvec(3,npwe)
123   ! G-vectors used to describe the two-point function (r.l.u.).
124 
125   real(dp),allocatable :: qibz(:,:)
126   ! qibz(3,nqibz)
127   ! q-points in reduced coordinates
128 
129   real(dp),allocatable :: qlwl(:,:)
130   ! qlwl(3,nqlwl)
131   ! q-points used for the long wave-length limit treatment.
132 
133   !real(gwp),allocatable :: epsm1_pole(:,:,:,:)
134   ! epsm1(npwe,npwe,nomega,nqibz)
135   ! Contains the two-point function $\epsilon_{G,Gp}(q,omega)$ in frequency and reciprocal space.
136 
137   complex(gwpc),allocatable :: epsm1(:,:,:,:)
138   ! epsm1(npwe,npwe,nomega,nqibz)
139   ! Contains the two-point function $\epsilon_{G,Gp}(q,omega)$ in frequency and reciprocal space.
140 
141   complex(dpc),allocatable :: omega(:)
142   ! omega(nomega)
143   ! Frequencies used both along the real and the imaginary axis.
144 
145   type(hscr_t) :: Hscr
146   ! The header reported in the _SCR of _SUSC file.
147   ! This object contains information on the susceptibility or the inverse dielectric matrix
148   ! as stored in the external file. These quantities do *NOT* correspond to the quantities
149   ! used during the GW calculation since some parameters might differ, actually they might be smaller.
150   ! For example, the number of G-vectors used can be smaller than the number of G"s stored on file.
151 
152  end type Epsilonm1_results
153 
154  public :: em1results_free                ! Free memory
155  public :: em1results_print               ! Print basic info
156  public :: Epsm1_symmetrizer              ! Symmetrize two-point function at a q-point in the BZ.
157  public :: Epsm1_symmetrizer_inplace      ! In-place version of the above
158  public :: init_Er_from_file              ! Initialize the object from file
159  public :: mkdump_Er                      ! Dump the object to a file.
160  public :: get_epsm1
161  public :: decompose_epsm1
162  public :: make_epsm1_driver              !  Calculate the inverse symmetrical dielectric matrix starting from chi0
163  public :: mkem1_q0                       ! construct the microscopic dielectric matrix for q-->0
164  public :: screen_mdielf                  ! Calculates W_{G,G'}(q,w) for a given q-point in the BZ using a model dielectric function.
165  public :: rpa_symepsm1

m_screening/Epsm1_symmetrizer [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  Epsm1_symmetrizer

FUNCTION

  Symmetrize the inverse dielectric matrix, namely calculate epsilon^{-1} at a generic
  q-point in the BZ starting from the knowledge of the matrix at a q-point in the IBZ.
  The procedure is quite generic and can be used for every two-point function which has
  the same symmetry as the crystal.

INPUTS

  nomega=Number of frequencies required. All frequencies from 1 up to nomega are symmetrized.
  npwc=Number of G vectors in symmetrized matrix, has to be smaller than Er%npwe.
  remove_exchange=If .TRUE., return e^{-1}-1 namely remove the exchange part.
  Er<Epsilonm1_results>=Data structure containing the inverse dielectric matrix.
  Gsph<gsphere_t>=data related to the G-sphere
    %grottb
    %phmSGt
  Qmesh<kmesh_t>=Structure defining the q-mesh used for Er.
    %nbz=Number of q-points in the BZ
    %tab(nbz)=Index of the symmetric q-point in the IBZ, for each point in the BZ
    %tabo(nbz)=The operation that rotates q_ibz onto \pm q_bz (depending on tabi)
    %tabi(nbz)=-1 if time-reversal has to be considered, 1 otherwise
  iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.

OUTPUT

  epsm1_qbz(npwc,npwc,nomega)=The inverse dielectric matrix at the q-point defined by iq_bz.
   Exchange part can be subtracted out.

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 \tilde\espilon^{-1}

    If q_bz = S q_ibz + G0:

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

    If time-reversal symmetry can be used then:

      $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau} \epsilon^{-1}_{-S^{-1}(G1+Go), -S^{-1}(G2+G0)}^*(q)

TODO

  Symmetrization can be skipped if iq_bz correspond to a point in the IBZ

SOURCE

500 subroutine Epsm1_symmetrizer(iq_bz,nomega,npwc,Er,Gsph,Qmesh,remove_exchange,epsm1_qbz)
501 
502 !Arguments ------------------------------------
503 !scalars
504  integer,intent(in) :: iq_bz,nomega,npwc
505  logical,intent(in) :: remove_exchange
506  type(Epsilonm1_results),intent(in) :: Er
507  type(gsphere_t),target,intent(in) :: Gsph
508  type(kmesh_t),intent(in) :: Qmesh
509 !arrays
510  complex(gwpc),intent(out) :: epsm1_qbz(npwc,npwc,nomega)
511 
512 !Local variables-------------------------------
513 !scalars
514  integer :: iw,ii,jj,iq_ibz,itim_q,isym_q,iq_loc,sg1,sg2
515  complex(gwpc) :: phmsg1t,phmsg2t_star
516 !arrays
517  real(dp) :: qbz(3)
518 
519 ! *********************************************************************
520 
521  ABI_CHECK(Er%nomega >= nomega, 'Too many frequencies required')
522  ABI_CHECK(Er%npwe >= npwc, 'Too many G-vectors required')
523 
524  ! Get iq_ibz, and symmetries from iq_ibz.
525  call Qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym_q, itim_q)
526 
527  ! If out-of-memory, only Er%espm1(:,:,:,1) has been allocated and filled.
528  iq_loc=iq_ibz; if (Er%mqmem==0) iq_loc=1
529 
530  ! MG: grottb is a 1-1 mapping, hence we can collapse the loops (false sharing is not an issue here).
531  !grottb => Gsph%rottb (1:npwc,itim_q,isym_q)
532  !phmSgt => Gsph%phmSGt(1:npwc,isym_q)
533 
534 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(sg2,sg1,phmsg1t,phmsg2t_star)
535  do iw=1,nomega
536    do jj=1,npwc
537      sg2 = Gsph%rottb(jj, itim_q, isym_q)
538      phmsg2t_star = CONJG(Gsph%phmSGt(jj, isym_q))
539      do ii=1,npwc
540        sg1 = Gsph%rottb(ii,itim_q,isym_q)
541        phmsg1t = Gsph%phmSGt(ii,isym_q)
542        epsm1_qbz(sg1,sg2,iw) = Er%epsm1(ii,jj,iw,iq_loc) * phmsg1t * phmsg2t_star
543        !epsm1_qbz(sg1,sg2,iw) = Er%epsm1(ii,jj,iw,iq_loc) * phmSgt(ii) * CONJG(phmSgt(jj))
544      end do
545    end do
546  end do
547  !
548  ! === Account for time-reversal ===
549  !epsm1_qbz(:,:,iw)=TRANSPOSE(epsm1_qbz(:,:,iw))
550  if (itim_q==2) then
551 !$OMP PARALLEL DO IF (nomega > 1)
552    do iw=1,nomega
553      call sqmat_itranspose(npwc,epsm1_qbz(:,:,iw))
554    end do
555  end if
556 
557  if (remove_exchange) then
558    ! === Subtract the exchange contribution ===
559    ! If it's a pole screening, the exchange contribution is already removed
560 !$OMP PARALLEL DO IF (nomega > 1)
561    do iw=1,nomega
562      do ii=1,npwc
563        epsm1_qbz(ii,ii,iw)=epsm1_qbz(ii,ii,iw)-CMPLX(1.0_gwp,0.0_gwp)
564      end do
565    end do
566  endif
567 
568 end subroutine Epsm1_symmetrizer

m_screening/Epsm1_symmetrizer_inplace [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  Epsm1_symmetrizer_inplace

FUNCTION

  Same function as Epsm1_symmetrizer, ecxept now the array Ep%epsm1 is modified inplace
  thorugh an auxiliary work array of dimension (npwc,npwc)

INPUTS

  nomega=Number of frequencies required. All frequencies from 1 up to nomega are symmetrized.
  npwc=Number of G vectors in symmetrized matrix, has to be smaller than Er%npwe.
  remove_exchange=If .TRUE., return e^{-1}-1 namely remove the exchange part.
  Er<Epsilonm1_results>=Data structure containing the inverse dielectric matrix.
  Gsph<gsphere_t>=data related to the G-sphere
  Er<Epsilonm1_results>=Data structure containing the inverse dielectric matrix.
  Gsph<gsphere_t>=data related to the G-sphere
    %grottb
    %phmSGt
  Qmesh<kmesh_t>=Structure defining the q-mesh used for Er.
    %nbz=Number of q-points in the BZ
    %tab(nbz)=Index of the symmetric q-point in the IBZ, for each point in the BZ
    %tabo(nbz)=The operation that rotates q_ibz onto \pm q_bz (depending on tabi)
    %tabi(nbz)=-1 if time-reversal has to be considered, 1 otherwise
  iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.

OUTPUT

  Er%epsm1(npwc,npwc,nomega,iq_loc) symmetrised

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 \tilde\espilon^{-1}
    If q_bz=Sq_ibz+G0:

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

    If time-reversal symmetry can be used then :
    $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau}\epsilon^{-1}_{-S^{-1}(G1+Go),-S^{-1}(G2+G0)}^*(q)

TODO

  Symmetrization can be skipped if iq_bz correspond to a point in the IBZ

SOURCE

620 subroutine Epsm1_symmetrizer_inplace(iq_bz,nomega,npwc,Er,Gsph,Qmesh,remove_exchange)
621 
622 !Arguments ------------------------------------
623 !scalars
624  integer,intent(in) :: iq_bz,nomega,npwc
625  logical,intent(in) :: remove_exchange
626  type(Epsilonm1_results) :: Er
627  type(gsphere_t),target,intent(in) :: Gsph
628  type(kmesh_t),intent(in) :: Qmesh
629 
630 !Local variables-------------------------------
631 !scalars
632  integer :: iw,ii,jj,iq_ibz,itim_q,isym_q,iq_loc,sg1,sg2
633 !arrays
634  real(dp) :: qbz(3)
635  complex(gwpc) :: phmsg1t,phmsg2t_star
636  complex(gwpc),allocatable :: work(:,:)
637 
638 ! *********************************************************************
639 
640  ABI_CHECK(Er%nomega>=nomega,'Too many frequencies required')
641  ABI_CHECK(Er%npwe  >=npwc , 'Too many G-vectors required')
642 
643  ABI_MALLOC(work,(npwc,npwc))
644 
645  ! * Get iq_ibz, and symmetries from iq_ibz.
646  call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q)
647 
648  ! If out-of-memory, only Er%espm1(:,:,:,1) has been allocated and filled.
649  iq_loc=iq_ibz; if (Er%mqmem==0) iq_loc=1
650 
651 !$OMP PARALLEL DO PRIVATE(sg2,sg1,phmsg1t,phmsg2t_star) IF (nomega > 1)
652  do iw=1,nomega
653    do jj=1,npwc
654      sg2 = Gsph%rottb(jj,itim_q,isym_q)
655      phmsg2t_star = CONJG(Gsph%phmSGt(jj,isym_q))
656      do ii=1,npwc
657        sg1 = Gsph%rottb(ii,itim_q,isym_q)
658        phmsg1t = Gsph%phmSGt(ii,isym_q)
659        work(sg1,sg2) = Er%epsm1(ii,jj,iw,iq_loc) * phmsg1t * phmsg2t_star
660        !work(grottb(ii),grottb(jj))=Er%epsm1(ii,jj,iw,iq_loc)*phmSgt(ii)*CONJG(phmSgt(jj))
661      end do
662    end do
663    Er%epsm1(:,:,iw,iq_loc) = work(:,:)
664  end do
665  !
666  ! === Account for time-reversal ===
667  !Er%epsm1(:,:,iw,iq_loc)=TRANSPOSE(Er%epsm1(:,:,iw,iq_loc))
668  if (itim_q==2) then
669 !$OMP PARALLEL DO IF (nomega > 1)
670    do iw=1,nomega
671      call sqmat_itranspose(npwc,Er%epsm1(:,:,iw,iq_loc))
672    end do
673  end if
674 
675  ! === Subtract the exchange contribution ===
676  if (remove_exchange) then
677 !$OMP PARALLEL DO IF (nomega > 1)
678    do iw=1,nomega
679      do ii=1,npwc
680        Er%epsm1(ii,ii,iw,iq_loc)=Er%epsm1(ii,ii,iw,iq_loc)-1.0_gwp
681      end do
682    end do
683  endif
684 
685  ABI_FREE(work)
686 
687 end subroutine Epsm1_symmetrizer_inplace

m_screening/get_epsm1 [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  get_epsm1

FUNCTION

  Work in progress but the main is idea is as follows:

  Return the symmetrized inverse dielectric matrix.
  This method implements both in-core and the out-of-core solution
  In the later, epsilon^-1 or chi0 are read from file.
  It is possible to specify options to retrieve (RPA |TDDDT, [TESTCHARGE|TESTPARTICLE]).
  All dimensions are already initialized in the Er% object, this method
  should act as a wrapper around rdscr and make_epsm1_driver. A better
  implementation will be done in the following once the coding of file handlers is completed.

INPUTS

  Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
  iqibzA[optional]=Index of the q-point to be read from file (only for out-of-memory solutions)
  iomode=option definig the file format.
  option_test
  comm=MPI communicator.

OUTPUT

  Er%epsm1

TODO

  Remove this routine. Now everything should be done with mkdump_Er

SOURCE

1108 subroutine get_epsm1(Er,Vcp,approx_type,option_test,iomode,comm,iqibzA)
1109 
1110 !Arguments ------------------------------------
1111 !scalars
1112  integer,intent(in) :: iomode,option_test,approx_type,comm
1113  integer,optional,intent(in) :: iqibzA
1114  type(vcoul_t),intent(in) :: Vcp
1115  type(Epsilonm1_results),intent(inout) :: Er
1116 
1117 !Local variables-------------------------------
1118 !scalars
1119  integer :: my_approx_type,my_option_test,ng,ierr
1120 ! *********************************************************************
1121 
1122  DBG_ENTER("COLL")
1123 
1124  my_approx_type = approx_type; my_option_test = option_test
1125 
1126  ! Vcp not yet used.
1127  ng = Vcp%ng
1128 
1129  select case (Er%mqmem)
1130  case (0)
1131    !  Out-of-core solution
1132    if (allocated(Er%epsm1))  then
1133      ABI_FREE(Er%epsm1)
1134    end if
1135    ABI_MALLOC_OR_DIE(Er%epsm1,(Er%npwe,Er%npwe,Er%nomega,1), ierr)
1136 
1137    ! FIXME there's a problem with SUSC files and MPI-IO
1138    !if (iomode == IO_MODE_MPI) then
1139    !  !write(std_out,*)"read_screening with iomode",iomode,"file: ",trim(er%fname)
1140    !  ABI_WARNING("SUSC files is buggy. Using Fortran IO")
1141    !  call read_screening(em1_ncname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm,iqiA=iqibzA)
1142    !else
1143    call read_screening(em1_ncname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm,iqiA=iqibzA)
1144    !end if
1145 
1146    if (Er%id == 4) then
1147      ! If q-slice of epsilon^-1 has been read then return
1148      !call em1results_print(Er)
1149      return
1150    else
1151      ABI_ERROR(sjoin('Wrong Er%ID', itoa(er%id)))
1152    end if
1153 
1154  case default
1155    ! In-core solution.
1156    ABI_ERROR("you should not be here")
1157  end select
1158 
1159  DBG_EXIT("COLL")
1160 
1161 end subroutine get_epsm1

m_screening/init_Er_from_file [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  init_Er_from_file

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

  fname=The name of the external file used to read the matrix.
  mqmem=0 for out-of-core solution, /=0 if entire matrix has to be stored in memory.
  npwe_asked=Number of G-vector to be used in the calculation, if <=0 use Max allowed number.
  comm=MPI communicator.

OUTPUT

  Er<Epsilonm1_results>=The structure initialized with basic dimensions and arrays.

SOURCE

711 subroutine init_Er_from_file(Er,fname,mqmem,npwe_asked,comm)
712 
713 !Arguments ------------------------------------
714  integer,intent(in) :: mqmem,npwe_asked,comm
715  character(len=*),intent(in) :: fname
716  type(Epsilonm1_results),intent(inout) :: Er
717 
718 !Local variables-------------------------------
719 !scalars
720  integer,parameter :: master=0
721  integer :: iw,fform,my_rank,unclassified
722  real(dp) :: re, im, tol
723  character(len=500) :: msg
724 
725 ! *********************************************************************
726 
727  DBG_ENTER("COLL")
728 
729  !@Epsilonm1_results
730  my_rank = xmpi_comm_rank(comm)
731 
732  ! Read header from file.
733  call wrtout(std_out,sjoin('init_Er_from_file- testing file: ',fname))
734  call Er%hscr%from_file(fname, fform, comm)
735 
736  ! Master echoes the header.
737  if (my_rank==master) call er%hscr%print()
738 
739  ! Generic Info
740  Er%ID         =0       ! Not yet initialized as epsm1 is calculated in mkdump_Er.F90
741  Er%fname      =fname
742  Er%fform      =fform
743  Er%Tordering=Er%Hscr%Tordering
744 
745 !TODO these quantitities should be checked and initiliazed in mkdump_Er
746 !BEGIN HARCODED
747  Er%nI       = 1
748  Er%nJ       = 1
749  Er%ikxc     = 0
750  Er%test_type=-1
751 
752  Er%Hscr%headform=HSCR_LATEST_HEADFORM   ! XG20090912
753 !END HARDCODED
754 
755  Er%nqibz=Er%Hscr%nqibz
756  Er%mqmem=mqmem ; if (mqmem/=0) Er%mqmem=Er%nqibz
757  ABI_MALLOC(Er%qibz,(3,Er%nqibz))
758  Er%qibz(:,:)=Er%Hscr%qibz(:,:)
759 
760  Er%nqlwl=Er%Hscr%nqlwl
761  ABI_MALLOC(Er%qlwl,(3,Er%nqlwl))
762  Er%qlwl(:,:)=Er%Hscr%qlwl(:,:)
763 
764  Er%nomega=Er%Hscr%nomega
765  ABI_MALLOC(Er%omega,(Er%nomega))
766  Er%omega(:)=Er%Hscr%omega(:)
767 
768  ! Count number of real, imaginary, and complex frequencies.
769  Er%nomega_r = 1
770  Er%nomega_i = 0
771  if (Er%nomega == 2) then
772    Er%nomega_i = 1
773  else
774    unclassified = 0
775    tol = tol6*Ha_eV
776    do iw = 2, Er%nomega
777      re =  REAL(Er%omega(iw))
778      im = AIMAG(Er%omega(iw))
779      if (re > tol .and. im < tol) then
780        Er%nomega_r = iw ! Real freqs are packed in the first locations.
781      else if (re < tol .and. im > tol) then
782        Er%nomega_i = Er%nomega_i + 1
783      else
784        unclassified = unclassified + 1
785      end if
786    end do
787    if (unclassified > 0) then
788      write(msg,'(3a,i6)')&
789 &      'Some complex frequencies are too small to qualify as real or imaginary.',ch10,&
790 &      'Number of unidentified frequencies = ', unclassified
791      ABI_WARNING(msg)
792    end if
793  end if
794 
795  ! Get G-vectors.
796  Er%npwe=Er%Hscr%npwe
797  if (npwe_asked>0) then
798    if (npwe_asked>Er%Hscr%npwe) then
799      write(msg,'(a,i8,2a,i8)')&
800 &     'Number of G-vectors saved on file is less than the value required = ',npwe_asked,ch10,&
801 &     'Calculation will proceed with Max available npwe = ',Er%Hscr%npwe
802      ABI_WARNING(msg)
803    else  ! Redefine the no. of G"s for W.
804      Er%npwe=npwe_asked
805    end if
806  end if
807 
808  ! pointer to Er%Hscr%gvec ?
809  ABI_MALLOC(Er%gvec,(3,Er%npwe))
810  Er%gvec=Er%Hscr%gvec(:,1:Er%npwe)
811 
812  DBG_EXIT("COLL")
813 
814 end subroutine init_Er_from_file

m_screening/lebedev_laikov_int [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  lebedev_laikov_int

FUNCTION

INPUTS

OUTPUT

SOURCE

2589 subroutine lebedev_laikov_int()
2590 
2591 !Arguments ------------------------------------
2592 
2593 !Local variables-------------------------------
2594 !scalars
2595  integer :: on,npts,ii,ll,mm,lmax,leb_idx !ierr,
2596  real(dp) :: accuracy
2597  complex(dpc) :: ang_int
2598 !arrays
2599  real(dp) :: cart_vpt(3) !,real_pars(0)
2600  real(dp),allocatable :: vx(:),vy(:),vz(:),ww(:)
2601  complex(dpc) :: tensor(3,3),cplx_pars(9)
2602  complex(dpc),allocatable :: ref_func(:),expd_func(:) !tmp_momenta(:)
2603 
2604 ! *************************************************************************
2605 
2606  ABI_ERROR("lebedev_laikov_int is still under development")
2607 
2608  !tensor=RESHAPE((/4.0,2.0,4.0,0.5,2.1,0.0,5.4,2.1,5.0/),(/3,3/))
2609  tensor=RESHAPE((/4.0,0.0,0.0,0.0,4.0,0.0,0.0,0.0,5.0/),(/3,3/))
2610  !tensor=RESHAPE((/1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0/),(/3,3/))
2611 
2612  npts=26
2613  ABI_MALLOC(vx,(npts))
2614  ABI_MALLOC(vy,(npts))
2615  ABI_MALLOC(vz,(npts))
2616  ABI_MALLOC(ww,(npts))
2617 
2618  !call LD0026(vx,vy,vz,ww,on)
2619 
2620  ang_int=czero
2621  do ii=1,npts
2622    cart_vpt = [vx(ii),vy(ii),vz(ii)]
2623    ang_int = ang_int + ww(ii)*DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt))
2624  end do
2625 
2626  !write(std_out,*)"quadratic form associated to tensor=",tensor
2627  write(std_out,*)"on ang_int",on,ang_int
2628 
2629  ABI_FREE(vx)
2630  ABI_FREE(vy)
2631  ABI_FREE(vz)
2632  ABI_FREE(ww)
2633 
2634  !call init_lebedev_gridset()
2635  cplx_pars = RESHAPE(tensor,(/9/)); accuracy=tol10
2636 
2637  ! This is the function to be expanded evaluated on the lebedev_laikov grid of index leb_idx
2638  leb_idx=3; npts=lebedev_npts(leb_idx)
2639  ABI_MALLOC(ref_func,(npts))
2640  do ii=1,npts
2641    !cart_vpt = Lgridset(leb_idx)%versor(:,ii)
2642    ref_func(ii) = one/DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt))
2643  end do
2644 
2645  ! Calculate the expansion in angular momenta of 1/{q.Tq}.
2646  ! Only even l-components contribute thanks to the parity of the integrand.
2647  ! tol6 seems to be an acceptable error, convergence wrt lmax is very slow even for simple tensors.
2648  ABI_MALLOC(expd_func,(npts))
2649  expd_func=czero
2650  lmax=10
2651  do ll=0,lmax,2
2652    !allocate(tmp_momenta(-ll:ll))
2653    do mm=-ll,ll
2654      ! MG: Commented becase it causes problems with the new version of abilint
2655      !call lebedev_quadrature(ylmstar_over_qTq,(/ll,mm/),real_pars,cplx_pars,ang_int,ierr,accuracy)
2656      write(std_out,*)ll,mm,ang_int
2657      !tmp_momenta(mm) = ang_int
2658      do ii=1,npts
2659        !cart_vpt = Lgridset(leb_idx)%versor(:,ii)
2660        expd_func(ii) = expd_func(ii) + four_pi*ang_int*ylmc(ll,mm,cart_vpt)
2661      end do
2662    end do
2663    !deallocate(tmp_momenta)
2664    write(std_out,*)"Error in angular expansion at l=",ll," is ",MAXVAL(ABS(expd_func-ref_func))
2665  end do
2666 
2667 !BEGINDEBUG
2668 ! do ii=1,npts
2669 !   write(777,*)ref_func(ii)
2670 !   write(778,*)expd_func(ii)
2671 ! end do
2672 !ENDDEBUG
2673 
2674  ABI_FREE(expd_func)
2675  ABI_FREE(ref_func)
2676 
2677  ABI_ERROR("Exiting from lebedev_laikov_int")
2678 
2679 end subroutine lebedev_laikov_int

m_screening/lwl_free [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 lwl_free

FUNCTION

INPUTS

OUTPUT

SOURCE

3333 subroutine lwl_free(lwl)
3334 
3335 !Arguments ------------------------------------
3336 !scalars
3337  type(lwl_t),intent(inout) :: lwl
3338 
3339 ! *************************************************************************
3340 
3341  if (allocated(lwl%head)) then
3342    ABI_FREE(lwl%head)
3343  end if
3344  if (allocated(lwl%lwing)) then
3345    ABI_FREE(lwl%lwing)
3346  end if
3347  if (allocated(lwl%uwing)) then
3348    ABI_FREE(lwl%uwing)
3349  end if
3350  if (allocated(lwl%body)) then
3351    ABI_FREE(lwl%body)
3352  end if
3353 
3354 end subroutine lwl_free

m_screening/lwl_init [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 lwl_init

FUNCTION

INPUTS

OUTPUT

SOURCE

3265 subroutine lwl_init(lwl, path, method, cryst, vcp, npwe, gvec, comm)
3266 
3267 !Arguments ------------------------------------
3268 !scalars
3269  integer,intent(in) :: comm,method,npwe
3270  character(len=*),intent(in) :: path
3271  type(crystal_t),intent(in) :: cryst
3272  type(vcoul_t),intent(in) :: vcp
3273  type(lwl_t),intent(out) :: lwl
3274 !arrays
3275  integer,intent(in) :: gvec(3,npwe)
3276 
3277 !Local variables-------------------------------
3278 !scalars
3279  integer,parameter :: master=0
3280  integer :: iomode,my_rank,nproc,unt
3281  character(len=500) :: msg
3282 !arrays
3283 
3284 ! *************************************************************************
3285 
3286  ABI_UNUSED((/cryst%natom, gvec(1,1), vcp%ng/))
3287 
3288  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
3289 
3290  lwl%fname = path
3291  lwl%method = method
3292  ABI_CHECK(any(method == [1,2,3]), sjoin("Wrong method:", itoa(method)))
3293 
3294  iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF
3295 
3296  ! Only master reads.
3297  if (my_rank == master) then
3298 
3299    select case (iomode)
3300    case (IO_MODE_FORTRAN)
3301      if (open_file(path, msg, newunit=unt, action="read", form="unformatted", status="old") /= 0) then
3302        ABI_ERROR(msg)
3303      end if
3304 
3305      close(unt)
3306 
3307    case default
3308      ABI_ERROR(sjoin("iomode:", itoa(iomode), "is not coded"))
3309    end select
3310  end if
3311 
3312  ! Broad cast data
3313  if (nproc > 1) then
3314  end if
3315 
3316 end subroutine lwl_init

m_screening/lwl_t [ Types ]

[ Top ] [ m_screening ] [ Types ]

NAME

  lwl_t

FUNCTION

SOURCE

215  type,public :: lwl_t
216 
217   integer :: npwe
218   ! Number of G vectors.
219 
220   integer :: nomega
221   ! Number of frequencies
222 
223   integer :: method
224   ! 1 = Only head
225   ! 2 = head + wings
226   ! 3 = head + wings + body corrections.
227 
228   character(len=fnlen) :: fname
229    ! Name of the file from which epsm1 is read.
230 
231    complex(dpc),allocatable :: head(:,:,:)
232    ! head(3,3,nomega)
233 
234    complex(dpc),allocatable :: lwing(:,:,:)
235    ! lwing(3,npwe,nomega)
236    ! Lower wings
237 
238    complex(dpc),allocatable :: uwing(:,:,:)
239    ! uwing(3,npwe,nomega)
240    ! Upper wings.
241 
242    complex(dpc),allocatable :: body(:,:,:)
243    ! uwing(npwe,npwe,nomega)
244    ! Body terms
245 
246  end type lwl_t
247 
248  public :: lwl_write
249  public :: lwl_init
250  !public :: lwl_from_file
251  public :: lwl_free

m_screening/lwl_write [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 lwl_write

FUNCTION

INPUTS

OUTPUT

SOURCE

3107 subroutine lwl_write(path, cryst, vcp, npwe, nomega, gvec, chi0, chi0_head, chi0_lwing, chi0_uwing, comm)
3108 
3109 !Arguments ------------------------------------
3110 !scalars
3111  integer,intent(in) :: npwe,nomega,comm
3112  character(len=*),intent(in) :: path
3113  type(crystal_t),intent(in) :: cryst
3114  type(vcoul_t),intent(in) :: Vcp
3115 !arrays
3116  integer,intent(in) :: gvec(3,npwe)
3117  complex(gwpc),intent(in) :: chi0(npwe,npwe,nomega)
3118  complex(dpc),intent(inout)  :: chi0_head(3,3,nomega),chi0_lwing(npwe,nomega,3),chi0_uwing(npwe,nomega,3)
3119 
3120 !Local variables-------------------------------
3121 !scalars
3122  integer,parameter :: master=0,prtvol=0
3123  integer :: iw,ii,iomode,unt,my_rank
3124  character(len=500) :: msg
3125  real(dp) :: length
3126  complex(dpc) :: wng(3),em1_00
3127 ! type(hscr_t),intent(out) :: hscr
3128 !arrays
3129  complex(dpc),allocatable :: wtest(:),eps_head(:,:,:)
3130 
3131 ! *************************************************************************
3132 
3133  !if (xmpi_comm_rank(comm) /= master) goto 100
3134  my_rank = xmpi_comm_rank(comm)
3135 
3136  ABI_MALLOC(wtest,(npwe))
3137 
3138  if (prtvol > 0 .and. my_rank == master) then
3139    call wrtout(std_out, "head of chi0")
3140    do iw=1,nomega
3141      call print_arr(chi0_head(:,:,iw),max_r=3,max_c=3)
3142    end do
3143 
3144    do iw=1,nomega
3145      call wrtout(std_out, "symmetrized e_00 via tensor")
3146      wng = MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT)
3147      write(std_out,*) one - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1)*Vcp%vcqlwl_sqrt(1,1)
3148 
3149      call wrtout(std_out, "symmetrized e_0G via tensor")
3150      do ii=1,npwe
3151        wng = chi0_uwing(ii,iw,:)
3152        wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1)
3153      end do
3154      call print_arr(wtest,max_r=9,unit=std_out)
3155 
3156      call wrtout(std_out, "symmetrized e_G0 via tensor")
3157      do ii=1,npwe
3158        wng = chi0_lwing(ii,iw,:)
3159        wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1)
3160      end do
3161      call print_arr(wtest,max_r=9,unit=std_out)
3162    end do
3163  end if
3164 
3165  ! Write chi0 data
3166  iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF
3167 
3168  if (my_rank == master) then
3169    if (iomode == IO_MODE_FORTRAN) then
3170      if (open_file(path,msg,newunit=unt,form="unformatted", action="write") /= 0) then
3171        ABI_ERROR(msg)
3172      end if
3173      !call hscr_io(er%hscr,fform,rdwr,unt,comm,master,iomode)
3174      do iw=1,nomega
3175        write(unt)chi0_head(:,:,iw)
3176      end do
3177      !do iw=1,nomega
3178      !  write(unt)chi0_lwing(:,iw,:)
3179      !end do
3180      !do iw=1,nomega
3181      !  write(unt)chi0_uwing(:,iw,:)
3182      !end do
3183 
3184    else
3185      ABI_ERROR(sjoin("iomode", itoa(iomode), "is not supported"))
3186    end if
3187  end if
3188 
3189  ABI_MALLOC(eps_head,(3,3,nomega))
3190  call mkem1_q0(npwe,1,1,nomega,cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm)
3191 
3192  if (my_rank == master) then
3193    if (iomode == IO_MODE_FORTRAN) then
3194      do iw=1,nomega
3195        write(unt)eps_head(:,:,iw)
3196      end do
3197      !do iw=1,nomega
3198      !  write(unt)chi0_lwing(:,iw,:)
3199      !end do
3200      !do iw=1,nomega
3201      !  write(unt)chi0_uwing(:,iw,:)
3202      !end do
3203    else
3204      ABI_ERROR(sjoin("iomode:", itoa(iomode), "is not supported"))
3205    end if
3206  end if
3207 
3208  ABI_FREE(eps_head)
3209 
3210  if (prtvol > 0 .and. my_rank == master) then
3211    length = normv(GW_Q0_DEFAULT,cryst%gmet,"G")
3212 
3213    do iw=1,nomega
3214      em1_00 = one / vdotw(GW_Q0_DEFAULT/length, MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT/length),cryst%gmet,"G")
3215      call wrtout(std_out, "e^1_{00} from tensor")
3216      write(std_out,*) em1_00
3217 
3218      call wrtout(std_out, "symmetrized e^-1_0G via tensor")
3219      do ii=1,npwe
3220        wng = chi0_uwing(ii,iw,:)
3221        wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G")
3222      end do
3223      wtest(1) = em1_00
3224      call print_arr(wtest,max_r=9,unit=std_out)
3225 
3226      call wrtout(std_out, "symmetrized e^-1_G0 via tensor")
3227      do ii=1,npwe
3228        wng = chi0_lwing(ii,iw,:)
3229        wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G")
3230      end do
3231      wtest(1) = em1_00
3232      call print_arr(wtest,max_r=9,unit=std_out)
3233    end do !iw
3234  end if
3235 
3236  ABI_FREE(wtest)
3237 
3238  if (my_rank == master) then
3239    if (iomode == IO_MODE_FORTRAN) then
3240      close(unt)
3241    else
3242      NCF_CHECK(nf90_close(unt))
3243    end if
3244  end if
3245 
3246 !100 call xmpi_barrier(comm)
3247 
3248 end subroutine lwl_write

m_screening/make_epsm1_driver [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 make_epsm1_driver

FUNCTION

  Driver routine to calculate the inverse symmetrical dielectric matrix starting
  from the irreducible polarizability. The routine considers a single q-point, and
  performs the following tasks:

  1) Calculate $\tilde\epsilon^{-1}$ using different approximations:
      * RPA
      * ALDA within TDDFT

  2) Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
     calculating these quantities for different small q-directions specified by the user
     (Not yet operative)

  3) Output the electron energy loss function and the macroscopic dielectric function with and
     without local field effects (only if non-zero real frequencies are available)

INPUTS

  iqibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
  dim_wing=Dimension of the wings (0 or 3 if q-->0)
  npwe=Number of G-vectors in chi0.
  nI,nJ=Number of rows/columns in chi0_ij (1,1 if collinear case)
  nomega=Number of frequencies.
  omega(nomega)=Frequencines in Hartree
  approx_type=Integer flag defining the type of approximation
   == 0 for RPA   ==
   == 1 for TDDFT ==
  option_test=Only for TDDFT:
   == 0 for TESTPARTICLE ==
   == 1 for TESTELECTRON ==
  Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
  nfftot=Total number of points in the FFT mesh.
  ngfft(18)=Info on the FFT mesh.
  nkxc=Dimension of the kernel in reciprocal space. 0 if kernel is not needed
  kxcg(nfftot,nkxc)=TDDFT kernel in reciprocal space on the FFT mesh. Used only if approx_type==1
  gvec(3,npwe)=G-vectors
  comm=MPI communicator.
  chi0_lwing(npwe*nI,nomega,dim_wing)=Lower wings of chi0 (only for q-->0)
  chi0_uwing(npwe*nJ,nomega,dim_wing)=Upper wings of chi0 (only for q-->0)
  chi0_head(dim_wing,dim_wing,nomega)=Head of of chi0 (only for q-->0)

OUTPUT

  spectra<spectra_t>Object containing e_macro(w) and EELS(w)

SIDE EFFECTS

  chi0(npwe*nI,npwe*nJ,nomega): in input the irreducible polarizability, in output
   the symmetrized inverse dielectric matrix.

SOURCE

1339 subroutine make_epsm1_driver(iqibz,dim_wing,npwe,nI,nJ,nomega,omega,&
1340   approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,chi0_head,&
1341   chi0_lwing,chi0_uwing,chi0,spectra,comm,&
1342   fxc_ADA,rhor) ! optional argument
1343 
1344 !Arguments ------------------------------------
1345 !scalars
1346  integer,intent(in) :: iqibz,nI,nJ,npwe,nomega,dim_wing,approx_type,option_test,nkxc,nfftot,comm
1347  real(dp),intent(in),optional :: rhor
1348  type(vcoul_t),target,intent(in) :: Vcp
1349  type(spectra_t),intent(out) :: Spectra
1350 !arrays
1351  integer,intent(in) :: ngfft(18),gvec(3,npwe)
1352  complex(gwpc),intent(in) :: kxcg(nfftot,nkxc)
1353  complex(dpc),intent(in) :: omega(nomega)
1354  complex(dpc),intent(inout) :: chi0_lwing(:,:,:)   !(npwe*nI,nomega,dim_wing)
1355  complex(dpc),intent(inout) :: chi0_uwing(:,:,:)   !(npwe*nJ,nomega,dim_wing)
1356  complex(dpc),intent(inout) :: chi0_head(:,:,:)   !(dim_wing,dim_wing,nomega)
1357  complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ,nomega)
1358  complex(gwpc),intent(in),optional :: fxc_ADA(npwe*nI,npwe*nJ)
1359 
1360 !Local variables-------------------------------
1361 !scalars
1362  integer,parameter :: master=0
1363  integer :: i1,i2,ig1,ig2,io,ierr,irank,my_nqlwl !iqlwl
1364  integer :: nor,my_rank,nprocs,comm_self,g1mg2_idx
1365  real(dp) :: ucvol
1366  logical :: is_qeq0,use_MPI
1367  character(len=500) :: msg
1368 !arrays
1369  integer :: omega_distrb(nomega)
1370  integer,allocatable :: istart(:),istop(:)
1371  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1372  real(dp),allocatable :: eelf(:,:),tmp_eelf(:)
1373  complex(dpc),allocatable :: epsm_lf(:,:),epsm_nlf(:,:),tmp_lf(:),tmp_nlf(:)
1374  complex(dpc),allocatable :: buffer_lwing(:,:),buffer_uwing(:,:)
1375  complex(gwpc),allocatable :: kxcg_mat(:,:)
1376 
1377 !bootstrap and LR
1378  integer :: istep,nstep
1379  real(dp) :: conv_err, alpha, Zr, qpg2(3), qpg2_nrm
1380  real(gwpc) :: chi00_head, fxc_head
1381  complex(gwpc),allocatable :: vfxc_boot(:,:), vfxc_boot0(:,:), vfxc_lr(:,:), vfxc_tmp(:,:), chi0_tmp(:,:), chi0_save(:,:,:)
1382  complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:)
1383 
1384 ! *************************************************************************
1385 
1386  DBG_ENTER("COLL")
1387 
1388  if (nI/=1.or.nJ/=1) then
1389    ABI_ERROR("nI or nJ=/1 not yet implemented")
1390  end if
1391 
1392  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1393 
1394  ! MG TODO We use comm_self for the inversion as the single precision version is not yet available
1395  comm_self = xmpi_comm_self
1396 
1397  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
1398 
1399  is_qeq0 = normv(Vcp%qibz(:,iqibz),gmet,'G') < GW_TOLQ0
1400 
1401  omega_distrb = my_rank
1402  use_MPI = .FALSE.
1403  use_MPI = nprocs >= nomega  ! Parallelism is not used
1404 
1405  if (use_MPI) then
1406    ! Initialize distribution table for frequencies.
1407    ABI_MALLOC(istart,(nprocs))
1408    ABI_MALLOC(istop,(nprocs))
1409    call xmpi_split_work2_i4b(nomega,nprocs,istart,istop)
1410    omega_distrb(:)=xmpi_undefined_rank
1411    do irank=0,nprocs-1
1412      i1 = istart(irank+1)
1413      i2 = istop (irank+1)
1414      if (i1<=i2) omega_distrb(i1:i2) = irank
1415    end do
1416    ABI_FREE(istart)
1417    ABI_FREE(istop)
1418  end if
1419 
1420  ! Initialize container for spectral results
1421  do nor=1,nomega
1422    if (ABS(AIMAG(omega(nor)))>1.e-3) EXIT
1423  end do
1424  nor=nor-1; if (nor==0) nor = 1 ! only imag !?
1425 
1426  if (dim_wing==3) then
1427    call wrtout(std_out,' Analyzing long wavelength limit for several q')
1428    call spectra_init(Spectra,nor,REAL(omega(1:nor)),Vcp%nqlwl,Vcp%qlwl)
1429    my_nqlwl = 1
1430    !my_nqlwl = dim_wing ! TODO
1431    !ABI_CHECK(dim_wing==SIZE(Vcp%vcqlwl_sqrt,DIM=2),"WRONG DIMS")
1432  else
1433    call spectra_init(Spectra,nor,REAL(omega(1:nor)),1,Vcp%qibz(:,iqibz))
1434    my_nqlwl = 1
1435  end if
1436  !
1437  ! NOTE: all processors have to perform this operation in order to have the
1438  !       epsm1 matrix when performing a sigma calculation starting with the file _SUS
1439  !
1440  ! Temporary arrays to store spectra.
1441  ABI_MALLOC(epsm_lf,(nomega,my_nqlwl))
1442  ABI_MALLOC(epsm_nlf,(nomega,my_nqlwl))
1443  ABI_MALLOC(eelf,(nomega,my_nqlwl))
1444  epsm_lf=czero; epsm_nlf=czero; eelf=zero
1445 
1446  ! Temporary arrays used to store output results.
1447  ABI_MALLOC(tmp_lf, (my_nqlwl))
1448  ABI_MALLOC(tmp_nlf, (my_nqlwl))
1449  ABI_MALLOC(tmp_eelf, (my_nqlwl))
1450 
1451  SELECT CASE (approx_type)
1452 
1453  CASE (0)
1454    ! RPA: \tepsilon = 1 - Vc^{1/2} chi0 Vc^{1/2}
1455    ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff.
1456    do io=1,nomega
1457      if (omega_distrb(io) == my_rank) then
1458        !write(std_out,*)"dim_wing",dim_wing
1459        call rpa_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),my_nqlwl,dim_wing, &
1460                          chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:), &
1461                          tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1462 
1463          ! Store results.
1464          epsm_lf(io,:) = tmp_lf
1465          epsm_nlf(io,:) = tmp_nlf
1466          eelf(io,:) = tmp_eelf
1467      end if
1468    end do ! nomega
1469 
1470  CASE (1)
1471    ! Vertex correction from Adiabatic TDDFT. chi_{G1,G2} = [\delta -\chi0 (vc+kxc)]^{-1}_{G1,G3} \chi0_{G3,G2}
1472    ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1473    ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1474 
1475    ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh.
1476    ABI_MALLOC_OR_DIE(kxcg_mat,(npwe,npwe), ierr)
1477 
1478    ierr=0
1479    do ig2=1,npwe
1480      do ig1=1,npwe
1481        g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft)
1482        if (g1mg2_idx>0) then
1483          kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1)
1484        else
1485          ierr=ierr+1
1486          kxcg_mat(ig1,ig2) = czero
1487        end if
1488      end do
1489    end do
1490 
1491    if (ierr/=0) then
1492      write(msg,'(a,i4,3a)')&
1493 &     'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1494 &     'Enlarge the FFT mesh to get rid of this problem. '
1495      ABI_WARNING(msg)
1496    end if
1497 
1498    !FIXME "recheck TDDFT code and parallel"
1499    ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1500    do io=1,nomega
1501      if (omega_distrb(io) == my_rank) then
1502        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),&
1503 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm)
1504 
1505        ! Store results.
1506        epsm_lf(io,:) = tmp_lf
1507        epsm_nlf(io,:) = tmp_nlf
1508        eelf(io,:) = tmp_eelf
1509      end if
1510    end do
1511 
1512    ABI_FREE(kxcg_mat)
1513    do io=1,nomega
1514      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1515      call wrtout(std_out,msg)
1516      call print_arr(chi0(:,:,io),mode_paral='PERS')
1517    end do
1518 
1519  CASE (2)
1520    ! ADA nonlocal vertex correction contained in fxc_ADA
1521    ABI_WARNING('Entered fxc_ADA branch: EXPERIMENTAL!')
1522    ! Test that argument was passed
1523    if (.not.present(fxc_ADA)) then
1524      ABI_ERROR('make_epsm1_driver was not called with optional argument fxc_ADA')
1525    end if
1526    ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1527 
1528    do io=1,nomega
1529      if (omega_distrb(io) == my_rank) then
1530        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),fxc_ADA,option_test,my_nqlwl,dim_wing,omega(io),&
1531 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm)
1532 
1533        ! Store results.
1534        epsm_lf(io,:) = tmp_lf
1535        epsm_nlf(io,:) = tmp_nlf
1536        eelf(io,:) = tmp_eelf
1537      end if
1538    end do
1539 
1540    do io=1,nomega
1541      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1542      call wrtout(std_out,msg)
1543      call print_arr(chi0(:,:,io),mode_paral='PERS')
1544    end do
1545 
1546  CASE (4)
1547    ! Bootstrap vertex kernel by Sharma [[cite:Sharma2011]]
1548    ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1549    ABI_MALLOC_OR_DIE(vfxc_boot0,(npwe*nI,npwe*nJ), ierr)
1550    ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
1551    ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1552 
1553    if (iqibz==1) then
1554      vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! Use Coulomb term for q-->0
1555    else
1556      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
1557    end if
1558 
1559    chi0_save = chi0 ! a copy of chi0 (ks)
1560    nstep = 50 ! max iteration steps
1561    alpha = 0.6 ! mixing
1562    chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1563    fxc_head = czero; vfxc_boot = czero; chi0_tmp = czero
1564    epsm_lf = czero; epsm_nlf = czero; eelf = zero
1565    write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head
1566    call wrtout(std_out,msg)
1567    ! loop
1568    conv_err = 0.1
1569    do istep=1, nstep
1570      chi0 = chi0_save
1571      io=1 ! for now only at omega=0
1572      call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),&
1573 &      chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1574      conv_err = smallest_real
1575      do ig2=1,npwe*nJ
1576        do ig1=1,npwe*nI
1577          conv_err= MAX(conv_err, ABS(chi0(ig1,ig2,1) - chi0_tmp(ig1,ig2)))
1578        end do
1579      end do
1580      write(msg,'(a,i4,a,f10.6)') ' => bootstrap itr# ', istep, ', eps^-1 max error: ', conv_err
1581      call wrtout(std_out,msg)
1582      write(msg,'(a,2f10.6)')  '    eps^-1(head):   ', chi0(1,1,1)
1583      call wrtout(std_out,msg)
1584      write(msg,'(a,2f10.6)')  '    v^-1*fxc(head): ', fxc_head
1585      call wrtout(std_out,msg)
1586      if (conv_err <= tol4) exit
1587      !
1588      chi0_tmp = chi0(:,:,1)
1589      vfxc_boot = chi0(:,:,1)/chi00_head
1590      if (istep > 1) then
1591        vfxc_boot = alpha*vfxc_boot0 + (one-alpha)*vfxc_boot
1592      end if
1593      vfxc_boot0 = vfxc_boot
1594      fxc_head = vfxc_boot(1,1)
1595      do ig1=1,npwe
1596        vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:)
1597      end do
1598    end do
1599    ! end loop
1600    if (istep <= nstep) then
1601      write(msg,'(a,i4,a)') ' => bootstrap fxc converged after ', istep, ' iterations'
1602      call wrtout(std_out,msg)
1603    else
1604      write(msg,'(a,i4,a)') ' -> bootstrap fxc not converged after ', nstep, ' iterations'
1605      ABI_WARNING(msg)
1606    end if
1607    !
1608    chi0 = chi0_save
1609    do io=1,nomega
1610      if (omega_distrb(io) == my_rank) then
1611        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1612 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1613        epsm_lf(io,:) = tmp_lf
1614        epsm_nlf(io,:) = tmp_nlf
1615        eelf(io,:) = tmp_eelf
1616      end if
1617    end do
1618 
1619    ABI_FREE(chi0_tmp)
1620    ABI_FREE(chi0_save)
1621    ABI_FREE(vfxc_boot)
1622    ABI_FREE(vfxc_boot0)
1623 
1624    do io=1,nomega
1625      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1626      call wrtout(std_out,msg)
1627      call print_arr(chi0(:,:,io),mode_paral='PERS')
1628    end do
1629 
1630  CASE(5)
1631    ! One-shot scalar bootstrap approximation
1632    ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1633    ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1634 
1635    if (iqibz==1) then
1636      vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! Use Coulomb term for q-->0
1637    else
1638      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
1639    end if
1640 
1641    chi0_save = chi0 ! a copy of chi0
1642    fxc_head = czero; vfxc_boot = czero;
1643    epsm_lf = czero; epsm_nlf = czero; eelf = zero
1644    chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1645    write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head
1646    call wrtout(std_out,msg)
1647 
1648    fxc_head = vc_sqrt(1)**2/chi00_head + vc_sqrt(1)**2/chi00_head - vc_sqrt(1)**2
1649    fxc_head = 0.5*fxc_head + 0.5*sqrt(fxc_head**2 - 4.0*vc_sqrt(1)**4/(chi00_head*chi00_head))
1650    vfxc_boot(1,1) = fxc_head
1651    write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head/vc_sqrt(1)**2
1652    call wrtout(std_out,msg)
1653 
1654    chi0 = chi0_save
1655    do io=1,nomega
1656      if (omega_distrb(io) == my_rank) then
1657        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1658 &         chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1659        epsm_lf(io,:) = tmp_lf
1660        epsm_nlf(io,:) = tmp_nlf
1661        eelf(io,:) = tmp_eelf
1662      end if
1663    end do
1664    write(msg,'(a,2f10.6)')  '    eps^-1(head):   ',chi0(1,1,1)
1665    call wrtout(std_out,msg)
1666 
1667    ABI_FREE(chi0_save)
1668    ABI_FREE(vfxc_boot)
1669 
1670    do io=1,nomega
1671      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1672      call wrtout(std_out,msg)
1673      call print_arr(chi0(:,:,io),mode_paral='PERS')
1674    end do
1675 
1676 CASE(6)
1677    ! RPA bootstrap by Rigamonti [[cite:Rigamonti2015]] and Berger [[cite:Berger2015]]
1678    ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1679    ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1680    ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
1681 
1682    if (iqibz==1) then
1683      vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! Use Coulomb term for q-->0
1684    else
1685      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
1686    end if
1687 
1688    chi0_save = chi0 ! a copy of chi0
1689    fxc_head = czero; vfxc_boot = czero;
1690    epsm_lf = czero; epsm_nlf = czero; eelf = zero
1691    !chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1692    write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head
1693    call wrtout(std_out,msg)
1694 
1695    io = 1 ! static
1696    call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),&
1697 &    chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1698    epsm_lf(1,:) = tmp_lf
1699 
1700    ! chi(RPA) = chi0 * (1 - chi0 * v_c)^-1
1701    chi0 = chi0_save
1702    do ig2=2,npwe
1703      do ig1=2,npwe
1704        chi0(ig1,ig2,io)=-vc_sqrt(ig1)*chi0(ig1,ig2,io)*vc_sqrt(ig2)
1705      end do
1706      chi0(ig2,ig2,io)=one+chi0(ig2,ig2,io)
1707    end do
1708    chi0(1,:,io) = czero; chi0(:,1,io) = czero; chi0(1,1,io) = one
1709    chi0_tmp = chi0(:,:,io)
1710    call xginv(chi0_tmp,npwe,comm=comm)
1711    chi0 = chi0_save
1712    chi0_tmp = MATMUL(chi0(:,:,io), chi0_tmp(:,:)) ! chi(RPA)
1713    do ig1=1,npwe
1714      chi0_tmp(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*chi0_tmp(ig1,:)
1715    end do
1716    !call xginv(chi0_tmp,npwe,comm=comm) ! chi(RPA)^-1
1717    !vfxc_boot = chi0_tmp/epsm_lf(1,1)
1718    !
1719    !vfxc_boot(1,1) = chi0_tmp(1,1)/epsm_lf(1,1)
1720    vfxc_boot(1,1) = one/chi0_tmp(1,1)/epsm_lf(1,1)
1721    !@WC: alternatively:
1722    !chi00_head = chi0(1,1,io)*vc_sqrt(1)**2
1723    !vfxc_boot(1,1) = one/chi00_head/epsm_lf(1,1)
1724    fxc_head = vfxc_boot(1,1)
1725    do ig1=1,npwe
1726      vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:)
1727    end do
1728    write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head
1729    call wrtout(std_out,msg)
1730 
1731    chi0 = chi0_save
1732    do io=1,nomega
1733      if (omega_distrb(io) == my_rank) then
1734        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1735 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1736        epsm_lf(io,:) = tmp_lf
1737        epsm_nlf(io,:) = tmp_nlf
1738        eelf(io,:) = tmp_eelf
1739      end if
1740    end do
1741    write(msg,'(a,2f10.6)')  '    eps^-1(head):   ',chi0(1,1,1)
1742    call wrtout(std_out,msg)
1743 
1744    ABI_FREE(chi0_save)
1745    ABI_FREE(vfxc_boot)
1746    ABI_FREE(chi0_tmp)
1747 
1748    do io=1,nomega
1749      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1750      call wrtout(std_out,msg)
1751      call print_arr(chi0(:,:,io),mode_paral='PERS')
1752    end do
1753 
1754  CASE (7)
1755    ! LR+ALDA hybrid vertex kernel by Tal
1756    ! First ALDA
1757    ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1758    ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1759    ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh.
1760    ABI_MALLOC_OR_DIE(kxcg_mat,(npwe,npwe), ierr)
1761    ierr=0
1762    do ig2=1,npwe
1763      do ig1=1,npwe
1764        g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft)
1765        if (g1mg2_idx>0) then
1766          kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1)
1767        else
1768          ierr=ierr+1
1769          kxcg_mat(ig1,ig2) = czero
1770        end if
1771      end do
1772    end do
1773    if (ierr/=0) then
1774      write(msg,'(a,i4,3a)')&
1775 &     'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1776 &     'Enlarge the FFT mesh to get rid of this problem. '
1777      ABI_WARNING(msg)
1778    end if
1779    !FIXME "recheck TDDFT code and parallel"
1780    ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1781 
1782    ! Now LR: (1-Z)*chi0^-1
1783    ABI_MALLOC_OR_DIE(vfxc_lr,(npwe*nI,npwe*nJ), ierr)
1784    ABI_MALLOC_OR_DIE(vfxc_tmp,(npwe*nI,npwe*nJ), ierr)
1785    ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
1786 
1787    if (iqibz==1) then
1788      vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! Use Coulomb term for q-->0
1789    else
1790      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
1791    end if
1792 
1793    Zr = 0.78
1794    chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1795    fxc_head = czero; vfxc_lr = czero; vfxc_tmp = czero
1796    epsm_lf = czero; epsm_nlf = czero; eelf = zero
1797    write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head
1798    call wrtout(std_out,msg)
1799    !
1800    chi0_tmp = chi0(:,:,1)
1801    call xginv(chi0_tmp,npwe,comm=comm)
1802    vfxc_lr = (one-Zr)*chi0_tmp(:,:)
1803    write(msg,'(a)') ' Constructing LR+ALDA fxc kernel'
1804    call wrtout(std_out,msg)
1805    !
1806    do ig1=1,npwe
1807       do ig2=1,npwe
1808        qpg2 = Vcp%qibz(:,iqibz) + gvec(:,ig1)
1809        qpg2_nrm = normv(qpg2,gmet,"G")
1810        qpg2 =  Vcp%qibz(:,iqibz) + gvec(:,ig2)
1811        qpg2_nrm = SQRT(qpg2_nrm * normv(qpg2,gmet,"G"))
1812        vfxc_tmp(ig1,ig2) = vfxc_lr(ig1,ig2)*exp(-(qpg2_nrm/k_thfermi(rhor))**2) + &
1813 &        kxcg_mat(ig1,ig2)*(one - exp(-(qpg2_nrm/k_thfermi(rhor))**2))
1814        !write(std_out,*) ig1, qpg2_nrm, k_thfermi(rhor), vfxc_lr(ig1,ig1), kxcg_mat(ig1,ig1), vfxc_tmp(ig1,ig1)
1815       end do
1816    end do
1817    !
1818    vfxc_lr = vfxc_tmp
1819 
1820    do io=1,nomega
1821      if (omega_distrb(io) == my_rank) then
1822        call atddft_hyb_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_lr,kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),&
1823 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1824        epsm_lf(io,:) = tmp_lf
1825        epsm_nlf(io,:) = tmp_nlf
1826        eelf(io,:) = tmp_eelf
1827      end if
1828    end do
1829 
1830    ABI_FREE(kxcg_mat)
1831    ABI_FREE(chi0_tmp)
1832    ABI_FREE(vfxc_lr)
1833    ABI_FREE(vfxc_tmp)
1834 
1835    do io=1,nomega
1836      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1837      call wrtout(std_out,msg)
1838      call print_arr(chi0(:,:,io),mode_paral='PERS')
1839    end do
1840 
1841  CASE DEFAULT
1842    ABI_BUG(sjoin('Wrong approx_type:',itoa(approx_type)))
1843  END SELECT
1844 
1845  if (use_MPI) then
1846    ! Collect results on each node.
1847    ABI_MALLOC(buffer_lwing, (size(chi0_lwing,dim=1), size(chi0_lwing, dim=3)))
1848    ABI_MALLOC(buffer_uwing, (size(chi0_uwing,dim=1), size(chi0_uwing, dim=3)))
1849 
1850    do io=1,nomega
1851      if (omega_distrb(io)/=my_rank) then
1852        ! Zero arrays.
1853        chi0(:,:,io) = czero_gw
1854        if(dim_wing>0) then
1855           chi0_lwing(:,io,:) = zero
1856           chi0_uwing(:,io,:) = zero
1857           chi0_head(:,:,io)  = czero
1858        endif
1859        epsm_lf(io,:) = czero
1860        epsm_nlf(io,:) = czero
1861        eelf(io,:) = zero
1862      end if
1863 
1864      call xmpi_sum(chi0(:,:,io), comm,ierr)
1865 
1866      if (dim_wing>0) then
1867         ! Build contiguous arrays
1868         buffer_lwing = chi0_lwing(:,io,:)
1869         buffer_uwing = chi0_uwing(:,io,:)
1870         call xmpi_sum(buffer_lwing,comm,ierr)
1871         call xmpi_sum(buffer_uwing,comm,ierr)
1872         chi0_lwing(:,io,:) = buffer_lwing
1873         chi0_uwing(:,io,:) = buffer_uwing
1874         if (size(chi0_head(:,:,io))/= zero) then
1875           call xmpi_sum(chi0_head(:,:,io),comm,ierr)
1876         end if
1877      end if
1878 
1879    end do ! iw
1880 
1881    call xmpi_sum(epsm_lf, comm,ierr )
1882    call xmpi_sum(epsm_nlf,comm,ierr)
1883    call xmpi_sum(eelf,    comm,ierr)
1884    ABI_FREE(buffer_lwing)
1885    ABI_FREE(buffer_uwing)
1886  end if
1887 
1888  ! Save results in Spectra%, mind the slicing.
1889  Spectra%emacro_nlf(:,:) = epsm_nlf(1:nor,:)
1890  Spectra%emacro_lf (:,:) = epsm_lf (1:nor,:)
1891  Spectra%eelf      (:,:) = eelf    (1:nor,:)
1892 
1893  ABI_FREE(epsm_lf)
1894  ABI_FREE(epsm_nlf)
1895  ABI_FREE(eelf)
1896  ABI_FREE(tmp_lf)
1897  ABI_FREE(tmp_nlf)
1898  ABI_FREE(tmp_eelf)
1899 
1900  DBG_EXIT("COLL")
1901 
1902 end subroutine make_epsm1_driver

m_screening/mdielf_bechstedt [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  mdielf_bechstedt

FUNCTION

  Calculates the model dielectric function for the homogeneous system
  as proposed by F. Bechstedt, in Solid State Commun. 84, 765 1992.

INPUTS

  eps_inf=Dielectric constant of the material
  qnrm=The modulus of the q-point.
  rhor=The local value of the density

SOURCE

2816 elemental function mdielf_bechstedt(eps_inf,qnrm,rhor) result(mdielf)
2817 
2818 !Arguments ------------------------------------
2819 !scalars
2820  real(dp),intent(in) :: eps_inf,qnrm,rhor
2821  real(dp) :: mdielf
2822 
2823 ! *************************************************************************
2824 
2825  mdielf = one + &
2826 &  one / ( one/(eps_inf-one) + (qnrm/k_thfermi(rhor))**2 + (three*qnrm**4)/(four*k_fermi(rhor)**2 * k_thfermi(rhor)**2) )
2827 
2828 end function mdielf_bechstedt

m_screening/mkdump_Er [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  mkdump_Er

FUNCTION

  Dump the content of an Epsilonm1_results data type on file.

INPUTS

  id_required=Identifier of the matrix to be calculated
  Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
  ngfft(18)=Info on the FFT mesh.
  nfftot=Total number of point on the FFT mesh.
  gvec(3,npwe)=Reduced coordinates of plane waves for the response functions
  npwe=Number of plane waves.
  comm=MPI communicator.

OUTPUT

SOURCE

 839 subroutine mkdump_Er(Er,Vcp,npwe,gvec,nkxc,kxcg,id_required,approx_type,&
 840 &                    ikxc_required,option_test,fname_dump,iomode,&
 841 &                    nfftot,ngfft,comm,fxc_ADA)
 842 
 843 !Arguments ------------------------------------
 844 !scalars
 845  integer,intent(in) :: id_required,approx_type,option_test,ikxc_required,nkxc
 846  integer,intent(in) :: iomode,nfftot,npwe,comm
 847  type(Epsilonm1_results),intent(inout) :: Er
 848  type(vcoul_t),intent(in) :: Vcp
 849  character(len=*),intent(in) :: fname_dump
 850 !arrays
 851  integer,intent(in) :: ngfft(18),gvec(3,npwe)
 852  complex(gwpc),intent(in) :: kxcg(nfftot,nkxc)
 853  complex(gwpc),intent(in), optional :: fxc_ADA(npwe*Er%nI,npwe*Er%nJ,Er%nqibz)
 854 
 855 !Local variables-------------------------------
 856 !scalars
 857  integer,parameter :: master=0
 858  integer :: dim_wing,iqibz,is_qeq0,mqmem_,npwe_asked
 859  integer :: unt_dump,fform,rdwr,ierr
 860  integer :: my_rank,comm_self
 861  real(dp) :: ucvol
 862  character(len=500) :: msg
 863  character(len=fnlen) :: ofname
 864  character(len=nctk_slen) :: in_varname,out_varname
 865  type(hscr_t) :: Hscr_cp
 866  type(spectra_t) :: spectra
 867 !arrays
 868  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
 869  complex(gwpc),allocatable :: epsm1(:,:,:)
 870  complex(dpc),allocatable :: dummy_lwing(:,:,:),dummy_uwing(:,:,:),dummy_head(:,:,:)
 871 
 872 ! *********************************************************************
 873 
 874  DBG_ENTER("COLL")
 875 
 876  ABI_CHECK(id_required==4,'Value of id_required not coded')
 877  ABI_CHECK(npwe==Er%npwe,"mismatch in npwe")
 878 
 879  my_rank = xmpi_comm_rank(comm); comm_self = xmpi_comm_self
 880  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
 881 
 882  ! if (Er%ID/=0) call reset_Epsilonm1(Er)
 883  Er%ID=id_required
 884 
 885  ofname = fname_dump
 886  in_varname = ncname_from_id(er%hscr%id)
 887  out_varname = ncname_from_id(id_required)
 888 
 889  write(std_out,*)'Er%ID: ',Er%ID,', Er%Hscr%ID: ',Er%Hscr%ID
 890 
 891  if (Er%ID==Er%Hscr%ID) then
 892    ! === The two-point function we are asking for is already stored on file ===
 893    ! * According to mqmem either read and store the entire matrix in memory or do nothing.
 894 
 895    if (Er%mqmem>0) then
 896      ! In-core solution.
 897      write(msg,'(a,f12.1,a)')' Memory needed for Er%epsm1 = ',two*gwpc*npwe**2*Er%nomega*Er%nqibz*b2Mb,' [Mb] <<< MEM'
 898      call wrtout(std_out,msg)
 899      ABI_MALLOC_OR_DIE(Er%epsm1,(npwe,npwe,Er%nomega,Er%nqibz), ierr)
 900 
 901      if (iomode == IO_MODE_MPI) then
 902        !call wrtout(std_out, "read_screening with MPI_IO")
 903        ABI_WARNING("SUSC files is buggy. Using Fortran IO")
 904        call read_screening(in_varname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm)
 905      else
 906        call read_screening(in_varname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm)
 907      end if
 908    else
 909      ! Out-of-core solution ===
 910      ABI_COMMENT("mqmem==0 => allocating a single q-slice of (W|chi0) (slower but less memory).")
 911      continue
 912    end if
 913 
 914    return
 915 
 916  else
 917    ! === The matrix stored on file do not correspond to the quantity required ===
 918    ! * Presently only the transformation chi0 => e^-1 is coded
 919    ! * According to Er%mqmem either calculate e^-1 dumping the result to a file
 920    !   for a subsequent use or calculate e^-1 keeping everything in memory.
 921 
 922    if (Er%mqmem==0) then
 923      ! === Open file and write the header for the SCR file ===
 924      ! * For the moment only master works.
 925 
 926      if (my_rank==master) then
 927        if (iomode == IO_MODE_ETSF) then
 928           ofname = nctk_ncify(ofname)
 929           NCF_CHECK(nctk_open_create(unt_dump, ofname, xmpi_comm_self))
 930        else
 931          if (open_file(ofname,msg,newunit=unt_dump,form="unformatted",status="unknown",action="write") /= 0) then
 932            ABI_ERROR(msg)
 933          end if
 934        end if
 935        call wrtout(std_out,sjoin('mkdump_Er: calculating and writing epsilon^-1 matrix on file: ',ofname))
 936 
 937        ! Update the entries in the header that have been modified.
 938        ! TODO, write function to return title, just for info
 939        call Er%Hscr%copy(Hscr_cp)
 940        Hscr_cp%ID = id_required
 941        Hscr_cp%ikxc = ikxc_required
 942        Hscr_cp%test_type = option_test
 943        Hscr_cp%titles(1)  = 'SCR file: epsilon^-1'
 944        Hscr_cp%titles(2)  = 'TESTPARTICLE'
 945        ! Treat the case in which a smaller matrix is used.
 946        Hscr_cp%npwe = npwe
 947 
 948        rdwr=2; fform=Hscr_cp%fform
 949        call hscr_io(hscr_cp,fform,rdwr,unt_dump,comm_self,master,iomode)
 950        call Hscr_cp%free()
 951 
 952        ABI_MALLOC_OR_DIE(epsm1,(npwe,npwe,Er%nomega), ierr)
 953 
 954        do iqibz=1,Er%nqibz
 955          is_qeq0=0
 956          if (normv(Er%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
 957          ! FIXME there's a problem with SUSC files and MPI-IO
 958          !if (iomode == IO_MODE_MPI) then
 959          !  ABI_WARNING("SUSC files is buggy. Using Fortran IO")
 960          !  call read_screening(in_varname,Er%fname,npwe,1,Er%nomega,epsm1,IO_MODE_FORTRAN,comm_self,iqiA=iqibz)
 961          !else
 962          call read_screening(in_varname,Er%fname,npwe,1,Er%nomega,epsm1,iomode,comm_self,iqiA=iqibz)
 963          !end if
 964 
 965          dim_wing=0; if (is_qeq0==1) dim_wing=3
 966          ABI_MALLOC(dummy_lwing,(npwe*Er%nI,Er%nomega,dim_wing))
 967          ABI_MALLOC(dummy_uwing,(npwe*Er%nJ,Er%nomega,dim_wing))
 968          ABI_MALLOC(dummy_head,(dim_wing,dim_wing,Er%nomega))
 969 
 970          if (approx_type<2 .or. approx_type>3) then
 971            ABI_WARNING('Entering out-of core RPA or Kxc branch')
 972            call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,&
 973 &                    approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
 974 &                    dummy_lwing,dummy_uwing,epsm1,spectra,comm_self)
 975          else
 976            ABI_WARNING('Entering out-of core fxc_ADA branch')
 977            call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,&
 978 &                    approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
 979 &                    dummy_lwing,dummy_uwing,epsm1,spectra,comm_self,fxc_ADA(:,:,iqibz))
 980          end if
 981 
 982          ABI_FREE(dummy_head)
 983          ABI_FREE(dummy_uwing)
 984          ABI_FREE(dummy_lwing)
 985 
 986          if (is_qeq0==1) then
 987            call spectra%repr(msg)
 988            call wrtout([std_out, ab_out], msg)
 989          end if
 990          call spectra%free()
 991 
 992          call write_screening(out_varname,unt_dump,iomode,npwe,Er%nomega,iqibz,epsm1)
 993        end do
 994 
 995        if (iomode == IO_MODE_ETSF) then
 996          NCF_CHECK(nf90_close(unt_dump))
 997        else
 998          close(unt_dump)
 999        endif
1000 
1001        ABI_FREE(epsm1)
1002      end if !master
1003 
1004      ! Master broadcasts ofname.
1005      ! NOTE: A synchronization is required here, else the other procs start to read the
1006      ! SCR file before it is written by the master. xmpi_bcast will synch the procs.
1007      call xmpi_bcast(ofname,  master, comm, ierr)
1008 
1009      ! Now Er% "belongs" to the file "ofname", thus
1010      ! each proc has to destroy and re-initialize the object.
1011      call em1results_free(Er)
1012 
1013      mqmem_=Er%mqmem; npwe_asked=npwe
1014      call init_Er_from_file(Er,ofname,mqmem_,npwe_asked,comm)
1015 
1016      !Now Er% has been reinitialized and ready-to-use.
1017      Er%id = id_required
1018      call em1results_print(Er)
1019    else
1020      ! ========================
1021      ! === In-core solution ===
1022      ! ========================
1023      ABI_MALLOC_OR_DIE(Er%epsm1,(npwe,npwe,Er%nomega,Er%nqibz), ierr)
1024 
1025      ! FIXME there's a problem with SUSC files and MPI-IO
1026      !if (iomode == IO_MODE_MPI) then
1027      !  !call wrtout(std_out, "read_screening with MPI_IO")
1028      !  ABI_WARNING("SUSC files is buggy. Using Fortran IO")
1029      !  call read_screening(in_varname,Er%fname,npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm)
1030      !else
1031      call read_screening(in_varname,Er%fname,npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm)
1032      !end if
1033 
1034      do iqibz=1,Er%nqibz
1035        is_qeq0=0; if (normv(Er%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1036 
1037        dim_wing=0; if (is_qeq0==1) dim_wing=3 ! FIXME
1038        ABI_MALLOC(dummy_lwing,(npwe*Er%nI,Er%nomega,dim_wing))
1039        ABI_MALLOC(dummy_uwing,(npwe*Er%nJ,Er%nomega,dim_wing))
1040        ABI_MALLOC(dummy_head,(dim_wing,dim_wing,Er%nomega))
1041 
1042        if (approx_type<2 .or. approx_type>3) then
1043          ABI_WARNING('Entering in-core RPA and Kxc branch')
1044          call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,&
1045 &                  approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1046 &                  dummy_lwing,dummy_uwing,Er%epsm1(:,:,:,iqibz),spectra,comm)
1047        else
1048          ABI_WARNING('Entering in-core fxc_ADA branch')
1049          call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,&
1050 &                  approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1051 &                  dummy_lwing,dummy_uwing,Er%epsm1(:,:,:,iqibz),spectra,comm,fxc_ADA=fxc_ADA(:,:,iqibz))
1052        end if
1053 
1054        ABI_FREE(dummy_lwing)
1055        ABI_FREE(dummy_uwing)
1056        ABI_FREE(dummy_head)
1057 
1058        if (is_qeq0==1) then
1059          call spectra%repr(msg)
1060          call wrtout([std_out, ab_out], msg)
1061        end if
1062 
1063        call spectra%free()
1064      end do
1065 
1066      Er%id = id_required
1067      call em1results_print(Er)
1068    end if
1069  end if
1070 
1071  DBG_EXIT("COLL")
1072 
1073 end subroutine mkdump_Er

m_screening/mkem1_q0 [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 mkem1_q0

FUNCTION

   This routine construct the microscopic dieletric matrix for q-->0 starting from the heads, wings and the body
   of the irreducible polarizability. Afterwards it calculates the symmetrized inverse dieletric matrix
   via a block wise inversion thus obtaining the heads and the wings of e^{-1} that can be
   used to describe the non-analytic behavior for q-->0.

INPUTS

 npwe=Number of Gs used to describe chi0
 nomega=Number of frequencies in chi0.
 n1,n2=Factors used to define the same of the chi0 matrix (1,1 if collinear, the typical case)
 Cryst<crystal_t>=Structure describing the crystal structure.
 Vcp<vcoul_t>=datatypes gathering info on the Coulomb term
 gvec(3,npwe)=G-vector for chi0 in reduced coordinates.
 comm=MPI communicator

OUTPUT

 eps_head(3,3,nomega)=The macroscopic dieletric tensor in reduced coordinates.
   The dieletric matrix along versor \hat q can be obtained with
     e(\hat q) = \hat q.eps_head \hat q if all quantities are given in Cartesian coordinates.

SIDE EFFECTS

 chi0(npwe*n1,npwe*n2,nomega)= Input: polarizability. output: inverse dieletric matrix (only the body is used)
 chi0_lwing(npwe*n1,nomega,3)
 chi0_uwing(npwe*n2,nomega,3)  Input:  the lower and upper wings of the polarizability
                               Output: the "lower" and "upper" wings of the inverse dieletric matrix. See notes below.
 chi0_head(3,3,nomega)= Input: the polarizability tensor in Cartesian coordinates.
                        Ouput: The "head" of the inverse dieletric matrix. See notes below.

NOTES

  Matrix inversion in block form.

         1  n-1
  M =  | c  u^t| 1     ==>   M^{-1} =  |  1/k          -u^t A^{-1}/k                    |
       | v  A  | n-1                   | -A^{-1} v/k    A^{-1} + (A^{-1}v u^t A^{-1})/k |

                             where k = c - u^t A^{-1} v

  Let q be a versor in reciprocal space, the symmetrized dielectric matrix with bare coulomb interaction
  can be written as

  \tilde\epsilon = | q.Eq      q.Wl(G2) |  where   E_ij = \delta_ij -4\pi chi0_head_ij
                   | q.Wl(G1)  B(G1,G2  |          Wl(G1) = -4\pi chi0_lwing(G1)
                                                   Wu(G2) = -4\pi chi0_uwing(G1)
  therefore, in Cartesian coordinates, we have:

  1) e^{-1}_{00}(q) = [ q_i q_j (E_{ij} - \sum_{GG'} Wu_i(G)a B_{GG'}^{-1} Wl_j(G')) ]^{-1} = 1/(q.Lq)

  2) e^{-1}_{0G'}(q) = -e^{-1}_{00}(q) [ \sum_{iG} q_i Wu_i(G)a B_{GG'}^{-1} ] = (q.Su) /(q.Lq)

  3) e^{-1}_{G0}(q)  = -e^{-1}_{00}(q) [ \sum_{iG'} q_i B_{GG'}^{-1} Wl_i(G') ] = (q.Sl) /(q.Lq)

  4) e^{-1}_{GG'}(q) = B_{GG'}^{-1} +
     [ \sum_{ij} q_i q_j ( \sum_T B^{-1}_{GT}^{-1} Wl_i(T)) (\sum_T' Wu_j(T') B^{-1}_{T'G'}^{-1} ] / (q.Lq)

  where Su(G,3) and Sl(G,3) are the "upper" and "lower" wings of the inverse dielectric matrix and
  L is the inverse dielectric tensor. Similar equations hold even if vectors and tensors are given in terms
  of the reciprocal lattice vectors provided that the metric tensor is taken into account.
  The main difference is in the expression for the tensor as only one metric tensor can be
  absorbed in the scalar product, the second metric multiplies one of the wings.

  *) The present implementation assumes that no cutoff technique is used in the Coulomb term.

  *) Once the tensor in know it is possible to average the quadratic form on the sphere exactly.
  In Cartesian coordinates one obtains.

    \dfrac{1}{4\pi} \int v.Tv d\Omega = Trace(T)/3

  For the inverse dielectric matrix we have to resort to a numerical integration

SOURCE

2472 subroutine mkem1_q0(npwe,n1,n2,nomega,Cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm)
2473 
2474 !Arguments ------------------------------------
2475 !scalars
2476  integer,intent(in) :: npwe,nomega,n1,n2,comm
2477  type(crystal_t),intent(in) :: Cryst
2478  type(vcoul_t),intent(in) :: Vcp
2479 !arrays
2480  integer,intent(in) :: gvec(3,npwe)
2481  complex(gwpc),intent(in) :: chi0(npwe*n1,npwe*n2,nomega)
2482  complex(dpc),intent(inout) :: chi0_lwing(npwe*n1,nomega,3)
2483  complex(dpc),intent(inout) :: chi0_uwing(npwe*n2,nomega,3)
2484  complex(dpc),intent(inout) :: chi0_head(3,3,nomega)
2485  complex(dpc),intent(out) :: eps_head(3,3,nomega)
2486 
2487 !Local variables ------------------------------
2488 !scalars
2489  integer :: iw,ig,ig1,ig2,idir,jdir
2490 !arrays
2491  real(dp),allocatable :: modg_inv(:)
2492  complex(dpc),allocatable :: eps_lwing(:,:),eps_uwing(:,:),eps_body(:,:),cvec(:)
2493 
2494 !************************************************************************
2495 
2496  ABI_CHECK(npwe/=1,"npwe must be >1")
2497  ABI_UNUSED(comm)
2498 
2499  ! Precompute 1/|G|.
2500  ABI_MALLOC(modg_inv,(npwe-1))
2501  do ig=1,npwe-1
2502    modg_inv(ig) = one/normv(gvec(:,ig+1),Cryst%gmet,'G')
2503  end do
2504 
2505  ABI_MALLOC(eps_uwing,((npwe-1)*n1,3))
2506  ABI_MALLOC(eps_lwing,((npwe-1)*n2,3))
2507  ABI_MALLOC(eps_body,(npwe-1,npwe-1))
2508  ABI_MALLOC(cvec,(npwe-1))
2509 
2510  do iw=1,nomega
2511    !
2512    ! Head and wings of the symmetrized epsilon.
2513    eps_head(:,:,iw) = -four_pi*chi0_head(:,:,iw)
2514    do idir=1,3
2515      eps_head(idir,idir,iw) = one + eps_head(idir,idir,iw)
2516      eps_lwing(:,idir) = -four_pi * modg_inv * chi0_lwing(2:,iw,idir)
2517      eps_uwing(:,idir) = -four_pi * modg_inv * chi0_uwing(2:,iw,idir)
2518      !eps_lwing(:,idir) = -chi0_lwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1)
2519      !eps_uwing(:,idir) = -chi0_uwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1)
2520    end do
2521 
2522    write(std_out,*)" espilon head"
2523    call print_arr(eps_head(:,:,iw))
2524    !
2525    ! Construct the body of the symmetrized epsilon then invert it.
2526    do ig2=1,npwe-1
2527      do ig1=1,npwe-1
2528        eps_body(ig1,ig2) = -four_pi * modg_inv(ig1)*chi0(ig1+1,ig2+1,iw )*modg_inv(ig2)
2529        !eps_body(ig1,ig2) = -Vcp%vcqlwl_sqrt(ig1+1,1)*chi0(ig1+1,ig2+1,iw)* Vcp%vcqlwl_sqrt(ig2+1,1)
2530      end do
2531      eps_body(ig2,ig2) = one + eps_body(ig2,ig2)
2532    end do
2533 
2534    call xginv(eps_body,npwe-1)
2535    !
2536    ! Overwrite chi0_head and chi0_wings with the head and the wings of the inverse dielectric matrix.
2537    do jdir=1,3
2538      !
2539      ! Head.
2540      cvec=czero
2541      do idir=1,3
2542        cvec = cvec + two_pi*Cryst%gmet(jdir,idir)*MATMUL(eps_body,eps_lwing(:,idir)) ! as we work in reciprocal coords.
2543      end do
2544      !cvec = MATMUL(eps_body,eps_lwing(:,jdir))
2545      do idir=1,3
2546        chi0_head(idir,jdir,iw) = eps_head(idir,jdir,iw) - xdotu(npwe-1,eps_uwing(:,idir),1,cvec,1)
2547      end do
2548      !
2549      ! Now the wings.
2550      chi0_uwing(2:,iw,jdir) = -MATMUL(eps_uwing(:,jdir),eps_body)
2551      chi0_lwing(2:,iw,jdir) = -MATMUL(eps_body,eps_lwing(:,jdir))
2552      !
2553    end do !jdir
2554 
2555    call wrtout(std_out, "espilon^1 head after block inversion")
2556    call print_arr(chi0_head(:,:,iw))
2557    !
2558    ! Change the body but do not add the corrections due to the head and the wings.
2559    ! since they can be obtained on the fly from eps_body and the wings of eps^{-1}.
2560    !%chi0(2:,2:,iw) = eps_body
2561  end do !iw
2562 
2563  ABI_FREE(modg_inv)
2564  ABI_FREE(cvec)
2565  ABI_FREE(eps_lwing)
2566  ABI_FREE(eps_uwing)
2567  ABI_FREE(eps_body)
2568 
2569  RETURN
2570  ABI_UNUSED(Vcp%ng)
2571 
2572 end subroutine mkem1_q0

m_screening/rpa_symepsm1 [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 rpa_symepsm1

FUNCTION

  Calculate RPA $\tilde\epsilon^{-1}$

  Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space
  calculating these quantities for different small q-directions specified by the user
  (Not yet operative)

INPUTS

  iqibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated
  Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction
  npwe=Number of G-vectors in chi0.
  nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case)
  dim_wing=Dimension of the wings (0 or 3 if q-->0)
  chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
  chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0)
  chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0)
  comm=MPI communicator.

SIDE EFFECTS

  chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability,
                         in output the symmetrized inverse dielectric matrix.

SOURCE

1935 subroutine rpa_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0,my_nqlwl,dim_wing,chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
1936 
1937 !Arguments ------------------------------------
1938 !scalars
1939  integer,intent(in) :: iqibz,nI,nJ,npwe,dim_wing,my_nqlwl,comm
1940  type(vcoul_t),target,intent(in) :: Vcp
1941 !arrays
1942  complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ)
1943  complex(dpc),intent(inout) :: chi0_lwing(:,:) !(npwe*nI,dim_wing)
1944  complex(dpc),intent(inout) :: chi0_uwing(:,:) !(npwe*nJ,dim_wing)
1945  complex(dpc),intent(inout) :: chi0_head(:,:) !(dim_wing,dim_wing)
1946  real(dp),intent(out) :: eelf(my_nqlwl)
1947  complex(dpc),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
1948 
1949 !Local variables-------------------------------
1950 !scalars
1951  integer,parameter :: master=0,prtvol=0
1952  integer :: ig1,ig2,iqlwl,my_rank,nprocs
1953  real(dp) :: ucvol
1954  logical :: is_qeq0
1955  !character(len=500) :: msg
1956 !arrays
1957  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1958  complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:)
1959  complex(gwpc),allocatable :: chi0_save(:,:)
1960 
1961 ! *************************************************************************
1962 
1963  ABI_UNUSED(chi0_head(1,1))
1964 
1965  if (nI/=1.or.nJ/=1) then
1966    ABI_ERROR("nI or nJ=/1 not yet implemented")
1967  end if
1968 
1969  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1970 
1971  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
1972 
1973  is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0)
1974  if (is_qeq0) then
1975    ABI_CHECK(iqibz==1, "q is 0 but iq_ibz /= 1")
1976  end if
1977 
1978  if (my_nqlwl>1) then
1979    ABI_MALLOC(chi0_save,(npwe*nI,npwe*nJ))
1980    chi0_save = chi0
1981  end if
1982  !
1983  ! Symmetrized RPA epsilon: 1 - Vc^{1/2} chi0 Vc^{1/2}
1984  ! vc_sqrt contains vc^{1/2}(q, G)
1985  !
1986  ! Loop over small q"s (if any) to treat the nonanalytical behavior.
1987  do iqlwl=my_nqlwl,1,-1
1988 
1989    if (my_nqlwl>1) then
1990      chi0(:,:) = chi0_save           ! restore pristine polarizability
1991      chi0(:,1) = chi0_lwing(:,iqlwl) ! change the wings
1992      chi0(1,:) = chi0_uwing(:,iqlwl)
1993    end if
1994 
1995    if (iqibz==1) then
1996      vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl)  ! Use Coulomb term for q-->0
1997    else
1998      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
1999    end if
2000 
2001    do ig2=1,npwe*nJ
2002      do ig1=1,npwe*nI
2003        chi0(ig1,ig2) = -vc_sqrt(ig1) * chi0(ig1,ig2) * vc_sqrt(ig2)
2004      end do
2005      chi0(ig2,ig2) = one + chi0(ig2,ig2)
2006    end do
2007    ! chi0, now contains \tepsilon.
2008 
2009    epsm_nlf(iqlwl)=chi0(1,1)
2010 
2011    if (prtvol > 0) then
2012      call wrtout(std_out,' Symmetrical epsilon(G,G'') ')
2013      call print_arr(chi0, unit=std_out)
2014    end if
2015    !
2016    ! === Invert tepsilon and calculate macroscopic dielectric constant ===
2017    ! * epsm_lf(w) = 1 / epsm1(G=0,Gp=0,w).
2018    ! * Since G=Gp=0, there is no difference btw symmetrical and not symmetrical.
2019 
2020    call xginv(chi0,npwe,comm=comm)
2021 
2022    epsm_lf(iqlwl) = one/chi0(1,1)
2023    eelf(iqlwl) = -AIMAG(chi0(1,1))
2024 
2025    if (prtvol > 0) then
2026      call wrtout(std_out," Symmetrical epsilon^-1(G,G'')")
2027      call print_arr(chi0, unit=std_out)
2028    end if
2029    !
2030    ! Save wings of e^-1 overwriting input values.
2031    if (dim_wing>0.and..FALSE.) then
2032      chi0_lwing(:,iqlwl) = chi0(:,1)
2033      chi0_uwing(:,iqlwl) = chi0(1,:)
2034    end if
2035 
2036  end do !iqlwl
2037 
2038  ABI_SFREE(chi0_save)
2039 
2040 end subroutine rpa_symepsm1

m_screening/screen_mdielf [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  screen_mdielf

FUNCTION

  Calculates W_{G,G'}(q,w) for a given q-point in the BZ using a model dielectric function.

INPUTS

  iq_bz=The index of the q-point in the BZ where W(q) is calculated.
  npw=Number of plane waves for W
  nomega=Number of frequency points.
  model_type=Flag defining the model.
  eps_inf=Dielectric constant of the material.
  Cryst<crystal_t>=Info on the unit cell
  Qmesh<kmesh_t>=Info on the set of q-points.
  Vcp<vcoul_t datatype>= containing information on the cutoff technique
  Gsph<Gsphere>=The G-sphere for W.
  nspden=Number of spin density components of the density.
  nfft=Number of FFT points on the dense FFT mesh
  ngfft(18)=contain all needed information about 3D FFT.
  rhor(nfft,nspden)=Electron density in real space (The PAW AE term is included)
  which= Set it to "EM1" if the symmetrized inverse dielectric matrix is wanted.
   By default the routines returns W.
  comm=MPI communicator.

OUTPUT

  w_qbz(npw,npw,nomega)

NOTES

   W_{G1,G2} =  1/2 {
     v(q+G1) \int em1(|q+G1|,r) e^{-i(G1-G2).r} dr  +
     v(q+G2) \int em1(|q+G2|,r) e^{-i(G1-G2).r} dr } / \Omega

SOURCE

2868 subroutine screen_mdielf(iq_bz,npw,nomega,model_type,eps_inf,Cryst,Qmesh,Vcp,Gsph,nspden,nfft,ngfft,rhor,which,w_qbz,comm)
2869 
2870 !Arguments ------------------------------------
2871 !scalars
2872  integer,intent(in) :: npw,nomega,nfft,nspden,iq_bz,comm,model_type
2873  real(dp),intent(in) :: eps_inf
2874  character(len=*),intent(in) :: which
2875  type(kmesh_t),intent(in) :: Qmesh
2876  type(crystal_t),intent(in) :: Cryst
2877  type(vcoul_t),target,intent(in) :: Vcp
2878  type(gsphere_t),intent(in) :: Gsph
2879 !arrays
2880  integer,intent(in) :: ngfft(18)
2881  real(dp),intent(in) :: rhor(nfft,nspden)
2882  complex(gwpc),intent(out) :: w_qbz(npw,npw,nomega)
2883 
2884 !Local variables-------------------------------
2885 !scalars
2886  integer,parameter :: tim_fourdp0=0,paral_kgb0=0,cplex1=1
2887  integer :: my_gstart,my_gstop,iq_ibz,ig,itim_q,isym_q
2888  integer :: ig1,ig2,g1mg2_fft,iw,ii,ierr,nprocs,isg,ifft !,row ,col
2889  real(dp) :: qpg2_nrm
2890  complex(dpc) :: ph_mqbzt
2891  logical :: is_qeq0,isirred
2892  !character(len=500) :: msg
2893  type(MPI_type) :: MPI_enreg_seq
2894 !arrays
2895  integer :: umklp(3)
2896  integer,allocatable :: igfft(:),g1mg2(:,:)
2897  real(dp) :: qpg2(3),qpt_bz(3)
2898  real(dp),allocatable :: em1_qpg2r(:),fofg(:,:)
2899  complex(gwpc),ABI_CONTIGUOUS pointer :: vc_sqrt_ibz(:)
2900  complex(gwpc),allocatable :: vc_qbz(:),ctmp(:,:)
2901  logical,allocatable :: mask(:)
2902 
2903 ! *************************************************************************
2904 
2905  ABI_CHECK(nomega==1,"screen_mdielf does not support nomega>1")
2906 
2907 !Fake MPI_type for the sequential part.
2908  call initmpi_seq(MPI_enreg_seq)
2909  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all')
2910 
2911  nprocs = xmpi_comm_size(comm)
2912  call xmpi_split_work(npw,comm,my_gstart,my_gstop)
2913 
2914  call qmesh%get_bz_item(iq_bz,qpt_bz,iq_ibz,isym_q,itim_q,ph_mqbzt,umklp,isirred)
2915 
2916  !if (itim_q/=1.or.isym_q/=1.or.ANY(umklp/=0) ) then
2917  !  ABI_ERROR("Bug in mdielf_bechstedt")
2918  !end if
2919  !
2920  ! Symmetrize Vc in the full BZ.
2921  is_qeq0 = (normv(qpt_bz,Cryst%gmet,'G')<GW_TOLQ0) ! Check if q==0
2922  if (is_qeq0) then
2923    vc_sqrt_ibz => Vcp%vcqlwl_sqrt(:,1)  ! Use Coulomb term for q-->0, only first Q is used, shall we average if nqlwl>1?
2924  else
2925    vc_sqrt_ibz => Vcp%vc_sqrt(:,iq_ibz)
2926  end if
2927 
2928  ABI_MALLOC(vc_qbz,(npw))
2929  do ig=1,npw
2930    isg = Gsph%rottb(ig,itim_q,isym_q)
2931    vc_qbz(isg) = vc_sqrt_ibz(ig)**2
2932  end do
2933 
2934  ABI_MALLOC(igfft,(npw))
2935  ABI_MALLOC(g1mg2,(3,npw))
2936  ABI_MALLOC(fofg,(2,nfft))
2937  ABI_MALLOC(em1_qpg2r,(nfft))
2938  ABI_MALLOC(mask,(npw))
2939 
2940  w_qbz=czero
2941  do ig2=my_gstart,my_gstop
2942    !
2943    ! Compute the index of G-G2 wave in the FFT grid.
2944    do ii=1,npw
2945      g1mg2(:,ii) = Gsph%gvec(:,ii) - Gsph%gvec(:,ig2)
2946    end do
2947    call kgindex(igfft,g1mg2,mask,MPI_enreg_seq,ngfft,npw)
2948 
2949    ! TODO can use zero-padding FFT to speed up the transform.
2950    !call sphereboundary(gbound,istwfk1,g1mg2,mgfft,npw)
2951 
2952    ! Evaluate em1_qpg2r = \int em1(|q+G2|,r) e^{-i(G1-G2).r} dr }.
2953    qpg2 = qpt_bz + Gsph%gvec(:,ig2)
2954    qpg2_nrm = normv(qpg2,Cryst%gmet,"G")
2955 
2956    do iw=1,nomega
2957      !
2958      select case (model_type)
2959      case (1)
2960        do ifft=1,nfft
2961          em1_qpg2r(ifft) = one / mdielf_bechstedt(eps_inf,qpg2_nrm,rhor(ifft,1))
2962        end do
2963      case default
2964        ABI_ERROR(sjoin("Unknown model_type:",itoa(model_type)))
2965      end select
2966 
2967      call fourdp(cplex1,fofg,em1_qpg2r,-1,MPI_enreg_seq,nfft,1,ngfft,tim_fourdp0)
2968      !
2969      ! Here, unlike the other parts of the code, the unsymmetrized e^{-1} is used.
2970      do ig1=1,npw
2971        g1mg2_fft = igfft(ig1)
2972        w_qbz(ig1,ig2,iw) = DCMPLX(fofg(1,g1mg2_fft), fofg(2,g1mg2_fft)) * vc_qbz(ig2) !/ Cryst%ucvol
2973      end do
2974    end do ! iw
2975  end do ! ig2
2976 
2977  ABI_FREE(em1_qpg2r)
2978  ABI_FREE(fofg)
2979  ABI_FREE(igfft)
2980  ABI_FREE(g1mg2)
2981  ABI_FREE(mask)
2982  !
2983  ! W = 1/2 * (A + A^H)
2984  ! The MPI sum is done inside the loop to avoid problems with the size of the packet.
2985  ABI_MALLOC_OR_DIE(ctmp,(npw,npw), ierr)
2986 
2987  do iw=1,nomega
2988    !ctmp = TRANSPOSE(CONJG(w_qbz(:,:,iw)))
2989    ctmp = GWPC_CONJG(w_qbz(:,:,iw))
2990    call sqmat_itranspose(npw,ctmp)
2991    w_qbz(:,:,iw) = half * (ctmp + w_qbz(:,:,iw))
2992    call xmpi_sum(w_qbz(:,:,iw),comm,ierr)
2993  end do
2994  !
2995  ! Calculate the symmetrized Em1. W = vc(G1)^{1/2} \tilde Em1 vc(G2)^{1/2} -------------------------
2996  if (toupper(which)=="EM1") then
2997    do ig=1,npw
2998      isg = Gsph%rottb(ig,itim_q,isym_q)
2999      vc_qbz(isg) = vc_sqrt_ibz(ig)  ! Workspace storing vc*{1/2}(q_BZ,G).
3000    end do
3001 
3002    do ig2=1,npw
3003      do ig1=1,npw
3004        ctmp(ig1,ig2) =  one / (vc_qbz(ig1) * vc_qbz(ig2))
3005      end do
3006    end do
3007 
3008    do iw=1,nomega
3009      w_qbz(:,:,iw) = w_qbz(:,:,iw) * ctmp(:,:)
3010    end do
3011  end if
3012 
3013  call destroy_mpi_enreg(MPI_enreg_seq)
3014 
3015  ABI_FREE(vc_qbz)
3016  if (allocated(ctmp)) then
3017    ABI_FREE(ctmp)
3018  end if
3019 
3020 end subroutine screen_mdielf

m_screening/ylmstar_over_qTq [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  ylmstar_over_qTq

FUNCTION

  Return Ylm(q)^*/(q.Tq) where q is a versor in Cartesian coordinates.
  and Ylm is a complex spherical Harmonics whose index (l,m) are
  passed via int_pars(1:2). T is a tensore in Cartesian coordinates
  passed via cplx_pars(1:9).

INPUTS

  cart_vers(3)=Cartesian components of the versor
  int_pars(1:2)=(l,m) indeces in Ylm. l>=0 and  m \in [-l,l]
  cplx_pars(1:9)=Tensor T in Cartesian coordinates.
  real_pars=Not used.

OUTPUT

  Value of Ylm(q)^*/(q.Tq)

SOURCE

2705 function ylmstar_over_qTq(cart_vers,int_pars,real_pars,cplx_pars)
2706 
2707 !Arguments ------------------------------------
2708 !scalars
2709  real(dp),intent(in) :: cart_vers(3)
2710  integer,intent(in) :: int_pars(:)
2711  real(dp),intent(in) :: real_pars(:)
2712  complex(dpc),intent(in) :: cplx_pars(:)
2713  complex(dpc) :: ylmstar_over_qTq
2714 !arrays
2715 
2716 !Local variables-------------------------------
2717 !scalars
2718  integer :: ll,mm
2719 !arrays
2720  complex(dpc) :: tensor(3,3)
2721 
2722 ! *************************************************************************
2723 
2724  tensor = RESHAPE(cplx_pars(1:9),(/3,3/))
2725  ll = int_pars(1) ! ll starts from zero.
2726  mm = int_pars(2) ! m \in [-l,l]
2727 
2728  ylmstar_over_qTq = CONJG(ylmc(ll,mm,cart_vers))/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers))
2729 
2730  RETURN
2731  ABI_UNUSED(real_pars(1))
2732 
2733 end function ylmstar_over_qTq

m_screening/ylmstar_wtq_over_qTq [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  ylmstar_wtq_over_qTq

FUNCTION

  Return Ylm(q)^* weight(q)/(q.Tq) where q is a versor in Cartesian coordinates.
  Ylm is a complex spherical Harmonics whose index (l,m) are
  passed via int_pars(1:2). T is a tensor in Cartesian coordinates
  passed via cplx_pars(1:9). weight(q) is the weighting function giving
  the length of the vector parallel to versor q that connects the origin
  of the lattice to one of the boundaries of the small cell centered at Gamma

INPUTS

  cart_vers(3)=Cartesian components of the versor
  int_pars(1:2)=(l,m) indeces in Ylm. l>=0 and  m \in [-l,l]
  cplx_pars(1:9)=Tensor T in Cartesian coordinates.
  real_pars(1:9)=The Cartesian vectors defining the small box centered around gamma point
    when referred to this vectors the points in the box are given by {(x,y,z) | x,y,z \in [-1,1]}.

OUTPUT

  Value of Ylm(q)^* weigh(q)/(q.Tq)

SOURCE

2762 function ylmstar_wtq_over_qTq(cart_vers,int_pars,real_pars,cplx_pars)
2763 
2764 !Arguments ------------------------------------
2765 !scalars
2766  real(dp),intent(in) :: cart_vers(3)
2767  integer,intent(in) :: int_pars(:)
2768  real(dp),intent(in) :: real_pars(:)
2769  complex(dpc),intent(in) :: cplx_pars(:)
2770  complex(dpc) :: ylmstar_wtq_over_qTq
2771 !arrays
2772 
2773 !Local variables-------------------------------
2774 !scalars
2775  integer :: ll,mm
2776  real(dp) :: wtq
2777 !arrays
2778  real(dp) :: gprimd(3,3),rprimd(3,3),red_vers(3)
2779  complex(dpc) :: tensor(3,3)
2780 
2781 ! *************************************************************************
2782 
2783  ABI_ERROR("Work in progress")
2784  ! box_len has to be tested
2785 
2786  gprimd = RESHAPE(real_pars(1:9),(/3,3/))
2787  red_vers = MATMUL(rprimd,cart_vers)
2788  wtq = box_len(red_vers,gprimd)
2789 
2790  tensor = RESHAPE(cplx_pars(1:9),(/3,3/))
2791  ll = int_pars(1) ! true ll i.e. not shifted
2792  mm = int_pars(2)
2793 
2794  ylmstar_wtq_over_qTq = CONJG(ylmc(ll,mm,cart_vers))*wtq/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers))
2795 
2796 end function ylmstar_wtq_over_qTq