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

Correct some circular dep. in EPW

This commit is contained in:
Samuel Ponce 2018-08-29 17:28:35 +01:00
parent e2f6df6f6d
commit 72a8f3c7d6
7 changed files with 539 additions and 849 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
EPW/src/Makefile
View File

@ -22,12 +22,12 @@ EPWOBJS = \
epwcom.o \
constants_epw.o \
io_epw.o \
io_eliashberg.o \
elph2.o \
eliashbergcom.o \
superconductivity_iso.o \
superconductivity_aniso.o \
superconductivity.o \
superconductivity_aniso.o \
superconductivity_iso.o \
io_eliashberg.o \
transportcom.o \
transport.o \
transport_iter.o \

4
EPW/src/eliashberg.f90
View File

@ -15,12 +15,12 @@
USE io_global, ONLY : stdout
USE epwcom, ONLY : liso, fila2f, gap_edge, lreal, limag, laniso
USE eliashbergcom, ONLY : gap0
USE superconductivity, ONLY : eliashberg_init, evaluate_a2f_lambda, estimate_tc_gap, &
USE superconductivity, ONLY : eliashberg_init, estimate_tc_gap, &
deallocate_eliashberg
USE io_eliashberg, ONLY : read_a2f, read_frequencies, read_eigenvalues, read_ephmat, &
read_kqmap
USE superconductivity_iso, ONLY : eliashberg_iso_iaxis, eliashberg_iso_raxis
USE superconductivity_aniso, ONLY : eliashberg_aniso_iaxis
USE superconductivity_aniso, ONLY : eliashberg_aniso_iaxis, evaluate_a2f_lambda
!
IMPLICIT NONE
!

104
EPW/src/io_eliashberg.f90
View File

@ -364,6 +364,110 @@
END SUBROUTINE eliashberg_write_raxis
!
!-----------------------------------------------------------------------
SUBROUTINE eliashberg_write_cont_raxis( itemp, cname )
!-----------------------------------------------------------------------
!
!
! This routine writes to files results from the solutions of the Eliashberg
! equations on the real-axis
!
USE kinds, ONLY : DP
USE io_epw, ONLY : iufilgap
USE io_files, ONLY : prefix
USE control_flags, ONLY : iverbosity
USE epwcom, ONLY : nqstep, fsthick, laniso, liso
USE eliashbergcom, ONLY : nsw, estemp, ws, gap, Agap, &
Delta, Znorm, ADelta, AZnorm, &
nkfs, nbndfs, ef0, ekfs
USE constants_epw, ONLY : kelvin2eV
!
IMPLICIT NONE
!
INTEGER :: iw, itemp, ik, ibnd
REAL(DP) :: temp
LOGICAL :: lgap
CHARACTER(len=256) :: name1, cname
!
temp = estemp(itemp) / kelvin2eV
!
IF ( laniso ) THEN
IF ( iverbosity .eq. 2 ) THEN
IF ( temp .lt. 10.d0 ) THEN
WRITE(name1,'(a,a1,a4,a9,f4.2)') TRIM(prefix), '.', cname, '_aniso_00', temp
ELSEIF ( temp .ge. 10.d0 .AND. temp .lt. 100.d0 ) THEN
WRITE(name1,'(a,a1,a4,a8,f5.2)') TRIM(prefix), '.', cname, '_aniso_0', temp
ELSEIF ( temp .ge. 100.d0 ) THEN
WRITE(name1,'(a,a1,a4,a7,f6.2)') TRIM(prefix), '.', cname, '_aniso_', temp
ENDIF
OPEN(iufilgap, file=name1, form='formatted')
WRITE(iufilgap,'(6a20)') '# w [eV]', 'Enk-Ef [eV]', 'Re[Znorm(w)]', 'Im[Znorm(w)]',&
'Re[Delta(w)] [eV]', 'Im[Delta(w)] [eV]'
ENDIF
!
DO ik = 1, nkfs
DO ibnd = 1, nbndfs
IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
lgap = .true.
! DO iw = 1, nsw
DO iw = 1, nsw-1 ! FG: this change is to prevent segfault in ws(iw+1) and ADelta(*,*,iw+1)
IF ( lgap .AND. iw .lt. nqstep .AND. real(ADelta(ibnd,ik,iw)) .gt. 0.d0 &
.AND. real(ADelta(ibnd,ik,iw+1)) .gt. 0.d0 &
.AND. ( ws(iw) - real(ADelta(ibnd,ik,iw)) )*( ws(iw+1) - real(ADelta(ibnd,ik,iw+1)) ) .lt. 0.d0 ) THEN
Agap(ibnd,ik,itemp) = ( ( real(ADelta(ibnd,ik,iw)) - ws(iw) ) * ws(iw+1) &
- ( real(ADelta(ibnd,ik,iw+1)) - ws(iw+1) ) * ws(iw) ) &
/ ( ( real(ADelta(ibnd,ik,iw)) - ws(iw) ) - ( real(ADelta(ibnd,ik,iw+1)) - ws(iw+1) ) )
lgap = .false.
ENDIF
IF ( iverbosity .eq. 2 ) THEN
WRITE(iufilgap,'(6ES20.10)') ws(iw), ekfs(ibnd,ik)-ef0, &
real(AZnorm(ibnd,ik,iw)), aimag(AZnorm(ibnd,ik,iw)), &
real(ADelta(ibnd,ik,iw)), aimag(ADelta(ibnd,ik,iw))
ENDIF
ENDDO ! iw
IF ( lgap ) &
Agap(ibnd,ik,itemp) = real(ADelta(ibnd,ik,1))
ENDIF
ENDDO ! ibnd
ENDDO ! ik
IF ( iverbosity .eq. 2 ) CLOSE(iufilgap)
!
CALL gap_distribution_FS ( itemp, cname )
!
ENDIF
!
! isotropic case
! SP: Only write isotropic for laniso if user really wants that
IF ( ( laniso .AND. iverbosity .eq. 2 ) .OR. liso ) THEN
IF ( temp .lt. 10.d0 ) THEN
WRITE(name1,'(a,a1,a4,a7,f4.2)') TRIM(prefix), '.', cname, '_iso_00', temp
ELSEIF ( temp .ge. 10.d0 .AND. temp .lt. 100.d0 ) THEN
WRITE(name1,'(a,a1,a4,a6,f5.2)') TRIM(prefix), '.', cname, '_iso_0', temp
ELSEIF ( temp .ge. 100.d0 ) THEN
WRITE(name1,'(a,a1,a4,a5,f6.2)') TRIM(prefix), '.', cname, '_iso_', temp
ENDIF
OPEN(iufilgap, file=name1, form='formatted')
WRITE(iufilgap,'(5a20)') 'w [eV]', 'Re[Znorm(w)]', 'Im[Znorm(w)]', 'Re[Delta(w)] [eV]', 'Im[Delta(w)] [eV]'
lgap = .true.
! DO iw = 1, nsw
DO iw = 1, nsw-1 ! this change is to prevent segfault in Delta(iw+1) and ws(iw+1)
IF ( lgap .AND. iw .lt. nqstep .AND. real(Delta(iw)) .gt. 0.d0 .AND. real(Delta(iw+1)) .gt. 0.d0 .AND. &
( ws(iw) - real(Delta(iw)) )*( ws(iw+1) - real(Delta(iw+1)) ) .lt. 0.d0 ) THEN
gap(itemp) = ( ( real(Delta(iw)) - ws(iw) ) * ws(iw+1) - ( real(Delta(iw+1)) - ws(iw+1) ) * ws(iw) ) &
/ ( ( real(Delta(iw)) - ws(iw) ) - ( real(Delta(iw+1)) - ws(iw+1) ) )
lgap = .false.
ENDIF
WRITE(iufilgap,'(5ES20.10)') ws(iw), real(Znorm(iw)), aimag(Znorm(iw)), &
real(Delta(iw)), aimag(Delta(iw))
ENDDO ! iw
CLOSE(iufilgap)
IF ( lgap ) &
gap(itemp) = real(Delta(1))
ENDIF
!
RETURN
!
END SUBROUTINE eliashberg_write_cont_raxis
!-----------------------------------------------------------------------
SUBROUTINE read_a2f
!-----------------------------------------------------------------------
!!

186
EPW/src/make.depend
View File

@ -7,75 +7,27 @@ allocate_epwq.o : epwcom.o
allocate_epwq.o : transportcom.o
bcast_epw_input.o : elph2.o
bcast_epw_input.o : epwcom.o
bloch2wan.o : constants_epw.o
bloch2wan.o : elph2.o
bloch2wan.o : epwcom.o
bloch2wan.o : io_epw.o
broyden.o : eliashbergcom.o
broyden.o : epwcom.o
close_epw.o : epwcom.o
close_epw.o : io_epw.o
create_mesh.o : eliashbergcom.o
create_mesh.o : elph2.o
create_mesh.o : epwcom.o
create_mesh.o : io_epw.o
createkmap.o : elph2.o
createkmap.o : epwcom.o
createkmap.o : kfold.o
deallocate_eliashberg.o : eliashbergcom.o
deallocate_eliashberg.o : elph2.o
deallocate_eliashberg.o : epwcom.o
deallocate_epw.o : elph2.o
deallocate_epw.o : epwcom.o
deallocate_epw.o : transportcom.o
dmebloch2wan.o : constants_epw.o
dmebloch2wan.o : elph2.o
dmewan2bloch.o : constants_epw.o
dmewan2bloch.o : elph2.o
dmewan2bloch.o : epwcom.o
dvqpsi_us3.o : elph2.o
dvqpsi_us_only3.o : elph2.o
dynbloch2wan.o : constants_epw.o
dynbloch2wan.o : elph2.o
dynbloch2wan.o : epwcom.o
dynwan2bloch.o : constants_epw.o
dynwan2bloch.o : elph2.o
dynwan2bloch.o : epwcom.o
eliashberg.o : eliashbergcom.o
eliashberg.o : epwcom.o
eliashberg_aniso_cont_raxis.o : constants_epw.o
eliashberg_aniso_cont_raxis.o : eliashbergcom.o
eliashberg_aniso_cont_raxis.o : elph2.o
eliashberg_aniso_cont_raxis.o : epwcom.o
eliashberg_aniso_iaxis.o : constants_epw.o
eliashberg_aniso_iaxis.o : eliashbergcom.o
eliashberg_aniso_iaxis.o : elph2.o
eliashberg_aniso_iaxis.o : epwcom.o
eliashberg_aniso_iaxis.o : io_epw.o
eliashberg_iso_cont_raxis.o : constants_epw.o
eliashberg_iso_cont_raxis.o : eliashbergcom.o
eliashberg_iso_cont_raxis.o : epwcom.o
eliashberg_iso_iaxis.o : constants_epw.o
eliashberg_iso_iaxis.o : eliashbergcom.o
eliashberg_iso_iaxis.o : epwcom.o
eliashberg_iso_raxis.o : constants_epw.o
eliashberg_iso_raxis.o : eliashbergcom.o
eliashberg_iso_raxis.o : epwcom.o
eliashberg_iso_raxis.o : io_epw.o
eliashberg_pp.o : constants_epw.o
eliashberg_pp.o : eliashbergcom.o
eliashberg_pp.o : epwcom.o
eliashberg_pp.o : io_epw.o
eliashberg_readfiles.o : constants_epw.o
eliashberg_readfiles.o : eliashbergcom.o
eliashberg_readfiles.o : elph2.o
eliashberg_readfiles.o : epwcom.o
eliashberg_readfiles.o : io_epw.o
eliashberg_setup.o : constants_epw.o
eliashberg_setup.o : eliashbergcom.o
eliashberg_setup.o : elph2.o
eliashberg_setup.o : epwcom.o
eliashberg_setup.o : io_epw.o
eliashberg_write.o : constants_epw.o
eliashberg_write.o : eliashbergcom.o
eliashberg_write.o : epwcom.o
eliashberg_write.o : io_epw.o
eliashberg.o : io_eliashberg.o
eliashberg.o : superconductivity.o
eliashberg.o : superconductivity_aniso.o
eliashberg.o : superconductivity_iso.o
elphel2_shuffle.o : constants_epw.o
elphel2_shuffle.o : elph2.o
elphon_shuffle.o : constants_epw.o
@ -84,42 +36,30 @@ elphon_shuffle_wrap.o : constants_epw.o
elphon_shuffle_wrap.o : elph2.o
elphon_shuffle_wrap.o : epwcom.o
elphon_shuffle_wrap.o : io_epw.o
ephbloch2wane.o : constants_epw.o
ephbloch2wane.o : io_epw.o
ephbloch2wanp.o : constants_epw.o
ephbloch2wanp.o : elph2.o
ephbloch2wanp.o : io_epw.o
ephwan2bloch.o : constants_epw.o
ephwan2bloch_mem.o : constants_epw.o
ephwan2blochp.o : constants_epw.o
ephwan2blochp.o : elph2.o
ephwan2blochp.o : epwcom.o
ephwan2blochp.o : io_epw.o
ephwan2blochp_mem.o : constants_epw.o
ephwan2blochp_mem.o : elph2.o
ephwan2blochp_mem.o : epwcom.o
ephwan2blochp_mem.o : io_epw.o
ephwann_shuffle.o : bloch2wan.o
ephwann_shuffle.o : constants_epw.o
ephwann_shuffle.o : elph2.o
ephwann_shuffle.o : epwcom.o
ephwann_shuffle.o : io_eliashberg.o
ephwann_shuffle.o : io_epw.o
ephwann_shuffle.o : transport.o
ephwann_shuffle.o : transport_iter.o
ephwann_shuffle.o : transportcom.o
ephwann_shuffle.o : wan2bloch.o
ephwann_shuffle.o : bloch2wan.o
ephwann_shuffle_mem-save.o : constants_epw.o
ephwann_shuffle_mem-save.o : elph2.o
ephwann_shuffle_mem-save.o : epwcom.o
ephwann_shuffle_mem-save.o : io_epw.o
ephwann_shuffle.o : wigner.o
ephwann_shuffle_mem.o : bloch2wan.o
ephwann_shuffle_mem.o : constants_epw.o
ephwann_shuffle_mem.o : elph2.o
ephwann_shuffle_mem.o : epwcom.o
ephwann_shuffle_mem.o : io_eliashberg.o
ephwann_shuffle_mem.o : io_epw.o
ephwann_shuffle_mem.o : transport.o
ephwann_shuffle_mem.o : transport_iter.o
ephwann_shuffle_mem.o : transportcom.o
ephwann_shuffle_mem.o : wan2bloch.o
ephwann_shuffle_mem.o : bloch2wan.o
ephwann_shuffle_mem.o : wigner.o
epw.o : elph2.o
epw.o : epwcom.o
epw_init.o : constants_epw.o
epw_init.o : elph2.o
epw_readin.o : constants_epw.o
epw_readin.o : elph2.o
@ -130,33 +70,29 @@ epw_setup.o : transportcom.o
epw_summary.o : epwcom.o
fermiwindow.o : elph2.o
fermiwindow.o : epwcom.o
gen_freqgrid.o : eliashbergcom.o
gen_freqgrid.o : epwcom.o
gmap_sym.o : constants_epw.o
hambloch2wan.o : constants_epw.o
hamwan2bloch.o : constants_epw.o
indabs.o : constants_epw.o
indabs.o : elph2.o
indabs.o : epwcom.o
indabs.o : io_epw.o
indabs.o : transportcom.o
io_eliashberg.o : constants_epw.o
io_eliashberg.o : eliashbergcom.o
io_eliashberg.o : elph2.o
io_eliashberg.o : epwcom.o
io_eliashberg.o : io_epw.o
io_eliashberg.o : superconductivity.o
io_scattering.o : constants_epw.o
io_scattering.o : elph2.o
io_scattering.o : epwcom.o
io_scattering.o : io_epw.o
io_scattering.o : transportcom.o
iterativebte.o : constants_epw.o
iterativebte.o : elph2.o
iterativebte.o : epwcom.o
iterativebte.o : io_epw.o
iterativebte.o : transportcom.o
kernels_aniso_iaxis.o : constants_epw.o
kernels_aniso_iaxis.o : eliashbergcom.o
kernels_aniso_iaxis.o : elph2.o
kernels_aniso_iaxis.o : epwcom.o
kernels_iso_iaxis.o : constants_epw.o
kernels_iso_iaxis.o : eliashbergcom.o
kernels_iso_iaxis.o : epwcom.o
kernels_raxis.o : constants_epw.o
kernels_raxis.o : eliashbergcom.o
kernels_raxis.o : epwcom.o
ktokpmq.o : elph2.o
ktokpmq.o : epwcom.o
load_rebal.o : constants_epw.o
load_rebal.o : elph2.o
load_rebal.o : epwcom.o
load_rebal.o : transportcom.o
loadkmesh.o : elph2.o
loadkmesh.o : epwcom.o
loadkmesh.o : io_epw.o
@ -190,6 +126,7 @@ readmat_shuffle2.o : elph2.o
readmat_shuffle2.o : epwcom.o
readmat_shuffle2.o : io_dyn_mat2.o
readmat_shuffle2.o : io_epw.o
readmat_shuffle2.o : wan2bloch.o
refold.o : io_epw.o
refold.o : kfold.o
rgd_blk_epw_fine_mem.o : constants_epw.o
@ -201,11 +138,6 @@ rotate_eigenm.o : elph2.o
rotate_epmat.o : constants_epw.o
rotate_epmat.o : elph2.o
rotate_epmat.o : epwcom.o
scattering_rate.o : constants_epw.o
scattering_rate.o : elph2.o
scattering_rate.o : epwcom.o
scattering_rate.o : io_epw.o
scattering_rate.o : transportcom.o
selfen_elec.o : constants_epw.o
selfen_elec.o : elph2.o
selfen_elec.o : epwcom.o
@ -214,7 +146,6 @@ selfen_elec.o : transportcom.o
selfen_phon.o : constants_epw.o
selfen_phon.o : elph2.o
selfen_phon.o : epwcom.o
selfen_phon.o : io_epw.o
selfen_pl.o : constants_epw.o
selfen_pl.o : elph2.o
selfen_pl.o : epwcom.o
@ -237,23 +168,42 @@ spectral_func_pl.o : constants_epw.o
spectral_func_pl.o : elph2.o
spectral_func_pl.o : epwcom.o
spectral_func_pl.o : io_epw.o
superconductivity.o : constants_epw.o
superconductivity.o : eliashbergcom.o
superconductivity.o : elph2.o
superconductivity.o : epwcom.o
superconductivity.o : io_epw.o
superconductivity_aniso.o : constants_epw.o
superconductivity_aniso.o : eliashbergcom.o
superconductivity_aniso.o : elph2.o
superconductivity_aniso.o : epwcom.o
superconductivity_aniso.o : io_eliashberg.o
superconductivity_aniso.o : io_epw.o
superconductivity_aniso.o : superconductivity.o
superconductivity_aniso.o : superconductivity_iso.o
superconductivity_iso.o : constants_epw.o
superconductivity_iso.o : eliashbergcom.o
superconductivity_iso.o : epwcom.o
superconductivity_iso.o : io_eliashberg.o
superconductivity_iso.o : io_epw.o
superconductivity_iso.o : superconductivity.o
system_mem_usage.o : io_epw.o
transport_coeffs.o : constants_epw.o
transport_coeffs.o : elph2.o
transport_coeffs.o : epwcom.o
transport_coeffs.o : io_epw.o
transport_coeffs.o : transportcom.o
vmebloch2wan.o : constants_epw.o
vmebloch2wan.o : elph2.o
vmewan2bloch.o : constants_epw.o
vmewan2bloch.o : elph2.o
vmewan2bloch.o : epwcom.o
transport.o : constants_epw.o
transport.o : elph2.o
transport.o : epwcom.o
transport.o : io_epw.o
transport.o : transportcom.o
transport_iter.o : constants_epw.o
transport_iter.o : elph2.o
transport_iter.o : epwcom.o
transport_iter.o : io_eliashberg.o
transport_iter.o : transportcom.o
wan2bloch.o : constants_epw.o
wan2bloch.o : elph2.o
wan2bloch.o : epwcom.o
wan2bloch.o : io_epw.o
wannierize.o : constants_epw.o
wannierize.o : epwcom.o
wannierize.o : io_epw.o
wannierize.o : wannierEPW.o
write_ephmat.o : constants_epw.o
write_ephmat.o : eliashbergcom.o
write_ephmat.o : elph2.o
write_ephmat.o : epwcom.o
write_ephmat.o : io_epw.o
wigner.o : constants_epw.o

725
EPW/src/superconductivity.f90
View File

@ -184,363 +184,6 @@
!
END SUBROUTINE eliashberg_init
!
!-----------------------------------------------------------------------
SUBROUTINE evaluate_a2f_lambda
!-----------------------------------------------------------------------
!
! computes the isotropic spectral function a2F(w), total lambda, and
! distribution of lambda
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE io_epw, ONLY : iua2ffil, iudosfil, iufillambda, iufillambdaFS
USE io_files, ONLY : prefix
USE phcom, ONLY : nmodes
USE cell_base, ONLY : bg
USE control_flags, ONLY : iverbosity
USE elph2, ONLY : nqtotf, wqf, wf
USE epwcom, ONLY : fsthick, eps_acustic, nqstep, degaussq, delta_qsmear, nqsmear, &
degaussw, nkf1, nkf2, nkf3
USE eliashbergcom, ONLY : nkfs, nbndfs, g2, ixkqf, ixqfs, nqfs, w0g, ekfs, ef0, dosef, wsph, &
wkfs, dwsph, a2f_iso, ixkff
USE constants_epw, ONLY : ryd2ev
USE io_global, ONLY : ionode_id
USE mp_global, ONLY : inter_pool_comm, my_pool_id, npool
USE mp_world, ONLY : mpime
USE mp, ONLY : mp_bcast, mp_barrier, mp_sum
USE superconductivity_aniso, ONLY : lambdar_aniso_ver1
!
IMPLICIT NONE
!
INTEGER :: ik, iq, iq0, iwph, ibnd, jbnd, imode, lower_bnd, upper_bnd, &
ismear, ibin, nbin, nbink, i, j, k
REAL(DP) :: weight, weightq, l_sum, lambda_eph, lambda_max(npool), &
sigma, dbin, dbink, x1, x2, x3
REAL(DP), ALLOCATABLE :: a2f(:,:), phdos(:,:), l_a2f(:), lambda_k(:,:), &
lambda_k_bin(:), lambda_pairs(:), a2f_modeproj(:,:), phdos_modeproj(:,:)
REAL(DP), EXTERNAL :: w0gauss
CHARACTER (len=256) :: name1
!
! This is only a quick fix since the subroutine was written for parallel execution - FG June 2014
#if ! defined(__MPI)
npool = 1
my_pool_id = 0
#endif
!
! degaussq is read from the input file in meV and converted to Ryd in epw_readin.f90
! go from Ryd to eV
degaussq = degaussq * ryd2ev
delta_qsmear = delta_qsmear * ryd2ev
!
CALL fkbounds( nkfs, lower_bnd, upper_bnd )
!
IF ( .not. ALLOCATED(a2f_iso) ) ALLOCATE( a2f_iso(nqstep) )
IF ( .not. ALLOCATED(a2f) ) ALLOCATE( a2f(nqstep,nqsmear) )
IF ( .not. ALLOCATED(a2f_modeproj) ) ALLOCATE( a2f_modeproj(nmodes,nqstep) )
a2f_iso(:) = 0.d0
a2f(:,:) = 0.d0
a2f_modeproj(:,:) = 0.d0
!
! RM - the 0 index in k is required when printing out values of lambda_k
! When the k-point is outside the Fermi shell, ixkff(ik)=0
IF ( .not. ALLOCATED(lambda_k) ) ALLOCATE(lambda_k(0:nkfs,nbndfs))
lambda_k(:,:) = 0.d0
!
l_sum = 0.d0
lambda_max(:) = 0.d0
DO ismear = 1, nqsmear
sigma = degaussq + (ismear-1) * delta_qsmear
DO ik = lower_bnd, upper_bnd
DO ibnd = 1, nbndfs
IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
DO iq = 1, nqfs(ik)
! iq0 - index of q-point on the full q-mesh
iq0 = ixqfs(ik,iq)
DO jbnd = 1, nbndfs
IF ( abs( ekfs(jbnd,ixkqf(ik,iq0)) - ef0 ) .lt. fsthick ) THEN
weight = wkfs(ik) * wqf(iq) * w0g(ibnd,ik) * w0g(jbnd,ixkqf(ik,iq0))
lambda_eph = 0.d0
DO imode = 1, nmodes
IF ( wf(imode,iq0) .gt. eps_acustic ) THEN
IF ( ismear .eq. 1 ) THEN
lambda_eph = lambda_eph + g2(ik,iq,ibnd,jbnd,imode) / wf(imode,iq0)
ENDIF
DO iwph = 1, nqstep
weightq = w0gauss( ( wsph(iwph) - wf(imode,iq0) ) / sigma, 0 ) / sigma
a2f(iwph,ismear) = a2f(iwph,ismear) + weight * weightq * g2(ik,iq,ibnd,jbnd,imode)
IF ( ismear .eq. 1 ) THEN
a2f_modeproj(imode,iwph) = a2f_modeproj(imode,iwph) +&
weight * weightq * g2(ik,iq,ibnd,jbnd,imode)
ENDIF
ENDDO ! iwph
ENDIF ! wf
ENDDO ! imode
IF ( ismear .eq. 1 .AND. lambda_eph .gt. 0.d0 ) THEN
l_sum = l_sum + weight * lambda_eph
weight = wqf(iq) * w0g(jbnd,ixkqf(ik,iq0))
lambda_k(ik,ibnd) = lambda_k(ik,ibnd) + weight * lambda_eph
IF ( lambda_eph .gt. lambda_max(my_pool_id+1) ) THEN
lambda_max(my_pool_id+1) = lambda_eph
ENDIF
ENDIF
ENDIF ! ekq
ENDDO ! jbnd
ENDDO ! iq
ENDIF ! ekk
ENDDO ! ibnd
ENDDO ! ik
ENDDO ! ismear
!
a2f(:,:) = 0.5d0 * a2f(:,:) / dosef
a2f_modeproj(:,:) = 0.5d0 * a2f_modeproj(:,:) / dosef
l_sum = l_sum / dosef
lambda_k(:,:) = 2.d0 * lambda_k(:,:)
lambda_max(:) = 2.d0 * dosef * lambda_max(:)
!
! collect contributions from all pools (sum over k-points)
CALL mp_sum( l_sum, inter_pool_comm )
CALL mp_sum( a2f, inter_pool_comm )
CALL mp_sum( a2f_modeproj, inter_pool_comm )
CALL mp_sum( lambda_max, inter_pool_comm )
CALL mp_sum( lambda_k, inter_pool_comm )
CALL mp_barrier(inter_pool_comm)
!
IF ( mpime .eq. ionode_id ) THEN
!
OPEN( unit = iua2ffil, file = TRIM(prefix)//".a2f", form = 'formatted')
OPEN( unit = iudosfil, file = TRIM(prefix)//".phdos", form = 'formatted')
!
IF ( .not. ALLOCATED(phdos) ) ALLOCATE( phdos(nqstep,nqsmear) )
IF ( .not. ALLOCATED(phdos_modeproj) ) ALLOCATE( phdos_modeproj(nmodes,nqstep) )
phdos(:,:) = 0.d0
phdos_modeproj(:,:) = 0.d0
!
DO ismear = 1, nqsmear
sigma = degaussq + (ismear-1) * delta_qsmear
DO iq = 1, nqtotf
DO imode = 1, nmodes
IF ( wf(imode,iq) .gt. eps_acustic ) THEN
DO iwph = 1, nqstep
weightq = w0gauss( ( wsph(iwph) - wf(imode,iq)) / sigma, 0 ) / sigma
phdos(iwph,ismear) = phdos(iwph,ismear) + wqf(iq) * weightq
IF ( ismear .eq. 1 ) THEN
phdos_modeproj(imode,iwph) = phdos_modeproj(imode,iwph) + wqf(iq) * weightq
ENDIF
ENDDO ! iwph
ENDIF ! wf
ENDDO ! imode
ENDDO ! iq
ENDDO ! ismear
!
IF ( .not. ALLOCATED(l_a2f) ) ALLOCATE( l_a2f(nqsmear) )
l_a2f(:) = 0.d0
!
DO ismear = 1, nqsmear
DO iwph = 1, nqstep
l_a2f(ismear) = l_a2f(ismear) + a2f(iwph,ismear) / wsph(iwph)
! wsph in meV (from eV) and phdos in states/meV (from states/eV)
IF (ismear .eq. nqsmear) WRITE (iua2ffil,'(f12.7,15f12.7)') wsph(iwph)*1000.d0, a2f(iwph,:)
IF (ismear .eq. nqsmear) WRITE (iudosfil,'(f12.7,15f15.7)') wsph(iwph)*1000.d0, phdos(iwph,:)/1000.d0
ENDDO
l_a2f(ismear) = 2.d0 * l_a2f(ismear) * dwsph
ENDDO
!
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)*1000.d0,ismear=1,nqsmear )
WRITE(iua2ffil,'(" Electron smearing (eV)", f12.7)') degaussw
WRITE(iua2ffil,'(" Fermi window (eV)", f12.7)') fsthick
WRITE(iua2ffil,'(" Summed el-ph coupling ", f12.7)') l_sum
CLOSE(iua2ffil)
CLOSE(iudosfil)
!
a2f_iso(:) = a2f(:,1)
OPEN( unit = iua2ffil, file = TRIM(prefix)//".a2f_iso", form = 'formatted')
OPEN( unit = iudosfil, file = TRIM(prefix)//".phdos_proj", form = 'formatted')
DO iwph = 1, nqstep
! wsph in meV (from eV) and phdos in states/meV (from states/eV)
WRITE(iua2ffil,'(f12.7,100f12.7)') wsph(iwph)*1000.d0, a2f_iso(iwph), a2f_modeproj(:,iwph)
WRITE(iudosfil,'(f12.7,100f15.7)') wsph(iwph)*1000.d0, phdos(iwph,1)/1000.d0, phdos_modeproj(:,iwph)/1000.d0
ENDDO
WRITE(iua2ffil,'(a,f18.7,a,f18.7)') 'lambda_int = ', l_a2f(1), ' lambda_sum = ',l_sum
CLOSE(iua2ffil)
CLOSE(iudosfil)
!
IF ( ALLOCATED(phdos) ) DEALLOCATE( phdos )
IF ( ALLOCATED(phdos_modeproj) ) DEALLOCATE( phdos_modeproj )
IF ( ALLOCATED(l_a2f) ) DEALLOCATE( l_a2f )
!
ENDIF
!
CALL mp_bcast( a2f_iso, ionode_id, inter_pool_comm )
CALL mp_barrier(inter_pool_comm)
!
IF ( ALLOCATED(a2f) ) DEALLOCATE( a2f )
IF ( ALLOCATED(a2f_modeproj) ) DEALLOCATE( a2f_modeproj )
!
nbink = int( 1.25d0 * maxval(lambda_k(:,:)) / 0.005d0 )
dbink = 1.25d0 * maxval(lambda_k(:,:)) / dble(nbink)
IF ( .not. ALLOCATED(lambda_k_bin) ) ALLOCATE ( lambda_k_bin(nbink) )
lambda_k_bin(:) = 0.d0
!
!SP : Should be initialized
nbin = 0
dbin = 0.0_DP
!
IF ( iverbosity == 2 ) THEN
nbin = int( 1.25d0 * maxval(lambda_max(:)) / 0.005d0 )
dbin = 1.25d0 * maxval(lambda_max(:)) / dble(nbin)
IF ( .not. ALLOCATED(lambda_pairs) ) ALLOCATE ( lambda_pairs(nbin) )
lambda_pairs(:) = 0.d0
ENDIF
!
WRITE(stdout,'(5x,a13,f21.7,a18,f21.7)') 'lambda_max = ', maxval(lambda_max(:)), ' lambda_k_max = ', maxval(lambda_k(:,:))
WRITE(stdout,'(a)') ' '
!
lambda_k(:,:) = 0.d0
DO ik = lower_bnd, upper_bnd
DO ibnd = 1, nbndfs
IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
DO iq = 1, nqfs(ik)
! iq0 - index of q-point on the full q-mesh
iq0 = ixqfs(ik,iq)
DO jbnd = 1, nbndfs
IF ( abs( ekfs(jbnd,ixkqf(ik,iq0)) - ef0 ) .lt. fsthick ) THEN
weight = wqf(iq) * w0g(jbnd,ixkqf(ik,iq0)) / dosef
CALL lambdar_aniso_ver1( ik, iq, ibnd, jbnd, 0.d0, lambda_eph )
lambda_k(ik,ibnd) = lambda_k(ik,ibnd) + weight * lambda_eph
IF ( iverbosity == 2 ) THEN
DO ibin = 1, nbin
sigma = 1.d0 * dbin
weight = w0gauss( ( lambda_eph - dble(ibin) * dbin ) / sigma, 0 ) / sigma
lambda_pairs(ibin) = lambda_pairs(ibin) + weight
ENDDO
ENDIF
ENDIF
ENDDO ! jbnd
ENDDO ! iq
DO ibin = 1, nbink
sigma = 1.d0 * dbink
weight = w0gauss( ( lambda_k(ik,ibnd) - dble(ibin) * dbink ) / sigma, 0 ) / sigma
lambda_k_bin(ibin) = lambda_k_bin(ibin) + weight
ENDDO
ENDIF
ENDDO ! ibnd
ENDDO ! ik
!
! collect contributions from all pools
CALL mp_sum( lambda_k, inter_pool_comm )
IF ( iverbosity .eq. 2 ) THEN
CALL mp_sum( lambda_pairs, inter_pool_comm )
ENDIF
CALL mp_sum( lambda_k_bin, inter_pool_comm )
CALL mp_barrier(inter_pool_comm)
!
IF ( mpime .eq. ionode_id ) THEN
!
! SP: Produced if user really wants it
IF ( iverbosity .eq. 2 ) THEN
OPEN(unit = iufillambda, file = TRIM(prefix)//".lambda_aniso", form = 'formatted')
WRITE(iufillambda,'(2a12,2a7)') '# enk-e0[eV]',' lambda_nk','# kpt','# band'
DO ik = 1, nkfs
DO ibnd = 1, nbndfs
IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
WRITE(iufillambda,'(2f12.7,2i7)') ekfs(ibnd,ik) - ef0, lambda_k(ik,ibnd), ik, ibnd
ENDIF
ENDDO
ENDDO
CLOSE(iufillambda)
ENDIF
!
OPEN(unit = iufillambda, file = TRIM(prefix)//".lambda_k_pairs", form = 'formatted')
WRITE(iufillambda,'(a12,a30)') '# lambda_nk',' \rho(lambda_nk) scaled to 1'
DO ibin = 1, nbink
WRITE(iufillambda,'(2f21.7)') dbink*dble(ibin), lambda_k_bin(ibin)/maxval(lambda_k_bin(:))
ENDDO
CLOSE(iufillambda)
!
! SP: Produced if user really wants it
IF ( iverbosity == 2 ) THEN
OPEN( unit = iufillambda, file = TRIM(prefix)//".lambda_pairs", form = 'formatted')
WRITE(iufillambda,'(a12,a30)') "# lambda_nk,n'k'", " \rho(lambda_nk,n'k') scaled to 1"
DO ibin = 1, nbin
WRITE(iufillambda,'(2f21.7)') dbin*dble(ibin), lambda_pairs(ibin)/maxval(lambda_pairs(:))
ENDDO
CLOSE(iufillambda)
ENDIF
!
! SP & RM: .cube file for VESTA plotting (only if iverbosity = 2)
!
! RM - If the k-point is outside the Fermi shell,
! ixkff(ik)=0 and lambda_k(0,ibnd) = 0.0
!
IF ( iverbosity .eq. 2 ) THEN
!
DO ibnd = 1, nbndfs
!
IF ( ibnd < 10 ) THEN
WRITE(name1,'(a,a8,i1,a5)') TRIM(prefix),'.lambda_', ibnd, '.cube'
ELSEIF ( ibnd < 100 ) THEN
WRITE(name1,'(a,a8,i2,a5)') TRIM(prefix),'.lambda_', ibnd, '.cube'
ELSEIF( ibnd < 1000 ) THEN
WRITE(name1,'(a,a8,i3,a5)') TRIM(prefix),'.lambda_', ibnd, '.cube'
ELSE
CALL errore( 'eliashberg_setup', 'Too many bands ',1)
ENDIF
!
OPEN(iufillambdaFS, file=name1, form='formatted')
WRITE(iufillambdaFS,*) 'Cubfile created from EPW calculation'
WRITE(iufillambdaFS,*) 'lambda'
WRITE(iufillambdaFS,'(i5,3f12.6)') 1, 0.0d0, 0.0d0, 0.0d0
WRITE(iufillambdaFS,'(i5,3f12.6)') nkf1, (bg(i,1)/DBLE(nkf1),i=1,3)
WRITE(iufillambdaFS,'(i5,3f12.6)') nkf2, (bg(i,2)/DBLE(nkf2),i=1,3)
WRITE(iufillambdaFS,'(i5,3f12.6)') nkf3, (bg(i,3)/DBLE(nkf3),i=1,3)
WRITE(iufillambdaFS,'(i5,4f12.6)') 1, 1.0d0, 0.0d0, 0.0d0, 0.0d0
WRITE(iufillambdaFS,'(6f12.6)') ( lambda_k(ixkff(ik),ibnd), ik=1,nkf1*nkf2*nkf3 )
CLOSE(iufillambdaFS)
ENDDO
!
ENDIF
!
! SP & RM : Write on file the lambda close to the Fermi surface along with
! Cartesian coordinate, band index, energy distance from Fermi level
! and lambda value.
!
OPEN(unit = iufillambdaFS, file = TRIM(prefix)//".lambda_FS", form='formatted')
WRITE(iufillambdaFS,'(a75)') '# k-point Band Enk-Ef [eV] lambda'
DO i = 1, nkf1
DO j = 1, nkf2
DO k = 1, nkf3
ik = k + (j-1)*nkf3 + (i-1)*nkf2*nkf3
IF ( ixkff(ik) .gt. 0 ) THEN
DO ibnd = 1, nbndfs
! SP: Here take a 0.2 eV interval around the FS.
IF ( abs( ekfs(ibnd,ixkff(ik)) - ef0 ) .lt. fsthick ) THEN
!IF ( abs( ekfs(ibnd,ixkff(ik)) - ef0 ) .lt. 0.2 ) THEN
x1 = bg(1,1)*(i-1)/nkf1+bg(1,2)*(j-1)/nkf2+bg(1,3)*(k-1)/nkf3
x2 = bg(2,1)*(i-1)/nkf1+bg(2,2)*(j-1)/nkf2+bg(2,3)*(k-1)/nkf3
x3 = bg(3,1)*(i-1)/nkf1+bg(3,2)*(j-1)/nkf2+bg(3,3)*(k-1)/nkf3
WRITE(iufillambdaFS,'(3f12.6,i8,f12.6,f24.15)') x1, x2, x3, ibnd, &
ekfs(ibnd,ixkff(ik))-ef0, lambda_k(ixkff(ik),ibnd)
ENDIF
ENDDO ! ibnd
ENDIF
ENDDO ! k
ENDDO ! j
ENDDO ! i
CLOSE(iufillambdaFS)
ENDIF
CALL mp_barrier(inter_pool_comm)
!
IF ( ALLOCATED(lambda_k) ) DEALLOCATE(lambda_k)
IF ( ALLOCATED(lambda_pairs) ) DEALLOCATE(lambda_pairs)
IF ( ALLOCATED(lambda_k_bin) ) DEALLOCATE(lambda_k_bin)
!
RETURN
!
END SUBROUTINE evaluate_a2f_lambda
!
!-----------------------------------------------------------------------
SUBROUTINE estimate_tc_gap
!-----------------------------------------------------------------------
@ -715,267 +358,6 @@
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE eliashberg_iso_raxis
!-----------------------------------------------------------------------
!!
!! This routine is the driver of the self-consistent cycle for the isotropic
!! Eliashberg equations on the real-axis.
!!
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE io_epw, ONLY : iufilgap
USE io_files, ONLY : prefix
USE epwcom, ONLY : nsiter, nstemp, broyden_beta, broyden_ndim
USE eliashbergcom, ONLY : nsw, Delta, Deltap, gap, estemp
USE constants_epw, ONLY : kelvin2eV, ci
USE mp, ONLY : mp_bcast, mp_barrier, mp_sum
!
IMPLICIT NONE
!
! Local variables
INTEGER :: itemp, iter
REAL(DP) :: tcpu, rdeltaout(nsw), rdeltain(nsw), cdeltaout(nsw), cdeltain(nsw)
REAL(DP), EXTERNAL :: get_clock
LOGICAL :: conv
CHARACTER (len=256) :: filgap
!
CALL start_clock( 'iso_raxis' )
!
WRITE(stdout,'(5x,a)') 'Solve isotropic Eliashberg equations on real-axis'
!
CALL gen_freqgrid_raxis
!
DO itemp = 1, nstemp ! loop over temperature
!
WRITE(stdout,'(a)') ' '
WRITE(stdout,'(5x,a,i3,a,f8.4,a,a,i3,a)') 'temp(', itemp, ') = ', estemp(itemp)/kelvin2eV, ' K '
WRITE(stdout,'(a)') ' '
iter = 1
conv = .false.
DO WHILE ( .not. conv .AND. iter .le. nsiter )
CALL integrate_eliashberg_iso_raxis( itemp, iter, conv )
rdeltain(:) = real(Deltap(:))
cdeltain(:) = aimag(Deltap(:))
rdeltaout(:) = real(Delta(:))
cdeltaout(:) = aimag(Delta(:))
CALL mix_broyden ( nsw, rdeltaout, rdeltain, broyden_beta, iter, broyden_ndim, conv )
CALL mix_broyden2( nsw, cdeltaout, cdeltain, broyden_beta, iter, broyden_ndim, conv )
Deltap(:) = rdeltain(:) + ci * cdeltain(:)
iter = iter + 1
ENDDO ! iter
WRITE(stdout,'(5x,a,i3,a,f8.4,a,a,i3,a,f10.6,a,a,f10.6,a)') &
'temp(', itemp, ') = ', estemp(itemp)/kelvin2eV, ' K ', &
' gap_edge(', itemp, ') = ', gap(itemp), ' eV ', &
' Re[Delta(1)] = ', real(Delta(1)), ' eV '
WRITE(stdout,'(a)') ' '
tcpu = get_clock( 'iso_raxis' )
WRITE( stdout,'(5x,a,i3,a,f8.1,a)') 'itemp = ', itemp, ' total cpu time :', tcpu, ' secs'
!
IF ( conv ) THEN
WRITE(stdout,'(a)') ' '
CALL print_clock( 'iso_raxis' )
WRITE(stdout,'(a)') ' '
ELSEIF ( .not. conv .AND. (iter-1) .eq. nsiter ) THEN
CALL deallocate_eliashberg
WRITE(stdout,'(a)') ' '
CALL stop_clock( 'iso_raxis' )
CALL print_clock( 'iso_raxis' )
CALL errore('integrate_eliashberg_iso_raxis','converged was not reached',1)
RETURN
ENDIF
!
ENDDO ! itemp
filgap = TRIM(prefix) // '.gap'
OPEN(iufilgap, file=filgap, status='unknown')
DO itemp = 1, nstemp ! loop over temperature
WRITE(iufilgap,'(2f12.6)') estemp(itemp)/kelvin2eV, gap(itemp)
ENDDO
CLOSE(iufilgap)
!
CALL stop_clock( 'iso_raxis' )
!
RETURN
!
END SUBROUTINE eliashberg_iso_raxis
!
!-----------------------------------------------------------------------
SUBROUTINE integrate_eliashberg_iso_raxis( itemp, iter, conv )
!-----------------------------------------------------------------------
!!
!! This routine solves the isotropic Eliashberg equations on the real-axis
!!
!
USE kinds, ONLY : DP
USE io_epw, ONLY : iufilker, iufilgap
USE io_global, ONLY : stdout
USE io_files, ONLY : prefix
USE epwcom, ONLY : nswfc, nqstep, nsiter, muc, conv_thr_raxis, &
kerwrite, kerread, nstemp
USE eliashbergcom, ONLY : nsw, estemp, ws, dws, gap0, gap, fdwp, Kp, Km, &
Delta, Deltap, Znorm
USE constants_epw, ONLY : kelvin2eV, ci
USE superconductivity_iso, ONLY : kernel_raxis
!
IMPLICIT NONE
!
INTEGER, INTENT (in) :: itemp
!! Counter on temperature
INTEGER, INTENT(in) :: iter
!! Counter on iteration steps
LOGICAL, INTENT(inout) :: conv
!! True if the calculation is converged
!
! Local variables
INTEGER :: iw, iwp
REAL(DP) :: dstep, a, b, c, d, absdelta, reldelta, errdelta, temp
REAL(DP), ALLOCATABLE :: wesqrt(:), desqrt(:)
REAL(DP) :: eps=1.0d-6
COMPLEX(DP) :: kernelp, kernelm, esqrt
COMPLEX(DP), ALLOCATABLE, SAVE :: Deltaold(:)
REAL(DP), EXTERNAL :: wgauss
LOGICAL :: lgap
CHARACTER(len=256) :: name1, name2
!
IF ( .not. ALLOCATED(wesqrt) ) ALLOCATE( wesqrt(nsw) )
IF ( .not. ALLOCATED(desqrt) ) ALLOCATE( desqrt(nsw) )
!
IF ( iter .eq. 1 ) THEN
IF ( .not. ALLOCATED(gap) ) ALLOCATE( gap(nstemp) )
IF ( .not. ALLOCATED(Delta) ) ALLOCATE( Delta(nsw) )
IF ( .not. ALLOCATED(Deltap) ) ALLOCATE( Deltap(nsw) )
IF ( .not. ALLOCATED(Znorm) ) ALLOCATE( Znorm(nsw) )
gap(itemp) = 0.d0
Deltap(:) = (0.d0, 0.d0)
Deltap(:) = gap0
IF ( .not. ALLOCATED(fdwp) ) ALLOCATE( fdwp(nsw) )
IF ( .not. ALLOCATED(Kp) ) ALLOCATE( Kp(nsw,nsw) )
IF ( .not. ALLOCATED(Km) ) ALLOCATE( Km(nsw,nsw) )
ENDIF
Delta(:) = (0.d0, 0.d0)
Znorm(:) = (0.d0, 0.d0)
!
temp = estemp(itemp) / kelvin2eV
IF ( temp .lt. 10.d0 ) THEN
WRITE(name2,'(a,a7,f4.2)') TRIM(prefix),'.ker_00', temp
ELSEIF ( temp .ge. 10.d0 ) THEN
WRITE(name2,'(a,a6,f5.2)') TRIM(prefix),'.ker_0', temp
ELSEIF ( temp .ge. 100.d0 ) THEN
WRITE(name2,'(a,a5,f6.2)') TRIM(prefix),'.ker_', temp
ENDIF
OPEN(iufilker, file=name2, form='unformatted')
!
IF ( iter .eq. 1 ) THEN
IF ( .not. ALLOCATED(Deltaold) ) ALLOCATE( Deltaold(nsw) )
Deltaold(:) = gap0
ENDIF
absdelta = 0.d0
reldelta = 0.d0
DO iw = 1, nsw ! loop over omega
DO iwp = 1, nsw ! loop over omega_prime
IF ( iter .eq. 1 ) THEN
IF ( iw .eq. 1 ) THEN
IF ( ABS(estemp(itemp)) < eps ) THEN
fdwp(iwp) = 0.d0
ELSE
fdwp(iwp) = wgauss( -ws(iwp) / estemp(itemp), -99 )
ENDIF
ENDIF
!
! read the kernels from file if they were calculated before otherwise calculate them
IF ( kerread ) THEN
READ(iufilker) a, b, c, d
Kp(iw,iwp) = a + ci*b
Km(iw,iwp) = c + ci*d
ENDIF
IF ( kerwrite ) THEN
CALL kernel_raxis( iw, iwp, itemp, kernelp, kernelm )
Kp(iw,iwp) = kernelp
Km(iw,iwp) = kernelm
WRITE(iufilker) real(Kp(iw,iwp)), aimag(Kp(iw,iwp)), &
real(Km(iw,iwp)), aimag(Km(iw,iwp))
ENDIF
ENDIF
!
! this step is performed at each iter step only for iw=1 since it is independ of w(iw)
IF ( iw .eq. 1 ) THEN
esqrt = 1.d0 / sqrt( ws(iwp)**2.d0 - Deltap(iwp)**2.d0 )
wesqrt(iwp) = real( ws(iwp) * esqrt )
desqrt(iwp) = real( Deltap(iwp) * esqrt )
ENDIF
!
! end points contribute only half ( trapezoidal integration rule )
IF ( (iwp .eq. 1) .OR. (iwp .eq. nsw) ) THEN
dstep = 0.5d0 * dws(iwp)
! boundary points contribute half from left and half from right side
ELSEIF ( iwp .eq. nswfc ) THEN
dstep = 0.5d0 * ( dws(iwp) + dws(iwp+1) )
ELSE
dstep = dws(iwp)
ENDIF
Znorm(iw) = Znorm(iw) + dstep * wesqrt(iwp) * Km(iw,iwp)
Delta(iw) = Delta(iw) + dstep * desqrt(iwp) &
* ( Kp(iw,iwp) - muc*( 1.d0 - 2.d0*fdwp(iwp) ) )
ENDDO ! iwp
Znorm(iw) = 1.d0 - Znorm(iw) / ws(iw)
Delta(iw) = Delta(iw) / Znorm(iw)
reldelta = reldelta + abs( Delta(iw) - Deltaold(iw) ) * dws(iw)
absdelta = absdelta + abs( Delta(iw) ) * dws(iw)
ENDDO ! iw
CLOSE(iufilker)
errdelta = reldelta / absdelta
Deltaold(:) = Delta(:)
!
WRITE(stdout,'(5x,a,i6,a,ES20.10,a,ES20.10,a,ES20.10)') 'iter = ', iter, ' error = ', errdelta, &
' Re[Znorm(1)] = ', real(Znorm(1)), ' Re[Delta(1)] = ', real(Delta(1))
!
IF ( errdelta .lt. conv_thr_raxis) conv = .true.
IF ( errdelta .lt. conv_thr_raxis .OR. iter .eq. nsiter ) THEN
IF ( temp .lt. 10.d0 ) THEN
WRITE(name1,'(a,a8,f4.2)') TRIM(prefix),'.gapr_00', temp
ELSEIF ( temp .ge. 10.d0 ) THEN
WRITE(name1,'(a,a7,f5.2)') TRIM(prefix),'.gapr_0', temp
ELSEIF ( temp .ge. 100.d0 ) THEN
WRITE(name1,'(a,a6,f6.2)') TRIM(prefix),'.gapr_', temp
ENDIF
OPEN(iufilgap, file=name1, form='formatted')
!
WRITE(iufilgap,'(5a18)') 'w', 'Re[Znorm(w)]', 'Im[Znorm(w)]', 'Re[Delta(w)]', 'Im[Delta(w)]'
lgap = .true.
! DO iw = 1, nsw
DO iw = 1, nsw-1 ! this change is to prevent segfault in Delta(iw+1) and ws(iw+1)
IF ( lgap .AND. iw .lt. (nqstep) .AND. real(Delta(iw)) .gt. 0.d0 .AND. real(Delta(iw+1)) .gt. 0.d0 .AND. &
( ws(iw) - real(Delta(iw)) )*( ws(iw+1) - real(Delta(iw+1)) ) .lt. 0.d0 ) THEN
gap(itemp) = ( ( real(Delta(iw)) - ws(iw) ) * ws(iw+1) - ( real(Delta(iw+1)) - ws(iw+1) ) * ws(iw) ) &
/ ( ( real(Delta(iw)) - ws(iw) ) - ( real(Delta(iw+1)) - ws(iw+1) ) )
lgap = .false.
ENDIF
WRITE(iufilgap,'(5ES20.10)') ws(iw), real(Znorm(iw)), aimag(Znorm(iw)), &
real(Delta(iw)), aimag(Delta(iw))
ENDDO
CLOSE(iufilgap)
IF ( lgap ) &
gap(itemp) = real(Delta(1))
gap0 = gap(itemp)
ENDIF
!
IF( ALLOCATED(wesqrt) ) DEALLOCATE(wesqrt)
IF( ALLOCATED(desqrt) ) DEALLOCATE(desqrt)
!
IF ( conv .OR. iter .eq. nsiter ) THEN
IF( ALLOCATED(Deltaold) ) DEALLOCATE(Deltaold)
ENDIF
IF ( .not. conv .AND. iter .eq. nsiter ) THEN
WRITE(stdout,'(5x,a,i6)') 'Convergence was not reached in nsiter = ', iter
CALL errore('integrate_eliashberg_iso_raxis','increase nsiter or reduce conv_thr_raxis',-1)
ENDIF
!
RETURN
!
END SUBROUTINE integrate_eliashberg_iso_raxis
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE dos_quasiparticle( itemp )
!-----------------------------------------------------------------------
!!
@ -1117,116 +499,9 @@
RETURN
!
END SUBROUTINE free_energy
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE eliashberg_write_cont_raxis( itemp, cname )
!-----------------------------------------------------------------------
!
!
! This routine writes to files results from the solutions of the Eliashberg
! equations on the real-axis
!
USE kinds, ONLY : DP
USE io_epw, ONLY : iufilgap
USE io_files, ONLY : prefix
USE control_flags, ONLY : iverbosity
USE epwcom, ONLY : nqstep, fsthick, laniso, liso
USE eliashbergcom, ONLY : nsw, estemp, ws, gap, Agap, &
Delta, Znorm, ADelta, AZnorm, &
nkfs, nbndfs, ef0, ekfs
USE constants_epw, ONLY : kelvin2eV
USE io_eliashberg, ONLY : gap_distribution_FS
!
IMPLICIT NONE
!
INTEGER :: iw, itemp, ik, ibnd
REAL(DP) :: temp
LOGICAL :: lgap
CHARACTER(len=256) :: name1, cname
!
temp = estemp(itemp) / kelvin2eV
!
IF ( laniso ) THEN
IF ( iverbosity .eq. 2 ) THEN
IF ( temp .lt. 10.d0 ) THEN
WRITE(name1,'(a,a1,a4,a9,f4.2)') TRIM(prefix), '.', cname, '_aniso_00', temp
ELSEIF ( temp .ge. 10.d0 .AND. temp .lt. 100.d0 ) THEN
WRITE(name1,'(a,a1,a4,a8,f5.2)') TRIM(prefix), '.', cname, '_aniso_0', temp
ELSEIF ( temp .ge. 100.d0 ) THEN
WRITE(name1,'(a,a1,a4,a7,f6.2)') TRIM(prefix), '.', cname, '_aniso_', temp
ENDIF
OPEN(iufilgap, file=name1, form='formatted')
WRITE(iufilgap,'(6a20)') '# w [eV]', 'Enk-Ef [eV]', 'Re[Znorm(w)]', 'Im[Znorm(w)]',&
'Re[Delta(w)] [eV]', 'Im[Delta(w)] [eV]'
ENDIF
!
DO ik = 1, nkfs
DO ibnd = 1, nbndfs
IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
lgap = .true.
! DO iw = 1, nsw
DO iw = 1, nsw-1 ! FG: this change is to prevent segfault in ws(iw+1) and ADelta(*,*,iw+1)
IF ( lgap .AND. iw .lt. nqstep .AND. real(ADelta(ibnd,ik,iw)) .gt. 0.d0 &
.AND. real(ADelta(ibnd,ik,iw+1)) .gt. 0.d0 &
.AND. ( ws(iw) - real(ADelta(ibnd,ik,iw)) )*( ws(iw+1) - real(ADelta(ibnd,ik,iw+1)) ) .lt. 0.d0 ) THEN
Agap(ibnd,ik,itemp) = ( ( real(ADelta(ibnd,ik,iw)) - ws(iw) ) * ws(iw+1) &
- ( real(ADelta(ibnd,ik,iw+1)) - ws(iw+1) ) * ws(iw) ) &
/ ( ( real(ADelta(ibnd,ik,iw)) - ws(iw) ) - ( real(ADelta(ibnd,ik,iw+1)) - ws(iw+1) ) )
lgap = .false.
ENDIF
IF ( iverbosity .eq. 2 ) THEN
WRITE(iufilgap,'(6ES20.10)') ws(iw), ekfs(ibnd,ik)-ef0, &
real(AZnorm(ibnd,ik,iw)), aimag(AZnorm(ibnd,ik,iw)), &
real(ADelta(ibnd,ik,iw)), aimag(ADelta(ibnd,ik,iw))
ENDIF
ENDDO ! iw
IF ( lgap ) &
Agap(ibnd,ik,itemp) = real(ADelta(ibnd,ik,1))
ENDIF
ENDDO ! ibnd
ENDDO ! ik
IF ( iverbosity .eq. 2 ) CLOSE(iufilgap)
!
CALL gap_distribution_FS ( itemp, cname )
!
ENDIF
!
! isotropic case
! SP: Only write isotropic for laniso if user really wants that
IF ( ( laniso .AND. iverbosity .eq. 2 ) .OR. liso ) THEN
IF ( temp .lt. 10.d0 ) THEN
WRITE(name1,'(a,a1,a4,a7,f4.2)') TRIM(prefix), '.', cname, '_iso_00', temp
ELSEIF ( temp .ge. 10.d0 .AND. temp .lt. 100.d0 ) THEN
WRITE(name1,'(a,a1,a4,a6,f5.2)') TRIM(prefix), '.', cname, '_iso_0', temp
ELSEIF ( temp .ge. 100.d0 ) THEN
WRITE(name1,'(a,a1,a4,a5,f6.2)') TRIM(prefix), '.', cname, '_iso_', temp
ENDIF
OPEN(iufilgap, file=name1, form='formatted')
WRITE(iufilgap,'(5a20)') 'w [eV]', 'Re[Znorm(w)]', 'Im[Znorm(w)]', 'Re[Delta(w)] [eV]', 'Im[Delta(w)] [eV]'
lgap = .true.
! DO iw = 1, nsw
DO iw = 1, nsw-1 ! this change is to prevent segfault in Delta(iw+1) and ws(iw+1)
IF ( lgap .AND. iw .lt. nqstep .AND. real(Delta(iw)) .gt. 0.d0 .AND. real(Delta(iw+1)) .gt. 0.d0 .AND. &
( ws(iw) - real(Delta(iw)) )*( ws(iw+1) - real(Delta(iw+1)) ) .lt. 0.d0 ) THEN
gap(itemp) = ( ( real(Delta(iw)) - ws(iw) ) * ws(iw+1) - ( real(Delta(iw+1)) - ws(iw+1) ) * ws(iw) ) &
/ ( ( real(Delta(iw)) - ws(iw) ) - ( real(Delta(iw+1)) - ws(iw+1) ) )
lgap = .false.
ENDIF
WRITE(iufilgap,'(5ES20.10)') ws(iw), real(Znorm(iw)), aimag(Znorm(iw)), &
real(Delta(iw)), aimag(Delta(iw))
ENDDO ! iw
CLOSE(iufilgap)
IF ( lgap ) &
gap(itemp) = real(Delta(1))
ENDIF
!
RETURN
!
END SUBROUTINE eliashberg_write_cont_raxis
!
!-----------------------------------------------------------------------
SUBROUTINE eliashberg_memlt_aniso_iaxis( itemp )
!-----------------------------------------------------------------------
!!

361
EPW/src/superconductivity_aniso.f90
View File

@ -1397,5 +1397,366 @@
RETURN
!
END SUBROUTINE evaluate_a2fij
!
!-----------------------------------------------------------------------
SUBROUTINE evaluate_a2f_lambda
!-----------------------------------------------------------------------
!
! computes the isotropic spectral function a2F(w), total lambda, and
! distribution of lambda
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE io_epw, ONLY : iua2ffil, iudosfil, iufillambda, iufillambdaFS
USE io_files, ONLY : prefix
USE phcom, ONLY : nmodes
USE cell_base, ONLY : bg
USE control_flags, ONLY : iverbosity
USE elph2, ONLY : nqtotf, wqf, wf
USE epwcom, ONLY : fsthick, eps_acustic, nqstep, degaussq, delta_qsmear, nqsmear, &
degaussw, nkf1, nkf2, nkf3
USE eliashbergcom, ONLY : nkfs, nbndfs, g2, ixkqf, ixqfs, nqfs, w0g, ekfs, ef0, dosef, wsph, &
wkfs, dwsph, a2f_iso, ixkff
USE constants_epw, ONLY : ryd2ev
USE io_global, ONLY : ionode_id
USE mp_global, ONLY : inter_pool_comm, my_pool_id, npool
USE mp_world, ONLY : mpime
USE mp, ONLY : mp_bcast, mp_barrier, mp_sum
!
IMPLICIT NONE
!
INTEGER :: ik, iq, iq0, iwph, ibnd, jbnd, imode, lower_bnd, upper_bnd, &
ismear, ibin, nbin, nbink, i, j, k
REAL(DP) :: weight, weightq, l_sum, lambda_eph, lambda_max(npool), &
sigma, dbin, dbink, x1, x2, x3
REAL(DP), ALLOCATABLE :: a2f(:,:), phdos(:,:), l_a2f(:), lambda_k(:,:), &
lambda_k_bin(:), lambda_pairs(:), a2f_modeproj(:,:), phdos_modeproj(:,:)
REAL(DP), EXTERNAL :: w0gauss
CHARACTER (len=256) :: name1
!
! This is only a quick fix since the subroutine was written for parallel execution - FG June 2014
#if ! defined(__MPI)
npool = 1
my_pool_id = 0
#endif
!
! degaussq is read from the input file in meV and converted to Ryd in epw_readin.f90
! go from Ryd to eV
degaussq = degaussq * ryd2ev
delta_qsmear = delta_qsmear * ryd2ev
!
CALL fkbounds( nkfs, lower_bnd, upper_bnd )
!
IF ( .not. ALLOCATED(a2f_iso) ) ALLOCATE( a2f_iso(nqstep) )
IF ( .not. ALLOCATED(a2f) ) ALLOCATE( a2f(nqstep,nqsmear) )
IF ( .not. ALLOCATED(a2f_modeproj) ) ALLOCATE( a2f_modeproj(nmodes,nqstep) )
a2f_iso(:) = 0.d0
a2f(:,:) = 0.d0
a2f_modeproj(:,:) = 0.d0
!
! RM - the 0 index in k is required when printing out values of lambda_k
! When the k-point is outside the Fermi shell, ixkff(ik)=0
IF ( .not. ALLOCATED(lambda_k) ) ALLOCATE(lambda_k(0:nkfs,nbndfs))
lambda_k(:,:) = 0.d0
!
l_sum = 0.d0
lambda_max(:) = 0.d0
DO ismear = 1, nqsmear
sigma = degaussq + (ismear-1) * delta_qsmear
DO ik = lower_bnd, upper_bnd
DO ibnd = 1, nbndfs
IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
DO iq = 1, nqfs(ik)
! iq0 - index of q-point on the full q-mesh
iq0 = ixqfs(ik,iq)
DO jbnd = 1, nbndfs
IF ( abs( ekfs(jbnd,ixkqf(ik,iq0)) - ef0 ) .lt. fsthick ) THEN
weight = wkfs(ik) * wqf(iq) * w0g(ibnd,ik) * w0g(jbnd,ixkqf(ik,iq0))
lambda_eph = 0.d0
DO imode = 1, nmodes
IF ( wf(imode,iq0) .gt. eps_acustic ) THEN
IF ( ismear .eq. 1 ) THEN
lambda_eph = lambda_eph + g2(ik,iq,ibnd,jbnd,imode) / wf(imode,iq0)
ENDIF
DO iwph = 1, nqstep
weightq = w0gauss( ( wsph(iwph) - wf(imode,iq0) ) / sigma, 0 ) / sigma
a2f(iwph,ismear) = a2f(iwph,ismear) + weight * weightq * g2(ik,iq,ibnd,jbnd,imode)
IF ( ismear .eq. 1 ) THEN
a2f_modeproj(imode,iwph) = a2f_modeproj(imode,iwph) +&
weight * weightq * g2(ik,iq,ibnd,jbnd,imode)
ENDIF
ENDDO ! iwph
ENDIF ! wf
ENDDO ! imode
IF ( ismear .eq. 1 .AND. lambda_eph .gt. 0.d0 ) THEN
l_sum = l_sum + weight * lambda_eph
weight = wqf(iq) * w0g(jbnd,ixkqf(ik,iq0))
lambda_k(ik,ibnd) = lambda_k(ik,ibnd) + weight * lambda_eph
IF ( lambda_eph .gt. lambda_max(my_pool_id+1) ) THEN
lambda_max(my_pool_id+1) = lambda_eph
ENDIF
ENDIF
ENDIF ! ekq
ENDDO ! jbnd
ENDDO ! iq
ENDIF ! ekk
ENDDO ! ibnd
ENDDO ! ik
ENDDO ! ismear
!
a2f(:,:) = 0.5d0 * a2f(:,:) / dosef
a2f_modeproj(:,:) = 0.5d0 * a2f_modeproj(:,:) / dosef
l_sum = l_sum / dosef
lambda_k(:,:) = 2.d0 * lambda_k(:,:)
lambda_max(:) = 2.d0 * dosef * lambda_max(:)
!
! collect contributions from all pools (sum over k-points)
CALL mp_sum( l_sum, inter_pool_comm )
CALL mp_sum( a2f, inter_pool_comm )
CALL mp_sum( a2f_modeproj, inter_pool_comm )
CALL mp_sum( lambda_max, inter_pool_comm )
CALL mp_sum( lambda_k, inter_pool_comm )
CALL mp_barrier(inter_pool_comm)
!
IF ( mpime .eq. ionode_id ) THEN
!
OPEN( unit = iua2ffil, file = TRIM(prefix)//".a2f", form = 'formatted')
OPEN( unit = iudosfil, file = TRIM(prefix)//".phdos", form = 'formatted')
!
IF ( .not. ALLOCATED(phdos) ) ALLOCATE( phdos(nqstep,nqsmear) )
IF ( .not. ALLOCATED(phdos_modeproj) ) ALLOCATE( phdos_modeproj(nmodes,nqstep) )
phdos(:,:) = 0.d0
phdos_modeproj(:,:) = 0.d0
!
DO ismear = 1, nqsmear
sigma = degaussq + (ismear-1) * delta_qsmear
DO iq = 1, nqtotf
DO imode = 1, nmodes
IF ( wf(imode,iq) .gt. eps_acustic ) THEN
DO iwph = 1, nqstep
weightq = w0gauss( ( wsph(iwph) - wf(imode,iq)) / sigma, 0 ) / sigma
phdos(iwph,ismear) = phdos(iwph,ismear) + wqf(iq) * weightq
IF ( ismear .eq. 1 ) THEN
phdos_modeproj(imode,iwph) = phdos_modeproj(imode,iwph) + wqf(iq) * weightq
ENDIF
ENDDO ! iwph
ENDIF ! wf
ENDDO ! imode
ENDDO ! iq
ENDDO ! ismear
!
IF ( .not. ALLOCATED(l_a2f) ) ALLOCATE( l_a2f(nqsmear) )
l_a2f(:) = 0.d0
!
DO ismear = 1, nqsmear
DO iwph = 1, nqstep
l_a2f(ismear) = l_a2f(ismear) + a2f(iwph,ismear) / wsph(iwph)
! wsph in meV (from eV) and phdos in states/meV (from states/eV)
IF (ismear .eq. nqsmear) WRITE (iua2ffil,'(f12.7,15f12.7)') wsph(iwph)*1000.d0, a2f(iwph,:)
IF (ismear .eq. nqsmear) WRITE (iudosfil,'(f12.7,15f15.7)') wsph(iwph)*1000.d0, phdos(iwph,:)/1000.d0
ENDDO
l_a2f(ismear) = 2.d0 * l_a2f(ismear) * dwsph
ENDDO
!
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)*1000.d0,ismear=1,nqsmear )
WRITE(iua2ffil,'(" Electron smearing (eV)", f12.7)') degaussw
WRITE(iua2ffil,'(" Fermi window (eV)", f12.7)') fsthick
WRITE(iua2ffil,'(" Summed el-ph coupling ", f12.7)') l_sum
CLOSE(iua2ffil)
CLOSE(iudosfil)
!
a2f_iso(:) = a2f(:,1)
OPEN( unit = iua2ffil, file = TRIM(prefix)//".a2f_iso", form = 'formatted')
OPEN( unit = iudosfil, file = TRIM(prefix)//".phdos_proj", form = 'formatted')
DO iwph = 1, nqstep
! wsph in meV (from eV) and phdos in states/meV (from states/eV)
WRITE(iua2ffil,'(f12.7,100f12.7)') wsph(iwph)*1000.d0, a2f_iso(iwph), a2f_modeproj(:,iwph)
WRITE(iudosfil,'(f12.7,100f15.7)') wsph(iwph)*1000.d0, phdos(iwph,1)/1000.d0, phdos_modeproj(:,iwph)/1000.d0
ENDDO
WRITE(iua2ffil,'(a,f18.7,a,f18.7)') 'lambda_int = ', l_a2f(1), ' lambda_sum = ',l_sum
CLOSE(iua2ffil)
CLOSE(iudosfil)
!
IF ( ALLOCATED(phdos) ) DEALLOCATE( phdos )
IF ( ALLOCATED(phdos_modeproj) ) DEALLOCATE( phdos_modeproj )
IF ( ALLOCATED(l_a2f) ) DEALLOCATE( l_a2f )
!
ENDIF
!
CALL mp_bcast( a2f_iso, ionode_id, inter_pool_comm )
CALL mp_barrier(inter_pool_comm)
!
IF ( ALLOCATED(a2f) ) DEALLOCATE( a2f )
IF ( ALLOCATED(a2f_modeproj) ) DEALLOCATE( a2f_modeproj )
!
nbink = int( 1.25d0 * maxval(lambda_k(:,:)) / 0.005d0 )
dbink = 1.25d0 * maxval(lambda_k(:,:)) / dble(nbink)
IF ( .not. ALLOCATED(lambda_k_bin) ) ALLOCATE ( lambda_k_bin(nbink) )
lambda_k_bin(:) = 0.d0
!
!SP : Should be initialized
nbin = 0
dbin = 0.0_DP
!
IF ( iverbosity == 2 ) THEN
nbin = int( 1.25d0 * maxval(lambda_max(:)) / 0.005d0 )
dbin = 1.25d0 * maxval(lambda_max(:)) / dble(nbin)
IF ( .not. ALLOCATED(lambda_pairs) ) ALLOCATE ( lambda_pairs(nbin) )
lambda_pairs(:) = 0.d0
ENDIF
!
WRITE(stdout,'(5x,a13,f21.7,a18,f21.7)') 'lambda_max = ', maxval(lambda_max(:)), ' lambda_k_max = ', maxval(lambda_k(:,:))
WRITE(stdout,'(a)') ' '
!
lambda_k(:,:) = 0.d0
DO ik = lower_bnd, upper_bnd
DO ibnd = 1, nbndfs
IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
DO iq = 1, nqfs(ik)
! iq0 - index of q-point on the full q-mesh
iq0 = ixqfs(ik,iq)
DO jbnd = 1, nbndfs
IF ( abs( ekfs(jbnd,ixkqf(ik,iq0)) - ef0 ) .lt. fsthick ) THEN
weight = wqf(iq) * w0g(jbnd,ixkqf(ik,iq0)) / dosef
CALL lambdar_aniso_ver1( ik, iq, ibnd, jbnd, 0.d0, lambda_eph )
lambda_k(ik,ibnd) = lambda_k(ik,ibnd) + weight * lambda_eph
IF ( iverbosity == 2 ) THEN
DO ibin = 1, nbin
sigma = 1.d0 * dbin
weight = w0gauss( ( lambda_eph - dble(ibin) * dbin ) / sigma, 0 ) / sigma
lambda_pairs(ibin) = lambda_pairs(ibin) + weight
ENDDO
ENDIF
ENDIF
ENDDO ! jbnd
ENDDO ! iq
DO ibin = 1, nbink
sigma = 1.d0 * dbink
weight = w0gauss( ( lambda_k(ik,ibnd) - dble(ibin) * dbink ) / sigma, 0 ) / sigma
lambda_k_bin(ibin) = lambda_k_bin(ibin) + weight
ENDDO
ENDIF
ENDDO ! ibnd
ENDDO ! ik
!
! collect contributions from all pools
CALL mp_sum( lambda_k, inter_pool_comm )
IF ( iverbosity .eq. 2 ) THEN
CALL mp_sum( lambda_pairs, inter_pool_comm )
ENDIF
CALL mp_sum( lambda_k_bin, inter_pool_comm )
CALL mp_barrier(inter_pool_comm)
!
IF ( mpime .eq. ionode_id ) THEN
!
! SP: Produced if user really wants it
IF ( iverbosity .eq. 2 ) THEN
OPEN(unit = iufillambda, file = TRIM(prefix)//".lambda_aniso", form = 'formatted')
WRITE(iufillambda,'(2a12,2a7)') '# enk-e0[eV]',' lambda_nk','# kpt','# band'
DO ik = 1, nkfs
DO ibnd = 1, nbndfs
IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
WRITE(iufillambda,'(2f12.7,2i7)') ekfs(ibnd,ik) - ef0, lambda_k(ik,ibnd), ik, ibnd
ENDIF
ENDDO
ENDDO
CLOSE(iufillambda)
ENDIF
!
OPEN(unit = iufillambda, file = TRIM(prefix)//".lambda_k_pairs", form = 'formatted')
WRITE(iufillambda,'(a12,a30)') '# lambda_nk',' \rho(lambda_nk) scaled to 1'
DO ibin = 1, nbink
WRITE(iufillambda,'(2f21.7)') dbink*dble(ibin), lambda_k_bin(ibin)/maxval(lambda_k_bin(:))
ENDDO
CLOSE(iufillambda)
!
! SP: Produced if user really wants it
IF ( iverbosity == 2 ) THEN
OPEN( unit = iufillambda, file = TRIM(prefix)//".lambda_pairs", form = 'formatted')
WRITE(iufillambda,'(a12,a30)') "# lambda_nk,n'k'", " \rho(lambda_nk,n'k') scaled to 1"
DO ibin = 1, nbin
WRITE(iufillambda,'(2f21.7)') dbin*dble(ibin), lambda_pairs(ibin)/maxval(lambda_pairs(:))
ENDDO
CLOSE(iufillambda)
ENDIF
!
! SP & RM: .cube file for VESTA plotting (only if iverbosity = 2)
!
! RM - If the k-point is outside the Fermi shell,
! ixkff(ik)=0 and lambda_k(0,ibnd) = 0.0
!
IF ( iverbosity .eq. 2 ) THEN
!
DO ibnd = 1, nbndfs
!
IF ( ibnd < 10 ) THEN
WRITE(name1,'(a,a8,i1,a5)') TRIM(prefix),'.lambda_', ibnd, '.cube'
ELSEIF ( ibnd < 100 ) THEN
WRITE(name1,'(a,a8,i2,a5)') TRIM(prefix),'.lambda_', ibnd, '.cube'
ELSEIF( ibnd < 1000 ) THEN
WRITE(name1,'(a,a8,i3,a5)') TRIM(prefix),'.lambda_', ibnd, '.cube'
ELSE
CALL errore( 'eliashberg_setup', 'Too many bands ',1)
ENDIF
!
OPEN(iufillambdaFS, file=name1, form='formatted')
WRITE(iufillambdaFS,*) 'Cubfile created from EPW calculation'
WRITE(iufillambdaFS,*) 'lambda'
WRITE(iufillambdaFS,'(i5,3f12.6)') 1, 0.0d0, 0.0d0, 0.0d0
WRITE(iufillambdaFS,'(i5,3f12.6)') nkf1, (bg(i,1)/DBLE(nkf1),i=1,3)
WRITE(iufillambdaFS,'(i5,3f12.6)') nkf2, (bg(i,2)/DBLE(nkf2),i=1,3)
WRITE(iufillambdaFS,'(i5,3f12.6)') nkf3, (bg(i,3)/DBLE(nkf3),i=1,3)
WRITE(iufillambdaFS,'(i5,4f12.6)') 1, 1.0d0, 0.0d0, 0.0d0, 0.0d0
WRITE(iufillambdaFS,'(6f12.6)') ( lambda_k(ixkff(ik),ibnd), ik=1,nkf1*nkf2*nkf3 )
CLOSE(iufillambdaFS)
ENDDO
!
ENDIF
!
! SP & RM : Write on file the lambda close to the Fermi surface along with
! Cartesian coordinate, band index, energy distance from Fermi level
! and lambda value.
!
OPEN(unit = iufillambdaFS, file = TRIM(prefix)//".lambda_FS", form='formatted')
WRITE(iufillambdaFS,'(a75)') '# k-point Band Enk-Ef [eV] lambda'
DO i = 1, nkf1
DO j = 1, nkf2
DO k = 1, nkf3
ik = k + (j-1)*nkf3 + (i-1)*nkf2*nkf3
IF ( ixkff(ik) .gt. 0 ) THEN
DO ibnd = 1, nbndfs
! SP: Here take a 0.2 eV interval around the FS.
IF ( abs( ekfs(ibnd,ixkff(ik)) - ef0 ) .lt. fsthick ) THEN
!IF ( abs( ekfs(ibnd,ixkff(ik)) - ef0 ) .lt. 0.2 ) THEN
x1 = bg(1,1)*(i-1)/nkf1+bg(1,2)*(j-1)/nkf2+bg(1,3)*(k-1)/nkf3
x2 = bg(2,1)*(i-1)/nkf1+bg(2,2)*(j-1)/nkf2+bg(2,3)*(k-1)/nkf3
x3 = bg(3,1)*(i-1)/nkf1+bg(3,2)*(j-1)/nkf2+bg(3,3)*(k-1)/nkf3
WRITE(iufillambdaFS,'(3f12.6,i8,f12.6,f24.15)') x1, x2, x3, ibnd, &
ekfs(ibnd,ixkff(ik))-ef0, lambda_k(ixkff(ik),ibnd)
ENDIF
ENDDO ! ibnd
ENDIF
ENDDO ! k
ENDDO ! j
ENDDO ! i
CLOSE(iufillambdaFS)
ENDIF
CALL mp_barrier(inter_pool_comm)
!
IF ( ALLOCATED(lambda_k) ) DEALLOCATE(lambda_k)
IF ( ALLOCATED(lambda_pairs) ) DEALLOCATE(lambda_pairs)
IF ( ALLOCATED(lambda_k_bin) ) DEALLOCATE(lambda_k_bin)
!
RETURN
!
END SUBROUTINE evaluate_a2f_lambda
!
!
END MODULE superconductivity_aniso

2
test-suite/userconfig.tmp
View File

@ -64,7 +64,7 @@ tolerance = ( (1.0e-6, 5.0e-3, 'e1'),
(1.0e-4, 1.0e-5, 'qdir'),
(5.0e-2, 5.0e-3, 'diel'),
(5.0e-2, 5.0e-3, 'born'),
(2.0e+0, 8.0e-2, 'phfreq'),
(2.0e+0, 2.5e-1, 'phfreq'),
(1.0e-5, 1.0e-5, 'q1'),
(1.0e-5, 1.0e-5, 'dos1'),
(1.0e-3, 5.0e-3, 'e2'),