TABLE OF CONTENTS
ABINIT/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.
COPYRIGHT
Copyright (C) 2008-2018 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 .
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.
PARENTS
sigma,wfk_analyze
CHILDREN
SOURCE
96 #if defined HAVE_CONFIG_H 97 #include "config.h" 98 #endif 99 100 #include "abi_common.h" 101 102 subroutine classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,& 103 & Cryst,BSt,Pawtab,Pawrad,Pawang,Psps,tolsym,BSym,& 104 & EDIFF_TOL) ! optional 105 106 use defs_basis 107 use defs_datatypes 108 use m_profiling_abi 109 use m_esymm 110 use m_errors 111 112 use m_numeric_tools, only : get_trace 113 use m_blas, only : xdotc, xdotu, xcopy 114 use m_fft_mesh, only : rotate_FFT_mesh, calc_ceigr 115 use m_crystal, only : crystal_t 116 use m_pawang, only : pawang_type 117 use m_pawrad, only : pawrad_type 118 use m_pawtab, only : pawtab_type,pawtab_get_lsize 119 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_print, pawfgrtab_free 120 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy 121 use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free 122 use m_wfd, only : wfd_get_ur, wfd_t, wfd_ug2cprj, wfd_change_ngfft, wfd_paw_get_aeur 123 124 !This section has been created automatically by the script Abilint (TD). 125 !Do not modify the following lines by hand. 126 #undef ABI_FUNC 127 #define ABI_FUNC 'classify_bands' 128 use interfaces_32_util 129 use interfaces_65_paw 130 use interfaces_69_wfdesc, except_this_one => classify_bands 131 !End of the abilint section 132 133 implicit none 134 135 !Arguments ------------------------------------ 136 !scalars 137 integer,intent(in) :: ik_ibz,spin,first_band,last_band 138 real(dp),intent(in) :: tolsym 139 real(dp),intent(in),optional :: EDIFF_TOL 140 logical,intent(in) :: use_paw_aeur 141 type(crystal_t),intent(in) :: Cryst 142 type(pawang_type),intent(in) :: Pawang 143 type(pseudopotential_type),intent(in) :: Psps 144 type(wfd_t),intent(inout) :: Wfd 145 type(ebands_t),target,intent(in) :: BSt 146 type(esymm_t),intent(out) :: BSym 147 !arrays 148 integer,intent(in) :: ngfftf(18) 149 type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 150 type(Pawrad_type),intent(inout) :: Pawrad(Cryst%ntypat*Wfd%usepaw) 151 152 !Local variables------------------------------- 153 !scalars 154 integer,parameter :: nspinor1=1 155 integer :: dim_degs,ib1,ib2,ib_stop,ib_start,iclass,idg,sym_idx 156 integer :: ir,isym,isym_class,tr_isym,jb1,jb2 157 integer :: nr1,nr2,nr3,nsym_class,nfft,cplex !,ifgd,nfgd,ifft_sph 158 integer :: ii,jj,lmax 159 integer :: optcut,optgr0,optgr1,optgr2,optrad 160 real(dp) :: EDIFF_TOL_,arg,fft_fact 161 complex(dpc) :: exp_mikg0t,exp_ikg0t,cmat_ab 162 logical :: iscompatibleFFT,found,only_trace 163 character(len=500) :: msg 164 !arrays 165 integer :: g0(3),toinv(Cryst%nsym),trial(3,3) 166 integer,pointer :: Rm1_rmt(:) 167 integer,target,allocatable :: irottb(:,:) 168 integer,allocatable :: tmp_sym(:,:,:),l_size_atm(:) 169 real(dp) :: kpt(3),kpg0(3),omat(2) 170 real(dp),pointer :: ene_k(:) 171 real(dp),pointer :: zarot(:,:,:,:) 172 complex(dpc),allocatable :: eig0r(:,:),tr_emig0r(:,:) 173 complex(gwpc),allocatable :: ur1(:),ur2(:),ur2_rot(:) 174 type(pawcprj_type),allocatable :: Cprj_b1(:,:),Cprj_b2(:,:),Cprj_b2rot(:,:) 175 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 176 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 177 178 ! ************************************************************************* 179 180 DBG_ENTER("COLL") 181 ! 182 ! Consistency check on input. 183 ABI_CHECK(Wfd%nspinor==1,'nspinor/=1 not coded') 184 ! 185 ! By default all bands are included 186 !first_band=1; last_band=Wfd%nband(ik_ibz,spin) 187 ABI_CHECK(first_band==1,"first_band/=1 not coded") 188 ABI_CHECK(last_band<=Wfd%nband(ik_ibz,spin),"last_band cannot be > nband_k") 189 190 EDIFF_TOL_=0.005/Ha_eV; if (PRESENT(EDIFF_TOL)) EDIFF_TOL_=ABS(EDIFF_TOL) 191 192 call wfd_change_ngfft(Wfd,Cryst,Psps,ngfftf) 193 ! 194 ! === Get index of the rotated FFT points === 195 ! * FFT mesh in real space _must_ be compatible with symmetries. 196 nr1=Wfd%ngfft(1) 197 nr2=Wfd%ngfft(2) 198 nr3=Wfd%ngfft(3) 199 nfft=Wfd%nfft ! No FFT parallelism 200 201 ABI_MALLOC(irottb,(nfft,Cryst%nsym)) 202 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,irottb,iscompatibleFFT) 203 204 if (.not.iscompatibleFFT) then 205 write(msg,'(3a)')& 206 & ' For symmetry analysis, the real space FFT mesh must be compatible with the symmetries of the space group',ch10,& 207 & ' classify_bands will return. Action: change the input variable ngfftf ' 208 MSG_WARNING(msg) 209 Bsym%err_status=1 210 Bsym%err_msg= msg 211 RETURN 212 end if 213 ! 214 ! only_trace=if .TRUE. only the trace of a single matrix per class is calculated (standard procedure if 215 ! only the symmetry of bands is required). If .FALSE. all the matrices for each irreducible representation 216 ! are calculated and stored in BSym 217 only_trace=.FALSE. 218 ! 219 ! ========================================== 220 ! ==== Analyse k-point symmetries first ==== 221 ! ========================================== 222 ! * The analysis is done here so that we already know if there is a problem. 223 kpt=Wfd%kibz(:,ik_ibz) 224 ! 225 !----Initialize the Bsym structure for this k-point and spin----! 226 ! * NOTE that all the degenerate states should be included! No check is done. 227 228 ene_k => BSt%eig(first_band:,ik_ibz,spin) ! Select a slice of eigenvalues 229 230 call esymm_init(Bsym,kpt,Cryst,only_trace,Wfd%nspinor,first_band,last_band,EDIFF_TOL_,ene_k,tolsym) 231 !Bsym%degs_bounds = Bsym%degs_bounds + (first_band -1) 232 233 if (Bsym%err_status/=0) then 234 write(msg,'(a,i0,a)')" esymm_init returned err_status= ",Bsym%err_status,& 235 & " Band classifications cannot be performed." 236 MSG_WARNING(msg) 237 RETURN 238 end if 239 240 do ii=1,Cryst%nsym 241 call mati3inv(Cryst%symrel(:,:,ii),trial) 242 trial=transpose(trial) 243 found=.FALSE. 244 do jj=1,Cryst%nsym 245 if (ALL(trial==Cryst%symrel(:,:,jj))) then 246 toinv(ii)=jj 247 !toinv(jj)=ii 248 found=.TRUE.; EXIT 249 end if 250 end do 251 if (.not.found) then 252 MSG_ERROR("inverse not found! ") 253 end if 254 end do 255 256 nullify(zarot) 257 258 if (Wfd%usepaw==1) then ! Allocate cprj_k and cprj_krot to store a set of bands for a single (K,SPIN). 259 ABI_DT_MALLOC(Cprj_b1 ,(Cryst%natom,Wfd%nspinor)) 260 call pawcprj_alloc(Cprj_b1, 0,Wfd%nlmn_atm) 261 ABI_DT_MALLOC(Cprj_b2 ,(Cryst%natom,Wfd%nspinor)) 262 call pawcprj_alloc(Cprj_b2, 0,Wfd%nlmn_atm) 263 ABI_DT_MALLOC(Cprj_b2rot,(Cryst%natom,Wfd%nspinor)) 264 call pawcprj_alloc(Cprj_b2rot,0,Wfd%nlmn_atm) 265 266 !zarot => Pawang%zarot 267 lmax = Pawang%l_max-1 268 ABI_MALLOC(zarot,(2*lmax+1,2*lmax+1,lmax+1,Cryst%nsym)) 269 zarot = Pawang%zarot 270 271 ABI_MALLOC(tmp_sym,(3,3,Cryst%nsym)) 272 do isym=1,Cryst%nsym 273 tmp_sym(:,:,isym) = Cryst%symrel(:,:,isym) 274 !tmp_sym(:,:,isym) = Cryst%symrel(:,:,toinv(isym)) 275 !tmp_sym(:,:,isym) = transpose(Cryst%symrel(:,:,isym)) 276 !tmp_sym(:,:,isym) = Cryst%symrec(:,:,isym) 277 !tmp_sym(:,:,isym) = TRANSPOSE(Cryst%symrec(:,:,isym)) 278 end do 279 !% call setsymrhoij(Cryst%rprimd,lmax,Cryst%nsym,3,Cryst%gprimd,tmp_sym,zarot) 280 !call setsymrhoij(Cryst%gprimd,lmax,Cryst%nsym,1,Cryst%rprimd,tmp_sym,zarot) 281 ABI_FREE(tmp_sym) 282 zarot = Pawang%zarot 283 284 cplex=1 285 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 286 ABI_DT_MALLOC(Pawfgrtab,(Cryst%natom)) 287 call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Wfd%nspden,Cryst%typat) 288 ABI_FREE(l_size_atm) 289 290 optcut=1 ! use rpaw to construct local_pawfgrtab 291 optgr0=0; optgr1=0; optgr2=0 ! dont need gY terms locally 292 optrad=1 ! do store r-R 293 294 295 call nhatgrid(Cryst%atindx1,Cryst%gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,Wfd%ngfft,Cryst%ntypat,& 296 & optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 297 298 !call pawfgrtab_print(Pawfgrtab,unit=std_out,Wfd%prtvol=10) 299 300 ABI_DT_MALLOC(Paw_onsite,(Cryst%natom)) 301 302 if (use_paw_aeur) then 303 MSG_WARNING("Using AE wavefunction for rotation in real space!") 304 call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,& 305 & Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab) 306 end if 307 end if 308 ! 309 ! =============================================== 310 ! ==== Calculate the representation matrices ==== 311 ! =============================================== 312 fft_fact=one/nfft 313 ABI_MALLOC(ur1,(nfft)) 314 ABI_MALLOC(ur2,(nfft)) 315 ABI_MALLOC(ur2_rot,(nfft)) 316 ! 317 ! * Precalculate eig0r = e^{iG0.r} on the FFT mesh. 318 ABI_MALLOC(eig0r,(nfft,Bsym%nsym_gk)) 319 320 do isym=1,Bsym%nsym_gk 321 g0=Bsym%g0(:,isym) 322 call calc_ceigr(g0,nfft,nspinor1,Wfd%ngfft,eig0r(:,isym)) 323 end do 324 325 if (Bsym%can_use_tr) then 326 ABI_MALLOC(tr_emig0r,(nfft,Bsym%nsym_trgk)) 327 do isym=1,Bsym%nsym_trgk 328 g0=Bsym%tr_g0(:,isym) 329 call calc_ceigr(-g0,nfft,nspinor1,Wfd%ngfft,tr_emig0r(:,isym)) 330 end do 331 end if 332 ! 333 do idg=1,Bsym%ndegs ! Loop over the set of degenerate states. 334 ib_start=Bsym%degs_bounds(1,idg) 335 ib_stop =Bsym%degs_bounds(2,idg) 336 dim_degs=Bsym%degs_dim(idg) 337 338 do ib1=ib_start,ib_stop ! First band index in the degenerate set. 339 jb1=ib1-ib_start+1 340 341 ! debugging: use AE wave on dense FFT mesh. 342 if (Wfd%usepaw==1..and.use_paw_aeur) then 343 call wfd_paw_get_aeur(Wfd,ib1,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur1) 344 else 345 call wfd_get_ur(Wfd,ib1,ik_ibz,spin,ur1) 346 if (Wfd%usepaw==1) then 347 call wfd_ug2cprj(Wfd,ib1,ik_ibz,spin,1,0,Cryst%natom,Cryst,Cprj_b1,sorted=.FALSE.) 348 end if 349 end if 350 351 do ib2=ib_start,ib_stop ! Second band index in the degenerate set. 352 353 if (Bsym%only_trace.and.ib1/=ib2) CYCLE ! Only the diagonal is needed. 354 355 if (ib2==ib1) then 356 call xcopy(nfft,ur1,1,ur2,1) 357 if (Wfd%usepaw==1) then 358 call pawcprj_copy(Cprj_b1,Cprj_b2) 359 end if 360 else 361 ! 362 ! debugging: use AE wave on dense FFT mesh. 363 if (Wfd%usepaw==1.and.use_paw_aeur) then 364 call wfd_paw_get_aeur(Wfd,ib2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur2) 365 else 366 call wfd_get_ur(Wfd,ib2,ik_ibz,spin,ur2) 367 if (Wfd%usepaw==1) then 368 call wfd_ug2cprj(Wfd,ib2,ik_ibz,spin,1,0,Cryst%natom,Cryst,Cprj_b2,sorted=.FALSE.) 369 end if 370 end if 371 end if 372 ! 373 ! =================================================== 374 ! ==== Loop over the classes of the little group ==== 375 ! =================================================== 376 sym_idx=0 377 do iclass=1,Bsym%nclass 378 nsym_class=Bsym%nelements(iclass) 379 ! 380 do isym_class=1,nsym_class ! Loop over elements in each class. 381 sym_idx=sym_idx+1 382 if (Bsym%only_trace.and.isym_class/=1) CYCLE ! Do it once if only the character is required. 383 384 isym=Bsym%sgk2symrec(sym_idx) 385 Rm1_rmt => irottb(:,isym) 386 ! 387 ! Classify states according to the irreps of the little group of k. 388 kpg0= kpt + Bsym%g0(:,sym_idx) 389 390 arg=-two_pi * DOT_PRODUCT(kpg0,Cryst%tnons(:,isym)) 391 392 if (ABS(arg) > tol6) then 393 exp_mikg0t=DCMPLX(DCOS(arg),DSIN(arg)) 394 else 395 exp_mikg0t=cone 396 end if 397 398 !if (Wfd%usepaw==1) then 399 !end if 400 ! 401 ! === Rotate the right wave function and apply the phase === 402 ! * Note that the k-point is the same within a lattice vector. 403 do ir=1,nfft 404 ur2_rot(ir)=ur2(Rm1_rmt(ir))*eig0r(ir,sym_idx) 405 end do 406 407 ! * The matrix element on the FFT mesh. 408 cmat_ab = xdotc(nfft,ur1,1,ur2_rot,1)*fft_fact*exp_mikg0t 409 410 if (Wfd%usepaw==1.and..not.use_paw_aeur) then ! Add the on-site contribution. 411 call rotate_cprj(kpt,isym,Wfd%nspinor,1,Cryst%natom,Cryst%nsym,Cryst%typat,Cryst%indsym,Cprj_b2,Cprj_b2rot) 412 413 omat = paw_phirotphj(Wfd%nspinor,Cryst%natom,Cryst%typat,& 414 & zarot(:,:,:,isym),Pawtab,Psps,Cprj_b1,Cprj_b2rot) 415 416 cmat_ab = cmat_ab + DCMPLX(omat(1),omat(2)) !* exp_mikg0t 417 end if 418 ! 419 jb2=ib2-ib_start+1 420 Bsym%Calc_irreps(idg)%mat(jb1,jb2,sym_idx)=cmat_ab 421 422 end do !isym_class 423 end do !iclass 424 ! 425 ! ========================================================= 426 ! ==== Loop over the symmetries such that -Sk = k + G0 ==== 427 ! ========================================================= 428 ! <-k,a| S |k b> = e^{i(k+G0).t} \int e^{-ig0.r} u_a u_b(R^{1}(r-t)) 429 if (Bsym%can_use_tr) then 430 431 do tr_isym=1,Bsym%nsym_trgk 432 433 isym=Bsym%tr_sgk2symrec(tr_isym) 434 Rm1_rmt => irottb(:,isym) 435 436 kpg0= kpt + Bsym%tr_g0(:,tr_isym) 437 arg= two_pi * DOT_PRODUCT(kpg0,Cryst%tnons(:,isym)) 438 439 if (ABS(arg) > tol6) then 440 exp_ikg0t=DCMPLX(DCOS(arg),DSIN(arg)) 441 else 442 exp_ikg0t=cone 443 end if 444 ! 445 ! === Rotate the right wave function and apply the phase === 446 ! * Note that the k-point is the same within a lattice vector. 447 do ir=1,nfft 448 ur2_rot(ir)=ur2(Rm1_rmt(ir)) * tr_emig0r(ir,tr_isym) 449 end do 450 ! 451 ! * The matrix element on the FFT mesh. 452 cmat_ab = xdotu(nfft,ur1,1,ur2_rot,1)*fft_fact*exp_ikg0t 453 454 if (Wfd%usepaw==1.and..not.use_paw_aeur) then ! Add the on-site contribution. ! TODO rechek this part. 455 call rotate_cprj(kpt,isym,Wfd%nspinor,1,Cryst%natom,Cryst%nsym,Cryst%typat,Cryst%indsym,Cprj_b2,Cprj_b2rot) 456 omat = paw_phirotphj(Wfd%nspinor,Cryst%natom,Cryst%typat,& 457 & zarot(:,:,:,isym),Pawtab,Psps,Cprj_b1,Cprj_b2rot,conjg_left=.TRUE.) 458 cmat_ab = cmat_ab + DCMPLX(omat(1),omat(2)) !* exp_ikg0t 459 end if 460 ! 461 jb2=ib2-ib_start+1 462 Bsym%trCalc_irreps(idg)%mat(jb1,jb2,tr_isym)=cmat_ab 463 end do ! tr_isym 464 end if 465 466 end do !ib2 467 end do !ib1 468 ! 469 ! === Calculate the trace for each class === 470 if (Bsym%only_trace) then ! TODO this is valid if only trace. 471 MSG_ERROR("Have to reconstruct missing traces") 472 else 473 do isym=1,Bsym%nsym_gk 474 Bsym%Calc_irreps(idg)%trace(isym) = get_trace( Bsym%Calc_irreps(idg)%mat(:,:,isym) ) 475 end do 476 if (Bsym%can_use_tr) then 477 do tr_isym=1,Bsym%nsym_trgk 478 Bsym%trCalc_irreps(idg)%trace(tr_isym) = get_trace( Bsym%trCalc_irreps(idg)%mat(:,:,tr_isym) ) 479 end do 480 end if 481 end if 482 483 end do ! idg 484 485 call esymm_finalize(Bsym,Wfd%prtvol) 486 487 call esymm_print(Bsym,unit=std_out,prtvol=Wfd%prtvol) 488 call esymm_print(Bsym,unit=ab_out ,prtvol=Wfd%prtvol) 489 ! 490 ! =================== 491 ! === Free memory === 492 ! =================== 493 ABI_FREE(irottb) 494 ABI_FREE(ur1) 495 ABI_FREE(ur2) 496 ABI_FREE(ur2_rot) 497 ABI_FREE(eig0r) 498 if (Bsym%can_use_tr) then 499 ABI_FREE(tr_emig0r) 500 end if 501 502 if (Wfd%usepaw==1) then 503 call pawcprj_free(Cprj_b1) 504 ABI_DT_FREE(Cprj_b1) 505 call pawcprj_free(Cprj_b2) 506 ABI_DT_FREE(Cprj_b2) 507 call pawcprj_free(Cprj_b2rot) 508 ABI_DT_FREE(Cprj_b2rot) 509 ABI_FREE(zarot) 510 call pawfgrtab_free(Pawfgrtab) 511 ABI_DT_FREE(Pawfgrtab) 512 call paw_pwaves_lmn_free(Paw_onsite) 513 ABI_DT_FREE(Paw_onsite) 514 end if 515 516 DBG_EXIT("COLL") 517 518 end subroutine classify_bands
ABINIT/paw_phirotphj [ 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.
PARENTS
SOURCE
653 function paw_phirotphj(nspinor,natom,typat,zarot_isym,Pawtab,Psps,Cprj_b1,Cprj_b2,conjg_left) result(omat) 654 655 use defs_basis 656 use defs_datatypes 657 use m_errors 658 659 use m_pawtab, only : pawtab_type 660 use m_pawcprj, only : pawcprj_type 661 662 !This section has been created automatically by the script Abilint (TD). 663 !Do not modify the following lines by hand. 664 #undef ABI_FUNC 665 #define ABI_FUNC 'paw_phirotphj' 666 !End of the abilint section 667 668 implicit none 669 670 !Arguments ------------------------------------ 671 !scalars 672 integer,intent(in) :: nspinor,natom 673 logical,optional,intent(in) :: conjg_left 674 type(pseudopotential_type),intent(in) :: Psps 675 !arrays 676 integer,intent(in) :: typat(natom) 677 real(dp),intent(in) :: zarot_isym(:,:,:) 678 real(dp) :: omat(2) 679 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat) 680 type(pawcprj_type),intent(in) :: Cprj_b1(natom,nspinor),Cprj_b2(natom,nspinor) 681 682 !Local variables------------------------------- 683 !scalars 684 integer :: iat,il,ilmn,ilpm,im,itypat,jl,jlmn,jlpm,jm,k0lmn,klmn,nlmn 685 real(dp) :: dmimj,fij,im_p,re_p,sij 686 logical :: do_conjg_left 687 688 ! ************************************************************************* 689 690 do_conjg_left = .FALSE.; if (PRESENT(conjg_left)) do_conjg_left = conjg_left 691 692 if (nspinor/=1) then 693 MSG_ERROR("nspinor/=1 not yet coded") 694 end if 695 696 ! === Rotate PAW projections === 697 ! * zarot_isym is the rotation matrix of real spherical harmonics associated to symrec(:,:,isym). 698 ! * zarot_isym multiply harmonics as row vectors, we need R^{-1} but we read R and invert m,mp in the equation below 699 omat=zero 700 701 do iat=1,natom 702 itypat=typat(iat) 703 nlmn=Pawtab(itypat)%lmn_size 704 705 do jlmn=1,nlmn 706 k0lmn=jlmn*(jlmn-1)/2 707 jl=Psps%indlmn(1,jlmn,itypat) 708 jm=Psps%indlmn(2,jlmn,itypat) 709 jlpm=1+jl+jm 710 711 do ilmn=1,jlmn 712 il=Psps%indlmn(1,ilmn,itypat) 713 im=Psps%indlmn(2,ilmn,itypat) 714 if (il/=jl.or.im/=jm) CYCLE ! Selection rule on l and m. 715 ilpm=1+il+im 716 717 klmn=k0lmn+ilmn 718 sij=Pawtab(itypat)%sij(klmn) !; if (ABS(sij)<tol14) CYCLE 719 720 ! Here we get the matrix associated to R^{-1}. 721 dmimj=zarot_isym(ilpm,jlpm,jl+1) 722 723 if (do_conjg_left) then ! take the complex conjugate of the left cprj. 724 re_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) & 725 & -Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) & 726 & +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) & 727 & -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) 728 729 im_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) & 730 & +Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) & 731 & -Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) & 732 & -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) 733 else 734 re_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) & 735 & +Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) & 736 & +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) & 737 & +Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) 738 739 im_p= Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) & 740 & -Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) & 741 & +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) & 742 & -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) 743 end if 744 ! 745 ! * Accumulate the atom-centered contributions. 746 fij = Pawtab(itypat)%dltij(klmn)/two 747 omat(1)= omat(1) + fij*sij*re_p*dmimj 748 omat(2)= omat(2) + fij*sij*im_p*dmimj 749 750 end do !ilmn 751 end do !jlmn 752 end do !iat 753 754 end function paw_phirotphj
ABINIT/rotate_cprj [ 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.
PARENTS
classify_bands
CHILDREN
SOURCE
557 subroutine rotate_cprj(kpoint,isym,nspinor,nbnds,natom,nsym,typat,indsym,Cprj_in,Cprj_out) 558 559 use defs_basis 560 use m_profiling_abi 561 use m_errors 562 563 use m_pawcprj, only : pawcprj_type, pawcprj_copy 564 565 !This section has been created automatically by the script Abilint (TD). 566 !Do not modify the following lines by hand. 567 #undef ABI_FUNC 568 #define ABI_FUNC 'rotate_cprj' 569 !End of the abilint section 570 571 implicit none 572 573 !Arguments ------------------------------------ 574 !scalars 575 integer,intent(in) :: nbnds,nspinor,natom,isym,nsym 576 !arrays 577 integer,intent(in) :: typat(natom),indsym(4,nsym,natom) 578 real(dp),intent(in) :: kpoint(3) 579 type(pawcprj_type),intent(in) :: Cprj_in(natom,nspinor*nbnds) 580 type(pawcprj_type),intent(out) :: Cprj_out(natom,nspinor*nbnds) 581 582 !Local variables------------------------------- 583 !scalars 584 integer :: iat,iband,itypat,iat_sym 585 real(dp) :: kdotr0 586 !arrays 587 integer :: r0(3) 588 real(dp) :: phase_kr0(2) 589 590 ! ************************************************************************* 591 592 !call pawcprj_copy(cprj_in,cprj_out) 593 !RETURN 594 595 do iat=1,natom 596 itypat=typat(iat) 597 ! 598 ! The index of the symmetric atom. 599 ! * R^{-1} (xred(:,iat)-tnons) = xred(:,iat_sym) + r0. 600 ! * phase_kr0 takes into account the case in which rotated atom is in another unit cell. 601 iat_sym=indsym(4,isym,iat); r0=indsym(1:3,isym,iat) 602 603 kdotr0 = two_pi*DOT_PRODUCT(kpoint,r0) 604 phase_kr0(1) = DCOS(kdotr0) 605 phase_kr0(2) = DSIN(kdotr0) 606 607 !phase_kr0 = (/one,zero/) 608 609 do iband=1,nspinor*nbnds 610 Cprj_out(iat,iband)%cp(1,:)= Cprj_in(iat_sym,iband)%cp(1,:)*phase_kr0(1) & 611 & -Cprj_in(iat_sym,iband)%cp(2,:)*phase_kr0(2) 612 613 Cprj_out(iat,iband)%cp(2,:)= Cprj_in(iat_sym,iband)%cp(1,:)*phase_kr0(2) & 614 & +Cprj_in(iat_sym,iband)%cp(2,:)*phase_kr0(1) 615 end do 616 end do ! iat 617 618 end subroutine rotate_cprj