DIAGFMAT Subroutine

    SUBROUTINE DIAGFMAT(NSPINS, NSBASIS, NELS, FMAT, DMAT, HFES, WORK, ECORE, ECUR)
        INTEGER NSPINS, NSBASIS, NELS(NSPINS)
        real(dp) HFES(NSBASIS, NSPINS)
        HElement_t(dp) FMAT(NSBASIS, NSBASIS, NSPINS)
        HElement_t(dp) DMAT(NSBASIS, NSBASIS, NSPINS)
!         real(dp) TMAT(NSBASIS*NSPINS,NSBASIS*NSPINS)
        real(dp) ECUR, ECORE
        HElement_t(dp) WORK(3 * NSBASIS), RWORK(3 * NSBASIS)

        INTEGER(int32) :: INFO
        HElement_t(dp) TOT, TOT2
        INTEGER ISPN, I, J, K
        character(*), parameter :: this_routine = 'DIAGFMAT'

!.. First calculate the HF energy double counting contrib
        TOT = 0.0_dp
        DO ISPN = 1, NSPINS
        DO J = 1, NSBASIS
        DO K = 1, NSBASIS
            TOT = TOT + DMAT(J, K, ISPN) * (FMAT(J, K, ISPN) - GetTMATEl((J - 1) * NSPINS + ISPN,(K - 1) * NSPINS + ISPN))
        END DO
        END DO
        END DO

!.. Now diagonalize the Fock matrix
        DO ISPN = 1, NSPINS
        IF (HElement_t_size == 1) THEN
            CALL DSYEV('V', 'U', NSBASIS, FMAT(1, 1, ISPN), NSBASIS, HFES(1, ISPN), WORK, 3 * NSBASIS, INFO)
        ELSE
            CALL ZHEEV('V', 'U', NSBASIS, FMAT(1, 1, ISPN), NSBASIS, HFES(1, ISPN), WORK, 3 * NSBASIS, RWORK, INFO)
        END IF
!.. eigenvector N is in FMAT(i,N,ISPN), where i is the component of the
!.. vector
        IF (INFO /= 0) THEN
            WRITE(stdout, *) 'DYSEV error: ', INFO
            call stop_all(this_routine, "DYSEV error")
        END IF
        IF (HFLogLevel > 0) CALL WRITE_HEMATRIX("EIGVEC", NSBASIS, NSBASIS, FMAT(1, 1, ISPN))
        IF (HFLogLevel > 0) CALL NECI_WRITE_MATRIX("EIGVAL", 1, NSBASIS, HFES(1, ISPN))
        END DO
!.. now calculate the sum of the occupied Fock orbitals
        TOT2 = 0.0_dp
        DO ISPN = 1, NSPINS
        DO I = 1, NELS(ISPN)
            TOT2 = TOT2 + (HFES(I, ISPN))
        END DO
        END DO
!.. subtract out the double counting term, and add in the core energy
        ECUR = (TOT2) - (TOT) / 2.0_dp + ECORE
    END subroutine diagfmat