subroutine pick_random_4ind(csf_i, elec_1, elec_2, orb_1, orb_2, elecs, orbs, pgen)
type(CSF_Info_t), intent(in) :: csf_i
integer, intent(in) :: elec_1, elec_2, orb_1, orb_2
integer, intent(out) :: elecs(2), orbs(2)
real(dp), intent(out) :: pgen
real(dp) :: r
elecs = 0
orbs = 0
pgen = 1.0_dp
select case (csf_i%stepvector(elec_1))
case (1)
elecs(1) = 2 * elec_1 - 1
case (2)
elecs(1) = 2 * elec_1
case (3)
r = genrand_real2_dSFMT()
if (r < 0.5_dp) then
elecs(1) = 2 * elec_1 - 1
else
elecs(1) = 2 * elec_1
end if
pgen = 0.5_dp
end select
select case (csf_i%stepvector(elec_2))
case (1)
elecs(2) = 2 * elec_2 - 1
case (2)
elecs(2) = 2 * elec_2
case (3)
r = genrand_real2_dSFMT()
if (r < 0.5_dp) then
elecs(2) = 2 * elec_2 - 1
else
elecs(2) = 2 * elec_2
end if
pgen = pgen * 0.5_dp
end select
select case (csf_i%stepvector(orb_1))
case (0)
if (same_spin(elecs(1), elecs(2))) then
if (is_beta(elecs(1))) then
orbs(1) = 2 * orb_1 - 1
else
orbs(1) = 2 * orb_1
end if
else
r = genrand_real2_dSFMT()
if (r < 0.5_dp) then
orbs(1) = 2 * orb_1 - 1
else
orbs(1) = 2 * orb_1
end if
pgen = 0.5_dp * pgen
end if
case (1)
orbs(1) = 2 * orb_1
case (2)
orbs(1) = 2 * orb_1 - 1
end select
! i think i can restrict the last one or?..
if (same_spin(elecs(1), elecs(2))) then
if (is_beta(elecs(1))) then
if (csf_i%stepvector(orb_2) == 1) then
pgen = 0.0_dp
elecs = 0
orbs = 0
return
else
orbs(2) = 2 * orb_2 - 1
end if
else
if (csf_i%stepvector(orb_2) == 2) then
pgen = 0.0_dp
elecs = 0
orbs = 0
return
else
orbs(2) = 2 * orb_2
end if
end if
else
if (is_beta(orbs(1))) then
if (csf_i%stepvector(orb_2) == 2) then
pgen = 0.0_dp
elecs = 0
orbs = 0
return
else
orbs(2) = 2 * orb_2
end if
else
if (csf_i%stepvector(orb_2) == 1) then
pgen = 0.0_dp
elecs = 0
orbs = 0
return
else
orbs(2) = 2 * orb_2 - 1
end if
end if
end if
end subroutine pick_random_4ind