TABLE OF CONTENTS


ABINIT/fsumrule [ Functions ]

[ Top ] [ 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