create_all_rdm_contribs Subroutine

public subroutine create_all_rdm_contribs(rdm_ind, x0, x1, rdm_ind_out, rdm_mat)


Type IntentOptional Attributes Name
integer(kind=int_rdm), intent(in) :: rdm_ind
real(kind=dp) :: x0
real(kind=dp) :: x1
integer(kind=int_rdm), intent(out), allocatable :: rdm_ind_out(:)
real(kind=dp), intent(out), allocatable :: rdm_mat(:)


Source Code

    subroutine create_all_rdm_contribs(rdm_ind, x0, x1, rdm_ind_out, rdm_mat)
        ! I also need a function which creates me the remaining
        ! index combination of not directly sampled excitation in the
        ! GUGA excitation generation
        integer(int_rdm), intent(in) :: rdm_ind
        real(dp) :: x0, x1
        integer(int_rdm), intent(out), allocatable :: rdm_ind_out(:)
        real(dp), intent(out), allocatable :: rdm_mat(:)
        integer :: ex_lvl, ex_typ, i, j, k, l

        ex_lvl = extract_excit_lvl_rdm(rdm_ind)
        ex_typ = extract_excit_type_rdm(rdm_ind)

        ASSERT(ex_lvl == 1 .or. ex_lvl == 2)
        ASSERT(ex_typ /= excit_type%invalid)

        if (ex_lvl == 1) then
            ! for singles we 'only' have the original index
            allocate(rdm_ind_out(1), source=rdm_ind)
            ! and the element is only the x0 contrib
            allocate(rdm_mat(1), source=x0)

        else if (ex_lvl == 2) then

            select case (ex_typ)

                ! we only take stuff from stochastic double excitations here..
                ! so this excludes some excitation types.. assert this below!
            case (excit_type%single_overlap_R_to_L, &
                  excit_type%single_overlap_L_to_R, &
                  excit_type%fullstop_lowering, &
                  excit_type%fullstop_raising, &
                  excit_type%fullstart_raising, &
                  excit_type%fullstart_lowering, &

                ! here we also only have the ij <-> kl conjugation which
                ! can (and should for now) be handled outside of this function!
                allocate(rdm_ind_out(1), source=rdm_ind)
                ! also x1 must be 0 in these cases
                allocate(rdm_mat(1), source=x0)

                ! group everything together where x0 is 0:
                ! the conjugated elements of those are actually dealt with
                ! in the single excitations (here they would actually be 0
                ! since we enforce a spin-flip in the overlap range!
            case (excit_type%fullstop_L_to_R, &
                  excit_type%fullstop_R_to_L, &
                  excit_type%fullstart_L_to_R, &
                  excit_type%fullstart_R_to_L, &

                allocate(rdm_ind_out(1), source=rdm_ind)
                allocate(rdm_mat(1), source=x1)

            case (excit_type%double_lowering, &

                ! here we need to do something
                call extract_2_rdm_ind(rdm_ind, i, j, k, l)
                ! a switch of j <-> l will induce a sign change between x0 and x1!
                allocate(rdm_ind_out(2), source=0_int_rdm)
                allocate(rdm_mat(2), source=0.0_dp)

                rdm_ind_out(1) = rdm_ind
                rdm_mat(1) = x0 + generator_sign(i, j, k, l) * x1

                rdm_ind_out(2) = contract_2_rdm_ind(i, l, k, j, &
                                                    excit_lvl=2, excit_typ=ex_typ)

                rdm_mat(2) = x0 + generator_sign(i, l, k, j) * x1

            case (excit_type%double_L_to_R_to_L, &
                  excit_type%double_R_to_L_to_R, &
                  excit_type%double_L_to_R, &

                ! these cases contain a non-overlap excitaion if no spin-flip
                ! in the overlap range happened (i hope for now that a
                ! zero x0 element indicates such a spin-flip.. although
                ! it should be fine, since when there is a spin-flip it
                ! is 0 definitely and even if there was, x0 is still 0 and this
                ! is all what is encoded for this flipped version!

                allocate(rdm_ind_out(2), source=0_int_rdm)
                allocate(rdm_mat(2), source=0.0_dp)

                ! the first is the original
                rdm_ind_out(1) = rdm_ind
                rdm_mat(1) = x0 + x1

                call extract_2_rdm_ind(rdm_ind, i, j, k, l)
                ! the second one does hopefully not need a type encoded
                ! in the rdm! since this is not available here!!
                rdm_ind_out(2) = contract_2_rdm_ind(i, l, k, j, excit_lvl=2)
                rdm_mat(2) = -2.0_dp * x0

            case default
                ! should not be here for now!

            end select
            ! but should actually not be here now!
            ! other wise no contributions
        end if

    end subroutine create_all_rdm_contribs