InitOrbitalSeparation Subroutine

    subroutine InitOrbitalSeparation()

! This subroutine is called if the SEPARATEOCCVIRT keyword is present in the input, it sets up SymLabelList2_rot so that the first
! NoOcc orbitals are the HF occupied, and the rest the virtual.  Within this separation, orbitals are ordered in symmetry
! groups.
! This means that two iterations of the rotate orbs routine will be performed, the first treats the occupied orbitals and the second
! the virtual.

        integer :: i, j, ierr, SymCurr, Symi
        integer(TagIntType) :: SymVirtOrbsTag, SymOccOrbsTag
        integer :: lo, hi
        integer, allocatable :: SymVirtOrbs(:), SymOccOrbs(:)
        character(len=*), parameter :: this_routine = 'InitOrbitalSeparation'

        allocate(SymLabelCounts2_rot(2, 16), stat=ierr)
        call LogMemAlloc('SymLabelCounts2_rot', 2 * 16, 4, this_routine, SymLabelCounts2_rotTag, ierr)
        SymLabelCounts2_rot(:, :) = 0
        ! first 8 refer to the occupied, and the second to the virtual.

        allocate(LabVirtOrbs(NoOrbs - NoOcc), stat=ierr)
        call LogMemAlloc('LabVirtOrbs', (NoOrbs - NoOcc), 4, this_routine, LabVirtOrbsTag, ierr)
        LabVirtOrbs(:) = 0
        allocate(LabOccOrbs(NoOcc), stat=ierr)
        call LogMemAlloc('LabOccOrbs', (NoOcc), 4, this_routine, LabOccOrbsTag, ierr)
        LabOccOrbs(:) = 0
        allocate(SymVirtOrbs(NoOrbs - NoOcc), stat=ierr)
        call LogMemAlloc('SymVirtOrbs', (NoOrbs - NoOcc), 4, this_routine, SymVirtOrbsTag, ierr)
        SymVirtOrbs(:) = 0
        allocate(SymOccOrbs(NoOcc), stat=ierr)
        call LogMemAlloc('SymOccOrbs', (NoOcc), 4, this_routine, SymOccOrbsTag, ierr)
        SymOccOrbs(:) = 0

! First fill SymLabelList2_rot.

! 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 (tStoreSpinOrbs) then
                LabOccOrbs(i) = BRR(i)
                SymOccOrbs(i) = int(G1(LabOccOrbs(i))%sym%S)
            else
                LabOccOrbs(i) = (BRR(2 * i)) / 2
                SymOccOrbs(i) = int(G1(LabOccOrbs(i) * 2)%sym%S)
            end if
        end do

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

        do i = 1, NoOrbs - NoOcc
            if (tStoreSpinOrbs) then
                LabVirtOrbs(i) = BRR(i + NEl)
                SymVirtOrbs(i) = int(G1(LabVirtOrbs(i))%sym%S)
            else
                LabVirtOrbs(i) = (BRR((2 * i) + NEl)) / 2
                SymVirtOrbs(i) = int(G1(LabVirtOrbs(i) * 2)%sym%S)
            end if
        end do

        call sort(SymVirtOrbs, LabVirtOrbs)

! SymLabelList2_rot is then filled with the symmetry ordered occupied then virtual arrays.
        do i = 1, NoOcc
            SymLabelList2_rot(i) = LabOccOrbs(i)
        end do
        j = 0
        do i = NoOcc + 1, NoOrbs
            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.
            SymLabelCounts2_rot(1, 1) = 1
            SymLabelCounts2_rot(1, 9) = NoOcc + 1
            SymLabelCounts2_rot(2, 1) = NoOcc
            SymLabelCounts2_rot(2, 9) = NoOrbs - NoOcc
        else
            ! otherwise we run through the occupied orbitals, counting the number with each symmetry
            ! and noting where in SymLabelList2_rot each symmetry block starts.
            SymCurr = 0
            SymLabelCounts2_rot(1, 1) = 1
            do i = 1, NoOcc
                if (tStoreSpinOrbs) then
                    Symi = int(G1(SymLabelList2_rot(i))%sym%S)
                else
                    Symi = int(G1(SymLabelList2_rot(i) * 2)%sym%S)
                end if
                SymLabelCounts2_rot(2, (Symi + 1)) = SymLabelCounts2_rot(2, (Symi + 1)) + 1
                if (Symi > SymCurr) then
                    SymLabelCounts2_rot(1, (Symi + 1)) = i
                    SymCurr = Symi
                end if
            end do
            ! the same is then done for the virtuals.
            SymCurr = 0
            SymLabelCounts2_rot(1, 9) = NoOcc + 1
            do i = NoOcc + 1, NoOrbs
                if (tStoreSpinOrbs) then
                    Symi = int(G1(SymLabelList2_rot(i))%sym%S)
                else
                    Symi = int(G1(SymLabelList2_rot(i) * 2)%sym%S)
                end if
                SymLabelCounts2_rot(2, (Symi + 9)) = SymLabelCounts2_rot(2, (Symi + 9)) + 1
                if (Symi > SymCurr) then
                    SymLabelCounts2_rot(1, (Symi + 9)) = i
                    SymCurr = Symi
                end if
            end do
        end if

        ! Go through each symmetry group, making sure the orbital pairs are ordered lowest to highest.
        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

! Deallocate the arrays just used in this routine.
        deallocate(LabOccOrbs)
        call LogMemDealloc(this_routine, LabOccOrbsTag)
        deallocate(LabVirtOrbs)
        call LogMemDealloc(this_routine, LabVirtOrbsTag)
        deallocate(SymOccOrbs)
        call LogMemDealloc(this_routine, SymOccOrbsTag)
        deallocate(SymVirtOrbs)
        call LogMemDealloc(this_routine, SymVirtOrbsTag)

    end subroutine InitOrbitalSeparation