! ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino ! ! This file is distributed under the terms of the GNU General Public ! License. See the file `LICENSE' in the root directory of the ! present distribution, or http://www.gnu.org/copyleft.gpl.txt . ! !----------------------------------------------------------------------- SUBROUTINE indabs_fly (iq) !----------------------------------------------------------------------- ! ! Calculate the phonon-assisted indirect absorption spectra ! Will be able to USE either the dipole matrix elements ! or the more correct velocity matrix elements ! 07/2010 JDN and Emmanouil Kioupakis ! !----------------------------------------------------------------------- #include "f_defs.h" USE kinds, ONLY : DP USE cell_base, ONLY : at, bg USE io_global, ONLY : stdout USE io_files, ONLY : prefix, find_free_unit USE phcom, ONLY : nmodes USE epwcom, ONLY : nbndsub, lrepmatf, iunepmatf, fsthick, & eptemp, ngaussw, degaussw, iuetf, & wmin, wmax, nw, nbndskip, a2f, epf_mem, & nsmear, delta_smear, eig_read, vme USE pwcom, ONLY : nelec, ef, isk, nbnd USE el_phon, ONLY : epf17, ibndmax, ibndmin, etf, & etfq, wkf, xqf, wqf, nksf, nxqf, & nksqf, wf, nkstotf, xkf, xqf, dmef, vmef, spectra #ifdef __PARA USE mp, ONLY : mp_barrier,mp_sum USE io_global, ONLY : ionode_id USE mp_global, ONLY : me_pool,inter_pool_comm,my_pool_id, mpime #endif ! implicit none ! real(kind=DP), PARAMETER :: eps = 20.d0 / 13.6058d0 / 8065.5d0, & ryd2mev = 13605.8, one = 1.d0, ryd2ev = 13.6058, & two = 2.d0, zero = 0.d0, pi = 3.14159265358979 complex(kind = 8), PARAMETER :: ci = (0.d0, 1.d0), cone = (1.d0, 0.d0), czero=(0.d0,0.d0) integer :: ik, ikk, ikq, ibnd, jbnd, mbnd, imode, nrec, & iomega, nomega, ipol, iq, ismear, iuaugerfil complex(kind=DP) epf (ibndmax-ibndmin+1, ibndmax-ibndmin+1), s1(3), s2(3) real(kind=DP) :: ekk, ekq, wq,ef0, wgkk, wgkq, & weight, wgauss, degaussw0, pfac, w0gauss, & nqv, cfac, delta_omega, omega, wg(nbndsub, nksf) ! real(kind=DP), ALLOCATABLE :: spectra_q(:,:,:) real(kind=DP), EXTERNAL :: efermig, efermit ! ! ! CALL start_clock('ABS SPECTRA') ! WRITE (6,*) "integrating :", iq, nxqf CALL flush(6) ! ! Fermi level and corresponding DOS ! ! Note that the weights of k+q points must be set to zero here ! CALL iweights(nksf, wkf, nbndsub, nelec, etf, ef0, wg, 0, isk) WRITE (6, 101) ef0 * ryd2ev ! ! if 'fine' Fermi level differs by more than 0.5 eV, there is probably ! something wrong with the Wannier functions IF (abs(ef0 - ef) * ryd2eV .gt. 0.5 .and. (.not.eig_read) ) & CALL errore ('selfen_phon', 'Something wrong with Ef, check MLWFs', 1) ! ! this should be an input PARAMETER nomega = 100 ! IF (iq.eq.1) THEN CALL epw_dos(ef0) ALLOCATE ( spectra(3, nomega, nsmear) ) spectra = zero WRITE(6,'(/5x,a)') REPEAT('=',67) WRITE(6,'(5x,"Indirect absorption spectra")') WRITE(6,'(5x,a/)') REPEAT('=',67) ! IF (nbndskip.gt.0) THEN ! ! here we take into account that we may skip bands when we wannierize ! (spin-unpolarized) ! nelec = nelec - two * nbndskip WRITE(stdout, '(/5x,"Skipping the first ",i4," bands:")' ) nbndskip WRITE(stdout, '(/5x,"The Fermi level will be determined with ",f9.5," electrons")' ) nelec ENDIF ENDIF ! IF (.not. ALLOCATEd(spectra_q)) ALLOCATE(spectra_q(3,nomega,nsmear)) spectra_q = zero ! ! DO iomega = 1, nomega WRITE (6,'(5x,a,2i5)') " iomega of nomega", iomega, nomega CALL flush(6) ! 50 meV delta_omega = 25 / ryd2mev omega = float(iomega) * delta_omega ! ! DO ik = 1, nksqf ! parallelized over these wavevectors ! ikk = 2 * ik - 1 ikq = ikk + 1 ! ! we read the hamiltonian eigenvalues (those at k+q depend on q!) ! IF (epf_mem) THEN etf (ibndmin:ibndmax, ikk) = etfq (ibndmin:ibndmax, ikk, 1) etf (ibndmin:ibndmax, ikq) = etfq (ibndmin:ibndmax, ikq, 1) ELSE nrec = ikk CALL davcio ( etf (ibndmin:ibndmax, ikk), ibndmax-ibndmin+1, iuetf, nrec, - 1) nrec = ikq CALL davcio ( etf (ibndmin:ibndmax, ikq), ibndmax-ibndmin+1, iuetf, nrec, - 1) ENDIF ! ! DO imode = 1, nmodes ! ! the phonon frequency wq = wf (imode, iq) ! ! we read the e-p matrix from disk / memory ! IF (epf_mem) THEN epf(:,:) = epf17 ( ik, 1, :, :, imode) ELSE nrec = (imode-1) * nksqf + ik CALL dasmio ( epf, ibndmax-ibndmin+1, lrepmatf, iunepmatf, nrec, -1) ENDIF ! ! DO ibnd = 1, ibndmax-ibndmin+1 ! ! the fermi occupation for k ekk = etf (ibndmin-1+ibnd, ikk) - ef0 wgkk = 0.d0 IF (ekk.lt. 0.02) wgkk = 1.d0 ! DO jbnd = 1, ibndmax-ibndmin+1 ! ! the fermi occupation for k+q ekq = etf (ibndmin-1+jbnd, ikq) - ef0 wgkq = 0.d0 IF (ekq.lt. 0.02) wgkq = 1.d0 ! ! ! P = (nqv + 0.5 +/- 0.5) * (fi(k) - fj(k+q) ) ! pfac = ( nqv + 0.5 +/- 0.5 )*( wgkk - wgkq ) ! IF (abs(wgkk - wgkq).lt.0.01) cycle ! DO ismear = 1, nsmear degaussw0 = (ismear-1) * delta_smear + degaussw ! nqv = 0.d0 IF (wq .gt. 3/ryd2mev) THEN !nqv = 1.0 / ( exp(wq/(degaussw00)) -1) nqv = wgauss( -wq/degaussw0, -99) nqv = nqv / ( one - two * nqv ) ENDIF pfac = nqv * ( wgkk - wgkq ) ! C = 4pi^2 e^2 / (nr) cm^2 ! internal units in Ry, ! C : e^2 = 2, m = 1/2, c= 1/(137.0359895*2) ! C = 4 * pi^2 * 2 *2^2 / (nr * 137*2) ! C = 16 pi^2/137 * 1/nr = 1.1523518/nr ! for Si, nr ~ 3.55 cfac = 0.324606 ! s1(:) = (0.d0, 0.d0) DO mbnd = 1, ibndmax-ibndmin+1 IF (vme) THEN s1(:) = epf(mbnd, jbnd) * 0.5 * vmef(:,ibnd, mbnd, ikk) / & ( etf(ibndmin-1+mbnd, ikk) - etf(ibndmin-1+ibnd, ikk) - omega + ci * degaussw0) ELSE s1(:) = epf(mbnd, jbnd) * dmef(:,ibnd, mbnd, ikk) / & ( etf(ibndmin-1+mbnd, ikk) - etf(ibndmin-1+ibnd, ikk) - omega + ci * degaussw0) ENDIF ENDDO ! s2(:) = (0.d0, 0.d0) DO mbnd = 1, ibndmax-ibndmin+1 IF (vme) THEN s2(:) = epf(ibnd, mbnd) * 0.5* vmef(:,mbnd, jbnd, ikq) / & ( etf(ibndmin-1+mbnd, ikq) - etf(ibndmin-1+ibnd, ikk) - wq + ci * degaussw0) ELSE s2(:) = epf(ibnd, mbnd) * dmef(:,mbnd, jbnd, ikq) / & ( etf(ibndmin-1+mbnd, ikq) - etf(ibndmin-1+ibnd, ikk) - wq + ci * degaussw0) ENDIF ! s2(:) = epf(ibnd, mbnd) * dmef(:,mbnd, jbnd, ikq) / ( etf(ibndmin-1+mbnd, ikq) - etf(ibndmin-1+ibnd, ikk) + wq + ci * degaussw0) ENDDO ! weight = w0gauss( ( ekq - ekk - omega - wq) / degaussw0, 0) / degaussw0 DO ipol = 1, 3 IF (wq .gt. 3/ryd2mev) THEN spectra_q(ipol,iomega, ismear) = spectra_q(ipol,iomega, ismear) + & 2 * cfac / omega * pfac * weight * abs( ( s1(ipol) + s2(ipol) )**2) ENDIF ENDDO ! ENDDO ! loop on smears ! ENDDO ! jbnd ENDDO ! ibnd ! ENDDO ! loop on q-modes ! ENDDO ! loop on k ! ENDDO ! omega #ifdef __PARA ! ! collect contributions from all pools (sum over k-points and q-points) ! this finishes the integral over the BZ (k,q) ! ! spectra for a given q CALL mp_sum(spectra_q,inter_pool_comm) ! #endif ! ! sum over q's, should fix the problem with parallelization split spectra(:,:,:) = spectra(:,:,:) + spectra_q(:,:,:) ! CALL stop_clock('ABS SPECTRA') ! ! print this out on the final iq IF (iq.eq.nxqf) THEN WRITE(6,'(5x,a/)') REPEAT('-',67) DO iomega = 1, nomega ! omega = iomega * delta_omega WRITE (6, '(5x,f8.4,3f24.8)') omega*ryd2ev, spectra(:, iomega, 1) ENDDO WRITE(6,'(5x,a/)') REPEAT('-',67) ! ! print out data to USEful file #ifdef __PARA IF (mpime.eq.ionode_id) THEN #endif iuaugerfil = find_free_unit() ! open (unit = iuaugerfil, file = trim(prefix)//".abs", form = 'formatted') DO ismear = 1, nsmear WRITE(iuaugerfil,'(5x,a,a/)') "# ", REPEAT('-',67) WRITE(iuaugerfil,*) "# Numerical broadening (meV)", & ryd2mev*((ismear-1) * delta_smear + eptemp) WRITE(iuaugerfil,'(a,5x,a,5x,a,5x,a)') "# E", "alphax", "alphay", "alphaz" WRITE(iuaugerfil,'(5x,a,a/)') "# ", REPEAT('-',67) DO iomega = 1, nomega omega = iomega * delta_omega WRITE (iuaugerfil, '(5x,f8.4,3f24.8)') omega*ryd2ev, spectra(:, iomega, ismear) ENDDO ENDDO close(iuaugerfil) ! #ifdef __PARA ENDIF #endif ! ENDIF ! 100 format(5x,'Gaussian Broadening: ',f7.3,' Ry, ngauss=',i4) 101 format(5x,'Ef=',f10.6,' eV') ! END SUBROUTINE indabs_fly !------------------------------------------------- SUBROUTINE epw_dos (ef) !------------------------------------------------- #include "f_defs.h" USE kinds, ONLY : DP USE io_global, ONLY : stdout USE io_files, ONLY : prefix, find_free_unit USE epwcom, ONLY : nbndsub, epf_mem, degaussw, iuetf USE pwcom, ONLY : nelec, isk, nbnd USE el_phon, ONLY : ibndmax, ibndmin, etf, nksqf, etfq, & nksf, wkf #ifdef __PARA USE mp, ONLY : mp_barrier,mp_sum, mp_min, mp_max USE io_global, ONLY : ionode_id, ionode USE mp_global, ONLY : me_pool,inter_pool_comm,my_pool_id, mpime #endif ! implicit none ! real(kind=DP), PARAMETER :: eps = 20.d0 / 13.6058d0 / 8065.5d0, & ryd2mev = 13605.8, one = 1.d0, ryd2ev = 13.6058, & two = 2.d0, zero = 0.d0, pi = 3.14159265358979 complex(kind = 8), PARAMETER :: ci = (0.d0, 1.d0), cone = (1.d0, 0.d0), czero=(0.d0,0.d0) integer :: ik, ikk, ibnd, jbnd, idos, ndos, nrec, iudosfil real(kind=DP) :: ebnd, ebndmin, ebndmax, en, w0gauss ! real(kind=DP) :: DDOT, tmp, wgauss, ef real(kind=DP), ALLOCATABLE :: fdos(:), focc(:) ! ebndmin = 1000 ebndmax = -1000 DO ik = 1, nksf DO ibnd = 1, nbndsub ebnd = etf (ibnd, ik) ! ibndmin = min(ibnd,ibndmin) ibndmax = max(ibnd,ibndmax) ebndmin = min(ebnd,ebndmin) ebndmax = max(ebnd,ebndmax) ! ENDDO ENDDO ! #ifdef __PARA tmp = float (ibndmin) CALL mp_min(tmp,inter_pool_comm) ibndmin = nint (tmp) CALL mp_min(ebndmin,inter_pool_comm) ! tmp = float (ibndmax) CALL mp_max(tmp, inter_pool_comm) ibndmax = nint (tmp) CALL mp_max(ebndmax,inter_pool_comm) #endif ! ndos = int( (ebndmax-ebndmin)/(degaussw/8.d0)) ALLOCATE (fdos(ndos) ) ALLOCATE (focc(ndos) ) fdos = zero focc = zero DO idos = 1, ndos en = ebndmin + (idos-1) * (degaussw/8.d0) focc(idos) = wgauss( -(en-ef)/degaussw, -99) focc(idos) = 0.d0 IF ((en-ef).lt. 0.02) focc(idos) = 1.d0 DO ik = 1, nksqf ikk = 2 * ik - 1 ! IF (epf_mem) THEN etf (ibndmin:ibndmax, ikk) = etfq (ibndmin:ibndmax, ikk, 1) ELSE nrec = ikk CALL davcio ( etf (ibndmin:ibndmax, ikk), ibndmax-ibndmin+1, iuetf, nrec, - 1) ENDIF ! DO ibnd = 1, nbndsub ! fdos(idos) = fdos(idos) + wkf(ikk) * & w0gauss( (etf(ibnd, ikk)-en) / degaussw, 0) / degaussw ENDDO ENDDO ENDDO #ifdef __PARA CALL mp_sum(fdos,inter_pool_comm) IF (ionode) THEN #endif iudosfil = find_free_unit() OPEN (unit = iudosfil, file = trim(prefix)//".dos", form = 'formatted') WRITE (iudosfil,*) "# Energy (eV) DOS" DO idos = 1, ndos WRITE (iudosfil,'(5x,f10.5,f14.8,f8.4)') (ebndmin + (idos-1) * (degaussw/8.d0))*ryd2ev, fdos(idos), focc(idos) ENDDO CLOSE (iudosfil) #ifdef __PARA ENDIF #endif WRITE (6,'(5x,a)') "Electronic DOS calculated" !------------------------- END SUBROUTINE epw_dos !-------------------------