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-2018 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 .

PARENTS

SOURCE

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

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.

PARENTS

      m_screening

CHILDREN

SOURCE

2168 subroutine atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,option_test,my_nqlwl,dim_wing,omega,&
2169 & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm)
2170 
2171 
2172 !This section has been created automatically by the script Abilint (TD).
2173 !Do not modify the following lines by hand.
2174 #undef ABI_FUNC
2175 #define ABI_FUNC 'atddft_symepsm1'
2176 !End of the abilint section
2177 
2178  implicit none
2179 
2180 !Arguments ------------------------------------
2181 !scalars
2182  integer,intent(in) :: iqibz,nI,nJ,npwe,dim_wing,my_nqlwl
2183  integer,intent(in) :: option_test,comm
2184  type(vcoul_t),target,intent(in) :: Vcp
2185 !arrays
2186  complex(gwpc),intent(in) :: kxcg_mat(npwe,npwe)
2187  complex(dpc),intent(in) :: omega
2188  complex(dpc),intent(inout) :: chi0_lwing(npwe*nI,dim_wing)
2189  complex(dpc),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing)
2190  complex(dpc),intent(inout) :: chi0_head(dim_wing,dim_wing)
2191  complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2192  real(dp),intent(out) :: eelf(my_nqlwl)
2193  complex(dpc),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2194 
2195 !Local variables-------------------------------
2196 !scalars
2197  integer,parameter :: master=0,prtvol=0
2198  integer :: ig1,ig2,my_rank,nprocs,ierr
2199  real(dp) :: ucvol
2200  logical :: is_qeq0
2201  character(len=500) :: msg
2202 !arrays
2203  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2204  complex(gwpc),allocatable :: chitmp(:,:)
2205  complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:)
2206 
2207 ! *************************************************************************
2208 
2209  ABI_UNUSED(chi0_head(1,1))
2210  ABI_UNUSED(chi0_lwing(1,1))
2211  ABI_UNUSED(chi0_uwing(1,1))
2212 
2213  if (nI/=1.or.nJ/=1) then
2214    MSG_ERROR("nI or nJ=/1 not yet implemented")
2215  end if
2216 
2217  ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
2218  ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded")
2219 
2220  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2221 
2222  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2223 
2224  is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0)
2225 
2226  if (iqibz==1) then
2227    !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl)  ! Use Coulomb term for q-->0
2228    vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! TODO add treatment of non-Analytic behavior
2229  else
2230    vc_sqrt => Vcp%vc_sqrt(:,iqibz)
2231  end if
2232 
2233  write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwpc*b2Mb," Mb"
2234  ABI_STAT_MALLOC(chitmp,(npwe,npwe), ierr)
2235  ABI_CHECK(ierr==0, msg)
2236  !
2237  ! * Calculate chi0*fxc.
2238  chitmp = MATMUL(chi0,kxcg_mat)
2239  ! * First calculate the NLF contribution
2240  do ig1=1,npwe
2241    do ig2=1,npwe
2242      chitmp(ig1,ig2)=-chitmp(ig1,ig2)
2243    end do
2244    chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2245  end do
2246 
2247  call xginv(chitmp,npwe,comm=comm)
2248 
2249  chitmp = MATMUL(chitmp,chi0)
2250  !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2251  chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1)
2252  chitmp(1,1)=chitmp(1,1)+one
2253 
2254  epsm_nlf(1)=chitmp(1,1)
2255 
2256  chitmp = MATMUL(chi0,kxcg_mat)
2257  ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp.
2258  do ig1=1,npwe
2259    do ig2=1,npwe
2260      chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2
2261    end do
2262    chitmp(ig1,ig1)=chitmp(ig1,ig1)+one
2263  end do
2264 
2265  ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0.
2266  call xginv(chitmp,npwe,comm=comm)
2267  chitmp=MATMUL(chitmp,chi0)
2268 
2269  ! * Save result, now chi0 contains chi.
2270  chi0=chitmp
2271 
2272  select case (option_test)
2273  case (0)
2274    ! Symmetrized TESTPARTICLE epsilon^-1
2275    call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi','COLL')
2276    do ig1=1,npwe
2277      chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)
2278      chi0(ig1,ig1)=one+chi0(ig1,ig1)
2279    end do
2280 
2281  case (1)
2282    ! Symmetrized TESTELECTRON epsilon^-1
2283    call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + (Vc + fxc)*chi',"COLL")
2284    chitmp=MATMUL(kxcg_mat,chi0)
2285 
2286    ! Perform hermitianization, only valid along the imaginary axis.
2287    if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All")
2288 
2289    do ig1=1,npwe
2290      chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:)
2291      chi0(ig1,ig1)=one+chi0(ig1,ig1)
2292    end do
2293 
2294  case default
2295    MSG_BUG(sjoin('Wrong option_test:',itoa(option_test)))
2296  end select
2297 
2298  ABI_FREE(chitmp)
2299  !
2300  ! === chi0 now contains symmetrical epsm1 ===
2301  ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w).
2302  epsm_lf(1) =  one/chi0(1,1)
2303  eelf   (1) = -AIMAG(chi0(1,1))
2304 
2305  if (prtvol > 0) then
2306    write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]'
2307    call wrtout(std_out,msg,'COLL')
2308    call print_arr(chi0,unit=std_out)
2309  end if
2310 
2311 end subroutine atddft_symepsm1

m_screening/chi_free [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 chi_free

FUNCTION

INPUTS

OUTPUT

PARENTS

      screening

CHILDREN

SOURCE

3081 subroutine chi_free(chi)
3082 
3083 
3084 !This section has been created automatically by the script Abilint (TD).
3085 !Do not modify the following lines by hand.
3086 #undef ABI_FUNC
3087 #define ABI_FUNC 'chi_free'
3088 !End of the abilint section
3089 
3090  implicit none
3091 
3092 !Arguments ------------------------------------
3093 !scalars
3094  type(chi_t),intent(inout) :: chi
3095 
3096 ! *************************************************************************
3097 
3098  if (allocated(chi%mat)) then
3099    ABI_FREE(chi%mat)
3100  end if
3101  if (allocated(chi%head)) then
3102    ABI_FREE(chi%head)
3103  end if
3104  if (allocated(chi%lwing)) then
3105    ABI_FREE(chi%lwing)
3106  end if
3107  if (allocated(chi%uwing)) then
3108    ABI_FREE(chi%uwing)
3109  end if
3110 
3111 end subroutine chi_free

m_screening/chi_new [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 chi_new

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

3035 type(chi_t) function chi_new(npwe, nomega) result(chi)
3036 
3037 
3038 !This section has been created automatically by the script Abilint (TD).
3039 !Do not modify the following lines by hand.
3040 #undef ABI_FUNC
3041 #define ABI_FUNC 'chi_new'
3042 !End of the abilint section
3043 
3044  implicit none
3045 
3046 !Arguments ------------------------------------
3047  integer,intent(in) :: npwe,nomega
3048 
3049 ! *************************************************************************
3050 
3051  chi%nomega = nomega; chi%npwe = npwe
3052 
3053  ABI_MALLOC(chi%mat, (npwe,npwe,nomega))
3054 
3055  ABI_MALLOC(chi%head, (3,3,nomega))
3056  ABI_MALLOC(chi%lwing, (npwe, nomega,3))
3057  ABI_MALLOC(chi%uwing, (npwe, nomega,3))
3058 
3059 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

183  type,public :: chi_t
184 
185    integer :: npwe
186    ! Number of G vectors.
187 
188    integer :: nomega
189    ! Number of frequencies
190 
191    complex(gwpc),allocatable :: mat(:,:,:)
192    ! mat(npwe, npwe, nomega)
193 
194    complex(dpc),allocatable :: head(:,:,:)
195    ! head(3,3,nomega)
196 
197    complex(dpc),allocatable :: lwing(:,:,:)
198    ! lwing(3,npwe,nomega)
199    ! Lower wings
200 
201    complex(dpc),allocatable :: uwing(:,:,:)
202    ! uwing(3,npwe,nomega)
203    ! Upper wings.
204 
205  end type chi_t
206 
207  public :: chi_new    ! Create new object (allocate memory)
208  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

PARENTS

      mrgscr

CHILDREN

SOURCE

1303 subroutine decompose_epsm1(Er,iqibz,eigs)
1304 
1305 
1306 !This section has been created automatically by the script Abilint (TD).
1307 !Do not modify the following lines by hand.
1308 #undef ABI_FUNC
1309 #define ABI_FUNC 'decompose_epsm1'
1310 !End of the abilint section
1311 
1312  implicit none
1313 
1314 !Arguments ------------------------------------
1315 !scalars
1316  integer,intent(in) :: iqibz
1317  type(Epsilonm1_results),intent(in) :: Er
1318 !arrays
1319  complex(dpc),intent(out) :: eigs(Er%npwe,Er%nomega)
1320 
1321 !Local variables-------------------------------
1322 !scalars
1323  integer :: info,lwork,iw,negw,ig1,ig2,idx,sdim,npwe,ierr
1324  character(len=500) :: msg
1325 !arrays
1326  real(dp),allocatable :: ww(:),rwork(:)
1327  complex(dpc),allocatable :: work(:),Adpp(:),eigvec(:,:),Afull(:,:),vs(:,:),wwc(:)
1328  logical,allocatable :: bwork(:)
1329  logical :: sortcplx !BUG in abilint
1330 
1331 ! *********************************************************************
1332 
1333  ABI_CHECK(Er%mqmem/=0,'mqmem==0 not implemented')
1334 
1335  npwe = Er%npwe
1336 
1337  do iw=1,Er%nomega
1338 
1339    if (ABS(REAL(Er%omega(iw)))>0.00001) then
1340      ! Eigenvalues for a generic complex matrix
1341      lwork=4*2*npwe
1342      ABI_MALLOC(wwc,(npwe))
1343      ABI_MALLOC(work,(lwork))
1344      ABI_MALLOC(rwork,(npwe))
1345      ABI_MALLOC(bwork,(npwe))
1346      ABI_MALLOC(vs,(npwe,npwe))
1347      ABI_MALLOC(Afull,(npwe,npwe))
1348 
1349      Afull=Er%epsm1(:,:,iw,iqibz)
1350 
1351      !for the moment no sort, maybe here I should sort using the real part?
1352      call ZGEES('V','N',sortcplx,npwe,Afull,npwe,sdim,wwc,vs,npwe,work,lwork,rwork,bwork,info)
1353      if (info/=0) then
1354        MSG_ERROR(sjoin("ZGEES returned info:",itoa(info)))
1355      end if
1356 
1357      eigs(:,iw)=wwc(:)
1358 
1359      ABI_FREE(wwc)
1360      ABI_FREE(work)
1361      ABI_FREE(rwork)
1362      ABI_FREE(bwork)
1363      ABI_FREE(vs)
1364      ABI_FREE(Afull)
1365 
1366    else
1367      ! Hermitian version.
1368      lwork=2*npwe-1
1369      ABI_MALLOC(ww,(npwe))
1370      ABI_MALLOC(work,(lwork))
1371      ABI_MALLOC(rwork,(3*npwe-2))
1372      ABI_MALLOC(eigvec,(npwe,npwe))
1373      ABI_STAT_MALLOC(Adpp,(npwe*(npwe+1)/2), ierr)
1374      ABI_CHECK(ierr==0, 'out of memory in Adpp')
1375 
1376      idx=0 ! Pack the matrix
1377      do ig2=1,npwe
1378        do ig1=1,ig2
1379          idx=idx+1
1380          Adpp(idx)=Er%epsm1(ig1,ig2,iw,iqibz)
1381        end do
1382      end do
1383 
1384      ! For the moment we require also the eigenvectors.
1385      call ZHPEV('V','U',npwe,Adpp,ww,eigvec,npwe,work,rwork,info)
1386      if (info/=0) then
1387        MSG_ERROR(sjoin('ZHPEV returned info=', itoa(info)))
1388      end if
1389 
1390      negw=(COUNT((REAL(ww)<tol6)))
1391      if (negw/=0) then
1392        write(msg,'(a,i5,a,i3,a,f8.4)')&
1393         'Found negative eigenvalues. No. ',negw,' at iqibz= ',iqibz,' minval= ',MINVAL(REAL(ww))
1394        MSG_WARNING(msg)
1395      end if
1396 
1397      eigs(:,iw)=ww(:)
1398 
1399      ABI_FREE(ww)
1400      ABI_FREE(work)
1401      ABI_FREE(rwork)
1402      ABI_FREE(eigvec)
1403      ABI_FREE(Adpp)
1404    end if
1405  end do !iw
1406 
1407 ! contains
1408 ! function sortcplx(carg) result(res)
1409 !  implicit none
1410 !  complex(dpc),intent(in) :: carg
1411 !  logical :: res
1412 !  res=.TRUE.
1413 ! end function sortcplx
1414 
1415 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

PARENTS

      m_screening,mrgscr,sigma

CHILDREN

SOURCE

283 subroutine em1results_free(Er)
284 
285 
286 !This section has been created automatically by the script Abilint (TD).
287 !Do not modify the following lines by hand.
288 #undef ABI_FUNC
289 #define ABI_FUNC 'em1results_free'
290 !End of the abilint section
291 
292  implicit none
293 
294 !Arguments ------------------------------------
295 !scalars
296  type(Epsilonm1_results),intent(inout) :: Er
297 ! *************************************************************************
298 
299  !@Epsilonm1_results
300  !integer
301  if (allocated(Er%gvec)) then
302    ABI_FREE(Er%gvec)
303  end if
304 
305  !real
306  if (allocated(Er%qibz)) then
307    ABI_FREE(Er%qibz)
308  end if
309  if (allocated(Er%qlwl)) then
310    ABI_FREE(Er%qlwl)
311  end if
312 
313  !complex
314  if (allocated(Er%epsm1)) then
315    ABI_FREE(Er%epsm1)
316  end if
317  if (allocated(Er%omega)) then
318    ABI_FREE(Er%omega)
319  end if
320 
321  !datatypes
322  call hscr_free(Er%Hscr)
323 
324 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.

PARENTS

      m_screening,mrgscr

CHILDREN

SOURCE

353 subroutine em1results_print(Er,unit,prtvol,mode_paral)
354 
355 
356 !This section has been created automatically by the script Abilint (TD).
357 !Do not modify the following lines by hand.
358 #undef ABI_FUNC
359 #define ABI_FUNC 'em1results_print'
360 !End of the abilint section
361 
362  implicit none
363 
364 !Arguments ------------------------------------
365  integer,optional,intent(in) :: unit,prtvol
366  character(len=4),optional,intent(in) :: mode_paral
367  type(Epsilonm1_results),intent(in) :: Er
368 
369 !Local variables-------------------------------
370  integer :: iw,iqibz,iqlwl,unt,my_prtvol
371  character(len=50) :: rfname,rforder,rfapprox,rftest,kxcname
372  character(len=500) :: msg
373  character(len=4) :: mode
374 ! *************************************************************************
375 
376  unt = std_out; if (present(unit)) unt = unit
377  my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol
378  mode   ='COLL'; if (present(mode_paral)) mode = mode_paral
379 
380  ! === chi0 or \epsilon^{-1} ? ===
381  SELECT CASE (Er%ID)
382  CASE (0)
383    rfname='Undefined'
384  CASE (1)
385    rfname='Irreducible Polarizability'
386  CASE (2)
387    rfname='Polarizability'
388  CASE (3)
389    rfname='Symmetrical Dielectric Matrix'
390  CASE (4)
391    rfname='Symmetrical Inverse Dielectric Matrix'
392  CASE DEFAULT
393    MSG_BUG(sjoin('Wrong Er%ID:',itoa(Er%ID)))
394  END SELECT
395 
396  ! For chi, \espilon or \epsilon^{-1}, define the approximation.
397  rfapprox='None'
398  if (Er%ID>=2.or.Er%ID<=4) then
399    if (Er%ikxc==0) then
400      rfapprox='RPA'
401    else if (Er%ikxc>0) then
402      rfapprox='Static TDDFT'
403    else
404      rfapprox='TDDFT'
405    end if
406  end if
407 
408  ! === If TDDFT and \epsilon^{-1}, define the type ===
409  rftest='None'
410 ! if (Er%ID==0) then
411 !  if (Er%test_type==0) then
412 !   rftest='TEST-PARTICLE'
413 !  else if (Er%test_type==1) then
414 !   rftest='TEST-ELECTRON'
415 !  else
416 !   write(msg,'(4a,i3)')ch10,&
417 !&   ' em1results_print : BUG - ',ch10,&
418 !&   ' Wrong value of Er%test_type = ',Er%test_type
419 !   MSG_ERROR(msg)
420 !  end if
421 ! end if
422 
423  ! === Define time-ordering ===
424  rforder='Undefined'
425  if (Er%Tordering==1) then
426    rforder='Time-Ordered'
427  else if (Er%Tordering==2) then
428    rforder='Advanced'
429  else if (Er%Tordering==3) then
430    rforder='Retarded'
431  else
432    MSG_BUG(sjoin('Wrong er%tordering= ',itoa(Er%Tordering)))
433  end if
434 
435  kxcname='None'
436  if (Er%ikxc/=0) then
437    !TODO Add function to retrieve kxc name
438    MSG_ERROR('Add function to retrieve kxc name')
439    kxcname='XXXXX'
440  end if
441 
442  write(msg,'(6a,5(3a))')ch10,&
443 &  ' ==== Info on the Response Function ==== ',ch10,&
444 &  '  Associated File ................  ',TRIM(Er%fname),ch10,&
445 &  '  Response Function Type .......... ',TRIM(rfname),ch10,&
446 &  '  Type of Approximation ........... ',TRIM(rfapprox),ch10,&
447 &  '  XC kernel used .................. ',TRIM(kxcname),ch10,&
448 &  '  Type of probing particle ........ ',TRIM(rftest),ch10,&
449 &  '  Time-Ordering ................... ',TRIM(rforder),ch10
450  call wrtout(unt,msg,mode)
451  write(msg,'(a,2i4,a,3(a,i4,a),a,3i4,2a,i4,a)')&
452 &  '  Number of components ............ ',Er%nI,Er%nJ,ch10,&
453 &  '  Number of q-points in the IBZ ... ',Er%nqibz,ch10,&
454 &  '  Number of q-points for q-->0 .... ',Er%nqlwl,ch10,&
455 &  '  Number of G-vectors ............. ',Er%npwe,ch10,&
456 &  '  Number of frequencies ........... ',Er%nomega,Er%nomega_r,Er%nomega_i,ch10,&
457 &  '  Value of mqmem .................. ',Er%mqmem,ch10
458  call wrtout(unt,msg,mode)
459 
460  if (Er%nqlwl/=0) then
461    write(msg,'(a,i3)')' q-points for long wavelength limit: ',Er%nqlwl
462    call wrtout(unt,msg,mode)
463    do iqlwl=1,Er%nqlwl
464      write(msg,'(1x,i5,a,3es16.8)')iqlwl,') ',Er%qlwl(:,iqlwl)
465      call wrtout(unt,msg,mode)
466    end do
467  end if
468 
469  if (my_prtvol>0) then ! Print out head and wings in the long-wavelength limit.
470    ! TODO add additional stuff.
471    write(msg,'(a,i4)')' Calculated Frequencies: ',Er%nomega
472    call wrtout(unt,msg,mode)
473    do iw=1,Er%nomega
474      write(msg,'(i4,es14.6)')iw,Er%omega(iw)*Ha_eV
475      call wrtout(unt,msg,mode)
476    end do
477 
478    write(msg,'(a,i4)')' Calculated q-points: ',Er%nqibz
479    call wrtout(unt,msg,mode)
480    do iqibz=1,Er%nqibz
481      write(msg,'(1x,i4,a,3es16.8)')iqibz,') ',Er%qibz(:,iqibz)
482      call wrtout(unt,msg,mode)
483    end do
484  end if ! my_prtvol>0
485 
486 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

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

PARENTS

      calc_sigc_me,cohsex_me

CHILDREN

SOURCE

544 subroutine Epsm1_symmetrizer(iq_bz,nomega,npwc,Er,Gsph,Qmesh,remove_exchange,epsm1_qbz)
545 
546 
547 !This section has been created automatically by the script Abilint (TD).
548 !Do not modify the following lines by hand.
549 #undef ABI_FUNC
550 #define ABI_FUNC 'Epsm1_symmetrizer'
551 !End of the abilint section
552 
553  implicit none
554 
555 !Arguments ------------------------------------
556 !scalars
557  integer,intent(in) :: iq_bz,nomega,npwc
558  logical,intent(in) :: remove_exchange
559  type(Epsilonm1_results),intent(in) :: Er
560  type(gsphere_t),target,intent(in) :: Gsph
561  type(kmesh_t),intent(in) :: Qmesh
562 !arrays
563  complex(gwpc),intent(out) :: epsm1_qbz(npwc,npwc,nomega)
564 
565 !Local variables-------------------------------
566 !scalars
567  integer :: iw,ii,jj,iq_ibz,itim_q,isym_q,iq_loc,sg1,sg2
568  complex(gwpc) :: phmsg1t,phmsg2t_star
569 !arrays
570  real(dp) :: qbz(3)
571 
572 ! *********************************************************************
573 
574  ABI_CHECK(Er%nomega>=nomega,'Too many frequencies required')
575  ABI_CHECK(Er%npwe  >=npwc , 'Too many G-vectors required')
576 
577  ! * Get iq_ibz, and symmetries from iq_ibz.
578  call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q)
579 
580  ! If out-of-memory, only Er%espm1(:,:,:,1) has been allocated and filled.
581  iq_loc=iq_ibz; if (Er%mqmem==0) iq_loc=1
582 
583  ! MG: grottb is a 1-1 mapping, hence we can collapse the loops (false sharing is not an issue here).
584  !grottb => Gsph%rottb (1:npwc,itim_q,isym_q)
585  !phmSgt => Gsph%phmSGt(1:npwc,isym_q)
586 
587 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(sg2,sg1,phmsg1t,phmsg2t_star)
588  do iw=1,nomega
589    do jj=1,npwc
590      sg2 = Gsph%rottb(jj,itim_q,isym_q)
591      phmsg2t_star = CONJG(Gsph%phmSGt(jj,isym_q))
592      do ii=1,npwc
593        sg1 = Gsph%rottb(ii,itim_q,isym_q)
594        phmsg1t = Gsph%phmSGt(ii,isym_q)
595        epsm1_qbz(sg1,sg2,iw) = Er%epsm1(ii,jj,iw,iq_loc) * phmsg1t * phmsg2t_star
596        !epsm1_qbz(sg1,sg2,iw) = Er%epsm1(ii,jj,iw,iq_loc) * phmSgt(ii) * CONJG(phmSgt(jj))
597      end do
598    end do
599  end do
600  !
601  ! === Account for time-reversal ===
602  !epsm1_qbz(:,:,iw)=TRANSPOSE(epsm1_qbz(:,:,iw))
603  if (itim_q==2) then
604 !$OMP PARALLEL DO IF (nomega > 1)
605    do iw=1,nomega
606      call sqmat_itranspose(npwc,epsm1_qbz(:,:,iw))
607    end do
608  end if
609 
610  if (remove_exchange) then
611    ! === Subtract the exchange contribution ===
612    ! If it's a pole screening, the exchange contribution is already removed
613 !$OMP PARALLEL DO IF (nomega > 1)
614    do iw=1,nomega
615      do ii=1,npwc
616        epsm1_qbz(ii,ii,iw)=epsm1_qbz(ii,ii,iw)-CMPLX(1.0_gwp,0.0_gwp)
617      end do
618    end do
619  endif
620 
621 end subroutine Epsm1_symmetrizer

m_screening/Epsm1_symmetrizer_inplace [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  Epsm1_symmetrizer_inplace

FUNCTION

  Same function as Epsm1_symmetrizer, ecept 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

PARENTS

      calc_sigc_me

CHILDREN

SOURCE

678 subroutine Epsm1_symmetrizer_inplace(iq_bz,nomega,npwc,Er,Gsph,Qmesh,remove_exchange)
679 
680 
681 !This section has been created automatically by the script Abilint (TD).
682 !Do not modify the following lines by hand.
683 #undef ABI_FUNC
684 #define ABI_FUNC 'Epsm1_symmetrizer_inplace'
685 !End of the abilint section
686 
687  implicit none
688 
689 !Arguments ------------------------------------
690 !scalars
691  integer,intent(in) :: iq_bz,nomega,npwc
692  logical,intent(in) :: remove_exchange
693  type(Epsilonm1_results) :: Er
694  type(gsphere_t),target,intent(in) :: Gsph
695  type(kmesh_t),intent(in) :: Qmesh
696 
697 !Local variables-------------------------------
698 !scalars
699  integer :: iw,ii,jj,iq_ibz,itim_q,isym_q,iq_loc,sg1,sg2
700 !arrays
701  real(dp) :: qbz(3)
702  complex(gwpc) :: phmsg1t,phmsg2t_star
703  complex(gwpc),allocatable :: work(:,:)
704 
705 ! *********************************************************************
706 
707  ABI_CHECK(Er%nomega>=nomega,'Too many frequencies required')
708  ABI_CHECK(Er%npwe  >=npwc , 'Too many G-vectors required')
709 
710  ABI_MALLOC(work,(npwc,npwc))
711 
712  ! * Get iq_ibz, and symmetries from iq_ibz.
713  call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q)
714 
715  ! If out-of-memory, only Er%espm1(:,:,:,1) has been allocated and filled.
716  iq_loc=iq_ibz; if (Er%mqmem==0) iq_loc=1
717 
718 !$OMP PARALLEL DO PRIVATE(sg2,sg1,phmsg1t,phmsg2t_star) IF (nomega > 1)
719  do iw=1,nomega
720    do jj=1,npwc
721      sg2 = Gsph%rottb(jj,itim_q,isym_q)
722      phmsg2t_star = CONJG(Gsph%phmSGt(jj,isym_q))
723      do ii=1,npwc
724        sg1 = Gsph%rottb(ii,itim_q,isym_q)
725        phmsg1t = Gsph%phmSGt(ii,isym_q)
726        work(sg1,sg2) = Er%epsm1(ii,jj,iw,iq_loc) * phmsg1t * phmsg2t_star
727        !work(grottb(ii),grottb(jj))=Er%epsm1(ii,jj,iw,iq_loc)*phmSgt(ii)*CONJG(phmSgt(jj))
728      end do
729    end do
730    Er%epsm1(:,:,iw,iq_loc) = work(:,:)
731  end do
732  !
733  ! === Account for time-reversal ===
734  !Er%epsm1(:,:,iw,iq_loc)=TRANSPOSE(Er%epsm1(:,:,iw,iq_loc))
735  if (itim_q==2) then
736 !$OMP PARALLEL DO IF (nomega > 1)
737    do iw=1,nomega
738      call sqmat_itranspose(npwc,Er%epsm1(:,:,iw,iq_loc))
739    end do
740  end if
741 
742  ! === Subtract the exchange contribution ===
743  if (remove_exchange) then
744 !$OMP PARALLEL DO IF (nomega > 1)
745    do iw=1,nomega
746      do ii=1,npwc
747        Er%epsm1(ii,ii,iw,iq_loc)=Er%epsm1(ii,ii,iw,iq_loc)-1.0_gwp
748      end do
749    end do
750  endif
751 
752  ABI_FREE(work)
753 
754 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

PARENTS

      calc_sigc_me,cohsex_me

CHILDREN

SOURCE

1217 subroutine get_epsm1(Er,Vcp,approx_type,option_test,iomode,comm,iqibzA)
1218 
1219 
1220 !This section has been created automatically by the script Abilint (TD).
1221 !Do not modify the following lines by hand.
1222 #undef ABI_FUNC
1223 #define ABI_FUNC 'get_epsm1'
1224 !End of the abilint section
1225 
1226  implicit none
1227 
1228 !Arguments ------------------------------------
1229 !scalars
1230  integer,intent(in) :: iomode,option_test,approx_type,comm
1231  integer,optional,intent(in) :: iqibzA
1232  type(vcoul_t),intent(in) :: Vcp
1233  type(Epsilonm1_results),intent(inout) :: Er
1234 
1235 !Local variables-------------------------------
1236 !scalars
1237  integer :: my_approx_type,my_option_test,ng,ierr
1238 ! *********************************************************************
1239 
1240  DBG_ENTER("COLL")
1241 
1242  my_approx_type = approx_type; my_option_test = option_test
1243 
1244  ! Vcp not yet used.
1245  ng = Vcp%ng
1246 
1247  select case (Er%mqmem)
1248  case (0)
1249    !  Out-of-core solution
1250    if (allocated(Er%epsm1))  then
1251      ABI_FREE(Er%epsm1)
1252    end if
1253    ABI_STAT_MALLOC(Er%epsm1,(Er%npwe,Er%npwe,Er%nomega,1), ierr)
1254    ABI_CHECK(ierr==0, 'Out-of-memory in Er%epsm1 (out-of-core)')
1255 
1256    ! FIXME there's a problem with SUSC files and MPI-IO
1257    !if (iomode == IO_MODE_MPI) then
1258    !  !write(std_out,*)"read_screening with iomode",iomode,"file: ",trim(er%fname)
1259    !  MSG_WARNING("SUSC files is buggy. Using Fortran IO")
1260    !  call read_screening(em1_ncname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm,iqiA=iqibzA)
1261    !else
1262    call read_screening(em1_ncname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm,iqiA=iqibzA)
1263    !end if
1264 
1265    if (Er%id == 4) then
1266      ! If q-slice of epsilon^-1 has been read then return
1267      !call em1results_print(Er)
1268      return
1269    else
1270      MSG_ERROR(sjoin('Wrong Er%ID', itoa(er%id)))
1271    end if
1272 
1273  case default
1274    ! In-core solution.
1275    MSG_ERROR("you should not be here")
1276  end select
1277 
1278  DBG_EXIT("COLL")
1279 
1280 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.

PARENTS

      m_screening,mrgscr,setup_sigma

CHILDREN

SOURCE

783 subroutine init_Er_from_file(Er,fname,mqmem,npwe_asked,comm)
784 
785 
786 !This section has been created automatically by the script Abilint (TD).
787 !Do not modify the following lines by hand.
788 #undef ABI_FUNC
789 #define ABI_FUNC 'init_Er_from_file'
790 !End of the abilint section
791 
792  implicit none
793 
794 !Arguments ------------------------------------
795  integer,intent(in) :: mqmem,npwe_asked,comm
796  character(len=*),intent(in) :: fname
797  type(Epsilonm1_results),intent(inout) :: Er
798 
799 !Local variables-------------------------------
800 !scalars
801  integer,parameter :: master=0
802  integer :: iw,fform,my_rank,unclassified
803  real(dp) :: re, im, tol
804  character(len=500) :: msg
805 
806 ! *********************************************************************
807 
808  DBG_ENTER("COLL")
809 
810  !@Epsilonm1_results
811  my_rank = xmpi_comm_rank(comm)
812 
813  ! Read header from file.
814  call wrtout(std_out,sjoin('init_Er_from_file- testing file: ',fname),'COLL')
815  call hscr_from_file(Er%hscr, fname, fform, comm)
816 
817  ! Master echoes the header.
818  if (my_rank==master) call hscr_print(er%hscr)
819 
820  ! Generic Info
821  Er%ID         =0       ! Not yet initialized as epsm1 is calculated in mkdump_Er.F90
822  Er%fname      =fname
823  Er%fform      =fform
824  Er%Tordering=Er%Hscr%Tordering
825 
826 !TODO these quantitities should be checked and initiliazed in mkdump_Er
827 !BEGIN HARCODED
828  Er%nI       = 1
829  Er%nJ       = 1
830  Er%ikxc     = 0
831  Er%test_type=-1
832 
833  Er%Hscr%headform=HSCR_LATEST_HEADFORM   ! XG20090912
834 !END HARDCODED
835 
836  Er%nqibz=Er%Hscr%nqibz
837  Er%mqmem=mqmem ; if (mqmem/=0) Er%mqmem=Er%nqibz
838  ABI_MALLOC(Er%qibz,(3,Er%nqibz))
839  Er%qibz(:,:)=Er%Hscr%qibz(:,:)
840 
841  Er%nqlwl=Er%Hscr%nqlwl
842  ABI_MALLOC(Er%qlwl,(3,Er%nqlwl))
843  Er%qlwl(:,:)=Er%Hscr%qlwl(:,:)
844 
845  Er%nomega=Er%Hscr%nomega
846  ABI_MALLOC(Er%omega,(Er%nomega))
847  Er%omega(:)=Er%Hscr%omega(:)
848 
849  ! Count number of real, imaginary, and complex frequencies.
850  Er%nomega_r = 1
851  Er%nomega_i = 0
852  if (Er%nomega == 2) then
853    Er%nomega_i = 1
854  else
855    unclassified = 0
856    tol = tol6*Ha_eV
857    do iw = 2, Er%nomega
858      re =  REAL(Er%omega(iw))
859      im = AIMAG(Er%omega(iw))
860      if (re > tol .and. im < tol) then
861        Er%nomega_r = iw ! Real freqs are packed in the first locations.
862      else if (re < tol .and. im > tol) then
863        Er%nomega_i = Er%nomega_i + 1
864      else
865        unclassified = unclassified + 1
866      end if
867    end do
868    if (unclassified > 0) then
869      write(msg,'(3a,i6)')&
870 &      'Some complex frequencies are too small to qualify as real or imaginary.',ch10,&
871 &      'Number of unidentified frequencies = ', unclassified
872      MSG_WARNING(msg)
873    end if
874  end if
875 
876  ! Get G-vectors.
877  Er%npwe=Er%Hscr%npwe
878  if (npwe_asked>0) then
879    if (npwe_asked>Er%Hscr%npwe) then
880      write(msg,'(a,i8,2a,i8)')&
881 &     'Number of G-vectors saved on file is less than the value required = ',npwe_asked,ch10,&
882 &     'Calculation will proceed with Max available npwe = ',Er%Hscr%npwe
883      MSG_WARNING(msg)
884    else  ! Redefine the no. of G"s for W.
885      Er%npwe=npwe_asked
886    end if
887  end if
888 
889  ! pointer to Er%Hscr%gvec ?
890  ABI_MALLOC(Er%gvec,(3,Er%npwe))
891  Er%gvec=Er%Hscr%gvec(:,1:Er%npwe)
892 
893  DBG_EXIT("COLL")
894 
895 end subroutine init_Er_from_file

m_screening/lebedev_laikov_int [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

  lebedev_laikov_int

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

2526 subroutine lebedev_laikov_int()
2527 
2528 
2529 !This section has been created automatically by the script Abilint (TD).
2530 !Do not modify the following lines by hand.
2531 #undef ABI_FUNC
2532 #define ABI_FUNC 'lebedev_laikov_int'
2533 !End of the abilint section
2534 
2535  implicit none
2536 
2537 !Arguments ------------------------------------
2538 
2539 !Local variables-------------------------------
2540 !scalars
2541  integer :: on,npts,ii,ll,mm,lmax,leb_idx !ierr,
2542  real(dp) :: accuracy
2543  complex(dpc) :: ang_int
2544 !arrays
2545  real(dp) :: cart_vpt(3) !,real_pars(0)
2546  real(dp),allocatable :: vx(:),vy(:),vz(:),ww(:)
2547  complex(dpc) :: tensor(3,3),cplx_pars(9)
2548  complex(dpc),allocatable :: ref_func(:),expd_func(:) !tmp_momenta(:)
2549 
2550 ! *************************************************************************
2551 
2552  MSG_ERROR("lebedev_laikov_int is still under development")
2553 
2554  !tensor=RESHAPE((/4.0,2.0,4.0,0.5,2.1,0.0,5.4,2.1,5.0/),(/3,3/))
2555  tensor=RESHAPE((/4.0,0.0,0.0,0.0,4.0,0.0,0.0,0.0,5.0/),(/3,3/))
2556  !tensor=RESHAPE((/1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0/),(/3,3/))
2557 
2558  npts=26
2559  ABI_MALLOC(vx,(npts))
2560  ABI_MALLOC(vy,(npts))
2561  ABI_MALLOC(vz,(npts))
2562  ABI_MALLOC(ww,(npts))
2563 
2564  !call LD0026(vx,vy,vz,ww,on)
2565 
2566  ang_int=czero
2567  do ii=1,npts
2568    cart_vpt = [vx(ii),vy(ii),vz(ii)]
2569    ang_int = ang_int + ww(ii)*DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt))
2570  end do
2571 
2572  !write(std_out,*)"quadratic form associated to tensor=",tensor
2573  write(std_out,*)"on ang_int",on,ang_int
2574 
2575  ABI_FREE(vx)
2576  ABI_FREE(vy)
2577  ABI_FREE(vz)
2578  ABI_FREE(ww)
2579 
2580  !call init_lebedev_gridset()
2581  cplx_pars = RESHAPE(tensor,(/9/)); accuracy=tol10
2582 
2583  ! This is the function to be expanded evaluated on the lebedev_laikov grid of index leb_idx
2584  leb_idx=3; npts=lebedev_npts(leb_idx)
2585  ABI_MALLOC(ref_func,(npts))
2586  do ii=1,npts
2587    !cart_vpt = Lgridset(leb_idx)%versor(:,ii)
2588    ref_func(ii) = one/DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt))
2589  end do
2590 
2591  ! Calculate the expansion in angular momenta of 1/{q.Tq}.
2592  ! Only even l-components contribute thanks to the parity of the integrand.
2593  ! tol6 seems to be an acceptable error, convergence wrt lmax is very slow even for simple tensors.
2594  ABI_MALLOC(expd_func,(npts))
2595  expd_func=czero
2596  lmax=10
2597  do ll=0,lmax,2
2598    !allocate(tmp_momenta(-ll:ll))
2599    do mm=-ll,ll
2600      ! MG: Commented becase it causes problems with the new version of abilint
2601      !call lebedev_quadrature(ylmstar_over_qTq,(/ll,mm/),real_pars,cplx_pars,ang_int,ierr,accuracy)
2602      write(std_out,*)ll,mm,ang_int
2603      !tmp_momenta(mm) = ang_int
2604      do ii=1,npts
2605        !cart_vpt = Lgridset(leb_idx)%versor(:,ii)
2606        expd_func(ii) = expd_func(ii) + four_pi*ang_int*ylmc(ll,mm,cart_vpt)
2607      end do
2608    end do
2609    !deallocate(tmp_momenta)
2610    write(std_out,*)"Error in angular expansion at l=",ll," is ",MAXVAL(ABS(expd_func-ref_func))
2611  end do
2612 
2613 !BEGINDEBUG
2614 ! do ii=1,npts
2615 !   write(777,*)ref_func(ii)
2616 !   write(778,*)expd_func(ii)
2617 ! end do
2618 !ENDDEBUG
2619 
2620  ABI_FREE(expd_func)
2621  ABI_FREE(ref_func)
2622 
2623  MSG_ERROR("Exiting from lebedev_laikov_int")
2624 
2625 end subroutine lebedev_laikov_int

m_screening/lwl_free [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 lwl_free

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

3387 subroutine lwl_free(lwl)
3388 
3389 
3390 !This section has been created automatically by the script Abilint (TD).
3391 !Do not modify the following lines by hand.
3392 #undef ABI_FUNC
3393 #define ABI_FUNC 'lwl_free'
3394 !End of the abilint section
3395 
3396  implicit none
3397 
3398 !Arguments ------------------------------------
3399 !scalars
3400  type(lwl_t),intent(inout) :: lwl
3401 
3402 ! *************************************************************************
3403 
3404  if (allocated(lwl%head)) then
3405    ABI_FREE(lwl%head)
3406  end if
3407  if (allocated(lwl%lwing)) then
3408    ABI_FREE(lwl%lwing)
3409  end if
3410  if (allocated(lwl%uwing)) then
3411    ABI_FREE(lwl%uwing)
3412  end if
3413  if (allocated(lwl%body)) then
3414    ABI_FREE(lwl%body)
3415  end if
3416 
3417 end subroutine lwl_free

m_screening/lwl_init [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 lwl_init

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

3306 subroutine lwl_init(lwl, path, method, cryst, vcp, npwe, gvec, comm)
3307 
3308 
3309 !This section has been created automatically by the script Abilint (TD).
3310 !Do not modify the following lines by hand.
3311 #undef ABI_FUNC
3312 #define ABI_FUNC 'lwl_init'
3313 !End of the abilint section
3314 
3315  implicit none
3316 
3317 !Arguments ------------------------------------
3318 !scalars
3319  integer,intent(in) :: comm,method,npwe
3320  character(len=*),intent(in) :: path
3321  type(crystal_t),intent(in) :: cryst
3322  type(vcoul_t),intent(in) :: vcp
3323  type(lwl_t),intent(out) :: lwl
3324 !arrays
3325  integer,intent(in) :: gvec(3,npwe)
3326 
3327 !Local variables-------------------------------
3328 !scalars
3329  integer,parameter :: master=0
3330  integer :: iomode,my_rank,nproc,unt
3331  character(len=500) :: msg
3332 !arrays
3333 
3334 ! *************************************************************************
3335 
3336  ABI_UNUSED((/cryst%natom, gvec(1,1), vcp%ng/))
3337 
3338  my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
3339 
3340  lwl%fname = path
3341  lwl%method = method
3342  ABI_CHECK(any(method == [1,2,3]), sjoin("Wrong method:", itoa(method)))
3343 
3344  iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF
3345 
3346  ! Only master reads.
3347  if (my_rank == master) then
3348 
3349    select case (iomode)
3350    case (IO_MODE_FORTRAN)
3351      if (open_file(path, msg, newunit=unt, action="read", form="unformatted", status="old") /= 0) then
3352        MSG_ERROR(msg)
3353      end if
3354 
3355      close(unt)
3356 
3357    case default
3358      MSG_ERROR(sjoin("iomode:", itoa(iomode), "is not coded"))
3359    end select
3360  end if
3361 
3362  ! Broad cast data
3363  if (nproc > 1) then
3364  end if
3365 
3366 end subroutine lwl_init

m_screening/lwl_t [ Types ]

[ Top ] [ m_screening ] [ Types ]

NAME

  lwl_t

FUNCTION

SOURCE

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

m_screening/lwl_write [ Functions ]

[ Top ] [ m_screening ] [ Functions ]

NAME

 lwl_write

FUNCTION

INPUTS

OUTPUT

PARENTS

      screening

CHILDREN

SOURCE

3133 subroutine lwl_write(path, cryst, vcp, npwe, nomega, gvec, chi0, chi0_head, chi0_lwing, chi0_uwing, comm)
3134 
3135 
3136 !This section has been created automatically by the script Abilint (TD).
3137 !Do not modify the following lines by hand.
3138 #undef ABI_FUNC
3139 #define ABI_FUNC 'lwl_write'
3140 !End of the abilint section
3141 
3142  implicit none
3143 
3144 !Arguments ------------------------------------
3145 !scalars
3146  integer,intent(in) :: npwe,nomega,comm
3147  character(len=*),intent(in) :: path
3148  type(crystal_t),intent(in) :: cryst
3149  type(vcoul_t),intent(in) :: Vcp
3150 !arrays
3151  integer,intent(in) :: gvec(3,npwe)
3152  complex(gwpc),intent(in) :: chi0(npwe,npwe,nomega)
3153  complex(dpc),intent(inout)  :: chi0_head(3,3,nomega),chi0_lwing(npwe,nomega,3),chi0_uwing(npwe,nomega,3)
3154 
3155 !Local variables-------------------------------
3156 !scalars
3157  integer,parameter :: master=0,prtvol=0
3158  integer :: iw,ii,iomode,unt,my_rank
3159  character(len=500) :: msg
3160  real(dp) :: length
3161  complex(dpc) :: wng(3),em1_00
3162 ! type(hscr_t),intent(out) :: hscr
3163 !arrays
3164  complex(dpc),allocatable :: wtest(:),eps_head(:,:,:)
3165 
3166 ! *************************************************************************
3167 
3168  !if (xmpi_comm_rank(comm) /= master) goto 100
3169  my_rank = xmpi_comm_rank(comm)
3170 
3171  ABI_MALLOC(wtest,(npwe))
3172 
3173  if (prtvol > 0 .and. my_rank == master) then
3174    call wrtout(std_out, "head of chi0")
3175    do iw=1,nomega
3176      call print_arr(chi0_head(:,:,iw),max_r=3,max_c=3)
3177    end do
3178 
3179    do iw=1,nomega
3180      call wrtout(std_out, "symmetrized e_00 via tensor")
3181      wng = MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT)
3182      write(std_out,*) one - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1)*Vcp%vcqlwl_sqrt(1,1)
3183 
3184      call wrtout(std_out, "symmetrized e_0G via tensor")
3185      do ii=1,npwe
3186        wng = chi0_uwing(ii,iw,:)
3187        wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1)
3188      end do
3189      call print_arr(wtest,max_r=9,unit=std_out)
3190 
3191      call wrtout(std_out, "symmetrized e_G0 via tensor")
3192      do ii=1,npwe
3193        wng = chi0_lwing(ii,iw,:)
3194        wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1)
3195      end do
3196      call print_arr(wtest,max_r=9,unit=std_out)
3197    end do
3198  end if
3199 
3200  ! Write chi0 data
3201  iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF
3202 
3203  if (my_rank == master) then
3204    if (iomode == IO_MODE_FORTRAN) then
3205      if (open_file(path,msg,newunit=unt,form="unformatted", action="write") /= 0) then
3206        MSG_ERROR(msg)
3207      end if
3208      !call hscr_io(er%hscr,fform,rdwr,unt,comm,master,iomode)
3209      do iw=1,nomega
3210        write(unt)chi0_head(:,:,iw)
3211      end do
3212      !do iw=1,nomega
3213      !  write(unt)chi0_lwing(:,iw,:)
3214      !end do
3215      !do iw=1,nomega
3216      !  write(unt)chi0_uwing(:,iw,:)
3217      !end do
3218 
3219    else
3220      MSG_ERROR(sjoin("iomode", itoa(iomode), "is not supported"))
3221    end if
3222  end if
3223 
3224  ABI_MALLOC(eps_head,(3,3,nomega))
3225  call mkem1_q0(npwe,1,1,nomega,cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm)
3226 
3227  if (my_rank == master) then
3228    if (iomode == IO_MODE_FORTRAN) then
3229      do iw=1,nomega
3230        write(unt)eps_head(:,:,iw)
3231      end do
3232      !do iw=1,nomega
3233      !  write(unt)chi0_lwing(:,iw,:)
3234      !end do
3235      !do iw=1,nomega
3236      !  write(unt)chi0_uwing(:,iw,:)
3237      !end do
3238    else
3239      MSG_ERROR(sjoin("iomode:", itoa(iomode), "is not supported"))
3240    end if
3241  end if
3242 
3243  ABI_FREE(eps_head)
3244 
3245  if (prtvol > 0 .and. my_rank == master) then
3246    length = normv(GW_Q0_DEFAULT,cryst%gmet,"G")
3247 
3248    do iw=1,nomega
3249      em1_00 = one / vdotw(GW_Q0_DEFAULT/length, MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT/length),cryst%gmet,"G")
3250      call wrtout(std_out, "e^1_{00} from tensor")
3251      write(std_out,*) em1_00
3252 
3253      call wrtout(std_out, "symmetrized e^-1_0G via tensor")
3254      do ii=1,npwe
3255        wng = chi0_uwing(ii,iw,:)
3256        wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G")
3257      end do
3258      wtest(1) = em1_00
3259      call print_arr(wtest,max_r=9,unit=std_out)
3260 
3261      call wrtout(std_out, "symmetrized e^-1_G0 via tensor")
3262      do ii=1,npwe
3263        wng = chi0_lwing(ii,iw,:)
3264        wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G")
3265      end do
3266      wtest(1) = em1_00
3267      call print_arr(wtest,max_r=9,unit=std_out)
3268    end do !iw
3269  end if
3270 
3271  ABI_FREE(wtest)
3272 
3273  if (my_rank == master) then
3274    if (iomode == IO_MODE_FORTRAN) then
3275      close(unt)
3276    else
3277 #ifdef HAVE_NETCDF
3278      NCF_CHECK(nf90_close(unt))
3279 #endif
3280    end if
3281  end if
3282 
3283 !100 call xmpi_barrier(comm)
3284 
3285 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.

PARENTS

      m_screening,screening

CHILDREN

SOURCE

1478 subroutine make_epsm1_driver(iqibz,dim_wing,npwe,nI,nJ,nomega,omega,&
1479   approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,chi0_head,&
1480   chi0_lwing,chi0_uwing,chi0,spectra,comm,&
1481   fxc_ADA) ! optional argument
1482 
1483 
1484 !This section has been created automatically by the script Abilint (TD).
1485 !Do not modify the following lines by hand.
1486 #undef ABI_FUNC
1487 #define ABI_FUNC 'make_epsm1_driver'
1488 !End of the abilint section
1489 
1490  implicit none
1491 
1492 !Arguments ------------------------------------
1493 !scalars
1494  integer,intent(in) :: iqibz,nI,nJ,npwe,nomega,dim_wing,approx_type,option_test,nkxc,nfftot,comm
1495  type(vcoul_t),target,intent(in) :: Vcp
1496  type(spectra_t),intent(out) :: Spectra
1497 !arrays
1498  integer,intent(in) :: ngfft(18),gvec(3,npwe)
1499  complex(gwpc),intent(in) :: kxcg(nfftot,nkxc)
1500  complex(dpc),intent(in) :: omega(nomega)
1501  complex(dpc),intent(inout) :: chi0_lwing(:,:,:)   !(npwe*nI,nomega,dim_wing)
1502  complex(dpc),intent(inout) :: chi0_uwing(:,:,:)   !(npwe*nJ,nomega,dim_wing)
1503  complex(dpc),intent(inout) :: chi0_head(:,:,:)   !(dim_wing,dim_wing,nomega)
1504  complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ,nomega)
1505  complex(gwpc),intent(in),optional :: fxc_ADA(npwe*nI,npwe*nJ)
1506 
1507 !Local variables-------------------------------
1508 !scalars
1509  integer,parameter :: master=0
1510  integer :: i1,i2,ig1,ig2,io,ierr,irank,my_nqlwl !iqlwl
1511  integer :: nor,my_rank,nprocs,comm_self,g1mg2_idx
1512  real(dp) :: ucvol
1513  logical :: is_qeq0,use_MPI
1514  character(len=500) :: msg
1515 !arrays
1516  integer :: omega_distrb(nomega)
1517  integer,allocatable :: istart(:),istop(:)
1518  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1519  real(dp),allocatable :: eelf(:,:),tmp_eelf(:)
1520  complex(dpc),allocatable :: epsm_lf(:,:),epsm_nlf(:,:),tmp_lf(:),tmp_nlf(:)
1521  complex(dpc),allocatable :: buffer_lwing(:,:),buffer_uwing(:,:)
1522  complex(gwpc),allocatable :: kxcg_mat(:,:)
1523 
1524 !bootstrap @WC
1525  integer :: istep,nstep
1526  logical :: converged
1527  real(dp) :: conv_err
1528  real(gwpc) :: chi00_head, fxc_head
1529  !real(gwpc) :: chi00_head, chi00rpa_head, fxc_head
1530  complex(gwpc),allocatable :: vfxc_boot(:,:), chi0_tmp(:,:), chi0_save(:,:,:)
1531  !complex(gwpc),allocatable :: fxc_lrc(:,:), vfxc_boot(:,:), chi0_tmp(:,:), chi0_save(:,:,:)
1532  complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:)
1533 
1534 ! *************************************************************************
1535 
1536  DBG_ENTER("COLL")
1537 
1538  if (nI/=1.or.nJ/=1) then
1539    MSG_ERROR("nI or nJ=/1 not yet implemented")
1540  end if
1541 
1542  nprocs  = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1543 
1544  ! MG TODO We use comm_self for the inversion as the single precision version is not yet available
1545  comm_self = xmpi_comm_self
1546 
1547  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
1548 
1549  is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0)
1550 
1551  omega_distrb = my_rank
1552  use_MPI = .FALSE.
1553  use_MPI = (nprocs>=nomega)  ! Parallelism is not used
1554 
1555  if (use_MPI) then
1556    ! * Initialize distribution table for frequencies.
1557    ABI_MALLOC(istart,(nprocs))
1558    ABI_MALLOC(istop,(nprocs))
1559    call xmpi_split_work2_i4b(nomega,nprocs,istart,istop,msg,ierr)
1560    omega_distrb(:)=xmpi_undefined_rank
1561    do irank=0,nprocs-1
1562      i1 = istart(irank+1)
1563      i2 = istop (irank+1)
1564      if (i1<=i2) omega_distrb(i1:i2) = irank
1565    end do
1566    ABI_FREE(istart)
1567    ABI_FREE(istop)
1568  end if
1569 
1570  ! Initialize container for spectral results
1571  do nor=1,nomega
1572    if (ABS(AIMAG(omega(nor)))>1.e-3) EXIT
1573  end do
1574  nor=nor-1; if (nor==0) nor = 1 ! only imag !?
1575 
1576  if (dim_wing==3) then
1577    call wrtout(std_out,' Analyzing long wavelength limit for several q','COLL')
1578    call spectra_init(Spectra,nor,REAL(omega(1:nor)),Vcp%nqlwl,Vcp%qlwl)
1579    my_nqlwl = 1
1580    !my_nqlwl = dim_wing ! TODO
1581    !ABI_CHECK(dim_wing==SIZE(Vcp%vcqlwl_sqrt,DIM=2),"WRONG DIMS")
1582  else
1583    call spectra_init(Spectra,nor,REAL(omega(1:nor)),1,Vcp%qibz(:,iqibz))
1584    my_nqlwl = 1
1585  end if
1586  !
1587  ! NOTE: all processors have to perform this operation in order to have the
1588  !       epsm1 matrix when performing a sigma calculation starting with the file _SUS
1589  !
1590  ! Temporary arrays to store spectra.
1591  ABI_MALLOC(epsm_lf,(nomega,my_nqlwl))
1592  ABI_MALLOC(epsm_nlf,(nomega,my_nqlwl))
1593  ABI_MALLOC(eelf,(nomega,my_nqlwl))
1594  epsm_lf=czero; epsm_nlf=czero; eelf=zero
1595 
1596  ! Temporary arrays used to store output results.
1597  ABI_MALLOC(tmp_lf, (my_nqlwl))
1598  ABI_MALLOC(tmp_nlf, (my_nqlwl))
1599  ABI_MALLOC(tmp_eelf, (my_nqlwl))
1600 
1601  SELECT CASE (approx_type)
1602 
1603  CASE (0)
1604    ! * RPA: \tepsilon=1 - Vc^{1/2} chi0 Vc^{1/2}
1605    ! * vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff.
1606    do io=1,nomega
1607      if (omega_distrb(io) == my_rank) then
1608        !write(std_out,*)"dim_wing",dim_wing
1609        call rpa_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),my_nqlwl,dim_wing,&
1610 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),&
1611 &        tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1612 
1613          ! Store results.
1614          epsm_lf(io,:) = tmp_lf
1615          epsm_nlf(io,:) = tmp_nlf
1616          eelf(io,:) = tmp_eelf
1617      end if
1618    end do ! nomega
1619 
1620  CASE (1)
1621    ! Vertex correction from Adiabatic TDDFT. chi_{G1,G2} = [\delta -\chi0 (vc+kxc)]^{-1}_{G1,G3} \chi0_{G3,G2}
1622    ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1623    ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1624 
1625    ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh.
1626    ABI_STAT_MALLOC(kxcg_mat,(npwe,npwe), ierr)
1627    ABI_CHECK(ierr==0, "out-of-memory kxcg_mat")
1628 
1629    ierr=0
1630    do ig2=1,npwe
1631      do ig1=1,npwe
1632        g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft)
1633        if (g1mg2_idx>0) then
1634          kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1)
1635        else
1636          ierr=ierr+1
1637          kxcg_mat(ig1,ig2) = czero
1638        end if
1639      end do
1640    end do
1641 
1642    if (ierr/=0) then
1643      write(msg,'(a,i4,3a)')&
1644 &     'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,&
1645 &     'Enlarge the FFT mesh to get rid of this problem. '
1646      MSG_WARNING(msg)
1647    end if
1648 
1649    !FIXME "recheck TDDFT code and parallel"
1650    ABI_CHECK(nkxc==1,"nkxc/=1 not coded")
1651    do io=1,nomega
1652      if (omega_distrb(io) == my_rank) then
1653        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),&
1654 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm)
1655 
1656        ! Store results.
1657        epsm_lf(io,:) = tmp_lf
1658        epsm_nlf(io,:) = tmp_nlf
1659        eelf(io,:) = tmp_eelf
1660      end if
1661    end do
1662 
1663    ABI_FREE(kxcg_mat)
1664    do io=1,nomega
1665      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1666      call wrtout(std_out,msg,'COLL')
1667      call print_arr(chi0(:,:,io),mode_paral='PERS')
1668    end do
1669 
1670  CASE (2)
1671    ! ADA nonlocal vertex correction contained in fxc_ADA
1672    MSG_WARNING('Entered fxc_ADA branch: EXPERIMENTAL!')
1673    ! Test that argument was passed
1674    if (.not.present(fxc_ADA)) then
1675      MSG_ERROR('make_epsm1_driver was not called with optional argument fxc_ADA')
1676    end if
1677    ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded")
1678 
1679    do io=1,nomega
1680      if (omega_distrb(io) == my_rank) then
1681        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),fxc_ADA,option_test,my_nqlwl,dim_wing,omega(io),&
1682 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm)
1683 
1684        ! Store results.
1685        epsm_lf(io,:) = tmp_lf
1686        epsm_nlf(io,:) = tmp_nlf
1687        eelf(io,:) = tmp_eelf
1688      end if
1689    end do
1690 
1691    do io=1,nomega
1692      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1693      call wrtout(std_out,msg,'COLL')
1694      call print_arr(chi0(:,:,io),mode_paral='PERS')
1695    end do
1696 
1697  CASE (4)
1698    !@WC bootstrap vertex correction, Sharma et al. PRL 107, 196401 (2011) [[cite:Sharma2011]] 
1699    ABI_STAT_MALLOC(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1700    ABI_CHECK(ierr==0, "out-of-memory in vfxc_boot")
1701    ABI_STAT_MALLOC(chi0_tmp,(npwe*nI,npwe*nJ), ierr)
1702    ABI_CHECK(ierr==0, "out-of-memory in chi0_tmp")
1703    ABI_STAT_MALLOC(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1704    ABI_CHECK(ierr==0, "out-of-memory in chi0_save")
1705 
1706    if (iqibz==1) then
1707      vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! Use Coulomb term for q-->0
1708    else
1709      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
1710    end if
1711 
1712    chi0_save = chi0 ! a copy of chi0 (ks)
1713    nstep = 50 ! max iteration steps
1714    chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1715    fxc_head = czero; vfxc_boot = czero; chi0_tmp = czero
1716    epsm_lf = czero; epsm_nlf = czero; eelf = zero
1717    write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head
1718    call wrtout(std_out,msg,'COLL')
1719 
1720    do istep=1,nstep
1721      chi0 = chi0_save
1722      do io=1,1 ! static
1723        !if (omega_distrb(io) == my_rank) then
1724        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),&
1725 &       chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1726        epsm_lf(io,:) = tmp_lf
1727        epsm_nlf(io,:) = tmp_nlf
1728        eelf(io,:) = tmp_eelf
1729        !end if
1730      end do
1731 
1732      conv_err = smallest_real
1733      do ig2=1,npwe*nJ
1734        do ig1=1,npwe*nI
1735          conv_err= MAX(conv_err, ABS(chi0(ig1,ig2,1) - chi0_tmp(ig1,ig2)))
1736        end do
1737      end do
1738      converged = (conv_err <= tol4)
1739      write(msg,'(a,i4,a,f10.6)') ' => bootstrap itr# ', istep, ', eps^-1 max error: ', conv_err
1740      call wrtout(std_out,msg,'COLL')
1741      write(msg,'(a,2f10.6)')  '    eps^-1(head):   ', chi0(1,1,1)
1742      call wrtout(std_out,msg,'COLL')
1743      write(msg,'(a,2f10.6)')  '    v^-1*fxc(head): ', fxc_head
1744      call wrtout(std_out,msg,'COLL')
1745 
1746      if (converged) then
1747        write(msg,'(a,i4,a)') ' => bootstrap fxc converged after ', istep, ' iterations'
1748        call wrtout(std_out,msg,'COLL')
1749        chi0 = chi0_save
1750        do io=1,nomega
1751          if (omega_distrb(io) == my_rank) then
1752            call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1753 &           chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1754            epsm_lf(io,:) = tmp_lf
1755            epsm_nlf(io,:) = tmp_nlf
1756            eelf(io,:) = tmp_eelf
1757          end if
1758        end do
1759        write(msg, '(a,2f10.6)') ' ->   eps^-1(head): ', chi0(1,1,1)
1760        call wrtout(std_out,msg,'COLL')
1761        write(msg,'(a,2f10.6)')  ' -> v^-1*fxc(head): ', fxc_head
1762        call wrtout(std_out,msg,'COLL')
1763        exit
1764      else if (istep < nstep) then
1765        chi0_tmp = chi0(:,:,1)
1766        vfxc_boot = chi0(:,:,1)/chi00_head ! full G vectors
1767        !vfxc_boot = czero; vfxc_boot(1,1) = chi0(1,1,1)/chi00_head ! head only
1768        fxc_head = vfxc_boot(1,1)
1769        do ig1=1,npwe
1770          vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:)
1771        end do
1772      else
1773        write(msg,'(a,i4,a)') ' -> bootstrap fxc not converged after ', nstep, ' iterations'
1774        MSG_WARNING(msg)
1775        ! proceed to calculate the dielectric function even fxc is not converged
1776        chi0 = chi0_save
1777        do io=1,nomega
1778          if (omega_distrb(io) == my_rank) then
1779            call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1780 &            chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1781            epsm_lf(io,:) = tmp_lf
1782            epsm_nlf(io,:) = tmp_nlf
1783            eelf(io,:) = tmp_eelf
1784          end if
1785        end do
1786      end if
1787    end do
1788 
1789    ABI_FREE(chi0_tmp)
1790    ABI_FREE(chi0_save)
1791    ABI_FREE(vfxc_boot)
1792 
1793    do io=1,nomega
1794      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1795      call wrtout(std_out,msg,'COLL')
1796      call print_arr(chi0(:,:,io),mode_paral='PERS')
1797    end do
1798 
1799  CASE(5)
1800    !@WC: one-shot scalar bootstrap approximation
1801    ABI_STAT_MALLOC(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1802    ABI_CHECK(ierr==0, "out-of-memory in vfxc_boot")
1803    ABI_STAT_MALLOC(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1804    ABI_CHECK(ierr==0, "out-of-memory in chi0_save")
1805 
1806    if (iqibz==1) then
1807      vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! Use Coulomb term for q-->0
1808    else
1809      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
1810    end if
1811 
1812    chi0_save = chi0 ! a copy of chi0
1813    fxc_head = czero; vfxc_boot = czero;
1814    epsm_lf = czero; epsm_nlf = czero; eelf = zero
1815    chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1816    write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head
1817    call wrtout(std_out,msg,'COLL')
1818 
1819    fxc_head = vc_sqrt(1)**2/chi00_head + vc_sqrt(1)**2/chi00_head - vc_sqrt(1)**2
1820    fxc_head = 0.5*fxc_head + 0.5*sqrt(fxc_head**2 - 4.0*vc_sqrt(1)**4/(chi00_head*chi00_head))
1821    vfxc_boot(1,1) = fxc_head
1822    write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head/vc_sqrt(1)**2
1823    call wrtout(std_out,msg,'COLL')
1824 
1825    chi0 = chi0_save
1826    do io=1,nomega
1827      if (omega_distrb(io) == my_rank) then
1828        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1829 &         chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1830        epsm_lf(io,:) = tmp_lf
1831        epsm_nlf(io,:) = tmp_nlf
1832        eelf(io,:) = tmp_eelf
1833      end if
1834    end do
1835    write(msg,'(a,2f10.6)')  '    eps^-1(head):   ',chi0(1,1,1)
1836    call wrtout(std_out,msg,'COLL')
1837 
1838    ABI_FREE(chi0_save)
1839    ABI_FREE(vfxc_boot)
1840 
1841    do io=1,nomega
1842      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1843      call wrtout(std_out,msg,'COLL')
1844      call print_arr(chi0(:,:,io),mode_paral='PERS')
1845    end do
1846 
1847 CASE(6)
1848    !@WC: RPA bootstrap by Rigamonti et al. (PRL 114, 146402) [[cite:Rigamonti2015]]
1849    !@WC: and Berger (PRL 115, 137402) [[cite:Berger2015]] 
1850    ABI_STAT_MALLOC(vfxc_boot,(npwe*nI,npwe*nJ), ierr)
1851    ABI_CHECK(ierr==0, "out-of-memory in vfxc_boot")
1852    ABI_STAT_MALLOC(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr)
1853    ABI_CHECK(ierr==0, "out-of-memory in chi0_save")
1854 
1855    if (iqibz==1) then
1856      vc_sqrt => Vcp%vcqlwl_sqrt(:,1)  ! Use Coulomb term for q-->0
1857    else
1858      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
1859    end if
1860 
1861    chi0_save = chi0 ! a copy of chi0
1862    fxc_head = czero; vfxc_boot = czero;
1863    epsm_lf = czero; epsm_nlf = czero; eelf = zero
1864    chi00_head = chi0(1,1,1)*vc_sqrt(1)**2
1865    write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head
1866    call wrtout(std_out,msg,'COLL')
1867 
1868    ! static
1869    io = 1
1870    call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),&
1871 &    chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1872    epsm_lf(1,:) = tmp_lf
1873 
1874    vfxc_boot(1,1) = 1.0/(chi00_head * epsm_lf(1,1))
1875    fxc_head = vfxc_boot(1,1)
1876    do ig1=1,npwe
1877      vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:)
1878    end do
1879    write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head
1880    call wrtout(std_out,msg,'COLL')
1881 
1882    chi0 = chi0_save
1883    do io=1,nomega
1884      if (omega_distrb(io) == my_rank) then
1885        call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),&
1886 &        chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self)
1887        epsm_lf(io,:) = tmp_lf
1888        epsm_nlf(io,:) = tmp_nlf
1889        eelf(io,:) = tmp_eelf
1890      end if
1891    end do
1892    write(msg,'(a,2f10.6)')  '    eps^-1(head):   ',chi0(1,1,1)
1893    call wrtout(std_out,msg,'COLL')
1894 
1895    ABI_FREE(chi0_save)
1896    ABI_FREE(vfxc_boot)
1897 
1898    do io=1,nomega
1899      write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]'
1900      call wrtout(std_out,msg,'COLL')
1901      call print_arr(chi0(:,:,io),mode_paral='PERS')
1902    end do
1903 
1904  CASE DEFAULT
1905    MSG_BUG(sjoin('Wrong approx_type:',itoa(approx_type)))
1906  END SELECT
1907 
1908  if (use_MPI) then
1909    ! Collect results on each node.
1910    ABI_MALLOC(buffer_lwing, (size(chi0_lwing,dim=1), size(chi0_lwing, dim=3)))
1911    ABI_MALLOC(buffer_uwing, (size(chi0_uwing,dim=1), size(chi0_uwing, dim=3)))
1912 
1913    do io=1,nomega
1914      if (omega_distrb(io)/=my_rank) then
1915        ! Zero arrays.
1916        chi0(:,:,io) = czero_gw
1917        if(dim_wing>0) then
1918           chi0_lwing(:,io,:) = zero
1919           chi0_uwing(:,io,:) = zero
1920           chi0_head(:,:,io)  = czero
1921        endif
1922        epsm_lf(io,:) = czero
1923        epsm_nlf(io,:) = czero
1924        eelf(io,:) = zero
1925      end if
1926 
1927      call xmpi_sum(chi0(:,:,io), comm,ierr)
1928 
1929      if (dim_wing>0) then
1930         ! Build contiguous arrays
1931         buffer_lwing = chi0_lwing(:,io,:)
1932         buffer_uwing = chi0_uwing(:,io,:)
1933         call xmpi_sum(buffer_lwing,comm,ierr)
1934         call xmpi_sum(buffer_uwing,comm,ierr)
1935         chi0_lwing(:,io,:) = buffer_lwing
1936         chi0_uwing(:,io,:) = buffer_uwing
1937         if (size(chi0_head(:,:,io))/= zero) then
1938           call xmpi_sum(chi0_head(:,:,io),comm,ierr)
1939         end if
1940      end if
1941 
1942    end do ! iw
1943 
1944    call xmpi_sum(epsm_lf, comm,ierr )
1945    call xmpi_sum(epsm_nlf,comm,ierr)
1946    call xmpi_sum(eelf,    comm,ierr)
1947    ABI_FREE(buffer_lwing)
1948    ABI_FREE(buffer_uwing)
1949  end if
1950 
1951  ! Save results in Spectra%, mind the slicing.
1952  Spectra%emacro_nlf(:,:) = epsm_nlf(1:nor,:)
1953  Spectra%emacro_lf (:,:) = epsm_lf (1:nor,:)
1954  Spectra%eelf      (:,:) = eelf    (1:nor,:)
1955 
1956  ABI_FREE(epsm_lf)
1957  ABI_FREE(epsm_nlf)
1958  ABI_FREE(eelf)
1959  ABI_FREE(tmp_lf)
1960  ABI_FREE(tmp_nlf)
1961  ABI_FREE(tmp_eelf)
1962 
1963  DBG_EXIT("COLL")
1964 
1965 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

PARENTS

SOURCE

2786 elemental function mdielf_bechstedt(eps_inf,qnrm,rhor) result(mdielf)
2787 
2788 
2789 !This section has been created automatically by the script Abilint (TD).
2790 !Do not modify the following lines by hand.
2791 #undef ABI_FUNC
2792 #define ABI_FUNC 'mdielf_bechstedt'
2793 !End of the abilint section
2794 
2795  implicit none
2796 
2797 !Arguments ------------------------------------
2798 !scalars
2799  real(dp),intent(in) :: eps_inf,qnrm,rhor
2800  real(dp) :: mdielf
2801 
2802 ! *************************************************************************
2803 
2804  mdielf = one + &
2805 &  one / ( one/(eps_inf-one) + (qnrm/k_thfermi(rhor))**2 + (three*qnrm**4)/(four*k_fermi(rhor)**2 * k_thfermi(rhor)**2) )
2806 
2807 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

PARENTS

      mrgscr,sigma

CHILDREN

SOURCE

 925 subroutine mkdump_Er(Er,Vcp,npwe,gvec,nkxc,kxcg,id_required,approx_type,&
 926 &                    ikxc_required,option_test,fname_dump,iomode,&
 927 &                    nfftot,ngfft,comm,fxc_ADA)
 928 
 929 
 930 !This section has been created automatically by the script Abilint (TD).
 931 !Do not modify the following lines by hand.
 932 #undef ABI_FUNC
 933 #define ABI_FUNC 'mkdump_Er'
 934 !End of the abilint section
 935 
 936  implicit none
 937 
 938 !Arguments ------------------------------------
 939 !scalars
 940  integer,intent(in) :: id_required,approx_type,option_test,ikxc_required,nkxc
 941  integer,intent(in) :: iomode,nfftot,npwe,comm
 942  type(Epsilonm1_results),intent(inout) :: Er
 943  type(vcoul_t),intent(in) :: Vcp
 944  character(len=*),intent(in) :: fname_dump
 945 !arrays
 946  integer,intent(in) :: ngfft(18),gvec(3,npwe)
 947  complex(gwpc),intent(in) :: kxcg(nfftot,nkxc)
 948  complex(gwpc),intent(in), optional :: fxc_ADA(npwe*Er%nI,npwe*Er%nJ,Er%nqibz)
 949 
 950 !Local variables-------------------------------
 951 !scalars
 952  integer,parameter :: master=0
 953  integer :: dim_wing,iqibz,is_qeq0,mqmem_,npwe_asked
 954  integer :: unt_dump,fform,rdwr,ierr
 955  integer :: my_rank,comm_self
 956  real(dp) :: ucvol
 957  character(len=500) :: msg
 958  character(len=fnlen) :: ofname
 959  character(len=nctk_slen) :: in_varname,out_varname
 960  type(hscr_t) :: Hscr_cp
 961  type(spectra_t) :: spectra
 962 !arrays
 963  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
 964  complex(gwpc),allocatable :: epsm1(:,:,:)
 965  complex(dpc),allocatable :: dummy_lwing(:,:,:),dummy_uwing(:,:,:),dummy_head(:,:,:)
 966 
 967 ! *********************************************************************
 968 
 969  DBG_ENTER("COLL")
 970 
 971  ABI_CHECK(id_required==4,'Value of id_required not coded')
 972  ABI_CHECK(npwe==Er%npwe,"mismatch in npwe")
 973 
 974  my_rank = xmpi_comm_rank(comm); comm_self = xmpi_comm_self
 975  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
 976 
 977  ! if (Er%ID/=0) call reset_Epsilonm1(Er)
 978  Er%ID=id_required
 979 
 980  ofname = fname_dump
 981  in_varname = ncname_from_id(er%hscr%id)
 982  out_varname = ncname_from_id(id_required)
 983 
 984  write(std_out,*)'Er%ID: ',Er%ID,', Er%Hscr%ID: ',Er%Hscr%ID
 985 
 986  if (Er%ID==Er%Hscr%ID) then
 987    ! === The two-point function we are asking for is already stored on file ===
 988    ! * According to mqmem either read and store the entire matrix in memory or do nothing.
 989 
 990    if (Er%mqmem>0) then
 991      ! In-core solution.
 992      write(msg,'(a,f12.1,a)')' Memory needed for Er%epsm1 = ',two*gwpc*npwe**2*Er%nomega*Er%nqibz*b2Mb,' [Mb]'
 993      call wrtout(std_out,msg,'PERS')
 994      ABI_STAT_MALLOC(Er%epsm1,(npwe,npwe,Er%nomega,Er%nqibz), ierr)
 995      ABI_CHECK(ierr==0, "Out-of-memory in Er%epsm1 (in-core)")
 996 
 997      if (iomode == IO_MODE_MPI) then
 998        !call wrtout(std_out,sjoin(ABI_FUNC, " read_screening with MPI_IO"))
 999        MSG_WARNING("SUSC files is buggy. Using Fortran IO")
1000        call read_screening(in_varname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm)
1001      else
1002        call read_screening(in_varname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm)
1003      end if
1004    else
1005      ! Out-of-core solution ===
1006      MSG_COMMENT("mqmem==0 => allocating a single q-slice of (W|chi0) (slower but less memory).")
1007      continue
1008    end if
1009 
1010    return
1011 
1012  else
1013    ! === The matrix stored on file do not correspond to the quantity required ===
1014    ! * Presently only the transformation chi0 => e^-1 is coded
1015    ! * According to Er%mqmem either calculate e^-1 dumping the result to a file
1016    !   for a subsequent use or calculate e^-1 keeping everything in memory.
1017 
1018    if (Er%mqmem==0) then
1019      ! === Open file and write the header for the SCR file ===
1020      ! * For the moment only master works.
1021 
1022      if (my_rank==master) then
1023        if (iomode == IO_MODE_ETSF) then
1024 #ifdef HAVE_NETCDF
1025           ofname = nctk_ncify(ofname)
1026           NCF_CHECK(nctk_open_create(unt_dump, ofname, xmpi_comm_self))
1027 #endif
1028        else
1029          if (open_file(ofname,msg,newunit=unt_dump,form="unformatted",status="unknown",action="write") /= 0) then
1030            MSG_ERROR(msg)
1031          end if
1032        end if
1033        call wrtout(std_out,sjoin('mkdump_Er: calculating and writing epsilon^-1 matrix on file: ',ofname),'COLL')
1034 
1035        ! Update the entries in the header that have been modified.
1036        ! TODO, write function to return title, just for info
1037        call hscr_copy(Er%Hscr,Hscr_cp)
1038        Hscr_cp%ID = id_required
1039        Hscr_cp%ikxc = ikxc_required
1040        Hscr_cp%test_type = option_test
1041        Hscr_cp%titles(1)  = 'SCR file: epsilon^-1'
1042        Hscr_cp%titles(2)  = 'TESTPARTICLE'
1043        ! Treat the case in which a smaller matrix is used.
1044        Hscr_cp%npwe = npwe
1045 
1046        rdwr=2; fform=Hscr_cp%fform
1047        call hscr_io(hscr_cp,fform,rdwr,unt_dump,comm_self,master,iomode)
1048        call hscr_free(Hscr_cp)
1049 
1050        ABI_STAT_MALLOC(epsm1,(npwe,npwe,Er%nomega), ierr)
1051        ABI_CHECK(ierr==0, "out of memory in epsm1")
1052 
1053        do iqibz=1,Er%nqibz
1054          is_qeq0=0
1055          if (normv(Er%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1056          ! FIXME there's a problem with SUSC files and MPI-IO
1057          !if (iomode == IO_MODE_MPI) then
1058          !  MSG_WARNING("SUSC files is buggy. Using Fortran IO")
1059          !  call read_screening(in_varname,Er%fname,npwe,1,Er%nomega,epsm1,IO_MODE_FORTRAN,comm_self,iqiA=iqibz)
1060          !else
1061          call read_screening(in_varname,Er%fname,npwe,1,Er%nomega,epsm1,iomode,comm_self,iqiA=iqibz)
1062          !end if
1063 
1064          dim_wing=0; if (is_qeq0==1) dim_wing=3
1065          ABI_MALLOC(dummy_lwing,(npwe*Er%nI,Er%nomega,dim_wing))
1066          ABI_MALLOC(dummy_uwing,(npwe*Er%nJ,Er%nomega,dim_wing))
1067          ABI_MALLOC(dummy_head,(dim_wing,dim_wing,Er%nomega))
1068 
1069          if (approx_type<2 .or. approx_type>3) then ! bootstrap
1070            MSG_WARNING('Entering out-of core RPA or Kxc branch')
1071            call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,&
1072 &                    approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1073 &                    dummy_lwing,dummy_uwing,epsm1,spectra,comm_self)
1074          else
1075            MSG_WARNING('Entering out-of core fxc_ADA branch')
1076            call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,&
1077 &                    approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1078 &                    dummy_lwing,dummy_uwing,epsm1,spectra,comm_self,fxc_ADA(:,:,iqibz))
1079          end if
1080 
1081          ABI_FREE(dummy_head)
1082          ABI_FREE(dummy_uwing)
1083          ABI_FREE(dummy_lwing)
1084 
1085          if (is_qeq0==1) then
1086            call spectra_repr(spectra,msg)
1087            call wrtout(std_out,msg,'COLL')
1088            call wrtout(ab_out,msg,'COLL')
1089          end if
1090          call spectra_free(spectra)
1091 
1092          call write_screening(out_varname,unt_dump,iomode,npwe,Er%nomega,iqibz,epsm1)
1093        end do
1094 
1095        if (iomode == IO_MODE_ETSF) then
1096 #ifdef HAVE_NETCDF
1097          NCF_CHECK(nf90_close(unt_dump))
1098 #endif
1099        else
1100          close(unt_dump)
1101        endif
1102 
1103        ABI_FREE(epsm1)
1104      end if !master
1105 
1106      ! Master broadcasts ofname.
1107      ! NOTE: A synchronization is required here, else the other procs start to read the
1108      ! SCR file before it is written by the master. xmpi_bcast will synch the procs.
1109      call xmpi_bcast(ofname,  master, comm, ierr)
1110 
1111      ! Now Er% "belongs" to the file "ofname", thus
1112      ! each proc has to destroy and re-initialize the object.
1113      call em1results_free(Er)
1114 
1115      mqmem_=Er%mqmem; npwe_asked=npwe
1116      call init_Er_from_file(Er,ofname,mqmem_,npwe_asked,comm)
1117 
1118      !Now Er% has been reinitialized and ready-to-use.
1119      Er%id = id_required
1120      call em1results_print(Er)
1121    else
1122      ! ========================
1123      ! === In-core solution ===
1124      ! ========================
1125      ABI_STAT_MALLOC(Er%epsm1,(npwe,npwe,Er%nomega,Er%nqibz), ierr)
1126      ABI_CHECK(ierr==0, 'Out-of-memory in Er%epsm1 (in-core)')
1127 
1128      ! FIXME there's a problem with SUSC files and MPI-IO
1129      !if (iomode == IO_MODE_MPI) then
1130      !  !call wrtout(std_out,sjoin(ABI_FUNC, "read_screening with MPI_IO"))
1131      !  MSG_WARNING("SUSC files is buggy. Using Fortran IO")
1132      !  call read_screening(in_varname,Er%fname,npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm)
1133      !else
1134      call read_screening(in_varname,Er%fname,npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm)
1135      !end if
1136 
1137      do iqibz=1,Er%nqibz
1138        is_qeq0=0; if (normv(Er%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1
1139 
1140        dim_wing=0; if (is_qeq0==1) dim_wing=3 ! FIXME
1141        ABI_MALLOC(dummy_lwing,(npwe*Er%nI,Er%nomega,dim_wing))
1142        ABI_MALLOC(dummy_uwing,(npwe*Er%nJ,Er%nomega,dim_wing))
1143        ABI_MALLOC(dummy_head,(dim_wing,dim_wing,Er%nomega))
1144 
1145        if (approx_type<2 .or. approx_type>3) then
1146          MSG_WARNING('Entering in-core RPA and Kxc branch')
1147          call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,&
1148 &                  approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1149 &                  dummy_lwing,dummy_uwing,Er%epsm1(:,:,:,iqibz),spectra,comm)
1150        else
1151          MSG_WARNING('Entering in-core fxc_ADA branch')
1152          call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,&
1153 &                  approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,&
1154 &                  dummy_lwing,dummy_uwing,Er%epsm1(:,:,:,iqibz),spectra,comm,fxc_ADA=fxc_ADA(:,:,iqibz))
1155        end if
1156 
1157        ABI_FREE(dummy_lwing)
1158        ABI_FREE(dummy_uwing)
1159        ABI_FREE(dummy_head)
1160 
1161        if (is_qeq0==1) then
1162          call spectra_repr(spectra,msg)
1163          call wrtout(std_out,msg,'COLL')
1164          call wrtout(ab_out,msg,'COLL')
1165        end if
1166 
1167        call spectra_free(spectra)
1168      end do
1169 
1170      Er%id = id_required
1171      call em1results_print(Er)
1172    end if
1173  end if
1174 
1175  DBG_EXIT("COLL")
1176 
1177 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

PARENTS

      m_screening

CHILDREN

SOURCE

2396 subroutine mkem1_q0(npwe,n1,n2,nomega,Cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm)
2397 
2398 
2399 !This section has been created automatically by the script Abilint (TD).
2400 !Do not modify the following lines by hand.
2401 #undef ABI_FUNC
2402 #define ABI_FUNC 'mkem1_q0'
2403 !End of the abilint section
2404 
2405  implicit none
2406 
2407 !Arguments ------------------------------------
2408 !scalars
2409  integer,intent(in) :: npwe,nomega,n1,n2,comm
2410  type(crystal_t),intent(in) :: Cryst
2411  type(vcoul_t),intent(in) :: Vcp
2412 !arrays
2413  integer,intent(in) :: gvec(3,npwe)
2414  complex(gwpc),intent(in) :: chi0(npwe*n1,npwe*n2,nomega)
2415  complex(dpc),intent(inout) :: chi0_lwing(npwe*n1,nomega,3)
2416  complex(dpc),intent(inout) :: chi0_uwing(npwe*n2,nomega,3)
2417  complex(dpc),intent(inout) :: chi0_head(3,3,nomega)
2418  complex(dpc),intent(out) :: eps_head(3,3,nomega)
2419 
2420 !Local variables ------------------------------
2421 !scalars
2422  integer :: iw,ig,ig1,ig2,idir,jdir
2423 !arrays
2424  real(dp),allocatable :: modg_inv(:)
2425  complex(dpc),allocatable :: eps_lwing(:,:),eps_uwing(:,:),eps_body(:,:),cvec(:)
2426 
2427 !************************************************************************
2428 
2429  ABI_CHECK(npwe/=1,"npwe must be >1")
2430  ABI_UNUSED(comm)
2431 
2432  ! Precompute 1/|G|.
2433  ABI_MALLOC(modg_inv,(npwe-1))
2434  do ig=1,npwe-1
2435    modg_inv(ig) = one/normv(gvec(:,ig+1),Cryst%gmet,'G')
2436  end do
2437 
2438  ABI_MALLOC(eps_uwing,((npwe-1)*n1,3))
2439  ABI_MALLOC(eps_lwing,((npwe-1)*n2,3))
2440  ABI_MALLOC(eps_body,(npwe-1,npwe-1))
2441  ABI_MALLOC(cvec,(npwe-1))
2442 
2443  do iw=1,nomega
2444    !
2445    ! Head and wings of the symmetrized epsilon.
2446    eps_head(:,:,iw) = -four_pi*chi0_head(:,:,iw)
2447    do idir=1,3
2448      eps_head(idir,idir,iw) = one + eps_head(idir,idir,iw)
2449      eps_lwing(:,idir) = -four_pi * modg_inv * chi0_lwing(2:,iw,idir)
2450      eps_uwing(:,idir) = -four_pi * modg_inv * chi0_uwing(2:,iw,idir)
2451      !eps_lwing(:,idir) = -chi0_lwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1)
2452      !eps_uwing(:,idir) = -chi0_uwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1)
2453    end do
2454 
2455    write(std_out,*)" espilon head"
2456    call print_arr(eps_head(:,:,iw))
2457    !
2458    ! Construct the body of the symmetrized epsilon then invert it.
2459    do ig2=1,npwe-1
2460      do ig1=1,npwe-1
2461        eps_body(ig1,ig2) = -four_pi * modg_inv(ig1)*chi0(ig1+1,ig2+1,iw )*modg_inv(ig2)
2462        !eps_body(ig1,ig2) = -Vcp%vcqlwl_sqrt(ig1+1,1)*chi0(ig1+1,ig2+1,iw)* Vcp%vcqlwl_sqrt(ig2+1,1)
2463      end do
2464      eps_body(ig2,ig2) = one + eps_body(ig2,ig2)
2465    end do
2466 
2467    call xginv(eps_body,npwe-1)
2468    !
2469    ! Overwrite chi0_head and chi0_wings with the head and the wings of the inverse dielectric matrix.
2470    do jdir=1,3
2471      !
2472      ! Head.
2473      cvec=czero
2474      do idir=1,3
2475        cvec = cvec + two_pi*Cryst%gmet(jdir,idir)*MATMUL(eps_body,eps_lwing(:,idir)) ! as we work in reciprocal coords.
2476      end do
2477      !cvec = MATMUL(eps_body,eps_lwing(:,jdir))
2478      do idir=1,3
2479        chi0_head(idir,jdir,iw) = eps_head(idir,jdir,iw) - xdotu(npwe-1,eps_uwing(:,idir),1,cvec,1)
2480      end do
2481      !
2482      ! Now the wings.
2483      chi0_uwing(2:,iw,jdir) = -MATMUL(eps_uwing(:,jdir),eps_body)
2484      chi0_lwing(2:,iw,jdir) = -MATMUL(eps_body,eps_lwing(:,jdir))
2485      !
2486    end do !jdir
2487 
2488    call wrtout(std_out, "espilon^1 head after block inversion", "COLL")
2489    call print_arr(chi0_head(:,:,iw))
2490    !
2491    ! Change the body but do not add the corrections due to the head and the wings.
2492    ! since they can be obtained on the fly from eps_body and the wings of eps^{-1}.
2493    !%chi0(2:,2:,iw) = eps_body
2494  end do !iw
2495 
2496  ABI_FREE(modg_inv)
2497  ABI_FREE(cvec)
2498  ABI_FREE(eps_lwing)
2499  ABI_FREE(eps_uwing)
2500  ABI_FREE(eps_body)
2501 
2502  RETURN
2503  ABI_UNUSED(Vcp%ng)
2504 
2505 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.

OUTPUT

SIDE EFFECTS

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

PARENTS

      m_screening

CHILDREN

SOURCE

2005 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)
2006 
2007 
2008 !This section has been created automatically by the script Abilint (TD).
2009 !Do not modify the following lines by hand.
2010 #undef ABI_FUNC
2011 #define ABI_FUNC 'rpa_symepsm1'
2012 !End of the abilint section
2013 
2014  implicit none
2015 
2016 !Arguments ------------------------------------
2017 !scalars
2018  integer,intent(in) :: iqibz,nI,nJ,npwe,dim_wing,my_nqlwl,comm
2019  type(vcoul_t),target,intent(in) :: Vcp
2020 !arrays
2021  complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ)
2022  complex(dpc),intent(inout) :: chi0_lwing(:,:) !(npwe*nI,dim_wing)
2023  complex(dpc),intent(inout) :: chi0_uwing(:,:) !(npwe*nJ,dim_wing)
2024  complex(dpc),intent(inout) :: chi0_head(:,:) !(dim_wing,dim_wing)
2025  real(dp),intent(out) :: eelf(my_nqlwl)
2026  complex(dpc),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl)
2027 
2028 !Local variables-------------------------------
2029 !scalars
2030  integer,parameter :: master=0,prtvol=0
2031  integer :: ig1,ig2,iqlwl,my_rank,nprocs
2032  real(dp) :: ucvol
2033  logical :: is_qeq0
2034  !character(len=500) :: msg
2035 !arrays
2036  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2037  complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:)
2038  complex(gwpc),allocatable :: chi0_save(:,:)
2039 
2040 ! *************************************************************************
2041 
2042  ABI_UNUSED(chi0_head(1,1))
2043 
2044  if (nI/=1.or.nJ/=1) then
2045    MSG_ERROR("nI or nJ=/1 not yet implemented")
2046  end if
2047 
2048  nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2049 
2050  call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol)
2051 
2052  is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0)
2053  if (is_qeq0) then
2054    ABI_CHECK(iqibz==1,"q is 0 but iq_ibz /= 1")
2055  end if
2056  !
2057  if (my_nqlwl>1) then
2058    ABI_MALLOC(chi0_save,(npwe*nI,npwe*nJ))
2059    chi0_save = chi0
2060  end if
2061  !
2062  ! Symmetrized RPA epsilon = 1 - Vc^{1/2} chi0 Vc^{1/2}
2063  !   * vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff.
2064  !
2065  ! * Loop over small q"s (if any) to treat the nonanalytical behavior.
2066  do iqlwl=my_nqlwl,1,-1
2067    !
2068    if (my_nqlwl>1) then
2069      chi0(:,:) = chi0_save           ! restore pristine polarizability
2070      chi0(:,1) = chi0_lwing(:,iqlwl) ! change the wings
2071      chi0(1,:) = chi0_uwing(:,iqlwl)
2072    end if
2073 
2074    if (iqibz==1) then
2075      vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl)  ! Use Coulomb term for q-->0
2076    else
2077      vc_sqrt => Vcp%vc_sqrt(:,iqibz)
2078    end if
2079 
2080    do ig2=1,npwe*nJ
2081      do ig1=1,npwe*nI
2082        chi0(ig1,ig2)=-vc_sqrt(ig1)*chi0(ig1,ig2)*vc_sqrt(ig2)
2083      end do
2084      chi0(ig2,ig2)=one+chi0(ig2,ig2)
2085    end do
2086 
2087    epsm_nlf(iqlwl)=chi0(1,1) ! * chi0, now contains \tepsilon.
2088 
2089    if (prtvol > 0) then
2090      call wrtout(std_out,' Symmetrical epsilon(G,G'') ','COLL')
2091      call print_arr(chi0, unit=std_out)
2092    end if
2093    !
2094    ! === Invert tepsilon and calculate macroscopic dielectric constant ===
2095    ! * epsm_lf(w)=1/epsm1(G=0,Gp=0,w).
2096    ! * Since G=Gp=0 there is no difference btw symmetrical and not symmetrical.
2097    !
2098    call xginv(chi0,npwe,comm=comm)
2099 
2100    epsm_lf(iqlwl) = one/chi0(1,1)
2101    eelf(iqlwl) = -AIMAG(chi0(1,1))
2102 
2103    if (prtvol > 0) then
2104      call wrtout(std_out," Symmetrical epsilon^-1(G,G'')",'COLL')
2105      call print_arr(chi0, unit=std_out)
2106    end if
2107    !
2108    ! Save wings of e^-1 overwriting input values.
2109    if (dim_wing>0.and..FALSE.) then
2110      chi0_lwing(:,iqlwl) = chi0(:,1)
2111      chi0_uwing(:,iqlwl) = chi0(1,:)
2112    end if
2113 
2114  end do !iqlwl
2115 
2116  if (allocated(chi0_save))  then
2117    ABI_FREE(chi0_save)
2118  end if
2119 
2120 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

PARENTS

      m_screen

CHILDREN

SOURCE

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

PARENTS

SOURCE

2653 function ylmstar_over_qTq(cart_vers,int_pars,real_pars,cplx_pars)
2654 
2655 
2656 !This section has been created automatically by the script Abilint (TD).
2657 !Do not modify the following lines by hand.
2658 #undef ABI_FUNC
2659 #define ABI_FUNC 'ylmstar_over_qTq'
2660 !End of the abilint section
2661 
2662  implicit none
2663 
2664 !Arguments ------------------------------------
2665 !scalars
2666  real(dp),intent(in) :: cart_vers(3)
2667  integer,intent(in) :: int_pars(:)
2668  real(dp),intent(in) :: real_pars(:)
2669  complex(dpc),intent(in) :: cplx_pars(:)
2670  complex(dpc) :: ylmstar_over_qTq
2671 !arrays
2672 
2673 !Local variables-------------------------------
2674 !scalars
2675  integer :: ll,mm
2676 !arrays
2677  complex(dpc) :: tensor(3,3)
2678 
2679 ! *************************************************************************
2680 
2681  tensor = RESHAPE(cplx_pars(1:9),(/3,3/))
2682  ll = int_pars(1) ! ll starts from zero.
2683  mm = int_pars(2) ! m \in [-l,l]
2684 
2685  ylmstar_over_qTq = CONJG(ylmc(ll,mm,cart_vers))/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers))
2686 
2687  RETURN
2688  ABI_UNUSED(real_pars(1))
2689 
2690 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)

PARENTS

SOURCE

2721 function ylmstar_wtq_over_qTq(cart_vers,int_pars,real_pars,cplx_pars)
2722 
2723 
2724 !This section has been created automatically by the script Abilint (TD).
2725 !Do not modify the following lines by hand.
2726 #undef ABI_FUNC
2727 #define ABI_FUNC 'ylmstar_wtq_over_qTq'
2728 !End of the abilint section
2729 
2730  implicit none
2731 
2732 !Arguments ------------------------------------
2733 !scalars
2734  real(dp),intent(in) :: cart_vers(3)
2735  integer,intent(in) :: int_pars(:)
2736  real(dp),intent(in) :: real_pars(:)
2737  complex(dpc),intent(in) :: cplx_pars(:)
2738  complex(dpc) :: ylmstar_wtq_over_qTq
2739 !arrays
2740 
2741 !Local variables-------------------------------
2742 !scalars
2743  integer :: ll,mm
2744  real(dp) :: wtq
2745 !arrays
2746  real(dp) :: gprimd(3,3),rprimd(3,3),red_vers(3)
2747  complex(dpc) :: tensor(3,3)
2748 
2749 ! *************************************************************************
2750 
2751  MSG_ERROR("Work in progress")
2752  ! box_len has to be tested
2753 
2754  gprimd = RESHAPE(real_pars(1:9),(/3,3/))
2755  red_vers = MATMUL(rprimd,cart_vers)
2756  wtq = box_len(red_vers,gprimd)
2757 
2758  tensor = RESHAPE(cplx_pars(1:9),(/3,3/))
2759  ll = int_pars(1) ! true ll i.e. not shifted
2760  mm = int_pars(2)
2761 
2762  ylmstar_wtq_over_qTq = CONJG(ylmc(ll,mm,cart_vers))*wtq/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers))
2763 
2764 end function ylmstar_wtq_over_qTq