Commit be93a2b5 authored by sponce's avatar sponce

All unused variables have been removed.

CMPLX (A, B) is replace by by CMPLX (A, B, kind=DP) as it was leading to loss of precision.



git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@12676 c92efa57-630b-4861-b058-cf58834340f0
parent 2d2f2b64
......@@ -29,7 +29,6 @@
USE pwcom, ONLY : nkstot
#endif
!
USE ions_base, ONLY : nat
USE pwcom, ONLY : nbnd, ngm, doublegrid, nks
USE kinds, ONLY : DP
USE modes, ONLY : nmodes, nirr, npert, u
......
......@@ -598,8 +598,8 @@
! bring epmatq in the mode representation of iq_first,
! and then in the cartesian representation of iq
!
CALL rotate_eigenm ( iq_first, iq, nqc, isym, nsym, s, invs, irt, &
rtau, xq, isq, cz1, cz2 )
CALL rotate_eigenm ( iq_first, nqc, isym, s, invs, irt, &
rtau, xq, cz1, cz2 )
!
CALL rotate_epmat ( cz1, cz2, xq, nqc, lwin, lwinq )
!DBSP
......@@ -641,8 +641,8 @@
! bring epmatq in the mode representation of iq_first,
! and then in the cartesian representation of iq
!
CALL rotate_eigenm ( iq_first, iq, nqc, isym, nsym, s, invs, irt, &
rtau, xq, isq, cz1, cz2 )
CALL rotate_eigenm ( iq_first, nqc, isym, s, invs, irt, &
rtau, xq, cz1, cz2 )
!
CALL rotate_epmat ( cz1, cz2, xq, nqc, lwin, lwinq )
!
......
......@@ -188,7 +188,7 @@
! Transform of position matrix elements
! PRB 74 195118 (2006)
IF (vme) CALL vmebloch2wan &
( nbnd, nbndsub, nks, nks, nkstot, xk, cu, nrr_k, irvec, wslen )
( nbnd, nbndsub, nks, nkstot, xk, cu, nrr_k, irvec, wslen )
!
! Electron-Phonon vertex (Bloch el and Bloch ph -> Wannier el and Bloch ph)
!
......@@ -555,7 +555,7 @@
!
IF (lpolar) THEN
!
CALL compute_bmn_para2( nbndsub, nkstot, cufkk, cufkq, bmatf )
CALL compute_bmn_para2( nbndsub, cufkk, cufkq, bmatf )
!
IF ( (abs(xxq(1)).gt.eps) .or. (abs(xxq(2)).gt.eps) .or. (abs(xxq(3)).gt.eps) ) THEN
CALL cryst_to_cart (1, xxq, bg, 1)
......@@ -763,7 +763,7 @@
!
IF (lpolar) THEN
!
CALL compute_bmn_para2( nbndsub, nkstot, cufkk, cufkq, bmatf )
CALL compute_bmn_para2( nbndsub, cufkk, cufkq, bmatf )
!
IF ( (abs(xxq(1)).gt.eps) .or. (abs(xxq(2)).gt.eps) .or. (abs(xxq(3)).gt.eps) ) THEN
CALL cryst_to_cart (1, xxq, bg, 1)
......
......@@ -28,7 +28,7 @@
PUBLIC :: epwdata, iundmedata, iunvmedata, iunksdata, iudyn, iukgmap, iuepb,&
iunepmatf, iurecover, iufilfreq, iufilegnv, iufileph, iufilkqmap, &
iufilikmap, iueig, iunepmatwp, iunepmatwe, iunkf, iunqf, &
iufileig
iufileig, iukmap
PUBLIC :: iuwinfil, iun_plot, iuukk, iuprojfil !, iummn
!
! Output of physically relevant quantities (60-100)
......@@ -92,6 +92,7 @@
INTEGER :: iunkf = 117 ! The unit with the fine k-point mesh in crystal coord.
INTEGER :: iunqf = 118 ! The unit with the fine q-point mesh in crystal coord.
INTEGER :: iufileig = 119 ! The unit with eigenenergies [band.eig]
INTEGER :: iukmap = 120 ! Unit for the k-point map generation
!
! Output quantites related to Wannier (201-250)
......
......@@ -495,7 +495,6 @@ END SUBROUTINE setup_nnkp
SUBROUTINE scan_file_to (iun_nnkp,keyword,found)
!-----------------------------------------------------------------------
!
USE io_global, ONLY : stdout
IMPLICIT NONE
CHARACTER(len=*), intent(in) :: keyword
INTEGER, intent(in) :: iun_nnkp
......@@ -620,17 +619,16 @@ END SUBROUTINE run_wannier
!-----------------------------------------------------------------------
SUBROUTINE compute_amn_para
!-----------------------------------------------------------------------
! adapted from compute_amn in pw2wannier90.f90
! parallelization on k-points has been added
! 10/2008 Jesse Noffsinger UC Berkeley
!
!! adapted from compute_amn in pw2wannier90.f90
!! parallelization on k-points has been added
!! 10/2008 Jesse Noffsinger UC Berkeley
!!
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
USE klist, ONLY : xk, nks, igk_k, ngk
USE klist, ONLY : xk, nks, igk_k
USE wvfct, ONLY : nbnd, npw, npwx, g2kin
USE gvecw, ONLY : ecutwfc
USE wavefunctions_module, ONLY : evc
USE units_ph, ONLY : lrwfc, iuwfc
USE gvect, ONLY : g, ngm
USE cell_base, ONLY : tpiba2
USE uspp, ONLY : nkb, vkb
......@@ -871,7 +869,6 @@ SUBROUTINE compute_mmn_para
logical :: any_uspp, exst
integer :: nkq, nkq_abs, ipool, ipol, istart, iend
integer :: ik_g, ikp_g, ind0 !, iummn
INTEGER :: ibnd_n, ibnd_m
!
any_uspp = ANY( upf(:)%tvanp )
!
......@@ -1043,7 +1040,7 @@ SUBROUTINE compute_mmn_para
DO na = 1, nat
!
arg = DOT_PRODUCT( dxk(:,ind), tau(:,na) ) * twopi
phase1 = CMPLX ( COS(arg), -SIN(arg) )
phase1 = CMPLX ( COS(arg), -SIN(arg), kind=DP )
!
IF ( ityp(na) == nt ) THEN
DO jh = 1, nh(nt)
......@@ -1226,7 +1223,6 @@ SUBROUTINE compute_pmn_para
!
USE io_global, ONLY : stdout
#ifdef __PARA
USE mp_global, ONLY : intra_pool_comm
USE mp, ONLY : mp_sum
#endif
USE kinds, ONLY : DP
......
......@@ -9,22 +9,20 @@
!--------------------------------------------------------------
subroutine readgmap ( nkstot, ngxx, ng0vec, g0vec_all_r, lower_bnd)
!--------------------------------------------------------------
!
! read map of G vectors G -> G-G_0 for a given q point
! (this is used for the folding of k+q into the first BZ)
!
!
!!
!! read map of G vectors G -> G-G_0 for a given q point
!! (this is used for the folding of k+q into the first BZ)
!!
!!
!--------------------------------------------------------------
#ifdef __PARA
USE mp_global,ONLY : my_pool_id, inter_pool_comm
USE mp_global,ONLY : inter_pool_comm
USE mp, ONLY : mp_bcast,mp_max,mp_min
USE mp_world, ONLY : mpime
use io_global,ONLY : ionode_id
#endif
USE kinds, ONLY : DP
use io_epw, ONLY : iukgmap
use wvfct, ONLY : npwx
use io_epw, ONLY : iukgmap, iukmap
use pwcom, ONLY : nks
use elph2, ONLY : shift, gmap, igk_k_all, ngk_all
USE io_files, ONLY : prefix
......@@ -32,13 +30,17 @@
!
! variables for folding of k+q grid
!
INTEGER :: ng0vec, ngxx, iukmap, nkstot, lower_bnd
! kpoint blocks (k points of all pools)
! number of inequivalent such translations
! bound for the allocation of the array gmap
! unit with map k --> k+q+G
REAL(kind=DP) :: g0vec_all_r(3,125)
! G-vectors needed to fold the k+q grid into the k grid, cartesian coord.
INTEGER, INTENT (in) :: nkstot
!! Total number of k and q points
INTEGER, INTENT (inout) :: ngxx
!! kpoint blocks (k points of all pools)
INTEGER, INTENT (out) :: ng0vec
!!
INTEGER, INTENT (in) :: lower_bnd
!! Lower bound for the k-parallellization
!
REAL(kind=DP), INTENT (out) :: g0vec_all_r(3,125)
!! G-vectors needed to fold the k+q grid into the k grid, cartesian coord.
!
! work variables
!
......@@ -69,7 +71,7 @@
ngxx = 0
DO ik = 1, nks
!
IF (maxval(igk_k_all(1:ngk_all(ik+lower_bnd-1),ik+lower_bnd-1)).gt.ngxx)&
IF (maxval(igk_k_all(1:ngk_all(ik+lower_bnd-1),ik+lower_bnd-1)) > ngxx)&
& ngxx = maxval(igk_k_all(1:ngk_all(ik+lower_bnd-1),ik+lower_bnd-1))
!
ENDDO
......@@ -104,7 +106,6 @@
! 'fake' readin is because the gmap appears *after* the
! wrong kmap.
!
iukmap = 97
open ( unit = iukmap, file = trim(prefix)//'.kmap', form = 'formatted',status='old',iostat=ios)
IF (ios /= 0) call errore ('readgmap', 'error opening kmap file', iukmap)
DO ik = 1, nkstot
......
......@@ -10,10 +10,10 @@
subroutine readmat_shuffle2 ( iq_irr, nqc_irr, nq, iq_first, sxq, imq, isq,&
invs, s, irt, rtau)
!-----------------------------------------------------------------------
!
! read dynamical matrix for the q points
! iq_first, iq_first+1, ... iq_first+nq-1
!
!!
!! read dynamical matrix for the q points
!! iq_first, iq_first+1, ... iq_first+nq-1
!!
!-----------------------------------------------------------------------
USE kinds, ONLY : DP
use io_files, ONLY : prefix
......@@ -31,10 +31,6 @@
USE io_global, ONLY : ionode, stdout
USE constants_epw, ONLY : cone, czero, twopi, rydcm1
USE io_epw, ONLY : iudyn
#ifdef __PARA
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
#endif
!
implicit none
!
......@@ -42,7 +38,7 @@
integer :: isym, invs (48), m1,m2,m3, s(3, 3, 48), sna, irt(48,nat), jsym, &
neig, iwork( 5*nmodes ), info, ifail( nmodes), isq (48), nsq, sym_sgq(48)
!
logical :: timerev, eqvect_strict
logical :: eqvect_strict
! True if we are using time reversal
!
real(kind=DP) :: arg, rwork( 7*nmodes ), aq(3), saq(3), raq(3)
......@@ -423,6 +419,8 @@
saq = xq
call cryst_to_cart (1, aq, at, - 1)
CALL cryst_to_cart (1, saq, at, -1)
! Initialize isym
isym = 1
DO jsym=1, nsq
ism1 = invs (sym_sgq(jsym))
raq = 0.d0
......
......@@ -24,8 +24,7 @@
USE noncollin_module,ONLY : npol
#ifdef __PARA
USE mp_global,ONLY : nproc_pool, me_pool
USE mp_global,ONLY : npool,my_pool_id
USE mp_world, ONLY : mpime
USE mp_global,ONLY : npool
#endif
!
implicit none
......
......@@ -9,16 +9,16 @@
!-----------------------------------------------------------------------
SUBROUTINE rgd_blk (nr1,nr2,nr3,nat,dyn,q,tau,epsil,zeu,bg,omega,sign)
!-----------------------------------------------------------------------
! This is adapted from QE PH/rigid.f90
!
! compute the rigid-ion (long-range) term for q
! The long-range term used here, to be added to or subtracted from the
! dynamical matrices, is exactly the same of the formula introduced
! in:
! X. Gonze et al, PRB 50. 13035 (1994) . Only the G-space term is
! implemented: the Ewald parameter alpha must be large enough to
! have negligible r-space contribution
!
!! This is adapted from QE PH/rigid.f90
!!
!! compute the rigid-ion (long-range) term for q
!! The long-range term used here, to be added to or subtracted from the
!! dynamical matrices, is exactly the same of the formula introduced
!! in:
!! X. Gonze et al, PRB 50. 13035 (1994) . Only the G-space term is
!! implemented: the Ewald parameter alpha must be large enough to
!! have negligible r-space contribution
!!
USE kinds, only: dp
USE constants_epw, only: pi, fpi, e2
implicit none
......@@ -164,7 +164,7 @@ SUBROUTINE rgd_blk_epw(q,uq,epmat,nmodes,epsil,zeu,bmat,sign,aaa)
!
USE kinds, only: dp
USE pwcom, ONLY : at, bg, omega, alat
USE ions_base, ONLY : amass, tau, nat, ntyp => nsp, ityp
USE ions_base, ONLY : tau, nat
USE constants_epw, only: twopi, fpi, e2, ci, czero, cone, two, ryd2mev
!
implicit none
......@@ -189,7 +189,7 @@ SUBROUTINE rgd_blk_epw(q,uq,epmat,nmodes,epsil,zeu,bmat,sign,aaa)
qeq, &! <q+G| epsil | q+G>
arg, zaq, &
g1, g2, g3, gmax, alph, geg
integer :: na, nb, ipol, im, m1, m2, m3
integer :: na, ipol, im, m1, m2, m3
complex(dp) :: fac, facqd, facq, aaa(nmodes)
!
IF (abs(sign) /= 1.0) &
......@@ -244,8 +244,8 @@ SUBROUTINE rgd_blk_epw2(nr1,nr2,nr3,q,uq,epmat,nmodes,epsil,zeu,bmat,sign)
! to be added or subtracted from the vertex
!
USE kinds, only: dp
USE pwcom, ONLY : at, bg, omega, alat
USE ions_base, ONLY : amass, tau, nat, ntyp => nsp, ityp
USE pwcom, ONLY : bg, omega, alat
USE ions_base, ONLY : tau, nat
USE constants_epw, only: twopi, fpi, e2, ci, czero, cone, two, ryd2mev
!
implicit none
......@@ -270,8 +270,8 @@ SUBROUTINE rgd_blk_epw2(nr1,nr2,nr3,q,uq,epmat,nmodes,epsil,zeu,bmat,sign)
real(DP) :: &
qeq, &! <q+G| epsil | q+G>
arg, zaq, &
g1, g2, g3, gmax, alph, geg, aaa1
integer :: na, nb, ipol, im, m1,m2,m3, nrx1,nrx2,nrx3
g1, g2, g3, gmax, alph, geg
integer :: na, ipol, im, m1,m2,m3, nrx1,nrx2,nrx3
complex(dp) :: fac, facqd, facq, aaa(nmodes)
!
IF (abs(sign) /= 1.0) &
......@@ -344,21 +344,20 @@ END SUBROUTINE rgd_blk_epw2
!
!
!-----------------------------------------------------------------------------
SUBROUTINE dyndia_epw (nmodes, xq, dyn, uq, wq)
SUBROUTINE dyndia_epw (nmodes, dyn, uq, wq)
!-----------------------------------------------------------------------------
!
!!
!! diagonalize dyn mat --> uq, wq
!!
USE kinds, ONLY : DP
USE pwcom, ONLY : at, bg, celldm, omega
USE ions_base, ONLY : amass, tau, nat, ntyp => nsp, ityp
USE ions_base, ONLY : amass, nat, ityp
USE constants_epw, ONLY : twopi, ci, czero, rydcm1
USE io_global, ONLY : stdout
!
implicit none
!
! input variables
integer :: nmodes
! number of branches
real(kind=DP) :: xq(3)
complex(kind=DP) :: dyn (nmodes, nmodes)
! dynamical matrix in bloch representation (Cartesian coordinates)
!
......@@ -371,11 +370,7 @@ SUBROUTINE dyndia_epw (nmodes, xq, dyn, uq, wq)
complex(kind=DP) :: cwork( 2*nmodes ), u1(nmodes,nmodes)
integer :: neig, info, ifail( nmodes ), iwork( 5*nmodes )
real(kind=DP) :: w1(nmodes), rwork( 7*nmodes ), massfac
integer :: imode, jmode, na, nb, i, mu, nu
!
! ------------------------------------------------------------------
! diagonalize dyn mat --> uq, wq
! ------------------------------------------------------------------
integer :: imode, jmode, na, nb, mu, nu
!
! divide by the square root of masses
!
......@@ -426,18 +421,17 @@ END SUBROUTINE dyndia_epw
!
!
!-----------------------------------------------------------------------
SUBROUTINE compute_bmn_para2 (nbnd, nkstot, cufkk, cufkq, bmatf)
SUBROUTINE compute_bmn_para2 (nbnd, cufkk, cufkq, bmatf)
!-----------------------------------------------------------------------
!
! calculates overlap U_k+q U_k^\dagger
!
!!
!! calculates overlap U_k+q U_k^\dagger
!!
USE kinds, ONLY : DP
USE constants_epw, ONLY : twopi, ci, czero, cone
USE mp_global, ONLY : my_pool_id
implicit none
!
! input variables
integer :: nbnd, nkstot
integer :: nbnd
! number of bands (possibly in the optimal subspace)
! number of kpoint blocks, total
complex(kind=DP) :: cufkk (nbnd, nbnd), cufkq (nbnd, nbnd)
......@@ -448,9 +442,6 @@ SUBROUTINE compute_bmn_para2 (nbnd, nkstot, cufkk, cufkq, bmatf)
complex(kind=DP) :: bmatf (nbnd, nbnd)
! overlap wfcs in Bloch representation, fine grid
!
! work variables
complex(kind=DP) :: btmp( nbnd, nbnd)
!
! every pool works with its own subset of k points on the fine grid
!
bmatf = czero
......@@ -468,29 +459,30 @@ END SUBROUTINE compute_bmn_para2
!-----------------------------------------------------------------------
SUBROUTINE compute_bmn_para3 (nbnd, nbndsub, nks, cuk, cukq, bmat)
!-----------------------------------------------------------------------
!
! calculates <u_mk+q|e^iGr|u_nk> in the approximation q+G->0
!
!!
!! calculates <u_mk+q|e^iGr|u_nk> in the approximation q+G->0
!!
USE kinds, ONLY : DP
USE constants_epw, ONLY : twopi, ci, czero, cone
USE mp_global, ONLY : my_pool_id
USE io_global, ONLY : stdout
implicit none
!
! input variables
integer :: nbnd, nbndsub, nks
! number of bands
! number of kpoint blocks, per pool
complex(kind=DP) :: cuk (nbnd,nbndsub,nks), cukq (nbnd,nbndsub,nks)
! rotation matrix U(k)
! rotation matrix U(k+q)
!
! output variables
complex(kind=DP) :: bmat (nbnd, nbnd, nks)
! overlap wfcs in Bloch representation, fine grid
INTEGER, INTENT (in) :: nbnd
!! number of bands
INTEGER, INTENT (in) :: nks
!! number of kpoint blocks, per pool
INTEGER, INTENT (in) :: nbndsub
!! Number of band on the subspace of Wannier
!
COMPLEX(kind=DP), INTENT (in) :: cuk (nbnd,nbndsub,nks)
!! rotation matrix U(k)
COMPLEX(kind=DP), INTENT (in) :: cukq (nbnd,nbndsub,nks)
!! rotation matrix U(k+q)
COMPLEX(kind=DP), INTENT (out) :: bmat (nbnd, nbnd, nks)
!! overlap wfcs in Bloch representation, fine grid
!
! work variables
integer :: ik,ibnd,jbnd
integer :: ik
!
! every pool works with its own subset of k points on the fine grid
!
......
......@@ -8,36 +8,36 @@
!
!
!--------------------------------------------------------------------------
subroutine rotate_eigenm ( iq_first, iq, nqc, isym, nsym, s, invs, irt, &
rtau, xq, isq, cz1, cz2)
subroutine rotate_eigenm ( iq_first, nqc, isym, s, invs, irt, &
rtau, xq, cz1, cz2)
!--------------------------------------------------------------------------
!
! Here:
!
! 1). determine the unitary matrix which gives the transformed eigenmodes
! from the first q in the star to the present q
!
! 2). rotate eigenmodes from the first q in the star (cz1) to the current
! q (cz2)
!
! In rotate_epmat.f90:
!
! 3). rotate the electron-phonon matrix from the cartesian representation
! of the first qpoint of the star to the eigenmode representation
! (using cz1).
!
! 4). rotate the electron-phonon matrix from the eigenmode representation
! to the cartesian representation of the qpoint iq in the star (with cz2).
!
! Step 3 requires using the rotated eigenmodes to set the phases
! (the diagonalizers of the rotated dyn mat will not
! work because of the gages in deg. subspaces and the phases).
!
! W.r.t. the standard method of q2qstar_ph.f90, here I construct the
! unitary matrix Gamma defined in Maradudin & Vosko, RMP 1968.
! In this way all rotations are just zgemm operations and we throw away
! all nasty indices.
!
!!
!! Here:
!!
!! 1). determine the unitary matrix which gives the transformed eigenmodes
!! from the first q in the star to the present q
!!
!! 2). rotate eigenmodes from the first q in the star (cz1) to the current
!! q (cz2)
!!
!! In rotate_epmat.f90:
!!
!! 3). rotate the electron-phonon matrix from the cartesian representation
!! of the first qpoint of the star to the eigenmode representation
!! (using cz1).
!!
!! 4). rotate the electron-phonon matrix from the eigenmode representation
!! to the cartesian representation of the qpoint iq in the star (with cz2).
!!
!! Step 3 requires using the rotated eigenmodes to set the phases
!! (the diagonalizers of the rotated dyn mat will not
!! work because of the gages in deg. subspaces and the phases).
!!
!! W.r.t. the standard method of q2qstar_ph.f90, here I construct the
!! unitary matrix Gamma defined in Maradudin & Vosko, RMP 1968.
!! In this way all rotations are just zgemm operations and we throw away
!! all nasty indices.
!!
!
!--------------------------------------------------------------------------
USE kinds, ONLY : DP
......@@ -50,7 +50,7 @@
use ions_base, ONLY : nat, amass, nat, ityp
implicit none
!
integer :: iq_first, nqc, isym, isq(48), nsym, iq
integer :: iq_first, nqc, isym
! originating q point of the star
! total number of qpoints
! fractional translation
......
......@@ -26,7 +26,6 @@
USE constants_epw, ONLY : cone, czero, ryd2mev
USE pwcom, ONLY : nbnd, nks
USE ions_base, ONLY : amass, ityp
USE modes, ONLY : nirr, npert
USE phcom, ONLY : nq1, nq2, nq3
implicit none
!
......
......@@ -34,23 +34,49 @@
USE io_epw, ONLY : iunepmatf, linewidth_elself
USE phcom, ONLY : nmodes
USE epwcom, ONLY : nbndsub, lrepmatf, &
fsthick, eptemp, ngaussw, degaussw, &
etf_mem, eps_acustic, efermi_read, fermi_energy
fsthick, eptemp, ngaussw, degaussw, &
etf_mem, eps_acustic, efermi_read, fermi_energy
USE pwcom, ONLY : ef !, nelec, isk
USE elph2, ONLY : etf, ibndmin, ibndmax, nkqf, &
epf17, wkf, nkf, nqtotf, wf, wqf, xkf, nkqtotf, &
sigmar_all, sigmai_all, sigmai_mode, zi_all, efnew, nqf
nkf, epf17, wkf, nqtotf, wf, wqf, xkf, nkqtotf, &
sigmar_all, sigmai_all, sigmai_mode, zi_all, efnew
USE control_flags, ONLY : iverbosity
USE constants_epw, ONLY : ryd2mev, one, ryd2ev, two, zero, pi, ci, eps6
#ifdef __PARA
USE mp, ONLY : mp_barrier, mp_sum
USE mp_global, ONLY : me_pool, inter_pool_comm, my_pool_id
USE mp_global, ONLY : inter_pool_comm
#endif
implicit none
!
INTEGER, INTENT (in) :: iq
!! Current q-point index
!
! Local variables
INTEGER :: n
!! Integer for the degenerate average over eigenstates
integer :: ik, ikk, ikq, ibnd, jbnd, imode, nrec, iq, fermicount
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 :: imode
!! Counter on mode
INTEGER :: nrec
!! Record index for reading the e-f matrix
INTEGER :: fermicount
!! Number of states on the Fermi surface
INTEGER :: nksqtotf
!! Total number of k+q points
INTEGER :: lower_bnd
!! Lower bounds index after k or q paral
INTEGER :: upper_bnd
!! Upper bounds index after k or q paral
!
REAL(kind=DP) :: tmp
!! Temporary variable to store real part of Sigma for the degenerate average
REAL(kind=DP) :: tmp2
......@@ -65,17 +91,54 @@
!! Temporary array to store the imag-part of Sigma
REAL(kind=DP) :: zi_tmp(ibndmax-ibndmin+1)
!! Temporary array to store the Z
REAL(kind=DP) :: g2, ekk, ekq, wq, ef0, wgq, wgkq, weight, dosef, &
w0g1, w0g2, inv_wq, inv_eptemp0, g2_tmp,&
inv_degaussw
REAL(kind=DP), external :: efermig, dos_ef, wgauss, w0gauss
complex(kind=DP) epf (ibndmax-ibndmin+1, ibndmax-ibndmin+1)
!
! variables for collecting data from all pools in parallel case
!
integer :: nksqtotf, lower_bnd, upper_bnd
REAL(kind=DP), ALLOCATABLE :: xkf_all(:,:), etf_all(:,:)
REAL(kind=DP) :: g2
!! Electron-phonon matrix elements squared in Ry^2
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) :: wq
!! Phonon frequency on the fine grid
REAL(kind=DP) :: ef0
!! Fermi energy level
REAL(kind=DP) :: wgq
!! Bose occupation factor $n_{q\nu}(T)$
REAL(kind=DP) :: wgkq
!! Fermi-Dirac occupation factor $f_{nk+q}(T)$
REAL(kind=DP) :: weight
!! Self-energy factor
!!$$ N_q \Re( \frac{f_{mk+q}(T) + n_{q\nu}(T)}{ \varepsilon_{nk} - \varepsilon_{mk+q} + \omega_{q\nu} - i\delta }) $$
!!$$ + N_q \Re( \frac{1- f_{mk+q}(T) + n_{q\nu}(T)}{ \varepsilon_{nk} - \varepsilon_{mk+q} - \omega_{q\nu} - i\delta }) $$
REAL(kind=DP) :: dosef
!! Density of state N(Ef)
REAL(kind=DP) :: w0g1
!! Dirac delta for the imaginary part of $\Sigma$
REAL(kind=DP) :: w0g2
!! Dirac delta for the imaginary part of $\Sigma$
REAL(kind=DP) :: inv_wq
!! $frac{1}{2\omega_{q\nu}}$ defined for efficiency reasons
REAL(kind=DP) :: inv_eptemp0
!! Inverse of temperature define for efficiency reasons
REAL(kind=DP) :: g2_tmp
!! If the phonon frequency is too small discart g
REAL(kind=DP) :: inv_degaussw
!! Inverse of the smearing for efficiency reasons
!REAL(kind=DP), external :: efermig
!! Function to compute the Fermi energy
REAL(kind=DP), external :: dos_ef
!! Function to compute the Density of States at the Fermi level
REAL(kind=DP), external :: wgauss
!! Fermi-Dirac distribution function (when -99)
REAL(kind=DP), external :: w0gauss
!! This function computes the derivative of the Fermi-Dirac function
!! It is therefore an approximation for a delta function
REAL(kind=DP), ALLOCATABLE :: xkf_all(:,:)
!! Collect k-point coordinate from all pools in parallel case
REAL(kind=DP), ALLOCATABLE :: etf_all(:,:)
!! Collect eigenenergies from all pools in parallel case
!
COMPLEX(kind=DP) epf (ibndmax-ibndmin+1, ibndmax-ibndmin+1)
!! Electron-phonon matrix element on the fine grid.
!
! SP: Define the inverse so that we can efficiently multiply instead of
! dividing
......@@ -403,7 +466,6 @@
100 FORMAT(5x,'Gaussian Broadening: ',f10.6,' eV, ngauss=',i4)
101 FORMAT(5x,'DOS =',f10.6,' states/spin/eV/Unit Cell at Ef=',f10.6,' eV')
102 FORMAT(5x,'E( ',i3,' )=',f9.4,' eV Re[Sigma]=',f15.6,' meV Im[Sigma]=',f15.6,' meV Z=',f15.6,' lam=',f15.6)
103 FORMAT(5x,'k( ',i7,' )=',f9.4,' eV Re[Sigma]=',f15.6,' meV Im[Sigma]=',f15.6,' meV Z=',f15.6)
!
RETURN
!
......@@ -413,15 +475,15 @@
!-----------------------------------------------------------------------
SUBROUTINE selfen_elec_k ( ik )
!-----------------------------------------------------------------------
!
! Compute the imaginary part of the electron self energy due to electron-
! phonon interaction in the Migdal approximation. This corresponds to
! the electron 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
!
!!
!! Compute the imaginary part of the electron self energy due to electron-
!! phonon interaction in the Migdal approximation. This corresponds to
!! the electron 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
......@@ -432,13 +494,13 @@
etf_mem, eps_acustic, efermi_read, fermi_energy
USE pwcom, ONLY : ef !, nelec, isk
USE elph2, ONLY : etf, ibndmin, ibndmax, nkqf, etf_k, &
epf17, wkf, nkf, nqtotf, wf, wqf, xkf, nkqtotf, &
epf17, wkf, nqtotf, wf, wqf, xkf, nkqtotf, &
sigmar_all, sigmai_all, sigmai_mode, zi_all, efnew, nqf
USE constants_epw, ONLY : ryd2mev, one, ryd2ev, two, zero, pi, ci, eps6
USE control_flags, ONLY : iverbosity