!--------------------------------------------------------------------- subroutine wanout !--------------------------------------------------------------------- ! ! Here we collect information for Wannier/Fourier ! interpolation of the electron-phonon matrix elements ! ! Parallelizazion on G-vectors and k-vectors on April 26 2006 ! No ultrasoft now ! Works for or Intel v9, IBM SP4, and Opteron with pgi ! ! Feliciano Giustino, University of California at Berkeley ! !--------------------------------------------------------------------- ! #include "f_defs.h" ! USE para, ONLY : kunit #ifdef __PARA use para USE mp_global, ONLY : nproc, my_pool_id, nproc_pool, & intra_image_comm, inter_pool_comm, me_pool, & root_pool, mpime USE mp, ONLY : mp_barrier, mp_bcast #endif USE kinds, only : DP use io_global, only : stdout, ionode_id USE wavefunctions_module, ONLY: evc USE io_files, ONLY: iunigk use wvfct, only : npwx use pwcom use phcom use el_phon USE wvfct, ONLY : et USE ions_base, ONLY : atm, tau, nat, ityp, ntyp => nsp ! implicit none ! ! work variables ! real(kind=DP) :: et_g( nbnd, nkstot ) integer :: ik, ikk, ikq, ibnd, jbnd, ir, ig, jg, ng, & ios, nkl, nkr, ik0, iks complex(kind=DP) , allocatable :: aux1 (:), aux2 (:), aux2_all(:) complex(kind=DP) :: ZDOTC ! ! the output files ! integer :: & iuwandat, & ! unit for wannier input file iukptskc, & ! unit for k- grid points iukptsqc, & ! unit for k+q- grid points iubandkc, & ! unit for energy eigenvalues on k- grid iubandqc, & ! unit for energy eigenvalues on k+q- grid iuwfwank, & ! unit for wavefunctions on k- grid iuwfwanq ! unit for wavefunctions on k+q- grid ! character(len=30) wfnname_kk, wfnname_kq, filbandk, & filbandq, filkptsk, filkptsq real(kind=DP), parameter :: ryd2ev = 13.60569172 real(kind = 8), parameter :: bohr2ang = 0.5291772108 complex(kind = 8), parameter :: czero = (0.d0, 0.d0), cone = (1.d0, 0.d0) integer :: ia, ix, it real(kind=DP) :: xktmp(3) real(kind=DP), allocatable :: rat(:,:,:) integer, allocatable :: natom(:) ! ! work variables ! integer :: ig0, ik1, igkq_tmp (npwx), imap, ishift ! allocate (rat( 3, nat, ntyp )) allocate (natom( ntyp )) allocate (aux1( nrxxs), aux2( nrxxs)) #ifdef __PARA allocate ( aux2_all( nrx1s*nrx2s*nrx3s) ) if (tshuffle.and.nproc_pool>1) call errore & ('wanout', 'only one proc per pool in shuffle mode', 1) #endif if (tshuffle.and.lgamma) call errore & ('wanout','shuffle mode cannot be used with q = 0',1) ! ! ... Open files for wannier code ! ! the file for the el-ph matrix is opened outside ! (in elphon) because of the loop on irreps ! ! here the files written by the first proc of every pool iuwfwank = 55 iuwfwanq = 56 ! ! here the files written by the first proc of the first pool iuwandat = 50 iubandkc = 51 iubandqc = 52 iukptskc = 53 iukptsqc = 54 ! ik0 = 0 ! #ifdef __PARA ! ! output on /dev/null for silent nodes ! if (me.ne.1) then iuwfwank = stdout iuwfwanq = stdout endif if ((me.ne.1).or.(mypool.ne.1)) then iubandkc = stdout iubandqc = stdout iukptskc = stdout iukptsqc = stdout iuwandat = stdout endif ! ! number of kpoint blocks, kpoints per pool and reminder ! npool = nproc/nproc_pool nkbl = nkstot / kunit nkl = kunit * ( nkbl / npool ) nkr = ( nkstot - nkl * npool ) / kunit ! ! the reminder goes to the first nkr pools ! (0...nkr-1) ! if ( my_pool_id < nkr ) nkl = nkl + kunit ! ! the index of the first k point in this pool ! iks = nkl * my_pool_id + 1 if ( my_pool_id >= nkr ) iks = iks + nkr * kunit ! ! the index of the first k point block in this pool - 1 ! (I will need the index of ik, not ikk) ! ik0 = ( iks - 1 ) / kunit ! #else if (lgamma) then kunit = 1 else kunit = 2 endif nkbl = nkstot / kunit #endif ! #ifdef __PARA if ((me.eq.1).and.(mypool.eq.1)) then #endif open (unit = iuwandat, file = 'wannier.data', form = 'formatted') open (unit = iukptskc, file = 'kpts_kc.dat' , form = 'formatted') open (unit = iukptsqc, file = 'kpts_qc.dat' , form = 'formatted') open (unit = iubandkc, file = 'BAND_kc.dat' , form = 'formatted') open (unit = iubandqc, file = 'BAND_qc.dat' , form = 'formatted') #ifdef __PARA endif #endif ! ! build input file WANNIER.DATA for wannier code ! natom = 0 do ia = 1, nat natom(ityp(ia)) = natom(ityp(ia)) + 1 rat(:,natom(ityp(ia)),ityp(ia)) = tau(:,ia) call cryst_to_cart (1,rat(:,natom(ityp(ia)),ityp(ia)),bg,-1) end do ! write (iuwandat,*) 'Crystal basis vectors, in unit of alat = ', & alat*bohr2ang, '(Ang)' write (iuwandat,'(3f18.10)') (at(ix,1),ix=1,3) write (iuwandat,'(3f18.10)') (at(ix,2),ix=1,3) write (iuwandat,'(3f18.10)') (at(ix,3),ix=1,3) ! write (iuwandat,*) 'Atomic positions ( lattice coordinates )' do it = 1, ntyp write (iuwandat,'(a,i3)') atm(it), natom(it) do ia = 1, natom(it) write (iuwandat,'(3f18.10)') (rat(ix,ia,it),ix=1,3) enddo end do ! write (iuwandat,*) 'Atomic positions ( From PW )' do ia = 1, nat write (iuwandat, 3) ia,(tau(ix,ia),ix=1,3) end do ! write (iuwandat,*) 'ecut= ', ecutwfc * ryd2ev, '(eV)' write (iuwandat,*) 'Grid ( nr1s, nr2s, nr3s, ngm_g ) = ' ! global number of G-vectors (from module pwcom) write (iuwandat,*) nr1s,nr2s,nr3s,ngm_g ! ! build input files KPTS.DAT and BAND.DAT for wannier code ! (eigenvalues are printed out in eV) ! note that the first processor has already a copy of all kpoints ! while for the eigenenergies we need to collect them across the pools ! et_g = et #ifdef __PARA CALL poolrecover( et_g, nbnd, nkstot, nks ) #endif ! do ik = 1, nkbl if (lgamma) then ikk = ik ikq = ik npwq = npw else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! xktmp = xk(:,ikk) call cryst_to_cart (1, xktmp, at, -1) write(iukptskc,'(4f20.10)') (xktmp(ix),ix=1,3), wk(ikk) ! ! this is to print out directly the k+q rather than the folded ones xktmp = xk(:,ikk) xktmp = xktmp + xqq call cryst_to_cart (1, xktmp, at, -1) ! write(iukptsqc,'(4f20.10)') (xktmp(ix),ix=1,3), wk(ikk) ! do ibnd = 1, nbnd write(iubandkc,*) ibnd, ik, ryd2ev*et_g( ibnd, ikk ) write(iubandqc,*) ibnd, ik, ryd2ev*et_g( ibnd, ikq ) enddo ! enddo ! #ifdef __PARA if ((me.eq.1).and.(mypool.eq.1)) then #endif close (iuwandat) close (iubandkc) close (iubandqc) close (iukptskc) close (iukptsqc) #ifdef __PARA endif #endif ! ! loop over the k-points for writing wfs ! if (nksq.gt.1) rewind (unit = iunigk) do ik = 1, nksq ! if (nksq.gt.1) then read (iunigk, err = 100, iostat = ios) npw, igk 100 call errore ('elphel', 'reading igk', abs (ios) ) endif ! ! ik = counter of k-points with vector k ! ikk= index of k-point with vector k ! ikq= index of k-point with vector k+q ! k and k+q are alternated if q!=0, are the same if q=0 ! if (lgamma) then ikk = ik ikq = ik npwq = npw else ikk = 2 * ik - 1 ikq = ikk + 1 endif if (lsda) current_spin = isk (ikk) if (.not.lgamma.and.nksq.gt.1) then read (iunigk, err = 200, iostat = ios) npwq, igkq 200 call errore ('wanout', 'reading igkq', abs (ios) ) endif ! ! read unperturbed wavefuctions psi(k) and psi(k+q) ! if (nksq.gt.1) then if (lgamma) then call davcio (evc, lrwfc, iuwfc, ikk, - 1) else call davcio (evc, lrwfc, iuwfc, ikk, - 1) call davcio (evq, lrwfc, iuwfc, ikq, - 1) endif endif ! ! ! ---------------------------------------------------------------- ! Set the gauge for the eigenstates: unitary transform and phases ! ---------------------------------------------------------------- ! ! With this option, different compilers and different machines ! should always give the same wavefunctions. ! if ( tphases ) then ! if (ik.eq.1) allocate ( umat(nbnd,nbnd,nks) ) if (lgamma) then call setphases ( kunit, ikk, npw, umat(:,:,ikk)) else call setphases ( kunit, ikk, npw, umat(:,:,ikk)) call setphases ( kunit, ikq, npwq,umat(:,:,ikq)) endif ! endif ! ! -------------------------------------------------- ! build input files WFN[ik].[ispin] for wannier code ! -------------------------------------------------- ! write(wfnname_kk, 1) ik + ik0, current_spin write(wfnname_kq, 2) ik + ik0, current_spin ! #ifdef __PARA if (me.eq.1) then #endif open (unit = iuwfwank, file = wfnname_kk, form = 'unformatted') open (unit = iuwfwanq, file = wfnname_kq, form = 'unformatted') #ifdef __PARA endif #endif ! write (iuwfwank) ik + ik0, nbnd write (iuwfwanq) ik + ik0, nbnd ! do ibnd = 1, nbnd ! aux2(:) = (0.d0, 0.d0) do ig = 1, npw ! if (tphases) then ! ! rotate before writing ! do jbnd = 1, nbnd aux2 (nls (igk (ig) ) ) = aux2 (nls (igk (ig) ) ) & + umat (jbnd, ibnd, ikk) * evc (ig,jbnd) enddo ! else ! ! or leave states unchanged ! aux2 (nls (igk (ig) ) ) = evc (ig, ibnd) endif ! enddo call cft3s (aux2, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2) ! #ifdef __PARA if (nproc_pool.gt.1) then ! aux2_all(:) = (0.d0, 0.d0) call cgather_sym( aux2, aux2_all) ! write(iuwfwank) ( aux2_all(ir), ir = 1, nrx1s*nrx2s*nrx3s ) ! else #endif write(iuwfwank) ( aux2(ir), ir = 1, nrx1s*nrx2s*nrx3s ) #ifdef __PARA endif #endif if ((.not.tshuffle).and.(.not.lgamma)) then ! aux2(:) = (0.d0, 0.d0) do ig = 1, npwq ! if (tphases) then ! ! rotate before writing ! do jbnd = 1, nbnd aux2 (nls (igkq (ig) ) ) = aux2 (nls (igkq (ig) ) ) & + umat (jbnd, ibnd, ikq) * evq (ig,jbnd) enddo ! else ! ! or leave states unchanged ! aux2 (nls (igkq (ig) ) ) = evq (ig, ibnd) endif ! enddo call cft3s (aux2, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2) #ifdef __PARA if (nproc_pool.gt.1) then ! aux2_all(:) = (0.d0, 0.d0) call cgather_sym( aux2, aux2_all) ! write(iuwfwanq) ( aux2_all(ir), ir = 1, nrx1s*nrx2s*nrx3s ) ! else #endif write(iuwfwanq) ( aux2(ir), ir = 1, nrx1s*nrx2s*nrx3s ) #ifdef __PARA endif #endif endif end do #ifdef __PARA if (me.eq.1) then #endif close (iuwfwank) close (iuwfwanq) #ifdef __PARA endif #endif enddo ! deallocate (aux1, aux2) #ifdef __PARA deallocate (aux2_all) #endif ! 1 format('kc/WFN',i5.5,'.',i1) 2 format('qc/WFN',i5.5,'.',i1) 3 format(' tau( ',i3,' )',3f10.6) ! return end subroutine wanout !------------------------------------------------------------ !