sum_hop_transcorr_factor_vec Function

    real(dp) function sum_hop_transcorr_factor_vec(r1, r2, r3, r4)
        ! function to perform the summation over the four hopping
        ! transcorrelation factors in the 2-body term of the hopping
        ! transcorrelated real-space hubbard model
        integer, intent(in) :: r1(3), r2(3), r3(3), r4(3)
#ifdef DEBUG_
        character(*), parameter :: this_routine = "sum_hop_transcorr_factor_vec"
#endif
        integer i, m_vec(3)
        integer :: ind_1(3), ind_2(3), ind_3(3), ind_4(3)

        ASSERT(associated(lat))

        sum_hop_transcorr_factor_vec = 0.0_dp

        if (allocated(hop_transcorr_factor_cached_vec)) then
            do i = 1, lat%get_nsites()
                m_vec = lat%get_r_vec(i)

                ind_1 = r1 - m_vec
                ind_2 = r2 - m_vec
                ind_3 = m_vec - r3
                ind_4 = m_vec - r4

                sum_hop_transcorr_factor_vec = sum_hop_transcorr_factor_vec + &
                                               hop_transcorr_factor_cached_vec(ind_1(1), ind_1(2), ind_1(3)) * &
                                               hop_transcorr_factor_cached_vec(ind_2(1), ind_2(2), ind_2(3)) * &
                                               hop_transcorr_factor_cached_vec_m(ind_3(1), ind_3(2), ind_3(3)) * &
                                               hop_transcorr_factor_cached_vec_m(ind_4(1), ind_4(2), ind_4(3))
            end do
        else

            do i = 1, lat%get_nsites()
                ! i need a routine to give me the real-space coordinates/vector
                ! and i also need to store that now!
                m_vec = lat%get_r_vec(i)
                sum_hop_transcorr_factor_vec = sum_hop_transcorr_factor_vec + &
                                               hop_transcorr_factor(trans_corr_param, r1 - m_vec) * &
                                               hop_transcorr_factor(trans_corr_param, r2 - m_vec) * &
                                               hop_transcorr_factor(-trans_corr_param, m_vec - r3) * &
                                               hop_transcorr_factor(-trans_corr_param, m_vec - r4)

            end do
        end if

    end function sum_hop_transcorr_factor_vec