EPA method.

This commit is contained in:
Samsonidze Georgy (CR/RTC2.2-NA) 2018-06-13 13:13:46 -07:00
parent dae010abfc
commit b77613de7a
8 changed files with 851 additions and 11 deletions

View File

@ -192,7 +192,7 @@ QEMODS = ../../Modules/libqemod.a ../../KS_Solvers/Davidson/libdavid.a ../../KS_
TLDEPS= bindir libs mods libdavid libcg lrmods
all : tldeps libs-ph ph.x dynmat.x matdyn.x q2r.x q2trans.x q2trans_fd.x lambda.x fqha.x q2qstar.x \
alpha2f.x
alpha2f.x epa.x
libs-ph : libph.a libphaux.a
@ -242,6 +242,11 @@ alpha2f.x : alpha2f.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
$(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS) $(LIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
epa.x : epa.o libph.a $(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS)
$(LD) $(LDFLAGS) -o $@ epa.o libph.a \
$(PWOBJS) $(LRMODS) $(QEMODS) $(LIBOBJS) $(LIBS)
- ( cd ../../bin ; ln -fs ../PHonon/PH/$@ . )
#fqha.o :
# $(MPIF90) $(FFLAGS_NOOPT) -c fqha.f90

View File

@ -72,7 +72,7 @@ SUBROUTINE check_initial_status(auxdyn)
USE io_files, ONLY : tmp_dir, postfix
USE lsda_mod, ONLY : nspin
USE scf, ONLY : rho
USE disp, ONLY : nqs, x_q, comp_iq, nq1, nq2, nq3, &
USE disp, ONLY : nqs, x_q, wq, comp_iq, nq1, nq2, nq3, &
done_iq, lgamma_iq
USE qpoint, ONLY : xq
USE control_lr, ONLY : lgamma
@ -138,8 +138,10 @@ SUBROUTINE check_initial_status(auxdyn)
nqs = 1
last_q = 1
ALLOCATE(x_q(3,1))
ALLOCATE(wq(1))
ALLOCATE(lgamma_iq(1))
x_q(:,1)=xq(:)
wq(1)=1.0d0
lgamma_iq(1)=lgamma
!
END IF

View File

@ -33,7 +33,7 @@ SUBROUTINE do_phonon(auxdyn)
USE disp, ONLY : nqs
USE control_ph, ONLY : epsil, trans, qplot, only_init, &
only_wfc, rec_code, where_rec
USE el_phon, ONLY : elph, elph_mat, elph_simple
USE el_phon, ONLY : elph, elph_mat, elph_simple, elph_epa
!
! YAMBO >
USE YAMBO, ONLY : elph_yambo
@ -116,6 +116,8 @@ SUBROUTINE do_phonon(auxdyn)
CALL elphsum_wannier(iq)
ELSEIF( elph_simple ) THEN
CALL elphsum_simple()
ELSEIF( elph_epa ) THEN
CALL elphfil_epa(iq)
ELSEIF( elph_yambo ) THEN
CALL elph_yambo_eval_and_IO()
ELSEIF(elph_tetra == 1) THEN

View File

@ -12,7 +12,7 @@ MODULE el_phon
!
SAVE
!
LOGICAL :: elph, elph_mat, elph_simple
LOGICAL :: elph, elph_mat, elph_simple, elph_epa
INTEGER :: elph_nbnd_min, elph_nbnd_max
INTEGER :: el_ph_ngauss, el_ph_nsigma
INTEGER :: iunwfcwann, lrwfcr

View File

@ -1139,6 +1139,300 @@ SUBROUTINE elphsum_simple
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
!
!-----------------------------------------------------------------------
USE cell_base, ONLY : ibrav, alat, omega, tpiba, at, bg
USE disp, ONLY : nq1, nq2, nq3, nqs, x_q, wq, lgamma_iq
USE dynmat, ONLY : dyn, w2
USE el_phon, ONLY : el_ph_mat, done_elph
USE fft_base, ONLY : dfftp, dffts, dfftb
USE gvect, ONLY : ngm_g, ecutrho
USE io_global, ONLY : ionode, ionode_id
USE ions_base, ONLY : nat, nsp, atm, ityp, tau
USE kinds, ONLY : DP
USE klist, ONLY : xk, wk, nelec, nks, nkstot, ngk
USE lsda_mod, ONLY : nspin, isk
USE modes, ONLY : nirr, nmodes, npert, npertx, u, t, tmq, &
name_rap_mode, num_rap_mode
USE lr_symm_base, ONLY : irgq, nsymq, irotmq, rtau, gi, gimq, &
minus_q, invsymq
USE mp, ONLY : mp_bcast, mp_sum
USE mp_images, ONLY : intra_image_comm
USE mp_pools, ONLY : npool, intra_pool_comm
USE qpoint, ONLY : nksq, nksqtot, ikks, ikqs, eigqts
USE start_k, ONLY : nk1, nk2, nk3, k1, k2, k3
USE symm_base, ONLY : s, invs, ftau, nrot, nsym, nsym_ns, &
nsym_na, ft, sr, sname, t_rev, irt, time_reversal, &
invsym, nofrac, allfrac, nosym, nosym_evc, no_t_rev
USE wvfct, ONLY : nbnd, et, wg
USE gvecw, ONLY : ecutwfc
USE io_files, ONLY : prefix
IMPLICIT NONE
INTEGER, INTENT(IN) :: iq
INTEGER :: iuelph, ios, irr, ii, jj, kk, ll
character :: cdate*9, ctime*9, sdate*32, stime*32, &
stitle*32, myaccess*10, mystatus*7
CHARACTER(LEN=80) :: filelph
REAL(DP), ALLOCATABLE :: xk_collect(:,:), wk_collect(:)
REAL(DP), ALLOCATABLE :: et_collect(:,:), wg_collect(:,:)
INTEGER, ALLOCATABLE :: ngk_collect(:)
INTEGER, ALLOCATABLE :: ikks_collect(:), ikqs_collect(:)
COMPLEX(DP), ALLOCATABLE :: el_ph_mat_collect(:,:,:,:)
INTEGER, EXTERNAL :: find_free_unit, atomic_number
filelph = TRIM(prefix) // '.epa.k'
DO irr = 1, nirr
IF (.NOT. done_elph(irr)) RETURN
ENDDO
IF (iq .EQ. 1) THEN
myaccess = 'sequential'
mystatus = 'replace'
ELSE
myaccess = 'append'
mystatus = 'old'
ENDIF
IF (ionode) THEN
iuelph = find_free_unit()
OPEN(unit = iuelph, file = TRIM(filelph), form = 'unformatted', &
access = myaccess, status = mystatus, iostat = ios)
ELSE
iuelph = 0
ENDIF
CALL mp_bcast(ios, ionode_id, intra_image_comm)
CALL errore('elphfil_epa', 'opening file ' // filelph, ABS(ios))
IF (iq .EQ. 1) THEN
CALL date_and_tim(cdate, ctime)
WRITE(sdate, '(A2,"-",A3,"-",A4,21X)') cdate(1:2), cdate(3:5), cdate(6:9)
WRITE(stime, '(A8,24X)') ctime(1:8)
WRITE(stitle, '("EPA-Complex",21X)')
CALL cryst_to_cart(nqs, x_q, at, -1)
! write header
IF (ionode) THEN
WRITE(iuelph) stitle, sdate, stime
WRITE(iuelph) ibrav, nat, nsp, nrot, nsym, nsym_ns, nsym_na, &
ngm_g, nspin, nbnd, nmodes, nqs
WRITE(iuelph) nq1, nq2, nq3, nk1, nk2, nk3, k1, k2, k3
WRITE(iuelph) time_reversal, invsym, nofrac, allfrac, nosym, &
nosym_evc, no_t_rev
WRITE(iuelph) alat, omega, tpiba, nelec, ecutrho, ecutwfc
WRITE(iuelph) dfftp%nr1, dfftp%nr2, dfftp%nr3
WRITE(iuelph) dffts%nr1, dffts%nr2, dffts%nr3
WRITE(iuelph) dfftb%nr1, dfftb%nr2, dfftb%nr3
WRITE(iuelph) ((at(ii, jj), ii = 1, 3), jj = 1, 3)
WRITE(iuelph) ((bg(ii, jj), ii = 1, 3), jj = 1, 3)
WRITE(iuelph) (atomic_number(atm(ii)), ii = 1, nsp)
WRITE(iuelph) (ityp(ii), ii = 1, nat)
WRITE(iuelph) ((tau(ii, jj), ii = 1, 3), jj = 1, nat)
WRITE(iuelph) ((x_q(ii, jj), ii = 1, 3), jj = 1, nqs)
WRITE(iuelph) (wq(ii), ii = 1, nqs)
WRITE(iuelph) (lgamma_iq(ii), ii = 1, nqs)
ENDIF
CALL cryst_to_cart(nqs, x_q, bg, 1)
ENDIF
! collect data for current q-point
ALLOCATE(xk_collect(3, nkstot))
ALLOCATE(wk_collect(nkstot))
ALLOCATE(et_collect(nbnd, nkstot))
ALLOCATE(wg_collect(nbnd, nkstot))
ALLOCATE(ngk_collect(nkstot))
ALLOCATE(ikks_collect(nksqtot))
ALLOCATE(ikqs_collect(nksqtot))
ALLOCATE(el_ph_mat_collect(nbnd, nbnd, nksqtot, nmodes))
IF (npool > 1) THEN
CALL poolcollect(3, nks, xk, nkstot, xk_collect)
CALL poolcollect(1, nks, wk, nkstot, wk_collect)
CALL poolcollect(nbnd, nks, et, nkstot, et_collect)
CALL poolcollect(nbnd, nks, wg, nkstot, wg_collect)
CALL ipoolcollect(1, nks, ngk, nkstot, ngk_collect)
CALL jpoolcollect(1, nksq, ikks, nksqtot, ikks_collect)
CALL jpoolcollect(1, nksq, ikqs, nksqtot, ikqs_collect)
CALL el_ph_collect(nmodes, el_ph_mat, el_ph_mat_collect, nksqtot, nksq)
ELSE
xk_collect(1:3, 1:nks) = xk(1:3, 1:nks)
wk_collect(1:nks) = wk(1:nks)
et_collect(1:nbnd, 1:nks) = et(1:nbnd, 1:nks)
wg_collect(1:nbnd, 1:nks) = wg(1:nbnd, 1:nks)
ngk_collect(1:nks) = ngk(1:nks)
ikks_collect(1:nksq) = ikks(1:nksq)
ikqs_collect(1:nksq) = ikqs(1:nksq)
el_ph_mat_collect(1:nbnd, 1:nbnd, 1:nksq, 1:nmodes) = &
el_ph_mat(1:nbnd, 1:nbnd, 1:nksq, 1:nmodes)
ENDIF
CALL cryst_to_cart(nkstot, xk_collect, at, -1)
! write data for current q-point
IF (ionode) THEN
WRITE(iuelph) nsymq, irotmq, nirr, npertx, nkstot, nksqtot
WRITE(iuelph) minus_q, invsymq
WRITE(iuelph) (irgq(ii), ii = 1, 48)
WRITE(iuelph) (npert(ii), ii = 1, nmodes)
WRITE(iuelph) (((rtau(ii, jj, kk), ii = 1, 3), jj = 1, 48), &
kk = 1, nat)
WRITE(iuelph) ((gi(ii, jj), ii = 1, 3), jj = 1, 48)
WRITE(iuelph) (gimq(ii), ii = 1, 3)
WRITE(iuelph) ((u(ii, jj), ii = 1, nmodes), jj = 1, nmodes)
WRITE(iuelph) ((((t(ii, jj, kk, ll), ii = 1, npertx), &
jj = 1, npertx), kk = 1, 48), ll = 1, nmodes)
WRITE(iuelph) (((tmq(ii, jj, kk), ii = 1, npertx), &
jj = 1, npertx), kk = 1, nmodes)
WRITE(iuelph) (name_rap_mode(ii), ii = 1, nmodes)
WRITE(iuelph) (num_rap_mode(ii), ii = 1, nmodes)
WRITE(iuelph) (((s(ii, jj, kk), ii = 1, 3), jj = 1, 3), kk = 1, 48)
WRITE(iuelph) (invs(ii), ii = 1, 48)
WRITE(iuelph) ((ftau(ii, jj), ii = 1, 3), jj = 1, 48)
WRITE(iuelph) ((ft(ii, jj), ii = 1, 3), jj = 1, 48)
WRITE(iuelph) (((sr(ii, jj, kk), ii = 1, 3), jj = 1, 3), kk = 1, 48)
WRITE(iuelph) (sname(ii), ii = 1, 48)
WRITE(iuelph) (t_rev(ii), ii = 1, 48)
WRITE(iuelph) ((irt(ii, jj), ii = 1, 48), jj = 1, nat)
WRITE(iuelph) ((xk_collect(ii, jj), ii = 1, 3), jj = 1, nkstot)
WRITE(iuelph) (wk_collect(ii), ii = 1, nkstot)
WRITE(iuelph) ((et_collect(ii, jj), ii = 1, nbnd), jj = 1, nkstot)
WRITE(iuelph) ((wg_collect(ii, jj), ii = 1, nbnd), jj = 1, nkstot)
WRITE(iuelph) (isk(ii), ii = 1, nkstot)
WRITE(iuelph) (ngk_collect(ii), ii = 1, nkstot)
WRITE(iuelph) (ikks_collect(ii), ii = 1, nksqtot)
WRITE(iuelph) (ikqs_collect(ii), ii = 1, nksqtot)
WRITE(iuelph) (eigqts(ii), ii = 1, nat)
WRITE(iuelph) (w2(ii), ii = 1, nmodes)
WRITE(iuelph) ((dyn(ii, jj), ii = 1, nmodes), jj = 1, nmodes)
WRITE(iuelph) ((((el_ph_mat_collect(ii, jj, kk, ll), ii = 1, nbnd), &
jj = 1, nbnd), kk = 1, nksqtot), ll = 1, nmodes)
CLOSE (unit = iuelph, status = 'keep')
ENDIF
CALL cryst_to_cart(nkstot, xk_collect, bg, 1)
DEALLOCATE(xk_collect)
DEALLOCATE(wk_collect)
DEALLOCATE(et_collect)
DEALLOCATE(wg_collect)
DEALLOCATE(ngk_collect)
DEALLOCATE(ikks_collect)
DEALLOCATE(ikqs_collect)
DEALLOCATE(el_ph_mat_collect)
RETURN
END SUBROUTINE elphfil_epa
!----------------------------------------------------------------------------
SUBROUTINE ipoolcollect( length, nks, f_in, nkstot, f_out )
!----------------------------------------------------------------------------
!
! ... as poolcollect, for an integer vector
!
USE mp_pools, ONLY : my_pool_id, npool, kunit, &
inter_pool_comm, intra_pool_comm
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length, nks, nkstot
! first dimension of arrays
! number of k-points per pool
! total number of k-points
INTEGER, INTENT(IN) :: f_in (length,nks)
! pool-distributed function
INTEGER, INTENT(OUT) :: f_out(length,nkstot)
! pool-collected function
!
INTEGER :: nbase, rest, nks1
!
nks1 = kunit * ( nkstot / kunit / npool )
!
rest = ( nkstot - nks1 * npool ) / kunit
!
IF ( ( my_pool_id + 1 ) <= rest ) nks1 = nks1 + kunit
!
IF (nks1.ne.nks) &
call errore('ipoolcollect','inconsistent number of k-points',1)
!
! ... calculates nbase = the position in the list of the first point that
! ... belong to this npool - 1
!
nbase = nks * my_pool_id
!
IF ( ( my_pool_id + 1 ) > rest ) nbase = nbase + rest * kunit
!
! copy the original points in the correct position of the list
!
f_out=0
f_out(:,nbase+1:nbase+nks) = f_in(:,1:nks)
!
CALL mp_sum( f_out, inter_pool_comm )
!
RETURN
!
END SUBROUTINE ipoolcollect
!----------------------------------------------------------------------------
SUBROUTINE jpoolcollect( length, nks, f_in, nkstot, f_out )
!----------------------------------------------------------------------------
!
! ... as ipoolcollect, without kunit and with an index shift
!
USE mp_pools, ONLY : my_pool_id, npool, kunit, &
inter_pool_comm, intra_pool_comm
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length, nks, nkstot
! first dimension of arrays
! number of k-points per pool
! total number of k-points
INTEGER, INTENT(IN) :: f_in (length,nks)
! pool-distributed function
INTEGER, INTENT(OUT) :: f_out(length,nkstot)
! pool-collected function
!
INTEGER :: nbase, rest, nks1
!
nks1 = ( nkstot / npool )
!
rest = ( nkstot - nks1 * npool )
!
IF ( ( my_pool_id + 1 ) <= rest ) nks1 = nks1 + 1
!
IF (nks1.ne.nks) &
call errore('jpoolcollect','inconsistent number of k-points',1)
!
! ... calculates nbase = the position in the list of the first point that
! ... belong to this npool - 1
!
nbase = nks * my_pool_id
!
IF ( ( my_pool_id + 1 ) > rest ) nbase = nbase + rest
!
! copy the original points in the correct position of the list
!
f_out=0
f_out(:,nbase+1:nbase+nks) = f_in(:,1:nks) + nbase * kunit
!
CALL mp_sum( f_out, inter_pool_comm )
!
RETURN
!
END SUBROUTINE jpoolcollect
!-----------------------------------------------------------------------
FUNCTION dos_ef (ngauss, degauss, ef, et, wk, nks, nbnd)
!-----------------------------------------------------------------------

524
PHonon/PH/epa.f90 Normal file
View File

@ -0,0 +1,524 @@
!-----------------------------------------------------------------------
!
! 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
!
! Adv. Energy Mater. 2018, 1800246
! doi:10.1002/aenm.201800246
! https://doi.org/10.1002/aenm.201800246
!
!-----------------------------------------------------------------------
program epa
use kinds, only : dp
use constants, only : ry_to_cmm1, rytoev, degspin, eps12
implicit none
integer, parameter :: nwin = 2
integer, parameter :: uni = 101
integer, parameter :: uno = 102
real(dp), parameter :: epsw2 = (20_dp / ry_to_cmm1)**2
logical :: minus_q
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), &
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
character(len=32) :: s1, s2, s3, s4, s5
real(dp), allocatable :: gdist(:,:,:)
real(dp), allocatable :: gepair(:,:,:)
real(dp), allocatable :: gavg(:,:,:,:)
real(dp), allocatable :: wavg(:)
real(dp), allocatable :: gtot(:,:,:)
integer, allocatable :: gnum(:,:,:)
real(dp), allocatable :: egrid(:,:,:)
real(dp), allocatable :: x_q(:,:)
real(dp), allocatable :: wq(:)
real(dp), allocatable :: w2(:)
complex(dp), allocatable :: dyn(:,:)
real(dp), allocatable :: wk(:)
real(dp), allocatable :: et(:,:)
integer, allocatable :: ikks(:)
integer, allocatable :: ikqs(:)
integer, allocatable :: irt(:,:)
real(dp), allocatable :: rtau(:,:,:)
complex(dp), pointer :: u(:,:)
complex(dp), allocatable :: el_ph_mat(:,:,:,:)
complex(dp), allocatable :: el_ph_sum(:,:)
write(6, '("Reading standard input")')
read(5, '(a)') fni
read(5, '(a)') fno
read(5, '(a)') job
read(5, *) ev, dev, nev, bvmin, bvmax
read(5, *) ec, dec, nec, bcmin, bcmax
read(5, *) gmax, ndist, nepmax
write(6, '(" input file name: ", a)') trim(fni)
write(6, '(" output file name: ", a)') trim(fno)
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)), &
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)), &
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 = ", &
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
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)
read(uni) s1, s2, s3
write(6, '(" title: ", a, ", date: ", a, ", time: ", a)') &
trim(s1), trim(s2), trim(s3)
read(uni) ii, nat, ii, ii, ii, ii, ii, &
ii, nspin, nbnd, nmodes, nqs
write(s1, '(i8)') nqs
write(s2, '(i8)') nbnd
write(s3, '(i8)') nspin
write(s4, '(i8)') nmodes
write(6, '(" nqs = ", a, " nbnd = ", a, " nspin = ", a, &
" nmodes = ", a)') trim(adjustl(s1)), trim(adjustl(s2)), &
trim(adjustl(s3)), trim(adjustl(s4))
if (nspin .eq. 1) then
wspin = 1.0d0 / degspin
else
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
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
if (ijob .eq. 1) then
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))
elseif (ijob .eq. 3) then
allocate(gdist(2, ndist, nwin))
allocate(gepair(4, nepmax, nwin))
allocate(egrid(2, nwin, 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
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)
enddo
enddo
endif
egrid(:,:,:) = egrid(:,:,:) / rytoev
endif
allocate(x_q(3, nqs))
allocate(wq(nqs))
allocate(w2(nmodes))
allocate(dyn(nmodes, nmodes))
allocate(irt(48, nat))
allocate(rtau(3, 48, nat))
allocate(u(nmodes, nmodes))
allocate(el_ph_sum(nmodes, nmodes))
read(uni)
read(uni)
read(uni)
read(uni)
read(uni)
read(uni)
read(uni) ((at(ii, jj), ii = 1, 3), jj = 1, 3)
read(uni) ((bg(ii, jj), ii = 1, 3), jj = 1, 3)
read(uni)
read(uni)
read(uni)
read(uni) ((x_q(ii, jj), ii = 1, 3), jj = 1, nqs)
read(uni) (wq(ii), ii = 1, nqs)
read(uni)
call cryst_to_cart(nqs, x_q, bg, 1)
if (ijob .eq. 1 .or. ijob .eq. 2) then
gavg(:,:,:,:) = 0.0d0
wavg(:) = 0.0d0
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
enddo
enddo
nepair(:) = 0
endif
if (ijob .eq. 2) gnum(:,:,:) = 0
do iq = 1, nqs
xq(:) = x_q(:, iq)
read(uni) nsymq, irotmq, ii, ii, nkstot, nksqtot
write(s1, '(i8)') iq
write(s2, '(i8)') nkstot
write(s3, '(i8)') nksqtot
write(6, '(" iq = ", a, " nkstot = ", a, " nksqtot = ", a)') &
trim(adjustl(s1)), trim(adjustl(s2)), trim(adjustl(s3))
allocate(wk(nkstot))
allocate(et(nbnd, nkstot))
allocate(ikks(nksqtot))
allocate(ikqs(nksqtot))
allocate(el_ph_mat(nbnd, nbnd, nksqtot, nmodes))
read(uni) minus_q
read(uni)
read(uni)
read(uni) (((rtau(ii, jj, kk), ii = 1, 3), jj = 1, 48), kk = 1, nat)
read(uni)
read(uni)
read(uni) ((u(ii, jj), ii = 1, nmodes), jj = 1, nmodes)
read(uni)
read(uni)
read(uni)
read(uni)
read(uni) (((s(ii, jj, kk), ii = 1, 3), jj = 1, 3), kk = 1, 48)
read(uni) (invs(ii), ii = 1, 48)
read(uni)
read(uni)
read(uni)
read(uni)
read(uni)
read(uni) ((irt(ii, jj), ii = 1, 48), jj = 1, nat)
read(uni)
read(uni) (wk(ii), ii = 1, nkstot)
read(uni) ((et(ii, jj), ii = 1, nbnd), jj = 1, nkstot)
read(uni)
read(uni)
read(uni)
read(uni) (ikks(ii), ii = 1, nksqtot)
read(uni) (ikqs(ii), ii = 1, nksqtot)
read(uni)
read(uni) (w2(ii), ii = 1, nmodes)
read(uni) ((dyn(ii, jj), ii = 1, nmodes), jj = 1, nmodes)
read(uni) ((((el_ph_mat(ii, jj, kk, ll), ii = 1, nbnd), &
jj = 1, nbnd), kk = 1, nksqtot), ll = 1, nmodes)
if (ijob .eq. 1) then
do ibnd = 1, nbnd
do jbnd = 1, nbnd
el_ph_sum(:,:) = (0.0d0, 0.0d0)
do ik = 1, nksqtot
ikk = ikks(ik)
ikq = ikqs(ik)
weight = wq(iq) * wk(ikk) * wspin
if (ibnd .eq. 1 .and. jbnd .eq. 1) gtot(1, 1, 1) = &
gtot(1, 1, 1) + weight
do jpert = 1, nmodes
do ipert = 1, nmodes
el_ph_sum(ipert, jpert) = el_ph_sum(ipert, jpert) &
+ conjg(el_ph_mat(jbnd, ibnd, ik, ipert)) &
* el_ph_mat(jbnd, ibnd, ik, jpert) * weight
enddo
enddo
enddo
call symdyn_munu_new(el_ph_sum, u, xq, s, invs, rtau, irt, &
at, bg, nsymq, nat, irotmq, minus_q)
do nu = 1, nmodes
if (w2(nu) > epsw2) then
factor = rytoev**2 / (2.0d0 * sqrt(w2(nu)))
else
factor = 0.0d0
endif
gbuf = 0.0d0
do mu = 1, nmodes
do vu = 1, nmodes
gbuf = gbuf + dble(conjg(dyn(mu, nu)) * &
el_ph_sum(mu, vu) * dyn(vu, nu))
enddo
enddo
gavg(nu, jbnd, ibnd, 1) = gavg(nu, jbnd, ibnd, 1) + gbuf * factor
enddo
enddo
enddo
elseif (ijob .eq. 2) then
do ii = 1, nwin
do jj = 1, nemax
do kk = 1, nemax
el_ph_sum(:,:) = (0.0d0, 0.0d0)
do ik = 1, nksqtot
ikk = ikks(ik)
ikq = ikqs(ik)
weight = wq(iq) * wk(ikk) * wspin
do ibnd = bpair(1, ii), bpair(2, ii)
if (et(ibnd, ikk) .gt. egrid(1, ii, jj) .and. &
et(ibnd, ikk) .le. egrid(2, ii, jj)) then
do jbnd = bpair(1, ii), bpair(2, ii)
if (et(jbnd, ikq) .gt. egrid(1, ii, kk) .and. &
et(jbnd, ikq) .le. egrid(2, ii, kk)) then
gnum(kk, jj, ii) = gnum(kk, jj, ii) + 1
gtot(kk, jj, ii) = gtot(kk, jj, ii) + weight
do jpert = 1, nmodes
do ipert = 1, nmodes
el_ph_sum(ipert, jpert) = el_ph_sum(ipert, jpert) &
+ conjg(el_ph_mat(jbnd, ibnd, ik, ipert)) &
* el_ph_mat(jbnd, ibnd, ik, jpert) * weight
enddo
enddo
endif
enddo
endif
enddo
enddo
call symdyn_munu_new(el_ph_sum, u, xq, s, invs, rtau, irt, at, &
bg, nsymq, nat, irotmq, minus_q)
do nu = 1, nmodes
if (w2(nu) > epsw2) then
factor = rytoev**2 / (2.0d0 * sqrt(w2(nu)))
else
factor = 0.0d0
endif
gbuf = 0.0d0
do mu = 1, nmodes
do vu = 1, nmodes
gbuf = gbuf + dble(conjg(dyn(mu, nu)) * &
el_ph_sum(mu, vu) * dyn(vu, nu))
enddo
enddo
gavg(nu, kk, jj, ii) = gavg(nu, kk, jj, ii) + gbuf * factor
enddo
enddo
enddo
enddo
elseif (ijob .eq. 3) then
do ii = 1, nwin
do ik = 1, nksqtot
ikk = ikks(ik)
ikq = ikqs(ik)
weight = wq(iq) * wk(ikk) * wspin
do ibnd = bpair(1, ii), bpair(2, ii)
if (et(ibnd, ikk) .gt. egrid(1, ii, 1) .and. &
et(ibnd, ikk) .le. egrid(2, ii, 1)) then
do jbnd = bpair(1, ii), bpair(2, ii)
if (et(jbnd, ikq) .gt. egrid(1, ii, 1) .and. &
et(jbnd, ikq) .le. egrid(2, ii, 1)) then
do jpert = 1, nmodes
do ipert = 1, nmodes
el_ph_sum(ipert, jpert) = &
conjg(el_ph_mat(jbnd, ibnd, ik, ipert)) &
* el_ph_mat(jbnd, ibnd, ik, jpert)
enddo
enddo
call symdyn_munu_new(el_ph_sum, u, xq, s, invs, &
rtau, irt, at, bg, nsymq, nat, irotmq, minus_q)
gsum = 0.0d0
do nu = 1, nmodes
if (w2(nu) > epsw2) then
factor = rytoev**2 / (2.0d0 * sqrt(w2(nu)))
else
factor = 0.0d0
endif
gbuf = 0.0d0
do mu = 1, nmodes
do vu = 1, nmodes
gbuf = gbuf + dble(conjg(dyn(mu, nu)) * &
el_ph_sum(mu, vu) * dyn(vu, nu))
enddo
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
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
endif
enddo
endif
enddo
enddo
enddo
endif
if (ijob .eq. 1 .or. ijob .eq. 2) then
weight = wq(iq)
wtot = wtot + weight
do nu = 1, nmodes
if (w2(nu) > epsw2) then
wavg(nu) = wavg(nu) + weight * sqrt(w2(nu))
endif
enddo
endif
deallocate(wk)
deallocate(et)
deallocate(ikks)
deallocate(ikqs)
deallocate(el_ph_mat)
enddo
close(uni, status = 'keep')
if (ijob .eq. 1) then
do ibnd = 1, nbnd
do jbnd = 1, nbnd
do nu = 1, nmodes
gavg(nu, jbnd, ibnd, 1) = gavg(nu, jbnd, ibnd, 1) / gtot(1, 1, 1)
enddo
enddo
enddo
elseif (ijob .eq. 2) then
do ii = 1, nwin
do jj = 1, nemax
do kk = 1, nemax
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)
enddo
endif
enddo
enddo
enddo
endif
if (ijob .eq. 1 .or. ijob .eq. 2) then
do nu = 1, nmodes
wavg(nu) = wavg(nu) * ry_to_cmm1 / wtot
enddo
endif
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
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), &
sum(gavg(:, kk, jj, ii)), gnum(kk, jj, ii), gtot(kk, jj, ii), &
ii = 1, nwin)
enddo
enddo
endif
open(uno, file = fno, form = 'formatted', status = 'replace')
write(6, '("Writing file ", a)') trim(fno)
if (ijob .eq. 1) then
write(uno, '(" mode <w> (cm^-1)")')
write(uno, '("--------------------------")')
do nu = 1, nmodes
write(uno, '(i8, e18.8)') nu, wavg(nu)
enddo
write(uno, *)
write(uno, '(" band band mode <|g|^2> (eV^2)")')
write(uno, '("------------------------------------------")')
do ibnd = 1, nbnd
do jbnd = 1, nbnd
do nu = 1, nmodes
write(uno, '(3i8, e18.8)') ibnd, jbnd, nu, gavg(nu, jbnd, ibnd, 1)
enddo
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)
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)
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)
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)
enddo
enddo
endif
close(uno, status = 'keep')
if (ijob .eq. 1 .or. ijob .eq. 2) then
deallocate(gavg)
deallocate(gtot)
deallocate(wavg)
endif
if (ijob .eq. 2) deallocate(gnum)
if (ijob .eq. 3) then
deallocate(gdist)
deallocate(gepair)
endif
if (ijob .eq. 2 .or. ijob .eq. 3) deallocate(egrid)
deallocate(x_q)
deallocate(wq)
deallocate(w2)
deallocate(dyn)
deallocate(irt)
deallocate(rtau)
deallocate(u)
deallocate(el_ph_sum)
end program epa

View File

@ -66,7 +66,7 @@ SUBROUTINE phq_readin()
USE freq_ph, ONLY : fpol, fiu, nfs
USE cryst_ph, ONLY : magnetic_sym
USE ph_restart, ONLY : ph_readfile
USE el_phon, ONLY : elph,elph_mat,elph_simple,elph_nbnd_min, elph_nbnd_max, &
USE el_phon, ONLY : elph,elph_mat,elph_simple,elph_epa,elph_nbnd_min, elph_nbnd_max, &
el_ph_sigma, el_ph_nsigma, el_ph_ngauss,auxdvscf
USE dfile_star, ONLY : drho_star, dvscf_star
@ -375,20 +375,29 @@ SUBROUTINE phq_readin()
elph=.true.
elph_mat=.false.
elph_simple=.true.
elph_epa=.false.
CASE( 'epa' )
elph=.true.
elph_mat=.false.
elph_simple=.false.
elph_epa=.true.
CASE( 'Wannier' )
elph=.true.
elph_mat=.true.
elph_simple=.false.
elph_epa=.false.
auxdvscf=trim(fildvscf)
CASE( 'interpolated' )
elph=.true.
elph_mat=.false.
elph_simple=.false.
elph_epa=.false.
! YAMBO >
CASE( 'yambo' )
elph=.true.
elph_mat=.false.
elph_simple=.false.
elph_epa=.false.
elph_yambo=.true.
nogg=.true.
auxdvscf=trim(fildvscf)
@ -396,6 +405,7 @@ SUBROUTINE phq_readin()
elph=.false.
elph_mat=.false.
elph_simple=.false.
elph_epa=.false.
elph_yambo=.false.
dvscf_yambo=.true.
nogg=.true.
@ -423,6 +433,7 @@ SUBROUTINE phq_readin()
elph=.false.
elph_mat=.false.
elph_simple=.false.
elph_epa=.false.
END SELECT
! YAMBO >
IF (.not.elph_yambo) then

View File

@ -11,7 +11,7 @@ SUBROUTINE q_points ( )
USE kinds, only : dp
USE io_global, ONLY : stdout, ionode, ionode_id
USE disp, ONLY : nq1, nq2, nq3, x_q, nqs, lgamma_iq
USE disp, ONLY : nq1, nq2, nq3, x_q, nqs, lgamma_iq, wq
USE output, ONLY : fildyn
USE symm_base, ONLY : nsym, s, time_reversal, t_rev, invs
USE cell_base, ONLY : at, bg
@ -25,7 +25,7 @@ SUBROUTINE q_points ( )
integer :: i, iq, ierr, iudyn = 26
logical :: exist_gamma, check, skip_equivalence=.FALSE.
logical, external :: check_q_points_sym
real(DP), allocatable :: xq(:,:), wq(:)
real(DP), allocatable :: xq(:,:), w_q(:)
INTEGER :: nqmax
!
@ -37,20 +37,22 @@ SUBROUTINE q_points ( )
nqmax= nq1 * nq2 * nq3
allocate (wq(nqmax))
allocate (w_q(nqmax))
allocate (xq(3,nqmax))
if(lshift_q) then
call kpoint_grid( nsym, time_reversal, skip_equivalence, s, t_rev, bg, nqmax,&
& 1,1,1, nq1,nq2,nq3, nqs, xq, wq )
& 1,1,1, nq1,nq2,nq3, nqs, xq, w_q )
else
call kpoint_grid( nsym, time_reversal, skip_equivalence, s, t_rev, bg, nqmax,&
& 0,0,0, nq1,nq2,nq3, nqs, xq, wq )
& 0,0,0, nq1,nq2,nq3, nqs, xq, w_q )
end if
allocate(wq(nqs))
allocate(x_q(3,nqs))
allocate(lgamma_iq(nqs))
wq(:)=w_q(1:nqs)
x_q(:,:)=xq(:,1:nqs)
deallocate (xq)
deallocate (wq)
deallocate (w_q)
!
! Check if the Gamma point is one of the points and put
! it in the first position (it should already be the first)