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