Created documentation for the EPA code, renamed some internal variables to improve readability.

This commit is contained in:
Georgy Samsonidze 2020-01-27 08:39:52 -08:00
parent b77613de7a
commit 32cd912137
3 changed files with 164 additions and 119 deletions

View File

@ -263,6 +263,14 @@ input_description -distribution {Quantum Espresso} -package PWscf -program ph.x
from the electron-phonon interactions
using the optimized tetrahedron method.
}
opt -val 'epa' {
Electron-phonon coupling matrix elements are written
to file prefix.epa.k for further processing by program
epa.x which implements electron-phonon averaged (EPA)
approximation as described in G. Samsonidze & B. Kozinsky,
Adv. Energy Mater. 2018, 1800246 doi:10.1002/aenm.201800246
arXiv:1511.08115
}
info {
For metals only, requires gaussian smearing.

View File

@ -1143,13 +1143,9 @@ END SUBROUTINE elphsum_simple
SUBROUTINE elphfil_epa(iq)
!-----------------------------------------------------------------------
!
! EPA (electron-phonon-averaged) approximation
! Writes electron-phonon matrix elements to a file
! Written by Georgy Samsonidze on 2015-01-28
!
! Adv. Energy Mater. 2018, 1800246
! doi:10.1002/aenm.201800246
! https://doi.org/10.1002/aenm.201800246
! which is subsequently processed by the epa code
! Original routine written by Georgy Samsonidze
!
!-----------------------------------------------------------------------
USE cell_base, ONLY : ibrav, alat, omega, tpiba, at, bg

View File

@ -1,22 +1,64 @@
!-----------------------------------------------------------------------
!
! EPA (electron-phonon-averaged) approximation
! Averages electron-phonon matrix elements over k and q wavevectors
! in the bins of incident and scattered electron energy grids
! Written by Georgy Samsonidze from 2015-01-28 to 2015-05-11
! Copyright (C) 2018 Quantum ESPRESSO group
! 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 .
!
! Adv. Energy Mater. 2018, 1800246
! doi:10.1002/aenm.201800246
! https://doi.org/10.1002/aenm.201800246
!
!-----------------------------------------------------------------------
!----------------------------------------------------------------------------
program epa
!----------------------------------------------------------------------------
!
!! epa.x:
!! reads electron-phonon coupling matrix elements produced by the
!! phonon code with electron_phonon = 'epa', and makes a transformation
!! from momentum to energy space according to electron-phonon averaged
!! (EPA) approximation as described in G. Samsonidze & B. Kozinsky,
!! Adv. Energy Mater. 2018, 1800246 doi:10.1002/aenm.201800246
!! arXiv:1511.08115
!!
!! Input data:
!!
!! fin
!! fout
!! job
!! edgev, stepv, nbinv, ivmin, ivmax
!! edgec, stepc, nbinc, icmin, icmax
!! stepg, nbing, ngmax
!!
!! fin : input file produced by the phonon code
!! fout : output file used by the BoltzTraP code
!! job : job type determines how to average the matrix elements
!! - 'egrid' is to average within bins of regular energy
!! grids (standard procedure)
!! - 'bpair' is to average over individual pairs of bands
!! (obsolete procedure)
!! - 'ghist' is to construct the histogram for plotting
!! a distribution of squared matrix elements
!! edgev : grid edge of the valence energy grid (in eV), set to
!! the valence band maximum
!! stepv : grid step of the valence energy grid (in eV), set to
!! a small negative value
!! nbinv : number of bins in the valence energy grid
!! ivmin/max : range of valence bands (0/0 for all valence bands)
!! edgec : grid edge of the conduction energy grid (in eV), set to
!! the conduction band minimum
!! stepc : grid step of the conduction energy grid (in eV), set to
!! a small positive value
!! nbinc : number of bins in the conduction energy grid
!! icmin/max : range of conduction bands (0/0 for all conduction bands)
!! stepg : grid step for the histogram (in eV^2)
!! nbing : number of bins for the histogram
!! ngmax : maximum number of elements in one bin of the histogram
!!
!! For more details on how to set up the energy grids please refer to
!! online documentation at https://github.com/mir-group/EPA
!
use kinds, only : dp
use constants, only : ry_to_cmm1, rytoev, degspin, eps12
implicit none
integer, parameter :: nwin = 2
integer, parameter :: ngrid = 2
integer, parameter :: uni = 101
integer, parameter :: uno = 102
real(dp), parameter :: epsw2 = (20_dp / ry_to_cmm1)**2
@ -25,17 +67,17 @@ program epa
integer :: nmodes, nqs, nspin, nbnd, nkstot, nksqtot, &
nat, nsymq, irotmq, iq, ik, ikk, ikq, ibnd, jbnd, &
nu, mu, vu, ipert, jpert, ii, jj, kk, ll, ijob, &
nev, bvmin, bvmax, nec, bcmin, bcmax, ndist, nepmax, &
nemin, nemax, ne(nwin), nepair(nwin), bpair(2, nwin), &
nbinv, ivmin, ivmax, nbinc, icmin, icmax, nbing, ngmax, &
nbinmax, nbin(ngrid), nhist(ngrid), bpair(2, ngrid), &
s(3, 3, 48), invs(48)
real(dp) :: wtot, weight, factor, gbuf, gsum, &
gmax, gstep, ev, dev, ec, dec, wspin, &
ee(nwin), de(nwin), xq(3), at(3, 3), bg(3, 3)
character(len=256) :: fni, fno, job, fmt
stepg, edgev, stepv, edgec, stepc, wspin, &
edge(ngrid), step(ngrid), xq(3), at(3, 3), bg(3, 3)
character(len=256) :: fin, fout, job, fmt
character(len=32) :: s1, s2, s3, s4, s5
real(dp), allocatable :: ghist(:,:,:)
real(dp), allocatable :: gdist(:,:,:)
real(dp), allocatable :: gepair(:,:,:)
real(dp), allocatable :: gavg(:,:,:,:)
real(dp), allocatable :: wavg(:)
real(dp), allocatable :: gtot(:,:,:)
@ -56,53 +98,53 @@ program epa
complex(dp), allocatable :: el_ph_sum(:,:)
write(6, '("Reading standard input")')
read(5, '(a)') fni
read(5, '(a)') fno
read(5, '(a)') fin
read(5, '(a)') fout
read(5, '(a)') job
read(5, *) ev, dev, nev, bvmin, bvmax
read(5, *) ec, dec, nec, bcmin, bcmax
read(5, *) gmax, ndist, nepmax
read(5, *) edgev, stepv, nbinv, ivmin, ivmax
read(5, *) edgec, stepc, nbinc, icmin, icmax
read(5, *) stepg, nbing, ngmax
write(6, '(" input file name: ", a)') trim(fni)
write(6, '(" output file name: ", a)') trim(fno)
write(6, '(" input file name: ", a)') trim(fin)
write(6, '(" output file name: ", a)') trim(fout)
write(6, '(" job type: ", a)') trim(job)
write(s1, '(f16.8)') ev
write(s2, '(f16.8)') dev
write(s3, '(i8)') nev
write(s4, '(i8)') bvmin
write(s5, '(i8)') bvmax
write(6, '(" ev = ", a, " eV dev = ", a, " eV nev = ", &
a, " bvmin = ", a, " bvmax = ", a)') trim(adjustl(s1)), &
write(s1, '(f16.8)') edgev
write(s2, '(f16.8)') stepv
write(s3, '(i8)') nbinv
write(s4, '(i8)') ivmin
write(s5, '(i8)') ivmax
write(6, '(" edgev = ", a, " eV stepv = ", a, " eV nbinv = ", &
a, " ivmin = ", a, " ivmax = ", a)') trim(adjustl(s1)), &
trim(adjustl(s2)), trim(adjustl(s3)), trim(adjustl(s4)), &
trim(adjustl(s5))
write(s1, '(f16.8)') ec
write(s2, '(f16.8)') dec
write(s3, '(i8)') nec
write(s4, '(i8)') bcmin
write(s5, '(i8)') bcmax
write(6, '(" ec = ", a, " eV dec = ", a, " eV nec = ", &
a, " bcmin = ", a, " bcmax = ", a)') trim(adjustl(s1)), &
write(s1, '(f16.8)') edgec
write(s2, '(f16.8)') stepc
write(s3, '(i8)') nbinc
write(s4, '(i8)') icmin
write(s5, '(i8)') icmax
write(6, '(" edgec = ", a, " eV stepc = ", a, " eV nbinc = ", &
a, " icmin = ", a, " icmax = ", a)') trim(adjustl(s1)), &
trim(adjustl(s2)), trim(adjustl(s3)), trim(adjustl(s4)), &
trim(adjustl(s5))
write(s1, '(f16.8)') gmax
write(s2, '(i8)') ndist
write(s3, '(i8)') nepmax
write(6, '(" gmax = ", a, " eV^2 ndist = ", a, " nepmax = ", &
write(s1, '(f16.8)') stepg
write(s2, '(i8)') nbing
write(s3, '(i8)') ngmax
write(6, '(" stepg = ", a, " eV^2 nbing = ", a, " ngmax = ", &
a)') trim(adjustl(s1)), trim(adjustl(s2)), trim(adjustl(s3))
if (trim(job) .eq. 'bpair') then
ijob = 1
elseif (trim(job) .eq. 'egrid') then
ijob = 2
elseif (trim(job) .eq. 'gdist') then
elseif (trim(job) .eq. 'ghist') then
ijob = 3
else
write(0, '("Error: wrong job type")')
stop 1
endif
open(uni, file = fni, form = 'unformatted', status = 'old')
write(6, '("Reading file ", a)') trim(fni)
open(uni, file = fin, form = 'unformatted', status = 'old')
write(6, '("Reading file ", a)') trim(fin)
read(uni) s1, s2, s3
write(6, '(" title: ", a, ", date: ", a, ", time: ", a)') &
trim(s1), trim(s2), trim(s3)
@ -122,13 +164,12 @@ program epa
wspin = 1.0d0
endif
ee = (/ev, ec/)
de = (/dev, dec/)
ne = (/nev, nec/)
nemin = minval(ne)
nemax = maxval(ne)
bpair = reshape((/bvmin, bvmax, bcmin, bcmax/), shape(bpair))
do ii = 1, nwin
edge = (/edgev, edgec/)
step = (/stepv, stepc/)
nbin = (/nbinv, nbinc/)
nbinmax = maxval(nbin)
bpair = reshape((/ivmin, ivmax, icmin, icmax/), shape(bpair))
do ii = 1, ngrid
if (bpair(1, ii) .lt. 1 .or. bpair(1, ii) .gt. nbnd) bpair(1, ii) = 1
if (bpair(2, ii) .lt. 1 .or. bpair(2, ii) .gt. nbnd) bpair(2, ii) = nbnd
enddo
@ -136,28 +177,28 @@ program epa
allocate(gavg(nmodes, nbnd, nbnd, 1))
allocate(gtot(1, 1, 1))
elseif (ijob .eq. 2) then
allocate(gavg(nmodes, nemax, nemax, nwin))
allocate(gtot(nemax, nemax, nwin))
allocate(gnum(nemax, nemax, nwin))
allocate(egrid(2, nwin, nemax))
allocate(gavg(nmodes, nbinmax, nbinmax, ngrid))
allocate(gtot(nbinmax, nbinmax, ngrid))
allocate(gnum(nbinmax, nbinmax, ngrid))
allocate(egrid(2, ngrid, nbinmax))
elseif (ijob .eq. 3) then
allocate(gdist(2, ndist, nwin))
allocate(gepair(4, nepmax, nwin))
allocate(egrid(2, nwin, 1))
allocate(ghist(2, nbing, ngrid))
allocate(gdist(4, ngmax, ngrid))
allocate(egrid(2, ngrid, 1))
endif
if (ijob .eq. 1 .or. ijob .eq. 2) then
allocate(wavg(nmodes))
endif
if (ijob .eq. 2 .or. ijob .eq. 3) then
do ii = 1, nwin
egrid(ii, ii, 1) = ee(ii) + de(ii)
egrid(nwin + 1 - ii, ii, 1) = sum(ee) / nwin
do ii = 1, ngrid
egrid(ii, ii, 1) = edge(ii) + step(ii)
egrid(ngrid + 1 - ii, ii, 1) = sum(edge) / ngrid
enddo
if (ijob .eq. 2 .and. nemax .gt. 1) then
do jj = 2, nemax
do ii = 1, nwin
egrid(ii, ii, jj) = ee(ii) + de(ii) * jj
egrid(nwin + 1 - ii, ii, jj) = ee(ii) + de(ii) * (jj - 1)
if (ijob .eq. 2 .and. nbinmax .gt. 1) then
do jj = 2, nbinmax
do ii = 1, ngrid
egrid(ii, ii, jj) = edge(ii) + step(ii) * jj
egrid(ngrid + 1 - ii, ii, jj) = edge(ii) + step(ii) * (jj - 1)
enddo
enddo
endif
@ -194,14 +235,13 @@ program epa
gtot(:,:,:) = 0.0d0
wtot = 0.0d0
elseif (ijob .eq. 3) then
gstep = gmax / ndist
do ii = 1, nwin
do jj = 1, ndist
gdist(1, jj, ii) = (jj - 1) * gstep
gdist(2, jj, ii) = 0.0d0
do ii = 1, ngrid
do jj = 1, nbing
ghist(1, jj, ii) = (jj - 1) * stepg
ghist(2, jj, ii) = 0.0d0
enddo
enddo
nepair(:) = 0
nhist(:) = 0
endif
if (ijob .eq. 2) gnum(:,:,:) = 0
do iq = 1, nqs
@ -288,9 +328,9 @@ program epa
enddo
enddo
elseif (ijob .eq. 2) then
do ii = 1, nwin
do jj = 1, nemax
do kk = 1, nemax
do ii = 1, ngrid
do jj = 1, nbin(ii)
do kk = 1, nbin(ii)
el_ph_sum(:,:) = (0.0d0, 0.0d0)
do ik = 1, nksqtot
ikk = ikks(ik)
@ -337,7 +377,7 @@ program epa
enddo
enddo
elseif (ijob .eq. 3) then
do ii = 1, nwin
do ii = 1, ngrid
do ik = 1, nksqtot
ikk = ikks(ik)
ikq = ikqs(ik)
@ -373,15 +413,15 @@ program epa
enddo
gsum = gsum + gbuf * factor
enddo
jj = nint(abs(gsum) / gstep)
if (jj .ge. 1 .and. jj .le. ndist) then
gdist(2, jj, ii) = gdist(2, jj, ii) + weight / gstep
jj = nint(abs(gsum) / stepg)
if (jj .ge. 1 .and. jj .le. nbing) then
ghist(2, jj, ii) = ghist(2, jj, ii) + weight / stepg
endif
nepair(ii) = nepair(ii) + 1
gepair(1, nepair(ii), ii) = et(ibnd, ikk) * rytoev - ee(ii)
gepair(2, nepair(ii), ii) = et(jbnd, ikq) * rytoev - ee(ii)
gepair(3, nepair(ii), ii) = gsum
gepair(4, nepair(ii), ii) = weight
nhist(ii) = nhist(ii) + 1
gdist(1, nhist(ii), ii) = et(ibnd, ikk) * rytoev - edge(ii)
gdist(2, nhist(ii), ii) = et(jbnd, ikq) * rytoev - edge(ii)
gdist(3, nhist(ii), ii) = gsum
gdist(4, nhist(ii), ii) = weight
endif
enddo
endif
@ -417,9 +457,9 @@ program epa
enddo
enddo
elseif (ijob .eq. 2) then
do ii = 1, nwin
do jj = 1, nemax
do kk = 1, nemax
do ii = 1, ngrid
do jj = 1, nbin(ii)
do kk = 1, nbin(ii)
if (gtot(kk, jj, ii) .gt. eps12) then
do nu = 1, nmodes
gavg(nu, kk, jj, ii) = gavg(nu, kk, jj, ii) / gtot(kk, jj, ii)
@ -439,20 +479,20 @@ program epa
if (ijob .eq. 2) then
write(fmt, '("(", i0, "(5x, ""e"", a, "" (eV)"", 4x, ""e"", a, ""''", &
" (eV)"", 3x, ""<|g"", a, ""|^2> (eV^2)"", 2x, ""count"", a, 5x,", &
" ""weight"", a))")') nwin
" ""weight"", a))")') ngrid
write(6, fmt) ('v', ll = 1, 5), ('c', ll = 1, 5)
write(fmt, '("(", i0, "(2f12.6, e18.8, i8, f12.6))")') nwin
do jj = 1, nemin
do kk = 1, nemin
write(6, fmt) ((jj - 0.5d0) * de(ii), (kk - 0.5d0) * de(ii), &
write(fmt, '("(", i0, "(2f12.6, e18.8, i8, f12.6))")') ngrid
do jj = 1, nbinmax
do kk = 1, nbinmax
write(6, fmt) ((jj - 0.5d0) * step(ii), (kk - 0.5d0) * step(ii), &
sum(gavg(:, kk, jj, ii)), gnum(kk, jj, ii), gtot(kk, jj, ii), &
ii = 1, nwin)
ii = 1, ngrid)
enddo
enddo
endif
open(uno, file = fno, form = 'formatted', status = 'replace')
write(6, '("Writing file ", a)') trim(fno)
open(uno, file = fout, form = 'formatted', status = 'replace')
write(6, '("Writing file ", a)') trim(fout)
if (ijob .eq. 1) then
write(uno, '(" mode <w> (cm^-1)")')
write(uno, '("--------------------------")')
@ -470,31 +510,31 @@ program epa
enddo
enddo
elseif (ijob .eq. 2) then
write(uno, '(2i8)') nwin, nmodes
do ii = 1, nwin
write(uno, '(2f14.8, i8)') ee(ii), de(ii), ne(ii)
write(uno, '(2i8)') ngrid, nmodes
do ii = 1, ngrid
write(uno, '(2f14.8, i8)') edge(ii), step(ii), nbin(ii)
enddo
write(fmt, '("(", i0, "e18.8)")') nmodes
write(uno, fmt) (wavg(nu), nu = 1, nmodes)
write(fmt, '("(3i8, ", i0, "e18.8)")') nmodes
do ii = 1, nwin
do jj = 1, ne(ii)
do kk = 1, ne(ii)
do ii = 1, ngrid
do jj = 1, nbin(ii)
do kk = 1, nbin(ii)
write(uno, fmt) ii, jj, kk, (gavg(nu, kk, jj, ii), &
nu = 1, nmodes)
enddo
enddo
enddo
elseif (ijob .eq. 3) then
write(uno, '(i8)') ndist
do jj = 1, ndist
write(uno, '(f14.8, 2e18.8)') gdist(1, jj, 1), gdist(2, jj, 1), &
gdist(2, jj, 2)
write(uno, '(i8)') nbing
do jj = 1, nbing
write(uno, '(f14.8, 2e18.8)') ghist(1, jj, 1), ghist(2, jj, 1), &
ghist(2, jj, 2)
enddo
write(uno, '(i8)') (nepair(ii), ii = 1, nwin)
do ii = 1, nwin
do jj = 1, nepair(ii)
write(uno, '(2f14.8, 2e18.8)') (gepair(kk, jj, ii), kk = 1, 4)
write(uno, '(i8)') (nhist(ii), ii = 1, ngrid)
do ii = 1, ngrid
do jj = 1, nhist(ii)
write(uno, '(2f14.8, 2e18.8)') (gdist(kk, jj, ii), kk = 1, 4)
enddo
enddo
endif
@ -507,8 +547,8 @@ program epa
endif
if (ijob .eq. 2) deallocate(gnum)
if (ijob .eq. 3) then
deallocate(ghist)
deallocate(gdist)
deallocate(gepair)
endif
if (ijob .eq. 2 .or. ijob .eq. 3) deallocate(egrid)
deallocate(x_q)
@ -521,4 +561,5 @@ program epa
deallocate(el_ph_sum)
end program epa
!
!----------------------------------------------------------------------------