Commit 42387d94 authored by sponce's avatar sponce

Adding some screening capabilities

Also cumulant post-processing.
Courtesy of C. Verdi and F. Caruso.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@13723 c92efa57-630b-4861-b058-cf58834340f0
parent ef5b32c4
......@@ -104,6 +104,7 @@ set_ndnmbr.o \
setphases_wrap.o \
sgama2.o \
sort.o \
spec_cumulant.o \
spectral_func.o \
spectral_func_ph.o \
spectral_func_pl.o \
......
......@@ -25,7 +25,7 @@
USE noncollin_module, ONLY : noncolin, npol
USE wavefunctions_module, ONLY: evc
USE spin_orb, ONLY : lspinorb
USE control_lr, ONLY : lgamma, nbnd_occ
USE control_lr, ONLY : nbnd_occ
USE phcom, ONLY : evq, dpsi, vlocq,&
dmuxc, npertx
USE phus, ONLY : int1, int1_nc, int2, int2_so, &
......@@ -54,19 +54,8 @@
!
! allocate space for the quantities needed in EPW
!
IF (lgamma) THEN
!
! q=0 : evq is a pointers to evc
!
evq => evc
ELSE
!
! q!=0 : evq is ALLOCATEd and calculated at point k+q
!
ALLOCATE (evq ( npwx*npol, nbnd))
ENDIF
!
ALLOCATE ( dpsi ( npwx*npol, nbnd))
ALLOCATE (evq ( npwx*npol, nbnd))
ALLOCATE (dpsi ( npwx*npol, nbnd))
ALLOCATE (vlocq ( ngm, ntyp))
! SP: nrxx is not used in QE 5 ==> tg_nnr is the maximum among nnr
! This SHOULD have the same dim as nrxx had.
......
......@@ -36,8 +36,9 @@
muc, mp_mesh_q, mp_mesh_k, max_memlt, lunif, &
lreal, lpolar, lpade, liso, limag, laniso, &
specfun_el, specfun_ph, lifc, asr_typ, &
lscreen, scr_typ, fermi_diff, smear_rpa, &
rand_q, rand_nq, rand_nk, rand_k, pwc, phonselfen, &
parallel_q, parallel_k, specfun_pl, &
parallel_q, parallel_k, specfun_pl, cumulant, bnd_cum, &
nw_specfun, nw, nswi, nswfc, nswc, nstemp, nsmear, &
wsfc, wscut, write_wfn, wmin_specfun, wmin, &
wmax_specfun, wmax, wepexst, wannierize, &
......@@ -52,7 +53,6 @@
USE mp_world, ONLY : world_comm
USE io_files, ONLY : prefix, tmp_dir
USE qpoint, ONLY : xq
USE control_lr, ONLY : lgamma
USE io_global, ONLY : meta_ionode_id
USE control_flags, ONLY : iverbosity
USE ions_base, ONLY : amass
......@@ -61,7 +61,6 @@
!
! logicals
!
CALL mp_bcast (lgamma, meta_ionode_id, world_comm)
CALL mp_bcast (epsil, meta_ionode_id, world_comm)
CALL mp_bcast (trans, meta_ionode_id, world_comm)
CALL mp_bcast (zue, meta_ionode_id, world_comm)
......@@ -108,6 +107,8 @@
CALL mp_bcast (laniso, meta_ionode_id, world_comm) !
CALL mp_bcast (lpolar, meta_ionode_id, world_comm) !
CALL mp_bcast (lifc, meta_ionode_id, world_comm)
CALL mp_bcast (lscreen, meta_ionode_id, world_comm)
CALL mp_bcast (cumulant, meta_ionode_id, world_comm)
CALL mp_bcast (lunif, meta_ionode_id, world_comm) !
CALL mp_bcast (kerwrite, meta_ionode_id, world_comm) !
CALL mp_bcast (kerread, meta_ionode_id, world_comm) !
......@@ -158,6 +159,8 @@
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 (scr_typ, meta_ionode_id, world_comm)
CALL mp_bcast (bnd_cum, meta_ionode_id, world_comm)
!
! real*8
!
......@@ -190,6 +193,8 @@
CALL mp_bcast (nel, meta_ionode_id, world_comm)
CALL mp_bcast (meff, meta_ionode_id, world_comm)
CALL mp_bcast (epsiHEG, meta_ionode_id, world_comm)
CALL mp_bcast (fermi_diff, meta_ionode_id, world_comm)
CALL mp_bcast (smear_rpa, meta_ionode_id, world_comm)
!
! characters
!
......
......@@ -26,8 +26,8 @@
REAL(DP), PARAMETER :: two = 2.d0
REAL(DP), PARAMETER :: zero = 0.d0
REAL(DP), PARAMETER :: e2 = 2.0_DP ! the square of the electron charge
COMPLEX(DP), PARAMETER :: ci = (0.d0, 1.d0)
COMPLEX(DP), PARAMETER :: cone = (1.d0, 0.d0)
COMPLEX(DP), PARAMETER :: ci = (0.d0, 1.d0)
COMPLEX(DP), PARAMETER :: cone = (1.d0, 0.d0)
COMPLEX(DP), PARAMETER :: czero = (0.d0, 0.d0)
!
! Unit conversion factors
......@@ -37,7 +37,8 @@
REAL(DP), PARAMETER :: bohr = 0.52917721092d0
REAL(DP), PARAMETER :: ryd2mev = 13605.6981d0
REAL(DP), PARAMETER :: ryd2ev = 13.6056981d0
REAL(DP), PARAMETER :: rydcm1 = 13.6056981d0 * 8065.541d0
REAL(DP), PARAMETER :: ha2ev = 2.d0*ryd2ev
REAL(DP), PARAMETER :: rydcm1 = ryd2ev * 8065.541d0
REAL(DP), PARAMETER :: bohr2ang = 0.52917721092d0
REAL(DP), PARAMETER :: ev2cmm1 = 8065.541d0
REAL(DP), PARAMETER :: kelvin2eV= 8.6173427909d-05
......
......@@ -29,11 +29,11 @@
USE phus, ONLY : int1, int1_nc, int2, int4, int4_nc
USE lr_symm_base, ONLY : rtau
USE noncollin_module, ONLY : m_loc
USE control_lr, ONLY : lgamma, nbnd_occ
USE control_lr, ONLY : nbnd_occ
USE becmod, ONLY : becp, deallocate_bec_type
USE elph2, ONLY : el_ph_mat, epf17, epsi, etf,&
etq, et_all, wf, wkf, wqf, wslen,&
xkq, xk_all, zstar, xkf, xqf, epmatwp
xkq, xk_all, zstar, xkf, xqf, epmatwp, eps_rpa
USE epwcom, ONLY : epbread, epwread
USE modes, ONLY : t, npert, u, name_rap_mode, num_rap_mode
USE qpoint, ONLY : eigqts, igkq
......@@ -59,16 +59,12 @@
IF(ALLOCATED(xk_all)) DEALLOCATE (xk_all)
IF(ALLOCATED(et_all)) DEALLOCATE (et_all)
IF(ALLOCATED(wslen)) DEALLOCATE (wslen)
IF(ALLOCATED(eps_rpa)) DEALLOCATE (eps_rpa)
!
ELSE
!
IF (lgamma) THEN
IF(ASSOCIATED(evq)) NULLIFY(evq)
IF(ASSOCIATED(igkq)) NULLIFY(igkq)
ELSE
IF(ASSOCIATED(evq)) DEALLOCATE(evq)
IF(ASSOCIATED(igkq)) DEALLOCATE(igkq)
END IF
IF(ASSOCIATED(evq)) DEALLOCATE(evq)
IF(ASSOCIATED(igkq)) DEALLOCATE(igkq)
!
IF(ALLOCATED(dvpsi)) DEALLOCATE (dvpsi)
IF(ALLOCATED(dpsi)) DEALLOCATE ( dpsi)
......@@ -139,6 +135,7 @@
IF(ALLOCATED(xk_all)) DEALLOCATE (xk_all)
IF(ALLOCATED(et_all)) DEALLOCATE (et_all)
IF(ALLOCATED(wslen)) DEALLOCATE (wslen)
IF(ALLOCATED(eps_rpa)) DEALLOCATE (eps_rpa)
ENDIF ! epwread .and. .not. epbread
!
END SUBROUTINE deallocate_epw
......@@ -29,11 +29,12 @@
umatq_all(:,:,:), &! the rotation matrix for the unique setting of the wfs gauge -- for all k+q points
dynq (:,:,:), &! dynamical matrix for every q (nmode, nmodes, nqtot)
epmatq (:,:,:,:,:), &! e-p matrix for every q (nbnd, nbnd, nks, nmodes, nqtot)
epf17 (:, :, :, :), &! full ep matrix in bloch rep stored in mem (nkqtotf, nbnd, nbnd, nmodes)-nbnd inside wndw
epf17 (:, :, :, :), &! full ep matrix in bloch rep stored in mem (nkqtotf, nbnd, nbnd, nmodes)-nbnd inside wndw
dmec(:,:,:,:), &! dipole matrix elements on the coarse mesh (ipol, nbnd, nbnd, nks)
dmef(:,:,:,:), &! dipole matrix elements on the fine mesh (ipol, nbnd, nbnd, nks)
vmef(:,:,:,:), &! velocity matrix elements on the fine mesh (ipol, nbnd, nbnd, nks)
bmat(:,:,:,:) ! overlap U_k+q U_k^\dagger on the coarse mesh (nbnd, nbnd, nks, nqtot)
bmat(:,:,:,:), &! overlap U_k+q U_k^\dagger on the coarse mesh (nbnd, nbnd, nks, nqtot)
eps_rpa(:) ! screening
REAL(KIND=DP), ALLOCATABLE ::&
a_all(:,:), &! electronic spectral function du to electron-phonon interaction
xk_all(:,:), &! full k point grid, coarse (3, nkstot)
......
......@@ -45,7 +45,6 @@
USE start_k, ONLY : nk1, nk2, nk3
USE phcom, ONLY : dpsi, dvpsi, evq, nq1, nq3, nq2
USE qpoint, ONLY : igkq
USE control_lr, ONLY : lgamma
USE qpoint, ONLY : xq
USE modes, ONLY : nmodes
USE lr_symm_base, ONLY : minus_q, rtau, gi, gimq, irotmq, nsymq, invsymq
......@@ -188,8 +187,6 @@
!
CALL start_clock ( 'elphon_wrap' )
!
IF (lgamma) CALL errore('elphon_shuffle_wrap','EPW does not support Gamma only calculation ',1)
!
! READ qpoint list from stdin
!
IF (meta_ionode) READ(5,*) nqc_irr
......
......@@ -32,8 +32,8 @@
USE ions_base, ONLY : nat, amass, ityp, tau
USE phcom, ONLY : nq1, nq2, nq3, nmodes
USE epwcom, ONLY : nbndsub, lrepmatf, fsthick, epwread, longrange, &
epwwrite, ngaussw, degaussw, lpolar, lifc, &
nbndskip, parallel_k, parallel_q, etf_mem, &
epwwrite, ngaussw, degaussw, lpolar, lifc, lscreen, &
nbndskip, parallel_k, parallel_q, etf_mem, scr_typ, &
elecselfen, phonselfen, nest_fn, a2f, specfun_ph, &
vme, eig_read, ephwrite, nkf1, nkf2, nkf3, &
efermi_read, fermi_energy, specfun_el, band_plot, &
......@@ -51,7 +51,7 @@
wkf, dynq, nqtotf, nkqf, epf17, nkf, nqf, et_ks, &
ibndmin, ibndmax, lambda_all, dmec, dmef, vmef, &
sigmai_all, sigmai_mode, gamma_all, epsi, zstar, &
efnew, ifc, sigmar_all, zi_all, nkqtotf
efnew, ifc, sigmar_all, zi_all, nkqtotf, eps_rpa
#if defined(__NAG)
USE f90_unix_io, ONLY : flush
#endif
......@@ -463,7 +463,7 @@
etf_ks ( nbndsub, nkqf), &
epmatf( nbndsub, nbndsub, nmodes), cufkk ( nbndsub, nbndsub), &
cufkq ( nbndsub, nbndsub), uf ( nmodes, nmodes), &
bmatf( nbndsub, nbndsub) )
bmatf( nbndsub, nbndsub), eps_rpa( nmodes) )
!
! Need to be initialized
epmatf(:,:,:) = czero
......@@ -485,7 +485,7 @@
! SP: Create a look-up table for the exponential of the factor.
! This can only work with homogeneous fine grids.
IF ( (nkf1 >0) .AND. (nkf2 > 0) .AND. (nkf3 > 0) .AND. &
(nqf1 >0) .AND. (nqf2 > 0) .AND. (nqf3 > 0) .AND. .NOT. mp_mesh_k ) THEN
(nqf1 >0) .AND. (nqf2 > 0) .AND. (nqf3 > 0) .AND. .NOT. mp_mesh_k .AND. .NOT. lscreen ) THEN
! Make a check
IF ((nqf1>nkf1) .or. (nqf2>nkf2) .or. (nqf3>nkf3)) &
CALL errore('The fine q-grid cannot be larger than the fine k-grid',1)
......@@ -770,6 +770,11 @@
! number of k points with a band on the Fermi surface
fermicount = 0
!
IF (lscreen) THEN
IF (scr_typ == 0) CALL rpa_epsilon (xxq, wf(:,iq), nmodes, epsi, eps_rpa)
IF (scr_typ == 1) CALL tf_epsilon (xxq, nmodes, epsi, eps_rpa)
ENDIF
!
! this is a loop over k blocks in the pool
! (size of the local k-set)
DO ik = 1, nkf
......@@ -919,7 +924,11 @@
DO jbnd = ibndmin, ibndmax
DO ibnd = ibndmin, ibndmax
!
epf17(ibnd-ibndmin+1,jbnd-ibndmin+1,:,ik) = epmatf(ibnd,jbnd,:)
IF (lscreen) THEN
epf17(ibnd-ibndmin+1,jbnd-ibndmin+1,:,ik) = epmatf(ibnd,jbnd,:) / eps_rpa(:)
ELSE
epf17(ibnd-ibndmin+1,jbnd-ibndmin+1,:,ik) = epmatf(ibnd,jbnd,:)
ENDIF
!
ENDDO
ENDDO
......
......@@ -20,14 +20,14 @@
!! @Note
!! 8/14/08 lnscf is unnecessary, as is nqs,iq_start
!!
USE io_global, ONLY : stdout
USE io_global, ONLY : stdout, ionode
USE mp, ONLY : mp_bcast, mp_barrier
USE mp_world, ONLY : mpime
USE mp_global, ONLY : mp_startup, ionode_id, mp_global_end
USE control_flags, ONLY : gamma_only
USE control_epw, ONLY : wannierize
USE global_version, ONLY : version_number
USE epwcom, ONLY : filukk, eliashberg, ep_coupling, epwread, epbread
USE epwcom, ONLY : filukk, eliashberg, ep_coupling, epwread, epbread, cumulant
USE environment, ONLY : environment_start
USE elph2, ONLY : elph
! Flag to perform an electron-phonon calculation. If .true.
......@@ -161,6 +161,10 @@ write(stdout,'(a)') "
CALL close_epw()
!
ENDIF
!
IF ( cumulant .and. ionode ) THEN
CALL spec_cumulant()
ENDIF
!
IF ( eliashberg ) THEN
CALL eliashberg_eqs()
......
......@@ -48,7 +48,8 @@
fermi_energy, efermi_read, max_memlt, fila2f, &
ep_coupling, nw_specfun, wmax_specfun, &
wmin_specfun, laniso, lpolar, lifc, asr_typ, &
proj, write_wfn, iswitch, ntempxx, &
lscreen, scr_typ, fermi_diff, smear_rpa, &
cumulant, bnd_cum, proj, write_wfn, iswitch, ntempxx, &
liso, lacon, lpade, etf_mem, epbwrite, &
nsiter, conv_thr_racon, specfun_el, specfun_ph, &
pwc, nswc, nswfc, nswi, filukq, filukk, &
......@@ -66,7 +67,6 @@
USE mp_world, ONLY : world_comm
USE partial, ONLY : atomo, nat_todo
USE constants, ONLY : AMU_RY
USE control_lr, ONLY : lgamma
USE mp_global, ONLY : my_pool_id, me_pool
USE io_global, ONLY : meta_ionode, meta_ionode_id, ionode, ionode_id, stdout
#if defined(__NAG)
......@@ -115,6 +115,7 @@
broyden_beta, broyden_ndim, nstemp, tempsmin, tempsmax, temps, &
conv_thr_raxis, conv_thr_iaxis, conv_thr_racon, &
gap_edge, nsiter, muc, lreal, limag, lpade, lacon, liso, laniso, lpolar,&
lscreen, scr_typ, fermi_diff, smear_rpa, cumulant, bnd_cum, &
lifc, asr_typ, lunif, kerwrite, kerread, imag_read, eliashberg, &
ep_coupling, fila2f, max_memlt, efermi_read, fermi_energy, &
specfun_el, specfun_ph, wmin_specfun, wmax_specfun, nw_specfun, &
......@@ -231,9 +232,18 @@
! delta_approx : if .true. the double delta approximation is used to compute the phonon self-energy
!
! added by CV & SP
! lpolar : if .true. enable the correct Wannier interpolation in the case of polar material.
! lifc : if .true. reads interatomic force constants produced by q2r.x for phonon interpolation
! lpolar : if .true. enable the correct Wannier interpolation in the case of polar material.
! lifc : if .true. reads interatomic force constants produced by q2r.x for phonon interpolation
! asr_typ : select type of ASR if lifc=.true. (as in matdyn); otherwise it is the usual simple sum rule
! lscreen : if .true. the e-ph matrix elements are screened by the RPA or TF dielectric function
! scr_typ : if 0 calculates the Lindhard screening, if 1 the Thomas-Fermi screening
! fermi_diff : difference between Fermi energy and band edge (in eV)
! smear_rpa : smearing for the calculation of the Lindhard function (in eV)
! cumulant : if .true. calculates the electron spectral function using the cumulant expansion method
! (can be used as independent postprocessing by setting ep_coupling=.false.)
! bnd_cum : band index for which the cumulant calculation is done
! (for more than one band, perform multiple calculations and add the results together)
!
!
! Added by SP
!
......@@ -263,8 +273,8 @@
! When etf_mem == 2, an additional loop is done on mode for the fine grid interpolation
! part. This reduces the memory further by a factor "nmodes".
! plselfen : Calculate the electron-plasmon self-energy.
! nel : Carrier concentration
! meff : Density of state effective mass
! nel : Fractional number of electrons in the unit cell
! meff : Density of state effective mass (in unit of the electron mass)
! epsiHEG : Dielectric constant at zero doping
!
CHARACTER (LEN=80) :: input_file
......@@ -400,10 +410,16 @@
lpolar = .false.
lifc = .false.
asr_typ = 'simple'
kerwrite= .false.
kerread = .false.
imag_read = .false.
eliashberg = .false.
lscreen = .false.
scr_typ = 0
fermi_diff = 1.d0
smear_rpa = 0.05d0
cumulant = .false.
bnd_cum = 1
kerwrite = .false.
kerread = .false.
imag_read = .false.
eliashberg = .false.
ep_coupling = .true.
nswfc = 0
nswc = 0
......@@ -442,9 +458,9 @@
longrange = .false.
shortrange = .false.
prtgkk = .false.
nel = 0.01d0
meff = 12.d0
epsiHEG = 0.25d0
nel = 0.0d0
meff = 1.d0
epsiHEG = 1.d0
!
! reading the namelist inputepw
!
......@@ -568,6 +584,12 @@
&'Error: kmaps has to be true for a restart run. ',1)
IF ( etf_mem == 2 .AND. parallel_q) CALL errore('epw_init',&
&'Error: Memory optimized version and q-parallelization not implemented. ',1)
IF ( lscreen .AND. parallel_q) CALL errore('epw_init',&
&'Error: lscreen can only be used with parallel_k ',1)
IF ( lscreen .AND. etf_mem == 2) CALL errore('epw_init',&
&'Error: lscreen not implemented with etf_mem=2 ',1)
IF ( cumulant .AND. parallel_q ) CALL errore('epw_init',&
&'Error: Cumulant and parallel_q is not implemented ',1)
!#ifndef __MPI
! IF ( etf_mem == 2 ) CALL errore('epw_init','Error: etf_mem == 2 only works with MPI.',1)
!#endif
......@@ -586,7 +608,7 @@
! fermi_energy read from the input file
! from eV to Ryd
IF ( efermi_read ) THEN
fermi_energy = fermi_energy / ryd2ev
fermi_energy = fermi_energy / ryd2ev
ENDIF
! eptemp : temperature for the electronic Fermi occupations in the e-p calculation (units of Kelvin)
! 1 K in eV = 8.6173423e-5
......@@ -633,7 +655,6 @@
!
xq(:) = 0.d0
!
lgamma = .false.
tmp_dir = trim(outdir)
dvscf_dir = trim(dvscf_dir)//'/'
!
......
......@@ -39,7 +39,7 @@
USE control_ph, ONLY : lgamma_gamma, search_sym, start_irr, &
last_irr, niter_ph, alpha_mix, &
flmixdpot, reduce_io, u_from_file
USE control_lr, ONLY : lgamma, alpha_pv, nbnd_occ
USE control_lr, ONLY : alpha_pv, nbnd_occ
USE gamma_gamma, ONLY : has_equivalent, asr, nasr, n_diff_sites, &
equiv_atoms, n_equiv_atoms, with_symmetry
USE partial, ONLY : &
......@@ -303,15 +303,6 @@
! ELSE
! xk_col(:,1:nks) = xk(:,1:nks)
! ENDIF
!
! The small group of q may be known. At a given q it is calculated
! by set_nscf, at gamma it coincides with the point group and we
! take nsymq=nsym
!
IF (lgamma.AND.modenum==0) THEN
nsymq=nsym
minus_q=.TRUE.
ENDIF
!
! If the code arrives here and nsymq is still 0 the small group of q has
! not been calculated by set_nscf because this is a recover run.
......@@ -322,8 +313,8 @@
!
!
IF (modenum > 0) THEN
search_sym=.FALSE.
minus_q = .FALSE.
search_sym=.FALSE.
minus_q = .FALSE.
ENDIF
!
! allocate and calculate rtau, the bravais lattice vector associated
......
......@@ -56,6 +56,10 @@
!! Maximum number of wannier functions
INTEGER :: etf_mem
!! If 0, all in memory. If 1, less is stored in memory (read files).
INTEGER :: scr_typ
!! If 0 calculates the Lindhard screening, if 1 the Thomas-Fermi screening
INTEGER :: bnd_cum
!! band index for which the cumulant calculation is done
!
! Superconductivity
INTEGER :: nswfc
......@@ -145,11 +149,15 @@
!
! Plasmon
REAL (KIND=DP) :: nel
!! Carrier concentration
!! fractional number of electrons in the unit cell
REAL (KIND=DP) :: meff
!! Density-of-state effective mass
!! Density-of-state effective mass (in unit of the electron mass)
REAL (KIND=DP) :: epsiHEG
!! Dielectric constant at zero doping.
REAL (KIND=DP) :: fermi_diff
!! difference between Fermi energy and band edge (in eV)
REAL (KIND=DP) :: smear_rpa
!! smearing for the calculation of the Lindhard function (in eV)
!
!LOGICAL :: tphases
!! tphases: if .TRUE. set absolute reference for unitary gauge of the eigenvectors
......@@ -210,8 +218,12 @@
! band_plot : if .true. write filrs to plot band structure and phonon dispersion
LOGICAL :: lpolar
!! if .true. enable the correct Wannier interpolation in the case of polar material.
LOGICAL :: lscreen
!! if .true. the e-ph matrix elements are screened by the RPA or TF dielectric function
LOGICAL :: lifc
!! if .true. reads interatomic force constants produced by q2r.x for phonon interpolation
LOGICAL :: cumulant
!! if .true. calculates the electron spectral function using the cumulant expansion method
LOGICAL :: delta_approx
!! if .true. the double delta approximation is used for the phonon self energy
LOGICAL :: ep_coupling
......
......@@ -24,7 +24,7 @@
PUBLIC :: lambda_phself, linewidth_phself, linewidth_elself, iospectral, &
iua2ffil, iudosfil, iufillambda, iuqdos, iufe, iufilker, &
iufilgap, iospectral_sup, iua2ftrfil, iufilgapFS, iufillambdaFS, &
iuwanep, iuwane, iunukk, iudvscf
iospectral_cum, iuwanep, iuwane, iunukk, iudvscf
PUBLIC :: epwdata, iundmedata, iunvmedata, iunksdata, iudyn, iukgmap, iuepb,&
iufilfreq, iufilegnv, iufileph, iufilkqmap, &
iufilikmap, iueig, iunepmatwp, iunepmatwe, iunkf, iunqf, &
......@@ -56,6 +56,8 @@
INTEGER :: iua2ftrfil = 72 ! Eliashberg transport a2f function [.a2f_tr]
INTEGER :: iufilgapFS = 73 ! Eliashberg superconducting gap on FS with k-points
INTEGER :: iufillambdaFS = 74 ! Electron-phonon coupling strength on FS with k-points
INTEGER :: iospectral_cum = 75 ! Electronic spectral function with the cumulant method
! [specfun_cum##.elself]
!DBSP : iukgmap was 96. Should be the same as set_kplusq.f90.
INTEGER :: iunukk = 77 ! Unit with rotation matrix U(k) from wannier code
INTEGER :: iudvscf = 80 ! Unit for the dvscf_q file
......
......@@ -20,7 +20,7 @@ SUBROUTINE loadqmesh_para
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE epwcom, ONLY : filqf, nqf1, nqf2, nqf3, &
rand_q, rand_nq, mp_mesh_q, system_2d
rand_q, rand_nq, mp_mesh_q, system_2d, lscreen
USE elph2, ONLY : xqf, wqf, nqf, nqtotf
USE cell_base, ONLY : at, bg
USE symm_base, ONLY : s, t_rev, time_reversal, set_sym_bl, nrot
......@@ -40,6 +40,7 @@ SUBROUTINE loadqmesh_para
IF (filqf .ne. '') THEN ! load from file (crystal coordinates)
!
WRITE(stdout, *) ' Using q-mesh file: ', trim(filqf)
IF (lscreen) WRITE(stdout, *) ' WARNING: if lscreen=.true., q-mesh needs to be [-0.5:0.5] (crystal)'
OPEN( unit = iunqf, file = filqf, status = 'old', form = 'formatted',err=100, iostat=ios)
100 CALL errore('loadkmesh_para','opening file '//filqf,abs(ios))
READ(iunqf, *) nqtotf
......@@ -55,7 +56,8 @@ SUBROUTINE loadqmesh_para
!
ELSEIF ( (nqf1.ne.0) .and. (nqf2.ne.0) .and. (nqf3.ne.0) ) THEN ! generate grid
IF (mp_mesh_q) THEN
! get size of the mp_mesh in the irr wedge
IF (lscreen) CALL errore ('If lscreen=.true. do not use mp_mesh_q',1)
! get size of the mp_mesh in the irr wedge
WRITE (stdout, '(a,3i4)') ' Using uniform MP q-mesh: ', nqf1, nqf2, nqf3
call set_sym_bl ( )
!
......@@ -91,6 +93,7 @@ SUBROUTINE loadqmesh_para
ENDDO
ENDDO
ENDDO
IF (lscreen) xqf_(:,:) = xqf_(:,:) - 0.5d0
!
ENDIF
!
......@@ -110,9 +113,11 @@ SUBROUTINE loadqmesh_para
!
IF ( system_2d ) THEN
CALL random_number(xqf_(1:2,iq))
IF (lscreen) xqf_(1:2,iq) = xqf_(1:2,iq) - 0.5d0
xqf_(3,iq) = 0.d0
ELSE
CALL random_number(xqf_(:,iq))
IF (lscreen) xqf_(:,iq) = xqf_(:,iq) - 0.5d0
ENDIF
!
ENDDO
......@@ -192,7 +197,7 @@ SUBROUTINE loadqmesh_serial
USE mp_world, ONLY : mpime
USE io_global, ONLY : stdout
USE epwcom, ONLY : filqf, nqf1, nqf2, nqf3, &
rand_q, rand_nq, mp_mesh_q, system_2d
rand_q, rand_nq, mp_mesh_q, system_2d, lscreen
USE elph2, ONLY : xqf, wqf, nqtotf, nqf
USE cell_base, ONLY : at, bg
USE symm_base, ONLY : s, t_rev, time_reversal, set_sym_bl, nrot
......@@ -207,6 +212,7 @@ SUBROUTINE loadqmesh_serial
! Each pool gets its own copy from the action=read statement
!
WRITE (stdout, *) ' Using q-mesh file: ', trim(filqf)
IF (lscreen) WRITE(stdout, *) ' WARNING: if lscreen=.true., q-mesh needs to be [-0.5:0.5] (crystal)'
OPEN ( unit = iunqf, file = filqf, status = 'old', form = 'formatted', err=100, iostat=ios)
100 CALL errore('loadqmesh_serial','opening file '//filqf,abs(ios))
READ(iunqf, *) nqtotf
......@@ -221,6 +227,7 @@ SUBROUTINE loadqmesh_serial
!
ELSEIF ( (nqf1.ne.0) .and. (nqf2.ne.0) .and. (nqf3.ne.0) ) THEN ! generate grid
IF (mp_mesh_q) THEN
IF (lscreen) CALL errore ('If lscreen=.true. do not use mp_mesh_q',1)
! get size of the mp_mesh in the irr wedge
WRITE (stdout, '(a,3i4)') ' Using uniform q-mesh: ', nqf1, nqf2, nqf3
call set_sym_bl ( )
......@@ -255,6 +262,7 @@ SUBROUTINE loadqmesh_serial
ENDDO
ENDDO
ENDDO
IF (lscreen) xqf(:,:) = xqf(:,:) - 0.5d0
!
!WRITE(stdout,'(a)') ' '
!DO iq = 1, nqtotf
......@@ -278,9 +286,11 @@ SUBROUTINE loadqmesh_serial
!
IF ( system_2d ) THEN
CALL random_number(xqf(1:2,iq))
IF (lscreen) xqf(1:2,iq) = xqf(1:2,iq) - 0.5d0
xqf(3,iq) = 0.d0
ELSE
CALL random_number(xqf(:,iq))
IF (lscreen) xqf(:,iq) = xqf(:,iq) - 0.5d0
ENDIF
!
!WRITE(stdout,'(4f12.7)') xqf(:,iq), wqf(iq)
......
......@@ -127,11 +127,14 @@
m_loc, nqs, lrigid, epsi_, zstar_, lraman, dchi_dtau)
ALLOCATE (dyn(3,3,nat,nat) )
IF (ionode) THEN
DO nt=1, ntyp
IF (amass(nt) <= 0.0d0) amass(nt)=amass2(nt)
END DO
DO nt=1, ntyp
IF (amass(nt) <= 0.0d0) amass(nt)=amass2(nt)
END DO
END IF
!
IF (lpolar .and. .not. lrigid) CALL errore('readmat_shuffle2', &
&'You set lpolar = .true. but did not put epsil = true in the PH calculation at Gamma. ',1)
!
IF (lrigid) THEN
WRITE (6,'(8x,a)') 'Read dielectric tensor and effective charges'
zstar = zstar_
......@@ -347,7 +350,7 @@
read (iudyn,'(a)') line
read (iudyn,'(a)') line
!
IF (lpolar) THEN
!IF (lpolar) THEN
IF (line(6:15).EQ.'Dielectric') THEN
read (iudyn,'(a)') line
read (iudyn,*) ((epsi(i,j), j=1,3), i=1,3)
......@@ -374,10 +377,11 @@
ENDDO
!
ELSE
call errore('readmat_shuffle2','You set lpolar = .true. but did not put epsil = true in the PH calculation at Gamma. ',1)
IF (lpolar) CALL errore('readmat_shuffle2',&
'You set lpolar = .true. but did not put epsil = true in the PH calculation at Gamma. ',1)
ENDIF
!
ENDIF
!ENDIF
!
ENDIF
!
......
......@@ -442,6 +442,209 @@ SUBROUTINE rgd_blk_epw_fine(nq1,nq2,nq3,q,uq,epmat,nmodes,epsil,zeu,bmat,signe)
!
END SUBROUTINE rgd_blk_epw_fine
!
!-----------------------------------------------------------------------------
SUBROUTINE rpa_epsilon (q, w, nmodes, epsil, eps_rpa)
!-----------------------------------------------------------------------------
!
! Compute the Lindhard dielectric function for the homogeneous electron gas
!
USE kinds, ONLY : dp
USE cell_base, ONLY : at, bg, omega, alat
USE constants_epw, ONLY : pi, twopi, ha2ev, cone
USE epwcom, ONLY : meff, fermi_diff, nel, smear_rpa
USE io_global, ONLY : stdout
!
implicit none
!
INTEGER, INTENT (in) :: nmodes
!! Number of phonon modes
!
REAL (kind=DP), INTENT (inout) :: q(3)
!! q vector (in crystal coordinates
REAL (kind=DP), INTENT (inout) :: w(nmodes)
!! phonon frequencies associated with q
REAL (kind=DP), INTENT (in) :: epsil(3,3)
!! dielectric constant tensor
!
COMPLEX (kind=DP), INTENT (out) :: eps_rpa(nmodes)
!! electronic screening
!
! Working variable
INTEGER :: im
!! Mode counter
!
REAL (kind=DP) :: n
!! Electron density in atomic units
REAL (kind=DP) :: rs
!! Prefactor for the dielectric screening
REAL (kind=DP) :: EF
!! Fermi-level in eV
REAL (kind=DP) :: kF
!! Fermi wavevector
REAL (kind=DP) :: pref
!! Prefactor for the dielectric function
REAL (kind=DP) :: eta
!! Broadening for the dielectric function
REAL (kind=DP) :: q2
!! q-point square
REAL (kind=DP) :: qm
!! Internal units for Hedin's formula
REAL (kind=DP) :: eps_ave
!! Average dielectric function (semiconductor/insulator)
!
COMPLEX (kind=DP) :: u
!! Complex frequency argument
!
LOGICAL, SAVE :: first_call=.true.
!! Logical for first_call the routine