find_nat_orb_occ_numbers Subroutine

    subroutine find_nat_orb_occ_numbers(rdm, irdm)

        use LoggingData, only: tPrintRODump
        use MemoryManager, only: LogMemAlloc
        use Parallel_neci, only: iProcIndex
        use rdm_data, only: one_rdm_t
        use RotateOrbsMod, only: FourIndInts, FourIndIntsTag
        use SystemData, only: tROHF, nbasis, G1, ARR, BRR

        ! Diagonalises the 1-RDM (rdm%matrix) so that, after this routine,
        ! rdm%matrix holds the eigenfunctions of the 1-RDM (the matrix
        ! transforming the MO's into the NOs). This also gets the NO
        ! occupation numbers (evaluse) and correlation entropy.

        type(one_rdm_t), intent(inout) :: rdm
        integer, intent(in) :: irdm

        integer :: ierr
        real(dp) :: SumDiag
        character(len=*), parameter :: t_r = 'find_nat_orb_occ_numbers'

        ierr = 0
        if (iProcIndex == 0) then

            ! Diagonalises the 1-RDM. rdm%matrix goes in as the 1-RDM, comes out
            ! as the eigenvector of the 1-RDM (the matrix transforming the MO's
            ! into the NOs).
            call DiagRDM(rdm, SumDiag)

            ! Writes out the NO occupation numbers and evectors to files.
            call write_evales_and_transform_mat(rdm, irdm, SumDiag)

            if (tPrintRODump .and. tROHF) then
                write(stdout, *) 'ROFCIDUMP not implemented for ROHF. Skip generation of ROFCIDUMP file.'
            else if (tPrintRODump) then
                write(stdout, "(A,F10.5,A)") "This will require at least ", (real(NoOrbs, dp)**4) * 8 / 10**9, "Gb to be available on head node"
                allocate(FourIndInts(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
                call LogMemAlloc('FourIndInts', (NoOrbs**4), 8, t_r, &
                                 FourIndIntsTag, ierr)
                if (ierr /= 0) then
                    write(stdout, *) "ierr: ", ierr
                    call Stop_All(t_r, 'Problem allocating FourIndInts array,')
                end if

                ! Then, transform2ElInts.
                write(stdout, *) ''
                write(stdout, *) 'Transforming the four index integrals'
                call Transform2ElIntsMemSave_RDM(rdm%matrix, rdm%sym_list_no)

                write(stdout, *) 'Re-calculating the fock matrix'
                call CalcFOCKMatrix_RDM(rdm)

                write(stdout, *) 'Refilling the UMAT and TMAT2D'

                ! The ROFCIDUMP is also printed out in here.
                call RefillUMATandTMAT2D_RDM(rdm%matrix, rdm%sym_list_no)

                call neci_flush(stdout)

                call writebasis(stdout, G1, nbasis, ARR, BRR)
            end if
        end if

    end subroutine find_nat_orb_occ_numbers