generate_cas Subroutine

    subroutine generate_cas(occ_orbs, virt_orbs, ilut_list, space_size)

        ! In: occ_orbs - The number of electrons in the CAS space.
        ! In: virt_orbs - The number of virtual spin orbitals in the CAS space.
        ! In/Out: ilut_list - List of determinants generated.
        ! In/Out: space_size - Number of determinants in the generated space.
        !             If ilut_list is not empty on input and you want to keep
        !             the states already in it, then on input space_size should
        !             be equal to the number of states to be kept in ilut_list,
        !             and new states will be added in from space_size+1.
        !             Otherwise, space_size must equal 0 on input.
        !             On output space_size will equal the total number of
        !             generated plus what space_size was on input.

        integer, intent(in) :: occ_orbs, virt_orbs
        integer(n_int), intent(inout) :: ilut_list(0:, :)
        integer, intent(inout) :: space_size

        type(BasisFN) :: CASSym
        integer(n_int) :: ilut(0:NIfTot)
        integer :: iCASDet
        integer :: num_active_orbs, nCASDet, i, j, comp, ierr
        integer, allocatable :: CASBrr(:), CASRef(:)
        integer(n_int) :: cas_bitmask(0:NIfD), cas_not_bitmask(0:NIfD)
        integer, pointer :: CASDets(:, :) => null()
        integer(TagIntType) :: CASDetsTag
        character(len=*), parameter :: t_r = "generate_cas"

        ! Start by adding the HF state.
        call add_state_to_space(ilutHF, ilut_list, space_size)

        ! This option should be true. It tells the subroutine gndts to only consider states
        ! with an Ms value in the correct spin subspace.
        if (.not. tSpn) call stop_all("generate_cas", "tSpn is not set to true.")

        ! The total number of orbitals in the active space:
        num_active_orbs = occ_orbs + virt_orbs
        allocate(CASBrr(1:num_active_orbs))
        allocate(CASRef(1:occ_orbs))
        do i = 1, num_active_orbs
            ! Run through the cas space, and create an array which will map these orbtials to the
            ! orbitals they actually represent.
            ! i.e. CASBRR(1) will store the lowest energy orbital in the CAS space and
            ! CASBRR(num_active_orbs) will store the highest energy orbital in the CAS space.
            CASBrr(i) = BRR(i + (nel - occ_orbs))
        end do

        ! Create a bit mask which has 1's in the bits which represent active orbitals and 0's in
        ! all other orbitals.
        cas_bitmask = 0_n_int
        do i = 1, num_active_orbs
            set_orb(cas_bitmask, CASBrr(i))
        end do
        ! Create a bit mask which has 0's in the bits which represent active orbitals and 1's in
        ! all other orbitals.
        cas_not_bitmask = not(cas_bitmask)

        ! CASRef holds the part of the HF determinant in the active space.
        CASRef = CasBRR(1:occ_orbs)
        call sort(CasRef)
        call GetSym(CASRef, occ_orbs, G1, nBasisMax, CASSym)

        ! First, generate all excitations so we know how many there are, to allocate the arrays.
        call gndts(occ_orbs, num_active_orbs, CASBrr, nBasisMax, CASDets, &
                   .true., G1, tSpn, LMS, .true., CASSym, nCASDet, iCASDet)

        if (nCASDet == 0) call stop_all("generate_cas", "No CAS determinants found.")

        ! Now allocate the array CASDets. CASDets(:,i) will store the orbitals in the active space
        ! which are occupied in the i'th determinant generated.
        allocate(CASDets(occ_orbs, nCASDet), stat=ierr)
        call LogMemAlloc("CASDets", occ_orbs * nCASDet, 8, t_r, CASDetsTag, ierr)
        CASDets(:, :) = 0

        ! Now fill up CASDets...
        call gndts(occ_orbs, num_active_orbs, CASBrr, nBasisMax, CASDets, &
                   .false., G1, tSpn, LMS, .true., CASSym, nCASDet, iCASDet)

        do i = 1, nCASDet
            ! First, create the bitstring representing this state:
            ! Start from the HF determinant and apply cas_not_bitmask to clear all active space
            ! orbitals.
            ilut(0:NIfD) = iand(ilutHF(0:NIfD), cas_not_bitmask)
            ! Then loop through the occupied orbitals in the active space, stored in CASDets(:,i),
            ! and set the corresponding bits.
            do j = 1, occ_orbs
                set_orb(ilut, CASDets(j, i))
            end do

            ! The function gndts always outputs the reference state. This is already included, so
            ! we want to cycle when we get to this state to ignore it.
            ! comp will be 0 if ilut and ilutHF are the same.
            comp = DetBitLT(ilut, ilutHF, NIfD)
            if (comp == 0) cycle

            ! Now that we have fully generated the determinant, add it to the main list.
            call add_state_to_space(ilut, ilut_list, space_size)

        end do

        deallocate(CASBrr)
        deallocate(CASRef)
        deallocate(CASDets, stat=ierr)
        call LogMemDealloc(t_r, CASDetsTag, ierr)

    end subroutine generate_cas