scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Further cleaning/beautification.

This commit is contained in:
Samuel Ponce 2019-11-20 15:31:37 +01:00
parent 049bef5119
commit 5408f9167a
18 changed files with 1277 additions and 1275 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
EPW/src/Makefile
View File

@ -30,6 +30,7 @@ division.o \
rigid_epw.o \
io_epw.o \
io_transport.o \
io_selfen.o \
wigner.o \
wan2bloch.o \
cum_mod.o \
@ -40,7 +41,6 @@ kfold.o \
dynmat_asr.o \
io_eliashberg.o \
utilities.o \
plot.o \
eliashbergcom.o \
supercond.o \
supercond_iso.o \

6
EPW/src/bcast_epw_input.f90
View File

@ -45,8 +45,8 @@
tempsmin, tempsmax, temps, delta_approx, title, &
scattering, scattering_serta, scattering_0rta, &
int_mob, scissor, carrier, ncarrier, &
restart, restart_freq, prtgkk, nel, meff, epsiheg,&
scatread, restart, restart_freq, restart_filq, &
restart, restart_step, prtgkk, nel, meff, epsiheg,&
scatread, restart, restart_step, restart_filq, &
lphase, omegamin, omegamax, omegastep, n_r, &
mob_maxiter, use_ws, epmatkqread, selecqread, &
scdm_sigma, assume_metal
@ -170,7 +170,7 @@
CALL mp_bcast(nstemp , meta_ionode_id, world_comm)
CALL mp_bcast(nsiter , meta_ionode_id, world_comm)
CALL mp_bcast(nw_specfun , meta_ionode_id, world_comm)
CALL mp_bcast(restart_freq, meta_ionode_id, world_comm)
CALL mp_bcast(restart_step, meta_ionode_id, world_comm)
CALL mp_bcast(scr_typ , meta_ionode_id, world_comm)
CALL mp_bcast(bnd_cum , meta_ionode_id, world_comm)
CALL mp_bcast(mob_maxiter , meta_ionode_id, world_comm)

14
EPW/src/bloch2wan.f90
View File

@ -224,13 +224,13 @@
ENDDO
!
! RMDB
DO ir = 1, nrr
DO jbnd = 1, nbndsub
DO ibnd = 1, nbndsub
WRITE(iudecayH, '(5I5,6F12.6)') irvec(:, ir), ibnd, jbnd, chw(ibnd, jbnd, ir) * ryd2ev
ENDDO
ENDDO
ENDDO
!DO ir = 1, nrr
! DO jbnd = 1, nbndsub
! DO ibnd = 1, nbndsub
! WRITE(iudecayH, '(5I5,6F12.6)') irvec(:, ir), ibnd, jbnd, chw(ibnd, jbnd, ir) * ryd2ev
! ENDDO
! ENDDO
!ENDDO
!
CLOSE(iudecayH)
!

4
EPW/src/elphon_shuffle_wrap.f90
View File

@ -66,7 +66,7 @@
USE phus, ONLY : int1, int1_nc, int2, int2_so, int4, int4_nc, int5, &
int5_so, alphap
USE kfold, ONLY : shift, createkmap_pw2, createkmap
USE low_lvl, ONLY : set_ndnmbr, eqvect_strict, read_modes
USE low_lvl, ONLY : set_ndnmbr, eqvect_strict, read_disp_pattern
USE io_epw, ONLY : read_ifc, readdvscf
USE poolgathering, ONLY : poolgather
USE rigid_epw, ONLY : compute_umn_c
@ -465,7 +465,7 @@
IF (.NOT. exst) CALL errore('elphon_shuffle_wrap', &
'cannot open file for reading or writing', ierr)
CALL iotk_open_read(iunpun, FILE = TRIM(filename), binary = .FALSE., ierr = ierr)
CALL read_modes(iunpun, iq_irr, ierr)
CALL read_disp_pattern(iunpun, iq_irr, ierr)
IF (ierr /= 0) CALL errore('elphon_shuffle_wrap', ' Problem with modes file', 1)
IF (meta_ionode) CALL iotk_close_read(iunpun)
ENDIF

21
EPW/src/ephwann_shuffle.f90
View File

@ -33,7 +33,7 @@
iterative_bte, longrange, scatread, nqf1, prtgkk, &
nqf2, nqf3, mp_mesh_k, restart, plselfen, &
specfun_pl, lindabs, use_ws, epbread, &
epmatkqread, selecqread, restart_freq, nsmear, &
epmatkqread, selecqread, restart_step, nsmear, &
nqc1, nqc2, nqc3, nkc1, nkc2, nkc3, assume_metal
USE control_flags, ONLY : iverbosity
USE noncollin_module, ONLY : noncolin
@ -66,10 +66,10 @@
USE io_eliashberg, ONLY : write_ephmat, count_kpoints, kmesh_fine, kqmap_fine
USE transport, ONLY : transport_coeffs, scattering_rate_q
USE grid, ONLY : qwindow
USE printing, ONLY : print_gkk
USE printing, ONLY : print_gkk, plot_band
USE io_epw, ONLY : rwepmatw, epw_read, epw_write
USE io_transport, ONLY : electron_read, tau_read, iter_open, print_ibte, &
iter_merge, spectral_read
USE io_transport, ONLY : tau_read, iter_open, print_ibte, iter_merge
USE io_selfen, ONLY : selfen_el_read, spectral_read
USE transport_iter,ONLY : iter_restart
USE close_epw, ONLY : iter_close
USE division, ONLY : fkbounds
@ -80,12 +80,13 @@
USE low_lvl, ONLY : system_mem_usage, mem_size
USE utilities, ONLY : compute_dos, broadening, fermicarrier, fermiwindow
USE grid, ONLY : loadqmesh_serial, loadkmesh_para, load_rebal
USE selfen, ONLY : selfen_phon_q, selfen_elec_q, selfen_pl_q
USE spectral_func, ONLY : spectral_func_el_q, spectral_func_ph_q, spectral_func_pl_q
USE selfen, ONLY : selfen_phon_q, selfen_elec_q, selfen_pl_q, &
nesting_fn_q
USE spectral_func, ONLY : spectral_func_el_q, spectral_func_ph_q, a2f_main, &
spectral_func_pl_q
USE io_epw, ONLY : read_ifc
USE rigid_epw, ONLY : rpa_epsilon, tf_epsilon, compute_umn_f, rgd_blk_epw_fine
USE indabs, ONLY : indabs_main, renorm_eig
USE plot, ONLY : nesting_fn_q, a2f_main, plot_band
#if defined(__MPI)
USE parallel_include, ONLY : MPI_MODE_RDONLY, MPI_INFO_NULL, MPI_OFFSET_KIND, &
MPI_OFFSET
@ -994,7 +995,7 @@
! Restart in SERTA case or self-energy (electron or plasmon) case
IF (restart) THEN
IF (elecselfen .OR. plselfen) THEN
CALL electron_read(iq_restart, totq, nktotf, sigmar_all, sigmai_all, zi_all)
CALL selfen_el_read(iq_restart, totq, nktotf, sigmar_all, sigmai_all, zi_all)
ENDIF
IF (specfun_el .OR. specfun_pl) THEN
CALL spectral_read(iq_restart, totq, nktotf, esigmar_all, esigmai_all)
@ -1104,7 +1105,7 @@
! elecselfen = true as nothing happen during the calculation otherwise.
!
IF (.NOT. phonselfen) THEN
IF (MOD(iqq, restart_freq) == 0) THEN
IF (MOD(iqq, restart_step) == 0) THEN
WRITE(stdout, '(5x, a, i10, a, i10)' ) 'Progression iq (fine) = ', iqq, '/', totq
ENDIF
ENDIF
@ -1322,7 +1323,7 @@
!
ENDDO ! end loop over k points
!
IF (MOD(iqq, restart_freq) == 0 .AND. adapt_smearing) THEN
IF (MOD(iqq, restart_step) == 0 .AND. adapt_smearing) THEN
! Min non-zero value
valmin(:) = zero
valmin(my_pool_id + 1) = 100.0d0

19
EPW/src/ephwann_shuffle_mem.f90
View File

@ -37,7 +37,7 @@
iterative_bte, longrange, scatread, nqf1, prtgkk, &
nqf2, nqf3, mp_mesh_k, restart, plselfen, &
specfun_pl, lindabs, use_ws, &
epmatkqread, selecqread, restart_freq, nsmear, &
epmatkqread, selecqread, restart_step, nsmear, &
nkc1, nkc2, nkc3, nqc1, nqc2, nqc3, assume_metal
USE control_flags, ONLY : iverbosity
USE noncollin_module, ONLY : noncolin
@ -70,10 +70,10 @@
USE io_eliashberg, ONLY : write_ephmat, count_kpoints, kmesh_fine, kqmap_fine
USE transport, ONLY : transport_coeffs, scattering_rate_q
USE grid, ONLY : qwindow
USE printing, ONLY : print_gkk
USE printing, ONLY : print_gkk, plot_band
USE io_epw, ONLY : rwepmatw, epw_read, epw_write
USE io_transport, ONLY : electron_read, tau_read, iter_open, print_ibte, &
iter_merge, spectral_read
USE io_transport, ONLY : tau_read, iter_open, print_ibte, iter_merge
USE io_selfen, ONLY : selen_el_read, spectral_read
USE transport_iter,ONLY : iter_restart
USE close_epw, ONLY : iter_close
USE division, ONLY : fkbounds
@ -84,11 +84,12 @@
USE low_lvl, ONLY : system_mem_usage, mem_size
USE utilities, ONLY : compute_dos, fermicarrier, fermiwindow
USE grid, ONLY : loadqmesh_serial, loadkmesh_para, load_rebal
USE selfen, ONLY : selfen_phon_q, selfen_elec_q, selfen_pl_q
USE spectral_func, ONLY : spectral_func_el_q, spectral_func_ph_q, spectral_func_pl_q
USE selfen, ONLY : selfen_phon_q, selfen_elec_q, selfen_pl_q, &
nesting_fn_q
USE spectral_func, ONLY : spectral_func_el_q, spectral_func_ph_q, a2f_main, &
spectral_func_pl_q
USE rigid_epw, ONLY : rpa_epsilon, tf_epsilon, compute_umn_f, rgd_blk_epw_fine_mem
USE indabs, ONLY : indabs_main, renorm_eig
USE plot, ONLY : nesting_fn_q, a2f_main, plot_band
#if defined(__MPI)
USE parallel_include, ONLY : MPI_MODE_RDONLY, MPI_INFO_NULL, MPI_OFFSET_KIND, &
MPI_OFFSET
@ -965,7 +966,7 @@
! Restart in SERTA case or self-energy (electron or plasmon) case
IF (restart) THEN
IF (elecselfen .OR. plselfen) THEN
CALL electron_read(iq_restart, totq, nktotf, sigmar_all, sigmai_all, zi_all)
CALL selfen_el_read(iq_restart, totq, nktotf, sigmar_all, sigmai_all, zi_all)
ENDIF
IF (specfun_el .OR. specfun_pl) THEN
CALL spectral_read(iq_restart, totq, nktotf, esigmar_all, esigmai_all)
@ -1079,7 +1080,7 @@
! elecselfen = true as nothing happen during the calculation otherwise.
!
IF (.NOT. phonselfen) THEN
IF (MOD(iqq, restart_freq) == 0) THEN
IF (MOD(iqq, restart_step) == 0) THEN
WRITE(stdout, '(5x,a,i10,a,i10)' ) 'Progression iq (fine) = ', iqq, '/', totq
ENDIF
ENDIF

8
EPW/src/epw_readin.f90
View File

@ -51,7 +51,7 @@
nsiter, conv_thr_racon, specfun_el, specfun_ph, nbndskip, &
system_2d, delta_approx, title, int_mob, scissor, &
iterative_bte, scattering, selecqread, epmatkqread, &
ncarrier, carrier, scattering_serta, restart, restart_freq,&
ncarrier, carrier, scattering_serta, restart, restart_step,&
scattering_0rta, longrange, shortrange, scatread, use_ws, &
restart_filq, prtgkk, nel, meff, epsiheg, lphase, &
omegamin, omegamax, omegastep, n_r, lindabs, mob_maxiter, &
@ -130,7 +130,7 @@
specfun_el, specfun_ph, wmin_specfun, wmax_specfun, nw_specfun, &
delta_approx, scattering, int_mob, scissor, ncarrier, carrier, &
iterative_bte, scattering_serta, scattering_0rta, longrange, shortrange,&
scatread, restart, restart_freq, restart_filq, prtgkk, nel, meff, &
scatread, restart, restart_step, restart_filq, prtgkk, nel, meff, &
epsiheg, lphase, omegamin, omegamax, omegastep, n_r, lindabs, &
mob_maxiter, auto_projections, scdm_proj, scdm_entanglement, scdm_mu, &
scdm_sigma, assume_metal
@ -266,7 +266,7 @@
! specfun_ph : if .TRUE. calculate phonon spectral function due to e-p interaction
! specfun_pl : if .TRUE. calculate plason spectral function
! restart : if .TRUE. a run can be restarted from the interpolation level
! restart_freq : Create a restart point every restart_freq q/k-points
! restart_step : Create a restart point every restart_step q/k-points
! restart_filq : Use to merge different q-grid scattering rates (name of the file)
! scattering : if .TRUE. scattering rates are calculated
! scattering_serta: if .TRUE. scattering rates are calculated using self-energy relaxation-time-approx
@ -355,7 +355,7 @@
epwread = .FALSE.
epwwrite = .TRUE.
restart = .FALSE.
restart_freq = 100
restart_step = 100
wannierize = .FALSE.
write_wfn = .FALSE.
kmaps = .FALSE.

4
EPW/src/epwcom.f90
View File

@ -230,8 +230,8 @@
!! nr. of iterations used in broyden mixing scheme
INTEGER :: nw_specfun
!! nr. of bins for frequency in electron spectral function due to e-p interaction
INTEGER :: restart_freq
!! Create a restart point during the interpolation part every restart_freq q/k-points.
INTEGER :: restart_step
!! Create a restart point during the interpolation part every restart_step q/k-points.
!
REAL(KIND = DP) :: degaussw
!! smearing width for Fermi surface average in e-ph coupling after wann interp

8
EPW/src/grid.f90
View File

@ -1289,7 +1289,7 @@
USE mp, ONLY : mp_sum, mp_bcast
USE constants_epw, ONLY : twopi, ci, zero, eps6, ryd2ev, czero
USE epwcom, ONLY : nbndsub, fsthick, use_ws, mp_mesh_k, nkf1, nkf2, &
nkf3, iterative_bte, restart_freq, scissor
nkf3, iterative_bte, restart_step, scissor
USE noncollin_module, ONLY : noncolin
USE pwcom, ONLY : ef, nelec
USE cell_base, ONLY : bg
@ -1512,7 +1512,7 @@
totq = totq + 1
selecq(totq) = iq
!
IF (MOD(totq, restart_freq) == 0) THEN
IF (MOD(totq, restart_step) == 0) THEN
WRITE(stdout, '(5x,a,i15,i15)')'Number selected, total', totq, iq
ENDIF
ENDIF
@ -1602,7 +1602,7 @@
IF (SUM(found(:)) > 0) THEN
totq = totq + 1
selecq(totq) = iq
IF (MOD(totq, restart_freq) == 0) THEN
IF (MOD(totq, restart_step) == 0) THEN
WRITE(stdout, '(5x,a,i12,i12)') 'Number selected, total', totq, iq
ENDIF
ENDIF
@ -1637,7 +1637,7 @@
IF (SUM(found(:)) > 0) THEN
totq = totq + 1
selecq(totq) = iq
IF (MOD(totq, restart_freq) == 0) THEN
IF (MOD(totq, restart_step) == 0) THEN
WRITE(stdout, '(5x,a,i12,i12)')'Number selected, total', totq, iq
ENDIF
ENDIF

459
EPW/src/io_selfen.f90 Normal file
View File

@ -0,0 +1,459 @@
!
! Copyright (C) 2016-2019 Samuel Ponce', Roxana Margine, 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 .
!
!
!----------------------------------------------------------------------
MODULE io_selfen
!----------------------------------------------------------------------
!!
!! This module contains various writing or reading routines related to self-energies.
!! Most of them are for restart purposes.
!!
IMPLICIT NONE
!
CONTAINS
!
!----------------------------------------------------------------------------
SUBROUTINE selfen_el_write(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
!----------------------------------------------------------------------------
!!
!! Write self-energy
!!
USE kinds, ONLY : DP
USE elph2, ONLY : lower_bnd, upper_bnd, nbndfst
USE io_var, ONLY : iufilsigma_all
USE io_files, ONLY : diropn
USE constants_epw, ONLY : zero
USE mp, ONLY : mp_barrier
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: iqq
!! Current q-point
INTEGER, INTENT(in) :: totq
!! Total number of q-points
INTEGER, INTENT(in) :: nktotf
!! Total number of k-points
REAL(KIND = DP), INTENT(inout) :: sigmar_all(nbndfst, nktotf)
!! Real part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(inout) :: sigmai_all(nbndfst, nktotf)
!! Imaginary part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(inout) :: zi_all(nbndfst, nktotf)
!! Z parameter of electron-phonon self-energy accross all pools
!
! Local variables
LOGICAL :: exst
!! Does the file exist
INTEGER :: i
!! Iterative index
INTEGER :: ik
!! K-point index
INTEGER :: ibnd
!! Local band index
INTEGER :: lsigma_all
!! Length of the vector
REAL(KIND = DP) :: aux(3 * nbndfst * nktotf + 2)
!! Vector to store the array
!
IF (mpime == ionode_id) THEN
!
lsigma_all = 3 * nbndfst * nktotf + 2
! First element is the current q-point
aux(1) = REAL(iqq - 1, KIND = DP) ! we need to start at the next q
! Second element is the total number of q-points
aux(2) = REAL(totq, KIND = DP)
!
i = 2
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
aux(i) = sigmar_all(ibnd, ik)
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
aux(i) = sigmai_all(ibnd, ik)
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
aux(i) = zi_all(ibnd, ik)
ENDDO
ENDDO
CALL diropn(iufilsigma_all, 'sigma_restart', lsigma_all, exst)
CALL davcio(aux, lsigma_all, iufilsigma_all, 1, +1)
CLOSE(iufilsigma_all)
ENDIF
!
! Make everythin 0 except the range of k-points we are working on
IF (lower_bnd > 1) THEN
sigmar_all(:, 1:lower_bnd - 1) = zero
sigmai_all(:, 1:lower_bnd - 1) = zero
zi_all(:, 1:lower_bnd - 1) = zero
ENDIF
IF (upper_bnd < nktotf) THEN
sigmar_all(:, upper_bnd + 1:nktotf) = zero
sigmai_all(:, upper_bnd + 1:nktotf) = zero
zi_all(:, upper_bnd + 1:nktotf) = zero
ENDIF
!
!----------------------------------------------------------------------------
END SUBROUTINE selfen_el_write
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE selfen_el_read(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
!----------------------------------------------------------------------------
!!
!! Self-energy reading
!!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE elph2, ONLY : lower_bnd, upper_bnd, nbndfst
USE io_var, ONLY : iufilsigma_all
USE io_files, ONLY : prefix, tmp_dir, diropn
USE constants_epw, ONLY : zero
USE mp, ONLY : mp_barrier, mp_bcast
USE mp_world, ONLY : mpime, world_comm
USE io_global, ONLY : ionode_id
!
IMPLICIT NONE
!
INTEGER, INTENT(inout) :: iqq
!! Current q-point
INTEGER, INTENT(in) :: totq
!! Total number of q-points
INTEGER, INTENT(in) :: nktotf
!! Total number of k-points
REAL(KIND = DP), INTENT(out) :: sigmar_all(nbndfst, nktotf)
!! Real part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(out) :: sigmai_all(nbndfst, nktotf)
!! Imaginary part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(out) :: zi_all(nbndfst, nktotf)
!! Z parameter of electron-phonon self-energy accross all pools
!
! Local variables
LOGICAL :: exst
!! Does the file exist
INTEGER :: i
!! Iterative index
INTEGER :: ik
!! K-point index
INTEGER :: ibnd
!! Local band index
INTEGER :: lsigma_all
!! Length of the vector
INTEGER :: nqtotf_read
!! Total number of q-point read
REAL(KIND = DP) :: aux(3 * nbndfst * nktotf + 2)
!! Vector to store the array
!
CHARACTER(LEN = 256) :: name1
!
IF (mpime == ionode_id) THEN
!
! First inquire if the file exists
#if defined(__MPI)
name1 = TRIM(tmp_dir) // TRIM(prefix) // '.sigma_restart1'
#else
name1 = TRIM(tmp_dir) // TRIM(prefix) // '.sigma_restart'
#endif
INQUIRE(FILE = name1, EXIST = exst)
!
IF (exst) THEN ! read the file
!
lsigma_all = 3 * nbndfst * nktotf + 2
CALL diropn(iufilsigma_all, 'sigma_restart', lsigma_all, exst)
CALL davcio(aux, lsigma_all, iufilsigma_all, 1, -1)
!
! First element is the iteration number
iqq = INT(aux(1))
iqq = iqq + 1 ! we need to start at the next q
nqtotf_read = INT(aux(2))
IF (nqtotf_read /= totq) CALL errore('selfen_el_read', &
&'Error: The current total number of q-point is not the same as the read one. ', 1)
!
i = 2
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
sigmar_all(ibnd, ik) = aux(i)
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
sigmai_all(ibnd, ik) = aux(i)
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
zi_all(ibnd, ik) = aux(i)
ENDDO
ENDDO
CLOSE(iufilsigma_all)
ENDIF
ENDIF
!
CALL mp_bcast(exst, ionode_id, world_comm)
!
IF (exst) THEN
CALL mp_bcast(iqq, ionode_id, world_comm)
CALL mp_bcast(sigmar_all, ionode_id, world_comm)
CALL mp_bcast(sigmai_all, ionode_id, world_comm)
CALL mp_bcast(zi_all, ionode_id, world_comm)
!
! Make everythin 0 except the range of k-points we are working on
IF (lower_bnd > 1) THEN
sigmar_all(:, 1:lower_bnd - 1) = zero
sigmai_all(:, 1:lower_bnd - 1) = zero
zi_all(:, 1:lower_bnd - 1) = zero
ENDIF
IF (upper_bnd < nktotf) THEN
sigmar_all(:, upper_bnd + 1:nktotf) = zero
sigmai_all(:, upper_bnd + 1:nktotf) = zero
zi_all(:, upper_bnd + 1:nktotf) = zero
ENDIF
!
WRITE(stdout, '(a,i10,a,i10)' ) ' Restart from: ', iqq,'/', totq
ENDIF
!
!----------------------------------------------------------------------------
END SUBROUTINE selfen_el_read
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE spectral_write(iqq, totq, nktotf, esigmar_all, esigmai_all)
!----------------------------------------------------------------------------
!!
!! Write self-energy
!!
USE kinds, ONLY : DP
USE elph2, ONLY : lower_bnd, upper_bnd, nbndfst
USE epwcom, ONLY : wmin_specfun, wmax_specfun, nw_specfun
USE io_var, ONLY : iufilesigma_all
USE io_files, ONLY : diropn
USE constants_epw, ONLY : zero
USE mp, ONLY : mp_barrier
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: iqq
!! Current q-point
INTEGER, INTENT(in) :: totq
!! Total number of q-points
INTEGER, INTENT(in) :: nktotf
!! Total number of k-points
REAL(KIND = DP), INTENT(inout) :: esigmar_all(nbndfst, nktotf, nw_specfun)
!! Real part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(inout) :: esigmai_all(nbndfst, nktotf, nw_specfun)
!! Imaginary part of the electron-phonon self-energy accross all pools
!
! Local variables
LOGICAL :: exst
!! Does the file exist
INTEGER :: i
!! Iterative index
INTEGER :: ik
!! K-point index
INTEGER :: ibnd
!! Local band index
INTEGER :: iw
!! Counter on the frequency
INTEGER :: lesigma_all
!! Length of the vector
REAL(KIND = DP) :: dw
!! Frequency intervals
REAL(KIND = DP) :: ww(nw_specfun)
!! Current frequency
REAL(KIND = DP) :: aux(2 * nbndfst * nktotf * nw_specfun + 2)
!! Vector to store the array
!
IF (mpime == ionode_id) THEN
!
! energy range and spacing for spectral function
!
dw = (wmax_specfun - wmin_specfun) / DBLE(nw_specfun - 1.d0)
DO iw = 1, nw_specfun
ww(iw) = wmin_specfun + DBLE(iw - 1) * dw
ENDDO
!
lesigma_all = 2 * nbndfst * nktotf * nw_specfun + 2
! First element is the current q-point
aux(1) = REAL(iqq - 1, KIND = DP) ! we need to start at the next q
! Second element is the total number of q-points
aux(2) = REAL(totq, KIND = DP)
!
i = 2
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO iw = 1, nw_specfun
i = i + 1
aux(i) = esigmar_all(ibnd, ik, iw)
ENDDO
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO iw = 1, nw_specfun
i = i + 1
aux(i) = esigmai_all(ibnd, ik, iw)
ENDDO
ENDDO
ENDDO
CALL diropn(iufilesigma_all, 'esigma_restart', lesigma_all, exst)
CALL davcio(aux, lesigma_all, iufilesigma_all, 1, +1)
CLOSE(iufilesigma_all)
ENDIF
!
! Make everythin 0 except the range of k-points we are working on
IF (lower_bnd > 1) THEN
esigmar_all(:, 1:lower_bnd - 1, :) = zero
esigmai_all(:, 1:lower_bnd - 1, :) = zero
ENDIF
IF (upper_bnd < nktotf) THEN
esigmar_all(:, upper_bnd + 1:nktotf, :) = zero
esigmai_all(:, upper_bnd + 1:nktotf, :) = zero
ENDIF
!
!----------------------------------------------------------------------------
END SUBROUTINE spectral_write
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE spectral_read(iqq, totq, nktotf, esigmar_all, esigmai_all)
!----------------------------------------------------------------------------
!!
!! Self-energy reading
!!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE elph2, ONLY : lower_bnd, upper_bnd, nbndfst
USE epwcom, ONLY : wmin_specfun, wmax_specfun, nw_specfun
USE io_var, ONLY : iufilesigma_all
USE io_files, ONLY : prefix, tmp_dir, diropn
USE constants_epw, ONLY : zero
USE mp, ONLY : mp_barrier, mp_bcast
USE mp_world, ONLY : mpime, world_comm
USE io_global, ONLY : ionode_id
!
IMPLICIT NONE
!
INTEGER, INTENT(inout) :: iqq
!! Current q-point
INTEGER, INTENT(in) :: totq
!! Total number of q-points
INTEGER, INTENT(in) :: nktotf
!! Total number of k-points
REAL(KIND = DP), INTENT(out) :: esigmar_all(nbndfst, nktotf, nw_specfun)
!! Real part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(out) :: esigmai_all(nbndfst, nktotf, nw_specfun)
!! Imaginary part of the electron-phonon self-energy accross all pools
!
! Local variables
LOGICAL :: exst
!! Does the file exist
INTEGER :: i
!! Iterative index
INTEGER :: ik
!! K-point index
INTEGER :: ibnd
!! Local band index
INTEGER :: iw
!! Counter on the frequency
INTEGER :: lesigma_all
!! Length of the vector
INTEGER :: nqtotf_read
!! Total number of q-point read
REAL(KIND = DP) :: dw
!! Frequency intervals
REAL(KIND = DP) :: ww(nw_specfun)
!! Current frequency
REAL(KIND = DP) :: aux(2 * nbndfst * nktotf * nw_specfun + 2)
!! Vector to store the array
!
CHARACTER(LEN = 256) :: name1
!
!
IF (mpime == ionode_id) THEN
!
! First inquire if the file exists
#if defined(__MPI)
name1 = TRIM(tmp_dir) // TRIM(prefix) // '.esigma_restart1'
#else
name1 = TRIM(tmp_dir) // TRIM(prefix) // '.esigma_restart'
#endif
INQUIRE(FILE = name1, EXIST = exst)
!
IF (exst) THEN ! read the file
!
lesigma_all = 2 * nbndfst * nktotf * nw_specfun + 2
CALL diropn(iufilesigma_all, 'esigma_restart', lesigma_all, exst)
CALL davcio(aux, lesigma_all, iufilesigma_all, 1, -1)
!
! First element is the iteration number
iqq = INT(aux(1))
iqq = iqq + 1 ! we need to start at the next q
nqtotf_read = INT(aux(2))
IF (nqtotf_read /= totq) CALL errore('electron_read',&
&'Error: The current total number of q-point is not the same as the read one. ', 1)
!
i = 2
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO iw = 1, nw_specfun
i = i + 1
esigmar_all(ibnd, ik, iw) = aux(i)
ENDDO
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO iw = 1, nw_specfun
i = i + 1
esigmai_all(ibnd, ik, iw) = aux(i)
ENDDO
ENDDO
ENDDO
CLOSE(iufilesigma_all)
ENDIF
ENDIF
!
CALL mp_bcast(exst, ionode_id, world_comm)
!
IF (exst) THEN
CALL mp_bcast(iqq, ionode_id, world_comm)
CALL mp_bcast(esigmar_all, ionode_id, world_comm)
CALL mp_bcast(esigmai_all, ionode_id, world_comm)
!
! Make everythin 0 except the range of k-points we are working on
IF (lower_bnd > 1) THEN
esigmar_all(:, 1:lower_bnd - 1, :) = zero
esigmai_all(:, 1:lower_bnd - 1, :) = zero
ENDIF
IF (upper_bnd < nktotf) THEN
esigmar_all(:, upper_bnd + 1:nktotf, :) = zero
esigmai_all(:, upper_bnd + 1:nktotf, :) = zero
ENDIF
!
WRITE(stdout, '(a,i10,a,i10)' ) ' Restart from: ', iqq,'/', totq
ENDIF
!
!----------------------------------------------------------------------------
END SUBROUTINE spectral_read
!----------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
END MODULE io_selfen
!------------------------------------------------------------------------------

439
EPW/src/io_transport.f90
View File

@ -209,7 +209,7 @@
WRITE(stdout, '(/5x,a)') REPEAT('=',67)
WRITE(stdout, '(5x,"Scattering rate for IBTE")')
WRITE(stdout, '(5x,a/)') REPEAT('=',67)
WRITE(stdout, '(5x,"Restart and restart_freq inputs deactivated (restart point at every q-points).")')
WRITE(stdout, '(5x,"Restart and restart_step inputs deactivated (restart point at every q-points).")')
WRITE(stdout, '(5x,"No intermediate mobility will be shown.")')
!
IF (fsthick < 1.d3) THEN
@ -1602,443 +1602,6 @@
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE electron_write(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
!----------------------------------------------------------------------------
!!
!! Write self-energy
!!
USE kinds, ONLY : DP
USE elph2, ONLY : lower_bnd, upper_bnd, nbndfst
USE io_var, ONLY : iufilsigma_all
USE io_files, ONLY : diropn
USE constants_epw, ONLY : zero
USE mp, ONLY : mp_barrier
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: iqq
!! Current q-point
INTEGER, INTENT(in) :: totq
!! Total number of q-points
INTEGER, INTENT(in) :: nktotf
!! Total number of k-points
REAL(KIND = DP), INTENT(inout) :: sigmar_all(nbndfst, nktotf)
!! Real part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(inout) :: sigmai_all(nbndfst, nktotf)
!! Imaginary part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(inout) :: zi_all(nbndfst, nktotf)
!! Z parameter of electron-phonon self-energy accross all pools
!
! Local variables
LOGICAL :: exst
!! Does the file exist
INTEGER :: i
!! Iterative index
INTEGER :: ik
!! K-point index
INTEGER :: ibnd
!! Local band index
INTEGER :: lsigma_all
!! Length of the vector
REAL(KIND = DP) :: aux(3 * nbndfst * nktotf + 2)
!! Vector to store the array
!
IF (mpime == ionode_id) THEN
!
lsigma_all = 3 * nbndfst * nktotf + 2
! First element is the current q-point
aux(1) = REAL(iqq - 1, KIND = DP) ! we need to start at the next q
! Second element is the total number of q-points
aux(2) = REAL(totq, KIND = DP)
!
i = 2
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
aux(i) = sigmar_all(ibnd, ik)
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
aux(i) = sigmai_all(ibnd, ik)
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
aux(i) = zi_all(ibnd, ik)
ENDDO
ENDDO
CALL diropn(iufilsigma_all, 'sigma_restart', lsigma_all, exst)
CALL davcio(aux, lsigma_all, iufilsigma_all, 1, +1)
CLOSE(iufilsigma_all)
ENDIF
!
! Make everythin 0 except the range of k-points we are working on
IF (lower_bnd > 1) THEN
sigmar_all(:, 1:lower_bnd - 1) = zero
sigmai_all(:, 1:lower_bnd - 1) = zero
zi_all(:, 1:lower_bnd - 1) = zero
ENDIF
IF (upper_bnd < nktotf) THEN
sigmar_all(:, upper_bnd + 1:nktotf) = zero
sigmai_all(:, upper_bnd + 1:nktotf) = zero
zi_all(:, upper_bnd + 1:nktotf) = zero
ENDIF
!
!----------------------------------------------------------------------------
END SUBROUTINE electron_write
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE electron_read(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
!----------------------------------------------------------------------------
!!
!! Self-energy reading
!!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE elph2, ONLY : lower_bnd, upper_bnd, nbndfst
USE io_var, ONLY : iufilsigma_all
USE io_files, ONLY : prefix, tmp_dir, diropn
USE constants_epw, ONLY : zero
USE mp, ONLY : mp_barrier, mp_bcast
USE mp_world, ONLY : mpime, world_comm
USE io_global, ONLY : ionode_id
!
IMPLICIT NONE
!
INTEGER, INTENT(inout) :: iqq
!! Current q-point
INTEGER, INTENT(in) :: totq
!! Total number of q-points
INTEGER, INTENT(in) :: nktotf
!! Total number of k-points
REAL(KIND = DP), INTENT(out) :: sigmar_all(nbndfst, nktotf)
!! Real part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(out) :: sigmai_all(nbndfst, nktotf)
!! Imaginary part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(out) :: zi_all(nbndfst, nktotf)
!! Z parameter of electron-phonon self-energy accross all pools
!
! Local variables
LOGICAL :: exst
!! Does the file exist
INTEGER :: i
!! Iterative index
INTEGER :: ik
!! K-point index
INTEGER :: ibnd
!! Local band index
INTEGER :: lsigma_all
!! Length of the vector
INTEGER :: nqtotf_read
!! Total number of q-point read
REAL(KIND = DP) :: aux(3 * nbndfst * nktotf + 2)
!! Vector to store the array
!
CHARACTER(LEN = 256) :: name1
!
IF (mpime == ionode_id) THEN
!
! First inquire if the file exists
#if defined(__MPI)
name1 = TRIM(tmp_dir) // TRIM(prefix) // '.sigma_restart1'
#else
name1 = TRIM(tmp_dir) // TRIM(prefix) // '.sigma_restart'
#endif
INQUIRE(FILE = name1, EXIST = exst)
!
IF (exst) THEN ! read the file
!
lsigma_all = 3 * nbndfst * nktotf + 2
CALL diropn(iufilsigma_all, 'sigma_restart', lsigma_all, exst)
CALL davcio(aux, lsigma_all, iufilsigma_all, 1, -1)
!
! First element is the iteration number
iqq = INT(aux(1))
iqq = iqq + 1 ! we need to start at the next q
nqtotf_read = INT(aux(2))
IF (nqtotf_read /= totq) CALL errore('electron_read',&
&'Error: The current total number of q-point is not the same as the read one. ', 1)
!
i = 2
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
sigmar_all(ibnd, ik) = aux(i)
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
sigmai_all(ibnd, ik) = aux(i)
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
i = i + 1
zi_all(ibnd, ik) = aux(i)
ENDDO
ENDDO
CLOSE(iufilsigma_all)
ENDIF
ENDIF
!
CALL mp_bcast(exst, ionode_id, world_comm)
!
IF (exst) THEN
CALL mp_bcast(iqq, ionode_id, world_comm)
CALL mp_bcast(sigmar_all, ionode_id, world_comm)
CALL mp_bcast(sigmai_all, ionode_id, world_comm)
CALL mp_bcast(zi_all, ionode_id, world_comm)
!
! Make everythin 0 except the range of k-points we are working on
IF (lower_bnd > 1) THEN
sigmar_all(:, 1:lower_bnd - 1) = zero
sigmai_all(:, 1:lower_bnd - 1) = zero
zi_all(:, 1:lower_bnd - 1) = zero
ENDIF
IF (upper_bnd < nktotf) THEN
sigmar_all(:, upper_bnd + 1:nktotf) = zero
sigmai_all(:, upper_bnd + 1:nktotf) = zero
zi_all(:, upper_bnd + 1:nktotf) = zero
ENDIF
!
WRITE(stdout, '(a,i10,a,i10)' ) ' Restart from: ', iqq,'/', totq
ENDIF
!
!----------------------------------------------------------------------------
END SUBROUTINE electron_read
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE spectral_write(iqq, totq, nktotf, esigmar_all, esigmai_all)
!----------------------------------------------------------------------------
!!
!! Write self-energy
!!
USE kinds, ONLY : DP
USE elph2, ONLY : lower_bnd, upper_bnd, nbndfst
USE epwcom, ONLY : wmin_specfun, wmax_specfun, nw_specfun
USE io_var, ONLY : iufilesigma_all
USE io_files, ONLY : diropn
USE constants_epw, ONLY : zero
USE mp, ONLY : mp_barrier
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: iqq
!! Current q-point
INTEGER, INTENT(in) :: totq
!! Total number of q-points
INTEGER, INTENT(in) :: nktotf
!! Total number of k-points
REAL(KIND = DP), INTENT(inout) :: esigmar_all(nbndfst, nktotf, nw_specfun)
!! Real part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(inout) :: esigmai_all(nbndfst, nktotf, nw_specfun)
!! Imaginary part of the electron-phonon self-energy accross all pools
!
! Local variables
LOGICAL :: exst
!! Does the file exist
INTEGER :: i
!! Iterative index
INTEGER :: ik
!! K-point index
INTEGER :: ibnd
!! Local band index
INTEGER :: iw
!! Counter on the frequency
INTEGER :: lesigma_all
!! Length of the vector
REAL(KIND = DP) :: dw
!! Frequency intervals
REAL(KIND = DP) :: ww(nw_specfun)
!! Current frequency
REAL(KIND = DP) :: aux(2 * nbndfst * nktotf * nw_specfun + 2)
!! Vector to store the array
!
IF (mpime == ionode_id) THEN
!
! energy range and spacing for spectral function
!
dw = (wmax_specfun - wmin_specfun) / DBLE(nw_specfun - 1.d0)
DO iw = 1, nw_specfun
ww(iw) = wmin_specfun + DBLE(iw - 1) * dw
ENDDO
!
lesigma_all = 2 * nbndfst * nktotf * nw_specfun + 2
! First element is the current q-point
aux(1) = REAL(iqq - 1, KIND = DP) ! we need to start at the next q
! Second element is the total number of q-points
aux(2) = REAL(totq, KIND = DP)
!
i = 2
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO iw = 1, nw_specfun
i = i + 1
aux(i) = esigmar_all(ibnd, ik, iw)
ENDDO
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO iw = 1, nw_specfun
i = i + 1
aux(i) = esigmai_all(ibnd, ik, iw)
ENDDO
ENDDO
ENDDO
CALL diropn(iufilesigma_all, 'esigma_restart', lesigma_all, exst)
CALL davcio(aux, lesigma_all, iufilesigma_all, 1, +1)
CLOSE(iufilesigma_all)
ENDIF
!
! Make everythin 0 except the range of k-points we are working on
IF (lower_bnd > 1) THEN
esigmar_all(:, 1:lower_bnd - 1, :) = zero
esigmai_all(:, 1:lower_bnd - 1, :) = zero
ENDIF
IF (upper_bnd < nktotf) THEN
esigmar_all(:, upper_bnd + 1:nktotf, :) = zero
esigmai_all(:, upper_bnd + 1:nktotf, :) = zero
ENDIF
!
!----------------------------------------------------------------------------
END SUBROUTINE spectral_write
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE spectral_read(iqq, totq, nktotf, esigmar_all, esigmai_all)
!----------------------------------------------------------------------------
!!
!! Self-energy reading
!!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE elph2, ONLY : lower_bnd, upper_bnd, nbndfst
USE epwcom, ONLY : wmin_specfun, wmax_specfun, nw_specfun
USE io_var, ONLY : iufilesigma_all
USE io_files, ONLY : prefix, tmp_dir, diropn
USE constants_epw, ONLY : zero
USE mp, ONLY : mp_barrier, mp_bcast
USE mp_world, ONLY : mpime, world_comm
USE io_global, ONLY : ionode_id
!
IMPLICIT NONE
!
INTEGER, INTENT(inout) :: iqq
!! Current q-point
INTEGER, INTENT(in) :: totq
!! Total number of q-points
INTEGER, INTENT(in) :: nktotf
!! Total number of k-points
REAL(KIND = DP), INTENT(out) :: esigmar_all(nbndfst, nktotf, nw_specfun)
!! Real part of the electron-phonon self-energy accross all pools
REAL(KIND = DP), INTENT(out) :: esigmai_all(nbndfst, nktotf, nw_specfun)
!! Imaginary part of the electron-phonon self-energy accross all pools
!
! Local variables
LOGICAL :: exst
!! Does the file exist
INTEGER :: i
!! Iterative index
INTEGER :: ik
!! K-point index
INTEGER :: ibnd
!! Local band index
INTEGER :: iw
!! Counter on the frequency
INTEGER :: lesigma_all
!! Length of the vector
INTEGER :: nqtotf_read
!! Total number of q-point read
REAL(KIND = DP) :: dw
!! Frequency intervals
REAL(KIND = DP) :: ww(nw_specfun)
!! Current frequency
REAL(KIND = DP) :: aux(2 * nbndfst * nktotf * nw_specfun + 2)
!! Vector to store the array
!
CHARACTER(LEN = 256) :: name1
!
!
IF (mpime == ionode_id) THEN
!
! First inquire if the file exists
#if defined(__MPI)
name1 = TRIM(tmp_dir) // TRIM(prefix) // '.esigma_restart1'
#else
name1 = TRIM(tmp_dir) // TRIM(prefix) // '.esigma_restart'
#endif
INQUIRE(FILE = name1, EXIST = exst)
!
IF (exst) THEN ! read the file
!
lesigma_all = 2 * nbndfst * nktotf * nw_specfun + 2
CALL diropn(iufilesigma_all, 'esigma_restart', lesigma_all, exst)
CALL davcio(aux, lesigma_all, iufilesigma_all, 1, -1)
!
! First element is the iteration number
iqq = INT(aux(1))
iqq = iqq + 1 ! we need to start at the next q
nqtotf_read = INT(aux(2))
IF (nqtotf_read /= totq) CALL errore('electron_read',&
&'Error: The current total number of q-point is not the same as the read one. ', 1)
!
i = 2
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO iw = 1, nw_specfun
i = i + 1
esigmar_all(ibnd, ik, iw) = aux(i)
ENDDO
ENDDO
ENDDO
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO iw = 1, nw_specfun
i = i + 1
esigmai_all(ibnd, ik, iw) = aux(i)
ENDDO
ENDDO
ENDDO
CLOSE(iufilesigma_all)
ENDIF
ENDIF
!
CALL mp_bcast(exst, ionode_id, world_comm)
!
IF (exst) THEN
CALL mp_bcast(iqq, ionode_id, world_comm)
CALL mp_bcast(esigmar_all, ionode_id, world_comm)
CALL mp_bcast(esigmai_all, ionode_id, world_comm)
!
! Make everythin 0 except the range of k-points we are working on
IF (lower_bnd > 1) THEN
esigmar_all(:, 1:lower_bnd - 1, :) = zero
esigmai_all(:, 1:lower_bnd - 1, :) = zero
ENDIF
IF (upper_bnd < nktotf) THEN
esigmar_all(:, upper_bnd + 1:nktotf, :) = zero
esigmai_all(:, upper_bnd + 1:nktotf, :) = zero
ENDIF
!
WRITE(stdout, '(a,i10,a,i10)' ) ' Restart from: ', iqq,'/', totq
ENDIF
!
!----------------------------------------------------------------------------
END SUBROUTINE spectral_read
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE tau_write(iqq, totq, nktotf, second)
!----------------------------------------------------------------------------
USE kinds, ONLY : DP

10
EPW/src/low_lvl.f90
View File

@ -956,7 +956,7 @@
!----------------------------------------------------------------------
!
!---------------------------------------------------------------------------
SUBROUTINE read_modes(iunpun, current_iq, ierr)
SUBROUTINE read_disp_pattern(iunpun, current_iq, ierr)
!---------------------------------------------------------------------------
!!
!! This routine reads the displacement patterns.
@ -995,7 +995,7 @@
CALL iotk_scan_dat(iunpun, "QPOINT_NUMBER", iq)
ENDIF
CALL mp_bcast(iq, meta_ionode_id, world_comm)
IF (iq /= current_iq) CALL errore('read_modes', ' Problems with current_iq', 1)
IF (iq /= current_iq) CALL errore('read_disp_pattern', ' Problems with current_iq', 1)
!
IF (meta_ionode) THEN
!
@ -1029,7 +1029,7 @@
RETURN
!
!---------------------------------------------------------------------------
END SUBROUTINE read_modes
END SUBROUTINE read_disp_pattern
!---------------------------------------------------------------------------
!
!------------------------------------------------------------
@ -1183,7 +1183,7 @@
!
! This is only a quick fix since the routine was written for parallel
! execution - FG June 2014
#if ! defined(__MPI)
#if defined(__MPI)
my_pool_id = 0
#endif
!
@ -1250,7 +1250,7 @@
!
! This is only a quick fix since the routine was written for parallel
! execution - FG June 2014
#if ! defined(__MPI)
#if defined(__MPI)
my_pool_id = 0
#endif
!

696
EPW/src/plot.f90
View File

@ -1,696 +0,0 @@
!
! Copyright (C) 2016-2019 Samuel Ponce', Roxana Margine, Feliciano Giustino
! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, 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 .
!
!----------------------------------------------------------------------
MODULE plot
!----------------------------------------------------------------------
!!
!! This module contains routine to plot data as well as DOS-like quantities.
!!
IMPLICIT NONE
!
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE a2f_main()
!-----------------------------------------------------------------------
!!
!! Compute the Eliasberg spectral function
!! in the Migdal approximation.
!!
!! If the q-points are not on a uniform grid (i.e. a line)
!! the function will not be correct
!!
!! 02/2009 works in serial on ionode at the moment. can be parallelized
!! 03/2009 added transport spectral function -- this involves a v_k dot v_kq term
!! in the quantities coming from selfen_phon.f90. Not fully implemented
!! 10/2009 the code is transitioning to 'on-the-fly' phonon selfenergies
!! and this routine is not currently functional
!! 10/2015 RM: added calcution of Tc based on Allen-Dynes formula
!! 09/2019 SP: Cleaning
!!
!
USE kinds, ONLY : DP
USE phcom, ONLY : nmodes
USE cell_base, ONLY : omega
USE epwcom, ONLY : degaussq, delta_qsmear, nqsmear, nqstep, nsmear, eps_acustic, &
delta_smear, degaussw, fsthick, nc
USE elph2, ONLY : nqtotf, wf, wqf, lambda_all, lambda_v_all
USE constants_epw, ONLY : ryd2mev, ryd2ev, kelvin2eV, one, two, zero, kelvin2Ry, pi
USE mp, ONLY : mp_barrier, mp_sum
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
USE io_global, ONLY : stdout
USE io_var, ONLY : iua2ffil, iudosfil, iua2ftrfil, iures
USE io_files, ONLY : prefix
!
IMPLICIT NONE
!
CHARACTER(LEN = 256) :: fila2f
!! File name for Eliashberg spectral function
CHARACTER(LEN = 256) :: fila2ftr
!! File name for transport Eliashberg spectral function
CHARACTER(LEN = 256) :: fildos
!! File name for phonon density of states
CHARACTER(LEN = 256) :: filres
!! File name for resistivity
!
INTEGER :: imode
!! Counter on mode
INTEGER :: iq
!! Counter on the q-point index
INTEGER :: iw
!! Counter on the frequency
INTEGER :: ismear
!! Counter on smearing values (phonons)
INTEGER :: isig
!! Counter on smearing values (electrons)
INTEGER :: i
!! Counter on mu
INTEGER :: itemp
!! Counter on temperature
INTEGER :: ierr
!! Error status
!
REAL(KIND = DP) :: weight
!! Factor in a2f
REAL(KIND = DP) :: temp
!! Temperature
REAL(KIND = DP) :: n
!! Carrier density
REAL(KIND = DP) :: be
!! Bose-Einstein distribution
REAL(KIND = DP) :: prefact
!! Prefactor in resistivity
REAL(KIND = DP) :: lambda_tot
!! Total e-ph coupling strength (summation)
REAL(KIND = DP) :: lambda_tr_tot
!! Total transport e-ph coupling strength (summation)
REAL(KIND = DP) :: degaussq0
!! Phonon smearing
REAL(KIND = DP) :: inv_degaussq0
!! Inverse of the smearing for efficiency reasons
REAL(KIND = DP) :: a2f_tmp
!! Temporary variable for Eliashberg spectral function
REAL(KIND = DP) :: a2f_tr_tmp
!! Temporary variable for transport Eliashberg spectral function
REAL(KIND = DP) :: om_max
!! max phonon frequency increased by 10%
REAL(KIND = DP) :: dw
!! Frequency intervals
REAL(KIND = DP) :: w0
!! Current frequency w(imode, iq)
REAL(KIND = DP) :: l
!! Temporary variable for e-ph coupling strength
REAL(KIND = DP) :: l_tr
!! Temporary variable for transport e-ph coupling strength
REAL(KIND = DP) :: tc
!! Critical temperature
REAL(KIND = DP) :: mu
!! Coulomb pseudopotential
REAL(KIND = DP), EXTERNAL :: w0gauss
!! The derivative of wgauss: an approximation to the delta function
REAL(KIND = DP) :: ww(nqstep)
!! Current frequency
REAL(KIND = DP), ALLOCATABLE :: a2f_(:, :)
!! Eliashberg spectral function for different ismear
REAL(KIND = DP), ALLOCATABLE :: a2f_tr(:, :)
!! Transport Eliashberg spectral function for different ismear
REAL(KIND = DP), ALLOCATABLE :: l_a2f(:)
!! total e-ph coupling strength (a2f_ integration) for different ismear
REAL(KIND = DP), ALLOCATABLE :: l_a2f_tr(:)
!! total transport e-ph coupling strength (a2f_tr integration) for different ismear
REAL(KIND = DP), ALLOCATABLE :: dosph(:, :)
!! Phonon density of states for different for different ismear
REAL(KIND = DP), ALLOCATABLE :: logavg(:)
!! logavg phonon frequency for different ismear
REAL(KIND = DP), ALLOCATABLE :: rho(:, :)
!! Resistivity for different for different ismear
!
CALL start_clock('a2F')
IF (mpime == ionode_id) THEN
!
ALLOCATE(a2f_(nqstep, nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating a2f_', 1)
ALLOCATE(a2f_tr(nqstep, nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating a2f_tr', 1)
ALLOCATE(dosph(nqstep, nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating dosph', 1)
ALLOCATE(l_a2f(nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating l_a2f', 1)
ALLOCATE(l_a2f_tr(nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating l_a2f_tr', 1)
ALLOCATE(logavg(nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating logavg', 1)
! The resitivity is computed for temperature between 0K-1000K by step of 10
! This is hardcoded and needs to be changed here if one wants to modify it
ALLOCATE(rho(100, nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating rho', 1)
!
DO isig = 1, nsmear
!
IF (isig < 10) THEN
WRITE(fila2f, '(a, a6, i1)') TRIM(prefix), '.a2f.0', isig
WRITE(fila2ftr, '(a, a9, i1)') TRIM(prefix), '.a2f_tr.0', isig
WRITE(filres, '(a, a6, i1)') TRIM(prefix), '.res.0', isig
WRITE(fildos, '(a, a8, i1)') TRIM(prefix), '.phdos.0', isig
ELSE
WRITE(fila2f, '(a, a5, i2)') TRIM(prefix), '.a2f.', isig
WRITE(fila2ftr, '(a, a8, i2)') TRIM(prefix), '.a2f_tr.', isig
WRITE(filres, '(a, a5, i2)') TRIM(prefix), '.res.', isig
WRITE(fildos, '(a, a7, i2)') TRIM(prefix), '.phdos.', isig
ENDIF
OPEN(UNIT = iua2ffil, FILE = fila2f, FORM = 'formatted')
OPEN(UNIT = iua2ftrfil, FILE = fila2ftr, FORM = 'formatted')
OPEN(UNIT = iures, FILE = filres, FORM = 'formatted')
OPEN(UNIT = iudosfil, FILE = fildos, FORM = 'formatted')
!
WRITE(stdout, '(/5x, a)') REPEAT('=',67)
WRITE(stdout, '(5x, "Eliashberg Spectral Function in the Migdal Approximation")')
WRITE(stdout, '(5x, a/)') REPEAT('=',67)
!
om_max = 1.1d0 * MAXVAL(wf(:, :)) ! increase by 10%
dw = om_max / DBLE(nqstep)
DO iw = 1, nqstep !
ww(iw) = DBLE(iw) * dw
ENDDO
!
lambda_tot = zero
l_a2f(:) = zero
a2f_(:, :) = zero
lambda_tr_tot = zero
l_a2f_tr(:) = zero
a2f_tr(:, :) = zero
dosph(:, :) = zero
logavg(:) = zero
!
DO ismear = 1, nqsmear
!
degaussq0 = degaussq + (ismear - 1) * delta_qsmear
inv_degaussq0 = one / degaussq0
!
DO iw = 1, nqstep ! loop over points on the a2F(w) graph
!
DO iq = 1, nqtotf ! loop over q-points
DO imode = 1, nmodes ! loop over modes
w0 = wf(imode, iq)
!
IF (w0 > eps_acustic) THEN
!
l = lambda_all(imode, iq, isig)
IF (lambda_all(imode, iq, isig) < 0.d0) l = zero ! sanity check
!
a2f_tmp = wqf(iq) * w0 * l / two
!
weight = w0gauss((ww(iw) - w0) * inv_degaussq0, 0) * inv_degaussq0
a2f_(iw, ismear) = a2f_(iw, ismear) + a2f_tmp * weight
dosph(iw, ismear) = dosph(iw, ismear) + wqf(iq) * weight
!
l_tr = lambda_v_all(imode, iq, isig)
IF (lambda_v_all(imode, iq, isig) < 0.d0) l_tr = zero !sanity check
!
a2f_tr_tmp = wqf(iq) * w0 * l_tr / two
!
a2f_tr(iw, ismear) = a2f_tr(iw, ismear) + a2f_tr_tmp * weight
!
ENDIF
ENDDO
ENDDO
!
! output a2f
!
IF (ismear == nqsmear) WRITE(iua2ffil, '(f12.7, 15f12.7)') ww(iw) * ryd2mev, a2f_(iw, :)
IF (ismear == nqsmear) WRITE(iua2ftrfil, '(f12.7, 15f12.7)') ww(iw) * ryd2mev, a2f_tr(iw, :)
IF (ismear == nqsmear) WRITE(iudosfil, '(f12.7, 15f12.7)') ww(iw) * ryd2mev, dosph(iw, :) / ryd2mev
!
! do the integral 2 int (a2F(w)/w dw)
!
l_a2f(ismear) = l_a2f(ismear) + two * a2f_(iw, ismear) / ww(iw) * dw
l_a2f_tr(ismear) = l_a2f_tr(ismear) + two * a2f_tr(iw, ismear) / ww(iw) * dw
logavg(ismear) = logavg(ismear) + two * a2f_(iw, ismear) * LOG(ww(iw)) / ww(iw) * dw
!
ENDDO
!
logavg(ismear) = EXP(logavg(ismear) / l_a2f(ismear))
!
ENDDO
!
DO iq = 1, nqtotf ! loop over q-points
DO imode = 1, nmodes ! loop over modes
IF (lambda_all(imode, iq, isig) > 0.d0 .AND. wf(imode, iq) > eps_acustic ) &
lambda_tot = lambda_tot + wqf(iq) * lambda_all(imode, iq, isig)
IF (lambda_v_all(imode, iq, isig) > 0.d0 .AND. wf(imode, iq) > eps_acustic) &
lambda_tr_tot = lambda_tr_tot + wqf(iq) * lambda_v_all(imode, iq, isig)
ENDDO
ENDDO
WRITE(stdout, '(5x, a, f12.7)') "lambda : ", lambda_tot
WRITE(stdout, '(5x, a, f12.7)') "lambda_tr : ", lambda_tr_tot
WRITE(stdout, '(a)') " "
!
!
! Allen-Dynes estimate of Tc for ismear = 1
!
WRITE(stdout, '(5x, a, f12.7, a)') "Estimated Allen-Dynes Tc"
WRITE(stdout, '(a)') " "
WRITE(stdout, '(5x, a, f12.7, a, f12.7)') "logavg = ", logavg(1), " l_a2f = ", l_a2f(1)
DO i = 1, 6
!
mu = 0.1d0 + 0.02d0 * DBLE(i - 1)
tc = logavg(1) / 1.2d0 * EXP(-1.04d0 * (1.d0 + l_a2f(1)) / (l_a2f(1) - mu * ( 1.d0 + 0.62d0 * l_a2f(1))))
! tc in K
!
tc = tc * ryd2ev / kelvin2eV
!SP: IF Tc is too big, it is not physical
IF (tc < 1000.0) THEN
WRITE(stdout, '(5x, a, f6.2, a, f22.12, a)') "mu = ", mu, " Tc = ", tc, " K"
ENDIF
!
ENDDO
!
rho(:, :) = zero
! Now compute the Resistivity of Metal using the Ziman formula
! rho(T,smearing) = 4 * pi * me/(n * e**2 * kb * T) int dw hbar w a2F_tr(w,smearing) n(w,T)(1+n(w,T))
! n is the number of electron per unit volume and n(w,T) is the Bose-Einstein distribution
! Usually this means "the number of electrons that contribute to the mobility" and so it is typically 8 (full shell)
! but not always. You might want to check this.
!
n = nc / omega
WRITE(iures, '(a)') '# Temperature [K] &
Resistivity [micro Ohm cm] for different Phonon smearing (meV) '
WRITE(iures, '("# ", 15f12.7)') ((degaussq + (ismear - 1) * delta_qsmear) * ryd2mev, ismear = 1, nqsmear)
DO ismear = 1, nqsmear
DO itemp = 1, 100 ! Per step of 10K
temp = itemp * 10.d0 * kelvin2Ry
! omega is the volume of the primitive cell in a.u.
!
prefact = 4.d0 * pi / (temp * n)
DO iw = 1, nqstep ! loop over points on the a2F(w)
!
be = one / (EXP(ww(iw) / temp) - one)
! Perform the integral with rectangle.
rho(itemp, ismear) = rho(itemp, ismear) + prefact * ww(iw) * a2f_tr(iw, ismear) * be * (1.d0 + be) * dw
!
ENDDO
! From a.u. to micro Ohm cm
! Conductivity 1 a.u. = 2.2999241E6 S/m
! Now to go from Ohm*m to micro Ohm cm we need to multiply by 1E8
rho(itemp, ismear) = rho(itemp, ismear) * 1E8 / 2.2999241E6
IF (ismear == nqsmear) WRITE (iures, '(i8, 15f12.7)') itemp * 10, rho(itemp, :)
ENDDO
ENDDO
CLOSE(iures)
!
WRITE(iua2ffil, *) "Integrated el-ph coupling"
WRITE(iua2ffil, '(" # ", 15f12.7)') l_a2f(:)
WRITE(iua2ffil, *) "Phonon smearing (meV)"
WRITE(iua2ffil, '(" # ", 15f12.7)') ((degaussq + (ismear - 1) * delta_qsmear) * ryd2mev, ismear = 1, nqsmear)
WRITE(iua2ffil, '(" Electron smearing (eV)", f12.7)') ((isig - 1) * delta_smear + degaussw) * ryd2ev
WRITE(iua2ffil, '(" Fermi window (eV)", f12.7)') fsthick * ryd2ev
WRITE(iua2ffil, '(" Summed el-ph coupling ", f12.7)') lambda_tot
CLOSE(iua2ffil)
!
WRITE(iua2ftrfil, *) "Integrated el-ph coupling"
WRITE(iua2ftrfil, '(" # ", 15f12.7)') l_a2f_tr(:)
WRITE(iua2ftrfil, *) "Phonon smearing (meV)"
WRITE(iua2ftrfil, '(" # ", 15f12.7)') ((degaussq + (ismear - 1) * delta_qsmear) * ryd2mev, ismear = 1, nqsmear)
WRITE(iua2ftrfil, '(" Electron smearing (eV)", f12.7)') ((isig - 1) * delta_smear + degaussw) * ryd2ev
WRITE(iua2ftrfil, '(" Fermi window (eV)", f12.7)') fsthick * ryd2ev
WRITE(iua2ftrfil, '(" Summed el-ph coupling ", f12.7)') lambda_tot
CLOSE(iua2ftrfil)
!
CLOSE(iudosfil)
!
ENDDO ! isig
!
DEALLOCATE(l_a2f, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating l_a2f', 1)
DEALLOCATE(l_a2f_tr, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating l_a2f_tr', 1)
DEALLOCATE(a2f_, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating a2f', 1)
DEALLOCATE(a2f_tr, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating a2f_tr', 1)
DEALLOCATE(rho, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating rho', 1)
DEALLOCATE(dosph, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating dosph', 1)
DEALLOCATE(logavg, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating logavg', 1)
!
ENDIF
!
CALL stop_clock('a2F')
CALL print_clock('a2F')
!
RETURN
!
!-----------------------------------------------------------------------
END SUBROUTINE a2f_main
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE nesting_fn_q(iqq, iq)
!-----------------------------------------------------------------------
!!
!! Compute the imaginary part of the phonon self energy due to electron-
!! phonon interaction in the Migdal approximation. This corresponds to
!! the phonon linewidth (half width). The phonon frequency is taken into
!! account in the energy selection rule.
!!
!! Use matrix elements, electronic eigenvalues and phonon frequencies
!! from ep-wannier interpolation.
!!
!-----------------------------------------------------------------------
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE epwcom, ONLY : nbndsub, fsthick, eptemp, ngaussw, degaussw, &
nsmear, delta_smear, efermi_read, fermi_energy
USE pwcom, ONLY : ef
USE elph2, ONLY : ibndmin, etf, wkf, xqf, wqf, nkqf, nktotf, &
nkf, xqf, nbndfst, efnew
USE constants_epw, ONLY : ryd2ev, zero, one, two
USE mp, ONLY : mp_barrier, mp_sum
USE mp_global, ONLY : inter_pool_comm
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: iqq
!! Current q-point index from selecq
INTEGER, INTENT(in) :: iq
!! Current q-point index
!
! Local variables
INTEGER :: ik
!! Counter on the k-point index
INTEGER :: ikk
!! k-point index
INTEGER :: ikq
!! q-point index
INTEGER :: ibnd
!! Counter on bands
INTEGER :: jbnd
!! Counter on bands
INTEGER :: fermicount
!! Number of states on the Fermi surface
INTEGER :: ismear
!! Counter on smearing values
!
REAL(KIND = DP) :: ekk
!! Eigen energy on the fine grid relative to the Fermi level
REAL(KIND = DP) :: ekq
!! Eigen energy of k+q on the fine grid relative to the Fermi level
REAL(KIND = DP) :: ef0
!! Fermi energy level
REAL(KIND = DP) :: weight
!! Imaginary part of the phonhon self-energy factor, sans e-ph matrix elements
REAL(KIND = DP) :: dosef
!! Density of state N(Ef)
REAL(KIND = DP) :: w0g1
!! Dirac delta at k for the imaginary part of $\Sigma$
REAL(KIND = DP) :: w0g2
!! Dirac delta at k+q for the imaginary part of $\Sigma$
REAL(KIND = DP) :: degaussw0
!! degaussw0 = (ismear-1) * delta_smear + degaussw
REAL(KIND = DP) :: inv_degaussw0
!! Inverse degaussw0 for efficiency reasons
REAL(KIND = DP) :: gamma
!! Nesting function
REAL(KIND = DP) :: dos_ef
!! Function returning the density of states at the Fermi level
REAL(KIND = DP) :: w0gauss
!! This function computes the derivative of the Fermi-Dirac function
!! It is therefore an approximation for a delta function
!
IF (iqq == 1) THEN
WRITE(stdout, '(/5x, a)') REPEAT('=', 67)
WRITE(stdout, '(5x, "Nesting Function in the double delta approx")')
WRITE(stdout, '(5x, a/)') REPEAT('=', 67)
!
IF (fsthick < 1.d3) WRITE(stdout, '(/5x, a, f10.6, a)' ) &
'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
WRITE(stdout, '(/5x, a, f10.6, a)' ) 'Golden Rule strictly enforced with T = ', eptemp * ryd2ev, ' eV'
ENDIF
!
! SP: The Gamma function needs to be put to 0 for each q
gamma = zero
!
! Here we loop on smearing values
DO ismear = 1, nsmear
!
degaussw0 = (ismear - 1) * delta_smear + degaussw
inv_degaussw0 = one / degaussw0
!
! Fermi level and corresponding DOS
!
! Note that the weights of k+q points must be set to zero here
! no spin-polarized calculation here
IF (efermi_read) THEN
ef0 = fermi_energy
ELSE
ef0 = efnew
ENDIF
!
dosef = dos_ef(ngaussw, degaussw0, ef0, etf, wkf, nkqf, nbndsub)
! N(Ef) in the equation for lambda is the DOS per spin
dosef = dosef / two
!
IF (iqq == 1) THEN
WRITE(stdout, 100) degaussw0 * ryd2ev, ngaussw
WRITE(stdout, 101) dosef / ryd2ev, ef0 * ryd2ev
ENDIF
!
!
CALL start_clock('nesting')
!
fermicount = 0
DO ik = 1, nkf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! here we must have ef, not ef0, to be consistent with ephwann_shuffle
IF ((MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) .AND. &
(MINVAL(ABS(etf(:, ikq) - ef)) < fsthick)) then
!
fermicount = fermicount + 1
!
DO ibnd = 1, nbndfst
!
ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0
w0g1 = w0gauss(ekk * inv_degaussw0, 0) * inv_degaussw0
!
DO jbnd = 1, nbndfst
!
ekq = etf(ibndmin - 1 + jbnd, ikq) - ef0
w0g2 = w0gauss(ekq *inv_degaussw0, 0) * inv_degaussw0
!
! = k-point weight * [f(E_k) - f(E_k+q)]/ [E_k+q - E_k -w_q +id]
! This is the imaginary part of the phonon self-energy, sans the matrix elements
!
! weight = wkf (ikk) * (wgkk - wgkq) * &
! aimag ( cone / ( ekq - ekk - ci * degaussw ) )
!
! the below expression is positive-definite, but also an approximation
! which neglects some fine features
!
weight = wkf(ikk) * w0g1 * w0g2
!
gamma = gamma + weight
!
ENDDO ! jbnd
ENDDO ! ibnd
ENDIF ! endif fsthick
ENDDO ! loop on k
!
! collect contributions from all pools (sum over k-points)
! this finishes the integral over the BZ (k)
!
CALL mp_sum(gamma, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
!
WRITE(stdout, '(/5x, "iq = ",i5," coord.: ", 3f9.5, " wt: ", f9.5)') iq, xqf(:, iq) , wqf(iq)
WRITE(stdout, '(5x, a)') REPEAT('-', 67)
!
WRITE(stdout, 102) gamma
WRITE(stdout, '(5x,a/)') REPEAT('-', 67)
!
WRITE(stdout, '(/5x, a, i8, a, i8/)') &
'Number of (k,k+q) pairs on the Fermi surface: ', fermicount, ' out of ', nktotf
!
CALL stop_clock('nesting')
ENDDO !smears
!
100 FORMAT(5x, 'Gaussian Broadening: ', f7.3,' eV, ngauss=', i4)
101 FORMAT(5x, 'DOS =', f10.6, ' states/spin/eV/Unit Cell at Ef=', f10.6, ' eV')
102 FORMAT(5x, 'Nesting function (q)=', E15.6, ' [Adimensional]')
!
!-----------------------------------------------------------------------
END SUBROUTINE nesting_fn_q
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE plot_band()
!-----------------------------------------------------------------------
!!
!! This routine writes output files for phonon dispersion and band structure
!! SP : Modified so that it works with the current plotband.x of QE 5
!!
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg
USE phcom, ONLY : nmodes
USE epwcom, ONLY : nbndsub, filqf, filkf
USE elph2, ONLY : etf, nkf, nqtotf, wf, xkf, xqf, nkqtotf, nktotf
USE constants_epw, ONLY : ryd2mev, ryd2ev, zero
USE io_var, ONLY : iufilfreq, iufileig
USE elph2, ONLY : nkqf
USE io_global, ONLY : ionode_id
USE mp, ONLY : mp_barrier, mp_sum
USE mp_global, ONLY : inter_pool_comm, my_pool_id
USE poolgathering, ONLY : poolgather2
!
IMPLICIT NONE
!
! Local variables
INTEGER :: ik
!! Global k-point index
INTEGER :: ikk
!! Index for the k-point
INTEGER :: ikq
!! Index for the q-point
INTEGER :: ibnd
!! Band index
INTEGER :: imode
!! Mode index
INTEGER :: iq
!! Global q-point index
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: dist
!! Distance from G-point
REAL(KIND = DP) :: dprev
!! Previous distance
REAL(KIND = DP) :: dcurr
!! Current distance
REAL(KIND = DP), ALLOCATABLE :: xkf_all(:, :)
!! K-points on the full k grid (all pools)
REAL(KIND = DP), ALLOCATABLE :: etf_all(:, :)
!! Eigenenergies on the full k grid (all pools)
!
IF (filqf /= ' ') THEN
!
IF (my_pool_id == ionode_id) THEN
!
OPEN(iufilfreq, FILE = "phband.freq", FORM = 'formatted')
WRITE(iufilfreq, '(" &plot nbnd=", i4, ", nks=", i6, " /")') nmodes, nqtotf
!
! crystal to cartesian coordinates
CALL cryst_to_cart(nqtotf, xqf, bg, 1)
!
dist = zero
dprev = zero
dcurr = zero
DO iq = 1, nqtotf
!
IF (iq /= 1) THEN
dist = DSQRT((xqf(1, iq) - xqf(1, iq - 1)) * (xqf(1, iq) - xqf(1, iq - 1)) &
+ (xqf(2, iq) - xqf(2, iq - 1)) * (xqf(2, iq) - xqf(2, iq - 1)) &
+ (xqf(3, iq) - xqf(3, iq - 1)) * (xqf(3, iq) - xqf(3, iq - 1)))
ELSE
dist = zero
ENDIF
dcurr = dprev + dist
dprev = dcurr
WRITE(iufilfreq, '(10x, 3f10.6)') xqf(:, iq)
WRITE(iufilfreq, '(1000f14.4)') (wf(imode, iq) * ryd2mev, imode = 1, nmodes)
!
ENDDO
CLOSE(iufilfreq)
!
! back from cartesian to crystal coordinates
CALL cryst_to_cart(nqtotf, xqf, at, -1)
!
ENDIF
ENDIF ! filqf
!
IF (filkf /= ' ') THEN
!
DO ik = 1, nkf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
ENDDO
!
ALLOCATE(xkf_all(3, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('plot_band', 'Error allocating xkf_all', 1)
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('plot_band', 'Error allocating etf_all', 1)
!
#if defined(__MPI)
CALL poolgather2(3, nkqtotf, nkqf, xkf, xkf_all)
CALL poolgather2(nbndsub, nkqtotf, nkqf, etf, etf_all)
CALL mp_barrier(inter_pool_comm)
#else
!
xkf_all = xkf
etf_all = etf
#endif
!
IF (my_pool_id == ionode_id) THEN
!
OPEN(iufileig, FILE = "band.eig", FORM = 'formatted')
WRITE(iufileig, '(" &plot nbnd=", i4, ", nks=", i6, " /")') nbndsub, nktotf
!
! crystal to cartesian coordinates
CALL cryst_to_cart(nkqtotf, xkf_all, bg, 1)
!
dist = zero
dprev = zero
dcurr = zero
DO ik = 1, nktotf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
IF (ikk /= 1) THEN
dist = DSQRT((xkf_all(1, ikk) - xkf_all(1, ikk - 2)) * (xkf_all(1, ikk) - xkf_all(1, ikk - 2)) &
+ (xkf_all(2, ikk) - xkf_all(2, ikk - 2)) * (xkf_all(2, ikk) - xkf_all(2, ikk - 2)) &
+ (xkf_all(3, ikk) - xkf_all(3, ikk - 2)) * (xkf_all(3, ikk) - xkf_all(3, ikk - 2)))
ELSE
dist = 0.d0
ENDIF
dcurr = dprev + dist
dprev = dcurr
WRITE(iufileig, '(10x, 3f10.6)') xkf_all(:, ikk)
WRITE(iufileig, '(1000f20.12)') (etf_all(ibnd, ikk) * ryd2ev, ibnd = 1, nbndsub)
!
ENDDO
CLOSE(iufileig)
!
! back from cartesian to crystal coordinates
CALL cryst_to_cart(nkqtotf, xkf_all, at, -1)
!
ENDIF
CALL mp_barrier(inter_pool_comm)
!
DEALLOCATE(xkf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('plot_band', 'Error deallocating xkf_all', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('plot_band', 'Error deallocating etf_all', 1)
!
ENDIF ! filkf
!
RETURN
!
!----------------------------------------------------------------------------
END SUBROUTINE plot_band
!----------------------------------------------------------------------------
!------------------------------------------------------------------------------
END MODULE plot
!------------------------------------------------------------------------------

175
EPW/src/printing.f90
View File

@ -313,7 +313,7 @@
IF (PRESENT(max_mob)) THEN
max_mob(:) = zero
ENDIF
CALL prtheader()
CALL prtheader_mob()
! compute conductivity
DO itemp = 1, nstemp
carrier_density = 0.0
@ -504,7 +504,7 @@
!! Carrier density [nb of carrier per unit cell]
REAL(KIND = DP) :: fnk
!! Fermi-Dirac occupation function
REAL(KIND = DP) :: Fi_check(3)
REAL(KIND = DP) :: fi_check(3)
!! Sum rule on population
REAL(KIND = DP), EXTERNAL :: wgauss
!! Compute the approximate theta function. Here computes Fermi-Dirac
@ -516,7 +516,7 @@
IF (PRESENT(max_mob)) THEN
max_mob(:) = zero
ENDIF
CALL prtheader()
CALL prtheader_mob()
DO itemp = 1, nstemp
carrier_density = 0.0
etemp = transp_temp(itemp)
@ -610,7 +610,7 @@
END SUBROUTINE compute_sigma
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
SUBROUTINE prtmob(itemp, sigma, carrier_density, Fi_check, ef0, etemp, max_mob)
SUBROUTINE prtmob(itemp, sigma, carrier_density, fi_check, ef0, etemp, max_mob)
!-----------------------------------------------------------------------
!!
!! This routine print the mobility (or conducrtivity for metals) in a
@ -632,7 +632,7 @@
!! Conductivity tensor
REAL(KIND = DP), INTENT(in) :: carrier_density
!! Carrier density in a.u.
REAL(KIND = DP), INTENT(in) :: Fi_check(3)
REAL(KIND = DP), INTENT(in) :: fi_check(3)
!! Integrated population vector
REAL(KIND = DP), INTENT(in) :: ef0
!! Fermi-level
@ -658,10 +658,10 @@
IF (ABS(nden) < eps80) CALL errore('prtmob', 'The carrier density is 0', 1)
mobility(:, :) = mobility(:, :) / (electron_SI * carrier_density * inv_cell) * (bohr2ang * ang2cm)**3
WRITE(stdout, '(5x, 1f8.3, 1f9.4, 1E14.5, 1E14.5, 3E16.6)') etemp * ryd2ev / kelvin2eV, ef0 * ryd2ev, &
nden, SUM(Fi_check(:)), mobility(1, 1), mobility(1, 2), mobility(1, 3)
nden, SUM(fi_check(:)), mobility(1, 1), mobility(1, 2), mobility(1, 3)
ELSE
WRITE(stdout, '(5x, 1f8.3, 1f9.4, 1E14.5, 1E14.5, 3E16.6)') etemp * ryd2ev / kelvin2eV, ef0 * ryd2ev, &
dos(itemp), SUM(Fi_check(:)), mobility(1, 1), mobility(1, 2), mobility(1, 3)
dos(itemp), SUM(fi_check(:)), mobility(1, 1), mobility(1, 2), mobility(1, 3)
ENDIF
WRITE(stdout, '(50x, 3E16.6)') mobility(2, 1), mobility(2, 2), mobility(2, 3)
WRITE(stdout, '(50x, 3E16.6)') mobility(3, 1), mobility(3, 2), mobility(3, 3)
@ -676,7 +676,7 @@
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE prtheader()
SUBROUTINE prtheader_mob()
!-----------------------------------------------------------------------
!!
!! This routine print a header for mobility calculation
@ -704,7 +704,7 @@
RETURN
!
!-----------------------------------------------------------------------
END SUBROUTINE prtheader
END SUBROUTINE prtheader_mob
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
@ -852,6 +852,163 @@
END SUBROUTINE print_clock_epw
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE plot_band()
!-----------------------------------------------------------------------
!!
!! This routine writes output files for phonon dispersion and band structure
!! SP : Modified so that it works with the current plotband.x of QE 5
!!
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg
USE phcom, ONLY : nmodes
USE epwcom, ONLY : nbndsub, filqf, filkf
USE elph2, ONLY : etf, nkf, nqtotf, wf, xkf, xqf, nkqtotf, nktotf
USE constants_epw, ONLY : ryd2mev, ryd2ev, zero
USE io_var, ONLY : iufilfreq, iufileig
USE elph2, ONLY : nkqf
USE io_global, ONLY : ionode_id
USE mp, ONLY : mp_barrier, mp_sum
USE mp_global, ONLY : inter_pool_comm, my_pool_id
USE poolgathering, ONLY : poolgather2
!
IMPLICIT NONE
!
! Local variables
INTEGER :: ik
!! Global k-point index
INTEGER :: ikk
!! Index for the k-point
INTEGER :: ikq
!! Index for the q-point
INTEGER :: ibnd
!! Band index
INTEGER :: imode
!! Mode index
INTEGER :: iq
!! Global q-point index
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: dist
!! Distance from G-point
REAL(KIND = DP) :: dprev
!! Previous distance
REAL(KIND = DP) :: dcurr
!! Current distance
REAL(KIND = DP), ALLOCATABLE :: xkf_all(:, :)
!! K-points on the full k grid (all pools)
REAL(KIND = DP), ALLOCATABLE :: etf_all(:, :)
!! Eigenenergies on the full k grid (all pools)
!
IF (filqf /= ' ') THEN
!
IF (my_pool_id == ionode_id) THEN
!
OPEN(iufilfreq, FILE = "phband.freq", FORM = 'formatted')
WRITE(iufilfreq, '(" &plot nbnd=", i4, ", nks=", i6, " /")') nmodes, nqtotf
!
! crystal to cartesian coordinates
CALL cryst_to_cart(nqtotf, xqf, bg, 1)
!
dist = zero
dprev = zero
dcurr = zero
DO iq = 1, nqtotf
!
IF (iq /= 1) THEN
dist = DSQRT((xqf(1, iq) - xqf(1, iq - 1)) * (xqf(1, iq) - xqf(1, iq - 1)) &
+ (xqf(2, iq) - xqf(2, iq - 1)) * (xqf(2, iq) - xqf(2, iq - 1)) &
+ (xqf(3, iq) - xqf(3, iq - 1)) * (xqf(3, iq) - xqf(3, iq - 1)))
ELSE
dist = zero
ENDIF
dcurr = dprev + dist
dprev = dcurr
WRITE(iufilfreq, '(10x, 3f10.6)') xqf(:, iq)
WRITE(iufilfreq, '(1000f14.4)') (wf(imode, iq) * ryd2mev, imode = 1, nmodes)
!
ENDDO
CLOSE(iufilfreq)
!
! back from cartesian to crystal coordinates
CALL cryst_to_cart(nqtotf, xqf, at, -1)
!
ENDIF
ENDIF ! filqf
!
IF (filkf /= ' ') THEN
!
DO ik = 1, nkf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
ENDDO
!
ALLOCATE(xkf_all(3, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('plot_band', 'Error allocating xkf_all', 1)
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('plot_band', 'Error allocating etf_all', 1)
!
#if defined(__MPI)
CALL poolgather2(3, nkqtotf, nkqf, xkf, xkf_all)
CALL poolgather2(nbndsub, nkqtotf, nkqf, etf, etf_all)
CALL mp_barrier(inter_pool_comm)
#else
!
xkf_all = xkf
etf_all = etf
#endif
!
IF (my_pool_id == ionode_id) THEN
!
OPEN(iufileig, FILE = "band.eig", FORM = 'formatted')
WRITE(iufileig, '(" &plot nbnd=", i4, ", nks=", i6, " /")') nbndsub, nktotf
!
! crystal to cartesian coordinates
CALL cryst_to_cart(nkqtotf, xkf_all, bg, 1)
!
dist = zero
dprev = zero
dcurr = zero
DO ik = 1, nktotf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
IF (ikk /= 1) THEN
dist = DSQRT((xkf_all(1, ikk) - xkf_all(1, ikk - 2)) * (xkf_all(1, ikk) - xkf_all(1, ikk - 2)) &
+ (xkf_all(2, ikk) - xkf_all(2, ikk - 2)) * (xkf_all(2, ikk) - xkf_all(2, ikk - 2)) &
+ (xkf_all(3, ikk) - xkf_all(3, ikk - 2)) * (xkf_all(3, ikk) - xkf_all(3, ikk - 2)))
ELSE
dist = 0.d0
ENDIF
dcurr = dprev + dist
dprev = dcurr
WRITE(iufileig, '(10x, 3f10.6)') xkf_all(:, ikk)
WRITE(iufileig, '(1000f20.12)') (etf_all(ibnd, ikk) * ryd2ev, ibnd = 1, nbndsub)
!
ENDDO
CLOSE(iufileig)
!
! back from cartesian to crystal coordinates
CALL cryst_to_cart(nkqtotf, xkf_all, at, -1)
!
ENDIF
CALL mp_barrier(inter_pool_comm)
!
DEALLOCATE(xkf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('plot_band', 'Error deallocating xkf_all', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('plot_band', 'Error deallocating etf_all', 1)
!
ENDIF ! filkf
!
RETURN
!
!----------------------------------------------------------------------------
END SUBROUTINE plot_band
!----------------------------------------------------------------------------
!-------------------------------------------------------------------------
END MODULE printing
!-------------------------------------------------------------------------

197
EPW/src/selfen.f90
View File

@ -43,7 +43,7 @@
USE io_var, ONLY : linewidth_elself
USE phcom, ONLY : nmodes
USE epwcom, ONLY : nbndsub, shortrange, fsthick, eptemp, ngaussw, degaussw, &
eps_acustic, efermi_read, fermi_energy, restart, restart_freq
eps_acustic, efermi_read, fermi_energy, restart, restart_step
USE pwcom, ONLY : ef
USE elph2, ONLY : etf, ibndmin, ibndmax, nkqf, xqf, eta, nbndfst, &
nkf, epf17, wf, wqf, xkf, nkqtotf, adapt_smearing, &
@ -55,7 +55,7 @@
USE mp_global, ONLY : inter_pool_comm
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
USE io_transport, ONLY : electron_write
USE io_selfen, ONLY : selfen_el_write
USE poolgathering, ONLY : poolgather2
!
IMPLICIT NONE
@ -346,14 +346,14 @@
!
! Creation of a restart point
IF (restart) THEN
IF (MOD(iqq, restart_freq) == 0) THEN
IF (MOD(iqq, restart_step) == 0) THEN
WRITE(stdout, '(5x, a, i10)' ) 'Creation of a restart point at ', iqq
CALL mp_sum(sigmar_all, inter_pool_comm)
CALL mp_sum(sigmai_all, inter_pool_comm)
CALL mp_sum(zi_all, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
CALL electron_write(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
CALL selfen_el_write(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
ENDIF
ENDIF
ENDIF ! in case of restart, do not do the first one
@ -941,7 +941,7 @@
USE io_var, ONLY : linewidth_elself
USE epwcom, ONLY : nbndsub, fsthick, eptemp, ngaussw, efermi_read, &
fermi_energy, degaussw, nel, meff, epsiheg, &
restart, restart_freq
restart, restart_step
USE pwcom, ONLY : ef
USE elph2, ONLY : etf, ibndmin, ibndmax, nkqf, xqf, dmef, adapt_smearing, &
nkf, wqf, xkf, nkqtotf, efnew, nbndfst, nktotf, &
@ -952,7 +952,7 @@
USE cell_base, ONLY : omega, alat, bg
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
USE io_transport, ONLY : electron_write
USE io_selfen, ONLY : selfen_el_write
USE poolgathering, ONLY : poolgather2
!
IMPLICIT NONE
@ -1278,14 +1278,14 @@
!
! Creation of a restart point
IF (restart) THEN
IF (MOD(iqq, restart_freq) == 0) THEN
IF (MOD(iqq, restart_step) == 0) THEN
WRITE(stdout, '(5x, a, i10)' ) 'Creation of a restart point at ', iqq
CALL mp_sum(sigmar_all, inter_pool_comm)
CALL mp_sum(sigmai_all, inter_pool_comm)
CALL mp_sum(zi_all, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
CALL electron_write(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
CALL selfen_el_write(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
ENDIF
ENDIF
ENDIF ! in case of restart, do not do the first one
@ -1501,6 +1501,187 @@
END FUNCTION dos_ef_seq
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE nesting_fn_q(iqq, iq)
!-----------------------------------------------------------------------
!!
!! Compute the imaginary part of the phonon self energy due to electron-
!! phonon interaction in the Migdal approximation. This corresponds to
!! the phonon linewidth (half width). The phonon frequency is taken into
!! account in the energy selection rule.
!!
!! Use matrix elements, electronic eigenvalues and phonon frequencies
!! from ep-wannier interpolation.
!!
!-----------------------------------------------------------------------
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE epwcom, ONLY : nbndsub, fsthick, eptemp, ngaussw, degaussw, &
nsmear, delta_smear, efermi_read, fermi_energy
USE pwcom, ONLY : ef
USE elph2, ONLY : ibndmin, etf, wkf, xqf, wqf, nkqf, nktotf, &
nkf, xqf, nbndfst, efnew
USE constants_epw, ONLY : ryd2ev, zero, one, two
USE mp, ONLY : mp_barrier, mp_sum
USE mp_global, ONLY : inter_pool_comm
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: iqq
!! Current q-point index from selecq
INTEGER, INTENT(in) :: iq
!! Current q-point index
!
! Local variables
INTEGER :: ik
!! Counter on the k-point index
INTEGER :: ikk
!! k-point index
INTEGER :: ikq
!! q-point index
INTEGER :: ibnd
!! Counter on bands
INTEGER :: jbnd
!! Counter on bands
INTEGER :: fermicount
!! Number of states on the Fermi surface
INTEGER :: ismear
!! Counter on smearing values
!
REAL(KIND = DP) :: ekk
!! Eigen energy on the fine grid relative to the Fermi level
REAL(KIND = DP) :: ekq
!! Eigen energy of k+q on the fine grid relative to the Fermi level
REAL(KIND = DP) :: ef0
!! Fermi energy level
REAL(KIND = DP) :: weight
!! Imaginary part of the phonhon self-energy factor, sans e-ph matrix elements
REAL(KIND = DP) :: dosef
!! Density of state N(Ef)
REAL(KIND = DP) :: w0g1
!! Dirac delta at k for the imaginary part of $\Sigma$
REAL(KIND = DP) :: w0g2
!! Dirac delta at k+q for the imaginary part of $\Sigma$
REAL(KIND = DP) :: degaussw0
!! degaussw0 = (ismear-1) * delta_smear + degaussw
REAL(KIND = DP) :: inv_degaussw0
!! Inverse degaussw0 for efficiency reasons
REAL(KIND = DP) :: gamma
!! Nesting function
REAL(KIND = DP) :: dos_ef
!! Function returning the density of states at the Fermi level
REAL(KIND = DP) :: w0gauss
!! This function computes the derivative of the Fermi-Dirac function
!! It is therefore an approximation for a delta function
!
IF (iqq == 1) THEN
WRITE(stdout, '(/5x, a)') REPEAT('=', 67)
WRITE(stdout, '(5x, "Nesting Function in the double delta approx")')
WRITE(stdout, '(5x, a/)') REPEAT('=', 67)
!
IF (fsthick < 1.d3) WRITE(stdout, '(/5x, a, f10.6, a)' ) &
'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
WRITE(stdout, '(/5x, a, f10.6, a)' ) 'Golden Rule strictly enforced with T = ', eptemp * ryd2ev, ' eV'
ENDIF
!
! SP: The Gamma function needs to be put to 0 for each q
gamma = zero
!
! Here we loop on smearing values
DO ismear = 1, nsmear
!
degaussw0 = (ismear - 1) * delta_smear + degaussw
inv_degaussw0 = one / degaussw0
!
! Fermi level and corresponding DOS
!
! Note that the weights of k+q points must be set to zero here
! no spin-polarized calculation here
IF (efermi_read) THEN
ef0 = fermi_energy
ELSE
ef0 = efnew
ENDIF
!
dosef = dos_ef(ngaussw, degaussw0, ef0, etf, wkf, nkqf, nbndsub)
! N(Ef) in the equation for lambda is the DOS per spin
dosef = dosef / two
!
IF (iqq == 1) THEN
WRITE(stdout, 100) degaussw0 * ryd2ev, ngaussw
WRITE(stdout, 101) dosef / ryd2ev, ef0 * ryd2ev
ENDIF
!
!
CALL start_clock('nesting')
!
fermicount = 0
DO ik = 1, nkf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! here we must have ef, not ef0, to be consistent with ephwann_shuffle
IF ((MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) .AND. &
(MINVAL(ABS(etf(:, ikq) - ef)) < fsthick)) then
!
fermicount = fermicount + 1
!
DO ibnd = 1, nbndfst
!
ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0
w0g1 = w0gauss(ekk * inv_degaussw0, 0) * inv_degaussw0
!
DO jbnd = 1, nbndfst
!
ekq = etf(ibndmin - 1 + jbnd, ikq) - ef0
w0g2 = w0gauss(ekq *inv_degaussw0, 0) * inv_degaussw0
!
! = k-point weight * [f(E_k) - f(E_k+q)]/ [E_k+q - E_k -w_q +id]
! This is the imaginary part of the phonon self-energy, sans the matrix elements
!
! weight = wkf (ikk) * (wgkk - wgkq) * &
! aimag ( cone / ( ekq - ekk - ci * degaussw ) )
!
! the below expression is positive-definite, but also an approximation
! which neglects some fine features
!
weight = wkf(ikk) * w0g1 * w0g2
!
gamma = gamma + weight
!
ENDDO ! jbnd
ENDDO ! ibnd
ENDIF ! endif fsthick
ENDDO ! loop on k
!
! collect contributions from all pools (sum over k-points)
! this finishes the integral over the BZ (k)
!
CALL mp_sum(gamma, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
!
WRITE(stdout, '(/5x, "iq = ",i5," coord.: ", 3f9.5, " wt: ", f9.5)') iq, xqf(:, iq) , wqf(iq)
WRITE(stdout, '(5x, a)') REPEAT('-', 67)
!
WRITE(stdout, 102) gamma
WRITE(stdout, '(5x,a/)') REPEAT('-', 67)
!
WRITE(stdout, '(/5x, a, i8, a, i8/)') &
'Number of (k,k+q) pairs on the Fermi surface: ', fermicount, ' out of ', nktotf
!
CALL stop_clock('nesting')
ENDDO !smears
!
100 FORMAT(5x, 'Gaussian Broadening: ', f7.3,' eV, ngauss=', i4)
101 FORMAT(5x, 'DOS =', f10.6, ' states/spin/eV/Unit Cell at Ef=', f10.6, ' eV')
102 FORMAT(5x, 'Nesting function (q)=', E15.6, ' [Adimensional]')
!
!-----------------------------------------------------------------------
END SUBROUTINE nesting_fn_q
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
END MODULE selfen
!-----------------------------------------------------------------------

344
EPW/src/spectral_func.f90
View File

@ -40,7 +40,7 @@
USE epwcom, ONLY : nbndsub, eps_acustic, fsthick, eptemp, ngaussw, &
degaussw, wmin_specfun, wmax_specfun, nw_specfun, &
shortrange, efermi_read, fermi_energy, restart, &
restart_freq
restart_step
USE pwcom, ONLY : ef
USE elph2, ONLY : etf, ibndmin, ibndmax, nkqf, xqf, nktotf, efnew, &
epf17, wkf, nkf, wf, wqf, xkf, nkqtotf, adapt_smearing, &
@ -291,7 +291,7 @@
!
! Creation of a restart point
IF (restart) THEN
IF (MOD(iqq, restart_freq) == 0) THEN
IF (MOD(iqq, restart_step) == 0) THEN
WRITE(stdout, '(5x, a, i10)' ) 'Creation of a restart point at ', iqq
CALL mp_sum(esigmar_all, inter_pool_comm)
CALL mp_sum(esigmai_all, inter_pool_comm)
@ -807,7 +807,7 @@
USE io_var, ONLY : iospectral_sup, iospectral
USE epwcom, ONLY : nbndsub, fsthick, eptemp, ngaussw, degaussw, nw_specfun, &
wmin_specfun, wmax_specfun, efermi_read, fermi_energy, &
nel, meff, epsiheg, restart, restart_freq
nel, meff, epsiheg, restart, restart_step
USE pwcom, ONLY : nelec, ef
USE elph2, ONLY : etf, ibndmin, ibndmax, nkqf, nbndfst, wkf, nkf, wqf, xkf, &
nkqtotf, xqf, dmef, esigmar_all, esigmai_all, a_all, &
@ -1139,7 +1139,7 @@
!
! Creation of a restart point
IF (restart) THEN
IF (MOD(iqq, restart_freq) == 0) THEN
IF (MOD(iqq, restart_step) == 0) THEN
WRITE(stdout, '(5x, a, i10)' ) 'Creation of a restart point at ', iqq
CALL mp_sum(esigmar_all, inter_pool_comm)
CALL mp_sum(esigmai_all, inter_pool_comm)
@ -1287,6 +1287,342 @@
END SUBROUTINE spectral_func_pl_q
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE a2f_main()
!-----------------------------------------------------------------------
!!
!! Compute the Eliasberg spectral function
!! in the Migdal approximation.
!!
!! If the q-points are not on a uniform grid (i.e. a line)
!! the function will not be correct
!!
!! 02/2009 works in serial on ionode at the moment. can be parallelized
!! 03/2009 added transport spectral function -- this involves a v_k dot v_kq term
!! in the quantities coming from selfen_phon.f90. Not fully implemented
!! 10/2009 the code is transitioning to 'on-the-fly' phonon selfenergies
!! and this routine is not currently functional
!! 10/2015 RM: added calcution of Tc based on Allen-Dynes formula
!! 09/2019 SP: Cleaning
!!
!
USE kinds, ONLY : DP
USE phcom, ONLY : nmodes
USE cell_base, ONLY : omega
USE epwcom, ONLY : degaussq, delta_qsmear, nqsmear, nqstep, nsmear, eps_acustic, &
delta_smear, degaussw, fsthick, nc
USE elph2, ONLY : nqtotf, wf, wqf, lambda_all, lambda_v_all
USE constants_epw, ONLY : ryd2mev, ryd2ev, kelvin2eV, one, two, zero, kelvin2Ry, pi
USE mp, ONLY : mp_barrier, mp_sum
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
USE io_global, ONLY : stdout
USE io_var, ONLY : iua2ffil, iudosfil, iua2ftrfil, iures
USE io_files, ONLY : prefix
!
IMPLICIT NONE
!
CHARACTER(LEN = 256) :: fila2f
!! File name for Eliashberg spectral function
CHARACTER(LEN = 256) :: fila2ftr
!! File name for transport Eliashberg spectral function
CHARACTER(LEN = 256) :: fildos
!! File name for phonon density of states
CHARACTER(LEN = 256) :: filres
!! File name for resistivity
!
INTEGER :: imode
!! Counter on mode
INTEGER :: iq
!! Counter on the q-point index
INTEGER :: iw
!! Counter on the frequency
INTEGER :: ismear
!! Counter on smearing values (phonons)
INTEGER :: isig
!! Counter on smearing values (electrons)
INTEGER :: i
!! Counter on mu
INTEGER :: itemp
!! Counter on temperature
INTEGER :: ierr
!! Error status
!
REAL(KIND = DP) :: weight
!! Factor in a2f
REAL(KIND = DP) :: temp
!! Temperature
REAL(KIND = DP) :: n
!! Carrier density
REAL(KIND = DP) :: be
!! Bose-Einstein distribution
REAL(KIND = DP) :: prefact
!! Prefactor in resistivity
REAL(KIND = DP) :: lambda_tot
!! Total e-ph coupling strength (summation)
REAL(KIND = DP) :: lambda_tr_tot
!! Total transport e-ph coupling strength (summation)
REAL(KIND = DP) :: degaussq0
!! Phonon smearing
REAL(KIND = DP) :: inv_degaussq0
!! Inverse of the smearing for efficiency reasons
REAL(KIND = DP) :: a2f_tmp
!! Temporary variable for Eliashberg spectral function
REAL(KIND = DP) :: a2f_tr_tmp
!! Temporary variable for transport Eliashberg spectral function
REAL(KIND = DP) :: om_max
!! max phonon frequency increased by 10%
REAL(KIND = DP) :: dw
!! Frequency intervals
REAL(KIND = DP) :: w0
!! Current frequency w(imode, iq)
REAL(KIND = DP) :: l
!! Temporary variable for e-ph coupling strength
REAL(KIND = DP) :: l_tr
!! Temporary variable for transport e-ph coupling strength
REAL(KIND = DP) :: tc
!! Critical temperature
REAL(KIND = DP) :: mu
!! Coulomb pseudopotential
REAL(KIND = DP), EXTERNAL :: w0gauss
!! The derivative of wgauss: an approximation to the delta function
REAL(KIND = DP) :: ww(nqstep)
!! Current frequency
REAL(KIND = DP), ALLOCATABLE :: a2f_(:, :)
!! Eliashberg spectral function for different ismear
REAL(KIND = DP), ALLOCATABLE :: a2f_tr(:, :)
!! Transport Eliashberg spectral function for different ismear
REAL(KIND = DP), ALLOCATABLE :: l_a2f(:)
!! total e-ph coupling strength (a2f_ integration) for different ismear
REAL(KIND = DP), ALLOCATABLE :: l_a2f_tr(:)
!! total transport e-ph coupling strength (a2f_tr integration) for different ismear
REAL(KIND = DP), ALLOCATABLE :: dosph(:, :)
!! Phonon density of states for different for different ismear
REAL(KIND = DP), ALLOCATABLE :: logavg(:)
!! logavg phonon frequency for different ismear
REAL(KIND = DP), ALLOCATABLE :: rho(:, :)
!! Resistivity for different for different ismear
!
CALL start_clock('a2F')
IF (mpime == ionode_id) THEN
!
ALLOCATE(a2f_(nqstep, nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating a2f_', 1)
ALLOCATE(a2f_tr(nqstep, nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating a2f_tr', 1)
ALLOCATE(dosph(nqstep, nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating dosph', 1)
ALLOCATE(l_a2f(nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating l_a2f', 1)
ALLOCATE(l_a2f_tr(nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating l_a2f_tr', 1)
ALLOCATE(logavg(nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating logavg', 1)
! The resitivity is computed for temperature between 0K-1000K by step of 10
! This is hardcoded and needs to be changed here if one wants to modify it
ALLOCATE(rho(100, nqsmear), STAT = ierr)
IF (ierr /= 0) CALL errore('a2f_main', 'Error allocating rho', 1)
!
DO isig = 1, nsmear
!
IF (isig < 10) THEN
WRITE(fila2f, '(a, a6, i1)') TRIM(prefix), '.a2f.0', isig
WRITE(fila2ftr, '(a, a9, i1)') TRIM(prefix), '.a2f_tr.0', isig
WRITE(filres, '(a, a6, i1)') TRIM(prefix), '.res.0', isig
WRITE(fildos, '(a, a8, i1)') TRIM(prefix), '.phdos.0', isig
ELSE
WRITE(fila2f, '(a, a5, i2)') TRIM(prefix), '.a2f.', isig
WRITE(fila2ftr, '(a, a8, i2)') TRIM(prefix), '.a2f_tr.', isig
WRITE(filres, '(a, a5, i2)') TRIM(prefix), '.res.', isig
WRITE(fildos, '(a, a7, i2)') TRIM(prefix), '.phdos.', isig
ENDIF
OPEN(UNIT = iua2ffil, FILE = fila2f, FORM = 'formatted')
OPEN(UNIT = iua2ftrfil, FILE = fila2ftr, FORM = 'formatted')
OPEN(UNIT = iures, FILE = filres, FORM = 'formatted')
OPEN(UNIT = iudosfil, FILE = fildos, FORM = 'formatted')
!
WRITE(stdout, '(/5x, a)') REPEAT('=',67)
WRITE(stdout, '(5x, "Eliashberg Spectral Function in the Migdal Approximation")')
WRITE(stdout, '(5x, a/)') REPEAT('=',67)
!
om_max = 1.1d0 * MAXVAL(wf(:, :)) ! increase by 10%
dw = om_max / DBLE(nqstep)
DO iw = 1, nqstep !
ww(iw) = DBLE(iw) * dw
ENDDO
!
lambda_tot = zero
l_a2f(:) = zero
a2f_(:, :) = zero
lambda_tr_tot = zero
l_a2f_tr(:) = zero
a2f_tr(:, :) = zero
dosph(:, :) = zero
logavg(:) = zero
!
DO ismear = 1, nqsmear
!
degaussq0 = degaussq + (ismear - 1) * delta_qsmear
inv_degaussq0 = one / degaussq0
!
DO iw = 1, nqstep ! loop over points on the a2F(w) graph
!
DO iq = 1, nqtotf ! loop over q-points
DO imode = 1, nmodes ! loop over modes
w0 = wf(imode, iq)
!
IF (w0 > eps_acustic) THEN
!
l = lambda_all(imode, iq, isig)
IF (lambda_all(imode, iq, isig) < 0.d0) l = zero ! sanity check
!
a2f_tmp = wqf(iq) * w0 * l / two
!
weight = w0gauss((ww(iw) - w0) * inv_degaussq0, 0) * inv_degaussq0
a2f_(iw, ismear) = a2f_(iw, ismear) + a2f_tmp * weight
dosph(iw, ismear) = dosph(iw, ismear) + wqf(iq) * weight
!
l_tr = lambda_v_all(imode, iq, isig)
IF (lambda_v_all(imode, iq, isig) < 0.d0) l_tr = zero !sanity check
!
a2f_tr_tmp = wqf(iq) * w0 * l_tr / two
!
a2f_tr(iw, ismear) = a2f_tr(iw, ismear) + a2f_tr_tmp * weight
!
ENDIF
ENDDO
ENDDO
!
! output a2f
!
IF (ismear == nqsmear) WRITE(iua2ffil, '(f12.7, 15f12.7)') ww(iw) * ryd2mev, a2f_(iw, :)
IF (ismear == nqsmear) WRITE(iua2ftrfil, '(f12.7, 15f12.7)') ww(iw) * ryd2mev, a2f_tr(iw, :)
IF (ismear == nqsmear) WRITE(iudosfil, '(f12.7, 15f12.7)') ww(iw) * ryd2mev, dosph(iw, :) / ryd2mev
!
! do the integral 2 int (a2F(w)/w dw)
!
l_a2f(ismear) = l_a2f(ismear) + two * a2f_(iw, ismear) / ww(iw) * dw
l_a2f_tr(ismear) = l_a2f_tr(ismear) + two * a2f_tr(iw, ismear) / ww(iw) * dw
logavg(ismear) = logavg(ismear) + two * a2f_(iw, ismear) * LOG(ww(iw)) / ww(iw) * dw
!
ENDDO
!
logavg(ismear) = EXP(logavg(ismear) / l_a2f(ismear))
!
ENDDO
!
DO iq = 1, nqtotf ! loop over q-points
DO imode = 1, nmodes ! loop over modes
IF (lambda_all(imode, iq, isig) > 0.d0 .AND. wf(imode, iq) > eps_acustic ) &
lambda_tot = lambda_tot + wqf(iq) * lambda_all(imode, iq, isig)
IF (lambda_v_all(imode, iq, isig) > 0.d0 .AND. wf(imode, iq) > eps_acustic) &
lambda_tr_tot = lambda_tr_tot + wqf(iq) * lambda_v_all(imode, iq, isig)
ENDDO
ENDDO
WRITE(stdout, '(5x, a, f12.7)') "lambda : ", lambda_tot
WRITE(stdout, '(5x, a, f12.7)') "lambda_tr : ", lambda_tr_tot
WRITE(stdout, '(a)') " "
!
!
! Allen-Dynes estimate of Tc for ismear = 1
!
WRITE(stdout, '(5x, a, f12.7, a)') "Estimated Allen-Dynes Tc"
WRITE(stdout, '(a)') " "
WRITE(stdout, '(5x, a, f12.7, a, f12.7)') "logavg = ", logavg(1), " l_a2f = ", l_a2f(1)
DO i = 1, 6
!
mu = 0.1d0 + 0.02d0 * DBLE(i - 1)
tc = logavg(1) / 1.2d0 * EXP(-1.04d0 * (1.d0 + l_a2f(1)) / (l_a2f(1) - mu * ( 1.d0 + 0.62d0 * l_a2f(1))))
! tc in K
!
tc = tc * ryd2ev / kelvin2eV
!SP: IF Tc is too big, it is not physical
IF (tc < 1000.0) THEN
WRITE(stdout, '(5x, a, f6.2, a, f22.12, a)') "mu = ", mu, " Tc = ", tc, " K"
ENDIF
!
ENDDO
!
rho(:, :) = zero
! Now compute the Resistivity of Metal using the Ziman formula
! rho(T,smearing) = 4 * pi * me/(n * e**2 * kb * T) int dw hbar w a2F_tr(w,smearing) n(w,T)(1+n(w,T))
! n is the number of electron per unit volume and n(w,T) is the Bose-Einstein distribution
! Usually this means "the number of electrons that contribute to the mobility" and so it is typically 8 (full shell)
! but not always. You might want to check this.
!
n = nc / omega
WRITE(iures, '(a)') '# Temperature [K] &
Resistivity [micro Ohm cm] for different Phonon smearing (meV) '
WRITE(iures, '("# ", 15f12.7)') ((degaussq + (ismear - 1) * delta_qsmear) * ryd2mev, ismear = 1, nqsmear)
DO ismear = 1, nqsmear
DO itemp = 1, 100 ! Per step of 10K
temp = itemp * 10.d0 * kelvin2Ry
! omega is the volume of the primitive cell in a.u.
!
prefact = 4.d0 * pi / (temp * n)
DO iw = 1, nqstep ! loop over points on the a2F(w)
!
be = one / (EXP(ww(iw) / temp) - one)
! Perform the integral with rectangle.
rho(itemp, ismear) = rho(itemp, ismear) + prefact * ww(iw) * a2f_tr(iw, ismear) * be * (1.d0 + be) * dw
!
ENDDO
! From a.u. to micro Ohm cm
! Conductivity 1 a.u. = 2.2999241E6 S/m
! Now to go from Ohm*m to micro Ohm cm we need to multiply by 1E8
rho(itemp, ismear) = rho(itemp, ismear) * 1E8 / 2.2999241E6
IF (ismear == nqsmear) WRITE (iures, '(i8, 15f12.7)') itemp * 10, rho(itemp, :)
ENDDO
ENDDO
CLOSE(iures)
!
WRITE(iua2ffil, *) "Integrated el-ph coupling"
WRITE(iua2ffil, '(" # ", 15f12.7)') l_a2f(:)
WRITE(iua2ffil, *) "Phonon smearing (meV)"
WRITE(iua2ffil, '(" # ", 15f12.7)') ((degaussq + (ismear - 1) * delta_qsmear) * ryd2mev, ismear = 1, nqsmear)
WRITE(iua2ffil, '(" Electron smearing (eV)", f12.7)') ((isig - 1) * delta_smear + degaussw) * ryd2ev
WRITE(iua2ffil, '(" Fermi window (eV)", f12.7)') fsthick * ryd2ev
WRITE(iua2ffil, '(" Summed el-ph coupling ", f12.7)') lambda_tot
CLOSE(iua2ffil)
!
WRITE(iua2ftrfil, *) "Integrated el-ph coupling"
WRITE(iua2ftrfil, '(" # ", 15f12.7)') l_a2f_tr(:)
WRITE(iua2ftrfil, *) "Phonon smearing (meV)"
WRITE(iua2ftrfil, '(" # ", 15f12.7)') ((degaussq + (ismear - 1) * delta_qsmear) * ryd2mev, ismear = 1, nqsmear)
WRITE(iua2ftrfil, '(" Electron smearing (eV)", f12.7)') ((isig - 1) * delta_smear + degaussw) * ryd2ev
WRITE(iua2ftrfil, '(" Fermi window (eV)", f12.7)') fsthick * ryd2ev
WRITE(iua2ftrfil, '(" Summed el-ph coupling ", f12.7)') lambda_tot
CLOSE(iua2ftrfil)
!
CLOSE(iudosfil)
!
ENDDO ! isig
!
DEALLOCATE(l_a2f, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating l_a2f', 1)
DEALLOCATE(l_a2f_tr, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating l_a2f_tr', 1)
DEALLOCATE(a2f_, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating a2f', 1)
DEALLOCATE(a2f_tr, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating a2f_tr', 1)
DEALLOCATE(rho, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating rho', 1)
DEALLOCATE(dosph, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating dosph', 1)
DEALLOCATE(logavg, STAT = ierr)
IF (ierr /= 0) CALL errore('eliashberg_a2f', 'Error deallocating logavg', 1)
!
ENDIF
!
CALL stop_clock('a2F')
CALL print_clock('a2F')
!
RETURN
!
!-----------------------------------------------------------------------
END SUBROUTINE a2f_main
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
END MODULE spectral_func
!-----------------------------------------------------------------------

2
EPW/src/supercond.f90
View File

@ -287,7 +287,7 @@
!! anisotropic e-ph coupling strength
!
! This is only a quick fix since the routine was written for parallel execution - FG June 2014
#if ! defined(__MPI)
#if defined(__MPI)
npool = 1
my_pool_id = 0
#endif

144
EPW/src/transport.f90
View File

@ -28,7 +28,7 @@
USE phcom, ONLY : nmodes
USE epwcom, ONLY : nbndsub, fsthick, eps_acustic, degaussw, restart, &
nstemp, scattering_serta, scattering_0rta, shortrange, &
restart_freq, restart_filq, vme, assume_metal
restart_step, restart_filq, vme, assume_metal
USE pwcom, ONLY : ef
USE elph2, ONLY : ibndmin, etf, nkqf, nkf, dmef, vmef, wf, wqf, &
epf17, nkqtotf, inv_tau_all, inv_tau_allcb, &
@ -373,7 +373,7 @@
!
! Creation of a restart point
IF (restart) THEN
IF (MOD(iqq, restart_freq) == 0) THEN
IF (MOD(iqq, restart_step) == 0) THEN
WRITE(stdout, '(a)' ) ' Creation of a restart point'
!
! The mp_sum will aggreage the results on each k-points.
@ -623,11 +623,11 @@
!! Number of points in the BZ corresponding to a point in IBZ
INTEGER :: ierr
!! Error status
INTEGER :: BZtoIBZ_tmp(nkf1 * nkf2 * nkf3)
INTEGER :: bztoibz_tmp(nkf1 * nkf2 * nkf3)
!! Temporary mapping
INTEGER :: BZtoIBZ(nkf1 * nkf2 * nkf3)
INTEGER :: bztoibz(nkf1 * nkf2 * nkf3)
!! BZ to IBZ mapping
INTEGER(SIK2) :: s_BZtoIBZ(nkf1 * nkf2 * nkf3)
INTEGER(SIK2) :: s_bztoibz(nkf1 * nkf2 * nkf3)
!! symmetry
REAL(KIND = DP) :: ekk
!! Energy relative to Fermi level: $$\varepsilon_{n\mathbf{k}}-\varepsilon_F$$
@ -657,11 +657,11 @@
!! Mobility along the zz axis after diagonalization [cm^2/Vs]
REAL(KIND = DP) :: vkk(3, nbndfst)
!! Electron velocity vector for a band.
REAL(KIND = DP) :: Sigma(9, nstemp)
REAL(KIND = DP) :: sigma(9, nstemp)
!! Conductivity matrix in vector form
REAL(KIND = DP) :: SigmaZ(9, nstemp)
REAL(KIND = DP) :: sigmaZ(9, nstemp)
!! Conductivity matrix in vector form with Znk
REAL(KIND = DP) :: Sigma_m(3, 3, nstemp)
REAL(KIND = DP) :: sigma_m(3, 3, nstemp)
!! Conductivity matrix
REAL(KIND = DP) :: sigma_up(3, 3)
!! Conductivity matrix in upper-triangle
@ -669,7 +669,7 @@
!! Eigenvalues from the diagonalized conductivity matrix
REAL(KIND = DP) :: sigma_vect(3, 3)
!! Eigenvectors from the diagonalized conductivity matrix
REAL(KIND = DP) :: Znk
REAL(KIND = DP) :: znk
!! Real Znk from \lambda_nk (called zi_allvb or zi_allcb)
REAL(KIND = DP) :: tdf_sigma(9)
!! Temporary file
@ -926,21 +926,21 @@
!
! SP - Uncomment to use symmetries on velocities
IF (mp_mesh_k) THEN
BZtoIBZ(:) = 0
s_BZtoIBZ(:) = 0
bztoibz(:) = 0
s_bztoibz(:) = 0
!
CALL set_sym_bl()
! What we get from this call is BZtoIBZ
CALL kpoint_grid_epw(nrot, time_reversal, .FALSE., s, t_rev, nkf1, nkf2, nkf3, BZtoIBZ, s_BZtoIBZ)
! What we get from this call is bztoibz
CALL kpoint_grid_epw(nrot, time_reversal, .FALSE., s, t_rev, nkf1, nkf2, nkf3, bztoibz, s_bztoibz)
!
IF (iterative_bte) THEN
! Now we have to remap the points because the IBZ k-points have been
! changed to to load balancing.
BZtoIBZ_tmp(:) = 0
bztoibz_tmp(:) = 0
DO ikbz = 1, nkf1 * nkf2 * nkf3
BZtoIBZ_tmp(ikbz) = map_rebal(BZtoIBZ(ikbz))
bztoibz_tmp(ikbz) = map_rebal(bztoibz(ikbz))
ENDDO
BZtoIBZ(:) = BZtoIBZ_tmp(:)
bztoibz(:) = bztoibz_tmp(:)
ENDIF
!
ENDIF
@ -967,8 +967,8 @@
! 7 = (3,1) = zx, 8 = (3,2) = zy, 9 = (3,3) = zz
! this can be reduced to 6 if we take into account symmetry xy=yx, ...
tdf_sigma(:) = zero
Sigma(:, :) = zero
SigmaZ(:, :) = zero
sigma(:, :) = zero
sigmaZ(:, :) = zero
!
ENDIF
!
@ -1006,10 +1006,10 @@
DO ikbz = 1, nkf1 * nkf2 * nkf3
! If the k-point from the full BZ is related by a symmetry operation
! to the current k-point, then take it.
IF (BZtoIBZ(ikbz) == ik + lower_bnd - 1) THEN
IF (bztoibz(ikbz) == ik + lower_bnd - 1) THEN
nb = nb + 1
! Transform the symmetry matrix from Crystal to cartesian
sa(:, :) = DBLE(s(:, :, s_BZtoIBZ(ikbz)))
sa(:, :) = DBLE(s(:, :, s_bztoibz(ikbz)))
sb = MATMUL(bg, sa)
sr(:, :) = MATMUL(at, TRANSPOSE(sb))
CALL DGEMV('n', 3, 3, 1.d0, sr, 3, vk_cart(:), 1, 0.d0, v_rot(:), 1)
@ -1055,13 +1055,13 @@
! (-df_nk/dE_nk) = (f_nk)*(1-f_nk)/ (k_B T)
dfnk = w0gauss(ekk / etemp, -99) / etemp
! electrical conductivity
Sigma(:, itemp) = Sigma(:, itemp) + dfnk * tdf_sigma(:) * tau
sigma(:, itemp) = Sigma(:, itemp) + dfnk * tdf_sigma(:) * tau
!
! Now do the same but with Znk multiplied
! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} )
Znk = one / (one + zi_allvb(itemp, ibnd, ik + lower_bnd - 1))
znk = one / (one + zi_allvb(itemp, ibnd, ik + lower_bnd - 1))
tau = one / (Znk * inv_tau_all(itemp, ibnd, ik + lower_bnd - 1))
SigmaZ(:, itemp) = SigmaZ(:, itemp) + dfnk * tdf_sigma(:) * tau
sigmaz(:, itemp) = sigmaz(:, itemp) + dfnk * tdf_sigma(:) * tau
!
ENDIF
ENDDO ! ibnd
@ -1070,8 +1070,8 @@
!
! The k points are distributed among pools: here we collect them
!
CALL mp_sum(Sigma(:, itemp), world_comm)
CALL mp_sum(SigmaZ(:, itemp), world_comm)
CALL mp_sum(sigma(:, itemp), world_comm)
CALL mp_sum(sigmaz(:, itemp), world_comm)
!
ENDDO ! nstemp
!
@ -1092,10 +1092,10 @@
!
DO itemp = 1, nstemp
etemp = transp_temp(itemp)
! Sigma in units of 1/(a.u.) is converted to 1/(Ohm * m)
! sigma in units of 1/(a.u.) is converted to 1/(Ohm * m)
IF (mpime == meta_ionode_id) THEN
WRITE(iufilsigma, '(11E16.8)') ef0(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, &
conv_factor1 * Sigma(:, itemp) * inv_cell
conv_factor1 * sigma(:, itemp) * inv_cell
ENDIF
carrier_density = 0.0
!
@ -1120,15 +1120,15 @@
! 4 = (2,1) = yx, 5 = (2,2) = yy, 6 = (2,3) = yz
! 7 = (3,1) = zx, 8 = (3,2) = zy, 9 = (3,3) = zz
sigma_up(:, :) = zero
sigma_up(1, 1) = Sigma(1, itemp)
sigma_up(1, 2) = Sigma(2, itemp)
sigma_up(1, 3) = Sigma(3, itemp)
sigma_up(2, 1) = Sigma(4, itemp)
sigma_up(2, 2) = Sigma(5, itemp)
sigma_up(2, 3) = Sigma(6, itemp)
sigma_up(3, 1) = Sigma(7, itemp)
sigma_up(3, 2) = Sigma(8, itemp)
sigma_up(3, 3) = Sigma(9, itemp)
sigma_up(1, 1) = sigma(1, itemp)
sigma_up(1, 2) = sigma(2, itemp)
sigma_up(1, 3) = sigma(3, itemp)
sigma_up(2, 1) = sigma(4, itemp)
sigma_up(2, 2) = sigma(5, itemp)
sigma_up(2, 3) = sigma(6, itemp)
sigma_up(3, 1) = sigma(7, itemp)
sigma_up(3, 2) = sigma(8, itemp)
sigma_up(3, 3) = sigma(9, itemp)
!
CALL rdiagh(3, sigma_up, 3, sigma_eig, sigma_vect)
!
@ -1184,8 +1184,8 @@
etemp = transp_temp(itemp)
IF (itemp == 1) THEN
tdf_sigma(:) = zero
Sigma(:, :) = zero
SigmaZ(:, :) = zero
sigma(:, :) = zero
sigmaZ(:, :) = zero
ENDIF
DO ik = 1, nkf
ikk = 2 * ik - 1
@ -1214,10 +1214,10 @@
DO ikbz = 1, nkf1 * nkf2 * nkf3
! If the k-point from the full BZ is related by a symmetry operation
! to the current k-point, then take it.
IF (BZtoIBZ(ikbz) == ik + lower_bnd - 1) THEN
IF (bztoibz(ikbz) == ik + lower_bnd - 1) THEN
nb = nb + 1
! Transform the symmetry matrix from Crystal to cartesian
sa(:, :) = DBLE(s(:, :, s_BZtoIBZ(ikbz)))
sa(:, :) = DBLE(s(:, :, s_bztoibz(ikbz)))
sb = MATMUL(bg, sa)
sr(:, :) = MATMUL(at, TRANSPOSE(sb))
CALL DGEMV('n', 3, 3, 1.d0, sr, 3, vk_cart(:), 1, 0.d0, v_rot(:), 1)
@ -1259,13 +1259,13 @@
ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0(itemp)
tau = one / inv_tau_all(itemp, ibnd, ik + lower_bnd - 1)
dfnk = w0gauss(ekk / etemp, -99) / etemp
Sigma(:, itemp) = Sigma(:, itemp) + dfnk * tdf_sigma(:) * tau
sigma(:, itemp) = Sigma(:, itemp) + dfnk * tdf_sigma(:) * tau
!
! Now do the same but with Znk multiplied
! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} )
Znk = one / (one + zi_allvb(itemp, ibnd, ik + lower_bnd - 1))
znk = one / (one + zi_allvb(itemp, ibnd, ik + lower_bnd - 1))
tau = one / (Znk * inv_tau_all(itemp, ibnd, ik + lower_bnd - 1))
SigmaZ(:, itemp) = SigmaZ(:, itemp) + dfnk * tdf_sigma(:) * tau
sigmaz(:, itemp) = sigmaz(:, itemp) + dfnk * tdf_sigma(:) * tau
ENDIF
ENDDO
ELSE ! In this case we have 2 Fermi levels
@ -1290,10 +1290,10 @@
DO ikbz = 1, nkf1 * nkf2 * nkf3
! If the k-point from the full BZ is related by a symmetry operation
! to the current k-point, then take it.
IF (BZtoIBZ(ikbz) == ik + lower_bnd - 1) THEN
IF (bztoibz(ikbz) == ik + lower_bnd - 1) THEN
nb = nb + 1
! Transform the symmetry matrix from Crystal to cartesian
sa(:, :) = DBLE(s(:, :, s_BZtoIBZ(ikbz)))
sa(:, :) = DBLE(s(:, :, s_bztoibz(ikbz)))
sb = MATMUL(bg, sa)
sr(:, :) = MATMUL(at, TRANSPOSE(sb))
CALL DGEMV('n', 3, 3, 1.d0, sr, 3, vk_cart(:), 1, 0.d0, v_rot(:), 1)
@ -1335,11 +1335,11 @@
ekk = etf(ibndmin - 1 + ibnd, ikk) - efcb(itemp)
tau = one / inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1)
dfnk = w0gauss(ekk / etemp, -99) / etemp
Sigma(:, itemp) = Sigma(:, itemp) + dfnk * tdf_sigma(:) * tau
sigma(:, itemp) = sigma(:, itemp) + dfnk * tdf_sigma(:) * tau
!
! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} )
Znk = one / (one + zi_allcb(itemp, ibnd, ik + lower_bnd - 1))
tau = one / (Znk * inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1))
znk = one / (one + zi_allcb(itemp, ibnd, ik + lower_bnd - 1))
tau = one / (znk * inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1))
ij = 0
DO j = 1, 3
DO i = 1, 3
@ -1347,14 +1347,14 @@
tdf_sigma(ij) = vkk(i, ibnd) * vkk(j, ibnd) * tau
ENDDO
ENDDO
SigmaZ(:, itemp) = SigmaZ(:, itemp) + wkf(ikk) * dfnk * tdf_sigma(:)
sigmaZ(:, itemp) = sigmaZ(:, itemp) + wkf(ikk) * dfnk * tdf_sigma(:)
ENDIF
ENDDO ! ibnd
ENDIF ! etcb
ENDIF ! endif fsthick
ENDDO ! end loop on k
CALL mp_sum(Sigma(:, itemp), world_comm)
CALL mp_sum(SigmaZ(:, itemp), world_comm)
CALL mp_sum(sigma(:, itemp), world_comm)
CALL mp_sum(sigmaZ(:, itemp), world_comm)
!
ENDDO ! nstemp
IF (mpime == meta_ionode_id) THEN
@ -1373,13 +1373,13 @@
DO itemp = 1, nstemp
etemp = transp_temp(itemp)
IF (mpime == meta_ionode_id) THEN
! Sigma in units of 1/(a.u.) is converted to 1/(Ohm * m)
! sigma in units of 1/(a.u.) is converted to 1/(Ohm * m)
IF (ABS(efcb(itemp)) < eps4) THEN
WRITE(iufilsigma, '(11E16.8)') ef0(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, &
conv_factor1 * Sigma(:, itemp) * inv_cell
conv_factor1 * sigma(:, itemp) * inv_cell
ELSE
WRITE(iufilsigma, '(11E16.8)') efcb(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, &
conv_factor1 * Sigma(:, itemp) * inv_cell
conv_factor1 * sigma(:, itemp) * inv_cell
ENDIF
ENDIF
carrier_density = 0.0
@ -1411,15 +1411,15 @@
! 4 = (2,1) = yx, 5 = (2,2) = yy, 6 = (2,3) = yz
! 7 = (3,1) = zx, 8 = (3,2) = zy, 9 = (3,3) = zz
sigma_up(:, :) = zero
sigma_up(1, 1) = Sigma(1, itemp)
sigma_up(1, 2) = Sigma(2, itemp)
sigma_up(1, 3) = Sigma(3, itemp)
sigma_up(2, 1) = Sigma(4, itemp)
sigma_up(2, 2) = Sigma(5, itemp)
sigma_up(2, 3) = Sigma(6, itemp)
sigma_up(3, 1) = Sigma(7, itemp)
sigma_up(3, 2) = Sigma(8, itemp)
sigma_up(3, 3) = Sigma(9, itemp)
sigma_up(1, 1) = sigma(1, itemp)
sigma_up(1, 2) = sigma(2, itemp)
sigma_up(1, 3) = sigma(3, itemp)
sigma_up(2, 1) = sigma(4, itemp)
sigma_up(2, 2) = sigma(5, itemp)
sigma_up(2, 3) = sigma(6, itemp)
sigma_up(3, 1) = sigma(7, itemp)
sigma_up(3, 2) = sigma(8, itemp)
sigma_up(3, 3) = sigma(9, itemp)
!
CALL rdiagh(3, sigma_up, 3, sigma_eig, sigma_vect)
!
@ -1446,15 +1446,15 @@
! ENDIF
! Now do the mobility with Znk factor ----------------------------------------------------------
sigma_up(:, :) = zero
sigma_up(1, 1) = SigmaZ(1, itemp)
sigma_up(1, 2) = SigmaZ(2, itemp)
sigma_up(1, 3) = SigmaZ(3, itemp)
sigma_up(2, 1) = SigmaZ(4, itemp)
sigma_up(2, 2) = SigmaZ(5, itemp)
sigma_up(2, 3) = SigmaZ(6, itemp)
sigma_up(3, 1) = SigmaZ(7, itemp)
sigma_up(3, 2) = SigmaZ(8, itemp)
sigma_up(3, 3) = SigmaZ(9, itemp)
sigma_up(1, 1) = sigmaZ(1, itemp)
sigma_up(1, 2) = sigmaZ(2, itemp)
sigma_up(1, 3) = sigmaZ(3, itemp)
sigma_up(2, 1) = sigmaZ(4, itemp)
sigma_up(2, 2) = sigmaZ(5, itemp)
sigma_up(2, 3) = sigmaZ(6, itemp)
sigma_up(3, 1) = sigmaZ(7, itemp)
sigma_up(3, 2) = sigmaZ(8, itemp)
sigma_up(3, 3) = sigmaZ(9, itemp)
CALL rdiagh(3, sigma_up, 3, sigma_eig, sigma_vect)
mobility_xx = (sigma_eig(1) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ)
mobility_yy = (sigma_eig(2) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ)