TABLE OF CONTENTS


ABINIT/m_classify_bands [ Modules ]

[ Top ] [ 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-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 .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_classify_bands
29 
30  use defs_basis
31  use defs_datatypes
32  use m_abicore
33  use m_esymm
34  use m_errors
35 
36  use m_numeric_tools,  only : get_trace
37  use m_symtk,          only : mati3inv
38  use m_hide_blas,      only : xdotc, xdotu, xcopy
39  use m_fft_mesh,       only : rotate_FFT_mesh, calc_ceigr
40  use m_crystal,        only : crystal_t
41  use m_pawang,         only : pawang_type
42  use m_pawrad,         only : pawrad_type
43  use m_pawtab,         only : pawtab_type, pawtab_get_lsize
44  use m_pawfgrtab,      only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_print, pawfgrtab_free
45  use m_pawcprj,        only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
46  use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
47  use m_paw_sphharm,    only : setsym_ylm
48  use m_paw_nhat,       only : nhatgrid
49  use m_wfd,            only : wfd_get_ur, wfd_t, wfd_ug2cprj, wfd_change_ngfft, wfd_paw_get_aeur
50 
51  implicit none
52 
53  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.

PARENTS

      sigma,wfk_analyze

CHILDREN

SOURCE

151 subroutine classify_bands(Wfd,use_paw_aeur,first_band,last_band,ik_ibz,spin,ngfftf,&
152 & Cryst,BSt,Pawtab,Pawrad,Pawang,Psps,tolsym,BSym,&
153 & EDIFF_TOL) ! optional
154 
155 
156 !This section has been created automatically by the script Abilint (TD).
157 !Do not modify the following lines by hand.
158 #undef ABI_FUNC
159 #define ABI_FUNC 'classify_bands'
160 !End of the abilint section
161 
162  implicit none
163 
164 !Arguments ------------------------------------
165 !scalars
166  integer,intent(in) :: ik_ibz,spin,first_band,last_band
167  real(dp),intent(in) :: tolsym
168  real(dp),intent(in),optional :: EDIFF_TOL
169  logical,intent(in) :: use_paw_aeur
170  type(crystal_t),intent(in) :: Cryst
171  type(pawang_type),intent(in) :: Pawang
172  type(pseudopotential_type),intent(in) :: Psps
173  type(wfd_t),intent(inout) :: Wfd
174  type(ebands_t),target,intent(in) :: BSt
175  type(esymm_t),intent(out) :: BSym
176 !arrays
177  integer,intent(in) :: ngfftf(18)
178  type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw)
179  type(Pawrad_type),intent(inout) :: Pawrad(Cryst%ntypat*Wfd%usepaw)
180 
181 !Local variables-------------------------------
182 !scalars
183  integer,parameter :: nspinor1=1
184  integer :: dim_degs,ib1,ib2,ib_stop,ib_start,iclass,idg,sym_idx
185  integer :: ir,isym,isym_class,tr_isym,jb1,jb2
186  integer :: nr1,nr2,nr3,nsym_class,nfft,cplex !,ifgd,nfgd,ifft_sph
187  integer :: ii,jj,lmax
188  integer :: optcut,optgr0,optgr1,optgr2,optrad
189  real(dp) :: EDIFF_TOL_,arg,fft_fact
190  complex(dpc) :: exp_mikg0t,exp_ikg0t,cmat_ab
191  logical :: iscompatibleFFT,found,only_trace
192  character(len=500) :: msg
193 !arrays
194  integer :: g0(3),toinv(Cryst%nsym),trial(3,3)
195  integer,pointer :: Rm1_rmt(:)
196  integer,target,allocatable :: irottb(:,:)
197  integer,allocatable :: tmp_sym(:,:,:),l_size_atm(:)
198  real(dp) :: kpt(3),kpg0(3),omat(2)
199  real(dp),pointer :: ene_k(:)
200  real(dp),pointer :: zarot(:,:,:,:)
201  complex(dpc),allocatable :: eig0r(:,:),tr_emig0r(:,:)
202  complex(gwpc),allocatable :: ur1(:),ur2(:),ur2_rot(:)
203  type(pawcprj_type),allocatable :: Cprj_b1(:,:),Cprj_b2(:,:),Cprj_b2rot(:,:)
204  type(Pawfgrtab_type),allocatable :: Pawfgrtab(:)
205  type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
206 
207 ! *************************************************************************
208 
209  DBG_ENTER("COLL")
210  !
211  ! Consistency check on input.
212  ABI_CHECK(Wfd%nspinor==1,'nspinor/=1 not coded')
213  !
214  ! By default all bands are included
215  !first_band=1; last_band=Wfd%nband(ik_ibz,spin)
216  ABI_CHECK(first_band==1,"first_band/=1 not coded")
217  ABI_CHECK(last_band<=Wfd%nband(ik_ibz,spin),"last_band cannot be > nband_k")
218 
219  EDIFF_TOL_=0.005/Ha_eV; if (PRESENT(EDIFF_TOL)) EDIFF_TOL_=ABS(EDIFF_TOL)
220 
221  call wfd_change_ngfft(Wfd,Cryst,Psps,ngfftf)
222  !
223  ! === Get index of the rotated FFT points ===
224  ! * FFT mesh in real space _must_ be compatible with symmetries.
225  nr1=Wfd%ngfft(1)
226  nr2=Wfd%ngfft(2)
227  nr3=Wfd%ngfft(3)
228  nfft=Wfd%nfft ! No FFT parallelism
229 
230  ABI_MALLOC(irottb,(nfft,Cryst%nsym))
231  call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,irottb,iscompatibleFFT)
232 
233  if (.not.iscompatibleFFT) then
234    write(msg,'(3a)')&
235 &    ' For symmetry analysis, the real space FFT mesh must be compatible with the symmetries of the space group',ch10,&
236 &    ' classify_bands will return. Action: change the input variable ngfftf '
237    MSG_WARNING(msg)
238    Bsym%err_status=1
239    Bsym%err_msg= msg
240    RETURN
241  end if
242  !
243  ! only_trace=if .TRUE. only the trace of a single matrix per class is calculated (standard procedure if
244  ! only the symmetry of bands is required). If .FALSE. all the matrices for each irreducible representation
245  ! are calculated and stored in BSym
246  only_trace=.FALSE.
247  !
248  ! ==========================================
249  ! ==== Analyse k-point symmetries first ====
250  ! ==========================================
251  ! * The analysis is done here so that we already know if there is a problem.
252  kpt=Wfd%kibz(:,ik_ibz)
253  !
254  !----Initialize the Bsym structure for this k-point and spin----!
255  ! * NOTE that all the degenerate states should be included! No check is done.
256 
257  ene_k => BSt%eig(first_band:,ik_ibz,spin) ! Select a slice of eigenvalues
258 
259  call esymm_init(Bsym,kpt,Cryst,only_trace,Wfd%nspinor,first_band,last_band,EDIFF_TOL_,ene_k,tolsym)
260  !Bsym%degs_bounds = Bsym%degs_bounds + (first_band -1)
261 
262  if (Bsym%err_status/=0) then
263    write(msg,'(a,i0,a)')" esymm_init returned err_status= ",Bsym%err_status,&
264 &    " Band classifications cannot be performed."
265    MSG_WARNING(msg)
266    RETURN
267  end if
268 
269  do ii=1,Cryst%nsym
270    call mati3inv(Cryst%symrel(:,:,ii),trial)
271    trial=transpose(trial)
272    found=.FALSE.
273    do jj=1,Cryst%nsym
274      if (ALL(trial==Cryst%symrel(:,:,jj))) then
275        toinv(ii)=jj
276        !toinv(jj)=ii
277        found=.TRUE.; EXIT
278      end if
279    end do
280    if (.not.found) then
281      MSG_ERROR("inverse not found! ")
282    end if
283  end do
284 
285  nullify(zarot)
286 
287  if (Wfd%usepaw==1) then ! Allocate cprj_k and cprj_krot to store a set of bands for a single (K,SPIN).
288    ABI_DT_MALLOC(Cprj_b1   ,(Cryst%natom,Wfd%nspinor))
289    call pawcprj_alloc(Cprj_b1,   0,Wfd%nlmn_atm)
290    ABI_DT_MALLOC(Cprj_b2   ,(Cryst%natom,Wfd%nspinor))
291    call pawcprj_alloc(Cprj_b2,   0,Wfd%nlmn_atm)
292    ABI_DT_MALLOC(Cprj_b2rot,(Cryst%natom,Wfd%nspinor))
293    call pawcprj_alloc(Cprj_b2rot,0,Wfd%nlmn_atm)
294 
295    !zarot => Pawang%zarot
296    lmax = Pawang%l_max-1
297    ABI_MALLOC(zarot,(2*lmax+1,2*lmax+1,lmax+1,Cryst%nsym))
298    zarot = Pawang%zarot
299 
300    ABI_MALLOC(tmp_sym,(3,3,Cryst%nsym))
301    do isym=1,Cryst%nsym
302      tmp_sym(:,:,isym) = Cryst%symrel(:,:,isym)
303      !tmp_sym(:,:,isym) = Cryst%symrel(:,:,toinv(isym))
304      !tmp_sym(:,:,isym) = transpose(Cryst%symrel(:,:,isym))
305      !tmp_sym(:,:,isym) = Cryst%symrec(:,:,isym)
306      !tmp_sym(:,:,isym) = TRANSPOSE(Cryst%symrec(:,:,isym))
307    end do
308    !% call setsym_ylm(Cryst%rprimd,lmax,Cryst%nsym,3,Cryst%gprimd,tmp_sym,zarot)
309    !call setsym_ylm(Cryst%gprimd,lmax,Cryst%nsym,1,Cryst%rprimd,tmp_sym,zarot)
310    ABI_FREE(tmp_sym)
311    zarot = Pawang%zarot
312 
313    cplex=1
314    call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat)
315    ABI_DT_MALLOC(Pawfgrtab,(Cryst%natom))
316    call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Wfd%nspden,Cryst%typat)
317    ABI_FREE(l_size_atm)
318 
319    optcut=1                     ! use rpaw to construct local_pawfgrtab
320    optgr0=0; optgr1=0; optgr2=0 ! dont need gY terms locally
321    optrad=1                     ! do store r-R
322 
323 
324    call nhatgrid(Cryst%atindx1,Cryst%gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,Wfd%ngfft,Cryst%ntypat,&
325 &    optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred)
326 
327    !call pawfgrtab_print(Pawfgrtab,unit=std_out,Wfd%prtvol=10)
328 
329    ABI_DT_MALLOC(Paw_onsite,(Cryst%natom))
330 
331    if (use_paw_aeur) then
332      MSG_WARNING("Using AE wavefunction for rotation in real space!")
333      call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,&
334 &                             Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab)
335    end if
336  end if
337  !
338  ! ===============================================
339  ! ==== Calculate the representation matrices ====
340  ! ===============================================
341  fft_fact=one/nfft
342  ABI_MALLOC(ur1,(nfft))
343  ABI_MALLOC(ur2,(nfft))
344  ABI_MALLOC(ur2_rot,(nfft))
345  !
346  ! * Precalculate eig0r = e^{iG0.r} on the FFT mesh.
347  ABI_MALLOC(eig0r,(nfft,Bsym%nsym_gk))
348 
349  do isym=1,Bsym%nsym_gk
350    g0=Bsym%g0(:,isym)
351    call calc_ceigr(g0,nfft,nspinor1,Wfd%ngfft,eig0r(:,isym))
352  end do
353 
354  if (Bsym%can_use_tr) then
355    ABI_MALLOC(tr_emig0r,(nfft,Bsym%nsym_trgk))
356    do isym=1,Bsym%nsym_trgk
357      g0=Bsym%tr_g0(:,isym)
358      call calc_ceigr(-g0,nfft,nspinor1,Wfd%ngfft,tr_emig0r(:,isym))
359    end do
360  end if
361  !
362  do idg=1,Bsym%ndegs ! Loop over the set of degenerate states.
363    ib_start=Bsym%degs_bounds(1,idg)
364    ib_stop =Bsym%degs_bounds(2,idg)
365    dim_degs=Bsym%degs_dim(idg)
366 
367    do ib1=ib_start,ib_stop ! First band index in the degenerate set.
368      jb1=ib1-ib_start+1
369 
370      ! debugging: use AE wave on dense FFT mesh.
371      if (Wfd%usepaw==1..and.use_paw_aeur) then
372        call wfd_paw_get_aeur(Wfd,ib1,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur1)
373      else
374        call wfd_get_ur(Wfd,ib1,ik_ibz,spin,ur1)
375        if (Wfd%usepaw==1) then
376          call wfd_ug2cprj(Wfd,ib1,ik_ibz,spin,1,0,Cryst%natom,Cryst,Cprj_b1,sorted=.FALSE.)
377        end if
378      end if
379 
380      do ib2=ib_start,ib_stop ! Second band index in the degenerate set.
381 
382        if (Bsym%only_trace.and.ib1/=ib2) CYCLE ! Only the diagonal is needed.
383 
384        if (ib2==ib1) then
385          call xcopy(nfft,ur1,1,ur2,1)
386          if (Wfd%usepaw==1) then
387            call pawcprj_copy(Cprj_b1,Cprj_b2)
388          end if
389        else
390          !
391          ! debugging: use AE wave on dense FFT mesh.
392          if (Wfd%usepaw==1.and.use_paw_aeur) then
393            call wfd_paw_get_aeur(Wfd,ib2,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur2)
394          else
395            call wfd_get_ur(Wfd,ib2,ik_ibz,spin,ur2)
396            if (Wfd%usepaw==1) then
397              call wfd_ug2cprj(Wfd,ib2,ik_ibz,spin,1,0,Cryst%natom,Cryst,Cprj_b2,sorted=.FALSE.)
398            end if
399          end if
400        end if
401        !
402        ! ===================================================
403        ! ==== Loop over the classes of the little group ====
404        ! ===================================================
405        sym_idx=0
406        do iclass=1,Bsym%nclass
407          nsym_class=Bsym%nelements(iclass)
408          !
409          do isym_class=1,nsym_class ! Loop over elements in each class.
410            sym_idx=sym_idx+1
411            if (Bsym%only_trace.and.isym_class/=1) CYCLE ! Do it once if only the character is required.
412 
413            isym=Bsym%sgk2symrec(sym_idx)
414            Rm1_rmt => irottb(:,isym)
415            !
416            ! Classify states according to the irreps of the little group of k.
417            kpg0= kpt + Bsym%g0(:,sym_idx)
418 
419            arg=-two_pi * DOT_PRODUCT(kpg0,Cryst%tnons(:,isym))
420 
421            if (ABS(arg) > tol6) then
422              exp_mikg0t=DCMPLX(DCOS(arg),DSIN(arg))
423            else
424              exp_mikg0t=cone
425            end if
426 
427            !if (Wfd%usepaw==1) then
428            !end if
429            !
430            ! === Rotate the right wave function and apply the phase ===
431            ! * Note that the k-point is the same within a lattice vector.
432            do ir=1,nfft
433              ur2_rot(ir)=ur2(Rm1_rmt(ir))*eig0r(ir,sym_idx)
434            end do
435 
436            ! * The matrix element on the FFT mesh.
437            cmat_ab = xdotc(nfft,ur1,1,ur2_rot,1)*fft_fact*exp_mikg0t
438 
439            if (Wfd%usepaw==1.and..not.use_paw_aeur) then ! Add the on-site contribution.
440              call rotate_cprj(kpt,isym,Wfd%nspinor,1,Cryst%natom,Cryst%nsym,Cryst%typat,Cryst%indsym,Cprj_b2,Cprj_b2rot)
441 
442              omat = paw_phirotphj(Wfd%nspinor,Cryst%natom,Cryst%typat,&
443 &              zarot(:,:,:,isym),Pawtab,Psps,Cprj_b1,Cprj_b2rot)
444 
445              cmat_ab = cmat_ab + DCMPLX(omat(1),omat(2)) !* exp_mikg0t
446            end if
447            !
448            jb2=ib2-ib_start+1
449            Bsym%Calc_irreps(idg)%mat(jb1,jb2,sym_idx)=cmat_ab
450 
451          end do !isym_class
452        end do !iclass
453        !
454        ! =========================================================
455        ! ==== Loop over the symmetries such that -Sk = k + G0 ====
456        ! =========================================================
457        ! <-k,a| S |k b>  = e^{i(k+G0).t} \int e^{-ig0.r} u_a u_b(R^{1}(r-t))
458        if (Bsym%can_use_tr) then
459 
460          do tr_isym=1,Bsym%nsym_trgk
461 
462            isym=Bsym%tr_sgk2symrec(tr_isym)
463            Rm1_rmt => irottb(:,isym)
464 
465            kpg0= kpt + Bsym%tr_g0(:,tr_isym)
466            arg= two_pi * DOT_PRODUCT(kpg0,Cryst%tnons(:,isym))
467 
468            if (ABS(arg) > tol6) then
469              exp_ikg0t=DCMPLX(DCOS(arg),DSIN(arg))
470            else
471              exp_ikg0t=cone
472            end if
473            !
474            ! === Rotate the right wave function and apply the phase ===
475            ! * Note that the k-point is the same within a lattice vector.
476            do ir=1,nfft
477              ur2_rot(ir)=ur2(Rm1_rmt(ir)) * tr_emig0r(ir,tr_isym)
478            end do
479            !
480            ! * The matrix element on the FFT mesh.
481            cmat_ab = xdotu(nfft,ur1,1,ur2_rot,1)*fft_fact*exp_ikg0t
482 
483            if (Wfd%usepaw==1.and..not.use_paw_aeur) then ! Add the on-site contribution. ! TODO rechek this part.
484                call rotate_cprj(kpt,isym,Wfd%nspinor,1,Cryst%natom,Cryst%nsym,Cryst%typat,Cryst%indsym,Cprj_b2,Cprj_b2rot)
485                omat = paw_phirotphj(Wfd%nspinor,Cryst%natom,Cryst%typat,&
486 &                zarot(:,:,:,isym),Pawtab,Psps,Cprj_b1,Cprj_b2rot,conjg_left=.TRUE.)
487              cmat_ab = cmat_ab + DCMPLX(omat(1),omat(2)) !* exp_ikg0t
488            end if
489            !
490            jb2=ib2-ib_start+1
491            Bsym%trCalc_irreps(idg)%mat(jb1,jb2,tr_isym)=cmat_ab
492          end do ! tr_isym
493        end if
494 
495      end do !ib2
496    end do !ib1
497    !
498    ! === Calculate the trace for each class ===
499    if (Bsym%only_trace) then ! TODO this is valid if only trace.
500      MSG_ERROR("Have to reconstruct missing traces")
501    else
502      do isym=1,Bsym%nsym_gk
503        Bsym%Calc_irreps(idg)%trace(isym) = get_trace( Bsym%Calc_irreps(idg)%mat(:,:,isym) )
504      end do
505      if (Bsym%can_use_tr) then
506        do tr_isym=1,Bsym%nsym_trgk
507          Bsym%trCalc_irreps(idg)%trace(tr_isym) = get_trace( Bsym%trCalc_irreps(idg)%mat(:,:,tr_isym) )
508        end do
509      end if
510    end if
511 
512  end do ! idg
513 
514  call esymm_finalize(Bsym,Wfd%prtvol)
515 
516  call esymm_print(Bsym,unit=std_out,prtvol=Wfd%prtvol)
517  call esymm_print(Bsym,unit=ab_out ,prtvol=Wfd%prtvol)
518  !
519  ! ===================
520  ! === Free memory ===
521  ! ===================
522  ABI_FREE(irottb)
523  ABI_FREE(ur1)
524  ABI_FREE(ur2)
525  ABI_FREE(ur2_rot)
526  ABI_FREE(eig0r)
527  if (Bsym%can_use_tr)  then
528    ABI_FREE(tr_emig0r)
529  end if
530 
531  if (Wfd%usepaw==1) then
532    call pawcprj_free(Cprj_b1)
533    ABI_DT_FREE(Cprj_b1)
534    call pawcprj_free(Cprj_b2)
535    ABI_DT_FREE(Cprj_b2)
536    call pawcprj_free(Cprj_b2rot)
537    ABI_DT_FREE(Cprj_b2rot)
538    ABI_FREE(zarot)
539    call pawfgrtab_free(Pawfgrtab)
540    ABI_DT_FREE(Pawfgrtab)
541    call paw_pwaves_lmn_free(Paw_onsite)
542    ABI_DT_FREE(Paw_onsite)
543  end if
544 
545  DBG_EXIT("COLL")
546 
547 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.

PARENTS

SOURCE

677 function paw_phirotphj(nspinor,natom,typat,zarot_isym,Pawtab,Psps,Cprj_b1,Cprj_b2,conjg_left) result(omat)
678 
679 
680 !This section has been created automatically by the script Abilint (TD).
681 !Do not modify the following lines by hand.
682 #undef ABI_FUNC
683 #define ABI_FUNC 'paw_phirotphj'
684 !End of the abilint section
685 
686  implicit none
687 
688 !Arguments ------------------------------------
689 !scalars
690  integer,intent(in) :: nspinor,natom
691  logical,optional,intent(in) :: conjg_left
692  type(pseudopotential_type),intent(in) :: Psps
693 !arrays
694  integer,intent(in) :: typat(natom)
695  real(dp),intent(in) :: zarot_isym(:,:,:)
696  real(dp) :: omat(2)
697  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
698  type(pawcprj_type),intent(in) :: Cprj_b1(natom,nspinor),Cprj_b2(natom,nspinor)
699 
700 !Local variables-------------------------------
701 !scalars
702  integer :: iat,il,ilmn,ilpm,im,itypat,jl,jlmn,jlpm,jm,k0lmn,klmn,nlmn
703  real(dp) :: dmimj,fij,im_p,re_p,sij
704  logical :: do_conjg_left
705 
706 ! *************************************************************************
707 
708  do_conjg_left = .FALSE.; if (PRESENT(conjg_left)) do_conjg_left = conjg_left
709 
710  if (nspinor/=1) then
711    MSG_ERROR("nspinor/=1 not yet coded")
712  end if
713 
714  ! === Rotate PAW projections ===
715  ! * zarot_isym is the rotation matrix of real spherical harmonics associated to symrec(:,:,isym).
716  ! * zarot_isym multiply harmonics as row vectors, we need R^{-1} but we read R and invert m,mp in the equation below
717  omat=zero
718 
719  do iat=1,natom
720    itypat=typat(iat)
721    nlmn=Pawtab(itypat)%lmn_size
722 
723    do jlmn=1,nlmn
724      k0lmn=jlmn*(jlmn-1)/2
725      jl=Psps%indlmn(1,jlmn,itypat)
726      jm=Psps%indlmn(2,jlmn,itypat)
727      jlpm=1+jl+jm
728 
729      do ilmn=1,jlmn
730        il=Psps%indlmn(1,ilmn,itypat)
731        im=Psps%indlmn(2,ilmn,itypat)
732        if (il/=jl.or.im/=jm) CYCLE ! Selection rule on l and m.
733        ilpm=1+il+im
734 
735        klmn=k0lmn+ilmn
736        sij=Pawtab(itypat)%sij(klmn) !; if (ABS(sij)<tol14) CYCLE
737 
738        ! Here we get the matrix associated to R^{-1}.
739        dmimj=zarot_isym(ilpm,jlpm,jl+1)
740 
741        if (do_conjg_left) then  ! take the complex conjugate of the left cprj.
742          re_p=  Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) &
743 &              -Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) &
744 &              +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) &
745 &              -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn)
746 
747          im_p=  Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) &
748 &              +Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) &
749 &              -Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) &
750 &              -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn)
751        else
752          re_p=  Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) &
753 &              +Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) &
754 &              +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn) &
755 &              +Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn)
756 
757          im_p=  Cprj_b1(iat,1)%cp(1,ilmn) * Cprj_b2(iat,1)%cp(2,jlmn) &
758 &              -Cprj_b1(iat,1)%cp(2,ilmn) * Cprj_b2(iat,1)%cp(1,jlmn) &
759 &              +Cprj_b1(iat,1)%cp(1,jlmn) * Cprj_b2(iat,1)%cp(2,ilmn) &
760 &              -Cprj_b1(iat,1)%cp(2,jlmn) * Cprj_b2(iat,1)%cp(1,ilmn)
761        end if
762        !
763        ! * Accumulate the atom-centered contributions.
764        fij = Pawtab(itypat)%dltij(klmn)/two
765        omat(1)= omat(1) + fij*sij*re_p*dmimj
766        omat(2)= omat(2) + fij*sij*im_p*dmimj
767 
768      end do !ilmn
769    end do !jlmn
770  end do !iat
771 
772 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.

PARENTS

      classify_bands

CHILDREN

SOURCE

586 subroutine rotate_cprj(kpoint,isym,nspinor,nbnds,natom,nsym,typat,indsym,Cprj_in,Cprj_out)
587 
588 
589 !This section has been created automatically by the script Abilint (TD).
590 !Do not modify the following lines by hand.
591 #undef ABI_FUNC
592 #define ABI_FUNC 'rotate_cprj'
593 !End of the abilint section
594 
595  implicit none
596 
597 !Arguments ------------------------------------
598 !scalars
599  integer,intent(in) :: nbnds,nspinor,natom,isym,nsym
600 !arrays
601  integer,intent(in) :: typat(natom),indsym(4,nsym,natom)
602  real(dp),intent(in) :: kpoint(3)
603  type(pawcprj_type),intent(in) :: Cprj_in(natom,nspinor*nbnds)
604  type(pawcprj_type),intent(out) :: Cprj_out(natom,nspinor*nbnds)
605 
606 !Local variables-------------------------------
607 !scalars
608  integer :: iat,iband,itypat,iat_sym
609  real(dp) :: kdotr0
610 !arrays
611  integer :: r0(3)
612  real(dp) :: phase_kr0(2)
613 
614 ! *************************************************************************
615 
616  !call pawcprj_copy(cprj_in,cprj_out)
617  !RETURN
618 
619  do iat=1,natom
620    itypat=typat(iat)
621    !
622    ! The index of the symmetric atom.
623    ! * R^{-1} (xred(:,iat)-tnons) = xred(:,iat_sym) + r0.
624    ! * phase_kr0 takes into account the case in which rotated atom is in another unit cell.
625    iat_sym=indsym(4,isym,iat); r0=indsym(1:3,isym,iat)
626 
627    kdotr0 = two_pi*DOT_PRODUCT(kpoint,r0)
628    phase_kr0(1) = DCOS(kdotr0)
629    phase_kr0(2) = DSIN(kdotr0)
630 
631    !phase_kr0 = (/one,zero/)
632 
633    do iband=1,nspinor*nbnds
634      Cprj_out(iat,iband)%cp(1,:)=  Cprj_in(iat_sym,iband)%cp(1,:)*phase_kr0(1) &
635 &                                 -Cprj_in(iat_sym,iband)%cp(2,:)*phase_kr0(2)
636 
637      Cprj_out(iat,iband)%cp(2,:)=  Cprj_in(iat_sym,iband)%cp(1,:)*phase_kr0(2) &
638 &                                 +Cprj_in(iat_sym,iband)%cp(2,:)*phase_kr0(1)
639    end do
640  end do ! iat
641 
642 end subroutine rotate_cprj