TABLE OF CONTENTS
ABINIT/fsumrule [ Functions ]
NAME
fsumrule
FUNCTION
This routine performs a check of the f-sum rule: \int_0^\infty d\omega \omega Im \epsilon_GGp(q,\omega) = = \frac{1}{2} \pi \omega_p^2 \frac{\rho(G-Gp)}{\rho(0)} versor(q+G) \dot versor(q+Gp) for q = G = Gp = 0, it reads: \int_0^\infty d\omega \omega Im \epsilon_00(q=0,\omega) = \pi \omega_p^2 / 2 check only the second one: calculate the integral to evaluate an omega_plasma^eff to compare with omega_plasma
COPYRIGHT
Copyright (C) 2007-2018 ABINIT group (MG,VO,LR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
nomega=number of real frequencies omega(nomega)= real frequencies eps(nomega)= function on the frequency grid in the present implementation must be the imaginary part of epsilon_00 (q=0) method=method used to perform the integration 0= naive integration 1=simpson rule
OUTPUT
Only check.
NOTES
Inspired to a similar routine of the dp code.
PARENTS
CHILDREN
simspon_int,wrtout
SOURCE
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 subroutine fsumrule(nomega,omega,eps,omegaplasma,method) 54 55 use m_profiling_abi 56 57 use defs_basis 58 use m_errors 59 60 !This section has been created automatically by the script Abilint (TD). 61 !Do not modify the following lines by hand. 62 #undef ABI_FUNC 63 #define ABI_FUNC 'fsumrule' 64 use interfaces_14_hidewrite 65 !End of the abilint section 66 67 implicit none 68 69 !Arguments ------------------------------------ 70 !scalars 71 integer,intent(in) :: method,nomega 72 real(dp),intent(in) :: omegaplasma 73 !arrays 74 real(dp),intent(in) :: eps(nomega),omega(nomega) 75 76 !Local variables------------------------------- 77 !scalars 78 integer :: ii,ip 79 real(dp) :: acc,check,domega,omegaplasmaeff,ww,wwp 80 character(len=500) :: msg 81 !arrays 82 real(dp) :: fs(nomega),fsint(nomega) 83 84 ! ************************************************************************* 85 86 !Check whether the frequency grid is linear or not 87 domega = (omega(nomega) - omega(1))/(nomega-1) 88 do ii =2,nomega 89 if (ABS(domega-(omega(ii)-omega(ii-1))) > 0.001) then 90 msg = 'Check cannot be performed since frequency step is not constant' 91 MSG_WARNING(msg) 92 RETURN 93 end if 94 end do 95 96 !Check whether omega(1) is small or not 97 if (omega(1) > 0.1/Ha_eV) then 98 msg='Check cannot be performed since first frequency on the grid > 0.1 eV' 99 MSG_WARNING(msg) 100 RETURN 101 end if 102 103 !If eps(nomega) is not 0 warn 104 if (eps(nomega) > 0.1) then 105 write(msg,'(a,f8.4,3a,f8.2,2a)')& 106 & ' Im epsilon for omega = ',omega(nomega)*Ha_eV,' eV',ch10,& 107 & ' is not yet zero, epsilon_2 = ',eps(nomega),ch10,& 108 & ' f-sum rule test could give wrong results' 109 MSG_WARNING(msg) 110 end if 111 112 !Check f-sum rule using naive integration 113 SELECT CASE (method) 114 115 CASE (0) 116 acc=zero 117 do ip=1,nomega 118 wwp=omega(ip) 119 acc=acc+wwp*e2(ip) 120 end do 121 acc=domega*acc 122 123 omegaplasmaeff=SQRT(acc*2/pi) 124 125 ! Perform Kramers-Kronig using Simpson integration 126 ! Simpson O(1/N^4), from NumRec in C p 134 NumRec in Fortran p 128 127 CASE (1) 128 do ip=1,nomega 129 wwp=omega(ip) 130 fs(ip)=wwp*e2(ip) 131 end do 132 133 call simspon_int(nomega,domega,fsint) 134 omegaplasmaeff=SQRT(fsint(nomega)*2/pi) 135 136 CASE DEFAULT 137 write(msg,'(a,i3)')' wrong value for method ',method 138 MSG_BUG(msg) 139 END SELECT 140 141 check=ABS((omegaplasmaeff-omegaplasma))/omegaplasma 142 143 !Write data 144 write(msg,'(2a,f6.2,3a,f6.2,3a,f6.2,2a)'),ch10,& 145 & '( omega_plasma = ',omegaplasma*Ha_eV, ' [eV] )',ch10,& 146 & '( omega_plasma^eff = 'omegaplasmaeff*Ha_eV,' [eV] )',ch10,& 147 & '( the f-sum rule is verified within ',check*100,'% )',ch10 148 call wrtout(std_out,msg,'COLL') 149 150 end subroutine fsumrule