TABLE OF CONTENTS

ABINIT/kramerskronig [ Functions ]

[ Top ] [ Functions ]

NAME

 kramerskronig

FUNCTION

 check or apply the Kramers Kronig relation:
  Re \epsilon(\omega) = 1 + \frac{2}{\pi}
  \int_0^\infty d\omega' frac{\omega'}{\omega'^2 - \omega^2} Im \epsilon(\omega')

COPYRIGHT

  Copyright (C) 2007-2017 ABINIT group (Valerio Olevano, Lucia Reining, Francesco Sottile, MG)
  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 (both real and imaginary part)
   real part can be used to check whether the K-K relation is satisfied or not
  method=method used to perform the integration
   0= naive integration
   1=simpson rule 
  only_check= if /=0 the real part of eps is checked against the imaginary part,
                a final report in written but the array eps is not modified 
              if ==0 the real part of eps is overwritten using the
              results obtained using the Kramers-Kronig relation

OUTPUT

 Inspired to check_kramerskronig of the DP code

PARENTS

      linear_optics_paw

CHILDREN

      simpson_int,wrtout

SOURCE

 42 #if defined HAVE_CONFIG_H
 43 #include "config.h"
 44 #endif
 45 
 46 #include "abi_common.h"
 47 
 48 subroutine kramerskronig(nomega,omega,eps,method,only_check)
 49 
 50  use defs_basis
 51  use m_errors
 52 
 53  use m_numeric_tools,   only : simpson_int
 54 
 55 !This section has been created automatically by the script Abilint (TD).
 56 !Do not modify the following lines by hand.
 57 #undef ABI_FUNC
 58 #define ABI_FUNC 'kramerskronig'
 59  use interfaces_14_hidewrite
 60 !End of the abilint section
 61 
 62  implicit none
 63 
 64 !Arguments ------------------------------------
 65 !scalars
 66  integer,intent(in) :: method,nomega,only_check
 67 !arrays
 68  real(dp),intent(in) :: omega(nomega)
 69  complex,intent(inout) :: eps(nomega)
 70 
 71 !Local variables-------------------------------
 72 !scalars
 73  integer,save :: enough=0
 74  integer :: ii,ip
 75  real(dp) :: acc,domega,eav,kkdif,kkrms,ww,wwp
 76  character(len=500) :: msg
 77 !arrays
 78  real(dp) :: e1kk(nomega),intkk(nomega),kk(nomega)
 79 
 80 ! *************************************************************************
 81 
 82 !Check whether the frequency grid is linear or not 
 83  domega = (omega(nomega) - omega(1)) / (nomega-1)
 84  do ii=2,nomega
 85    if (ABS(domega-(omega(ii)-omega(ii-1))) > 0.001) then 
 86      if (only_check/=1) then
 87        msg="check cannot be performed since the frequency step is not constant"
 88        MSG_WARNING(msg)
 89        RETURN
 90      else 
 91        msg=' Cannot perform integration since frequency step is not constant'
 92        MSG_ERROR(msg)
 93      end if 
 94    end if
 95  end do
 96 
 97 !Check whether omega(1) is small or not
 98  if (omega(1) > 0.1/Ha_eV) then 
 99    if (only_check/=1) then
100      msg=' Check cannot be performed since first frequency on the grid > 0.1 eV'
101      MSG_WARNING(msg)
102      RETURN
103    else 
104      msg=' Cannot perform integration since first frequency on the grid > 0.1 eV'
105      MSG_ERROR(msg)
106    end if
107  end if 
108 
109 !If eps(nomega) is not 0 warn
110  if (AIMAG(eps(nomega)) > 0.1 .and. enough<50) then
111    enough=enough+1
112    write(msg,'(a,f8.4,3a,f8.2,2a)')&
113 &   'Im epsilon for omega = ',omega(nomega)*Ha_eV,' eV',ch10,&
114 &   'is not yet zero, epsilon_2 = ',AIMAG(eps(nomega)),ch10,&
115 &   'Kramers Kronig could give wrong results'
116    MSG_WARNING(msg)
117    if (enough==50) then 
118      write(msg,'(3a)')' sufficient number of WARNINGS-',ch10,' stop writing '
119      call wrtout(std_out,msg,'COLL')
120    end if 
121  end if
122  
123  
124 !Perform Kramers-Kronig using naive integration
125  select case (method)
126  case (0)
127 
128    do ii=1,nomega
129      ww = omega(ii)
130      acc = 0.0_dp
131      do ip=1,nomega
132        if (ip == ii) CYCLE
133        wwp = omega(ip)
134        acc = acc + wwp/(wwp**2-ww**2) *AIMAG(eps(ip))
135      end do
136      e1kk(ii) = one + two/pi*domega* acc
137    end do
138    
139 !    Perform Kramers-Kronig using Simpson integration    
140 !    Simpson O(1/N^4), from NumRec in C p 134  NumRec in Fortran p 128
141  case (1)
142 
143    kk=zero
144 
145    do ii=1,nomega
146      ww=omega(ii)
147      do ip=1,nomega
148        if (ip == ii) CYCLE
149        wwp = omega(ip)
150        kk(ip) = wwp/(wwp**2-ww**2) *AIMAG(eps(ip))
151      end do
152      call simpson_int(nomega,domega,kk,intkk)
153      e1kk(ii) = one + two/pi * intkk(nomega)
154    end do
155 
156  case default
157    write(msg,'(a,i0)')' Wrong value for method ',method
158    MSG_BUG(msg)
159  end select
160 
161 !at this point real part is in e1kk, need to put it into eps
162  do ii=1,nomega
163    eps(ii)=CMPLX(e1kk(ii),AIMAG(eps(ii)))
164  end do 
165 
166 !Verify Kramers-Kronig
167  eav   = zero
168  kkdif = zero
169  kkrms = zero
170 
171  do ii=1,nomega
172    kkdif = kkdif + ABS(REAL(eps(ii)) - e1kk(ii))
173    kkrms = kkrms + (REAL(eps(ii)) - e1kk(ii))*(REAL(eps(ii)) - e1kk(ii))
174    eav = eav + ABS(REAL(eps(ii)))
175  end do
176 
177  eav = eav/nomega
178  kkdif = (kkdif/nomega) / eav
179  kkrms = (kkrms/nomega) / (eav*eav)
180 
181  kk = ABS(REAL(eps(1)) - e1kk(1)) / REAL(eps(1))
182 
183 !Write data
184  write(msg,'(a,f7.2,a)')' Kramers-Kronig transform is verified within ',MAXVAL(kk)*100,"%"
185  call wrtout(std_out,msg,'COLL')
186 
187 end subroutine kramerskronig