gasci_util.F90 Source File

!> This module contains functions for GAS that are not bound
!>  to a specific GAS excitation generator.
module gasci_util
    use constants, only: n_int, dp, int64
    use SystemData, only: nEl, nBasis
    use gasci, only: GASSpec_t, LocalGASSpec_t, CumulGASSpec_t, GAS_specification
    use gasci_supergroup_index, only: get_supergroups
    use orb_idx_mod, only: SpinProj_t, calc_spin_raw, sum, alpha, beta
    use sort_mod, only: sort
    use excitation_types, only: Excite_1_t, Excite_2_t, excite, get_last_tgt, &
        set_last_tgt, UNKNOWN, canonicalize
    use util_mod, only: lex_leq, cumsum, operator(.div.), near_zero, binary_search_first_ge, &
        operator(.isclose.), custom_findloc, choose_i64, swap
    use dSFMT_interface, only: genrand_real2_dSFMT
    use DetBitOps, only: ilut_lt, ilut_gt, EncodeBitDet
    use bit_rep_data, only: NIfTot, NIfD
    use sltcnd_mod, only: sltcnd_excit
    use sets_mod, only: disjoint, is_sorted
    use growing_buffers, only: buffer_int_2D_t, buffer_int_1D_t
    use bit_reps, only: decode_bit_det
    implicit none
    private
    public :: get_available_singles, get_available_doubles, &
        gen_all_excits, get_n_SDs


    public :: get_cumulative_list, draw_from_cum_list

    public :: write_GAS_info, t_output_GAS_sizes

    interface get_cumulative_list
            module procedure get_cumulative_list_Excite_1_t
            module procedure get_cumulative_list_Excite_2_t
    end interface

    logical :: t_output_GAS_sizes = .false.

contains


    !>  @brief
    !>  Get all single excitated determinants from det_I that are allowed under GAS constraints.
    pure function get_available_singles(GAS_spec, det_I) result(singles_exc_list)
        class(GASSpec_t), intent(in) :: GAS_spec
        integer, intent(in) :: det_I(:)
        integer, allocatable :: singles_exc_list(:, :)
            !! Dimension is (nEl, n_configurations)
        character(*), parameter :: this_routine = 'get_available_singles'

        integer, allocatable :: possible_holes(:)
        integer :: i, j, src1, tgt1
        type(SpinProj_t) :: m_s_1
        type(buffer_int_2D_t) :: buffer


#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (GAS_spec%contains_conf(det_I))) then
            call stop_all (this_routine, "Assert fail: GAS_spec%contains_conf(det_I)")
        end if
    end block
#endif

        call buffer%init(size(det_I))

        do i = 1, size(det_I)
            src1 = det_I(i)
            m_s_1 = calc_spin_raw(src1)
            possible_holes = GAS_spec%get_possible_holes(&
                        det_I, add_holes=det_I(i:i), &
                        excess=-m_s_1)
            do j = 1, size(possible_holes)
                tgt1 = possible_holes(j)
                call buffer%push_back(excite(det_I, Excite_1_t(src1, tgt1)))
            end do
        end do
        call buffer%dump_reset(singles_exc_list)

    

        ! merge_sort
        block
            integer, dimension(:, :), allocatable :: tmp
            integer :: current_size, left, mid, right
            integer, parameter :: along = 2, &
                bubblesort_size = 20

            associate(X => singles_exc_list)
                ! Determine good number of starting splits
                block
                    integer :: n_splits
                    n_splits = 1
                    do while (size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0) > bubblesort_size)
                        n_splits = n_splits + 1
                    end do
                    current_size = size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0)
                end block

                ! Reuse this variable as tmp for swap in bubble_sort
                ! and for merge in merge_sort.
                block
                    integer :: max_buffer_size, n_merges
                    n_merges = ceiling(log(real(size(X, along)) / real(current_size)) / log(2.0))
                    max_buffer_size = current_size * merge(2**(n_merges - 1), 1, n_merges >= 1)
                    allocate(tmp(size(X, 1), max_buffer_size))
                end block

                ! Bubble sort bins of size `bubblesort_size`.
                do left = lbound(X, along), ubound(X, along), current_size
                    right = min(left + bubblesort_size - 1, ubound(X, along))
                        ! bubblesort
    block
      integer :: n, i
      associate(V => X(:, left : right))
        do n = ubound(V, 2), lbound(V, 2) + 1, -1
          do i = lbound(V, 2), ubound(V, 2) - 1
            if (.not.   lex_leq(V(:, i), V(:, i + 1))) then
                  ! swap
    block
        associate(tmp => tmp(:, 1))
            tmp = V(:, i)
            V(:, i) = V(:, i + 1)
            V(:, i + 1) = tmp
        end associate
    end block
            end if
          end do
        end do
      end associate
    end block
                end do

                do while (current_size < size(X, along))
                    do left = lbound(X, along), ubound(X, along), 2 * current_size
                        right = min(left + 2 * current_size - 1, ubound(X, along))
                        mid = min(left + current_size, right) - 1
                        tmp(:, : mid - left + 1) = X(:, left : mid)
                            ! merge
    block
        integer :: i, j, k

        associate(A => tmp(:, : mid - left + 1), B => X(:, mid + 1 : right), C => X(:, left : right))

            if (size(A) + size(B) > size(C)) then
                error stop
            end if

            i = lbound(A, 2)
            j = lbound(B, 2)
            do k = lbound(C, 2), ubound(C, 2)
                if (i <= ubound(A, 2) .and. j <= ubound(B, 2)) then
                    if (  lex_leq(A(:, i), B(:, j))) then
                        C(:, k) = A(:, i)
                        i = i + 1
                    else
                        C(:, k) = B(:, j)
                        j = j + 1
                    end if
                else if (i <= ubound(A, 2)) then
                    C(:, k) = A(:, i)
                    i = i + 1
                else if (j <= ubound(B, 2)) then
                    C(:, k) = B(:, j)
                    j = j + 1
                end if
            end do
        end associate
    end block
                    end do
                    current_size = 2 * current_size
                end do
            end associate
        end block
    end function

    !>  @brief
    !>  Get all double excitated determinants from det_I that are allowed under GAS constraints.
    pure function get_available_doubles(GAS_spec, det_I) result(doubles_exc_list)
        class(GASSpec_t), intent(in) :: GAS_spec
        integer, intent(in) :: det_I(:)
        integer, allocatable :: doubles_exc_list(:, :)
        character(*), parameter :: this_routine = 'get_available_doubles'

        integer, allocatable :: first_pick_possible_holes(:), second_pick_possible_holes(:), deleted(:)
        integer :: i, j, k, l, src1, src2, tgt1, tgt2
        type(SpinProj_t) :: m_s_1
        type(buffer_int_2D_t) :: buffer

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (GAS_spec%contains_conf(det_I))) then
            call stop_all (this_routine, "Assert fail: GAS_spec%contains_conf(det_I)")
        end if
    end block
#endif

        call buffer%init(size(det_I))

        do i = 1, size(det_I)
            do j = i + 1, size(det_I)
                src1 = det_I(i)
                src2 = det_I(j)
                deleted = det_I([i, j])
                first_pick_possible_holes = GAS_spec%get_possible_holes(det_I, &
                                        add_holes=deleted, excess=-sum(calc_spin_raw(deleted)), &
                                        n_total=2)
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (disjoint(first_pick_possible_holes, det_I))) then
            call stop_all (this_routine, "Assert fail: disjoint(first_pick_possible_holes, det_I)")
        end if
    end block
#endif
                do k = 1, size(first_pick_possible_holes)
                    tgt1 = first_pick_possible_holes(k)
                    m_s_1 = calc_spin_raw(tgt1)
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (any(m_s_1 == [alpha, beta]))) then
            call stop_all (this_routine, "Assert fail: any(m_s_1 == [alpha, beta])")
        end if
    end block
#endif

                    second_pick_possible_holes = GAS_spec%get_possible_holes(&
                            det_I, add_holes=deleted, &
                            add_particles=[tgt1], &
                            n_total=1, excess=m_s_1 - sum(calc_spin_raw(deleted)))

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (disjoint(second_pick_possible_holes, [tgt1]))) then
            call stop_all (this_routine, "Assert fail: disjoint(second_pick_possible_holes, [tgt1])")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (disjoint(second_pick_possible_holes, det_I))) then
            call stop_all (this_routine, "Assert fail: disjoint(second_pick_possible_holes, det_I)")
        end if
    end block
#endif

                    do l = 1, size(second_pick_possible_holes)
                        tgt2 = second_pick_possible_holes(l)
                        call buffer%push_back(excite(det_I, canonicalize(Excite_2_t(src1, tgt1, src2, tgt2))))
                    end do
                end do
            end do
        end do

        call buffer%dump_reset(doubles_exc_list)

    

        ! merge_sort
        block
            integer, dimension(:, :), allocatable :: tmp
            integer :: current_size, left, mid, right
            integer, parameter :: along = 2, &
                bubblesort_size = 20

            associate(X => doubles_exc_list)
                ! Determine good number of starting splits
                block
                    integer :: n_splits
                    n_splits = 1
                    do while (size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0) > bubblesort_size)
                        n_splits = n_splits + 1
                    end do
                    current_size = size(X, along) / n_splits + merge(0, 1, mod(size(X, along), n_splits) == 0)
                end block

                ! Reuse this variable as tmp for swap in bubble_sort
                ! and for merge in merge_sort.
                block
                    integer :: max_buffer_size, n_merges
                    n_merges = ceiling(log(real(size(X, along)) / real(current_size)) / log(2.0))
                    max_buffer_size = current_size * merge(2**(n_merges - 1), 1, n_merges >= 1)
                    allocate(tmp(size(X, 1), max_buffer_size))
                end block

                ! Bubble sort bins of size `bubblesort_size`.
                do left = lbound(X, along), ubound(X, along), current_size
                    right = min(left + bubblesort_size - 1, ubound(X, along))
                        ! bubblesort
    block
      integer :: n, i
      associate(V => X(:, left : right))
        do n = ubound(V, 2), lbound(V, 2) + 1, -1
          do i = lbound(V, 2), ubound(V, 2) - 1
            if (.not.   lex_leq(V(:, i), V(:, i + 1))) then
                  ! swap
    block
        associate(tmp => tmp(:, 1))
            tmp = V(:, i)
            V(:, i) = V(:, i + 1)
            V(:, i + 1) = tmp
        end associate
    end block
            end if
          end do
        end do
      end associate
    end block
                end do

                do while (current_size < size(X, along))
                    do left = lbound(X, along), ubound(X, along), 2 * current_size
                        right = min(left + 2 * current_size - 1, ubound(X, along))
                        mid = min(left + current_size, right) - 1
                        tmp(:, : mid - left + 1) = X(:, left : mid)
                            ! merge
    block
        integer :: i, j, k

        associate(A => tmp(:, : mid - left + 1), B => X(:, mid + 1 : right), C => X(:, left : right))

            if (size(A) + size(B) > size(C)) then
                error stop
            end if

            i = lbound(A, 2)
            j = lbound(B, 2)
            do k = lbound(C, 2), ubound(C, 2)
                if (i <= ubound(A, 2) .and. j <= ubound(B, 2)) then
                    if (  lex_leq(A(:, i), B(:, j))) then
                        C(:, k) = A(:, i)
                        i = i + 1
                    else
                        C(:, k) = B(:, j)
                        j = j + 1
                    end if
                else if (i <= ubound(A, 2)) then
                    C(:, k) = A(:, i)
                    i = i + 1
                else if (j <= ubound(B, 2)) then
                    C(:, k) = B(:, j)
                    j = j + 1
                end if
            end do
        end associate
    end block
                    end do
                    current_size = 2 * current_size
                end do
            end associate
        end block

        remove_double_appearances : block
            integer, allocatable :: tmp_buffer(:, :)
            allocate(tmp_buffer(size(doubles_exc_list, 1), size(doubles_exc_list, 2)))
            j = 1
            tmp_buffer(:, j) = doubles_exc_list(:, 1)
            do i = 2, size(doubles_exc_list, 2)
                if (any(doubles_exc_list(:, i - 1) /= doubles_exc_list(:, i))) then
                    j = j + 1
                    tmp_buffer(:, j) = doubles_exc_list(:, i)
                end if
            end do
            doubles_exc_list = tmp_buffer(:, : j)
        end block remove_double_appearances
    end function


    !>  @brief
    !>  Get all excitated determinants from
    !>  det_I that are allowed under GAS constraints.
    !>
    !>  @param[in] GAS_spec GAS specification
    !>  @param[in] nI Starting determinant
    !>  @param[out] n_excits Number of determinants
    !>  @param[out] det_list Allocatable array of determinants in ilut format
    !>  @param[in] ic Optional input for excitation level (ic=1 => singles, ic=2 => doubles)
    !>      If ommited generate all.
    subroutine gen_all_excits(GAS_spec, nI, n_excits, det_list, ic)
        class(GASSpec_t), intent(in) :: GAS_spec
        integer, intent(in) :: nI(:)
        integer, intent(out) :: n_excits
        integer(n_int), intent(out), allocatable :: det_list(:,:)
        integer, intent(in), optional :: ic

        integer, allocatable :: singles(:, :), doubles(:, :)
        integer :: i, j

        if (present(ic)) then
            select case(ic)
            case(1)
                singles = get_available_singles(GAS_spec, nI)
                allocate(doubles(0, 0))
            case(2)
                allocate(singles(0, 0))
                doubles = get_available_doubles(GAS_spec, nI)
            end select
        else
            singles = get_available_singles(GAS_spec, nI)
            doubles = get_available_doubles(GAS_spec, nI)
        end if

        n_excits = size(singles, 2) + size(doubles, 2)
        allocate(det_list(0:niftot, n_excits))
        j = 1
        do i = 1, size(singles, 2)
            call EncodeBitDet(singles(:, i), det_list(:, j))
            j = j + 1
        end do

        do i = 1, size(doubles, 2)
            call EncodeBitDet(doubles(:, i), det_list(:, j))
            j = j + 1
        end do

        call sort(det_list, ilut_lt, ilut_gt)
    end subroutine gen_all_excits



    !>  @brief
    !>  Build up a cumulative list of matrix elements.
    !>
    !>  @details
    !>  Calculate the matrix elements for the possible excitations from det_I
    !>  to the possible holes using the incomplete defined excitation.
    !>
    !>  @param[in] det_I, Reference determinant in "nI-format".
    !>  @param[in] incomplete_exc, An excitation where the last target is unknown.
    !>  @param[in] possible_holes, Possible holes for the last target.
    function get_cumulative_list_Excite_1_t(det_I, incomplete_exc, possible_holes) result(cSum)
        integer, intent(in) :: det_I(:)
        type(Excite_1_t), intent(in) :: incomplete_exc
        integer, intent(in) :: possible_holes(:)
        real(dp) :: cSum(size(possible_holes))
        character(*), parameter :: this_routine = 'get_cumulative_list_Excite_1_t'

        real(dp) :: previous
        type(Excite_1_t) :: exc
        integer :: i

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (get_last_tgt(exc) == UNKNOWN)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion get_last_tgt(exc) == UNKNOWN"
            write(stderr, *) "failed in /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/gasci_uti&
                &l.fpp:224"
            call stop_all (this_routine, "Assert fail: get_last_tgt(exc) == UNKNOWN")
        end if
    end block
#endif
        exc = incomplete_exc

        ! build the cumulative list of matrix elements <src|H|tgt>
        previous = 0.0_dp
        do i = 1, size(possible_holes)
            call set_last_tgt(exc, possible_holes(i))
            cSum(i) = abs(sltcnd_excit(det_I, canonicalize(exc), .false.)) + previous
            previous = cSum(i)
        end do

        ! Normalize
        if (near_zero(cSum(size(cSum)))) then
            cSum(:) = 0.0_dp
        else
            cSum(:) = cSum(:) / cSum(size(cSum))
        end if
    end function get_cumulative_list_Excite_1_t
    !>  @brief
    !>  Build up a cumulative list of matrix elements.
    !>
    !>  @details
    !>  Calculate the matrix elements for the possible excitations from det_I
    !>  to the possible holes using the incomplete defined excitation.
    !>
    !>  @param[in] det_I, Reference determinant in "nI-format".
    !>  @param[in] incomplete_exc, An excitation where the last target is unknown.
    !>  @param[in] possible_holes, Possible holes for the last target.
    function get_cumulative_list_Excite_2_t(det_I, incomplete_exc, possible_holes) result(cSum)
        integer, intent(in) :: det_I(:)
        type(Excite_2_t), intent(in) :: incomplete_exc
        integer, intent(in) :: possible_holes(:)
        real(dp) :: cSum(size(possible_holes))
        character(*), parameter :: this_routine = 'get_cumulative_list_Excite_2_t'

        real(dp) :: previous
        type(Excite_2_t) :: exc
        integer :: i

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        use constants, only: stderr
        if (.not. (get_last_tgt(exc) == UNKNOWN)) then
            write(stderr, *) ""
            write(stderr, *) "Assertion get_last_tgt(exc) == UNKNOWN"
            write(stderr, *) "failed in /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/gasci_uti&
                &l.fpp:224"
            call stop_all (this_routine, "Assert fail: get_last_tgt(exc) == UNKNOWN")
        end if
    end block
#endif
        exc = incomplete_exc

        ! build the cumulative list of matrix elements <src|H|tgt>
        previous = 0.0_dp
        do i = 1, size(possible_holes)
            call set_last_tgt(exc, possible_holes(i))
            cSum(i) = abs(sltcnd_excit(det_I, canonicalize(exc), .false.)) + previous
            previous = cSum(i)
        end do

        ! Normalize
        if (near_zero(cSum(size(cSum)))) then
            cSum(:) = 0.0_dp
        else
            cSum(:) = cSum(:) / cSum(size(cSum))
        end if
    end function get_cumulative_list_Excite_2_t


    !>  @brief
    !>  Draw from a cumulative list.
    subroutine draw_from_cum_list(c_sum, idx, pgen)
        real(dp), intent(in) :: c_sum(:)
        integer, intent(out) :: idx
        real(dp), intent(out) :: pgen
        character(*), parameter :: this_routine = 'draw_from_cum_list'

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. ((c_sum(size(c_sum)) .isclose. 0.0_dp)             .or. (c_sum(size(c_sum)) .isclose. 1.0_dp))) then
            call stop_all (this_routine, "Assert fail: (c_sum(size(c_sum)) .isclose. 0.0_dp)             .or. (c_sum(size(c_sum))&
                & .isclose. 1.0_dp)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(c_sum))) then
            call stop_all (this_routine, "Assert fail: is_sorted(c_sum)")
        end if
    end block
#endif

        ! there might not be such an excitation
        if (c_sum(size(c_sum)) > 0) then
            ! find the index of the target hole in the cumulative list
            idx = binary_search_first_ge(c_sum, genrand_real2_dSFMT())
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (1 <= idx .and. idx <= size(c_sum))) then
            call stop_all (this_routine, "Assert fail: 1 <= idx .and. idx <= size(c_sum)")
        end if
    end block
#endif

            ! adjust pgen with the probability for picking tgt from the cumulative list
            if (idx == 1) then
                pgen = c_sum(1)
            else
                pgen = c_sum(idx) - c_sum(idx - 1)
            end if
        else
            idx = 0
            pgen = 0.0_dp
        end if
    end subroutine

    pure function get_alpha_supergroups(sg, n_orbs, S_z) result(sg_alpha)
        !! Return the possible supergroups/distributions for alpha electrons.
        !!
        !! If `sg_alpha` and `sg_beta` are the distributions of alpha/beta electrons among the
        !! GAS spaces. Then we have `sg(:) = sg_alpha(:) + sg_beta(:)`.
        !! We want to generate all possible `sg_alpha` such that
        !! The GAS constraints are still fullfilled and the total
        !! spin projection is maintained. (`sum(sg_alpha) + sum(sg_beta) == 2*S_z`)
        integer, intent(in) :: sg(:)
            !! The overall supergroup
        integer, intent(in) :: n_orbs(size(sg))
            !! The number of spatial orbitals per GAS space
        type(SpinProj_t), intent(in) :: S_z
            !! The Spin projection
        integer, allocatable :: sg_alpha(:, :)
            !! All possible distributions of alpha electrons among the
            !!      GAS spaces.
        integer :: N, i, j

        N = sum(sg)
        if (mod(S_z%val + N, 2) == 0) then
        block
            type(LocalGASSpec_t) :: sg_alpha_constraint
            integer :: N_alpha
            N_alpha = (N + S_z%val) .div. 2
            sg_alpha_constraint = LocalGASSpec_t(&
                                n_min=max(sg - n_orbs, 0), &
                                n_max=min(sg, n_orbs), &
                                spat_GAS_orbs=[((j, i = 1, n_orbs(j)), j = 1, size(n_orbs))])
            sg_alpha = get_supergroups(sg_alpha_constraint, N_alpha)
        end block
        else
            allocate(sg_alpha(0, 0))
        end if
    end function

    elemental function get_n_SDs(GAS_spec, N, S_z) result(n_SDs)
        !! Return the number of Slater-determinants.
        class(GASSpec_t), intent(in) :: GAS_spec
            !! GAS specification.
        integer, intent(in) :: N
            !! The number of particles
        type(SpinProj_t), intent(in) :: S_z
            !! Spin projection
        integer(int64) :: n_SDs
        integer, allocatable :: supergroups(:, :), alpha_supergroups(:, :)
        integer :: i, j, iGAS

        supergroups = get_supergroups(GAS_spec, N)

        block
            integer :: N_alpha(GAS_spec%nGAS()), N_beta(GAS_spec%nGAS()), n_spat_orbs(GAS_spec%nGAS())
                !! These are the number of alpha/beta electrons and the number of
                !!  spatial orbitals per GAS space
            n_spat_orbs = GAS_spec%GAS_size() .div. 2

            n_SDs = 0_int64
            do j = 1, size(supergroups, 2)
                alpha_supergroups = get_alpha_supergroups(supergroups(:, j), n_spat_orbs, S_z)
                do i = 1, size(alpha_supergroups, 2)
                    N_alpha = alpha_supergroups(:, i)
                    N_beta = supergroups(:, j) - N_alpha(:)
                    ! For a given supergroup and S_z in each GAS space
                    !   (coming from the distribution of alpha electrons)
                    !   the number of possible configurations in each GAS space is calculated.
                    ! The overall number is just the product over the GAS spaces.
                    n_SDs = n_SDs + product([(choose_i64(n_spat_orbs(iGAS), N_alpha(iGAS)) &
                                            * choose_i64(n_spat_orbs(iGAS), N_beta(iGAS)), iGAS = 1, GAS_spec%nGAS())])
                end do
            end do
        end block
    end function

    subroutine write_GAS_info(GAS_spec, N, S_z, iunit)
        !! Write info about the GAS constraints to `iunit`
        !!
        !! The routine especially compares the CAS and GAS Hilbert space sizes.
        class(GASSpec_t), intent(in) :: GAS_spec
            !! GAS constraints.
        integer, intent(in) :: N
            !! The particle number.
        type(SpinProj_t), intent(in) :: S_z
            !! The total spin projection.
        integer, intent(in) :: iunit

        call GAS_spec%write_to(iunit)

#ifdef WARNING_WORKAROUND_
        associate(N => N); end associate
        associate(S_z => S_z); end associate
#endif

        if (t_output_GAS_sizes) then
        block
            integer :: n_alpha, n_beta, n_spat_orbs
            integer(int64) :: size_CAS, size_GAS
            N_alpha = (N + S_z%val) .div. 2
            N_beta = N - N_alpha
            n_spat_orbs = GAS_spec%n_spin_orbs() .div. 2
            size_CAS = choose_i64(n_spat_orbs, N_alpha) * choose_i64(n_spat_orbs, N_beta)
            write(iunit, '(A, 1x, I0)') 'The size of the CAS space is:', size_CAS
            size_GAS = get_n_SDs(GAS_spec, nEl, S_z)
            write(iunit, '(A, 1x, I0)') 'The size of the GAS space is:', size_GAS
            write(iunit, '(A, 1x, E10.5)') 'The fraction of the GAS space is:', real(size_GAS, dp) / real(size_CAS, dp)
        end block
        end if
    end subroutine
end module