Diagonalizehij Subroutine

    subroutine Diagonalizehij()

! This routine takes the original <i|h|j> matrix and diagonalises it.  The resulting coefficients from this process
! are then the rotation coefficients to be applied to the four index integrals etc.
! This eliminates the <i|h|j> elements from the single excitations, and leaves only coulomb and exchange terms.
! In order to maintain the same HF energy, only the virtual elements are diagonalised, within symmetry blocks.

        integer :: i, j, Sym, ierr, NoSymBlock, WorkSize, WorkCheck, SymStartInd
        integer(TagIntType) WorkTag, DiagTMAT2DBlockTag, TMAT2DSymBlockTag
        real(dp), allocatable :: TMAT2DSymBlock(:, :), DiagTMAT2DBlock(:), Work(:)
        character(len=*), parameter :: this_routine = 'Diagonalizehij'

        write(stdout, *) 'The original coefficient matrix'
        do i = 1, NoOrbs
            do j = 1, NoOrbs
                write(stdout, '(F20.10)', advance='no') CoeffT1(j, i)
            end do
            write(stdout, *) ''
        end do

        write(stdout, *) 'The original TMAT2D matrix'
        do i = 1, NoOrbs
            do j = 1, NoOrbs
                write(stdout, '(F20.10)', advance='no') TMAT2DTemp(j, i)
            end do
            write(stdout, *) ''
        end do
        TMAT2DRot(:, :) = 0.0_dp
        DiagTMAT2Dfull(:) = 0.0_dp

! Now need to pick out symmetry blocks, from the virtual orbitals and diagonalize them.

! Take first symmetry, (0) and find the number of virtual orbitals with this symmetry.  If this is greater than 1,
! take the block, diagonlize it, and put it into TMAT2DRot.

        Sym = 0
        WorkSize = -1
        do while (Sym <= 7)

            NoSymBlock = SymLabelCounts2_rot(2, Sym + 9)

            SymStartInd = SymLabelCounts2_rot(1, Sym + 9) - 1
            ! This is one less than the index that the symmetry starts, so that when we run through i = 1,..., we can
            ! start at SymStartInd+i.

            if (NoSymBlock > 1) then
                allocate(TMAT2DSymBlock(NoSymBlock, NoSymBlock), stat=ierr)
                call LogMemAlloc('TMAT2DSymBlock', NoSymBlock**2, 8, this_routine, TMAT2DSymBlockTag, ierr)
                allocate(DiagTMAT2DBlock(NoSymBlock), stat=ierr)
                call LogMemAlloc('DiagTMAT2DBlock', NoSymBlock, 8, this_routine, DiagTMAT2DBlockTag, ierr)

                WorkCheck = 3 * NoSymBlock + 1
                WorkSize = WorkCheck
                allocate(Work(WorkSize), stat=ierr)
                call LogMemAlloc('Work', WorkSize, 8, this_routine, WorkTag, ierr)

                do j = 1, NoSymBlock
                    do i = 1, NoSymBlock
                        TMAT2DSymBlock(i, j) = TMAT2DTemp(SymStartInd + i, SymStartInd + j)
                    end do
                end do

                write(stdout, *) '*****'
                write(stdout, *) 'Symmetry ', Sym, ' has ', NoSymBlock, ' orbitals .'
                write(stdout, *) 'The TMAT2D for this symmetry block is '
                do i = 1, NoSymBlock
                    do j = 1, NoSymBlock
                        write(stdout, '(F20.10)', advance='no') TMAT2DSymBlock(j, i)
                    end do
                    write(stdout, *) ''
                end do

                call DSYEV('V', 'U', NoSymBlock, TMAT2DSymBlock, NoSymBlock, DiagTMAT2Dblock, Work, WorkSize, ierr)
                ! TMAT2DSymBlock goes in as the original TMAT2DSymBlock, comes out as the eigenvectors (Coefficients).
                ! TMAT2DBlock comes out as the eigenvalues in ascending order.
                if (ierr /= 0) then
                    write(stdout, *) 'Problem with symmetry, ', Sym, ' of TMAT2D'
                    call neci_flush(stdout)
                    call Stop_All(this_routine, "Diagonalization of TMAT2DSymBlock failed...")
                end if

                write(stdout, *) 'After diagonalization, the e-vectors (diagonal elements) of this matrix are,'
                do i = 1, NoSymBlock
                    write(stdout, '(F20.10)', advance='no') DiagTMAT2Dblock(i)
                end do
                write(stdout, *) ''
                write(stdout, *) 'These go from orbital,', SymStartInd + 1, ' to ', SymStartInd + NoSymBlock

                do i = 1, NoSymBlock
                    DiagTMAT2Dfull(SymStartInd + i - NoOcc) = DiagTMAT2DBlock(i)
                end do

                ! CAREFUL if eigenvalues are put in ascending order, this may not be correct, with the labelling system.
                ! may be better to just take coefficients and transform TMAT2DRot in transform2elints.
                ! a check that comes out as diagonal is a check of this routine anyway.

                write(stdout, *) 'The eigenvectors (coefficients) for symmtry block ', Sym
                do i = 1, NoSymBlock
                    do j = 1, NoSymBlock
                        write(stdout, '(F20.10)', advance='no') TMAT2DSymBlock(j, i)
                    end do
                    write(stdout, *) ''
                end do

                ! Directly fill the coefficient matrix with the eigenvectors from the diagonalization.
                do j = 1, NoSymBlock
                    do i = 1, NoSymBlock
                        CoeffT1(SymStartInd + i, SymStartInd + j) = TMAT2DSymBlock(i, j)
                    end do
                end do

                deallocate(Work)
                call LogMemDealloc(this_routine, WorkTag)

                deallocate(DiagTMAT2DBlock)
                call LogMemDealloc(this_routine, DiagTMAT2DBlockTag)

                deallocate(TMAT2DSymBlock)
                call LogMemDealloc(this_routine, TMAT2DSymBlockTag)
            else if (NoSymBlock == 1) then
                DiagTMAT2Dfull(SymStartInd + 1 - NoOcc) = TMAT2DTemp(SymStartInd + 1, SymStartInd + 1)
                write(stdout, *) '*****'
                write(stdout, *) 'Symmetry ', Sym, ' has only one orbital.'
                write(stdout, *) 'Copying diagonal element,', SymStartInd + 1, 'to DiagTMAT2Dfull'
            end if

            Sym = Sym + 1
        end do

        write(stdout, *) '*****'
        write(stdout, *) 'The final coefficient matrix'
        do i = 1, NoOrbs
            do j = 1, NoOrbs
                write(stdout, '(F20.10)', advance='no') CoeffT1(j, i)
            end do
            write(stdout, *) ''
        end do

        write(stdout, *) '*****'
        write(stdout, *) 'The diagonal elements of TMAT2D'
        do i = 1, (NoOrbs - NoOcc)
            write(stdout, *) DiagTMAT2Dfull(i)
        end do

    end subroutine Diagonalizehij