TABLE OF CONTENTS
ABINIT/calc_optical_mels [ Functions ]
NAME
calc_optical_mels
FUNCTION
Calculate all optical matrix elements in the BZ.
COPYRIGHT
Copyright (C) 2009-2018 ABINIT group (L.Reining, V.Olevano, F.Sottile, S.Albrecht, G.Onida, M.Giantomassi) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
lomo_spin(Wfd%nsppol)=Index of the lomo band for the different spins. lomo_min,max_band=minimum and max band index to be calculated. nkbz=Number of points in the full Brillouin zone. inclvkb=if different from 0, [Vnl,r] is included in the calculation of the matrix element of the velocity operator No meaning for PAW (except for LDA+U) qpt(3) Kmesh<kmesh_t>=Info on the k-point sampling for wave functions. Cryst<crystal_t>=Structure defining the crystalline structure. KS_Bst<ebands_t> Pawtab(Cryst%ntypat*usepaw)<pawtab_type>=PAW tabulated starting data Psps <pseudopotential_type>=variables related to pseudopotentials. Hur(Cryst%natom*usepaw)<pawhur_t>=Only for PAW and LDA+U, quantities used to evaluate the commutator [H_u,r]. Wfd<wfd_t>=Handler for the wavefunctions. nkibz=Number of irreducible k-points. nsppol=Number of independent spin polarizations. usepaw=1 for PAW, 0 otherwise. nspinor=Number of spinorial components. comm=MPI Communicator.
OUTPUT
opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,nsppol)=Matrix elements <c k|e^{+iqr}|v k>
PARENTS
m_exc_spectra,m_haydock
CHILDREN
get_bz_item,matrginv,pawcprj_alloc,pawcprj_free,vkbr_free,vkbr_init wfd_distribute_bbp,wfd_get_cprj,wrtout,xmpi_barrier,xmpi_sum
SOURCE
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 subroutine calc_optical_mels(Wfd,Kmesh,KS_Bst,Cryst,Psps,Pawtab,Hur,& 54 & inclvkb,lomo_spin,lomo_min,max_band,nkbz,qpoint,opt_cvk) 55 56 use defs_basis 57 use m_errors 58 use m_profiling_abi 59 use m_xmpi 60 61 use defs_datatypes, only : ebands_t, pseudopotential_type 62 use m_abilasi, only : matrginv 63 use m_bz_mesh, only : kmesh_t, get_BZ_item 64 use m_crystal, only : crystal_t 65 use m_vkbr, only : vkbr_t, vkbr_free, vkbr_init, nc_ihr_comm 66 use m_wfd, only : wfd_t, wfd_get_cprj, wfd_distribute_bbp 67 use m_pawtab, only : pawtab_type 68 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 69 use m_pawhr, only : pawhur_t, paw_ihr 70 71 !This section has been created automatically by the script Abilint (TD). 72 !Do not modify the following lines by hand. 73 #undef ABI_FUNC 74 #define ABI_FUNC 'calc_optical_mels' 75 use interfaces_14_hidewrite 76 !End of the abilint section 77 78 implicit none 79 80 !Arguments ------------------------------------ 81 !scalars 82 integer,intent(in) :: nkbz,inclvkb,lomo_min,max_band 83 type(kmesh_t),intent(in) :: Kmesh 84 type(crystal_t),intent(in) :: Cryst 85 type(pseudopotential_type),intent(in) :: Psps 86 type(ebands_t),intent(in) :: KS_Bst 87 type(wfd_t),target,intent(inout) :: Wfd 88 !arrays 89 integer,intent(in) :: lomo_spin(Wfd%nsppol) 90 real(dp),intent(in) :: qpoint(3) 91 complex(dpc),intent(out) :: opt_cvk(lomo_min:max_band,lomo_min:max_band,nkbz,Wfd%nsppol) 92 type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 93 type(pawhur_t),intent(in) :: Hur(Cryst%natom*Wfd%usepaw) 94 95 !Local variables ------------------------------ 96 !scalars 97 integer :: nsppol,usepaw,nspinor,comm,spin,npw_k,istwf_k,my_nbbp 98 integer :: ik_bz,ik_ibz,itim_k,isym_k,ib_c,ib_v,ierr,my_rank 99 real(dp) :: ediff 100 complex(dpc) :: emcvk 101 type(vkbr_t) :: vkbr 102 !arrays 103 integer,allocatable :: bbp_distrb(:,:) 104 integer,ABI_CONTIGUOUS pointer :: kg_k(:,:) 105 real(dp) :: mat_dp(3,3),qrot(3),b1(3),b2(3),b3(3),kbz(3) 106 complex(dpc),allocatable :: ir_kibz(:,:,:,:,:) 107 complex(gwpc), ABI_CONTIGUOUS pointer :: ug_c(:),ug_v(:) 108 complex(gwpc) :: ihrc(3,Wfd%nspinor**2) 109 logical :: bbp_mask(Wfd%mband,Wfd%mband) 110 type(pawcprj_type),allocatable :: Cp_v(:,:),Cp_c(:,:) 111 112 !************************************************************************ 113 114 call wrtout(std_out," Calculating optical matrix elements in the IBZ","COLL") 115 ABI_CHECK(Wfd%nspinor==1,"nspinor==2 not coded") 116 117 comm = Wfd%comm 118 my_rank = Wfd%my_rank 119 120 nsppol = Wfd%nsppol 121 nspinor = Wfd%nspinor 122 usepaw = Wfd%usepaw 123 124 if (usepaw==1) then 125 ABI_DT_MALLOC(Cp_v,(Wfd%natom,nspinor)) 126 call pawcprj_alloc(Cp_v,0,Wfd%nlmn_atm) 127 ABI_DT_MALLOC(Cp_c,(Wfd%natom,nspinor)) 128 call pawcprj_alloc(Cp_c,0,Wfd%nlmn_atm) 129 end if 130 131 if (inclvkb==1.and.usepaw==0) then 132 MSG_ERROR("inclvkb==1 not coded,using inclvkb==2") 133 end if 134 ! 135 ! Calculate the matrix elements of ir in the IBZ. 136 ABI_MALLOC(ir_kibz,(3,lomo_min:max_band,lomo_min:max_band,Wfd%nkibz,nsppol)) 137 ir_kibz=czero 138 139 ABI_MALLOC(bbp_distrb, (Wfd%mband,Wfd%mband)) 140 141 do spin=1,nsppol 142 do ik_ibz=1,Wfd%nkibz 143 ! 144 ! Distribute the (b,b') entries. 145 bbp_mask=.FALSE.; bbp_mask(lomo_spin(spin):max_band,lomo_spin(spin):max_band)=.TRUE. 146 call wfd_distribute_bbp(Wfd,ik_ibz,spin,"All",my_nbbp,bbp_distrb,bbp_mask=bbp_mask) 147 if (ALL(bbp_distrb/=my_rank)) CYCLE 148 149 istwf_k = Wfd%istwfk(ik_ibz) 150 ABI_CHECK(istwf_k==1,"istwf_k/=1 not coded") ! KB stuff is missing. 151 npw_k = Wfd%npwarr(ik_ibz) 152 kg_k => Wfd%Kdata(ik_ibz)%kg_k 153 154 if (inclvkb/=0.and.usepaw==0) then ! Prepare term i <n,k|[Vnl,r]|n"k> 155 call vkbr_init(vkbr,Cryst,Psps,inclvkb,istwf_k,npw_k,Kmesh%ibz(:,ik_ibz),kg_k) 156 end if 157 158 ! Note: spinorial case is not coded therefore we work with ihrc(:,1). 159 ! TODO: The lower triangle can be Reconstructed by symmetry. 160 do ib_v=lomo_spin(spin),max_band ! Loop over bands 161 if ( ALL(bbp_distrb(ib_v,:)/=my_rank) ) CYCLE 162 ug_v => Wfd%Wave(ib_v,ik_ibz,spin)%ug 163 if (usepaw==1) then 164 call wfd_get_cprj(Wfd,ib_v,ik_ibz,spin,Cryst,Cp_v,sorted=.FALSE.) 165 end if 166 167 do ib_c=lomo_spin(spin),max_band 168 if (bbp_distrb(ib_v,ib_c)/=my_rank) CYCLE 169 ug_c => Wfd%Wave(ib_c,ik_ibz,spin)%ug 170 171 if (usepaw==0) then 172 ! Calculate matrix elements of i[H,r] for NC pseudopotentials. 173 ihrc = nc_ihr_comm(vkbr,cryst,psps,npw_k,nspinor,istwf_k,inclvkb,Kmesh%ibz(:,ik_ibz),ug_c,ug_v,kg_k) 174 175 else 176 ! Matrix elements of i[H,r] for PAW. 177 call wfd_get_cprj(Wfd,ib_c,ik_ibz,spin,Cryst,Cp_c,sorted=.FALSE.) 178 179 ihrc = paw_ihr(spin,nspinor,npw_k,istwf_k,Kmesh%ibz(:,ik_ibz),Cryst,Pawtab,ug_c,ug_v,kg_k,Cp_c,Cp_v,HUr) 180 end if 181 ! 182 ! Save matrix elements of i*r in the IBZ 183 ediff = KS_Bst%eig(ib_c,ik_ibz,spin) - KS_BSt%eig(ib_v,ik_ibz,spin) 184 if (ABS(ediff)<tol16) ediff=tol6 ! Treat a possible degeneracy between v and c. 185 ir_kibz(:,ib_c,ib_v,ik_ibz,spin) = ihrc(:,1)/ediff 186 187 end do !ib_c 188 end do !ib_v 189 190 call vkbr_free(vkbr) 191 end do !spin 192 end do !ik_ibz 193 194 ! Collect results on each node. 195 call xmpi_sum(ir_kibz,comm,ierr) 196 197 ABI_FREE(bbp_distrb) 198 199 if (usepaw==1) then 200 call pawcprj_free(Cp_v) 201 ABI_DT_FREE(Cp_v) 202 call pawcprj_free(Cp_c) 203 ABI_DT_FREE(Cp_c) 204 end if 205 ! 206 ! ====================================================== 207 ! ==== Calculate Fcv(kBZ) in the full Brilouin zone ==== 208 ! ====================================================== 209 ! 210 ! Symmetrization of the matrix elements of the position operator. 211 ! <Sk b|r|Sk b'> = R <k b|r|k b'> + \tau \delta_{bb'} 212 ! where S is one of the symrec operations in reciprocal space, R is the 213 ! corresponding operation in real space, \tau being the associated fractional translations. 214 ! 215 ! q.Mcv( Sk) = S^{-1}q. Mcv(k) 216 ! q.Mcv(-Sk) = -S^{-1}q. CONJG(Mcv(k)) if time-reversal is used. 217 218 b1=Cryst%gprimd(:,1)*two_pi 219 b2=Cryst%gprimd(:,2)*two_pi 220 b3=Cryst%gprimd(:,3)*two_pi 221 222 opt_cvk = czero 223 do spin=1,nsppol 224 do ik_bz=1,nkbz 225 ! 226 ! * Get ik_ibz, and symmetries index from ik_bz. 227 call get_BZ_item(Kmesh,ik_bz,kbz,ik_ibz,isym_k,itim_k) 228 229 mat_dp = DBLE(Cryst%symrec(:,:,isym_k)) 230 call matrginv(mat_dp,3,3) ! Invert 231 qrot = (3-2*itim_k) * MATMUL(mat_dp,qpoint) 232 233 do ib_v=lomo_spin(spin),max_band ! Loops over the bands C and V start 234 do ib_c=lomo_spin(spin),max_band 235 !if (ib_c==ib_v) CYCLE 236 emcvk = pdtqrc(qrot,ir_kibz(:,ib_c,ib_v,ik_ibz,spin),b1,b2,b3) 237 if (itim_k==2) emcvk = CONJG(emcvk) 238 opt_cvk(ib_c,ib_v,ik_bz,spin) = emcvk 239 end do !ib_c 240 end do !ib_v 241 242 end do !ik_bz 243 end do !spin 244 245 ABI_FREE(ir_kibz) 246 247 call xmpi_barrier(comm) 248 249 contains
ABINIT/pdtqrc [ Functions ]
NAME
pdtqrc
FUNCTION
Calculate the dot product of a real vector with a complex vector, where each is in terms of b1-b3
INPUTS
OUTPUT
PARENTS
SOURCE
267 pure function pdtqrc(R,C,b1,b2,b3) 268 269 270 !This section has been created automatically by the script Abilint (TD). 271 !Do not modify the following lines by hand. 272 #undef ABI_FUNC 273 #define ABI_FUNC 'pdtqrc' 274 !End of the abilint section 275 276 implicit none 277 278 !Arguments ------------------------------------ 279 !arrays 280 real(dp),intent(in) :: R(3),b1(3),b2(3),b3(3) 281 complex(dpc),intent(in) :: C(3) 282 complex(dpc) :: pdtqrc 283 284 !Local variables ------------------------------ 285 !scalars 286 integer :: ii 287 288 !************************************************************************ 289 290 pdtqrc=czero 291 do ii=1,3 292 pdtqrc = pdtqrc + (R(1)*b1(ii)+R(2)*b2(ii)+R(3)*b3(ii)) * & 293 & (C(1)*b1(ii)+C(2)*b2(ii)+C(3)*b3(ii)) 294 end do 295 296 end function pdtqrc