! !********************************************************************* !********************************************************************* !**** **** !**** EPHWANN **** !**** **** !********************************************************************* !********************************************************************* ! ! Wannier interpolation of electron-phonon vertex ! ! Scalar implementation Feb 2006 ! Parallel version May 2006 ! Disentenglement Oct 2006 ! ! No ultrasoft now ! No spin polarization ! ! Feliciano Giustino, University of California at Berkeley ! !==================================================================== ! ! subroutines and modules for the electron-phonon interpolation ! ! * ../PH/ephwann.f90 ! * ../PH/elphel2.f90 ! * ../PH/elphsum3_strict.f90 ! * ../PH/setphases.f90 ! * ../PH/elph.f90 ! * ../PW/refold.f90 ! * ../PW/set_kplusq.f90 ! * ../pwtools/q2r2_strict.f90 ! * ../pwtools/phinterp.f90 ! * ../pwtools/phinterpmod.f90 ! * ../pwtools/alpha2F.f90 ! * ../pwtools/unfold.f90 ! !----------------------------------------------------------------------- subroutine ephwann (imode0) !----------------------------------------------------------------------- ! #include "f_defs.h" ! USE kinds, only : DP USE io_global, ONLY : stdout USE ions_base, ONLY : nat use pwcom use phcom use el_phon use control_flags, ONLY : iverbosity use io_files, only: prefix USE para, ONLY : kunit #ifdef __PARA use para USE io_global, ONLY : ionode_id USE mp_global, ONLY : nproc, my_pool_id, nproc_pool, & intra_image_comm, inter_pool_comm, me_pool, root_pool, & mpime, intra_pool_comm USE mp, ONLY : mp_barrier, mp_bcast #endif implicit none ! real(kind=DP), parameter :: ryd2ev = 13.60569172, & bohr2ang = 0.5291772108, & twopi = 6.28318530717959,& eps = 1.0e-4 ! complex(kind=DP), parameter :: ci = (0.d0,1.d0), & czero = (0.d0, 0.d0), & cone = (1.d0, 0.d0) ! complex(kind=DP), allocatable :: &! chs(:,:,:), &! Hamiltonian in smooth Bloch basis, coarse mesh epmats (:,:,:), &! e-p matrix in smooth Bloch basis, coarse mesh chw (:,:,:), &! Hamiltonian in wannier basis epmatw (:,:,:), &! e-p matrix in wannier basis cuf(:,:), &! Rotation matrix, fine mesh cufkk(:,:), &! the same, for points k cufkq(:,:), &! the same, for points k+q chf(:,:), &! Hamiltonian in smooth Bloch basis, fine mesh epmatf (:,:), &! e-p matrix in smooth Bloch basis, fine mesh eptmp (:,:) ! e-p matrix, temporary ! integer :: ik, ikk, ikq, ibnd, jbnd, imode0, ipol, nsize, rest, & ik0, ind, ix, jx, ir, mbnd, irkk, irkq, irk, pbnd, i ! integer, parameter :: & iunkkf = 70, &! unit with fine k point mesh (crys coord) iunukk = 71, &! unit with rotation matrix U(k) from wannier code iunukq = 72, &! unit with rotation matrix U(k+q) from wannier code iunkmap = 73 ! unit with map {k+q} --> {k} for shuffle mode ! ! variables for lapack ZHPEVX ! integer :: neig, info integer, allocatable :: ifail(:), iwork(:) real(kind=DP), allocatable :: w(:), rwork(:) complex(kind=DP), allocatable :: champ(:), cwork(:), cz(:,:) ! ! work variables ! real(kind=DP) :: rdotk, rnorm, rvec(3), tmp, xxq(3), dirc(3,3), & ebnd, ebndmin, ebndmax real(kind=DP), allocatable :: xkf_ (:,:), wkf_ (:), aux(:,:) complex(kind=DP) :: cfac logical :: tfermikk, tfermikq, exst integer :: fermicount, nrec character (len=256) :: filint ! call start_clock ( 'ephwann' ) ! if (imode0.eq.0) then ! !==================================================================== write(stdout, '(/5x,a)' ) 'EPHWANN: load fine k mesh and distribute' !==================================================================== ! call start_clock ( 'load data' ) ! #ifdef __PARA if (mpime.eq.ionode_id) then #endif ! ! read fine mesh ! open ( unit = iunkkf, file = filkf, status = 'old', form = 'formatted') read(iunkkf, *) nkstotf write( stdout, '(14x,"Size of k point mesh for interpolation: ",i10)' ) nkstotf ! #ifdef __PARA endif ! call mp_bcast (nkstotf, ionode_id, inter_pool_comm) call mp_bcast (nkstotf, root_pool, intra_pool_comm) ! #endif ! ! temporary storage of k points and weights ! if (lgamma) then allocate( xkf_ (3, nkstotf), wkf_ (nkstotf) ) else allocate( xkf_ (3, 2*nkstotf), wkf_ (2*nkstotf) ) endif ! #ifdef __PARA if (mpime.eq.ionode_id) then #endif if (lgamma) then ! do ik = 1, nkstotf read (iunkkf, *) (xkf_ (ipol, ik), ipol=1,3), wkf_ (ik) enddo ! else ! ! bring q in crystal coordinates ! xxq = xq call cryst_to_cart (1,xxq,at,-1) ! ! read k and generate k+q in crystal coordinates ! do ik = 1, nkstotf ! ikk = 2 * ik - 1 ikq = ikk + 1 ! read (iunkkf, *) (xkf_ (ipol, ikk ), ipol=1,3), wkf_ (ikk) xkf_ (:, ikq) = xkf_ (:, ikk) + xxq(:) wkf_ ( ikq ) = 0.d0 ! enddo ! ! redefine nkstotf to include the k+q points ! nkstotf = 2 * nkstotf ! endif close ( iunkkf ) ! ! check the sum of the kpoint weights: ! must be 2 for spin-unpolarized calculations ! if (abs(sum( wkf_ )-2.d0).gt.eps) call errore & ('ephwann','wrong weights in kpoint grid',1) ! #ifdef __PARA endif ! call mp_bcast (nkstotf, ionode_id, inter_pool_comm) call mp_bcast (nkstotf, root_pool, intra_pool_comm) ! ! scatter the k points of the fine mesh across the pools ! nksf = kunit * ( nkstotf / kunit / npool ) rest = ( nkstotf - nksf * npool ) / kunit if ( mypool .le. rest ) nksf = nksf + kunit nksqf = nksf / kunit ! nsize = 3 call poolscatter (nsize, nkstotf, xkf_ , nksf, xkf_ ) nsize = 1 call poolscatter (nsize, nkstotf, wkf_ , nksf, wkf_ ) ! write( stdout, '(14x,"max number of k points per pool:",7x,i10)' ) nksf ! nkbl = nkstot / kunit #else if (lgamma) then kunit = 1 else kunit = 2 endif nkbl = nkstot / kunit nksqf = nkbl #endif ! ! now copy the kpoints and weights into smaller arrays ! and deallocate the bigger ones ! allocate ( xkf ( 3, nksf ), wkf( nksf ) ) xkf = xkf_ ( :, 1:nksf ) wkf = wkf_ ( 1:nksf ) deallocate ( xkf_ , wkf_ ) ! !========================================================================= write(stdout, '(14x,a)' ) 'load U on coarse mesh and distribute' !========================================================================= ! ! modify default (in openfilq.f90 we did not know nbnd yet) if (nbndsub.eq.0) nbndsub = nbnd ! if (nbndsub.ne.nbnd) & write(stdout, '(14x,a)' ) 'band disentanglement is used' ! allocate ( cu ( nbnd, nbndsub, nkstot) ) ! #ifdef __PARA if (mpime.eq.ionode_id) then #endif ! ! first proc read rotation matrix (coarse mesh) from file ! open ( unit = iunukk, file = filukk, status = 'old', form = 'formatted') if ((.not.tshuffle).and.(.not.lgamma)) & open ( unit = iunukq, file = filukq, status = 'old', form = 'formatted') ! do ik = 1, nkbl ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! do ibnd = 1, nbnd do jbnd = 1, nbndsub read(iunukk, *) cu (ibnd, jbnd, ikk) if ((.not.tshuffle).and.(.not.lgamma)) & read(iunukq, *) cu (ibnd, jbnd, ikq) enddo enddo enddo ! close ( iunukk ) if ((.not.tshuffle).and.(.not.lgamma)) close ( iunukq ) ! ! in shuffle mode generate U(k+q) through the map ! if (tshuffle.and.(.not.lgamma)) then ! ! first read map {k+q}->{k} ! open ( unit = iunkmap, file = 'kmap.dat', status = 'old', form = 'formatted') do ik = 1, nkbl read(iunkmap,*) ik0, kmap( ik0 ) enddo close ( iunkmap ) ! ! then generate the matrix for the q-displaced mesh ! do ik = 1, nkbl ikq = 2*ik ikk = 2*kmap(ik)-1 cu (:, :, ikq) = cu (:, :, ikk ) enddo ! endif ! #ifdef __PARA endif ! ! scatter the rotation matrix acroos the pools ! allocate ( aux( 2*nbnd*nbndsub, nkstot) ) ! ! first pack the array... ! nsize = 2 * nbnd * nbndsub do ik = 1, nkstot ind = 0 do ibnd = 1, nbnd do jbnd = 1, nbndsub ind = ind + 1 aux( ind, ik) = real ( cu (ibnd, jbnd, ik) ) ind = ind + 1 aux( ind, ik) = aimag ( cu (ibnd, jbnd, ik) ) if (ind.gt.nsize) call errore ('ephwann','problem with packing of cu',1) enddo enddo enddo ! ! ... then scatter ! call poolscatter (nsize, nkstot, aux, nks, aux ) ! ! ... finally unpack ! do ik = 1, nks ind = 0 do ibnd = 1, nbnd do jbnd = 1, nbndsub ind = ind + 1 cu (ibnd, jbnd, ik) = aux( 2*ind-1, ik) + ci * aux( 2*ind, ik) enddo enddo enddo ! deallocate ( aux ) ! #endif ! call stop_clock ( 'load data' ) ! !========================================================================= write(stdout, '(14x,a)' ) 'generate wigner-seitz points' !========================================================================= ! ! the positions of the unit cells in real space where the wannier ! functions will be sitting ! allocate ( irvec (3, 2*nk1*nk2*nk3), ndegen (2*nk1*nk2*nk3), & wslen(2*nk1*nk2*nk3) ) ! call wigner_seitz ( nk1, nk2, nk3, irvec, nrr, ndegen, wslen) ! !==================================================================== write(stdout, '(14x,a)' ) 'interpolation of the Hamiltonian' !==================================================================== ! call start_clock ( 'Ham: step 1' ) ! !---------------------------------------------------------- ! STEP 1: rotation to optimally smooth Bloch states !---------------------------------------------------------- ! ! H~_k = U_k^\dagger H_k U_k ! H~_k+q = U_k+q^\dagger H_k+q U_k+q ! ! H~_k is chs( nbndsub, nbndsub, 2*ik-1 ) ! H~_k+q is chs( nbndsub, nbndsub, 2*ik ) ! allocate ( chs ( nbndsub, nbndsub, nks) ) ! do ik = 1, nks ! do jbnd = 1, nbndsub do ibnd = 1, jbnd ! chs ( ibnd, jbnd, ik ) = czero ! do mbnd = 1, nbnd chs ( ibnd, jbnd, ik) = chs ( ibnd, jbnd, ik) + & conjg( cu ( mbnd, ibnd, ik)) * et ( mbnd, ik) * cu ( mbnd, jbnd, ik) enddo ! if (ibnd .lt. jbnd ) & chs (jbnd , ibnd , ik) = conjg (chs ( ibnd, jbnd, ik)) ! enddo enddo enddo ! call stop_clock ( 'Ham: step 1' ) ! !---------------------------------------------------------- ! STEP 2: Fourier transform to go into Wannier basis !---------------------------------------------------------- ! ! H_k (R) = (1/nk) sum_k e^{-ikR } H~_k(k ) ! H_k+q (R) = (1/nk) sum_{k+q} e^{-i(k+q)R} H~_q(k+q) ! ! chw (nbndsub, nbndsub, 2 * ir -1 ) is H_k (R) ! chw (nbndsub, nbndsub, 2 * ir ) is H_k+q (R) ! ! because of the periodic gauge, H_k (R) and H_k+q (R) ! should be the same when using shuffe mode (commensurate q) ! call start_clock ( 'Ham: step 2' ) ! allocate ( chw (nbndsub, nbndsub, 2 * nrr ) ) ! ! bring xk in crystal coordinates ! call cryst_to_cart (nks, xk, at, -1) ! chw ( :, :, :) = czero ! do ir = 1, nrr ! ! this convention is used both for q = 0 and q /= 0 irkk = 2 * ir - 1 irkq = irkk + 1 ! do ik = 1, nksq ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! rdotk = twopi * dot_product( xk ( :, ikk), float(irvec( :, ir) )) cfac = exp( -ci*rdotk ) / float(nkbl) chw ( :, :, irkk ) = chw ( :, :, irkk ) + cfac * chs ( :, :, ikk ) ! rdotk = twopi * dot_product( xk ( :, ikq), float(irvec( :, ir) )) cfac = exp( -ci*rdotk ) / float(nkbl) chw ( :, :, irkq ) = chw ( :, :, irkq ) + cfac * chs ( :, :, ikq ) ! enddo ! enddo #ifdef __PARA call poolreduce (2 * 2 * nrr * nbndsub * nbndsub , chw ) #endif ! ! bring xk back into cart coord ! call cryst_to_cart (nks, xk, bg, 1) ! deallocate ( chs ) ! if (iverbosity.eq.1) then ! ! check spatial decay of Hamiltonian in Wannier basis ! the unit in r-space is angstrom, and I am plotting ! the matrix for the first mode only ! #ifdef __PARA if (mpime.eq.ionode_id) then #endif if (imode0.eq.0) then write(300, '(/3x,a/)') & 'Spatial decay of Hamiltonian in Wannier basis' do ir =1,nrr ! irkk = 2 * ir - 1 ! tmp = maxval ( abs( chw (:,:,irkk)) ) write(300, *) wslen(ir) * celldm (1) * bohr2ang, tmp ! enddo endif #ifdef __PARA endif call mp_barrier() #endif endif ! call stop_clock ( 'Ham: step 2' ) ! ! Now inverse Fourier transform to fine k and k+q meshes ! and diagonalize. Every pool works with its own subset ! of k points on the fine grid. ! allocate ( chf ( nbndsub, nbndsub), cuf ( nbndsub, nbndsub), etf ( nbndsub, nksf) ) ! ! temporary stuff for lapack zhpevx ! allocate ( ifail( nbndsub ), iwork( 5*nbndsub ), w( nbndsub ), & rwork( 7*nbndsub ), champ( nbndsub*(nbndsub+1)/2 ), & cwork( 2*nbndsub ), cz( nbndsub, nbndsub) ) ! call start_clock ( 'Ham: step 3' ) ! do ik = 1, nksf ! !---------------------------------------------------------- ! STEP 3: inverse Fourier transform to fine k and k+q meshes !---------------------------------------------------------- ! ! H~_k' = sum_R 1/ndegen(R) e^{-ik'R } H_k(R) ! H~_k'+q = sum_R 1/ndegen(R) e^{-i(k'+q)R} H_k+q(R) ! ! H~_k is chf ( nbndsub, nbndsub, 2*ik-1 ) ! H~_k+q is chf ( nbndsub, nbndsub, 2*ik ) ! chf (:,:) = czero ! do ir = 1, nrr ! irkk = 2 * ir - 1 irkq = irkk + 1 ! if (lgamma) then irk = irkk else if (2*(ik/2).eq.ik-1) then ! this is a k point irk = irkk else ! this is a k+q point irk = irkq endif endif ! rdotk = twopi * dot_product( xkf ( :, ik), float(irvec( :, ir) )) cfac = exp( ci*rdotk ) / float( ndegen(ir) ) chf (:,:) = chf (:,:) + cfac * chw (:,:, irk ) ! enddo ! !--------------------------------------------------------------------- ! STEP 4: diagonalize smooth Hamiltonian on k points of the fine grid !--------------------------------------------------------------------- ! ! champ: complex hamiltonian packed (upper triangular part for zhpevx) ! do jbnd = 1, nbndsub do ibnd = 1, jbnd champ (ibnd + (jbnd - 1) * jbnd/2 ) = chf ( ibnd, jbnd) enddo enddo ! call zhpevx ('V', 'A', 'U', nbndsub, champ , 0.0, 0.0, & 0, 0,-1.0, neig, w, cz, nbndsub, cwork, & rwork, iwork, ifail, info) ! cuf = conjg( transpose ( cz ) ) ! ! eigenvalues in Ry (mind when comparing with wannier code) ! etf (:, ik) = w(:) ! ! temporary storage to file ! call dasmio ( cuf, nbndsub, lrcuf, iuncuf, ik, +1) ! enddo ! call stop_clock ( 'Ham: step 3' ) ! deallocate ( chw, chf, cuf ) deallocate ( ifail, iwork, w, rwork, champ, cwork, cz ) ! !------------------------------------------------------------------------- ! endif imode0.eq.0 endif ! write(stdout,'(14x,a,i3)') & !================================================================================== 'interp e-p matrix elements for pert # ', imode0 + 1 !================================================================================== ! !---------------------------------------------------------- ! STEP 1: rotation to optimally smooth Bloch states !---------------------------------------------------------- ! ! g~ = U_k+q^\dagger g U_k ! ! g is el_ph_mat (ibnd, jbnd, ik, imode0+1) ! g~ is epmats (ibnd, jbnd, ik) ! call start_clock ( 'ep: step 1' ) ! allocate ( epmats (nbndsub, nbndsub, nksq), eptmp(nbndsub, nbnd) ) ! do ik = 1, nksq ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! ! the two zgemm calls perform the following ops: ! epmats = [ cu(ikq)^\dagger * el_ph_mat ] * cu(ikk) ! [here we have a size-reduction from nbnd*nbnd to nbndsub*nbndsub] ! call zgemm ('c', 'n', nbndsub, nbnd, nbnd, cone, cu(:,:,ikq), & nbnd, el_ph_mat(:,:,ik,imode0+1), nbnd, czero, eptmp, nbndsub) call zgemm ('n', 'n', nbndsub, nbndsub, nbnd, cone, eptmp, & nbndsub, cu(:,:,ikk), nbnd, czero, epmats(:,:,ik), nbndsub) ! enddo ! call stop_clock ( 'ep: step 1' ) ! !---------------------------------------------------------------------- ! STEP 2: Fourier transform to obtain matrix elements in wannier basis !---------------------------------------------------------------------- ! ! g (R) = (1/nkc) sum_k e^{-ikR} g~(k) ! ! epmatw (nbndsub,nbndsub,ir) is g(R) ! call start_clock ( 'ep: step 2' ) ! allocate ( epmatw (nbndsub, nbndsub, nrr) ) epmatw (:, :, :) = czero ! ! bring xk in crystal coordinates ! call cryst_to_cart (nks, xk, at, -1) ! do ir = 1, nrr ! do ik = 1, nksq ! if (lgamma) then ikk = ik else ikk = 2 * ik - 1 endif ! rdotk = twopi * dot_product( xk ( :, ikk), float(irvec( :, ir) )) cfac = exp( -ci*rdotk ) / float(nkbl) epmatw ( :, :, ir) = epmatw ( :, :, ir) + cfac * epmats ( :, :, ik) ! enddo ! enddo ! #ifdef __PARA call poolreduce (2 * nrr * nbndsub * nbndsub, epmatw ) #endif ! ! bring xk back into cart coord ! call cryst_to_cart (nks, xk, bg, 1) ! deallocate ( epmats ) ! !@ if (iverbosity.eq.1) then ! ! check spatial decay of matrix elements in Wannier basis ! the unit in r-space is angstrom, and I am plotting ! the matrix for the first mode only ! #ifdef __PARA if (mpime.eq.ionode_id) then #endif if (imode0.eq.0) then !@ write(301, '(/3x,a/)') & !@ 'Spatial decay of e-p matrix elements in Wannier basis' do ir =1,nrr ! tmp = maxval ( abs(epmatw(:,:,ir)) ) write(301, *) wslen(ir) * celldm (1) * bohr2ang, tmp ! enddo endif #ifdef __PARA endif #endif !@ endif ! call stop_clock ( 'ep: step 2' ) ! allocate ( epmatf ( nbndsub, nbndsub), & cufkk ( nbndsub, nbndsub), & cufkq ( nbndsub, nbndsub)) ! call start_clock ( 'ep: step 3' ) ! if (imode0.eq.0) then ! ! the following is to slim down the electron-phonon matrix ! for read/write by keeping only the bands within fsthick from ! the Fermi level. We always rotate the full el_ph_mat but ! we only write el_ph_mat(ibnmin:ibndmax,ibndmin:ibndmax) ! and subsequently read in elphsum3 (and elphsum3_strict) ! ibndmin = 100000 ibndmax = 0 ebndmin = 1.d8 ebndmax = -1.d8 ! do ik = 1, nksf do ibnd = 1, nbndsub ebnd = etf (ibnd, ik) ! if ( abs(ebnd - ef) .lt. fsthick ) then ibndmin = min(ibnd,ibndmin) ibndmax = max(ibnd,ibndmax) ebndmin = min(ebnd,ebndmin) ebndmax = max(ebnd,ebndmax) endif ! enddo enddo ! #ifdef __PARA tmp = float (ibndmin) call poolextreme( tmp, -1 ) ibndmin = nint (tmp) call poolextreme( ebndmin, -1 ) ! tmp = float (ibndmax) call poolextreme( tmp, +1 ) ibndmax = nint (tmp) call poolextreme( ebndmax, +1 ) #endif ! if (iverbosity.eq.1) then write(stdout,'(/14x,a,i5,2x,a,f9.3)') 'ibndmin = ', ibndmin, 'ebndmin = ', ebndmin write(stdout,'(14x,a,i5,2x,a,f9.3/)') 'ibndmax = ', ibndmax, 'ebndmax = ', ebndmax endif ! ! now 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) ! endif ! ! counter for the number of k points with a band falling on the Fermi surface fermicount = 0 ! do ik = 1, nksqf ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! ! perform interpolation only for those pair of states both within ! a Fermi shell of size fsthick (should give a significant speedup) ! tfermikk = .false. tfermikq = .false. ! do ibnd = 1, nbndsub if ( abs(etf (ibnd, ikk) - ef) .lt. fsthick ) tfermikk = .true. if ( abs(etf (ibnd, ikq) - ef) .lt. fsthick ) tfermikq = .true. enddo ! if (tfermikk.and.tfermikq) then ! fermicount = fermicount + 1 ! !---------------------------------------------------------- ! STEP 3: inverse Fourier transform of g to fine k mesh !---------------------------------------------------------- ! ! g~ (k') = sum_R 1/ndegen(R) e^{-ik'R} g (R) ! ! g~(k') is epmatf (nbndsub, nbndsub, ik ) ! every pool works with its own subset of k points on the fine grid ! epmatf ( :, :) = czero ! do ir = 1, nrr ! ! note xkf is assumed to be already in cryst coord ! rdotk = twopi * dot_product ( xkf(:, ikk), float(irvec(:, ir)) ) cfac = exp( ci*rdotk ) / float( ndegen(ir) ) ! epmatf (:, :) = epmatf (:, :) + cfac * epmatw ( :, :, ir) ! enddo ! !---------------------------------------------------------- ! STEP 4: un-rotate to Bloch space, fine grid !---------------------------------------------------------- ! ! g (k') = U_q^\dagger (k') g~ (k') U_k (k') ! ! read cufkk and cufkq from file ! call dasmio ( cufkk, nbndsub, lrcuf, iuncuf, ikk, -1) call dasmio ( cufkq, nbndsub, lrcuf, iuncuf, ikq, -1) ! ! the two zgemm calls perform the following ops: ! epmatf = [ cufkq * epmatf ] * cufkk^\dagger ! call zgemm ('n', 'n', nbndsub, nbndsub, nbndsub, cone, cufkq, & nbndsub, epmatf, nbndsub, czero, eptmp, nbndsub) call zgemm ('n', 'c', nbndsub, nbndsub, nbndsub, cone, eptmp, & nbndsub, cufkk, nbndsub, czero, epmatf, nbndsub) ! ! write epmatf to file, mind the outer loop on modes! ! nrec = imode0 * nksqf + ik call dasmio ( epmatf(ibndmin:ibndmax,ibndmin:ibndmax), & ibndmax-ibndmin+1, lrepmatf, iunepmatf, nrec, +1) ! endif ! enddo ! call stop_clock ( 'ep: step 3' ) ! if (imode0+1.eq.3*nat) then write( stdout, '(/14x,a,i5,a,i5)' ) & 'Number of (k,k+q) pairs on the Fermi surface: ', & fermicount, ' out of ', nksqf endif ! deallocate ( epmatw, epmatf, cufkk, cufkq, eptmp ) ! !============================================================== ! call stop_clock ( 'ephwann' ) ! return end subroutine ephwann ! !-------------------------------------------------------------- subroutine dasmio ( mat, nsize, lrec, iun, nrec, iop ) !-------------------------------------------------------------- ! ! A simple wrapper to the davcio routine to read/write square ! matrices instead of vectors (Direc Access Square Matrix I/O) ! by @ FG ! USE kinds, only : DP implicit none integer :: nsize, lrec, iun, nrec, iop, i complex(kind=DP):: mat(nsize,nsize), aux ( nsize*nsize ) ! if ( iop .eq. -1 ) then ! ! read matrix ! call davcio ( aux, lrec, iun, nrec, -1 ) do i = 1, nsize * nsize mat ( 1 + (i-1)/nsize , i - nsize*((i-1)/nsize) ) = aux (i) enddo ! elseif ( iop .eq. 1 ) then ! ! write matrix ! do i = 1, nsize * nsize aux (i) = mat ( 1 + (i-1)/nsize , i - nsize*((i-1)/nsize) ) enddo call davcio ( aux, lrec, iun, nrec, +1 ) ! else ! call errore ('dasmio','iop not permitted',1) ! endif ! return end subroutine dasmio !