subroutine create_cum_list_tJ_model(ilutI, src, neighbors, cum_arr, cum_sum, &
ic_list, tgt, cpt)
integer(n_int), intent(in) :: ilutI(0:NIfTot)
integer, intent(in) :: src, neighbors(:)
real(dp), intent(out), allocatable :: cum_arr(:)
real(dp), intent(out) :: cum_sum
integer, intent(out), allocatable :: ic_list(:)
integer, intent(in), optional :: tgt
real(dp), intent(out), optional :: cpt
#ifdef DEBUG_
character(*), parameter :: this_routine = "create_cum_list_tJ_model"
#endif
integer :: i, nI(nel), temp_ex(2, maxExcit)
integer, allocatable :: single_excits(:)
integer, allocatable :: spin_flips(:)
real(dp) :: elem
logical :: t_single, t_flip, t_single_possible, t_flip_possible
ASSERT(IsOcc(ilutI, src))
ASSERT(allocated(exchange_matrix))
allocate(cum_arr(size(neighbors)))
allocate(ic_list(size(neighbors)))
allocate(single_excits(size(neighbors)))
allocate(spin_flips(size(neighbors)))
cum_arr = 0
cum_sum = 0.0_dp
ic_list = 0
call decode_bit_det(nI, ilutI)
temp_ex(1, 1) = src
if (is_beta(src)) then
single_excits = 2 * neighbors - 1
spin_flips = 2 * neighbors
! fill in the corresponding alpha orbital
temp_ex(2, 1) = src + 1
else
single_excits = 2 * neighbors
spin_flips = 2 * neighbors - 1
! fill in the beta orbital
temp_ex(2, 1) = src - 1
end if
if (present(tgt)) then
t_single = .false.
t_flip = .false.
! find the probability of choosing orbital target
if (is_beta(src) .eqv. is_beta(tgt)) then
! then it was definetly a single excitation
t_single = .true.
else
t_flip = .true.
end if
ASSERT(present(cpt))
cpt = 0.0_dp
do i = 1, ubound(neighbors, 1)
elem = 0.0_dp
t_single_possible = .false.
t_flip_possible = .false.
if (IsNotOcc(ilutI, single_excits(i)) .and. &
IsNotOcc(ilutI, spin_flips(i))) then
! just to be sure use the tmat, so both orbitals are
! definetly connected
elem = abs(get_offdiag_helement_tJ(nI, [src, single_excits(i)], .false.))
! elem = abs(GetTMatEl(src, single_excits(i)))
t_single_possible = .true.
else if (IsOcc(ilutI, spin_flips(i)) .and. &
IsNotOcc(ilutI, single_excits(i))) then
temp_ex(1, 2) = spin_flips(i)
temp_ex(2, 2) = single_excits(i)
elem = abs(get_offdiag_helement_heisenberg(nI, temp_ex, .false.))
! elem = abs(get_heisenberg_exchange(src, spin_flips(i)))
t_flip_possible = .true.
end if
cum_sum = cum_sum + elem
if (t_single .and. t_single_possible .and. tgt == single_excits(i)) then
cpt = elem
else if (t_flip .and. t_flip_possible .and. tgt == spin_flips(i)) then
cpt = elem
end if
end do
if (cum_sum < EPS) then
cpt = 0.0_dp
else
cpt = cpt / cum_sum
end if
else
! create the list depending on the possible excitations
do i = 1, ubound(neighbors, 1)
elem = 0.0_dp
if (IsNotOcc(ilutI, single_excits(i)) .and. &
IsNotOcc(ilutI, spin_flips(i))) then
! then the orbital is empty an we can do a hopping
! reuse the hubbard matrix elements..
! or just the -t element?
! since we only need absolute value of matrix element
! it would be better to just use -t .. anyway.. keep it
! general
elem = abs(get_offdiag_helement_tJ(nI, [src, single_excits(i)], .false.))
! elem = abs(GetTMatEl(src, single_excits(i)))
ic_list(i) = 1
else if (IsOcc(IlutI, spin_flips(i)) .and. &
IsNotOcc(IlutI, single_excits(i))) then
! then we can do a spin flip
temp_ex(1, 2) = spin_flips(i)
temp_ex(2, 2) = single_excits(i)
elem = abs(get_offdiag_helement_heisenberg(nI, temp_ex, .false.))
! elem = abs(get_heisenberg_exchange(src, spin_flips(i)))
ic_list(i) = 2
else
! if the spin-parallel is occupied, no exciation
! possible, and also prohibit double occupancies
elem = 0.0_dp
end if
cum_sum = cum_sum + elem
cum_arr(i) = cum_sum
end do
end if
end subroutine create_cum_list_tJ_model