1!=============================================================================== 2! 3! Utilities: 4! 5! (1) epsomega Originally By gsm Last Modified 11/18/2010 (gsm) 6! 7! Parallel code that plots frequency dependence of Plasmon-Pole, Retarded, 8! Advanced epsilon for a given q, G, G` vectors. Input parameters are read 9! from file epsomega.inp in the working directory. 10! 11! Note that epsPP constructed with full-frequency epsmat will be 12! different from epsPP 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 epsomega 27 28 use global_m 29 use input_utils_m 30 use inversion_m 31 use read_matrix_m 32 use scalapack_m 33 use epsread_hdf5_m 34#ifdef HDF5 35 use hdf5 36#endif 37 38 implicit none 39 40 logical :: sflag 41 integer :: i,j,k,iq,mq,ig,igp,jg,jgp,iunit,ierr 42 integer :: icurr,icol,irow,icolm,irowm,icols,irows 43 integer :: iout,nFFTgridpts,FFTgrid(3) 44 real(DP) :: omega 45 46 character*256, parameter :: fninp = "epsomega.inp" 47 character*256 :: fneps,fnrho,fngpp,fnffr,fnffa 48 real(DP) :: q(3) 49 integer :: g(3),gp(3) 50 51 character :: infile*80 52 53 type(gspace) :: gvec 54 integer :: freq_dep,nFreq,nfreq_imag,ii,nq,ng,nmtx,kgrid(3) 55 real(DP) :: dDeltaFreq,dBrdning,ecuts,delta,qvec(3) 56 real(DP), allocatable :: qpt(:,:) 57 real(DP), allocatable :: ekin(:) 58 real(DP), allocatable :: dFreqGrid(:) 59 integer, allocatable :: isrtx(:) 60 integer, allocatable :: isorti(:) 61 integer, allocatable :: nmtx_of_q(:) 62 SCALAR, allocatable :: eps(:,:) 63 complex(DPC), allocatable :: dFreqBrd(:) 64 complex(DPC), allocatable :: epsPP(:,:,:) 65 complex(DPC), allocatable :: epsR(:,:,:) 66 complex(DPC), allocatable :: epsA(:,:,:) 67 complex(DPC), allocatable :: epsAux(:,:) 68 SCALAR, allocatable :: rhog(:) 69 real(DP) :: bvec(3,3),bdot(3,3),blat,celvol,reccelvol 70 integer :: nvecs,nspin 71 72 integer :: nproc_para,num_gvec,gx,gy,gz,qgrid(3),error 73 real(DP) :: rho0 74 integer, allocatable :: gvec_rho(:,:) 75 SCALAR, allocatable :: xcdum(:,:) 76 77 real(DP) :: wp2,qg(3),qgp(3),qgqg,qgqgp,lambda,phi 78 SCALAR :: Omega2,wtilde2,epsggp,I_epsggp 79 complex(DPC) :: wtilde2_temp 80 81 type (scalapack) :: scal 82 83 call peinfo_init() 84 85!----------------------------- 86! read input file 87 88 if (peinf%inode.eq.0) then 89 write(6,'(/,1x,"reading",1x,a,1x,"file",/)')trim(fninp) 90 91 call open_file(55,file=trim(fninp),form='formatted',status='old') 92 93 read(55,'(a)') fneps 94 read(55,'(a)') fnrho 95 read(55,*) (q(i),i=1,3) 96 read(55,*) (g(i),i=1,3) 97 read(55,*) (gp(i),i=1,3) 98 read(55,*) nFreq 99 read(55,*) dDeltaFreq 100 read(55,*) dBrdning 101 read(55,'(a)') fngpp 102 read(55,'(a)') fnffr 103 read(55,'(a)') fnffa 104 105 call close_file(55) 106 107 write(6,'(2a)') " eps file = ", trim(fneps) 108 write(6,'(2a)') " rho file = ", trim(fnrho) 109 write(6,'(a, 3f15.12)') " q vector = ", q(1:3) 110 write(6,'(a, 3i4)') " G vector = ", g(1:3) 111 write(6,'(a, 3i4)') " G' vector = ", gp(1:3) 112 write(6,'(a, i6)') " nFreq = ", nFreq 113 write(6,'(a, f7.3)') " dDeltaFreq = ", dDeltaFreq 114 write(6,'(a, f7.3)') " dBrdning = ", dBrdning 115 write(6,'(2a)') " GPP file = ", trim(fngpp) 116 write(6,'(2a)') " FFR file = ", trim(fnffr) 117 write(6,'(2a,/)') " FFA file = ", trim(fnffa) 118 endif 119 120#ifdef MPI 121 call MPI_Bcast(fneps,256,MPI_CHARACTER,0,MPI_COMM_WORLD,mpierr) 122 call MPI_Bcast(fnrho,256,MPI_CHARACTER,0,MPI_COMM_WORLD,mpierr) 123 call MPI_Bcast(q,3,MPI_REAL_DP,0,MPI_COMM_WORLD,mpierr) 124 call MPI_Bcast(g,3,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 125 call MPI_Bcast(gp,3,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 126 call MPI_Bcast(fngpp,256,MPI_CHARACTER,0,MPI_COMM_WORLD,mpierr) 127 call MPI_Bcast(fnffr,256,MPI_CHARACTER,0,MPI_COMM_WORLD,mpierr) 128 call MPI_Bcast(fnffa,256,MPI_CHARACTER,0,MPI_COMM_WORLD,mpierr) 129#endif 130 131!----------------------------- 132! read eps file 133 134 if (peinf%inode.eq.0) then 135 write(6,'(1x,"reading",1x,a,1x,"file",/)')trim(fneps) 136 137#ifdef HDF5 138 139 ! JRD: There is no equivalent on read_matrix_d 140 ! yet. So, the beginings of HDF5 support are here. 141 ! But, it is not complete 142 call die("No Support for HDF5 Yet. epsinvomega routine does support HDF5.") 143 144 call h5open_f(error) 145 146 infile = trim(fneps) 147 148 call read_eps_grid_sizes_hdf5(ng, nq, ecuts, nfreq, nfreq_imag, nmtx, qgrid, freq_dep, infile) 149 SAFE_ALLOCATE(dFreqGrid,(nFreq)) 150 SAFE_ALLOCATE(dFreqBrd,(nFreq)) 151 SAFE_ALLOCATE(qpt, (3,nq)) 152 SAFE_ALLOCATE(nmtx_of_q, (nq)) 153 SAFE_ALLOCATE(gvec%components, (3,ng)) 154 155 call read_eps_qgrid_hdf5(nq,qpt,nmtx_of_q,infile) 156 if (freq_dep .ne. 0) then 157 call read_eps_freqgrid_hdf5(nFreq,dFreqGrid,dFreqBrd,infile) 158 endif 159 160#else 161 162 iunit=12 163 call open_file(unit=iunit,file=trim(fneps),form='unformatted',status='old') 164 165 read(iunit) 166 read(iunit) freq_dep,ii 167 if (freq_dep.ne.0) then 168 nFreq=ii 169 endif 170 read(iunit) (kgrid(i),i=1,3) 171 SAFE_ALLOCATE(dFreqGrid,(nFreq)) 172 SAFE_ALLOCATE(dFreqBrd,(nFreq)) 173 if (freq_dep.ne.0) then 174 read(iunit) (dFreqGrid(i),i=1,nFreq),(dFreqBrd(i),i=1,nFreq) 175 if (nFreq.gt.1) dDeltaFreq=dFreqGrid(2)-dFreqGrid(1) 176 dBrdning=IMAG(dFreqBrd(1)) 177 else 178 do i=1,nFreq 179 dFreqGrid(i)=dDeltaFreq*dble(i-1) 180 dFreqBrd(i)=dBrdning 181 enddo 182 read(iunit) 183 endif 184 read(iunit) 185 read(iunit) 186 read(iunit) ecuts 187 read(iunit) nq 188 read(iunit) ng 189 190 rewind(iunit) 191 192 SAFE_ALLOCATE(qpt, (3,nq)) 193 SAFE_ALLOCATE(gvec%components, (3,ng)) 194 195 read(iunit) 196 read(iunit) 197 read(iunit) 198 read(iunit) 199 read(iunit) 200 read(iunit) 201 read(iunit) 202 read(iunit) nq,((qpt(j,i),j=1,3),i=1,nq) 203 read(iunit) ng,((gvec%components(j,i),j=1,3),i=1,ng) 204 205#endif 206 207 ig=0 208 igp=0 209 do i=1,ng 210 if (gvec%components(1,i).eq.g(1).and.gvec%components(2,i).eq.g(2).and. & 211 gvec%components(3,i).eq.g(3)) ig=i 212 if (gvec%components(1,i).eq.gp(1).and.gvec%components(2,i).eq.gp(2).and. & 213 gvec%components(3,i).eq.gp(3)) igp=i 214 enddo 215 216 mq=-1 217 do iq=1,nq 218 if (abs(q(1)-qpt(1,iq)).lt.TOL_Zero.and. & 219 abs(q(2)-qpt(2,iq)).lt.TOL_Zero.and. & 220 abs(q(3)-qpt(3,iq)).lt.TOL_Zero) then 221 mq=iq-1 222 exit 223 endif 224 enddo 225 226 SAFE_DEALLOCATE(qpt) 227 228 write(6,'(a,i6)') " omega num = ", nFreq 229 endif 230 231#ifdef MPI 232 call MPI_Bcast(freq_dep,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 233 call MPI_Bcast(nFreq,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 234 call MPI_Bcast(dDeltaFreq,1,MPI_REAL_DP,0,MPI_COMM_WORLD,mpierr) 235 call MPI_Bcast(dBrdning,1,MPI_REAL_DP,0,MPI_COMM_WORLD,mpierr) 236 call MPI_Bcast(mq,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 237 call MPI_Bcast(ig,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 238 call MPI_Bcast(igp,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 239 call MPI_Bcast(ng,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 240 if (peinf%inode.ne.0) then 241 SAFE_ALLOCATE(dFreqGrid,(nFreq)) 242 SAFE_ALLOCATE(dFreqBrd,(nFreq)) 243 SAFE_ALLOCATE(gvec%components, (3,ng)) 244 endif 245 call MPI_Bcast(dFreqGrid,nFreq,MPI_REAL_DP,0,MPI_COMM_WORLD,mpierr) 246 call MPI_Bcast(dFreqBrd,nFreq,MPI_COMPLEX_DPC,0,MPI_COMM_WORLD,mpierr) 247 call MPI_Bcast(gvec%components,3*ng,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 248#endif 249 250 if (mq.eq.-1) then 251 call die("cannot find q vector in file " // trim(fneps)) 252 endif 253 if (ig.eq.0) then 254 call die("cannot find G vector in file " // trim(fneps)) 255 endif 256 if (igp.eq.0) then 257 call die("cannot find G' vector in file " // trim(fneps)) 258 endif 259 260 if (peinf%inode.eq.0) then 261 262#ifdef HDF5 263 264 SAFE_ALLOCATE(isrtx, (ng)) 265 SAFE_ALLOCATE(isorti, (ng)) 266 SAFE_ALLOCATE(ekin, (ng)) 267 call read_eps_gvecsofq_hdf5(ng,isrtx,isorti,ekin,mq+1,infile) 268 nmtx = nmtx_of_q(mq+1) 269 SAFE_DEALLOCATE(ekin) 270 SAFE_DEALLOCATE(isorti) 271 SAFE_DEALLOCATE(nmtx_of_q) 272 273#else 274 275 do iq=1,mq 276 read(iunit) ng,nmtx 277 read(iunit) 278 read(iunit) 279 if (freq_dep.eq.0) then 280 do j=1,nmtx 281 read(iunit) 282 enddo 283 else 284 do j=1,nmtx 285 do i=1,nmtx 286 read(iunit) 287#ifdef CPLX 288 if(freq_dep.ne.3) then 289 read(iunit) 290 endif 291#endif 292 enddo 293 enddo 294 endif 295 enddo 296 297 SAFE_ALLOCATE(isrtx, (ng)) 298 SAFE_ALLOCATE(isorti, (ng)) 299 isrtx(:)=0 300 isorti(:)=0 301 302 read(iunit) ng,nmtx,(isrtx(i),isorti(i),i=1,ng) 303 read(iunit) 304 read(iunit) 305 306 ig=isorti(ig) 307 igp=isorti(igp) 308 309 SAFE_DEALLOCATE(isorti) 310 311#endif 312 313 endif 314 315#ifdef MPI 316 call MPI_Bcast(nmtx,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 317 call MPI_Bcast(ig,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 318 call MPI_Bcast(igp,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 319 if (peinf%inode.ne.0) then 320 SAFE_ALLOCATE(isrtx, (ng)) 321 endif 322 call MPI_Bcast(isrtx,ng,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 323#endif 324 325 if (ig.eq.0) then 326 call die("cannot find G vector in file " // trim(fneps)) 327 endif 328 if (igp.eq.0) then 329 call die("cannot find G' vector in file " // trim(fneps)) 330 endif 331 332 call blacs_setup(scal, nmtx, .true.) 333 334!----------------------------- 335! read distributed matrices 336 337 SAFE_ALLOCATE(eps, (scal%npr,scal%npc)) 338 SAFE_ALLOCATE(epsPP, (nFreq,scal%npr,scal%npc)) 339 if (freq_dep.ne.0) then 340 SAFE_ALLOCATE(epsR, (nFreq,scal%npr,scal%npc)) 341 SAFE_ALLOCATE(epsA, (nFreq,scal%npr,scal%npc)) 342 endif 343 SAFE_ALLOCATE(epsAux, (scal%npr,scal%npc)) 344 345 if (freq_dep.eq.0) then 346 call read_matrix_d(scal,eps,nmtx,iunit) 347 else 348 peinf%rank_mtxel=0 349 call read_matrix_f(scal, nFreq, nFreq, epsR, nmtx, 1, iunit, advanced=epsA) 350#ifndef CPLX 351 do icolm=1,scal%npc 352 do irowm=1,scal%npr 353 do i=1,nFreq 354 epsA(i,irowm,icolm)=CONJG(epsR(i,irowm,icolm)) 355 enddo 356 enddo 357 enddo 358#else 359 if (freq_dep.eq.3) then 360 do icolm=1,scal%npc 361 do irowm=1,scal%npr 362 do i=1,nFreq 363 epsA(i,irowm,icolm)=CONJG(epsR(i,irowm,icolm)) 364 enddo 365 enddo 366 enddo 367 endif 368#endif 369 do icolm=1,scal%npc 370 do irowm=1,scal%npr 371#ifdef CPLX 372 eps(irowm,icolm)=epsR(1,irowm,icolm) 373#else 374 eps(irowm,icolm)=dble(epsR(1,irowm,icolm)) 375#endif 376 enddo 377 enddo 378 endif 379 380 if (peinf%inode.eq.0) then 381 call close_file(iunit) 382 endif 383 384 if (peinf%inode.eq.0) then 385 write(6,'(a,i6)') " omega num = ", nFreq 386 write(6,'(a,f7.3)') " omega step = ", dDeltaFreq 387 write(6,'(a,f7.3,/)') " omega brd = ", dBrdning 388 endif 389 390!----------------------------- 391! read rho file 392 393 SAFE_ALLOCATE(rhog, (ng)) 394 rhog(:)=ZERO 395 396 if (peinf%inode.eq.0) then 397 write(6,'(1x,"reading",1x,a,1x,"file",/)')trim(fnrho) 398 399 iunit=95 400 call open_file(unit=iunit,file=trim(fnrho),form='unformatted',status='old') 401 402 rho0=0.0d0 403 404 read(iunit) 405 read(iunit) nspin,nvecs 406 read(iunit) (FFTgrid(i),i=1,3) 407 endif 408! (gsm) this is ugly but the only way I see to quickly fix the mess 409! with gvec structure and gvec_index subroutine 410#ifdef MPI 411 call MPI_Bcast(FFTgrid,3,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 412#endif 413 gvec%ng=ng 414 gvec%FFTgrid=FFTgrid 415 call gvec_index(gvec) 416 nFFTgridpts=gvec%nFFTgridpts 417 if (peinf%inode.eq.0) then 418 read(iunit) celvol 419 read(iunit) reccelvol, blat, ((bvec(i,j),i=1,3),j=1,3), ((bdot(i,j),i=1,3),j=1,3) 420 read(iunit) 421 read(iunit) 422 read(iunit) 423 SAFE_ALLOCATE(gvec_rho, (3, nvecs)) 424 SAFE_ALLOCATE(xcdum, (nvecs, nspin)) 425 read(iunit) nproc_para 426 jg = 1 427 do i=1,nproc_para 428 read(iunit) num_gvec 429 read(iunit) ((gvec_rho(k, j), k = 1, 3), j = jg, jg + num_gvec - 1) 430 jg = jg + num_gvec 431 enddo 432 read(iunit) nproc_para 433 jg = 1 434 do i=1,nproc_para 435 read(iunit) num_gvec 436 read(iunit) ((xcdum(j, k), j = jg, jg + num_gvec - 1), k = 1, nspin) 437 jg = jg + num_gvec 438 enddo 439 ierr=0 440 do j=1,nvecs 441 if (gvec_rho(1,j).eq.0.and.gvec_rho(2,j).eq.0.and.gvec_rho(3,j).eq.0) then 442 do k=1,nspin 443 rho0=rho0+dble(xcdum(j,k)) 444 enddo 445 endif 446 iout=((gvec_rho(1,j)+FFTgrid(1)/2)*FFTgrid(2)+gvec_rho(2,j)+FFTgrid(2)/2)* & 447 FFTgrid(3)+gvec_rho(3,j)+FFTgrid(3)/2+1 448 if (iout.ge.1.and.iout.le.nFFTgridpts) then 449 iout=gvec%index_vec(iout) 450 if (iout.ge.1.and.iout.le.ng) then 451 rhog(iout)=ZERO 452 do k=1,nspin 453 rhog(iout)=rhog(iout) + xcdum(j,k) 454 enddo 455 else 456 ierr=1 457 endif 458 else 459 ierr=1 460 endif 461 enddo 462 SAFE_DEALLOCATE(gvec_rho) 463 SAFE_DEALLOCATE(xcdum) 464 465 call close_file(iunit) 466 467 write(6,'(5x,"cel vol =",e20.12)') celvol 468 write(6,'(5x,"rho(0) =",f10.3,/)') rho0 469 endif 470 471#ifdef MPI 472 call MPI_Bcast(bdot,9,MPI_REAL_DP,0,MPI_COMM_WORLD,mpierr) 473 call MPI_Bcast(celvol,1,MPI_REAL_DP,0,MPI_COMM_WORLD,mpierr) 474 call MPI_Bcast(FFTgrid,3,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 475 call MPI_Bcast(ierr,1,MPI_INTEGER,0,MPI_COMM_WORLD,mpierr) 476 call MPI_Bcast(rho0,1,MPI_REAL_DP,0,MPI_COMM_WORLD,mpierr) 477 call MPI_Bcast(rhog,ng,MPI_SCALAR,0,MPI_COMM_WORLD,mpierr) 478#endif 479 480 if (ierr.ne.0) then 481 call die("unknown G vector in file " // trim(fnrho)) 482 endif 483 if (abs(rho0).le.TOL_Zero) then 484 call die("cannot find rho(0) in file " // trim(fnrho)) 485 endif 486 487!----------------------------- 488! construct generalized plasmon pole model 489 490 if (peinf%inode.eq.0) then 491 write(6,'(1x,"constructing GPP model",/)') 492 endif 493 494 epsPP(:,:,:)=(0.0d0,0.0d0) 495 wp2=ryd*ryd*16.0d0*PI_D*rho0/celvol 496 sflag=.false. 497 irows=0 498 icols=0 499 icurr=0 500 do jgp=1,nmtx 501 icol=mod(int(((jgp-1)/scal%nbl)+TOL_Small),scal%npcol) 502 do jg=1,nmtx 503 irow=mod(int(((jg-1)/scal%nbl)+TOL_Small),scal%nprow) 504 if (irow.eq.scal%myprow.and.icol.eq.scal%mypcol) then 505 icurr=icurr+1 506 icolm=int((icurr-1)/scal%npr+TOL_Small)+1 507 irowm=mod((icurr-1),scal%npr)+1 508 if (jg.eq.ig.and.jgp.eq.igp) then 509 sflag=.true. 510 irows=irowm 511 icols=icolm 512 endif 513 qg(:)=q(:)+dble(gvec%components(:,isrtx(jg))) 514 qgp(:)=q(:)+dble(gvec%components(:,isrtx(jgp))) 515 qgqg=dot_product(qg,matmul(bdot,qg)) 516 qgqgp=dot_product(qg,matmul(bdot,qgp)) 517 if (abs(qgqg).lt.TOL_Zero) cycle 518 gx=gvec%components(1,isrtx(jg))-gvec%components(1,isrtx(jgp)) 519 gy=gvec%components(2,isrtx(jg))-gvec%components(2,isrtx(jgp)) 520 gz=gvec%components(3,isrtx(jg))-gvec%components(3,isrtx(jgp)) 521 iout=((gx+FFTgrid(1)/2)*FFTgrid(2)+gy+FFTgrid(2)/2)* & 522 FFTgrid(3)+gz+FFTgrid(3)/2+1 523 if (iout.lt.1.or.iout.gt.nFFTgridpts) cycle 524 iout=gvec%index_vec(iout) 525 if (iout.lt.1.or.iout.gt.ng) cycle 526 Omega2=wp2*qgqgp/qgqg*rhog(iout)/rho0 527 epsggp=eps(irowm,icolm) 528 if (jg.eq.jgp) then 529 delta=1.0d0 530 else 531 delta=0.0d0 532 endif 533 I_epsggp=delta-epsggp 534 if (abs(I_epsggp).lt.TOL_Small) cycle 535#ifdef CPLX 536! Complex GPP [PRB 40, 3162 (1989)] 537 wtilde2_temp = Omega2 / I_epsggp 538 lambda = abs(wtilde2_temp) 539 if (lambda .lt. TOL_Small) cycle 540 phi = atan2(IMAG(wtilde2_temp), dble(wtilde2_temp)) 541 if (abs(cos(phi)) .lt. TOL_Small) cycle 542 wtilde2 = lambda / cos(phi) 543 Omega2 = Omega2 * CMPLX(1.0d0, -tan(phi)) 544#else 545! Real GPP [PRB 34, 5390 (1986)] 546 wtilde2 = Omega2 / I_epsggp 547 if (abs(wtilde2) .lt. TOL_Small) cycle 548#endif 549 do i=1,nFreq 550 omega=dFreqGrid(i) 551 epsPP(i,irowm,icolm)=delta+Omega2/(omega**2-wtilde2-CMPLX(0D0,dBrdning)) 552 enddo 553 endif 554 enddo 555 enddo 556 557 if (peinf%inode.eq.0) then 558 write(6,'(5x,"plasma frequency =",f10.3," eV",/)') sqrt(wp2) 559 endif 560 561!----------------------------- 562! invert matrices 563 564 if (peinf%inode.eq.0) then 565 write(6,'(1x,"inverting matrices",/)') 566 endif 567 568#if defined USESCALAPACK 569 do i=1,nFreq 570 epsAux(:,:) = epsPP(i,:,:) 571 call zinvert_with_scalapack(nmtx, scal, epsAux) 572 epsPP(i,:,:) = epsAux(:,:) 573 enddo 574 if(freq_dep .ne. 0) then 575 do i=1,nFreq 576 epsAux(:,:) = epsR(i,:,:) 577 call zinvert_with_scalapack(nmtx, scal, epsAux) 578 epsR(i,:,:) = epsAux(:,:) 579 epsAux(:,:) = epsA(i,:,:) 580 call zinvert_with_scalapack(nmtx, scal, epsAux) 581 epsA(i,:,:) = epsAux(:,:) 582 enddo 583 endif 584#else 585 do i=1,nFreq 586 epsAux(:,:) = epsPP(i,:,:) 587 call zinvert_serial(nmtx, epsAux) 588 epsPP(i,:,:) = epsAux(:,:) 589 enddo 590 if(freq_dep .ne. 0) then 591 do i=1,nFreq 592 epsAux(:,:) = epsR(i,:,:) 593 call zinvert_serial(nmtx, epsAux) 594 epsR(i,:,:) = epsAux(:,:) 595 epsAux(:,:) = epsA(i,:,:) 596 call zinvert_serial(nmtx, epsAux) 597 epsA(i,:,:) = epsAux(:,:) 598 enddo 599 endif 600#endif 601 602!----------------------------- 603! write generalized plasmon pole file 604 605 if (sflag) then 606 write(6,'(1x,"writing",1x,a,1x,"file",/)')trim(fngpp) 607 608 call open_file(unit=7,file=trim(fngpp),form='formatted',status='replace') 609 610 do i=1,nFreq 611 omega=dFreqGrid(i) 612 write(7,100)omega,epsPP(i,irows,icols) 613 enddo 614 615 call close_file(7) 616 endif 617 618!----------------------------- 619! write full frequency files 620 621 if (sflag) then 622 if (freq_dep.ne.0) then 623 write(6,'(1x,"writing",1x,a,1x,"and",1x,a,1x,"files",/)') & 624 trim(fnffr),trim(fnffa) 625 626 call open_file(unit=8,file=trim(fnffr),form='formatted',status='replace') 627 call open_file(unit=9,file=trim(fnffa),form='formatted',status='replace') 628 629 do i=1,nFreq 630 omega=dFreqGrid(i) 631 write(8,100)omega,epsR(i,irows,icols) 632 write(9,100)omega,epsA(i,irows,icols) 633 enddo 634 635 call close_file(8) 636 call close_file(9) 637 endif 638 endif 639 640!----------------------------- 641! deallocate and finish 642 643 SAFE_DEALLOCATE_P(gvec%components) 644 SAFE_DEALLOCATE(isrtx) 645 SAFE_DEALLOCATE_P(gvec%index_vec) 646 SAFE_DEALLOCATE(eps) 647 SAFE_DEALLOCATE(epsPP) 648 SAFE_DEALLOCATE(dFreqGrid) 649 SAFE_DEALLOCATE(dFreqBrd) 650 if (freq_dep.ne.0) then 651 SAFE_DEALLOCATE(epsR) 652 SAFE_DEALLOCATE(epsA) 653 endif 654 SAFE_DEALLOCATE(epsAux) 655 SAFE_DEALLOCATE(rhog) 656 657#ifdef MPI 658 call MPI_Finalize(mpierr) 659#endif 660 661100 format(3f25.15) 662 663end program epsomega 664