SetupNatOrbLabels Subroutine

    subroutine SetupNatOrbLabels()

        use MemoryManager, only: TagIntType

        integer :: x, i, j, ierr, NoOcc
        integer :: StartFill01, StartFill02, Symi, SymCurr, Prev, EndFill01, EndFill02
        character(len=*), parameter :: t_r = 'SetupNatOrbLabels'
        integer, allocatable :: LabVirtOrbs(:), LabOccOrbs(:), SymVirtOrbs(:), SymOccOrbs(:)
        integer(TagIntType) :: LabVirtOrbsTag, LabOccOrbsTag, SymVirtOrbsTag, SymOccOrbsTag
        integer :: lo, hi

        ! The earlier test should pick this up, if it crashes here, will want
        ! to put in an earlier test so that we don't get all the way to this
        ! stage.
        if ((LMS /= 0) .and. (.not. tStoreSpinOrbs)) then
            call stop_all("FindNatOrbs", "Open shell system, and UMAT is not being stored as spin orbitals.")
        end if

        ! We now need two slightly different sets of orbital labels for the
        ! case of spin orbitals and spatial orbitals. When using spin orbitals
        ! we want all the beta spin followed by all the alpha spin. Then we
        ! want two values for the number of occupied orbitals to allow for high
        ! spin cases. With spatial, it is equivalent to just keeping the beta
        ! spin.
        if (tStoreSpinOrbs) then
            NoSpinCyc = 2
        else
            NoSpinCyc = 1
        end if

        do x = 1, NoSpinCyc
            if (.not. tSeparateOccVirt) then
                NoOcc = 0
            else
                if (x == 1) then
                    if (tStoreSpinOrbs) then
                        NoOcc = nOccBeta
                    else
                        NoOcc = NEl / 2
                    end if
                end if
                if (x == 2) NoOcc = nOccAlpha
            end if

            if (tSeparateOccVirt) then
                allocate(LabOccOrbs(NoOcc), stat=ierr)
                call LogMemAlloc('LabOccOrbs', (NoOcc), 4, t_r, LabOccOrbsTag, ierr)
                if (ierr /= 0) call stop_all(t_r, "Mem allocation for LabOccOrbs failed.")
                LabOccOrbs(:) = 0
                allocate(SymOccOrbs(NoOcc), stat=ierr)
                call LogMemAlloc('SymOccOrbs', (NoOcc), 4, t_r, SymOccOrbsTag, ierr)
                if (ierr /= 0) call stop_all(t_r, "Mem allocation for SymOccOrbs failed.")
                SymOccOrbs(:) = 0
            end if

            allocate(LabVirtOrbs(SpatOrbs - NoOcc), stat=ierr)
            call LogMemAlloc('LabVirtOrbs', (SpatOrbs - NoOcc), 4, t_r, LabVirtOrbsTag, ierr)
            if (ierr /= 0) call stop_all(t_r, "Mem allocation for LabVirtOrbs failed.")
            LabVirtOrbs(:) = 0
            allocate(SymVirtOrbs(SpatOrbs - NoOcc), stat=ierr)
            call LogMemAlloc('SymVirtOrbs', (SpatOrbs - NoOcc), 4, t_r, SymVirtOrbsTag, ierr)
            if (ierr /= 0) call stop_all(t_r, "Mem allocation for SymVirtOrbs failed.")
            SymVirtOrbs(:) = 0

            ! First fill SymLabelList2_rot.

            ! Brr has the orbital numbers in order of energy...
            ! i.e Brr(2) = the orbital index with the second lowest energy.

            ! This picks out the NoOcc lowest energy orbitals from BRR as these
            ! will be the occupied. These are then ordered according to
            ! symmetry, and the same done to the virtual.
            do i = 1, NoOcc
                if (x == 1) then
                    if (tStoreSpinOrbs) then
                        LabOccOrbs(i) = BRR(2 * i) - 1
                        SymOccOrbs(i) = int(G1(LabOccOrbs(i))%sym%S, 4)
                    else
                        LabOccOrbs(i) = BRR(2 * i) / 2
                        SymOccOrbs(i) = int(G1(LabOccOrbs(i) * 2)%sym%S, 4)
                    end if
                else if (x == 2) then
                    LabOccOrbs(i) = BRR(2 * i)
                    SymOccOrbs(i) = int(G1(LabOccOrbs(i))%sym%S, 4)
                end if
            end do

            ! Sorts LabOrbs according to the order of SymOccOrbs (i.e. in
            ! terms of symmetry).
            if (tSeparateOccVirt) call sort(SymOccOrbs, LabOccOrbs)

            ! Same for the virtual.
            do i = 1, SpatOrbs - NoOcc
                if (x == 1) then
                    if (tStoreSpinOrbs) then
                        if (tSeparateOccVirt) then
                            LabVirtOrbs(i) = BRR((2 * i) + (NoOcc * 2)) - 1
                            SymVirtOrbs(i) = int(G1(LabVirtOrbs(i))%sym%S, 4)
                        else
                            LabVirtOrbs(i) = BRR((2 * i)) - 1
                            SymVirtOrbs(i) = int(G1(LabVirtOrbs(i))%sym%S, 4)
                        end if
                    else
                        if (tSeparateOccVirt) then
                            LabVirtOrbs(i) = BRR((2 * i) + (NoOcc * 2)) / 2
                            SymVirtOrbs(i) = int(G1(LabVirtOrbs(i) * 2)%sym%S, 4)
                        else
                            LabVirtOrbs(i) = BRR((2 * i)) / 2
                            SymVirtOrbs(i) = int(G1(LabVirtOrbs(i) * 2)%sym%S, 4)
                        end if
                    end if
                else if (x == 2) then
                    if (tSeparateOccVirt) then
                        LabVirtOrbs(i) = BRR((2 * i) + (NoOcc * 2))
                        SymVirtOrbs(i) = int(G1(LabVirtOrbs(i))%sym%S, 4)
                    else
                        LabVirtOrbs(i) = BRR((2 * i))
                        SymVirtOrbs(i) = int(G1(LabVirtOrbs(i))%sym%S, 4)
                    end if
                end if
            end do

            ! Sorts LabOrbs according to the order of SymOccOrbs (i.e. in
            ! terms of symmetry).
            call sort(SymVirtOrbs, LabVirtOrbs)

            ! SymLabelList2_rot is then filled with the symmetry ordered
            ! occupied then virtual arrays for each spin.
            if (x == 1) then
                StartFill01 = 1
                StartFill02 = NoOcc + 1
                EndFill01 = NoOcc
                EndFill02 = SpatOrbs
            else if (x == 2) then
                StartFill01 = SpatOrbs + 1
                StartFill02 = SpatOrbs + NoOcc + 1
                EndFill01 = SpatOrbs + NoOcc
                EndFill02 = NoOrbs
            end if

            j = 0
            do i = StartFill01, EndFill01
                j = j + 1
                SymLabelList2_rot(i) = LabOccOrbs(j)
            end do
            j = 0
            do i = StartFill02, EndFill02
                j = j + 1
                SymLabelList2_rot(i) = LabVirtOrbs(j)
            end do

            ! Second fill SymLabelCounts2_rot.
            ! - the first 8 places of SymLabelCounts2_rot(1,:) and
            !     SymLabelCounts2_rot(2,:) refer to the occupied orbitals
            ! - and the second 8 to the virtuals.

            if (lNoSymmetry) then
                ! If we are ignoring symmetry, all orbitals essentially have
                ! symmetry 0.
                if (x == 1) then
                    SymLabelCounts2_rot(1, 1) = 1
                    SymLabelCounts2_rot(1, 9) = NoOcc + 1
                    SymLabelCounts2_rot(2, 1) = NoOcc
                    SymLabelCounts2_rot(2, 9) = SpatOrbs - NoOcc
                else if (x == 2) then
                    SymLabelCounts2_rot(1, 17) = 1
                    SymLabelCounts2_rot(1, 25) = NoOcc + 1
                    SymLabelCounts2_rot(2, 17) = NoOcc
                    SymLabelCounts2_rot(2, 25) = SpatOrbs - NoOcc
                end if

            else
                ! Otherwise we run through the occupied orbitals, counting the
                ! number with each symmetry and noting where in
                ! SymLabelList2_rot each symmetry block starts.
                if (x == 1) then
                    StartFill01 = 1
                    StartFill02 = 9
                    Prev = 0
                else if (x == 2) then
                    StartFill01 = 17
                    StartFill02 = 25
                    Prev = SpatOrbs
                end if
                SymCurr = 0
                SymLabelCounts2_rot(1, StartFill01) = 1 + Prev
                do i = 1, NoOcc
                    if (tStoreSpinOrbs) then
                        Symi = int(G1(SymLabelList2_rot(i + Prev))%sym%S, 4)
                    else
                        Symi = int(G1((SymLabelList2_rot(i + Prev) * 2))%sym%S, 4)
                    end if
                    SymLabelCounts2_rot(2, (Symi + StartFill01)) = SymLabelCounts2_rot(2, (Symi + StartFill01)) + 1
                    if (Symi /= SymCurr) then
                        SymLabelCounts2_rot(1, (Symi + StartFill01)) = i + Prev
                        SymCurr = Symi
                    end if
                end do
                ! The same is then done for the virtuals.
                SymCurr = 0
                SymLabelCounts2_rot(1, StartFill02) = NoOcc + 1 + Prev
                do i = NoOcc + 1, SpatOrbs
                    if (tStoreSpinOrbs) then
                        Symi = int(G1(SymLabelList2_rot(i + Prev))%sym%S, 4)
                    else
                        Symi = int(G1((SymLabelList2_rot(i + Prev) * 2))%sym%S, 4)
                    end if
                    SymLabelCounts2_rot(2, (Symi + StartFill02)) = SymLabelCounts2_rot(2, (Symi + StartFill02)) + 1
                    if (Symi /= SymCurr) then
                        SymLabelCounts2_rot(1, (Symi + StartFill02)) = i + Prev
                        SymCurr = Symi
                    end if
                end do
            end if

            ! Go through each symmetry group, making sure the orbital pairs
            ! are ordered lowest to highest.
            if (x == 1) then
                do i = 1, 16
                    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
            else if (x == 2) then
                do i = 17, 32
                    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
            end if

            ! Deallocate the arrays just used in this routine.
            if (tSeparateOccVirt) then
                deallocate(SymOccOrbs)
                call LogMemDealloc(t_r, SymOccOrbsTag)

                deallocate(LabOccOrbs)
                call LogMemDealloc(t_r, LabOccOrbsTag)
            end if

            deallocate(SymVirtOrbs)
            call LogMemDealloc(t_r, SymVirtOrbsTag)

            deallocate(LabVirtOrbs)
            call LogMemDealloc(t_r, LabVirtOrbsTag)

        end do

        do i = 1, NoOrbs
            SymLabelListInv_rot(SymLabelList2_rot(i)) = i
        end do

        do i = 1, NoOrbs
            SymLabelList3_rot(i) = SymLabelList2_rot(i)
        end do

        if (.not. tSeparateOccVirt) then
            ! Basically we treat all the orbitals as virtuals and set NoOcc
            ! to zero in each routine.
            tRotateVirtOnly = .true.
        end if

    end subroutine SetupNatOrbLabels