symm_hf.F90 23.8 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
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 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
      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
137
        rotkpt = matmul( rrot(:,:,iop), kpts%bkf(:,nk) )
138 139 140 141 142

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

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

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

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

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

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

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

            !check if rotkpt is identical to bk(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
256
            IF( maxval( abs( rotkpt - kpts%bkf(:,ikpt) ) ) .le. 1E-06) THEN
257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276
              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
277 278
      DO i=1,hybdat%nbands(nk)
        DO j=i+1,hybdat%nbands(nk)
279 280 281 282 283 284
          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
285
      DO i=1,hybdat%ne_eig(nk)
286 287 288 289 290 291 292 293 294 295 296 297 298 299 300
        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
301 302
        DO iband = 1,hybdat%nbands(nk)/5+1
          WRITE(6,'(5i5)')degenerat(iband*5-4:min(iband*5,hybdat%nbands(nk)))
303 304 305 306 307 308 309 310 311
        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
312
         CALL read_cmt(cmt,nk)
313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
#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
336
     &                    hybdat%bas1,hybdat%bas2)
337 338 339 340 341 342


#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
343
        DO iband=1,hybdat%nbands(nk)
344 345
          IF(degenerat(iband) .ge. 1) THEN
            ndb = degenerat(iband)
Daniel Wortmann's avatar
Daniel Wortmann committed
346
            DO i= 1,hybdat%nbands(nk)!iband,iband+ndb-1
347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368
              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
369
          rotkpt   = matmul(rrot(:,:,iop),kpts%bkf(:,nk))
370
          rotkpt   = modulo1(rotkpt,kpts%nkpt3)
Daniel Wortmann's avatar
Daniel Wortmann committed
371
          IF( maxval( abs( rotkpt - kpts%bkf(:,nk) ) ) .gt. 1e-08) THEN
372 373 374 375
            STOP 'symm: error in psym'
          END IF
#endif
          ic = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
376
          DO i=1,hybdat%nbands(nk)
377 378 379 380 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
            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
414
        DO iband = 1,hybdat%nbands(nk)
415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435
          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
436
        DO iband1 = 1,hybdat%nbands(nk)
437 438 439 440
          ndb1 = degenerat(iband1)
          IF( ndb1 .ge. 1 ) THEN
            ic1 = ic1 + 1
            ic2 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
441
            DO iband2 = 1,hybdat%nbands(nk)
442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463
              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
464 465 466
      
           CALL read_cmt(cmt,nk)
        !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
467 468 469 470 471 472 473 474 475 476 477 478
        
        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
479 480 481
     &                        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)
482 483 484 485 486
              END DO
            END DO
          END DO
        END DO
      
Daniel Wortmann's avatar
Daniel Wortmann committed
487
        ALLOCATE( wavefolap(hybdat%nbands(nk),hybdat%nbands(nk)),carr(hybrid%maxindx),stat=ok )
488 489 490 491 492 493 494 495 496 497 498
        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
499
                DO iband1 = 1,hybdat%nbands(nk)
500 501 502 503 504 505 506 507 508 509 510 511 512 513
                  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
514
        DO iband1 = 1,hybdat%nbands(nk)
515 516 517 518 519 520 521 522 523
          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
524
        DO iband1 = 1,hybdat%nbands(nk)
525 526 527 528
          ndb1 = degenerat(iband1)
          IF( ndb1 .eq. 0 ) CYCLE
          ic1 = ic1 + 1
          ic2 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
529
          DO iband2 = 1,hybdat%nbands(nk)
530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549
            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
550
      DO iband1 = 1,hybdat%nbands(nk)
551 552 553 554 555 556 557 558
        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
559
        DO iband2 = 1,hybdat%nbands(nk)
560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588
          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
589 590
      IF( hybdat%lmaxcd .gt. atoms%lmaxd ) STOP &
     & 'symm_hf: The very impropable case that hybdat%lmaxcd > atoms%lmaxd occurs'
591

Daniel Wortmann's avatar
Daniel Wortmann committed
592
      ALLOCATE(rep_c(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,nsymop,atoms%nat)&
593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609
     &        , 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
610 611
            rotkpt = matmul( rrot(:,:,isym), kpts%bkf(:,nk) )
            g      = nint(rotkpt - kpts%bkf(:,nk))
612 613 614 615 616

            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
617
     &         hybrid%d_wgn2(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,isym) * cdum
618 619 620 621 622 623 624 625 626 627
          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
628
      USE m_types
629 630 631 632 633 634 635 636 637 638 639 640 641
      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
642 643
      INTEGER               ::  neqvkpt(kpts%nkptf),list(kpts%nkptf),parent(kpts%nkptf),&
     &                          symop(kpts%nkptf)
644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664
      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
665
        rotkpt = matmul( rrot(:,:,iop), kpts%bkf(:,nk) )
666 667 668 669 670

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

        !check if rotkpt is identical to bk(:,nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
671
        IF( maxval( abs( rotkpt - kpts%bkf(:,nk) ) ) .le. 1E-07) THEN
672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692
          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
693
      DO ikpt=1,kpts%nkptf
694 695 696
        list(ikpt) = ikpt-1
      END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
697
      DO ikpt=2,kpts%nkptf
698 699
        DO iop=1,nsymop

Daniel Wortmann's avatar
Daniel Wortmann committed
700
          rotkpt = matmul( rrot(:,:,psym(iop)), kpts%bkf(:,ikpt) )
701 702 703 704 705 706

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

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