TABLE OF CONTENTS
- ABINIT/m_screening
- m_screening/atddft_hyb_symepsm1
- m_screening/atddft_symepsm1
- m_screening/chi_free
- m_screening/chi_new
- m_screening/chi_t
- m_screening/decompose_epsm1
- m_screening/em1results_free
- m_screening/em1results_print
- m_screening/epsilonm1_results
- m_screening/Epsm1_symmetrizer
- m_screening/Epsm1_symmetrizer_inplace
- m_screening/get_epsm1
- m_screening/init_Er_from_file
- m_screening/lebedev_laikov_int
- m_screening/lwl_free
- m_screening/lwl_init
- m_screening/lwl_t
- m_screening/lwl_write
- m_screening/make_epsm1_driver
- m_screening/mdielf_bechstedt
- m_screening/mkdump_Er
- m_screening/mkem1_q0
- m_screening/rpa_symepsm1
- m_screening/screen_mdielf
- m_screening/ylmstar_over_qTq
- m_screening/ylmstar_wtq_over_qTq
ABINIT/m_screening [ 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-2024 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_screening 24 25 use defs_basis 26 use m_abicore 27 use m_hide_blas 28 use m_linalg_interfaces 29 use m_xmpi 30 use m_errors 31 use m_copy 32 use m_splines 33 use m_lebedev 34 use m_spectra 35 use m_nctk 36 use m_distribfft 37 use netcdf 38 39 use defs_abitypes, only : MPI_type 40 use m_gwdefs, only : GW_TOLQ0, czero_gw, GW_Q0_DEFAULT 41 use m_fstrings, only : toupper, endswith, sjoin, itoa 42 use m_io_tools, only : open_file 43 use m_numeric_tools, only : print_arr, hermitianize 44 use m_special_funcs, only : k_fermi, k_thfermi 45 use m_geometry, only : normv, vdotw, metric 46 use m_hide_lapack, only : xginv 47 use m_crystal, only : crystal_t 48 use m_bz_mesh, only : kmesh_t, box_len 49 use m_fft_mesh, only : g2ifft 50 use m_fftcore, only : kgindex 51 use m_fft, only : fourdp 52 use m_gsphere, only : gsphere_t 53 use m_vcoul, only : vcoul_t 54 use m_io_screening, only : hscr_io, read_screening, write_screening, & 55 HSCR_LATEST_HEADFORM, hscr_t, ncname_from_id, em1_ncname 56 use m_paw_sphharm, only : ylmc 57 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 58 59 implicit none 60 61 private
m_screening/atddft_hyb_symepsm1 [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
atddft_hyb_symepsm1
FUNCTION
Calculate $\tilde\epsilon^{-1}$ using ALDA within TDDFT Based on atddft_symepsm1, modified for the LR+ALDA hybrid kernel Output the electron energy loss function and the macroscopic dielectric function with and without local field effects (only if non-zero real frequencies are available)
INPUTS
iqibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated npwe=Number of G-vectors in chi0. nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case) dim_wing=Dimension of the wings (0 or 3 if q-->0) kxcg_mat_sr=Short-range fxc kernel used in the TE epsilon^-1 option_test=Only for TDDFT: == 0 for TESTPARTICLE == == 1 for TESTELECTRON == Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction %nqibz=Number of q-points. %qibz(3,nqibz)=q-points in the IBZ. comm=MPI communicator. chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0) chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0) chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
OUTPUT
SIDE EFFECTS
chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output the symmetrized inverse dielectric matrix.
SOURCE
2258 subroutine atddft_hyb_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,kxcg_mat_sr,option_test,my_nqlwl,dim_wing,omega,& 2259 & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm) 2260 2261 !Arguments ------------------------------------ 2262 !scalars 2263 integer,intent(in) :: iqibz,nI,nJ,npwe,dim_wing,my_nqlwl 2264 integer,intent(in) :: option_test,comm 2265 type(vcoul_t),target,intent(in) :: Vcp 2266 !arrays 2267 complex(gwpc),intent(in) :: kxcg_mat(npwe,npwe), kxcg_mat_sr(npwe,npwe) 2268 complex(dpc),intent(in) :: omega 2269 complex(dpc),intent(inout) :: chi0_lwing(npwe*nI,dim_wing) 2270 complex(dpc),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing) 2271 complex(dpc),intent(inout) :: chi0_head(dim_wing,dim_wing) 2272 complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ) 2273 real(dp),intent(out) :: eelf(my_nqlwl) 2274 complex(dpc),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl) 2275 2276 !Local variables------------------------------- 2277 !scalars 2278 integer,parameter :: master=0,prtvol=0 2279 integer :: ig1,ig2,my_rank,nprocs,ierr 2280 real(dp) :: ucvol, Zr 2281 logical :: is_qeq0 2282 character(len=500) :: msg 2283 !arrays 2284 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 2285 complex(gwpc),allocatable :: chitmp(:,:) 2286 complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:) 2287 2288 ! ************************************************************************* 2289 2290 ABI_UNUSED(chi0_head(1,1)) 2291 ABI_UNUSED(chi0_lwing(1,1)) 2292 ABI_UNUSED(chi0_uwing(1,1)) 2293 2294 if (nI/=1.or.nJ/=1) then 2295 ABI_ERROR("nI or nJ=/1 not yet implemented") 2296 end if 2297 2298 ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded") 2299 ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded") 2300 2301 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2302 2303 call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol) 2304 2305 is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) 2306 2307 if (iqibz==1) then 2308 !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0 2309 vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! TODO add treatment of non-Analytic behavior 2310 else 2311 vc_sqrt => Vcp%vc_sqrt(:,iqibz) 2312 end if 2313 2314 write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwpc*b2Mb," Mb" 2315 ABI_MALLOC_OR_DIE(chitmp,(npwe,npwe), ierr) 2316 ! 2317 ! * Calculate chi0*fxc. 2318 chitmp = MATMUL(chi0,kxcg_mat) 2319 ! * First calculate the NLF contribution 2320 do ig1=1,npwe 2321 do ig2=1,npwe 2322 chitmp(ig1,ig2)=-chitmp(ig1,ig2) 2323 end do 2324 chitmp(ig1,ig1)=chitmp(ig1,ig1)+one 2325 end do 2326 2327 call xginv(chitmp,npwe,comm=comm) 2328 2329 chitmp = MATMUL(chitmp,chi0) 2330 !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All") 2331 chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1) 2332 chitmp(1,1)=chitmp(1,1)+one 2333 2334 epsm_nlf(1)=chitmp(1,1) 2335 2336 chitmp = MATMUL(chi0,kxcg_mat) 2337 ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp. 2338 do ig1=1,npwe 2339 do ig2=1,npwe 2340 chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2 2341 end do 2342 chitmp(ig1,ig1)=chitmp(ig1,ig1)+one 2343 end do 2344 2345 ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0. 2346 call xginv(chitmp,npwe,comm=comm) 2347 chitmp=MATMUL(chitmp,chi0) 2348 2349 ! * Save result, now chi0 contains chi. 2350 chi0=chitmp 2351 2352 select case (option_test) 2353 case (0) 2354 ! Symmetrized TESTPARTICLE epsilon^-1 2355 call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi') 2356 do ig1=1,npwe 2357 chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:) 2358 chi0(ig1,ig1)=one+chi0(ig1,ig1) 2359 end do 2360 2361 case (1) 2362 ! Symmetrized TESTELECTRON epsilon^-1 2363 call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + Vc*chi + Zr*Kxc_sr*chi') 2364 chitmp=MATMUL(kxcg_mat_sr,chi0) 2365 Zr = 0.78 2366 2367 ! Perform hermitianization, only valid along the imaginary axis. 2368 if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All") 2369 2370 do ig1=1,npwe 2371 chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:)*Zr 2372 chi0(ig1,ig1)=one+chi0(ig1,ig1) 2373 end do 2374 2375 case default 2376 ABI_BUG(sjoin('Wrong option_test:',itoa(option_test))) 2377 end select 2378 2379 ABI_FREE(chitmp) 2380 ! 2381 ! === chi0 now contains symmetrical epsm1 === 2382 ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w). 2383 epsm_lf(1) = one/chi0(1,1) 2384 eelf (1) = -AIMAG(chi0(1,1)) 2385 2386 if (prtvol > 0) then 2387 write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]' 2388 call wrtout(std_out,msg) 2389 call print_arr(chi0,unit=std_out) 2390 end if 2391 2392 end subroutine atddft_hyb_symepsm1
m_screening/atddft_symepsm1 [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
atddft_symepsm1
FUNCTION
Calculate $\tilde\epsilon^{-1}$ using ALDA within TDDFT 2) Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space calculating these quantities for different small q-directions specified by the user (Not yet operative) Output the electron energy loss function and the macroscopic dielectric function with and without local field effects (only if non-zero real frequencies are available)
INPUTS
iqibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated npwe=Number of G-vectors in chi0. nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case) dim_wing=Dimension of the wings (0 or 3 if q-->0) option_test=Only for TDDFT: == 0 for TESTPARTICLE == == 1 for TESTELECTRON == Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction %nqibz=Number of q-points. %qibz(3,nqibz)=q-points in the IBZ. comm=MPI communicator. chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0) chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0) chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0)
OUTPUT
SIDE EFFECTS
chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output the symmetrized inverse dielectric matrix.
SOURCE
2083 subroutine atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0,kxcg_mat,option_test,my_nqlwl,dim_wing,omega,& 2084 & chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm) 2085 2086 !Arguments ------------------------------------ 2087 !scalars 2088 integer,intent(in) :: iqibz,nI,nJ,npwe,dim_wing,my_nqlwl 2089 integer,intent(in) :: option_test,comm 2090 type(vcoul_t),target,intent(in) :: Vcp 2091 !arrays 2092 complex(gwpc),intent(in) :: kxcg_mat(npwe,npwe) 2093 complex(dpc),intent(in) :: omega 2094 complex(dpc),intent(inout) :: chi0_lwing(npwe*nI,dim_wing) 2095 complex(dpc),intent(inout) :: chi0_uwing(npwe*nJ,dim_wing) 2096 complex(dpc),intent(inout) :: chi0_head(dim_wing,dim_wing) 2097 complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ) 2098 real(dp),intent(out) :: eelf(my_nqlwl) 2099 complex(dpc),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl) 2100 2101 !Local variables------------------------------- 2102 !scalars 2103 integer,parameter :: master=0,prtvol=0 2104 integer :: ig1,ig2,my_rank,nprocs,ierr 2105 real(dp) :: ucvol 2106 logical :: is_qeq0 2107 character(len=500) :: msg 2108 !arrays 2109 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 2110 complex(gwpc),allocatable :: chitmp(:,:) 2111 complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:) 2112 2113 ! ************************************************************************* 2114 2115 ABI_UNUSED(chi0_head(1,1)) 2116 ABI_UNUSED(chi0_lwing(1,1)) 2117 ABI_UNUSED(chi0_uwing(1,1)) 2118 2119 if (nI/=1.or.nJ/=1) then 2120 ABI_ERROR("nI or nJ=/1 not yet implemented") 2121 end if 2122 2123 ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded") 2124 ABI_CHECK(my_nqlwl==1,"my_nqlwl/=1 not coded") 2125 2126 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2127 2128 call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol) 2129 2130 is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) 2131 2132 if (iqibz==1) then 2133 !%vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0 2134 vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! TODO add treatment of non-Analytic behavior 2135 else 2136 vc_sqrt => Vcp%vc_sqrt(:,iqibz) 2137 end if 2138 2139 write(msg,'(a,f8.2,a)')" chitmp requires: ",npwe**2*gwpc*b2Mb," Mb" 2140 ABI_MALLOC_OR_DIE(chitmp, (npwe,npwe), ierr) 2141 ! 2142 ! Calculate chi0*fxc. 2143 chitmp = MATMUL(chi0,kxcg_mat) 2144 ! First calculate the NLF contribution 2145 do ig1=1,npwe 2146 do ig2=1,npwe 2147 chitmp(ig1,ig2)=-chitmp(ig1,ig2) 2148 end do 2149 chitmp(ig1,ig1)=chitmp(ig1,ig1)+one 2150 end do 2151 2152 call xginv(chitmp,npwe,comm=comm) 2153 2154 chitmp = MATMUL(chitmp,chi0) 2155 !if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All") 2156 chitmp(1,1)=-vc_sqrt(1)*chitmp(1,1)*vc_sqrt(1) 2157 chitmp(1,1)=chitmp(1,1)+one 2158 2159 epsm_nlf(1)=chitmp(1,1) 2160 2161 chitmp = MATMUL(chi0,kxcg_mat) 2162 ! * Calculate (1-chi0*Vc-chi0*Kxc) and put it in chitmp. 2163 do ig1=1,npwe 2164 do ig2=1,npwe 2165 chitmp(ig1,ig2)=-chitmp(ig1,ig2)-chi0(ig1,ig2)*vc_sqrt(ig2)**2 2166 end do 2167 chitmp(ig1,ig1)=chitmp(ig1,ig1)+one 2168 end do 2169 2170 ! * Invert (1-chi0*Vc-chi0*Kxc) and Multiply by chi0. 2171 call xginv(chitmp,npwe,comm=comm) 2172 chitmp=MATMUL(chitmp,chi0) 2173 2174 ! * Save result, now chi0 contains chi. 2175 chi0=chitmp 2176 2177 select case (option_test) 2178 case (0) 2179 ! Symmetrized TESTPARTICLE epsilon^-1 2180 call wrtout(std_out,' Calculating TESTPARTICLE epsilon^-1(G,G") = 1 + Vc*chi') 2181 do ig1=1,npwe 2182 chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:) 2183 chi0(ig1,ig1)=one+chi0(ig1,ig1) 2184 end do 2185 2186 case (1) 2187 ! Symmetrized TESTELECTRON epsilon^-1 2188 call wrtout(std_out,' Calculating TESTELECTRON epsilon^-1(G,G") = 1 + (Vc + fxc)*chi') 2189 chitmp=MATMUL(kxcg_mat,chi0) 2190 2191 ! Perform hermitianization, only valid along the imaginary axis. 2192 if (.not. ABS(REAL(omega))> tol3) call hermitianize(chitmp,"All") 2193 2194 do ig1=1,npwe 2195 chi0(ig1,:)=(vc_sqrt(ig1)*vc_sqrt(:))*chi0(ig1,:)+chitmp(ig1,:) 2196 chi0(ig1,ig1)=one+chi0(ig1,ig1) 2197 end do 2198 2199 case default 2200 ABI_BUG(sjoin('Wrong option_test:',itoa(option_test))) 2201 end select 2202 2203 ABI_FREE(chitmp) 2204 ! 2205 ! === chi0 now contains symmetrical epsm1 === 2206 ! * Calculate macroscopic dielectric constant epsm_lf(w)=1/epsm1(G=0,Gp=0,w). 2207 epsm_lf(1) = one/chi0(1,1) 2208 eelf (1) = -AIMAG(chi0(1,1)) 2209 2210 if (prtvol > 0) then 2211 write(msg,'(a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at omega',omega*Ha_eV,' [eV]' 2212 call wrtout(std_out,msg) 2213 call print_arr(chi0,unit=std_out) 2214 end if 2215 2216 end subroutine atddft_symepsm1
m_screening/chi_free [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
chi_free
FUNCTION
INPUTS
OUTPUT
SOURCE
3069 subroutine chi_free(chi) 3070 3071 !Arguments ------------------------------------ 3072 !scalars 3073 type(chi_t),intent(inout) :: chi 3074 3075 ! ************************************************************************* 3076 3077 if (allocated(chi%mat)) then 3078 ABI_FREE(chi%mat) 3079 end if 3080 if (allocated(chi%head)) then 3081 ABI_FREE(chi%head) 3082 end if 3083 if (allocated(chi%lwing)) then 3084 ABI_FREE(chi%lwing) 3085 end if 3086 if (allocated(chi%uwing)) then 3087 ABI_FREE(chi%uwing) 3088 end if 3089 3090 end subroutine chi_free
m_screening/chi_new [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
chi_new
FUNCTION
INPUTS
OUTPUT
SOURCE
3037 type(chi_t) function chi_new(npwe, nomega) result(chi) 3038 3039 !Arguments ------------------------------------ 3040 integer,intent(in) :: npwe,nomega 3041 3042 ! ************************************************************************* 3043 3044 chi%nomega = nomega; chi%npwe = npwe 3045 3046 ABI_MALLOC(chi%mat, (npwe,npwe,nomega)) 3047 3048 ABI_MALLOC(chi%head, (3,3,nomega)) 3049 ABI_MALLOC(chi%lwing, (npwe, nomega,3)) 3050 ABI_MALLOC(chi%uwing, (npwe, nomega,3)) 3051 3052 end function chi_new
m_screening/chi_t [ Types ]
[ Top ] [ m_screening ] [ Types ]
NAME
chi_t
FUNCTION
This object contains the head and the wings of the polarizability These quantities are used to treat the q-->0 limit
SOURCE
179 type,public :: chi_t 180 181 integer :: npwe 182 ! Number of G vectors. 183 184 integer :: nomega 185 ! Number of frequencies 186 187 complex(gwpc),allocatable :: mat(:,:,:) 188 ! mat(npwe, npwe, nomega) 189 190 complex(dpc),allocatable :: head(:,:,:) 191 ! head(3,3,nomega) 192 193 complex(dpc),allocatable :: lwing(:,:,:) 194 ! lwing(3,npwe,nomega) 195 ! Lower wings 196 197 complex(dpc),allocatable :: uwing(:,:,:) 198 ! uwing(3,npwe,nomega) 199 ! Upper wings. 200 201 end type chi_t 202 203 public :: chi_new ! Create new object (allocate memory) 204 public :: chi_free ! Free memory.
m_screening/decompose_epsm1 [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
decompose_epsm1
FUNCTION
Decompose the complex symmetrized dielectric
INPUTS
OUTPUT
SOURCE
1179 subroutine decompose_epsm1(Er,iqibz,eigs) 1180 1181 !Arguments ------------------------------------ 1182 !scalars 1183 integer,intent(in) :: iqibz 1184 type(Epsilonm1_results),intent(in) :: Er 1185 !arrays 1186 complex(dpc),intent(out) :: eigs(Er%npwe,Er%nomega) 1187 1188 !Local variables------------------------------- 1189 !scalars 1190 integer :: info,lwork,iw,negw,ig1,ig2,idx,sdim,npwe,ierr 1191 character(len=500) :: msg 1192 !arrays 1193 real(dp),allocatable :: ww(:),rwork(:) 1194 complex(dpc),allocatable :: work(:),Adpp(:),eigvec(:,:),Afull(:,:),vs(:,:),wwc(:) 1195 logical,allocatable :: bwork(:) 1196 logical :: sortcplx !BUG in abilint 1197 1198 ! ********************************************************************* 1199 1200 ABI_CHECK(Er%mqmem/=0,'mqmem==0 not implemented') 1201 1202 npwe = Er%npwe 1203 1204 do iw=1,Er%nomega 1205 1206 if (ABS(REAL(Er%omega(iw)))>0.00001) then 1207 ! Eigenvalues for a generic complex matrix 1208 lwork=4*2*npwe 1209 ABI_MALLOC(wwc,(npwe)) 1210 ABI_MALLOC(work,(lwork)) 1211 ABI_MALLOC(rwork,(npwe)) 1212 ABI_MALLOC(bwork,(npwe)) 1213 ABI_MALLOC(vs,(npwe,npwe)) 1214 ABI_MALLOC(Afull,(npwe,npwe)) 1215 1216 Afull=Er%epsm1(:,:,iw,iqibz) 1217 1218 !for the moment no sort, maybe here I should sort using the real part? 1219 call ZGEES('V','N',sortcplx,npwe,Afull,npwe,sdim,wwc,vs,npwe,work,lwork,rwork,bwork,info) 1220 if (info/=0) then 1221 ABI_ERROR(sjoin("ZGEES returned info:",itoa(info))) 1222 end if 1223 1224 eigs(:,iw)=wwc(:) 1225 1226 ABI_FREE(wwc) 1227 ABI_FREE(work) 1228 ABI_FREE(rwork) 1229 ABI_FREE(bwork) 1230 ABI_FREE(vs) 1231 ABI_FREE(Afull) 1232 1233 else 1234 ! Hermitian version. 1235 lwork=2*npwe-1 1236 ABI_MALLOC(ww,(npwe)) 1237 ABI_MALLOC(work,(lwork)) 1238 ABI_MALLOC(rwork,(3*npwe-2)) 1239 ABI_MALLOC(eigvec,(npwe,npwe)) 1240 ABI_MALLOC_OR_DIE(Adpp,(npwe*(npwe+1)/2), ierr) 1241 1242 idx=0 ! Pack the matrix 1243 do ig2=1,npwe 1244 do ig1=1,ig2 1245 idx=idx+1 1246 Adpp(idx)=Er%epsm1(ig1,ig2,iw,iqibz) 1247 end do 1248 end do 1249 1250 ! For the moment we require also the eigenvectors. 1251 call ZHPEV('V','U',npwe,Adpp,ww,eigvec,npwe,work,rwork,info) 1252 if (info/=0) then 1253 ABI_ERROR(sjoin('ZHPEV returned info=', itoa(info))) 1254 end if 1255 1256 negw=(COUNT((REAL(ww)<tol6))) 1257 if (negw/=0) then 1258 write(msg,'(a,i5,a,i3,a,f8.4)')& 1259 'Found negative eigenvalues. No. ',negw,' at iqibz= ',iqibz,' minval= ',MINVAL(REAL(ww)) 1260 ABI_WARNING(msg) 1261 end if 1262 1263 eigs(:,iw)=ww(:) 1264 1265 ABI_FREE(ww) 1266 ABI_FREE(work) 1267 ABI_FREE(rwork) 1268 ABI_FREE(eigvec) 1269 ABI_FREE(Adpp) 1270 end if 1271 end do !iw 1272 1273 ! contains 1274 ! function sortcplx(carg) result(res) 1275 ! implicit none 1276 ! complex(dpc),intent(in) :: carg 1277 ! logical :: res 1278 ! res=.TRUE. 1279 ! end function sortcplx 1280 1281 end subroutine decompose_epsm1
m_screening/em1results_free [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
em1results_free
FUNCTION
Deallocate all the pointers in Er that result to be associated. Perform also a cleaning of the Header.
INPUTS
OUTPUT
SOURCE
274 subroutine em1results_free(Er) 275 276 !Arguments ------------------------------------ 277 !scalars 278 type(Epsilonm1_results),intent(inout) :: Er 279 ! ************************************************************************* 280 281 !@Epsilonm1_results 282 !integer 283 ABI_SFREE(Er%gvec) 284 285 !real 286 ABI_SFREE(Er%qibz) 287 ABI_SFREE(Er%qlwl) 288 289 !complex 290 ABI_SFREE(Er%epsm1) 291 ABI_SFREE(Er%omega) 292 293 !datatypes 294 call Er%Hscr%free() 295 296 end subroutine em1results_free
m_screening/em1results_print [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
em1results_print
FUNCTION
Print the basic dimensions and the most important quantities reported in the Epsilonm1_results data type.
INPUTS
Er<Epsilonm1_results>=The data type. unit[optional]=the unit number for output. prtvol[optional]=verbosity level. mode_paral[optional]=either COLL or PERS.
OUTPUT
Only printing.
SOURCE
320 subroutine em1results_print(Er,unit,prtvol,mode_paral) 321 322 !Arguments ------------------------------------ 323 integer,optional,intent(in) :: unit,prtvol 324 character(len=4),optional,intent(in) :: mode_paral 325 type(Epsilonm1_results),intent(in) :: Er 326 327 !Local variables------------------------------- 328 integer :: iw,iqibz,iqlwl,unt,my_prtvol 329 character(len=50) :: rfname,rforder,rfapprox,rftest,kxcname 330 character(len=500) :: msg 331 character(len=4) :: mode 332 ! ************************************************************************* 333 334 unt = std_out; if (present(unit)) unt = unit 335 my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol 336 mode ='COLL'; if (present(mode_paral)) mode = mode_paral 337 338 ! === chi0 or \epsilon^{-1} ? === 339 SELECT CASE (Er%ID) 340 CASE (0) 341 rfname = 'Undefined' 342 CASE (1) 343 rfname = 'Irreducible Polarizability' 344 CASE (2) 345 rfname = 'Polarizability' 346 CASE (3) 347 rfname = 'Symmetrical Dielectric Matrix' 348 CASE (4) 349 rfname = 'Symmetrical Inverse Dielectric Matrix' 350 CASE DEFAULT 351 ABI_BUG(sjoin('Wrong Er%ID:',itoa(Er%ID))) 352 END SELECT 353 354 ! For chi, \espilon or \epsilon^{-1}, define the approximation. 355 rfapprox='None' 356 if (Er%ID>=2.or.Er%ID<=4) then 357 if (Er%ikxc==0) then 358 rfapprox='RPA' 359 else if (Er%ikxc>0) then 360 rfapprox='Static TDDFT' 361 else 362 rfapprox='TDDFT' 363 end if 364 end if 365 366 ! === If TDDFT and \epsilon^{-1}, define the type === 367 rftest='None' 368 ! if (Er%ID==0) then 369 ! if (Er%test_type==0) then 370 ! rftest='TEST-PARTICLE' 371 ! else if (Er%test_type==1) then 372 ! rftest='TEST-ELECTRON' 373 ! else 374 ! write(msg,'(4a,i3)')ch10,& 375 !& ' em1results_print : BUG - ',ch10,& 376 !& ' Wrong value of Er%test_type = ',Er%test_type 377 ! ABI_ERROR(msg) 378 ! end if 379 ! end if 380 381 ! === Define time-ordering === 382 rforder='Undefined' 383 if (Er%Tordering==1) then 384 rforder='Time-Ordered' 385 else if (Er%Tordering==2) then 386 rforder='Advanced' 387 else if (Er%Tordering==3) then 388 rforder='Retarded' 389 else 390 ABI_BUG(sjoin('Wrong er%tordering= ',itoa(Er%Tordering))) 391 end if 392 393 kxcname='None' 394 if (Er%ikxc/=0) then 395 !TODO Add function to retrieve kxc name 396 ABI_ERROR('Add function to retrieve kxc name') 397 kxcname='XXXXX' 398 end if 399 400 write(msg,'(6a,5(3a))')ch10,& 401 & ' ==== Info on the Response Function ==== ',ch10,& 402 & ' Associated File ................ ',TRIM(Er%fname),ch10,& 403 & ' Response Function Type .......... ',TRIM(rfname),ch10,& 404 & ' Type of Approximation ........... ',TRIM(rfapprox),ch10,& 405 & ' XC kernel used .................. ',TRIM(kxcname),ch10,& 406 & ' Type of probing particle ........ ',TRIM(rftest),ch10,& 407 & ' Time-Ordering ................... ',TRIM(rforder),ch10 408 call wrtout(unt,msg,mode) 409 write(msg,'(a,2i4,a,3(a,i4,a),a,3i4,2a,i4,a)')& 410 ' Number of components ............ ',Er%nI,Er%nJ,ch10,& 411 ' Number of q-points in the IBZ ... ',Er%nqibz,ch10,& 412 ' Number of q-points for q-->0 .... ',Er%nqlwl,ch10,& 413 ' Number of G-vectors ............. ',Er%npwe,ch10,& 414 ' Number of frequencies ........... ',Er%nomega,Er%nomega_r,Er%nomega_i,ch10,& 415 ' Value of mqmem .................. ',Er%mqmem,ch10 416 call wrtout(unt,msg,mode) 417 418 if (Er%nqlwl/=0) then 419 write(msg,'(a,i3)')' q-points for long wavelength limit: ',Er%nqlwl 420 call wrtout(unt,msg,mode) 421 do iqlwl=1,Er%nqlwl 422 write(msg,'(1x,i5,a,3es16.8)')iqlwl,') ',Er%qlwl(:,iqlwl) 423 call wrtout(unt,msg,mode) 424 end do 425 end if 426 427 if (my_prtvol>0) then 428 ! Print out head and wings in the long-wavelength limit. 429 ! TODO add additional stuff. 430 write(msg,'(a,i4)')' Calculated Frequencies: ',Er%nomega 431 call wrtout(unt,msg,mode) 432 do iw=1,Er%nomega 433 write(msg,'(i4,es14.6)')iw,Er%omega(iw)*Ha_eV 434 call wrtout(unt,msg,mode) 435 end do 436 437 write(msg,'(a,i4)')' Calculated q-points: ',Er%nqibz 438 call wrtout(unt,msg,mode) 439 do iqibz=1,Er%nqibz 440 write(msg,'(1x,i4,a,3es16.8)')iqibz,') ',Er%qibz(:,iqibz) 441 call wrtout(unt,msg,mode) 442 end do 443 end if ! my_prtvol>0 444 445 end subroutine em1results_print
m_screening/epsilonm1_results [ Types ]
[ Top ] [ m_screening ] [ Types ]
NAME
epsilonm1_results
FUNCTION
For the GW part of ABINIT, the epsilonm1_results structured datatype gather the results of screening: the inverse dielectric matrix, and the omega matrices.
SOURCE
76 type,public :: Epsilonm1_results 77 78 integer :: id 79 ! Matrix identifier: O if not yet defined, 1 for chi0, 80 ! 2 for chi, 3 for epsilon, 4 for espilon^{-1}, 5 for W. 81 82 integer :: ikxc 83 ! Kxc kernel used, 0 for None (RPA), >0 for static TDDFT (=ixc), <0 for TDDFT 84 85 integer :: fform 86 ! File format: 1002 for SCR|SUSC files. 87 88 integer :: mqmem 89 ! =0 for out-of-core solution, =nqibz if entire matrix is stored in memory. 90 91 integer :: nI,nJ 92 ! Number of components (rows,columns) in chi|eps^-1. (1,1) if collinear. 93 94 integer :: nqibz 95 ! Number of q-points in the IBZ used. 96 97 integer :: nqlwl 98 ! Number of point used for the treatment of the long wave-length limit. 99 100 integer :: nomega 101 ! Total number of frequencies. 102 103 integer :: nomega_i 104 ! Number of purely imaginary frequencies used. 105 106 integer :: nomega_r 107 ! Number of real frequencies used. 108 109 integer :: npwe 110 ! Number of G vectors. 111 112 integer :: test_type 113 ! 0 for None, 1 for TEST-PARTICLE, 2 for TEST-ELECTRON (only for TDDFT) 114 115 integer :: tordering 116 ! 0 if not defined, 1 for Time-Ordered, 2 for Advanced, 3 for Retarded. 117 118 character(len=fnlen) :: fname 119 ! Name of the file from which epsm1 is read. 120 121 integer,allocatable :: gvec(:,:) 122 ! gvec(3,npwe) 123 ! G-vectors used to describe the two-point function (r.l.u.). 124 125 real(dp),allocatable :: qibz(:,:) 126 ! qibz(3,nqibz) 127 ! q-points in reduced coordinates 128 129 real(dp),allocatable :: qlwl(:,:) 130 ! qlwl(3,nqlwl) 131 ! q-points used for the long wave-length limit treatment. 132 133 !real(gwp),allocatable :: epsm1_pole(:,:,:,:) 134 ! epsm1(npwe,npwe,nomega,nqibz) 135 ! Contains the two-point function $\epsilon_{G,Gp}(q,omega)$ in frequency and reciprocal space. 136 137 complex(gwpc),allocatable :: epsm1(:,:,:,:) 138 ! epsm1(npwe,npwe,nomega,nqibz) 139 ! Contains the two-point function $\epsilon_{G,Gp}(q,omega)$ in frequency and reciprocal space. 140 141 complex(dpc),allocatable :: omega(:) 142 ! omega(nomega) 143 ! Frequencies used both along the real and the imaginary axis. 144 145 type(hscr_t) :: Hscr 146 ! The header reported in the _SCR of _SUSC file. 147 ! This object contains information on the susceptibility or the inverse dielectric matrix 148 ! as stored in the external file. These quantities do *NOT* correspond to the quantities 149 ! used during the GW calculation since some parameters might differ, actually they might be smaller. 150 ! For example, the number of G-vectors used can be smaller than the number of G"s stored on file. 151 152 end type Epsilonm1_results 153 154 public :: em1results_free ! Free memory 155 public :: em1results_print ! Print basic info 156 public :: Epsm1_symmetrizer ! Symmetrize two-point function at a q-point in the BZ. 157 public :: Epsm1_symmetrizer_inplace ! In-place version of the above 158 public :: init_Er_from_file ! Initialize the object from file 159 public :: mkdump_Er ! Dump the object to a file. 160 public :: get_epsm1 161 public :: decompose_epsm1 162 public :: make_epsm1_driver ! Calculate the inverse symmetrical dielectric matrix starting from chi0 163 public :: mkem1_q0 ! construct the microscopic dielectric matrix for q-->0 164 public :: screen_mdielf ! Calculates W_{G,G'}(q,w) for a given q-point in the BZ using a model dielectric function. 165 public :: rpa_symepsm1
m_screening/Epsm1_symmetrizer [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
Epsm1_symmetrizer
FUNCTION
Symmetrize the inverse dielectric matrix, namely calculate epsilon^{-1} at a generic q-point in the BZ starting from the knowledge of the matrix at a q-point in the IBZ. The procedure is quite generic and can be used for every two-point function which has the same symmetry as the crystal.
INPUTS
nomega=Number of frequencies required. All frequencies from 1 up to nomega are symmetrized. npwc=Number of G vectors in symmetrized matrix, has to be smaller than Er%npwe. remove_exchange=If .TRUE., return e^{-1}-1 namely remove the exchange part. Er<Epsilonm1_results>=Data structure containing the inverse dielectric matrix. Gsph<gsphere_t>=data related to the G-sphere %grottb %phmSGt Qmesh<kmesh_t>=Structure defining the q-mesh used for Er. %nbz=Number of q-points in the BZ %tab(nbz)=Index of the symmetric q-point in the IBZ, for each point in the BZ %tabo(nbz)=The operation that rotates q_ibz onto \pm q_bz (depending on tabi) %tabi(nbz)=-1 if time-reversal has to be considered, 1 otherwise iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.
OUTPUT
epsm1_qbz(npwc,npwc,nomega)=The inverse dielectric matrix at the q-point defined by iq_bz. Exchange part can be subtracted out.
NOTES
In the present implementation we are not considering a possible umklapp vector G0 in the expression Sq = q+G0. Treating this case would require some changes in the G-sphere since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required to reconstruct the BZ. * Remember the symmetry properties of \tilde\espilon^{-1} If q_bz = S q_ibz + G0: $\epsilon^{-1}_{SG1-G0, SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau} \epsilon^{-1}_{G1, G2)}(q) If time-reversal symmetry can be used then: $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau} \epsilon^{-1}_{-S^{-1}(G1+Go), -S^{-1}(G2+G0)}^*(q)
TODO
Symmetrization can be skipped if iq_bz correspond to a point in the IBZ
SOURCE
500 subroutine Epsm1_symmetrizer(iq_bz,nomega,npwc,Er,Gsph,Qmesh,remove_exchange,epsm1_qbz) 501 502 !Arguments ------------------------------------ 503 !scalars 504 integer,intent(in) :: iq_bz,nomega,npwc 505 logical,intent(in) :: remove_exchange 506 type(Epsilonm1_results),intent(in) :: Er 507 type(gsphere_t),target,intent(in) :: Gsph 508 type(kmesh_t),intent(in) :: Qmesh 509 !arrays 510 complex(gwpc),intent(out) :: epsm1_qbz(npwc,npwc,nomega) 511 512 !Local variables------------------------------- 513 !scalars 514 integer :: iw,ii,jj,iq_ibz,itim_q,isym_q,iq_loc,sg1,sg2 515 complex(gwpc) :: phmsg1t,phmsg2t_star 516 !arrays 517 real(dp) :: qbz(3) 518 519 ! ********************************************************************* 520 521 ABI_CHECK(Er%nomega >= nomega, 'Too many frequencies required') 522 ABI_CHECK(Er%npwe >= npwc, 'Too many G-vectors required') 523 524 ! Get iq_ibz, and symmetries from iq_ibz. 525 call Qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym_q, itim_q) 526 527 ! If out-of-memory, only Er%espm1(:,:,:,1) has been allocated and filled. 528 iq_loc=iq_ibz; if (Er%mqmem==0) iq_loc=1 529 530 ! MG: grottb is a 1-1 mapping, hence we can collapse the loops (false sharing is not an issue here). 531 !grottb => Gsph%rottb (1:npwc,itim_q,isym_q) 532 !phmSgt => Gsph%phmSGt(1:npwc,isym_q) 533 534 !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(sg2,sg1,phmsg1t,phmsg2t_star) 535 do iw=1,nomega 536 do jj=1,npwc 537 sg2 = Gsph%rottb(jj, itim_q, isym_q) 538 phmsg2t_star = CONJG(Gsph%phmSGt(jj, isym_q)) 539 do ii=1,npwc 540 sg1 = Gsph%rottb(ii,itim_q,isym_q) 541 phmsg1t = Gsph%phmSGt(ii,isym_q) 542 epsm1_qbz(sg1,sg2,iw) = Er%epsm1(ii,jj,iw,iq_loc) * phmsg1t * phmsg2t_star 543 !epsm1_qbz(sg1,sg2,iw) = Er%epsm1(ii,jj,iw,iq_loc) * phmSgt(ii) * CONJG(phmSgt(jj)) 544 end do 545 end do 546 end do 547 ! 548 ! === Account for time-reversal === 549 !epsm1_qbz(:,:,iw)=TRANSPOSE(epsm1_qbz(:,:,iw)) 550 if (itim_q==2) then 551 !$OMP PARALLEL DO IF (nomega > 1) 552 do iw=1,nomega 553 call sqmat_itranspose(npwc,epsm1_qbz(:,:,iw)) 554 end do 555 end if 556 557 if (remove_exchange) then 558 ! === Subtract the exchange contribution === 559 ! If it's a pole screening, the exchange contribution is already removed 560 !$OMP PARALLEL DO IF (nomega > 1) 561 do iw=1,nomega 562 do ii=1,npwc 563 epsm1_qbz(ii,ii,iw)=epsm1_qbz(ii,ii,iw)-CMPLX(1.0_gwp,0.0_gwp) 564 end do 565 end do 566 endif 567 568 end subroutine Epsm1_symmetrizer
m_screening/Epsm1_symmetrizer_inplace [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
Epsm1_symmetrizer_inplace
FUNCTION
Same function as Epsm1_symmetrizer, ecxept now the array Ep%epsm1 is modified inplace thorugh an auxiliary work array of dimension (npwc,npwc)
INPUTS
nomega=Number of frequencies required. All frequencies from 1 up to nomega are symmetrized. npwc=Number of G vectors in symmetrized matrix, has to be smaller than Er%npwe. remove_exchange=If .TRUE., return e^{-1}-1 namely remove the exchange part. Er<Epsilonm1_results>=Data structure containing the inverse dielectric matrix. Gsph<gsphere_t>=data related to the G-sphere Er<Epsilonm1_results>=Data structure containing the inverse dielectric matrix. Gsph<gsphere_t>=data related to the G-sphere %grottb %phmSGt Qmesh<kmesh_t>=Structure defining the q-mesh used for Er. %nbz=Number of q-points in the BZ %tab(nbz)=Index of the symmetric q-point in the IBZ, for each point in the BZ %tabo(nbz)=The operation that rotates q_ibz onto \pm q_bz (depending on tabi) %tabi(nbz)=-1 if time-reversal has to be considered, 1 otherwise iq_bz=Index of the q-point in the BZ where epsilon^-1 is required.
OUTPUT
Er%epsm1(npwc,npwc,nomega,iq_loc) symmetrised
NOTES
In the present implementation we are not considering a possible umklapp vector G0 in the expression Sq = q+G0. Treating this case would require some changes in the G-sphere since we have to consider G-G0. The code however stops in sigma if a nonzero G0 is required to reconstruct the BZ. * Remember the symmetry properties of \tilde\espilon^{-1} If q_bz=Sq_ibz+G0: $\epsilon^{-1}_{SG1-G0,SG2-G0}(q_bz) = e^{+iS(G2-G1).\tau}\epsilon^{-1}_{G1,G2)}(q) If time-reversal symmetry can be used then : $\epsilon^{-1}_{G1,G2}(-q_bz) = e^{+i(G1-G2).\tau}\epsilon^{-1}_{-S^{-1}(G1+Go),-S^{-1}(G2+G0)}^*(q)
TODO
Symmetrization can be skipped if iq_bz correspond to a point in the IBZ
SOURCE
620 subroutine Epsm1_symmetrizer_inplace(iq_bz,nomega,npwc,Er,Gsph,Qmesh,remove_exchange) 621 622 !Arguments ------------------------------------ 623 !scalars 624 integer,intent(in) :: iq_bz,nomega,npwc 625 logical,intent(in) :: remove_exchange 626 type(Epsilonm1_results) :: Er 627 type(gsphere_t),target,intent(in) :: Gsph 628 type(kmesh_t),intent(in) :: Qmesh 629 630 !Local variables------------------------------- 631 !scalars 632 integer :: iw,ii,jj,iq_ibz,itim_q,isym_q,iq_loc,sg1,sg2 633 !arrays 634 real(dp) :: qbz(3) 635 complex(gwpc) :: phmsg1t,phmsg2t_star 636 complex(gwpc),allocatable :: work(:,:) 637 638 ! ********************************************************************* 639 640 ABI_CHECK(Er%nomega>=nomega,'Too many frequencies required') 641 ABI_CHECK(Er%npwe >=npwc , 'Too many G-vectors required') 642 643 ABI_MALLOC(work,(npwc,npwc)) 644 645 ! * Get iq_ibz, and symmetries from iq_ibz. 646 call qmesh%get_BZ_item(iq_bz,qbz,iq_ibz,isym_q,itim_q) 647 648 ! If out-of-memory, only Er%espm1(:,:,:,1) has been allocated and filled. 649 iq_loc=iq_ibz; if (Er%mqmem==0) iq_loc=1 650 651 !$OMP PARALLEL DO PRIVATE(sg2,sg1,phmsg1t,phmsg2t_star) IF (nomega > 1) 652 do iw=1,nomega 653 do jj=1,npwc 654 sg2 = Gsph%rottb(jj,itim_q,isym_q) 655 phmsg2t_star = CONJG(Gsph%phmSGt(jj,isym_q)) 656 do ii=1,npwc 657 sg1 = Gsph%rottb(ii,itim_q,isym_q) 658 phmsg1t = Gsph%phmSGt(ii,isym_q) 659 work(sg1,sg2) = Er%epsm1(ii,jj,iw,iq_loc) * phmsg1t * phmsg2t_star 660 !work(grottb(ii),grottb(jj))=Er%epsm1(ii,jj,iw,iq_loc)*phmSgt(ii)*CONJG(phmSgt(jj)) 661 end do 662 end do 663 Er%epsm1(:,:,iw,iq_loc) = work(:,:) 664 end do 665 ! 666 ! === Account for time-reversal === 667 !Er%epsm1(:,:,iw,iq_loc)=TRANSPOSE(Er%epsm1(:,:,iw,iq_loc)) 668 if (itim_q==2) then 669 !$OMP PARALLEL DO IF (nomega > 1) 670 do iw=1,nomega 671 call sqmat_itranspose(npwc,Er%epsm1(:,:,iw,iq_loc)) 672 end do 673 end if 674 675 ! === Subtract the exchange contribution === 676 if (remove_exchange) then 677 !$OMP PARALLEL DO IF (nomega > 1) 678 do iw=1,nomega 679 do ii=1,npwc 680 Er%epsm1(ii,ii,iw,iq_loc)=Er%epsm1(ii,ii,iw,iq_loc)-1.0_gwp 681 end do 682 end do 683 endif 684 685 ABI_FREE(work) 686 687 end subroutine Epsm1_symmetrizer_inplace
m_screening/get_epsm1 [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
get_epsm1
FUNCTION
Work in progress but the main is idea is as follows: Return the symmetrized inverse dielectric matrix. This method implements both in-core and the out-of-core solution In the later, epsilon^-1 or chi0 are read from file. It is possible to specify options to retrieve (RPA |TDDDT, [TESTCHARGE|TESTPARTICLE]). All dimensions are already initialized in the Er% object, this method should act as a wrapper around rdscr and make_epsm1_driver. A better implementation will be done in the following once the coding of file handlers is completed.
INPUTS
Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction iqibzA[optional]=Index of the q-point to be read from file (only for out-of-memory solutions) iomode=option definig the file format. option_test comm=MPI communicator.
OUTPUT
Er%epsm1
TODO
Remove this routine. Now everything should be done with mkdump_Er
SOURCE
1108 subroutine get_epsm1(Er,Vcp,approx_type,option_test,iomode,comm,iqibzA) 1109 1110 !Arguments ------------------------------------ 1111 !scalars 1112 integer,intent(in) :: iomode,option_test,approx_type,comm 1113 integer,optional,intent(in) :: iqibzA 1114 type(vcoul_t),intent(in) :: Vcp 1115 type(Epsilonm1_results),intent(inout) :: Er 1116 1117 !Local variables------------------------------- 1118 !scalars 1119 integer :: my_approx_type,my_option_test,ng,ierr 1120 ! ********************************************************************* 1121 1122 DBG_ENTER("COLL") 1123 1124 my_approx_type = approx_type; my_option_test = option_test 1125 1126 ! Vcp not yet used. 1127 ng = Vcp%ng 1128 1129 select case (Er%mqmem) 1130 case (0) 1131 ! Out-of-core solution 1132 if (allocated(Er%epsm1)) then 1133 ABI_FREE(Er%epsm1) 1134 end if 1135 ABI_MALLOC_OR_DIE(Er%epsm1,(Er%npwe,Er%npwe,Er%nomega,1), ierr) 1136 1137 ! FIXME there's a problem with SUSC files and MPI-IO 1138 !if (iomode == IO_MODE_MPI) then 1139 ! !write(std_out,*)"read_screening with iomode",iomode,"file: ",trim(er%fname) 1140 ! ABI_WARNING("SUSC files is buggy. Using Fortran IO") 1141 ! call read_screening(em1_ncname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm,iqiA=iqibzA) 1142 !else 1143 call read_screening(em1_ncname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm,iqiA=iqibzA) 1144 !end if 1145 1146 if (Er%id == 4) then 1147 ! If q-slice of epsilon^-1 has been read then return 1148 !call em1results_print(Er) 1149 return 1150 else 1151 ABI_ERROR(sjoin('Wrong Er%ID', itoa(er%id))) 1152 end if 1153 1154 case default 1155 ! In-core solution. 1156 ABI_ERROR("you should not be here") 1157 end select 1158 1159 DBG_EXIT("COLL") 1160 1161 end subroutine get_epsm1
m_screening/init_Er_from_file [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
init_Er_from_file
FUNCTION
Initialize basic dimensions and the important (small) arrays in an Epsilonm1_results data type starting from a file containing either epsilon^{-1} (_SCR) or chi0 (_SUSC).
INPUTS
fname=The name of the external file used to read the matrix. mqmem=0 for out-of-core solution, /=0 if entire matrix has to be stored in memory. npwe_asked=Number of G-vector to be used in the calculation, if <=0 use Max allowed number. comm=MPI communicator.
OUTPUT
Er<Epsilonm1_results>=The structure initialized with basic dimensions and arrays.
SOURCE
711 subroutine init_Er_from_file(Er,fname,mqmem,npwe_asked,comm) 712 713 !Arguments ------------------------------------ 714 integer,intent(in) :: mqmem,npwe_asked,comm 715 character(len=*),intent(in) :: fname 716 type(Epsilonm1_results),intent(inout) :: Er 717 718 !Local variables------------------------------- 719 !scalars 720 integer,parameter :: master=0 721 integer :: iw,fform,my_rank,unclassified 722 real(dp) :: re, im, tol 723 character(len=500) :: msg 724 725 ! ********************************************************************* 726 727 DBG_ENTER("COLL") 728 729 !@Epsilonm1_results 730 my_rank = xmpi_comm_rank(comm) 731 732 ! Read header from file. 733 call wrtout(std_out,sjoin('init_Er_from_file- testing file: ',fname)) 734 call Er%hscr%from_file(fname, fform, comm) 735 736 ! Master echoes the header. 737 if (my_rank==master) call er%hscr%print() 738 739 ! Generic Info 740 Er%ID =0 ! Not yet initialized as epsm1 is calculated in mkdump_Er.F90 741 Er%fname =fname 742 Er%fform =fform 743 Er%Tordering=Er%Hscr%Tordering 744 745 !TODO these quantitities should be checked and initiliazed in mkdump_Er 746 !BEGIN HARCODED 747 Er%nI = 1 748 Er%nJ = 1 749 Er%ikxc = 0 750 Er%test_type=-1 751 752 Er%Hscr%headform=HSCR_LATEST_HEADFORM ! XG20090912 753 !END HARDCODED 754 755 Er%nqibz=Er%Hscr%nqibz 756 Er%mqmem=mqmem ; if (mqmem/=0) Er%mqmem=Er%nqibz 757 ABI_MALLOC(Er%qibz,(3,Er%nqibz)) 758 Er%qibz(:,:)=Er%Hscr%qibz(:,:) 759 760 Er%nqlwl=Er%Hscr%nqlwl 761 ABI_MALLOC(Er%qlwl,(3,Er%nqlwl)) 762 Er%qlwl(:,:)=Er%Hscr%qlwl(:,:) 763 764 Er%nomega=Er%Hscr%nomega 765 ABI_MALLOC(Er%omega,(Er%nomega)) 766 Er%omega(:)=Er%Hscr%omega(:) 767 768 ! Count number of real, imaginary, and complex frequencies. 769 Er%nomega_r = 1 770 Er%nomega_i = 0 771 if (Er%nomega == 2) then 772 Er%nomega_i = 1 773 else 774 unclassified = 0 775 tol = tol6*Ha_eV 776 do iw = 2, Er%nomega 777 re = REAL(Er%omega(iw)) 778 im = AIMAG(Er%omega(iw)) 779 if (re > tol .and. im < tol) then 780 Er%nomega_r = iw ! Real freqs are packed in the first locations. 781 else if (re < tol .and. im > tol) then 782 Er%nomega_i = Er%nomega_i + 1 783 else 784 unclassified = unclassified + 1 785 end if 786 end do 787 if (unclassified > 0) then 788 write(msg,'(3a,i6)')& 789 & 'Some complex frequencies are too small to qualify as real or imaginary.',ch10,& 790 & 'Number of unidentified frequencies = ', unclassified 791 ABI_WARNING(msg) 792 end if 793 end if 794 795 ! Get G-vectors. 796 Er%npwe=Er%Hscr%npwe 797 if (npwe_asked>0) then 798 if (npwe_asked>Er%Hscr%npwe) then 799 write(msg,'(a,i8,2a,i8)')& 800 & 'Number of G-vectors saved on file is less than the value required = ',npwe_asked,ch10,& 801 & 'Calculation will proceed with Max available npwe = ',Er%Hscr%npwe 802 ABI_WARNING(msg) 803 else ! Redefine the no. of G"s for W. 804 Er%npwe=npwe_asked 805 end if 806 end if 807 808 ! pointer to Er%Hscr%gvec ? 809 ABI_MALLOC(Er%gvec,(3,Er%npwe)) 810 Er%gvec=Er%Hscr%gvec(:,1:Er%npwe) 811 812 DBG_EXIT("COLL") 813 814 end subroutine init_Er_from_file
m_screening/lebedev_laikov_int [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
lebedev_laikov_int
FUNCTION
INPUTS
OUTPUT
SOURCE
2589 subroutine lebedev_laikov_int() 2590 2591 !Arguments ------------------------------------ 2592 2593 !Local variables------------------------------- 2594 !scalars 2595 integer :: on,npts,ii,ll,mm,lmax,leb_idx !ierr, 2596 real(dp) :: accuracy 2597 complex(dpc) :: ang_int 2598 !arrays 2599 real(dp) :: cart_vpt(3) !,real_pars(0) 2600 real(dp),allocatable :: vx(:),vy(:),vz(:),ww(:) 2601 complex(dpc) :: tensor(3,3),cplx_pars(9) 2602 complex(dpc),allocatable :: ref_func(:),expd_func(:) !tmp_momenta(:) 2603 2604 ! ************************************************************************* 2605 2606 ABI_ERROR("lebedev_laikov_int is still under development") 2607 2608 !tensor=RESHAPE((/4.0,2.0,4.0,0.5,2.1,0.0,5.4,2.1,5.0/),(/3,3/)) 2609 tensor=RESHAPE((/4.0,0.0,0.0,0.0,4.0,0.0,0.0,0.0,5.0/),(/3,3/)) 2610 !tensor=RESHAPE((/1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0/),(/3,3/)) 2611 2612 npts=26 2613 ABI_MALLOC(vx,(npts)) 2614 ABI_MALLOC(vy,(npts)) 2615 ABI_MALLOC(vz,(npts)) 2616 ABI_MALLOC(ww,(npts)) 2617 2618 !call LD0026(vx,vy,vz,ww,on) 2619 2620 ang_int=czero 2621 do ii=1,npts 2622 cart_vpt = [vx(ii),vy(ii),vz(ii)] 2623 ang_int = ang_int + ww(ii)*DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt)) 2624 end do 2625 2626 !write(std_out,*)"quadratic form associated to tensor=",tensor 2627 write(std_out,*)"on ang_int",on,ang_int 2628 2629 ABI_FREE(vx) 2630 ABI_FREE(vy) 2631 ABI_FREE(vz) 2632 ABI_FREE(ww) 2633 2634 !call init_lebedev_gridset() 2635 cplx_pars = RESHAPE(tensor,(/9/)); accuracy=tol10 2636 2637 ! This is the function to be expanded evaluated on the lebedev_laikov grid of index leb_idx 2638 leb_idx=3; npts=lebedev_npts(leb_idx) 2639 ABI_MALLOC(ref_func,(npts)) 2640 do ii=1,npts 2641 !cart_vpt = Lgridset(leb_idx)%versor(:,ii) 2642 ref_func(ii) = one/DOT_PRODUCT(cart_vpt,MATMUL(tensor,cart_vpt)) 2643 end do 2644 2645 ! Calculate the expansion in angular momenta of 1/{q.Tq}. 2646 ! Only even l-components contribute thanks to the parity of the integrand. 2647 ! tol6 seems to be an acceptable error, convergence wrt lmax is very slow even for simple tensors. 2648 ABI_MALLOC(expd_func,(npts)) 2649 expd_func=czero 2650 lmax=10 2651 do ll=0,lmax,2 2652 !allocate(tmp_momenta(-ll:ll)) 2653 do mm=-ll,ll 2654 ! MG: Commented becase it causes problems with the new version of abilint 2655 !call lebedev_quadrature(ylmstar_over_qTq,(/ll,mm/),real_pars,cplx_pars,ang_int,ierr,accuracy) 2656 write(std_out,*)ll,mm,ang_int 2657 !tmp_momenta(mm) = ang_int 2658 do ii=1,npts 2659 !cart_vpt = Lgridset(leb_idx)%versor(:,ii) 2660 expd_func(ii) = expd_func(ii) + four_pi*ang_int*ylmc(ll,mm,cart_vpt) 2661 end do 2662 end do 2663 !deallocate(tmp_momenta) 2664 write(std_out,*)"Error in angular expansion at l=",ll," is ",MAXVAL(ABS(expd_func-ref_func)) 2665 end do 2666 2667 !BEGINDEBUG 2668 ! do ii=1,npts 2669 ! write(777,*)ref_func(ii) 2670 ! write(778,*)expd_func(ii) 2671 ! end do 2672 !ENDDEBUG 2673 2674 ABI_FREE(expd_func) 2675 ABI_FREE(ref_func) 2676 2677 ABI_ERROR("Exiting from lebedev_laikov_int") 2678 2679 end subroutine lebedev_laikov_int
m_screening/lwl_free [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
lwl_free
FUNCTION
INPUTS
OUTPUT
SOURCE
3333 subroutine lwl_free(lwl) 3334 3335 !Arguments ------------------------------------ 3336 !scalars 3337 type(lwl_t),intent(inout) :: lwl 3338 3339 ! ************************************************************************* 3340 3341 if (allocated(lwl%head)) then 3342 ABI_FREE(lwl%head) 3343 end if 3344 if (allocated(lwl%lwing)) then 3345 ABI_FREE(lwl%lwing) 3346 end if 3347 if (allocated(lwl%uwing)) then 3348 ABI_FREE(lwl%uwing) 3349 end if 3350 if (allocated(lwl%body)) then 3351 ABI_FREE(lwl%body) 3352 end if 3353 3354 end subroutine lwl_free
m_screening/lwl_init [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
lwl_init
FUNCTION
INPUTS
OUTPUT
SOURCE
3265 subroutine lwl_init(lwl, path, method, cryst, vcp, npwe, gvec, comm) 3266 3267 !Arguments ------------------------------------ 3268 !scalars 3269 integer,intent(in) :: comm,method,npwe 3270 character(len=*),intent(in) :: path 3271 type(crystal_t),intent(in) :: cryst 3272 type(vcoul_t),intent(in) :: vcp 3273 type(lwl_t),intent(out) :: lwl 3274 !arrays 3275 integer,intent(in) :: gvec(3,npwe) 3276 3277 !Local variables------------------------------- 3278 !scalars 3279 integer,parameter :: master=0 3280 integer :: iomode,my_rank,nproc,unt 3281 character(len=500) :: msg 3282 !arrays 3283 3284 ! ************************************************************************* 3285 3286 ABI_UNUSED((/cryst%natom, gvec(1,1), vcp%ng/)) 3287 3288 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 3289 3290 lwl%fname = path 3291 lwl%method = method 3292 ABI_CHECK(any(method == [1,2,3]), sjoin("Wrong method:", itoa(method))) 3293 3294 iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF 3295 3296 ! Only master reads. 3297 if (my_rank == master) then 3298 3299 select case (iomode) 3300 case (IO_MODE_FORTRAN) 3301 if (open_file(path, msg, newunit=unt, action="read", form="unformatted", status="old") /= 0) then 3302 ABI_ERROR(msg) 3303 end if 3304 3305 close(unt) 3306 3307 case default 3308 ABI_ERROR(sjoin("iomode:", itoa(iomode), "is not coded")) 3309 end select 3310 end if 3311 3312 ! Broad cast data 3313 if (nproc > 1) then 3314 end if 3315 3316 end subroutine lwl_init
m_screening/lwl_t [ Types ]
[ Top ] [ m_screening ] [ Types ]
NAME
lwl_t
FUNCTION
SOURCE
215 type,public :: lwl_t 216 217 integer :: npwe 218 ! Number of G vectors. 219 220 integer :: nomega 221 ! Number of frequencies 222 223 integer :: method 224 ! 1 = Only head 225 ! 2 = head + wings 226 ! 3 = head + wings + body corrections. 227 228 character(len=fnlen) :: fname 229 ! Name of the file from which epsm1 is read. 230 231 complex(dpc),allocatable :: head(:,:,:) 232 ! head(3,3,nomega) 233 234 complex(dpc),allocatable :: lwing(:,:,:) 235 ! lwing(3,npwe,nomega) 236 ! Lower wings 237 238 complex(dpc),allocatable :: uwing(:,:,:) 239 ! uwing(3,npwe,nomega) 240 ! Upper wings. 241 242 complex(dpc),allocatable :: body(:,:,:) 243 ! uwing(npwe,npwe,nomega) 244 ! Body terms 245 246 end type lwl_t 247 248 public :: lwl_write 249 public :: lwl_init 250 !public :: lwl_from_file 251 public :: lwl_free
m_screening/lwl_write [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
lwl_write
FUNCTION
INPUTS
OUTPUT
SOURCE
3107 subroutine lwl_write(path, cryst, vcp, npwe, nomega, gvec, chi0, chi0_head, chi0_lwing, chi0_uwing, comm) 3108 3109 !Arguments ------------------------------------ 3110 !scalars 3111 integer,intent(in) :: npwe,nomega,comm 3112 character(len=*),intent(in) :: path 3113 type(crystal_t),intent(in) :: cryst 3114 type(vcoul_t),intent(in) :: Vcp 3115 !arrays 3116 integer,intent(in) :: gvec(3,npwe) 3117 complex(gwpc),intent(in) :: chi0(npwe,npwe,nomega) 3118 complex(dpc),intent(inout) :: chi0_head(3,3,nomega),chi0_lwing(npwe,nomega,3),chi0_uwing(npwe,nomega,3) 3119 3120 !Local variables------------------------------- 3121 !scalars 3122 integer,parameter :: master=0,prtvol=0 3123 integer :: iw,ii,iomode,unt,my_rank 3124 character(len=500) :: msg 3125 real(dp) :: length 3126 complex(dpc) :: wng(3),em1_00 3127 ! type(hscr_t),intent(out) :: hscr 3128 !arrays 3129 complex(dpc),allocatable :: wtest(:),eps_head(:,:,:) 3130 3131 ! ************************************************************************* 3132 3133 !if (xmpi_comm_rank(comm) /= master) goto 100 3134 my_rank = xmpi_comm_rank(comm) 3135 3136 ABI_MALLOC(wtest,(npwe)) 3137 3138 if (prtvol > 0 .and. my_rank == master) then 3139 call wrtout(std_out, "head of chi0") 3140 do iw=1,nomega 3141 call print_arr(chi0_head(:,:,iw),max_r=3,max_c=3) 3142 end do 3143 3144 do iw=1,nomega 3145 call wrtout(std_out, "symmetrized e_00 via tensor") 3146 wng = MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT) 3147 write(std_out,*) one - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1)*Vcp%vcqlwl_sqrt(1,1) 3148 3149 call wrtout(std_out, "symmetrized e_0G via tensor") 3150 do ii=1,npwe 3151 wng = chi0_uwing(ii,iw,:) 3152 wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1) 3153 end do 3154 call print_arr(wtest,max_r=9,unit=std_out) 3155 3156 call wrtout(std_out, "symmetrized e_G0 via tensor") 3157 do ii=1,npwe 3158 wng = chi0_lwing(ii,iw,:) 3159 wtest(ii) = - vdotw(GW_Q0_DEFAULT,wng,cryst%gmet,"G") * Vcp%vcqlwl_sqrt(1,1) * Vcp%vcqlwl_sqrt(ii,1) 3160 end do 3161 call print_arr(wtest,max_r=9,unit=std_out) 3162 end do 3163 end if 3164 3165 ! Write chi0 data 3166 iomode = IO_MODE_FORTRAN; if (endswith(path, ".nc")) iomode = IO_MODE_ETSF 3167 3168 if (my_rank == master) then 3169 if (iomode == IO_MODE_FORTRAN) then 3170 if (open_file(path,msg,newunit=unt,form="unformatted", action="write") /= 0) then 3171 ABI_ERROR(msg) 3172 end if 3173 !call hscr_io(er%hscr,fform,rdwr,unt,comm,master,iomode) 3174 do iw=1,nomega 3175 write(unt)chi0_head(:,:,iw) 3176 end do 3177 !do iw=1,nomega 3178 ! write(unt)chi0_lwing(:,iw,:) 3179 !end do 3180 !do iw=1,nomega 3181 ! write(unt)chi0_uwing(:,iw,:) 3182 !end do 3183 3184 else 3185 ABI_ERROR(sjoin("iomode", itoa(iomode), "is not supported")) 3186 end if 3187 end if 3188 3189 ABI_MALLOC(eps_head,(3,3,nomega)) 3190 call mkem1_q0(npwe,1,1,nomega,cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm) 3191 3192 if (my_rank == master) then 3193 if (iomode == IO_MODE_FORTRAN) then 3194 do iw=1,nomega 3195 write(unt)eps_head(:,:,iw) 3196 end do 3197 !do iw=1,nomega 3198 ! write(unt)chi0_lwing(:,iw,:) 3199 !end do 3200 !do iw=1,nomega 3201 ! write(unt)chi0_uwing(:,iw,:) 3202 !end do 3203 else 3204 ABI_ERROR(sjoin("iomode:", itoa(iomode), "is not supported")) 3205 end if 3206 end if 3207 3208 ABI_FREE(eps_head) 3209 3210 if (prtvol > 0 .and. my_rank == master) then 3211 length = normv(GW_Q0_DEFAULT,cryst%gmet,"G") 3212 3213 do iw=1,nomega 3214 em1_00 = one / vdotw(GW_Q0_DEFAULT/length, MATMUL(chi0_head(:,:,iw),GW_Q0_DEFAULT/length),cryst%gmet,"G") 3215 call wrtout(std_out, "e^1_{00} from tensor") 3216 write(std_out,*) em1_00 3217 3218 call wrtout(std_out, "symmetrized e^-1_0G via tensor") 3219 do ii=1,npwe 3220 wng = chi0_uwing(ii,iw,:) 3221 wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G") 3222 end do 3223 wtest(1) = em1_00 3224 call print_arr(wtest,max_r=9,unit=std_out) 3225 3226 call wrtout(std_out, "symmetrized e^-1_G0 via tensor") 3227 do ii=1,npwe 3228 wng = chi0_lwing(ii,iw,:) 3229 wtest(ii) = em1_00*vdotw(GW_Q0_DEFAULT/length,wng,cryst%gmet,"G") 3230 end do 3231 wtest(1) = em1_00 3232 call print_arr(wtest,max_r=9,unit=std_out) 3233 end do !iw 3234 end if 3235 3236 ABI_FREE(wtest) 3237 3238 if (my_rank == master) then 3239 if (iomode == IO_MODE_FORTRAN) then 3240 close(unt) 3241 else 3242 NCF_CHECK(nf90_close(unt)) 3243 end if 3244 end if 3245 3246 !100 call xmpi_barrier(comm) 3247 3248 end subroutine lwl_write
m_screening/make_epsm1_driver [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
make_epsm1_driver
FUNCTION
Driver routine to calculate the inverse symmetrical dielectric matrix starting from the irreducible polarizability. The routine considers a single q-point, and performs the following tasks: 1) Calculate $\tilde\epsilon^{-1}$ using different approximations: * RPA * ALDA within TDDFT 2) Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space calculating these quantities for different small q-directions specified by the user (Not yet operative) 3) Output the electron energy loss function and the macroscopic dielectric function with and without local field effects (only if non-zero real frequencies are available)
INPUTS
iqibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated dim_wing=Dimension of the wings (0 or 3 if q-->0) npwe=Number of G-vectors in chi0. nI,nJ=Number of rows/columns in chi0_ij (1,1 if collinear case) nomega=Number of frequencies. omega(nomega)=Frequencines in Hartree approx_type=Integer flag defining the type of approximation == 0 for RPA == == 1 for TDDFT == option_test=Only for TDDFT: == 0 for TESTPARTICLE == == 1 for TESTELECTRON == Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction nfftot=Total number of points in the FFT mesh. ngfft(18)=Info on the FFT mesh. nkxc=Dimension of the kernel in reciprocal space. 0 if kernel is not needed kxcg(nfftot,nkxc)=TDDFT kernel in reciprocal space on the FFT mesh. Used only if approx_type==1 gvec(3,npwe)=G-vectors comm=MPI communicator. chi0_lwing(npwe*nI,nomega,dim_wing)=Lower wings of chi0 (only for q-->0) chi0_uwing(npwe*nJ,nomega,dim_wing)=Upper wings of chi0 (only for q-->0) chi0_head(dim_wing,dim_wing,nomega)=Head of of chi0 (only for q-->0)
OUTPUT
spectra<spectra_t>Object containing e_macro(w) and EELS(w)
SIDE EFFECTS
chi0(npwe*nI,npwe*nJ,nomega): in input the irreducible polarizability, in output the symmetrized inverse dielectric matrix.
SOURCE
1339 subroutine make_epsm1_driver(iqibz,dim_wing,npwe,nI,nJ,nomega,omega,& 1340 approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,chi0_head,& 1341 chi0_lwing,chi0_uwing,chi0,spectra,comm,& 1342 fxc_ADA,rhor) ! optional argument 1343 1344 !Arguments ------------------------------------ 1345 !scalars 1346 integer,intent(in) :: iqibz,nI,nJ,npwe,nomega,dim_wing,approx_type,option_test,nkxc,nfftot,comm 1347 real(dp),intent(in),optional :: rhor 1348 type(vcoul_t),target,intent(in) :: Vcp 1349 type(spectra_t),intent(out) :: Spectra 1350 !arrays 1351 integer,intent(in) :: ngfft(18),gvec(3,npwe) 1352 complex(gwpc),intent(in) :: kxcg(nfftot,nkxc) 1353 complex(dpc),intent(in) :: omega(nomega) 1354 complex(dpc),intent(inout) :: chi0_lwing(:,:,:) !(npwe*nI,nomega,dim_wing) 1355 complex(dpc),intent(inout) :: chi0_uwing(:,:,:) !(npwe*nJ,nomega,dim_wing) 1356 complex(dpc),intent(inout) :: chi0_head(:,:,:) !(dim_wing,dim_wing,nomega) 1357 complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ,nomega) 1358 complex(gwpc),intent(in),optional :: fxc_ADA(npwe*nI,npwe*nJ) 1359 1360 !Local variables------------------------------- 1361 !scalars 1362 integer,parameter :: master=0 1363 integer :: i1,i2,ig1,ig2,io,ierr,irank,my_nqlwl !iqlwl 1364 integer :: nor,my_rank,nprocs,comm_self,g1mg2_idx 1365 real(dp) :: ucvol 1366 logical :: is_qeq0,use_MPI 1367 character(len=500) :: msg 1368 !arrays 1369 integer :: omega_distrb(nomega) 1370 integer,allocatable :: istart(:),istop(:) 1371 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 1372 real(dp),allocatable :: eelf(:,:),tmp_eelf(:) 1373 complex(dpc),allocatable :: epsm_lf(:,:),epsm_nlf(:,:),tmp_lf(:),tmp_nlf(:) 1374 complex(dpc),allocatable :: buffer_lwing(:,:),buffer_uwing(:,:) 1375 complex(gwpc),allocatable :: kxcg_mat(:,:) 1376 1377 !bootstrap and LR 1378 integer :: istep,nstep 1379 real(dp) :: conv_err, alpha, Zr, qpg2(3), qpg2_nrm 1380 real(gwpc) :: chi00_head, fxc_head 1381 complex(gwpc),allocatable :: vfxc_boot(:,:), vfxc_boot0(:,:), vfxc_lr(:,:), vfxc_tmp(:,:), chi0_tmp(:,:), chi0_save(:,:,:) 1382 complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:) 1383 1384 ! ************************************************************************* 1385 1386 DBG_ENTER("COLL") 1387 1388 if (nI/=1.or.nJ/=1) then 1389 ABI_ERROR("nI or nJ=/1 not yet implemented") 1390 end if 1391 1392 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1393 1394 ! MG TODO We use comm_self for the inversion as the single precision version is not yet available 1395 comm_self = xmpi_comm_self 1396 1397 call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol) 1398 1399 is_qeq0 = normv(Vcp%qibz(:,iqibz),gmet,'G') < GW_TOLQ0 1400 1401 omega_distrb = my_rank 1402 use_MPI = .FALSE. 1403 use_MPI = nprocs >= nomega ! Parallelism is not used 1404 1405 if (use_MPI) then 1406 ! Initialize distribution table for frequencies. 1407 ABI_MALLOC(istart,(nprocs)) 1408 ABI_MALLOC(istop,(nprocs)) 1409 call xmpi_split_work2_i4b(nomega,nprocs,istart,istop) 1410 omega_distrb(:)=xmpi_undefined_rank 1411 do irank=0,nprocs-1 1412 i1 = istart(irank+1) 1413 i2 = istop (irank+1) 1414 if (i1<=i2) omega_distrb(i1:i2) = irank 1415 end do 1416 ABI_FREE(istart) 1417 ABI_FREE(istop) 1418 end if 1419 1420 ! Initialize container for spectral results 1421 do nor=1,nomega 1422 if (ABS(AIMAG(omega(nor)))>1.e-3) EXIT 1423 end do 1424 nor=nor-1; if (nor==0) nor = 1 ! only imag !? 1425 1426 if (dim_wing==3) then 1427 call wrtout(std_out,' Analyzing long wavelength limit for several q') 1428 call spectra_init(Spectra,nor,REAL(omega(1:nor)),Vcp%nqlwl,Vcp%qlwl) 1429 my_nqlwl = 1 1430 !my_nqlwl = dim_wing ! TODO 1431 !ABI_CHECK(dim_wing==SIZE(Vcp%vcqlwl_sqrt,DIM=2),"WRONG DIMS") 1432 else 1433 call spectra_init(Spectra,nor,REAL(omega(1:nor)),1,Vcp%qibz(:,iqibz)) 1434 my_nqlwl = 1 1435 end if 1436 ! 1437 ! NOTE: all processors have to perform this operation in order to have the 1438 ! epsm1 matrix when performing a sigma calculation starting with the file _SUS 1439 ! 1440 ! Temporary arrays to store spectra. 1441 ABI_MALLOC(epsm_lf,(nomega,my_nqlwl)) 1442 ABI_MALLOC(epsm_nlf,(nomega,my_nqlwl)) 1443 ABI_MALLOC(eelf,(nomega,my_nqlwl)) 1444 epsm_lf=czero; epsm_nlf=czero; eelf=zero 1445 1446 ! Temporary arrays used to store output results. 1447 ABI_MALLOC(tmp_lf, (my_nqlwl)) 1448 ABI_MALLOC(tmp_nlf, (my_nqlwl)) 1449 ABI_MALLOC(tmp_eelf, (my_nqlwl)) 1450 1451 SELECT CASE (approx_type) 1452 1453 CASE (0) 1454 ! RPA: \tepsilon = 1 - Vc^{1/2} chi0 Vc^{1/2} 1455 ! vc_sqrt contains vc^{1/2}(q,G), complex-valued to allow for a possible cutoff. 1456 do io=1,nomega 1457 if (omega_distrb(io) == my_rank) then 1458 !write(std_out,*)"dim_wing",dim_wing 1459 call rpa_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),my_nqlwl,dim_wing, & 1460 chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:), & 1461 tmp_lf,tmp_nlf,tmp_eelf,comm_self) 1462 1463 ! Store results. 1464 epsm_lf(io,:) = tmp_lf 1465 epsm_nlf(io,:) = tmp_nlf 1466 eelf(io,:) = tmp_eelf 1467 end if 1468 end do ! nomega 1469 1470 CASE (1) 1471 ! Vertex correction from Adiabatic TDDFT. chi_{G1,G2} = [\delta -\chi0 (vc+kxc)]^{-1}_{G1,G3} \chi0_{G3,G2} 1472 ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded") 1473 ABI_CHECK(nkxc==1,"nkxc/=1 not coded") 1474 1475 ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh. 1476 ABI_MALLOC_OR_DIE(kxcg_mat,(npwe,npwe), ierr) 1477 1478 ierr=0 1479 do ig2=1,npwe 1480 do ig1=1,npwe 1481 g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft) 1482 if (g1mg2_idx>0) then 1483 kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1) 1484 else 1485 ierr=ierr+1 1486 kxcg_mat(ig1,ig2) = czero 1487 end if 1488 end do 1489 end do 1490 1491 if (ierr/=0) then 1492 write(msg,'(a,i4,3a)')& 1493 & 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,& 1494 & 'Enlarge the FFT mesh to get rid of this problem. ' 1495 ABI_WARNING(msg) 1496 end if 1497 1498 !FIXME "recheck TDDFT code and parallel" 1499 ABI_CHECK(nkxc==1,"nkxc/=1 not coded") 1500 do io=1,nomega 1501 if (omega_distrb(io) == my_rank) then 1502 call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),& 1503 & chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm) 1504 1505 ! Store results. 1506 epsm_lf(io,:) = tmp_lf 1507 epsm_nlf(io,:) = tmp_nlf 1508 eelf(io,:) = tmp_eelf 1509 end if 1510 end do 1511 1512 ABI_FREE(kxcg_mat) 1513 do io=1,nomega 1514 write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]' 1515 call wrtout(std_out,msg) 1516 call print_arr(chi0(:,:,io),mode_paral='PERS') 1517 end do 1518 1519 CASE (2) 1520 ! ADA nonlocal vertex correction contained in fxc_ADA 1521 ABI_WARNING('Entered fxc_ADA branch: EXPERIMENTAL!') 1522 ! Test that argument was passed 1523 if (.not.present(fxc_ADA)) then 1524 ABI_ERROR('make_epsm1_driver was not called with optional argument fxc_ADA') 1525 end if 1526 ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded") 1527 1528 do io=1,nomega 1529 if (omega_distrb(io) == my_rank) then 1530 call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),fxc_ADA,option_test,my_nqlwl,dim_wing,omega(io),& 1531 & chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm) 1532 1533 ! Store results. 1534 epsm_lf(io,:) = tmp_lf 1535 epsm_nlf(io,:) = tmp_nlf 1536 eelf(io,:) = tmp_eelf 1537 end if 1538 end do 1539 1540 do io=1,nomega 1541 write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]' 1542 call wrtout(std_out,msg) 1543 call print_arr(chi0(:,:,io),mode_paral='PERS') 1544 end do 1545 1546 CASE (4) 1547 ! Bootstrap vertex kernel by Sharma [[cite:Sharma2011]] 1548 ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr) 1549 ABI_MALLOC_OR_DIE(vfxc_boot0,(npwe*nI,npwe*nJ), ierr) 1550 ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr) 1551 ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr) 1552 1553 if (iqibz==1) then 1554 vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0 1555 else 1556 vc_sqrt => Vcp%vc_sqrt(:,iqibz) 1557 end if 1558 1559 chi0_save = chi0 ! a copy of chi0 (ks) 1560 nstep = 50 ! max iteration steps 1561 alpha = 0.6 ! mixing 1562 chi00_head = chi0(1,1,1)*vc_sqrt(1)**2 1563 fxc_head = czero; vfxc_boot = czero; chi0_tmp = czero 1564 epsm_lf = czero; epsm_nlf = czero; eelf = zero 1565 write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head 1566 call wrtout(std_out,msg) 1567 ! loop 1568 conv_err = 0.1 1569 do istep=1, nstep 1570 chi0 = chi0_save 1571 io=1 ! for now only at omega=0 1572 call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),& 1573 & chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self) 1574 conv_err = smallest_real 1575 do ig2=1,npwe*nJ 1576 do ig1=1,npwe*nI 1577 conv_err= MAX(conv_err, ABS(chi0(ig1,ig2,1) - chi0_tmp(ig1,ig2))) 1578 end do 1579 end do 1580 write(msg,'(a,i4,a,f10.6)') ' => bootstrap itr# ', istep, ', eps^-1 max error: ', conv_err 1581 call wrtout(std_out,msg) 1582 write(msg,'(a,2f10.6)') ' eps^-1(head): ', chi0(1,1,1) 1583 call wrtout(std_out,msg) 1584 write(msg,'(a,2f10.6)') ' v^-1*fxc(head): ', fxc_head 1585 call wrtout(std_out,msg) 1586 if (conv_err <= tol4) exit 1587 ! 1588 chi0_tmp = chi0(:,:,1) 1589 vfxc_boot = chi0(:,:,1)/chi00_head 1590 if (istep > 1) then 1591 vfxc_boot = alpha*vfxc_boot0 + (one-alpha)*vfxc_boot 1592 end if 1593 vfxc_boot0 = vfxc_boot 1594 fxc_head = vfxc_boot(1,1) 1595 do ig1=1,npwe 1596 vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:) 1597 end do 1598 end do 1599 ! end loop 1600 if (istep <= nstep) then 1601 write(msg,'(a,i4,a)') ' => bootstrap fxc converged after ', istep, ' iterations' 1602 call wrtout(std_out,msg) 1603 else 1604 write(msg,'(a,i4,a)') ' -> bootstrap fxc not converged after ', nstep, ' iterations' 1605 ABI_WARNING(msg) 1606 end if 1607 ! 1608 chi0 = chi0_save 1609 do io=1,nomega 1610 if (omega_distrb(io) == my_rank) then 1611 call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),& 1612 & chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self) 1613 epsm_lf(io,:) = tmp_lf 1614 epsm_nlf(io,:) = tmp_nlf 1615 eelf(io,:) = tmp_eelf 1616 end if 1617 end do 1618 1619 ABI_FREE(chi0_tmp) 1620 ABI_FREE(chi0_save) 1621 ABI_FREE(vfxc_boot) 1622 ABI_FREE(vfxc_boot0) 1623 1624 do io=1,nomega 1625 write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]' 1626 call wrtout(std_out,msg) 1627 call print_arr(chi0(:,:,io),mode_paral='PERS') 1628 end do 1629 1630 CASE(5) 1631 ! One-shot scalar bootstrap approximation 1632 ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr) 1633 ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr) 1634 1635 if (iqibz==1) then 1636 vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0 1637 else 1638 vc_sqrt => Vcp%vc_sqrt(:,iqibz) 1639 end if 1640 1641 chi0_save = chi0 ! a copy of chi0 1642 fxc_head = czero; vfxc_boot = czero; 1643 epsm_lf = czero; epsm_nlf = czero; eelf = zero 1644 chi00_head = chi0(1,1,1)*vc_sqrt(1)**2 1645 write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head 1646 call wrtout(std_out,msg) 1647 1648 fxc_head = vc_sqrt(1)**2/chi00_head + vc_sqrt(1)**2/chi00_head - vc_sqrt(1)**2 1649 fxc_head = 0.5*fxc_head + 0.5*sqrt(fxc_head**2 - 4.0*vc_sqrt(1)**4/(chi00_head*chi00_head)) 1650 vfxc_boot(1,1) = fxc_head 1651 write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head/vc_sqrt(1)**2 1652 call wrtout(std_out,msg) 1653 1654 chi0 = chi0_save 1655 do io=1,nomega 1656 if (omega_distrb(io) == my_rank) then 1657 call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),& 1658 & chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self) 1659 epsm_lf(io,:) = tmp_lf 1660 epsm_nlf(io,:) = tmp_nlf 1661 eelf(io,:) = tmp_eelf 1662 end if 1663 end do 1664 write(msg,'(a,2f10.6)') ' eps^-1(head): ',chi0(1,1,1) 1665 call wrtout(std_out,msg) 1666 1667 ABI_FREE(chi0_save) 1668 ABI_FREE(vfxc_boot) 1669 1670 do io=1,nomega 1671 write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]' 1672 call wrtout(std_out,msg) 1673 call print_arr(chi0(:,:,io),mode_paral='PERS') 1674 end do 1675 1676 CASE(6) 1677 ! RPA bootstrap by Rigamonti [[cite:Rigamonti2015]] and Berger [[cite:Berger2015]] 1678 ABI_MALLOC_OR_DIE(vfxc_boot,(npwe*nI,npwe*nJ), ierr) 1679 ABI_MALLOC_OR_DIE(chi0_save,(npwe*nI,npwe*nJ,nomega), ierr) 1680 ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr) 1681 1682 if (iqibz==1) then 1683 vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0 1684 else 1685 vc_sqrt => Vcp%vc_sqrt(:,iqibz) 1686 end if 1687 1688 chi0_save = chi0 ! a copy of chi0 1689 fxc_head = czero; vfxc_boot = czero; 1690 epsm_lf = czero; epsm_nlf = czero; eelf = zero 1691 !chi00_head = chi0(1,1,1)*vc_sqrt(1)**2 1692 write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ',chi00_head 1693 call wrtout(std_out,msg) 1694 1695 io = 1 ! static 1696 call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,0,my_nqlwl,dim_wing,omega(io),& 1697 & chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self) 1698 epsm_lf(1,:) = tmp_lf 1699 1700 ! chi(RPA) = chi0 * (1 - chi0 * v_c)^-1 1701 chi0 = chi0_save 1702 do ig2=2,npwe 1703 do ig1=2,npwe 1704 chi0(ig1,ig2,io)=-vc_sqrt(ig1)*chi0(ig1,ig2,io)*vc_sqrt(ig2) 1705 end do 1706 chi0(ig2,ig2,io)=one+chi0(ig2,ig2,io) 1707 end do 1708 chi0(1,:,io) = czero; chi0(:,1,io) = czero; chi0(1,1,io) = one 1709 chi0_tmp = chi0(:,:,io) 1710 call xginv(chi0_tmp,npwe,comm=comm) 1711 chi0 = chi0_save 1712 chi0_tmp = MATMUL(chi0(:,:,io), chi0_tmp(:,:)) ! chi(RPA) 1713 do ig1=1,npwe 1714 chi0_tmp(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*chi0_tmp(ig1,:) 1715 end do 1716 !call xginv(chi0_tmp,npwe,comm=comm) ! chi(RPA)^-1 1717 !vfxc_boot = chi0_tmp/epsm_lf(1,1) 1718 ! 1719 !vfxc_boot(1,1) = chi0_tmp(1,1)/epsm_lf(1,1) 1720 vfxc_boot(1,1) = one/chi0_tmp(1,1)/epsm_lf(1,1) 1721 !@WC: alternatively: 1722 !chi00_head = chi0(1,1,io)*vc_sqrt(1)**2 1723 !vfxc_boot(1,1) = one/chi00_head/epsm_lf(1,1) 1724 fxc_head = vfxc_boot(1,1) 1725 do ig1=1,npwe 1726 vfxc_boot(ig1,:) = vc_sqrt(ig1)*vc_sqrt(:)*vfxc_boot(ig1,:) 1727 end do 1728 write(msg,'(a,2f10.6)') ' -> v^-1*fxc(head): ',fxc_head 1729 call wrtout(std_out,msg) 1730 1731 chi0 = chi0_save 1732 do io=1,nomega 1733 if (omega_distrb(io) == my_rank) then 1734 call atddft_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_boot,option_test,my_nqlwl,dim_wing,omega(io),& 1735 & chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self) 1736 epsm_lf(io,:) = tmp_lf 1737 epsm_nlf(io,:) = tmp_nlf 1738 eelf(io,:) = tmp_eelf 1739 end if 1740 end do 1741 write(msg,'(a,2f10.6)') ' eps^-1(head): ',chi0(1,1,1) 1742 call wrtout(std_out,msg) 1743 1744 ABI_FREE(chi0_save) 1745 ABI_FREE(vfxc_boot) 1746 ABI_FREE(chi0_tmp) 1747 1748 do io=1,nomega 1749 write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]' 1750 call wrtout(std_out,msg) 1751 call print_arr(chi0(:,:,io),mode_paral='PERS') 1752 end do 1753 1754 CASE (7) 1755 ! LR+ALDA hybrid vertex kernel by Tal 1756 ! First ALDA 1757 ABI_CHECK(Vcp%nqlwl==1,"nqlwl/=1 not coded") 1758 ABI_CHECK(nkxc==1,"nkxc/=1 not coded") 1759 ! Make kxcg_mat(G1,G2) = kxcg(G1-G2) from kxcg defined on the FFT mesh. 1760 ABI_MALLOC_OR_DIE(kxcg_mat,(npwe,npwe), ierr) 1761 ierr=0 1762 do ig2=1,npwe 1763 do ig1=1,npwe 1764 g1mg2_idx = g2ifft(gvec(:,ig1)-gvec(:,ig2),ngfft) 1765 if (g1mg2_idx>0) then 1766 kxcg_mat(ig1,ig2) = kxcg(g1mg2_idx,1) 1767 else 1768 ierr=ierr+1 1769 kxcg_mat(ig1,ig2) = czero 1770 end if 1771 end do 1772 end do 1773 if (ierr/=0) then 1774 write(msg,'(a,i4,3a)')& 1775 & 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,& 1776 & 'Enlarge the FFT mesh to get rid of this problem. ' 1777 ABI_WARNING(msg) 1778 end if 1779 !FIXME "recheck TDDFT code and parallel" 1780 ABI_CHECK(nkxc==1,"nkxc/=1 not coded") 1781 1782 ! Now LR: (1-Z)*chi0^-1 1783 ABI_MALLOC_OR_DIE(vfxc_lr,(npwe*nI,npwe*nJ), ierr) 1784 ABI_MALLOC_OR_DIE(vfxc_tmp,(npwe*nI,npwe*nJ), ierr) 1785 ABI_MALLOC_OR_DIE(chi0_tmp,(npwe*nI,npwe*nJ), ierr) 1786 1787 if (iqibz==1) then 1788 vc_sqrt => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0 1789 else 1790 vc_sqrt => Vcp%vc_sqrt(:,iqibz) 1791 end if 1792 1793 Zr = 0.78 1794 chi00_head = chi0(1,1,1)*vc_sqrt(1)**2 1795 fxc_head = czero; vfxc_lr = czero; vfxc_tmp = czero 1796 epsm_lf = czero; epsm_nlf = czero; eelf = zero 1797 write(msg,'(a,2f10.6)') ' -> chi0_dft(head): ', chi00_head 1798 call wrtout(std_out,msg) 1799 ! 1800 chi0_tmp = chi0(:,:,1) 1801 call xginv(chi0_tmp,npwe,comm=comm) 1802 vfxc_lr = (one-Zr)*chi0_tmp(:,:) 1803 write(msg,'(a)') ' Constructing LR+ALDA fxc kernel' 1804 call wrtout(std_out,msg) 1805 ! 1806 do ig1=1,npwe 1807 do ig2=1,npwe 1808 qpg2 = Vcp%qibz(:,iqibz) + gvec(:,ig1) 1809 qpg2_nrm = normv(qpg2,gmet,"G") 1810 qpg2 = Vcp%qibz(:,iqibz) + gvec(:,ig2) 1811 qpg2_nrm = SQRT(qpg2_nrm * normv(qpg2,gmet,"G")) 1812 vfxc_tmp(ig1,ig2) = vfxc_lr(ig1,ig2)*exp(-(qpg2_nrm/k_thfermi(rhor))**2) + & 1813 & kxcg_mat(ig1,ig2)*(one - exp(-(qpg2_nrm/k_thfermi(rhor))**2)) 1814 !write(std_out,*) ig1, qpg2_nrm, k_thfermi(rhor), vfxc_lr(ig1,ig1), kxcg_mat(ig1,ig1), vfxc_tmp(ig1,ig1) 1815 end do 1816 end do 1817 ! 1818 vfxc_lr = vfxc_tmp 1819 1820 do io=1,nomega 1821 if (omega_distrb(io) == my_rank) then 1822 call atddft_hyb_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0(:,:,io),vfxc_lr,kxcg_mat,option_test,my_nqlwl,dim_wing,omega(io),& 1823 & chi0_head(:,:,io),chi0_lwing(:,io,:),chi0_uwing(:,io,:),tmp_lf,tmp_nlf,tmp_eelf,comm_self) 1824 epsm_lf(io,:) = tmp_lf 1825 epsm_nlf(io,:) = tmp_nlf 1826 eelf(io,:) = tmp_eelf 1827 end if 1828 end do 1829 1830 ABI_FREE(kxcg_mat) 1831 ABI_FREE(chi0_tmp) 1832 ABI_FREE(vfxc_lr) 1833 ABI_FREE(vfxc_tmp) 1834 1835 do io=1,nomega 1836 write(msg,'(a,i4,a,2f9.4,a)')' Symmetrical epsilon^-1(G,G'') at the ',io,' th omega',omega(io)*Ha_eV,' [eV]' 1837 call wrtout(std_out,msg) 1838 call print_arr(chi0(:,:,io),mode_paral='PERS') 1839 end do 1840 1841 CASE DEFAULT 1842 ABI_BUG(sjoin('Wrong approx_type:',itoa(approx_type))) 1843 END SELECT 1844 1845 if (use_MPI) then 1846 ! Collect results on each node. 1847 ABI_MALLOC(buffer_lwing, (size(chi0_lwing,dim=1), size(chi0_lwing, dim=3))) 1848 ABI_MALLOC(buffer_uwing, (size(chi0_uwing,dim=1), size(chi0_uwing, dim=3))) 1849 1850 do io=1,nomega 1851 if (omega_distrb(io)/=my_rank) then 1852 ! Zero arrays. 1853 chi0(:,:,io) = czero_gw 1854 if(dim_wing>0) then 1855 chi0_lwing(:,io,:) = zero 1856 chi0_uwing(:,io,:) = zero 1857 chi0_head(:,:,io) = czero 1858 endif 1859 epsm_lf(io,:) = czero 1860 epsm_nlf(io,:) = czero 1861 eelf(io,:) = zero 1862 end if 1863 1864 call xmpi_sum(chi0(:,:,io), comm,ierr) 1865 1866 if (dim_wing>0) then 1867 ! Build contiguous arrays 1868 buffer_lwing = chi0_lwing(:,io,:) 1869 buffer_uwing = chi0_uwing(:,io,:) 1870 call xmpi_sum(buffer_lwing,comm,ierr) 1871 call xmpi_sum(buffer_uwing,comm,ierr) 1872 chi0_lwing(:,io,:) = buffer_lwing 1873 chi0_uwing(:,io,:) = buffer_uwing 1874 if (size(chi0_head(:,:,io))/= zero) then 1875 call xmpi_sum(chi0_head(:,:,io),comm,ierr) 1876 end if 1877 end if 1878 1879 end do ! iw 1880 1881 call xmpi_sum(epsm_lf, comm,ierr ) 1882 call xmpi_sum(epsm_nlf,comm,ierr) 1883 call xmpi_sum(eelf, comm,ierr) 1884 ABI_FREE(buffer_lwing) 1885 ABI_FREE(buffer_uwing) 1886 end if 1887 1888 ! Save results in Spectra%, mind the slicing. 1889 Spectra%emacro_nlf(:,:) = epsm_nlf(1:nor,:) 1890 Spectra%emacro_lf (:,:) = epsm_lf (1:nor,:) 1891 Spectra%eelf (:,:) = eelf (1:nor,:) 1892 1893 ABI_FREE(epsm_lf) 1894 ABI_FREE(epsm_nlf) 1895 ABI_FREE(eelf) 1896 ABI_FREE(tmp_lf) 1897 ABI_FREE(tmp_nlf) 1898 ABI_FREE(tmp_eelf) 1899 1900 DBG_EXIT("COLL") 1901 1902 end subroutine make_epsm1_driver
m_screening/mdielf_bechstedt [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
mdielf_bechstedt
FUNCTION
Calculates the model dielectric function for the homogeneous system as proposed by F. Bechstedt, in Solid State Commun. 84, 765 1992.
INPUTS
eps_inf=Dielectric constant of the material qnrm=The modulus of the q-point. rhor=The local value of the density
SOURCE
2816 elemental function mdielf_bechstedt(eps_inf,qnrm,rhor) result(mdielf) 2817 2818 !Arguments ------------------------------------ 2819 !scalars 2820 real(dp),intent(in) :: eps_inf,qnrm,rhor 2821 real(dp) :: mdielf 2822 2823 ! ************************************************************************* 2824 2825 mdielf = one + & 2826 & one / ( one/(eps_inf-one) + (qnrm/k_thfermi(rhor))**2 + (three*qnrm**4)/(four*k_fermi(rhor)**2 * k_thfermi(rhor)**2) ) 2827 2828 end function mdielf_bechstedt
m_screening/mkdump_Er [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
mkdump_Er
FUNCTION
Dump the content of an Epsilonm1_results data type on file.
INPUTS
id_required=Identifier of the matrix to be calculated Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction ngfft(18)=Info on the FFT mesh. nfftot=Total number of point on the FFT mesh. gvec(3,npwe)=Reduced coordinates of plane waves for the response functions npwe=Number of plane waves. comm=MPI communicator.
OUTPUT
SOURCE
839 subroutine mkdump_Er(Er,Vcp,npwe,gvec,nkxc,kxcg,id_required,approx_type,& 840 & ikxc_required,option_test,fname_dump,iomode,& 841 & nfftot,ngfft,comm,fxc_ADA) 842 843 !Arguments ------------------------------------ 844 !scalars 845 integer,intent(in) :: id_required,approx_type,option_test,ikxc_required,nkxc 846 integer,intent(in) :: iomode,nfftot,npwe,comm 847 type(Epsilonm1_results),intent(inout) :: Er 848 type(vcoul_t),intent(in) :: Vcp 849 character(len=*),intent(in) :: fname_dump 850 !arrays 851 integer,intent(in) :: ngfft(18),gvec(3,npwe) 852 complex(gwpc),intent(in) :: kxcg(nfftot,nkxc) 853 complex(gwpc),intent(in), optional :: fxc_ADA(npwe*Er%nI,npwe*Er%nJ,Er%nqibz) 854 855 !Local variables------------------------------- 856 !scalars 857 integer,parameter :: master=0 858 integer :: dim_wing,iqibz,is_qeq0,mqmem_,npwe_asked 859 integer :: unt_dump,fform,rdwr,ierr 860 integer :: my_rank,comm_self 861 real(dp) :: ucvol 862 character(len=500) :: msg 863 character(len=fnlen) :: ofname 864 character(len=nctk_slen) :: in_varname,out_varname 865 type(hscr_t) :: Hscr_cp 866 type(spectra_t) :: spectra 867 !arrays 868 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 869 complex(gwpc),allocatable :: epsm1(:,:,:) 870 complex(dpc),allocatable :: dummy_lwing(:,:,:),dummy_uwing(:,:,:),dummy_head(:,:,:) 871 872 ! ********************************************************************* 873 874 DBG_ENTER("COLL") 875 876 ABI_CHECK(id_required==4,'Value of id_required not coded') 877 ABI_CHECK(npwe==Er%npwe,"mismatch in npwe") 878 879 my_rank = xmpi_comm_rank(comm); comm_self = xmpi_comm_self 880 call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol) 881 882 ! if (Er%ID/=0) call reset_Epsilonm1(Er) 883 Er%ID=id_required 884 885 ofname = fname_dump 886 in_varname = ncname_from_id(er%hscr%id) 887 out_varname = ncname_from_id(id_required) 888 889 write(std_out,*)'Er%ID: ',Er%ID,', Er%Hscr%ID: ',Er%Hscr%ID 890 891 if (Er%ID==Er%Hscr%ID) then 892 ! === The two-point function we are asking for is already stored on file === 893 ! * According to mqmem either read and store the entire matrix in memory or do nothing. 894 895 if (Er%mqmem>0) then 896 ! In-core solution. 897 write(msg,'(a,f12.1,a)')' Memory needed for Er%epsm1 = ',two*gwpc*npwe**2*Er%nomega*Er%nqibz*b2Mb,' [Mb] <<< MEM' 898 call wrtout(std_out,msg) 899 ABI_MALLOC_OR_DIE(Er%epsm1,(npwe,npwe,Er%nomega,Er%nqibz), ierr) 900 901 if (iomode == IO_MODE_MPI) then 902 !call wrtout(std_out, "read_screening with MPI_IO") 903 ABI_WARNING("SUSC files is buggy. Using Fortran IO") 904 call read_screening(in_varname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm) 905 else 906 call read_screening(in_varname,Er%fname,Er%npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm) 907 end if 908 else 909 ! Out-of-core solution === 910 ABI_COMMENT("mqmem==0 => allocating a single q-slice of (W|chi0) (slower but less memory).") 911 continue 912 end if 913 914 return 915 916 else 917 ! === The matrix stored on file do not correspond to the quantity required === 918 ! * Presently only the transformation chi0 => e^-1 is coded 919 ! * According to Er%mqmem either calculate e^-1 dumping the result to a file 920 ! for a subsequent use or calculate e^-1 keeping everything in memory. 921 922 if (Er%mqmem==0) then 923 ! === Open file and write the header for the SCR file === 924 ! * For the moment only master works. 925 926 if (my_rank==master) then 927 if (iomode == IO_MODE_ETSF) then 928 ofname = nctk_ncify(ofname) 929 NCF_CHECK(nctk_open_create(unt_dump, ofname, xmpi_comm_self)) 930 else 931 if (open_file(ofname,msg,newunit=unt_dump,form="unformatted",status="unknown",action="write") /= 0) then 932 ABI_ERROR(msg) 933 end if 934 end if 935 call wrtout(std_out,sjoin('mkdump_Er: calculating and writing epsilon^-1 matrix on file: ',ofname)) 936 937 ! Update the entries in the header that have been modified. 938 ! TODO, write function to return title, just for info 939 call Er%Hscr%copy(Hscr_cp) 940 Hscr_cp%ID = id_required 941 Hscr_cp%ikxc = ikxc_required 942 Hscr_cp%test_type = option_test 943 Hscr_cp%titles(1) = 'SCR file: epsilon^-1' 944 Hscr_cp%titles(2) = 'TESTPARTICLE' 945 ! Treat the case in which a smaller matrix is used. 946 Hscr_cp%npwe = npwe 947 948 rdwr=2; fform=Hscr_cp%fform 949 call hscr_io(hscr_cp,fform,rdwr,unt_dump,comm_self,master,iomode) 950 call Hscr_cp%free() 951 952 ABI_MALLOC_OR_DIE(epsm1,(npwe,npwe,Er%nomega), ierr) 953 954 do iqibz=1,Er%nqibz 955 is_qeq0=0 956 if (normv(Er%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1 957 ! FIXME there's a problem with SUSC files and MPI-IO 958 !if (iomode == IO_MODE_MPI) then 959 ! ABI_WARNING("SUSC files is buggy. Using Fortran IO") 960 ! call read_screening(in_varname,Er%fname,npwe,1,Er%nomega,epsm1,IO_MODE_FORTRAN,comm_self,iqiA=iqibz) 961 !else 962 call read_screening(in_varname,Er%fname,npwe,1,Er%nomega,epsm1,iomode,comm_self,iqiA=iqibz) 963 !end if 964 965 dim_wing=0; if (is_qeq0==1) dim_wing=3 966 ABI_MALLOC(dummy_lwing,(npwe*Er%nI,Er%nomega,dim_wing)) 967 ABI_MALLOC(dummy_uwing,(npwe*Er%nJ,Er%nomega,dim_wing)) 968 ABI_MALLOC(dummy_head,(dim_wing,dim_wing,Er%nomega)) 969 970 if (approx_type<2 .or. approx_type>3) then 971 ABI_WARNING('Entering out-of core RPA or Kxc branch') 972 call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,& 973 & approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,& 974 & dummy_lwing,dummy_uwing,epsm1,spectra,comm_self) 975 else 976 ABI_WARNING('Entering out-of core fxc_ADA branch') 977 call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,& 978 & approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,& 979 & dummy_lwing,dummy_uwing,epsm1,spectra,comm_self,fxc_ADA(:,:,iqibz)) 980 end if 981 982 ABI_FREE(dummy_head) 983 ABI_FREE(dummy_uwing) 984 ABI_FREE(dummy_lwing) 985 986 if (is_qeq0==1) then 987 call spectra%repr(msg) 988 call wrtout([std_out, ab_out], msg) 989 end if 990 call spectra%free() 991 992 call write_screening(out_varname,unt_dump,iomode,npwe,Er%nomega,iqibz,epsm1) 993 end do 994 995 if (iomode == IO_MODE_ETSF) then 996 NCF_CHECK(nf90_close(unt_dump)) 997 else 998 close(unt_dump) 999 endif 1000 1001 ABI_FREE(epsm1) 1002 end if !master 1003 1004 ! Master broadcasts ofname. 1005 ! NOTE: A synchronization is required here, else the other procs start to read the 1006 ! SCR file before it is written by the master. xmpi_bcast will synch the procs. 1007 call xmpi_bcast(ofname, master, comm, ierr) 1008 1009 ! Now Er% "belongs" to the file "ofname", thus 1010 ! each proc has to destroy and re-initialize the object. 1011 call em1results_free(Er) 1012 1013 mqmem_=Er%mqmem; npwe_asked=npwe 1014 call init_Er_from_file(Er,ofname,mqmem_,npwe_asked,comm) 1015 1016 !Now Er% has been reinitialized and ready-to-use. 1017 Er%id = id_required 1018 call em1results_print(Er) 1019 else 1020 ! ======================== 1021 ! === In-core solution === 1022 ! ======================== 1023 ABI_MALLOC_OR_DIE(Er%epsm1,(npwe,npwe,Er%nomega,Er%nqibz), ierr) 1024 1025 ! FIXME there's a problem with SUSC files and MPI-IO 1026 !if (iomode == IO_MODE_MPI) then 1027 ! !call wrtout(std_out, "read_screening with MPI_IO") 1028 ! ABI_WARNING("SUSC files is buggy. Using Fortran IO") 1029 ! call read_screening(in_varname,Er%fname,npwe,Er%nqibz,Er%nomega,Er%epsm1,IO_MODE_FORTRAN,comm) 1030 !else 1031 call read_screening(in_varname,Er%fname,npwe,Er%nqibz,Er%nomega,Er%epsm1,iomode,comm) 1032 !end if 1033 1034 do iqibz=1,Er%nqibz 1035 is_qeq0=0; if (normv(Er%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) is_qeq0=1 1036 1037 dim_wing=0; if (is_qeq0==1) dim_wing=3 ! FIXME 1038 ABI_MALLOC(dummy_lwing,(npwe*Er%nI,Er%nomega,dim_wing)) 1039 ABI_MALLOC(dummy_uwing,(npwe*Er%nJ,Er%nomega,dim_wing)) 1040 ABI_MALLOC(dummy_head,(dim_wing,dim_wing,Er%nomega)) 1041 1042 if (approx_type<2 .or. approx_type>3) then 1043 ABI_WARNING('Entering in-core RPA and Kxc branch') 1044 call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,& 1045 & approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,& 1046 & dummy_lwing,dummy_uwing,Er%epsm1(:,:,:,iqibz),spectra,comm) 1047 else 1048 ABI_WARNING('Entering in-core fxc_ADA branch') 1049 call make_epsm1_driver(iqibz,dim_wing,npwe,Er%nI,Er%nJ,Er%nomega,Er%omega,& 1050 & approx_type,option_test,Vcp,nfftot,ngfft,nkxc,kxcg,gvec,dummy_head,& 1051 & dummy_lwing,dummy_uwing,Er%epsm1(:,:,:,iqibz),spectra,comm,fxc_ADA=fxc_ADA(:,:,iqibz)) 1052 end if 1053 1054 ABI_FREE(dummy_lwing) 1055 ABI_FREE(dummy_uwing) 1056 ABI_FREE(dummy_head) 1057 1058 if (is_qeq0==1) then 1059 call spectra%repr(msg) 1060 call wrtout([std_out, ab_out], msg) 1061 end if 1062 1063 call spectra%free() 1064 end do 1065 1066 Er%id = id_required 1067 call em1results_print(Er) 1068 end if 1069 end if 1070 1071 DBG_EXIT("COLL") 1072 1073 end subroutine mkdump_Er
m_screening/mkem1_q0 [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
mkem1_q0
FUNCTION
This routine construct the microscopic dieletric matrix for q-->0 starting from the heads, wings and the body of the irreducible polarizability. Afterwards it calculates the symmetrized inverse dieletric matrix via a block wise inversion thus obtaining the heads and the wings of e^{-1} that can be used to describe the non-analytic behavior for q-->0.
INPUTS
npwe=Number of Gs used to describe chi0 nomega=Number of frequencies in chi0. n1,n2=Factors used to define the same of the chi0 matrix (1,1 if collinear, the typical case) Cryst<crystal_t>=Structure describing the crystal structure. Vcp<vcoul_t>=datatypes gathering info on the Coulomb term gvec(3,npwe)=G-vector for chi0 in reduced coordinates. comm=MPI communicator
OUTPUT
eps_head(3,3,nomega)=The macroscopic dieletric tensor in reduced coordinates. The dieletric matrix along versor \hat q can be obtained with e(\hat q) = \hat q.eps_head \hat q if all quantities are given in Cartesian coordinates.
SIDE EFFECTS
chi0(npwe*n1,npwe*n2,nomega)= Input: polarizability. output: inverse dieletric matrix (only the body is used) chi0_lwing(npwe*n1,nomega,3) chi0_uwing(npwe*n2,nomega,3) Input: the lower and upper wings of the polarizability Output: the "lower" and "upper" wings of the inverse dieletric matrix. See notes below. chi0_head(3,3,nomega)= Input: the polarizability tensor in Cartesian coordinates. Ouput: The "head" of the inverse dieletric matrix. See notes below.
NOTES
Matrix inversion in block form. 1 n-1 M = | c u^t| 1 ==> M^{-1} = | 1/k -u^t A^{-1}/k | | v A | n-1 | -A^{-1} v/k A^{-1} + (A^{-1}v u^t A^{-1})/k | where k = c - u^t A^{-1} v Let q be a versor in reciprocal space, the symmetrized dielectric matrix with bare coulomb interaction can be written as \tilde\epsilon = | q.Eq q.Wl(G2) | where E_ij = \delta_ij -4\pi chi0_head_ij | q.Wl(G1) B(G1,G2 | Wl(G1) = -4\pi chi0_lwing(G1) Wu(G2) = -4\pi chi0_uwing(G1) therefore, in Cartesian coordinates, we have: 1) e^{-1}_{00}(q) = [ q_i q_j (E_{ij} - \sum_{GG'} Wu_i(G)a B_{GG'}^{-1} Wl_j(G')) ]^{-1} = 1/(q.Lq) 2) e^{-1}_{0G'}(q) = -e^{-1}_{00}(q) [ \sum_{iG} q_i Wu_i(G)a B_{GG'}^{-1} ] = (q.Su) /(q.Lq) 3) e^{-1}_{G0}(q) = -e^{-1}_{00}(q) [ \sum_{iG'} q_i B_{GG'}^{-1} Wl_i(G') ] = (q.Sl) /(q.Lq) 4) e^{-1}_{GG'}(q) = B_{GG'}^{-1} + [ \sum_{ij} q_i q_j ( \sum_T B^{-1}_{GT}^{-1} Wl_i(T)) (\sum_T' Wu_j(T') B^{-1}_{T'G'}^{-1} ] / (q.Lq) where Su(G,3) and Sl(G,3) are the "upper" and "lower" wings of the inverse dielectric matrix and L is the inverse dielectric tensor. Similar equations hold even if vectors and tensors are given in terms of the reciprocal lattice vectors provided that the metric tensor is taken into account. The main difference is in the expression for the tensor as only one metric tensor can be absorbed in the scalar product, the second metric multiplies one of the wings. *) The present implementation assumes that no cutoff technique is used in the Coulomb term. *) Once the tensor in know it is possible to average the quadratic form on the sphere exactly. In Cartesian coordinates one obtains. \dfrac{1}{4\pi} \int v.Tv d\Omega = Trace(T)/3 For the inverse dielectric matrix we have to resort to a numerical integration
SOURCE
2472 subroutine mkem1_q0(npwe,n1,n2,nomega,Cryst,Vcp,gvec,chi0_head,chi0_lwing,chi0_uwing,chi0,eps_head,comm) 2473 2474 !Arguments ------------------------------------ 2475 !scalars 2476 integer,intent(in) :: npwe,nomega,n1,n2,comm 2477 type(crystal_t),intent(in) :: Cryst 2478 type(vcoul_t),intent(in) :: Vcp 2479 !arrays 2480 integer,intent(in) :: gvec(3,npwe) 2481 complex(gwpc),intent(in) :: chi0(npwe*n1,npwe*n2,nomega) 2482 complex(dpc),intent(inout) :: chi0_lwing(npwe*n1,nomega,3) 2483 complex(dpc),intent(inout) :: chi0_uwing(npwe*n2,nomega,3) 2484 complex(dpc),intent(inout) :: chi0_head(3,3,nomega) 2485 complex(dpc),intent(out) :: eps_head(3,3,nomega) 2486 2487 !Local variables ------------------------------ 2488 !scalars 2489 integer :: iw,ig,ig1,ig2,idir,jdir 2490 !arrays 2491 real(dp),allocatable :: modg_inv(:) 2492 complex(dpc),allocatable :: eps_lwing(:,:),eps_uwing(:,:),eps_body(:,:),cvec(:) 2493 2494 !************************************************************************ 2495 2496 ABI_CHECK(npwe/=1,"npwe must be >1") 2497 ABI_UNUSED(comm) 2498 2499 ! Precompute 1/|G|. 2500 ABI_MALLOC(modg_inv,(npwe-1)) 2501 do ig=1,npwe-1 2502 modg_inv(ig) = one/normv(gvec(:,ig+1),Cryst%gmet,'G') 2503 end do 2504 2505 ABI_MALLOC(eps_uwing,((npwe-1)*n1,3)) 2506 ABI_MALLOC(eps_lwing,((npwe-1)*n2,3)) 2507 ABI_MALLOC(eps_body,(npwe-1,npwe-1)) 2508 ABI_MALLOC(cvec,(npwe-1)) 2509 2510 do iw=1,nomega 2511 ! 2512 ! Head and wings of the symmetrized epsilon. 2513 eps_head(:,:,iw) = -four_pi*chi0_head(:,:,iw) 2514 do idir=1,3 2515 eps_head(idir,idir,iw) = one + eps_head(idir,idir,iw) 2516 eps_lwing(:,idir) = -four_pi * modg_inv * chi0_lwing(2:,iw,idir) 2517 eps_uwing(:,idir) = -four_pi * modg_inv * chi0_uwing(2:,iw,idir) 2518 !eps_lwing(:,idir) = -chi0_lwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1) 2519 !eps_uwing(:,idir) = -chi0_uwing(2:,iw,idir) * SQRT(four_pi) * Vcp%vcqlwl_sqrt(2:npwe,1) 2520 end do 2521 2522 write(std_out,*)" espilon head" 2523 call print_arr(eps_head(:,:,iw)) 2524 ! 2525 ! Construct the body of the symmetrized epsilon then invert it. 2526 do ig2=1,npwe-1 2527 do ig1=1,npwe-1 2528 eps_body(ig1,ig2) = -four_pi * modg_inv(ig1)*chi0(ig1+1,ig2+1,iw )*modg_inv(ig2) 2529 !eps_body(ig1,ig2) = -Vcp%vcqlwl_sqrt(ig1+1,1)*chi0(ig1+1,ig2+1,iw)* Vcp%vcqlwl_sqrt(ig2+1,1) 2530 end do 2531 eps_body(ig2,ig2) = one + eps_body(ig2,ig2) 2532 end do 2533 2534 call xginv(eps_body,npwe-1) 2535 ! 2536 ! Overwrite chi0_head and chi0_wings with the head and the wings of the inverse dielectric matrix. 2537 do jdir=1,3 2538 ! 2539 ! Head. 2540 cvec=czero 2541 do idir=1,3 2542 cvec = cvec + two_pi*Cryst%gmet(jdir,idir)*MATMUL(eps_body,eps_lwing(:,idir)) ! as we work in reciprocal coords. 2543 end do 2544 !cvec = MATMUL(eps_body,eps_lwing(:,jdir)) 2545 do idir=1,3 2546 chi0_head(idir,jdir,iw) = eps_head(idir,jdir,iw) - xdotu(npwe-1,eps_uwing(:,idir),1,cvec,1) 2547 end do 2548 ! 2549 ! Now the wings. 2550 chi0_uwing(2:,iw,jdir) = -MATMUL(eps_uwing(:,jdir),eps_body) 2551 chi0_lwing(2:,iw,jdir) = -MATMUL(eps_body,eps_lwing(:,jdir)) 2552 ! 2553 end do !jdir 2554 2555 call wrtout(std_out, "espilon^1 head after block inversion") 2556 call print_arr(chi0_head(:,:,iw)) 2557 ! 2558 ! Change the body but do not add the corrections due to the head and the wings. 2559 ! since they can be obtained on the fly from eps_body and the wings of eps^{-1}. 2560 !%chi0(2:,2:,iw) = eps_body 2561 end do !iw 2562 2563 ABI_FREE(modg_inv) 2564 ABI_FREE(cvec) 2565 ABI_FREE(eps_lwing) 2566 ABI_FREE(eps_uwing) 2567 ABI_FREE(eps_body) 2568 2569 RETURN 2570 ABI_UNUSED(Vcp%ng) 2571 2572 end subroutine mkem1_q0
m_screening/rpa_symepsm1 [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
rpa_symepsm1
FUNCTION
Calculate RPA $\tilde\epsilon^{-1}$ Use a special treatment of non-Analytic behavior of heads and wings in reciprocal space calculating these quantities for different small q-directions specified by the user (Not yet operative)
INPUTS
iqibz=index of the q-point in the array Vcp%qibz where epsilon^-1 has to be calculated Vcp<vcoul_t>=Structure gathering data on the Coulombian interaction npwe=Number of G-vectors in chi0. nI,nJ=Number of rows/columns in chi0_ij (1,1 in collinear case) dim_wing=Dimension of the wings (0 or 3 if q-->0) chi0_head(dim_wing,dim_wing)=Head of of chi0 (only for q-->0) chi0_lwing(npwe*nI,dim_wing)=Lower wings of chi0 (only for q-->0) chi0_uwing(npwe*nJ,dim_wing)=Upper wings of chi0 (only for q-->0) comm=MPI communicator.
SIDE EFFECTS
chi0(npwe*nI,npwe*nJ): in input the irreducible polarizability, in output the symmetrized inverse dielectric matrix.
SOURCE
1935 subroutine rpa_symepsm1(iqibz,Vcp,npwe,nI,nJ,chi0,my_nqlwl,dim_wing,chi0_head,chi0_lwing,chi0_uwing,epsm_lf,epsm_nlf,eelf,comm) 1936 1937 !Arguments ------------------------------------ 1938 !scalars 1939 integer,intent(in) :: iqibz,nI,nJ,npwe,dim_wing,my_nqlwl,comm 1940 type(vcoul_t),target,intent(in) :: Vcp 1941 !arrays 1942 complex(gwpc),intent(inout) :: chi0(npwe*nI,npwe*nJ) 1943 complex(dpc),intent(inout) :: chi0_lwing(:,:) !(npwe*nI,dim_wing) 1944 complex(dpc),intent(inout) :: chi0_uwing(:,:) !(npwe*nJ,dim_wing) 1945 complex(dpc),intent(inout) :: chi0_head(:,:) !(dim_wing,dim_wing) 1946 real(dp),intent(out) :: eelf(my_nqlwl) 1947 complex(dpc),intent(out) :: epsm_lf(my_nqlwl),epsm_nlf(my_nqlwl) 1948 1949 !Local variables------------------------------- 1950 !scalars 1951 integer,parameter :: master=0,prtvol=0 1952 integer :: ig1,ig2,iqlwl,my_rank,nprocs 1953 real(dp) :: ucvol 1954 logical :: is_qeq0 1955 !character(len=500) :: msg 1956 !arrays 1957 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 1958 complex(gwpc), ABI_CONTIGUOUS pointer :: vc_sqrt(:) 1959 complex(gwpc),allocatable :: chi0_save(:,:) 1960 1961 ! ************************************************************************* 1962 1963 ABI_UNUSED(chi0_head(1,1)) 1964 1965 if (nI/=1.or.nJ/=1) then 1966 ABI_ERROR("nI or nJ=/1 not yet implemented") 1967 end if 1968 1969 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1970 1971 call metric(gmet,gprimd,-1,rmet,Vcp%rprimd,ucvol) 1972 1973 is_qeq0 = (normv(Vcp%qibz(:,iqibz),gmet,'G')<GW_TOLQ0) 1974 if (is_qeq0) then 1975 ABI_CHECK(iqibz==1, "q is 0 but iq_ibz /= 1") 1976 end if 1977 1978 if (my_nqlwl>1) then 1979 ABI_MALLOC(chi0_save,(npwe*nI,npwe*nJ)) 1980 chi0_save = chi0 1981 end if 1982 ! 1983 ! Symmetrized RPA epsilon: 1 - Vc^{1/2} chi0 Vc^{1/2} 1984 ! vc_sqrt contains vc^{1/2}(q, G) 1985 ! 1986 ! Loop over small q"s (if any) to treat the nonanalytical behavior. 1987 do iqlwl=my_nqlwl,1,-1 1988 1989 if (my_nqlwl>1) then 1990 chi0(:,:) = chi0_save ! restore pristine polarizability 1991 chi0(:,1) = chi0_lwing(:,iqlwl) ! change the wings 1992 chi0(1,:) = chi0_uwing(:,iqlwl) 1993 end if 1994 1995 if (iqibz==1) then 1996 vc_sqrt => Vcp%vcqlwl_sqrt(:,iqlwl) ! Use Coulomb term for q-->0 1997 else 1998 vc_sqrt => Vcp%vc_sqrt(:,iqibz) 1999 end if 2000 2001 do ig2=1,npwe*nJ 2002 do ig1=1,npwe*nI 2003 chi0(ig1,ig2) = -vc_sqrt(ig1) * chi0(ig1,ig2) * vc_sqrt(ig2) 2004 end do 2005 chi0(ig2,ig2) = one + chi0(ig2,ig2) 2006 end do 2007 ! chi0, now contains \tepsilon. 2008 2009 epsm_nlf(iqlwl)=chi0(1,1) 2010 2011 if (prtvol > 0) then 2012 call wrtout(std_out,' Symmetrical epsilon(G,G'') ') 2013 call print_arr(chi0, unit=std_out) 2014 end if 2015 ! 2016 ! === Invert tepsilon and calculate macroscopic dielectric constant === 2017 ! * epsm_lf(w) = 1 / epsm1(G=0,Gp=0,w). 2018 ! * Since G=Gp=0, there is no difference btw symmetrical and not symmetrical. 2019 2020 call xginv(chi0,npwe,comm=comm) 2021 2022 epsm_lf(iqlwl) = one/chi0(1,1) 2023 eelf(iqlwl) = -AIMAG(chi0(1,1)) 2024 2025 if (prtvol > 0) then 2026 call wrtout(std_out," Symmetrical epsilon^-1(G,G'')") 2027 call print_arr(chi0, unit=std_out) 2028 end if 2029 ! 2030 ! Save wings of e^-1 overwriting input values. 2031 if (dim_wing>0.and..FALSE.) then 2032 chi0_lwing(:,iqlwl) = chi0(:,1) 2033 chi0_uwing(:,iqlwl) = chi0(1,:) 2034 end if 2035 2036 end do !iqlwl 2037 2038 ABI_SFREE(chi0_save) 2039 2040 end subroutine rpa_symepsm1
m_screening/screen_mdielf [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
screen_mdielf
FUNCTION
Calculates W_{G,G'}(q,w) for a given q-point in the BZ using a model dielectric function.
INPUTS
iq_bz=The index of the q-point in the BZ where W(q) is calculated. npw=Number of plane waves for W nomega=Number of frequency points. model_type=Flag defining the model. eps_inf=Dielectric constant of the material. Cryst<crystal_t>=Info on the unit cell Qmesh<kmesh_t>=Info on the set of q-points. Vcp<vcoul_t datatype>= containing information on the cutoff technique Gsph<Gsphere>=The G-sphere for W. nspden=Number of spin density components of the density. nfft=Number of FFT points on the dense FFT mesh ngfft(18)=contain all needed information about 3D FFT. rhor(nfft,nspden)=Electron density in real space (The PAW AE term is included) which= Set it to "EM1" if the symmetrized inverse dielectric matrix is wanted. By default the routines returns W. comm=MPI communicator.
OUTPUT
w_qbz(npw,npw,nomega)
NOTES
W_{G1,G2} = 1/2 { v(q+G1) \int em1(|q+G1|,r) e^{-i(G1-G2).r} dr + v(q+G2) \int em1(|q+G2|,r) e^{-i(G1-G2).r} dr } / \Omega
SOURCE
2868 subroutine screen_mdielf(iq_bz,npw,nomega,model_type,eps_inf,Cryst,Qmesh,Vcp,Gsph,nspden,nfft,ngfft,rhor,which,w_qbz,comm) 2869 2870 !Arguments ------------------------------------ 2871 !scalars 2872 integer,intent(in) :: npw,nomega,nfft,nspden,iq_bz,comm,model_type 2873 real(dp),intent(in) :: eps_inf 2874 character(len=*),intent(in) :: which 2875 type(kmesh_t),intent(in) :: Qmesh 2876 type(crystal_t),intent(in) :: Cryst 2877 type(vcoul_t),target,intent(in) :: Vcp 2878 type(gsphere_t),intent(in) :: Gsph 2879 !arrays 2880 integer,intent(in) :: ngfft(18) 2881 real(dp),intent(in) :: rhor(nfft,nspden) 2882 complex(gwpc),intent(out) :: w_qbz(npw,npw,nomega) 2883 2884 !Local variables------------------------------- 2885 !scalars 2886 integer,parameter :: tim_fourdp0=0,paral_kgb0=0,cplex1=1 2887 integer :: my_gstart,my_gstop,iq_ibz,ig,itim_q,isym_q 2888 integer :: ig1,ig2,g1mg2_fft,iw,ii,ierr,nprocs,isg,ifft !,row ,col 2889 real(dp) :: qpg2_nrm 2890 complex(dpc) :: ph_mqbzt 2891 logical :: is_qeq0,isirred 2892 !character(len=500) :: msg 2893 type(MPI_type) :: MPI_enreg_seq 2894 !arrays 2895 integer :: umklp(3) 2896 integer,allocatable :: igfft(:),g1mg2(:,:) 2897 real(dp) :: qpg2(3),qpt_bz(3) 2898 real(dp),allocatable :: em1_qpg2r(:),fofg(:,:) 2899 complex(gwpc),ABI_CONTIGUOUS pointer :: vc_sqrt_ibz(:) 2900 complex(gwpc),allocatable :: vc_qbz(:),ctmp(:,:) 2901 logical,allocatable :: mask(:) 2902 2903 ! ************************************************************************* 2904 2905 ABI_CHECK(nomega==1,"screen_mdielf does not support nomega>1") 2906 2907 !Fake MPI_type for the sequential part. 2908 call initmpi_seq(MPI_enreg_seq) 2909 call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfft(2),ngfft(3),'all') 2910 2911 nprocs = xmpi_comm_size(comm) 2912 call xmpi_split_work(npw,comm,my_gstart,my_gstop) 2913 2914 call qmesh%get_bz_item(iq_bz,qpt_bz,iq_ibz,isym_q,itim_q,ph_mqbzt,umklp,isirred) 2915 2916 !if (itim_q/=1.or.isym_q/=1.or.ANY(umklp/=0) ) then 2917 ! ABI_ERROR("Bug in mdielf_bechstedt") 2918 !end if 2919 ! 2920 ! Symmetrize Vc in the full BZ. 2921 is_qeq0 = (normv(qpt_bz,Cryst%gmet,'G')<GW_TOLQ0) ! Check if q==0 2922 if (is_qeq0) then 2923 vc_sqrt_ibz => Vcp%vcqlwl_sqrt(:,1) ! Use Coulomb term for q-->0, only first Q is used, shall we average if nqlwl>1? 2924 else 2925 vc_sqrt_ibz => Vcp%vc_sqrt(:,iq_ibz) 2926 end if 2927 2928 ABI_MALLOC(vc_qbz,(npw)) 2929 do ig=1,npw 2930 isg = Gsph%rottb(ig,itim_q,isym_q) 2931 vc_qbz(isg) = vc_sqrt_ibz(ig)**2 2932 end do 2933 2934 ABI_MALLOC(igfft,(npw)) 2935 ABI_MALLOC(g1mg2,(3,npw)) 2936 ABI_MALLOC(fofg,(2,nfft)) 2937 ABI_MALLOC(em1_qpg2r,(nfft)) 2938 ABI_MALLOC(mask,(npw)) 2939 2940 w_qbz=czero 2941 do ig2=my_gstart,my_gstop 2942 ! 2943 ! Compute the index of G-G2 wave in the FFT grid. 2944 do ii=1,npw 2945 g1mg2(:,ii) = Gsph%gvec(:,ii) - Gsph%gvec(:,ig2) 2946 end do 2947 call kgindex(igfft,g1mg2,mask,MPI_enreg_seq,ngfft,npw) 2948 2949 ! TODO can use zero-padding FFT to speed up the transform. 2950 !call sphereboundary(gbound,istwfk1,g1mg2,mgfft,npw) 2951 2952 ! Evaluate em1_qpg2r = \int em1(|q+G2|,r) e^{-i(G1-G2).r} dr }. 2953 qpg2 = qpt_bz + Gsph%gvec(:,ig2) 2954 qpg2_nrm = normv(qpg2,Cryst%gmet,"G") 2955 2956 do iw=1,nomega 2957 ! 2958 select case (model_type) 2959 case (1) 2960 do ifft=1,nfft 2961 em1_qpg2r(ifft) = one / mdielf_bechstedt(eps_inf,qpg2_nrm,rhor(ifft,1)) 2962 end do 2963 case default 2964 ABI_ERROR(sjoin("Unknown model_type:",itoa(model_type))) 2965 end select 2966 2967 call fourdp(cplex1,fofg,em1_qpg2r,-1,MPI_enreg_seq,nfft,1,ngfft,tim_fourdp0) 2968 ! 2969 ! Here, unlike the other parts of the code, the unsymmetrized e^{-1} is used. 2970 do ig1=1,npw 2971 g1mg2_fft = igfft(ig1) 2972 w_qbz(ig1,ig2,iw) = DCMPLX(fofg(1,g1mg2_fft), fofg(2,g1mg2_fft)) * vc_qbz(ig2) !/ Cryst%ucvol 2973 end do 2974 end do ! iw 2975 end do ! ig2 2976 2977 ABI_FREE(em1_qpg2r) 2978 ABI_FREE(fofg) 2979 ABI_FREE(igfft) 2980 ABI_FREE(g1mg2) 2981 ABI_FREE(mask) 2982 ! 2983 ! W = 1/2 * (A + A^H) 2984 ! The MPI sum is done inside the loop to avoid problems with the size of the packet. 2985 ABI_MALLOC_OR_DIE(ctmp,(npw,npw), ierr) 2986 2987 do iw=1,nomega 2988 !ctmp = TRANSPOSE(CONJG(w_qbz(:,:,iw))) 2989 ctmp = GWPC_CONJG(w_qbz(:,:,iw)) 2990 call sqmat_itranspose(npw,ctmp) 2991 w_qbz(:,:,iw) = half * (ctmp + w_qbz(:,:,iw)) 2992 call xmpi_sum(w_qbz(:,:,iw),comm,ierr) 2993 end do 2994 ! 2995 ! Calculate the symmetrized Em1. W = vc(G1)^{1/2} \tilde Em1 vc(G2)^{1/2} ------------------------- 2996 if (toupper(which)=="EM1") then 2997 do ig=1,npw 2998 isg = Gsph%rottb(ig,itim_q,isym_q) 2999 vc_qbz(isg) = vc_sqrt_ibz(ig) ! Workspace storing vc*{1/2}(q_BZ,G). 3000 end do 3001 3002 do ig2=1,npw 3003 do ig1=1,npw 3004 ctmp(ig1,ig2) = one / (vc_qbz(ig1) * vc_qbz(ig2)) 3005 end do 3006 end do 3007 3008 do iw=1,nomega 3009 w_qbz(:,:,iw) = w_qbz(:,:,iw) * ctmp(:,:) 3010 end do 3011 end if 3012 3013 call destroy_mpi_enreg(MPI_enreg_seq) 3014 3015 ABI_FREE(vc_qbz) 3016 if (allocated(ctmp)) then 3017 ABI_FREE(ctmp) 3018 end if 3019 3020 end subroutine screen_mdielf
m_screening/ylmstar_over_qTq [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
ylmstar_over_qTq
FUNCTION
Return Ylm(q)^*/(q.Tq) where q is a versor in Cartesian coordinates. and Ylm is a complex spherical Harmonics whose index (l,m) are passed via int_pars(1:2). T is a tensore in Cartesian coordinates passed via cplx_pars(1:9).
INPUTS
cart_vers(3)=Cartesian components of the versor int_pars(1:2)=(l,m) indices in Ylm. l>=0 and m \in [-l,l] cplx_pars(1:9)=Tensor T in Cartesian coordinates. real_pars=Not used.
OUTPUT
Value of Ylm(q)^*/(q.Tq)
SOURCE
2705 function ylmstar_over_qTq(cart_vers,int_pars,real_pars,cplx_pars) 2706 2707 !Arguments ------------------------------------ 2708 !scalars 2709 real(dp),intent(in) :: cart_vers(3) 2710 integer,intent(in) :: int_pars(:) 2711 real(dp),intent(in) :: real_pars(:) 2712 complex(dpc),intent(in) :: cplx_pars(:) 2713 complex(dpc) :: ylmstar_over_qTq 2714 !arrays 2715 2716 !Local variables------------------------------- 2717 !scalars 2718 integer :: ll,mm 2719 !arrays 2720 complex(dpc) :: tensor(3,3) 2721 2722 ! ************************************************************************* 2723 2724 tensor = RESHAPE(cplx_pars(1:9),(/3,3/)) 2725 ll = int_pars(1) ! ll starts from zero. 2726 mm = int_pars(2) ! m \in [-l,l] 2727 2728 ylmstar_over_qTq = CONJG(ylmc(ll,mm,cart_vers))/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers)) 2729 2730 RETURN 2731 ABI_UNUSED(real_pars(1)) 2732 2733 end function ylmstar_over_qTq
m_screening/ylmstar_wtq_over_qTq [ Functions ]
[ Top ] [ m_screening ] [ Functions ]
NAME
ylmstar_wtq_over_qTq
FUNCTION
Return Ylm(q)^* weight(q)/(q.Tq) where q is a versor in Cartesian coordinates. Ylm is a complex spherical Harmonics whose index (l,m) are passed via int_pars(1:2). T is a tensor in Cartesian coordinates passed via cplx_pars(1:9). weight(q) is the weighting function giving the length of the vector parallel to versor q that connects the origin of the lattice to one of the boundaries of the small cell centered at Gamma
INPUTS
cart_vers(3)=Cartesian components of the versor int_pars(1:2)=(l,m) indices in Ylm. l>=0 and m \in [-l,l] cplx_pars(1:9)=Tensor T in Cartesian coordinates. real_pars(1:9)=The Cartesian vectors defining the small box centered around gamma point when referred to this vectors the points in the box are given by {(x,y,z) | x,y,z \in [-1,1]}.
OUTPUT
Value of Ylm(q)^* weigh(q)/(q.Tq)
SOURCE
2762 function ylmstar_wtq_over_qTq(cart_vers,int_pars,real_pars,cplx_pars) 2763 2764 !Arguments ------------------------------------ 2765 !scalars 2766 real(dp),intent(in) :: cart_vers(3) 2767 integer,intent(in) :: int_pars(:) 2768 real(dp),intent(in) :: real_pars(:) 2769 complex(dpc),intent(in) :: cplx_pars(:) 2770 complex(dpc) :: ylmstar_wtq_over_qTq 2771 !arrays 2772 2773 !Local variables------------------------------- 2774 !scalars 2775 integer :: ll,mm 2776 real(dp) :: wtq 2777 !arrays 2778 real(dp) :: gprimd(3,3),rprimd(3,3),red_vers(3) 2779 complex(dpc) :: tensor(3,3) 2780 2781 ! ************************************************************************* 2782 2783 ABI_ERROR("Work in progress") 2784 ! box_len has to be tested 2785 2786 gprimd = RESHAPE(real_pars(1:9),(/3,3/)) 2787 red_vers = MATMUL(rprimd,cart_vers) 2788 wtq = box_len(red_vers,gprimd) 2789 2790 tensor = RESHAPE(cplx_pars(1:9),(/3,3/)) 2791 ll = int_pars(1) ! true ll i.e. not shifted 2792 mm = int_pars(2) 2793 2794 ylmstar_wtq_over_qTq = CONJG(ylmc(ll,mm,cart_vers))*wtq/DOT_PRODUCT(cart_vers,MATMUL(tensor,cart_vers)) 2795 2796 end function ylmstar_wtq_over_qTq