subroutine calcLoweringSemiStart(ilut, csf_i, excitInfo, tempExcits, nExcits, &
plusWeight, minusWeight, zeroWeight)
! try at creating a reusable raising generator semi start
! just realize, for double excitations i have to save 2 matrix elements
! so hopefully a way is to use the imaginary matrix element storage
! for these kind of excitations! -> ask simon
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(in) :: excitInfo
integer(n_int), intent(inout) :: tempExcits(:, :)
integer, intent(inout) :: nExcits
real(dp), intent(in) :: plusWeight, minusWeight, zeroWeight
character(*), parameter :: this_routine = "calcLoweringSemiStart"
integer :: iEx, deltaB, s, cnt, en, fe
real(dp) :: tempWeight_0, tempWeight_1, tempWeight, bVal
integer(n_int) :: t(0:nifguga)
! at semi start with alike generator full stops only the deltaB=0
! are compatible, i hope this is correctly accounted by the weights..
s = excitInfo%secondStart
bVal = csf_i%B_real(s)
en = excitInfo%fullEnd
fe = excitInfo%firstEnd
! now if there is a LR(3) end it also restricts these type of
! exitations to pseudo doubles -> make distinction
! the second if is to only apply this to fullstop cases but also be
! able to use it fore general double excitations
if (csf_i%stepvector(en) == 3 .and. fe == en) then
! here only swiches to 0 branch are possible
select case (csf_i%stepvector(s))
case (1)
! only -1 branches can lead to 0 branch
! not yet asssured by weights that only -1 branch is
! arriving -> exclude wrong excitations
cnt = 0
do iEx = 1, nExcits
if (getDeltaB(tempExcits(:, iEx)) == -1) then
t = tempExcits(:, iEx)
! change 1 -> 0
clr_orb(t, 2 * s - 1)
call getDoubleMatrixElement(0, 1, -1, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight)
call update_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
call setDeltaB(0, t)
cnt = cnt + 1
tempExcits(:, cnt) = t
end if
end do
nExcits = cnt
case (2)
! only +1 can lead to 0 branch
! not yet correctly accounted in weights, so an arriving +1
! branch is ensured... todo
cnt = 0
do iEx = 1, nExcits
if (getDeltaB(tempExcits(:, iEx)) == +1) then
t = tempExcits(:, iEx)
! change 2-> 0
clr_orb(t, 2 * s)
call getDoubleMatrixElement(0, 2, 1, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight)
call update_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
call setDeltaB(0, t)
cnt = cnt + 1
tempExcits(:, cnt) = t
end if
end do
nExcits = cnt
case (3)
! has to be 3 in this generator combination.
ASSERT(isThree(ilut, s))
! switches to 0 branch always possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == -1) then
! change 3 -> 2
clr_orb(t, 2 * s - 1)
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight)
else
! change 3->1
clr_orb(t, 2 * s)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight)
end if
call update_matrix_element(t, tempWeight, 1)
call encode_matrix_element(t, 0.0_dp, 2)
call setDeltaB(0, t)
tempExcits(:, iEx) = t
end do
end select
else
select case (csf_i%stepvector(s))
case (1)
! always works if weights are fitting, check possibs.
! think i can do it generally for both... since its the same
! operations only... only need to check if atleast on of the
! relevant weights is nonzero
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
! change 1 -> 0
clr_orb(t, 2 * s - 1)
call setDeltaB(deltaB + 1, t)
call getDoubleMatrixElement(0, 1, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight_0, tempWeight_1)
call encode_matrix_element(t, tempWeight_1 * &
extract_matrix_element(t, 1), 2)
call update_matrix_element(t, tempWeight_0, 1)
tempExcits(:, iEx) = t
end do
case (2)
! should also work generally, since if a weight is zero the
! non compatible branch shouldnt even arrive here
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
! change 2 -> 0
clr_orb(t, 2 * s)
call setDeltaB(deltaB - 1, t)
call getDoubleMatrixElement(0, 2, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight_0, tempWeight_1)
call encode_matrix_element(t, tempWeight_1 * &
extract_matrix_element(t, 1), 2)
call update_matrix_element(t, tempWeight_0, 1)
tempExcits(:, iEx) = t
end do
case (3)
! has to be 3 in lowering case
! update: for a fulldouble excitation the 0-weight
! is not always > 0 dependending on the semistop and
! fullend values!
! so do it new!
if ((.not. near_zero(minusWeight)) .and. (.not. near_zero(plusWeight)) &
.and. (.not. near_zero(zeroWeight))) then
! all branches possible
! how to most efficiently do all that...
do iEx = 1, nExcits
t = tempExcits(:, iEx)
tempWeight = extract_matrix_element(t, 1)
deltaB = getDeltaB(t)
! first do 3 -> 1
clr_orb(t, 2 * s)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
tempWeight_0, tempWeight_1)
call update_matrix_element(t, tempWeight_0, 1)
call encode_matrix_element(t, tempWeight * tempWeight_1, 2)
call setDeltaB(deltaB - 1, t)
tempExcits(:, iEx) = t
! then do 3->2
nExcits = nExcits + 1
clr_orb(t, 2 * s - 1)
set_orb(t, 2 * s)
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
tempWeight_0, tempWeight_1)
call encode_matrix_element(t, tempWeight_0 * tempWeight, 1)
call encode_matrix_element(t, tempWeight_1 * tempWeight, 2)
call setDeltaB(deltaB + 1, t)
tempExcits(:, nExcits) = t
end do
else if ((.not. near_zero(minusWeight)) .and. (.not. near_zero(plusWeight)) &
.and. near_zero(zeroWeight)) then
! cont. 0 branches not possible
! here a arriving -1 branch can only become a -2 branch
! and a +1 branch a +2 branch
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == -1) then
! then there is a 3 -> 1 switch
clr_orb(t, 2 * s)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, x1_element=tempWeight_1)
call setDeltaB(deltaB - 1, t)
else
! there is a 3 -> 2 switch
clr_orb(t, 2 * s - 1)
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, x1_element=tempWeight_1)
call setDeltaB(deltaB + 1, t)
end if
call encode_matrix_element(t, extract_matrix_element(t, 1) * &
tempWeight_1, 2)
call encode_matrix_element(t, 0.0_dp, 1)
tempExcits(:, iEx) = t
end do
else if ((.not. near_zero(minusWeight)) .and. near_zero(plusWeight) &
.and. (.not. near_zero(zeroWeight))) then
! cont. + branches not possible
! so for an arriving -1 branch both -2 and 0 branch are
! possible but for a +1 branch only the 0 branch
do iEx = 1, nExcits
t = tempExcits(:, iEx)
! make all possible 3->1 excitation first
deltaB = getDeltaB(t)
clr_orb(t, 2 * s)
! todo: see how deltaB access correctly works...
! have to use new deltaB here i Think!!!
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight_0, tempWeight_1)
! need unchanged matrix element for branching
tempWeight = extract_matrix_element(t, 1)
call encode_matrix_element(t, tempWeight * &
tempWeight_1, 2)
call update_matrix_element(t, tempWeight_0, 1)
call setDeltaB(deltaB - 1, t)
tempExcits(:, iEx) = t
! and for +1 branch there is also a switch
if (deltaB == -1) then
! have to switch from previuosly switched 1 -> 2
set_orb(t, 2 * s)
clr_orb(t, 2 * s - 1)
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight_0, tempWeight_1)
call encode_matrix_element(t, tempWeight_0 * &
tempWeight, 1)
call encode_matrix_element(t, tempWeight_1 * &
tempWeight, 2)
call setDeltaB(0, t)
nExcits = nExcits + 1
tempExcits(:, nExcits) = t
end if
end do
else if ((.not. near_zero(minusWeight)) .and. near_zero(plusWeight) &
.and. near_zero(zeroWeight)) then
! only - branch possible
! so only a -1 arriving branch can cont. to a -2 branch
! hopefully the weights coming before handle that correctly
! so no +1 branch arrives here -> check!
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == +1) then
call print_indices(excitInfo)
call write_guga_list(stdout, tempExcits(:, 1:nExcits))
ASSERT(deltaB /= 1)
end if
! do the 3 -> 1 switch to -2 branch
clr_orb(t, 2 * s)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bval, excitInfo%order, x1_element=tempWeight_1)
call setDeltaB(deltaB - 1, t)
call encode_matrix_element(t, extract_matrix_element(t, 1) * &
tempWeight_1, 2)
call encode_matrix_element(t, 0.0_dp, 1)
tempExcits(:, iEx) = t
end do
else if (near_zero(minusWeight) .and. (.not. near_zero(plusWeight)) &
.and. (.not. near_zero(zeroWeight))) then
! cont. - branch not possible
! so both options are possible for an incoming +1 branch
! and only the 0 branch for the -1
do iEx = 1, nExcits
t = tempExcits(:, iEx)
! make all possible 3->2 excitation first
deltaB = getDeltaB(t)
clr_orb(t, 2 * s - 1)
! todo: see how deltaB access correctly works...
! have to use new deltaB here i Think!!!
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight_0, tempWeight_1)
! need unchanged matrix element for branching
tempWeight = extract_matrix_element(t, 1)
call encode_matrix_element(t, tempWeight_0 * &
tempWeight, 1)
call encode_matrix_element(t, tempWeight_1 * &
tempWeight, 2)
call setDeltaB(deltaB + 1, t)
tempExcits(:, iEx) = t
! and for +1 branch there is also a switch
if (deltaB == +1) then
! have to switch from previuosly switched 2 -> 1
clr_orb(t, 2 * s)
set_orb(t, 2 * s - 1)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight_0, tempWeight_1)
call encode_matrix_element(t, tempWeight_0 * &
tempWeight, 1)
call encode_matrix_element(t, tempWeight_1 * &
tempWeight, 2)
call setDeltaB(0, t)
nExcits = nExcits + 1
tempExcits(:, nExcits) = t
end if
end do
else if (near_zero(minusWeight) .and. (.not. near_zero(plusWeight)) &
.and. near_zero(zeroWeight)) then
! only + possible
! only an arriving +1 branch can go on.. so check and
! hope there is no -1 branch arriving
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == -1) then
call print_indices(excitInfo)
call write_guga_list(stdout, tempExcits(:, 1:nExcits))
ASSERT(deltaB /= -1)
end if
! do the 3 -> 2 switch to the +2 branch
clr_orb(t, 2 * s - 1)
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, x1_element=tempWeight_1)
call setDeltaB(deltaB + 1, t)
call encode_matrix_element(t, extract_matrix_element(t, 1) * &
tempWeight_1, 2)
call encode_matrix_element(t, 0.0_dp, 1)
tempExcits(:, iEx) = t
end do
else if (near_zero(minusWeight) .and. near_zero(plusWeight) &
.and. (.not. near_zero(zeroWeight))) then
! only 0 branch possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == -1) then
! -1 branch 3->2
clr_orb(t, 2 * s - 1)
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, tempWeight_0, &
tempWeight_1)
else
! +1 branch: 3->1
clr_orb(t, 2 * s)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, tempWeight_0, &
tempWeight_1)
end if
call setDeltaB(0, t)
call encode_matrix_element(t, tempWeight_1 * &
extract_matrix_element(t, 1), 2)
call update_matrix_element(t, tempWeight_0, 1)
tempExcits(:, iEx) = t
end do
else
! something went wrong -> no branches possible
call stop_all(this_routine, " no branches possible in raising semistart!")
end if
end select
end if
end subroutine calcLoweringSemiStart