subroutine doubleUpdate(ilut, csf_i, sO, excitInfo, weights, tempExcits, nExcits, &
negSwitches, posSwitches)
! for double excitation weights are usually always already calculated
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
integer, intent(in) :: sO
type(ExcitationInformation_t), intent(in) :: excitInfo
type(WeightObj_t), intent(in) :: weights
integer(n_int), intent(inout) :: tempExcits(:, :)
integer, intent(inout) :: nExcits
real(dp), intent(in) :: negSwitches(nSpatOrbs), posSwitches(nSpatOrbs)
character(*), parameter :: this_routine = "doubleUpdate"
real(dp) :: minusWeight, plusWeight, zeroWeight
integer :: gen1, gen2, iEx, deltaB
real(dp) :: bVal, tempWeight_0, tempWeight_1, order
integer(n_int) :: t(0:nifguga), s(0:nifguga)
ASSERT(isProperCSF_ilut(ilut))
ASSERT(sO > 0 .and. sO <= nSpatOrbs)
if (csf_i%Occ_int(sO) /= 1) then
! in this case no change in stepvector or matrix element
return
end if
! and for more readibility extract certain values:
gen1 = excitInfo%gen1
gen2 = excitInfo%gen2
bVal = csf_i%B_real(sO)
! stupid!! order only needed at semistarts and semistops! so just set
! it to 1.0 here
order = 1.0_dp
! else extract the relevant weights:
minusWeight = weights%proc%minus(negSwitches(sO), bVal, weights%dat)
plusWeight = weights%proc%plus(posSwitches(sO), bVal, weights%dat)
zeroWeight = weights%proc%zero(negSwitches(sO), posSwitches(sO), bVal, &
weights%dat)
! try new implementation with checking all 3 weights...
! ===================================================================
if (csf_i%stepvector(sO) == 1) then
if (minusWeight > 0.0_dp .and. plusWeight > 0.0_dp .and. zeroWeight > 0.0_dp) then
! all excitations are possible in this case.
! first do the on track
do iEx = 1, nExcits
t = tempExcits(:, iEx)
! need it two times for matrix elements
s = t
deltaB = getDeltaB(t)
! just calc matrix element for 1 -> 1 on track
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
! if deltaB != 2 do :
if (deltaB /= 2) then
! then do the 1 -> 2 switch
clr_orb(s, 2 * sO - 1)
set_orb(s, 2 * sO)
! change deltaB
call setDeltaB(deltaB + 2, s)
! calc matrix elements, only x1 matrix elements here
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call encode_matrix_element(s, 0.0_dp, 1)
call update_matrix_element(s, tempWeight_1, 2)
nExcits = nExcits + 1
tempExcits(:, nExcits) = s
end if
end do
else if (near_zero(plusWeight) .and. (.not. near_zero(minusWeight)) &
.and. (.not. near_zero(zeroWeight))) then
! all purely + branches are not possible
! do the on track possibles first
do iEx = 1, nExcits
t = tempExcits(:, iEx)
s = t
deltaB = getDeltaB(t)
! i think i never should have a +2 here then now...
ASSERT(deltaB /= 2)
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
! and do the 1->2 to 0 branch in case of dB = -2
if (deltaB == -2) then
! 1 -> 2 switch
clr_orb(s, 2 * sO - 1)
set_orb(s, 2 * sO)
! change deltaB
call setDeltaB(deltaB + 2, s)
! calc matrix elements, only x1 matrix elements here
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call encode_matrix_element(s, 0.0_dp, 1)
call update_matrix_element(s, tempWeight_1, 2)
nExcits = nExcits + 1
tempExcits(:, nExcits) = s
end if
end do
else if (near_zero(minusWeight) .and. (.not. near_zero(plusWeight)) &
.and. (.not. near_zero(zeroWeight))) then
! -2 branch not possible.. not sure if i ever come here with
! dB = -2 then.. -> check
do iEx = 1, nExcits
t = tempExcits(:, iEx)
s = t
deltaB = getDeltaB(t)
! do the on-track specifically if a -2 comes
if (deltaB == -2) then
clr_orb(t, 2 * sO - 1)
set_orb(t, 2 * sO)
call setDeltaB(0, t)
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call encode_matrix_element(t, 0.0_dp, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
else
! elsewise do the staying possib
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
! and if dB = 0, do the possible +2 switch too
if (deltaB == 0) then
clr_orb(s, 2 * sO - 1)
set_orb(s, 2 * sO)
call setDeltaB(2, s)
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call encode_matrix_element(s, 0.0_dp, 1)
call update_matrix_element(s, tempWeight_1, 2)
nExcits = nExcits + 1
tempExcits(:, nExcits) = s
end if
end if
end do
else if (near_zero(minusWeight) .and. near_zero(plusWeight) &
.and. (.not. near_zero(zeroWeight))) then
! only the 0 branches are possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
ASSERT(deltaB /= +2)
if (deltaB == -2) then
! do the 1->2 0 swicth
clr_orb(t, 2 * sO - 1)
set_orb(t, 2 * sO)
call setDeltaB(0, t)
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
else
! if 0 comes just get the matrix elements
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
end if
! and update matrix elements
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
end do
else if ((.not. near_zero(plusWeight)) .and. (.not. near_zero(minusWeight)) &
.and. near_zero(zeroWeight)) then
! the continuing 0-branches are not allowed
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == 0) then
! do the switch i a 0 branch comes
clr_orb(t, 2 * sO - 1)
set_orb(t, 2 * sO)
call setDeltaB(2, t)
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
else
! otherwise just get the matrix elemet
! x0-always zero in this case
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
end if
call encode_matrix_element(t, 0.0_dp, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
end do
else if (near_zero(zeroWeight) .and. near_zero(plusWeight) &
.and. (.not. near_zero(minusWeight))) then
! only -2 staying branch is possible... so i should just be here
! with a -2
do iEx = 1, nExcits
t = tempExcits(:, iEx)
ASSERT(getDeltaB(t) == -2)
call getDoubleMatrixElement(1, 1, -2, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
! x0 element should already be 0
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
end do
else if (near_zero(zeroWeight) .and. near_zero(zeroWeight) &
.and. (.not. near_zero(plusWeight))) then
! only the +2 cont. branches are possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
ASSERT(deltaB /= -2)
if (deltaB == 0) then
! do the switch to +2
clr_orb(t, 2 * sO - 1)
set_orb(t, 2 * sO)
call setDeltaB(2, t)
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
else
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
end if
call encode_matrix_element(t, 0.0_dp, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
end do
else
! should not be here... all weights 0#
call stop_all(this_routine, "should not be here! all 3 weights 0 in double update at sO=1")
end if
! also do the possibilities at a step=2
else if (csf_i%stepvector(sO) == 2) then
if (minusWeight > 0.0_dp .and. plusWeight > 0.0_dp .and. zeroWeight > 0.0_dp) then
! everything possible!
do iEx = 1, nExcits
t = tempExcits(:, iEx)
s = t
deltaB = getDeltaB(t)
! do on track first
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
! if /= -2 2 -> 1 branching also possible
if (deltaB /= -2) then
clr_orb(s, 2 * sO)
set_orb(s, 2 * sO - 1)
call setDeltaB(deltaB - 2, s)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call encode_matrix_element(s, 0.0_dp, 1)
call update_matrix_element(s, tempWeight_1, 2)
nExcits = nExcits + 1
tempExcits(:, nExcits) = s
end if
end do
else if (near_zero(plusWeight) .and. (.not. near_zero(minusWeight)) &
.and. (.not. near_zero(zeroWeight))) then
! +2 cont. not possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
if (deltaB == -2) then
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
else if (deltaB == 2) then
clr_orb(t, 2 * sO)
set_orb(t, 2 * sO - 1)
call setDeltaB(0, t)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
else
s = t
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
! then do switch
clr_orb(s, 2 * sO)
set_orb(s, 2 * sO - 1)
call setDeltaB(-2, s)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call encode_matrix_element(s, 0.0_dp, 1)
call update_matrix_element(s, tempWeight_1, 2)
nExcits = nExcits + 1
tempExcits(:, nExcits) = s
end if
end do
else if (near_zero(minusWeight) .and. (.not. near_zero(plusWeight)) &
.and. (.not. near_zero(zeroWeight))) then
! cont. -2 branches not possible
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
s = t
ASSERT(deltaB /= -2)
! do staying first
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
! then do possible switch
if (deltaB == 2) then
clr_orb(s, 2 * sO)
set_orb(s, 2 * sO - 1)
call setDeltaB(0, s)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call update_matrix_element(s, tempWeight_1, 2)
nExcits = nExcits + 1
tempExcits(:, nExcits) = s
end if
end do
else if (near_zero(minusWeight) .and. near_zero(plusWeight) &
.and. (.not. near_zero(zeroWeight))) then
! only cont. 0 branch valid
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
ASSERT(deltaB /= -2)
if (deltaB == 2) then
! do switch
clr_orb(t, 2 * sO)
set_orb(t, 2 * sO - 1)
call setDeltaB(0, t)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
else
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
end if
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
end do
else if ((.not. near_zero(plusWeight)) .and. (.not. near_zero(minusWeight)) &
.and. near_zero(zeroWeight)) then
! no 0 branch valid
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
! only switch if 0 branch arrives
if (deltaB == 0) then
! do swicht
clr_orb(t, 2 * sO)
set_orb(t, 2 * sO - 1)
call setDeltaB(-2, t)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call encode_matrix_element(t, 0.0_dp, 1)
else
! staying is possible
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
end if
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
end do
else if (near_zero(zeroWeight) .and. near_zero(plusWeight) &
.and. (.not. near_zero(minusWeight))) then
! only -2 branches valis
do iEx = 1, nExcits
t = tempExcits(:, iEx)
deltaB = getDeltaB(t)
ASSERT(deltaB /= 2)
if (deltaB == 0) then
clr_orb(t, 2 * sO)
set_orb(t, 2 * sO - 1)
call setDeltaB(-2, t)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call encode_matrix_element(t, 0.0_dp, 1)
else
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
end if
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
end do
else if (near_zero(zeroWeight) .and. near_zero(minusWeight) &
.and. (.not. near_zero(plusWeight))) then
! only +2 staying branch is possible -> assert that only that
! comes
do iEx = 1, nExcits
t = tempExcits(:, iEx)
ASSERT(getDeltaB(t) == 2)
call getDoubleMatrixElement(2, 2, 2, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
call update_matrix_element(t, tempWeight_1, 2)
tempExcits(:, iEx) = t
end do
else
! should not be here.. all 3 weights 0
call stop_all(this_routine, "should not be here. all 3 weights 0 at s=2 in!")
end if
end if
end subroutine doubleUpdate