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

Final cleaning and modularization.

This commit is contained in:
Samuel Ponce 2019-09-13 17:28:34 +01:00
parent 27878e7f31
commit 8392d54d2d
15 changed files with 418 additions and 357 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
EPW/src/Makefile
View File

@ -26,6 +26,7 @@ elph2.o \
io_epw.o \
low_lvl.o \
division.o \
rigid_epw.o \
wan2bloch.o \
wigner.o \
io_dyn_mat2.o \
@ -33,6 +34,7 @@ readmat.o \
io_scattering.o \
cum_mod.o \
close_epw.o \
poolgathering.o \
printing.o \
kfold.o \
io_eliashberg.o \
@ -40,11 +42,9 @@ eliashbergcom.o \
superconductivity.o \
superconductivity_iso.o \
superconductivity_aniso.o \
poolgathering.o \
grid.o \
transport.o \
transport_iter.o \
rigid_epw.o \
a2f.o \
adddvscf2.o \
bcast_epw_input.o \

6
EPW/src/bloch2wan.f90
View File

@ -578,10 +578,8 @@
!
CALL cryst_to_cart(nq, xk, bg, 1)
!
!
! check spatial decay of dynamical matrix in Wannier basis
! the unit in r-space is angstrom, and I am plotting
! the matrix for the first mode only
! check spatial decay of dynamical matrix in Wannier basis
! the unit in r-space is angstrom, and I am plotting the matrix for the first mode only
!
IF (mpime == ionode_id) THEN
OPEN(UNIT = iudecaydyn, FILE = 'decay.dynmat')

51
EPW/src/rotate_eigenm.f90
View File

@ -140,8 +140,7 @@
!
DO jmode = 1, nmodes
DO imode = 1, jmode
dynp(imode + (jmode - 1) * jmode / 2) = &
(dyn1(imode, jmode) + CONJG(dyn1(jmode, imode))) / 2.d0
dynp(imode + (jmode - 1) * jmode / 2) = (dyn1(imode, jmode) + CONJG(dyn1(jmode, imode))) / 2.d0
ENDDO
ENDDO
!
@ -150,8 +149,7 @@
!
IF (iverbosity == 1) THEN
!
! check the frequencies
!
! check the frequencies
DO nu = 1, nmodes
IF (w1(nu) > 0.d0) THEN
wtmp(nu) = SQRT(ABS(w1(nu)))
@ -160,25 +158,25 @@
ENDIF
ENDDO
WRITE(stdout, '(5x,"Frequencies of the matrix for the first q in the star (cm^{-1})")')
WRITE(stdout, '(6(2x,f10.5))' ) (wtmp(nu) * rydcm1, nu = 1, nmodes)
WRITE(stdout, '(6(2x,f10.5))') (wtmp(nu) * rydcm1, nu = 1, nmodes)
!
ENDIF
!
! here dyn1 is dynq(:,:,iq_first) after dividing by the masses
! (the true dyn matrix D_q)
! here dyn1 is dynq(:,:,iq_first) after dividing by the masses
! (the true dyn matrix D_q)
!
! -----------------------------------------------------------------------
! the matrix gamma (Maradudin & Vosko, RMP, eq. 2.37)
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! the matrix gamma (Maradudin & Vosko, RMP, eq. 2.37)
! -----------------------------------------------------------------------
!
! I have built the matrix by following step-by-step q2star_ph.f90 and
! rotate_and_add_dyn.f90
! I have built the matrix by following step-by-step q2star_ph.f90 and
! rotate_and_add_dyn.f90
!
ism1 = invs(isym)
!
! the symmetry matrix in cartesian coordinates
! (so that we avoid going back and forth with the dynmat)
! note the presence of both at and bg in the transform!
! the symmetry matrix in cartesian coordinates
! (so that we avoid going back and forth with the dynmat)
! note the presence of both at and bg in the transform!
!
scart = DBLE(s(:, :, ism1))
scart = MATMUL(MATMUL(bg, scart), TRANSPOSE(at))
@ -208,22 +206,21 @@
!
! D_{Sq} = gamma * D_q * gamma^\dagger (Maradudin & Vosko, RMP, eq. 3.5)
!
CALL zgemm ('n', 'n', nmodes, nmodes, nmodes, cone, gamma, &
CALL ZGEMM('n', 'n', nmodes, nmodes, nmodes, cone, gamma, &
nmodes, dyn1 , nmodes, czero , dyn2, nmodes)
CALL zgemm ('n', 'c', nmodes, nmodes, nmodes, cone, dyn2, &
CALL ZGEMM('n', 'c', nmodes, nmodes, nmodes, cone, dyn2, &
nmodes, gamma, nmodes, czero, dyn1, nmodes)
!
DO jmode = 1, nmodes
DO imode = 1, jmode
dynp(imode + (jmode - 1) * jmode / 2) = &
(dyn1(imode, jmode) + CONJG(dyn1(jmode, imode))) / 2.d0
dynp(imode + (jmode - 1) * jmode / 2) = (dyn1(imode, jmode) + CONJG(dyn1(jmode, imode))) / 2.d0
ENDDO
ENDDO
!
CALL zhpevx('V', 'A', 'U', nmodes, dynp, 0.0, 0.0, &
CALL ZHPEVX('V', 'A', 'U', nmodes, dynp, 0.0, 0.0, &
0, 0, -1.0, neig, w2, cz2, nmodes, cwork, rwork, iwork, ifail, info)
!
! check the frequencies
! Check the frequencies
!
DO nu = 1, nmodes
IF (w2(nu) > 0.d0) THEN
@ -238,9 +235,9 @@
ENDIF
!
!
! -----------------------------------------------------------------------
! transformation of the eigenvectors: e_{Sq} = gamma * e_q (MV eq. 2.36)
! -----------------------------------------------------------------------
! -----------------------------------------------------------------------
! transformation of the eigenvectors: e_{Sq} = gamma * e_q (MV eq. 2.36)
! -----------------------------------------------------------------------
!
CALL ZGEMM('n', 'n', nmodes, nmodes, nmodes, cone, gamma, &
nmodes, cz1 , nmodes, czero , cz2, nmodes)
@ -256,7 +253,7 @@
ENDDO
!
! the rotated matrix and the one read from file
IF (iverbosity == 1) WRITE(stdout,'(2f15.10)') dyn2 - dyn1
IF (iverbosity == 1) WRITE(stdout, '(2f15.10)') dyn2 - dyn1
!
! here I have checked that the matrix rotated with gamma
! is perfectly equal to the one read from file for this q in the star
@ -273,7 +270,7 @@
!
IF (iverbosity == 1) THEN
!
! a simple check on the frequencies
! A simple check on the frequencies
!
DO nu = 1, nmodes
IF (w2(nu) > 0.d0 ) then
@ -283,7 +280,7 @@
ENDIF
ENDDO
WRITE(stdout, '(5x,"Frequencies of the matrix for the current q in the star (cm^{-1})")')
WRITE(stdout, '(6(2x,f10.5))' ) (wtmp(nu) * rydcm1, nu = 1, nmodes)
WRITE(stdout, '(6(2x,f10.5))') (wtmp(nu) * rydcm1, nu = 1, nmodes)
!
ENDIF
!

17
EPW/src/rotate_epmat.f90
View File

@ -30,6 +30,12 @@
!
IMPLICIT NONE
!
LOGICAL, INTENT(in) :: lwin(nbnd, nks)
!! Bands at k within outer energy window
LOGICAL, INTENT(in) :: lwinq(nbnd, nks)
!! Bands at k+q within outer energy window
LOGICAL, INTENT(in) :: exband(nbnd)
!! Bands excluded from the calculation of overlap and projection matrices
INTEGER, INTENT(in) :: iq
!! Current qpoint
REAL(KIND = DP), INTENT(in) :: xq(3)
@ -38,12 +44,6 @@
!! eigenvectors for the first q in the star
COMPLEX(KIND = DP), INTENT(inout) :: cz2(nmodes, nmodes)
!! Rotated eigenvectors for the current q in the star
LOGICAL, INTENT(in) :: lwin(nbnd, nks)
!! Bands at k within outer energy window
LOGICAL, INTENT(in) :: lwinq(nbnd, nks)
!! Bands at k+q within outer energy window
LOGICAL, INTENT(in) :: exband(nbnd)
!! Bands excluded from the calculation of overlap and projection matrices
!
! Local variables
INTEGER :: mu
@ -167,7 +167,7 @@
! bring e-p matrix from the cartesian representation of the
! first q in the star to the corresponding eigenmode representation
!
CALL zgemv('t', nmodes, nmodes, cone, cz1, nmodes, &
CALL ZGEMV('t', nmodes, nmodes, cone, cz1, nmodes, &
epmatq_opt(ibnd, jbnd, ik, :), 1, czero, eptmp, 1)
!
IF (lpolar) THEN
@ -179,7 +179,7 @@
!
! rotate epmat in the cartesian representation for this q in the star
!
CALL zgemv('t', nmodes, nmodes, cone, cz2, nmodes, &
CALL ZGEMV('t', nmodes, nmodes, cone, cz2, nmodes, &
eptmp, 1, czero, epmatq(ibnd, jbnd, ik, :, iq), 1)
ENDDO
ENDDO
@ -188,4 +188,3 @@
!---------------------------------------------------------------------------
END SUBROUTINE rotate_epmat
!---------------------------------------------------------------------------

232
EPW/src/selfen.f90
View File

@ -72,7 +72,6 @@
!! Total number of q-points from the selecq.fmt grid.
!
! Local variables
!
INTEGER :: n
!! Integer for the degenerate average over eigenstates
INTEGER :: ik
@ -89,6 +88,8 @@
!! Counter on mode
INTEGER :: fermicount
!! Number of states on the Fermi surface
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: tmp
!! Temporary variable to store real part of Sigma for the degenerate average
REAL(KIND = DP) :: tmp2
@ -194,9 +195,9 @@
!
IF (iqq == 1) THEN
!
WRITE(stdout,'(/5x,a)') REPEAT('=',67)
WRITE(stdout,'(5x,"Electron (Imaginary) Self-Energy in the Migdal Approximation")')
WRITE(stdout,'(5x,a/)') REPEAT('=',67)
WRITE(stdout, '(/5x,a)') REPEAT('=',67)
WRITE(stdout, '(5x,"Electron (Imaginary) Self-Energy in the Migdal Approximation")')
WRITE(stdout, '(5x,a/)') REPEAT('=',67)
!
IF (fsthick < 1.d3) WRITE(stdout, '(/5x,a,f10.6,a)' ) 'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
WRITE(stdout, '(/5x,a,f10.6,a)') 'Golden Rule strictly enforced with T = ', eptemp * ryd2ev, ' eV'
@ -217,7 +218,7 @@
!
IF ((iqq == 1) .AND. .NOT. adapt_smearing) THEN
WRITE (stdout, 100) degaussw * ryd2ev, ngaussw
WRITE (stdout,'(a)') ' '
WRITE (stdout, '(a)') ' '
ENDIF
!
IF (restart) THEN
@ -240,92 +241,88 @@
!
fermicount = 0
DO ik = 1, nkf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! here we must have ef, not ef0, to be consistent with ephwann_shuffle
! (but in this case they are the same)
!
IF ((MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) .AND. &
(MINVAL(ABS(etf(:, ikq) - ef)) < fsthick)) THEN
!
fermicount = fermicount + 1
DO imode = 1, nmodes
DO ibnd = 1, nbndfst
!
! the energy of the electron at k (relative to Ef)
ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0
!
eta_tmp = eta2(ibnd, imode, ik)
inv_eta_tmp = inv_eta(ibnd, imode, ik)
!
DO jbnd = 1, nbndfst
!
! the fermi occupation for k+q
ekq = etf(ibndmin - 1 + jbnd, ikq) - ef0
wgkq = wgauss(-ekq * inv_eptemp, -99)
!
! here we take into account the zero-point SQRT(hbar/2M\omega)
! with hbar = 1 and M already contained in the eigenmodes
! g2 is Ry^2, wkf must already account for the spin factor
!
! SP: Shortrange is disabled for efficiency reasons
!IF (shortrange .AND. ( ABS(xqf (1, iq))> eps8 .OR. ABS(xqf (2, iq))> eps8 &
! .OR. ABS(xqf (3, iq))> eps8 )) THEN
! ! SP: The abs has to be removed. Indeed the epf17 can be a pure imaginary
! ! number, in which case its square will be a negative number.
! g2 = REAL((epf17 (jbnd, ibnd, imode, ik)**two) * inv_wq(imode) * g2_tmp(imode))
!ELSE
! g2 = (ABS(epf17 (jbnd, ibnd, imode, ik))**two) * inv_wq(imode) * g2_tmp(imode)
!ENDIF
!
g2 = (ABS(epf17(jbnd, ibnd, imode, ik))**two) * inv_wq(imode) * g2_tmp(imode)
!
! There is a sign error for wq in Eq. 9 of Comp. Phys. Comm. 181, 2140 (2010). - RM
! The sign was corrected according to Eq. (7.282) page 489 from Mahan's book
! (Many-Particle Physics, 3rd edition)
!
weight = wqf(iq) * REAL(((wgkq + wgq(imode)) / (ekk - (ekq - wq(imode)) - ci * eta_tmp) + &
(one - wgkq + wgq(imode)) / (ekk - (ekq + wq(imode)) - ci * eta_tmp)))
!
sigmar_all(ibnd, ik + lower_bnd - 1) = sigmar_all(ibnd, ik + lower_bnd - 1) + g2 * weight
!
! Logical implementation
! weight = wqf(iq) * aimag ( &
! ( ( wgkq + wgq ) / ( ekk - ( ekq - wq ) - ci * degaussw ) + &
! ( one - wgkq + wgq ) / ( ekk - ( ekq + wq ) - ci * degaussw ) ) )
!
! Delta implementation
w0g1 = w0gauss((ekk - ekq + wq(imode)) * inv_eta_tmp, 0) * inv_eta_tmp
w0g2 = w0gauss((ekk - ekq - wq(imode)) * inv_eta_tmp, 0) * inv_eta_tmp
weight = pi * wqf(iq) * ((wgkq + wgq(imode)) * w0g1 + (one - wgkq + wgq(imode)) * w0g2)
!
sigmai_all(ibnd, ik + lower_bnd - 1) = sigmai_all(ibnd, ik + lower_bnd - 1) + g2 * weight
!
! Mode-resolved
IF (iverbosity == 3) THEN
sigmai_mode(ibnd, imode, ik + lower_bnd - 1) = sigmai_mode(ibnd, imode, ik + lower_bnd - 1) + g2 * weight
ENDIF
!
! Z FACTOR: -\frac{\partial\Re\Sigma}{\partial\omega}
!
weight = wqf(iq) * &
(( wgkq + wgq(imode)) * ((ekk - (ekq - wq(imode)))**two - eta_tmp**two) / &
((ekk - (ekq - wq(imode)))**two + eta_tmp**two)**two + &
(one - wgkq + wgq(imode)) * ((ekk - (ekq + wq(imode)))**two - eta_tmp**two) / &
((ekk - (ekq + wq(imode)))**two + eta_tmp**two)**two )
!
zi_all(ibnd, ik + lower_bnd - 1) = zi_all(ibnd, ik + lower_bnd - 1) + g2 * weight
!
ENDDO !jbnd
!
ENDDO !ibnd
!
ENDDO !imode
!
ENDIF ! endif fsthick
!
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! here we must have ef, not ef0, to be consistent with ephwann_shuffle
! (but in this case they are the same)
!
IF ((MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) .AND. &
(MINVAL(ABS(etf(:, ikq) - ef)) < fsthick)) THEN
!
fermicount = fermicount + 1
DO imode = 1, nmodes
DO ibnd = 1, nbndfst
!
! the energy of the electron at k (relative to Ef)
ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0
!
eta_tmp = eta2(ibnd, imode, ik)
inv_eta_tmp = inv_eta(ibnd, imode, ik)
!
DO jbnd = 1, nbndfst
!
! the fermi occupation for k+q
ekq = etf(ibndmin - 1 + jbnd, ikq) - ef0
wgkq = wgauss(-ekq * inv_eptemp, -99)
!
! here we take into account the zero-point SQRT(hbar/2M\omega)
! with hbar = 1 and M already contained in the eigenmodes
! g2 is Ry^2, wkf must already account for the spin factor
!
! SP: Shortrange is disabled for efficiency reasons
!IF (shortrange .AND. ( ABS(xqf (1, iq))> eps8 .OR. ABS(xqf (2, iq))> eps8 &
! .OR. ABS(xqf (3, iq))> eps8 )) THEN
! ! SP: The abs has to be removed. Indeed the epf17 can be a pure imaginary
! ! number, in which case its square will be a negative number.
! g2 = REAL((epf17 (jbnd, ibnd, imode, ik)**two) * inv_wq(imode) * g2_tmp(imode))
!ELSE
! g2 = (ABS(epf17 (jbnd, ibnd, imode, ik))**two) * inv_wq(imode) * g2_tmp(imode)
!ENDIF
!
g2 = (ABS(epf17(jbnd, ibnd, imode, ik))**two) * inv_wq(imode) * g2_tmp(imode)
!
! There is a sign error for wq in Eq. 9 of Comp. Phys. Comm. 181, 2140 (2010). - RM
! The sign was corrected according to Eq. (7.282) page 489 from Mahan's book
! (Many-Particle Physics, 3rd edition)
!
weight = wqf(iq) * REAL(((wgkq + wgq(imode)) / (ekk - (ekq - wq(imode)) - ci * eta_tmp) + &
(one - wgkq + wgq(imode)) / (ekk - (ekq + wq(imode)) - ci * eta_tmp)))
!
sigmar_all(ibnd, ik + lower_bnd - 1) = sigmar_all(ibnd, ik + lower_bnd - 1) + g2 * weight
!
! Logical implementation
! weight = wqf(iq) * aimag ( &
! ( ( wgkq + wgq ) / ( ekk - ( ekq - wq ) - ci * degaussw ) + &
! ( one - wgkq + wgq ) / ( ekk - ( ekq + wq ) - ci * degaussw ) ) )
!
! Delta implementation
w0g1 = w0gauss((ekk - ekq + wq(imode)) * inv_eta_tmp, 0) * inv_eta_tmp
w0g2 = w0gauss((ekk - ekq - wq(imode)) * inv_eta_tmp, 0) * inv_eta_tmp
weight = pi * wqf(iq) * ((wgkq + wgq(imode)) * w0g1 + (one - wgkq + wgq(imode)) * w0g2)
!
sigmai_all(ibnd, ik + lower_bnd - 1) = sigmai_all(ibnd, ik + lower_bnd - 1) + g2 * weight
!
! Mode-resolved
IF (iverbosity == 3) THEN
sigmai_mode(ibnd, imode, ik + lower_bnd - 1) = sigmai_mode(ibnd, imode, ik + lower_bnd - 1) + g2 * weight
ENDIF
!
! Z FACTOR: -\frac{\partial\Re\Sigma}{\partial\omega}
!
weight = wqf(iq) * &
(( wgkq + wgq(imode)) * ((ekk - (ekq - wq(imode)))**two - eta_tmp**two) / &
((ekk - (ekq - wq(imode)))**two + eta_tmp**two)**two + &
(one - wgkq + wgq(imode)) * ((ekk - (ekq + wq(imode)))**two - eta_tmp**two) / &
((ekk - (ekq + wq(imode)))**two + eta_tmp**two)**two )
!
zi_all(ibnd, ik + lower_bnd - 1) = zi_all(ibnd, ik + lower_bnd - 1) + g2 * weight
!
ENDDO !jbnd
ENDDO !ibnd
ENDDO !imode
ENDIF ! endif fsthick
ENDDO ! end loop on k
!
! Creation of a restart point
@ -346,8 +343,10 @@
!
IF (iqq == totq) THEN
!
ALLOCATE(xkf_all(3, nkqtotf))
ALLOCATE(etf_all(nbndsub, nkqtotf))
ALLOCATE(xkf_all(3, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error allocating xkf_all', 1)
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error allocating etf_all', 1)
xkf_all(:, :) = zero
etf_all(:, :) = zero
!
@ -440,21 +439,21 @@
ryd2mev * sigmai_all(ibnd,ik), zi_all(ibnd, ik), one / zi_all(ibnd, ik) - one
IF (iverbosity == 3) THEN
DO imode = 1, nmodes
WRITE(linewidth_elself,'(i9,2x)', ADVANCE = 'no') ik
WRITE(linewidth_elself,'(i9,2x)', ADVANCE = 'no') ibndmin - 1 + ibnd
WRITE(linewidth_elself,'(E22.14,2x)', ADVANCE = 'no') ryd2ev * ekk
WRITE(linewidth_elself,'(i9,2x)', ADVANCE = 'no') imode
WRITE(linewidth_elself,'(E22.14,2x)') ryd2mev * sigmai_mode(ibnd, imode, ik)
WRITE(linewidth_elself, '(i9,2x)', ADVANCE = 'no') ik
WRITE(linewidth_elself, '(i9,2x)', ADVANCE = 'no') ibndmin - 1 + ibnd
WRITE(linewidth_elself, '(E22.14,2x)', ADVANCE = 'no') ryd2ev * ekk
WRITE(linewidth_elself, '(i9,2x)', ADVANCE = 'no') imode
WRITE(linewidth_elself, '(E22.14,2x)') ryd2mev * sigmai_mode(ibnd, imode, ik)
ENDDO
ELSE
WRITE(linewidth_elself,'(i9,2x)', ADVANCE = 'no') ik
WRITE(linewidth_elself,'(i9,2x)', ADVANCE = 'no') ibndmin - 1 + ibnd
WRITE(linewidth_elself,'(E22.14,2x)', ADVANCE = 'no') ryd2ev * ekk
WRITE(linewidth_elself,'(E22.14,2x)') ryd2mev * sigmai_all(ibnd, ik)
WRITE(linewidth_elself, '(i9,2x)', ADVANCE = 'no') ik
WRITE(linewidth_elself, '(i9,2x)', ADVANCE = 'no') ibndmin - 1 + ibnd
WRITE(linewidth_elself, '(E22.14,2x)', ADVANCE = 'no') ryd2ev * ekk
WRITE(linewidth_elself, '(E22.14,2x)') ryd2mev * sigmai_all(ibnd, ik)
ENDIF
!
ENDDO
WRITE(stdout,'(5x,a/)') REPEAT('-',67)
WRITE(stdout, '(5x,a/)') REPEAT('-',67)
!
ENDDO
ENDIF
@ -471,7 +470,7 @@
! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} )
!zi_all (ibnd,ik) = one / ( one + zi_all (ibnd,ik) )
!
WRITE(stdout,'(2i9,5f12.4)') ik, ibndmin - 1 + ibnd, ryd2ev * ekk, ryd2mev * sigmar_all(ibnd, ik), &
WRITE(stdout, '(2i9,5f12.4)') ik, ibndmin - 1 + ibnd, ryd2ev * ekk, ryd2mev * sigmar_all(ibnd, ik), &
ryd2mev * sigmai_all(ibnd, ik), zi_all(ibnd, ik), one / zi_all(ibnd, ik) - one
!
ENDDO
@ -481,8 +480,10 @@
ENDDO
!
CLOSE(linewidth_elself)
DEALLOCATE(xkf_all)
DEALLOCATE(etf_all)
DEALLOCATE(xkf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error deallocating xkf_all', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error deallocating etf_all', 1)
!
ENDIF
!
@ -830,8 +831,8 @@
! collect contributions from all pools (sum over k-points)
! this finishes the integral over the BZ (k)
!
CALL mp_sum(gamma,inter_pool_comm)
CALL mp_sum(gamma_v,inter_pool_comm)
CALL mp_sum(gamma, inter_pool_comm)
CALL mp_sum(gamma_v, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
!
@ -965,6 +966,8 @@
!!
INTEGER :: n
!!
INTEGER :: ierr
!! Error status
!! Integer for the degenerate average over eigenstates
REAL(KIND = DP) :: tmp
!! Temporary variable to store real part of Sigma for the degenerate average
@ -1213,8 +1216,10 @@
!
IF (iqq == totq) THEN
!
ALLOCATE(xkf_all(3, nkqtotf))
ALLOCATE(etf_all(nbndsub, nkqtotf))
ALLOCATE(xkf_all(3, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_pl_q', 'Error allocating xkf_all', 1)
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_pl_q', 'Error allocating etf_all', 1)
xkf_all(:, :) = zero
etf_all(:, :) = zero
!
@ -1308,8 +1313,10 @@
!
CLOSE(linewidth_elself)
!
DEALLOCATE(xkf_all)
DEALLOCATE(etf_all)
DEALLOCATE(xkf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_pl_q', 'Error deallocating xkf_all', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_pl_q', 'Error deallocating etf_all', 1)
!
ENDIF
!
@ -1417,4 +1424,3 @@
!-----------------------------------------------------------------------
END MODULE selfen
!-----------------------------------------------------------------------

30
EPW/src/setphases_wrap.f90
View File

@ -7,7 +7,8 @@
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!
SUBROUTINE setphases_wrap
!---------------------------------------------------------------------
SUBROUTINE setphases_wrap()
!---------------------------------------------------------------------
!!
!! This is the wrapper which is used to set the phases of the wavefunctions
@ -17,7 +18,6 @@
!!
!! The phases for all wavefunctions are now stored in umat_all
!!
!!
!---------------------------------------------------------------------
USE kinds, ONLY : DP
use klist, only : nkstot
@ -35,19 +35,21 @@
!! Band-index
INTEGER :: jbnd
!! Band index
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: zero_vect(3)
!! Real vector
!
IF (nproc_pool>1) call errore &
('setphases_wrap', 'only one proc per pool', 1)
IF (nproc_pool>1) CALL errore('setphases_wrap', 'only one proc per pool', 1)
!
ALLOCATE(umat_all(nbnd, nbnd, nkstot))
ALLOCATE(umat(nbnd, nbnd, nks))
ALLOCATE(umat_all(nbnd, nbnd, nkstot), STAT = ierr)
IF (ierr /= 0) CALL errore('setphases_wrap', 'Error allocating umat_all', 1)
ALLOCATE(umat(nbnd, nbnd, nks), STAT = ierr)
IF (ierr /= 0) CALL errore('setphases_wrap', 'Error allocating umat', 1)
umat_all = (0.d0, 0.d0)
zero_vect = 0.d0
!
!
WRITE(stdout,'(5x,a)') 'No wavefunction gauge setting applied'
WRITE(stdout, '(5x,a)') 'No wavefunction gauge setting applied'
!
IF (ionode) THEN
DO ik = 1, nkstot
@ -78,14 +80,6 @@
!
CALL mp_sum(umat_all, inter_pool_comm)
!
!IF (iverbosity == 1) then
! WRITE (stdout,* ) "Phase setting matrices:"
! DO ik = 1, nkstot
! DO ibnd = 1, nbnd
! WRITE(stdout, '(8f8.5)') umat_all(:,ibnd, ik)
! ENDDO
! WRITE(stdout,*)
! ENDDO
!ENDIF
!
!---------------------------------------------------------------------
END SUBROUTINE setphases_wrap
!---------------------------------------------------------------------

47
EPW/src/spectral_func.f90
View File

@ -62,7 +62,6 @@
!! Total number of q-point in window
!
! Local variables
!
INTEGER :: iw
!! Counter on the frequency
INTEGER :: ik
@ -83,6 +82,8 @@
!! Lower bounds index after k or q paral
INTEGER :: upper_bnd
!! Upper bounds index after k or q paral
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: g2
!! Electron-phonon matrix elements squared in Ry^2
REAL(KIND = DP) :: ekk
@ -264,8 +265,10 @@
!
IF (iqq == totq) THEN
!
ALLOCATE(xkf_all(3, nkqtotf))
ALLOCATE(etf_all(nbndsub, nkqtotf))
ALLOCATE(xkf_all(3, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_func_q', 'Error allocating xkf_all', 1)
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_func_q', 'Error allocating etf_all', 1)
xkf_all(:, :) = zero
etf_all(:, :) = zero
!
@ -392,8 +395,10 @@
!
IF (me_pool == 0) CLOSE(iospectral_sup)
!
DEALLOCATE(xkf_all)
DEALLOCATE(etf_all)
DEALLOCATE(xkf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_func_q', 'Error deallocating xkf_all', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_func_q', 'Error deallocating etf_all', 1)
!
ENDIF
!
@ -448,6 +453,7 @@
INTEGER, INTENT(in) :: totq
!! Total q-points in selecq window
!
! Local variables
INTEGER :: ik
!! Counter on the k-point index
INTEGER :: ikk
@ -794,6 +800,8 @@
!!
INTEGER :: upper_bnd
!!
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: g2
!!
REAL(KIND = DP) :: ekk
@ -869,14 +877,14 @@
dw = (wmax_specfun - wmin_specfun) / DBLE(nw_specfun - 1)
!
IF (iqq == 1) THEN
!
WRITE(stdout, '(/5x,a)') REPEAT('=',67)
WRITE(stdout, '(5x,"Electron Spectral Function in the Migdal Approximation")')
WRITE(stdout, '(5x,a/)') REPEAT('=',67)
!
IF (fsthick < 1.d3 ) WRITE(stdout, '(/5x,a,f10.6,a)' ) 'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
WRITE(stdout, '(/5x,a,f10.6,a)' ) 'Golden Rule strictly enforced with T = ', eptemp * ryd2ev, ' eV'
!
!
WRITE(stdout, '(/5x,a)') REPEAT('=',67)
WRITE(stdout, '(5x,"Electron Spectral Function in the Migdal Approximation")')
WRITE(stdout, '(5x,a/)') REPEAT('=',67)
!
IF (fsthick < 1.d3 ) WRITE(stdout, '(/5x,a,f10.6,a)' ) 'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
WRITE(stdout, '(/5x,a,f10.6,a)' ) 'Golden Rule strictly enforced with T = ', eptemp * ryd2ev, ' eV'
!
ENDIF
!
! Fermi level and corresponding DOS
@ -1026,8 +1034,10 @@
!
IF (iqq == totq) THEN
!
ALLOCATE(xkf_all(3, nkqtotf))
ALLOCATE(etf_all(nbndsub, nkqtotf))
ALLOCATE(xkf_all(3, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_func_pl_q', 'Error allocating xkf_all', 1)
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_func_pl_q', 'Error allocating etf_all', 1)
xkf_all(:, :) = zero
etf_all(:, :) = zero
!
@ -1150,8 +1160,10 @@
!
IF (me_pool == 0) CLOSE(iospectral_sup)
!
DEALLOCATE(xkf_all)
DEALLOCATE(etf_all)
DEALLOCATE(xkf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_func_pl_q', 'Error deallocating xkf_all', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_func_pl_q', 'Error deallocating etf_all', 1)
ENDIF
!
100 FORMAT(5x, 'Gaussian Broadening: ', f10.6, ' eV, ngauss=', i4)
@ -1166,4 +1178,3 @@
!-----------------------------------------------------------------------
END MODULE spectral_func
!-----------------------------------------------------------------------

8
EPW/src/star_q2.f90
View File

@ -71,7 +71,7 @@
!! a zero vector: used in eqvect
REAL(KIND = DP), PARAMETER :: accep = 1.e-5_dp
!! Tolerence
!
zero(:) = 0.d0
!
! go to crystal coordinates
@ -112,9 +112,9 @@
saq(:, nq) = raq(:)
sym_smallq(nq) = isym
DO i = 1, 3
sxq(i, nq) = bg(i, 1) * saq(1, nq) &
+ bg(i, 2) * saq(2, nq) &
+ bg(i, 3) * saq(3, nq)
sxq(i, nq) = bg(i, 1) * saq(1, nq) &
+ bg(i, 2) * saq(2, nq) &
+ bg(i, 3) * saq(3, nq)
ENDDO
ENDIF
ENDDO

19
EPW/src/stop_epw.f90
View File

@ -11,7 +11,7 @@
! Modified from stop_ph
!
!-----------------------------------------------------------------------
SUBROUTINE stop_epw
SUBROUTINE stop_epw()
!-----------------------------------------------------------------------
!!
!! Close all files and synchronize processes before stopping.
@ -26,7 +26,7 @@
!
IMPLICIT NONE
!
CALL print_clock_epw
CALL print_clock_epw()
!
WRITE(stdout, '(a)') " "
WRITE(stdout, '(a)') " Please consider citing: "
@ -40,31 +40,30 @@
!
! Eliashberg superconductivity
IF (eliashberg) THEN
WRITE(stdout, '(a)') " eliashberg :: E. R. Margine and F. Giustino, Phys. Rev. B 87, 024505 (2013) "
WRITE(stdout, '(a)') " eliashberg :: E. R. Margine and F. Giustino, Phys. Rev. B 87, 024505 (2013) "
ENDIF
!
! Plasmons
IF (plselfen .OR. specfun_pl) THEN
WRITE(stdout, '(a)') " plselfen or specfun_pl :: F. Caruso, C. Verdi, S. Ponce and F. Giustino, Phys. Rev. B 97, 165113 (2018) "
WRITE(stdout, '(a)') &
" plselfen or specfun_pl :: F. Caruso, C. Verdi, S. Ponce and F. Giustino, Phys. Rev. B 97, 165113 (2018) "
ENDIF
!
! Transport module
IF (scattering) THEN
WRITE(stdout, '(a)') " scattering :: S. Ponce, E. R. Margine and F. Giustino, Phys. Rev. B 97, 121201 (2018) "
WRITE(stdout, '(a)') " scattering :: S. Ponce, E. R. Margine and F. Giustino, Phys. Rev. B 97, 121201 (2018) "
ENDIF
IF (iterative_bte) THEN
WRITE(stdout, '(a)') " iterative_bte :: S. Ponce, E. R. Margine and F. Giustino, Phys. Rev. B 97, 121201 (2018) "
WRITE(stdout, '(a)') " iterative_bte :: F. Macheda and N. Bonini, Phys. Rev. B 98, 201201 (2018) "
WRITE(stdout, '(a)') " iterative_bte :: S. Ponce, E. R. Margine and F. Giustino, Phys. Rev. B 97, 121201 (2018) "
WRITE(stdout, '(a)') " iterative_bte :: F. Macheda and N. Bonini, Phys. Rev. B 98, 201201 (2018) "
ENDIF
!
! Improvements
IF (adapt_smearing) THEN
WRITE(stdout, '(a)') " adapt_smearing :: F. Macheda and N. Bonini, Phys. Rev. B 98, 201201 (2018) "
WRITE(stdout, '(a)') " adapt_smearing :: F. Macheda and N. Bonini, Phys. Rev. B 98, 201201 (2018) "
ENDIF
!
CALL mp_barrier(inter_pool_comm)
CALL mp_end(inter_pool_comm)
!
CALL mp_global_end()
!
STOP

91
EPW/src/transport.f90
View File

@ -20,21 +20,21 @@
SUBROUTINE scattering_rate_q(iqq, iq, totq, ef0, efcb, first_cycle)
!-----------------------------------------------------------------------
!!
!! This routine computes the scattering rate (inv_tau)
!! This routine computes the scattering rate (inv_tau)
!!
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE phcom, ONLY : nmodes
USE epwcom, ONLY : nbndsub, fsthick, eps_acustic, degaussw, &
nstemp, scattering_serta, scattering_0rta, shortrange,&
restart, restart_freq, restart_filq, vme
USE epwcom, ONLY : nbndsub, fsthick, eps_acustic, degaussw, restart, &
nstemp, scattering_serta, scattering_0rta, shortrange, &
restart_freq, restart_filq, vme
USE pwcom, ONLY : ef
USE elph2, ONLY : ibndmax, ibndmin, etf, nkqf, nkf, dmef, vmef, wf, wqf, &
epf17, nqtotf, nkqtotf, inv_tau_all, inv_tau_allcb, &
epf17, nqtotf, nkqtotf, inv_tau_all, inv_tau_allcb, &
xqf, zi_allvb, zi_allcb, nbndfst, nktotf, transp_temp, &
lower_bnd
USE constants_epw, ONLY : zero, one, two, pi, ryd2mev, kelvin2eV, ryd2ev, &
USE constants_epw, ONLY : zero, one, two, pi, ryd2mev, kelvin2eV, ryd2ev, &
eps6, eps8, eps4
USE mp, ONLY : mp_barrier, mp_sum
USE mp_global, ONLY : world_comm
@ -75,7 +75,8 @@
!! Index over temperature range
INTEGER :: nqtotf_new
!! Number of q-point in the new dataset
!
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: tmp
!! Temporary variable to store real part of Sigma for the degenerate average
REAL(KIND = DP) :: tmp2
@ -125,7 +126,6 @@
!! Temporary array to store the zi
REAL(KIND = DP), ALLOCATABLE :: inv_tau_all_new (:, :, :)
!! New scattering rates to be merged
!
REAL(KIND = DP), ALLOCATABLE :: etf_all(:, :)
!! Eigen-energies on the fine grid collected from all pools in parallel case
REAL(KIND = DP), EXTERNAL :: DDOT
@ -145,7 +145,7 @@
WRITE(stdout, '(5x,"Scattering rate")')
WRITE(stdout, '(5x,a/)') REPEAT('=',67)
!
IF (fsthick < 1.d3 ) &
IF (fsthick < 1.d3) &
WRITE(stdout, '(/5x,a,f10.6,a)' ) 'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
WRITE(stdout, '(5x,a,f10.6,a)' ) 'This is computed with respect to the fine Fermi level ',ef * ryd2ev, ' eV'
WRITE(stdout, '(5x,a,f10.6,a,f10.6,a)' ) 'Only states between ',(ef - fsthick) * ryd2ev, ' eV and ', &
@ -189,7 +189,7 @@
!
IF (ABS(vkk(1, ibnd)**2 + vkk(2, ibnd)**2 + vkk(3, ibnd)**2 ) > eps4) &
vel_factor(ibnd, jbnd) = DDOT(3, vkk(:, ibnd), 1, vkq(:, jbnd), 1) / &
DDOT(3, vkk(:, ibnd), 1, vkk(:, ibnd), 1)
DDOT(3, vkk(:, ibnd), 1, vkk(:, ibnd), 1)
ENDDO
ENDDO
ELSE
@ -205,7 +205,7 @@
!
IF (ABS(vkk(1, ibnd)**2 + vkk(2, ibnd)**2 + vkk(3, ibnd)**2 ) > eps4) &
vel_factor(ibnd, jbnd) = DDOT(3, vkk(:, ibnd), 1, vkq(:, jbnd), 1) / &
DDOT(3, vkk(:, ibnd), 1, vkk(:, ibnd), 1)
DDOT(3, vkk(:, ibnd), 1, vkk(:, ibnd), 1)
ENDDO
ENDDO
ENDIF
@ -272,8 +272,7 @@
! [1 - f(E_k+q) + n(w_q)] * delta[E_k - E_k+q - w_q] }
!
! DBSP Just to try
trans_prob = pi * wqf(iq) * g2 * &
((fmkq + wgq) * w0g1 + (one - fmkq + wgq) * w0g2)
trans_prob = pi * wqf(iq) * g2 * ((fmkq + wgq) * w0g1 + (one - fmkq + wgq) * w0g2)
!
IF (scattering_serta) THEN
! energy relaxation time approximation
@ -403,7 +402,8 @@
!
! The total number of k points
!
ALLOCATE(etf_all(nbndsub, nkqtotf))
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('scattering_rate_q', 'Error allocating etf_all', 1)
!
CALL mp_sum(inv_tau_all, world_comm)
IF (ABS(efcb(1)) > eps4) CALL mp_sum(inv_tau_allcb, world_comm)
@ -429,14 +429,16 @@
! In case we read another q-file, merge the scattering here
IF (restart_filq /= '') THEN
!
ALLOCATE(inv_tau_all_new(nstemp, nbndfst, nktotf))
ALLOCATE(inv_tau_all_new(nstemp, nbndfst, nktotf), STAT = ierr)
IF (ierr /= 0) CALL errore('scattering_rate_q', 'Error allocating inv_tau_all_new', 1)
inv_tau_all_new(:, :, :) = zero
!
CALL merge_read(nktotf, nqtotf_new, inv_tau_all_new)
!
inv_tau_all(:, :, :) = (inv_tau_all(:, :, :) * totq &
+ inv_tau_all_new(:, :, :) * nqtotf_new) / (totq + nqtotf_new)
DEALLOCATE(inv_tau_all_new)
DEALLOCATE(inv_tau_all_new, STAT = ierr)
IF (ierr /= 0) CALL errore('scattering_rate_q', 'Error deallocating inv_tau_all_new', 1)
!
WRITE(stdout, '(a)' ) ' '
WRITE(stdout, '(a,i10,a)' ) ' Merge scattering for a total of ', totq + nqtotf_new, ' q-points'
@ -515,7 +517,8 @@
!
ENDDO !nstemp
!
DEALLOCATE(etf_all)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('scattering_rate_q', 'Error deallocating etf_all', 1)
!
! Creation of a restart point at the end
IF (restart) THEN
@ -550,7 +553,7 @@
SUBROUTINE transport_coeffs(ef0, efcb)
!-----------------------------------------------------------------------
!!
!! This SUBROUTINE computes the transport coefficients
!! This routine computes the transport coefficients
!! SP - June 2018 - Update for symmetries in velocities when using homogeneous grids.
!! This is currently commented out since it is ONLY needed if we are
!! interested in the off-diagonal mobility_\alpha\beta terms.
@ -562,20 +565,17 @@
USE cell_base, ONLY : alat, at, omega
USE io_files, ONLY : prefix
USE io_epw, ONLY : iufilsigma
USE epwcom, ONLY : nbndsub, fsthick, &
system_2d, nstemp, &
int_mob, ncarrier, scatread, &
iterative_bte, vme
USE epwcom, ONLY : nbndsub, fsthick, system_2d, nstemp, &
int_mob, ncarrier, scatread, iterative_bte, vme
USE pwcom, ONLY : ef
USE elph2, ONLY : ibndmax, ibndmin, etf, nkf, wkf, dmef, vmef, &
USE elph2, ONLY : ibndmax, ibndmin, etf, nkf, wkf, dmef, vmef, &
inv_tau_all, nkqtotf, inv_tau_allcb, transp_temp, &
zi_allvb, zi_allcb, map_rebal, nbndfst, nktotf
USE constants_epw, ONLY : zero, one, bohr2ang, ryd2ev, electron_SI, &
USE constants_epw, ONLY : zero, one, bohr2ang, ryd2ev, electron_SI, &
kelvin2eV, hbar, Ang2m, hbarJ, ang2cm, czero
USE mp, ONLY : mp_sum, mp_bcast
USE mp_global, ONLY : world_comm
USE mp_world, ONLY : mpime
! SP - Uncomment to use symmetries on velocities
USE symm_base, ONLY : s, t_rev, time_reversal, set_sym_bl, nrot
USE io_global, ONLY : ionode_id
USE cell_base, ONLY : bg
@ -622,13 +622,14 @@
!! k-point index that run on the full BZ
INTEGER :: nb
!! Number of points in the BZ corresponding to a point in IBZ
INTEGER :: ierr
!! Error status
INTEGER :: BZtoIBZ_tmp(nkf1 * nkf2 * nkf3)
!! Temporary mapping
INTEGER :: BZtoIBZ(nkf1 * nkf2 * nkf3)
!! BZ to IBZ mapping
INTEGER(SIK2) :: s_BZtoIBZ(nkf1 * nkf2 * nkf3)
!! symmetry
!
REAL(KIND = DP) :: ekk
!! Energy relative to Fermi level: $$\varepsilon_{n\mathbf{k}}-\varepsilon_F$$
REAL(KIND = DP) :: dfnk
@ -712,7 +713,6 @@
inv_cell = 1.0d0 / omega
! for 2d system need to divide by area (vacuum in z-direction)
IF (system_2d ) inv_cell = inv_cell * at(3, 3) * alat
!
! We can read the scattering rate from files.
IF (scatread) THEN
@ -727,32 +727,39 @@
! Lets gather the velocities from all pools
#if defined(__MPI)
IF (vme) THEN
ALLOCATE(vmef_all(3, nbndsub, nbndsub, nkqtotf))
ALLOCATE(vmef_all(3, nbndsub, nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error allocating vmef_all', 1)
vmef_all(:, :, :, :) = czero
CALL poolgatherc4(3, nbndsub, nbndsub, nkqtotf, 2 * nkf, vmef, vmef_all)
ELSE
ALLOCATE(dmef_all(3, nbndsub, nbndsub, nkqtotf))
ALLOCATE(dmef_all(3, nbndsub, nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error allocating dmef_all', 1)
dmef_all(:, :, :, :) = czero
CALL poolgatherc4(3, nbndsub, nbndsub, nkqtotf, 2 * nkf, dmef, dmef_all)
ENDIF
ALLOCATE(wkf_all(nkqtotf))
ALLOCATE(wkf_all(nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error allocating wkf_all', 1)
wkf_all(:) = zero
CALL poolgather2(1, nkqtotf, 2 * nkf, wkf, wkf_all)
#else
IF (vme) THEN
ALLOCATE(vmef_all(3, nbndsub, nbndsub, nkqtotf))
ALLOCATE(vmef_all(3, nbndsub, nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error allocating vmef_all', 1)
vmef_all = vmef
ELSE
ALLOCATE(dmef_all(3, nbndsub, nbndsub, nkqtotf))
ALLOCATE(dmef_all(3, nbndsub, nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error allocating dmef_all', 1)
dmef_all = dmef
ENDIF
#endif
ALLOCATE(tdf_sigma_m(3, 3, nbndfst, nkqtotf))
ALLOCATE(tdf_sigma_m(3, 3, nbndfst, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error allocating tdf_sigma_m', 1)
tdf_sigma_m(:, :, :, :) = zero
!
! In this case, the sum over q has already been done. It should therefore be ok
! to do the mobility in sequential. Each cpu does the same thing below
ALLOCATE(etf_all(nbndsub, nktotf))
ALLOCATE(etf_all(nbndsub, nktotf), STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error allocating etf_all', 1)
!
CALL scattering_read(etemp, ef0(itemp), etf_all, inv_tau_all)
!
@ -905,12 +912,16 @@
ENDIF ! int_mob .OR. (ncarrier > 1E5)
!
IF (vme) THEN
DEALLOCATE(vmef_all)
DEALLOCATE(vmef_all, STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error deallocating vmef_all', 1)
ELSE
DEALLOCATE(dmef_all)
DEALLOCATE(dmef_all, STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error deallocating dmef_all', 1)
ENDIF
DEALLOCATE(tdf_sigma_m)
DEALLOCATE(etf_all)
DEALLOCATE(tdf_sigma_m, STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error deallocating tdf_sigma_m', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('transport_coeffs', 'Error deallocating etf_all', 1)
ENDDO ! itemp
!
ELSE ! Case without reading the scattering rates from files.
@ -1161,7 +1172,6 @@
! WRITE(stdout,'(45x, 1E18.6, a)') mobility_yy, ' y-axis [Z]'
! WRITE(stdout,'(45x, 1E18.6, a)') mobility_zz, ' z-axis [Z]'
! WRITE(stdout,'(45x, 1E18.6, a)') mobility, ' avg [Z]'
!
ENDDO ! nstemp
!
@ -1464,7 +1474,6 @@
! WRITE(stdout,'(45x, 1E18.6, a)') mobility_yy, ' y-axis [Z]'
! WRITE(stdout,'(45x, 1E18.6, a)') mobility_zz, ' z-axis [Z]'
! WRITE(stdout,'(45x, 1E18.6, a)') mobility, ' avg [Z]'
!
ENDDO ! nstemp
WRITE(stdout,'(5x)')

125
EPW/src/transport_iter.f90
View File

@ -26,7 +26,6 @@
!! The fine k-point and q-point grid have to be commensurate.
!! The k-point grid uses crystal symmetry to decrease computational cost.
!!
!-----------------------------------------------------------------------
USE kinds, ONLY : DP, sgl
USE io_global, ONLY : stdout
USE cell_base, ONLY : alat, at, omega, bg
@ -35,12 +34,13 @@
system_2d, int_mob, ncarrier, restart, restart_freq, &
mp_mesh_k, nkf1, nkf2, nkf3, vme, broyden_beta
USE pwcom, ONLY : ef
USE elph2, ONLY : ibndmax, ibndmin, etf, nkqf, nkf, wkf, dmef, vmef, &
wf, xkf, epf17, nqtotf, nkqtotf, nbndfst, nktotf, &
map_rebal, xqf, wqf, nqf, transp_temp, mobilityh_save, &
mobilityel_save, lower_bnd, ixkqf_tr, s_BZtoIBZ_full
USE constants_epw, ONLY : zero, one, two, pi, kelvin2eV, ryd2ev, &
electron_SI, bohr2ang, ang2cm, hbarJ, eps6, eps8, eps10, &
USE elph2, ONLY : ibndmax, ibndmin, etf, nkqf, nkf, wkf, dmef, vmef, &
wf, xkf, epf17, nqtotf, nkqtotf, nbndfst, nktotf, &
map_rebal, xqf, wqf, nqf, transp_temp, &
mobilityel_save, lower_bnd, ixkqf_tr, s_BZtoIBZ_full,&
mobilityh_save
USE constants_epw, ONLY : zero, one, two, pi, kelvin2eV, ryd2ev, eps10, &
electron_SI, bohr2ang, ang2cm, hbarJ, eps6, eps8, &
eps2, eps4, eps20, eps80, eps160, hbar, cm2m, byte2Mb
USE mp, ONLY : mp_barrier, mp_sum, mp_bcast
USE mp_global, ONLY : inter_pool_comm, world_comm, my_pool_id
@ -94,8 +94,6 @@
!! Is the current k-point a special k-point
LOGICAL :: special_map(nind)
!! Special k-point map
INTEGER :: ierr
!! Error flag
INTEGER :: ind
!! Index for sparse matrix
INTEGER :: iter
@ -182,6 +180,8 @@
!! Symmetry mapping
INTEGER :: max_size
!! Max size of the arrays
INTEGER :: ierr
!! Error status
INTEGER, ALLOCATABLE :: sp_map(:, :)
!! Mapping for special points
INTEGER, ALLOCATABLE :: counter_map(:)
@ -194,7 +194,6 @@
!! Inverse of k_inside_fsthick_id_nrot_inv
INTEGER, ALLOCATABLE :: xkf_sp(:, :)
!! Special k-points
!
REAL(KIND = DP) :: tau
!! Relaxation time
REAL(KIND = DP) :: ekk
@ -330,21 +329,23 @@
!
IF (iverbosity == 4) THEN
! Array size reporting
WRITE(stdout,'(5x,a)') 'Big array size reporting [Mb]'
WRITE(stdout,'(5x,a)') '-- ibte --'
WRITE(stdout,'(5x,a,f12.6)') 'BZtoIBZ : ', nkf1 * nkf2 * nkf3 * byte2Mb / 2 !INTEGER(4)
WRITE(stdout,'(5x,a,f12.6)') 's_BZtoIBZ : ', nkf1 * nkf2 * nkf3 * byte2Mb / 2 / 4 !INTEGER(SIK2)
WRITE(stdout,'(5x,a,f12.6)') 'F_SERTA : ', (3 * (nbndfst) * (nktotf) * nstemp) * byte2Mb!REAL(8)
WRITE(stdout,'(5x,a,f12.6)') 'BZtoIBZ_mat : ', nrot * (nktotf) * byte2Mb / 2 !INTEGER(4)
WRITE(stdout,'(5x,a,f12.6)') 'xkf_bz : ', 3 * nkf1 * nkf2 * nkf3 * byte2Mb / 2 !REAL(4)
WRITE(stdout, '(5x,a)') 'Big array size reporting [Mb]'
WRITE(stdout, '(5x,a)') '-- ibte --'
WRITE(stdout, '(5x,a,f12.6)') 'BZtoIBZ : ', nkf1 * nkf2 * nkf3 * byte2Mb / 2 !INTEGER(4)
WRITE(stdout, '(5x,a,f12.6)') 's_BZtoIBZ : ', nkf1 * nkf2 * nkf3 * byte2Mb / 2 / 4 !INTEGER(SIK2)
WRITE(stdout, '(5x,a,f12.6)') 'F_SERTA : ', (3 * (nbndfst) * (nktotf) * nstemp) * byte2Mb!REAL(8)
WRITE(stdout, '(5x,a,f12.6)') 'BZtoIBZ_mat : ', nrot * (nktotf) * byte2Mb / 2 !INTEGER(4)
WRITE(stdout, '(5x,a,f12.6)') 'xkf_bz : ', 3 * nkf1 * nkf2 * nkf3 * byte2Mb / 2 !REAL(4)
rws(:, :) = zero
CALL wsinit(rws, nrwsx, nrws, bg)
ENDIF
!
! Deal with symmetries
IF (mp_mesh_k) THEN
ALLOCATE(ixkqf_tr(nind))
ALLOCATE(s_BZtoIBZ_full(nind))
ALLOCATE(ixkqf_tr(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('ibte', 'Error allocating ixkqf_tr', 1)
ALLOCATE(s_BZtoIBZ_full(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('ibte', 'Error allocating s_BZtoIBZ_full', 1)
! For a given k-point in the IBZ gives the k-point index
! of all the k-point in the full BZ that are connected to the current
! one by symmetry. nrot is the max number of symmetry
@ -411,7 +412,7 @@
IF (SUM(ABS(F_SERTA(:, ibnd, ik, itemp))) > eps160) THEN
xkk = xkf_all(:, 2 * ik - 1)
CALL cryst_to_cart(1, xkk, bg, 1)
WRITE(stdout,'(3i8,4f12.6,3E14.5)') itemp, ik, ibnd, xkk, ekk, F_SERTA(:, ibnd, ik, itemp)
WRITE(stdout, '(3i8,4f12.6,3E14.5)') itemp, ik, ibnd, xkk, ekk, F_SERTA(:, ibnd, ik, itemp)
ENDIF
ENDIF ! iverbosity 4
ENDIF
@ -551,10 +552,13 @@
ENDDO ! end of while loop
!
! Deallocate
DEALLOCATE(xkf_sp)
DEALLOCATE(xkf_sp, STAT = ierr)
IF (ierr /= 0) CALL errore('ibte', 'Error deallocating xkf_sp', 1)
IF (mp_mesh_k) THEN
DEALLOCATE(ixkqf_tr)
DEALLOCATE(s_BZtoIBZ_full)
DEALLOCATE(ixkqf_tr, STAT = ierr)
IF (ierr /= 0) CALL errore('ibte', 'Error deallocating ixkqf_tr', 1)
DEALLOCATE(s_BZtoIBZ_full, STAT = ierr)
IF (ierr /= 0) CALL errore('ibte', 'Error deallocating s_BZtoIBZ_full', 1)
ENDIF
!
RETURN
@ -619,10 +623,8 @@
!! Fermi level for the temperature itemp for cb band
!
! Local variables
!
CHARACTER(LEN = 256) :: filint
!! Name of the file to write/read
!
INTEGER :: ierr
!! Error status
INTEGER :: ik
@ -682,7 +684,6 @@
INTEGER(KIND = 8) :: upper_bnd
!! end for current CPU
#endif
!
REAL(KIND = DP) :: dum1
!! Dummy variable
REAL(KIND = DP), ALLOCATABLE :: trans_prob(:)
@ -764,7 +765,8 @@
! Allocate the local size
nind = upper_bnd - lower_bnd + 1
WRITE(stdout, '(5x,a,i10)') 'Number of elements per core ', nind
ALLOCATE(trans_prob(nind))
ALLOCATE(trans_prob(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating trans_prob', 1)
trans_prob(:) = 0.0d0
!
! Open file containing trans_prob
@ -795,11 +797,16 @@
CALL MPI_FILE_OPEN(world_comm, 'sparset', MPI_MODE_RDONLY, MPI_INFO_NULL, iunsparset, ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'error in MPI_FILE_OPEN sparset', 1)
!
ALLOCATE(sparse_q(nind))
ALLOCATE(sparse_k(nind))
ALLOCATE(sparse_i(nind))
ALLOCATE(sparse_j(nind))
ALLOCATE(sparse_t(nind))
ALLOCATE(sparse_q(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparse_q', 1)
ALLOCATE(sparse_k(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparse_k', 1)
ALLOCATE(sparse_i(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparse_i', 1)
ALLOCATE(sparse_j(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparse_j', 1)
ALLOCATE(sparse_t(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparse_t', 1)
sparse_q(:) = 0.0d0
sparse_k(:) = 0.0d0
sparse_i(:) = 0.0d0
@ -845,12 +852,18 @@
IF (ierr /= 0) CALL errore('iter_restart', 'error in MPI_FILE_CLOSE', 1)
CALL MPI_FILE_CLOSE(iunsparset, ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'error in MPI_FILE_CLOSE', 1)
DEALLOCATE(trans_prob)
DEALLOCATE(sparse_q)
DEALLOCATE(sparse_k)
DEALLOCATE(sparse_i)
DEALLOCATE(sparse_j)
DEALLOCATE(sparse_t)
DEALLOCATE(trans_prob, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating trans_prob', 1)
DEALLOCATE(sparse_q, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparse_q', 1)
DEALLOCATE(sparse_k, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparse_k', 1)
DEALLOCATE(sparse_i, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparse_i', 1)
DEALLOCATE(sparse_j, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparse_j', 1)
DEALLOCATE(sparse_t, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparse_t', 1)
!
ENDIF
! Electrons
@ -860,7 +873,8 @@
! Allocate the local size
nind = upper_bnd - lower_bnd + 1
WRITE(stdout, '(5x,a,i10)') 'Number of elements per core ', nind
ALLOCATE(trans_probcb(nind))
ALLOCATE(trans_probcb(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating trans_probcb', 1)
trans_probcb(:) = 0.0d0
!
! Open file containing trans_prob
@ -891,11 +905,16 @@
CALL MPI_FILE_OPEN(world_comm, 'sparsetcb', MPI_MODE_RDONLY, MPI_INFO_NULL, iunsparsetcb, ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'error in MPI_FILE_OPEN sparsetcb', 1)
!
ALLOCATE(sparsecb_q(nind))
ALLOCATE(sparsecb_k(nind))
ALLOCATE(sparsecb_i(nind))
ALLOCATE(sparsecb_j(nind))
ALLOCATE(sparsecb_t(nind))
ALLOCATE(sparsecb_q(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparsecb_q', 1)
ALLOCATE(sparsecb_k(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparsecb_k', 1)
ALLOCATE(sparsecb_i(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparsecb_i', 1)
ALLOCATE(sparsecb_j(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparsecb_j', 1)
ALLOCATE(sparsecb_t(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error allocating sparsecb_t', 1)
sparsecb_q(:) = 0.0d0
sparsecb_k(:) = 0.0d0
sparsecb_i(:) = 0.0d0
@ -940,12 +959,18 @@
IF (ierr /= 0) CALL errore('iter_restart', 'error in MPI_FILE_CLOSE', 1)
CALL MPI_FILE_CLOSE(iunsparsetcb, ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'error in MPI_FILE_CLOSE', 1)
DEALLOCATE(trans_probcb)
DEALLOCATE(sparsecb_q)
DEALLOCATE(sparsecb_k)
DEALLOCATE(sparsecb_i)
DEALLOCATE(sparsecb_j)
DEALLOCATE(sparsecb_t)
DEALLOCATE(trans_probcb, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating trans_probcb', 1)
DEALLOCATE(sparsecb_q, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparsecb_q', 1)
DEALLOCATE(sparsecb_k, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparsecb_k', 1)
DEALLOCATE(sparsecb_i, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparsecb_i', 1)
DEALLOCATE(sparsecb_j, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparsecb_j', 1)
DEALLOCATE(sparsecb_t, STAT = ierr)
IF (ierr /= 0) CALL errore('iter_restart', 'Error deallocating sparsecb_t', 1)
!
ENDIF
#endif

22
EPW/src/wan2bloch.f90
View File

@ -54,10 +54,8 @@
!! number of WS points
INTEGER, INTENT(in) :: dims
!! dims = nbndsub if use_ws or 1 otherwise
!
REAL(KIND = DP), INTENT(out) :: eig (nbnd)
!! interpolated hamiltonian eigenvalues for this kpoint
!
COMPLEX(KIND = DP), INTENT(in) :: cfac(nrr, dims, dims)
!! Exponential factor
COMPLEX(KIND = DP), INTENT(in) :: chw(nbnd, nbnd, nrr)
@ -153,7 +151,7 @@
ENDDO
ENDDO
ELSE
CALL zgemv('n', nbnd**2, nrr, cone, chw, nbnd**2, cfac(:, 1, 1), 1, cone, chf, 1)
CALL ZGEMV('n', nbnd**2, nrr, cone, chw, nbnd**2, cfac(:, 1, 1), 1, cone, chf, 1)
ENDIF
!
!---------------------------------------------------------------------
@ -327,7 +325,6 @@
!! coordinates of phononic WS points
INTEGER, INTENT(in) :: ndegen_q(nrr_q, nat, nat)
!! degeneracy of WS points
!
REAL(KIND = DP), INTENT(in) :: xxq(3)
!! kpoint coordinates for the interpolation
REAL(KIND = DP), INTENT(out) :: eig(nmodes)
@ -959,7 +956,7 @@
ENDDO
ENDDO
ELSE
CALL zgemv('n', 3 * (nbnd**2), nrr, cone, cdmew(:, :, :, :), 3 * (nbnd**2), cfac(:, 1, 1), 1, cone, cdmef(:, :, :), 1)
CALL ZGEMV('n', 3 * (nbnd**2), nrr, cone, cdmew(:, :, :, :), 3 * (nbnd**2), cfac(:, 1, 1), 1, cone, cdmef(:, :, :), 1)
ENDIF
!
!----------------------------------------------------------
@ -1130,7 +1127,7 @@
ENDDO
ENDDO
ELSE
CALL zgemv('n', 3 * (nbnd**2), nrr, cone, cvmew(:, :, :, :), 3 * (nbnd**2), cfac(:, 1, 1), 1, cone, cvmef(:, :, :), 1)
CALL ZGEMV('n', 3 * (nbnd**2), nrr, cone, cvmew(:, :, :, :), 3 * (nbnd**2), cfac(:, 1, 1), 1, cone, cvmef(:, :, :), 1)
ENDIF
!
! k-derivative of the Hamiltonian in the Wannier gauge
@ -1274,7 +1271,7 @@
ENDDO
ENDDO
!
CALL zhpevx('V', 'A', 'U', deg_dim(ideg), champ , zero, zero, &
CALL ZHPEVX('V', 'A', 'U', deg_dim(ideg), champ , zero, zero, &
0, 0, -one, neig, w, cz, deg_dim(ideg), cwork, rwork, iwork, ifail, info)
!
vmef_deg(ipol, :, :) = zero
@ -1372,7 +1369,7 @@
!! interpolated velocity matrix elements in Bloch basis, fine mesh, in the degenerate subspaces
!
! Local variables
LOGICAL, SAVE :: first=.TRUE.
LOGICAL, SAVE :: first = .TRUE.
!! First entrance [used when lifc == .TRUE.]
LOGICAL :: duplicates
!! Returns if the modes contains degeneracices for that q-point.
@ -1666,7 +1663,7 @@
ENDDO
ENDDO
!
CALL zhpevx('V', 'A', 'U', deg_dim(ideg), champ , zero, zero, &
CALL ZHPEVX('V', 'A', 'U', deg_dim(ideg), champ , zero, zero, &
0, 0, -one, neig, w, cz, deg_dim(ideg), cwork, rwork, iwork, ifail, info)
!
vmef_deg(ipol, :, :) = zero
@ -1738,8 +1735,6 @@
!
IMPLICIT NONE
!
! input variables
!
INTEGER, INTENT(in) :: nmodes
!! Total number of modes
INTEGER, INTENT(in) :: nrr_g
@ -1764,7 +1759,6 @@
!! e-p matrix in Bloch representation, fine grid
!
! Local variables
!
INTEGER :: ir
!! Real space WS index
INTEGER :: iw
@ -2077,7 +2071,6 @@
!! Counter on the number of Wannier functions
INTEGER :: imode
!! Counter on phonon modes
!
COMPLEX(KIND = DP) :: eptmp(nbnd, nbnd)
!! Temporary variable
CALL start_clock('ephW2B')
@ -2196,7 +2189,7 @@
ENDDO
ENDDO
ELSE
CALL zgemv('n', nbnd**2, nrr, cone, epmatw(:, :, :), nbnd**2, cfac(:, 1, 1), 1, cone, epmatf(:, :), 1)
CALL ZGEMV('n', nbnd**2, nrr, cone, epmatw(:, :, :), nbnd**2, cfac(:, 1, 1), 1, cone, epmatf(:, :), 1)
ENDIF
!
!----------------------------------------------------------
@ -2305,7 +2298,6 @@
!
REAL(KIND = DP) :: rdotk
!! Exponential for the FT
!
COMPLEX(KIND = DP) :: cfac(nrr_g, dims, dims)
!! Factor for the FT
COMPLEX(KIND = DP), ALLOCATABLE :: epmatw(:, :, :)

14
EPW/src/wannierEPW.f90
View File

@ -13,10 +13,16 @@
MODULE wannierEPW
!-----------------------------------------------------------------------------------------
!!
!!
!! Interface between Wannier and EPW
!!
USE kinds, ONLY : DP
!
CHARACTER(LEN = 15) :: wan_mode
!! running mode
CHARACTER(LEN = 256) :: seedname2
!! prepended to file names in wannier90. For implementation of wannier_lib
CHARACTER(LEN = 3), ALLOCATABLE :: atsym(:)
!! atomic symbols. atsym(nat)
LOGICAL :: logwann
!!
LOGICAL :: write_unk
@ -39,12 +45,6 @@
!! lwindow(num_bands,iknum)
LOGICAL, ALLOCATABLE :: excluded_band(:)
!! list of bands to exclude from the calculation of WFs
CHARACTER(LEN = 15) :: wan_mode
!! running mode
CHARACTER(LEN = 256) :: seedname2
!! prepended to file names in wannier90. For implementation of wannier_lib
CHARACTER(LEN = 3), ALLOCATABLE :: atsym(:)
!! atomic symbols. atsym(nat)
INTEGER :: nnb
!! number of b-vectors
INTEGER :: iun_nnkp

60
EPW/src/wannierize.f90
View File

@ -7,10 +7,10 @@
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!--------------------------------------------------------------------
SUBROUTINE wann_run
SUBROUTINE wann_run()
!---------------------------------------------------------------------
!!
!! This is the SUBROUTINE which controls the w90 run. Primarily,
!! This is the routine which controls the w90 run. Primarily,
!! we get the phases to remove degeneracies in the wfs, and
!! call pw2wan90epw
!!
@ -30,7 +30,8 @@
! Local variables
INTEGER :: num_kpts
!! number of k-points in wannierization
!
INTEGER :: ierr
!! Error status
!
CALL start_clock('WANNIER')
!
@ -40,9 +41,10 @@
num_kpts = mp_grid(1) * mp_grid(2) * mp_grid(3)
!
IF (num_kpts /= nkstot) CALL errore('wannierize', 'inconsistent nscf and elph k-grids', 1)
IF (nbnd < n_wannier ) CALL errore('wannierize', 'Must have as many or more bands than Wannier functions', 1)
IF (nbnd < n_wannier) CALL errore('wannierize', 'Must have as many or more bands than Wannier functions', 1)
!
ALLOCATE(kpt_latt(3, num_kpts))
ALLOCATE(kpt_latt(3, num_kpts), STAT = ierr)
IF (ierr /= 0) CALL errore('wann_run', 'Error allocating kpt_latt', 1)
!
WRITE(stdout, '(5x,a)') REPEAT("-",67)
WRITE(stdout, '(a, i2,a,i2,a,i2,a)') " Wannierization on ", nkc1, " x ", nkc2, " x ", nkc3 , " electronic grid"
@ -53,15 +55,16 @@
!
! write the short input file for the wannier90 code
!
CALL write_winfil
CALL write_winfil()
!
! run the wannier90 code to create MLWFs
!
CALL pw2wan90epw
CALL pw2wan90epw()
!
! project the Wannier functions onto energy space
!
DEALLOCATE(kpt_latt)
DEALLOCATE(kpt_latt, STAT = ierr)
IF (ierr /= 0) CALL errore('wann_run', 'Error deallocating kpt_latt', 1)
!
WRITE(stdout, '(5x,a)') REPEAT("-",67)
CALL print_clock('WANNIER')
@ -72,15 +75,13 @@
!------------------------------------------------------------
!
!------------------------------------------------------------
SUBROUTINE write_winfil
SUBROUTINE write_winfil()
!------------------------------------------------------------
!!
!!
!! This SUBROUTINE writes the prefix.win file which wannier90.x
!! This routine writes the prefix.win file which wannier90.x
!! needs to run. Primarily it contains information about the
!! windows used for the disentanglement, and the initial projections.
!! JN - 10/2008 projections now in elph.in file
!------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_files, ONLY : prefix
@ -95,10 +96,11 @@
!
IMPLICIT NONE
!
! Local variables
LOGICAL :: random
!! Random
INTEGER :: i
!!
!! Band index
REAL(KIND = DP) :: et_tmp(nbnd, nkstot)
!! eigenvalues on full coarse k-mesh
!
@ -151,17 +153,17 @@
!
CLOSE(iuwinfil)
!
ENDIF
ENDIF ! meta_ionode
!
!------------------------------------------------------------
END SUBROUTINE write_winfil
!------------------------------------------------------------
!
!------------------------------------------------------------
SUBROUTINE proj_w90
SUBROUTINE proj_w90()
!------------------------------------------------------------
!!
!! This SUBROUTINE computes the energy projections of
!! This routine computes the energy projections of
!! the computed Wannier functions
!! 07/2010 Needs work. Right now this sub is nearly worthless
!
@ -196,6 +198,8 @@
!!
INTEGER :: iwann
!!
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: dE
!!
REAL(KIND = DP) :: sigma
@ -226,14 +230,19 @@
ne = INT((dis_win_max - dis_win_min + 1) / dE)
IF (ne < 1) CALL errore('proj_wan', 'Problem with disentanglement window', 1)
!
ALLOCATE(proj_wf(n_wannier, ne + 1))
ALLOCATE(proj_wf(n_wannier, ne + 1), STAT = ierr)
IF (ierr /= 0) CALL errore('proj_w90', 'Error allocating proj_wf', 1)
proj_wf = 0.d0
!
ALLOCATE(cu(nbnd, n_wannier, nks))
ALLOCATE(cuq(nbnd, n_wannier, nks))
ALLOCATE(xkq(3, nks))
ALLOCATE(cu(nbnd, n_wannier, nks), STAT = ierr)
IF (ierr /= 0) CALL errore('proj_w90', 'Error allocating cu', 1)
ALLOCATE(cuq(nbnd, n_wannier, nks), STAT = ierr)
IF (ierr /= 0) CALL errore('proj_w90', 'Error allocating cuq', 1)
ALLOCATE(xkq(3, nks), STAT = ierr)
IF (ierr /= 0) CALL errore('proj_w90', 'Error allocating xkq', 1)
CALL loadumat(nbnd, n_wannier, nks, nkstot, xxq, cu, cuq, lwin, lwinq, exband)
DEALLOCATE(xkq)
DEALLOCATE(xkq, STAT = ierr)
IF (ierr /= 0) CALL errore('proj_w90', 'Error deallocating xkq', 1)
!
DO iwann = 1, n_wannier
DO ie = 1, ne
@ -263,9 +272,12 @@
!
CLOSE(iuprojfil)
ENDIF
DEALLOCATE(proj_wf)
DEALLOCATE(cu)
DEALLOCATE(cuq)
DEALLOCATE(proj_wf, STAT = ierr)
IF (ierr /= 0) CALL errore('proj_w90', 'Error deallocating proj_wf', 1)
DEALLOCATE(cu, STAT = ierr)
IF (ierr /= 0) CALL errore('proj_w90', 'Error deallocating cu', 1)
DEALLOCATE(cuq, STAT = ierr)
IF (ierr /= 0) CALL errore('proj_w90', 'Error deallocating cuq', 1)
!
!------------------------------------------------------------
END SUBROUTINE proj_w90

49
EPW/src/wigner.f90
View File

@ -110,6 +110,8 @@
!! maximum number of WS vectors for the phonons
INTEGER :: nrr_g
!! maximum number of WS vectors for the electron-phonon
INTEGER :: ierr
!! Error status
INTEGER :: irvec_kk (3, 20 * nkc1 * nkc2 * nkc3)
!! local INTEGER components of the ir-th Wigner-Seitz grid point
!! in the basis of the lattice vectors for electrons
@ -142,15 +144,24 @@
CALL wigner_seitzkq(nqc1, nqc2, nqc3, irvec_qq, ndegen_qq, wslen_qq, nrr_q, tau, dims2)
CALL wigner_seitzg(nqc1, nqc2, nqc3, irvec_gg, ndegen_gg, wslen_gg, nrr_g, w_centers, tau, dims, dims2)
!
ALLOCATE(irvec_k(3, nrr_k))
ALLOCATE(irvec_q(3, nrr_q))
ALLOCATE(irvec_g(3, nrr_g))
ALLOCATE(ndegen_k(nrr_k, dims, dims))
ALLOCATE(ndegen_q(nrr_q, dims2, dims2))
ALLOCATE(ndegen_g(nrr_g, dims2, dims, dims))
ALLOCATE(wslen_k(nrr_k))
ALLOCATE(wslen_q(nrr_q))
ALLOCATE(wslen_g(nrr_g))
ALLOCATE(irvec_k(3, nrr_k), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating irvec_k', 1)
ALLOCATE(irvec_q(3, nrr_q), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating irvec_q', 1)
ALLOCATE(irvec_g(3, nrr_g), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating irvec_g', 1)
ALLOCATE(ndegen_k(nrr_k, dims, dims), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating ndegen_k', 1)
ALLOCATE(ndegen_q(nrr_q, dims2, dims2), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating ndegen_q', 1)
ALLOCATE(ndegen_g(nrr_g, dims2, dims, dims), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating ndegen_g', 1)
ALLOCATE(wslen_k(nrr_k), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating wslen_k', 1)
ALLOCATE(wslen_q(nrr_q), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating wslen_q', 1)
ALLOCATE(wslen_g(nrr_g), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating wslen_g', 1)
!
! Create vectors with correct size.
DO ir = 1, nrr_k
@ -224,6 +235,8 @@
!! Index of sorting
INTEGER :: nind
!! The metric tensor
INTEGER :: ierr
!! Error status
INTEGER :: nrr_tmp(dims, dims)
!! Temporary WS matrix
INTEGER :: irvec_tmp(3, 20 * nc1 * nc2 * nc3, dims, dims)
@ -245,7 +258,8 @@
IF (nind < 125) THEN
nind = 125
ENDIF
ALLOCATE(ind(nind))
ALLOCATE(ind(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitzkq', 'Error allocating ind', 1)
!
DO ipol = 1, 3
DO jpol = 1, 3
@ -363,7 +377,7 @@
tot = tot + 1.d0 / DBLE(ndegen_tmp(i, iw, iw2))
ENDDO
!
IF (ABS(tot - DBLE(nc1*nc2*nc3)) > eps6) call errore &
IF (ABS(tot - DBLE(nc1*nc2*nc3)) > eps6) CALL errore &
('wigner_seitzkq',' weights do not add up to nc1*nc2*nc3', 1)
IF (ABS(tot - tot2) > eps6) CALL errore &
('wigner_seitzkq', ' weigths of pair of atoms is not equal to global weights', 1)
@ -385,7 +399,8 @@
!
CALL cryst_to_cart(dims, shift(:, :), at, 1)
!
DEALLOCATE(ind)
DEALLOCATE(ind, STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitzkq', 'Error deallocating ind', 1)
!
!-----------------------------------------------------------------------------
END SUBROUTINE wigner_seitzkq
@ -454,6 +469,8 @@
!! Index of sorting
INTEGER :: nind
!! The metric tensor
INTEGER :: ierr
!! Error status
INTEGER :: nrr_tmp(dims2, dims, dims)
!! Temporary array that contains the max number of WS vectors
!! for a pair of atoms.
@ -476,7 +493,8 @@
IF (nind < 125) THEN
nind = 125
ENDIF
ALLOCATE(ind(nind))
ALLOCATE(ind(nind), STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitzg', 'Error allocating ind', 1)
!
DO ipol = 1, 3
DO jpol = 1, 3
@ -537,7 +555,7 @@
found = .FALSE.
i = 1
mindist = dist(1)
DO WHILE (ABS(dist(i) - mindist) < eps6 .AND. i < 125 )
DO WHILE (ABS(dist(i) - mindist) < eps6 .AND. i < 125)
IF (ind(i) == 63) found = .TRUE.
i = i + 1
ENDDO
@ -633,7 +651,8 @@
CALL cryst_to_cart(dims2, tau(:, :), at, 1)
CALL cryst_to_cart(dims, w_centers(:, :), at, 1)
!
DEALLOCATE(ind)
DEALLOCATE(ind, STAT = ierr)
IF (ierr /= 0) CALL errore('wigner_seitzg', 'Error deallocating ind', 1)
!
!------------------------------------------------------------------------------------------
END SUBROUTINE wigner_seitzg