subroutine create_random_spin_orbs(ilut, csf_i, excitInfo, elecs, orbs, pgen)
! a subroutine to create random spin-orbitals from chosen
! spatial orbital of a GUGA excitation.
! this is needed in the crude back-spawn approximation to
! create determinant-like excitations
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(in) :: excitInfo
integer, intent(out) :: elecs(2)
integer, intent(out) :: orbs(2)
real(dp), intent(out) :: pgen
character(*), parameter :: this_routine = "create_random_spin_orbs"
integer :: elec_1, orb_1, elec_2, orb_2, d_elec, s_elec, d_orb, s_orb
real(dp) :: r
! i essentially need electron and empty spin-orbitals to
! create the excitation. the validity of the CSF will be checked
! afterwards and then maybe thrown away..
elecs = 0
orbs = 0
pgen = 1.0_dp
if (excitInfo%typ == excit_type%single) then
! the case of a single excitation
elec_1 = excitInfo%j
orb_1 = excitInfo%i
select case (csf_i%stepvector(elec_1))
case (0)
call stop_all(this_routine, "empty elec index")
case (1)
! 1 corresponds to a beta orbital
elecs(1) = 2 * elec_1 - 1
orbs(1) = 2 * orb_1 - 1
case (2)
! 2 corresponds to alpha orbitals
elecs(1) = 2 * elec_1
orbs(1) = 2 * orb_1
case (3)
select case (csf_i%stepvector(orb_1))
case (0)
! now we can have two possibilities
r = genrand_real2_dSFMT()
if (r < 0.5_dp) then
elecs(1) = 2 * elec_1 - 1
orbs(1) = 2 * orb_1 - 1
else
elecs(1) = 2 * elec_1
orbs(1) = 2 * orb_1
end if
pgen = 0.5_dp
case (1)
elecs(1) = 2 * elec_1
orbs(1) = 2 * orb_1
case (2)
elecs(1) = 2 * elec_1 - 1
orbs(1) = 2 * orb_1 - 1
end select
end select
else
elec_1 = excitInfo%j
elec_2 = excitInfo%l
orb_1 = excitInfo%i
orb_2 = excitInfo%k
ASSERT(csf_i%stepvector(elec_1) /= 0)
ASSERT(csf_i%stepvector(elec_2) /= 0)
ASSERT(csf_i%stepvector(orb_1) /= 3)
ASSERT(csf_i%stepvector(orb_2) /= 3)
select case (excitInfo%typ)
case (excit_type%single_overlap_L_to_R, &
excit_type%fullstop_lowering, &
excit_type%fullstart_raising)
! here i know the orbital indices are identical
ASSERT(orb_1 == orb_2)
orbs(1) = 2 * orb_1 - 1
orbs(2) = 2 * orb_1
! write a general function which gives me valid spin-orbs
! for GUGA CSFs (mayb as an input the 'neutral' number 3 here.c.
! maybe later..
if ((csf_i%stepvector(elec_1) == csf_i%stepvector(elec_2)) &
.and. csf_i%Occ_int(elec_1) == 1) then
! in this case no crude excitation is possible
pgen = 0.0_dp
elecs = 0
orbs = 0
return
end if
select case (csf_i%stepvector(elec_1))
case (0)
call stop_all(this_routine, "something wrong happened")
case (1)
elecs(1) = 2 * elec_1 - 1
! then the second must be of 'opposite' spin
elecs(2) = 2 * elec_2
case (2)
elecs(1) = 2 * elec_1
elecs(2) = 2 * elec_2 - 1
case (3)
if (csf_i%stepvector(elec_2) == 3) then
! here we have a choice
r = genrand_real2_dSFMT()
if (r < 0.5_dp) then
elecs(1) = 2 * elec_1 - 1
elecs(2) = 2 * elec_2
else
elecs(1) = 2 * elec_1
elecs(2) = 2 * elec_2 - 1
end if
pgen = 0.5_dp
else if (csf_i%stepvector(elec_2) == 1) then
elecs(2) = 2 * elec_2 - 1
elecs(1) = 2 * elec_1
else if (csf_i%stepvector(elec_2) == 2) then
elecs(2) = 2 * elec_2
elecs(1) = 2 * elec_1 - 1
end if
end select
case (excit_type%single_overlap_R_to_L, &
excit_type%fullstop_raising, &
excit_type%fullstart_lowering)
! here i know the electron indices are the same
ASSERT(elec_1 == elec_2)
elecs(1) = 2 * elec_1 - 1
elecs(2) = 2 * elec_1
if ((csf_i%stepvector(orb_1) == csf_i%stepvector(orb_2)) &
.and. csf_i%Occ_int(orb_1) == 1) then
pgen = 0
elecs = 0
orbs = 0
return
end if
select case (csf_i%stepvector(orb_1))
case (3)
call stop_all(this_routine, "something went wrong")
case (1)
orbs(1) = 2 * orb_1
orbs(2) = 2 * orb_2 - 1
case (2)
orbs(1) = 2 * orb_1 - 1
orbs(2) = 2 * orb_2
case (0)
if (csf_i%stepvector(orb_2) == 0) then
r = genrand_real2_dSFMT()
if (r < 0.5_dp) then
orbs(1) = 2 * orb_1 - 1
orbs(2) = 2 * orb_2
else
orbs(1) = 2 * orb_1
orbs(2) = 2 * orb_2 - 1
end if
pgen = 0.5_dp
else if (csf_i%stepvector(orb_2) == 1) then
orbs(2) = 2 * orb_2
orbs(1) = 2 * orb_1 - 1
else if (csf_i%stepvector(orb_2) == 2) then
orbs(2) = 2 * orb_2 - 1
orbs(1) = 2 * orb_1
end if
end select
! do the same as above, just for two hole indices!
case (excit_type%fullstop_L_to_R, &
excit_type%fullstop_R_to_L, &
excit_type%fullStart_L_to_R, &
excit_type%fullstart_R_to_L)
! here i know one electron and hole index are the same
ASSERT(elec_1 /= elec_2)
ASSERT(orb_1 /= orb_2)
! this case is very restrictive..
if (elec_1 == orb_1) then
s_elec = elec_1
s_orb = orb_1
d_elec = elec_2
d_orb = orb_2
else if (elec_1 == orb_2) then
s_elec = elec_1
s_orb = orb_1
d_elec = elec_2
d_orb = orb_1
else if (elec_2 == orb_1) then
s_elec = elec_2
s_orb = orb_1
d_elec = elec_1
d_orb = orb_2
else if (elec_2 == orb_2) then
s_elec = elec_2
s_orb = orb_2
d_elec = elec_1
d_orb = orb_1
#ifdef DEBUG_
else
call stop_all(this_routine, "something went wrong")
#endif
end if
ASSERT(csf_i%Occ_int(s_elec) == 1)
if (csf_i%stepvector(s_elec) == 1) then
if (IsOcc(ilut, 2 * d_elec) .and. IsNotOcc(ilut, 2 * d_orb - 1)) then
! this is the only case this is possible
elecs(1) = 2 * s_elec - 1
elecs(2) = 2 * d_elec
orbs(1) = 2 * s_orb
orbs(2) = 2 * d_orb - 1
else
pgen = 0.0_dp
end if
else if (csf_i%stepvector(s_elec) == 2) then
if (IsOcc(ilut, 2 * d_elec - 1) .and. IsNotOcc(ilut, 2 * d_orb)) then
elecs(1) = 2 * s_elec
elecs(2) = 2 * d_elec - 1
orbs(1) = 2 * s_orb - 1
orbs(2) = 2 * d_orb
else
pgen = 0.0_dp
end if
end if
case (excit_type%fullstart_stop_mixed)
! full start full stop mixed
ASSERT(elec_1 /= elec_2)
ASSERT(orb_1 /= orb_2)
ASSERT(csf_i%Occ_int(elec_1) == 1)
ASSERT(csf_i%Occ_int(elec_2) == 1)
if (csf_i%stepvector(elec_1) == csf_i%stepvector(elec_2)) then
pgen = 0.0_dp
else if (csf_i%stepvector(elec_1) == 1) then
elecs(1) = 2 * elec_1 - 1
elecs(2) = 2 * elec_2
orbs(1) = 2 * elec_2 - 1
orbs(2) = 2 * elec_1
else if (csf_i%stepvector(elec_1) == 2) then
elecs(1) = 2 * elec_1
elecs(2) = 2 * elec_2 - 1
orbs(1) = 2 * elec_2
orbs(2) = 2 * elec_1 - 1
end if
case (excit_type%fullstart_stop_alike)
! full-start full-stop alike
ASSERT(elec_1 == elec_2)
ASSERT(orb_1 == orb_2)
elecs(1) = 2 * elec_1 - 1
elecs(2) = 2 * elec_1
orbs(1) = 2 * orb_1 - 1
orbs(2) = 2 * orb_1
case default
! now the general 4-index double excitations..
! this can be nasty again..
call pick_random_4ind(csf_i, elec_1, elec_2, orb_1, orb_2, elecs, orbs, pgen)
end select
end if
end subroutine create_random_spin_orbs