TABLE OF CONTENTS
ABINIT/check_fsumrule [ 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
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
n=Number of frequencies. o(n)=Frequency mesh. e2(n)=imaginary part of epsilon_00 omegaplasma=Drude plasma frequency.
OUTPUT
Only checking.
PARENTS
m_exc_spectra
CHILDREN
wrtout
SOURCE
194 #if defined HAVE_CONFIG_H 195 #include "config.h" 196 #endif 197 198 #include "abi_common.h" 199 200 subroutine check_fsumrule(n,o,e2,omegaplasma) 201 202 use defs_basis 203 use m_errors 204 use m_profiling_abi 205 206 use m_numeric_tools, only : simpson_cplx 207 208 !This section has been created automatically by the script Abilint (TD). 209 !Do not modify the following lines by hand. 210 #undef ABI_FUNC 211 #define ABI_FUNC 'check_fsumrule' 212 use interfaces_14_hidewrite 213 !End of the abilint section 214 215 implicit none 216 217 !Arguments ------------------------------------ 218 !scalars 219 integer,intent(in) :: n 220 real(dp),intent(in) :: omegaplasma 221 !arrays 222 real(dp),intent(in) :: o(n),e2(n) 223 224 !Local variables ------------------------------ 225 !scalars 226 integer :: ii,ip 227 real(dp) :: omegap,domega,integral,omegaplasmaeff,fsumrule 228 character(len=500) :: msg 229 !arrays 230 complex(dpc) :: intg(n) 231 232 !************************************************************************ 233 234 ! calculate domega step and verify 235 domega = (o(n) - o(1)) / (n-1) 236 237 do ii=2,n 238 if (domega-(o(ii)-o(ii-1)) > tol3) then 239 MSG_WARNING("Frequency mesh not linear. Returning") 240 return 241 end if 242 end do 243 244 if (o(1) > 0.1/Ha_eV) then 245 MSG_WARNING("First frequency is not zero. Returning") 246 return 247 end if 248 249 if (e2(n) > 0.1) then 250 write(msg,'(a,f12.6,3a,f12.6,2a)')& 251 & ' Im epsilon for omega= ',o(n)*Ha_eV,' eV ',ch10,& 252 & ' is not yet zero, epsilon_2= ',e2(n),ch10,& 253 & ' f-sum rule test could give wrong results.' 254 MSG_WARNING(msg) 255 end if 256 257 ! integrate to obtain f-sum rule 258 integral=zero 259 do ip=1,n 260 omegap=o(ip) 261 integral = integral + omegap * e2(ip) 262 end do 263 integral = domega * integral 264 265 !integrate with simpson to obtain f-sum rule 266 do ip = 1, n 267 omegap = o(ip) 268 intg(ip) = omegap * e2(ip) 269 end do 270 271 integral = real(simpson_cplx(n,domega,intg)) 272 if(integral < 0) then 273 MSG_ERROR("The integral of the imaginary of dielectric function is negative !!!") 274 else 275 omegaplasmaeff = sqrt(integral*two/pi) 276 end if 277 278 fsumrule = abs((omegaplasmaeff - omegaplasma)) / omegaplasma 279 280 ! write data 281 write(msg,'(3(a,f6.2,2a))')& 282 & " omega_plasma = ",omegaplasma*Ha_eV, " [eV]",ch10,& 283 & " omega_plasma^eff = ",omegaplasmaeff*Ha_eV," [eV]",ch10,& 284 & " the f-sum rule is verified within ",fsumrule*100,"%",ch10 285 call wrtout(std_out,msg,"COLL") 286 287 end subroutine check_fsumrule
ABINIT/check_kramerskronig [ 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)
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
n=Number of frequency points. eps(n)=Dielectric function. o(n)=Frequency mesh.
OUTPUT
Only checking.
PARENTS
m_exc_spectra
CHILDREN
wrtout
SOURCE
34 #if defined HAVE_CONFIG_H 35 #include "config.h" 36 #endif 37 38 #include "abi_common.h" 39 40 subroutine check_kramerskronig(n,o,eps) 41 42 use defs_basis 43 use m_errors 44 use m_profiling_abi 45 46 use m_numeric_tools, only : simpson_cplx 47 48 !This section has been created automatically by the script Abilint (TD). 49 !Do not modify the following lines by hand. 50 #undef ABI_FUNC 51 #define ABI_FUNC 'check_kramerskronig' 52 use interfaces_14_hidewrite 53 !End of the abilint section 54 55 implicit none 56 57 !Arguments ------------------------------------ 58 !scalars 59 integer,intent(in) :: n 60 !arrays 61 complex(dpc),intent(in) :: eps(n) 62 real(dp),intent(in) :: o(n) 63 64 !Local variables ------------------------------ 65 !scalars 66 integer :: ii,ip 67 real(dp) :: omega,omegap,domega,kk,kkrms,eav 68 complex(dpc) c 69 character(len=500) :: msg 70 !arrays 71 real(dp) :: e1kk(n) 72 complex(dpc) :: intg(n) 73 74 !************************************************************************ 75 ! init jmb 76 e1kk=zero 77 intg=(zero,zero) 78 79 ! calculate domega step and verify all 80 domega = (o(n) - o(1)) / (n-1) 81 82 do ii=2,n 83 if (domega-(o(ii)-o(ii-1)) > tol3) then 84 MSG_WARNING("Frequency mesh not linear. Returning") 85 return 86 end if 87 end do 88 89 if(o(1) > 0.1/Ha_eV) then 90 MSG_WARNING("First frequency is not zero. Returning") 91 return 92 end if 93 94 if (aimag(eps(n)) > 0.1) then 95 write(msg,'(a,f12.6,3a,f12.6,2a)')& 96 & ' Im epsilon for omega= ',o(n)*Ha_eV,'eV',ch10,& 97 & ' is not yet zero, epsilon_2= ',aimag(eps(n)),ch10,& 98 & ' Kramers Kronig test could give wrong results. ' 99 MSG_WARNING(msg) 100 end if 101 102 ! Fill array for kramers kronig. 103 do ii=1,n 104 omega=o(ii) 105 c = (0.0,0.0) 106 do ip=1,n 107 if(ip == ii) cycle 108 omegap = o(ip) 109 c = c + omegap / (omegap**2-omega**2) * aimag(eps(ip)) 110 end do 111 e1kk(ii) = one + two/pi * domega*real(c) 112 end do 113 114 !perform kramers kronig with simpson integration 115 do ii=1,n 116 omega=o(ii) 117 do ip=1,n 118 if (ip==ii) cycle 119 omegap = o(ip) 120 intg(ip) = omegap / (omegap**2 - omega**2) * aimag(eps(ip)) 121 end do 122 c = simpson_cplx(n,domega,intg) 123 e1kk(ii) = one + two/pi * real(c) 124 end do 125 126 !verify kramers kronig 127 eav=zero; kk=zero; kkrms=zero 128 do ii=1,n 129 kk = kk + abs(real(eps(ii)) - e1kk(ii)) 130 kkrms = kkrms +(real(eps(ii)) - e1kk(ii))*(real(eps(ii)) - e1kk(ii)) 131 eav = eav + abs(real(eps(ii))) 132 end do 133 134 eav = eav/n 135 kk = (kk/n)/eav 136 kkrms = (kkrms/n) / (eav*eav) 137 138 kk = abs(real(eps(1)) - e1kk(1)) / real(eps(1)) 139 140 ! write data 141 write(msg,'(a,f7.2,a)')" The Kramers-Kronig is verified within ",100*kk,"%" 142 call wrtout(std_out,msg,"COLL") 143 144 ! write(std_out,'("# Kramers Kronig calculation of epsilon1")') 145 ! write(std_out,'("# omega epsilon1 epsilon1kk")') 146 ! do ii=1,n 147 ! write(std_out,'(f7.3,2e15.7)') o(ii)*Ha_eV, real(eps(ii)), e1kk(ii) 148 ! end do 149 150 end subroutine check_kramerskronig