Branch data Line data Source code
1 : : !!****m* ABINIT/m_exc_spectra
2 : : !! NAME
3 : : !! m_exc_spectra
4 : : !!
5 : : !! FUNCTION
6 : : !! Routines to compute the macroscopic dielectric function in the Bethe-Salpeter code.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2009-2026 ABINIT and EXC groups (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi, Y. Gillet)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : MODULE m_exc_spectra
23 : :
24 : : use, intrinsic :: iso_c_binding
25 : : use defs_basis
26 : : use m_bs_defs
27 : : use m_abicore
28 : : use m_xmpi
29 : : use m_errors
30 : : use netcdf
31 : : use m_nctk
32 : : use m_ebands
33 : : use m_hdr
34 : :
35 : : use defs_datatypes, only : pseudopotential_type
36 : : use m_io_tools, only : open_file
37 : : use m_fstrings, only : toupper, strcat, sjoin, int2char4
38 : : use m_numeric_tools, only : simpson_int, simpson_cplx
39 : : use m_hide_blas, only : xdotu,xdotc
40 : : use m_special_funcs, only : gaussian
41 : : use m_crystal, only : crystal_t
42 : : use m_bz_mesh, only : kmesh_t
43 : : use m_eprenorms, only : eprenorms_t, renorm_bst
44 : : use m_pawtab, only : pawtab_type
45 : : use m_paw_hr, only : pawhur_t
46 : : use m_wfd, only : wfdgw_t
47 : : !use m_bse_io, only : exc_amplitude
48 : : use m_wfd_optic, only : calc_optical_mels
49 : :
50 : : implicit none
51 : :
52 : : private
53 : :
54 : : public :: build_spectra ! Driver routine for the computation of optical spectra.
55 : : public :: exc_write_data ! This routine drives the writing of the files produced by the Bethe-Salpeter code.
56 : : public :: exc_eps_rpa ! Build epsilon within RPA and GW.
57 : : public :: mdfs_ncwrite ! Writes the MDF.nc file with the final results.
58 : : public :: exc_write_tensor ! Write of complex dielectric tensor
59 : : !public :: exc_eps_resonant ! Build the macroscopic dielectric function with excitonic effects.
60 : : !!***
61 : :
62 : : contains
63 : :
64 : : !!****f* m_exc_spectra/build_spectra
65 : : !! NAME
66 : : !! build_spectra
67 : : !!
68 : : !! FUNCTION
69 : : !! Driver routine for the computation of optical spectra.
70 : : !!
71 : : !! INPUTS
72 : : !! usepaw=1 for PAW calculations, 0 otherwise.
73 : : !! drude_plsmf=Drude plasma frequency.
74 : : !! Bsp<excparam>=Data type gathering the paramenters used for the Bethe-Salpeter calculation.
75 : : !! inclvkb=If different from 0, [Vnl,r] is included in the calculation of the matrix elements of the velocity operator.
76 : : !! BS_files<excfiles>=filenames used in the Bethe-Salpeter part.
77 : : !! Kmesh<kmesh_t>=the k-point sampling for the wave functions.
78 : : !! Cryst<crystal_t>=Structure defining the crystalline structure.
79 : : !! KS_BSt=The KS energies.
80 : : !! QP_BSt=The QP energies.
81 : : !! Psps <pseudopotential_type>=variables related to pseudopotentials.
82 : : !! Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data
83 : : !! Hur(Cryst%natom*usepaw)<pawhur_t>=Only for PAW and DFT+U, quantities used to evaluate the commutator [H_u,r].
84 : : !! Wfd<wfdgw_t>=Handler for the wavefunctions.
85 : : !! nsppol=Number of independent spin polarizations.
86 : : !! nspinor=Number of spinorial components.
87 : : !! comm=MPI communicator.
88 : : !!
89 : : !! OUTPUT
90 : : !! No output. The routine calls specialized routines where the computation and the output of the spectra is done.
91 : : !!
92 : : !! SOURCE
93 : :
94 : 7 : subroutine build_spectra(BSp,BS_files,Cryst,Kmesh,KS_BSt,QP_BSt,Psps,Pawtab,Wfd,Hur,drude_plsmf,comm,Epren)
95 : :
96 : : !Arguments ------------------------------------
97 : : !scalars
98 : : integer,intent(in) :: comm
99 : : real(dp),intent(in) :: drude_plsmf
100 : : type(excparam),intent(in) :: BSp
101 : : type(excfiles),intent(in) :: BS_files
102 : : type(pseudopotential_type),intent(in) :: Psps
103 : : type(kmesh_t),intent(in) :: Kmesh
104 : : type(crystal_t),intent(in) :: Cryst
105 : : type(ebands_t),intent(in) :: KS_BSt,QP_BSt
106 : : type(wfdgw_t),intent(inout) :: Wfd
107 : : type(eprenorms_t),optional,intent(in) :: Epren
108 : : !arrays
109 : : type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
110 : : type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw)
111 : :
112 : : !Local variables ------------------------------
113 : : !scalars
114 : : integer :: my_rank,master,iq,io,nsppol,lomo_min,max_band,ncid
115 : : integer :: itemp,ntemp
116 : : logical :: do_ep_renorm
117 : : real(dp) :: omegaev
118 : : complex(dp) :: ks_avg,gw_avg,exc_avg
119 : : character(len=4) :: ts
120 : : character(len=fnlen) :: path,prefix,filbseig, ost_fname
121 : : !character(len=500) :: msg
122 : 7 : type(ebands_t) :: EPBSt, EP_QPBSt
123 : : !arrays
124 : 7 : real(dp),allocatable :: dos_exc(:),dos_gw(:),dos_ks(:)
125 : 7 : complex(dp),allocatable :: eps_rpanlf(:,:),eps_gwnlf(:,:), eps_exc(:,:),opt_cvk(:,:,:,:,:)
126 : : !************************************************************************
127 : :
128 : 7 : my_rank = Wfd%my_rank
129 : 7 : master = Wfd%master
130 : 7 : nsppol = Wfd%nsppol
131 : :
132 : 7 : do_ep_renorm = .False.
133 : 7 : ntemp = 1
134 [ + + + + ]: 7 : if (BSp%do_ep_renorm .and. PRESENT(Epren)) then
135 : 1 : do_ep_renorm = .True.
136 : 1 : ntemp = Epren%ntemp
137 : : end if
138 : :
139 : : ! =====================================================
140 : : ! === Calculate fcv(k)=<c k s|e^{-iqr}|v k s> in BZ ===
141 : : ! =====================================================
142 : 7 : lomo_min=Bsp%lomo_min; max_band=Bsp%nbnds
143 [ - + + - : 49 : ABI_MALLOC(opt_cvk,(lomo_min:max_band,lomo_min:max_band,BSp%nkbz,nsppol,BSp%nq))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + ]
144 : :
145 [ + + ]: 49 : do iq=1,BSp%nq
146 : : call calc_optical_mels(Wfd,Kmesh,KS_BSt,Cryst,Psps,Pawtab,Hur,BSp%inclvkb,Bsp%lomo_spin,lomo_min,max_band,&
147 : 49 : BSp%nkbz,BSp%q(:,iq),opt_cvk(:,:,:,:,iq))
148 : : end do
149 : : !
150 : : ! ============================
151 : : ! ==== Make EPS EXCITONIC ====
152 : : ! ============================
153 [ + - ]: 7 : if (my_rank==master) then ! Only master works.
154 : :
155 [ - + + - : 28 : ABI_MALLOC(eps_exc,(BSp%nomega,BSp%nq))
- + + - +
- + - - +
- + ]
156 [ - + + - : 21 : ABI_MALLOC(dos_exc,(BSp%nomega))
+ - + - -
+ - + ]
157 [ - + + - : 21 : ABI_MALLOC(eps_rpanlf,(BSp%nomega,BSp%nq))
- + + - -
+ - + ]
158 [ - + + - : 14 : ABI_MALLOC(dos_ks,(BSp%nomega))
- + - + ]
159 [ - + + - : 21 : ABI_MALLOC(eps_gwnlf ,(BSp%nomega,BSp%nq))
- + + - -
+ - + ]
160 [ - + + - : 14 : ABI_MALLOC(dos_gw,(BSp%nomega))
- + - + ]
161 : :
162 [ + + ]: 23 : do itemp = 1, ntemp
163 : 16 : call int2char4(itemp,ts)
164 : :
165 [ + + ]: 16 : if(do_ep_renorm) then
166 : 10 : prefix = TRIM("_T") // ts
167 : : else
168 : 6 : prefix = ""
169 : : end if
170 : :
171 [ + - ]: 16 : ost_fname = strcat(BS_files%out_basename,prefix,"_EXC_OST")
172 : :
173 : : !TODO for RPA
174 : 16 : call KS_BST%copy(EPBSt)
175 : 16 : call QP_BST%copy(EP_QPBSt)
176 : :
177 [ - + ]: 16 : if (BS_files%in_eig /= BSE_NOFILE) then
178 [ # # ]: 0 : filbseig = strcat(BS_files%in_eig,prefix)
179 : : else
180 [ + - ]: 16 : filbseig = strcat(BS_files%out_eig,prefix)
181 : : end if
182 : :
183 [ + + ]: 16 : if(do_ep_renorm) then
184 : : ! No scissor with KSBST
185 : 10 : call renorm_bst(Epren, EPBSt, Cryst, itemp, do_lifetime=.TRUE.,do_check=.TRUE.)
186 : :
187 : 10 : call renorm_bst(Epren, EP_QPBSt, Cryst, itemp, do_lifetime=.TRUE.,do_check=.FALSE.)
188 : : end if
189 : :
190 : :
191 [ + + ]: 16 : if (BSp%use_coupling==0) then
192 : : call exc_eps_resonant(BSp,filbseig,ost_fname,lomo_min,max_band,BSp%nkbz,nsppol,opt_cvk,&
193 : 15 : Cryst%ucvol,BSp%nomega,BSp%omega,eps_exc,dos_exc,elph_lifetime=do_ep_renorm)
194 : : else
195 : : call exc_eps_coupling(Bsp,BS_files,lomo_min,max_band,BSp%nkbz,nsppol,opt_cvk,&
196 : 1 : Cryst%ucvol,BSp%nomega,BSp%omega,eps_exc,dos_exc)
197 : : end if
198 : : !
199 : : ! =======================================================
200 : : ! === Make EPS RPA and GW without local-field effects ===
201 : : ! =======================================================
202 : 16 : call wrtout(std_out," Calculating RPA NLF and QP NLF epsilon")
203 : :
204 : : call exc_eps_rpa(BSp%nbnds,BSp%lomo_spin,Bsp%lomo_min,BSp%homo_spin,Kmesh,EPBSt,BSp%nq,nsppol,opt_cvk,&
205 : 16 : Cryst%ucvol,BSp%broad,BSp%nomega,BSp%omega,eps_rpanlf,dos_ks)
206 : :
207 : : call exc_eps_rpa(BSp%nbnds,BSp%lomo_spin,Bsp%lomo_min,BSp%homo_spin,Kmesh,EP_QPBSt,BSp%nq,nsppol,opt_cvk,&
208 : 16 : Cryst%ucvol,Bsp%broad,BSp%nomega,BSp%omega,eps_gwnlf,dos_gw)
209 : : !
210 : : ! =========================
211 : : ! === Write out Epsilon ===
212 : : ! =========================
213 : : !this is just for the automatic tests, It will be removed when fldiff
214 : : !will be able to compare two optical spectral
215 : 16 : write(ab_out,*)" "
216 : 16 : write(ab_out,*)"Macroscopic dielectric function:"
217 : 16 : write(ab_out,*)"omega [eV] <KS_RPA_nlf> <GW_RPA_nlf> <BSE> "
218 [ + + ]: 176 : do io=1,MIN(10,BSp%nomega)
219 : 160 : omegaev = REAL(BSp%omega(io))*Ha_eV
220 [ + + ]: 1120 : ks_avg = SUM( eps_rpanlf(io,:)) / Bsp%nq
221 [ + + ]: 1120 : gw_avg = SUM( eps_gwnlf (io,:)) / Bsp%nq
222 [ + + ]: 1120 : exc_avg = SUM( eps_exc (io,:)) / Bsp%nq
223 : 176 : write(ab_out,'(7f9.4)')omegaev,ks_avg,gw_avg,exc_avg
224 : : end do
225 : 16 : write(ab_out,*)" "
226 : :
227 : : !
228 : : ! Master node writes final results on file.
229 : 16 : call exc_write_data(BSp,BS_files,"RPA_NLF_MDF",eps_rpanlf,prefix=prefix,dos=dos_ks)
230 : 16 : call exc_write_data(BSp,BS_files,"GW_NLF_MDF",eps_gwnlf,prefix=prefix,dos=dos_gw)
231 : 16 : call exc_write_data(BSp,BS_files,"EXC_MDF",eps_exc,prefix=prefix,dos=dos_exc)
232 : :
233 : 16 : call wrtout(std_out," Checking Kramers Kronig on Excitonic Macroscopic Epsilon")
234 [ + - + + : 14240 : call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_exc(:,1))
- + ]
235 : :
236 : 16 : call wrtout(std_out," Checking Kramers Kronig on RPA NLF Macroscopic Epsilon")
237 [ + + - + ]: 14240 : call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_rpanlf(:,1))
238 : :
239 : 16 : call wrtout(std_out," Checking Kramers Kronig on GW NLF Macroscopic Epsilon")
240 [ + + - + ]: 14240 : call check_kramerskronig(BSp%nomega,REAL(BSp%omega),eps_gwnlf(:,1))
241 : :
242 : 16 : call wrtout(std_out," Checking f-sum rule on Excitonic Macroscopic Epsilon")
243 : :
244 [ + - ]: 16 : if (BSp%exchange_term>0) then
245 : 16 : ABI_COMMENT(' f-sum rule should be checked without LF')
246 : : end if
247 [ + + + - : 28464 : call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_exc(:,1)),drude_plsmf)
+ + - + -
+ ]
248 : :
249 : 16 : call wrtout(std_out," Checking f-sum rule on RPA NLF Macroscopic Epsilon")
250 [ + + + + : 28464 : call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_rpanlf(:,1)),drude_plsmf)
- + - + ]
251 : :
252 : 16 : call wrtout(std_out," Checking f-sum rule on GW NLF Macroscopic Epsilon")
253 [ + + + + : 28464 : call check_fsumrule(BSp%nomega,REAL(BSp%omega),AIMAG(eps_gwnlf(:,1)),drude_plsmf)
- + - + ]
254 : :
255 [ + - ]: 16 : path = strcat(BS_files%out_basename, strcat(prefix,"_MDF.nc"))
256 [ - + ]: 16 : NCF_CHECK_MSG(nctk_open_create(ncid, path, xmpi_comm_self), sjoin("Creating MDF file:", path))
257 [ - + ]: 16 : NCF_CHECK(cryst%ncwrite(ncid))
258 [ - + ]: 16 : NCF_CHECK(QP_BSt%ncwrite(ncid))
259 : : ! Write dielectric functions.
260 : 16 : call mdfs_ncwrite(ncid, Bsp, eps_exc,eps_rpanlf,eps_gwnlf)
261 [ - + ]: 16 : NCF_CHECK(nf90_close(ncid))
262 : :
263 : : !TODO
264 : 16 : call EPBSt%free()
265 : 23 : call EP_QPBSt%free()
266 : : end do
267 : :
268 [ - + ]: 7 : ABI_FREE(eps_rpanlf)
269 [ - + ]: 7 : ABI_FREE(eps_gwnlf)
270 [ - + ]: 7 : ABI_FREE(eps_exc)
271 [ - + ]: 7 : ABI_FREE(dos_exc)
272 [ - + ]: 7 : ABI_FREE(dos_ks)
273 [ - + ]: 7 : ABI_FREE(dos_gw)
274 : : end if ! my_rank==master
275 : :
276 [ - + ]: 7 : ABI_FREE(opt_cvk)
277 : :
278 : 7 : call xmpi_barrier(comm)
279 : :
280 [ - + - + : 7 : end subroutine build_spectra
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
281 : : !!***
282 : :
283 : : !----------------------------------------------------------------------
284 : :
285 : : !!****f* m_exc_spectra/exc_write_data
286 : : !! NAME
287 : : !! exc_write_data
288 : : !!
289 : : !! FUNCTION
290 : : !! This routine drives the writing of the files produced by the Bethe-Salpeter code.
291 : : !!
292 : : !! INPUTS
293 : : !! BSp<excparam>=Bethe-Salpeter Parameters.
294 : : !! what= "EXC_MDF"
295 : : !! "RPA_NLF_MDF"
296 : : !! "GW_NLF_MDF"
297 : : !! [dos(nomega)]
298 : : !!
299 : : !! OUTPUT
300 : : !! Only writing.
301 : : !!
302 : : !! SIDE EFFECTS
303 : : !! eps(BSp%nomega,BSp%nq) = Macroscopic dielectric function to be written.
304 : : !!
305 : : !! SOURCE
306 : :
307 : 132 : subroutine exc_write_data(BSp,BS_files,what,eps,prefix,dos)
308 : :
309 : : !Arguments ------------------------------------
310 : : !scalars
311 : : character(len=*),intent(in) :: what
312 : : type(excparam),intent(in) :: BSp
313 : : type(excfiles),intent(in) :: BS_files
314 : : character(len=*),optional,intent(in) :: prefix
315 : : !arrays
316 : : real(dp),optional,intent(in) :: dos(BSp%nomega)
317 : : complex(dp),intent(in) :: eps(BSp%nomega,BSp%nq)
318 : :
319 : : !Local variables ------------------------------
320 : : !scalars
321 : : integer :: io,iq,funt
322 : : real(dp) :: omegaev,step
323 : : !real(dp),parameter :: SMALL=5.0d-99
324 : : !arrays
325 : 264 : real(dp) :: int_dos(BSp%nomega)
326 : 264 : real(dp) :: tmp_eps(2,BSp%nq)
327 : : character(len=500) :: lf_type,block_type,wgg_type,frm,str_type,msg
328 : : character(len=fnlen) :: fname
329 : : !************************************************************************
330 : :
331 [ + - ]: 132 : if (PRESENT(prefix)) then
332 [ + - ]: 132 : fname = strcat(BS_files%out_basename,prefix,'_',toupper(what))
333 : : else
334 [ # # ]: 0 : fname = strcat(BS_files%out_basename,'_',toupper(what))
335 : : end if
336 : :
337 [ - + ]: 132 : if (open_file(fname,msg,newunit=funt,form="formatted", action="write") /= 0) then
338 : 0 : ABI_ERROR(msg)
339 : : end if
340 : :
341 : 308 : select case (toupper(what))
342 : : case ("EXC_MDF")
343 [ + - ]: 44 : call wrtout(ab_out," Writing EXC Macroscopic dielectric function to file: "//trim(fname))
344 : :
345 : 44 : write(funt,'("# Macroscopic dielectric function obtained with the BS equation.")')
346 : :
347 : 44 : lf_type = 'WITHOUT LOCAL FIELD EFFECTS'
348 [ + - ]: 44 : if (BSp%exchange_term>0) lf_type='LOCAL FIELD EFFECTS INCLUDED'
349 : 44 : call bsp%calctype2str(str_type)
350 : 44 : write(funt,'("# ",a," " ,a)') TRIM(str_type), TRIM(lf_type)
351 : :
352 : 44 : block_type = 'RESONANT-ONLY calculation'
353 [ + + ]: 44 : if (BSp%use_coupling>0) block_type = 'RESONANT+COUPLING calculation'
354 : 44 : write(funt,'("# ",a)') TRIM(block_type)
355 : :
356 [ + - ]: 44 : if (BSp%use_coulomb_term) then
357 : 44 : wgg_type = "Coulomb term constructed with full W(G1,G2)"
358 [ + + ]: 44 : if ( BSp%use_diagonal_Wgg ) wgg_type = "Coulomb term constructed with diagonal approximation W(G1,G1)"
359 : 44 : write(funt,'("# ",a)') TRIM(wgg_type)
360 : : end if
361 : :
362 : 44 : write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
363 : :
364 : : case ("RPA_NLF_MDF")
365 [ + - ]: 44 : call wrtout(ab_out," Writing KS-RPA macroscopic dielectric function without local fields to file: "//trim(fname))
366 : 44 : write(funt,'("# RPA macroscopic dielectric function without local fields")')
367 : :
368 : : case ("GW_NLF_MDF")
369 [ + - ]: 44 : call wrtout(ab_out," Writing GW-RPA macroscopic dielectric function without local fields to file: "//trim(fname))
370 : :
371 : 44 : write(funt,'("# GW Macroscopic dielectric function without local field effects ")')
372 : 44 : write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
373 : :
374 : : case default
375 [ + + + - : 132 : ABI_ERROR("Unknown value for what: "//trim(what))
- - ]
376 : : end select
377 : :
378 : : ! Paramaters common to the different calculations.
379 [ + + ]: 132 : if (BSp%algorithm /= BSE_ALGO_HAYDOCK) then
380 : 48 : write(funt,'(a,i0)')"# nstates included in the diagonalization = ",BSp%nstates
381 : : end if
382 : :
383 [ + + ]: 132 : if (BSp%algorithm == BSE_ALGO_HAYDOCK) then
384 : 84 : write(funt,'(a,2f7.4)')'# Tolerance = ',BSp%haydock_tol
385 : : end if
386 : :
387 : 132 : write(funt,'(a,i0)')"# npweps = ",BSp%npweps
388 : 132 : write(funt,'(a,i0)')"# npwwfn = ",BSp%npwwfn
389 : 132 : write(funt,'(a,i0)')"# nbands = ",BSp%nbnds
390 : 132 : write(funt,'(a,i0)')"# loband = ",BSp%lomo_spin(1)
391 [ + + ]: 132 : if (Bsp%nsppol==2) write(funt,'(a,i0)')"# loband(spin=2) = ",BSp%lomo_spin(2)
392 : 132 : write(funt,'(a,i0)')"# nkibz = ",BSp%nkibz
393 : 132 : write(funt,'(a,i0)')"# nkbz = ",BSp%nkbz
394 : 132 : write(funt,'(a,f7.4,a)')'# Lorentzian broadening = ',BSp%broad*Ha_eV,' [eV]'
395 : :
396 : : ! Write the list of q-points.
397 : 132 : write(funt,'(a)')"# List of q-points for the optical limit:"
398 [ + + ]: 924 : do iq=1,BSp%nq
399 : 924 : write(funt,'(a,3(f9.6,","),a)')'# q = ',BSp%q(:,iq),' [Reduced coords] '
400 : : end do
401 : :
402 : : ! Write spectra.
403 [ + + ]: 132 : if (.not.PRESENT(dos)) then
404 : 28 : write(funt,'(a)')"# omega [eV] RE(eps(q=1)) IM(eps(q=1) RE(eps(q=2) ) ... "
405 : 28 : write(frm,*)'(f7.3,',2*BSp%nq,'(1x,f9.4))'
406 [ + + ]: 13349 : do io=1,BSp%nomega
407 : 13321 : omegaev = DBLE(BSp%omega(io))*Ha_eV
408 [ + + ]: 93247 : tmp_eps(1,:) = REAL (eps(io,:))
409 [ + + ]: 93247 : tmp_eps(2,:) = AIMAG(eps(io,:))
410 : : !where (ABS(tmp_eps) < SMALL) ! this to improve the portability of the automatic tests.
411 : : ! tmp_eps = zero
412 : : !end where
413 [ + - + + ]: 93275 : write(funt,frm) omegaev,(tmp_eps(:,iq), iq=1,BSp%nq)
414 : : end do
415 : :
416 : : else
417 : 104 : write(funt,'(a)')"# omega [eV] RE(eps(q=1)) IM(eps(q=1) RE(eps(q=2) ) ... DOS IDOS"
418 : 104 : step = DBLE(BSp%omega(2) - BSp%omega(1))
419 [ - + ]: 104 : if ( ABS( step - DBLE((BSp%omega(BSp%nomega) - BSp%omega(BSp%nomega-1)))) > tol6 ) then
420 : 0 : ABI_WARNING("Frequency mesh must be linear for using simpson_int")
421 : : end if
422 : 104 : call simpson_int(Bsp%nomega,step,dos,int_dos)
423 : : !write(frm,*)'(f7.3,',2*BSp%nq,'es12.4,2es12.4)'
424 : 104 : write(frm,*)'(f7.3,',2*BSp%nq,'(1x,f9.4,1x,f9.4,1x,f9.4))'
425 [ + + ]: 69418 : do io=1,BSp%nomega
426 : 69314 : omegaev = DBLE(BSp%omega(io))*Ha_eV
427 [ + + ]: 485198 : tmp_eps(1,:) = REAL (eps(io,:))
428 [ + + ]: 485198 : tmp_eps(2,:) = AIMAG(eps(io,:))
429 : : !where (ABS(tmp_eps) < SMALL) ! this to improve the portability of the automatic tests.
430 : : ! tmp_eps = zero
431 : : !end where
432 : : !write(funt,frm) omegaev,(eps(io,iq), iq=1,BSp%nq), dos(io), int_dos(io)
433 [ + - + + ]: 485302 : write(funt,frm) omegaev,(tmp_eps(:,iq), iq=1,BSp%nq), dos(io), int_dos(io)
434 : : end do
435 : : end if
436 : :
437 : 132 : close(funt)
438 : :
439 : 132 : end subroutine exc_write_data
440 : : !!***
441 : :
442 : : !----------------------------------------------------------------------
443 : :
444 : : !!****f* m_exc_spectra/exc_eps_rpa
445 : : !! NAME
446 : : !! exc_eps_rpa
447 : : !!
448 : : !! FUNCTION
449 : : !! Build epsilon within RPA and GW.
450 : : !!
451 : : !! INPUTS
452 : : !! nkbz=Number of points in the BZ
453 : : !! nbnds=Number of bands
454 : : !! lomo_spin(nsppol)
455 : : !! lomo_min=Lowest occupied state
456 : : !! homo=Number of occupied states.
457 : : !! homo_spin(nsppol)
458 : : !! nsppol=Number of independent spin polarizations.
459 : : !! nomega=Number of frequencies
460 : : !! omega(nomega)=Frequency mesh.
461 : : !! ucvol=Unit cell volume.
462 : : !! broad=Broadening used for the DOS.
463 : : !! opt_cvk(nbnds,nbnds,nkbz)=Matrix elements <b k|e^{-iqr}|b" k> for a given q in the full BZ.
464 : : !!
465 : : !! OUTPUT
466 : : !! eps_rpa(nomega)=RPA spectrum without local-field effects.
467 : : !! dos(nomega)=The DOS.
468 : : !!
469 : : !! SOURCE
470 : :
471 : 96 : subroutine exc_eps_rpa(nbnds,lomo_spin,lomo_min,homo_spin,Kmesh,Bst,nq,nsppol,opt_cvk,ucvol,broad,nomega,omega,eps_rpa,dos)
472 : :
473 : : !Arguments ------------------------------------
474 : : !scalars
475 : : integer,intent(in) :: nbnds,lomo_min,nsppol,nomega,nq
476 : : real(dp),intent(in) :: ucvol,broad
477 : : type(kmesh_t),intent(in) :: Kmesh
478 : : type(ebands_t),intent(in) :: BSt
479 : : !arrays
480 : : integer,intent(in) :: lomo_spin(nsppol),homo_spin(nsppol)
481 : : real(dp),intent(out) :: dos(nomega)
482 : : complex(dp),intent(in) :: omega(nomega)
483 : : complex(dp),intent(in) :: opt_cvk(lomo_min:nbnds,lomo_min:nbnds,Kmesh%nbz,nsppol,nq)
484 : : complex(dp),intent(out) :: eps_rpa(nomega,nq)
485 : :
486 : : !Local variables ------------------------------
487 : : !scalars
488 : : integer :: iw,ib_v,ib_c,ik_bz,ik_ibz,spin,iq
489 : : real(dp) :: fact,arg,ediff
490 : : real(dp) :: linewidth
491 : : complex(dp) :: ctemp
492 : : logical :: do_linewidth
493 : : !************************************************************************
494 : :
495 : : ! TODO: four_pi comes from the bare Coulomb term hence the
496 : : ! present implementation is not compatible with the cutoff technique.
497 : 96 : fact=four_pi/(ucvol*Kmesh%nbz)
498 [ + + ]: 96 : if (nsppol==1) fact=two*fact ! two accounts for the occupation factors.
499 : :
500 [ + + + + : 403158 : eps_rpa=czero; dos=zero
+ + ]
501 : :
502 : 96 : do_linewidth = .FALSE.
503 : 96 : do_linewidth = allocated(BSt%linewidth)
504 : :
505 : : !write(std_out,*)nsppol,Kmesh%nbz,lomo_min,homo,nbnds
506 : : !
507 : : ! Sum over all QP transitions.
508 [ + + ]: 194 : do spin=1,nsppol
509 [ + + ]: 3314 : do ik_bz=1,Kmesh%nbz
510 : 3120 : ik_ibz = Kmesh%tab(ik_bz)
511 [ + + ]: 12466 : do ib_v=lomo_spin(spin),homo_spin(spin)
512 [ + + ]: 42192 : do ib_c=homo_spin(spin)+1,nbnds
513 : : !
514 : : ! TODO here energies are always assumed to be real.
515 : 29824 : ediff = BSt%eig(ib_c,ik_ibz,spin) - BSt%eig(ib_v,ik_ibz,spin)
516 : :
517 [ + + ]: 39072 : if(do_linewidth) then
518 : 3840 : linewidth = BSt%linewidth(1,ib_c,ik_ibz,spin) + BSt%linewidth(1,ib_v,ik_ibz,spin)
519 [ + + ]: 26880 : do iq=1,nq
520 : 23040 : ctemp = opt_cvk(ib_c,ib_v,ik_bz,spin,iq)
521 [ + + ]: 23089920 : do iw=1,nomega
522 : : eps_rpa(iw,iq) = eps_rpa(iw,iq) + ctemp * CONJG(ctemp) *&
523 : 23086080 : (one/(ediff-j_dpc*linewidth-omega(iw)) + one/(ediff+j_dpc*linewidth+omega(iw)))
524 : : end do
525 : : end do
526 : : ! The JDOS at q=0
527 : : !if (ediff*Ha_eV < 0.3) then
528 : : ! write(std_out,*)"Small transition ",ik_ibz,ib_v,ib_c
529 : : !end if
530 : :
531 [ + + ]: 3847680 : do iw=1,nomega
532 : 3843840 : arg = DBLE(omega(iw)) - ediff
533 : 3847680 : dos(iw) = dos(iw) + gaussian(arg, linewidth)
534 : : end do
535 : : else
536 [ + + ]: 181888 : do iq=1,nq
537 : 155904 : ctemp = opt_cvk(ib_c,ib_v,ik_bz,spin,iq)
538 [ + + ]: 40669696 : do iw=1,nomega
539 : : eps_rpa(iw,iq) = eps_rpa(iw,iq) + ctemp * CONJG(ctemp) *&
540 : 40643712 : (one/(ediff-omega(iw)) + one/(ediff+omega(iw)))
541 : : end do
542 : : end do
543 : : !
544 : : ! The JDOS at q=0
545 : : !if (ediff*Ha_eV < 0.3) then
546 : : ! write(std_out,*)"Small transition ",ik_ibz,ib_v,ib_c
547 : : !end if
548 : :
549 [ + + ]: 6773952 : do iw=1,nomega
550 : 6747968 : arg = DBLE(omega(iw)) - ediff
551 : 6773952 : dos(iw) = dos(iw) + gaussian(arg, broad)
552 : : end do
553 : : end if
554 : : !
555 : : end do !ib_c
556 : : end do !ib_v
557 : : end do !ik_bz
558 : : end do !spin
559 : :
560 [ + + ]: 57594 : dos = dos/Kmesh%nbz
561 [ + + + + ]: 345660 : eps_rpa = cone + fact*eps_rpa
562 : :
563 : 96 : end subroutine exc_eps_rpa
564 : : !!***
565 : :
566 : : !----------------------------------------------------------------------
567 : :
568 : : !!****f* m_exc_spectra/exc_eps_resonant
569 : : !! NAME
570 : : !! exc_eps_resonant
571 : : !!
572 : : !! FUNCTION
573 : : !! This routine builds the macroscopic dielectric function with excitonic effects.
574 : : !!
575 : : !! INPUTS
576 : : !! Bsp
577 : : !! lomo_min,max_band
578 : : !! nkbz=Number of points in the BZ
579 : : !! nsppol=Number of independent polarizations.
580 : : !! nomega=Number of frequencies
581 : : !! omega(nomega)=frequency mesh (complex shift is already included)
582 : : !! ucvol=Volume of the unit cell.
583 : : !! opt_cvk(lomo_min:max_band,mib:max_band,nkbz,nsppol,Bsp%nq)=Matrix elements <b k|e^{-iqr}|b" k> for a given q in the full BZ.
584 : : !!
585 : : !! OUTPUT
586 : : !! eps_exc(nomega,Bsp%nq)=Macroscopic dielectric function with excitonic effects.
587 : : !! dos_exc(nomega)=The DOS of the excitonic Hamiltonian
588 : : !!
589 : : !! SOURCE
590 : :
591 : 15 : subroutine exc_eps_resonant(Bsp,filbseig,ost_fname,lomo_min,max_band,nkbz,nsppol,opt_cvk,&
592 : 15 : ucvol,nomega,omega,eps_exc,dos_exc,elph_lifetime)
593 : :
594 : : !Arguments ------------------------------------
595 : : !scalars
596 : : integer,intent(in) :: lomo_min,max_band,nkbz,nomega,nsppol
597 : : real(dp),intent(in) :: ucvol
598 : : type(excparam),intent(in) :: BSp
599 : : character(len=fnlen),intent(in) :: filbseig,ost_fname
600 : : logical,optional,intent(in) :: elph_lifetime
601 : : !arrays
602 : : real(dp),intent(out) :: dos_exc(nomega)
603 : : complex(dp),intent(in) :: opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol,BSp%nq),omega(nomega)
604 : : complex(dp),intent(out) :: eps_exc(nomega,BSp%nq)
605 : :
606 : : !Local variables ------------------------------
607 : : !scalars
608 : : integer :: ll,it,iw,ib_v,ib_c,ik_bz,neig_read,eig_unt,exc_size,iq !fform,
609 : : integer :: spin,spad,hsize_read,nstates,ost_unt
610 : : logical :: do_ep_lifetime, file_do_lifetime
611 : : real(dp) :: fact,arg
612 : : complex(dp) :: dotprod
613 : : character(len=500) :: msg,frm,errmsg
614 : : !arrays
615 : 15 : real(dp),allocatable :: exc_ene(:)
616 : 15 : complex(dp) :: ctemp(BSp%nq),dtemp(BSp%nq)
617 : 15 : complex(dp),allocatable :: ostrength(:,:),exc_ene_cplx(:),exc_state(:),exc_state2(:)
618 : : !************************************************************************
619 : :
620 : 15 : call wrtout(std_out," Calculating excitonic epsilon with antiresonant")
621 : :
622 [ + + ]: 15 : if (nsppol==2) then
623 : 1 : ABI_WARNING("nsppol==2 still under development")
624 : : end if
625 : :
626 [ + + ]: 31 : exc_size = SUM(BSp%nreh)
627 : 15 : nstates = BSp%nstates
628 : :
629 : 15 : do_ep_lifetime = .FALSE.
630 [ + - ]: 15 : if (PRESENT(elph_lifetime)) do_ep_lifetime = elph_lifetime
631 : :
632 [ + + + - : 31 : if (ANY(Bsp%nreh/=Bsp%nreh(1))) then
- + ]
633 : 0 : write(msg,'(a,2(i0,1x))')"BSE does not support different number of transitions for the two spin channels. nreh: ",Bsp%nreh
634 : 0 : ABI_WARNING(msg)
635 : : end if
636 : : !
637 : : ! TODO:
638 : : ! four_pi comes from the bare Coulomb term hence the
639 : : ! present implementation is not compatible with the cutoff technique.
640 [ + + ]: 15 : fact=four_pi/(ucvol*nkbz); if (nsppol==1) fact=two*fact ! two to account for the occupation numbers.
641 : :
642 [ + - ]: 15 : call wrtout(std_out," Reading excitonic eigenstates from file: "//TRIM(filbseig))
643 [ - + ]: 15 : if (open_file(filbseig,msg,newunit=eig_unt,form="unformatted",status="old",action="read") /= 0) then
644 : 0 : ABI_ERROR(msg)
645 : : end if
646 : :
647 [ + - ]: 15 : read(eig_unt, err=10, iomsg=errmsg) file_do_lifetime
648 : :
649 [ + + - + ]: 15 : if(do_ep_lifetime .and. .not. file_do_lifetime) then
650 : 0 : ABI_ERROR("Cannot do lifetime as the data is not present in the file !")
651 : : end if
652 : :
653 [ + - ]: 15 : read(eig_unt, err=10, iomsg=errmsg) hsize_read,neig_read
654 : :
655 [ - + ]: 15 : if (hsize_read /= exc_size) then
656 : 0 : write(msg,'(2(a,i0))')" Wrong size of the Hamiltonian: read: ",hsize_read," expected= ",exc_size
657 : 0 : ABI_ERROR(msg)
658 : : end if
659 : :
660 [ - + ]: 15 : if (neig_read /= nstates) then
661 : 0 : write(msg,'(2(a,i0))')" Wrong number of eigenstates: read: ",neig_read," expected= ",nstates
662 : 0 : ABI_ERROR(msg)
663 : : end if
664 : : !
665 : : ! Read eigenvalues, ignore possibly small imaginary part.
666 [ - + + - : 45 : ABI_MALLOC(exc_ene_cplx,(neig_read))
+ - + - -
+ - + ]
667 [ + - ]: 15 : read(eig_unt, err=10, iomsg=errmsg) exc_ene_cplx
668 : :
669 [ - + + - : 45 : ABI_MALLOC(exc_ene,(neig_read))
+ - + - -
+ - + ]
670 [ + - - + : 2258 : exc_ene = DBLE(exc_ene_cplx)
- - + + ]
671 : : !ABI_FREE(exc_ene_cplx)
672 : : !
673 : : ! Calculate oscillator strength.
674 [ - + + - : 45 : ABI_MALLOC(exc_state,(exc_size))
+ - + - -
+ - + ]
675 [ - + + - : 60 : ABI_MALLOC(ostrength,(neig_read,BSp%nq))
- + + - +
- + - - +
- + ]
676 : :
677 [ + + ]: 15 : if (do_ep_lifetime) then
678 [ - + + - : 20 : ABI_MALLOC(exc_state2,(exc_size))
- + - + ]
679 [ + + ]: 970 : do ll=1,neig_read ! Loop over excitonic eigenstates reported on file.
680 [ + - ]: 960 : read(eig_unt, err=10, iomsg=errmsg) exc_state(:) ! Righteigenvector
681 [ + - ]: 960 : read(eig_unt, err=10, iomsg=errmsg) exc_state2(:) ! Lefteigenvector
682 : :
683 : : ! Here assuming that eigenvectors are such as Xl_i' Xr_j = delta_ij
684 : : ! Otherwise, I need to invert the overlap matrix !
685 : :
686 : : ! Rescale the vectors so that they are "normalized" with respect to the other one !
687 : 960 : dotprod = xdotc(exc_size,exc_state2(:),1,exc_state(:),1)
688 [ + + ]: 93120 : exc_state2(:) = exc_state2(:)/CONJG(dotprod)
689 : :
690 [ + + ]: 6720 : ctemp(:) = czero
691 [ + + ]: 6720 : dtemp(:) = czero
692 [ + + ]: 1920 : do spin=1,nsppol
693 : 960 : spad=(spin-1)*BSp%nreh(1) ! Loop over spin channels.
694 [ + + ]: 94080 : do it=1,BSp%nreh(spin) ! Loop over resonant transition t = (k,v,c,s)
695 : 92160 : ik_bz = Bsp%Trans(it,spin)%k
696 : 92160 : ib_v = Bsp%Trans(it,spin)%v
697 : 92160 : ib_c = Bsp%Trans(it,spin)%c
698 [ + + ]: 646080 : do iq=1,BSp%nq
699 : 552960 : ctemp(iq) = ctemp(iq) + CONJG(opt_cvk(ib_c,ib_v,ik_bz,spin,iq)) * exc_state(it+spad)
700 : 645120 : dtemp(iq) = dtemp(iq) + CONJG(opt_cvk(ib_c,ib_v,ik_bz,spin,iq)) * exc_state2(it+spad)
701 : : end do
702 : : end do ! it
703 : : end do
704 [ + + ]: 6730 : ostrength(ll,:) = ctemp(:)*CONJG(dtemp(:))
705 : : end do ! ll
706 [ - + ]: 10 : ABI_FREE(exc_state2)
707 : : else
708 [ + + ]: 1273 : do ll=1,neig_read ! Loop over excitonic eigenstates reported on file.
709 [ + - ]: 1268 : read(eig_unt, err=10, iomsg=errmsg) exc_state(:)
710 [ - + - - ]: 1268 : if(file_do_lifetime) read(eig_unt, err=10, iomsg=errmsg)
711 : :
712 [ + + ]: 8876 : ctemp(:) = czero
713 [ + + ]: 3304 : do spin=1,nsppol
714 : 2036 : spad=(spin-1)*BSp%nreh(1) ! Loop over spin channels.
715 [ + + ]: 757480 : do it=1,BSp%nreh(spin) ! Loop over resonant transition t = (k,v,c,s)
716 : 754176 : ik_bz = Bsp%Trans(it,spin)%k
717 : 754176 : ib_v = Bsp%Trans(it,spin)%v
718 : 754176 : ib_c = Bsp%Trans(it,spin)%c
719 [ + + ]: 5281268 : do iq=1,BSp%nq
720 : 5279232 : ctemp(iq) = ctemp(iq) + CONJG(opt_cvk(ib_c,ib_v,ik_bz,spin,iq)) * exc_state(it+spad)
721 : : end do
722 : : end do ! it
723 : : end do
724 [ + + ]: 8881 : ostrength(ll,:) = ctemp(:)*CONJG(ctemp(:))
725 : : end do ! ll
726 : : end if
727 : :
728 : :
729 [ + - ]: 15 : close(eig_unt, err=10, iomsg=errmsg)
730 [ - + ]: 15 : ABI_FREE(exc_state)
731 : :
732 [ + + ]: 15 : if(do_ep_lifetime) then
733 [ + + + + ]: 60130 : eps_exc = one
734 [ + + ]: 970 : do ll=1,neig_read ! Sum over all excitonic eigenstates read from file.
735 [ + + ]: 6730 : do iq=1,BSp%nq
736 [ + + ]: 5772480 : do iw=1,nomega
737 : : eps_exc(iw,iq) = eps_exc(iw,iq) + &
738 : 5771520 : fact * ostrength(ll,iq) * (one/(exc_ene_cplx(ll) - omega(iw)) - one/(-DCONJG(exc_ene_cplx(ll)) - omega(iw)))
739 : : end do
740 : : end do !ll
741 : : end do !iw
742 : : else
743 [ + + + + ]: 24713 : eps_exc = one
744 [ + + ]: 1273 : do ll=1,neig_read ! Sum over all excitonic eigenstates read from file.
745 [ + + ]: 8881 : do iq=1,BSp%nq
746 [ + + ]: 1458164 : do iw=1,nomega
747 : : eps_exc(iw,iq) = eps_exc(iw,iq) + &
748 : 1456896 : fact * ostrength(ll,iq) * (one/(exc_ene(ll) - omega(iw)) - one/(-exc_ene(ll) - omega(iw)))
749 : : end do
750 : : end do !ll
751 : : end do !iw
752 : : end if
753 : :
754 : : ! The excitonic DOS.
755 [ + + ]: 14138 : dos_exc=zero
756 [ + + ]: 2243 : do ll=1,neig_read ! Sum over the calculate excitonic eigenstates.
757 [ + + ]: 1204751 : do iw=1,nomega
758 : 1202508 : arg = ( DBLE(omega(iw)) - exc_ene(ll))
759 [ + + ]: 1204736 : if(do_ep_lifetime) then
760 : 960960 : dos_exc(iw) = dos_exc(iw) + gaussian(arg, AIMAG(exc_ene_cplx(ll)))
761 : : else
762 : 241548 : dos_exc(iw) = dos_exc(iw) + gaussian(arg, Bsp%broad)
763 : : end if
764 : : end do
765 : : end do
766 : :
767 : : ! Write oscillator strengths to file.
768 [ - + ]: 15 : if (open_file(ost_fname,msg,newunit=ost_unt,form="formatted",action="write") /= 0) then
769 : 0 : ABI_ERROR(msg)
770 : : end if
771 : :
772 : 15 : write(ost_unt,'("# Oscillator strengths of the excitonic states for the different q-polarizations.")')
773 : :
774 : : ! Write the list of q-points.
775 : 15 : write(ost_unt,*)"# List of q-points for the optical limit"
776 [ + + ]: 105 : do iq=1,BSp%nq
777 : 105 : write(ost_unt,'(a,3(f9.6,","),a)')'# q = ',BSp%q(:,iq),' [Reduced coords] '
778 : : end do
779 : :
780 : 15 : write(ost_unt,*)"# E_lambda [eV] ostrength(q=1) ostrength(q=2) .... "
781 : 15 : write(frm,*)'(f8.4,',BSp%nq,'es12.4)'
782 [ + + ]: 2243 : do ll=1,neig_read
783 : 2243 : write(ost_unt,frm)exc_ene(ll)*Ha_eV,(ostrength(ll,iq), iq=1,BSp%nq)
784 : : end do
785 : :
786 : 15 : close(ost_unt)
787 : :
788 [ - + ]: 15 : ABI_FREE(ostrength)
789 [ - + ]: 15 : ABI_FREE(exc_ene)
790 [ - + ]: 15 : ABI_FREE(exc_ene_cplx)
791 : :
792 : : !call exc_amplitude(Bsp,filbseig,1,(/(ll,ll=1,10)/),"TEST_AMPLITUDE")
793 : : !call exc_amplitude(Bsp,filbseig,1,(/30/),"TEST_AMPLITUDE")
794 : 15 : return
795 : :
796 : : ! Handler IO-error
797 : : 10 continue
798 : 0 : ABI_ERROR(errmsg)
799 : :
800 [ + - - - : 15 : end subroutine exc_eps_resonant
- - - - -
- - - ]
801 : : !!***
802 : :
803 : : !----------------------------------------------------------------------
804 : :
805 : : !!****f* m_exc_spectra/exc_eps_coupling
806 : : !! NAME
807 : : !! exc_eps_coupling
808 : : !!
809 : : !! FUNCTION
810 : : !! Make epsilon EXCITONIC with full COUPLING.
811 : : !!
812 : : !! INPUTS
813 : : !! Bsp
814 : : !! nkbz=Number of points in the BZ
815 : : !! lomo_min,max_band
816 : : !! nomega=Number of frequencies
817 : : !! omega(nomega)=frequency mesh.
818 : : !! nsppol=Number of independent spin polarizations.
819 : : !! ucvol=Unit cell volume.
820 : : !! BS_files<excfiles>File names used in the Bethe-Salpeter code.
821 : : !! opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol)=Matrix elements <b k|e^{-iqr}|b" k> for a given q in the full BZ.
822 : : !!
823 : : !! OUTPUT
824 : : !! eps_exc(nomega)=Macroscopic dielectric function with excitonic effects calculated including the COUPLING.
825 : : !! dos_exc(nomega)=The DOS of the excitonic Hamiltonian
826 : : !!
827 : : !! SOURCE
828 : :
829 : 1 : subroutine exc_eps_coupling(Bsp,BS_files,lomo_min,max_band,nkbz,nsppol,opt_cvk,ucvol,nomega,omega,eps_exc,dos_exc)
830 : :
831 : : !Arguments ------------------------------------
832 : : !scalars
833 : : integer,intent(in) :: lomo_min,max_band,nkbz,nomega,nsppol
834 : : real(dp),intent(in) :: ucvol
835 : : type(excfiles),intent(in) :: BS_files
836 : : type(excparam),intent(in) :: BSp
837 : : !arrays
838 : : real(dp),intent(out) :: dos_exc(nomega)
839 : : complex(dp),intent(in) :: opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol,BSp%nq),omega(nomega)
840 : : complex(dp),intent(out) :: eps_exc(nomega,BSp%nq)
841 : :
842 : : !Local variables ------------------------------
843 : : !scalars
844 : : integer :: mi,it,ii,ib_v,ib_c,ik_bz,exc_size_read,nstates_read,eig_unt !,fform
845 : : integer :: exc_size,iq,spin,tr_idx,tar_idx,nstates,iw,ll,ierr
846 : : real(dp) :: fact,arg
847 : : complex(dp) :: eps,fam,famp
848 : : character(len=500) :: msg,errmsg
849 : : character(len=fnlen) :: filbseig
850 : : logical :: do_lifetime
851 : : !arrays
852 : 1 : complex(dp),allocatable :: Ami(:),exc_ene(:),Sm1mi(:)
853 : 1 : complex(dp),allocatable :: msfap(:,:),fa(:,:),fap(:,:)
854 : : !************************************************************************
855 : :
856 : 1 : call wrtout(std_out," Calculating absorption strength with full coupling")
857 : :
858 [ - + ]: 1 : if (nsppol==2) then
859 : 0 : ABI_WARNING("nsppol==2 is still under development")
860 : : end if
861 : :
862 : : ! Rank of the entire excitonic Hamiltonian including the coupling block.
863 [ + + - + : 2 : exc_size = 2*SUM(BSp%nreh); if (nsppol==2) exc_size = 2*(SUM(BSp%nreh) + BSp%nreh(2))
- - ]
864 : 1 : nstates = BSp%nstates
865 : :
866 : : ! TODO: four_pi comes from the bare Coulomb term hence the
867 : : ! present implementation is not compatible with the cutoff technique.
868 : : ! factor two is due to the occupation factors.
869 [ + - ]: 1 : fact=four_pi/(ucvol*nkbz); if (nsppol==1) fact=two*fact
870 : :
871 [ - + ]: 1 : if (BS_files%in_eig /= BSE_NOFILE) then
872 : 0 : filbseig = BS_files%in_eig
873 : : else
874 : 1 : filbseig = BS_files%out_eig
875 : : end if
876 : :
877 [ + - ]: 1 : call wrtout(std_out," Reading excitonic eigenstates from file: "//trim(filbseig))
878 [ - + ]: 1 : if (open_file(filbseig,msg,newunit=eig_unt,form="unformatted", status="old", action="read") /= 0) then
879 : 0 : ABI_ERROR(msg)
880 : : end if
881 : :
882 [ + - ]: 1 : read(eig_unt, err=10, iomsg=errmsg) do_lifetime
883 : :
884 [ - + ]: 1 : if (do_lifetime) then
885 : 0 : ABI_CHECK(.not. do_lifetime, "Finite lifetime with coupling is not supported yet !")
886 : : end if
887 : :
888 [ + - ]: 1 : read(eig_unt, err=10, iomsg=errmsg) exc_size_read, nstates_read
889 [ - + ]: 1 : ABI_CHECK(exc_size_read==exc_size,"wrong file")
890 [ - + ]: 1 : ABI_CHECK(nstates_read==nstates,"Partial diago not supported yet")
891 : : !
892 : : ! Read eigenvalues
893 [ - + + - : 3 : ABI_MALLOC(exc_ene,(nstates))
+ - + - -
+ - + ]
894 [ + - ]: 1 : read(eig_unt, err=10, iomsg=errmsg) exc_ene(:)
895 : :
896 [ - + + - : 4 : ABI_MALLOC(fa,(nstates,BSp%nq))
- + + - +
- + - - +
- + ]
897 [ - + + - : 3 : ABI_MALLOC(fap,(nstates,BSp%nq))
- + + - -
+ - + ]
898 [ - + + - : 3 : ABI_MALLOC_OR_DIE(Ami,(exc_size), ierr)
+ - + - +
- - + -
+ ]
899 : :
900 [ + + ]: 769 : do mi=1,nstates ! Loop on excitonic eigenvalues mi
901 [ + - ]: 768 : read(eig_unt, err=10, iomsg=errmsg) Ami(:)
902 : :
903 [ + + ]: 5377 : do iq=1,BSp%nq
904 : : fam = czero
905 : : famp = czero
906 [ + + ]: 9216 : do spin=1,nsppol
907 [ + + ]: 1778688 : do it=1,BSp%nreh(spin) ! Loop over transition t = (k,v,c)
908 : 1769472 : ik_bz = Bsp%Trans(it,spin)%k
909 : 1769472 : ib_v = Bsp%Trans(it,spin)%v
910 : 1769472 : ib_c = Bsp%Trans(it,spin)%c
911 : 1769472 : tr_idx = it + (spin-1)*Bsp%nreh(1)
912 [ + - ]: 1769472 : if (nsppol==1) then
913 : 1769472 : tar_idx = it + Bsp%nreh(1)
914 : : else
915 [ # # # # ]: 0 : if (spin==1) tar_idx = it + SUM(Bsp%nreh)
916 [ # # ]: 0 : if (spin==2) tar_idx = it + 2*Bsp%nreh(1)+Bsp%nreh(2)
917 : : end if
918 : :
919 : : fam = fam + CONJG(opt_cvk(ib_c,ib_v,ik_bz,spin,iq)) * Ami(tr_idx) &
920 : 1769472 : + CONJG(opt_cvk(ib_v,ib_c,ik_bz,spin,iq)) * Ami(tar_idx)
921 : :
922 : : famp = famp - opt_cvk(ib_c,ib_v,ik_bz,spin,iq) * CONJG(Ami(tr_idx)) &
923 : 1774080 : + opt_cvk(ib_v,ib_c,ik_bz,spin,iq) * CONJG(Ami(tar_idx))
924 : : end do
925 : : end do
926 : : ! Save results.
927 : 4608 : fa (mi,iq) = fam
928 : 5376 : fap(mi,iq) = famp
929 : : end do
930 : : end do ! mi
931 : :
932 [ - + ]: 1 : ABI_FREE(Ami)
933 : :
934 : : ! Read O{-1} and sum over the eigenstates.
935 [ - + + - : 3 : ABI_MALLOC(msfap,(nstates,BSp%nq))
- + + - -
+ - + ]
936 [ - + + - : 2 : ABI_MALLOC(Sm1mi,(nstates))
- + - + ]
937 : :
938 [ + + ]: 769 : do mi=1,nstates
939 [ + - ]: 768 : read(eig_unt, err=10, iomsg=errmsg) Sm1mi
940 [ + + ]: 590592 : Sm1mi = DCONJG(Sm1mi) ! This gives the row since O^{-1} is Hermitian.
941 [ + + ]: 5377 : do iq=1,BSp%nq
942 : 5376 : msfap(mi,iq) = xdotu(exc_size,Sm1mi,1,fap(:,iq),1)
943 : : end do
944 : : end do
945 : :
946 [ - + ]: 1 : ABI_FREE(Sm1mi)
947 : :
948 [ + - ]: 1 : close(eig_unt, err=10, iomsg=errmsg)
949 : :
950 : : ! === Calculate excitonic epsilon with coupling ===
951 [ + + ]: 7 : do iq=1,BSp%nq
952 : : !
953 [ + + ]: 613 : do ii=1,nomega
954 : : eps = czero
955 [ + + ]: 466014 : do mi=1,nstates ! sum over all exciton eigenstates
956 : 466014 : eps = eps - fa(mi,iq) * msfap(mi,iq) / (exc_ene(mi) - omega(ii))
957 : : end do
958 : 612 : eps_exc(ii,iq) = one + fact * eps
959 : : end do
960 : : end do
961 : :
962 [ - + ]: 1 : ABI_FREE(fa)
963 [ - + ]: 1 : ABI_FREE(msfap)
964 [ - + ]: 1 : ABI_FREE(fap)
965 : :
966 : : ! The excitonic DOS.
967 [ + + ]: 102 : dos_exc=zero
968 [ + + ]: 769 : do ll=1,nstates ! Sum over the calculate excitonic eigenstates.
969 [ + + ]: 78337 : do iw=1,nomega
970 : 77568 : arg = DBLE(omega(iw) - exc_ene(ll))
971 : 78336 : dos_exc(iw) = dos_exc(iw) + gaussian(arg, Bsp%broad)
972 : : end do
973 : : end do
974 : :
975 [ - + ]: 1 : ABI_FREE(exc_ene)
976 : :
977 : : return
978 : :
979 : : 10 continue
980 : 0 : ABI_ERROR(errmsg)
981 : :
982 [ # # # # : 0 : end subroutine exc_eps_coupling
# # # # #
# # # ]
983 : : !!***
984 : :
985 : : !----------------------------------------------------------------------
986 : :
987 : : !!****f* m_exc_spectra/exc_write_tensor
988 : : !! NAME
989 : : !! exc_write_tensor
990 : : !!
991 : : !! FUNCTION
992 : : !! This routine drives the writing of complex dielectric tensor
993 : : !!
994 : : !! INPUTS
995 : : !! BSp<excparam>=Bethe-Salpeter Parameters.
996 : : !! what= "EXC_TSR_CART" or "EXC_TSR_RED"
997 : : !! "RPA_NLF_TSR_CART" or "RPA_NLF_TSR_RED"
998 : : !! "GW_NLF_TSR_CART" or "GW_NLF_TSR_RED"
999 : : !!
1000 : : !! OUTPUT
1001 : : !! Only writing.
1002 : : !!
1003 : : !! SIDE EFFECTS
1004 : : !! tensor(BSp%nomega,6) = Complex dielectric tensor to be written
1005 : : !!
1006 : : !! SOURCE
1007 : :
1008 : 168 : subroutine exc_write_tensor(BSp,BS_files,what,tensor)
1009 : :
1010 : : !Arguments ------------------------------------
1011 : : !scalars
1012 : : character(len=*),intent(in) :: what
1013 : : type(excparam),intent(in) :: BSp
1014 : : type(excfiles),intent(in) :: BS_files
1015 : : !arrays
1016 : : complex(dp),intent(in) :: tensor(BSp%nomega,6)
1017 : :
1018 : : !Local variables ------------------------------
1019 : : !scalars
1020 : : integer :: io,iq,funt
1021 : : real(dp) :: omegaev
1022 : : !arrays
1023 : : character(len=500) :: lf_type,block_type,wgg_type,frm,str_type, msg
1024 : : character(len=fnlen) :: fname
1025 : : !************************************************************************
1026 : :
1027 [ + - ]: 168 : fname = strcat(BS_files%out_basename,'_',toupper(what))
1028 [ - + ]: 168 : if (open_file(fname,msg,newunit=funt,form="formatted", action="write") /= 0) then
1029 : 0 : ABI_ERROR(msg)
1030 : : end if
1031 : :
1032 : 364 : select case (toupper(what))
1033 : : case ("EXC_TSR_CART")
1034 : 28 : write(funt,'("# Complex dielectric tensor (cart. coord.) obtained with the BS equation.")')
1035 : :
1036 : 28 : lf_type = 'WITHOUT LOCAL FIELD EFFECTS'
1037 [ + - ]: 28 : if (BSp%exchange_term>0) lf_type='LOCAL FIELD EFFECTS INCLUDED'
1038 : 28 : call bsp%calctype2str(str_type)
1039 : 28 : write(funt,'("# ",a," " ,a)') TRIM(str_type), TRIM(lf_type)
1040 : :
1041 : 28 : block_type = 'RESONANT-ONLY calculation'
1042 [ - + ]: 28 : if (BSp%use_coupling>0) block_type = 'RESONANT+COUPLING calculation'
1043 : 28 : write(funt,'("# ",a)') TRIM(block_type)
1044 : :
1045 [ + - ]: 28 : if (BSp%use_coulomb_term) then
1046 : 28 : wgg_type = "Coulomb term constructed with full W(G1,G2)"
1047 [ + + ]: 28 : if ( BSp%use_diagonal_Wgg ) wgg_type = "Coulomb term constructed with diagonal approximation W(G1,G1)"
1048 : 28 : write(funt,'("# ",a)') TRIM(wgg_type)
1049 : : end if
1050 : :
1051 : 28 : write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
1052 : :
1053 : : case ("EXC_TSR_RED")
1054 : 28 : write(funt,'("# Complex dielectric tensor (red. coord.) obtained with the BS equation.")')
1055 : :
1056 : 28 : lf_type = 'WITHOUT LOCAL FIELD EFFECTS'
1057 [ + - ]: 28 : if (BSp%exchange_term>0) lf_type='LOCAL FIELD EFFECTS INCLUDED'
1058 : 28 : call bsp%calctype2str(str_type)
1059 : 28 : write(funt,'("# ",a," " ,a)') TRIM(str_type), TRIM(lf_type)
1060 : :
1061 : 28 : block_type = 'RESONANT-ONLY calculation'
1062 [ - + ]: 28 : if (BSp%use_coupling>0) block_type = 'RESONANT+COUPLING calculation'
1063 : 28 : write(funt,'("# ",a)') TRIM(block_type)
1064 : :
1065 [ + - ]: 28 : if (BSp%use_coulomb_term) then
1066 : 28 : wgg_type = "Coulomb term constructed with full W(G1,G2)"
1067 [ + + ]: 28 : if ( BSp%use_diagonal_Wgg ) wgg_type = "Coulomb term constructed with diagonal approximation W(G1,G1)"
1068 : 28 : write(funt,'("# ",a)') TRIM(wgg_type)
1069 : : end if
1070 : :
1071 : 28 : write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
1072 : :
1073 : : case ("RPA_NLF_TSR_CART")
1074 : 28 : write(funt,'("# RPA complex dielectric tensor (cart. coord.) without local fields")')
1075 : :
1076 : : case ("RPA_NLF_TSR_RED")
1077 : 28 : write(funt,'("# RPA complex dielectric tensor (red. coord.) without local fields")')
1078 : :
1079 : : case ("GW_NLF_TSR_CART")
1080 : 28 : write(funt,'("# GW complex dielectric tensor (cart. coord.) without local field effects ")')
1081 : 28 : write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
1082 : :
1083 : : case ("GW_NLF_TSR_RED")
1084 : 28 : write(funt,'("# GW complex dielectric tensor (red. coord.) without local field effects ")')
1085 : 28 : write(funt,'(a,f7.4,a)')'# Scissor operator energy = ',BSp%mbpt_sciss*Ha_eV,' [eV]'
1086 : :
1087 : : case default
1088 [ + + + + : 168 : ABI_ERROR("Unknown value for what: "//TRIM(what))
+ + - -
- ]
1089 : : end select
1090 : :
1091 : : ! Paramaters common to the different calculations.
1092 [ - + ]: 168 : if (BSp%algorithm /= BSE_ALGO_HAYDOCK) then
1093 : 0 : write(funt,'(a,i0)')"# nstates included in the diagonalization = ",BSp%nstates
1094 : : end if
1095 : :
1096 [ + - ]: 168 : if (BSp%algorithm == BSE_ALGO_HAYDOCK) then
1097 : 168 : write(funt,'(a,2f7.4)')'# Tolerance = ',BSp%haydock_tol
1098 : : end if
1099 : :
1100 : 168 : write(funt,'(a,i0)')"# npweps = ",BSp%npweps
1101 : 168 : write(funt,'(a,i0)')"# npwwfn = ",BSp%npwwfn
1102 : 168 : write(funt,'(a,i0)')"# nbands = ",BSp%nbnds
1103 : 168 : write(funt,'(a,i0)')"# loband = ",BSp%lomo_spin(1)
1104 [ - + ]: 168 : if (Bsp%nsppol==2) write(funt,'(a,i0)')"# loband(spin=2) = ",BSp%lomo_spin(2)
1105 : 168 : write(funt,'(a,i0)')"# nkibz = ",BSp%nkibz
1106 : 168 : write(funt,'(a,i0)')"# nkbz = ",BSp%nkbz
1107 : 168 : write(funt,'(a,f7.4,a)')'# Lorentzian broadening = ',BSp%broad*Ha_eV,' [eV]'
1108 : :
1109 : : ! Write tensor.
1110 : 168 : write(funt,'(3a)') "# omega [eV] RE(eps_11) IM(eps_11) RE(eps_22)", &
1111 : 168 : "IM(eps_22) RE(eps_33) IM(eps_33) RE(eps_12) IM(eps_12)", &
1112 : 336 : "RE(eps_13) IM(eps_13) RE(eps_23) IM(eps_23))"
1113 : 168 : write(frm,*) '(f7.3,12es14.6)'
1114 [ + + ]: 80094 : do io=1,BSp%nomega
1115 : 79926 : omegaev = DBLE(BSp%omega(io))*Ha_eV
1116 : 80094 : write(funt,frm) omegaev,(tensor(io,iq), iq=1,6)
1117 : : end do
1118 : :
1119 : 168 : close(funt)
1120 : :
1121 : 168 : end subroutine exc_write_tensor
1122 : : !!***
1123 : :
1124 : : !----------------------------------------------------------------------
1125 : :
1126 : : !!****f* m_exc_spectra/mdfs_ncwrite
1127 : : !! NAME
1128 : : !! mdfs_ncwrite
1129 : : !!
1130 : : !! FUNCTION
1131 : : !! Writes the MDF.nc file with the final results.
1132 : : !!
1133 : : !! INPUTS
1134 : : !! ncid =NC file handle
1135 : : !! Bsp<excparam>=Data type gathering the paramenters used for the Bethe-Salpeter calculation.
1136 : : !! eps_exc = Excitonic MDF
1137 : : !! eps_rpanlf = KS-RPA MDF without local-field effects.
1138 : : !! eps_gwnlf = GW-RPA MDF without local-field effects.
1139 : : !!
1140 : : !! OUTPUT
1141 : : !! Only writing.
1142 : : !!
1143 : : !! SOURCE
1144 : :
1145 : 44 : subroutine mdfs_ncwrite(ncid,Bsp,eps_exc,eps_rpanlf,eps_gwnlf)
1146 : :
1147 : : !Arguments ------------------------------------
1148 : : !scalars
1149 : : integer,intent(in) :: ncid
1150 : : type(excparam),intent(in) :: BSp
1151 : : !arrays
1152 : : complex(dp),target,intent(in) :: eps_exc(BSp%nomega,BSp%nq)
1153 : : complex(dp),target,intent(in) :: eps_rpanlf(BSp%nomega,BSp%nq)
1154 : : complex(dp),target,intent(in) :: eps_gwnlf(BSp%nomega,BSp%nq)
1155 : :
1156 : : !Local variables-------------------------------
1157 : : !scalars
1158 : : integer :: ncerr
1159 : 44 : real(dp), contiguous, pointer :: rvals(:,:,:)
1160 : : ! *************************************************************************
1161 : : ! =========================
1162 : : ! === Write the dimensions
1163 : : ! =========================
1164 : :
1165 : : ncerr = nctk_defnwrite_ivars(ncid, [character(len=nctk_slen) :: &
1166 : : "mdf_version", "nsppol", "npwwfn", "npweps", "nkibz", "nkbz",&
1167 : : "nkibz_iterp", "nkbz_interp", "wtype", "interp_mode"],&
1168 : : [1, Bsp%nsppol, Bsp%npwwfn, Bsp%npweps, Bsp%nkibz, Bsp%nkbz, &
1169 [ + + + + ]: 924 : Bsp%nkibz_interp, Bsp%nkbz_interp,Bsp%wtype, Bsp%interp_mode])
1170 [ - + ]: 44 : NCF_CHECK(ncerr)
1171 : :
1172 : : ncerr = nctk_defnwrite_dpvars(ncid, [character(len=nctk_slen) :: &
1173 : : "ecutwfn", "ecuteps", "mbpt_sciss", "broad", "eps_inf"],&
1174 [ + + + + ]: 484 : [Bsp%ecutwfn, Bsp%ecuteps, Bsp%mbpt_sciss, Bsp%broad, Bsp%eps_inf])
1175 [ - + ]: 44 : NCF_CHECK(ncerr)
1176 : :
1177 : : ncerr = nctk_def_dims(ncid, [nctkdim_t("two", 2), nctkdim_t("three", 3), nctkdim_t("number_of_qpoints", Bsp%nq),&
1178 [ + + ]: 264 : nctkdim_t("number_of_frequencies", Bsp%nomega), nctkdim_t("number_of_spins", bsp%nsppol)], defmode=.True.)
1179 [ - + ]: 44 : NCF_CHECK(ncerr)
1180 : :
1181 : : ! Define variables.
1182 : : !arrays
1183 : : ncerr = nctk_def_arrays(ncid, [&
1184 : : nctkarr_t('qpoints', "dp", 'three, number_of_qpoints'),&
1185 : : nctkarr_t('wmesh', "dp", 'number_of_frequencies'),&
1186 : : nctkarr_t('nreh', "i", "number_of_spins"),&
1187 : : nctkarr_t('lomo_spin', "i", "number_of_spins"),&
1188 : : nctkarr_t('humo_spin', "i", "number_of_spins"),&
1189 : : nctkarr_t('exc_mdf', "dp", 'two, number_of_frequencies, number_of_qpoints'),&
1190 : : nctkarr_t('rpanlf_mdf', "dp", 'two, number_of_frequencies, number_of_qpoints'),&
1191 [ + + ]: 396 : nctkarr_t('gwnlf_mdf', "dp", 'two, number_of_frequencies, number_of_qpoints')])
1192 [ - + ]: 44 : NCF_CHECK(ncerr)
1193 : :
1194 : : ! Write data.
1195 [ - + ]: 44 : NCF_CHECK(nctk_set_datamode(ncid))
1196 [ - + ]: 44 : NCF_CHECK(nf90_put_var(ncid, vid('qpoints'), Bsp%q))
1197 [ - + ]: 44 : NCF_CHECK(nf90_put_var(ncid, vid('nreh'), bsp%nreh))
1198 [ - + ]: 44 : NCF_CHECK(nf90_put_var(ncid, vid('lomo_spin'), bsp%lomo_spin))
1199 [ - + ]: 44 : NCF_CHECK(nf90_put_var(ncid, vid('humo_spin'), bsp%humo_spin))
1200 : :
1201 : : ! Write frequency in mesh in eV.
1202 [ + + - + : 27589 : NCF_CHECK(nf90_put_var(ncid, vid('wmesh'), REAL(Bsp%omega)*Ha_eV))
- - ]
1203 : :
1204 [ + + ]: 176 : call c_f_pointer(c_loc(eps_exc(1,1)), rvals, shape=[2, bsp%nomega, bsp%nq])
1205 [ - + ]: 44 : NCF_CHECK(nf90_put_var(ncid, vid('exc_mdf'), rvals))
1206 : :
1207 [ + + ]: 176 : call c_f_pointer(c_loc(eps_rpanlf(1,1)), rvals, shape=[2, bsp%nomega, bsp%nq])
1208 [ - + ]: 44 : NCF_CHECK(nf90_put_var(ncid, vid('rpanlf_mdf'), rvals))
1209 : :
1210 [ + + ]: 176 : call c_f_pointer(c_loc(eps_gwnlf(1,1)), rvals, shape=[2, bsp%nomega, bsp%nq])
1211 [ - + ]: 88 : NCF_CHECK(nf90_put_var(ncid, vid("gwnlf_mdf"), rvals))
1212 : :
1213 : : contains
1214 : 352 : integer function vid(vname)
1215 : : character(len=*),intent(in) :: vname
1216 [ + - ]: 352 : vid = nctk_idname(ncid, vname)
1217 : : end function vid
1218 : :
1219 : : end subroutine mdfs_ncwrite
1220 : : !!***
1221 : :
1222 : : !!****f* m_exc_spectra/check_kramerskronig
1223 : : !! NAME
1224 : : !! check_kramerskronig
1225 : : !!
1226 : : !! FUNCTION
1227 : : !! check Kramers Kronig
1228 : : !! \int_0^\infty d\omega' frac{\omega'}{\omega'^2 - \omega^2}
1229 : : !! Im \epsilon(\omega') = Re \epsilon(\omega)
1230 : : !!
1231 : : !! INPUTS
1232 : : !! n=Number of frequency points.
1233 : : !! eps(n)=Dielectric function.
1234 : : !! o(n)=Frequency mesh.
1235 : : !!
1236 : : !! OUTPUT
1237 : : !! Only checking.
1238 : : !!
1239 : : !! SOURCE
1240 : :
1241 : 48 : subroutine check_kramerskronig(n,o,eps)
1242 : :
1243 : : !Arguments ------------------------------------
1244 : : !scalars
1245 : : integer,intent(in) :: n
1246 : : !arrays
1247 : : complex(dp),intent(in) :: eps(n)
1248 : : real(dp),intent(in) :: o(n)
1249 : :
1250 : : !Local variables ------------------------------
1251 : : !scalars
1252 : : integer :: ii,ip
1253 : : real(dp) :: omega,omegap,domega,kk,kkrms,eav
1254 : : complex(dp) :: c
1255 : : character(len=500) :: msg
1256 : : !arrays
1257 : 96 : real(dp) :: e1kk(n)
1258 : 96 : complex(dp) :: intg(n)
1259 : : !************************************************************************
1260 : : ! init jmb
1261 [ + + ]: 42720 : e1kk=zero
1262 [ + + ]: 42720 : intg=(zero,zero)
1263 : :
1264 : : ! calculate domega step and verify all
1265 : 48 : domega = (o(n) - o(1)) / (n-1)
1266 : :
1267 [ + + ]: 42672 : do ii=2,n
1268 [ - + ]: 42672 : if (domega-(o(ii)-o(ii-1)) > tol3) then
1269 : 0 : ABI_WARNING("Frequency mesh not linear. Returning")
1270 : 0 : return
1271 : : end if
1272 : : end do
1273 : :
1274 [ - + ]: 48 : if(o(1) > 0.1/Ha_eV) then
1275 : 0 : ABI_WARNING("First frequency is not zero. Returning")
1276 : 0 : return
1277 : : end if
1278 : :
1279 [ + + ]: 48 : if (aimag(eps(n)) > 0.1) then
1280 : : write(msg,'(a,f12.6,3a,f12.6,2a)')&
1281 : 36 : ' Im epsilon for omega= ',o(n)*Ha_eV,'eV',ch10,&
1282 : 36 : ' is not yet zero, epsilon_2= ',aimag(eps(n)),ch10,&
1283 : 72 : ' Kramers Kronig test could give wrong results. '
1284 : 36 : ABI_WARNING(msg)
1285 : : end if
1286 : :
1287 : : ! Fill array for kramers kronig.
1288 [ + + ]: 42720 : do ii=1,n
1289 : 42672 : omega=o(ii)
1290 : 42672 : c = (0.0,0.0)
1291 [ + + ]: 45902664 : do ip=1,n
1292 [ + + ]: 45859992 : if(ip == ii) cycle
1293 : 45817320 : omegap = o(ip)
1294 : 45902664 : c = c + omegap / (omegap**2-omega**2) * aimag(eps(ip))
1295 : : end do
1296 : 42720 : e1kk(ii) = one + two/pi * domega*real(c)
1297 : : end do
1298 : :
1299 : : ! perform kramers kronig with simpson integration
1300 [ + + ]: 42720 : do ii=1,n
1301 : 42672 : omega=o(ii)
1302 [ + + ]: 45902664 : do ip=1,n
1303 [ + + ]: 45859992 : if (ip==ii) cycle
1304 : 45817320 : omegap = o(ip)
1305 : 45902664 : intg(ip) = omegap / (omegap**2 - omega**2) * aimag(eps(ip))
1306 : : end do
1307 : 42672 : c = simpson_cplx(n,domega,intg)
1308 : 42720 : e1kk(ii) = one + two/pi * real(c)
1309 : : end do
1310 : :
1311 : : !verify kramers kronig
1312 : : eav=zero; kk=zero; kkrms=zero
1313 : : do ii=1,n
1314 : : kk = kk + abs(real(eps(ii)) - e1kk(ii))
1315 : : kkrms = kkrms +(real(eps(ii)) - e1kk(ii))*(real(eps(ii)) - e1kk(ii))
1316 : : eav = eav + abs(real(eps(ii)))
1317 : : end do
1318 : :
1319 : 48 : eav = eav/n
1320 : 48 : kk = (kk/n)/eav
1321 : 48 : kkrms = (kkrms/n) / (eav*eav)
1322 : :
1323 : 48 : kk = abs(real(eps(1)) - e1kk(1)) / real(eps(1))
1324 : :
1325 : : ! write data
1326 : 48 : write(msg,'(a,f7.2,a)')" The Kramers-Kronig is verified within ",100*kk,"%"
1327 : 48 : call wrtout(std_out, msg)
1328 : :
1329 : : ! write(std_out,'("# Kramers Kronig calculation of epsilon1")')
1330 : : ! write(std_out,'("# omega epsilon1 epsilon1kk")')
1331 : : ! do ii=1,n
1332 : : ! write(std_out,'(f7.3,2e15.7)') o(ii)*Ha_eV, real(eps(ii)), e1kk(ii)
1333 : : ! end do
1334 : :
1335 : : end subroutine check_kramerskronig
1336 : : !!***
1337 : :
1338 : : !----------------------------------------------------------------------
1339 : :
1340 : : !!****f* m_exc_spectra/check_fsumrule
1341 : : !! NAME
1342 : : !! check_fsumrule
1343 : : !!
1344 : : !! FUNCTION
1345 : : !! check f-sum rule
1346 : : !! \int_0^\infty d\omega \omega Im \epsilon_GG'(q,\omega) =
1347 : : !! = \frac{1}{2} \pi \omega_p^2 \frac{\rho(G-G')}{\rho(0)}
1348 : : !! versor(q+G) \dot versor(q+G')
1349 : : !! for q = G = G' = 0, it reads:
1350 : : !! \int_0^\infty d\omega \omega Im \epsilon_00(q=0,\omega) =
1351 : : !! = \pi \omega_p^2 / 2
1352 : : !! calculate only the second one
1353 : : !! calculate the integral to evaluate an omega_plasma^eff to compare with omega_plasma
1354 : : !!
1355 : : !! INPUTS
1356 : : !! n=Number of frequencies.
1357 : : !! o(n)=Frequency mesh.
1358 : : !! e2(n)=imaginary part of epsilon_00
1359 : : !! omegaplasma=Drude plasma frequency.
1360 : : !!
1361 : : !! OUTPUT
1362 : : !! Only checking.
1363 : : !!
1364 : : !! SOURCE
1365 : :
1366 : 48 : subroutine check_fsumrule(n,o,e2,omegaplasma)
1367 : :
1368 : : !Arguments ------------------------------------
1369 : : !scalars
1370 : : integer,intent(in) :: n
1371 : : real(dp),intent(in) :: omegaplasma
1372 : : !arrays
1373 : : real(dp),intent(in) :: o(n),e2(n)
1374 : :
1375 : : !Local variables ------------------------------
1376 : : !scalars
1377 : : integer :: ii,ip
1378 : : real(dp) :: omegap,domega,integral,omegaplasmaeff,fsumrule
1379 : : character(len=500) :: msg
1380 : : !arrays
1381 : 48 : complex(dp) :: intg(n)
1382 : : !************************************************************************
1383 : :
1384 : : ! calculate domega step and verify
1385 : 48 : domega = (o(n) - o(1)) / (n-1)
1386 : :
1387 [ + + ]: 42672 : do ii=2,n
1388 [ - + ]: 42672 : if (domega-(o(ii)-o(ii-1)) > tol3) then
1389 : 0 : ABI_WARNING("Frequency mesh not linear. Returning")
1390 : 0 : return
1391 : : end if
1392 : : end do
1393 : :
1394 [ - + ]: 48 : if (o(1) > 0.1/Ha_eV) then
1395 : 0 : ABI_WARNING("First frequency is not zero. Returning")
1396 : 0 : return
1397 : : end if
1398 : :
1399 [ + + ]: 48 : if (e2(n) > 0.1) then
1400 : : write(msg,'(a,f12.6,3a,f12.6,2a)')&
1401 : 36 : ' Im epsilon for omega= ',o(n)*Ha_eV,' eV ',ch10,&
1402 : 36 : ' is not yet zero, epsilon_2= ',e2(n),ch10,&
1403 : 72 : ' f-sum rule test could give wrong results.'
1404 : 36 : ABI_WARNING(msg)
1405 : : end if
1406 : :
1407 : : ! integrate to obtain f-sum rule
1408 : 48 : integral=zero
1409 : 48 : do ip=1,n
1410 : : omegap=o(ip)
1411 : 48 : integral = integral + omegap * e2(ip)
1412 : : end do
1413 : 42720 : integral = domega * integral
1414 : :
1415 : : ! integrate with simpson to obtain f-sum rule
1416 [ + + ]: 42720 : do ip = 1, n
1417 : 42672 : omegap = o(ip)
1418 : 42720 : intg(ip) = omegap * e2(ip)
1419 : : end do
1420 : :
1421 : 48 : integral = real(simpson_cplx(n,domega,intg))
1422 [ - + ]: 48 : if (integral < 0) then
1423 : 0 : ABI_ERROR("The integral of the imaginary of dielectric function is negative !!!")
1424 : : else
1425 : 48 : omegaplasmaeff = sqrt(integral*two/pi)
1426 : : end if
1427 : 48 : fsumrule = abs((omegaplasmaeff - omegaplasma)) / omegaplasma
1428 : :
1429 : : ! write data
1430 : : write(msg,'(3(a,f6.2,2a))')&
1431 : 48 : " omega_plasma = ",omegaplasma*Ha_eV, " [eV]",ch10,&
1432 : 48 : " omega_plasma^eff = ",omegaplasmaeff*Ha_eV," [eV]",ch10,&
1433 : 96 : " the f-sum rule is verified within ",fsumrule*100,"%",ch10
1434 : 48 : call wrtout(std_out, msg)
1435 : :
1436 : : end subroutine check_fsumrule
1437 : : !!***
1438 : :
1439 : : END MODULE m_exc_spectra
|
