symm_hf.F90 22.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   This module generates the little group of k and the extended irr. !
!   BZ. Furthermore it calculates the irr. representation             !
!                                                                     !
!   P(R,T)\phi_n,k = \sum_{n'} rep_v(n',n) *\phi_n',k        !
!   where                                                             !
!         P  is an element of the little group of k                   !
!         n' runs over group of degenerat states belonging to n.      !
!                                                                     !
!                                             M.Betzinger (09/07)     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      MODULE m_symm_hf

#define irreps .false.

      CONTAINS

Daniel Wortmann's avatar
Daniel Wortmann committed
18
      SUBROUTINE symm_hf(kpts,nk,sym,&
Daniel Wortmann's avatar
Daniel Wortmann committed
19 20
     &                   dimension,hybdat,eig_irr,&
     &                   atoms,hybrid,cell,&
21 22 23 24 25 26 27 28 29 30
     &                   lapw,jsp,&
     &                   gpt,&
     &                   mpi,irank2,&
     &                   nsymop,psym,nkpt_EIBZ,n_q,parent,&
     &                   symop,degenerat,pointer_EIBZ,maxndb,nddb,&
     &                   nsest,indx_sest,rep_c ) 


      USE m_constants
      USE m_util   ,ONLY: modulo1,intgrf,intgrf_init
Daniel Wortmann's avatar
Daniel Wortmann committed
31
      USE m_olap   ,ONLY: wfolap_inv,wfolap_noinv,wfolap1,wfolap_init
32
      USE m_trafo  ,ONLY: waveftrafo_symm
Daniel Wortmann's avatar
Daniel Wortmann committed
33
      USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
34
      USE m_io_hybrid
35 36
      IMPLICIT NONE

Daniel Wortmann's avatar
Daniel Wortmann committed
37 38
      TYPE(t_hybdat),INTENT(IN)   :: hybdat

39 40
      TYPE(t_mpi),INTENT(IN)   :: mpi
      TYPE(t_dimension),INTENT(IN)   :: dimension
Daniel Wortmann's avatar
Daniel Wortmann committed
41 42
      TYPE(t_hybrid),INTENT(IN) :: hybrid
      TYPE(t_sym),INTENT(IN)    :: sym
43 44
      TYPE(t_cell),INTENT(IN)   :: cell
      TYPE(t_kpts),INTENT(IN)   :: kpts
Daniel Wortmann's avatar
Daniel Wortmann committed
45
      TYPE(t_atoms),INTENT(IN)  :: atoms
46 47 48
      TYPE(t_lapw),INTENT(IN)   :: lapw

!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
49
      INTEGER,INTENT(IN)              :: nk  
50 51 52 53 54 55 56 57
      INTEGER,INTENT(IN)              :: jsp
      INTEGER,INTENT(IN)              :: irank2
      INTEGER,INTENT(OUT)             :: nkpt_EIBZ
      INTEGER,INTENT(OUT)             :: nsymop
      INTEGER,INTENT(OUT)             :: maxndb,nddb

!     - arrays -
      INTEGER,INTENT(IN)              :: gpt(3,lapw%nv(jsp))
Daniel Wortmann's avatar
Daniel Wortmann committed
58 59
      INTEGER,INTENT(OUT)             :: parent(kpts%nkptf)
      INTEGER,INTENT(OUT)             :: symop(kpts%nkptf)
Daniel Wortmann's avatar
Daniel Wortmann committed
60 61
      INTEGER,INTENT(INOUT)           :: degenerat(hybdat%ne_eig(nk))
      INTEGER,INTENT(OUT)             :: nsest(hybdat%nbands(nk)), indx_sest(hybdat%nbands(nk),hybdat%nbands(nk))  
62 63 64
      INTEGER,ALLOCATABLE,INTENT(OUT) :: pointer_EIBZ(:)
      INTEGER,ALLOCATABLE,INTENT(OUT) :: psym(:),n_q(:)      

Daniel Wortmann's avatar
Daniel Wortmann committed
65
      REAL,INTENT(IN)                 :: eig_irr(dimension%neigd,kpts%nkpt)
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86
      COMPLEX,ALLOCATABLE,INTENT(OUT) :: rep_c(:,:,:,:,:)

!     - local scalars -
      INTEGER                         :: ikpt,ikpt1,iop,isym,iisym,m
      INTEGER                         :: itype,ieq,iatom,ratom
      INTEGER                         :: iband,iband1,iband2,iatom0
      INTEGER                         :: i,j,ic,ic1,ic2
      INTEGER                         :: irecl_cmt,irecl_z
      INTEGER                         :: ok
      INTEGER                         :: l,lm
      INTEGER                         :: n1,n2,nn
      INTEGER                         :: ndb,ndb1,ndb2
      INTEGER                         :: nrkpt

      REAL                            :: tolerance,pi

      COMPLEX                         :: cdum
      COMPLEX , PARAMETER             :: img = (0d0,1d0)

!     - local arrays -
      INTEGER                         :: rrot(3,3,sym%nsym)
Daniel Wortmann's avatar
Daniel Wortmann committed
87 88
      INTEGER                         :: neqvkpt(kpts%nkptf)
      INTEGER                         :: list(kpts%nkptf)
89 90 91 92 93 94 95 96 97 98 99 100
      INTEGER,ALLOCATABLE             :: help(:)
      
      REAL                            :: rotkpt(3),g(3)
      REAL   ,ALLOCATABLE             :: olapmt(:,:,:,:)

      COMPLEX                         :: cmt(dimension%neigd,hybrid%maxlmindx,atoms%nat)
      COMPLEX                         :: carr1(hybrid%maxlmindx,atoms%nat)
      COMPLEX,ALLOCATABLE             :: carr(:),wavefolap(:,:)
      COMPLEX,ALLOCATABLE             :: cmthlp(:,:,:)
      COMPLEX,ALLOCATABLE             :: cpwhlp(:,:)
      COMPLEX,ALLOCATABLE             :: trace(:,:)

Daniel Wortmann's avatar
Daniel Wortmann committed
101
      TYPE(t_mat)                     :: olappw,z
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
      COMPLEX,ALLOCATABLE             :: rep_d(:,:,:)
      LOGICAL,ALLOCATABLE             :: symequivalent(:,:)

      IF ( irank2 == 0 ) THEN
        WRITE(6,'(A)') new_line('n') // new_line('n') // '### subroutine: symm ###'
      END IF

      ! calculate rotations in reciprocal space
      DO i = 1,sym%nsym
        IF( i .le. sym%nop ) THEN
          rrot(:,:,i) = transpose(sym%mrot(:,:,sym%invtab(i)))
        ELSE
          rrot(:,:,i) = -rrot(:,:,i-sym%nop)
        END IF
      END DO


      ! determine little group of k., i.e. those symmetry operations
      ! which keep bk(:,nk) invariant
      ! nsymop :: number of such symmetry-operations
      ! psym   :: points to the symmetry-operation

      ic = 0
      ALLOCATE(psym(sym%nsym))
      
      DO iop=1,sym%nsym

Daniel Wortmann's avatar
Daniel Wortmann committed
129
        rotkpt = matmul( rrot(:,:,iop), kpts%bkf(:,nk) )
130 131 132 133 134

        !transfer rotkpt into BZ
        rotkpt = modulo1(rotkpt,kpts%nkpt3)

        !check if rotkpt is identical to bk(:,nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
135
        IF( maxval( abs( rotkpt - kpts%bkf(:,nk) ) ) .le. 1E-07) THEN
136 137 138 139 140 141 142 143 144
          ic = ic + 1
          psym(ic) = iop
        END IF
      END DO
      nsymop = ic
     

      IF ( irank2 == 0 ) THEN
        WRITE(6,'(A,i3)') ' nk',nk
Daniel Wortmann's avatar
Daniel Wortmann committed
145
        WRITE(6,'(A,3f10.5)') ' kpts%bkf(:,nk):',kpts%bkf(:,nk)
146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
        WRITE(6,'(A,i3)') ' Number of elements in the little group:',nsymop
      END IF

      ! reallocate psym
      ALLOCATE(help(ic))
      help = psym(1:ic)
      DEALLOCATE(psym)
      ALLOCATE(psym(ic))
      psym =help
      DEALLOCATE(help)

      ! determine extented irreducible BZ of k ( EIBZ(k) ), i.e.
      ! those k-points, which can generate the whole BZ by
      ! applying the symmetry operations of the little group of k

      neqvkpt = 0

Daniel Wortmann's avatar
Daniel Wortmann committed
163
      DO i=1,kpts%nkptf
164 165 166 167
        list(i) = i-1
      END DO

      symop = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
168
      DO ikpt=2,kpts%nkptf
169 170
        DO iop=1,nsymop

Daniel Wortmann's avatar
Daniel Wortmann committed
171
          rotkpt = matmul( rrot(:,:,psym(iop)), kpts%bkf(:,ikpt) )
172 173 174 175 176 177

          !transfer rotkpt into BZ
          rotkpt = modulo1(rotkpt,kpts%nkpt3)

          !determine number of rotkpt
          nrkpt = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
178 179
          DO ikpt1=1,kpts%nkptf
            IF ( maxval( abs( rotkpt - kpts%bkf(:,ikpt1) ) ) <= 1E-06 ) THEN
180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
              nrkpt = ikpt1
              EXIT
            END IF
          END DO
          IF( nrkpt .eq.0 ) STOP 'symm: Difference vector not found !'

          IF( list(nrkpt) .ne. 0) THEN
            list(nrkpt)   = 0
            neqvkpt(ikpt) = neqvkpt(ikpt) + 1
            parent(nrkpt) = ikpt
            symop(nrkpt)  = psym(iop)
          END IF
          IF ( all(list .eq. 0) ) EXIT

        END DO
      END DO

      ! for the Gamma-point holds:
      parent(1)  = 1
      neqvkpt(1) = 1

Daniel Wortmann's avatar
Daniel Wortmann committed
201

202 203 204

      ! determine number of members in the EIBZ(k)
      ic = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
205
      DO ikpt=1,kpts%nkptf
206 207 208 209 210 211
        IF(parent(ikpt) .eq. ikpt) ic = ic + 1
      END DO
      nkpt_EIBZ = ic

      ALLOCATE( pointer_EIBZ(nkpt_EIBZ) )
      ic = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
212
      DO ikpt=1,kpts%nkptf
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
        IF(parent(ikpt) .eq. ikpt) THEN
          ic = ic + 1
          pointer_EIBZ(ic) = ikpt
        END IF
      END DO

      IF ( irank2 == 0 ) THEN
        WRITE(6,'(A,i5)') ' Number of k-points in the EIBZ',nkpt_EIBZ
      END IF


      ! determine the factor n_q, that means the number of symmetrie operations of the little group of bk(:,nk)
      ! which keep q (in EIBZ) invariant

      ALLOCATE( n_q(nkpt_EIBZ) )

      ic  = 0
      n_q = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
231
      DO ikpt = 1,kpts%nkptf
232 233 234 235
        IF ( parent(ikpt) .eq. ikpt ) THEN
          ic = ic + 1
          DO iop=1,nsymop
            isym = psym(iop)
Daniel Wortmann's avatar
Daniel Wortmann committed
236
            rotkpt = matmul( rrot(:,:,isym), kpts%bkf(:,ikpt) )
237 238 239 240 241

            !transfer rotkpt into BZ
            rotkpt = modulo1(rotkpt,kpts%nkpt3)

            !check if rotkpt is identical to bk(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
242
            IF( maxval( abs( rotkpt - kpts%bkf(:,ikpt) ) ) .le. 1E-06) THEN
243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
              n_q(ic) = n_q(ic) + 1 
            END IF
          END DO
        END IF
      END DO
      IF( ic .ne. nkpt_EIBZ ) STOP 'symm: failure EIBZ'


      ! calculate degeneracy:
      ! degenerat(i) = 1 state i  is not degenerat,
      ! degenerat(i) = j state i has j-1 degenerat states at {i+1,...,i+j-1}
      ! degenerat(i) = 0 state i is degenerat

      tolerance = 1E-07 !0.00001

      degenerat = 1
      IF ( irank2 == 0 ) THEN
        WRITE(6,'(A,f10.8)') ' Tolerance for determining degenerate states=', tolerance
     END IF

Daniel Wortmann's avatar
Daniel Wortmann committed
263 264
      DO i=1,hybdat%nbands(nk)
        DO j=i+1,hybdat%nbands(nk)
265 266 267 268 269 270
          IF( abs(eig_irr(i,nk)-eig_irr(j,nk)) .le. tolerance) THEN
            degenerat(i) = degenerat(i) + 1
          END IF
        END DO
      END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
271
      DO i=1,hybdat%ne_eig(nk)
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286
        IF ( degenerat(i) .ne. 1 .or. degenerat(i) .ne. 0 ) THEN
          degenerat(i+1:i+degenerat(i)-1) = 0
        END IF
      END DO


      ! maximal number of degenerate bands -> maxndb
      maxndb = maxval(degenerat)

      ! number of different degenerate bands/states
      nddb   = count( degenerat .ge. 1)


      IF ( irank2 == 0 ) THEN
        WRITE(6,*) ' Degenerate states:'
Daniel Wortmann's avatar
Daniel Wortmann committed
287 288
        DO iband = 1,hybdat%nbands(nk)/5+1
          WRITE(6,'(5i5)')degenerat(iband*5-4:min(iband*5,hybdat%nbands(nk)))
289 290 291 292 293 294 295 296 297
        END DO
      END IF

      IF( irreps ) THEN
        ! calculate representation, i.e. the action of an element of
        ! the little group of k on \phi_n,k:
        ! P(R,T)\phi_n,k = \sum_{n'\in degenerat(n)} rep_v(n',n) *\phi_n',k
      
        ! read in cmt and z at current k-point (nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
298
         CALL read_cmt(cmt,nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
299
         call read_z(z,nk)
300 301 302 303 304

        ALLOCATE(rep_d(maxndb,nddb,nsymop),stat=ok )
        IF( ok .ne. 0) STOP 'symm: failure allocation rep_v'


Daniel Wortmann's avatar
Daniel Wortmann committed
305 306 307
        call olappw%alloc(z%l_real,lapw%nv(jsp),lapw%nv(jsp))
        ALLOCATE( olapmt(hybrid%maxindx,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype),stat=ok)
        IF( ok .ne. 0) STOP 'symm: failure allocation olapmt'
308 309 310 311 312

        olapmt = 0
        CALL wfolap_init (olappw,olapmt,gpt,&
     &                   atoms,hybrid,&
     &                    cell,&
Daniel Wortmann's avatar
Daniel Wortmann committed
313
     &                    hybdat%bas1,hybdat%bas2)
314 315 316 317 318 319 320 321 322 323 324




        ALLOCATE( cmthlp(hybrid%maxlmindx,atoms%nat,maxndb), cpwhlp(lapw%nv(jsp),maxndb),&
     &            stat= ok )
        IF( ok .ne. 0 ) STOP 'symm: failure allocation cmthlp/cpwhlp' 

        DO isym=1,nsymop 
          iop= psym(isym) 

Daniel Wortmann's avatar
Daniel Wortmann committed
325

326
          ic = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
327
          DO i=1,hybdat%nbands(nk)
328 329 330 331 332 333
            ndb = degenerat(i)
            IF ( ndb .ge. 1 ) THEN
              ic     = ic + 1
              cmthlp = 0
              cpwhlp = 0
              CALL waveftrafo_symm( &
Daniel Wortmann's avatar
Daniel Wortmann committed
334
     &                      cmthlp(:,:,:ndb),cpwhlp(:,:ndb),cmt,z%l_real,z%data_r,z%data_c,&
335 336 337 338 339 340 341 342 343
     &                      i,ndb,nk,iop,atoms,&
     &                      hybrid,kpts,&
     &                      sym,&
     &                      jsp,dimension,&
     &                      cell,gpt,lapw )

            
              DO iband = 1,ndb
                carr1 = cmt(iband+i-1,:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
344 345 346 347 348 349 350
                IF(z%l_real) THEN
                   rep_d(iband,ic,isym) = wfolap_inv( carr1,z%data_r(:lapw%nv(jsp),iband+i-1),cmthlp(:,:,iband),&
                        cpwhlp(:,iband), lapw%nv(jsp),lapw%nv(jsp),olappw%data_r,olapmt,atoms,hybrid)
                else
                   rep_d(iband,ic,isym) = wfolap_noinv( carr1,z%data_c(:lapw%nv(jsp),iband+i-1),cmthlp(:,:,iband),&
                        cpwhlp(:,iband), lapw%nv(jsp),lapw%nv(jsp),olappw%data_c,olapmt,atoms,hybrid)
                endif
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366
              END DO
            
            END IF
          END DO

        END DO

        DEALLOCATE( cmthlp,cpwhlp )


        ! calculate trace of irrecudible representation
        ALLOCATE( trace(sym%nsym,nddb), stat=ok )
        IF( ok .ne. 0) STOP 'symm: failure allocation trace'

        ic    = 0
        trace = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
367
        DO iband = 1,hybdat%nbands(nk)
368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388
          ndb = degenerat(iband)
          IF( ndb .ge. 1 ) THEN
            ic = ic + 1
            !calculate trace
            DO iop = 1,nsymop
              isym = psym(iop)
              DO i = 1,ndb
                trace(isym,ic) = trace(isym,ic) + rep_d(i,ic,iop)
              END DO
            END DO
          END IF
        END DO


        ! determine symmetry equivalent bands/irreducible representations by comparing the trace

        ALLOCATE( symequivalent(nddb,nddb), stat=ok )
        IF( ok .ne. 0) STOP 'symm: failure allocation symequivalent'

        ic1 = 0
        symequivalent = .false.
Daniel Wortmann's avatar
Daniel Wortmann committed
389
        DO iband1 = 1,hybdat%nbands(nk)
390 391 392 393
          ndb1 = degenerat(iband1)
          IF( ndb1 .ge. 1 ) THEN
            ic1 = ic1 + 1
            ic2 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
394
            DO iband2 = 1,hybdat%nbands(nk)
395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416
              ndb2 = degenerat(iband2)
              IF( ndb2 .ge. 1 ) THEN
                ic2 = ic2 + 1
                IF( ndb2 .eq. ndb1 ) THEN
                  ! note that this criterium is only valid for pure spatial rotations
                  ! if one combines spatial rotations with time reversal symmetry there
                  ! is no unique criteria to identify symequivalent state
                  ! however, also in the latter case the trace of the spatial rotations 
                  ! for two symmetry equivalent states must be equivalent
                  IF( all(abs(trace(:sym%nop,ic1)-trace(:sym%nop,ic2)) <= 1E-8))&
     &            THEN
                    symequivalent(ic2,ic1) = .true.
                  END IF
                END IF
              END IF
            END DO
          END IF
        END DO
      
        
      ELSE
        ! read in cmt and z at current k-point (nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
417 418 419
      
           CALL read_cmt(cmt,nk)
        !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
420 421 422 423 424 425 426 427 428 429 430 431
        
        IF( allocated(olapmt)) deallocate(olapmt)
        ALLOCATE( olapmt(hybrid%maxindx,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype),stat=ok)
        IF( ok .ne. 0) STOP 'symm: failure allocation olapmt'
        olapmt = 0

        DO itype = 1,atoms%ntype
          DO l = 0,atoms%lmax(itype)
            nn = hybrid%nindx(l,itype)
            DO n2 = 1,nn
              DO n1 = 1,nn
                olapmt(n1,n2,l,itype) = intgrf ( &
Daniel Wortmann's avatar
Daniel Wortmann committed
432 433 434
     &                        hybdat%bas1(:,n1,l,itype)*hybdat%bas1(:,n2,l,itype)&
     &                       +hybdat%bas2(:,n1,l,itype)*hybdat%bas2(:,n2,l,itype),&
     &                        atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,hybdat%gridf)
435 436 437 438 439
              END DO
            END DO
          END DO
        END DO
      
Daniel Wortmann's avatar
Daniel Wortmann committed
440
        ALLOCATE( wavefolap(hybdat%nbands(nk),hybdat%nbands(nk)),carr(hybrid%maxindx),stat=ok )
441 442 443 444 445 446 447 448 449 450 451
        IF( ok .ne. 0) STOP 'symm: failure allocation wfolap/maxindx'
        wavefolap = 0
        
        iatom  = 0
        DO itype = 1,atoms%ntype
          DO ieq = 1,atoms%neq(itype)
            iatom = iatom + 1
            lm = 0
            DO l = 0,atoms%lmax(itype)
              DO M = -l,l
                nn     = hybrid%nindx(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
452
                DO iband1 = 1,hybdat%nbands(nk)
453 454 455 456 457 458 459 460 461 462 463 464 465 466
                  carr(:nn) = matmul( olapmt(:nn,:nn,l,itype),&
     &                                cmt(iband1,lm+1:lm+nn,iatom)) 
                  DO iband2 = 1,iband1 
                    wavefolap(iband2,iband1)&
     &            = wavefolap(iband2,iband1)&
     &            + dot_product(cmt(iband2,lm+1:lm+nn,iatom),carr(:nn))
                  END DO
                END DO
                lm     = lm + nn
              END DO
            END DO
          END DO
        END DO
        
Daniel Wortmann's avatar
Daniel Wortmann committed
467
        DO iband1 = 1,hybdat%nbands(nk)
468 469 470 471 472 473 474 475 476
          DO iband2 = 1,iband1
            wavefolap(iband1,iband2) = conjg(wavefolap(iband2,iband1))
          END DO
        END DO
        
        ALLOCATE( symequivalent(nddb,nddb), stat=ok )
        IF( ok .ne. 0) STOP 'symm: failure allocation symequivalent'
        symequivalent = .false.
        ic1 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
477
        DO iband1 = 1,hybdat%nbands(nk)
478 479 480 481
          ndb1 = degenerat(iband1)
          IF( ndb1 .eq. 0 ) CYCLE
          ic1 = ic1 + 1
          ic2 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
482
          DO iband2 = 1,hybdat%nbands(nk)
483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502
            ndb2 = degenerat(iband2)
            IF( ndb2 .eq. 0 ) CYCLE
            ic2 = ic2 + 1
            IF( any( abs(wavefolap(iband1:iband1+ndb1-1,&
     &                             iband2:iband2+ndb2-1)) > 1E-9 ) )THEN
!     &          .and. ndb1 .eq. ndb2 ) THEN
              symequivalent(ic2,ic1) = .true.
            END IF
          END DO
        END DO
      END IF
      
      !
      ! generate index field which contain the band combinations (n1,n2),
      ! which are non zero 
      !   
      
      ic1       = 0
      indx_sest = 0
      nsest     = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
503
      DO iband1 = 1,hybdat%nbands(nk)
504 505 506 507 508 509 510 511
        ndb1 = degenerat(iband1)
        IF( ndb1 .ge. 1 ) ic1 = ic1 + 1
        i = 0
        DO WHILE ( degenerat(iband1-i) == 0 )
          i = i+1
        END DO
        ndb1 = degenerat(iband1-i)
        ic2 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
512
        DO iband2 = 1,hybdat%nbands(nk)
513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541
          ndb2 = degenerat(iband2)
          IF( ndb2 .ge. 1 ) ic2 = ic2 + 1
          i = 0
          DO WHILE ( degenerat(iband2-i) == 0 )
            i = i+1
          END DO
          ndb2 = degenerat(iband2-i)
          ! only upper triangular part
          IF( symequivalent(ic2,ic1) .and. iband2 .le. iband1 ) THEN
!            IF( ndb1 .ne. ndb2 ) STOP 'symm_hf: failure symequivalent'
            nsest(iband1)                   = nsest(iband1) + 1
            indx_sest(nsest(iband1),iband1) = iband2
          END IF
        END DO
      END DO

      !
      ! calculate representations for core states
      ! (note that for a core state, these are proportional to the Wigner D matrices)
      !
      ! Definition of the Wigner rotation matrices
      !
      !                     -1                l
      ! P(R) Y  (r)  = Y  (R  r)  =  sum     D    (R)  Y   (r)
      !       lm        lm              m'    m m'      lm'
      !

      pi = pimach()

Daniel Wortmann's avatar
Daniel Wortmann committed
542 543
      IF( hybdat%lmaxcd .gt. atoms%lmaxd ) STOP &
     & 'symm_hf: The very impropable case that hybdat%lmaxcd > atoms%lmaxd occurs'
544

Daniel Wortmann's avatar
Daniel Wortmann committed
545
      ALLOCATE(rep_c(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,nsymop,atoms%nat)&
546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562
     &        , stat=ok )
      IF( ok .ne. 0) STOP 'symm_hf: failure allocation rep_c'

      iatom  = 0
      iatom0 = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
          DO iop = 1,nsymop
            isym = psym(iop)
            IF( isym .le. sym%nop) THEN
              iisym = isym
            ELSE
              iisym = isym - sym%nop
            END IF

            ratom  = hybrid%map(iatom,isym)
Daniel Wortmann's avatar
Daniel Wortmann committed
563 564
            rotkpt = matmul( rrot(:,:,isym), kpts%bkf(:,nk) )
            g      = nint(rotkpt - kpts%bkf(:,nk))
565 566 567 568 569

            cdum   = exp(-2*pi*img*dot_product(rotkpt,sym%tau(:,iisym)))* &
     &               exp( 2*pi*img*dot_product(g,atoms%taual(:,ratom)))

            rep_c(:,:,:,iop,iatom) = &
Daniel Wortmann's avatar
Daniel Wortmann committed
570
     &         hybrid%d_wgn2(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,isym) * cdum
571 572 573 574 575 576 577 578 579 580
          END DO
        END DO
        iatom0 = iatom0 + atoms%neq(itype)
      END DO

      END SUBROUTINE symm_hf

      INTEGER FUNCTION symm_hf_nkpt_EIBZ(kpts,nk,sym)

      USE m_util, ONLY: modulo1
Daniel Wortmann's avatar
Daniel Wortmann committed
581
      USE m_types
582 583 584 585 586 587 588 589 590 591 592 593 594
      IMPLICIT NONE
      TYPE(t_sym),INTENT(IN)   :: sym
      TYPE(t_kpts),INTENT(IN)   :: kpts

!     - scalar input -
      INTEGER, INTENT(IN)   :: nk
!     - array input -

!     - local scalars -
      INTEGER               ::  isym,ic,iop,ikpt,ikpt1
      INTEGER               ::  nsymop,nrkpt
!     - local arrays -
      INTEGER               ::  rrot(3,3,sym%nsym)
Daniel Wortmann's avatar
Daniel Wortmann committed
595 596
      INTEGER               ::  neqvkpt(kpts%nkptf),list(kpts%nkptf),parent(kpts%nkptf),&
     &                          symop(kpts%nkptf)
597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617
      INTEGER, ALLOCATABLE  ::  psym(:)!,help(:)
      REAL                  ::  rotkpt(3)

      ! calculate rotations in reciprocal space
      DO isym = 1,sym%nsym
        IF( isym .le. sym%nop ) THEN
          rrot(:,:,isym) = transpose(sym%mrot(:,:,sym%invtab(isym)))
        ELSE
          rrot(:,:,isym) = -rrot(:,:,isym-sym%nop)
        END IF
      END DO

      ! determine little group of k., i.e. those symmetry operations
      ! which keep bk(:,nk,nw) invariant
      ! nsymop :: number of such symmetry-operations
      ! psym   :: points to the symmetry-operation

      ic = 0
      ALLOCATE(psym(sym%nsym))

      DO iop=1,sym%nsym
Daniel Wortmann's avatar
Daniel Wortmann committed
618
        rotkpt = matmul( rrot(:,:,iop), kpts%bkf(:,nk) )
619 620 621 622 623

        !transfer rotkpt into BZ
        rotkpt = modulo1(rotkpt,kpts%nkpt3)

        !check if rotkpt is identical to bk(:,nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
624
        IF( maxval( abs( rotkpt - kpts%bkf(:,nk) ) ) .le. 1E-07) THEN
625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645
          ic = ic + 1
          psym(ic) = iop
        END IF
      END DO
      nsymop = ic

      ! reallocate psym
!       ALLOCATE(help(ic))
!       help = psym(1:ic)
!       DEALLOCATE(psym)
!       ALLOCATE(psym(ic))
!       psym = help
!       DEALLOCATE(help)

      ! determine extented irreducible BZ of k ( EIBZ(k) ), i.e.
      ! those k-points, which can generate the whole BZ by
      ! applying the symmetry operations of the little group of k

      neqvkpt = 0

!       list = (/ (ikpt-1, ikpt=1,nkpt) /)
Daniel Wortmann's avatar
Daniel Wortmann committed
646
      DO ikpt=1,kpts%nkptf
647 648 649
        list(ikpt) = ikpt-1
      END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
650
      DO ikpt=2,kpts%nkptf
651 652
        DO iop=1,nsymop

Daniel Wortmann's avatar
Daniel Wortmann committed
653
          rotkpt = matmul( rrot(:,:,psym(iop)), kpts%bkf(:,ikpt) )
654 655 656 657 658 659

          !transfer rotkpt into BZ
          rotkpt = modulo1(rotkpt,kpts%nkpt3)

          !determine number of rotkpt
          nrkpt = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
660 661
          DO ikpt1=1,kpts%nkptf
            IF ( maxval( abs( rotkpt - kpts%bkf(:,ikpt1) ) ) .le. 1E-06 ) THEN
662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684
              nrkpt = ikpt1
              EXIT
            END IF
          END DO
          IF( nrkpt .eq.0 ) STOP 'symm: Difference vector not found !'

          IF( list(nrkpt) .ne. 0) THEN
            list(nrkpt)   = 0
            neqvkpt(ikpt) = neqvkpt(ikpt) + 1
            parent(nrkpt) = ikpt
            symop(nrkpt)  = psym(iop)
          END IF
          IF ( all(list .eq. 0) ) EXIT

        END DO
      END DO

      ! for the Gamma-point holds:
      parent(1)  = 1
      neqvkpt(1) = 1

      ! determine number of members in the EIBZ(k)
      ic = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
685
      DO ikpt=1,kpts%nkptf
686 687 688 689 690 691 692
        IF(parent(ikpt) .eq. ikpt) ic = ic + 1
      END DO
      symm_hf_nkpt_EIBZ = ic

      END FUNCTION symm_hf_nkpt_EIBZ

      END MODULE m_symm_hf