#include "macros.h" ! all excitations with rank higher than max_excit_rank are definitely zero ! note that this excludes further excitations, which must be handled manually ! Excite_Further_t is for all ranks > max_excit_rank !> A module for representing different excitations. !> !> There is one abstract base class excitation_t that represents !> an arbitrary excitation. !> Possible excitations are the !> non trivial excitations Excite_{1,2,3}_t !> and the trivial excitations Excite_0_t and Excite_Further_t. !> !> The non trivial excitations can "know" only some indices !> and leave the rest UNKNOWN, which is an arbitrary integer constant. !> !> The procedures create_excitation, get_excitation, and get_bit_excitation !> can be used, to create excitations from nIs, or iluts at runtime. module excitation_types use orb_idx_mod, only: size, sum, calc_spin_raw use constants, only: dp, n_int, bits_n_int, maxExcit use bit_rep_data, only: nIfTot use util_mod, only: stop_all use SystemData, only: nEl, G1 use orb_idx_mod, only: SpinOrbIdx_t use sets_mod, only: disjoint, subset, is_sorted, special_union_complement, is_set, & operator(.in.) use DetBitOps, only: GetBitExcitation use excit_mod, only: GetExcitation use sort_mod, only: sort implicit none private public :: Excitation_t, UNKNOWN, defined, dyn_defined, get_last_tgt, set_last_tgt, & create_excitation, get_excitation, get_bit_excitation, & ilut_excite, excite, dyn_excite, dyn_nI_excite, & is_sorted, is_canonical, spin_allowed, canonicalize, occupation_allowed, make_canonical public :: Excite_0_t public :: Excite_1_t public :: Excite_2_t public :: Excite_3_t public :: Excite_Further_t !> Arbitrary non occuring (?!) orbital index. integer, parameter :: UNKNOWN = huge(UNKNOWN) !> Abstract base class for excitations. type, abstract :: Excitation_t end type type, extends(Excitation_t) :: Excite_0_t !! Represents the orbital indices of a 0-order excitation !! The array is sorted like: !! [srcs, tgts] integer :: val(2, 0) = UNKNOWN end type type, extends(Excitation_t) :: Excite_1_t !! Represents the orbital indices of a 1-order excitation !! The array is sorted like: !! [srcs, tgts] integer :: val(2, 1) = UNKNOWN end type type, extends(Excitation_t) :: Excite_2_t !! Represents the orbital indices of a 2-order excitation !! The array is sorted like: !! [srcs, tgts] integer :: val(2, 2) = UNKNOWN end type type, extends(Excitation_t) :: Excite_3_t !! Represents the orbital indices of a 3-order excitation !! The array is sorted like: !! [srcs, tgts] integer :: val(2, 3) = UNKNOWN end type type, extends(Excitation_t) :: Excite_Further_t !! Represents an excitation with so many different indices, it has to be zero integer :: val(2, 0) = UNKNOWN end type !> Additional constructors for the excitation types from integers instead !> of an integer array. !> !> The non trivial excitations !> are initialized by passing the respective integer arrays into !> the type. !> Alternatively one can use integer arguments to initialize. !> Omitted indices are set to UNKNOWN. !> !> \code{.unparsed} !> Excite_t([1, 2]) == Excite_t(src=1, tgt=2) !> ! If the target should be UNKNOWN, just omit it !> Excite_t(src=1) !> \endcode !> !> The signature is `(src_1, tgt_1, src_2, tgt_2, ...)`. !> depending on the actual type. interface Excite_1_t module procedure from_integer_Excite_1_t end interface !> Additional constructors for the excitation types from integers instead !> of an integer array. !> !> The non trivial excitations !> are initialized by passing the respective integer arrays into !> the type. !> Alternatively one can use integer arguments to initialize. !> Omitted indices are set to UNKNOWN. !> !> \code{.unparsed} !> Excite_t([1, 2]) == Excite_t(src=1, tgt=2) !> ! If the target should be UNKNOWN, just omit it !> Excite_t(src=1) !> \endcode !> !> The signature is `(src_1, tgt_1, src_2, tgt_2, ...)`. !> depending on the actual type. interface Excite_2_t module procedure from_integer_Excite_2_t end interface !> Additional constructors for the excitation types from integers instead !> of an integer array. !> !> The non trivial excitations !> are initialized by passing the respective integer arrays into !> the type. !> Alternatively one can use integer arguments to initialize. !> Omitted indices are set to UNKNOWN. !> !> \code{.unparsed} !> Excite_t([1, 2]) == Excite_t(src=1, tgt=2) !> ! If the target should be UNKNOWN, just omit it !> Excite_t(src=1) !> \endcode !> !> The signature is `(src_1, tgt_1, src_2, tgt_2, ...)`. !> depending on the actual type. interface Excite_3_t module procedure from_integer_Excite_3_t end interface ! workaround for pesky intel compiler errors ! this should be removed if IFORT_ works properly at some pointer in the future #ifdef IFORT_ interface Excite_0_t module procedure Excite_0_t_ctor end interface #endif !> Return true if all sources and targets are not UNKNOWN. interface defined module procedure defined_Excite_0_t module procedure defined_Excite_1_t module procedure defined_Excite_2_t module procedure defined_Excite_3_t end interface !> Return true if all sources and targets are not UNKNOWN. interface is_sorted module procedure is_sorted_Excite_0_t module procedure is_sorted_Excite_1_t module procedure is_sorted_Excite_2_t module procedure is_sorted_Excite_3_t end interface !> Return true if the excitation is canonical !> !> Canonical means: !> 1. that the excitation is defined, i.e. it has no UNKNOWN, !> 2. the sources and the targets are sets, i.e. they are unique and ordered, !> 3. the sources and the targets are disjoint, interface is_canonical module procedure is_canonical_Excite_0_t module procedure is_canonical_Excite_1_t module procedure is_canonical_Excite_2_t module procedure is_canonical_Excite_3_t end interface !> Return true if the excitation preserves the overall spin-projection interface spin_allowed module procedure spin_allowed_Excite_0_t module procedure spin_allowed_Excite_1_t module procedure spin_allowed_Excite_2_t module procedure spin_allowed_Excite_3_t end interface !> Return true if the excitation is allowed by occupation of the starting determinant !> !> The input excitation has to be canonical. interface occupation_allowed module procedure occupation_allowed_Excite_0_t module procedure occupation_allowed_Excite_1_t module procedure occupation_allowed_Excite_2_t module procedure occupation_allowed_Excite_3_t end interface !> Canonicalize an excitation !> !> Canonical means that the excitation is defined, i.e. it has no UNKNOWN, !> the sources and the targets are sets, i.e. they are unique and ordered, !> and the sources and the targets are disjoint. interface canonicalize module procedure canonicalize_Excite_0_t module procedure canonicalize_Excite_1_t module procedure canonicalize_Excite_2_t module procedure canonicalize_Excite_3_t end interface !> Canonicalize an excitation and count the necessary swaps !> !> Canonical means that the excitation is defined, i.e. it has no UNKNOWN, !> the sources and the targets are sets, i.e. they are unique and ordered, !> and the sources and the targets are disjoint. interface make_canonical module procedure make_canonical_Excite_0_t module procedure make_canonical_Excite_1_t module procedure make_canonical_Excite_2_t module procedure make_canonical_Excite_3_t end interface !> Get the last target of a non trivial excitation. interface get_last_tgt module procedure get_last_tgt_Excite_1_t module procedure get_last_tgt_Excite_2_t module procedure get_last_tgt_Excite_3_t end interface !> Set the last target of a non trivial excitation. interface set_last_tgt module procedure set_last_tgt_Excite_1_t module procedure set_last_tgt_Excite_2_t module procedure set_last_tgt_Excite_3_t end interface !> Perform the excitation on a given determinant. !> !> It is assumed that the excitations are non trivial. !> I.e. for single excitations source /= target !> and for double excitations the set of sources and targets has !> to be disjoint. interface excite module procedure excite_nI_Excite_0_t module procedure excite_nI_Excite_1_t module procedure excite_nI_Excite_2_t module procedure excite_nI_Excite_3_t module procedure excite_SpinOrbIdx_t_Excite_0_t module procedure excite_SpinOrbIdx_t_Excite_1_t module procedure excite_SpinOrbIdx_t_Excite_2_t module procedure excite_SpinOrbIdx_t_Excite_3_t end interface interface ilut_excite module procedure excite_Ilut_t_Excite_0_t module procedure excite_Ilut_t_Excite_1_t module procedure excite_Ilut_t_Excite_2_t module procedure excite_Ilut_t_Excite_3_t end interface interface get_excitation module procedure get_excitation_old, get_excitation_new end interface contains ! workaround for pesky intel compiler errors ! this should be removed if IFORT_ works properly at some pointer in the future #ifdef IFORT_ type(Excite_0_t) function Excite_0_t_ctor() result(this) integer :: tmpval(2, 0) tmpval = UNKNOWN this%val = tmpval end function #endif elemental function defined_Excite_0_t (exc) result(res) type(Excite_0_t), intent(in) :: exc logical :: res res = all(exc%val /= UNKNOWN) end function elemental function is_sorted_Excite_0_t (exc) result(res) type(Excite_0_t), intent(in) :: exc logical :: res res = is_sorted(exc%val(1, :)) .and. is_sorted(exc%val(2, :)) end function elemental function is_canonical_Excite_0_t (exc) result(res) type(Excite_0_t), intent(in) :: exc logical :: res res = .false. if (.not. defined(exc)) return if (.not. (is_set(exc%val(1, :)) .and. is_set(exc%val(2, :)))) return if (.not. disjoint(exc%val(1, :), exc%val(2, :))) return res = .true. end function elemental function spin_allowed_Excite_0_t (exc) result(res) type(Excite_0_t), intent(in) :: exc logical :: res res = sum(G1(exc%val(1, :))%Ms) == sum(G1(exc%val(2, :))%Ms) end function pure function occupation_allowed_Excite_0_t (nI, exc) result(res) integer, intent(in) :: nI(:) type(Excite_0_t), intent(in) :: exc logical :: res debug_function_name("occupation_allowed_Excite_0_t") #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif res = subset(exc%val(1, :), nI) .and. disjoint(exc%val(2, :), nI) end function elemental subroutine make_canonical_Excite_0_t (exc, even_swaps) type(Excite_0_t), intent(inout) :: exc logical, intent(out) :: even_swaps integer :: n_swaps(2) type(Excite_0_t) :: copy_exc routine_name("make_canonical") copy_exc = exc n_swaps(1) = 0 ! merge_sort block integer, dimension(:), allocatable :: tmp integer :: current_size, left, mid, right integer, parameter :: along = 1, & bubblesort_size = 20 associate(X => exc%val(1, :)) ! 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(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, 1), lbound(V, 1) + 1, -1 do i = lbound(V, 1), ubound(V, 1) - 1 if (.not. V(i) <= V(i + 1)) then ! swap block n_swaps(1) = n_swaps(1) + 1 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, 1) j = lbound(B, 1) do k = lbound(C, 1), ubound(C, 1) if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then if ( A(i) <= B(j)) then C(k) = A(i) i = i + 1 else C(k) = B(j) j = j + 1 n_swaps(1) = n_swaps(1) + 1 end if else if (i <= ubound(A, 1)) then C(k) = A(i) i = i + 1 else if (j <= ubound(B, 1)) then C(k) = B(j) j = j + 1 n_swaps(1) = n_swaps(1) + 1 end if end do end associate end block end do current_size = 2 * current_size end do end associate end block n_swaps(2) = 0 ! merge_sort block integer, dimension(:), allocatable :: tmp integer :: current_size, left, mid, right integer, parameter :: along = 1, & bubblesort_size = 20 associate(X => exc%val(2, :)) ! 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(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, 1), lbound(V, 1) + 1, -1 do i = lbound(V, 1), ubound(V, 1) - 1 if (.not. V(i) <= V(i + 1)) then ! swap block n_swaps(2) = n_swaps(2) + 1 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, 1) j = lbound(B, 1) do k = lbound(C, 1), ubound(C, 1) if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then if ( A(i) <= B(j)) then C(k) = A(i) i = i + 1 else C(k) = B(j) j = j + 1 n_swaps(2) = n_swaps(2) + 1 end if else if (i <= ubound(A, 1)) then C(k) = A(i) i = i + 1 else if (j <= ubound(B, 1)) then C(k) = B(j) j = j + 1 n_swaps(2) = n_swaps(2) + 1 end if end do end associate end block end do current_size = 2 * current_size end do end associate end block even_swaps = mod(sum(n_swaps), 2) == 0 #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif end subroutine elemental function canonicalize_Excite_0_t (exc) result(res) type(Excite_0_t), intent(in) :: exc type(Excite_0_t) :: res logical :: dummy res = exc call make_canonical(res, dummy) end function elemental function defined_Excite_1_t (exc) result(res) type(Excite_1_t), intent(in) :: exc logical :: res res = all(exc%val /= UNKNOWN) end function elemental function is_sorted_Excite_1_t (exc) result(res) type(Excite_1_t), intent(in) :: exc logical :: res res = is_sorted(exc%val(1, :)) .and. is_sorted(exc%val(2, :)) end function elemental function is_canonical_Excite_1_t (exc) result(res) type(Excite_1_t), intent(in) :: exc logical :: res res = .false. if (.not. defined(exc)) return if (.not. (is_set(exc%val(1, :)) .and. is_set(exc%val(2, :)))) return if (.not. disjoint(exc%val(1, :), exc%val(2, :))) return res = .true. end function elemental function spin_allowed_Excite_1_t (exc) result(res) type(Excite_1_t), intent(in) :: exc logical :: res res = sum(G1(exc%val(1, :))%Ms) == sum(G1(exc%val(2, :))%Ms) end function pure function occupation_allowed_Excite_1_t (nI, exc) result(res) integer, intent(in) :: nI(:) type(Excite_1_t), intent(in) :: exc logical :: res debug_function_name("occupation_allowed_Excite_1_t") #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif res = subset(exc%val(1, :), nI) .and. disjoint(exc%val(2, :), nI) end function elemental subroutine make_canonical_Excite_1_t (exc, even_swaps) type(Excite_1_t), intent(inout) :: exc logical, intent(out) :: even_swaps integer :: n_swaps(2) type(Excite_1_t) :: copy_exc routine_name("make_canonical") copy_exc = exc n_swaps(1) = 0 ! merge_sort block integer, dimension(:), allocatable :: tmp integer :: current_size, left, mid, right integer, parameter :: along = 1, & bubblesort_size = 20 associate(X => exc%val(1, :)) ! 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(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, 1), lbound(V, 1) + 1, -1 do i = lbound(V, 1), ubound(V, 1) - 1 if (.not. V(i) <= V(i + 1)) then ! swap block n_swaps(1) = n_swaps(1) + 1 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, 1) j = lbound(B, 1) do k = lbound(C, 1), ubound(C, 1) if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then if ( A(i) <= B(j)) then C(k) = A(i) i = i + 1 else C(k) = B(j) j = j + 1 n_swaps(1) = n_swaps(1) + 1 end if else if (i <= ubound(A, 1)) then C(k) = A(i) i = i + 1 else if (j <= ubound(B, 1)) then C(k) = B(j) j = j + 1 n_swaps(1) = n_swaps(1) + 1 end if end do end associate end block end do current_size = 2 * current_size end do end associate end block n_swaps(2) = 0 ! merge_sort block integer, dimension(:), allocatable :: tmp integer :: current_size, left, mid, right integer, parameter :: along = 1, & bubblesort_size = 20 associate(X => exc%val(2, :)) ! 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(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, 1), lbound(V, 1) + 1, -1 do i = lbound(V, 1), ubound(V, 1) - 1 if (.not. V(i) <= V(i + 1)) then ! swap block n_swaps(2) = n_swaps(2) + 1 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, 1) j = lbound(B, 1) do k = lbound(C, 1), ubound(C, 1) if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then if ( A(i) <= B(j)) then C(k) = A(i) i = i + 1 else C(k) = B(j) j = j + 1 n_swaps(2) = n_swaps(2) + 1 end if else if (i <= ubound(A, 1)) then C(k) = A(i) i = i + 1 else if (j <= ubound(B, 1)) then C(k) = B(j) j = j + 1 n_swaps(2) = n_swaps(2) + 1 end if end do end associate end block end do current_size = 2 * current_size end do end associate end block even_swaps = mod(sum(n_swaps), 2) == 0 #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif end subroutine elemental function canonicalize_Excite_1_t (exc) result(res) type(Excite_1_t), intent(in) :: exc type(Excite_1_t) :: res logical :: dummy res = exc call make_canonical(res, dummy) end function elemental function defined_Excite_2_t (exc) result(res) type(Excite_2_t), intent(in) :: exc logical :: res res = all(exc%val /= UNKNOWN) end function elemental function is_sorted_Excite_2_t (exc) result(res) type(Excite_2_t), intent(in) :: exc logical :: res res = is_sorted(exc%val(1, :)) .and. is_sorted(exc%val(2, :)) end function elemental function is_canonical_Excite_2_t (exc) result(res) type(Excite_2_t), intent(in) :: exc logical :: res res = .false. if (.not. defined(exc)) return if (.not. (is_set(exc%val(1, :)) .and. is_set(exc%val(2, :)))) return if (.not. disjoint(exc%val(1, :), exc%val(2, :))) return res = .true. end function elemental function spin_allowed_Excite_2_t (exc) result(res) type(Excite_2_t), intent(in) :: exc logical :: res res = sum(G1(exc%val(1, :))%Ms) == sum(G1(exc%val(2, :))%Ms) end function pure function occupation_allowed_Excite_2_t (nI, exc) result(res) integer, intent(in) :: nI(:) type(Excite_2_t), intent(in) :: exc logical :: res debug_function_name("occupation_allowed_Excite_2_t") #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif res = subset(exc%val(1, :), nI) .and. disjoint(exc%val(2, :), nI) end function elemental subroutine make_canonical_Excite_2_t (exc, even_swaps) type(Excite_2_t), intent(inout) :: exc logical, intent(out) :: even_swaps integer :: n_swaps(2) type(Excite_2_t) :: copy_exc routine_name("make_canonical") copy_exc = exc n_swaps(1) = 0 ! merge_sort block integer, dimension(:), allocatable :: tmp integer :: current_size, left, mid, right integer, parameter :: along = 1, & bubblesort_size = 20 associate(X => exc%val(1, :)) ! 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(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, 1), lbound(V, 1) + 1, -1 do i = lbound(V, 1), ubound(V, 1) - 1 if (.not. V(i) <= V(i + 1)) then ! swap block n_swaps(1) = n_swaps(1) + 1 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, 1) j = lbound(B, 1) do k = lbound(C, 1), ubound(C, 1) if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then if ( A(i) <= B(j)) then C(k) = A(i) i = i + 1 else C(k) = B(j) j = j + 1 n_swaps(1) = n_swaps(1) + 1 end if else if (i <= ubound(A, 1)) then C(k) = A(i) i = i + 1 else if (j <= ubound(B, 1)) then C(k) = B(j) j = j + 1 n_swaps(1) = n_swaps(1) + 1 end if end do end associate end block end do current_size = 2 * current_size end do end associate end block n_swaps(2) = 0 ! merge_sort block integer, dimension(:), allocatable :: tmp integer :: current_size, left, mid, right integer, parameter :: along = 1, & bubblesort_size = 20 associate(X => exc%val(2, :)) ! 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(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, 1), lbound(V, 1) + 1, -1 do i = lbound(V, 1), ubound(V, 1) - 1 if (.not. V(i) <= V(i + 1)) then ! swap block n_swaps(2) = n_swaps(2) + 1 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, 1) j = lbound(B, 1) do k = lbound(C, 1), ubound(C, 1) if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then if ( A(i) <= B(j)) then C(k) = A(i) i = i + 1 else C(k) = B(j) j = j + 1 n_swaps(2) = n_swaps(2) + 1 end if else if (i <= ubound(A, 1)) then C(k) = A(i) i = i + 1 else if (j <= ubound(B, 1)) then C(k) = B(j) j = j + 1 n_swaps(2) = n_swaps(2) + 1 end if end do end associate end block end do current_size = 2 * current_size end do end associate end block even_swaps = mod(sum(n_swaps), 2) == 0 #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif end subroutine elemental function canonicalize_Excite_2_t (exc) result(res) type(Excite_2_t), intent(in) :: exc type(Excite_2_t) :: res logical :: dummy res = exc call make_canonical(res, dummy) end function elemental function defined_Excite_3_t (exc) result(res) type(Excite_3_t), intent(in) :: exc logical :: res res = all(exc%val /= UNKNOWN) end function elemental function is_sorted_Excite_3_t (exc) result(res) type(Excite_3_t), intent(in) :: exc logical :: res res = is_sorted(exc%val(1, :)) .and. is_sorted(exc%val(2, :)) end function elemental function is_canonical_Excite_3_t (exc) result(res) type(Excite_3_t), intent(in) :: exc logical :: res res = .false. if (.not. defined(exc)) return if (.not. (is_set(exc%val(1, :)) .and. is_set(exc%val(2, :)))) return if (.not. disjoint(exc%val(1, :), exc%val(2, :))) return res = .true. end function elemental function spin_allowed_Excite_3_t (exc) result(res) type(Excite_3_t), intent(in) :: exc logical :: res res = sum(G1(exc%val(1, :))%Ms) == sum(G1(exc%val(2, :))%Ms) end function pure function occupation_allowed_Excite_3_t (nI, exc) result(res) integer, intent(in) :: nI(:) type(Excite_3_t), intent(in) :: exc logical :: res debug_function_name("occupation_allowed_Excite_3_t") #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif res = subset(exc%val(1, :), nI) .and. disjoint(exc%val(2, :), nI) end function elemental subroutine make_canonical_Excite_3_t (exc, even_swaps) type(Excite_3_t), intent(inout) :: exc logical, intent(out) :: even_swaps integer :: n_swaps(2) type(Excite_3_t) :: copy_exc routine_name("make_canonical") copy_exc = exc n_swaps(1) = 0 ! merge_sort block integer, dimension(:), allocatable :: tmp integer :: current_size, left, mid, right integer, parameter :: along = 1, & bubblesort_size = 20 associate(X => exc%val(1, :)) ! 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(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, 1), lbound(V, 1) + 1, -1 do i = lbound(V, 1), ubound(V, 1) - 1 if (.not. V(i) <= V(i + 1)) then ! swap block n_swaps(1) = n_swaps(1) + 1 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, 1) j = lbound(B, 1) do k = lbound(C, 1), ubound(C, 1) if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then if ( A(i) <= B(j)) then C(k) = A(i) i = i + 1 else C(k) = B(j) j = j + 1 n_swaps(1) = n_swaps(1) + 1 end if else if (i <= ubound(A, 1)) then C(k) = A(i) i = i + 1 else if (j <= ubound(B, 1)) then C(k) = B(j) j = j + 1 n_swaps(1) = n_swaps(1) + 1 end if end do end associate end block end do current_size = 2 * current_size end do end associate end block n_swaps(2) = 0 ! merge_sort block integer, dimension(:), allocatable :: tmp integer :: current_size, left, mid, right integer, parameter :: along = 1, & bubblesort_size = 20 associate(X => exc%val(2, :)) ! 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(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, 1), lbound(V, 1) + 1, -1 do i = lbound(V, 1), ubound(V, 1) - 1 if (.not. V(i) <= V(i + 1)) then ! swap block n_swaps(2) = n_swaps(2) + 1 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, 1) j = lbound(B, 1) do k = lbound(C, 1), ubound(C, 1) if (i <= ubound(A, 1) .and. j <= ubound(B, 1)) then if ( A(i) <= B(j)) then C(k) = A(i) i = i + 1 else C(k) = B(j) j = j + 1 n_swaps(2) = n_swaps(2) + 1 end if else if (i <= ubound(A, 1)) then C(k) = A(i) i = i + 1 else if (j <= ubound(B, 1)) then C(k) = B(j) j = j + 1 n_swaps(2) = n_swaps(2) + 1 end if end do end associate end block end do current_size = 2 * current_size end do end associate end block even_swaps = mod(sum(n_swaps), 2) == 0 #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif end subroutine elemental function canonicalize_Excite_3_t (exc) result(res) type(Excite_3_t), intent(in) :: exc type(Excite_3_t) :: res logical :: dummy res = exc call make_canonical(res, dummy) end function elemental function dyn_defined(exc) result(res) class(Excitation_t), intent(in) :: exc logical :: res routine_name("dyn_defined") select type (exc) type is (Excite_0_t) res = defined(exc) type is (Excite_1_t) res = defined(exc) type is (Excite_2_t) res = defined(exc) type is (Excite_3_t) res = defined(exc) class default call stop_all(this_routine, "Excitation type invalid.") end select end function pure function from_integer_Excite_1_t(src, tgt) result(res) integer, intent(in), optional :: src, tgt type(Excite_1_t) :: res ! The values default to UNKNOWN automatically. if (present(src)) res%val(1, 1) = src if (present(tgt)) res%val(2, 1) = tgt end function pure function from_integer_Excite_2_t(src1, tgt1, src2, tgt2) result(res) integer, intent(in), optional :: src1, tgt1, src2, tgt2 type(Excite_2_t) :: res ! The values default to UNKNOWN automatically. if (present(src1)) res%val(1, 1) = src1 if (present(tgt1)) res%val(2, 1) = tgt1 if (present(src2)) res%val(1, 2) = src2 if (present(tgt2)) res%val(2, 2) = tgt2 end function pure function from_integer_Excite_3_t(src1, tgt1, src2, tgt2, src3, tgt3) result(res) integer, intent(in), optional :: src1, tgt1, src2, tgt2, src3, tgt3 type(Excite_3_t) :: res ! The values default to UNKNOWN automatically. if (present(src1)) res%val(1, 1) = src1 if (present(tgt1)) res%val(2, 1) = tgt1 if (present(src2)) res%val(1, 2) = src2 if (present(tgt2)) res%val(2, 2) = tgt2 if (present(src2)) res%val(1, 3) = src3 if (present(tgt2)) res%val(2, 3) = tgt3 end function !> Create an excitation from an excitation matrix and excitation level IC pure function create_excitation(ic, ex) result(exc) !> The excitation level. (1=Excite_1_t, 2=Excite_2_t, ...) integer, intent(in) :: IC !> An excitation matrix as in the %val component of !> the excitation types. integer, intent(in), optional :: ex(2, ic) !> An excitation of type excitation_t. !> By using select type(exc) one can select the actual type at runtime !> **and** statically dispatch as much as possible at runtime. class(Excitation_t), allocatable :: exc #ifdef DEBUG_ character(*), parameter :: this_routine = 'create_excitation' #endif select case (IC) case (0) allocate(Excite_0_t :: exc) case (1) allocate(Excite_1_t :: exc) case (2) allocate(Excite_2_t :: exc) case (3) allocate(Excite_3_t :: exc) case (4 :) allocate(Excite_Further_t :: exc) case (:-1) #ifdef DEBUG_ call stop_all(this_routine, 'invalid IC < 0 passed.') #endif end select if (present(ex)) then select type (exc) type is (Excite_1_t) exc%val = ex type is (Excite_2_t) exc%val = ex type is (Excite_3_t) exc%val = ex end select end if end function !> Create an excitation from nI to nJ where the excitation level !> is already known. subroutine get_excitation_new(nI, nJ, IC, exc, tParity) !> Two Slater determinants in nI format. integer, intent(in) :: nI(nEl), nJ(nEl) !> The excitation level. (1=Excite_1_t, 2=Excite_2_t, ...) integer, intent(in) :: IC !> An excitation of type excitation_t. !> By using select type(exc) one can select the actual type at runtime !> **and** statically dispatch as much as possible at compile time. class(Excitation_t), allocatable, intent(out) :: exc !> The parity of the excitation. logical, intent(out) :: tParity exc = create_excitation(ic) ! The compiler has to statically know, what the type of exc is. select type (exc) type is (Excite_1_t) exc%val(1, 1) = 1 call GetExcitation(nI, nJ, nel, exc%val, tParity) type is (Excite_2_t) exc%val(1, 1) = 2 call GetExcitation(nI, nJ, nel, exc%val, tParity) type is (Excite_3_t) exc%val(1, 1) = 3 call GetExcitation(nI, nJ, nel, exc%val, tParity) end select end subroutine get_excitation_new subroutine get_excitation_old(nI, nJ, ic, exc, tParity) integer, intent(in) :: nI(nEl), nJ(nEl), IC integer, intent(out) :: exc(2, maxExcit) logical, intent(out) :: tParity character(*), parameter :: this_routine = 'get_excitation_old' #ifdef DEBUG_ block use util_mod, only: stop_all use constants, only: stderr if (.not. (ic .in. [1, 2, 3])) then write(stderr, *) "" write(stderr, *) "Assertion ic .in. [1, 2, 3]" write(stderr, *) "failed in /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/excitatio& &n_types.fpp:412" call stop_all (this_routine, "Assert fail: ic .in. [1, 2, 3]") end if end block #endif exc(1, 1) = ic call GetExcitation(nI, nJ, nel, exc, tParity) end subroutine !> Create canonical excitation from ilutI to ilutJ !> where the excitation level is already known. subroutine get_bit_excitation(ilutI, ilutJ, IC, exc, tParity) !> Two Slater determinants in bitmask format. integer(kind=n_int), intent(in) :: iLutI(0:NIfTot), iLutJ(0:NIfTot) !> The excitation level. (1=Excite_1_t, 2=Excite_2_t, ...) integer, intent(in) :: IC !> The parity of the excitation. class(Excitation_t), allocatable, intent(out) :: exc logical, intent(out) :: tParity routine_name("get_bit_excitation") exc = create_excitation(ic) tParity = .false. ! The compiler has to statically know, what the type of exc is. select type (exc) type is (Excite_0_t) continue type is (Excite_1_t) exc%val(1, 1) = IC call GetBitExcitation(iLutI, iLutJ, exc%val, tParity) exc = canonicalize(exc) type is (Excite_2_t) exc%val(1, 1) = IC call GetBitExcitation(iLutI, iLutJ, exc%val, tParity) exc = canonicalize(exc) type is (Excite_3_t) exc%val(1, 1) = IC call GetBitExcitation(iLutI, iLutJ, exc%val, tParity) exc = canonicalize(exc) type is (Excite_Further_t) continue class default call stop_all(this_routine, "Excitation type invalid.") end select end subroutine get_bit_excitation pure subroutine set_last_tgt_Excite_1_t (exc, tgt) type(Excite_1_t), intent(inout) :: exc integer, intent(in) :: tgt exc%val(2, size(exc%val, 2)) = tgt end subroutine set_last_tgt_Excite_1_t pure subroutine set_last_tgt_Excite_2_t (exc, tgt) type(Excite_2_t), intent(inout) :: exc integer, intent(in) :: tgt exc%val(2, size(exc%val, 2)) = tgt end subroutine set_last_tgt_Excite_2_t pure subroutine set_last_tgt_Excite_3_t (exc, tgt) type(Excite_3_t), intent(inout) :: exc integer, intent(in) :: tgt exc%val(2, size(exc%val, 2)) = tgt end subroutine set_last_tgt_Excite_3_t pure function get_last_tgt_Excite_1_t (exc) result(res) type(Excite_1_t), intent(in) :: exc integer :: res res = exc%val(2, size(exc%val, 2)) end function get_last_tgt_Excite_1_t pure function get_last_tgt_Excite_2_t (exc) result(res) type(Excite_2_t), intent(in) :: exc integer :: res res = exc%val(2, size(exc%val, 2)) end function get_last_tgt_Excite_2_t pure function get_last_tgt_Excite_3_t (exc) result(res) type(Excite_3_t), intent(in) :: exc integer :: res res = exc%val(2, size(exc%val, 2)) end function get_last_tgt_Excite_3_t pure function excite_nI_Excite_0_t(det_I, exc) result(res) integer, intent(in) :: det_I(:) type(Excite_0_t), intent(in) :: exc integer :: res(size(det_I)) #ifdef WARNING_WORKAROUND_ associate(exc => exc); end associate #endif res = det_I end function pure function excite_nI_Excite_1_t(det_I, exc) result(res) integer, intent(in) :: det_I(:) type(Excite_1_t), intent(in) :: exc integer :: res(size(det_I)) debug_function_name('excite_nI_Excite_1_t') #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif res = special_union_complement(det_I, exc%val(2, :), exc%val(1, :)) end function pure function excite_nI_Excite_2_t(det_I, exc) result(res) integer, intent(in) :: det_I(:) type(Excite_2_t), intent(in) :: exc integer :: res(size(det_I)) debug_function_name('excite_nI_Excite_2_t') #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif res = special_union_complement(det_I, exc%val(2, :), exc%val(1, :)) end function pure function excite_nI_Excite_3_t(det_I, exc) result(res) integer, intent(in) :: det_I(:) type(Excite_3_t), intent(in) :: exc integer :: res(size(det_I)) debug_function_name('excite_nI_Excite_3_t') #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (is_canonical(exc))) then call stop_all (this_routine, "Assert fail: is_canonical(exc)") end if end block #endif res = special_union_complement(det_I, exc%val(2, :), exc%val(1, :)) end function pure function excite_SpinOrbIdx_t_Excite_0_t(det_I, exc) result(res) type(SpinOrbIdx_t), intent(in) :: det_I type(Excite_0_t), intent(in) :: exc type(SpinOrbIdx_t) :: res res%idx = excite(det_I%idx, exc) end function pure function excite_SpinOrbIdx_t_Excite_1_t(det_I, exc) result(res) type(SpinOrbIdx_t), intent(in) :: det_I type(Excite_1_t), intent(in) :: exc type(SpinOrbIdx_t) :: res res%idx = excite(det_I%idx, exc) end function pure function excite_SpinOrbIdx_t_Excite_2_t(det_I, exc) result(res) type(SpinOrbIdx_t), intent(in) :: det_I type(Excite_2_t), intent(in) :: exc type(SpinOrbIdx_t) :: res res%idx = excite(det_I%idx, exc) end function pure function excite_SpinOrbIdx_t_Excite_3_t(det_I, exc) result(res) type(SpinOrbIdx_t), intent(in) :: det_I type(Excite_3_t), intent(in) :: exc type(SpinOrbIdx_t) :: res res%idx = excite(det_I%idx, exc) end function pure function excite_Ilut_t_Excite_0_t(ilut_I, exc) result(res) integer(n_int), intent(in) :: ilut_I(:) type(Excite_0_t), intent(in) :: exc integer(n_int) :: res(0:size(ilut_I) - 1) character(*), parameter :: this_routine = 'excite_Ilut_t_Excite_0_t' integer :: src(0), tgt(0), i src = exc%val(1, :) tgt = exc%val(2, :) #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (defined(exc))) then call stop_all (this_routine, "Assert fail: defined(exc)") end if end block #endif res = ilut_I do i = 1, 0 #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (all(src(i) /= tgt))) then call stop_all (this_routine, "Assert fail: all(src(i) /= tgt)") end if end block #endif clr_orb(res, src(i)) set_orb(res, tgt(i)) end do end function pure function excite_Ilut_t_Excite_1_t(ilut_I, exc) result(res) integer(n_int), intent(in) :: ilut_I(:) type(Excite_1_t), intent(in) :: exc integer(n_int) :: res(0:size(ilut_I) - 1) character(*), parameter :: this_routine = 'excite_Ilut_t_Excite_1_t' integer :: src(1), tgt(1), i src = exc%val(1, :) tgt = exc%val(2, :) #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (defined(exc))) then call stop_all (this_routine, "Assert fail: defined(exc)") end if end block #endif res = ilut_I do i = 1, 1 #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (all(src(i) /= tgt))) then call stop_all (this_routine, "Assert fail: all(src(i) /= tgt)") end if end block #endif clr_orb(res, src(i)) set_orb(res, tgt(i)) end do end function pure function excite_Ilut_t_Excite_2_t(ilut_I, exc) result(res) integer(n_int), intent(in) :: ilut_I(:) type(Excite_2_t), intent(in) :: exc integer(n_int) :: res(0:size(ilut_I) - 1) character(*), parameter :: this_routine = 'excite_Ilut_t_Excite_2_t' integer :: src(2), tgt(2), i src = exc%val(1, :) tgt = exc%val(2, :) #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (defined(exc))) then call stop_all (this_routine, "Assert fail: defined(exc)") end if end block #endif res = ilut_I do i = 1, 2 #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (all(src(i) /= tgt))) then call stop_all (this_routine, "Assert fail: all(src(i) /= tgt)") end if end block #endif clr_orb(res, src(i)) set_orb(res, tgt(i)) end do end function pure function excite_Ilut_t_Excite_3_t(ilut_I, exc) result(res) integer(n_int), intent(in) :: ilut_I(:) type(Excite_3_t), intent(in) :: exc integer(n_int) :: res(0:size(ilut_I) - 1) character(*), parameter :: this_routine = 'excite_Ilut_t_Excite_3_t' integer :: src(3), tgt(3), i src = exc%val(1, :) tgt = exc%val(2, :) #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (defined(exc))) then call stop_all (this_routine, "Assert fail: defined(exc)") end if end block #endif res = ilut_I do i = 1, 3 #ifdef DEBUG_ block use util_mod, only: stop_all if (.not. (all(src(i) /= tgt))) then call stop_all (this_routine, "Assert fail: all(src(i) /= tgt)") end if end block #endif clr_orb(res, src(i)) set_orb(res, tgt(i)) end do end function pure function dyn_excite(det_I, exc) result(res) type(SpinOrbIdx_t), intent(in) :: det_I class(Excitation_t), intent(in) :: exc type(SpinOrbIdx_t) :: res select type (exc) type is (Excite_1_t) res = excite(det_I, exc) type is (Excite_2_t) res = excite(det_I, exc) type is (Excite_3_t) res = excite(det_I, exc) end select end function pure function dyn_nI_excite(det_I, exc) result(res) integer, intent(in) :: det_I(:) class(Excitation_t), intent(in) :: exc integer :: res(size(det_I)) select type (exc) type is (Excite_0_t) res = excite(det_I, exc) type is (Excite_1_t) res = excite(det_I, exc) type is (Excite_2_t) res = excite(det_I, exc) type is (Excite_3_t) res = excite(det_I, exc) class default call stop_all("dyn_nI_excite", "Excitation type invalid.") end select end function end module