subroutine calcRaisingSemiStart(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 = "calcRaisingSemiStart"
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)
fe = excitInfo%firstEnd
en = excitInfo%fullEnd
! now if there is a LR(3) end it also restricts these type of
! exitations to pseudo doubles -> make distinction
! to use it for general double excitation i also have to check if
! i am dealing with a full stop. otherwise stepvalue doesnt matter
! at end
! could avoid that is Three(end) now since included in weights...
if (csf_i%stepvector(en) == 3 .and. en == fe) 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 assured correctly by weights that only -1 branches
! arrive...
cnt = 0
do iEx = 1, nExcits
if (getDeltaB(tempExcits(:, iEx)) == -1) then
t = tempExcits(:, iEx)
! change 1 -> 3
set_orb(t, 2 * s)
call getDoubleMatrixElement(3, 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
cnt = 0
do iEx = 1, nExcits
if (getDeltaB(tempExcits(:, iEx)) == +1) then
t = tempExcits(:, iEx)
! change 2-> 3
set_orb(t, 2 * s - 1)
call getDoubleMatrixElement(3, 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 (0)
! has to be 0 in this generator combination.
ASSERT(isZero(ilut, s))
! switches to 0 branch always possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == -1) then
! change 0->2
set_orb(t, 2 * s)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, &
excitInfo%order, tempWeight)
else
! change 0->1
set_orb(t, 2 * s - 1)
call getDoubleMatrixElement(1, 0, 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
! since no choice here do not have to check...
! hope this implementation is correct ...
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
! change 1 -> 3
set_orb(t, 2 * s)
call setDeltaB(deltaB + 1, t)
call getDoubleMatrixElement(3, 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 -> 3
set_orb(t, 2 * s - 1)
call setDeltaB(deltaB - 1, t)
call getDoubleMatrixElement(3, 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 (0)
! has to be 0 in raising 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
do iEx = 1, nExcits
t = tempExcits(:, iEx)
tempWeight = extract_matrix_element(t, 1)
deltaB = getDeltaB(t)
! first do 0 -> 1
set_orb(t, 2 * s - 1)
call getDoubleMatrixElement(1, 0, 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(:, iEx) = t
! then do 0->2
nExcits = nExcits + 1
clr_orb(t, 2 * s - 1)
set_orb(t, 2 * s)
call getDoubleMatrixElement(2, 0, 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
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == -1) then
! do the 0 > 1 switch to the -2 branch
set_orb(t, 2 * s - 1)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, x1_element=tempWeight_1)
call setDeltaB(deltaB - 1, t)
else
! do the 0 -> 2 switch to +2 branch
set_orb(t, 2 * s)
call getDoubleMatrixElement(2, 0, 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
! +2 excitations not possible
! only branching for -1 branch
! but excitation 0->1 everytime possible for both
do iEx = 1, nExcits
t = tempExcits(:, iEx)
! make all possible 0->1 excitation first
deltaB = getDeltaB(t)
set_orb(t, 2 * s - 1)
! todo: see how deltaB access correctly works...
! have to use new deltaB here i Think!!!
call getDoubleMatrixElement(1, 0, 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 1 -> 2
set_orb(t, 2 * s)
clr_orb(t, 2 * s - 1)
call getDoubleMatrixElement(2, 0, 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 ensure no +1 branch arrives...
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 0 > 1 switch to -2 branch
set_orb(t, 2 * s - 1)
call getDoubleMatrixElement(1, 0, 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
! -2 excitations not possible
! only branching for +1 branch
! 0->2 excitation works for both
do iEx = 1, nExcits
t = tempExcits(:, iEx)
! make all possible 0->2 excitation first
deltaB = getDeltaB(t)
set_orb(t, 2 * s)
! todo: see how deltaB access correctly works...
! have to use new deltaB here i Think!!!
call getDoubleMatrixElement(2, 0, 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, 0, 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
! so check that no -1 branch arrives
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 0 -> 2 switch to +2 branch
set_orb(t, 2 * s)
call getDoubleMatrixElement(2, 0, 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 0->2
set_orb(t, 2 * s)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, tempWeight_0, &
tempWeight_1)
else
! +1 branch: 0->1
set_orb(t, 2 * s - 1)
call getDoubleMatrixElement(1, 0, 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 lowering semistart!")
end if
end select
end if
end subroutine calcRaisingSemiStart