subroutine calcLoweringSemiStop(ilut, csf_i, excitInfo, tempExcits, nExcits, plusWeight, &
minusWeight, t_no_singles_opt)
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t) :: excitInfo
integer(n_int), intent(inout), allocatable :: tempExcits(:, :)
integer, intent(inout) :: nExcits
real(dp), intent(in) :: plusWeight, minusWeight
logical, intent(in), optional :: t_no_singles_opt
character(*), parameter :: this_routine = "calcLoweringSemiStop"
integer :: se, iEx, deltaB, st, ss
integer(n_int) :: t(0:nifguga), s(0:nifguga)
real(dp) :: tempWeight, tempWeight_0, tempWeight_1, bVal
logical :: t_no_singles
ASSERT(isProperCSF_ilut(ilut))
ASSERT(.not. isThree(ilut, excitInfo%firstEnd))
ASSERT(plusWeight >= 0.0_dp)
ASSERT(minusWeight >= 0.0_dp)
se = excitInfo%firstEnd
bVal = csf_i%B_real(se)
st = excitInfo%fullStart
ss = excitInfo%secondStart
if (present(t_no_singles_opt)) then
t_no_singles = t_no_singles_opt
else
t_no_singles = .false.
end if
! do it similar to semi-starts and check if its a pseudo excit
! TODO!! have to make additional if here, to check i comes from
! a full start, because i also want to use it for normal double
! excitations, and there it doesnt matter what the stepvector value
! at the fullstart is!!!
if (csf_i%stepvector(st) == 3 .and. ss == st) then
! only 0 branches in this case
! first do the non-branching possibs
select case (csf_i%stepvector(se))
case (1)
! 1 -> 3 switch
if (near_zero(plusWeight)) then
tempExcits = 0
nExcits = 0
return
end if
do iEx = 1, nExcits
t = tempExcits(:, iEx)
ASSERT(getDeltaB(t) == 0)
! also need to assert weight is non-zero as it should be
! for mixed full-starts or general double excitations
! i cant really guaruantee it comes here with
! a plusWeight > 0, since 0 branch is always a
! possibility -> so i have to remove excitations if the
! weight is 0
set_orb(t, 2 * se)
call getDoubleMatrixElement(3, 1, 0, excitInfo%gen1, &
excitInfo%gen2, bVal, 1.0_dp, tempWeight)
call update_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
call setDeltaB(1, t)
tempExcits(:, iEx) = t
end do
case (2)
! 2 -> 3 switch
if (near_zero(minusWeight)) then
tempExcits = 0
nExcits = 0
return
end if
do iEx = 1, nExcits
t = tempExcits(:, iEx)
ASSERT(getDeltaB(t) == 0)
set_orb(t, 2 * se - 1)
call getDoubleMatrixElement(3, 2, 0, excitInfo%gen1, &
excitInfo%gen2, bVal, 1.0_dp, tempWeight)
call update_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
call setDeltaB(-1, t)
tempExcits(:, iEx) = t
end do
case (0)
! zero case -> branching possibilities according to b an weights
if (csf_i%B_int(se) > 0 .and. plusWeight > 0.0_dp &
.and. minusWeight > 0.0_dp) then
! both excitations are possible then
do iEx = 1, nExcits
t = tempExcits(:, iEx)
ASSERT(getDeltaB(t) == 0)
! have to store t weight since i need it twice
tempWeight_1 = extract_matrix_element(t, 1)
! do 0->1 first
set_orb(t, 2 * se - 1)
call setDeltaB(-1, t)
call getDoubleMatrixElement(1, 0, 0, excitInfo%gen1, &
excitInfo%gen2, bVal, 1.0_dp, tempWeight)
call update_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
! then change to 0->2 branch
set_orb(t, 2 * se)
clr_orb(t, 2 * se - 1)
call setDeltaB(+1, t)
call getDoubleMatrixElement(2, 0, 0, excitInfo%gen1, &
excitInfo%gen2, bVal, 1.0_dp, tempWeight)
call encode_matrix_element(t, tempWeight * tempWeight_1, 1)
call encode_matrix_element(t, 0.0_dp, 2)
nExcits = nExcits + 1
tempExcits(:, nExcits) = t
end do
else if (csf_i%B_int(se) == 0 .or. near_zero(plusWeight)) then
! only -1 branch possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
ASSERT(getDeltaB(t) == 0)
ASSERT(minusWeight > 0.0_dp)
! 0 -> 1
set_orb(t, 2 * se - 1)
call setDeltaB(-1, t)
call getDoubleMatrixElement(1, 0, 0, excitInfo%gen1, &
excitInfo%gen2, bVal, 1.0_dp, tempWeight)
call update_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
end do
else if (near_zero(minusWeight) .and. csf_i%B_int(se) > 0) then
! only +1 branches possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
ASSERT(getDeltaB(t) == 0)
ASSERT(plusWeight > 0.0_dp)
! 0 -> 2
set_orb(t, 2 * se)
call setDeltaB(+1, t)
call getDoubleMatrixElement(2, 0, 0, excitInfo%gen1, &
excitInfo%gen2, bVal, 1.0_dp, tempWeight)
call update_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
end do
else if (near_zero(minusWeight) .and. csf_i%B_int(se) == 0) then
! in this case no excitaiton is possible due to b value todo
call stop_all(this_routine, "implement cancelled excitations")
else
! shouldnt be here...
call stop_all(this_routine, "somethin went wrong! shouldnt be here!")
end if
end select
else
! also +2, and -2 branches arriving possible!
! again the non-branching values first
select case (csf_i%stepvector(se))
case (1)
! 1 -> 3 for 0 and -2 branch
! have to check different weight combinations
! although excitations should only get there if the weight
! fits.... -> check that
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (t_no_singles) then
if (.not. near_zero(extract_matrix_element(t, 1))) then
ASSERT(DeltaB == 0)
call encode_matrix_element(t, 0.0_dp, 1)
call encode_matrix_element(t, 0.0_dp, 2)
end if
end if
! check if weights fit.
! do not need actually, since no choice in excitations...
! and if dealt correctly until now it should be fine.
! hopefully atleast
! do 1 -> 3
set_orb(t, 2 * se)
call setDeltaB(deltaB + 1, t)
call getDoubleMatrixElement(3, 1, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order1, tempWeight_0, tempWeight_1)
! after semi-stop i only need the sum of the matrix
! elements
tempWeight = extract_matrix_element(t, 1) * tempWeight_0 + &
extract_matrix_element(t, 2) * tempWeight_1
call encode_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
end do
case (2)
! 2 -> 3 for 0 and +2 branches
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (t_no_singles) then
if (.not. near_zero(extract_matrix_element(t, 1))) then
ASSERT(DeltaB == 0)
call encode_matrix_element(t, 0.0_dp, 1)
call encode_matrix_element(t, 0.0_dp, 2)
end if
end if
! do 2 -> 3
set_orb(t, 2 * se - 1)
call setDeltaB(deltaB - 1, t)
call getDoubleMatrixElement(3, 2, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order1, tempWeight_0, tempWeight_1)
tempWeight = extract_matrix_element(t, 1) * tempWeight_0 + &
extract_matrix_element(t, 2) * tempWeight_1
call encode_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
end do
case (0)
! 0 case -> have to check weights more thourougly
! when a certain +-2 excitation comes to a 0 semi-stop
! the ongoing weights should allow an excitation or
! otherwise the excitation shouldnt even have been created!
! for the 0 branch arriving i have to check if a branching
! is possible.. and have to do that outside of the do-loops
! to be more efficient
if (csf_i%B_int(se) > 0 .and. plusWeight > 0.0_dp &
.and. minusWeight > 0.0_dp) then
! all excitations for 0 branch possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (t_no_singles) then
if (.not. near_zero(extract_matrix_element(t, 1))) then
ASSERT(DeltaB == 0)
call encode_matrix_element(t, 0.0_dp, 1)
call encode_matrix_element(t, 0.0_dp, 2)
call encode_matrix_element(tempExcits(:, iex), 0.0_dp, 1)
call encode_matrix_element(tempExcits(:, iex), 0.0_dp, 2)
end if
end if
if (deltaB == 0) then
! do 0 -> 1 first
set_orb(t, 2 * se - 1)
call setDeltaB(-1, t)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
tempWeight = extract_matrix_element(t, 1) * tempWeight_0 + &
extract_matrix_element(t, 2) * tempWeight_1
call encode_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
! need fresh temp. ilut for matrix elements
s = tempExcits(:, iEx)
tempExcits(:, iEx) = t
! then do 0->2 also
set_orb(s, 2 * se)
call setDeltaB(+1, s)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
tempWeight = extract_matrix_element(s, 1) * tempWeight_0 + &
extract_matrix_element(s, 2) * tempWeight_1
call encode_matrix_element(s, tempWeight, 1)
call encode_matrix_element(s, 0.0_dp, 2)
nExcits = nExcits + 1
tempExcits(:, nExcits) = s
else if (deltaB == -2) then
! only 0 -> 2 branch
set_orb(t, 2 * se)
call setDeltaB(-1, t)
! not sure anymore how matrix elements get
! adressed. maybe with incoming/outgoing deltaB
! todo! check that!!!
! x0 matrix element is always 0 in this case
! only take out x1 element
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
x1_element=tempWeight_1)
tempWeight_1 = extract_matrix_element(t, 2) * tempWeight_1
call encode_matrix_element(t, tempWeight_1, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
else
! only 0 -> 1 branch
set_orb(t, 2 * se - 1)
call setDeltaB(deltaB - 1, t)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
x1_element=tempWeight_1)
tempWeight_1 = extract_matrix_element(t, 2) * tempWeight_1
call encode_matrix_element(t, tempWeight_1, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
end if
end do
else if (csf_i%B_int(se) == 0 .or. near_zero(plusWeight)) then
! only -1 branch when 0 branch arrives... the switch from
! +2 -> +1 branch shouldnt be affected, since i wouldn not
! arrive at semi.stop if 0 weight, and if b value would
! be so low it also wouldn be possible to have this kind
! of excitation at this point
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (t_no_singles) then
if (.not. near_zero(extract_matrix_element(t, 1))) then
ASSERT(DeltaB == 0)
call encode_matrix_element(t, 0.0_dp, 1)
call encode_matrix_element(t, 0.0_dp, 2)
end if
end if
ASSERT(deltaB /= 2)
if (deltaB == 0) then
! do 0 -> 1 switch
set_orb(t, 2 * se - 1)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
tempWeight = extract_matrix_element(t, 1) * tempWeight_0 + &
extract_matrix_element(t, 2) * tempWeight_1
else
! -2 branch
! do 0->2
set_orb(t, 2 * se)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
x1_element=tempWeight_1)
tempWeight = extract_matrix_element(t, 2) * tempWeight_1
end if
call setDeltaB(-1, t)
call encode_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
end do
else if (csf_i%B_int(se) > 0 .and. near_zero(minusWeight)) then
! only +1 branch possible afterwards
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (t_no_singles) then
if (.not. near_zero(extract_matrix_element(t, 1))) then
ASSERT(DeltaB == 0)
call encode_matrix_element(t, 0.0_dp, 1)
call encode_matrix_element(t, 0.0_dp, 2)
end if
end if
ASSERT(deltaB /= -2)
if (deltaB == 0) then
! do 0->2 +1 branch
set_orb(t, 2 * se)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
tempWeight = extract_matrix_element(t, 1) * tempWeight_0 + &
extract_matrix_element(t, 2) * tempWeight_1
else
! +2 branch
! do 0 -> 1
set_orb(t, 2 * se - 1)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
x1_element=tempWeight_1)
tempWeight = extract_matrix_element(t, 2) * tempWeight_1
end if
call setDeltaB(+1, t)
call encode_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
tempExcits(:, iEx) = t
end do
else if (csf_i%B_int(se) == 0 .and. near_zero(plusWeight)) then
! broken excitation due to b value restriction
! todo how to deal with that ...
call stop_all(this_routine, "broken excitation due to b value. todo!")
else
! shouldnt be here
call stop_all(this_routine, "something went wrong. shouldnt be here!")
end if
end select
end if
end subroutine calcLoweringSemiStop