symm_hf.F90 24.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   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

      SUBROUTINE symm_hf(kpts,nkpti,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 31 32 33 34 35
     &                   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
      USE m_olap   ,ONLY: wfolap,wfolap1,wfolap_init
      USE m_trafo  ,ONLY: waveftrafo_symm
    USE m_types
      IMPLICIT NONE

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

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

!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
48
      INTEGER,INTENT(IN)              :: nkpti  ,nk  
49 50 51 52 53 54 55 56 57 58
      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))
      INTEGER,INTENT(OUT)             :: parent(kpts%nkpt)
      INTEGER,INTENT(OUT)             :: symop(kpts%nkpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
59 60
      INTEGER,INTENT(INOUT)           :: degenerat(hybdat%ne_eig(nk))
      INTEGER,INTENT(OUT)             :: nsest(hybdat%nbands(nk)), indx_sest(hybdat%nbands(nk),hybdat%nbands(nk))  
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 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 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275
      INTEGER,ALLOCATABLE,INTENT(OUT) :: pointer_EIBZ(:)
      INTEGER,ALLOCATABLE,INTENT(OUT) :: psym(:),n_q(:)      

      REAL,INTENT(IN)                 :: eig_irr(dimension%neigd,nkpti)
      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)
      INTEGER                         :: neqvkpt(kpts%nkpt)
      INTEGER                         :: list(kpts%nkpt)
      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(:,:)

#if ( defined(CPP_INVERSION) )
      REAL                            :: carr2(lapw%nv(jsp))
      REAL                            :: z(dimension%nbasfcn,dimension%neigd)
      REAL,ALLOCATABLE                :: olappw(:,:)
#else
      COMPLEX                         :: carr2(lapw%nv(jsp))
      COMPLEX                         :: z(dimension%nbasfcn,dimension%neigd)
      COMPLEX,ALLOCATABLE             :: olappw(:,:)
#endif 
      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

        rotkpt = matmul( rrot(:,:,iop), kpts%bk(:,nk) )

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

        !check if rotkpt is identical to bk(:,nk)
        IF( maxval( abs( rotkpt - kpts%bk(:,nk) ) ) .le. 1E-07) THEN
          ic = ic + 1
          psym(ic) = iop
        END IF
      END DO
      nsymop = ic
     

      IF ( irank2 == 0 ) THEN
        WRITE(6,'(A,i3)') ' nk',nk
        WRITE(6,'(A,3f10.5)') ' kpts%bk(:,nk):',kpts%bk(:,nk)
        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

      DO i=1,kpts%nkpt
        list(i) = i-1
      END DO

      symop = 0
      DO ikpt=2,kpts%nkpt
        DO iop=1,nsymop

          rotkpt = matmul( rrot(:,:,psym(iop)), kpts%bk(:,ikpt) )

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

          !determine number of rotkpt
          nrkpt = 0
          DO ikpt1=1,kpts%nkpt
            IF ( maxval( abs( rotkpt - kpts%bk(:,ikpt1) ) ) <= 1E-06 ) THEN
              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

#ifdef CPP_DEBUG
      IF( sum(neqvkpt(:)) .ne. kpts%nkpt) THEN
        IF ( mpi%irank == 0 ) WRITE(6,'(A,i3)') ' Check EIBZ(k) by summing&&
     & up equivalent k-points:',sum(neqvkpt(:))
        STOP 'symm: neqvkpt not identical to nkpt'
      END IF
#endif

      ! determine number of members in the EIBZ(k)
      ic = 0
      DO ikpt=1,kpts%nkpt
        IF(parent(ikpt) .eq. ikpt) ic = ic + 1
      END DO
      nkpt_EIBZ = ic

      ALLOCATE( pointer_EIBZ(nkpt_EIBZ) )
      ic = 0
      DO ikpt=1,kpts%nkpt
        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
      DO ikpt = 1,kpts%nkpt
        IF ( parent(ikpt) .eq. ikpt ) THEN
          ic = ic + 1
          DO iop=1,nsymop
            isym = psym(iop)
            rotkpt = matmul( rrot(:,:,isym), kpts%bk(:,ikpt) )

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

            !check if rotkpt is identical to bk(:,ikpt)
            IF( maxval( abs( rotkpt - kpts%bk(:,ikpt) ) ) .le. 1E-06) THEN
              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
276 277
      DO i=1,hybdat%nbands(nk)
        DO j=i+1,hybdat%nbands(nk)
278 279 280 281 282 283
          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
284
      DO i=1,hybdat%ne_eig(nk)
285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
        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
300 301
        DO iband = 1,hybdat%nbands(nk)/5+1
          WRITE(6,'(5i5)')degenerat(iband*5-4:min(iband*5,hybdat%nbands(nk)))
302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339
        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)
        irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16
        OPEN(unit=777,file='cmt',form='unformatted',access='direct',&
     &       recl=irecl_cmt)
        READ(777,rec=nk) cmt
        CLOSE(777)

#ifdef CPP_INVERSION
        irecl_z   =  dimension%nbasfcn*dimension%neigd*8
#else
        irecl_z   =  dimension%nbasfcn*dimension%neigd*16
#endif
        OPEN(unit=778,file='z',form='unformatted',access='direct',&
     &      recl=irecl_z)
        READ(778,rec=nk) z
        CLOSE(778)

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


        ALLOCATE( olappw(lapw%nv(jsp),lapw%nv(jsp)),&
     &            olapmt(hybrid%maxindx,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype),stat=ok)
        IF( ok .ne. 0) STOP 'symm: failure allocation olappw/olapmt'

        olappw = 0
        olapmt = 0
        CALL wfolap_init (olappw,olapmt,gpt,&
     &                   atoms,hybrid,&
     &                    cell,&
Daniel Wortmann's avatar
Daniel Wortmann committed
340
     &                    hybdat%bas1,hybdat%bas2)
341 342 343 344 345 346


#ifdef CPP_DEBUG1
        ! check orthogonality of wavefunctions
        OPEN(677,file='ortho')
        WRITE(677,*) 'iteration',it,'k-point',nk
Daniel Wortmann's avatar
Daniel Wortmann committed
347
        DO iband=1,hybdat%nbands(nk)
348 349
          IF(degenerat(iband) .ge. 1) THEN
            ndb = degenerat(iband)
Daniel Wortmann's avatar
Daniel Wortmann committed
350
            DO i= 1,hybdat%nbands(nk)!iband,iband+ndb-1
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379
              WRITE(677,*) iband,i,&
     &          wfolap1( cmt(iband,:,:),z(:lapw%nv(jsp),iband),&
     &           cmt(i,:,:),z(:lapw%nv(jsp),i),&
     &          lapw%nv(jsp),lapw%nv(jsp),olappw,olapmt,&
     &          atoms%ntype,atoms%neq,atoms%nat,atoms%lmax,atoms%lmaxd,hybrid%nindx,hybrid%maxindx,hybrid%maxlmindx)


            END DO
          
          END IF
        END DO
        CLOSE(677)
#endif

        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) 

#ifdef CPP_DEBUG
          rotkpt   = matmul(rrot(:,:,iop),kpts%bk(:,nk))
          rotkpt   = modulo1(rotkpt,kpts%nkpt3)
          IF( maxval( abs( rotkpt - kpts%bk(:,nk) ) ) .gt. 1e-08) THEN
            STOP 'symm: error in psym'
          END IF
#endif
          ic = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
380
          DO i=1,hybdat%nbands(nk)
381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417
            ndb = degenerat(i)
            IF ( ndb .ge. 1 ) THEN
              ic     = ic + 1
              cmthlp = 0
              cpwhlp = 0
              STOP "WAVETRAVO_SYM real/complex"
              CALL waveftrafo_symm( &
     &                      cmthlp(:,:,:ndb),cpwhlp(:,:ndb),cmt,.false.,(/0.1/),z,&
     &                      i,ndb,nk,iop,atoms,&
     &                      hybrid,kpts,&
     &                      sym,&
     &                      jsp,dimension,&
     &                      cell,gpt,lapw )

            
              DO iband = 1,ndb
                carr1 = cmt(iband+i-1,:,:)
                carr2 = z(:lapw%nv(jsp),iband+i-1)
                rep_d(iband,ic,isym) &
     &        = wfolap( carr1,carr2,cmthlp(:,:,iband),cpwhlp(:,iband),&
     &                  lapw%nv(jsp),lapw%nv(jsp),olappw,olapmt,atoms,hybrid)
              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
418
        DO iband = 1,hybdat%nbands(nk)
419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439
          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
440
        DO iband1 = 1,hybdat%nbands(nk)
441 442 443 444
          ndb1 = degenerat(iband1)
          IF( ndb1 .ge. 1 ) THEN
            ic1 = ic1 + 1
            ic2 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
445
            DO iband2 = 1,hybdat%nbands(nk)
446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473
              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)
        irecl_cmt = dimension%neigd*hybrid%maxlmindx*atoms%nat*16
        OPEN(unit=777,file='cmt',form='unformatted',access='direct',&
     &       recl=irecl_cmt)
        READ(777,rec=nk) cmt
        CLOSE(777)
        
Daniel Wortmann's avatar
Daniel Wortmann committed
474
        CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
475 476 477 478 479 480 481 482 483 484 485 486
        
        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
487 488 489
     &                        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)
490 491 492 493 494
              END DO
            END DO
          END DO
        END DO
      
Daniel Wortmann's avatar
Daniel Wortmann committed
495
        ALLOCATE( wavefolap(hybdat%nbands(nk),hybdat%nbands(nk)),carr(hybrid%maxindx),stat=ok )
496 497 498 499 500 501 502 503 504 505 506
        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
507
                DO iband1 = 1,hybdat%nbands(nk)
508 509 510 511 512 513 514 515 516 517 518 519 520 521
                  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
522
        DO iband1 = 1,hybdat%nbands(nk)
523 524 525 526 527 528 529 530 531
          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
532
        DO iband1 = 1,hybdat%nbands(nk)
533 534 535 536
          ndb1 = degenerat(iband1)
          IF( ndb1 .eq. 0 ) CYCLE
          ic1 = ic1 + 1
          ic2 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
537
          DO iband2 = 1,hybdat%nbands(nk)
538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557
            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
558
      DO iband1 = 1,hybdat%nbands(nk)
559 560 561 562 563 564 565 566
        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
567
        DO iband2 = 1,hybdat%nbands(nk)
568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596
          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
597 598
      IF( hybdat%lmaxcd .gt. atoms%lmaxd ) STOP &
     & 'symm_hf: The very impropable case that hybdat%lmaxcd > atoms%lmaxd occurs'
599

Daniel Wortmann's avatar
Daniel Wortmann committed
600
      ALLOCATE(rep_c(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,nsymop,atoms%nat)&
601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624
     &        , 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)
            rotkpt = matmul( rrot(:,:,isym), kpts%bk(:,nk) )
            g      = nint(rotkpt - kpts%bk(:,nk))

            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
625
     &         sym%d_wgn(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,isym) * cdum
626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747
          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
    USE m_types
      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)
      INTEGER               ::  neqvkpt(kpts%nkpt),list(kpts%nkpt),parent(kpts%nkpt),&
     &                          symop(kpts%nkpt)
      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
        rotkpt = matmul( rrot(:,:,iop), kpts%bk(:,nk) )

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

        !check if rotkpt is identical to bk(:,nk)
        IF( maxval( abs( rotkpt - kpts%bk(:,nk) ) ) .le. 1E-07) THEN
          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) /)
      DO ikpt=1,kpts%nkpt
        list(ikpt) = ikpt-1
      END DO

      DO ikpt=2,kpts%nkpt
        DO iop=1,nsymop

          rotkpt = matmul( rrot(:,:,psym(iop)), kpts%bk(:,ikpt) )

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

          !determine number of rotkpt
          nrkpt = 0
          DO ikpt1=1,kpts%nkpt
            IF ( maxval( abs( rotkpt - kpts%bk(:,ikpt1) ) ) .le. 1E-06 ) THEN
              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
      DO ikpt=1,kpts%nkpt
        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