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
|
