setup_pchb_sampler_conditional_FullFull Subroutine

private subroutine setup_pchb_sampler_conditional_FullFull(this, rank_with_info, ij_weights)

Type Bound

PropVec_PCHB_FullFull_t

Arguments

Type IntentOptional Attributes Name
class(PropVec_PCHB_FullFull_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_FullFull(this, rank_with_info, ij_weights)
        class(PropVec_PCHB_FullFull_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
        real(dp), allocatable :: w_A(:), w_B(:), w_B_spat(:)
        integer, allocatable :: prop_vecs(:, :)

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

        ! initialize the mapping ab -> (a, b)
        abMax = fuse_symm_idx(nSpatOrbs, nSpatOrbs)
        ijMax = abMax

        this%unoccupied = UnoccupiedGetter_t(nBasis)

        call this%a_sampler%shared_alloc([int(ijMax), size(prop_vecs, 2)], nBasis, 'PCHB')
        call this%b_sampler%shared_alloc([nBasis, int(ijMax), size(prop_vecs, 2)], nBasis, '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_A(nBasis), w_B(nBasis), w_B_spat(nSpatOrbs), source=0._dp)
            allocate(ij_weights(nSpatOrbs, nSpatOrbs, size(prop_vecs, 2)), source=0._dp)
        else
            allocate(w_A(0), w_B(0), w_B_spat(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_A(:) = 0.0_dp
                    ij = fuse_symm_idx(i, j)
                    do a = 1, nSpatOrbs
                        w_B_spat(:) = 0.0_dp
                        do b = 1, nSpatOrbs
                            if (iProcIndex_intra == rank_with_info) then
                                associate(idx => ijab_Index_t([i, j, a, b]))
                                associate(exc => DistinctDouble_t(idx%normalize()))
                                if (this%indexer%is_allowed(exc, prop_vecs(:, i_sg)) &
                                        .and. symmetry_allowed(exc)) then
                                    w_B_spat(b) = this%get_weight(exc)
                                end if
                                end associate
                                end associate
                            end if
                        end do
                        if (iProcIndex_intra == root) then
                            w_B(1::2) = w_B_spat(:)
                            w_B(2::2) = w_B_spat(:)
                            w_B(2 * a - 1) = 0._dp
                        end if
                        call this%B_sampler%setup_entry(2 * a - 1, IJ, i_sg, root, w_B)
                        if (iProcIndex_intra == root) then
                            call swap(w_B(2 * a - 1), w_B(2 * a))
                        end if
                        call this%B_sampler%setup_entry(2 * a, IJ, i_sg, root, w_B)
                        if (iProcIndex_intra == root) then
                            w_A(2 * A - 1) = sum(w_B_spat)
                            w_A(2 * A) = w_A(2 * A - 1)
                        end if
                    end do
                    call this%A_sampler%setup_entry(IJ, i_sg, root, w_A)
                    if (iProcIndex_intra == root) then
                        ij_weights(i, j, i_sg) = sum(w_A)
                        ij_weights(j, i, i_sg) = ij_weights(i, j, i_sg)
                    end if
                end do
            end do
        end do
    end subroutine setup_pchb_sampler_conditional_FullFull