mirror of https://gitlab.com/QEF/q-e.git
cumulative module updated
This commit is contained in:
zx199323
parent
5589beee2e
commit
1f2f97aacf
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Reference in New Issue