get_offdiag_helement_heisenberg Function

    function get_offdiag_helement_heisenberg(nI, ex, tpar) result(hel)
        integer, intent(in) :: nI(nel), ex(2, 2)
        logical, intent(in) :: tpar
        HElement_t(dp) :: hel
#ifdef DEBUG_
        character(*), parameter :: this_routine = "get_offdiag_helement_heisenberg"
#endif
        integer(n_int) :: ilutI(0:NIfTot)
        integer :: src(2), tgt(2)
        real(dp) :: ni_z, nj_z, spin_fac

        src = get_src(ex)
        tgt = get_tgt(ex)

        ! oh and i have to check the occupancy of nI here or? otherwise this
        ! does not make any sense and is just based on the ex-matrix
        ! it is done in the same way in the slater condon routines, to not
        ! check occupancy.. do it in the debug mode atleast. .
        call EncodeBitDet(nI, ilutI)

        ASSERT(IsOcc(ilutI, src(1)))
        ASSERT(IsOcc(ilutI, src(2)))
        ASSERT(IsNotOcc(ilutI, tgt(1)))
        ASSERT(IsNotOcc(ilutI, tgt(2)))

        ! should i assert the same "same-spinness" here or just return
        ! zero is spins and orbitals do not fit?? i guess that would be
        ! better
        if (same_spin(src(1), src(2))) then
            hel = h_cast(0.0_dp)
        else
            ! i have to check if the orbitals fit.. ex is sorted here or?
            ! can i be sure about that?? check that!

            if (.not. (is_in_pair(src(1), tgt(1)) .and. is_in_pair(src(2), tgt(2)))) then
                hel = h_cast(0.0_dp)

            else
                ! this matrix access checks if the orbitals are connected
                hel = get_heisenberg_exchange(src(1), src(2))

            end if
        end if

        if (t_trans_corr_2body) then
            ! i just realise the sign of the spin densities around the
            ! involved orbitals depend on the spin to be flipped!
            ni_z = get_spin_density_neighbors(ilutI, gtid(src(1)))
            nj_z = get_spin_density_neighbors(ilutI, gtid(src(2)))

            ! beta is defined as negative ms!
            if (is_beta(src(1))) then
                spin_fac = ni_z - nj_z
            else
                spin_fac = nj_z - ni_z
            end if

            hel = hel * exp(2.0_dp * trans_corr_param_2body * (spin_fac - 1.0_dp))

        end if

        if (tpar) hel = -hel

    end function get_offdiag_helement_heisenberg