TABLE OF CONTENTS
ABINIT/msig [ Functions ]
NAME
msig
FUNCTION
This program computes the elements of the optical frequency dependent conductivity tensor and the conductivity along the three principal axes from the Kubo-Greenwood formula for PAW formalism
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (SMazevet,VR) 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/Infos/contributors .
INPUTS
fcti(npti)= conductivity, as calculated in conducti npti= number of points to calculate conductivity xi(npti)= energies where the conductivity is calculated
OUTPUT
no output, only files
NOTES
this program calculates the imaginary part of the conductivity (principal value)
+derived optical properties.
the calculation is performed on the same grid as the initial input
to calculate the principal value, a trapezoidale integration +taylor expansion to
third order is used (W.J. Thomson computer in physics vol 12 p94 1998)
two input files are needed inppv.dat (parameters) and sigma.dat (energy,sigma_1)
two output files ppsigma.dat (energy,sigma_1,sigma_2,epsilon_1,epsilon_2)
abs.dat (energy,nomega,komega,romega,absomega)
march 2002 s.mazevet
PARENTS
conducti_nc,conducti_paw
CHILDREN
intrpl
SOURCE
45 #if defined HAVE_CONFIG_H 46 #include "config.h" 47 #endif 48 49 #include "abi_common.h" 50 51 subroutine msig(fcti,npti,xi,filnam_out_sig) 52 53 use defs_basis 54 use m_profiling_abi 55 use m_errors 56 57 use m_io_tools, only : open_file 58 use m_splines, only : intrpl 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 'msig' 64 !End of the abilint section 65 66 implicit none 67 68 !Arguments ----------------------------------- 69 !scalars 70 integer,intent(in) :: npti 71 !arrays 72 real(dp),intent(in) :: fcti(npti),xi(npti) 73 character(len=fnlen),intent(in) :: filnam_out_sig 74 75 !Local variables------------------------------- 76 !scalars 77 integer,parameter :: npt=10000 78 integer :: ii,ip,npt1,npt2,eps_unt,abs_unt 79 real(dp),parameter :: del=0.001_dp,ohmtosec=9.d11 80 real(dp) :: dx,dx1,dx2,eps1,eps2,idel,komega,pole,refl,sigma2,xsum 81 character(len=500) :: msg 82 !arrays 83 real(dp),allocatable :: abso(:),fct(:),fct1(:),fct2(:),fct3(:),fct4(:),fct5(:) 84 real(dp),allocatable :: fctii(:),fp(:),fpp(:),fppp(:),nomega(:),ppsig(:) 85 real(dp),allocatable :: x1(:),x2(:) 86 87 ! ********************************************************************************* 88 !BEGIN EXECUTABLE SECTION 89 90 write(std_out,'(2a)')ch10,'Calculate the principal value and related optical properties' 91 write(std_out,'(a)')'following W.J. Thomson computer in physics vol 12 p94 1998 for ' 92 write(std_out,'(a)')'the principal value. S. Mazevet' 93 write(std_out,'(a)')'OPTIONS' 94 write(std_out,'(a)')'use default number of integration pts: npt=10000' 95 write(std_out,'(a)')'Use default value for delta interval: del=1e-3' 96 97 if (open_file(trim(filnam_out_sig)//'_eps',msg, newunit=eps_unt,status='replace',action="write") /= 0) then 98 MSG_ERROR(msg) 99 end if 100 write(eps_unt,'(a)')'#energy (eV),sigma_1(Ohm-1cm-1),sigma_2(Ohm-1cm-1),epsilon_1,epsilon_2' 101 102 if (open_file(trim(filnam_out_sig)//'_abs', msg, newunit=abs_unt, status='replace',action="write") /= 0) then 103 MSG_ERROR(msg) 104 end if 105 write(abs_unt,'(a)')'#energy(eV),nomega,komega,refl.,abso.(cm-1)' 106 107 ABI_ALLOCATE(fct,(npt)) 108 ABI_ALLOCATE(fct2,(npt)) 109 ABI_ALLOCATE(fct3,(npt)) 110 ABI_ALLOCATE(fct4,(npt)) 111 ABI_ALLOCATE(fct5,(npt)) 112 ABI_ALLOCATE(fp,(npt)) 113 ABI_ALLOCATE(fpp,(npt)) 114 ABI_ALLOCATE(fppp,(npt)) 115 ABI_ALLOCATE(x1,(npt)) 116 ABI_ALLOCATE(x2,(npt)) 117 ABI_ALLOCATE(fct1,(npt)) 118 ABI_ALLOCATE(ppsig,(npt)) 119 ABI_ALLOCATE(fctii,(npt)) 120 ABI_ALLOCATE(abso,(npt)) 121 ABI_ALLOCATE(nomega,(npt)) 122 123 if (npti > npt) then 124 write (std_out,*) 'msig: input npti is too large for hard coded npt array size = ', npt 125 MSG_ERROR("Aborting now") 126 end if 127 128 !loop on the initial energy grid 129 do ip=1,npti 130 131 ! adjust the interval before and after the pole to reflect range/npt interval 132 xsum=zero 133 dx=(xi(npti)-xi(1))/dble(npt-1) 134 pole=xi(ip) 135 npt1=int((pole-del)/dx) 136 dx1=zero 137 if(npt1/=1) dx1=(pole-del)/(npt1-1) 138 npt2=int((xi(npti)-pole-del)/dx) 139 dx2=(xi(npti)-pole-del)/(npt2-1) 140 141 ! for the moment skip the pp calculation when the pole if too close to the end of the range 142 if(npt1<=1.or.npt2<=1) then 143 144 xsum=zero 145 ppsig(ip)=zero 146 147 else 148 149 ! define the fct for which the pp calculation is needed using xi^2-pole^2 factorization 150 fctii(1:npti) = zero 151 fctii(1:npti)=fcti(1:npti)*pole/(xi(1:npti)+pole) 152 153 ! define the grid on each side of the pole x1 before x2 after 154 do ii=1,npt1 155 x1(ii)=dx1*dble(ii-1) 156 end do 157 do ii=1,npt2 158 x2(ii)=pole+del+dx2*dble(ii-1) 159 end do 160 161 ! interpolate the initial fct fii on the new grids x1 and x2 (cubic spline) 162 ! write(std_out,*) npti,npt1 163 164 ! MJV 6/12/2008: 165 ! for each use of fctii should ensure that npt1 npt2 etc... are less than 166 ! npt=len(fctii) 167 call intrpl(npti,xi,fctii,npt1,x1,fct4,fct1,fct5,1) 168 call intrpl(npti,xi,fctii,npt2,x2,fct3,fct2,fct5,1) 169 170 ! calculate the two integrals from 0-->pole-lamda and pole+lamda--> end range 171 ! trapezoidal integration 172 do ii=1,npt1 173 fct1(ii)=fct4(ii)/(x1(ii)-pole) 174 end do 175 do ii=1,npt2 176 fct2(ii)=fct3(ii)/(x2(ii)-pole) 177 end do 178 179 do ii=2,npt1-1 180 xsum=xsum+fct1(ii)*dx1 181 end do 182 do ii=2,npt2-1 183 xsum=xsum+fct2(ii)*dx2 184 end do 185 xsum=xsum+half*(fct1(1)+fct1(npt1))*dx1+half*(fct2(1)+fct2(npt2))*dx2 186 187 ! calculate the first and third derivative at the pole and add the taylor expansion 188 call intrpl(npti,xi,fctii,npti,xi,fct3,fct4,fct5,1) 189 call intrpl(npti,xi,fct4,1,(/pole/),fp,fpp,fppp,1) 190 191 idel=two*fp(1)*(del)+fppp(1)*(del**3)/nine 192 xsum=xsum+idel 193 194 end if 195 196 ! calculate the derivated optical quantities and output the value 197 sigma2=(-two/pi)*xsum 198 eps1=one-(four_pi*sigma2/(pole)) 199 eps2=four*fcti(ip)*pi/(pole) 200 201 ! A special treatment of the case where eps2 is very small compared to eps1 is needed 202 if(eps2**2 > eps1**2 * tol12)then 203 nomega(ip)=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2))) 204 komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2))) 205 abso(ip)=four_pi*fcti(ip)*ohmtosec*Ohmcm/nomega(ip)/(Sp_Lt_SI*100._dp) 206 else if(eps1>zero)then 207 nomega(ip)=sqrt(half*(eps1 + sqrt(eps1**2 + eps2**2))) 208 komega=half*abs(eps2/sqrt(eps1)) 209 abso(ip)=four_pi*fcti(ip)*ohmtosec*Ohmcm/nomega(ip)/(Sp_Lt_SI*100._dp) 210 else if(eps1<zero)then 211 nomega(ip)=half*abs(eps2/sqrt(-eps1)) 212 komega=sqrt(half*(-eps1 + sqrt(eps1**2 + eps2**2))) 213 abso(ip)=two*sqrt(-eps1)*pole*ohmtosec*Ohmcm/(Sp_Lt_SI*100._dp) 214 end if 215 216 refl=((one-nomega(ip))**2 + komega**2)/ & 217 & ((one+nomega(ip))**2 + komega**2) 218 219 write(eps_unt,'(5e18.10)') Ha_eV*pole,fcti(ip)*Ohmcm,sigma2*Ohmcm,eps1,eps2 220 write(abs_unt,'(5e18.10)') Ha_eV*pole,nomega(ip),komega,refl,abso(ip) 221 222 end do 223 224 close(eps_unt) 225 close(abs_unt) 226 227 ABI_DEALLOCATE(fct) 228 ABI_DEALLOCATE(x1) 229 ABI_DEALLOCATE(x2) 230 ABI_DEALLOCATE(fct2) 231 ABI_DEALLOCATE(fct3) 232 ABI_DEALLOCATE(fct4) 233 ABI_DEALLOCATE(fct5) 234 ABI_DEALLOCATE(fp) 235 ABI_DEALLOCATE(fpp) 236 ABI_DEALLOCATE(fppp) 237 ABI_DEALLOCATE(fct1) 238 ABI_DEALLOCATE(ppsig) 239 ABI_DEALLOCATE(fctii) 240 ABI_DEALLOCATE(abso) 241 ABI_DEALLOCATE(nomega) 242 243 end subroutine msig