! !----------------------------------------------------------------------- subroutine createkmap ( xq ) !----------------------------------------------------------------------- ! ! ! USE kinds, only : DP USE cell_base, ONLY : at, bg USE klist, ONLY : nkstot, nks, xk USE para, ONLY : kunit USE pwcom, ONLY : nk1, nk2, nk3 use io_global, only : stdout use klist, only : kmap use kfold USE control_flags, ONLY : iverbosity USE control_ph, ONLY : lgamma #ifdef __PARA use para, only : me use mp_global, only : my_pool_id, mpime use mp, only : mp_barrier #endif implicit none ! real(kind=DP) :: eps, xq (3) ! integer :: ik, j, nksqtot, ikk, ikq ! ! obsolete: ! ! variables for folding of the k+q mesh into the kmesh @FG ! there may be at most 9 different G_0 to refold k+q into k ! of the first BZ (including G_0=0) ! [actually there are 3^3=27 translations but for a fixed q only ! 2^3=8 out of these 27 are possible since q has a definite direction] ! ! new: ! ! I need all q-vectors int the same run, so I consider all the ! 27 possibilities by setting ng0vec = 27 in ../PW/set_kplusq.f90 ! real(kind=DP) :: xx, yy, zz integer :: i, k, n, g0vec(3), ig0, iukmap, ig1, ig2, ig3 logical :: in_the_list ! ! #ifdef __PARA if (mpime.eq.0) then #endif ! ! the first proc keeps a copy of all kpoints ! ! nksqtot = nkstot / kunit ! if ( .not. allocated (shift) ) allocate ( shift(nksqtot) ) ! ! Now fold k+q back into the k-grid for wannier interpolation. ! Since this is done before divide&impera, every pool has all the kpoints. @FG ! ! bring q in crystal coordinates and check commensuration ! ! loosy tolerance: not important since k+q is defined trhough nint() eps = 1.d-5 ! call cryst_to_cart (1,xq,at,-1) ! xx = xq(1)*nk1 yy = xq(2)*nk2 zz = xq(3)*nk3 in_the_list = abs(xx-nint(xx)).le.eps .and. & abs(yy-nint(yy)).le.eps .and. & abs(zz-nint(zz)).le.eps if (iverbosity.eq.1) & write(stdout,'(a,3i3)') ' q in integer coord:',nint(xx),nint(yy),nint(zz) if (.not.in_the_list) call errore('create_kmap','q-vec not commensurate',1) ! ! bring all the k's in crystal coordinates ! call cryst_to_cart ( nkstot, xk, at, -1) ! !@ ! the 27 possible translations !@ now we have 125 possibilities for LSCO, because the unit cell is far from cubic... sob! !@ ng0vec must be redefined ! ng0vec = 0 !@ do ig1 = -1,1 !@ do ig2 = -1,1 !@ do ig3 = -1,1 do ig1 = -2,2 do ig2 = -2,2 do ig3 = -2,2 ng0vec = ng0vec + 1 g0vec_all(1,ng0vec) = ig1 g0vec_all(2,ng0vec) = ig2 g0vec_all(3,ng0vec) = ig3 enddo enddo enddo ! do ik = 1, nksqtot ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = 2 * ik endif ! ! check that the k's are actually on a uniform mesh centered at gamma ! xx = xk(1, ikk)*nk1 yy = xk(2, ikk)*nk2 zz = xk(3, ikk)*nk3 if (iverbosity.eq.1) & write(stdout,'(a,i3,a,3i3)') 'ik = ',ik,', k in integer coord:',nint(xx),nint(yy),nint(zz) in_the_list = abs(xx-nint(xx)).le.eps .and. & abs(yy-nint(yy)).le.eps .and. & abs(zz-nint(zz)).le.eps if (.not.in_the_list) call errore('createkmap','is this a uniform k-mesh?',1) ! ! now add the phonon wavevector and check that k+q falls again on the k grid ! xk (:, ikq) = xk (:, ikk) + xq (:) ! xx = xk(1, ikq)*nk1 yy = xk(2, ikq)*nk2 zz = xk(3, ikq)*nk3 if (iverbosity.eq.1) & write(stdout,'(a,i3,a,3i3)') 'ik = ',ik,', k+q in integer coord:',nint(xx),nint(yy),nint(zz) in_the_list = abs(xx-nint(xx)).le.eps .and. & abs(yy-nint(yy)).le.eps .and. & abs(zz-nint(zz)).le.eps if (.not.in_the_list) call errore('createkmap','k+q does not fall on k-grid',1) ! ! find the index of this k+q in the k-grid ! i = mod ( nint ( xx + 2*nk1), nk1 ) j = mod ( nint ( yy + 2*nk2), nk2 ) k = mod ( nint ( zz + 2*nk3), nk3 ) if (iverbosity.eq.1) & write(stdout,'(a,i3,a,3i3)') 'ik = ',ik,', f-k+q in integer coord:',i,j,k n = i*nk2*nk3 + j*nk3 + k + 1 kmap( ik ) = n ! ! determine the G_0 such that k+q+G_0 belongs to the first BZ ! g0vec(1) = ( i - nint(xx) )/ nk1 g0vec(2) = ( j - nint(yy) )/ nk2 g0vec(3) = ( k - nint(zz) )/ nk3 if (iverbosity.eq.1) then write(stdout,'(a,i3,a,3i3)') 'ik = ',ik,', G_0 in integer coord:',g0vec(:) write(stdout,'(2i5)') ik, kmap(ik) endif ! ! ! now store the shift for this k+q point ! in_the_list = .false. ig0 = 0 do while ((ig0.le.ng0vec).and.(.not.in_the_list)) ig0 = ig0 + 1 in_the_list = (g0vec(1) .eq. g0vec_all(1,ig0)) .and. & (g0vec(2) .eq. g0vec_all(2,ig0)) .and. & (g0vec(3) .eq. g0vec_all(3,ig0)) enddo shift( ik ) = ig0 ! ! if (.not.in_the_list) call errore & ('createkmap','cannot find the folding vector in the list',1) if (iverbosity.eq.1) then write(stdout,'(a,i3,10x,a,i3)') 'ik = ',ik,'shift code:', shift(ik) write(stdout,*) endif ! ! obsolete: ! ! very important: now redefine k+q through the corresponding kpoint on the k mesh ! Note that this will require using the periodic gauge in the calculation of the ! electron-phonon matrix elements (factor e^iG_0*r if G_0 is the vector used for ! refolding) ! xk (:, ikq) = xk (:, 2 * n - 1) ! enddo ! ! bring everybody back to cartesian coordinates ! call cryst_to_cart ( 1, xq, bg, 1) call cryst_to_cart ( nkstot, xk, bg, 1) ! g0vec_all_r = float ( g0vec_all ) call cryst_to_cart ( ng0vec, g0vec_all_r, bg, 1) ! if (iverbosity.eq.1) then write(stdout,*) write(stdout,*) ' ik, kmap(ik), shift(ik)' do ik = 1, nksqtot write(stdout,'(3i4)') ik, kmap(ik), shift(ik) enddo write(stdout,*) ' There are ',ng0vec,'inequivalent folding G_0 vectors' write(stdout,*) ' crystal coordinates: ' do ig0 = 1,ng0vec write(stdout,'(a,i3,a,3i3)') 'g0vec( ',ig0,') = ',g0vec_all(:,ig0) enddo write(stdout,*) ' cartesian coordinates: ' do ig0 = 1,ng0vec write(stdout,'(3f9.5)') g0vec_all_r(:,ig0) enddo endif ! ! the unit with kmap(ik) and shift(ik) ! iukmap = 97 open ( unit = iukmap, file = 'kmap.dat', form = 'formatted') do ik = 1, nksqtot write (iukmap, '(3i6)') ik, kmap(ik), shift(ik) enddo close ( unit = iukmap ) ! #ifdef __PARA endif call mp_barrier() #endif ! return end subroutine createkmap