SetUpSymLabels_RDM Subroutine

    subroutine SetUpSymLabels_RDM()

        ! This routine just sets up the symmetry labels so that the orbitals
        ! are ordered according to symmetry (all beta then all alpha if spin orbs).

        use rdm_data, only: tOpenShell
        use RotateOrbsData, only: SymLabelList2_rot, SymLabelCounts2_rot, SymLabelListInv_rot
        use RotateOrbsData, only: NoOrbs, SpatOrbs, NoSymLabelCounts
        use sort_mod, only: sort
        use SystemData, only: G1, BRR, lNoSymmetry, tFixLz, iMaxLz
        use UMatCache, only: gtID
        use MemoryManager, only: LogMemAlloc, LogMemDealloc

        integer, allocatable :: SymOrbs_rot(:)
        integer :: SymOrbs_rotTag, ierr, i, j, SpatSym, LzSym
        integer :: lo, hi, Symi, SymCurr, Symi2, SymCurr2
        character(len=*), parameter :: t_r = 'SetUpSymLabels_RDM'

        ! This is only allocated temporarily to be used to order the orbitals by.
        allocate(SymOrbs_rot(NoOrbs), stat=ierr)
        call LogMemAlloc('SymOrbs_rot', NoOrbs, 4, t_r, SymOrbs_rotTag, ierr)
        if (ierr /= 0) call stop_all(t_r, "Mem allocation for SymOrbs_rot failed.")

        ! Now we want to put the spatial orbital index, followed by the symmetry.
        SymLabelList2_rot(:) = 0
        SymOrbs_rot(:) = 0

        ! *** STEP 1 *** Fill SymLabelList2_rot.
        ! Find the orbitals and order them in terms of symmetry.
        do i = 1, SpatOrbs
            if (tOpenShell) then
                ! For open shell systems, all alpha are followed by all beta.
                SymLabelList2_rot(i) = BRR(2 * i)
                SymLabelList2_rot(i + SpatOrbs) = BRR((2 * i) - 1)

                if (tFixLz) then
                    SpatSym = int(G1(BRR(2 * i))%sym%S)
                    LzSym = int(G1(BRR(2 * i))%Ml)
                    SymOrbs_rot(i) = (SpatSym * ((2 * iMaxLz) + 1)) + (LzSym + iMaxLz)

                    SpatSym = int(G1(BRR((2 * i) - 1))%sym%S)
                    LzSym = int(G1(BRR((2 * i) - 1))%Ml)
                    SymOrbs_rot(i + SpatOrbs) = (SpatSym * ((2 * iMaxLz) + 1)) + (LzSym + iMaxLz)
                else
                    SymOrbs_rot(i) = int(G1(BRR(2 * i))%sym%S)
                    SymOrbs_rot(i + SpatOrbs) = int(G1(BRR((2 * i) - 1))%sym%S)
                end if
            else
                SymLabelList2_rot(i) = gtID(BRR(2 * i))
                if (tFixLz) then
                    SpatSym = int(G1(BRR(2 * i))%sym%S)
                    LzSym = int(G1(BRR(2 * i))%Ml)
                    SymOrbs_rot(i) = (SpatSym * ((2 * iMaxLz) + 1)) + (LzSym + iMaxLz)
                else
                    SymOrbs_rot(i) = int(G1(BRR(2 * i))%sym%S)
                end if
                ! Orbital BRR(2*i) for i = 1 will be the beta orbital with the
                ! second lowest energy - want the spatial orbital index to go with this.
                ! G1 also in spin orbitals - get symmetry of this beta orbital, will
                ! be the same as the spatial orbital.
            end if
        end do

        call sort(SymOrbs_rot(1:SpatOrbs), SymLabelList2_rot(1:SpatOrbs))
        ! Sorts SymLabelList2_rot according to the order of SymOrbs_rot
        ! (i.e. in terms of symmetry).
        if (tOpenShell) &
            call sort(SymOrbs_rot(SpatOrbs + 1:nBasis), SymLabelList2_rot(SpatOrbs + 1:nBasis))
        ! Also do this for the beta set if spin orbitals.

        ! *** STEP 2 *** Fill SymLabelCounts2_rot_rot. This is like
        ! SymLabelCounts2_rot, but has a different ordering - BEWARE.
        ! SymLabelCounts(1,:) contains the position in SymLabelList2_rot
        ! where the symmetry index starts, SymLabelCounts(2,:) contains the
        ! number of orbitals in that symmetry index. Again if spin orbs, all
        ! alpha are followed by all beta - i.e. first 8 refer to alpha, second
        ! 8 to beta.

        if (lNoSymmetry) then
            ! If we are ignoring symmetry, all orbitals essentially have
            ! symmetry 0.
            SymLabelCounts2_rot(1, 1) = 1
            SymLabelCounts2_rot(2, 1) = SpatOrbs
            if (tOpenShell) then
                SymLabelCounts2_rot(1, 9) = SpatOrbs + 1
                SymLabelCounts2_rot(2, 9) = SpatOrbs
            end if
        else
            ! Otherwise we run through the occupied orbitals, counting the
            ! number with each symmetry (spatial and Lz) and noting where in
            ! SymLabelList2_rot each symmetry block starts.
            SymCurr = 0
            SymLabelCounts2_rot(1, 1) = 1
            if (tOpenShell) then
                SymCurr2 = 0
                SymLabelCounts2_rot(1, 9) = SpatOrbs + 1
            end if
            do i = 1, SpatOrbs
                if (tOpenShell) then
                    Symi = SymOrbs_rot(i)
                    Symi2 = SymOrbs_rot(i + SpatOrbs)
                else
                    Symi = SymOrbs_rot(i)
                end if
                SymLabelCounts2_rot(2, (Symi + 1)) = SymLabelCounts2_rot(2, (Symi + 1)) + 1
                if (Symi /= SymCurr) then
                    do j = SymCurr + 1, Symi
                        SymLabelCounts2_rot(1, (j + 1)) = i
                    end do
                    SymCurr = Symi
                end if
                if (tOpenShell) then
                    SymLabelCounts2_rot(2, (Symi2 + 9)) = SymLabelCounts2_rot(2, (Symi2 + 9)) + 1
                    if (Symi2 /= SymCurr2) then
                        do j = SymCurr2 + 1, Symi2
                            SymLabelCounts2_rot(1, (j + 9)) = i + SpatOrbs
                        end do
                        SymCurr2 = Symi2
                    end if
                end if
            end do
        end if

        ! Go through each symmetry group, making sure the orbitals are
        ! ordered lowest to highest within each symmetry.
        do i = 1, NoSymLabelCounts
            if (SymLabelCounts2_rot(2, i) /= 0) then
                lo = SymLabelCounts2_rot(1, i)
                hi = lo + SymLabelCounts2_rot(2, i) - 1
                call sort(SymLabelList2_rot(lo:hi))
            end if
        end do

        ! Construct the inverse matrix. While SymLabelList2_rot takes a
        ! position and tells us what orbital is in it, we also might need to
        ! take an orbital and find out what position to put its contribution in.
        do i = 1, NoOrbs
            SymLabelListInv_rot(SymLabelList2_rot(i)) = i
        end do

        ! Deallocate the arrays just used in this routine.
        deallocate(SymOrbs_rot)
        call LogMemDealloc(t_r, SymOrbs_rotTag)

    end subroutine SetUpSymLabels_RDM