diff --git a/EPW/src/star_q2.f90 b/EPW/src/star_q2.f90 index 17b297bee..6f2595de4 100644 --- a/EPW/src/star_q2.f90 +++ b/EPW/src/star_q2.f90 @@ -10,123 +10,141 @@ ! Adapted from QE ! !----------------------------------------------------------------------- - SUBROUTINE star_q2 (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, verbosity, sym_smallq ) + SUBROUTINE star_q2(xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, verbosity, sym_smallq) !----------------------------------------------------------------------- - ! - ! Generate the star of q vectors that are equivalent to the input one - ! NB: input s(:,:,1:nsym) must contain all crystal symmetries, - ! i.e. not those of the small-qroup of q only - ! + !! + !! Generate the star of q vectors that are equivalent to the input one + !! NB: input s(:,:,1:nsym) must contain all crystal symmetries, + !! i.e. not those of the small-qroup of q only + !! + !! SP - Sept. 2019 - Cleaning USE io_global, ONLY : stdout USE kinds, ONLY : DP + ! IMPLICIT NONE ! - REAL(DP), parameter :: accep=1.e-5_dp - ! + LOGICAL, INTENT(in) :: verbosity + !! if true prints several messages. + INTEGER, INTENT(in) :: nsym + !! nsym matrices of symmetry operations + INTEGER, INTENT(in) :: s(3, 3, 48) + !! Symmetry operations + INTEGER, INTENT(in) :: invs(48) + !! list of inverse operation indices + INTEGER, INTENT(out) :: nq + !! degeneracy of the star of q + INTEGER, INTENT(out) :: isq(48) + !! index of q in the star for a given sym + INTEGER, INTENT(out) :: imq + !! index of -q in the star (0 if not present) + REAL(KIND = DP), INTENT(in) :: xq(3) + !! q vector + REAL(KIND = DP), INTENT(in) :: at(3, 3) + !! direct lattice vectors + REAL(KIND = DP), INTENT(in) :: bg(3, 3) + !! reciprocal lattice vectors + REAL(KIND = DP), INTENT(in) :: + !! + REAL(KIND = DP), INTENT(out) :: sxq(3, 48) + !! list of vectors in the star of q + ! + ! Local variables + LOGICAL, EXTERNAL :: eqvect + !! function used to compare two vectors INTEGER :: sym_smallq(48) - ! - integer, intent(in) :: nsym, s (3, 3, 48), invs(48) - ! nsym matrices of symmetry operations - ! invs: list of inverse operation indices - REAL(DP), intent(in) :: xq (3), at (3, 3), bg (3, 3) - ! xq: q vector - ! at: direct lattice vectors - ! bg: reciprocal lattice vectors - ! - integer, intent(out) :: nq, isq (48), imq - ! nq : degeneracy of the star of q - ! isq : index of q in the star for a given sym - ! imq : index of -q in the star (0 if not present) - ! - REAL(DP), intent(out) :: sxq (3, 48) - ! list of vectors in the star of q - logical, intent(in) :: verbosity - ! if true prints several messages. - ! - INTEGER :: nsq (48), isym, ism1, iq, i - ! number of symmetry ops. of bravais lattice - ! counters on symmetry ops. - ! index of inverse of isym - ! counters - REAL(DP) :: saq (3, 48), aq (3), raq (3), zero (3) - ! auxiliary list of q (crystal coordinates) - ! input q in crystal coordinates - ! rotated q in crystal coordinates - ! coordinates of fractionary translations - ! a zero vector: used in eqvect - ! - logical, EXTERNAL :: eqvect - ! function used to compare two vectors - ! + !! + INTEGER :: nsq(48) + !! number of symmetry ops. of bravais lattice + INTEGER :: isym + !! counters on symmetry ops. + INTEGER :: ism1 + !! index of inverse of isym + INTEGER :: iq + !! q-counter + INTEGER :: i + !! Counter + REAL(KIND = DP) :: saq(3, 48) + !! auxiliary list of q (crystal coordinates) + REAL(KIND = DP) :: aq(3) + !! input q in crystal coordinates + REAL(KIND = DP) :: raq(3) + !! rotated q in crystal coordinates + REAL(KIND = DP) :: zero(3) + !! a zero vector: used in eqvect + REAL(KIND = DP), PARAMETER :: accep = 1.e-5_dp + !! Tolerence + zero(:) = 0.d0 ! ! go to crystal coordinates ! - do i = 1, 3 - aq(i) = xq(1) * at(1,i) + xq(2) * at(2,i) + xq(3) * at(3,i) - enddo + DO i = 1, 3 + aq(i) = xq(1) * at(1,i) + xq(2) * at(2,i) + xq(3) * at(3,i) + ENDDO ! ! create the list of rotated q ! - do i = 1, 48 - nsq (i) = 0 - isq (i) = 0 - enddo + DO i = 1, 48 + nsq(i) = 0 + isq(i) = 0 + ENDDO nq = 0 - do isym = 1, nsym - ism1 = invs (isym) - do i = 1, 3 - raq (i) = s (i, 1, ism1) * aq (1) & - + s (i, 2, ism1) * aq (2) & - + s (i, 3, ism1) * aq (3) - enddo - do i = 1, 3 - sxq (i, 48) = bg (i, 1) * raq (1) & - + bg (i, 2) * raq (2) & - + bg (i, 3) * raq (3) - enddo - do iq = 1, nq - if (eqvect (raq, saq (1, iq), zero, accep) ) then - isq (isym) = iq - nsq (iq) = nsq (iq) + 1 - endif - enddo - if (isq (isym) == 0) then - nq = nq + 1 - nsq (nq) = 1 - isq (isym) = nq - 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) - enddo - endif - enddo + DO isym = 1, nsym + ism1 = invs(isym) + DO i = 1, 3 + raq(i) = s(i, 1, ism1) * aq(1) & + + s(i, 2, ism1) * aq(2) & + + s(i, 3, ism1) * aq(3) + ENDDO + DO i = 1, 3 + sxq(i, 48) = bg(i, 1) * raq(1) & + + bg(i, 2) * raq(2) & + + bg(i, 3) * raq(3) + ENDDO + DO iq = 1, nq + IF (eqvect(raq, saq(1, iq), zero, accep)) THEN + isq(isym) = iq + nsq(iq) = nsq(iq) + 1 + ENDDIF + ENDDO + IF (isq(isym) == 0) THEN + nq = nq + 1 + nsq(nq) = 1 + isq(isym) = nq + 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) + ENDDO + ENDIF + ENDDO ! ! SP: Now determine which q in the star is obtained ! ! Set imq index if needed and check star degeneracy ! - raq (:) = - aq(:) + raq(:) = -aq(:) imq = 0 - do iq = 1, nq - if (eqvect (raq, saq (1, iq), zero, accep) ) imq = iq - if (nsq(iq)*nq /= nsym) call errore ('star_q', 'wrong degeneracy', iq) - enddo + DO iq = 1, nq + IF (eqvect(raq, saq(1, iq), zero, accep)) imq = iq + IF (nsq(iq) * nq /= nsym) CALL errore('star_q', 'wrong degeneracy', iq) + ENDDO ! ! writes star of q ! IF (verbosity) THEN - WRITE( stdout, * ) - WRITE( stdout, '(5x,a,i4)') 'Number of q in the star = ', nq - WRITE( stdout, '(5x,a)') 'List of q in the star:' - WRITE( stdout, '(7x,i4,3f14.9)') (iq, (sxq(i,iq), i = 1,3), iq = 1,nq) - if (imq == 0) then - WRITE( stdout, '(5x,a)') 'In addition there is the -q list: ' - WRITE( stdout, '(7x,i4,3f14.9)') (iq, (-sxq(i,iq), i = 1,3), iq = 1,nq) - endif + WRITE(stdout, *) + WRITE(stdout, '(5x,a,i4)') 'Number of q in the star = ', nq + WRITE(stdout, '(5x,a)') 'List of q in the star:' + WRITE(stdout, '(7x,i4,3f14.9)') (iq, (sxq(i, iq), i = 1, 3), iq = 1, nq) + IF (imq == 0) THEN + WRITE(stdout, '(5x,a)') 'In addition there is the -q list: ' + WRITE(stdout, '(7x,i4,3f14.9)') (iq, (-sxq(i, iq), i = 1, 3), iq = 1, nq) + ENDIF ENDIF - return -END SUBROUTINE star_q2 + RETURN + !----------------------------------------------------------------------- + END SUBROUTINE star_q2 + !----------------------------------------------------------------------- diff --git a/EPW/src/transport.f90 b/EPW/src/transport.f90 index af8bd6959..bffb26d84 100644 --- a/EPW/src/transport.f90 +++ b/EPW/src/transport.f90 @@ -32,17 +32,17 @@ 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, & - xqf, zi_allvb, zi_allcb + xqf, zi_allvb, zi_allcb, nbndfst, nktotf USE transportcom, ONLY : transp_temp, lower_bnd USE constants_epw, ONLY : zero, one, two, pi, ryd2mev, kelvin2eV, ryd2ev, & - eps6, eps8 + eps6, eps8, eps4 USE mp, ONLY : mp_barrier, mp_sum USE mp_global, ONLY : world_comm USE io_scattering, ONLY : scattering_write, tau_write, merge_read ! IMPLICIT NONE ! - LOGICAL, INTENT (INOUT) :: first_cycle + LOGICAL, INTENT(inout) :: first_cycle !! Use to determine weather this is the first cycle after restart INTEGER, INTENT(in) :: iqq !! Q-point index from the selected q @@ -112,11 +112,11 @@ !! Fermi-Dirac occupation function $$f_{m\mathbf{k+q}}$$ REAL(KIND = DP) :: trans_prob !! Transition probability function - REAL(KIND = DP) :: vkk(3,nbndfst) + REAL(KIND = DP) :: vkk(3, nbndfst) !! Electronic velocity $$v_{n\mathbf{k}}$$ - REAL(KIND = DP) :: vkq(3,nbndfst) + REAL(KIND = DP) :: vkq(3, nbndfst) !! Electronic velocity $$v_{m\mathbf{k+q}}$$ - REAL(KIND = DP) :: vel_factor(nbndfst,nbndfst) + REAL(KIND = DP) :: vel_factor(nbndfst, nbndfst) !! Velocity factor $$ 1 - \frac{(v_{nk} \cdot v_{mk+q})}{ |v_{nk}|^2} $$ REAL(KIND = DP) :: inv_tau_tmp(nbndfst) !! Temporary array to store the scattering rates @@ -135,26 +135,21 @@ !! Compute the approximate theta function. Here computes Fermi-Dirac REAL(KIND = DP), EXTERNAL :: w0gauss !! The derivative of wgauss: an approximation to the delta function - REAL(KIND = DP), PARAMETER :: eps = 1.d-4 - !! Tolerence parameter for the velocity ! - CALL start_clock ( 'SCAT' ) + CALL start_clock('SCAT') ! IF (iqq == 1) THEN ! - WRITE(stdout,'(/5x,a)') REPEAT('=',67) - WRITE(stdout,'(5x,"Scattering rate")') - WRITE(stdout,'(5x,a/)') REPEAT('=',67) + WRITE(stdout, '(/5x,a)') REPEAT('=',67) + WRITE(stdout, '(5x,"Scattering rate")') + 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)' ) '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 ',& - (ef+fsthick) * ryd2ev, ' eV will be included' - WRITE(stdout,'(5x,a/)') - ! - !IF (.NOT. ALLOCATED (inv_tau_all) ) ALLOCATE(inv_tau_all(nstemp,nbndfst,nktotf) ) - !inv_tau_all(:, :, :) = zero + WRITE(stdout, '(5x,a,f10.6,a,f10.6,a)' ) 'Only states between ',(ef - fsthick) * ryd2ev, ' eV and ', & + (ef + fsthick) * ryd2ev, ' eV will be included' + WRITE(stdout, '(5x,a/)') ! ENDIF ! @@ -167,11 +162,10 @@ ! etemp = transp_temp(itemp) ! - ! SP: Define the inverse so that we can efficiently multiply instead of - ! dividing + ! SP: Define the inverse so that we can efficiently multiply instead of dividing ! - inv_etemp = 1.0/etemp - inv_degaussw = 1.0/degaussw + inv_etemp = 1.0 / etemp + inv_degaussw = 1.0 / degaussw ! DO ik = 1, nkf ! @@ -185,32 +179,32 @@ DO ibnd = 1, nbndfst ! ! vkk(3,nbnd) - velocity for k - vkk(:,ibnd) = REAL (vmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = REAL(vmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ! DO jbnd = 1, nbndfst ! ! vkq(3,nbnd) - velocity for k + q - vkq(:,jbnd) = REAL (vmef (:, ibndmin-1+jbnd, ibndmin-1+jbnd, ikq ) ) + vkq(:, jbnd) = REAL(vmef(:, ibndmin - 1 + jbnd, ibndmin - 1 + jbnd, ikq)) ! - IF (ABS(vkk(1,ibnd)**2 + vkk(2,ibnd)**2 + vkk(3,ibnd)**2 ) > eps) & - vel_factor(ibnd,jbnd) = DDOT(3, vkk(:,ibnd), 1, vkq(:,jbnd), 1) / & - DDOT(3, vkk(:,ibnd), 1, vkk(:,ibnd), 1) + 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) ENDDO ENDDO ELSE DO ibnd = 1, nbndfst ! ! vkk(3,nbnd) - velocity for k - vkk(:,ibnd) = 2.0 * REAL (dmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = 2.0 * REAL(dmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ! DO jbnd = 1, nbndfst ! ! vkq(3,nbnd) - velocity for k + q - vkq(:,jbnd) = 2.0 * REAL (dmef (:, ibndmin-1+jbnd, ibndmin-1+jbnd, ikq ) ) + vkq(:, jbnd) = 2.0 * REAL(dmef(:, ibndmin - 1 + jbnd, ibndmin - 1 + jbnd, ikq)) ! - IF (ABS(vkk(1,ibnd)**2 + vkk(2,ibnd)**2 + vkk(3,ibnd)**2 ) > eps) & - vel_factor(ibnd,jbnd) = DDOT(3, vkk(:,ibnd), 1, vkq(:,jbnd), 1) / & - DDOT(3, vkk(:,ibnd), 1, vkk(:,ibnd), 1) + 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) ENDDO ENDDO ENDIF @@ -219,24 +213,22 @@ ! ! We are not consistent with ef from ephwann_shuffle but it should not ! matter if fstick is large enough. - !IF (( minval ( ABS(etf (:, ikk) - ef0(itemp)) ) < fsthick ) .AND. & - ! ( minval ( ABS(etf (:, ikq) - ef0(itemp)) ) < fsthick )) THEN - IF (( minval ( ABS(etf (:, ikk) - ef) ) < fsthick ) .AND. & - ( minval ( ABS(etf (:, ikq) - ef) ) < fsthick )) THEN + IF ((MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) .AND. & + (MINVAL(ABS(etf(:, ikq) - ef)) < fsthick)) THEN ! DO imode = 1, nmodes ! ! the phonon frequency and bose occupation - wq = wf (imode, iq) + wq = wf(imode, iq) ! ! SP : Avoid if statement in inner loops ! the coupling from Gamma acoustic phonons is negligible IF (wq > eps_acustic) THEN g2_tmp = 1.0 - wgq = wgauss( -wq*inv_etemp, -99) - wgq = wgq / ( one - two * wgq ) + wgq = wgauss(-wq * inv_etemp, -99) + wgq = wgq / (one - two * wgq) ! SP : Define the inverse for efficiency - inv_wq = 1.0/(two * wq) + inv_wq = 1.0 / (two * wq) ELSE g2_tmp = 0.0 wgq = 0.0 @@ -246,13 +238,13 @@ DO ibnd = 1, nbndfst ! ! energy at k (relative to Ef) - ekk = etf (ibndmin-1+ibnd, ikk) - ef0(itemp) + ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0(itemp) ! DO jbnd = 1, nbndfst ! ! energy and fermi occupation at k+q - ekq = etf (ibndmin-1+jbnd, ikq) - ef0(itemp) - fmkq = wgauss( -ekq*inv_etemp, -99) + ekq = etf(ibndmin - 1 + jbnd, ikq) - ef0(itemp) + fmkq = wgauss(-ekq * inv_etemp, -99) ! ! here we take into account the zero-point SQRT(hbar/2M\omega) ! with hbar = 1 and M already contained in the eigenmodes @@ -260,18 +252,18 @@ ! ! In case of q=\Gamma, then the short-range = the normal g. We therefore ! need to treat it like the normal g with ABS(g). - IF (shortrange .AND. ( ABS(xqf (1, iq))> eps8 .OR. ABS(xqf (2, iq))> eps8 & - .OR. ABS(xqf (3, iq))> eps8 )) THEN + 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*g2_tmp, KIND = DP ) + g2 = REAL((epf17(jbnd, ibnd, imode, ik)**two) * inv_wq * g2_tmp, KIND = DP) ELSE - g2 = (abs(epf17 (jbnd, ibnd, imode, ik))**two)*inv_wq*g2_tmp + g2 = (ABS(epf17(jbnd, ibnd, imode, ik))**two) * inv_wq * g2_tmp ENDIF ! ! delta[E_k - E_k+q + w_q] and delta[E_k - E_k+q - w_q] - w0g1 = w0gauss( (ekk-ekq+wq) * inv_degaussw, 0) * inv_degaussw - w0g2 = w0gauss( (ekk-ekq-wq) * inv_degaussw, 0) * inv_degaussw + w0g1 = w0gauss((ekk - ekq + wq) * inv_degaussw, 0) * inv_degaussw + w0g2 = w0gauss((ekk - ekq - wq) * inv_degaussw, 0) * inv_degaussw ! ! transition probability ! (2 pi/hbar) * (k+q-point weight) * g2 * @@ -279,42 +271,27 @@ ! [1 - f(E_k+q) + n(w_q)] * delta[E_k - E_k+q - w_q] } ! ! DBSP Just to try - !trans_prob = pi * g2 * & - ! ( (fmkq+wgq)*w0g1 + (one-fmkq+wgq)*w0g2 ) trans_prob = pi * wqf(iq) * g2 * & - ( (fmkq+wgq)*w0g1 + (one-fmkq+wgq)*w0g2 ) + ((fmkq + wgq) * w0g1 + (one - fmkq + wgq) * w0g2) ! - !if ((ik == 6) .AND. (ibnd == 4) .AND. jbnd==1 .AND. imode==6) then - ! print*,'wqf(iq) ',wqf(iq) - ! print*,'wq ',wq - ! print*,'inv_etemp ',inv_etemp - ! print*,'g2 ',g2 - ! print*,'fmkq ',fmkq - ! print*,'wgq ',wgq - ! print*,'w0g1 ',w0g1 - ! print*,'trans_prob ik, ibnd, jbnd, imode ', trans_prob, ik, ibnd, jbnd, imode - ! print*,'inv_tau',SUM(inv_tau_all(1,5:8,27)) - !end if - ! IF (scattering_serta) THEN ! energy relaxation time approximation - inv_tau_all(itemp,ibnd,ik+lower_bnd-1) = inv_tau_all(itemp,ibnd,ik+lower_bnd-1) + two * trans_prob - + inv_tau_all(itemp, ibnd, ik + lower_bnd - 1) = inv_tau_all(itemp, ibnd, ik + lower_bnd - 1) + two * trans_prob ELSEIF (scattering_0rta) THEN ! momentum relaxation time approximation - inv_tau_all(itemp,ibnd,ik+lower_bnd-1) = inv_tau_all(itemp,ibnd,ik+lower_bnd-1) & - + two * trans_prob * vel_factor(ibnd,jbnd) + inv_tau_all(itemp, ibnd, ik + lower_bnd - 1) = inv_tau_all(itemp, ibnd, ik + lower_bnd - 1) & + + two * trans_prob * vel_factor(ibnd, jbnd) ENDIF ! ! Z FACTOR: -\frac{\partial\Re\Sigma}{\partial\omega} ! weight = wqf(iq) * & - ( ( fmkq + wgq ) * ( (ekk - ( ekq - wq ))**two - degaussw**two ) / & - ( (ekk - ( ekq - wq ))**two + degaussw**two )**two + & - ( one - fmkq + wgq ) * ( (ekk - ( ekq + wq ))**two - degaussw**two ) / & - ( (ekk - ( ekq + wq ))**two + degaussw**two )**two ) + (( fmkq + wgq ) * ((ekk - (ekq - wq))**two - degaussw**two) / & + ((ekk - (ekq - wq))**two + degaussw**two)**two + & + (one - fmkq + wgq ) * ((ekk - (ekq + wq))**two - degaussw**two) / & + ((ekk - (ekq + wq))**two + degaussw**two)**two) ! - zi_allvb(itemp,ibnd,ik+lower_bnd-1) = zi_allvb(itemp,ibnd,ik+lower_bnd-1) + g2 * weight + zi_allvb(itemp, ibnd, ik + lower_bnd - 1) = zi_allvb(itemp, ibnd, ik + lower_bnd - 1) + g2 * weight ! ENDDO !jbnd ! @@ -322,18 +299,18 @@ ! ! In this case we are also computing the scattering rate for another Fermi level position ! This is used to compute both the electron and hole mobility at the same time. - IF (ABS(efcb(itemp)) > eps) THEN + IF (ABS(efcb(itemp)) > eps4) THEN ! DO ibnd = 1, nbndfst ! ! energy at k (relative to Ef) - ekk = etf (ibndmin-1+ibnd, ikk) - efcb(itemp) + ekk = etf(ibndmin - 1 + ibnd, ikk) - efcb(itemp) ! DO jbnd = 1, nbndfst ! ! energy and fermi occupation at k+q - ekq = etf (ibndmin-1+jbnd, ikq) - efcb(itemp) - fmkq = wgauss( -ekq*inv_etemp, -99) + ekq = etf(ibndmin - 1 + jbnd, ikq) - efcb(itemp) + fmkq = wgauss(-ekq * inv_etemp, -99) ! ! here we take into account the zero-point SQRT(hbar/2M\omega) ! with hbar = 1 and M already contained in the eigenmodes @@ -341,18 +318,18 @@ ! ! In case of q=\Gamma, then the short-range = the normal g. We therefore ! need to treat it like the normal g with ABS(g). - IF (shortrange .AND. ( ABS(xqf (1, iq))> eps8 .OR. ABS(xqf (2, iq))> eps8 & - .OR. ABS(xqf (3, iq))> eps8 )) THEN + 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*g2_tmp, KIND = DP) + g2 = REAL((epf17(jbnd, ibnd, imode, ik)**two) * inv_wq * g2_tmp, KIND = DP) ELSE - g2 = (abs(epf17 (jbnd, ibnd, imode, ik))**two)*inv_wq*g2_tmp + g2 = (ABS(epf17(jbnd, ibnd, imode, ik))**two) * inv_wq * g2_tmp ENDIF ! ! delta[E_k - E_k+q + w_q] and delta[E_k - E_k+q - w_q] - w0g1 = w0gauss( (ekk-ekq+wq) * inv_degaussw, 0) * inv_degaussw - w0g2 = w0gauss( (ekk-ekq-wq) * inv_degaussw, 0) * inv_degaussw + w0g1 = w0gauss((ekk - ekq + wq) * inv_degaussw, 0) * inv_degaussw + w0g2 = w0gauss((ekk - ekq - wq) * inv_degaussw, 0) * inv_degaussw ! ! transition probability ! (2 pi/hbar) * (k+q-point weight) * g2 * @@ -360,71 +337,63 @@ ! [1 - f(E_k+q) + n(w_q)] * delta[E_k - E_k+q - w_q] } ! trans_prob = pi * wqf(iq) * g2 * & - ( (fmkq+wgq)*w0g1 + (one-fmkq+wgq)*w0g2 ) + ((fmkq + wgq) * w0g1 + (one - fmkq + wgq) * w0g2) ! IF (scattering_serta) THEN ! energy relaxation time approximation - inv_tau_allcb(itemp,ibnd,ik+lower_bnd-1) = inv_tau_allcb(itemp,ibnd,ik+lower_bnd-1) + two * trans_prob + inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1) = inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1) + two * trans_prob ! ELSEIF (scattering_0rta) THEN ! momentum relaxation time approximation - inv_tau_allcb(itemp,ibnd,ik+lower_bnd-1) = inv_tau_allcb(itemp,ibnd,ik+lower_bnd-1) & - + two * trans_prob * vel_factor(ibnd,jbnd) + inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1) = inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1) & + + two * trans_prob * vel_factor(ibnd, jbnd) ENDIF ! ! Z FACTOR: -\frac{\partial\Re\Sigma}{\partial\omega} ! weight = wqf(iq) * & - ( ( fmkq + wgq ) * ( (ekk - ( ekq - wq ))**two - degaussw**two ) / & - ( (ekk - ( ekq - wq ))**two + degaussw**two )**two + & - ( one - fmkq + wgq ) * ( (ekk - ( ekq + wq ))**two - degaussw**two ) / & - ( (ekk - ( ekq + wq ))**two + degaussw**two )**two ) + (( fmkq + wgq) * ((ekk - (ekq - wq))**two - degaussw**two) / & + ((ekk - (ekq - wq))**two + degaussw**two)**two + & + (one - fmkq + wgq) * ((ekk - (ekq + wq))**two - degaussw**two) / & + ((ekk - (ekq + wq))**two + degaussw**two)**two) ! - zi_allcb(itemp,ibnd,ik+lower_bnd-1) = zi_allcb(itemp,ibnd,ik+lower_bnd-1) + g2 * weight + zi_allcb(itemp, ibnd, ik + lower_bnd - 1) = zi_allcb(itemp, ibnd, ik + lower_bnd - 1) + g2 * weight ! ENDDO !jbnd - ! ENDDO !ibnd - ! - ENDIF ! ABS(efcb) < eps - ! + ENDIF ! ABS(efcb) < eps4 ENDDO !imode - ! ENDIF ! endif fsthick - ! ENDDO ! end loop on k ENDDO ! itemp ! ! Creation of a restart point IF (restart) THEN - IF (MOD(iqq,restart_freq) == 0) THEN + IF (MOD(iqq, restart_freq) == 0) THEN WRITE(stdout, '(a)' ) ' Creation of a restart point' ! ! The mp_sum will aggreage the results on each k-points. - CALL mp_sum( inv_tau_all, world_comm ) - CALL mp_sum( zi_allvb, world_comm ) + CALL mp_sum(inv_tau_all, world_comm) + CALL mp_sum(zi_allvb, world_comm) ! - IF (ABS(efcb(1)) > eps) THEN + IF (ABS(efcb(1)) > eps4) THEN ! - CALL mp_sum( inv_tau_allcb, world_comm ) - CALL mp_sum( zi_allcb, world_comm ) + CALL mp_sum(inv_tau_allcb, world_comm) + CALL mp_sum(zi_allcb, world_comm) ! ENDIF ! - IF (ABS(efcb(1)) > eps) THEN - CALL tau_write(iqq,totq,nktotf,.TRUE.) + IF (ABS(efcb(1)) > eps4) THEN + CALL tau_write(iqq, totq, nktotf, .TRUE.) ELSE - CALL tau_write(iqq,totq,nktotf,.FALSE.) + CALL tau_write(iqq, totq, nktotf, .FALSE.) ENDIF ! ! Now show intermediate mobility with that amount of q-points - CALL transport_coeffs(ef0,efcb) - ! + CALL transport_coeffs(ef0, efcb) ENDIF ENDIF - ! ENDIF ! first_cycle - ! ! ! The k points are distributed among pools: here we collect them ! @@ -432,21 +401,19 @@ ! ! The total number of k points ! - ALLOCATE(etf_all ( nbndsub, nkqtotf )) + ALLOCATE(etf_all(nbndsub, nkqtotf)) ! - CALL mp_sum( inv_tau_all, world_comm ) - IF (ABS(efcb(1)) > eps) CALL mp_sum( inv_tau_allcb, world_comm ) - CALL mp_sum( zi_allvb, world_comm ) - IF (ABS(efcb(1)) > eps) CALL mp_sum( zi_allcb, world_comm ) - !print*,'zi_allvb SUM ',SUM(zi_allvb) - !print*,'inv_tau_all SUM ',SUM(inv_tau_all) + CALL mp_sum(inv_tau_all, world_comm) + IF (ABS(efcb(1)) > eps4) CALL mp_sum(inv_tau_allcb, world_comm) + CALL mp_sum(zi_allvb, world_comm) + IF (ABS(efcb(1)) > eps4) CALL mp_sum(zi_allcb, world_comm) ! #if defined(__MPI) ! - ! collect contributions from all pools (sum over k-points) + ! Collect contributions from all pools (sum over k-points) ! this finishes the integral over the BZ (k) ! - CALL poolgather2 ( nbndsub, nkqtotf, nkqf, etf, etf_all ) + CALL poolgather2(nbndsub, nkqtotf, nkqf, etf, etf_all) #else ! etf_all = etf @@ -455,7 +422,7 @@ DO itemp = 1, nstemp ! etemp = transp_temp(itemp) - WRITE(stdout, '(a,f8.3,a)' ) ' Temperature ',etemp * ryd2ev / kelvin2eV,' K' + WRITE(stdout, '(a,f8.3,a)' ) ' Temperature ', etemp * ryd2ev / kelvin2eV, ' K' ! ! In case we read another q-file, merge the scattering here IF (restart_filq /= '') THEN @@ -463,39 +430,39 @@ ALLOCATE(inv_tau_all_new(nstemp, nbndfst, nktotf)) inv_tau_all_new(:, :, :) = zero ! - CALL merge_read( nktotf, nqtotf_new, inv_tau_all_new ) + 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) + inv_tau_all(:, :, :) = (inv_tau_all(:, :, :) * totq & + + inv_tau_all_new(:, :, :) * nqtotf_new) / (totq + nqtotf_new) DEALLOCATE(inv_tau_all_new) ! WRITE(stdout, '(a)' ) ' ' - WRITE(stdout, '(a,i10,a)' ) ' Merge scattering for a total of ',totq+nqtotf_new,' q-points' + WRITE(stdout, '(a,i10,a)' ) ' Merge scattering for a total of ', totq + nqtotf_new, ' q-points' ! - CALL tau_write(iqq+nqtotf_new,totq+nqtotf_new,nktotf, .FALSE.) + CALL tau_write(iqq + nqtotf_new, totq + nqtotf_new, nktotf, .FALSE.) WRITE(stdout, '(a)' ) ' Write to restart file the sum' WRITE(stdout, '(a)' ) ' ' ! ! ENDIF ! Average over degenerate eigenstates: - WRITE(stdout,'(5x,"Average over degenerate eigenstates is performed")') + WRITE(stdout, '(5x,"Average over degenerate eigenstates is performed")') ! DO ik = 1, nktotf ikk = 2 * ik - 1 ikq = ikk + 1 ! DO ibnd = 1, nbndfst - ekk = etf_all (ibndmin-1+ibnd, ikk) + ekk = etf_all(ibndmin - 1 + ibnd, ikk) n = 0 tmp = 0.0_DP tmp2 = 0.0_DP DO jbnd = 1, nbndfst - ekk2 = etf_all (ibndmin-1+jbnd, ikk) - IF (ABS(ekk2-ekk) < eps6) THEN + ekk2 = etf_all(ibndmin - 1 + jbnd, ikk) + IF (ABS(ekk2 - ekk) < eps6) THEN n = n + 1 - tmp = tmp + inv_tau_all (itemp,jbnd,ik) - tmp2 = tmp2 + zi_allvb(itemp,jbnd,ik) + tmp = tmp + inv_tau_all(itemp, jbnd, ik) + tmp2 = tmp2 + zi_allvb(itemp, jbnd, ik) ENDIF ! ENDDO ! jbnd @@ -503,30 +470,30 @@ zi_tmp(ibnd) = tmp2 / FLOAT(n) ! ENDDO ! ibnd - inv_tau_all (itemp,:,ik) = inv_tau_tmp(:) - zi_allvb (itemp,:,ik) = zi_tmp(:) + inv_tau_all(itemp, :, ik) = inv_tau_tmp(:) + zi_allvb(itemp, :, ik) = zi_tmp(:) ! ENDDO ! nkqtotf ! - IF (ABS(efcb(itemp)) > eps) THEN + IF (ABS(efcb(itemp)) > eps4) THEN ! Average over degenerate eigenstates: - WRITE(stdout,'(5x,"Average over degenerate eigenstates in CB is performed")') + WRITE(stdout, '(5x,"Average over degenerate eigenstates in CB is performed")') ! DO ik = 1, nktotf ikk = 2 * ik - 1 ikq = ikk + 1 ! DO ibnd = 1, nbndfst - ekk = etf_all (ibndmin-1+ibnd, ikk) + ekk = etf_all(ibndmin - 1 + ibnd, ikk) n = 0 tmp = 0.0_DP tmp2 = 0.0_DP DO jbnd = 1, nbndfst - ekk2 = etf_all (ibndmin-1+jbnd, ikk) - IF (ABS(ekk2-ekk) < eps6) THEN + ekk2 = etf_all(ibndmin - 1 + jbnd, ikk) + IF (ABS(ekk2 - ekk) < eps6) THEN n = n + 1 - tmp = tmp + inv_tau_allcb (itemp,jbnd,ik) - tmp2 = tmp2 + zi_allcb (itemp,jbnd,ik) + tmp = tmp + inv_tau_allcb(itemp, jbnd, ik) + tmp2 = tmp2 + zi_allcb(itemp, jbnd, ik) ENDIF ! ENDDO ! jbnd @@ -534,15 +501,14 @@ zi_tmp(ibnd) = tmp2 / FLOAT(n) ! ENDDO ! ibnd - inv_tau_allcb (itemp,:,ik) = inv_tau_tmp(:) - zi_allcb (itemp,:,ik) = zi_tmp(:) + inv_tau_allcb(itemp, :, ik) = inv_tau_tmp(:) + zi_allcb(itemp, :, ik) = zi_tmp(:) ! ENDDO ! nkqtotf ENDIF ! ! Output scattering rates here after looping over all q-points ! (with their contributions summed in inv_tau_all, etc.) - !print*,'inv_tau_all(1,1,1) ',inv_tau_all(1,1,1) CALL scattering_write(itemp, etemp, ef0, etf_all) ! ENDDO !nstemp @@ -553,17 +519,17 @@ IF (restart) THEN WRITE(stdout, '(a)' ) ' Creation of the final restart point' ! - IF (ABS(efcb(1)) > eps) THEN - CALL tau_write(iqq,totq,nktotf,.TRUE.) + IF (ABS(efcb(1)) > eps4) THEN + CALL tau_write(iqq, totq, nktotf, .TRUE.) ELSE - CALL tau_write(iqq,totq,nktotf,.FALSE.) + CALL tau_write(iqq, totq, nktotf, .FALSE.) ENDIF ! ENDIF ! restart ! ENDIF ! iqq ! - CALL stop_clock ( 'SCAT' ) + CALL stop_clock ('SCAT') ! DBSP !write(stdout,*),'iqq ',iqq !print*,shape(inv_tau_all) @@ -574,11 +540,12 @@ ! RETURN ! + !----------------------------------------------------------------------- END SUBROUTINE scattering_rate_q !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- - SUBROUTINE transport_coeffs (ef0, efcb) + SUBROUTINE transport_coeffs(ef0, efcb) !----------------------------------------------------------------------- !! !! This SUBROUTINE computes the transport coefficients @@ -600,7 +567,7 @@ USE pwcom, ONLY : ef USE elph2, ONLY : ibndmax, ibndmin, etf, nkf, wkf, dmef, vmef, & inv_tau_all, nkqtotf, inv_tau_allcb, & - zi_allvb, zi_allcb, map_rebal + zi_allvb, zi_allcb, map_rebal, nbndfst, nktotf USE transportcom, ONLY : transp_temp USE constants_epw, ONLY : zero, one, bohr2ang, ryd2ev, electron_SI, & kelvin2eV, hbar, Ang2m, hbarJ, ang2cm, czero @@ -614,7 +581,7 @@ USE mp, ONLY : mp_bcast USE mp_global, ONLY : inter_pool_comm USE epwcom, ONLY : mp_mesh_k, nkf1, nkf2, nkf3 - USE constants_epw, ONLY : eps6 + USE constants_epw, ONLY : eps6, eps4 USE noncollin_module, ONLY : noncolin USE io_scattering, ONLY : scattering_read USE division, ONLY : fkbounds @@ -645,24 +612,21 @@ !! Lower bounds index after k or q paral INTEGER :: upper_bnd !! Upper bounds index after k or q paral - ! SP - Uncomment to use symmetries on velocities INTEGER :: nkqtotf_tmp !! Temporary k-q points. INTEGER :: ikbz !! k-point index that run on the full BZ INTEGER :: nb !! Number of points in the BZ corresponding to a point in IBZ - INTEGER :: BZtoIBZ_tmp(nkf1*nkf2*nkf3) + INTEGER :: BZtoIBZ_tmp(nkf1 * nkf2 * nkf3) !! Temporary mapping - INTEGER :: BZtoIBZ(nkf1*nkf2*nkf3) + INTEGER :: BZtoIBZ(nkf1 * nkf2 * nkf3) !! BZ to IBZ mapping - INTEGER :: s_BZtoIBZ(3,3,nkf1*nkf2*nkf3) + INTEGER :: s_BZtoIBZ(3, 3, nkf1 * nkf2 * nkf3) !! symmetry ! REAL(KIND = DP) :: ekk !! Energy relative to Fermi level: $$\varepsilon_{n\mathbf{k}}-\varepsilon_F$$ - ! REAL(KIND = DP) :: ef0(nstemp) - !! Fermi level for the temperature itemp REAL(KIND = DP) :: dfnk !! Derivative Fermi distribution $$-df_{nk}/dE_{nk}$$ REAL(KIND = DP) :: etemp @@ -687,26 +651,24 @@ !! Mobility along the yy axis after diagonalization [cm^2/Vs] REAL(KIND = DP) :: mobility_zz !! Mobility along the zz axis after diagonalization [cm^2/Vs] - REAL(KIND = DP) :: vkk(3,nbndfst) + REAL(KIND = DP) :: vkk(3, nbndfst) !! Electron velocity vector for a band. - REAL(KIND = DP) :: Sigma(9,nstemp) + REAL(KIND = DP) :: Sigma(9, nstemp) !! Conductivity matrix in vector form - REAL(KIND = DP) :: SigmaZ(9,nstemp) + REAL(KIND = DP) :: SigmaZ(9, nstemp) !! Conductivity matrix in vector form with Znk - REAL(KIND = DP) :: Sigma_m(3,3,nstemp) + REAL(KIND = DP) :: Sigma_m(3, 3, nstemp) !! Conductivity matrix - REAL(KIND = DP) :: sigma_up(3,3) + REAL(KIND = DP) :: sigma_up(3, 3) !! Conductivity matrix in upper-triangle REAL(KIND = DP) :: sigma_eig(3) !! Eigenvalues from the diagonalized conductivity matrix - REAL(KIND = DP) :: sigma_vect(3,3) + REAL(KIND = DP) :: sigma_vect(3, 3) !! Eigenvectors from the diagonalized conductivity matrix REAL(KIND = DP) :: Znk !! Real Znk from \lambda_nk (called zi_allvb or zi_allcb) REAL(KIND = DP) :: tdf_sigma(9) !! Temporary file - REAL(KIND = DP), PARAMETER :: eps = 1.d-4 - !! Tolerence REAL(KIND = DP), EXTERNAL :: wgauss !! Compute the approximate theta function. Here computes Fermi-Dirac REAL(KIND = DP), EXTERNAL :: w0gauss @@ -734,18 +696,23 @@ !! Rotated velocity by the symmetry operation REAL(KIND = DP) :: vk_cart(3) !! veloctiy in cartesian coordinate - REAL(KIND = DP) :: sa(3,3), sb(3,3), sr(3,3) + REAL(KIND = DP) :: sa(3, 3) + !! Rotation matrix + REAL(KIND = DP) :: sb(3, 3) + !! Rotation matrix + REAL(KIND = DP) :: sr(3, 3) + !! Rotation matrix + ! CALL start_clock('MOB') ! - inv_cell = 1.0d0/omega + 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 + IF (system_2d ) inv_cell = inv_cell * at(3, 3) * alat ! ! We can read the scattering rate from files. IF (scatread) THEN - conv_factor1 = electron_SI / ( hbar * bohr2ang * Ang2m ) + conv_factor1 = electron_SI / (hbar * bohr2ang * Ang2m) Sigma_m(:, :, :) = zero ! ! Compute the Fermi level @@ -788,37 +755,37 @@ ! This is hole mobility. ---------------------------------------------------- IF (int_mob .OR. (ncarrier < -1E5)) THEN IF (itemp == 1) THEN - WRITE(stdout,'(/5x,a)') REPEAT('=',67) - WRITE(stdout,'(5x,"Temp [K] Fermi [eV] Hole density [cm^-3] Hole mobility [cm^2/Vs]")') - WRITE(stdout,'(5x,a/)') REPEAT('=',67) + WRITE(stdout, '(/5x,a)') REPEAT('=',67) + WRITE(stdout, '(5x,"Temp [K] Fermi [eV] Hole density [cm^-3] Hole mobility [cm^2/Vs]")') + WRITE(stdout, '(5x,a/)') REPEAT('=',67) ENDIF ! DO ik = 1, nktotf - ikk=2 * ik - 1 + ikk = 2 * ik - 1 ! here we must have ef, not ef0, to be consistent with ephwann_shuffle - IF (minval ( ABS(etf_all (:, ik) - ef ) ) < fsthick) THEN + IF (MINVAL(ABS(etf_all(:, ik) - ef)) < fsthick) THEN DO ibnd = 1, nbndfst ! This selects only valence bands for hole conduction - IF (etf_all (ibndmin-1+ibnd, ik) < ef0(itemp) ) THEN + IF (etf_all(ibndmin - 1 + ibnd, ik) < ef0(itemp)) THEN IF (vme) THEN - vkk(:,ibnd) = REAL (vmef_all (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = REAL(vmef_all(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ELSE - vkk(:,ibnd) = 2.0 * REAL (dmef_all (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = 2.0 * REAL(dmef_all (:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ENDIF - tau = one / inv_tau_all(itemp,ibnd,ik) - ekk = etf_all (ibndmin-1+ibnd, ik) - ef0(itemp) + tau = one / inv_tau_all(itemp, ibnd, ik) + ekk = etf_all(ibndmin - 1 + ibnd, ik) - ef0(itemp) ! DO j = 1, 3 DO i = 1, 3 - tdf_sigma_m(i,j,ibnd,ik) = vkk(i,ibnd) * vkk(j,ibnd) * tau + tdf_sigma_m(i, j, ibnd, ik) = vkk(i, ibnd) * vkk(j, ibnd) * tau ENDDO ENDDO ! ! derivative Fermi distribution - dfnk = w0gauss( ekk / etemp, -99 ) / etemp + dfnk = w0gauss(ekk / etemp, -99) / etemp ! ! electrical conductivity matrix - Sigma_m(:,:,itemp) = Sigma_m(:,:,itemp) + wkf_all(ikk) * dfnk * tdf_sigma_m(:,:,ibnd,ik) + Sigma_m(:, :, itemp) = Sigma_m(:, :, itemp) + wkf_all(ikk) * dfnk * tdf_sigma_m(:, :, ibnd, ik) ! ENDIF ! valence bands ENDDO ! ibnd @@ -831,39 +798,39 @@ ikk = 2 * ik - 1 DO ibnd = 1, nbndfst ! This selects only valence bands for hole conduction - IF (etf_all (ibndmin-1+ibnd, ik) < ef0(itemp) ) THEN + IF (etf_all(ibndmin - 1 + ibnd, ik) < ef0(itemp)) THEN ! energy at k (relative to Ef) - ekk = etf_all (ibndmin-1+ibnd, ik) - ef0(itemp) - fnk = wgauss( -ekk / etemp, -99) + ekk = etf_all(ibndmin - 1 + ibnd, ik) - ef0(itemp) + fnk = wgauss(-ekk / etemp, -99) ! The wkf(ikk) already include a factor 2 - carrier_density = carrier_density + wkf_all(ikk) * (1.0d0 - fnk ) + carrier_density = carrier_density + wkf_all(ikk) * (1.0d0 - fnk) ENDIF ENDDO ENDDO ! ! Diagonalize the conductivity matrix - CALL rdiagh(3,Sigma_m(:,:,itemp),3,sigma_eig,sigma_vect) + CALL rdiagh(3, Sigma_m(:, :, itemp), 3, sigma_eig, sigma_vect) ! - mobility_xx = ( sigma_eig(1) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility_yy = ( sigma_eig(2) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility_zz = ( sigma_eig(3) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility = (mobility_xx+mobility_yy+mobility_zz)/3 + mobility_xx = (sigma_eig(1) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_yy = (sigma_eig(2) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_zz = (sigma_eig(3) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility = (mobility_xx + mobility_yy + mobility_zz) / 3 ! carrier_density in cm^-1 - carrier_density_prt = carrier_density * inv_cell * ( bohr2ang * ang2cm )**(-3) - WRITE(stdout,'(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev /kelvin2eV, & - ef0(itemp)*ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility_yy, ' y-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility_zz, ' z-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility, ' avg' + carrier_density_prt = carrier_density * inv_cell * (bohr2ang * ang2cm)**(-3) + WRITE(stdout, '(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev / kelvin2eV, & + ef0(itemp) * ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility_yy, ' y-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility_zz, ' z-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility, ' avg' ! ENDIF ! int_mob .OR. (ncarrier < -1E5) ! ! This is electron mobility. ---------------------------------------------------- IF (int_mob .OR. (ncarrier > 1E5)) THEN IF (itemp == 1) THEN - WRITE(stdout,'(/5x,a)') REPEAT('=',67) - WRITE(stdout,'(5x,"Temp [K] Fermi [eV] Electron density [cm^-3] Electron mobility [cm^2/Vs]")') - WRITE(stdout,'(5x,a/)') REPEAT('=',67) + WRITE(stdout, '(/5x,a)') REPEAT('=',67) + WRITE(stdout, '(5x,"Temp [K] Fermi [eV] Electron density [cm^-3] Electron mobility [cm^2/Vs]")') + WRITE(stdout, '(5x,a/)') REPEAT('=',67) ENDIF ! tdf_sigma_m(:, :, :, :) = zero @@ -871,29 +838,29 @@ DO ik = 1, nktotf ikk = 2 * ik - 1 ! here we must have ef, not ef0, to be consistent with ephwann_shuffle - IF (minval ( ABS(etf_all (:, ik) - ef ) ) < fsthick) THEN + IF (MINVAL(ABS(etf_all(:, ik) - ef)) < fsthick) THEN DO ibnd = 1, nbndfst ! This selects only conduction bands for electron conduction - IF (etf_all (ibndmin-1+ibnd, ik) > ef0(itemp) ) THEN + IF (etf_all(ibndmin - 1 + ibnd, ik) > ef0(itemp)) THEN IF (vme) THEN - vkk(:,ibnd) = REAL (vmef_all (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = REAL(vmef_all(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ELSE - vkk(:,ibnd) = 2.0 * REAL (dmef_all (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = 2.0 * REAL(dmef_all(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ENDIF - tau = one / inv_tau_all(itemp,ibnd,ik) - ekk = etf_all (ibndmin-1+ibnd, ik) - ef0(itemp) + tau = one / inv_tau_all(itemp, ibnd, ik) + ekk = etf_all(ibndmin - 1 + ibnd, ik) - ef0(itemp) ! DO j = 1, 3 DO i = 1, 3 - tdf_sigma_m(i,j,ibnd,ik) = vkk(i,ibnd) * vkk(j,ibnd) * tau + tdf_sigma_m(i, j, ibnd, ik) = vkk(i, ibnd) * vkk(j, ibnd) * tau ENDDO ENDDO ! ! derivative Fermi distribution - dfnk = w0gauss( ekk / etemp, -99 ) / etemp + dfnk = w0gauss(ekk / etemp, -99) / etemp ! ! electrical conductivity matrix - Sigma_m(:,:,itemp) = Sigma_m(:,:,itemp) + wkf_all(ikk) * dfnk * tdf_sigma_m(:,:,ibnd,ik) + Sigma_m(:, :, itemp) = Sigma_m(:, :, itemp) + wkf_all(ikk) * dfnk * tdf_sigma_m(:, :, ibnd, ik) ! ENDIF ! valence bands ENDDO ! ibnd @@ -906,10 +873,10 @@ ikk = 2 * ik - 1 DO ibnd = 1, nbndfst ! This selects only conduction bands for electron conduction - IF (etf_all (ibndmin-1+ibnd, ik) > ef0(itemp) ) THEN + IF (etf_all(ibndmin - 1 + ibnd, ik) > ef0(itemp)) THEN ! energy at k (relative to Ef) - ekk = etf_all (ibndmin-1+ibnd, ik) - ef0(itemp) - fnk = wgauss( -ekk / etemp, -99) + ekk = etf_all(ibndmin - 1 + ibnd, ik) - ef0(itemp) + fnk = wgauss(-ekk / etemp, -99) ! The wkf(ikk) already include a factor 2 carrier_density = carrier_density + wkf_all(ikk) * fnk ENDIF @@ -917,19 +884,19 @@ ENDDO ! ! Diagonalize the conductivity matrix - CALL rdiagh(3,Sigma_m(:,:,itemp),3,sigma_eig,sigma_vect) + CALL rdiagh(3, Sigma_m(:, :, itemp), 3, sigma_eig, sigma_vect) ! - mobility_xx = ( sigma_eig(1) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility_yy = ( sigma_eig(2) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility_zz = ( sigma_eig(3) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility = (mobility_xx+mobility_yy+mobility_zz)/3 + mobility_xx = (sigma_eig(1) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_yy = (sigma_eig(2) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_zz = (sigma_eig(3) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility = (mobility_xx + mobility_yy + mobility_zz) / 3 ! carrier_density in cm^-1 - carrier_density_prt = carrier_density * inv_cell * ( bohr2ang * ang2cm )**(-3) - WRITE(stdout,'(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev /kelvin2eV, & - ef0(itemp)*ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility_yy, ' y-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility_zz, ' z-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility, ' avg' + carrier_density_prt = carrier_density * inv_cell * (bohr2ang * ang2cm)**(-3) + WRITE(stdout, '(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev / kelvin2eV, & + ef0(itemp) * ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility_yy, ' y-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility_zz, ' z-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility, ' avg' ! ENDIF ! int_mob .OR. (ncarrier > 1E5) ! @@ -946,33 +913,23 @@ ! ! SP - Uncomment to use symmetries on velocities IF (mp_mesh_k) THEN - IF (mpime == ionode_id) THEN - ! - CALL set_sym_bl( ) - BZtoIBZ(:) = 0 - s_BZtoIBZ(:, :, :) = 0 - ! What we get from this call is BZtoIBZ - CALL kpoint_grid_epw ( nrot, time_reversal, .FALSE., s, t_rev, bg, nkf1*nkf2*nkf3, & - nkf1,nkf2,nkf3, nkqtotf_tmp, xkf_tmp, wkf_tmp,BZtoIBZ,s_BZtoIBZ) - ! - IF (iterative_bte) THEN - ! Now we have to remap the points because the IBZ k-points have been - ! changed to to load balancing. - BZtoIBZ_tmp(:) = 0 - DO ikbz = 1, nkf1*nkf2*nkf3 - BZtoIBZ_tmp(ikbz) = map_rebal( BZtoIBZ( ikbz ) ) - ENDDO - BZtoIBZ(:) = BZtoIBZ_tmp(:) - ENDIF - ! + BZtoIBZ(:) = 0 + s_BZtoIBZ(:) = 0 + ! + CALL set_sym_bl() + ! What we get from this call is BZtoIBZ + CALL kpoint_grid_epw(nrot, time_reversal, .FALSE., s, t_rev, bg, nkf1, nkf2, nkf3, BZtoIBZ, s_BZtoIBZ) + ! + IF (iterative_bte) THEN + ! Now we have to remap the points because the IBZ k-points have been + ! changed to to load balancing. + BZtoIBZ_tmp(:) = 0 + DO ikbz = 1, nkf1 * nkf2 * nkf3 + BZtoIBZ_tmp(ikbz) = map_rebal(BZtoIBZ(ikbz)) + ENDDO + BZtoIBZ(:) = BZtoIBZ_tmp(:) ENDIF - CALL mp_bcast( s_BZtoIBZ, ionode_id, inter_pool_comm ) - CALL mp_bcast( BZtoIBZ, ionode_id, inter_pool_comm ) - ! - !WRITE(stdout,*)'s ',s - !WRITE(stdout,*)'BZtoIBZ ',BZtoIBZ - !WRITE(stdout,*)'s_BZtoIBZ ',s_BZtoIBZ - ! + ! ENDIF ! ! This is hole mobility. In the case of intrinsic mobilities we can do both @@ -980,17 +937,13 @@ ! the case for doped mobilities. ! ! find the bounds of k-dependent arrays in the parallel case in each pool - CALL fkbounds( nktotf, lower_bnd, upper_bnd ) + CALL fkbounds(nktotf, lower_bnd, upper_bnd) ! IF (int_mob .OR. (ncarrier < -1E5)) THEN ! DO itemp = 1, nstemp ! etemp = transp_temp(itemp) - !DBSP - !write(stdout,*)'etemp ',etemp - !write(stdout,*)'inv_tau_all ', SUM(inv_tau_all(itemp,:,:)) - !write(stdout,*)'inv_tau_all ', SUM(inv_tau_all(:, :, :)) ! IF (itemp == 1) THEN ! @@ -1000,9 +953,9 @@ ! 4 = (2,1) = yx, 5 = (2,2) = yy, 6 = (2,3) = yz ! 7 = (3,1) = zx, 8 = (3,2) = zy, 9 = (3,3) = zz ! this can be reduced to 6 if we take into account symmetry xy=yx, ... - tdf_sigma(:) = zero - Sigma(:, :) = zero - SigmaZ(:, :) = zero + tdf_sigma(:) = zero + Sigma(:, :) = zero + SigmaZ(:, :) = zero ! ENDIF ! @@ -1011,66 +964,63 @@ ikk = 2 * ik - 1 ! ! here we must have ef, not ef0, to be consistent with ephwann_shuffle - IF (minval ( ABS(etf (:, ikk) - ef) ) < fsthick) THEN + IF (MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) THEN ! DO ibnd = 1, nbndfst ! ! This selects only valence bands for hole conduction - IF (etf (ibndmin-1+ibnd, ikk) < ef0(itemp)) THEN + IF (etf(ibndmin - 1 + ibnd, ikk) < ef0(itemp)) THEN ! ! vkk(3,nbnd) - velocity for k tdf_sigma(:) = zero IF (vme) THEN ! vmef is in units of Ryd * bohr - vkk(:,ibnd) = REAL (vmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = REAL(vmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ELSE ! v_(k,i) = 1/m = 2 * dmef (:, i,i,k) ! 1/m = 2 in Rydberg atomic units ! dmef is in units of 1/a.u. (where a.u. is bohr) ! v_(k,i) is in units of Ryd * a.u. - vkk(:,ibnd) = 2.0 * REAL (dmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = 2.0 * REAL(dmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ENDIF ! Use symmetries on k-point (from Homogeneous grid only) IF (mp_mesh_k) THEN ! - vk_cart(:) = vkk(:,ibnd) + vk_cart(:) = vkk(:, ibnd) ! ! Loop on full BZ nb = 0 - DO ikbz = 1, nkf1*nkf2*nkf3 + DO ikbz = 1, nkf1 * nkf2 * nkf3 ! If the k-point from the full BZ is related by a symmetry operation ! to the current k-point, then take it. - IF (BZtoIBZ(ikbz) == ik+lower_bnd-1) THEN + IF (BZtoIBZ(ikbz) == ik + lower_bnd - 1) THEN nb = nb + 1 ! Transform the symmetry matrix from Crystal to cartesian - sa (:, :) = dble ( s_BZtoIBZ(:,:,ikbz) ) - sb = matmul ( bg, sa ) - sr (:, :) = matmul ( at, transpose (sb) ) - CALL dgemv( 'n', 3, 3, 1.d0,& - sr, 3, vk_cart(:),1 ,0.d0 , v_rot(:), 1 ) + sa(:, :) = DBLE(s_BZtoIBZ(:, :, ikbz)) + sb = MATMUL(bg, sa) + sr(:, :) = MATMUL(at, TRANSPOSE(sb)) + CALL dgemv('n', 3, 3, 1.d0, sr, 3, vk_cart(:), 1, 0.d0, v_rot(:), 1) ij = 0 DO j = 1, 3 DO i = 1, 3 ij = ij + 1 ! The factor two in the weight at the end is to account for spin IF (noncolin) THEN - tdf_sigma(ij) = tdf_sigma(ij) + ( v_rot(i) * v_rot(j) ) * 1.0 / (nkf1*nkf2*nkf3) + tdf_sigma(ij) = tdf_sigma(ij) + (v_rot(i) * v_rot(j)) * 1.0 / (nkf1 * nkf2 * nkf3) ELSE - tdf_sigma(ij) = tdf_sigma(ij) + ( v_rot(i) * v_rot(j) ) * 2.0 / (nkf1*nkf2*nkf3) + tdf_sigma(ij) = tdf_sigma(ij) + (v_rot(i) * v_rot(j)) * 2.0 / (nkf1 * nkf2 * nkf3) ENDIF ENDDO ENDDO ENDIF ENDDO ! ikbz IF (noncolin) THEN - IF (ABS(nb*1.0/(nkf1*nkf2*nkf3) - wkf(ikk)) > eps6) THEN - CALL errore ('transport', & - &' The number of kpoint in the IBZ is not equal to the weight', 1) + IF (ABS(nb * 1.0 / (nkf1 * nkf2 * nkf3) - wkf(ikk)) > eps6) THEN + CALL errore('transport', ' The number of kpoint in the IBZ is not equal to the weight', 1) ENDIF ELSE - IF (ABS(nb*2.0/(nkf1*nkf2*nkf3) - wkf(ikk)) > eps6) THEN - CALL errore ('transport', & - &' The number of kpoint in the IBZ is not equal to the weight', 1) + IF (ABS(nb * 2.0 / (nkf1 * nkf2 * nkf3) - wkf(ikk)) > eps6) THEN + CALL errore('transport', ' The number of kpoint in the IBZ is not equal to the weight', 1) ENDIF ENDIF ! withtout symmetries @@ -1080,70 +1030,59 @@ DO j = 1, 3 DO i = 1, 3 ij = ij + 1 - tdf_sigma(ij) = vkk(i,ibnd) * vkk(j,ibnd) * wkf(ikk) + tdf_sigma(ij) = vkk(i, ibnd) * vkk(j, ibnd) * wkf(ikk) ENDDO ENDDO - ! ENDIF ! mp_mesh_k ! ! energy at k (relative to Ef) - ekk = etf (ibndmin-1+ibnd, ikk) - ef0(itemp) - tau = one / inv_tau_all(itemp,ibnd,ik+lower_bnd-1) + ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0(itemp) + tau = one / inv_tau_all(itemp, ibnd, ik + lower_bnd - 1) ! derivative Fermi distribution ! (-df_nk/dE_nk) = (f_nk)*(1-f_nk)/ (k_B T) - dfnk = w0gauss( ekk / etemp, -99 ) / etemp + dfnk = w0gauss(ekk / etemp, -99) / etemp ! electrical conductivity - Sigma(:,itemp) = Sigma(:,itemp) + dfnk * tdf_sigma(:) * tau + Sigma(:, itemp) = Sigma(:, itemp) + dfnk * tdf_sigma(:) * tau ! ! Now do the same but with Znk multiplied ! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} ) - Znk = one / ( one + zi_allvb (itemp,ibnd,ik+lower_bnd-1) ) - tau = one / ( Znk * inv_tau_all(itemp,ibnd,ik+lower_bnd-1) ) - SigmaZ(:,itemp) = SigmaZ(:,itemp) + dfnk * tdf_sigma(:) * tau - - !print*,'itemp ik ibnd ',itemp, ik, ibnd - !print*,'Sigma ',Sigma(:,itemp) - !print*,'SigmaZ ',SigmaZ(:,itemp) - !print*,'Znk ',Znk + Znk = one / (one + zi_allvb(itemp, ibnd, ik + lower_bnd - 1)) + tau = one / (Znk * inv_tau_all(itemp, ibnd, ik + lower_bnd - 1)) + SigmaZ(:, itemp) = SigmaZ(:, itemp) + dfnk * tdf_sigma(:) * tau ! ENDIF - ! ENDDO ! ibnd - ! ENDIF ! endif fsthick - ! ENDDO ! end loop on k ! ! The k points are distributed among pools: here we collect them ! - CALL mp_sum( Sigma(:,itemp), world_comm ) - CALL mp_sum( SigmaZ(:,itemp), world_comm ) - !DBSP - !write(stdout,*) 'ef0(itemp) ',ef0(itemp) + CALL mp_sum(Sigma(:, itemp), world_comm) + CALL mp_sum(SigmaZ(:, itemp), world_comm) ! ENDDO ! nstemp ! IF (mpime == meta_ionode_id) THEN filsigma = TRIM(prefix) // '_elcond_h' - OPEN(iufilsigma, file = filsigma, form = 'formatted') - WRITE(iufilsigma,'(a)') "# Electrical conductivity in 1/(Ohm * m)" - WRITE(iufilsigma,'(a)') "# Ef(eV) Temp(K) Sigma_xx Sigma_xy Sigma_xz" // & - " Sigma_yx Sigma_yy Sigma_yz " // & - " Sigma_xz Sigma_yz Sigma_zz" + OPEN(iufilsigma, FILE = filsigma, FORM = 'formatted') + WRITE(iufilsigma, '(a)') "# Electrical conductivity in 1/(Ohm * m)" + WRITE(iufilsigma, '(a)') "# Ef(eV) Temp(K) Sigma_xx Sigma_xy Sigma_xz" // & + " Sigma_yx Sigma_yy Sigma_yz " // & + " Sigma_xz Sigma_yz Sigma_zz" ENDIF ! - conv_factor1 = electron_SI / ( hbar * bohr2ang * Ang2m ) + conv_factor1 = electron_SI / (hbar * bohr2ang * Ang2m) ! - WRITE(stdout,'(/5x,a)') REPEAT('=',67) - WRITE(stdout,'(5x,"Temp [K] Fermi [eV] Hole density [cm^-3] Hole mobility [cm^2/Vs]")') - WRITE(stdout,'(5x,a/)') REPEAT('=',67) + WRITE(stdout, '(/5x,a)') REPEAT('=',67) + WRITE(stdout, '(5x,"Temp [K] Fermi [eV] Hole density [cm^-3] Hole mobility [cm^2/Vs]")') + WRITE(stdout, '(5x,a/)') REPEAT('=',67) ! DO itemp = 1, nstemp etemp = transp_temp(itemp) ! Sigma in units of 1/(a.u.) is converted to 1/(Ohm * m) IF (mpime == meta_ionode_id) THEN - WRITE(iufilsigma,'(11E16.8)') ef0(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, & - conv_factor1 * Sigma(:,itemp) * inv_cell + WRITE(iufilsigma, '(11E16.8)') ef0(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, & + conv_factor1 * Sigma(:, itemp) * inv_cell ENDIF carrier_density = 0.0 ! @@ -1151,49 +1090,49 @@ ikk = 2 * ik - 1 DO ibnd = 1, nbndfst ! This selects only valence bands for hole conduction - IF (etf (ibndmin-1+ibnd, ikk) < ef0(itemp)) THEN + IF (etf(ibndmin - 1 + ibnd, ikk) < ef0(itemp)) THEN ! energy at k (relative to Ef) - ekk = etf (ibndmin-1+ibnd, ikk) - ef0(itemp) - fnk = wgauss( -ekk / etemp, -99) + ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0(itemp) + fnk = wgauss(-ekk / etemp, -99) ! The wkf(ikk) already include a factor 2 - carrier_density = carrier_density + wkf(ikk) * (1.0d0 - fnk ) + carrier_density = carrier_density + wkf(ikk) * (1.0d0 - fnk) ENDIF ENDDO ENDDO ! - CALL mp_sum( carrier_density, world_comm ) + CALL mp_sum(carrier_density, world_comm) ! ! Diagonalize the conductivity matrix ! 1 = (1,1) = xx, 2 = (1,2) = xy, 3 = (1,3) = xz ! 4 = (2,1) = yx, 5 = (2,2) = yy, 6 = (2,3) = yz ! 7 = (3,1) = zx, 8 = (3,2) = zy, 9 = (3,3) = zz sigma_up(:, :) = zero - sigma_up(1,1) = Sigma(1,itemp) - sigma_up(1,2) = Sigma(2,itemp) - sigma_up(1,3) = Sigma(3,itemp) - sigma_up(2,1) = Sigma(4,itemp) - sigma_up(2,2) = Sigma(5,itemp) - sigma_up(2,3) = Sigma(6,itemp) - sigma_up(3,1) = Sigma(7,itemp) - sigma_up(3,2) = Sigma(8,itemp) - sigma_up(3,3) = Sigma(9,itemp) + sigma_up(1, 1) = Sigma(1, itemp) + sigma_up(1, 2) = Sigma(2, itemp) + sigma_up(1, 3) = Sigma(3, itemp) + sigma_up(2, 1) = Sigma(4, itemp) + sigma_up(2, 2) = Sigma(5, itemp) + sigma_up(2, 3) = Sigma(6, itemp) + sigma_up(3, 1) = Sigma(7, itemp) + sigma_up(3, 2) = Sigma(8, itemp) + sigma_up(3, 3) = Sigma(9, itemp) ! - CALL rdiagh(3,sigma_up,3,sigma_eig,sigma_vect) + CALL rdiagh(3, sigma_up, 3, sigma_eig, sigma_vect) ! !Sigma_diag = (Sigma(1,itemp)+Sigma(5,itemp)+Sigma(9,itemp))/3 !Sigma_offdiag = (Sigma(2,itemp)+Sigma(3,itemp)+Sigma(4,itemp)+& ! Sigma(6,itemp)+Sigma(7,itemp)+Sigma(8,itemp))/6 - mobility_xx = ( sigma_eig(1) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility_yy = ( sigma_eig(2) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility_zz = ( sigma_eig(3) * electron_SI * ( bohr2ang * ang2cm )**2) /( carrier_density * hbarJ) - mobility = (mobility_xx+mobility_yy+mobility_zz)/3 + mobility_xx = (sigma_eig(1) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_yy = (sigma_eig(2) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_zz = (sigma_eig(3) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility = (mobility_xx + mobility_yy + mobility_zz) / 3 ! carrier_density in cm^-1 - carrier_density_prt = carrier_density * inv_cell * ( bohr2ang * ang2cm )**(-3) - WRITE(stdout,'(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev /kelvin2eV, & - ef0(itemp)*ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility_yy, ' y-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility_zz, ' z-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility, ' avg' + carrier_density_prt = carrier_density * inv_cell * (bohr2ang * ang2cm)**(-3) + WRITE(stdout, '(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev / kelvin2eV, & + ef0(itemp) * ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility_yy, ' y-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility_zz, ' z-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility, ' avg' ! ! ! Now do Znk ---------------------------------------------------------- ! sigma_up(:, :) = zero @@ -1230,72 +1169,67 @@ ! IF (int_mob .OR. (ncarrier > 1E5)) THEN DO itemp = 1, nstemp - ! etemp = transp_temp(itemp) IF (itemp == 1) THEN - tdf_sigma(:) = zero - Sigma(:, :) = zero - SigmaZ(:, :) = zero + tdf_sigma(:) = zero + Sigma(:, :) = zero + SigmaZ(:, :) = zero ENDIF DO ik = 1, nkf ikk = 2 * ik - 1 - IF (minval ( ABS(etf (:, ikk) - ef) ) < fsthick) THEN - IF (ABS(efcb(itemp)) < eps) THEN + IF (MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) THEN + IF (ABS(efcb(itemp)) < eps4) THEN DO ibnd = 1, nbndfst ! This selects only cond bands for electron conduction - IF (etf (ibndmin-1+ibnd, ikk) > ef0(itemp)) THEN - ! vkk(3,nbnd) - velocity for k + IF (etf(ibndmin - 1 + ibnd, ikk) > ef0(itemp)) THEN tdf_sigma(:) = zero IF (vme) THEN ! vmef is in units of Ryd * bohr - vkk(:,ibnd) = REAL (vmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = REAL(vmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ELSE ! v_(k,i) = 1/m = 2 * dmef (:, i,i,k) ! 1/m = 2 in Rydberg atomic units ! dmef is in units of 1/a.u. (where a.u. is bohr) ! v_(k,i) is in units of Ryd * a.u. - vkk(:,ibnd) = 2.0 * REAL (dmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = 2.0 * REAL(dmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ENDIF IF (mp_mesh_k) THEN ! - vk_cart(:) = vkk(:,ibnd) + vk_cart(:) = vkk(:, ibnd) ! ! Loop on full BZ nb = 0 - DO ikbz = 1, nkf1*nkf2*nkf3 + DO ikbz = 1, nkf1 * nkf2 * nkf3 ! If the k-point from the full BZ is related by a symmetry operation ! to the current k-point, then take it. - IF (BZtoIBZ(ikbz) == ik+lower_bnd-1) THEN + IF (BZtoIBZ(ikbz) == ik + lower_bnd - 1) THEN nb = nb + 1 ! Transform the symmetry matrix from Crystal to cartesian - sa (:, :) = dble ( s_BZtoIBZ(:,:,ikbz) ) - sb = matmul ( bg, sa ) - sr (:, :) = matmul ( at, transpose (sb) ) - CALL dgemv( 'n', 3, 3, 1.d0,& - sr, 3, vk_cart(:),1 ,0.d0 , v_rot(:), 1 ) + sa(:, :) = DBLE(s_BZtoIBZ(:, :, ikbz)) + sb = MATMUL(bg, sa) + sr(:, :) = MATMUL(at, TRANSPOSE(sb)) + CALL dgemv( 'n', 3, 3, 1.d0, sr, 3, vk_cart(:), 1, 0.d0, v_rot(:), 1) ij = 0 DO j = 1, 3 DO i = 1, 3 ij = ij + 1 ! The factor two in the weight at the end is to account for spin IF (noncolin) THEN - tdf_sigma(ij) = tdf_sigma(ij) + ( v_rot(i) * v_rot(j) ) * 1.0 / (nkf1*nkf2*nkf3) + tdf_sigma(ij) = tdf_sigma(ij) + (v_rot(i) * v_rot(j)) * 1.0 / (nkf1 * nkf2 * nkf3) ELSE - tdf_sigma(ij) = tdf_sigma(ij) + ( v_rot(i) * v_rot(j) ) * 2.0 / (nkf1*nkf2*nkf3) + tdf_sigma(ij) = tdf_sigma(ij) + (v_rot(i) * v_rot(j)) * 2.0 / (nkf1 * nkf2 * nkf3) ENDIF ENDDO ENDDO ENDIF ENDDO ! ikbz IF (noncolin) THEN - IF (ABS(nb*1.0/(nkf1*nkf2*nkf3) - wkf(ikk)) > eps6) THEN - CALL errore ('transport', & - &' The number of kpoint in the IBZ is not equal to the weight', 1) + IF (ABS(nb * 1.0 / (nkf1 * nkf2 * nkf3) - wkf(ikk)) > eps6) THEN + CALL errore('transport', ' The number of kpoint in the IBZ is not equal to the weight', 1) ENDIF ELSE - IF (ABS(nb*2.0/(nkf1*nkf2*nkf3) - wkf(ikk)) > eps6) THEN - CALL errore ('transport', & - &' The number of kpoint in the IBZ is not equal to the weight', 1) + IF (ABS(nb * 2.0 / (nkf1 * nkf2 * nkf3) - wkf(ikk)) > eps6) THEN + CALL errore('transport', ' The number of kpoint in the IBZ is not equal to the weight', 1) ENDIF ENDIF ! withtout symmetries @@ -1305,77 +1239,73 @@ DO j = 1, 3 DO i = 1, 3 ij = ij + 1 - tdf_sigma(ij) = vkk(i,ibnd) * vkk(j,ibnd) * wkf(ikk) + tdf_sigma(ij) = vkk(i, ibnd) * vkk(j, ibnd) * wkf(ikk) ENDDO ENDDO - ! ENDIF ! mp_mesh_k ! - ekk = etf (ibndmin-1+ibnd, ikk) - ef0(itemp) - tau = one / inv_tau_all(itemp,ibnd,ik+lower_bnd-1) - dfnk = w0gauss( ekk / etemp, -99 ) / etemp - Sigma(:,itemp) = Sigma(:,itemp) + dfnk * tdf_sigma(:) * tau + ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0(itemp) + tau = one / inv_tau_all(itemp, ibnd, ik + lower_bnd - 1) + dfnk = w0gauss(ekk / etemp, -99) / etemp + Sigma(:, itemp) = Sigma(:, itemp) + dfnk * tdf_sigma(:) * tau ! ! Now do the same but with Znk multiplied ! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} ) - Znk = one / ( one + zi_allvb (itemp,ibnd,ik+lower_bnd-1) ) - tau = one / ( Znk * inv_tau_all(itemp,ibnd,ik+lower_bnd-1) ) - SigmaZ(:,itemp) = SigmaZ(:,itemp) + dfnk * tdf_sigma(:) * tau + Znk = one / (one + zi_allvb(itemp, ibnd, ik + lower_bnd - 1)) + tau = one / (Znk * inv_tau_all(itemp, ibnd, ik + lower_bnd - 1)) + SigmaZ(:, itemp) = SigmaZ(:, itemp) + dfnk * tdf_sigma(:) * tau ENDIF ENDDO ELSE ! In this case we have 2 Fermi levels DO ibnd = 1, nbndfst ! This selects only cond bands for hole conduction - IF (etf (ibndmin-1+ibnd, ikk) > efcb(itemp)) THEN + IF (etf(ibndmin - 1 + ibnd, ikk) > efcb(itemp)) THEN ! ! SP - Uncomment to use symmetries on velocities tdf_sigma(:) = zero IF (vme) THEN - vkk(:,ibnd) = REAL (vmef (:, ibndmin-1+ibnd,ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = REAL(vmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ELSE - vkk(:,ibnd) = 2.0 * REAL (dmef (:, ibndmin-1+ibnd,ibndmin-1+ibnd, ikk)) + vkk(:, ibnd) = 2.0 * REAL(dmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk)) ENDIF ! IF (mp_mesh_k) THEN ! - vk_cart(:) = vkk(:,ibnd) + vk_cart(:) = vkk(:, ibnd) ! ! Loop on full BZ nb = 0 - DO ikbz = 1, nkf1*nkf2*nkf3 + DO ikbz = 1, nkf1 * nkf2 * nkf3 ! If the k-point from the full BZ is related by a symmetry operation ! to the current k-point, then take it. - IF (BZtoIBZ(ikbz) == ik+lower_bnd-1) THEN + IF (BZtoIBZ(ikbz) == ik + lower_bnd - 1) THEN nb = nb + 1 ! Transform the symmetry matrix from Crystal to cartesian - sa (:, :) = dble ( s_BZtoIBZ(:,:,ikbz) ) - sb = matmul ( bg, sa ) - sr (:, :) = matmul ( at, transpose (sb) ) - CALL dgemv( 'n', 3, 3, 1.d0,& - sr, 3, vk_cart(:),1 ,0.d0 , v_rot(:), 1 ) + sa(:, :) = DBLE(s_BZtoIBZ(:, :, ikbz)) + sb = MATMUL(bg, sa) + sr(:, :) = MATMUL(at, TRANSPOSE(sb)) + CALL dgemv('n', 3, 3, 1.d0, sr, 3, vk_cart(:), 1, 0.d0, v_rot(:), 1) ij = 0 DO j = 1, 3 DO i = 1, 3 ij = ij + 1 ! The factor two in the weight at the end is to account for spin IF (noncolin) THEN - tdf_sigma(ij) = tdf_sigma(ij) + ( v_rot(i) * v_rot(j) ) * 1.0 / (nkf1*nkf2*nkf3) + tdf_sigma(ij) = tdf_sigma(ij) + (v_rot(i) * v_rot(j)) * 1.0 / (nkf1 * nkf2 * nkf3) ELSE - tdf_sigma(ij) = tdf_sigma(ij) + ( v_rot(i) * v_rot(j) ) * 2.0 / (nkf1*nkf2*nkf3) + tdf_sigma(ij) = tdf_sigma(ij) + (v_rot(i) * v_rot(j)) * 2.0 / (nkf1 * nkf2 * nkf3) ENDIF ENDDO ENDDO ENDIF ENDDO ! ikbz IF (noncolin) THEN - IF (ABS(nb*1.0/(nkf1*nkf2*nkf3) - wkf(ikk)) > eps6) THEN - CALL errore ('transport', & - &' The number of kpoint in the IBZ is not equal to the weight', 1) + IF (ABS(nb * 1.0 / (nkf1 * nkf2 * nkf3) - wkf(ikk)) > eps6) THEN + CALL errore('transport', ' The number of kpoint in the IBZ is not equal to the weight', 1) ENDIF ELSE - IF (ABS(nb*2.0/(nkf1*nkf2*nkf3) - wkf(ikk)) > eps6) THEN - CALL errore ('transport', & - &' The number of kpoint in the IBZ is not equal to the weight', 1) + IF (ABS(nb * 2.0 / (nkf1 * nkf2 * nkf3) - wkf(ikk)) > eps6) THEN + CALL errore ('transport', ' The number of kpoint in the IBZ is not equal to the weight', 1) ENDIF ENDIF ! withtout symmetries @@ -1385,77 +1315,59 @@ DO j = 1, 3 DO i = 1, 3 ij = ij + 1 - tdf_sigma(ij) = vkk(i,ibnd) * vkk(j,ibnd) * wkf(ikk) + tdf_sigma(ij) = vkk(i, ibnd) * vkk(j, ibnd) * wkf(ikk) ENDDO ENDDO ! ENDIF ! mp_mesh_k - ekk = etf (ibndmin-1+ibnd, ikk) - efcb(itemp) - tau = one / inv_tau_allcb(itemp,ibnd,ik+lower_bnd-1) - dfnk = w0gauss( ekk / etemp, -99 ) / etemp - Sigma(:,itemp) = Sigma(:,itemp) + dfnk * tdf_sigma(:) * tau - - - !IF (vme) THEN - ! vkk(:,ibnd) = REAL (vmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) - !ELSE - ! vkk(:,ibnd) = 2.0 * REAL (dmef (:, ibndmin-1+ibnd, ibndmin-1+ibnd, ikk)) - !ENDIF - !ekk = etf (ibndmin-1+ibnd, ikk) - efcb(itemp) - !tau = one / inv_tau_allcb(itemp,ibnd,ik+lower_bnd-1) - !ij = 0 - !DO j = 1, 3 - ! DO i = 1, 3 - ! ij = ij + 1 - ! tdf_sigma(ij) = vkk(i,ibnd) * vkk(j,ibnd) * tau - ! ENDDO - !ENDDO - !dfnk = w0gauss( ekk / etemp, -99 ) / etemp - !Sigma(:,itemp) = Sigma(:,itemp) + wkf(ikk) * dfnk * tdf_sigma(:) - ! + ekk = etf(ibndmin - 1 + ibnd, ikk) - efcb(itemp) + tau = one / inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1) + dfnk = w0gauss(ekk / etemp, -99) / etemp + Sigma(:, itemp) = Sigma(:, itemp) + dfnk * tdf_sigma(:) * tau + ! ! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} ) - Znk = one / ( one + zi_allcb (itemp,ibnd,ik+lower_bnd-1) ) - tau = one / ( Znk * inv_tau_allcb(itemp,ibnd,ik+lower_bnd-1) ) + Znk = one / (one + zi_allcb(itemp, ibnd, ik + lower_bnd - 1)) + tau = one / (Znk * inv_tau_allcb(itemp, ibnd, ik + lower_bnd - 1)) ij = 0 DO j = 1, 3 DO i = 1, 3 ij = ij + 1 - tdf_sigma(ij) = vkk(i,ibnd) * vkk(j,ibnd) * tau + tdf_sigma(ij) = vkk(i, ibnd) * vkk(j, ibnd) * tau ENDDO ENDDO - SigmaZ(:,itemp) = SigmaZ(:,itemp) + wkf(ikk) * dfnk * tdf_sigma(:) + SigmaZ(:, itemp) = SigmaZ(:, itemp) + wkf(ikk) * dfnk * tdf_sigma(:) ENDIF ENDDO ! ibnd ENDIF ! etcb ENDIF ! endif fsthick ENDDO ! end loop on k - CALL mp_sum( Sigma(:,itemp), world_comm ) - CALL mp_sum( SigmaZ(:,itemp), world_comm ) + CALL mp_sum(Sigma(:, itemp), world_comm) + CALL mp_sum(SigmaZ(:, itemp), world_comm) ! ENDDO ! nstemp IF (mpime == meta_ionode_id) THEN filsigma = TRIM(prefix) // '_elcond_e' - OPEN(iufilsigma, file = filsigma, form = 'formatted') - WRITE(iufilsigma,'(a)') "# Electrical conductivity in 1/(Ohm * m)" - WRITE(iufilsigma,'(a)') "# Ef(eV) Temp(K) Sigma_xx Sigma_xy Sigma_xz" // & - " Sigma_yx Sigma_yy Sigma_yz " // & - " Sigma_xz Sigma_yz Sigma_zz" + OPEN(iufilsigma, FILE = filsigma, FORM = 'formatted') + WRITE(iufilsigma, '(a)') "# Electrical conductivity in 1/(Ohm * m)" + WRITE(iufilsigma, '(a)') "# Ef(eV) Temp(K) Sigma_xx Sigma_xy Sigma_xz" // & + " Sigma_yx Sigma_yy Sigma_yz " // & + " Sigma_xz Sigma_yz Sigma_zz" ENDIF ! conv_factor1 = electron_SI / ( hbar * bohr2ang * Ang2m ) - WRITE(stdout,'(/5x,a)') REPEAT('=',67) - WRITE(stdout,'(5x,"Temp [K] Fermi [eV] Elec density [cm^-3] Elec mobility [cm^2/Vs]")') - WRITE(stdout,'(5x,a/)') REPEAT('=',67) + WRITE(stdout, '(/5x,a)') REPEAT('=',67) + WRITE(stdout, '(5x,"Temp [K] Fermi [eV] Elec density [cm^-3] Elec mobility [cm^2/Vs]")') + WRITE(stdout, '(5x,a/)') REPEAT('=',67) DO itemp = 1, nstemp etemp = transp_temp(itemp) IF (mpime == meta_ionode_id) THEN ! Sigma in units of 1/(a.u.) is converted to 1/(Ohm * m) - IF (ABS(efcb(itemp)) < eps) THEN - WRITE(iufilsigma,'(11E16.8)') ef0(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, & - conv_factor1 * Sigma(:,itemp) * inv_cell + IF (ABS(efcb(itemp)) < eps4) THEN + WRITE(iufilsigma, '(11E16.8)') ef0(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, & + conv_factor1 * Sigma(:, itemp) * inv_cell ELSE - WRITE(iufilsigma,'(11E16.8)') efcb(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, & - conv_factor1 * Sigma(:,itemp) * inv_cell + WRITE(iufilsigma, '(11E16.8)') efcb(itemp) * ryd2ev, etemp * ryd2ev / kelvin2eV, & + conv_factor1 * Sigma(:, itemp) * inv_cell ENDIF ENDIF carrier_density = 0.0 @@ -1464,85 +1376,81 @@ DO ibnd = 1, nbndfst ikk = 2 * ik - 1 ! This selects only conduction bands for electron conduction - IF (ABS(efcb(itemp)) < eps) THEN - IF (etf (ibndmin-1+ibnd, ikk) > ef0(itemp)) THEN - ekk = etf (ibndmin-1+ibnd, ikk) - ef0(itemp) + IF (ABS(efcb(itemp)) < eps4) THEN + IF (etf(ibndmin - 1 + ibnd, ikk) > ef0(itemp)) THEN + ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0(itemp) fnk = wgauss( -ekk / etemp, -99) ! The wkf(ikk) already include a factor 2 carrier_density = carrier_density + wkf(ikk) * fnk ENDIF ELSE - IF (etf (ibndmin-1+ibnd, ikk) > efcb(itemp)) THEN - ekk = etf (ibndmin-1+ibnd, ikk) - efcb(itemp) - fnk = wgauss( -ekk / etemp, -99) + IF (etf(ibndmin - 1 + ibnd, ikk) > efcb(itemp)) THEN + ekk = etf(ibndmin - 1 + ibnd, ikk) - efcb(itemp) + fnk = wgauss(-ekk / etemp, -99) ! The wkf(ikk) already include a factor 2 carrier_density = carrier_density + wkf(ikk) * fnk ENDIF ENDIF ENDDO ENDDO - CALL mp_sum( carrier_density, world_comm ) + CALL mp_sum(carrier_density, world_comm) ! Diagonalize the conductivity matrix ! 1 = (1,1) = xx, 2 = (1,2) = xy, 3 = (1,3) = xz ! 4 = (2,1) = yx, 5 = (2,2) = yy, 6 = (2,3) = yz ! 7 = (3,1) = zx, 8 = (3,2) = zy, 9 = (3,3) = zz sigma_up(:, :) = zero - sigma_up(1,1) = Sigma(1,itemp) - sigma_up(1,2) = Sigma(2,itemp) - sigma_up(1,3) = Sigma(3,itemp) - sigma_up(2,1) = Sigma(4,itemp) - sigma_up(2,2) = Sigma(5,itemp) - sigma_up(2,3) = Sigma(6,itemp) - sigma_up(3,1) = Sigma(7,itemp) - sigma_up(3,2) = Sigma(8,itemp) - sigma_up(3,3) = Sigma(9,itemp) + sigma_up(1, 1) = Sigma(1, itemp) + sigma_up(1, 2) = Sigma(2, itemp) + sigma_up(1, 3) = Sigma(3, itemp) + sigma_up(2, 1) = Sigma(4, itemp) + sigma_up(2, 2) = Sigma(5, itemp) + sigma_up(2, 3) = Sigma(6, itemp) + sigma_up(3, 1) = Sigma(7, itemp) + sigma_up(3, 2) = Sigma(8, itemp) + sigma_up(3, 3) = Sigma(9, itemp) ! - CALL rdiagh(3,sigma_up,3,sigma_eig,sigma_vect) + CALL rdiagh(3, sigma_up, 3, sigma_eig, sigma_vect) ! - !Sigma_diag = (Sigma(1,itemp)+Sigma(5,itemp)+Sigma(9,itemp))/3 - !Sigma_offdiag = (Sigma(2,itemp)+Sigma(3,itemp)+Sigma(4,itemp)+& - ! Sigma(6,itemp)+Sigma(7,itemp)+Sigma(8,itemp))/6 - mobility_xx = ( sigma_eig(1) * electron_SI * ( bohr2ang * ang2cm )**2) / ( carrier_density * hbarJ) - mobility_yy = ( sigma_eig(2) * electron_SI * ( bohr2ang * ang2cm )**2) / ( carrier_density * hbarJ) - mobility_zz = ( sigma_eig(3) * electron_SI * ( bohr2ang * ang2cm )**2) / ( carrier_density * hbarJ) - mobility = (mobility_xx+mobility_yy+mobility_zz)/3 + mobility_xx = (sigma_eig(1) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_yy = (sigma_eig(2) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_zz = (sigma_eig(3) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility = (mobility_xx + mobility_yy + mobility_zz) / 3 ! - - ! carrier_density in cm^-1 - carrier_density_prt = carrier_density * inv_cell * ( bohr2ang * ang2cm )**(-3) - IF (ABS(efcb(itemp)) < eps) THEN - WRITE(stdout,'(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev / kelvin2eV,& - ef0(itemp)*ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' + ! Carrier_density in cm^-1 + carrier_density_prt = carrier_density * inv_cell * (bohr2ang * ang2cm)**(-3) + IF (ABS(efcb(itemp)) < eps4) THEN + WRITE(stdout, '(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev / kelvin2eV,& + ef0(itemp) * ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' ELSE - WRITE(stdout,'(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev / kelvin2eV,& - efcb(itemp)*ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' + WRITE(stdout, '(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev / kelvin2eV,& + efcb(itemp) * ryd2ev, carrier_density_prt, mobility_xx, ' x-axis' ENDIF - WRITE(stdout,'(45x, 1E18.6, a)') mobility_yy, ' y-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility_zz, ' z-axis' - WRITE(stdout,'(45x, 1E18.6, a)') mobility, ' avg' + WRITE(stdout, '(45x, 1E18.6, a)') mobility_yy, ' y-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility_zz, ' z-axis' + WRITE(stdout, '(45x, 1E18.6, a)') mobility, ' avg' ! Issue warning if the material is anisotropic ! IF (Sigma_offdiag > 0.1*Sigma_diag) THEN ! WRITE(stdout,'(5x,a,1f10.5,a)') 'Warning: Sigma_offdiag = ',(Sigma_offdiag*100)/Sigma_diag, '% of Sigma_diag' ! ENDIF ! Now do the mobility with Znk factor ---------------------------------------------------------- sigma_up(:, :) = zero - sigma_up(1,1) = SigmaZ(1,itemp) - sigma_up(1,2) = SigmaZ(2,itemp) - sigma_up(1,3) = SigmaZ(3,itemp) - sigma_up(2,1) = SigmaZ(4,itemp) - sigma_up(2,2) = SigmaZ(5,itemp) - sigma_up(2,3) = SigmaZ(6,itemp) - sigma_up(3,1) = SigmaZ(7,itemp) - sigma_up(3,2) = SigmaZ(8,itemp) - sigma_up(3,3) = SigmaZ(9,itemp) - CALL rdiagh(3,sigma_up,3,sigma_eig,sigma_vect) - mobility_xx = ( sigma_eig(1) * electron_SI * ( bohr2ang * ang2cm )**2) / ( carrier_density * hbarJ) - mobility_yy = ( sigma_eig(2) * electron_SI * ( bohr2ang * ang2cm )**2) / ( carrier_density * hbarJ) - mobility_zz = ( sigma_eig(3) * electron_SI * ( bohr2ang * ang2cm )**2) / ( carrier_density * hbarJ) - mobility = (mobility_xx+mobility_yy+mobility_zz)/3 + sigma_up(1, 1) = SigmaZ(1, itemp) + sigma_up(1, 2) = SigmaZ(2, itemp) + sigma_up(1, 3) = SigmaZ(3, itemp) + sigma_up(2, 1) = SigmaZ(4, itemp) + sigma_up(2, 2) = SigmaZ(5, itemp) + sigma_up(2, 3) = SigmaZ(6, itemp) + sigma_up(3, 1) = SigmaZ(7, itemp) + sigma_up(3, 2) = SigmaZ(8, itemp) + sigma_up(3, 3) = SigmaZ(9, itemp) + CALL rdiagh(3, sigma_up, 3, sigma_eig, sigma_vect) + mobility_xx = (sigma_eig(1) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_yy = (sigma_eig(2) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility_zz = (sigma_eig(3) * electron_SI * (bohr2ang * ang2cm)**2) / (carrier_density * hbarJ) + mobility = (mobility_xx + mobility_yy + mobility_zz) / 3 ! ! DBSP - Z-factor - ! IF (ABS(efcb(itemp)) < eps) THEN + ! IF (ABS(efcb(itemp)) < eps4) THEN ! WRITE(stdout,'(5x, 1f8.3, 1f12.4, 1E19.6, 1E19.6, a)') etemp * ryd2ev / kelvin2eV,& ! ef0(itemp)*ryd2ev, carrier_density_prt, mobility_xx, ' x-axis [Z]' ! ELSE @@ -1567,13 +1475,14 @@ ! CALL stop_clock('MOB') ! - WRITE( stdout, * ) ' Total time so far' - CALL print_clock ('SCAT') - CALL print_clock ('MOB') - WRITE(stdout,'(5x)') + WRITE(stdout, * ) ' Total time so far' + CALL print_clock('SCAT') + CALL print_clock('MOB') + WRITE(stdout, '(5x)') ! RETURN ! + !-------------------------------------------------------------------------- END SUBROUTINE transport_coeffs !-------------------------------------------------------------------------- !