gasci.F90 Source File


Contents

Source Code


Source Code

#include "macros.h"


module gasci
    use constants, only: n_int, dp
    use SystemData, only: nBasis
    use util_mod, only: cumsum, stop_all, operator(.div.)
    use excitation_types, only: Excite_1_t, Excite_2_t
    use util_mod, only: lex_leq, cumsum, operator(.div.), near_zero, binary_search_first_ge, &
        operator(.isclose.), custom_findloc, EnumBase_t, lex_geq
    use sets_mod, only: disjoint, operator(.U.), is_sorted, operator(.complement.)
    use orb_idx_mod, only: SpinProj_t, calc_spin_raw, sum, alpha, beta
    use excitation_types, only: Excite_1_t, Excite_2_t, excite, get_last_tgt, set_last_tgt, UNKNOWN
    use sltcnd_mod, only: sltcnd_excit
    use dSFMT_interface, only: genrand_real2_dSFMT
    use bit_rep_data, only: nIfTot
    use bit_reps, only: decode_bit_det
    use growing_buffers, only: buffer_int_2D_t, buffer_int_1D_t

    better_implicit_none
    private
    public :: possible_GAS_exc_gen, &
        GAS_exc_gen, GAS_specification, GASSpec_t, &
        user_input_GAS_exc_gen, get_name, &
        LocalGASSpec_t, CumulGASSpec_t, FlexibleGASSpec_t

    type, extends(EnumBase_t) :: GAS_exc_gen_t
    end type

    type :: possible_GAS_exc_gen_t
        type(GAS_exc_gen_t) :: &
            DISCONNECTED = GAS_exc_gen_t(1), &
            ON_FLY_HEAT_BATH = GAS_exc_gen_t(2), &
            DISCARDING = GAS_exc_gen_t(3), &
            PCHB = GAS_exc_gen_t(4)
    end type

    type(possible_GAS_exc_gen_t), parameter :: possible_GAS_exc_gen = possible_GAS_exc_gen_t()

    type(GAS_exc_gen_t) :: GAS_exc_gen = possible_GAS_exc_gen%ON_FLY_HEAT_BATH
    type(GAS_exc_gen_t), allocatable :: user_input_GAS_exc_gen

    ! NOTE: At the current state of implementation `GASSpec_t` is a completely immutable
    ! datastructure to outside code after the constructor has been called.
    ! If possible keep it like this!

    type, abstract :: GASSpec_t
        !! Speficies the GAS spaces.
        private
        integer, allocatable :: GAS_table(:)
            !! GAS_table(i) returns the GAS space for the i-th spin orbital
        integer, allocatable  :: GAS_sizes(:)
            !! The number of spin orbitals per GAS space
        integer :: largest_GAS_size
            !! maxval(GAS_sizes)
        integer, allocatable :: splitted_orbitals(:, :)
            !! This is the preimage of `GASSpec_t%GAS_table.`
            !!
            !! An array that contains the spin orbitals per GAS space of dimension
            !! `splitted_orbitals(1 : maxval(GAS_sizes), 1 : nGAS)`.
            !! Only `splitted_orbitals(i, j), 1 <= i <= GAS_sizes(j)`
            !! is defined.
        logical :: lookup_is_connected
            !! These lookup variables stay valid, because the data structure is
            !!  immutable
        logical :: exchange_recoupling
    contains
        ! Nearly all member functions should be public.
        procedure(contains_supergroup_t), deferred :: contains_supergroup
        procedure(is_valid_t), deferred :: is_valid
        procedure(write_to_t), deferred :: write_to
        procedure(get_possible_spaces_t), deferred :: get_possible_spaces
        procedure :: get_possible_holes
        procedure :: contains_conf
        procedure :: contains_ilut
        procedure :: is_connected => get_is_connected
        procedure :: nGAS => get_nGAS
        procedure :: n_spin_orbs => get_nOrbs
        procedure :: max_GAS_size => get_max_GAS_size
        procedure :: recoupling

        generic :: GAS_size => get_GAS_size_i, get_GAS_size_idx, get_GAS_size_all
        procedure, private :: get_GAS_size_i
        procedure, private :: get_GAS_size_idx
        procedure, private :: get_GAS_size_all

        procedure :: get_iGAS
        procedure :: get_orb_idx
        procedure :: count_per_GAS
        generic :: is_allowed => is_allowed_single, is_allowed_double
        procedure, private :: is_allowed_single
        procedure, private :: is_allowed_double
    end type

    abstract interface
        logical pure function contains_supergroup_t(this, supergroup)
            import :: GASSpec_t
            implicit none
            class(GASSpec_t), intent(in) :: this
            integer, intent(in) :: supergroup(:)
        end function

        logical pure function is_valid_t(this, n_basis)
            import :: GASSpec_t
            implicit none
            class(GASSpec_t), intent(in) :: this
            integer, intent(in), optional :: n_basis
        end function

        subroutine write_to_t(this, iunit)
            !! Write a string representation of this GAS specification to iunit
            import :: GASSpec_t
            implicit none
            class(GASSpec_t), intent(in) :: this
            integer, intent(in) :: iunit
        end subroutine

        pure function get_possible_spaces_t(this, supergroup, add_holes, add_particles, n_total) result(spaces)
            !!  Return the GAS spaces, where one particle can be created.
            !!
            !!  The returned array can be empty (allocated, but size == 0).
            import :: GASSpec_t
            implicit none
            class(GASSpec_t), intent(in) :: this
                !!  Specification of GAS spaces.
            integer, intent(in) :: supergroup(size(this%GAS_sizes))
                !!  The particles per GAS space.
            integer, intent(in), optional :: add_holes(:)
                !!  An index of orbitals where particles should be deleted
                !!  before creating the new particle.
            integer, intent(in), optional :: add_particles(:)
                !!  Index of orbitals where particles should be created
                !! before creating the new particle.
            integer, intent(in), optional :: n_total
                !! The total number of particles that will be created.
                !! Defaults to one.  (Relevant for double excitations)
            integer, allocatable :: spaces(:)
        end function
    end interface

    type, extends(GASSpec_t) :: LocalGASSpec_t
        private
        !> The indices are:
        !> `min(1 : nGAS), max(1 : nGAS)`.
        !> `min(iGAS)` specifies the minimum particle number per GAS space.
        !> `max(iGAS)` specifies the maximum particle number per GAS space.
        integer, allocatable :: min(:), max(:)
    contains
        procedure :: contains_supergroup => Local_contains_supergroup
        procedure :: is_valid => Local_is_valid
        procedure :: write_to => Local_write_to
        procedure :: get_possible_spaces => Local_get_possible_spaces

        generic :: get_min => get_min_i, get_min_all
        procedure, private :: get_min_i, get_min_all
        generic :: get_max => get_max_i, get_max_all
        procedure, private :: get_max_i, get_max_all
    end type

    interface LocalGASSpec_t
        module procedure construct_LocalGASSpec_t
    end interface

    type, extends(GASSpec_t) :: CumulGASSpec_t
        private
        !> The indices are:
        !> `c_min(1 : nGAS), c_max(1 : nGAS)`.
        !> `c_min(iGAS)` specifies the cumulated minimum particle number per GAS space.
        !> `c_max(iGAS)` specifies the cumulated maximum particle number per GAS space.
        integer, allocatable :: c_min(:), c_max(:)
    contains
        procedure :: contains_supergroup => Cumul_contains_supergroup
        procedure :: is_valid => Cumul_is_valid
        procedure :: write_to => Cumul_write_to
        procedure :: get_possible_spaces => Cumul_get_possible_spaces

        generic :: get_cmin => get_cmin_i, get_cmin_all
        procedure, private :: get_cmin_i
        procedure, private :: get_cmin_all
        generic :: get_cmax => get_cmax_i, get_cmax_all
        procedure, private :: get_cmax_i
        procedure, private :: get_cmax_all
    end type

    interface CumulGASSpec_t
        module procedure construct_CumulGASSpec_t
    end interface

    type, extends(GASSpec_t) :: FlexibleGASSpec_t
        private
        !> List the allowed supergroups.
        !> The indices are: (nGAS, n_supergroups)
        integer, public, allocatable :: supergroups(:, :)
        !> The number of particles.
        integer :: N
    contains
        procedure :: contains_supergroup => Flexible_contains_supergroup
        procedure :: is_valid => Flexible_is_valid
        procedure :: write_to => Flexible_write_to
        procedure :: get_possible_spaces => Flexible_get_possible_spaces

        procedure :: N_particle => Flexible_N_particle
    end type

    interface FlexibleGASSpec_t
        module procedure construct_FlexibleGASSpec_t
    end interface

    class(GASSpec_t), allocatable :: GAS_specification

contains

    !> Returns the total number of GAS spaces.
    integer pure function get_nGAS(this)
        class(GASSpec_t), intent(in) :: this
        get_nGAS = size(this%GAS_sizes)
    end function

    !> Returns the size of the largest GAS space.
    integer pure function get_max_GAS_size(this)
        class(GASSpec_t), intent(in) :: this
        get_max_GAS_size = this%largest_GAS_size
    end function

    !>  Returns the size of the i-th GAS space in number of spin orbitals.
    integer pure function get_GAS_size_i(this, iGAS)
        class(GASSpec_t), intent(in) :: this
        integer, intent(in) :: iGAS
        get_GAS_size_i = this%GAS_sizes(iGAS)
    end function

    !>  Returns the sizes for GAS spaces specified in idx.
    pure function get_GAS_size_idx(this, idx) result(res)
        class(GASSpec_t), intent(in) :: this
        integer, intent(in) :: idx(:)
        integer :: res(size(idx))
        res(:) = this%GAS_sizes(idx)
    end function

    !>  Returns the sizes for all GAS spaces.
    pure function get_GAS_size_all(this) result(res)
        class(GASSpec_t), intent(in) :: this
        integer :: res(size(this%GAS_sizes))
        res(:) = this%GAS_sizes(:)
    end function

    !> Returns the GAS space for a given spin orbital index.
    integer elemental function get_iGAS(this, spin_orb_idx)
        class(GASSpec_t), intent(in) :: this
        integer, intent(in) :: spin_orb_idx
        get_iGAS = this%GAS_table(spin_orb_idx)
    end function

    integer elemental function get_nOrbs(this)
        class(GASSpec_t), intent(in) :: this
        get_nOrbs = size(this%GAS_table)
    end function

    !> Returns the i-th spin orbital in the iGAS GAS space.
    !>
    !> Can be seen as the preimage of get_iGAS (which is usually not injective).
    integer elemental function get_orb_idx(this, i, iGAS)
        class(GASSpec_t), intent(in) :: this
        integer, intent(in) :: i, iGAS
        character(*), parameter :: this_routine = 'get_orb_idx'

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (1 <= i .and. i <= this%GAS_size(iGAS))) then
            call stop_all (this_routine, "Assert fail: 1 <= i .and. i <= this%GAS_size(iGAS)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (1 <= iGAS .and. iGAS <= this%nGAS())) then
            call stop_all (this_routine, "Assert fail: 1 <= iGAS .and. iGAS <= this%nGAS()")
        end if
    end block
#endif

        get_orb_idx = this%splitted_orbitals(i, iGAS)
    end function


    logical pure function get_is_connected(this)
        !! Query if there are connected GAS spaces under the GAS specification.
        class(GASSpec_t), intent(in) :: this
            !!  Specification of GAS spaces.
        get_is_connected = this%lookup_is_connected
    end function


    !>  Query wether a determinant or CSF is contained in the GAS space.
    pure function contains_conf(this, nI) result(res)
        !> Specification of GAS spaces.
        class(GASSpec_t), intent(in) :: this
        !> An index of occupied spin orbitals.
        integer, intent(in) :: nI(:)
        logical :: res
        res = this%contains_supergroup(this%count_per_GAS(nI))
    end function

    logical elemental function recoupling(this)
        class(GASSpec_t), intent(in) :: this
        recoupling = this%exchange_recoupling
    end function

    !>  Query wether a determinant in bitmask format is contained in the GAS space.
    !>
    !>  The function in nI-format is faster!
    !>  It is **assumed** that the determinant is contained in the
    !>  Full CI space and obeys e.g. the Pauli principle.
    !>  The return value is not defined, if that is not the case!
    pure function contains_ilut(this, ilut) result(res)
        !>  Specification of GAS spaces.
        class(GASSpec_t), intent(in) :: this
        !>  An index of occupied spin orbitals.
        integer(n_int), intent(in) :: ilut(0 : nIfTot)
        logical :: res
        integer :: nI(sum(popcnt(ilut)))
        call decode_bit_det(nI, ilut)
        res = this%contains_conf(nI)
    end function

    !> Count the particles per GAS space. i.e. return the supergroup.
    pure function count_per_GAS(this, occupied) result(supergroup)
        class(GASSpec_t), intent(in) :: this
        integer, intent(in) :: occupied(:)

        integer :: supergroup(get_nGAS(this))

        integer :: iel, iGAS

        supergroup = 0
        do iel = 1, size(occupied)
            iGAS = this%get_iGAS(occupied(iel))
            supergroup(iGAS) = supergroup(iGAS) + 1
        end do
    end function

    pure function get_name(impl) result(res)
        type(GAS_exc_gen_t), intent(in) :: impl
        character(len=:), allocatable :: res
        if (impl == possible_GAS_exc_gen%DISCONNECTED) then
            res = 'Heat-bath on-the-fly GAS implementation for disconnected spaces'
        else if (impl == possible_GAS_exc_gen%ON_FLY_HEAT_BATH) then
            res = 'Heat-bath on-the-fly GAS implementation'
        else if (impl == possible_GAS_exc_gen%DISCARDING) then
            res = 'Discarding GAS implementation'
        else if (impl == possible_GAS_exc_gen%PCHB) then
            res = 'PCHB GAS implementation'
        end if
    end function

    !> Check if a single excitation is allowed.
    !>
    !> Is called once at initialization, so it does not have to be super fast.
    logical pure function is_allowed_single(this, exc, supergroup)
        class(GASSpec_t), intent(in) :: this
        type(Excite_1_t), intent(in) :: exc
        integer, intent(in) :: supergroup(:)

        integer :: excited_supergroup(size(supergroup))
        integer :: src_space, tgt_space

        src_space = this%get_iGAS(exc%val(1, 1))
        tgt_space = this%get_iGAS(exc%val(2, 1))

        if (src_space == tgt_space) then
            ! All electrons come from the same space and there are no restrictions
            ! regarding recoupling or GAS.
            is_allowed_single = .true.
        else
            ! Ensure that GAS specifications contain supergroup **after** excitation.
            excited_supergroup = supergroup
            excited_supergroup(src_space) = excited_supergroup(src_space) - 1
            excited_supergroup(tgt_space) = excited_supergroup(tgt_space) + 1

            is_allowed_single = this%contains_supergroup(excited_supergroup)
        end if
    end function


    !> Check if a double excitation is allowed.
    !>
    !> Is called once at initialization, so it does not have to be super fast.
    !> `recoupling` allows recoupling excitations that change the spin projection
    !> of individual GAS spaces.
    logical pure function is_allowed_double(this, exc, supergroup)
        class(GASSpec_t), intent(in) :: this
        type(Excite_2_t), intent(in) :: exc
        integer, intent(in) :: supergroup(:)

        integer :: excited_supergroup(size(supergroup))
        integer :: src_spaces(2), tgt_spaces(2)

        src_spaces = this%get_iGAS(exc%val(1, :))
        tgt_spaces = this%get_iGAS(exc%val(2, :))

        if (all(src_spaces == tgt_spaces) .and. src_spaces(1) == src_spaces(2)) then
            ! All electrons come from the same space and there are no restrictions
            ! regarding recoupling or GAS.
            is_allowed_double = .true.
        else
            ! Ensure that GAS specifications contain supergroup **after** excitation.
            excited_supergroup = supergroup
            excited_supergroup(src_spaces) = excited_supergroup(src_spaces) - 1
            excited_supergroup(tgt_spaces) = excited_supergroup(tgt_spaces) + 1

            is_allowed_double = this%contains_supergroup(excited_supergroup)

            if (is_allowed_double .and. .not. this%exchange_recoupling) then
                block
                    type(SpinProj_t) :: src_spins(2), tgt_spins(2)

                    src_spins = calc_spin_raw(exc%val(1, :))
                    tgt_spins = calc_spin_raw(exc%val(2, :))

                    if (src_spaces(1) > src_spaces(2)) then
    ! swap
    block
            type(SpinProj_t), allocatable :: tmp
            tmp = src_spins(1)
            src_spins(1) = src_spins(2)
            src_spins(2) = tmp
    end block
                    end if

                    if (tgt_spaces(1) > tgt_spaces(2)) then
    ! swap
    block
            type(SpinProj_t), allocatable :: tmp
            tmp = tgt_spins(1)
            tgt_spins(1) = tgt_spins(2)
            tgt_spins(2) = tmp
    end block
                    end if

                    is_allowed_double = all(src_spins == tgt_spins)
                end block
            end if
        end if
    end function

    !>  Return the possible holes where a particle can be created under GAS constraints.
    !>
    !>  This function uses `get_possible_spaces` to find possible GAS spaces
    !>  where a particle can be created and returns only unoccupied
    !>  sites of correct spin.
    !>
    !>  "Trivial" excitations are avoided. That means, that a site is only counted
    !>  as unoccupied if it was unoccupied in nI from the beginning on.
    !>  (A double excitation where a particle is deleted, but immediately
    !>  recreated would be such a trivial excitations.)
    pure function get_possible_holes(this, det_I, add_holes, add_particles, n_total, excess) result(possible_holes)
        !>  Specification of GAS spaces.
        class(GASSpec_t), intent(in) :: this
        !>  The starting determinant
        integer, intent(in) :: det_I(:)
        !>  An index of orbitals
        !>      where particles should be deleted before creating the new particle.
        integer, intent(in), optional :: add_holes(:)
        !>  An index of orbitals
        !>      where particles should be created before creating the new particle.
        integer, intent(in), optional :: add_particles(:)
        !>  The total number of particles
        !>      that will be created. Defaults to one.
        integer, intent(in), optional :: n_total
        !>  The current excess of spin projections.
        !>      If a beta electron was deleted, the excess is \(1 \cdot \alpha)\.
        type(SpinProj_t), intent(in), optional :: excess
        character(*), parameter :: this_routine = 'get_possible_holes'

        integer, allocatable :: possible_holes(:)

        integer, allocatable :: spaces(:)
        integer :: n_total_

if(present(n_total)) then
    n_total_ = n_total
else
    n_total_ = 1
endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (1 <= n_total_)) then
            call stop_all (this_routine, "Assert fail: 1 <= n_total_")
        end if
    end block
#endif
        if (present(excess)) then
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (abs(excess%val) <= n_total_)) then
            call stop_all (this_routine, "Assert fail: abs(excess%val) <= n_total_")
        end if
    end block
#endif
        end if


        ! Note, that non-present optional arguments can be passed
        ! into optional arguments without checking!
        spaces = this%get_possible_spaces(&
             this%count_per_GAS(det_I), add_holes=add_holes, &
             add_particles=add_particles, n_total=n_total_)

        if (size(spaces) == 0) then
            possible_holes = [integer::]
            return
        end if

        block
            integer :: i, iGAS, incr, curr_value, iGAS_min_val, idx_space
            integer :: L(size(spaces)), counter(size(spaces))
            integer, allocatable :: possible_values(:)
            type(SpinProj_t) :: m_s

            m_s = SpinProj_t(0)
            if (present(excess)) then
                if (abs(excess%val) == n_total_) m_s = -SpinProj_t(sign(1, excess%val))
            end if


            L(:) = this%GAS_size(spaces)
            if (m_s == beta) then
                allocate(possible_values(sum(L) .div. 2))
                counter(:) = 1
                incr = 2
            else if (m_s == alpha) then
                allocate(possible_values(sum(L) .div. 2))
                counter(:) = 2
                incr = 2
            else
                allocate(possible_values(sum(L)))
                counter(:) = 1
                incr = 1
            end if

            ! Here we merge the values from splitted_orbitals sortedly into
            ! possible_values
            i = 1
            do while (any(counter <= L))
                curr_value = huge(curr_value)
                ! foreach iGAS in spaces
                do idx_space = 1, size(spaces)
                    iGAS = spaces(idx_space)
                    if (counter(idx_space) <= L(idx_space)) then
                        if (this%get_orb_idx(counter(idx_space), iGAS) < curr_value) then
                            curr_value = this%get_orb_idx(counter(idx_space), iGAS)
                            iGAS_min_val = idx_space
                        end if
                    end if
                end do

                counter(iGAS_min_val) = counter(iGAS_min_val) + incr
                possible_values(i) = curr_value
                i = i + 1
            end do

            if (present(add_particles)) then
                possible_holes = possible_values .complement. (det_I .U. add_particles)
            else
                possible_holes = possible_values .complement. det_I
            end if
        end block
    end function

    !>  Constructor of LocalGASSpec_t
    pure function construct_LocalGASSpec_t(n_min, n_max, spat_GAS_orbs, recoupling) result(GAS_spec)
        !> Minimum particle number per GAS space.
        integer, intent(in) :: n_min(:)
        !> Maximum particle number per GAS space.
        integer, intent(in) :: n_max(:)
        !> GAS space for the i-th **spatial** orbital.
        integer, intent(in) :: spat_GAS_orbs(:)
        !> Exchange double excitations that recouple the spin are allowed
        logical, intent(in), optional :: recoupling
        logical :: recoupling_

        type(LocalGASSpec_t) :: GAS_spec
        character(*), parameter :: this_routine = 'construct_LocalGASSpec_t'

        integer :: n_spin_orbs, max_GAS_size
        integer, allocatable :: splitted_orbitals(:, :), GAS_table(:), GAS_sizes(:)
        integer :: i, iel, iGAS, nGAS

if(present(recoupling)) then
    recoupling_ = recoupling
else
    recoupling_ = .true.
endif

        nGAS = maxval(spat_GAS_orbs)
        GAS_sizes = 2 * frequency(spat_GAS_orbs)

        if (any(min(n_max, GAS_sizes) < n_min)) then
            call stop_all(this_routine, 'any(min(n_max, GAS_sizes) < n_min) violated')
        end if


        max_GAS_size = maxval(GAS_sizes)
        n_spin_orbs = sum(GAS_sizes)

        allocate(GAS_table(n_spin_orbs))
        GAS_table(1::2) = spat_GAS_orbs
        GAS_table(2::2) = spat_GAS_orbs

        allocate(splitted_orbitals(max_GAS_size, nGAS))
        block
            integer :: all_orbs(n_spin_orbs), splitted_sizes(nGAS)
            all_orbs = [(i, i = 1, n_spin_orbs)]
            splitted_sizes = 0
            do iel = 1, size(all_orbs)
                iGAS = GAS_table(iel)
                splitted_sizes(iGAS) = splitted_sizes(iGAS) + 1
                splitted_orbitals(splitted_sizes(iGAS), iGAS) = all_orbs(iel)
            end do
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (all(GAS_sizes == splitted_sizes))) then
            call stop_all (this_routine, "Assert fail: all(GAS_sizes == splitted_sizes)")
        end if
    end block
#endif
        end block


        ! This block and the additional declaration of new_n_max
        ! is just necessary because of shitty intel compilers (Ifort 18).
        ! If you can remove it, I am happy.
        block
            integer :: new_n_max(size(n_max))
            new_n_max = min(n_max, GAS_sizes)
            GAS_spec = LocalGASSpec_t(&
                    min=n_min, max=new_n_max, GAS_table=GAS_table, &
                    GAS_sizes=GAS_sizes, largest_GAS_size=max_GAS_size, &
                    splitted_orbitals=splitted_orbitals, &
                    lookup_is_connected=any(n_min(:) /= n_max(:)), &
                    exchange_recoupling=recoupling_)
        end block

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

        contains

        pure function frequency(N) result(res)
            integer, intent(in) :: N(:)
            integer, allocatable :: res(:)
            integer :: i
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (minval(N) == 1)) then
            call stop_all (this_routine, "Assert fail: minval(N) == 1")
        end if
    end block
#endif
            allocate(res(maxval(N)), source=0)
            do i = 1, size(N)
                res(N(i)) = res(N(i)) + 1
            end do
        end function
    end function

    !>  Query wether a supergroup is contained in the GAS space.
    pure function Local_contains_supergroup(this, supergroup) result(res)
        !> Specification of GAS spaces.
        class(LocalGASSpec_t), intent(in) :: this
        !> A supergroup.
        integer, intent(in) :: supergroup(:)
        logical :: res
        res = all(this%min(:) <= supergroup &
                  .and. supergroup <= min(this%max(:), this%GAS_sizes(:)))
    end function

    !> Check if the GAS specification is valid
    !>
    !> If the number of particles or the number of spin orbitals
    !> is provided, then the consistency with these numbers
    !> is checked as well.
    logical pure function Local_is_valid(this, n_basis)
        !>  Specification of GAS spaces with local constraints.
        class(LocalGASSpec_t), intent(in) :: this
        integer, intent(in), optional :: n_basis

        logical :: shapes_match, pauli_principle, n_orbs_correct

        associate(GAS_sizes => this%GAS_sizes, n_min => this%min, &
                  n_max => this%max, nGAS => this%nGAS())

            shapes_match = &
                all([size(GAS_sizes), size(n_min), size(n_max)] == nGAS) &
                .and. maxval(this%GAS_table) == nGAS

            pauli_principle = all(n_min(:) <= GAS_sizes)

            if (present(n_basis)) then
                n_orbs_correct = sum(GAS_sizes) == n_basis
            else
                n_orbs_correct = .true.
            end if
        end associate

        Local_is_valid = all([shapes_match, pauli_principle, n_orbs_correct])
    end function

    subroutine Local_write_to(this, iunit)
        class(LocalGASSpec_t), intent(in) :: this
        integer, intent(in) :: iunit
        integer :: iGAS, iorb

        write(iunit, '(A)') 'Local constraints'
        write(iunit, '(A)') 'n_i: number of spatial orbitals per i-th GAS space'
        write(iunit, '(A)') 'n_min_i: minimum of particle number per i-th GAS space'
        write(iunit, '(A)') 'n_max_i: maximum of particle number per i-th GAS space'
        write(iunit, '(A10, 1x, A10, 1x, A10)') 'n_i', 'n_min_i', 'n_max_i'
        write(iunit, '(A)') '--------------------------------'
        do iGAS = 1, this%nGAS()
            write(iunit, '(I10, 1x, I10, 1x, I10)') this%GAS_size(iGAS) .div. 2, this%get_min(iGAS), this%get_max(iGAS)
        end do
        write(iunit, '(A)') '--------------------------------'
        write(iunit, '(A)') 'The distribution of spatial orbitals to GAS spaces is given by:'
        do iorb = 1, this%n_spin_orbs(), 2
            write(iunit, '(I0, 1x)', advance='no') this%get_iGAS(iorb)
        end do
        write(iunit, *)

        if (any(this%GAS_sizes < this%max)) then
            write(iunit, '(A)') 'In at least one GAS space, the maximum allowed particle number by GAS constraints'
            write(iunit, '(A)') '   is larger than the particle number allowed by the Pauli principle.'
            write(iunit, '(A)') '   Was this intended when preparing your input?'
        end if

        if (.not. this%recoupling()) then
            write(iunit, '(A)') 'Double excitations with exchange are forbidden.'
        end if
    end subroutine


    !> Returns the minimum particle number for a given GAS space.
    integer elemental function get_min_i(this, iGAS)
        class(LocalGASSpec_t), intent(in) :: this
        integer, intent(in) :: iGAS
        get_min_i = this%min(iGAS)
    end function

    !> Returns the minimum particle number for all GAS spaces.
    pure function get_min_all(this) result(res)
        class(LocalGASSpec_t), intent(in) :: this
        integer, allocatable :: res(:)
        res = this%min(:)
    end function

    !> Returns the maximum particle number for a given GAS space.
    integer elemental function get_max_i(this, iGAS)
        class(LocalGASSpec_t), intent(in) :: this
        integer, intent(in) :: iGAS
        get_max_i = this%max(iGAS)
    end function

    !> Returns the maximum particle number for all GAS spaces.
    pure function get_max_all(this) result(res)
        class(LocalGASSpec_t), intent(in) :: this
        integer, allocatable :: res(:)
        res = this%max(:)
    end function

    pure function Local_get_possible_spaces(this, supergroup, add_holes, add_particles, n_total) result(spaces)
        class(LocalGASSpec_t), intent(in) :: this
        integer, intent(in) :: supergroup(size(this%GAS_sizes))
        integer, intent(in), optional :: add_holes(:), add_particles(:), n_total
        integer, allocatable :: spaces(:)

        integer :: n_total_, iGAS
        integer :: &
        !> pumber of particles per iGAS
            n_particle(this%nGAS()), &
        !> deficit per iGAS
            deficit(this%nGAS()), &
        !> vacant orbitals per iGAS
            vacant(this%nGAS())
        type(buffer_int_1D_t) :: space_buffer

if(present(n_total)) then
    n_total_ = n_total
else
    n_total_ = 1
endif

        block
            integer :: B(this%nGAS()), C(this%nGAS())
            if (present(add_holes)) then
                B = this%count_per_GAS(add_holes)
            else
                B = 0
            end if
            if (present(add_particles)) then
                C = this%count_per_GAS(add_particles)
            else
                C = 0
            end if
            n_particle = supergroup - B + C
        end block

        deficit(:) = this%get_min() - n_particle(:)
        vacant(:) = this%get_max() - n_particle(:)

        if (n_total_ < sum(deficit) .or. sum(vacant) < n_total_) then
            spaces = [integer::]
        else if (any(deficit == n_total_)) then
            spaces = [custom_findloc(deficit == n_total_, val=.true.)]
        else
            call space_buffer%init()
            do iGAS = 1, this%nGAS()
                if (vacant(iGAS) > 0) then
                    call space_buffer%push_back(iGAS)
                end if
            end do
            call space_buffer%dump_reset(spaces)
        end if
    end function


    !>  Constructor of CumulGASSpec_t
    pure function construct_CumulGASSpec_t(cn_min, cn_max, spat_GAS_orbs, recoupling) result(GAS_spec)
        !>  Cumulative minimum particle number.
        integer, intent(in) :: cn_min(:)
        !>  Cumulative maximum particle number.
        integer, intent(in) :: cn_max(:)
        !>  GAS space for the i-th **spatial** orbital.
        integer, intent(in) :: spat_GAS_orbs(:)
        !> Exchange double excitations that recouple the spin are allowed
        logical, intent(in), optional :: recoupling
        logical :: recoupling_

        type(CumulGASSpec_t) :: GAS_spec
        character(*), parameter :: this_routine = 'construct_CumulGASSpec_t'

        integer :: n_spin_orbs, max_GAS_size
        integer, allocatable :: splitted_orbitals(:, :), GAS_table(:), GAS_sizes(:)
        integer :: i, iel, iGAS, nGAS

if(present(recoupling)) then
    recoupling_ = recoupling
else
    recoupling_ = .true.
endif

        nGAS = maxval(spat_GAS_orbs)
        GAS_sizes = 2 * frequency(spat_GAS_orbs)

        max_GAS_size = maxval(GAS_sizes)
        n_spin_orbs = sum(GAS_sizes)

        allocate(GAS_table(n_spin_orbs))
        GAS_table(1::2) = spat_GAS_orbs
        GAS_table(2::2) = spat_GAS_orbs

        allocate(splitted_orbitals(max_GAS_size, nGAS))
        block
            integer :: all_orbs(n_spin_orbs), splitted_sizes(nGAS)
            all_orbs = [(i, i = 1, n_spin_orbs)]
            splitted_sizes = 0
            do iel = 1, size(all_orbs)
                iGAS = GAS_table(iel)
                splitted_sizes(iGAS) = splitted_sizes(iGAS) + 1
                splitted_orbitals(splitted_sizes(iGAS), iGAS) = all_orbs(iel)
            end do
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (all(GAS_sizes == splitted_sizes))) then
            call stop_all (this_routine, "Assert fail: all(GAS_sizes == splitted_sizes)")
        end if
    end block
#endif
        end block

        GAS_spec = CumulGASSpec_t(&
                c_min=cn_min, c_max=cn_max, GAS_table=GAS_table, &
                GAS_sizes=GAS_sizes, largest_GAS_size=max_GAS_size, &
                splitted_orbitals=splitted_orbitals, &
                lookup_is_connected=any(cn_min(:) /= cn_max(:)), &
                exchange_recoupling=recoupling_)

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

        contains

        pure function frequency(N) result(res)
            integer, intent(in) :: N(:)
            integer, allocatable :: res(:)
            integer :: i
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (minval(N) == 1)) then
            call stop_all (this_routine, "Assert fail: minval(N) == 1")
        end if
    end block
#endif
            allocate(res(maxval(N)), source=0)
            do i = 1, size(N)
                res(N(i)) = res(N(i)) + 1
            end do
        end function
    end function

    !> Query wether a supergroup is contained in the GAS space.
    pure function Cumul_contains_supergroup(this, supergroup) result(res)
        class(CumulGASSpec_t), intent(in) :: this
        integer, intent(in) :: supergroup(:)
        logical :: res
        associate(cumulated => cumsum(supergroup))
            res = all(this%c_min(:) <= cumulated &
                      .and. cumulated <= this%c_max(:) &
                      .and. supergroup <= this%GAS_sizes(:))
        end associate
    end function

    !> Check if the GAS specification is valid
    !>
    !> If the number of particles or the number of spin orbitals
    !> is provided, then the consistency with these numbers
    !> is checked as well.
    pure function Cumul_is_valid(this, n_basis) result(is_valid)
        !>  Specification of GAS spaces.
        class(CumulGASSpec_t), intent(in) :: this
        !>  The number of spin orbitals.
        integer, intent(in), optional :: n_basis
        logical :: is_valid

        logical :: shapes_match, pauli_principle, monotonic, &
            n_orbs_correct

        associate(GAS_sizes => this%GAS_size(), n_min => this%get_cmin(), &
                  n_max => this%get_cmax(), nGAS => this%nGAS())

            shapes_match = &
                all([size(GAS_sizes), size(n_min), size(n_max)] == nGAS) &
                .and. maxval(this%GAS_table) == nGAS

            pauli_principle = all(n_min(:) <= cumsum(GAS_sizes))

            if (nGAS >= 2) then
                monotonic = all([all(n_max(2:) >= n_max(: nGAS - 1)), &
                                 all(n_min(2:) >= n_min(: nGAS - 1))])
            else
                monotonic = .true.
            end if

            if (present(n_basis)) then
                n_orbs_correct = sum(GAS_sizes) == n_basis
            else
                n_orbs_correct = .true.
            end if
        end associate

        is_valid = all([shapes_match, pauli_principle, monotonic, n_orbs_correct])
    end function

    !> Write a string representation of this GAS specification to iunit
    subroutine Cumul_write_to(this, iunit)
        class(CumulGASSpec_t), intent(in) :: this
        integer, intent(in) :: iunit
        integer :: iGAS, iorb

        write(iunit, '(A)') 'Cumulative constraints'
        write(iunit, '(A)') 'n_i: number of spatial orbitals per i-th GAS space'
        write(iunit, '(A)') 'cn_min_i: minimum of cumulative particle number per i-th GAS space'
        write(iunit, '(A)') 'cn_max_i: maximum of cumulative particle number per i-th GAS space'
        write(iunit, '(A10, 1x, A10, 1x, A10)') 'n_i', 'cn_min_i', 'cn_max_i'
        write(iunit, '(A)') '--------------------------------'
        do iGAS = 1, this%nGAS()
            write(iunit, '(I10, 1x, I10, 1x, I10)') &
                this%GAS_size(iGAS) .div. 2, this%get_cmin(iGAS), this%get_cmax(iGAS)
        end do
        write(iunit, '(A)') '--------------------------------'
        write(iunit, '(A)') 'The distribution of spatial orbitals to GAS spaces is given by:'
        do iorb = 1, this%n_spin_orbs(), 2
            write(iunit, '(I0, 1x)', advance='no') this%get_iGAS(iorb)
        end do
        write(iunit, *)

        if (any(this%GAS_size() < (this%get_cmax() - eoshift(this%get_cmin(), -1)))) then
            write(iunit, '(A)') 'In at least one GAS space, the maximum allowed particle number by GAS constraints'
            write(iunit, '(A)') '   is larger than the particle number allowed by the Pauli principle.'
            write(iunit, '(A)') '   Was this intended when preparing your input?'
        end if

        if (.not. this%recoupling()) then
            write(iunit, '(A)') 'Double excitations with exchange are forbidden.'
        end if
    end subroutine


    !> Returns the minimum particle number for a given GAS space.
    integer elemental function get_cmin_i(this, iGAS)
        class(CumulGASSpec_t), intent(in) :: this
        integer, intent(in) :: iGAS
        get_cmin_i = this%c_min(iGAS)
    end function

    !> Returns the minimum particle number for all GAS spaces.
    pure function get_cmin_all(this) result(res)
        class(CumulGASSpec_t), intent(in) :: this
        integer, allocatable :: res(:)
        res = this%c_min(:)
    end function

    !> Returns the maximum particle number for a given GAS space.
    integer elemental function get_cmax_i(this, iGAS)
        class(CumulGASSpec_t), intent(in) :: this
        integer, intent(in) :: iGAS
        get_cmax_i = this%c_max(iGAS)
    end function

    !> Returns the maximum particle number for all GAS spaces.
    pure function get_cmax_all(this) result(res)
        class(CumulGASSpec_t), intent(in) :: this
        integer, allocatable :: res(:)
        res = this%c_max(:)
    end function

    pure function Cumul_get_possible_spaces(this, supergroup, add_holes, add_particles, n_total) result(spaces)
        class(CumulGASSpec_t), intent(in) :: this
        integer, intent(in) :: supergroup(size(this%GAS_sizes))
        integer, intent(in), optional :: add_holes(:), add_particles(:), n_total
        !> Lower and upper bound for spaces where a particle can be created.
        !> If no particle can be created, then spaces == 0 .
        integer, allocatable :: spaces(:)

        integer :: n_total_, iGAS, lower_bound, upper_bound
        type(buffer_int_1D_t) :: space_buffer

        integer :: &
        !> Cumulated number of particles per iGAS
            cum_n_particle(this%nGAS()), &
        !> Cumulated deficit per iGAS
            deficit(this%nGAS()), &
        !> Cumulated vacant orbitals per iGAS
            vacant(this%nGAS())

if(present(n_total)) then
    n_total_ = n_total
else
    n_total_ = 1
endif

        block
            integer :: B(this%nGAS()), C(this%nGAS())
            if (present(add_holes)) then
                B = this%count_per_GAS(add_holes)
            else
                B = 0
            end if
            if (present(add_particles)) then
                C = this%count_per_GAS(add_particles)
            else
                C = 0
            end if
            cum_n_particle = cumsum(supergroup - B + C)
        end block

        deficit(:) = this%get_cmin() - cum_n_particle(:)
        vacant(:) = this%get_cmax() - cum_n_particle(:)

        if (any(n_total_ < deficit) .or. all(vacant < n_total_)) then
            spaces = [integer :: ]
            return
        end if

        ! We assume that it is possible to create a particle at least in
        ! the last GAS space.
        ! Search from behind the first occurence where it is not possible
        ! anymore to create a particle.
        ! The lower bound is one GAS index above.
        lower_bound = custom_findloc(vacant <= 0, .true., back=.true.) + 1
        ! Find the first index, where a particle has to be created.
        upper_bound = custom_findloc(deficit, n_total_)

        if (lower_bound > upper_bound .or. lower_bound > this%nGAS()) then
            spaces = [integer :: ]
        else
        block
            integer :: new_deficit(size(deficit)), new_vacant(size(vacant)), jGAS
            logical :: allowed
            call space_buffer%init()
            do iGAS = lower_bound, upper_bound
                new_deficit = deficit
                new_deficit(iGAS :) = new_deficit(iGAS :) - 1
                new_vacant = vacant
                new_vacant(iGAS :) = new_vacant(iGAS :) - 1
                allowed = .true.
                do jGAS = 1, size(new_deficit) - 1
                    if (any(new_deficit(jGAS) > new_vacant(jGAS + 1 : ))) then
                        allowed = .false.
                        exit
                    end if
                end do
                if (allowed) call space_buffer%push_back(iGAS)
            end do
            call space_buffer%dump_reset(spaces)
        end block
        end if
    end function


    !>  Query wether a supergroup is contained in the GAS space.
    pure function Flexible_contains_supergroup(this, supergroup) result(res)
        !> Specification of GAS spaces.
        class(FlexibleGASSpec_t), intent(in) :: this
        !> A supergroup.
        integer, intent(in) :: supergroup(:)
        logical :: res
        integer :: i
        ! This function can be considerably sped up by applying
        ! the composition index trick from gasci_supergroup_index.fpp
        res = .false.
        do i = 1, size(this%supergroups, 2)
            if (all(supergroup == this%supergroups(:, i))) then
                res = .true.
                return
            end if
        end do
    end function


    !> Check if the GAS specification is valid
    !>
    !> If the number of particles or the number of spin orbitals
    !> is provided, then the consistency with these numbers
    !> is checked as well.
    logical pure function Flexible_is_valid(this, n_basis)
        !>  Specification of GAS spaces with local constraints.
        class(FlexibleGASSpec_t), intent(in) :: this
        integer, intent(in), optional :: n_basis

        logical :: shapes_match, pauli_principle, n_orbs_correct

        associate(GAS_sizes => this%GAS_sizes, nGAS => this%nGAS())

            shapes_match = &
                size(GAS_sizes) == nGAS .and. maxval(this%GAS_table) == nGAS

            pauli_principle = all(maxval(this%supergroups, dim=2) <= GAS_sizes)

            if (present(n_basis)) then
                n_orbs_correct = sum(GAS_sizes) == n_basis
            else
                n_orbs_correct = .true.
            end if
        end associate

        Flexible_is_valid = all([shapes_match, pauli_principle, n_orbs_correct])
    end function

    subroutine Flexible_write_to(this, iunit)
        class(FlexibleGASSpec_t), intent(in) :: this
        integer, intent(in) :: iunit
        integer :: i_sg

        write(iunit, '(A)') 'Flexible GAS constraints'
        if (.not. this%recoupling()) then
            write(iunit, '(A)') 'Double excitations with exchange are forbidden.'
        end if
        write(iunit, '(A)') 'The following supergroups are allowed'
        do i_sg = 1, size(this%supergroups, 2)
            write(iunit, *) this%supergroups(:, i_sg)
        end do
    end subroutine

    pure function Flexible_get_possible_spaces(this, supergroup, add_holes, add_particles, n_total) result(spaces)
        class(FlexibleGASSpec_t), intent(in) :: this
        integer, intent(in) :: supergroup(size(this%GAS_sizes))
        integer, intent(in), optional :: add_holes(:), add_particles(:), n_total
        integer, allocatable :: spaces(:)
        character(*), parameter :: this_routine = 'Flexible_get_possible_spaces'

        integer :: n_total_

if(present(n_total)) then
    n_total_ = n_total
else
    n_total_ = 1
endif
#ifdef WARNING_WORKAROUND_
        associate(supergroup => supergroup); end associate
        associate(add_holes => add_holes); end associate
        associate(add_particles => add_particles); end associate
#endif
        spaces = [-1]

        call stop_all(this_routine, "Has to be implemented.")
    end function

    !>  Constructor of FlexibleGASSpec_t
    pure function construct_FlexibleGASSpec_t(supergroups, spat_GAS_orbs, recoupling) result(GAS_spec)
        !> The allowed supergroups.
        integer, intent(in) :: supergroups(:, :)
        !> GAS space for the i-th **spatial** orbital.
        integer, intent(in) :: spat_GAS_orbs(:)
        !> Exchange double excitations that recouple the spin are allowed
        logical, intent(in), optional :: recoupling
        logical :: recoupling_

        type(FlexibleGASSpec_t) :: GAS_spec

        integer :: n_spin_orbs, max_GAS_size
        integer, allocatable :: splitted_orbitals(:, :), GAS_table(:), &
            GAS_sizes(:), sorted_supergroup(:, :)
        integer :: i, iel, iGAS, nGAS, N
        character(*), parameter :: this_routine = 'construct_FlexibleGASSpec_t'

if(present(recoupling)) then
    recoupling_ = recoupling
else
    recoupling_ = .true.
endif

        N = sum(supergroups(:, 1))
        if (any(N /= sum(supergroups, dim=1))) then
            call stop_all(&
                this_routine, 'Different particle numbers in the supergroups.')
        end if

        sorted_supergroup = supergroups
    

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

            associate(X => sorted_supergroup)
                ! 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_geq(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_geq(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


        nGAS = maxval(spat_GAS_orbs)
        GAS_sizes = 2 * frequency(spat_GAS_orbs)

        if (any(GAS_sizes < maxval(sorted_supergroup, dim=2))) then
            call stop_all(this_routine, 'Pauli violation: Too many particles per GAS space.')
        end if

        max_GAS_size = maxval(GAS_sizes)
        n_spin_orbs = sum(GAS_sizes)

        allocate(GAS_table(n_spin_orbs))
        GAS_table(1::2) = spat_GAS_orbs
        GAS_table(2::2) = spat_GAS_orbs

        allocate(splitted_orbitals(max_GAS_size, nGAS))
        block
            integer :: all_orbs(n_spin_orbs), splitted_sizes(nGAS)
            all_orbs = [(i, i = 1, n_spin_orbs)]
            splitted_sizes = 0
            do iel = 1, size(all_orbs)
                iGAS = GAS_table(iel)
                splitted_sizes(iGAS) = splitted_sizes(iGAS) + 1
                splitted_orbitals(splitted_sizes(iGAS), iGAS) = all_orbs(iel)
            end do
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (all(GAS_sizes == splitted_sizes))) then
            call stop_all (this_routine, "Assert fail: all(GAS_sizes == splitted_sizes)")
        end if
    end block
#endif
        end block


        GAS_spec = FlexibleGASSpec_t(&
                supergroups=sorted_supergroup, &
                N=N, &
                GAS_table=GAS_table, &
                GAS_sizes=GAS_sizes, largest_GAS_size=max_GAS_size, &
                splitted_orbitals=splitted_orbitals, &
                lookup_is_connected=size(sorted_supergroup, 2) > 1, &
                exchange_recoupling=recoupling_)

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

        contains

            pure function frequency(N) result(res)
                integer, intent(in) :: N(:)
                integer, allocatable :: res(:)
                integer :: i
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (minval(N) == 1)) then
            call stop_all (this_routine, "Assert fail: minval(N) == 1")
        end if
    end block
#endif
                allocate(res(maxval(N)), source=0)
                do i = 1, size(N)
                    res(N(i)) = res(N(i)) + 1
                end do
            end function
    end function

    elemental function Flexible_N_particle(this) result(res)
        class(FlexibleGASSpec_t), intent(in) :: this
        integer :: res
        res = this%N
    end function

end module gasci