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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   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 31 32
     &                   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
Daniel Wortmann's avatar
Daniel Wortmann committed
33
      USE m_types
34 35
      IMPLICIT NONE

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

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

!     - scalars -
Daniel Wortmann's avatar
Daniel Wortmann committed
48
      INTEGER,INTENT(IN)              :: nk  
49 50 51 52 53 54 55 56
      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
57 58
      INTEGER,INTENT(OUT)             :: parent(kpts%nkptf)
      INTEGER,INTENT(OUT)             :: symop(kpts%nkptf)
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
      INTEGER,ALLOCATABLE,INTENT(OUT) :: pointer_EIBZ(:)
      INTEGER,ALLOCATABLE,INTENT(OUT) :: psym(:),n_q(:)      

Daniel Wortmann's avatar
Daniel Wortmann committed
64
      REAL,INTENT(IN)                 :: eig_irr(dimension%neigd,kpts%nkpt)
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
      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
86 87
      INTEGER                         :: neqvkpt(kpts%nkptf)
      INTEGER                         :: list(kpts%nkptf)
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
      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

Daniel Wortmann's avatar
Daniel Wortmann committed
136
        rotkpt = matmul( rrot(:,:,iop), kpts%bkf(:,nk) )
137 138 139 140 141

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

        !check if rotkpt is identical to bk(:,nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
142
        IF( maxval( abs( rotkpt - kpts%bkf(:,nk) ) ) .le. 1E-07) THEN
143 144 145 146 147 148 149 150 151
          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
152
        WRITE(6,'(A,3f10.5)') ' kpts%bkf(:,nk):',kpts%bkf(:,nk)
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
        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
170
      DO i=1,kpts%nkptf
171 172 173 174
        list(i) = i-1
      END DO

      symop = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
175
      DO ikpt=2,kpts%nkptf
176 177
        DO iop=1,nsymop

Daniel Wortmann's avatar
Daniel Wortmann committed
178
          rotkpt = matmul( rrot(:,:,psym(iop)), kpts%bkf(:,ikpt) )
179 180 181 182 183 184

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

          !determine number of rotkpt
          nrkpt = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
185 186
          DO ikpt1=1,kpts%nkptf
            IF ( maxval( abs( rotkpt - kpts%bkf(:,ikpt1) ) ) <= 1E-06 ) THEN
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
              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
Daniel Wortmann's avatar
Daniel Wortmann committed
209
      IF( sum(neqvkpt(:)) .ne. kpts%nkptf) THEN
210 211 212 213 214 215 216 217
        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
Daniel Wortmann's avatar
Daniel Wortmann committed
218
      DO ikpt=1,kpts%nkptf
219 220 221 222 223 224
        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
225
      DO ikpt=1,kpts%nkptf
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
        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
244
      DO ikpt = 1,kpts%nkptf
245 246 247 248
        IF ( parent(ikpt) .eq. ikpt ) THEN
          ic = ic + 1
          DO iop=1,nsymop
            isym = psym(iop)
Daniel Wortmann's avatar
Daniel Wortmann committed
249
            rotkpt = matmul( rrot(:,:,isym), kpts%bkf(:,ikpt) )
250 251 252 253 254

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

            !check if rotkpt is identical to bk(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
255
            IF( maxval( abs( rotkpt - kpts%bkf(:,ikpt) ) ) .le. 1E-06) THEN
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275
              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
              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
Daniel Wortmann's avatar
Daniel Wortmann committed
373
          rotkpt   = matmul(rrot(:,:,iop),kpts%bkf(:,nk))
374
          rotkpt   = modulo1(rotkpt,kpts%nkpt3)
Daniel Wortmann's avatar
Daniel Wortmann committed
375
          IF( maxval( abs( rotkpt - kpts%bkf(:,nk) ) ) .gt. 1e-08) THEN
376 377 378 379
            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
     &        , 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
618 619
            rotkpt = matmul( rrot(:,:,isym), kpts%bkf(:,nk) )
            g      = nint(rotkpt - kpts%bkf(:,nk))
620 621 622 623 624

            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
          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
636
      USE m_types
637 638 639 640 641 642 643 644 645 646 647 648 649
      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
650 651
      INTEGER               ::  neqvkpt(kpts%nkptf),list(kpts%nkptf),parent(kpts%nkptf),&
     &                          symop(kpts%nkptf)
652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672
      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
673
        rotkpt = matmul( rrot(:,:,iop), kpts%bkf(:,nk) )
674 675 676 677 678

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

        !check if rotkpt is identical to bk(:,nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
679
        IF( maxval( abs( rotkpt - kpts%bkf(:,nk) ) ) .le. 1E-07) THEN
680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700
          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
701
      DO ikpt=1,kpts%nkptf
702 703 704
        list(ikpt) = ikpt-1
      END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
705
      DO ikpt=2,kpts%nkptf
706 707
        DO iop=1,nsymop

Daniel Wortmann's avatar
Daniel Wortmann committed
708
          rotkpt = matmul( rrot(:,:,psym(iop)), kpts%bkf(:,ikpt) )
709 710 711 712 713 714

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

          !determine number of rotkpt
          nrkpt = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
715 716
          DO ikpt1=1,kpts%nkptf
            IF ( maxval( abs( rotkpt - kpts%bkf(:,ikpt1) ) ) .le. 1E-06 ) THEN
717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739
              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
740
      DO ikpt=1,kpts%nkptf
741 742 743 744 745 746 747
        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