!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_grp_k
  USE m_juDFT

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: grp_k, euler

CONTAINS 

  SUBROUTINE grp_k(sym,mrot_k,cell,bk,nclass,nirr,char_table,&
                   grpname,irrname,su,sp_alph,sp_beta)

    !************************************************************
    !
    !   Determines the group of k, returns the number of classes, the name of the
    !   group, the irreducible representations and the character table.
    !   All the groups are not implemented yet, and the identification is not tested for 
    !   the groups. 
    !
    !  The subroutine works also with double groups, however, no double groups are 
    !  tabulated at the moment.
    !
    ! 
    !   Character tables are taken from 
    !   T. Inui, Y Tanabe, and Y. Onodera, "Group theory and its applications in physics",
    !   Springer (1996)
    !
    !   Jussi Enkovaara, Juelich 2004
    !**************************************************************

    !      USE m_mrot2su
    USE m_inv3
    USE m_constants
    USE m_socsym,    ONLY : soc_sym, cross
    USE m_types
    IMPLICIT NONE

    TYPE(t_sym),INTENT(IN)   :: sym
    TYPE(t_cell),INTENT(IN)   :: cell
   

    REAL,INTENT(IN)      :: bk(3)
    COMPLEX, INTENT(OUT) :: char_table(:,:)
    INTEGER, INTENT(OUT) :: mrot_k(:,:,:),nclass,nirr
    CHARACTER(LEN=7),  INTENT(OUT) :: grpname
    CHARACTER(LEN=5),  INTENT(OUT) :: irrname(:)
    COMPLEX, OPTIONAL, INTENT(OUT) :: su(:,:,:)
    REAL,    OPTIONAL, INTENT(IN)  :: sp_alph,sp_beta ! spin quant. axis

    ! locals
    INTEGER ::  nopk,c
    LOGICAL, ALLOCATABLE :: error(:)
    INTEGER :: mrot2(3,3,48)
    REAL    :: ktest(3),rtmp,rtmp2(4),kt(3)
    INTEGER :: mtest(3,3),mtmp(3,3),mtmpinv(3,3),munit(3,3)
    INTEGER :: n,n2 ,i,itmp
    INTEGER :: elem(48),belongs(48),members(48)
    LOGICAL :: soc,l_sorted
    ! COMPLEX :: sutmp(2,2), sutmpinv(2,2), su2(2,2,48)
    COMPLEX :: sutest(2,2)
    INTEGER :: d,det(48),rot(12),rot_type(48)
    REAL :: alpha,beta,gamma,theta(48),rax(4,48)
    REAL,PARAMETER :: eps=0.000001
    COMPLEX :: omega
    COMPLEX,PARAMETER :: one=CMPLX(1.0,0.0),zero=CMPLX(0.0,0.0)
    COMPLEX,PARAMETER :: imi=CMPLX(0.0,1.0)


    soc=.FALSE.
    IF (PRESENT(su)) soc=.TRUE.
    !sym%nop=SIZE(sym%mrot,3)

    ALLOCATE(error(sym%nop))
    error=.FALSE.
    ! Reduce the symmetry due to spin-orbit
    IF (soc.AND.PRESENT(sp_alph)) THEN
       CALL soc_sym(sym%nop,sym%mrot,sp_beta,sp_alph,cell%amat,error)!keep
    ENDIF

    ! determine the group of k
    nopk=0
    ksymloop: DO n=1,sym%nop
       ktest(1)=bk(1)-sym%mrot(1,1,n)*bk(1)-sym%mrot(2,1,n)*bk(2)-sym%mrot(3,1,n)*bk(3)
       ktest(2)=bk(2)-sym%mrot(1,2,n)*bk(1)-sym%mrot(2,2,n)*bk(2)-sym%mrot(3,2,n)*bk(3)
       ktest(3)=bk(3)-sym%mrot(1,3,n)*bk(1)-sym%mrot(2,3,n)*bk(2)-sym%mrot(3,3,n)*bk(3)
       IF (( ABS( ktest(1) - NINT(ktest(1)) ) < eps ) .AND.&
            &         ( ABS( ktest(2) - NINT(ktest(2)) ) < eps ) .AND.&
            &         ( ABS( ktest(3) - NINT(ktest(3)) ) < eps ) .AND.&
            &         (.NOT.error(n)) ) THEN
          nopk=nopk+1
          mrot_k(:,:,nopk)=sym%mrot(:,:,n)
          CYCLE ksymloop
       ENDIF
    ENDDO ksymloop

    DEALLOCATE(error)

    ! Determine the spin-rotations
    ! Double groups not used at the moment, the groups are classified 
    ! without the spin rotations
    !       IF (soc) CALL mrot2su(mrot_k(:,:,1:nopk),amat,su)

    ! identify the group
    ! first determine the classes
    members=1
    nclass=1
    elem(nclass)=1
    belongs(1)=1
    !       IF (soc) THEN 
    ! double group
    !         DO n=1,nopk        
    !            mrot_k(:,:,n+nopk)=mrot_k(:,:,n)
    !            su(:,:,n+nopk)=-su(:,:,n)
    !         ENDDO
    !         nopk=2*nopk
    !       ENDIF

    cloop: DO n=2,nopk
       DO n2=1,nopk
          mtmp=mrot_k(:,:,n2)
          CALL inv3(mrot_k(:,:,n2),mtmpinv,d)
          mtest=MATMUL(MATMUL(mtmp,mrot_k(:,:,n)),mtmpinv)
          !            IF (soc) THEN
          !               sutmp=su(:,:,n2)
          !               sutmpinv=CONJG(TRANSPOSE(sutmp))
          !               sutest=MATMUL(MATMUL(sutmp,su(:,:,n)),sutmpinv)
          !            ENDIF
          DO c=1,nclass
             !               IF (soc) THEN
             !                  IF (ALL((mtest-mrot_k(:,:,elem(c))).EQ.0).AND.
             !c     &   ALL(ABS(REAL(sutest-su(:,:,elem(c)))).LE.0.0001).AND.
             !     &   ALL(ABS(AIMAG(sutest-su(:,:,elem(c)))).LE.0.0001)) THEN
             !                     belongs(n)=c
             !                     CYCLE classloop
             !                  ENDIF
             !               ELSE
             IF (ALL((mtest-mrot_k(:,:,elem(c))).EQ.0)) THEN 
                belongs(n)=c
                members(c)=members(c)+1
                CYCLE cloop
             ENDIF
             !               ENDIF
          ENDDO
       ENDDO
       nclass=nclass+1
       elem(nclass)=n
       belongs(n)=nclass
    ENDDO cloop
    nirr=nclass
    IF (soc) nirr=2*nclass

    ! sort the classes according the rotation type
    mrot2(:,:,1:nopk)=mrot_k(:,:,1:nopk)
    DO c=1,nclass
       mrot_k(:,:,c)=mrot2(:,:,elem(c))
    ENDDO
    !        IF (soc) THEN
    !           su2(:,:,1:nopk)=su(:,:,1:nopk)
    !           DO c=1,nclass
    !              su(:,:,c)=su2(:,:,elem(c))
    !           ENDDO
    !        ENDIF



    ! identify the group
    rot=0
    rot_type=0
    ! rot_type: 2=two fold rot, 3=three fold rot, etc, 7=identity
    !           8=two fold improper rot, ...
    munit=0
    munit(1,1)=1
    munit(2,2)=1
    munit(3,3)=1
    ! determine the number of different rotations
    DO c=1,nclass
       det(c)=&
            &      mrot_k(1,1,c)*mrot_k(2,2,c)*mrot_k(3,3,c)+&
            &      mrot_k(1,2,c)*mrot_k(2,3,c)*mrot_k(3,1,c)+&
            &      mrot_k(2,1,c)*mrot_k(3,2,c)*mrot_k(1,3,c)-&
            &      mrot_k(1,3,c)*mrot_k(2,2,c)*mrot_k(3,1,c)-&
            &      mrot_k(2,3,c)*mrot_k(3,2,c)*mrot_k(1,1,c)-&
            &      mrot_k(2,1,c)*mrot_k(1,2,c)*mrot_k(3,3,c)
       mtest=det(c)*mrot_k(:,:,c)
       rotloop: DO i=1,6
          IF (ALL((mtest-munit).EQ.0)) THEN
             rot(i+(1-det(c))*3)=rot(i+(1-det(c))*3)+1
             IF (i.EQ.1) THEN
                rot_type(c)=7-(1-det(c))*3
             ELSE
                rot_type(c)=i-(1-det(c))*3
             ENDIF
             EXIT rotloop
          ENDIF
          mtest=MATMUL(det(c)*mrot_k(:,:,c),mtest)
       ENDDO rotloop
       IF (ANY((mtest-munit).NE.0)) THEN
          WRITE(oUnit,*) 'grp_k: Cannot find the order of rotation'
       ENDIF
       IF ((rot(5).GT.0).OR.(rot(11).GT.0)) THEN
          WRITE(oUnit,*) 'grp_k: 5 fold rotation found!'
       ENDIF
       CALL euler(c,sym,cell,alpha,beta,gamma)
       CALL rotaxis(alpha,beta,gamma,rax(1:4,c),theta(c))
       IF (soc) THEN
          su(1,1,c)=COS(beta/2.0)*EXP(CMPLX(0.0,-(alpha+gamma)/2.0))
          su(1,2,c)=-SIN(beta/2.0)*EXP(CMPLX(0.0,-(alpha-gamma)/2.0))
          su(2,1,c)=SIN(beta/2.0)*EXP(CMPLX(0.0,(alpha-gamma)/2.0))
          su(2,2,c)=COS(beta/2.0)*EXP(CMPLX(0.0,(alpha+gamma)/2.0))
       ENDIF
    ENDDO

    !<-- Sort the classes
    ! The group elements are sorted in the following way:
    ! First the proper rotations, then improper ones, with increasing rotation angle
    ! Rotations with the same angle are arranged with increasing magnitude of the rotation
    ! axis, e.g. (I, 90 deg around 001, 90 around 110, mirror, ...

    l_sorted=.FALSE.
    DO WHILE (.NOT.l_sorted)
       l_sorted=.TRUE.
       DO c=1,nclass-1
          IF (rot_type(c).LT.rot_type(c+1)) THEN
             mtest=mrot_k(:,:,c)
             mrot_k(:,:,c)=mrot_k(:,:,c+1)
             mrot_k(:,:,c+1)=mtest
             rtmp=theta(c)
             theta(c)=theta(c+1)
             theta(c+1)=rtmp
             rtmp2=rax(:,c)
             rax(:,c)=rax(:,c+1)
             rax(:,c+1)=rtmp2
             itmp=det(c)
             det(c)=det(c+1)
             det(c+1)=itmp
             itmp=rot_type(c)
             rot_type(c)=rot_type(c+1)
             rot_type(c+1)=itmp
             itmp=members(c)
             members(c)=members(c+1)
             members(c+1)=itmp
             IF (soc) THEN
                sutest=su(:,:,c)
                su(:,:,c)=su(:,:,c+1)
                su(:,:,c+1)=sutest
             ENDIF
             l_sorted=.FALSE.
          ENDIF
          IF ((rot_type(c).EQ.rot_type(c+1)).AND.&
               &            (theta(c).GT.theta(c+1))) THEN
             !              IF (theta(c)+(1-det(c))*2.0*pi.GT.
             !     &             theta(c+1)+(1-det(c+1))*2.0*pi) THEN
             mtest=mrot_k(:,:,c)
             mrot_k(:,:,c)=mrot_k(:,:,c+1)
             mrot_k(:,:,c+1)=mtest
             rtmp=theta(c)
             theta(c)=theta(c+1)
             theta(c+1)=rtmp
             rtmp2=rax(:,c)
             rax(:,c)=rax(:,c+1)
             rax(:,c+1)=rtmp2
             itmp=det(c)
             det(c)=det(c+1)
             det(c+1)=itmp
             itmp=members(c)
             members(c)=members(c+1)
             members(c+1)=itmp
             IF (soc) THEN
                sutest=su(:,:,c)
                su(:,:,c)=su(:,:,c+1)
                su(:,:,c+1)=sutest
             ENDIF
             l_sorted=.FALSE.
          ENDIF
          IF ((rot_type(c).EQ.rot_type(c+1)).AND.&
               &            (ABS(theta(c)-theta(c+1)).LT.0.0001).AND.    &
               &            (rax(4,c).GT.rax(4,c+1))) THEN
             !              IF ((ABS(theta(c)+(1-det(c))*2.0*pi-
             !     &           theta(c+1)-(1-det(c+1))*2.0*pi).LT.0.0001).AND.
             !     &           (rax(4,c).GT.rax(4,c+1))) THEN
             mtest=mrot_k(:,:,c)
             mrot_k(:,:,c)=mrot_k(:,:,c+1)
             mrot_k(:,:,c+1)=mtest
             rtmp2=rax(:,c)
             rax(:,c)=rax(:,c+1)
             rax(:,c+1)=rtmp2
             itmp=members(c)
             members(c)=members(c+1)
             members(c+1)=itmp
             IF (soc) THEN
                sutest=su(:,:,c)
                su(:,:,c)=su(:,:,c+1)
                su(:,:,c+1)=sutest
             ENDIF
             l_sorted=.FALSE.
          ENDIF
       ENDDO
    ENDDO
    !>


    WRITE(24,110) bk
    DO c=1,nclass
       IF (det(c).EQ.1) THEN
          WRITE(24,111) NINT(theta(c)*180/pi_const),(rax(1:3,c)),members(c)
       ELSE
          WRITE(24,112) NINT(theta(c)*180/pi_const),(rax(1:3,c)),members(c)
       ENDIF
    ENDDO
110 FORMAT('Symmetry operations of group of k=',3f6.3)
111 FORMAT(i3,1x,'degree proper rotation around ',3f6.2,3x,&
         &      i3,' members in class')
112 FORMAT(i3,1x,'degree improper rotation around ',3f6.2,3x,i3&
         &       ,' members in class')


    !<-- Character tables

    char_table=0.0
    char_table(1,:)=1.0
    grpname='Unknown'
    irrname='Unkno'   ! Only 5 characters wide, cf. ../sympsi.F.
    ! First look the number of classes, within groups with the same number of classes
    ! check the number of different rotations
    SELECT CASE(nclass)
    CASE(1)                           ! only C1
       grpname='C1'
       char_table(1,1)=1.0
       irrname(1)='Gam1'
       IF (soc) THEN
          char_table(2,1)=1.0
          irrname(2)='Gam2'
       ENDIF
    CASE(2)                           ! C-1, C2 and Cs
       IF (rot(2).GT.0) THEN
          grpname='C2'
          char_table(1,1:2)=(/1.0,  1.0/)
          char_table(2,1:2)=(/1.0, -1.0/)
          irrname(1)='Gam1'
          irrname(2)='Gam2'
       ELSE
          grpname='C1h'
          char_table(1,1:2)=(/1.0,  1.0/)
          char_table(2,1:2)=(/1.0, -1.0/)
          irrname(1)='Gam1'
          irrname(2)='Gam2'
          IF (soc) THEN
             char_table(3,1:2)=(/one, -imi/)
             char_table(4,1:2)=(/one,  imi/)
             irrname(3)='Gam3'
             irrname(4)='Gam4'
          ENDIF
       ENDIF
    CASE(3)                           ! C3, D3 and C3v
       IF (ANY(det(1:3).EQ.-1)) THEN
          grpname='C3v'
          char_table(1,1:3)=(/1.0,  1.0,  1.0/)
          char_table(2,1:3)=(/1.0,  1.0, -1.0/)
          char_table(3,1:3)=(/2.0, -1.0,  0.0/)
          irrname(1)='Lam1'
          irrname(2)='Lam2'
          irrname(3)='Lam3'
          IF (soc) THEN
             nirr=6
             char_table(4,1:3)=(/2.0,  1.0,  0.0/)
             char_table(5,1:3)=(/one, -one,  imi/)
             char_table(6,1:3)=(/one, -one, -imi/)
             irrname(4)='Lam6'
             irrname(5)='Lam4'
             irrname(6)='Lam5'
          ENDIF
       ELSE
          grpname='C3'
          omega=EXP(CMPLX(0.0,-2*pi_const/3.0))
          char_table(1,1:3)=(/1.0,  1.0,  1.0/)
          char_table(2,1:3)=(/one, omega, omega**2/)
          char_table(3,1:3)=(/one, omega**2, omega /)
          irrname(1)='Gam1'
          irrname(2)='Gam2'
          irrname(3)='Gam3'
       ENDIF
    CASE(4)                           ! C2h, D2, C2v, C4, S4, and T
       IF((rot(2).EQ.1).AND.(rot(8).EQ.2)) THEN
          grpname='C2v'
          char_table(1,1:4)=(/1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:4)=(/1.0,  1.0, -1.0, -1.0/)
          char_table(3,1:4)=(/1.0, -1.0, -1.0,  1.0/)
          char_table(4,1:4)=(/1.0, -1.0,  1.0, -1.0/)
          irrname(1)='Z1'
          irrname(2)='Z2'
          irrname(3)='Z3'
          irrname(4)='Z4'            
       ELSE IF (rot(3).GT.1) THEN
          grpname='T'
          omega=EXP(CMPLX(0.0,-2*pi_const/3.0))
          char_table(1,1:4)=(/1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:4)=(/one,  one, omega, omega**2/)
          char_table(3,1:4)=(/one,  one, omega**2,  omega/)
          char_table(4,1:4)=(/3.0, -1.0,  0.0,  0.0/)
          irrname(1)='Gam1'
          irrname(2)='Gam2'
          irrname(3)='Gam3'
          irrname(4)='Gam4'
       ELSE IF ((rot(4).GT.1).OR.(rot(10).GT.1)) THEN
          IF (rot(4).GT.1) grpname='C4'
          IF (rot(10).GT.1) grpname='S4'
          char_table(1,1:4)=(/1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:4)=(/1.0, -1.0,  1.0, -1.0/)
          char_table(3,1:4)=(/one, -imi, -one,  imi/)
          char_table(4,1:4)=(/one,  imi, -one, -imi/)
          irrname(1)='Gam1'
          irrname(2)='Gam2'
          irrname(3)='Gam3'
          irrname(4)='Gam4'
          IF (soc) THEN
             nirr=8
             omega=EXP(CMPLX(0.0,-pi_const/4.0))
             char_table(5,1:4)=(/one, omega, -imi, &
                  &                   -CONJG(omega)/)
             char_table(6,1:4)=(/one, CONJG(omega), imi, &
                  &                   -omega/)
             char_table(7,1:4)=(/one, -omega, -imi, &
                  &                   CONJG(omega)/)
             char_table(8,1:4)=(/one, -CONJG(omega), imi, &
                  &                   omega/)
             irrname(5)='Gam5'
             irrname(6)='Gam6'
             irrname(7)='Gam7'
             irrname(8)='Gam8'                   
          ENDIF
       ELSE IF ((rot(2).EQ.1).AND.(rot(8).EQ.1)) THEN
          grpname='ZKUS'
       ENDIF
    CASE(5)                           ! D4, C4v, D2d, O and Td
       IF (rot(3).GT.0) THEN
          grpname='Td'
          char_table(1,1:5)=(/1.0,  1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:5)=(/1.0, -1.0,  1.0, -1.0,  1.0/)
          char_table(3,1:5)=(/2.0,  0.0,  2.0,  0.0, -1.0/)
          char_table(4,1:5)=(/3.0,  1.0, -1.0, -1.0,  0.0/)
          char_table(5,1:5)=(/3.0, -1.0, -1.0,  1.0,  0.0/)
          irrname(1)='P1'
          irrname(2)='P2'
          irrname(3)='P3'
          irrname(4)='P5'
          irrname(5)='P4'
       ELSE IF (rot(10).GT.0) THEN
          grpname='D2d'
          char_table(1,1:5)=(/1.0,  1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:5)=(/1.0,  1.0, -1.0,  1.0, -1.0/)
          char_table(3,1:5)=(/1.0,  1.0,  1.0, -1.0, -1.0/)
          char_table(4,1:5)=(/1.0,  1.0, -1.0, -1.0,  1.0/)
          char_table(5,1:5)=(/2.0, -2.0,  0.0,  0.0,  0.0/)
          !               standard order: E, 2, 2_x, -4, m
          irrname(1)='W1'
          irrname(2)='W2'
          irrname(3)='W1`'
          irrname(4)='W2`'
          irrname(5)='W3'
       ELSE
          grpname='C4v'
          char_table(1,1:5)=(/1.0,  1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:5)=(/1.0,  1.0,  1.0, -1.0, -1.0/)
          char_table(3,1:5)=(/1.0, -1.0,  1.0,  1.0, -1.0/)
          char_table(4,1:5)=(/1.0, -1.0,  1.0, -1.0,  1.0/)
          char_table(5,1:5)=(/2.0,  0.0, -2.0,  0.0,  0.0/)
          irrname(1)='Del1'
          irrname(2)='Del1`'
          irrname(3)='Del2'
          irrname(4)='Del2`'
          irrname(5)='Del5'
       ENDIF
    CASE(6)                           ! C3i, D3d, C6, C3h, D6, C6v and D3h
       IF (rot(2).EQ.0) THEN
          IF (rot(7).EQ.0) THEN
             grpname='C3h'
             char_table(:,:) = 0.0
             WRITE(24,*) 'C3h: Character table missing in grp_k.F'
          ELSE
             grpname='C3i'
             omega=CMPLX(-0.5,SQRT(3./4.))
             char_table(1,1:6)=(/one, one, one, one, one, one/)
             char_table(2,1:6)=(/one, one, one,-one,-one,-one/)
             char_table(3,1:6)=(/one, omega, CONJG(omega),&
                  &                                one, omega, CONJG(omega)/)
             char_table(4,1:6)=(/one, omega, CONJG(omega),&
                  &                               -one,-omega,-CONJG(omega)/)
             char_table(5,1:6)=(/one, CONJG(omega), omega,&
                  &                                one, CONJG(omega), omega/)
             char_table(6,1:6)=(/one, CONJG(omega), omega,&
                  &                               -one,-CONJG(omega),-omega/)
             irrname(1)='Ag  '
             irrname(2)='Au  '
             irrname(3)='E1g '
             irrname(4)='E1u '
             irrname(5)='E2g '
             irrname(6)='E2u '
             IF (soc) THEN
                char_table(7,1:6)=(/one,-one, one, one,-one, one/)
                char_table(8,1:6)=(/one,-one, one,-one, one,-one/)
                char_table( 9,1:6)=(/one,-omega, CONJG(omega),&
                     &                                one,-omega, CONJG(omega)/)
                char_table(10,1:6)=(/one,-omega, CONJG(omega),&
                     &                               -one, omega,-CONJG(omega)/)
                char_table(11,1:6)=(/one,-CONJG(omega), omega,&
                     &                                one,-CONJG(omega), omega/)
                char_table(12,1:6)=(/one,-CONJG(omega), omega,&
                     &                               -one, CONJG(omega),-omega/)
                irrname(7)="Ag' "
                irrname(8)="Au' "
                irrname(9)="E1g'"
                irrname(10)="E1u'"
                irrname(11)="E2g'"
                irrname(12)="E2u'"
             ENDIF
          ENDIF
       ELSEIF (rot(2).GT.1) THEN
          grpname='D6'
          char_table(1,1:6)=(/1.0,  1.0,  1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:6)=(/1.0,  1.0,  1.0,  1.0, -1.0, -1.0/)
          char_table(3,1:6)=(/1.0, -1.0,  1.0, -1.0,  1.0, -1.0/)
          char_table(4,1:6)=(/1.0, -1.0,  1.0, -1.0, -1.0,  1.0/)
          char_table(5,1:6)=(/2.0,  1.0, -1.0, -2.0,  0.0,  0.0/)
          char_table(6,1:6)=(/2.0, -1.0, -1.0,  2.0,  0.0,  0.0/)
          irrname(1)='Gam1'
          irrname(2)='Gam2'
          irrname(3)='Gam3'
          irrname(4)='Gam4'
          irrname(5)='Gam6'
          irrname(6)='Gam5'
       ELSE
          IF (rot(6).EQ.0) THEN
             grpname='D3d'
             char_table(1,1:6)=(/1.0,  1.0,  1.0,  1.0,  1.0,  1.0/)
             char_table(2,1:6)=(/1.0,  1.0, -1.0,  1.0,  1.0, -1.0/)
             char_table(3,1:6)=(/2.0, -1.0,  0.0,  2.0, -1.0,  0.0/)
             char_table(4,1:6)=(/1.0,  1.0,  1.0, -1.0, -1.0, -1.0/)
             char_table(5,1:6)=(/1.0,  1.0, -1.0, -1.0, -1.0,  1.0/)
             char_table(6,1:6)=(/2.0, -1.0,  0.0, -2.0,  1.0,  0.0/)
             irrname(1)='A1g'
             irrname(2)='A2g'
             irrname(3)='Eg'
             irrname(4)='A1u'
             irrname(5)='A2u'
             irrname(6)='Eu'
          ELSEIF (rot(6).GT.1) THEN
             grpname='C6' 
             char_table(:,:) = 0.0
             WRITE(24,*) 'C6: Character table missing in grp_k.F'
          ELSE
             IF (members(3).GT.1) grpname='D3h'    ! maybe this works
             IF (members(3).EQ.1) grpname='C6v'
             char_table(1,1:6)=(/1.0,  1.0,  1.0,  1.0,  1.0,  1.0/)
             char_table(2,1:6)=(/1.0,  1.0,  1.0,  1.0, -1.0, -1.0/)
             char_table(3,1:6)=(/1.0, -1.0,  1.0, -1.0,  1.0, -1.0/)
             char_table(4,1:6)=(/1.0, -1.0,  1.0, -1.0, -1.0,  1.0/)
             char_table(5,1:6)=(/2.0,  1.0, -1.0, -2.0,  0.0,  0.0/)
             char_table(6,1:6)=(/2.0, -1.0, -1.0,  2.0,  0.0,  0.0/)
             irrname(1)='Gam1'
             irrname(2)='Gam2'
             irrname(3)='Gam3'
             irrname(4)='Gam4'
             irrname(5)='Gam6'
             irrname(6)='Gam5'
          ENDIF
       ENDIF
    CASE(8)                           ! Th, C4h and D2h
       IF (rot(3).GT.0) THEN
          grpname='Th'
          char_table(:,:) = 0.0
          WRITE(24,*) 'Th: Character table missing in grp_k.F'
       ELSE IF (rot(2).EQ.3) THEN
          grpname='D2h'
          char_table(:,:) = 0.0
          WRITE(24,*) 'D2h: Character table missing in grp_k.F'
       ELSE
          grpname='C4h'
          char_table(:,:) = 0.0
          WRITE(24,*) 'C4h: Character table missing in grp_k.F'
       ENDIF
    CASE(10)                           ! D4h and Oh
       IF (rot(3).GT.0) THEN
          grpname='Oh'
          char_table(1,1:10)= (/1.0,  1.0,  1.0,  1.0,  1.0,&
               &                                1.0,  1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:10)= (/1.0, -1.0,  1.0,  1.0, -1.0,&
               &                                1.0, -1.0,  1.0,  1.0, -1.0/)
          char_table(3,1:10)= (/2.0,  0.0, -1.0,  2.0,  0.0,&
               &                                2.0,  0.0, -1.0,  2.0,  0.0/)
          char_table(4,1:10)= (/3.0,  1.0,  0.0, -1.0, -1.0,&
               &                                3.0,  1.0,  0.0, -1.0, -1.0/)
          char_table(5,1:10)= (/3.0, -1.0,  0.0, -1.0,  1.0,&
               &                                3.0, -1.0,  0.0, -1.0,  1.0/)
          char_table(6,1:10)= (/1.0,  1.0,  1.0,  1.0,  1.0,&
               &                               -1.0, -1.0, -1.0, -1.0, -1.0/)
          char_table(7,1:10)= (/1.0, -1.0,  1.0,  1.0, -1.0,&
               &                               -1.0,  1.0, -1.0, -1.0,  1.0/)
          char_table(8,1:10)= (/2.0,  0.0, -1.0,  2.0,  0.0,&
               &                               -2.0,  0.0,  1.0, -2.0,  0.0/)
          char_table(9,1:10)= (/3.0,  1.0,  0.0, -1.0, -1.0,&
               &                               -3.0, -1.0,  0.0,  1.0,  1.0/)
          char_table(10,1:10)=(/3.0, -1.0,  0.0, -1.0,  1.0,&
               &                               -3.0,  1.0,  0.0,  1.0, -1.0/)

          irrname(1)='Gam1+'
          irrname(2)='Gam2+'
          irrname(3)='Gam3+'
          irrname(4)='Gam4+'
          irrname(5)='Gam5+'
          irrname(6)='Gam1-'
          irrname(7)='Gam2-'
          irrname(8)='Gam3-'
          irrname(9)='Gam4-'
          irrname(10)='Gam5-'
       ELSE
          grpname='D4h'
          char_table(1,1:10)= (/1.0,  1.0,  1.0,  1.0,  1.0,&
               &                                1.0,  1.0,  1.0,  1.0,  1.0/)
          char_table(2,1:10)= (/1.0,  1.0,  1.0, -1.0, -1.0,&
               &                                1.0,  1.0,  1.0, -1.0, -1.0/)
          char_table(3,1:10)= (/1.0, -1.0,  1.0,  1.0, -1.0,&
               &                                1.0, -1.0,  1.0,  1.0, -1.0/)
          char_table(4,1:10)= (/1.0, -1.0,  1.0, -1.0,  1.0,&
               &                                1.0, -1.0,  1.0, -1.0,  1.0/)
          char_table(5,1:10)= (/2.0,  0.0, -2.0,  0.0,  0.0,&
               &                                2.0,  0.0, -2.0,  0.0,  0.0/)
          char_table(6,1:10)= (/1.0,  1.0,  1.0,  1.0,  1.0,&
               &                               -1.0, -1.0, -1.0, -1.0, -1.0/)
          char_table(7,1:10)= (/1.0,  1.0,  1.0, -1.0, -1.0,&
               &                               -1.0, -1.0, -1.0,  1.0,  1.0/)
          char_table(8,1:10)= (/1.0, -1.0,  1.0,  1.0, -1.0,&
               &                               -1.0,  1.0, -1.0, -1.0,  1.0/)
          char_table(9,1:10)= (/1.0, -1.0,  1.0, -1.0,  1.0,&
               &                               -1.0,  1.0, -1.0,  1.0, -1.0/)
          char_table(10,1:10)=(/2.0,  0.0, -2.0,  0.0,  0.0,&
               &                               -2.0,  0.0,  2.0,  0.0,  0.0/)

          ! standard char-table assumes the order: E, 4, 2, 2_h, 2_h' where 4 and 2 
          ! rotate around the same axis and 2_h is perpendicular to this axis. Check:

          CALL cross( rax(1,2),rax(1,3),kt )
          IF (kt(1)**2+kt(2)**2+kt(3)**2 < eps) THEN
             !                  all ok
          ELSE 
             CALL cross( rax(1,2),rax(1,4),kt )
             IF (kt(1)**2+kt(2)**2+kt(3)**2 < eps) THEN      ! change 3 & 4
                CALL change_column(char_table,10,10,3,4)
                ! check, whether element 3 is perpendicular to element 2
                IF (DOT_PRODUCT(rax(1:3,2),rax(1:3,3)) < eps) THEN
                   !                     all ok
                ELSE             ! change 4 & 5
                   WRITE(*,*) DOT_PRODUCT(rax(1:3,2),rax(1:3,3))
                   CALL change_column(char_table,10,10,4,5)
                ENDIF
             ELSE
                CALL cross( rax(1,2),rax(1,5),kt )
                IF (kt(1)**2+kt(2)**2+kt(3)**2 < eps) THEN    ! change 3 & 5
                   CALL change_column(char_table,10,10,3,5)
                   IF (DOT_PRODUCT(rax(1:3,2),rax(1:3,3)) < eps) THEN ! change 4 & 5
                      CALL change_column(char_table,10,10,4,5)
                   ELSE             
                      !                     all ok
                   ENDIF
                ELSE
                   CALL juDFT_error("D4h",calledby="grp_k")
                ENDIF
             ENDIF
          ENDIF
          irrname(1)='A1g'
          irrname(2)='A2g'
          irrname(3)='B1g'
          irrname(4)='B2g'
          irrname(5)='Eg'
          irrname(6)='A1u'
          irrname(7)='A2u'
          irrname(8)='B1u'
          irrname(9)='B2u'
          irrname(10)='Eu'
       ENDIF
    CASE(12)                           ! D6h and C6h
       IF (rot(2).EQ.3) THEN
          grpname='D6h'
          char_table(:,:) = 0.0
          WRITE(24,*) 'D6h: Character table missing in grp_k.F'
       ELSE
          grpname='C6h'
          char_table(:,:) = 0.0
          WRITE(24,*) 'C6h: Character table missing in grp_k.F'
       ENDIF
    CASE DEFAULT
       WRITE(24,*) 'Group of k not identified'
    END SELECT

    !>

  END SUBROUTINE grp_k
  !-------------------------------------------------------------------------------

  !<--SUBROUTINE rotaxis(alpha,beta,gamma,rax,theta)
  SUBROUTINE rotaxis(alpha,beta,gamma,rax,theta)

    ! Determines the rotation axis based on the Euler angles
    IMPLICIT NONE
    REAL, INTENT(IN) :: alpha,beta,gamma
    REAL, INTENT(OUT) :: rax(4),theta 

    REAL :: costhe2,the2
    INTEGER :: i

    costhe2=COS(beta/2.0)*COS((alpha+gamma)/2.0)
    the2=ACOS(costhe2)
    rax=0.0
    IF (the2.GT.0.00001) THEN
       rax(1)=-SIN(beta/2.0)*SIN((alpha-gamma)/2.0)/SIN(the2)
       rax(2)= SIN(beta/2.0)*COS((alpha-gamma)/2.0)/SIN(the2)
       rax(3)= COS(beta/2.0)*SIN((alpha+gamma)/2.0)/SIN(the2)
    ENDIF
    DO i=1,3
       IF (ABS(rax(i)).GT.0.0001) THEN
          rax(i)=rax(i)/ABS(rax(i))
       ENDIF
    ENDDO
    rax(4)=rax(1)**2+rax(2)**2+rax(3)**2
    theta=the2*2.0
  END SUBROUTINE rotaxis
  !>

  !<--Subroutine euler
  SUBROUTINE euler(n,sym,cell,alpha,beta,gamma)

    ! determines the Euler angles corresponding the proper rotation part of mrot      

    USE m_constants, ONLY : pi_const
    USE m_inv3
    USE m_types
    IMPLICIT NONE
    INTEGER,INTENT(IN)         :: n
    TYPE(t_sym),INTENT(IN)     :: sym
    TYPE(t_cell),INTENT(IN)    :: cell
    
    REAL, INTENT(OUT) :: alpha,beta,gamma

    INTEGER :: det
    REAL :: mprop(3,3)      
    REAL :: sina,sinb,sinc,cosa,cosb,cosc
    REAL :: amatinv(3,3),detr

    CALL inv3(cell%amat,amatinv,detr)
    det= sym%mrot(1,1,n)*sym%mrot(2,2,n)*sym%mrot(3,3,n) +&
         &        sym%mrot(1,2,n)*sym%mrot(2,3,n)*sym%mrot(3,1,n) +&
         &        sym%mrot(2,1,n)*sym%mrot(3,2,n)*sym%mrot(1,3,n) -&
         &        sym%mrot(1,3,n)*sym%mrot(2,2,n)*sym%mrot(3,1,n) -&
         &        sym%mrot(2,3,n)*sym%mrot(3,2,n)*sym%mrot(1,1,n) -&
         &        sym%mrot(2,1,n)*sym%mrot(1,2,n)*sym%mrot(3,3,n)

    ! Take the proper rotation         
    mprop=REAL(det)*MATMUL(cell%amat,MATMUL(REAL(sym%mrot(:,:,n)),amatinv))
    ! Euler angles         
    cosb = mprop(3,3)
    sinb = 1.00 - cosb*cosb
    sinb = MAX(sinb,0.00)
    sinb = SQRT(sinb)
    !
    ! if beta = 0 or pi , only alpha+gamma or -gamma have a meaning:
    !
    IF ( ABS(sinb).LT.1.0e-5 ) THEN 
       beta = 0.0
       IF ( cosb.LT.0.0 ) beta = pi_const
       gamma = 0.0
       cosa = mprop(1,1)/cosb
       sina = mprop(1,2)/cosb
       IF ( ABS(sina).LT.1.0e-5 ) THEN
          alpha=0.0
          IF ( cosa.LT.0.0 ) alpha=alpha+pi_const
       ELSE
          alpha = 0.5*pi_const - ATAN(cosa/sina)
          IF ( sina.LT.0.0 ) alpha=alpha+pi_const
       ENDIF
    ELSE
       beta = 0.5*pi_const - ATAN(cosb/sinb)
       !     
       ! determine alpha and gamma from d13 d23 d32 d31
       !
       cosa = mprop(3,1)/sinb
       sina = mprop(3,2)/sinb
       cosc =-mprop(1,3)/sinb
       sinc = mprop(2,3)/sinb
       IF ( ABS(sina).LT.1.0e-5 ) THEN
          alpha=0.0
          IF ( cosa.LT.0.0 ) alpha=alpha+pi_const
       ELSE
          alpha = 0.5*pi_const - ATAN(cosa/sina)
          IF ( sina.LT.0.0 ) alpha=alpha+pi_const
       ENDIF
       IF ( ABS(sinc).LT.1.0e-5 ) THEN
          gamma = 0.0
          IF ( cosc.LT.0.0 ) gamma=gamma+pi_const
       ELSE
          gamma = 0.5*pi_const - ATAN(cosc/sinc)
          IF ( sinc.LT.0.0 ) gamma=gamma+pi_const
       ENDIF

    ENDIF

  END SUBROUTINE euler
  !>
  !-----------------------------------------------------------
  SUBROUTINE change_column( char_table,nx,ny,r1,r2 )

    INTEGER, INTENT (IN)   :: nx,ny,r1,r2
    COMPLEX, INTENT(INOUT) :: char_table(:,:)

    COMPLEX tmpc(nx)

    tmpc(1:nx) = char_table(1:nx,r1)
    char_table(1:nx,r1) = char_table(1:nx,r2)
    char_table(1:nx,r2) = tmpc(1:nx)
    tmpc(1:nx) = char_table(1:nx,r1+nx/2)
    char_table(1:nx,r1+nx/2) = char_table(1:nx,r2+nx/2)
    char_table(1:nx,r2+nx/2) = tmpc(1:nx)

  END SUBROUTINE change_column
END MODULE m_grp_k
