! !********************************************************************* !********************************************************************* !**** **** !**** EPHWANN **** !**** **** !********************************************************************* !********************************************************************* ! ! Wannier interpolation of electron-phonon vertex ! ! Scalar implementation Feb 2006 ! Parallel version May 2006 ! Disentenglement Oct 2006 ! Compact formalism Dec 2006 ! Phonon irreducible zone Mar 2007 ! ! No ultrasoft now ! No spin polarization ! ! Feliciano Giustino, University of California at Berkeley ! !----------------------------------------------------------------------- subroutine ephwann_shuffle ( nqc, xqc ) !----------------------------------------------------------------------- ! #include "f_defs.h" use kinds, only : DP use pwcom, only : nbnd, nks, nkstot, & nk1, nk2, nk3, et, xk, at, bg, ef, amass use phcom, only : dyn, nbndsub, nksq, lgamma, nq1, nq2, nq3, & nmodes, lrepmatf, iunepmatf, fsthick, w2, iuetf, lretf, epwread, epwwrite use control_flags, only : iverbosity use io_files, only: prefix use ions_base, only : ityp use io_global, only : stdout use el_phon #ifdef __PARA use para use mp, only : mp_barrier USE mp_global, only : mpime use io_global, only : ionode_id #endif implicit none ! complex(kind=DP), allocatable :: & chw (:,:,:), &! Hamiltonian in wannier basis epmatwe (:,:,:,:,:), &! e-p matrix in wannier basis - electrons epmatwp (:,:,:,:,:), &! e-p matrix in wannier basis - electrons and phonons epmatwef (:,:,:,:), &! e-p matrix in el wannier - fine Bloch phonon grid chf(:,:) ! Hamiltonian in smooth Bloch basis, fine mesh !@ real(kind=DP), allocatable :: & complex(kind=DP), allocatable :: & rdw (:,:,:) ! dynamical matrix in wannier basis (real) complex(kind=DP), allocatable :: & epmatf( :, :, :), &! e-p matrix in smooth Bloch basis, fine mesh cufkk ( :, :), &! Rotation matrix, fine mesh, points k cufkq ( :, :), &! the same, for points k+q uf ( :, :) ! Rotation matrix for phonons integer :: & nqc, &! number of qpoints in the coarse grid iunepwdata ! number of qpoints in the coarse grid real(kind=DP) :: & xqc (3, nqc) ! qpoint list, coarse mesh ! ! work variables ! integer :: iq, ik, ikk, ikq, ibnd, jbnd, imode, jmode,ir, na, nu, mu, & fermicount, nrec, indnew, indold, irk, irq, ipool logical :: exst character (len=256) :: filint real(kind=DP) :: xxq(3), xxk(3), xkk(3), xkq(3), ww(nmodes), tmp1, tmp2 real(kind=DP), parameter :: rydcm1 = 13.6058d0 * 8065.5d0 ! unit for I/O all quantities in Wannier representation iunepwdata = 79 ! if (nbndsub.ne.nbnd) & write(stdout, '(14x,a)' ) 'band disentanglement is used' ! allocate ( cu ( nbnd, nbndsub, nks), & irvec (3, 2*nk1*nk2*nk3), & ndegen_k (2*nk1*nk2*nk3), & ndegen_q (2*nq1*nq2*nq3), & wslen(2*nk1*nk2*nk3) ) ! call start_clock ( 'ephwann' ) ! call loadkmesh ! ! determine Wigner-Seitz points ! call wigner_seitz2 & ( nk1, nk2, nk3, nq1, nq2, nq3, nrr_k, nrr_q, irvec, wslen, ndegen_k, ndegen_q ) ! #ifdef __PARA ! this is to make sure that everybody allocates at the same time, mainly for debugging purposes call mp_barrier() #endif ! if ( epwread ) then ! ! read all quantities in Wannier representation from file ! in parallel case all pools read the same file ! write(6,'(/5x,"Reading Hamiltonian, Dynamical matrix and EP vertex in Wannier rep from file"/)') ! ! this sequential reading is incredibly slow as compared to letting the operating system care about it !#ifdef __PARA ! do ipool = 1, npool !#endif ! open(unit=iunepwdata,file='epwdata.fmt') read (iunepwdata,*) ef read (iunepwdata,*) nbndsub, nrr_k, nmodes, nrr_q ! allocate ( epmatwp ( nbndsub, nbndsub, nrr_k, nmodes, nrr_q), & chw ( nbndsub, nbndsub, nrr_k ), rdw( nmodes, nmodes, nrr_q ) ) ! do ibnd = 1, nbndsub do jbnd = 1, nbndsub do irk = 1, nrr_k read (iunepwdata,*) chw(ibnd,jbnd,irk) enddo enddo enddo ! do imode = 1, nmodes do jmode = 1, nmodes do irq = 1, nrr_q read (iunepwdata,*) rdw(imode,jmode,irq) enddo enddo enddo ! do ibnd = 1, nbndsub do jbnd = 1, nbndsub do irk = 1, nrr_k do imode = 1, nmodes do irq = 1, nrr_q read (iunepwdata,*) epmatwp(ibnd,jbnd,irk,imode,irq) enddo enddo enddo enddo enddo close(iunepwdata) ! !#ifdef __PARA ! endif ! call mp_barrier ! enddo !#endif else ! xxq = 0.d0 ! not used at this point call loadumat & ( nbnd, nbndsub, nks, nkstot, xxq, cu ) ! ! ------------------------------------------------------ ! Bloch to Wannier transform ! ------------------------------------------------------ ! allocate ( epmatwe ( nbndsub, nbndsub, nrr_k, nmodes, nqc), & epmatwp ( nbndsub, nbndsub, nrr_k, nmodes, nrr_q), & chw ( nbndsub, nbndsub, nrr_k ), & rdw ( nmodes, nmodes, nrr_q ) ) ! ! Hamiltonian ! call hambloch2wan & ( nbnd, nbndsub, nks, nksq, nkbl, lgamma, et, xk, cu, nrr_k, irvec, wslen, chw ) ! ! Dynamical Matrix ! call dynbloch2wan & ( nmodes, nqc, xqc, dynq, nrr_q, irvec, wslen, rdw ) ! ! Electron-Phonon vertex (Bloch el and Bloch ph -> Wannier el and Bloch ph) ! do iq = 1, nqc ! xxq = xqc (:, iq) ! ! we need the cu again for the k+q points, we generate the map here ! call loadumat ( nbnd, nbndsub, nks, nkstot, xxq, cu ) ! do imode = 1, nmodes ! call ephbloch2wane & ( nbnd, nbndsub, nks, nksq, nkbl, lgamma, xk, cu, & epmatq (:,:,:,imode,iq), nrr_k, irvec, wslen, epmatwe(:,:,:,imode,iq) ) ! enddo ! enddo ! ! Electron-Phonon vertex (Wannier el and Bloch ph -> Wannier el and Wannier ph) ! call ephbloch2wanp & ( nbndsub, nmodes, xqc, nqc, dynq, irvec, wslen, nrr_k, nrr_q, epmatwe, epmatwp ) ! endif if ( epwwrite ) then ! ! write all quantities in Wannier representation to file ! write(6,'(/5x,"Writing Hamiltonian, Dynamical matrix and EP vertex in Wannier rep to file"/)') ! #ifdef __PARA if (mpime.eq.ionode_id) then #endif ! open(unit=iunepwdata,file='epwdata.fmt') write (iunepwdata,*) ef write (iunepwdata,*) nbndsub, nrr_k, nmodes, nrr_q do ibnd = 1, nbndsub do jbnd = 1, nbndsub do irk = 1, nrr_k write (iunepwdata,*) chw(ibnd,jbnd,irk) enddo enddo enddo ! do imode = 1, nmodes do jmode = 1, nmodes do irq = 1, nrr_q write (iunepwdata,*) rdw(imode,jmode,irq) enddo enddo enddo ! do ibnd = 1, nbndsub do jbnd = 1, nbndsub do irk = 1, nrr_k do imode = 1, nmodes do irq = 1, nrr_q write (iunepwdata,*) epmatwp(ibnd,jbnd,irk,imode,irq) enddo enddo enddo enddo enddo close(iunepwdata) ! #ifdef __PARA endif #endif ! endif ! if ( allocated (epmatwe) ) deallocate (epmatwe) if ( allocated (epmatq) ) deallocate (epmatq) if ( allocated (cu) ) deallocate (cu) ! allocate ( epmatwef( nbndsub, nbndsub, nrr_k, nmodes), & wf ( nmodes, nxqf ), etf ( nbndsub, nksf), & epmatf( nbndsub, nbndsub, nmodes), cufkk ( nbndsub, nbndsub), & cufkq ( nbndsub, nbndsub), uf ( nmodes, nmodes) ) ! ! test only ! do iq = 1, nxqf ! xxq = xqf (:, iq) ! ! ------------------------------------------------------ ! dynamical matrix : Wannier -> Bloch ! ------------------------------------------------------ ! call dynwan2bloch & ( nmodes, nrr_q, irvec, ndegen_q, rdw, xxq, uf, w2 ) ! do nu = 1, nmodes if ( w2 (nu) .gt. 0.d0 ) then wf(nu,iq) = sqrt(abs( w2 (nu) )) else wf(nu,iq) = -sqrt(abs( w2 (nu) )) endif enddo ! write ( stdout, '(6(2x,f10.5))' ) (wf(nu,iq)*rydcm1, nu = 1, nmodes) ! enddo ! ! ------------------------------------------------------ ! hamiltonian : Wannier -> Bloch (preliminary) ! ------------------------------------------------------ ! ! we here perform a preliminary interpolation of the hamiltonian ! in order to determine the fermi window ibndmin:ibndmax for later use. ! We will interpolate again afterwards, for each k and k+q separately ! xxq = 0.d0 do ik = 1, nksf ! xxk = xkf (:, ik) ! if (.not.lgamma .and. 2*(ik/2).eq.ik ) then ! ! this is a k+q point : redefine as xkf (:, ik-1) + xxq ! call cryst_to_cart ( 1, xxq, at,-1 ) xxk = xkf (:, ik-1) + xxq call cryst_to_cart ( 1, xxq, bg, 1 ) ! endif ! call hamwan2bloch & ( nbndsub, nrr_k, irvec, ndegen_k, chw, xxk, cufkk, etf (:, ik)) ! enddo ! ! identify the bands within fsthick from the Fermi level ! (in shuffle mode this actually does not depend on q) ! call fermiwindow ! ! and open the .epf file[s] with the proper record lenght ! lrepmatf = 2 * (ibndmax-ibndmin+1) * (ibndmax-ibndmin+1) filint = trim(prefix)//'.epf' call diropn (iunepmatf, filint, lrepmatf, exst) ! lretf = (ibndmax-ibndmin+1) filint = trim(prefix)//'.etf' call diropn (iuetf, filint, lretf, exst) ! ! xqf must be in crystal coordinates ! do iq = 1, nxqf ! if (iverbosity .eq. 1) then write(6,'(/5x,"Interpolating iq = ",i6," out of ",i6)') iq, nxqf else ! ! sorry I had a compulsive-obsessive desire to try this... FG ! if (iq.eq.1) then write(6,'(/5x,"Interpolation progress: ")',advance='no') indold = 0 endif indnew = nint(float(iq)/float(nxqf)*43) if (indnew.ne.indold) write(6,'(a)',advance='no') '#' indold = indnew call flush(6) ! only on opteron endif ! xxq = xqf (:, iq) ! ! ------------------------------------------------------ ! dynamical matrix : Wannier -> Bloch ! ------------------------------------------------------ ! call dynwan2bloch & ( nmodes, nrr_q, irvec, ndegen_q, rdw, xxq, uf, w2 ) ! ! ...then take into account the mass factors and square-root the frequencies... ! do nu = 1, nmodes ! if ( w2 (nu) .gt. 0.d0 ) then wf(nu,iq) = sqrt(abs( w2 (nu) )) else wf(nu,iq) = -sqrt(abs( w2 (nu) )) endif ! do mu = 1, nmodes na = (mu - 1) / 3 + 1 uf (mu, nu) = uf (mu, nu) / sqrt(amass(ityp(na))) enddo enddo ! write ( stdout, '(6(2x,f10.5))' ) (wf(nu,iq)*rydcm1, nu = 1, nmodes) ! ! -------------------------------------------------------------- ! epmat : Wannier el and Wannier ph -> Wannier el and Bloch ph ! -------------------------------------------------------------- ! call ephwan2blochp & ( nmodes, xxq, irvec, ndegen_q, nrr_q, uf, epmatwp, epmatwef, nbndsub, nrr_k ) ! ! number of k points with a band on the Fermi surface fermicount = 0 ! do ik = 1, nksqf ! ! xkf is assumed to be in crys coord ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! xkk = xkf(:, ikk) xkq = xkk + xxq ! ! ------------------------------------------------------ ! hamiltonian : Wannier -> Bloch ! ------------------------------------------------------ ! call hamwan2bloch & ( nbndsub, nrr_k, irvec, ndegen_k, chw, xkk, cufkk, etf (:, ikk)) call hamwan2bloch & ( nbndsub, nrr_k, irvec, ndegen_k, chw, xkq, cufkq, etf (:, ikq)) ! ! we write to file the interpolated hamiltonian aigenvalues ! (the eigenvalues for ikk are written nxqf times, which is ! very redundant but easier to handle. If we have a big system ! this can be dropped) ! nrec = (iq-1) * nksf + ikk call davcio ( etf (ibndmin:ibndmax, ikk), ibndmax-ibndmin+1, iuetf, nrec, + 1) nrec = (iq-1) * nksf + ikq call davcio ( etf (ibndmin:ibndmax, ikq), ibndmax-ibndmin+1, iuetf, nrec, + 1) ! ! interpolate only when (k,k+q) both have at least one band ! within a Fermi shell of size fsthick ! if ( ( minval ( abs(etf (:, ikk) - ef) ) .lt. fsthick ) .and. & ( minval ( abs(etf (:, ikq) - ef) ) .lt. fsthick ) ) then ! fermicount = fermicount + 1 ! ! -------------------------------------------------------------- ! epmat : Wannier el and Bloch ph -> Bloch el and Bloch ph ! -------------------------------------------------------------- ! call ephwan2bloch & ( nbndsub, nrr_k, irvec, ndegen_k, epmatwef, xkk, cufkk, cufkq, epmatf, nmodes ) ! ! write epmatf to file ! do imode = 1, nmodes ! nrec = (iq-1) * nmodes * nksqf + (imode-1) * nksqf + ik call dasmio ( epmatf(ibndmin:ibndmax,ibndmin:ibndmax,imode), & ibndmax-ibndmin+1, lrepmatf, iunepmatf, nrec, +1) ! enddo ! endif ! enddo ! enddo ! write( stdout,*) write( stdout, '(/5x,a,i5,a,i5)' ) & 'Number of (k,k+q) pairs on the Fermi surface: ', & fermicount, ' out of ', nksqf ! call stop_clock ( 'ephwann' ) ! return end subroutine ephwann_shuffle !