From 1f2f97aacff69ceee3e800bcbf7db57d5c64b038 Mon Sep 17 00:00:00 2001 From: zx199323 Date: Fri, 3 Jul 2020 16:37:44 -0400 Subject: [PATCH] cumulative module updated --- EPW/src/cum_mod.f90 | 305 +++++++++++++++++++++++--------------------- 1 file changed, 157 insertions(+), 148 deletions(-) diff --git a/EPW/src/cum_mod.f90 b/EPW/src/cum_mod.f90 index ffc86acda..81225c053 100644 --- a/EPW/src/cum_mod.f90 +++ b/EPW/src/cum_mod.f90 @@ -35,16 +35,18 @@ !! !----------------------------------------------------------------------- USE kinds, ONLY : DP, i4b - USE constants_epw, ONLY : two, zero, ryd2ev, ryd2mev, ci + USE constants_epw, ONLY : kelvin2eV, two, zero, ryd2ev, ryd2mev, ci USE constants, ONLY : pi USE io_global, ONLY : stdout USE io_var, ONLY : iospectral_sup, iospectral_cum USE epwcom, ONLY : eptemp, wmin_specfun, wmax_specfun, nw_specfun, & - bnd_cum - USE elph2, ONLY : ibndmin, ibndmax + bnd_cum, nstemp + USE elph2, ONLY : ibndmin, ibndmax, gtemp ! IMPLICIT NONE ! + CHARACTER(LEN = 20) :: tp + CHARACTER(LEN = 256) :: filespecsup CHARACTER(LEN = 64) :: line !! Auxiliary string CHARACTER(LEN = 64) :: filespec @@ -67,6 +69,8 @@ !! Total number of k-points INTEGER :: i0 !! Energy index of Fermi level (w=0) + INTEGER :: itemp + !! Counter on temperatures INTEGER :: ierr !! Error status REAL(KIND = DP) :: dw @@ -107,154 +111,159 @@ WRITE(stdout, '(5x,a)') 'Warning: the routine is sequential but very fast.' WRITE(stdout, '(5x,a/)') REPEAT('=',75) ! - OPEN (UNIT = iospectral_sup, FILE = 'specfun_sup.elself', STATUS = 'old', IOSTAT = ios) - IF (ios /= 0) CALL errore('spectral_cumulant', 'opening file specfun_sup.elself', ABS(ios)) - ! - ! determine number of k points, ibndmin, ibndmax - DO im = 1, 6 - READ(iospectral_sup, '(a)') line - ENDDO - DO im = 1, maxrecs - READ (iospectral_sup, *, IOSTAT = ios) i1, i2 - IF (im == 1) ibndmin = i2 - IF (ios /= 0) EXIT - IF (im == maxrecs) CALL errore('spectral_cumulant', 'increase maxrecs', 1) - ENDDO - ! - REWIND(iospectral_sup) - ! - nk = i1 - ibndmax = i2 - WRITE(stdout, '(5x,a/)') "Read self-energy from file specfun_sup.elself" - WRITE(stdout, '(5x,a,i4,a,i4,a,i4,a,f12.6/)') "Check: nk = ", nk, & - ", ibndmin = ", ibndmin, ", ibndmax = ", ibndmax, " kbT (eV) = ", eptemp * ryd2ev - ! - ALLOCATE(ww(nw_specfun), STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating ww', 1) - ALLOCATE(ek(nk), STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating ek', 1) - ALLOCATE(sigmar(nk, nw_specfun), STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating sigmar', 1) - ALLOCATE(sigmai(nk, nw_specfun), STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating sigmai', 1) - ALLOCATE(a_mig(nw_specfun, nk), STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_mig', 1) - ALLOCATE(a_cw(nw_specfun), STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_cw', 1) - ALLOCATE(a_ct(nw_specfun), STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_ct', 1) - ALLOCATE(a_tmp(nw_specfun), STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_tmp', 1) - ! - ! read and store Kohn-Sham energy, energy grid, real and im sigma for designated band - DO im = 1,6 - READ(iospectral_sup, '(a)') line - ENDDO - DO ibnd = 1, ibndmax - ibndmin + 1 - DO ik = 1, nk - DO iw = 1, nw_specfun - READ(iospectral_sup,*) i1, i2, a1, a2, a3, a4 - IF (i2 == bnd_cum) THEN - ! ek, w read in eV; Sigma read in meV - ek(ik) = a1 / ryd2ev - ww(iw) = a2 / ryd2ev - sigmar(ik, iw) = a3 / ryd2mev ! / ( EXP(ww(iw)/eptemp )+1.d0 ) - sigmai(ik, iw) = a4 / ryd2mev ! / ( EXP(ww(iw)/eptemp )+1.d0 ) - ! spec func as in spectral_func.f90 - a_mig(iw, ik) = ABS(sigmai(ik, iw)) / pi / ((ww(iw) - ek(ik) - sigmar(ik, iw))**two + (sigmai(ik, iw) )**two) - ENDIF + DO itemp = 1, nstemp + WRITE(tp, "(f8.3)") gtemp(itemp) * ryd2ev / kelvin2eV + filespecsup = 'specfun_sup.elself.' // trim(adjustl(tp)) // 'K' + OPEN (UNIT = iospectral_sup, FILE = filespecsup, STATUS = 'old', IOSTAT = ios) + IF (ios /= 0) CALL errore('spectral_cumulant', 'opening file specfun_sup.elself', ABS(ios)) + ! + ! determine number of k points, ibndmin, ibndmax + DO im = 1, 6 + READ(iospectral_sup, '(a)') line + ENDDO + DO im = 1, maxrecs + READ (iospectral_sup, *, IOSTAT = ios) i1, i2 + IF (im == 1) ibndmin = i2 + IF (ios /= 0) EXIT + IF (im == maxrecs) CALL errore('spectral_cumulant', 'increase maxrecs', 1) + ENDDO + ! + REWIND(iospectral_sup) + ! + nk = i1 + ibndmax = i2 + WRITE(stdout, '(5x,a/)') "Read self-energy from file specfun_sup.elself" + WRITE(stdout, '(5x,a,i4,a,i4,a,i4,a,f12.6/)') "Check: nk = ", nk, & + ", ibndmin = ", ibndmin, ", ibndmax = ", ibndmax, " kbT (eV) = ", gtemp(itemp) * ryd2ev + ! + ALLOCATE(ww(nw_specfun), STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating ww', 1) + ALLOCATE(ek(nk), STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating ek', 1) + ALLOCATE(sigmar(nk, nw_specfun), STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating sigmar', 1) + ALLOCATE(sigmai(nk, nw_specfun), STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating sigmai', 1) + ALLOCATE(a_mig(nw_specfun, nk), STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_mig', 1) + ALLOCATE(a_cw(nw_specfun), STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_cw', 1) + ALLOCATE(a_ct(nw_specfun), STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_ct', 1) + ALLOCATE(a_tmp(nw_specfun), STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_tmp', 1) + ! + ! read and store Kohn-Sham energy, energy grid, real and im sigma for designated band + DO im = 1,6 + READ(iospectral_sup, '(a)') line + ENDDO + DO ibnd = 1, ibndmax - ibndmin + 1 + DO ik = 1, nk + DO iw = 1, nw_specfun + READ(iospectral_sup,*) i1, i2, a1, a2, a3, a4 + IF (i2 == bnd_cum) THEN + ! ek, w read in eV; Sigma read in meV + ek(ik) = a1 / ryd2ev + ww(iw) = a2 / ryd2ev + sigmar(ik, iw) = a3 / ryd2mev ! / ( EXP(ww(iw)/eptemp )+1.d0 ) + sigmai(ik, iw) = a4 / ryd2mev ! / ( EXP(ww(iw)/eptemp )+1.d0 ) + ! spec func as in spectral_func.f90 + a_mig(iw, ik) = ABS(sigmai(ik, iw)) / pi / ((ww(iw) - ek(ik) - sigmar(ik, iw))**two + (sigmai(ik, iw) )**two) + ENDIF + ENDDO ENDDO ENDDO - ENDDO - ! - CLOSE(iospectral_sup) - ! - ! open file for cumulant spectral function - IF (bnd_cum < 10) THEN - WRITE(filespec, '(a,i1,a)') 'specfun_cum', bnd_cum, '.elself' - ELSEIF (bnd_cum > 9 .AND. bnd_cum < 100) THEN - WRITE(filespec, '(a,i2,a)') 'specfun_cum', bnd_cum, '.elself' - ELSE - WRITE(filespec, '(a,i3,a)') 'specfun_cum', bnd_cum, '.elself' - ENDIF - OPEN(UNIT = iospectral_cum, FILE = filespec) - ! - WRITE(iospectral_cum, '(a)') '# k Energy [eV] A(k,w) [meV^-1] Z-factor ' - WRITE(iospectral_cum, '(a)') '# with convolutions | using FFT ' - WRITE(stdout, '(8x,a)') 'k Energy [eV] A(k,w) [meV^-1] Z-factor ' - WRITE(stdout, '(8x,a)') ' with convolutions | using FFT ' - ! - ! define index corresponding to omega=0 (Fermi level) - i0 = MINLOC(ABS(ww(:)), DIM = 1) - IF (ABS(ww(i0)) > dw) CALL errore('spectral_cumulant', 'w=0 needs to be included in [wmin:wmax]', 1) - a_cw = zero - a_ct = zero - a_tmp = zero - ! - DO ik = 1, nk - IF (ek(ik) < e_thresh) THEN - ! - ekk = ek(ik) - ! - ! Cumulant from convolutions in frequency space (ImSigma>0 in EPW) - CALL cumulant_conv(ekk, ww, sigmar(ik, :), -sigmai(ik, :), a_mig(:, ik), a_cw, zeta) - ! - ! Cumulant calculated in time domain + FFT (ImSigma>0 in EPW) - CALL cumulant_time(ekk, ww, sigmar(ik, :), -sigmai(ik, :), a_mig(:, ik), a_tmp) - ! - DO iw = 1, nw_specfun - ! - ! map the indices of the FFT frequency grid onto the original one - IF (iw >= i0) THEN - a_ct(iw) = a_tmp(iw - i0 + 1) - ELSE - a_ct(iw) = a_tmp(iw + nw_specfun - i0 + 1) - ENDIF - ! - ! write cumulant spectral function on file (in meV^-1, as in spectral_func.f90) - ! 3rd column: A_cum using convolutions; 4th column: A_cum using FFT - IF (iw == 1) THEN - WRITE(iospectral_cum, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,f8.4)') & - ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev, zeta - WRITE(stdout,'(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,f8.4)') & - ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev, zeta - ELSE - WRITE(iospectral_cum, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7)') & - ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev !/ ( EXP(ww(iw)/eptemp )+1.d0 ) - WRITE(stdout, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7)') & - ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev !/ ( EXP(ww(iw)/eptemp )+1.d0 ) - ! uncomment to multiply by Fermi occupation factor - ENDIF - ! - ENDDO - ! - WRITE(iospectral_cum, '(a)') ' ' - WRITE(stdout, '(a)') ' ' - ! - ENDIF ! only states below energy e_thresh ! - ENDDO ! main loop k - ! - WRITE(stdout, '(5x,a)') 'The file specfun_cum[BND].elself has been correctly written' - ! - CLOSE(iospectral_cum) - ! - DEALLOCATE(ww, STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating ww', 1) - DEALLOCATE(ek, STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating ek', 1) - DEALLOCATE(sigmar, STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating sigmar', 1) - DEALLOCATE(sigmai, STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating sigmai', 1) - DEALLOCATE(a_mig, STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_mig', 1) - DEALLOCATE(a_cw, STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_cw', 1) - DEALLOCATE(a_ct, STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_ct', 1) - DEALLOCATE(a_tmp, STAT = ierr) - IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_tmp', 1) + CLOSE(iospectral_sup) + ! + ! open file for cumulant spectral function + WRITE(tp, "(f8.3)") gtemp(itemp) * ryd2ev / kelvin2eV + IF (bnd_cum < 10) THEN + WRITE(filespec, '(a,i1,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K' + ELSEIF (bnd_cum > 9 .AND. bnd_cum < 100) THEN + WRITE(filespec, '(a,i2,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K' + ELSE + WRITE(filespec, '(a,i3,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K' + ENDIF + OPEN(UNIT = iospectral_cum, FILE = filespec) + ! + WRITE(iospectral_cum, '(a)') '# k Energy [eV] A(k,w) [meV^-1] Z-factor ' + WRITE(iospectral_cum, '(a)') '# with convolutions | using FFT ' + WRITE(stdout, '(8x,a)') 'k Energy [eV] A(k,w) [meV^-1] Z-factor ' + WRITE(stdout, '(8x,a)') ' with convolutions | using FFT ' + ! + ! define index corresponding to omega=0 (Fermi level) + i0 = MINLOC(ABS(ww(:)), DIM = 1) + IF (ABS(ww(i0)) > dw) CALL errore('spectral_cumulant', 'w=0 needs to be included in [wmin:wmax]', 1) + a_cw = zero + a_ct = zero + a_tmp = zero + ! + DO ik = 1, nk + IF (ek(ik) < e_thresh) THEN + ! + ekk = ek(ik) + ! + ! Cumulant from convolutions in frequency space (ImSigma>0 in EPW) + CALL cumulant_conv(ekk, ww, sigmar(ik, :), -sigmai(ik, :), a_mig(:, ik), a_cw, zeta) + ! + ! Cumulant calculated in time domain + FFT (ImSigma>0 in EPW) + CALL cumulant_time(ekk, ww, sigmar(ik, :), -sigmai(ik, :), a_mig(:, ik), a_tmp) + ! + DO iw = 1, nw_specfun + ! + ! map the indices of the FFT frequency grid onto the original one + IF (iw >= i0) THEN + a_ct(iw) = a_tmp(iw - i0 + 1) + ELSE + a_ct(iw) = a_tmp(iw + nw_specfun - i0 + 1) + ENDIF + ! + ! write cumulant spectral function on file (in meV^-1, as in spectral_func.f90) + ! 3rd column: A_cum using convolutions; 4th column: A_cum using FFT + IF (iw == 1) THEN + WRITE(iospectral_cum, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,f8.4)') & + ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev, zeta + WRITE(stdout,'(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,f8.4)') & + ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev, zeta + ELSE + WRITE(iospectral_cum, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7)') & + ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev !/ ( EXP(ww(iw)/eptemp )+1.d0 ) + WRITE(stdout, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7)') & + ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev !/ ( EXP(ww(iw)/eptemp )+1.d0 ) + ! uncomment to multiply by Fermi occupation factor + ENDIF + ! + ENDDO + ! + WRITE(iospectral_cum, '(a)') ' ' + WRITE(stdout, '(a)') ' ' + ! + ENDIF ! only states below energy e_thresh + ! + ENDDO ! main loop k + ! + WRITE(stdout, '(5x,a)') 'The file specfun_cum[BND].elself has been correctly written' + ! + CLOSE(iospectral_cum) + ! + DEALLOCATE(ww, STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating ww', 1) + DEALLOCATE(ek, STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating ek', 1) + DEALLOCATE(sigmar, STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating sigmar', 1) + DEALLOCATE(sigmai, STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating sigmai', 1) + DEALLOCATE(a_mig, STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_mig', 1) + DEALLOCATE(a_cw, STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_cw', 1) + DEALLOCATE(a_ct, STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_ct', 1) + DEALLOCATE(a_tmp, STAT = ierr) + IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_tmp', 1) + ENDDO !itemp ! !----------------------------------------------------------------------- END SUBROUTINE spectral_cumulant