TABLE OF CONTENTS
- ABINIT/m_classify_bands
- m_classify_bands/classify_bands
- m_classify_bands/paw_phirotphj
- m_classify_bands/rotate_cprj
ABINIT/m_classify_bands [ Modules ]
NAME
m_classify_bands
FUNCTION
Finds the irreducible representation associated to a set of degenerate bands at a given k-point and spin.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_classify_bands 24 25 use defs_basis 26 use m_abicore 27 use m_esymm 28 use m_errors 29 30 use defs_datatypes, only : pseudopotential_type, ebands_t 31 use m_numeric_tools, only : get_trace 32 use m_symtk, only : mati3inv 33 use m_hide_blas, only : xdotc, xdotu, xcopy 34 use m_fft_mesh, only : rotate_FFT_mesh, calc_ceigr 35 use m_crystal, only : crystal_t 36 use m_pawang, only : pawang_type 37 use m_pawrad, only : pawrad_type 38 use m_pawtab, only : pawtab_type, pawtab_get_lsize 39 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_print, pawfgrtab_free 40 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy 41 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 42 use m_paw_sphharm, only : setsym_ylm 43 use m_paw_nhat, only : nhatgrid 44 use m_wfd, only : wfd_t 45 46 implicit none 47 48 private
m_classify_bands/classify_bands [ Functions ]
[ Top ] [ m_classify_bands ] [ Functions ]
NAME
classify_bands
FUNCTION
This routine finds the irreducible representation associated to a set of degenerate bands at a given k-point and spin. The irreducible representation is obtained by rotating the set of degenerate wavefunctions using the symmetry operations in the little group of k. Two states are treated as degenerate if their energy differs by less than EDIFF_TOL.
INPUTS
Wfd(wfd_t)= structure gathering information on wave functions %nibz=Number of points in the IBZ. %nsppol=number of independent spin polarizations %usepaw=1 if PAW ik_ibz=The index of the k-point in the IBZ. spin=The spin index. ngfft(18)=Info on the FFT mesh to be used for evaluting u(r) and the rotated u(R^{1}(r-t)). ngfft must be compatible with the symmetries of the crystal and can differ from Wfd%ngfft. wfd_change_ngfft is called if ANY(Wfd%ngfft(1:3) =/ ngfft). Cryst<crystal_t>=Type gathering info on the crystal structure. %nsym=Number of operations in space group %ntypat=Number of type of atoms (onlu for PAW) %symrec(3,3,nsym)=Symmetry operations in reciprocal space (reduced coordinates) %tnons(3,nsym)=Fractional translations %typat(natom)=Type of each atom BSt<ebands_t>=Datatype with electronic energies. Pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data. Pawang <type(pawang_type)>=paw angular mesh and related data Psps<pseudopotential_type> %indlmn(6,lmnmax,ntypat)=array giving l,m,n,lm,ln,spin for i=lmn (for each atom type) Dtfil<datafiles_type>=variables related to files %unpaw tolsym=Tolerance for the symmetries (input variable) [EDIFF_TOL]= tolerance on the energy difference of two states (if not specified is set to 0.005 eV)
OUTPUT
BSym<Bands_Symmetries>=structure containing info on the little group of the k-point as well
as the character of the representation associated to each set of degenerate states
if BSym%isymmorphic the symmetry analysis cannot be performed, usually it means that
k is at zone border and there are non-symmorphic translations (see Notes)
NOTES
* Let M(R_t) the irreducible representation associated to the space group symmetry (R_t). * By convention M(R_t) multiplies wave functions as a row vector: $ R_t \psi_a(r) = \psi_a (R^{-1}(r-\tau)) = \sum_b M(R_t)_{ba} \psi_b $ Therefore, if R_t belongs to the little group of k (i.e. Sk=k+G0), one obtains: $ M_ab(R_t) = e^{-i(k+G0).\tau} \int e^{iG0.r} u_{ak}(r)^* u_{bk}(R^{-1}(r-\tau)) \,dr $. * The irreducible representation of the small _point_ group of k, M_ab(R), suffices to classify the degenerate eigenstates provided that particular conditions are fulfilled (see limitations below). The matrix is indeed given by: $ M_ab(R) = e^{+ik.\tau} M_ab(R_t) = e^{-iG0.\tau} \int e^{iG0.r} u_{ak}(r)^* u_{bk}(R^{-1}(r-\tau))\,dr $ The phase factor outside the integral should be zero since symmetry analysis at border zone in non-symmorphic space groups is not available. Anyway it is included in our expressions for the sake of consistency. * For PAW there is an additional onsite terms involving <phi_i|phi_j(R^{-1}(r-\tau)> and the pseudized version that can be evaluated using the rotation matrix for real spherical harmonis, zarot(mp,m,l,R). $ Y_{lm}(Rr)= \sum_{m'} zarot(m',m,ll,R) Y_{lm'}(r) $ $ M^{onsite}_ab(R_t) = sum_{c ij} <\tpsi_a| p_i^c> <p_j^{c'}|\tpsi_b\> \times [ <\phi_i^c|\phi_j^{c'}> - <\tphi_i^c|\tphi_j^{c'}> ]. $ $ [ <\phi_i^c|\phi_j^{c'}> - <\tphi_i^c|\tphi_j^{c'}> ] = s_{ij} D_{\mi\mj}^\lj(R^{-1}) $ where c' is the rotated atom i.e c' = R^{-1}( c-\tau) and D is the rotation matrix for real spherical harmonics. Remember that zarot(m',m,l,R)=zarot(m,m',l,R^{-1}) and $ Y^l_m(ISG) = sum_{m'} D_{m'm}(S) Y_{m'}^l(G) (-i)^l $ $ D_{m'm}^l (R) = D_{m,m'}^l (R^{-1}) $ * LIMITATIONS: The method does not work if k is at zone border and the little group of k contains a non-symmorphic fractional translation.
SOURCE
141 subroutine classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,& 142 & Cryst,BSt,Pawtab,Pawrad,Pawang,Psps,tolsym,BSym,& 143 & EDIFF_TOL) ! optional 144 145 !Arguments ------------------------------------ 146 !scalars 147 integer,intent(in) :: ik_ibz,spin,first_band,last_band 148 real(dp),intent(in) :: tolsym 149 real(dp),intent(in),optional :: EDIFF_TOL 150 logical,intent(in) :: use_paw_aeur 151 type(crystal_t),intent(in) :: Cryst 152 type(pawang_type),intent(in) :: Pawang 153 type(pseudopotential_type),intent(in) :: Psps 154 class(wfd_t),intent(inout) :: Wfd 155 type(ebands_t),target,intent(in) :: BSt 156 type(esymm_t),intent(out) :: BSym 157 !arrays 158 integer,intent(in) :: ngfftf(18) 159 type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 160 type(Pawrad_type),intent(inout) :: Pawrad(Cryst%ntypat*Wfd%usepaw) 161 162 !Local variables------------------------------- 163 !scalars 164 integer,parameter :: nspinor1=1 165 integer :: dim_degs,ib1,ib2,ib_stop,ib_start,iclass,idg,sym_idx 166 integer :: ir,isym,isym_class,tr_isym,jb1,jb2 167 integer :: nr1,nr2,nr3,nsym_class,nfft,cplex !,ifgd,nfgd,ifft_sph 168 integer :: ii,jj,lmax 169 integer :: optcut,optgr0,optgr1,optgr2,optrad 170 real(dp) :: EDIFF_TOL_,arg,fft_fact 171 complex(dpc) :: exp_mikg0t,exp_ikg0t,cmat_ab 172 logical :: iscompatibleFFT,found,only_trace 173 character(len=500) :: msg 174 !arrays 175 integer :: g0(3),toinv(Cryst%nsym),trial(3,3) 176 integer,pointer :: Rm1_rmt(:) 177 integer,target,allocatable :: irottb(:,:) 178 integer,allocatable :: tmp_sym(:,:,:),l_size_atm(:) 179 real(dp) :: kpt(3),kpg0(3),omat(2) 180 real(dp),pointer :: ene_k(:) 181 real(dp),pointer :: zarot(:,:,:,:) 182 complex(dpc),allocatable :: eig0r(:,:),tr_emig0r(:,:) 183 complex(gwpc),allocatable :: ur1(:),ur2(:),ur2_rot(:) 184 type(pawcprj_type),allocatable :: Cprj_b1(:,:),Cprj_b2(:,:),Cprj_b2rot(:,:) 185 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 186 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 187 188 ! ************************************************************************* 189 190 DBG_ENTER("COLL") 191 ! 192 ! Consistency check on input. 193 ABI_CHECK(Wfd%nspinor==1,'nspinor/=1 not coded') 194 ! 195 ! By default all bands are included 196 !first_band=1; last_band=Wfd%nband(ik_ibz,spin) 197 ABI_CHECK(first_band==1,"first_band/=1 not coded") 198 ABI_CHECK(last_band<=Wfd%nband(ik_ibz,spin),"last_band cannot be > nband_k") 199 200 EDIFF_TOL_=0.005/Ha_eV; if (PRESENT(EDIFF_TOL)) EDIFF_TOL_=ABS(EDIFF_TOL) 201 202 call wfd%change_ngfft(Cryst,Psps,ngfftf) 203 ! 204 ! === Get index of the rotated FFT points === 205 ! * FFT mesh in real space _must_ be compatible with symmetries. 206 nr1=Wfd%ngfft(1) 207 nr2=Wfd%ngfft(2) 208 nr3=Wfd%ngfft(3) 209 nfft=Wfd%nfft ! No FFT parallelism 210 211 ABI_MALLOC(irottb,(nfft,Cryst%nsym)) 212 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,irottb,iscompatibleFFT) 213 214 if (.not.iscompatibleFFT) then 215 write(msg,'(3a)')& 216 & ' For symmetry analysis, the real space FFT mesh must be compatible with the symmetries of the space group',ch10,& 217 & ' classify_bands will return. Action: change the input variable ngfftf ' 218 ABI_WARNING(msg) 219 Bsym%err_status=1 220 Bsym%err_msg= msg 221 RETURN 222 end if 223 ! 224 ! only_trace=if .TRUE. only the trace of a single matrix per class is calculated (standard procedure if 225 ! only the symmetry of bands is required). If .FALSE. all the matrices for each irreducible representation 226 ! are calculated and stored in BSym 227 only_trace=.FALSE. 228 ! 229 ! ========================================== 230 ! ==== Analyse k-point symmetries first ==== 231 ! ========================================== 232 ! * The analysis is done here so that we already know if there is a problem. 233 kpt=Wfd%kibz(:,ik_ibz) 234 ! 235 !----Initialize the Bsym structure for this k-point and spin----! 236 ! * NOTE that all the degenerate states should be included! No check is done. 237 238 ene_k => BSt%eig(first_band:,ik_ibz,spin) ! Select a slice of eigenvalues 239 240 call esymm_init(Bsym,kpt,Cryst,only_trace,Wfd%nspinor,first_band,last_band,EDIFF_TOL_,ene_k,tolsym) 241 !Bsym%degs_bounds = Bsym%degs_bounds + (first_band -1) 242 243 if (Bsym%err_status/=0) then 244 write(msg,'(a,i0,a)')" esymm_init returned err_status= ",Bsym%err_status,& 245 & " Band classifications cannot be performed." 246 ABI_WARNING(msg) 247 RETURN 248 end if 249 250 do ii=1,Cryst%nsym 251 call mati3inv(Cryst%symrel(:,:,ii),trial) 252 trial=transpose(trial) 253 found=.FALSE. 254 do jj=1,Cryst%nsym 255 if (ALL(trial==Cryst%symrel(:,:,jj))) then 256 toinv(ii)=jj 257 !toinv(jj)=ii 258 found=.TRUE.; EXIT 259 end if 260 end do 261 if (.not.found) then 262 ABI_ERROR("inverse not found! ") 263 end if 264 end do 265 266 nullify(zarot) 267 268 if (Wfd%usepaw==1) then ! Allocate cprj_k and cprj_krot to store a set of bands for a single (K,SPIN). 269 ABI_MALLOC(Cprj_b1 ,(Cryst%natom,Wfd%nspinor)) 270 call pawcprj_alloc(Cprj_b1, 0,Wfd%nlmn_atm) 271 ABI_MALLOC(Cprj_b2 ,(Cryst%natom,Wfd%nspinor)) 272 call pawcprj_alloc(Cprj_b2, 0,Wfd%nlmn_atm) 273 ABI_MALLOC(Cprj_b2rot,(Cryst%natom,Wfd%nspinor)) 274 call pawcprj_alloc(Cprj_b2rot,0,Wfd%nlmn_atm) 275 276 !zarot => Pawang%zarot 277 lmax = Pawang%l_max-1 278 ABI_MALLOC(zarot,(2*lmax+1,2*lmax+1,lmax+1,Cryst%nsym)) 279 zarot = Pawang%zarot 280 281 ABI_MALLOC(tmp_sym,(3,3,Cryst%nsym)) 282 do isym=1,Cryst%nsym 283 tmp_sym(:,:,isym) = Cryst%symrel(:,:,isym) 284 !tmp_sym(:,:,isym) = Cryst%symrel(:,:,toinv(isym)) 285 !tmp_sym(:,:,isym) = transpose(Cryst%symrel(:,:,isym)) 286 !tmp_sym(:,:,isym) = Cryst%symrec(:,:,isym) 287 !tmp_sym(:,:,isym) = TRANSPOSE(Cryst%symrec(:,:,isym)) 288 end do 289 !% call setsym_ylm(Cryst%rprimd,lmax,Cryst%nsym,3,Cryst%gprimd,tmp_sym,zarot) 290 !call setsym_ylm(Cryst%gprimd,lmax,Cryst%nsym,1,Cryst%rprimd,tmp_sym,zarot) 291 ABI_FREE(tmp_sym) 292 zarot = Pawang%zarot 293 294 cplex=1 295 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 296 ABI_MALLOC(Pawfgrtab,(Cryst%natom)) 297 call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Wfd%nspden,Cryst%typat) 298 ABI_FREE(l_size_atm) 299 300 optcut=1 ! use rpaw to construct local_pawfgrtab 301 optgr0=0; optgr1=0; optgr2=0 ! dont need gY terms locally 302 optrad=1 ! do store r-R 303 304 305 call nhatgrid(Cryst%atindx1,Cryst%gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,Wfd%ngfft,Cryst%ntypat,& 306 & optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 307 308 !call pawfgrtab_print(Pawfgrtab,unit=std_out,Wfd%prtvol=10) 309 310 ABI_MALLOC(Paw_onsite,(Cryst%natom)) 311 312 if (use_paw_aeur) then 313 ABI_WARNING("Using AE wavefunction for rotation in real space!") 314 call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,& 315 & Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab) 316 end if 317 end if 318 ! 319 ! =============================================== 320 ! ==== Calculate the representation matrices ==== 321 ! =============================================== 322 fft_fact=one/nfft 323 ABI_MALLOC(ur1,(nfft)) 324 ABI_MALLOC(ur2,(nfft)) 325 ABI_MALLOC(ur2_rot,(nfft)) 326 ! 327 ! * Precalculate eig0r = e^{iG0.r} on the FFT mesh. 328 ABI_MALLOC(eig0r,(nfft,Bsym%nsym_gk)) 329 330 do isym=1,Bsym%nsym_gk 331 g0=Bsym%g0(:,isym) 332 call calc_ceigr(g0,nfft,nspinor1,Wfd%ngfft,eig0r(:,isym)) 333 end do 334 335 if (Bsym%can_use_tr) then 336 ABI_MALLOC(tr_emig0r,(nfft,Bsym%nsym_trgk)) 337 do isym=1,Bsym%nsym_trgk 338 g0=Bsym%tr_g0(:,isym) 339 call calc_ceigr(-g0,nfft,nspinor1,Wfd%ngfft,tr_emig0r(:,isym)) 340 end do 341 end if 342 ! 343 do idg=1,Bsym%ndegs ! Loop over the set of degenerate states. 344 ib_start=Bsym%degs_bounds(1,idg) 345 ib_stop =Bsym%degs_bounds(2,idg) 346 dim_degs=Bsym%degs_dim(idg) 347 348 do ib1=ib_start,ib_stop ! First band index in the degenerate set. 349 jb1=ib1-ib_start+1 350 351 ! debugging: use AE wave on dense FFT mesh. 352 if (Wfd%usepaw==1..and.use_paw_aeur) then 353 call wfd%paw_get_aeur(ib1,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur1) 354 else 355 call wfd%get_ur(ib1,ik_ibz,spin,ur1) 356 if (Wfd%usepaw==1) then 357 call wfd%ug2cprj(ib1,ik_ibz,spin,1,0,Cryst%natom,Cryst,Cprj_b1,sorted=.FALSE.) 358 end if 359 end if 360 361 do ib2=ib_start,ib_stop ! Second band index in the degenerate set. 362 363 if (Bsym%only_trace.and.ib1/=ib2) CYCLE ! Only the diagonal is needed. 364 365 if (ib2==ib1) then 366 call xcopy(nfft,ur1,1,ur2,1) 367 if (Wfd%usepaw==1) then 368 call pawcprj_copy(Cprj_b1,Cprj_b2) 369 end if 370 else 371 ! 372 ! debugging: use AE wave on dense FFT mesh. 373 if (Wfd%usepaw==1.and.use_paw_aeur) then 374 call wfd%paw_get_aeur(ib2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur2) 375 else 376 call wfd%get_ur(ib2,ik_ibz,spin,ur2) 377 if (Wfd%usepaw==1) then 378 call wfd%ug2cprj(ib2,ik_ibz,spin,1,0,Cryst%natom,Cryst,Cprj_b2,sorted=.FALSE.) 379 end if 380 end if 381 end if 382 ! 383 ! =================================================== 384 ! ==== Loop over the classes of the little group ==== 385 ! =================================================== 386 sym_idx=0 387 do iclass=1,Bsym%nclass 388 nsym_class=Bsym%nelements(iclass) 389 ! 390 do isym_class=1,nsym_class ! Loop over elements in each class. 391 sym_idx=sym_idx+1 392 if (Bsym%only_trace.and.isym_class/=1) CYCLE ! Do it once if only the character is required. 393 394 isym=Bsym%sgk2symrec(sym_idx) 395 Rm1_rmt => irottb(:,isym) 396 ! 397 ! Classify states according to the irreps of the little group of k. 398 kpg0= kpt + Bsym%g0(:,sym_idx) 399 400 arg=-two_pi * DOT_PRODUCT(kpg0,Cryst%tnons(:,isym)) 401 402 if (ABS(arg) > tol6) then 403 exp_mikg0t=DCMPLX(DCOS(arg),DSIN(arg)) 404 else 405 exp_mikg0t=cone 406 end if 407 408 !if (Wfd%usepaw==1) then 409 !end if 410 ! 411 ! === Rotate the right wave function and apply the phase === 412 ! * Note that the k-point is the same within a lattice vector. 413 do ir=1,nfft 414 ur2_rot(ir)=ur2(Rm1_rmt(ir))*eig0r(ir,sym_idx) 415 end do 416 417 ! * The matrix element on the FFT mesh. 418 cmat_ab = xdotc(nfft,ur1,1,ur2_rot,1)*fft_fact*exp_mikg0t 419 420 if (Wfd%usepaw==1.and..not.use_paw_aeur) then ! Add the on-site contribution. 421 call rotate_cprj(kpt,isym,Wfd%nspinor,1,Cryst%natom,Cryst%nsym,Cryst%typat,Cryst%indsym,Cprj_b2,Cprj_b2rot) 422 423 omat = paw_phirotphj(Wfd%nspinor,Cryst%natom,Cryst%typat,& 424 & zarot(:,:,:,isym),Pawtab,Psps,Cprj_b1,Cprj_b2rot) 425 426 cmat_ab = cmat_ab + DCMPLX(omat(1),omat(2)) !* exp_mikg0t 427 end if 428 ! 429 jb2=ib2-ib_start+1 430 Bsym%Calc_irreps(idg)%mat(jb1,jb2,sym_idx)=cmat_ab 431 432 end do !isym_class 433 end do !iclass 434 ! 435 ! ========================================================= 436 ! ==== Loop over the symmetries such that -Sk = k + G0 ==== 437 ! ========================================================= 438 ! <-k,a| S |k b> = e^{i(k+G0).t} \int e^{-ig0.r} u_a u_b(R^{1}(r-t)) 439 if (Bsym%can_use_tr) then 440 441 do tr_isym=1,Bsym%nsym_trgk 442 443 isym=Bsym%tr_sgk2symrec(tr_isym) 444 Rm1_rmt => irottb(:,isym) 445 446 kpg0= kpt + Bsym%tr_g0(:,tr_isym) 447 arg= two_pi * DOT_PRODUCT(kpg0,Cryst%tnons(:,isym)) 448 449 if (ABS(arg) > tol6) then 450 exp_ikg0t=DCMPLX(DCOS(arg),DSIN(arg)) 451 else 452 exp_ikg0t=cone 453 end if 454 ! 455 ! === Rotate the right wave function and apply the phase === 456 ! * Note that the k-point is the same within a lattice vector. 457 do ir=1,nfft 458 ur2_rot(ir)=ur2(Rm1_rmt(ir)) * tr_emig0r(ir,tr_isym) 459 end do 460 ! 461 ! * The matrix element on the FFT mesh. 462 cmat_ab = xdotu(nfft,ur1,1,ur2_rot,1)*fft_fact*exp_ikg0t 463 464 if (Wfd%usepaw==1.and..not.use_paw_aeur) then ! Add the on-site contribution. ! TODO rechek this part. 465 call rotate_cprj(kpt,isym,Wfd%nspinor,1,Cryst%natom,Cryst%nsym,Cryst%typat,Cryst%indsym,Cprj_b2,Cprj_b2rot) 466 omat = paw_phirotphj(Wfd%nspinor,Cryst%natom,Cryst%typat,& 467 & zarot(:,:,:,isym),Pawtab,Psps,Cprj_b1,Cprj_b2rot,conjg_left=.TRUE.) 468 cmat_ab = cmat_ab + DCMPLX(omat(1),omat(2)) !* exp_ikg0t 469 end if 470 ! 471 jb2=ib2-ib_start+1 472 Bsym%trCalc_irreps(idg)%mat(jb1,jb2,tr_isym)=cmat_ab 473 end do ! tr_isym 474 end if 475 476 end do !ib2 477 end do !ib1 478 ! 479 ! === Calculate the trace for each class === 480 if (Bsym%only_trace) then ! TODO this is valid if only trace. 481 ABI_ERROR("Have to reconstruct missing traces") 482 else 483 do isym=1,Bsym%nsym_gk 484 Bsym%Calc_irreps(idg)%trace(isym) = get_trace( Bsym%Calc_irreps(idg)%mat(:,:,isym) ) 485 end do 486 if (Bsym%can_use_tr) then 487 do tr_isym=1,Bsym%nsym_trgk 488 Bsym%trCalc_irreps(idg)%trace(tr_isym) = get_trace( Bsym%trCalc_irreps(idg)%mat(:,:,tr_isym) ) 489 end do 490 end if 491 end if 492 493 end do ! idg 494 495 call esymm_finalize(Bsym,Wfd%prtvol) 496 497 call esymm_print(Bsym,unit=std_out,prtvol=Wfd%prtvol) 498 call esymm_print(Bsym,unit=ab_out ,prtvol=Wfd%prtvol) 499 ! 500 ! =================== 501 ! === Free memory === 502 ! =================== 503 ABI_FREE(irottb) 504 ABI_FREE(ur1) 505 ABI_FREE(ur2) 506 ABI_FREE(ur2_rot) 507 ABI_FREE(eig0r) 508 if (Bsym%can_use_tr) then 509 ABI_FREE(tr_emig0r) 510 end if 511 512 if (Wfd%usepaw==1) then 513 call pawcprj_free(Cprj_b1) 514 ABI_FREE(Cprj_b1) 515 call pawcprj_free(Cprj_b2) 516 ABI_FREE(Cprj_b2) 517 call pawcprj_free(Cprj_b2rot) 518 ABI_FREE(Cprj_b2rot) 519 ABI_FREE(zarot) 520 call pawfgrtab_free(Pawfgrtab) 521 ABI_FREE(Pawfgrtab) 522 call paw_pwaves_lmn_free(Paw_onsite) 523 ABI_FREE(Paw_onsite) 524 end if 525 526 DBG_EXIT("COLL") 527 528 end subroutine classify_bands
m_classify_bands/paw_phirotphj [ Functions ]
[ Top ] [ m_classify_bands ] [ Functions ]
NAME
paw_phirotphj
FUNCTION
This routine calculates
<\tPsi_1|\tprj_i> <\tprj_j|\tPsi_2> [ <\phi_i|\phi_j(R^{-1}r> - <\tphi_i|\tphi_j(R^{-1}r> ]
[ <\phi_i|\phi_j(R^{-1}r> - <\tphi_i|\tphi_j(R^{-1}r> ] = s_ij D_{mi,mi}^{li}(R)
INPUTS
nspinor=Number of spinorial components. natom=number of atoms typat(natom)=type of eahc atom zarot_isym Pawtab(ntypat)<Pawtab_type>=PAW tabulated starting data Psps<pseudopotential_type>=Info on pseudopotentials. Cprj_b1(natom,nspinor)<type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with all NL projectors at fixed k-point Cprj_b2(natom,nspinor)<type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with all NL projectors at fixed k-point [conjg_left]=.TRUE if the complex conjugate of the left wavefunctions has to be taken. Defaults to .FALSE.
OUTPUT
omat(2)=The onsite matrix element.
SOURCE
642 function paw_phirotphj(nspinor,natom,typat,zarot_isym,Pawtab,Psps,Cprj_b1,Cprj_b2,conjg_left) result(omat) 643 644 !Arguments ------------------------------------ 645 !scalars 646 integer,intent(in) :: nspinor,natom 647 logical,optional,intent(in) :: conjg_left 648 type(pseudopotential_type),intent(in) :: Psps 649 !arrays 650 integer,intent(in) :: typat(natom) 651 real(dp),intent(in) :: zarot_isym(:,:,:) 652 real(dp) :: omat(2) 653 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat) 654 type(pawcprj_type),intent(in) :: Cprj_b1(natom,nspinor),Cprj_b2(natom,nspinor) 655 656 !Local variables------------------------------- 657 !scalars 658 integer :: iat,il,ilmn,ilpm,im,itypat,jl,jlmn,jlpm,jm,k0lmn,klmn,nlmn 659 real(dp) :: dmimj,fij,im_p,re_p,sij 660 logical :: do_conjg_left 661 662 ! ************************************************************************* 663 664 do_conjg_left = .FALSE.; if (PRESENT(conjg_left)) do_conjg_left = conjg_left 665 666 if (nspinor/=1) then 667 ABI_ERROR("nspinor/=1 not yet coded") 668 end if 669 670 ! === Rotate PAW projections === 671 ! * zarot_isym is the rotation matrix of real spherical harmonics associated to symrec(:,:,isym). 672 ! * zarot_isym multiply harmonics as row vectors, we need R^{-1} but we read R and invert m,mp in the equation below 673 omat=zero 674 675 do iat=1,natom 676 itypat=typat(iat) 677 nlmn=Pawtab(itypat)%lmn_size 678 679 do jlmn=1,nlmn 680 k0lmn=jlmn*(jlmn-1)/2 681 jl=Psps%indlmn(1,jlmn,itypat) 682 jm=Psps%indlmn(2,jlmn,itypat) 683 jlpm=1+jl+jm 684 685 do ilmn=1,jlmn 686 il=Psps%indlmn(1,ilmn,itypat) 687 im=Psps%indlmn(2,ilmn,itypat) 688 if (il/=jl.or.im/=jm) CYCLE ! Selection rule on l and m. 689 ilpm=1+il+im 690 691 klmn=k0lmn+ilmn 692 sij=Pawtab(itypat)%sij(klmn) !; if (ABS(sij)<tol14) CYCLE 693 694 ! Here we get the matrix associated to R^{-1}. 695 dmimj=zarot_isym(ilpm,jlpm,jl+1) 696 697 if (do_conjg_left) then ! take the complex conjugate of the left cprj. 698 re_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) & 699 & -Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) & 700 & +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) & 701 & -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) 702 703 im_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) & 704 & +Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) & 705 & -Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) & 706 & -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) 707 else 708 re_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) & 709 & +Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) & 710 & +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) & 711 & +Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) 712 713 im_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) & 714 & -Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) & 715 & +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) & 716 & -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) 717 end if 718 ! 719 ! * Accumulate the atom-centered contributions. 720 fij = Pawtab(itypat)%dltij(klmn)/two 721 omat(1)= omat(1) + fij*sij*re_p*dmimj 722 omat(2)= omat(2) + fij*sij*im_p*dmimj 723 724 end do !ilmn 725 end do !jlmn 726 end do !iat 727 728 end function paw_phirotphj
m_classify_bands/rotate_cprj [ Functions ]
[ Top ] [ m_classify_bands ] [ Functions ]
NAME
rotate_cprj
FUNCTION
Rotate cprj matrix elements by applying the symmetry operation of index isym that preserves the given k-point within a reciprocal lattice vector.
INPUTS
isym=index of the symmetry in the symrec arrays that preserves the given k-point within a reciprocal lattice vector ntypat=number of types of atom. natom=number of atoms. Cryst<crystal_t>=Datatype gathering info on the unit cell. typat(natom)=type of each atom. nbnds=number of bands for this k-point ans spin Cprj_in(natom,nbnds)<type(pawcprj_type)>= projected input wave functions <Proj_i|Cnk> with all NL projectors at fixed k-point
OUTPUT
Cprj_out(natom,nbnds) <type(pawcprj_type)>= projection of the smooth PAW wave function onto projectors centered on equivalent sites of the crystal (non restricted to be in the firs unit cell) The equivalent site is defined according to the symmetry operation isym. Thus Cprj_out contains Cprj_out(at,b)=<p_j^{R^{-1}(L_{at}-\tau)} | \tpsi_b> if R is the isym operation with fractional translation \tau L_{at} is the position of the initial atom inside the first unit cell Note that atom a might be in a cell different from the initial one. No wrapping is done.
SOURCE
562 subroutine rotate_cprj(kpoint,isym,nspinor,nbnds,natom,nsym,typat,indsym,Cprj_in,Cprj_out) 563 564 !Arguments ------------------------------------ 565 !scalars 566 integer,intent(in) :: nbnds,nspinor,natom,isym,nsym 567 !arrays 568 integer,intent(in) :: typat(natom),indsym(4,nsym,natom) 569 real(dp),intent(in) :: kpoint(3) 570 type(pawcprj_type),intent(in) :: Cprj_in(natom,nspinor*nbnds) 571 type(pawcprj_type),intent(out) :: Cprj_out(natom,nspinor*nbnds) 572 573 !Local variables------------------------------- 574 !scalars 575 integer :: iat,iband,itypat,iat_sym 576 real(dp) :: kdotr0 577 !arrays 578 integer :: r0(3) 579 real(dp) :: phase_kr0(2) 580 581 ! ************************************************************************* 582 583 !call pawcprj_copy(cprj_in,cprj_out) 584 !RETURN 585 586 do iat=1,natom 587 itypat=typat(iat) 588 ! 589 ! The index of the symmetric atom. 590 ! * R^{-1} (xred(:,iat)-tnons) = xred(:,iat_sym) + r0. 591 ! * phase_kr0 takes into account the case in which rotated atom is in another unit cell. 592 iat_sym=indsym(4,isym,iat); r0=indsym(1:3,isym,iat) 593 594 kdotr0 = two_pi*DOT_PRODUCT(kpoint,r0) 595 phase_kr0(1) = DCOS(kdotr0) 596 phase_kr0(2) = DSIN(kdotr0) 597 598 !phase_kr0 = (/one,zero/) 599 600 do iband=1,nspinor*nbnds 601 Cprj_out(iat,iband)%cp(1,:)= Cprj_in(iat_sym,iband)%cp(1,:)*phase_kr0(1) & 602 & -Cprj_in(iat_sym,iband)%cp(2,:)*phase_kr0(2) 603 604 Cprj_out(iat,iband)%cp(2,:)= Cprj_in(iat_sym,iband)%cp(1,:)*phase_kr0(2) & 605 & +Cprj_in(iat_sym,iband)%cp(2,:)*phase_kr0(1) 606 end do 607 end do ! iat 608 609 end subroutine rotate_cprj