#include "macros.h" module gasci use constants, only: n_int, dp, int64 use SystemData, only: nBasis, tGUGA use util_mod, only: lex_leq, cumsum, operator(.div.), near_zero, binary_search_first_ge, & operator(.isclose.), custom_findloc, EnumBase_t, lex_geq, stop_all 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_3_t, excite, canonicalize 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 sort_mod, only: sort use DetBitOps, only: ilut_lt, ilut_gt, EncodeBitDet use growing_buffers, only: buffer_int_2D_t, buffer_int_1D_t use guga_bitRepOps, only: CSF_Info_t use guga_data, only: DistinctDouble_t, ijab_Index_t use composition_utils, only: composition_idx, composition_from_idx better_implicit_none private public :: possible_GAS_exc_gen, & GAS_exc_gen, GAS_specification, GASSpec_t, & get_name, LocalGASSpec_t, CumulGASSpec_t, FlexibleGASSpec_t, tGAS, tSD_GAS, tGUGA_GAS, & CAS_spec, finalize_GAS, P_specification, frequency, & GAS_SpecKind_t, GAS_SpecKind_vals_t, GAS_spec_kind_vals logical :: tGAS = .false. !! Are we restricting excitations 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), allocatable :: GAS_exc_gen type, extends(EnumBase_t) :: GAS_SpecKind_t end type type :: GAS_SpecKind_vals_t type(GAS_SpecKind_t) :: & LOCAL = GAS_SpecKind_t(1), & CUMULATIVE = GAS_SpecKind_t(2), & FLEXIBLE = GAS_SpecKind_t(3) end type type(GAS_SpecKind_vals_t), parameter :: GAS_spec_kind_vals = GAS_SpecKind_vals_t() ! 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 !! Do we do exchange recoupling? integer :: N_particle !! The particle number 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, private :: contains_conf_nI procedure, private :: contains_conf_csf generic :: contains_conf => contains_conf_nI generic :: contains_conf => contains_conf_csf procedure :: contains_ilut procedure :: is_connected => get_is_connected procedure :: nGAS => get_nGAS procedure :: max_GAS_size => get_max_GAS_size procedure :: recoupling procedure :: n_spin_orbs => get_nOrbs procedure :: nEl => get_nEl 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, private :: count_per_GAS_nI procedure, private :: count_per_GAS_csf generic :: count_per_GAS => count_per_GAS_nI, count_per_GAS_csf procedure, private :: is_allowed_single procedure, private :: is_allowed_double procedure, private :: is_allowed_triple procedure, private :: is_allowed_distinct_double generic :: is_allowed => is_allowed_single, is_allowed_double, is_allowed_triple, is_allowed_distinct_double procedure, private :: excite_single procedure, private :: excite_double procedure, private :: excite_triple generic :: excite => excite_single, excite_double, excite_triple procedure :: get_available_singles procedure :: get_available_doubles procedure :: gen_all_excits 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(:, :) integer(int64), public, allocatable :: composition_indices(:) 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 end type interface FlexibleGASSpec_t module procedure construct_FlexibleGASSpec_t end interface class(GASSpec_t), allocatable :: GAS_specification class(GASSpec_t), allocatable :: P_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 is contained in the GAS space. pure function contains_conf_nI(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 debug_function_name("contains_conf_nI") #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (size(nI) == this%nEl())) then call stop_all (this_routine, "Assert fail: size(nI) == this%nEl()") end if end block #endif res = this%contains_supergroup(this%count_per_GAS(nI)) end function !> Query wether a CSF is contained in the GAS space. pure function contains_conf_csf(this, csf_i) result(res) !> Specification of GAS spaces. class(GASSpec_t), intent(in) :: this type(CSF_Info_t), intent(in) :: csf_i logical :: res debug_function_name("contains_conf_csf") #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (sum(csf_i%occ_int) == this%nEl())) then call stop_all (this_routine, "Assert fail: sum(csf_i%occ_int) == this%nEl()") end if end block #endif res = this%contains_supergroup(this%count_per_GAS(csf_i)) end function logical elemental function recoupling(this) !! Query whether we do exchange recoupling. class(GASSpec_t), intent(in) :: this recoupling = this%exchange_recoupling end function integer elemental function get_nEl(this) !! Get the number of electrons. class(GASSpec_t), intent(in) :: this get_nEl = this%N_particle 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_nI(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 !> Count the particles per GAS space. i.e. return the supergroup. pure function count_per_GAS_csf(this, csf_i) result(supergroup) class(GASSpec_t), intent(in) :: this type(CSF_Info_t), intent(in) :: csf_i integer :: supergroup(this%nGAS()) integer :: i_orb, iGAS supergroup = 0 do i_orb = 1, size(csf_i%occ_int) iGAS = this%get_iGAS(i_orb * 2) supergroup(iGAS) = supergroup(iGAS) + csf_i%occ_int(i_orb) 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 !Query the supergroup after a single excitation pure function excite_single(this, exc, supergroup) result(excited_supergroup) class(GASSpec_t), intent(in) :: this type(Excite_1_t), intent(in) :: exc integer, intent(in) :: supergroup(:) integer :: src_space, tgt_space integer :: excited_supergroup(size(supergroup)) !get source and target space src_space = this%get_iGAS(exc%val(1, 1)) tgt_space = this%get_iGAS(exc%val(2, 1)) excited_supergroup = supergroup excited_supergroup(src_space) = excited_supergroup(src_space) - 1 excited_supergroup(tgt_space) = excited_supergroup(tgt_space) + 1 end function !Query the supergroup after a double excitation pure function excite_double(this, exc, supergroup) result(excited_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), i src_spaces(:) = this%get_iGAS(exc%val(1, :)) tgt_spaces(:) = this%get_iGAS(exc%val(2, :)) excited_supergroup = supergroup do i = 1, 2 !Remove particles form src_spaces excited_supergroup(src_spaces(i)) = excited_supergroup(src_spaces(i)) - 1 !Add particles to tgt_spaces excited_supergroup(tgt_spaces(i)) = excited_supergroup(tgt_spaces(i)) + 1 end do end function !Query the supergroup after a double excitation pure function excite_triple(this, exc, supergroup) result(excited_supergroup) class(GASSpec_t), intent(in) :: this type(Excite_3_t), intent(in) :: exc integer, intent(in) :: supergroup(:) integer :: excited_supergroup(size(supergroup)) integer :: src_spaces(3), tgt_spaces(3), i src_spaces(:) = this%get_iGAS(exc%val(1, :)) tgt_spaces(:) = this%get_iGAS(exc%val(2, :)) excited_supergroup = supergroup do i = 1, 3 !Remove particles form src_spaces excited_supergroup(src_spaces(i)) = excited_supergroup(src_spaces(i)) - 1 !Add particles to tgt_spaces excited_supergroup(tgt_spaces(i)) = excited_supergroup(tgt_spaces(i)) + 1 end do 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 ! Keep this for following if-statement 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 = this%excite_single(exc, supergroup) 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) ! Keep for if-statement src_spaces(1) = this%get_iGAS(exc%val(1, 1)) src_spaces(2) = this%get_iGAS(exc%val(1, 2)) 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 = this%excite_double(exc, supergroup) 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 !> Check if a triple 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_triple(this, exc, supergroup) class(GASSpec_t), intent(in) :: this type(Excite_3_t), intent(in) :: exc integer, intent(in) :: supergroup(:) integer :: excited_supergroup(size(supergroup)) integer :: src_spaces(3), tgt_spaces(3) ! Keep for if-statement 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) .and. src_spaces(2) == src_spaces(3)) then ! All electrons come from the same space and there are no restrictions ! regarding recoupling or GAS. is_allowed_triple = .true. else ! Ensure that GAS specifications contain supergroup **after** excitation. excited_supergroup = this%excite(exc, supergroup) is_allowed_triple = this%contains_supergroup(excited_supergroup) end if end function logical pure function is_allowed_distinct_double(this, exc, supergroup) class(GASSpec_t), intent(in) :: this type(DistinctDouble_t), intent(in) :: exc integer, intent(in) :: supergroup(:) integer :: i, j, a, b call exc%idx%extract_ijkl(a, i, b, j) is_allowed_distinct_double = this%is_allowed(Excite_2_t(2 * i, 2 * a, 2 * j, 2 * b), supergroup) 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 pure function construct_LocalGASSpec_t(n_min, n_max, spat_GAS_orbs, nEl, recoupling) result(GAS_spec) !! Constructor of LocalGASSpec_t integer, intent(in) :: n_min(:) !! Minimum particle number per GAS space. integer, intent(in) :: n_max(:) !! Maximum particle number per GAS space. integer, intent(in) :: spat_GAS_orbs(:) !! GAS space for the i-th **spatial** orbital. integer, intent(in) :: nEl !! Number of particles. logical, intent(in), optional :: recoupling !! Exchange double excitations that recouple the spin are allowed 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_, & n_particle=nEl) 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 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, & has_enough_room 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) has_enough_room = sum(this%get_min()) <= this%nEl() .and. this%nEl() <= sum(this%get_max()) 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, & has_enough_room]) 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 pure function construct_CumulGASSpec_t(cn_min, cn_max, spat_GAS_orbs, nEl, recoupling) result(GAS_spec) !! Constructor of CumulGASSpec_t integer, intent(in) :: cn_min(:) !! Cumulative minimum particle number. integer, intent(in) :: cn_max(:) !! Cumulative maximum particle number. integer, intent(in) :: spat_GAS_orbs(:) !! GAS space for the i-th **spatial** orbital. integer, intent(in) :: nEl !! Number of particles. logical, intent(in), optional :: recoupling !! Exchange double excitations that recouple the spin are allowed 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_, & n_particle=nEl) #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 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, has_enough_room 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 has_enough_room = this%get_cmin(this%nGAS()) <= this%nEl() & .and. this%nEl() <= this%get_cmax(this%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, & has_enough_room, 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 debug_function_name("Flexible_is_valid") logical :: shapes_match, pauli_principle, n_orbs_correct, & has_enough_room 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) has_enough_room = this%nEl() == sum(this%supergroups(:, 1)) #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (all(this%nEl() == sum(this%supergroups, dim=1)))) then call stop_all (this_routine, "Assert fail: all(this%nEl() == sum(this%supergroups, dim=1))") end if end block #endif 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, & has_enough_room]) 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(int64), allocatable :: composition_indices(:) 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 allocate(composition_indices(size(sorted_supergroup, dim=2))) block integer :: i_sg do i_sg = 1, size(composition_indices) composition_indices(i_sg) = composition_idx(sorted_supergroup(:, i_sg)) end do 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_particle=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_, & composition_indices=composition_indices) #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 end function pure logical function tSD_GAS() !! Doing GAS in Slater Determinant (SD) space tSD_GAS = tGAS .and. .not. tGUGA end function pure logical function tGUGA_GAS() !! Doing GAS in GUGA space tGUGA_GAS = tGAS .and. tGUGA end function type(LocalGASSpec_t) pure function CAS_spec(n_el, n_spat_orbs) integer, intent(in) :: n_el, n_spat_orbs integer :: i ! It does not matter if we use local or cumulative GAS ! constraints CAS_spec = LocalGASSpec_t(n_min=[n_el], n_max=[n_el], & spat_GAS_orbs=[(1, i = 1, n_spat_orbs)], nEl=n_el) end function subroutine finalize_GAS() if (allocated(GAS_specification)) deallocate(GAS_specification) if (allocated(GAS_exc_gen)) deallocate(GAS_exc_gen) end subroutine pure function get_available_singles(this, det_I) result(singles_exc_list) !> Get all single excitated determinants from det_I that are allowed under GAS constraints. class(GASSpec_t), intent(in) :: this 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. (this%contains_conf(det_I))) then call stop_all (this_routine, "Assert fail: this%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 = this%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(this, det_I) result(doubles_exc_list) class(GASSpec_t), intent(in) :: this 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. (this%contains_conf(det_I))) then call stop_all (this_routine, "Assert fail: this%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 = this%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 = this%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 subroutine gen_all_excits(this, nI, n_excits, det_list, ic) !! Get all excitated determinants from !! det_I that are allowed under GAS constraints. class(GASSpec_t), intent(in) :: this !! GAS specification integer, intent(in) :: nI(:) !! Starting determinant integer, intent(out) :: n_excits !! Number of determinants integer(n_int), intent(out), allocatable :: det_list(:,:) !! Allocatable array of determinants in ilut format integer, intent(in), optional :: ic !! Optional input for excitation level (ic=1 => singles, ic=2 => doubles) !! If ommited generate all. integer, allocatable :: singles(:, :), doubles(:, :) integer :: i, j if (present(ic)) then select case(ic) case(1) singles = this%get_available_singles(nI) allocate(doubles(0, 0)) case(2) allocate(singles(0, 0)) doubles = this%get_available_doubles(nI) end select else singles = this%get_available_singles(nI) doubles = this%get_available_doubles(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 pure function frequency(N) result(res) !!! Count the frequency of numbers in the array. integer, intent(in) :: N(:) integer, allocatable :: res(:) integer :: i debug_function_name("gasci::frequency") #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 module gasci