1!=============================================================================== 2! 3! Utilities: 4! 5! (1) epsinvomega Originally By gsm Last Modified 11/17/2010 (gsm) 6! 7! Serial code that plots frequency dependence of Plasmon-Pole, Retarded, 8! Advanced epsilon inverse for a given q, G, G` vectors. Input parameters 9! are read from file epsinvomega.inp in the working directory. 10! 11! Note that epsInvPP constructed with full-frequency epsmat will be 12! different from epsInvPP constructed with static epsmat because of 13! finite broadening used in full-frequency epsmat. 14! 15! FIXME: use WFN_RHO_VXC for reading RHO 16! 17! FIXME: [2013-01-15 gsm] compact form of GPP model doesn`t work when 18! wtilde2 < 0 which is often the case for G .ne. G`; replace with the 19! explicit expression from Hybertsen & Louie (see commented out code 20! in the GPP section of epsinvomega.f90) 21! 22!=============================================================================== 23 24#include "f_defs.h" 25 26program epsinvomega 27 28 use global_m 29 use epsread_hdf5_m 30#ifdef HDF5 31 use hdf5 32#endif 33 34 implicit none 35 36 integer :: i,j,k,iq,mq,ig,igp,iunit 37 complex(DPC) :: omega 38 39 character*256, parameter :: fninp = "epsinvomega.inp" 40 character*256 :: fneps,fnrho,fngpp,fnffr,fnffa 41 character :: infile*80 42 real(DP) :: q(3) 43 integer :: g(3),gp(3),gmgp(3) 44 45 integer :: freq_dep,nFreq_tmp,nFreq,nfreq_imag,ii,nq,ng,nmtx,kgrid(3) 46 real(DP) :: dDeltaFreq,dBrdning,ecuts,delta 47 real(DP), allocatable :: qpt(:,:) 48 real(DP), allocatable :: ekin(:) 49 real(DP), allocatable :: dFreqGrid(:) 50 integer, allocatable :: isrtx(:) 51 integer, allocatable :: isorti(:) 52 integer, allocatable :: nmtx_of_q(:) 53 SCALAR, allocatable :: eps(:) 54 complex(DPC), allocatable :: dFreqBrd(:) 55 complex(DPC), allocatable :: epsPP(:) 56 complex(DPC), allocatable :: epsR(:) 57 complex(DPC), allocatable :: epsA(:) 58 complex(DPC), allocatable :: epsRTemp(:,:) 59 complex(DPC), allocatable :: epsATemp(:,:) 60 61 real(DP) :: bvec(3,3),bdot(3,3),blat,celvol,reccelvol,ebind 62 integer :: nvecs,nspin,qgrid(3) 63 64 integer nproc_para,num_gvec, error 65 complex(DPC) :: epsStatic 66 real(DP) :: rho0 67 SCALAR :: rhogmgp 68 integer, allocatable :: gvec(:,:) 69 SCALAR, allocatable :: xcdum(:,:) 70 71 real(DP) :: wp2,qg(3),qgp(3),qgqg,qgqgp,lambda,phi 72 SCALAR :: Omega2,wtilde2,epsggp,I_epsggp,eps_static,eps_dynamic 73 complex(DPC) :: wtilde2_temp 74! real(DP) :: epsPPRe,epsPPIm 75! complex(DPC) :: wtilde,ampl 76 77!----------------------------- 78! read input file 79 80 write(6,'(/,1x,"reading",1x,a,1x,"file",/)')trim(fninp) 81 82 call open_file(55,file=trim(fninp),form='formatted',status='old') 83 84 read(55,'(a)') fneps 85 read(55,'(a)') fnrho 86 read(55,*) (q(i),i=1,3) 87 read(55,*) (g(i),i=1,3) 88 read(55,*) (gp(i),i=1,3) 89 read(55,*) nFreq 90 read(55,*) dDeltaFreq 91 read(55,*) dBrdning 92 read(55,*) ebind 93 read(55,'(a)') fngpp 94 read(55,'(a)') fnffr 95 read(55,'(a)') fnffa 96 97 call close_file(55) 98 99 write(6,'(2a)') " eps file = ", trim(fneps) 100 write(6,'(2a)') " rho file = ", trim(fnrho) 101 write(6,'(a, 3f15.12)') " q vector = ", q(1:3) 102 write(6,'(a, 3i4)') " G vector = ", g(1:3) 103 write(6,'(a, 3i4)') " G' vector = ", gp(1:3) 104 write(6,'(a, i6)') " nFreq = ", nFreq 105 write(6,'(a, f7.3)') " dDeltaFreq = ", dDeltaFreq 106 write(6,'(a, f7.3)') " dBrdning = ", dBrdning 107 write(6,'(a, f7.3)') " ebind = ", ebind 108 write(6,'(2a)') " GPP file = ", trim(fngpp) 109 write(6,'(2a)') " FFR file = ", trim(fnffr) 110 write(6,'(2a,/)') " FFA file = ", trim(fnffa) 111 112 gmgp(:)=g(:)-gp(:) 113 114!----------------------------- 115! read eps file 116 117 write(6,'(1x,"reading",1x,a,1x,"file",/)')trim(fneps) 118 119#ifdef HDF5 120 121 call h5open_f(error) 122 123 infile = trim(fneps) 124 125 call read_eps_grid_sizes_hdf5(ng, nq, ecuts, nfreq_tmp, nfreq_imag, nmtx, qgrid, freq_dep, infile) 126 if (freq_dep==2 .and. nFreq_tmp/=nFreq) nFreq = nFreq_tmp 127 128 SAFE_ALLOCATE(dFreqGrid,(nFreq)) 129 SAFE_ALLOCATE(dFreqBrd,(nFreq)) 130 SAFE_ALLOCATE(qpt, (3,nq)) 131 SAFE_ALLOCATE(nmtx_of_q, (nq)) 132 SAFE_ALLOCATE(gvec, (3,ng)) 133 134 call read_eps_qgrid_hdf5(nq,qpt,nmtx_of_q,infile) 135 if (freq_dep .eq. 2) then 136 call read_eps_freqgrid_hdf5(nFreq,dFreqGrid,dFreqBrd,infile) 137 else 138 do i=1,nFreq 139 dFreqGrid(i)=dDeltaFreq*dble(i-1) 140 dFreqBrd(i)=dBrdning 141 enddo 142 endif 143 144 call read_eps_old_gvecs_hdf5(ng,gvec,infile) 145 146#else 147 148 iunit=12 149 call open_file(unit=iunit,file=trim(fneps),form='unformatted',status='old') 150 151 read(iunit) 152 read(iunit) freq_dep,ii 153 if (freq_dep.eq.2) then 154 nFreq_tmp=ii 155 endif 156 if (freq_dep==2 .and. nFreq_tmp/=nFreq) nFreq = nFreq_tmp 157 read(iunit) (kgrid(i),i=1,3) 158 SAFE_ALLOCATE(dFreqGrid,(nFreq)) 159 SAFE_ALLOCATE(dFreqBrd,(nFreq)) 160 if (freq_dep.eq.2) then 161 read(iunit) (dFreqGrid(i),i=1,nFreq),(dFreqBrd(i),i=1,nFreq) 162 if (nFreq.gt.1) dDeltaFreq=dFreqGrid(2)-dFreqGrid(1) 163 dBrdning=IMAG(dFreqBrd(1)) 164 else 165 do i=1,nFreq 166 dFreqGrid(i)=dDeltaFreq*dble(i-1) 167 dFreqBrd(i)=dBrdning 168 enddo 169 read(iunit) 170 endif 171 read(iunit) 172 read(iunit) 173 read(iunit) ecuts 174 read(iunit) nq 175 read(iunit) ng 176 177 rewind(iunit) 178 179 SAFE_ALLOCATE(qpt, (3,nq)) 180 SAFE_ALLOCATE(gvec, (3,ng)) 181 182 read(iunit) 183 read(iunit) 184 read(iunit) 185 read(iunit) 186 read(iunit) 187 read(iunit) 188 read(iunit) 189 read(iunit) nq,((qpt(j,i),j=1,3),i=1,nq) 190 read(iunit) ng,((gvec(j,i),j=1,3),i=1,ng) 191 192#endif 193 194 ig=0 195 igp=0 196 do i=1,ng 197 if (gvec(1,i).eq.g(1).and.gvec(2,i).eq.g(2).and. & 198 gvec(3,i).eq.g(3)) ig=i 199 if (gvec(1,i).eq.gp(1).and.gvec(2,i).eq.gp(2).and. & 200 gvec(3,i).eq.gp(3)) igp=i 201 enddo 202 if (ig.eq.0) then 203 call die("cannot find G vector in file " // trim(fneps)) 204 endif 205 if (igp.eq.0) then 206 call die("cannot find G' vector in file " // trim(fneps)) 207 endif 208 209 mq=-1 210 do iq=1,nq 211 if (abs(q(1)-qpt(1,iq)).lt.TOL_Zero.and. & 212 abs(q(2)-qpt(2,iq)).lt.TOL_Zero.and. & 213 abs(q(3)-qpt(3,iq)).lt.TOL_Zero) then 214 mq=iq-1 215 exit 216 endif 217 enddo 218 if (mq.eq.-1) then 219 call die("cannot find q vector in file " // trim(fneps)) 220 endif 221 222#ifdef HDF5 223 224 SAFE_ALLOCATE(isrtx, (ng)) 225 SAFE_ALLOCATE(isorti, (ng)) 226 SAFE_ALLOCATE(ekin, (ng)) 227 call read_eps_gvecsofq_hdf5(ng,isrtx,isorti,ekin,mq+1,infile) 228 nmtx = nmtx_of_q(mq+1) 229 SAFE_DEALLOCATE(ekin) 230 SAFE_DEALLOCATE(nmtx_of_q) 231 232#else 233 234 do iq=1,mq 235 read(iunit) ng,nmtx 236 read(iunit) 237 read(iunit) 238 if (freq_dep.eq.0) then 239 do j=1,nmtx 240 read(iunit) 241 enddo 242 endif 243 if (freq_dep.eq.2) then 244 do j=1,nmtx 245 do i=1,nmtx 246 read(iunit) 247#ifdef CPLX 248 read(iunit) 249#endif 250 enddo 251 enddo 252 endif 253 enddo 254 255 SAFE_ALLOCATE(isrtx, (ng)) 256 SAFE_ALLOCATE(isorti, (ng)) 257 258 read(iunit) ng,nmtx,(isrtx(i),isorti(i),i=1,ng) 259 read(iunit) 260 read(iunit) 261 262#endif 263 264 ig=isorti(ig) 265 igp=isorti(igp) 266 if (ig.eq.0) then 267 call die("cannot find G vector in file " // trim(fneps)) 268 endif 269 if (igp.eq.0) then 270 call die("cannot find G' vector in file " // trim(fneps)) 271 endif 272 273 SAFE_ALLOCATE(epsPP, (nFreq)) 274 SAFE_ALLOCATE(epsR, (nFreq)) 275 SAFE_ALLOCATE(epsA, (nFreq)) 276 277#ifdef HDF5 278 279 if (freq_dep.eq.0) then 280 SAFE_ALLOCATE(eps, (nmtx)) 281 call read_eps_matrix_col_hdf5(eps,igp,nmtx,mq+1,1,infile) 282 epsStatic=eps(ig) 283 SAFE_DEALLOCATE(eps) 284 endif 285 if (freq_dep.eq.2) then 286 SAFE_ALLOCATE(epsRTemp, (nFreq,nmtx)) 287#ifdef CPLX 288 SAFE_ALLOCATE(epsATemp, (nFreq,nmtx)) 289 call read_eps_matrix_col_f_hdf5(epsRTemp,nFreq,igp,nmtx,mq+1,1,infile,advanced=epsATemp) 290 epsA(:)=epsATemp(:,ig) 291 SAFE_DEALLOCATE(epsATemp) 292#else 293 call read_eps_matrix_col_f_hdf5(epsRTemp,nFreq,igp,nmtx,mq+1,1,infile) 294#endif 295 epsR(:)=epsRTemp(:,ig) 296 SAFE_DEALLOCATE(epsRTemp) 297 epsStatic=epsR(1) 298 endif 299 300#else 301 302 if (freq_dep.eq.0) then 303 SAFE_ALLOCATE(eps, (nmtx)) 304 do j=1,nmtx 305 if (j.eq.igp) then 306 read(iunit) (eps(i),i=1,nmtx) 307 epsStatic=eps(ig) 308 else 309 read(iunit) 310 endif 311 enddo 312 SAFE_DEALLOCATE(eps) 313 endif 314 if (freq_dep.eq.2) then 315 do j=1,nmtx 316 do i=1,nmtx 317 if (j.eq.igp.and.i.eq.ig) then 318 read(iunit) (epsR(k),k=1,nFreq) 319 epsStatic=epsR(1) 320 else 321 read(iunit) 322 endif 323 enddo 324#ifdef CPLX 325 do i=1,nmtx 326 if (j.eq.igp.and.i.eq.ig) then 327 read(iunit) (epsA(k),k=1,nFreq) 328 else 329 read(iunit) 330 endif 331 enddo 332#endif 333 enddo 334#ifndef CPLX 335 do i=1,nFreq 336 epsA(i)=CONJG(epsR(i)) 337 enddo 338#endif 339 endif 340 341 call close_file(iunit) 342 343#endif 344 345 SAFE_DEALLOCATE(isrtx) 346 SAFE_DEALLOCATE(isorti) 347 348 SAFE_DEALLOCATE(qpt) 349 SAFE_DEALLOCATE(gvec) 350 351 write(6,'(a,i6)') " omega num = ", nFreq 352 write(6,'(a,f7.3)') " omega step = ", dDeltaFreq 353 write(6,'(a,f7.3,/)') " omega brd = ", dBrdning 354 355!----------------------------- 356! read rho file 357 358 write(6,'(1x,"reading",1x,a,1x,"file",/)')trim(fnrho) 359 360 iunit=95 361 call open_file(unit=iunit,file=trim(fnrho),form='unformatted',status='old') 362 363 rho0=0.0d0 364 rhogmgp=ZERO 365 366 read(iunit) 367 read(iunit) nspin,nvecs 368 read(iunit) 369 read(iunit) celvol 370 read(iunit) reccelvol, blat, ((bvec(i,j),i=1,3),j=1,3), ((bdot(i,j),i=1,3),j=1,3) 371 read(iunit) 372 read(iunit) 373 read(iunit) 374 SAFE_ALLOCATE(gvec, (3, nvecs)) 375 SAFE_ALLOCATE(xcdum, (nvecs, nspin)) 376 read(iunit) nproc_para 377 ig = 1 378 do i=1,nproc_para 379 read(iunit) num_gvec 380 read(iunit) ((gvec(k, j), k = 1, 3), j = ig, ig + num_gvec - 1) 381 ig = ig + num_gvec 382 enddo 383 read(iunit) nproc_para 384 ig = 1 385 do i=1,nproc_para 386 read(iunit) num_gvec 387 read(iunit) ((xcdum(j, k), j = ig, ig + num_gvec - 1), k = 1, nspin) 388 ig = ig + num_gvec 389 enddo 390 do j=1,nvecs 391 if (gvec(1,j).eq.0.and.gvec(2,j).eq.0.and.gvec(3,j).eq.0) then 392 do k=1,nspin 393 rho0=rho0+dble(xcdum(j,k)) 394 enddo 395 endif 396 if (gvec(1,j).eq.gmgp(1).and.gvec(2,j).eq.gmgp(2).and.gvec(3,j).eq.gmgp(3)) then 397 do k=1,nspin 398 rhogmgp=rhogmgp+xcdum(j,k) 399 enddo 400 endif 401 enddo 402 SAFE_DEALLOCATE(gvec) 403 SAFE_DEALLOCATE(xcdum) 404 405 if (abs(rho0).le.TOL_Zero) then 406 call die("cannot find rho(0) in file " // trim(fnrho)) 407 endif 408 if (abs(rhogmgp).le.TOL_Zero) then 409 call die("cannot find rho(G-G') in file " // trim(fnrho)) 410 endif 411 412 call close_file(iunit) 413 414 write(6,'(5x,"cel vol =",e20.12)') celvol 415 write(6,'(5x,"rho(0) =",f10.3,/)') rho0 416 417!----------------------------- 418! construct generalized plasmon pole model 419 420 write(6,'(1x,"constructing GPP model",/)') 421 422 wp2=ryd*ryd*16.0d0*PI_D*rho0/celvol 423 qg(:)=q(:)+dble(g(:)) 424 qgp(:)=q(:)+dble(gp(:)) 425 qgqg=dot_product(qg,matmul(bdot,qg)) 426 qgqgp=dot_product(qg,matmul(bdot,qgp)) 427 if (abs(qgqg) .lt. TOL_Zero) call die("GPP model diverges") 428 Omega2=wp2*qgqgp/qgqg*rhogmgp/rho0 429#ifdef CPLX 430 epsggp = epsStatic 431#else 432 epsggp = dble(epsStatic) 433#endif 434 if (all(g(1:3) .eq. gp(1:3))) then 435 delta = 1.0d0 436 else 437 delta = 0.0d0 438 endif 439 I_epsggp = delta - epsggp 440 if (abs(I_epsggp) .lt. TOL_Small) call die("GPP model diverges") 441#ifdef CPLX 442! Complex GPP [PRB 40, 3162 (1989)] 443 wtilde2_temp = Omega2 / I_epsggp 444 lambda = abs(wtilde2_temp) 445 if (lambda .lt. TOL_Small) call die("GPP model diverges") 446 phi = atan2(IMAG(wtilde2_temp), dble(wtilde2_temp)) 447 if (abs(cos(phi)) .lt. TOL_Small) call die("GPP model diverges") 448 wtilde2 = lambda / cos(phi) 449 Omega2 = Omega2 * CMPLX(1.0d0, -tan(phi)) 450#else 451! Real GPP [PRB 34, 5390 (1986)] 452 wtilde2 = Omega2 / I_epsggp 453 if (abs(wtilde2) .lt. TOL_Small) call die("GPP model diverges") 454#endif 455! wtilde=dble(sqrt(COMPLEXIFY(wtilde2))) 456! ampl=-0.5d0*PI_D*sqrt(COMPLEXIFY(Omega2/wtilde2)) 457 do i=1,nFreq 458 omega=dFreqGrid(i) 459! epsPPRe=delta+dble(Omega2/(omega**2-wtilde2)) 460! epsPPIm=dble(ampl)/sqrt(PI_D)/dBrdning* & 461! (exp(-(omega-wtilde)**2/dBrdning**2) & 462! -exp(-(omega+wtilde)**2/dBrdning**2)) 463! epsPP(i)=CMPLX(epsPPRe,epsPPIm) 464! Instead of using the above, we write real and imaginary in a compact 465! form by adding dbrdning in denominator. The imaginary part should be 466! a sharp (delta function) at the plasmon frequency. 467 epsPP(i)=delta+Omega2/(omega**2-wtilde2-CMPLX(0D0,dBrdning)) 468 enddo 469 470 write(6,'(5x,"plasma frequency =",f10.3," eV",/)') sqrt(wp2) 471 472!----------------------------- 473! write generalized plasmon pole file 474 475 write(6,'(1x,"writing",1x,a,1x,"file",/)')trim(fngpp) 476 477 call open_file(unit=7,file=trim(fngpp),form='formatted',status='replace') 478 479 do i=1,nFreq 480 omega=dFreqGrid(i) 481 write(7,100)omega,epsPP(i) 482 enddo 483 484 call close_file(7) 485 486!----------------------------- 487! write full frequency files 488 489 if (freq_dep.eq.2) then 490 write(6,'(1x,"writing",1x,a,1x,"and",1x,a,1x,"files",/)') & 491 trim(fnffr),trim(fnffa) 492 493 call open_file(unit=8,file=trim(fnffr),form='formatted',status='replace') 494 call open_file(unit=9,file=trim(fnffa),form='formatted',status='replace') 495 496 eps_static=epsR(1) 497 eps_dynamic=eps_static 498 499 do i=1,nFreq 500 omega=dFreqGrid(i)+dFreqBrd(i) 501 write(8,100)omega,epsR(i) 502 write(9,100)omega,epsA(i) 503 if (i .gt. 1) then 504 eps_dynamic=eps_dynamic-(2.0d0/PI_D)*(dFreqGrid(i)-dFreqGrid(i-1))* & 505 IMAG(epsR(i))*((1.0d0/omega)-1.0d0/(omega+ebind)) 506 endif 507 enddo 508 write(6,*) 'Static Head:', eps_static, 'Dynamic Head:', eps_dynamic 509 write(6,*) 510 511 call close_file(8) 512 call close_file(9) 513 endif 514 515!----------------------------- 516! deallocate and finish 517 518 SAFE_DEALLOCATE(epsPP) 519 SAFE_DEALLOCATE(epsR) 520 SAFE_DEALLOCATE(epsA) 521 SAFE_DEALLOCATE(dFreqGrid) 522 SAFE_DEALLOCATE(dFreqBrd) 523 524100 format(4f25.15) 525 526end program epsinvomega 527