TABLE OF CONTENTS
ABINIT/qijb_kk [ Functions ]
NAME
qijb_kk
FUNCTION
Routine which computes PAW onsite part of wavefunction overlap for Bloch functions at two k-points k and k+b. These quantities are used in PAW-based computations of polarization and magnetization.
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
dkvecs(3) :: $\Delta k$ input vector expibi(2,my_natom,3) :: phase factors at each atomic site for given k offset gprimd(3,3)=dimensioned primitive translations of reciprocal lattice lmn2max :: lmnmax*(lmnmax+1)/2 natom=number of atoms in unit cell ntypat=number of types of atoms in unit cell pawang <type(pawang_type)>=paw angular mesh and related data pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data typat=typat(natom) list of atom types
OUTPUT
calc_qijb(2,lmn2max,natom) :: PAW on-site overlaps of wavefunctions at neighboring
k point
SIDE EFFECTS
NOTES
this function computes the on-site data for the PAW version of <u_nk|u_mk+b>, that is, two Bloch vectors at two different k points.
PARENTS
initberry,overlap_k1k2_paw
CHILDREN
initylmr,sbf8,simp_gen
SOURCE
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 subroutine qijb_kk(calc_qijb,dkvecs,expibi,gprimd,lmn2max,natom,ntypat,& 54 & pawang,pawrad,pawtab,typat) 55 56 use m_profiling_abi 57 58 use defs_basis 59 use m_errors 60 61 use m_xmpi, only : xmpi_sum 62 use m_special_funcs, only : sbf8 63 use m_pawang, only : pawang_type 64 use m_pawrad, only : pawrad_type, simp_gen 65 use m_pawtab, only : pawtab_type 66 use m_paw_sphharm, only : initylmr 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'qijb_kk' 72 !End of the abilint section 73 74 implicit none 75 76 !Arguments--------------------------- 77 !scalars 78 integer,intent(in) :: lmn2max,natom,ntypat 79 type(pawang_type),intent(in) :: pawang 80 real(dp),intent(out) :: calc_qijb(2,lmn2max,natom) 81 !arrays 82 integer,intent(in) :: typat(natom) 83 real(dp),intent(in) :: dkvecs(3),expibi(2,natom),gprimd(3,3) 84 type(pawrad_type),intent(in) :: pawrad(ntypat) 85 type(pawtab_type),intent(in) :: pawtab(ntypat) 86 87 !Local variables--------------------------- 88 !scalars 89 integer :: iatom,ir,isel,itypat 90 integer :: klm,kln,klmn,lbess,lbesslm,lmin,lmax,mbess,mesh_size 91 integer :: ylmr_normchoice,ylmr_npts,ylmr_option 92 real(dp) :: arg,bessg,bnorm,intg,rterm 93 complex(dpc) :: cterm,etb,ifac 94 !arrays 95 real(dp) :: bb(3),bbn(3),bcart(3),ylmgr(1,1,0),ylmr_nrm(1) 96 real(dp),allocatable :: ff(:),j_bessel(:,:),ylmb(:),sb_out(:) 97 ! the following is (i)^L mod 4. 98 complex(dpc),dimension(0:3) :: il(0:3)=(/cone,j_dpc,-cone,-j_dpc/) 99 100 ! ************************************************************************* 101 102 calc_qijb(:,:,:) = zero 103 104 ylmr_normchoice = 0 ! input to initylmr are normalized 105 ylmr_npts = 1 ! only 1 point to compute in initylmr 106 ylmr_nrm(1) = one ! weight of normed point for initylmr 107 ylmr_option = 1 ! compute only ylm's in initylmr 108 109 ABI_ALLOCATE(sb_out, (pawang%l_size_max)) 110 111 do iatom = 1, natom 112 113 itypat = typat(iatom) 114 mesh_size = pawtab(itypat)%mesh_size 115 116 ABI_ALLOCATE(j_bessel,(mesh_size,pawang%l_size_max)) 117 ABI_ALLOCATE(ff,(mesh_size)) 118 ABI_ALLOCATE(ylmb,(pawang%l_size_max*pawang%l_size_max)) 119 120 ! here is exp(-i b.R) for current atom: recall storage in expibi 121 etb = cmplx(expibi(1,iatom),expibi(2,iatom)) 122 123 ! note the definition used for the k-dependence of the PAW basis functions: 124 !$|\phi_{i,k}\rangle = exp(-i k\cdot r)|\phi_i\rangle 125 ! see Umari, Gonze, and Pasquarello, PRB 69,235102 Eq. 23. Thus the k-vector on the 126 ! bra side enters as k, while on the ket side it enters as -k. 127 bb(:) = -dkvecs(:) 128 129 ! reference bb to cartesian axes 130 bcart(1:3)=MATMUL(gprimd(1:3,1:3),bb(1:3)) 131 132 ! bbn is b-hat (the unit vector in the b direction) 133 bnorm=dsqrt(dot_product(bcart,bcart)) 134 bbn(:) = bcart(:)/bnorm 135 136 ! as an argument to the bessel function, need 2pi*b*r = 1 so b is re-normed to two_pi 137 bnorm = two_pi*bnorm 138 do ir=1,mesh_size 139 arg=bnorm*pawrad(itypat)%rad(ir) 140 call sbf8(pawang%l_size_max,arg,sb_out) ! spherical bessel functions at each mesh point 141 j_bessel(ir,:) = sb_out 142 end do ! end loop over mesh 143 144 ! compute Y_LM(b) here 145 call initylmr(pawang%l_size_max,ylmr_normchoice,ylmr_npts,ylmr_nrm,ylmr_option,bbn,ylmb(:),ylmgr) 146 147 do klmn = 1, pawtab(itypat)%lmn2_size 148 klm =pawtab(itypat)%indklmn(1,klmn) 149 kln =pawtab(itypat)%indklmn(2,klmn) 150 lmin=pawtab(itypat)%indklmn(3,klmn) 151 lmax=pawtab(itypat)%indklmn(4,klmn) 152 do lbess = lmin, lmax, 2 ! only possible choices for L s.t. Gaunt integrals 153 ! will be non-zero 154 ifac = il(mod(lbess,4)) 155 do mbess = -lbess, lbess 156 lbesslm = lbess*lbess+lbess+mbess+1 157 isel=pawang%gntselect(lbesslm,klm) 158 if (isel > 0) then 159 bessg = pawang%realgnt(isel) 160 ff(1:mesh_size)=(pawtab(itypat)%phiphj(1:mesh_size,kln)& 161 & -pawtab(itypat)%tphitphj(1:mesh_size,kln))& 162 & *j_bessel(1:mesh_size,lbess+1) 163 call simp_gen(intg,ff,pawrad(itypat)) 164 rterm = four_pi*bessg*intg*ylmb(lbesslm) 165 cterm = etb*ifac*rterm 166 calc_qijb(1,klmn,iatom) = & 167 & calc_qijb(1,klmn,iatom) + dreal(cterm) 168 calc_qijb(2,klmn,iatom) = & 169 & calc_qijb(2,klmn,iatom) + dimag(cterm) 170 171 end if ! end selection on non-zero Gaunt factors 172 end do ! end loop on mbess = -lbess, lbess 173 end do ! end loop on lmin-lmax bessel l values 174 end do ! end loop on lmn2_size klmn basis pairs 175 176 ABI_DEALLOCATE(j_bessel) 177 ABI_DEALLOCATE(ff) 178 ABI_DEALLOCATE(ylmb) 179 end do ! end loop over atoms 180 181 ABI_DEALLOCATE(sb_out) 182 183 end subroutine qijb_kk