subroutine gen_excit_tJ_model(nI, ilutI, nJ, ilutJ, exFlag, ic, &
ex, tParity, pGen, hel, store, run)
integer, intent(in) :: nI(nel), exFlag
integer(n_int), intent(in) :: ilutI(0:NIfTot)
integer, intent(out) :: nJ(nel), ic, ex(2, maxExcit)
integer(n_int), intent(out) :: ilutJ(0:NifTot)
real(dp), intent(out) :: pGen
logical, intent(out) :: tParity
HElement_t(dp), intent(out) :: hel
type(excit_gen_store_type), intent(inout), target :: store
integer, intent(in), optional :: run
character(*), parameter :: this_routine = "gen_excit_tJ_model"
integer :: elec, src, id, ind, elec_2, tgt_1, tgt_2
real(dp) :: p_elec, p_orb, cum_sum, r, cum_sum_opp, cpt_opp
integer, allocatable :: neighbors(:), ic_list(:), tmp_ic_list(:)
real(dp), allocatable :: cum_arr(:), cum_arr_opp(:)
! the idea for the tJ excitation generator on a lattice is to
! still pick the electron at random and then check for its
! neighbors, if the neighboring site is empty, do a hop
! and if there is a electron of opposite spin do a spin-flip
! if it is occupied by an electron of same spin no excitation with
! this neighbor is possible
unused_var(exFlag)
unused_var(store)
unused_var(run)
hel = h_cast(0.0_dp)
! use the lattice type like in the real-space hubbard implementation
ASSERT(associated(lat))
#ifdef WARNING_WORKAROUND_
hel = 0.0_dp
if (present(run)) then
unused_var(run)
end if
#endif
unused_var(store)
unused_var(exflag)
elec = 1 + int(genrand_real2_dsfmt() * nel)
p_elec = 1.0_dp / real(nel, dp)
src = nI(elec)
id = gtid(src)
neighbors = lat%get_neighbors(id)
call create_cum_list_tJ_model(ilutI, src, neighbors, cum_arr, cum_sum, &
ic_list)
if (cum_sum < EPS) then
nJ(1) = 0
pgen = 0.0_dp
return
end if
r = genrand_real2_dsfmt() * cum_sum
ind = binary_search_first_ge(cum_arr, r)
! i just realised that for spin-flip excitation we have to take the
! opposite order of orbital picking into account too..
! because the same excitation could have been picked with the spin
! opposite electron in the neighborhood of the first electron too..
ic = ic_list(ind)
! for spin-flips i have to add a contribution down below!
if (ind == 1) then
p_orb = cum_arr(1) / cum_sum
else
p_orb = (cum_arr(ind) - cum_arr(ind - 1)) / cum_sum
end if
if (ic == 1) then
if (is_beta(src)) then
tgt_1 = 2 * neighbors(ind) - 1
else
tgt_1 = 2 * neighbors(ind)
end if
call make_single(nI, nJ, elec, tgt_1, ex, tParity)
ilutJ = make_ilutJ(ilutI, ex, 1)
else if (ic == 2) then
! here i have to recalc the contribution if i would have picked
! the electron in orbital spin_orb first
! but i made some assumptions about the order of the picked
! electrons and holes, which is not valid to do..
if (is_beta(src)) then
! need to the the index of electron 2 in nI
! the second electron must be alpha
elec_2 = find_elec_in_ni(nI, 2 * neighbors(ind))
! we need the orbital alpha of src
! and the beta of the second orbital
tgt_1 = get_alpha(src)
tgt_2 = 2 * neighbors(ind) - 1
else
! v.v here
elec_2 = find_elec_in_ni(nI, 2 * neighbors(ind) - 1)
tgt_1 = get_beta(src)
tgt_2 = 2 * neighbors(ind)
end if
! the idea is to target the spin-orbital of the other electron!
! (the first in this case!
call create_cum_list_tJ_model(ilutI, nI(elec_2), &
lat%get_neighbors(gtid(nI(elec_2))), cum_arr_opp, cum_sum_opp, &
tmp_ic_list, src, cpt_opp)
p_orb = p_orb + cpt_opp
call make_double(nI, nJ, elec, elec_2, tgt_1, tgt_2, ex, tParity)
ilutJ = make_ilutJ(ilutI, ex, 2)
else
! something went wrong..
call stop_all(this_routine, &
"something went wrong ic > 2!")
end if
pgen = p_elec * p_orb
end subroutine gen_excit_tJ_model