order_one_rdm Subroutine

    subroutine order_one_rdm(rdm)

        ! Here, if symmetry is kept, we are going to have to reorder the
        ! eigenvectors according to the size of the eigenvalues, while taking
        ! the orbital labels (and therefore symmetries) with them. This will
        ! be put back into MP2VDM from MP2VDMTemp.

        ! Want to reorder the eigenvalues from largest to smallest, taking the
        ! eigenvectors with them and the symmetry as well. If using spin
        ! orbitals, do this for the alpha spin and then the beta.

        ! The newly sorted SymLabelList2_rot and SymLabelListInv_rot will be
        ! stored in rdm%sym_list_no and rdm%sym_list_inv_no (as they are
        ! specific to this RDM).

        use MemoryManager, only: LogMemAlloc
        use rdm_data, only: tOpenShell, one_rdm_t
        use sort_mod, only: sort
        use SystemData, only: nbasis

        type(one_rdm_t), intent(inout) :: rdm

        integer :: spin, i, j, ierr, StartSort, EndSort
        character(len=*), parameter :: t_r = 'order_one_rdm'
        integer, allocatable :: SymLabelList_temp(:)
        real(dp), allocatable :: one_rdm_Temp(:, :), EvaluesTemp(:)
        integer :: Orb, New_Pos

        ! Temporary arrays.
        allocate(one_rdm_Temp(NoOrbs, NoOrbs), stat=ierr)
        allocate(SymLabelList_temp(NoOrbs), stat=ierr)
        allocate(EvaluesTemp(NoOrbs), stat=ierr)

        ! This is just a temporary array for some sorting.
        SymLabelList_temp = SymLabelList2_rot

        ! This object will contain the updated version of SymLabelList2_rot,
        ! specifically for this RDM, and ordered with the 1-RDM eigenvalues.
        rdm%sym_list_no = SymLabelList2_rot

        StartSort = 1
        EndSort = SpatOrbs

        ! Unfortunately this sort routine orders the orbitals in ascending
        ! order... which is not quite what we want.  Just remember this when
        ! printing out the Evalues.
        call sort(rdm%EValues(startSort:endSort), &
                  rdm%matrix(1:NoOrbs, startSort:endSort), &
                  rdm%sym_list_no(startSort:endSort))

        if (tOpenShell) then
            StartSort = SpatOrbs + 1
            EndSort = nBasis

            call sort(rdm%EValues(startSort:endSort), &
                      rdm%matrix(1:NoOrbs, startSort:endSort), &
                      rdm%sym_list_no(startSort:endSort))

        end if

        ! We now have the NO's ordered according to the size of their Evalues
        ! (occupation numbers).  This will have jumbled up their symmetries.
        ! Want to reorder the MO's to match this ordering (so that we only
        ! have one SymLabelList array).

        ! Need a new SymLabelListInv_rot too. This will be stored in
        ! rdm%sym_list_inv_no.
        rdm%sym_list_inv_no = 0
        do i = 1, NoOrbs
            rdm%sym_list_inv_no(rdm%sym_list_no(NoOrbs - i + 1)) = i
        end do

        one_rdm_Temp = rdm%matrix
        rdm%matrix = 0.0_dp

        do i = 1, NoOrbs
            do j = 1, NoOrbs

                ! In position j, the MO orbital Orb is currently there.
                Orb = SymLabelList_temp(j)

                ! Want to move it to the position the NO's are in.
                New_Pos = rdm%sym_list_inv_no(Orb)

                ! But we also want to reverse the order of everything...
                rdm%matrix(New_Pos, NoOrbs - i + 1) = one_rdm_Temp(j, i)
            end do
        end do

        SymLabelList_temp = rdm%sym_list_no
        EvaluesTemp = rdm%evalues
        do i = 1, NoOrbs
            rdm%sym_list_no(i) = SymLabelList_temp(NoOrbs - i + 1)
            rdm%evalues(i) = EvaluesTemp(NoOrbs - i + 1)
        end do

        deallocate(one_rdm_Temp)
        deallocate(SymLabelList_temp)
        deallocate(EvaluesTemp)

    end subroutine order_one_rdm