LCOV - code coverage report
Current view: top level - hybrid - symm_hf.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 208 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 3 0.0 %

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !   This module generates the little group of k and the extended irr. !
       3             : !   BZ. Furthermore it calculates the irr. representation             !
       4             : !                                                                     !
       5             : !   P(R,T)\phi_n,k = \sum_{n'} rep_v(n',n) *\phi_n',k        !
       6             : !   where                                                             !
       7             : !         P  is an element of the little group of k                   !
       8             : !         n' runs over group of degenerat states belonging to n.      !
       9             : !                                                                     !
      10             : !                                             M.Betzinger (09/07)     !
      11             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      12             : MODULE m_symm_hf
      13             : 
      14             : #define irreps .false.
      15             : 
      16             : CONTAINS
      17             : 
      18           0 :    SUBROUTINE symm_hf_init(sym, kpts, nk, nsymop, rrot, psym)
      19             : 
      20             :       USE m_types
      21             :       USE m_util, ONLY: modulo1
      22             : 
      23             :       IMPLICIT NONE
      24             : 
      25             :       TYPE(t_sym), INTENT(IN)    :: sym
      26             :       TYPE(t_kpts), INTENT(IN)    :: kpts
      27             :       INTEGER, INTENT(IN)    :: nk
      28             :       INTEGER, INTENT(OUT)   :: nsymop
      29             :       INTEGER, INTENT(INOUT) :: rrot(3, 3, sym%nsym)
      30             :       INTEGER, INTENT(INOUT) :: psym(sym%nsym) ! Note: psym is only filled up to index nsymop
      31             : 
      32             :       INTEGER :: i
      33             :       REAL    :: rotkpt(3)
      34             : 
      35             :       ! calculate rotations in reciprocal space
      36           0 :       DO i = 1, sym%nsym
      37           0 :          IF (i <= sym%nop) THEN
      38           0 :             rrot(:, :, i) = transpose(sym%mrot(:, :, sym%invtab(i)))
      39             :          ELSE
      40           0 :             rrot(:, :, i) = -rrot(:, :, i - sym%nop)
      41             :          END IF
      42             :       END DO
      43             : 
      44             :       ! determine little group of k., i.e. those symmetry operations
      45             :       ! which keep bk(:,nk) invariant
      46             :       ! nsymop :: number of such symmetry-operations
      47             :       ! psym   :: points to the symmetry-operation
      48             : 
      49           0 :       psym = 0
      50           0 :       nsymop = 0
      51           0 :       DO i = 1, sym%nsym
      52           0 :          rotkpt = matmul(rrot(:, :, i), kpts%bkf(:, nk))
      53             : 
      54             :          !transfer rotkpt into BZ
      55           0 :          rotkpt = modulo1(rotkpt, kpts%nkpt3)
      56             : 
      57             :          !check if rotkpt is identical to bk(:,nk)
      58           0 :          IF (maxval(abs(rotkpt - kpts%bkf(:, nk))) <= 1E-07) THEN
      59           0 :             nsymop = nsymop + 1
      60           0 :             psym(nsymop) = i
      61             :          END IF
      62             :       END DO
      63             : 
      64           0 :       WRITE (6, '(A,i3)') ' nk', nk
      65           0 :       WRITE (6, '(A,3f10.5)') ' kpts%bkf(:,nk):', kpts%bkf(:, nk)
      66           0 :       WRITE (6, '(A,i3)') ' Number of elements in the little group:', nsymop
      67             : 
      68           0 :    END SUBROUTINE symm_hf_init
      69             : 
      70           0 :    SUBROUTINE symm_hf(kpts, nk, sym, dimension, hybdat, eig_irr, atoms, hybrid, cell, &
      71           0 :                       lapw, jsp, mpi, rrot, nsymop, psym, nkpt_EIBZ, n_q, parent, &
      72           0 :                       pointer_EIBZ, nsest, indx_sest)
      73             : 
      74             :       USE m_constants
      75             :       USE m_types
      76             :       USE m_util, ONLY: modulo1, intgrf, intgrf_init
      77             :       USE m_olap, ONLY: wfolap_inv, wfolap_noinv, wfolap1, wfolap_init
      78             :       USE m_trafo, ONLY: waveftrafo_symm
      79             :       USE m_io_hybrid
      80             :       IMPLICIT NONE
      81             : 
      82             :       TYPE(t_hybdat), INTENT(IN)   :: hybdat
      83             : 
      84             :       TYPE(t_mpi), INTENT(IN)   :: mpi
      85             :       TYPE(t_dimension), INTENT(IN)   :: dimension
      86             :       TYPE(t_hybrid), INTENT(IN) :: hybrid
      87             :       TYPE(t_sym), INTENT(IN)    :: sym
      88             :       TYPE(t_cell), INTENT(IN)   :: cell
      89             :       TYPE(t_kpts), INTENT(IN)   :: kpts
      90             :       TYPE(t_atoms), INTENT(IN)  :: atoms
      91             :       TYPE(t_lapw), INTENT(IN)   :: lapw
      92             : 
      93             : !     - scalars -
      94             :       INTEGER, INTENT(IN)              :: nk
      95             :       INTEGER, INTENT(IN)              :: jsp
      96             :       INTEGER, INTENT(OUT)             :: nkpt_EIBZ
      97             :       INTEGER, INTENT(IN)              :: nsymop
      98             : 
      99             : !     - arrays -
     100             :       INTEGER, INTENT(IN)              :: rrot(3, 3, sym%nsym)
     101             :       INTEGER, INTENT(IN)              :: psym(sym%nsym)
     102             :       INTEGER, INTENT(OUT)             :: parent(kpts%nkptf)
     103             :       INTEGER, INTENT(OUT)             :: nsest(hybrid%nbands(nk)), indx_sest(hybrid%nbands(nk), hybrid%nbands(nk))
     104             :       INTEGER, ALLOCATABLE, INTENT(OUT) :: pointer_EIBZ(:)
     105             :       INTEGER, ALLOCATABLE, INTENT(OUT) :: n_q(:)
     106             : 
     107             :       REAL, INTENT(IN)                 :: eig_irr(dimension%neigd, kpts%nkpt)
     108             : 
     109             : !     - local scalars -
     110             :       INTEGER                         :: ikpt, ikpt1, iop, isym, iisym, m
     111             :       INTEGER                         :: itype, ieq, iatom, ratom
     112             :       INTEGER                         :: iband, iband1, iband2, iatom0
     113             :       INTEGER                         :: i, j, ic, ic1, ic2
     114             :       INTEGER                         :: irecl_cmt, irecl_z
     115             :       INTEGER                         :: ok
     116             :       INTEGER                         :: l, lm
     117             :       INTEGER                         :: n1, n2, nn
     118             :       INTEGER                         :: ndb, ndb1, ndb2
     119             :       INTEGER                         :: nrkpt
     120             :       INTEGER                         :: maxndb, nddb
     121             : 
     122             :       REAL                            :: tolerance, pi
     123             : 
     124             :       COMPLEX                         :: cdum
     125             :       COMPLEX, PARAMETER             :: img = (0.0, 1.0)
     126             : 
     127             : !     - local arrays -
     128           0 :       INTEGER                         :: neqvkpt(kpts%nkptf)
     129           0 :       INTEGER                         :: list(kpts%nkptf)
     130           0 :       INTEGER                         :: degenerat(hybrid%ne_eig(nk))
     131             :       INTEGER, ALLOCATABLE             :: help(:)
     132             : 
     133             :       REAL                            :: rotkpt(3), g(3)
     134           0 :       REAL, ALLOCATABLE             :: olapmt(:, :, :, :)
     135             : 
     136           0 :       COMPLEX                         :: cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
     137           0 :       COMPLEX                         :: carr1(hybrid%maxlmindx, atoms%nat)
     138           0 :       COMPLEX, ALLOCATABLE             :: carr(:), wavefolap(:, :)
     139           0 :       COMPLEX, ALLOCATABLE             :: cmthlp(:, :, :)
     140           0 :       COMPLEX, ALLOCATABLE             :: cpwhlp(:, :)
     141           0 :       COMPLEX, ALLOCATABLE             :: trace(:, :)
     142             : 
     143           0 :       TYPE(t_mat)                     :: olappw, z
     144           0 :       COMPLEX, ALLOCATABLE             :: rep_d(:, :, :)
     145           0 :       LOGICAL, ALLOCATABLE             :: symequivalent(:, :)
     146             : 
     147           0 :       WRITE (6, '(A)') new_line('n')//new_line('n')//'### subroutine: symm ###'
     148             : 
     149             :       ! determine extented irreducible BZ of k ( EIBZ(k) ), i.e.
     150             :       ! those k-points, which can generate the whole BZ by
     151             :       ! applying the symmetry operations of the little group of k
     152             : 
     153           0 :       neqvkpt = 0
     154             : 
     155           0 :       DO i = 1, kpts%nkptf
     156           0 :          list(i) = i - 1
     157             :       END DO
     158             : 
     159           0 :       DO ikpt = 2, kpts%nkptf
     160           0 :          DO iop = 1, nsymop
     161             : 
     162           0 :             rotkpt = matmul(rrot(:, :, psym(iop)), kpts%bkf(:, ikpt))
     163             : 
     164             :             !transfer rotkpt into BZ
     165           0 :             rotkpt = modulo1(rotkpt, kpts%nkpt3)
     166             : 
     167             :             !determine number of rotkpt
     168           0 :             nrkpt = 0
     169           0 :             DO ikpt1 = 1, kpts%nkptf
     170           0 :                IF (maxval(abs(rotkpt - kpts%bkf(:, ikpt1))) <= 1E-06) THEN
     171             :                   nrkpt = ikpt1
     172             :                   EXIT
     173             :                END IF
     174             :             END DO
     175           0 :             IF (nrkpt == 0) STOP 'symm: Difference vector not found !'
     176             : 
     177           0 :             IF (list(nrkpt) /= 0) THEN
     178           0 :                list(nrkpt) = 0
     179           0 :                neqvkpt(ikpt) = neqvkpt(ikpt) + 1
     180           0 :                parent(nrkpt) = ikpt
     181             :             END IF
     182           0 :             IF (all(list == 0)) EXIT
     183             : 
     184             :          END DO
     185             :       END DO
     186             : 
     187             :       ! for the Gamma-point holds:
     188           0 :       parent(1) = 1
     189           0 :       neqvkpt(1) = 1
     190             : 
     191             :       ! determine number of members in the EIBZ(k)
     192           0 :       ic = 0
     193           0 :       DO ikpt = 1, kpts%nkptf
     194           0 :          IF (parent(ikpt) == ikpt) ic = ic + 1
     195             :       END DO
     196           0 :       nkpt_EIBZ = ic
     197             : 
     198           0 :       ALLOCATE (pointer_EIBZ(nkpt_EIBZ))
     199           0 :       ic = 0
     200           0 :       DO ikpt = 1, kpts%nkptf
     201           0 :          IF (parent(ikpt) == ikpt) THEN
     202           0 :             ic = ic + 1
     203           0 :             pointer_EIBZ(ic) = ikpt
     204             :          END IF
     205             :       END DO
     206             : 
     207           0 :       WRITE (6, '(A,i5)') ' Number of k-points in the EIBZ', nkpt_EIBZ
     208             : 
     209             :       ! determine the factor n_q, that means the number of symmetrie operations of the little group of bk(:,nk)
     210             :       ! which keep q (in EIBZ) invariant
     211             : 
     212           0 :       ALLOCATE (n_q(nkpt_EIBZ))
     213             : 
     214           0 :       ic = 0
     215           0 :       n_q = 0
     216           0 :       DO ikpt = 1, kpts%nkptf
     217           0 :          IF (parent(ikpt) == ikpt) THEN
     218           0 :             ic = ic + 1
     219           0 :             DO iop = 1, nsymop
     220           0 :                isym = psym(iop)
     221           0 :                rotkpt = matmul(rrot(:, :, isym), kpts%bkf(:, ikpt))
     222             : 
     223             :                !transfer rotkpt into BZ
     224           0 :                rotkpt = modulo1(rotkpt, kpts%nkpt3)
     225             : 
     226             :                !check if rotkpt is identical to bk(:,ikpt)
     227           0 :                IF (maxval(abs(rotkpt - kpts%bkf(:, ikpt))) <= 1E-06) THEN
     228           0 :                   n_q(ic) = n_q(ic) + 1
     229             :                END IF
     230             :             END DO
     231             :          END IF
     232             :       END DO
     233           0 :       IF (ic /= nkpt_EIBZ) STOP 'symm: failure EIBZ'
     234             : 
     235             :       ! calculate degeneracy:
     236             :       ! degenerat(i) = 1 state i  is not degenerat,
     237             :       ! degenerat(i) = j state i has j-1 degenerat states at {i+1,...,i+j-1}
     238             :       ! degenerat(i) = 0 state i is degenerat
     239             : 
     240           0 :       tolerance = 1E-07 !0.00001
     241             : 
     242           0 :       degenerat = 1
     243             : 
     244           0 :       WRITE (6, '(A,f10.8)') ' Tolerance for determining degenerate states=', tolerance
     245             : 
     246           0 :       DO i = 1, hybrid%nbands(nk)
     247           0 :          DO j = i + 1, hybrid%nbands(nk)
     248           0 :             IF (abs(eig_irr(i, nk) - eig_irr(j, nk)) <= tolerance) THEN
     249           0 :                degenerat(i) = degenerat(i) + 1
     250             :             END IF
     251             :          END DO
     252             :       END DO
     253             : 
     254           0 :       DO i = 1, hybrid%ne_eig(nk)
     255           0 :          IF (degenerat(i) /= 1 .or. degenerat(i) /= 0) THEN
     256           0 :             degenerat(i + 1:i + degenerat(i) - 1) = 0
     257             :          END IF
     258             :       END DO
     259             : 
     260             :       ! maximal number of degenerate bands -> maxndb
     261             :       maxndb = maxval(degenerat)
     262             : 
     263             :       ! number of different degenerate bands/states
     264           0 :       nddb = count(degenerat >= 1)
     265             : 
     266           0 :       WRITE (6, *) ' Degenerate states:'
     267           0 :       DO iband = 1, hybrid%nbands(nk)/5 + 1
     268           0 :          WRITE (6, '(5i5)') degenerat(iband*5 - 4:min(iband*5, hybrid%nbands(nk)))
     269             :       END DO
     270             : 
     271             :       IF (irreps) THEN
     272             :          ! calculate representation, i.e. the action of an element of
     273             :          ! the little group of k on \phi_n,k:
     274             :          ! P(R,T)\phi_n,k = \sum_{n'\in degenerat(n)} rep_v(n',n) *\phi_n',k
     275             : 
     276             :          ! read in cmt and z at current k-point (nk)
     277             :          CALL read_cmt(cmt, nk)
     278             :          call read_z(z, nk)
     279             : 
     280             :          ALLOCATE (rep_d(maxndb, nddb, nsymop), stat=ok)
     281             :          IF (ok /= 0) STOP 'symm: failure allocation rep_v'
     282             : 
     283             :          call olappw%alloc(z%l_real, lapw%nv(jsp), lapw%nv(jsp))
     284             :          ALLOCATE (olapmt(hybrid%maxindx, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype), stat=ok)
     285             :          IF (ok /= 0) STOP 'symm: failure allocation olapmt'
     286             : 
     287             :          olapmt = 0
     288             :          CALL wfolap_init(olappw, olapmt, lapw%gvec(:, :, jsp), atoms, hybrid, cell, hybdat%bas1, hybdat%bas2)
     289             : 
     290             :          ALLOCATE (cmthlp(hybrid%maxlmindx, atoms%nat, maxndb), cpwhlp(lapw%nv(jsp), maxndb), stat=ok)
     291             :          IF (ok /= 0) STOP 'symm: failure allocation cmthlp/cpwhlp'
     292             : 
     293             :          DO isym = 1, nsymop
     294             :             iop = psym(isym)
     295             : 
     296             :             ic = 0
     297             :             DO i = 1, hybrid%nbands(nk)
     298             :                ndb = degenerat(i)
     299             :                IF (ndb >= 1) THEN
     300             :                   ic = ic + 1
     301             :                   cmthlp = 0
     302             :                   cpwhlp = 0
     303             : 
     304             :                   CALL waveftrafo_symm(cmthlp(:, :, :ndb), cpwhlp(:, :ndb), cmt, z%l_real, z%data_r, z%data_c, &
     305             :                                        i, ndb, nk, iop, atoms, hybrid, kpts, sym, jsp, dimension, cell, lapw)
     306             : 
     307             :                   DO iband = 1, ndb
     308             :                      carr1 = cmt(iband + i - 1, :, :)
     309             :                      IF (z%l_real) THEN
     310             :                         rep_d(iband, ic, isym) = wfolap_inv(carr1, z%data_r(:lapw%nv(jsp), iband + i - 1), cmthlp(:, :, iband), &
     311             :                                                             cpwhlp(:, iband), lapw%nv(jsp), lapw%nv(jsp), olappw%data_r, olapmt, atoms, hybrid)
     312             :                      else
     313             :                         rep_d(iband, ic, isym) = wfolap_noinv(carr1, z%data_c(:lapw%nv(jsp), iband + i - 1), cmthlp(:, :, iband), &
     314             :                                                               cpwhlp(:, iband), lapw%nv(jsp), lapw%nv(jsp), olappw%data_c, olapmt, atoms, hybrid)
     315             :                      endif
     316             :                   END DO
     317             : 
     318             :                END IF
     319             :             END DO
     320             : 
     321             :          END DO
     322             : 
     323             :          DEALLOCATE (cmthlp, cpwhlp)
     324             : 
     325             :          ! calculate trace of irrecudible representation
     326             :          ALLOCATE (trace(sym%nsym, nddb), stat=ok)
     327             :          IF (ok /= 0) STOP 'symm: failure allocation trace'
     328             : 
     329             :          ic = 0
     330             :          trace = 0
     331             :          DO iband = 1, hybrid%nbands(nk)
     332             :             ndb = degenerat(iband)
     333             :             IF (ndb >= 1) THEN
     334             :                ic = ic + 1
     335             :                !calculate trace
     336             :                DO iop = 1, nsymop
     337             :                   isym = psym(iop)
     338             :                   DO i = 1, ndb
     339             :                      trace(isym, ic) = trace(isym, ic) + rep_d(i, ic, iop)
     340             :                   END DO
     341             :                END DO
     342             :             END IF
     343             :          END DO
     344             : 
     345             :          ! determine symmetry equivalent bands/irreducible representations by comparing the trace
     346             : 
     347             :          ALLOCATE (symequivalent(nddb, nddb), stat=ok)
     348             :          IF (ok /= 0) STOP 'symm: failure allocation symequivalent'
     349             : 
     350             :          ic1 = 0
     351             :          symequivalent = .false.
     352             :          DO iband1 = 1, hybrid%nbands(nk)
     353             :             ndb1 = degenerat(iband1)
     354             :             IF (ndb1 >= 1) THEN
     355             :                ic1 = ic1 + 1
     356             :                ic2 = 0
     357             :                DO iband2 = 1, hybrid%nbands(nk)
     358             :                   ndb2 = degenerat(iband2)
     359             :                   IF (ndb2 >= 1) THEN
     360             :                      ic2 = ic2 + 1
     361             :                      IF (ndb2 == ndb1) THEN
     362             :                         ! note that this criterium is only valid for pure spatial rotations
     363             :                         ! if one combines spatial rotations with time reversal symmetry there
     364             :                         ! is no unique criteria to identify symequivalent state
     365             :                         ! however, also in the latter case the trace of the spatial rotations
     366             :                         ! for two symmetry equivalent states must be equivalent
     367             :                         IF (all(abs(trace(:sym%nop, ic1) - trace(:sym%nop, ic2)) <= 1E-8))&
     368             :            &            THEN
     369             :                            symequivalent(ic2, ic1) = .true.
     370             :                         END IF
     371             :                      END IF
     372             :                   END IF
     373             :                END DO
     374             :             END IF
     375             :          END DO
     376             : 
     377             :       ELSE
     378             :          ! read in cmt and z at current k-point (nk)
     379             : 
     380           0 :          CALL read_cmt(cmt, nk)
     381             :          !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
     382             : 
     383             :          IF (allocated(olapmt)) deallocate (olapmt)
     384           0 :          ALLOCATE (olapmt(hybrid%maxindx, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype), stat=ok)
     385           0 :          IF (ok /= 0) STOP 'symm: failure allocation olapmt'
     386           0 :          olapmt = 0
     387             : 
     388           0 :          DO itype = 1, atoms%ntype
     389           0 :             DO l = 0, atoms%lmax(itype)
     390           0 :                nn = hybrid%nindx(l, itype)
     391           0 :                DO n2 = 1, nn
     392           0 :                   DO n1 = 1, nn
     393             :                      olapmt(n1, n2, l, itype) = intgrf( &
     394             :           &                        hybdat%bas1(:, n1, l, itype)*hybdat%bas1(:, n2, l, itype)&
     395             :           &                       + hybdat%bas2(:, n1, l, itype)*hybdat%bas2(:, n2, l, itype),&
     396           0 :           &                        atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, hybdat%gridf)
     397             :                   END DO
     398             :                END DO
     399             :             END DO
     400             :          END DO
     401             : 
     402           0 :          ALLOCATE (wavefolap(hybrid%nbands(nk), hybrid%nbands(nk)), carr(hybrid%maxindx), stat=ok)
     403           0 :          IF (ok /= 0) STOP 'symm: failure allocation wfolap/maxindx'
     404           0 :          wavefolap = 0
     405             : 
     406           0 :          iatom = 0
     407           0 :          DO itype = 1, atoms%ntype
     408           0 :             DO ieq = 1, atoms%neq(itype)
     409           0 :                iatom = iatom + 1
     410           0 :                lm = 0
     411           0 :                DO l = 0, atoms%lmax(itype)
     412           0 :                   DO M = -l, l
     413           0 :                      nn = hybrid%nindx(l, itype)
     414           0 :                      DO iband1 = 1, hybrid%nbands(nk)
     415             :                         carr(:nn) = matmul(olapmt(:nn, :nn, l, itype),&
     416           0 :            &                                cmt(iband1, lm + 1:lm + nn, iatom))
     417           0 :                         DO iband2 = 1, iband1
     418             :                            wavefolap(iband2, iband1)&
     419             :             &            = wavefolap(iband2, iband1)&
     420           0 :             &            + dot_product(cmt(iband2, lm + 1:lm + nn, iatom), carr(:nn))
     421             :                         END DO
     422             :                      END DO
     423           0 :                      lm = lm + nn
     424             :                   END DO
     425             :                END DO
     426             :             END DO
     427             :          END DO
     428             : 
     429           0 :          DO iband1 = 1, hybrid%nbands(nk)
     430           0 :             DO iband2 = 1, iband1
     431           0 :                wavefolap(iband1, iband2) = conjg(wavefolap(iband2, iband1))
     432             :             END DO
     433             :          END DO
     434             : 
     435           0 :          ALLOCATE (symequivalent(nddb, nddb), stat=ok)
     436           0 :          IF (ok /= 0) STOP 'symm: failure allocation symequivalent'
     437           0 :          symequivalent = .false.
     438           0 :          ic1 = 0
     439           0 :          DO iband1 = 1, hybrid%nbands(nk)
     440           0 :             ndb1 = degenerat(iband1)
     441           0 :             IF (ndb1 == 0) CYCLE
     442           0 :             ic1 = ic1 + 1
     443           0 :             ic2 = 0
     444           0 :             DO iband2 = 1, hybrid%nbands(nk)
     445           0 :                ndb2 = degenerat(iband2)
     446           0 :                IF (ndb2 == 0) CYCLE
     447           0 :                ic2 = ic2 + 1
     448           0 :                IF (any(abs(wavefolap(iband1:iband1 + ndb1 - 1,&
     449           0 :         &                             iband2:iband2 + ndb2 - 1)) > 1E-9)) THEN
     450             : !     &          .and. ndb1 .eq. ndb2 ) THEN
     451           0 :                   symequivalent(ic2, ic1) = .true.
     452             :                END IF
     453             :             END DO
     454             :          END DO
     455             :       END IF
     456             : 
     457             :       !
     458             :       ! generate index field which contain the band combinations (n1,n2),
     459             :       ! which are non zero
     460             :       !
     461             : 
     462             :       ic1 = 0
     463           0 :       indx_sest = 0
     464           0 :       nsest = 0
     465           0 :       DO iband1 = 1, hybrid%nbands(nk)
     466           0 :          ndb1 = degenerat(iband1)
     467           0 :          IF (ndb1 >= 1) ic1 = ic1 + 1
     468           0 :          i = 0
     469             :          DO WHILE (degenerat(iband1 - i) == 0)
     470           0 :             i = i + 1
     471             :          END DO
     472             :          ndb1 = degenerat(iband1 - i)
     473             :          ic2 = 0
     474           0 :          DO iband2 = 1, hybrid%nbands(nk)
     475           0 :             ndb2 = degenerat(iband2)
     476           0 :             IF (ndb2 >= 1) ic2 = ic2 + 1
     477           0 :             i = 0
     478           0 :             DO WHILE (degenerat(iband2 - i) == 0)
     479             :                i = i + 1
     480             :             END DO
     481           0 :             ndb2 = degenerat(iband2 - i)
     482             :             ! only upper triangular part
     483           0 :             IF (symequivalent(ic2, ic1) .and. iband2 <= iband1) THEN
     484             : !            IF( ndb1 .ne. ndb2 ) STOP 'symm_hf: failure symequivalent'
     485           0 :                nsest(iband1) = nsest(iband1) + 1
     486           0 :                indx_sest(nsest(iband1), iband1) = iband2
     487             :             END IF
     488             :          END DO
     489             :       END DO
     490             : 
     491             :       !
     492             :       ! calculate representations for core states
     493             :       ! (note that for a core state, these are proportional to the Wigner D matrices)
     494             :       !
     495             :       ! Definition of the Wigner rotation matrices
     496             :       !
     497             :       !                     -1                l
     498             :       ! P(R) Y  (r)  = Y  (R  r)  =  sum     D    (R)  Y   (r)
     499             :       !       lm        lm              m'    m m'      lm'
     500             :       !
     501             : 
     502           0 :       pi = pimach()
     503             : 
     504           0 :       IF (hybdat%lmaxcd > atoms%lmaxd) STOP &
     505           0 :      & 'symm_hf: The very impropable case that hybdat%lmaxcd > atoms%lmaxd occurs'
     506             : 
     507           0 :       iatom = 0
     508           0 :       iatom0 = 0
     509           0 :       DO itype = 1, atoms%ntype
     510           0 :          DO ieq = 1, atoms%neq(itype)
     511           0 :             iatom = iatom + 1
     512           0 :             DO iop = 1, nsymop
     513           0 :                isym = psym(iop)
     514           0 :                IF (isym <= sym%nop) THEN
     515             :                   iisym = isym
     516             :                ELSE
     517             :                   iisym = isym - sym%nop
     518             :                END IF
     519             : 
     520           0 :                ratom = hybrid%map(iatom, isym)
     521           0 :                rotkpt = matmul(rrot(:, :, isym), kpts%bkf(:, nk))
     522             :                g = nint(rotkpt - kpts%bkf(:, nk))
     523             : 
     524             :                cdum = exp(-2*pi*img*dot_product(rotkpt, sym%tau(:, iisym)))* &
     525           0 :         &               exp(2*pi*img*dot_product(g, atoms%taual(:, ratom)))
     526             :             END DO
     527             :          END DO
     528           0 :          iatom0 = iatom0 + atoms%neq(itype)
     529             :       END DO
     530             : 
     531           0 :    END SUBROUTINE symm_hf
     532             : 
     533           0 :    INTEGER FUNCTION symm_hf_nkpt_EIBZ(kpts, nk, sym)
     534             : 
     535             :       USE m_util, ONLY: modulo1
     536             :       USE m_types
     537             :       IMPLICIT NONE
     538             :       TYPE(t_sym), INTENT(IN)   :: sym
     539             :       TYPE(t_kpts), INTENT(IN)   :: kpts
     540             : 
     541             : !     - scalar input -
     542             :       INTEGER, INTENT(IN)   :: nk
     543             : !     - array input -
     544             : 
     545             : !     - local scalars -
     546             :       INTEGER               ::  isym, ic, iop, ikpt, ikpt1
     547             :       INTEGER               ::  nsymop, nrkpt
     548             : !     - local arrays -
     549           0 :       INTEGER               ::  rrot(3, 3, sym%nsym)
     550           0 :       INTEGER               ::  neqvkpt(kpts%nkptf), list(kpts%nkptf), parent(kpts%nkptf),&
     551           0 :      &                          symop(kpts%nkptf)
     552           0 :       INTEGER, ALLOCATABLE  ::  psym(:)!,help(:)
     553             :       REAL                  ::  rotkpt(3)
     554             : 
     555             :       ! calculate rotations in reciprocal space
     556           0 :       DO isym = 1, sym%nsym
     557           0 :          IF (isym <= sym%nop) THEN
     558           0 :             rrot(:, :, isym) = transpose(sym%mrot(:, :, sym%invtab(isym)))
     559             :          ELSE
     560           0 :             rrot(:, :, isym) = -rrot(:, :, isym - sym%nop)
     561             :          END IF
     562             :       END DO
     563             : 
     564             :       ! determine little group of k., i.e. those symmetry operations
     565             :       ! which keep bk(:,nk,nw) invariant
     566             :       ! nsymop :: number of such symmetry-operations
     567             :       ! psym   :: points to the symmetry-operation
     568             : 
     569           0 :       ic = 0
     570           0 :       ALLOCATE (psym(sym%nsym))
     571             : 
     572           0 :       DO iop = 1, sym%nsym
     573           0 :          rotkpt = matmul(rrot(:, :, iop), kpts%bkf(:, nk))
     574             : 
     575             :          !transfer rotkpt into BZ
     576           0 :          rotkpt = modulo1(rotkpt, kpts%nkpt3)
     577             : 
     578             :          !check if rotkpt is identical to bk(:,nk)
     579           0 :          IF (maxval(abs(rotkpt - kpts%bkf(:, nk))) <= 1E-07) THEN
     580           0 :             ic = ic + 1
     581           0 :             psym(ic) = iop
     582             :          END IF
     583             :       END DO
     584             :       nsymop = ic
     585             : 
     586             :       ! reallocate psym
     587             : !       ALLOCATE(help(ic))
     588             : !       help = psym(1:ic)
     589             : !       DEALLOCATE(psym)
     590             : !       ALLOCATE(psym(ic))
     591             : !       psym = help
     592             : !       DEALLOCATE(help)
     593             : 
     594             :       ! determine extented irreducible BZ of k ( EIBZ(k) ), i.e.
     595             :       ! those k-points, which can generate the whole BZ by
     596             :       ! applying the symmetry operations of the little group of k
     597             : 
     598           0 :       neqvkpt = 0
     599             : 
     600             : !       list = (/ (ikpt-1, ikpt=1,nkpt) /)
     601           0 :       DO ikpt = 1, kpts%nkptf
     602           0 :          list(ikpt) = ikpt - 1
     603             :       END DO
     604             : 
     605           0 :       DO ikpt = 2, kpts%nkptf
     606           0 :          DO iop = 1, nsymop
     607             : 
     608           0 :             rotkpt = matmul(rrot(:, :, psym(iop)), kpts%bkf(:, ikpt))
     609             : 
     610             :             !transfer rotkpt into BZ
     611           0 :             rotkpt = modulo1(rotkpt, kpts%nkpt3)
     612             : 
     613             :             !determine number of rotkpt
     614           0 :             nrkpt = 0
     615           0 :             DO ikpt1 = 1, kpts%nkptf
     616           0 :                IF (maxval(abs(rotkpt - kpts%bkf(:, ikpt1))) <= 1E-06) THEN
     617             :                   nrkpt = ikpt1
     618             :                   EXIT
     619             :                END IF
     620             :             END DO
     621           0 :             IF (nrkpt == 0) STOP 'symm: Difference vector not found !'
     622             : 
     623           0 :             IF (list(nrkpt) /= 0) THEN
     624           0 :                list(nrkpt) = 0
     625           0 :                neqvkpt(ikpt) = neqvkpt(ikpt) + 1
     626           0 :                parent(nrkpt) = ikpt
     627           0 :                symop(nrkpt) = psym(iop)
     628             :             END IF
     629           0 :             IF (all(list == 0)) EXIT
     630             : 
     631             :          END DO
     632             :       END DO
     633             : 
     634             :       ! for the Gamma-point holds:
     635           0 :       parent(1) = 1
     636           0 :       neqvkpt(1) = 1
     637             : 
     638             :       ! determine number of members in the EIBZ(k)
     639           0 :       ic = 0
     640           0 :       DO ikpt = 1, kpts%nkptf
     641           0 :          IF (parent(ikpt) == ikpt) ic = ic + 1
     642             :       END DO
     643           0 :       symm_hf_nkpt_EIBZ = ic
     644             : 
     645           0 :    END FUNCTION symm_hf_nkpt_EIBZ
     646             : 
     647             : END MODULE m_symm_hf

Generated by: LCOV version 1.13