subroutine pickOrbs_sym_uniform_mol_single(ilut, nI, csf_i, excitInfo, pgen)
! new implementation to pick single orbitals, more similar to the
! other neci implementations
! with this new looping over other orbitals it will probably also
! be easier to include the neccesarry switch conditions!
! this also applies for double excitations!! -> think about that !
integer(n_int), intent(in) :: ilut(0:nifguga)
integer, intent(in) :: nI(nel)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(out) :: excitInfo
real(dp), intent(out) :: pgen
character(*), parameter :: this_routine = "pickOrbs_sym_uniform_mol_single"
integer :: elec, cc_i, ierr, nOrb, orb_ind, orb_i, orb_a
real(dp), allocatable :: cum_arr(:)
real(dp) :: cum_sum, r, elec_factor
unused_var(ilut)
! first pick completely random from electrons only!
elec = 1 + floor(genrand_real2_dSFMT() * nEl)
! have to adjust pgen if it is a doubly occupied orbital afterwards
! -> since twice the chance to pick that orbital then!
! pick associated "spin orbital"
orb_i = nI(elec)
! get the symmetry index:
! since there is no spin restriction here have to consider both
! again
cc_i = ClassCountInd(1, SpinOrbSymLabel(orb_i), G1(orb_i)%Ml)
! get the number of orbitals in this symmetry sector
nOrb = OrbClassCount(cc_i)
allocate(cum_arr(nOrb), stat=ierr)
! actually only need one of the symmetry list, but then just have to
! only check if one of the two spin-orbitals is empty
! now have to have specific picking routines depending on the
! stepvalue of the already picked orbital i
! TODO: hm -> for the singles the matrix element gets calculated
! exactly! -> thats NOT EASY in the guga case!
! at least not as easy as in the determinant case!!
! write an email to simon and ask ali if that makes any sense then
! eg. to use the one particle elements as an approximation..
select case (csf_i%stepvector(gtID(orb_i)))
! der stepvalue sagt mir auch, ob es ein alpha oder beta
! elektron war..
case (1)
elec_factor = 1.0_dp
call gen_cum_list_guga_single_1(nI, csf_i, orb_i, cc_i, cum_arr)
case (2)
! to do
elec_factor = 1.0_dp
call gen_cum_list_guga_single_2(nI, csf_i, orb_i, cc_i, cum_arr)
case (3)
! adjust pgen, the chance to pick a doubly occupied with
! spinorbitals is twice as high..
elec_factor = 2.0_dp
call gen_cum_list_guga_single_3(nI, csf_i, orb_i, cc_i, cum_arr)
case default
call stop_all(this_routine, "should not have picked empty orbital")
end select
! assign the spatial orbital:
orb_i = gtID(orb_i)
! get the orbital
cum_sum = cum_arr(nOrb)
if (near_zero(cum_sum)) then
orb_a = 0
excitInfo%valid = .false.
return
else
r = genrand_real2_dSFMT() * cum_sum
orb_ind = binary_search_first_ge(cum_arr, r)
orb_a = sym_label_list_spat(SymLabelCounts2(1, cc_i) + orb_ind - 1)
if (orb_ind == 1) then
pgen = cum_arr(1) / cum_sum
else
pgen = (cum_arr(orb_ind) - cum_arr(orb_ind - 1)) / cum_sum
end if
end if
deallocate(cum_arr)
ASSERT(orb_a /= orb_i)
! assign excitInfo and calc. final pgen
pgen = pgen * elec_factor / real(nEl, dp)
if (orb_a < orb_i) then
! raising generator
excitInfo = assign_excitInfo_values_single(gen_type%R, orb_a, orb_i, orb_a, orb_i)
else
! lowering generator
excitInfo = assign_excitInfo_values_single(gen_type%L, orb_a, orb_i, orb_i, orb_a)
end if
end subroutine pickOrbs_sym_uniform_mol_single