diff --git a/EPW/src/Makefile b/EPW/src/Makefile index da90fd3bb..29385c947 100644 --- a/EPW/src/Makefile +++ b/EPW/src/Makefile @@ -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 \ diff --git a/EPW/src/eliashberg.f90 b/EPW/src/eliashberg.f90 index dee7d7bf1..c03bbfcba 100644 --- a/EPW/src/eliashberg.f90 +++ b/EPW/src/eliashberg.f90 @@ -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 ! diff --git a/EPW/src/io_eliashberg.f90 b/EPW/src/io_eliashberg.f90 index f4bc826a6..8c8d85c92 100644 --- a/EPW/src/io_eliashberg.f90 +++ b/EPW/src/io_eliashberg.f90 @@ -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 !----------------------------------------------------------------------- !! diff --git a/EPW/src/make.depend b/EPW/src/make.depend index 91beaa12c..507dd42aa 100644 --- a/EPW/src/make.depend +++ b/EPW/src/make.depend @@ -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 diff --git a/EPW/src/superconductivity.f90 b/EPW/src/superconductivity.f90 index ff79e8f11..139ff3091 100644 --- a/EPW/src/superconductivity.f90 +++ b/EPW/src/superconductivity.f90 @@ -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 ) !----------------------------------------------------------------------- !! diff --git a/EPW/src/superconductivity_aniso.f90 b/EPW/src/superconductivity_aniso.f90 index 04081557f..d2c61ed05 100644 --- a/EPW/src/superconductivity_aniso.f90 +++ b/EPW/src/superconductivity_aniso.f90 @@ -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 diff --git a/test-suite/userconfig.tmp b/test-suite/userconfig.tmp index 2206b0db1..fc8b0f0a5 100644 --- a/test-suite/userconfig.tmp +++ b/test-suite/userconfig.tmp @@ -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'),