setup_pchb_sampler_conditional_FastFast Subroutine

private subroutine setup_pchb_sampler_conditional_FastFast(this, rank_with_info, ij_weights)

Type Bound

PropVec_PCHB_FastFast_t

Arguments

Type IntentOptional Attributes Name
class(PropVec_PCHB_FastFast_t), intent(inout) :: this
integer, intent(in) :: rank_with_info
real(kind=dp), intent(out), allocatable :: ij_weights(:,:,:)

Source Code

    subroutine setup_pchb_sampler_conditional_FastFast(this, rank_with_info, ij_weights)
        class(PropVec_PCHB_FastFast_t), intent(inout) :: this
        integer, intent(in) :: rank_with_info
        real(dp), allocatable, intent(out) :: ij_weights(:, :, :)
        integer :: i, j, ij, a, b, ab, i_sg
        integer(int64) :: ijMax, abMax, memCost
        integer, allocatable :: prop_vecs(:, :)
        real(dp), allocatable :: w(:)

        ! possible prop_vecs
        prop_vecs = this%indexer%get_prop_vecs()

        ! initialize the mapping ab -> (a, b)
        abMax = fuse_symm_idx(nSpatOrbs, nSpatOrbs)
        ijMax = abMax
        allocate(this%tgtOrbs(2, 0:abMax), source=0)
        do a = 1, nSpatOrbs
            do b = 1, a
                ab = fuse_symm_idx(a, b)
                this%tgtOrbs(1, ab) = b
                this%tgtOrbs(2, ab) = a
            end do
        end do


        calculate_memory_demand: block
            integer(int64) :: number_of_fused_indices, bytes_per_sampler, n_prop_vec
            number_of_fused_indices = int(ijMax, int64)
            bytes_per_sampler = (int(abMax, int64) * 3_int64 * 8_int64)
            n_prop_vec = size(prop_vecs, 2, kind=int64)
            memCost = number_of_fused_indices * bytes_per_sampler * n_prop_vec
        endblock calculate_memory_demand

        write(stdout, *) "Excitation generator requires", real(memCost, dp) / 2.0_dp**30, "GB of memory"
        write(stdout, *) "The number of prop_vecs is", size(prop_vecs, 2)
        write(stdout, *) "Generating samplers for PCHB excitation generator"
        write(stdout, *) "Depending on the number of prop_vecs this can take up to 10min."
        call this%alias_sampler%shared_alloc([int(ijMax), size(prop_vecs, 2)], int(abMax), 'PCHB')
        call this%new_info_table(int(ijMax, int64), int(abMax, int64))

        block
            type(ByteSize_t) :: memory_demand
            memory_demand = this%get_memory_demand()
            write(stdout, *) "The number of prop_vecs is", size(prop_vecs, 2)
            write(stdout, *) "The excitation generator requires roughly: ", &
                memory_demand%to_human_str(binary=.true.), ' or ', memory_demand%to_human_str()
        end block


        if (iProcIndex_intra == rank_with_info) then
            allocate(w(abMax), source=0._dp)
            allocate(ij_weights(nSpatOrbs, nSpatOrbs, size(prop_vecs, 2)), source=0._dp)
        else
            allocate(w(0))
            allocate(ij_weights(0, 0, 0))
        end if

        ! The info for double excitations is independent of the super group
        do i = 1, nSpatOrbs
            do j = i, nSpatOrbs
                ij = fuse_symm_idx(i, j)
                do a = 1, nSpatOrbs
                    do b = a, nSpatOrbs
                        ab = fuse_symm_idx(a, b)
                        call this%set_info_entry(ij, ab, get_info(ijab_Index_t([i, j, a, b])))
                    end do
                end do
            end do
        end do

        do i_sg = 1, size(prop_vecs, 2)
            do i = 1, nSpatOrbs
                do j = i, nSpatOrbs
                    w = 0.0_dp
                    ij = fuse_symm_idx(i, j)
                    do a = 1, nSpatOrbs
                        do b = a, nSpatOrbs
                            ab = fuse_symm_idx(a, b)
                            if (iProcIndex_intra == rank_with_info) then
                                associate(exc => DistinctDouble_t(ijab_Index_t([i, j, a, b])))
                                if (this%indexer%is_allowed(exc, prop_vecs(:, i_sg)) .and. symmetry_allowed(exc)) then
                                    w(ab) = this%get_weight(exc)
                                end if
                                end associate
                            end if
                        end do
                    end do
                    if (iProcIndex_intra == root) then
                        ij_weights(i, j, i_sg) = sum(w)
                        ij_weights(j, i, i_sg) = ij_weights(i, j, i_sg)
                    end if
                    call this%alias_sampler%setup_entry(ij, i_sg, rank_with_info, w)
                end do
            end do
        end do
    end subroutine setup_pchb_sampler_conditional_FastFast