SUBROUTINE symham (n, hamil)
USE Numeric_Kinds_Module
USE FileUnits_Module
USE SymGroup_Module
USE CommonInfo_Module
IMPLICIT NONE
!==============================================================================
!    I N P U T / O U T P U T
INTEGER  :: n 
REAL(dp) :: hamil (n, n)
!==============================================================================
!    I N T E R N A L S
LOGICAL :: debug,endpt
CHARACTER(len=2) :: label2,l2
INTEGER :: gamma, irrep, j, i, iw, isymop, isum, k, nlast, index 
INTEGER :: fac,nb, ne
!==============================================================================
!    A L L O C A T A B L E
INTEGER,ALLOCATABLE  :: iblock(:,:),n1(:)
REAL(dp),ALLOCATABLE :: ttmp(:,:), tblock(:,:), norm(:)
!==============================================================================
! Allocate necessary arrays
 ALLOCATE( tblock(n,n), ttmp(n,n) )
 ALLOCATE( iblock(n,n) )
 ALLOCATE( n1(8), norm(n) )

 debug = .false.

!==============================================================================
! Zero out i block and t block
 DO i = 1, n
   DO j = 1, n
     iblock (i, j) = 0
     tblock (i, j) = 0.d0
   ENDDO
 ENDDO

!==============================================================================
! Begin symmetrization
 nlast = 0
 CALL symop (irrep, 1, i, gamma, n, isymop, n1, 0, fac)

! Loop over irred. reps
 DO irrep = 1, nirrep ! DO1
   nchi_ir(irrep) = n1(irrep)
   iw = 0                   ! Counting index
   nb = nlast + 1           ! Lower grid index of irred rep
   ne = nlast + n1 (irrep)  ! Upper grid index of irred rep
   DO j = nb, ne  ! DO2
     isum = 0
     DO while (isum.eq.0) !DO3
       ! Loop over symmetry operations for a fixed j , nlast , iw
       DO isymop = 1, nsymop !DO4
         index = j - nlast + iw ! Counts from 1 to dim of irrep
         CALL symop (irrep, index, i, gamma, n, isymop, n1, 1, fac)
         iblock (i, j) = gamma + iblock (i, j)
         ! For 2D representations?
         IF (fac.ne.0.and.i.ne.1) THEN
           iblock (n + 2 - i, j) = iblock (n + 2 - i, j) + fac * gamma
         ELSEIF (i.eq.1) THEN
           iblock (i, j) = iblock (i, j) + fac * gamma
         ENDIF !
       ENDDO !ENDDO4
      
       isum = 0
       DO k = 1, n !D05
         isum = isum + iblock (k, j) * iblock (k, j)
       ENDDO !ENDDO5

       IF (isum.eq.0) iw = iw + 1

     ENDDO !ENDDO3 DO while !
     IF (j.ne.nlast + 1) THEN
       CALL ischmidt (iblock (1, nlast + 1), n, j - nlast - 1, iblock (1, j) )
     ENDIF
   ENDDO !ENDDO2

   nlast = nlast + n1 (irrep)
   DO i = nb, ne
     DO j = i, ne
       isum = 0
       DO k = 1, n
         isum = isum + iblock (k, i) * iblock (k, j)
       ENDDO
     ENDDO
   ENDDO

 ENDDO !ENDDO1

 DO i = 1, n
   j = i
   isum = 0
   DO k = 1, n
     isum = isum + iblock (k, i) * iblock (k, j)
   ENDDO
   IF (i.eq.j) THEN
     IF (isum.eq.0) THEN
       WRITE(Output_Unit, * ) 'error ,i,j,isum=', i, j, isum
     ENDIF
     norm (i) = isum
     norm (i) = 1.d0 / sqrt (norm (i) )
   ELSE
     IF (isum.ne.0) THEN
       WRITE(Output_Unit, * ) 'error', i, j, isum
     ENDIF !
   ENDIF
 ENDDO

 IF (debug) THEN
   WRITE(Output_Unit, * ) 'iblock matrix'
   DO i = 1, min (n, 25)
   WRITE(Output_Unit, 10) (iblock (i, j), j = 1, min (n, 25) )
   ENDDO
 ENDIF
 10 FORMAT(1x,25i3)

 DO i = 1, n
   DO j = 1, n
     tblock (i, j) = 0.d0
     DO k = 1, n
       tblock (i, j) = tblock (i, j) + hamil (i, k) * iblock (k, j)
     ENDDO
   ENDDO
 ENDDO

 DO i = 1, n
   DO j = 1, n
     hamil (i, j) = tblock (i, j) * norm (j)
   ENDDO
 ENDDO


 DO i = 1, n
   DO j = 1, n
     tblock (i, j) = 0.d0
     DO k = 1, n
       tblock (i, j) = tblock (i, j) + iblock (k, i) * hamil (k, j)
     ENDDO
   ENDDO
 ENDDO 


 DO i = 1, n
   DO j = 1, n
     hamil (i, j) = tblock (i, j) * norm (i)
   ENDDO
 ENDDO

 DO j = 1, n
   DO i=1,n
    ttmp(i,j) = iblock(i,j)*norm(j)
   ENDDO
 ENDDO

 OPEN(binout_unit0,file=TRIM(BinOutdir)//'symmat.bin',form='unformatted')
 WRITE(binout_unit0) ttmp
 CLOSE(binout_unit0)

 nlast = 1
 DO irrep = 1, nirrep
   DO i = 1, n
     DO j = 1, n
       IF (i.ge.nlast.and.i.lt.nlast + n1 (irrep) ) THEN
         IF (j.lt.nlast.or.j.ge.nlast + n1 (irrep) ) THEN
           IF (abs (hamil (i, j) ) .lt.1.d-10)THEN
             hamil (i, j) = 0.d0
           ELSE
             WRITE(* , * ) 'irrep,i,j=', irrep, i, j
             WRITE(* , * ) 'Error hamil=', hamil (i, j),abs(hamil (i, j)), &
                                     1.d-10,abs (hamil (i, j) ) .lt.1.d-10
             STOP 'symham'
           ENDIF
         ENDIF
       ENDIF
     ENDDO
   ENDDO
   nlast = nlast + n1 (irrep)
 ENDDO

 IF (debug) THEN
   WRITE(Output_Unit, * ) 'Hamiltonian matrix'
   CALL mxout (hamil, n, n)
 ENDIF
 20 FORMAT(1x,24f5.0)

! DO i = 1, n
!   DO j = 1, n
!     tblock (i, j) = iblock (i, j) * norm (j)
!   ENDDO
! ENDDO

 DEALLOCATE(iblock)
 DEALLOCATE(tblock, ttmp, norm)

RETURN
END SUBROUTINE symham