subroutine createStochasticExcitation_single(ilut, nI, csf_i, exc, pgen)
! calculate one possible single excitation and the corresponding
! probabilistic weight and hamilton matrix element for a given CSF
! store matrix element in ilut for now... maybe change that later
integer(n_int), intent(in) :: ilut(0:nifguga)
integer, intent(in) :: nI(nel)
type(CSF_Info_t), intent(in) :: csf_i
integer(n_int), intent(out) :: exc(0:nifguga)
real(dp), intent(out) :: pgen
character(*), parameter :: this_routine = "createStochasticExcitation_single"
type(ExcitationInformation_t) :: excitInfo
real(dp) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs), &
branch_pgen, orb_pgen, temp_pgen
HElement_t(dp) :: integral, mat_ele
type(WeightObj_t) :: weights
integer :: iO, st, en, step, i, j, gen, deltaB, step2
integer(n_int) :: ilutI(0:niftot), ilutJ(0:niftot)
ASSERT(isProperCSF_ilut(ilut))
! first pick possible orbitals:
! todo: combine that with integrals elements... otherwise it does not
! make too much sense
! but since there is a sum of contribution, which can lead to the
! same excitaiton shouldn't i check if the whole sum is zero, and not
! only the one-particle matrix element.
! this pickOrbitals picker is the only difference in the symmetric
! and non-symmetric excitation generator
! so in the initialize function point a general orbital picker to
! the specific one, depending on the input..
call pickOrbitals_single(ilut, nI, csf_i, excitInfo, orb_pgen)
if (.not. excitInfo%valid) then
! if no valid indices were picked, return 0 excitation and return
exc = 0_n_int
pgen = 0.0_dp
return
end if
if (t_guga_back_spawn) then
! do smth like this:
! if i find to increase the excit-lvl with the chosen
! orbitals and the current CSF is a non-initiator ->
! perform a crude excitation
if (increase_ex_levl(csf_i, excitInfo) .and. .not. is_init_guga) then
if (t_guga_back_spawn_trunc) then
! a 2 indicated we want to cancel excit-lvl increasing
! excitations.
pgen = 0.0_dp
exc = 0_n_int
return
end if
call create_crude_guga_single(ilut, nI, csf_i, exc, branch_pgen, excitInfo)
! there is also this routine I already wrote:
! I should combine those two as they do the same job
pgen = orb_pgen * branch_pgen
return
end if
end if
! do the crude approximation here for now..
if (tgen_guga_crude .and. .not. tgen_guga_mixed) then
call create_crude_single(ilut, csf_i, exc, branch_pgen, excitInfo)
if (near_zero(branch_pgen)) then
exc = 0_n_int
pgen = 0.0_dp
return
end if
call convert_ilut_toNECI(ilut, ilutI)
call convert_ilut_toNECI(exc, ilutJ)
call calc_guga_matrix_element(ilutI, csf_i, ilutJ, CSF_Info_t(ilutJ), excitInfo, mat_ele, .true.)
if (near_zero(mat_ele)) then
exc = 0
pgen = 0.0_dp
return
end if
call encode_matrix_element(exc, 0.0_dp, 2)
call encode_matrix_element(exc, mat_ele, 1)
pgen = orb_pgen * branch_pgen
return
end if
! have to have these short variable names or otherwise compilation
! fails due to line-length restrictions with is Three(), etc. macros
st = excitInfo%fullStart
en = excitInfo%fullEnd
! reimplement it from scratch
i = excitInfo%i
j = excitInfo%j
! first the "normal" contribution
! not sure if i have to subtract that element or not...
integral = getTmatEl(2 * i, 2 * j)
! ! then calculate the remaing switche given indices
call calcRemainingSwitches_excitInfo_single(csf_i, excitInfo, posSwitches, negSwitches)
! intitialize the weights
weights = init_singleWeight(csf_i, excitInfo%fullEnd)
! create the start randomly(if multiple possibilities)
! create the start in such a way to use it for double excitations too
call createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
negSwitches, exc, branch_pgen)
! can it be zero here? maybe due to matrix element issues...
if (near_zero(branch_pgen)) then
pgen = 0.0_dp
exc = 0_n_int
return
end if
! update at last...
gen = excitInfo%currentGen
! then do the stochastic updates..
do iO = st + 1, en - 1
! need the ongoing deltaB value to access the multFactor table in
! the same way as single and double excitations..
deltaB = getDeltaB(exc)
call singleStochasticUpdate(ilut, csf_i, iO, excitInfo, weights, posSwitches, &
negSwitches, exc, temp_pgen)
branch_pgen = branch_pgen * temp_pgen
! also get the double contribution during this loop
! depends on both stepvalues...
if (t_trunc_guga_pgen .or. &
(t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
if (branch_pgen < trunc_guga_pgen) then
pgen = 0.0_dp
exc = 0_n_int
return
end if
end if
! how to modifiy the values?
! one of it is just additional
if (.not. (treal .or. t_new_real_space_hubbard .or. t_mixed_hubbard &
.or. t_tJ_model)) then
integral = integral + get_umat_el(i, iO, j, iO) * csf_i%Occ_real(iO)
! but the r_k part confuses me a bit ...
step = csf_i%stepvector(iO)
step2 = getStepvalue(exc, iO)
integral = integral + get_umat_el(i, iO, iO, j) * &
getDoubleContribution(step2, step, deltaB, gen, csf_i%B_real(iO))
end if
end do
! the end step should be easy in this case. since due to the
! correct use of probabilistic weights only valid excitations should
! come to this point
call singleStochasticEnd(csf_i, excitInfo, exc)
! maybe but a check here if the matrix element anyhow turned out zero
if (near_zero(extract_matrix_element(exc, 1))) then
pgen = 0.0_dp
exc = 0_n_int
return
end if
pgen = orb_pgen * branch_pgen
! do all the integral calulation afterwards..
! since it depends on the created excitation.
! updates integral variable:
! should i intermediately but a if-statement for the real-space
! hubbard model here? or should i write a more efficient
! single-excitation creator? i guess i should..
! and also for the matrix element calculation maybe..
if (.not. (treal .or. t_new_real_space_hubbard .or. &
t_heisenberg_model .or. t_tJ_model .or. t_mixed_hubbard)) then
! TODO(@Oskar): Don't calculate here
call calc_integral_contribution_single(csf_i, CSF_Info_t(exc), i, j, st, en, integral)
end if
if (tFillingStochRDMOnFly) then
! if we want to do 'fast' GUGA RDMs we need to store the
! rdm index and the x0 (for singles here) coupling coefficient
! as part of the ilut(0:nifguga). this also necessitates
! a 'longer' nifguga (+3 i think for rdm_index and x0 and x1..)
! with an accompanying change to niftot.. (this could get a bit
! messy in the rest of the code..)
call encode_stochastic_rdm_info(GugaBits, exc, rdm_ind= &
contract_1_rdm_ind(i, j, excit_lvl=1, excit_typ=excit_type%single), &
x0=extract_matrix_element(exc, 1), x1=0.0_dp)
end if
call encode_matrix_element(exc, 0.0_dp, 2)
call update_matrix_element(exc, integral, 1)
if (near_zero(extract_matrix_element(exc, 1))) then
pgen = 0.0_dp
exc = 0_n_int
end if
! store the most recent excitation information
global_excitInfo = excitInfo
end subroutine createStochasticExcitation_single