HistNatOrbEvalues Subroutine

    subroutine HistNatOrbEvalues()

        use LoggingData, only: tTruncRODump
        use RotateOrbsData, only: CoeffT1, EvaluesTrunc

        integer :: i, k, a, b, NoOcc, io1, io2
        real(dp) :: EvalueEnergies(1:NoOrbs), OrbEnergies(1:NoOrbs)
        real(dp) :: SumEvalues

        io1 = get_free_unit()
        NoOcc = NEl / 2 ! Is this correct in all cases?!

        open(io1, file='EVALUES-PLOTRAT', status='unknown')
        if (tStoreSpinOrbs) then
            k = 0
            do i = 1, SpatOrbs
                k = k + 2
                write(io1, '(F20.10,ES20.10)') real(k - 1, dp) / real(NoOrbs, dp), Evalues(i)
                write(io1, '(F20.10,ES20.10)') real(k, dp) / real(NoOrbs, dp), Evalues(SpatOrbs + i)
            end do
        else if (tRotateOccOnly) then
            k = 0
            do i = 1, NoOcc
                k = k + 1
                write(io1, '(F20.10,ES20.10)') real(k, dp) / real(NoOcc, dp), Evalues(i)
            end do
        else if (tRotateVirtOnly) then
            k = NoOcc
            do i = NoOcc + 1, NoOrbs
                k = k + 1
                write(io1, '(F20.10,ES20.10)') real(k - NoOcc, dp) / real(NoOrbs - NoOcc, dp), Evalues(i)
            end do
        else
            k = 0
            do i = 1, SpatOrbs
                k = k + 1
                write(io1, '(F20.10,ES20.10)') real(k, dp) / real(NoOrbs, dp), Evalues(i)
            end do
        end if
        close(io1)

        ! Want to write out the eigenvectors in order of the energy of the new
        ! orbitals - so that we can see the occupations of the type of orbital.
        ! For now, keep this separate to the transformation of ARR - even
        ! though it is equivalent.

        OrbEnergies(:) = 0.0_dp
        EvalueEnergies(:) = 0.0_dp
        SumEvalues = 0.0_dp
        do i = 1, NoOrbs
            if (tStoreSpinOrbs) then
                SumEvalues = SumEvalues + Evalues(i)
            else
                SumEvalues = SumEvalues + (2 * Evalues(i))
            end if

            if (tTruncRODump) then
                EvalueEnergies(i) = EvaluesTrunc(i)
            else
                EvalueEnergies(i) = Evalues(i)
            end if

            ! We are only interested in the diagonal elements.
            do a = 1, NoOrbs
                b = SymLabelList2_rot(a)
                if (tStoreSpinOrbs) then
                    OrbEnergies(i) = OrbEnergies(i) + (CoeffT1(a, i) * ARR(b, 2) * CoeffT1(a, i))
                else
                    OrbEnergies(i) = OrbEnergies(i) + (CoeffT1(a, i) * ARR(2 * b, 2) * CoeffT1(a, i))
                end if
            end do
        end do

        ! ARR(:,1) - ordered by energy, ARR(:,2) - ordered by spin-orbital
        ! index.

        call sort(orbEnergies(1:noOrbs), EvalueEnergies(1:noOrbs))

        io2 = get_free_unit()
        open(io2, file='EVALUES-ENERGY', status='unknown')
        do i = 1, NoOrbs
            write(io2, *) OrbEnergies(NoOrbs - i + 1), EvalueEnergies(NoOrbs - i + 1)
        end do
        write(io2, *) 'The sum of the occupation numbers (eigenvalues) = ', SumEvalues
        write(io2, *) 'The number of electrons = ', NEl
        call neci_flush(io2)
        close(io2)
        call neci_flush(stdout)

        call PrintOccTable()

    end subroutine HistNatOrbEvalues