calc_single_overlap_mixed_ex Subroutine

private pure subroutine calc_single_overlap_mixed_ex(csf_i, csf_j, excitInfo, mat_ele, t_calc_full, rdm_ind, rdm_mat)

Arguments

Type IntentOptional Attributes Name
type(CSF_Info_t), intent(in) :: csf_i
type(CSF_Info_t), intent(in) :: csf_j
type(ExcitationInformation_t), intent(in) :: excitInfo
real(kind=dp), intent(out) :: mat_ele
logical, intent(in), optional :: t_calc_full
integer(kind=int_rdm), intent(out), optional, allocatable :: rdm_ind(:)
real(kind=dp), intent(out), optional, allocatable :: rdm_mat(:)

Contents


Source Code

    pure subroutine calc_single_overlap_mixed_ex(csf_i, csf_j, excitInfo, mat_ele, &
                                            t_calc_full, rdm_ind, rdm_mat)
        ! routine to exactly calculate the matrix element between 2 CSFs
        ! connected by a single overlap excitation with mixed generators
        type(CSF_Info_t), intent(in) :: csf_i, csf_j
        type(ExcitationInformation_t), intent(in) :: excitInfo
        HElement_t(dp), intent(out) :: mat_ele
        logical, intent(in), optional :: t_calc_full
        integer(int_rdm), intent(out), allocatable, optional :: rdm_ind(:)
        real(dp), intent(out), allocatable, optional :: rdm_mat(:)
        character(*), parameter :: this_routine = "calc_single_overlap_mixed_ex"

        real(dp) :: bVal, temp_mat, guga_mat
        HElement_t(dp) :: umat
        integer :: i, db, gen2, step1, step2, delta_b(nSpatOrbs)
        logical :: t_calc_full_

        ! in the case of rdm calculation, i know that this type of exitation
        ! only has one (or two, with switches 2-body integrals..??)
        ! rdm-contribution..

        ASSERT(present(rdm_ind) .eqv. present(rdm_mat))
        def_default(t_calc_full_, t_calc_full, .true.)

        ! set some defaults in case of early exit
        mat_ele = h_cast(0.0_dp)
        delta_b = csf_i%B_int - csf_j%B_int

        associate (ii => excitInfo%i, jj => excitInfo%j, kk => excitInfo%k, &
                   ll => excitInfo%l, st => excitInfo%fullStart, &
                   ss => excitInfo%secondStart, en => excitInfo%fullEnd, &
                   gen => excitInfo%firstGen, fe => excitInfo%firstEnd, &
                   typ => excitInfo%typ)

            if (present(rdm_ind) .and. present(rdm_mat)) then
                ! i am not sure yet if I will use symmetries in the RDM
                ! calculation (some are also left out in the SD based implo..
                ! so for now sample both combinations
                allocate(rdm_ind(1), source=0_int_rdm)
                allocate(rdm_mat(1), source=0.0_dp)
                rdm_ind(1) = contract_2_rdm_ind(ii, jj, kk, ll)
            end if

            if (any(abs(delta_b) > 1)) return

            if (t_calc_full_) then
                if (typ == excit_type%single_overlap_L_to_R) then

                    umat = (get_umat_el(fe, ss, st, en) + &
                            get_umat_el(ss, fe, en, st)) / 2.0_dp

                else if (typ == excit_type%single_overlap_R_to_L) then
                    umat = (get_umat_el(st, en, fe, ss) + &
                            get_umat_el(en, st, ss, fe)) / 2.0_dp
                else
                    call stop_all(this_routine, "shouldnt be here!")
                end if
            else
                umat = h_cast(1.0_dp)
            end if

            ! for the hamiltonian matrix element i can exit here if umat is 0
            ! but for the rdm-contribution i need to calc. the GUGA element
            if (t_calc_full .and. near_zero(umat)) return

            guga_mat = 1.0_dp

            ! i have to do the start specifically, due to the access of the
            ! single matrix elements
            step1 = csf_i%stepvector(st)
            step2 = csf_j%stepvector(st)
            bVal = csf_i%b_real(st)
            db = delta_b(st)

            guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, gen, bVal)

            if (near_zero(guga_mat)) return

            do i = st + 1, ss - 1

                step1 = csf_i%stepvector(i)
                step2 = csf_j%stepvector(i)
                bVal = csf_i%b_real(i)
                db = delta_b(i - 1)

                guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, gen, bVal)

                if (near_zero(guga_mat)) return

            end do

            step1 = csf_i%stepvector(ss)
            step2 = csf_j%stepvector(ss)
            bVal = csf_i%b_real(ss)
            db = delta_b(ss - 1)
            gen2 = -gen

            call getDoubleMatrixElement(step2, step1, db, gen, gen2, bVal, 1.0_dp, temp_mat)

            guga_mat = guga_mat * temp_mat

            if (near_zero(guga_mat)) return

            do i = ss + 1, en

                step1 = csf_i%stepvector(i)
                step2 = csf_j%stepvector(i)
                bVal = csf_i%b_real(i)
                db = delta_b(i - 1)

                guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, gen2, bVal)

                if (near_zero(guga_mat)) return

            end do

        end associate

        mat_ele = guga_mat * umat

        ! both coupling coeffs are the same
        if (present(rdm_mat)) rdm_mat = guga_mat

    end subroutine calc_single_overlap_mixed_ex