From 8392d54d2d1d210aded5477ec188acbb48950e90 Mon Sep 17 00:00:00 2001 From: Samuel Ponce Date: Fri, 13 Sep 2019 17:28:34 +0100 Subject: [PATCH] Final cleaning and modularization. --- EPW/src/Makefile | 4 +- EPW/src/bloch2wan.f90 | 6 +- EPW/src/rotate_eigenm.f90 | 51 ++++---- EPW/src/rotate_epmat.f90 | 17 ++- EPW/src/selfen.f90 | 232 +++++++++++++++++++------------------ EPW/src/setphases_wrap.f90 | 30 ++--- EPW/src/spectral_func.f90 | 47 +++++--- EPW/src/star_q2.f90 | 8 +- EPW/src/stop_epw.f90 | 19 ++- EPW/src/transport.f90 | 91 ++++++++------- EPW/src/transport_iter.f90 | 125 ++++++++++++-------- EPW/src/wan2bloch.f90 | 22 ++-- EPW/src/wannierEPW.f90 | 14 +-- EPW/src/wannierize.f90 | 60 ++++++---- EPW/src/wigner.f90 | 49 +++++--- 15 files changed, 418 insertions(+), 357 deletions(-) diff --git a/EPW/src/Makefile b/EPW/src/Makefile index 88b3cb030..0b413483a 100644 --- a/EPW/src/Makefile +++ b/EPW/src/Makefile @@ -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 \ diff --git a/EPW/src/bloch2wan.f90 b/EPW/src/bloch2wan.f90 index 46f26c753..79454d5c4 100644 --- a/EPW/src/bloch2wan.f90 +++ b/EPW/src/bloch2wan.f90 @@ -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') diff --git a/EPW/src/rotate_eigenm.f90 b/EPW/src/rotate_eigenm.f90 index b0d403624..dac8942b8 100644 --- a/EPW/src/rotate_eigenm.f90 +++ b/EPW/src/rotate_eigenm.f90 @@ -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 ! diff --git a/EPW/src/rotate_epmat.f90 b/EPW/src/rotate_epmat.f90 index ef0c8c6c5..86e905423 100644 --- a/EPW/src/rotate_epmat.f90 +++ b/EPW/src/rotate_epmat.f90 @@ -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 !--------------------------------------------------------------------------- - diff --git a/EPW/src/selfen.f90 b/EPW/src/selfen.f90 index 710c8a070..7699fdbb6 100644 --- a/EPW/src/selfen.f90 +++ b/EPW/src/selfen.f90 @@ -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 !----------------------------------------------------------------------- - diff --git a/EPW/src/setphases_wrap.f90 b/EPW/src/setphases_wrap.f90 index e56e20325..878af74a4 100644 --- a/EPW/src/setphases_wrap.f90 +++ b/EPW/src/setphases_wrap.f90 @@ -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 + !--------------------------------------------------------------------- diff --git a/EPW/src/spectral_func.f90 b/EPW/src/spectral_func.f90 index ed20036d1..7e3d410c5 100644 --- a/EPW/src/spectral_func.f90 +++ b/EPW/src/spectral_func.f90 @@ -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 !----------------------------------------------------------------------- - diff --git a/EPW/src/star_q2.f90 b/EPW/src/star_q2.f90 index a9f7dcfc9..0c2db22bf 100644 --- a/EPW/src/star_q2.f90 +++ b/EPW/src/star_q2.f90 @@ -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 diff --git a/EPW/src/stop_epw.f90 b/EPW/src/stop_epw.f90 index 8c027736b..f0760c7f9 100644 --- a/EPW/src/stop_epw.f90 +++ b/EPW/src/stop_epw.f90 @@ -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 diff --git a/EPW/src/transport.f90 b/EPW/src/transport.f90 index 34a450a2c..7ac03e86a 100644 --- a/EPW/src/transport.f90 +++ b/EPW/src/transport.f90 @@ -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)') diff --git a/EPW/src/transport_iter.f90 b/EPW/src/transport_iter.f90 index d124e50b3..2c3a230f8 100644 --- a/EPW/src/transport_iter.f90 +++ b/EPW/src/transport_iter.f90 @@ -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 diff --git a/EPW/src/wan2bloch.f90 b/EPW/src/wan2bloch.f90 index 95188a727..5638fafb3 100644 --- a/EPW/src/wan2bloch.f90 +++ b/EPW/src/wan2bloch.f90 @@ -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(:, :, :) diff --git a/EPW/src/wannierEPW.f90 b/EPW/src/wannierEPW.f90 index fb72e6c0e..e3dd05efe 100644 --- a/EPW/src/wannierEPW.f90 +++ b/EPW/src/wannierEPW.f90 @@ -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 diff --git a/EPW/src/wannierize.f90 b/EPW/src/wannierize.f90 index 483cb746f..93a67fe7a 100644 --- a/EPW/src/wannierize.f90 +++ b/EPW/src/wannierize.f90 @@ -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 diff --git a/EPW/src/wigner.f90 b/EPW/src/wigner.f90 index 442180b0a..68476586e 100644 --- a/EPW/src/wigner.f90 +++ b/EPW/src/wigner.f90 @@ -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