function get_diag_helement_heisenberg(nI) result(hel)
integer, intent(in) :: nI(nel)
HElement_t(dp) :: hel
#ifdef DEBUG_
character(*), parameter :: this_routine = "get_diag_helement_heisenberg"
#endif
integer :: i, j, src, flip
integer, allocatable :: spin_neighbors(:)
integer(n_int) :: ilut(0:NIfTot)
! here i have to sum over all the nearest neigbhor pairs and
! find the spin alignment, if it is parallel or not..
! in the heisenberg case the contribution form n_i n_j form
! nearest neighbors is always the same for every determinant nI
! but for the tJ model with less than half-filling this quantitiy is
! not a constant.. so i should differentiate in here between tJ
! and heisenberg models ..
! and what about double counting? should i loop over all the
! orbitals and their neighbors or must i differentiate between already
! counted contributions?
! it is easier to check occupancy in the ilut format
ASSERT(associated(lat))
call EncodeBitDet(nI, ilut)
hel = h_cast(0.0_dp)
do i = 1, nel
src = ni(i)
spin_neighbors = lat%get_spinorb_neighbors(src)
if (is_beta(src)) then
flip = +1
else
flip = -1
end if
do j = 1, size(spin_neighbors)
if (IsOcc(ilut, spin_neighbors(j))) then
! then it is same spin
! but i really think that we have to take the
! occupancy into account
! and in the tJ model this cancels for parallel spin
if (t_heisenberg_model) then
hel = hel + h_cast(exchange_j / 4.0_dp)
end if
else if (IsOcc(ilut, spin_neighbors(j) + flip)) then
if (t_heisenberg_model) then
hel = hel - h_cast(exchange_j / 4.0_dp)
else
hel = hel - h_cast(exchange_j / 2.0_dp)
end if
! it can be empty too, then there is no contribution
end if
end do
end do
! i think i would double count if i do not divide by 2..
! i hope i am right..
hel = hel / 2.0_dp
end function get_diag_helement_heisenberg