subroutine calcFullStartRaising(ilut, csf_i, excitInfo, excitations, nExcits, &
posSwitches, negSwitches)
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(inout) :: excitInfo
integer(n_int), intent(out), allocatable :: excitations(:, :)
integer, intent(out) :: nExcits
real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
character(*), parameter :: this_routine = "calcFullStartRaising"
integer :: nMax, ierr, iOrb, start, ende, semi, gen
integer(n_int) :: t(0:nifguga)
real(dp) :: tempWeight, minusWeight, plusWeight, bVal, nOpen
type(WeightObj_t) :: weights
integer(n_int), allocatable :: tempExcits(:, :)
ASSERT(isZero(ilut, excitInfo%fullStart))
ASSERT(isProperCSF_ilut(ilut))
! create the fullStart
start = excitInfo%fullStart
ende = excitInfo%fullEnd
semi = excitInfo%firstEnd
gen = excitInfo%lastGen
bVal = csf_i%B_real(semi)
! first have to allocate both arrays for the determinant and weight list
! worse then worst case for single excitations are 2^|i-j| excitations
! for a given CSF -> for now use that to allocate the arrays first
! update only need the number of open orbitals between i and j, and
! some additional room if 0/3 at start
! use already provided open orbital counting function.
! nMax = 2**(ende - start)
nMax = 2 + 2**count_open_orbs_ij(csf_i, start, ende)
allocate(tempExcits(0:nifguga, nMax), stat=ierr)
t = ilut
! have to make the switch and store the first matrix element in it too
! additionally also already calculate the sign coming from the
! pseudo double excitation which only depends on the number of open
! orbitals in the overlap region
nOpen = real(count_open_orbs_ij(csf_i, start, semi - 1), dp)
! set 0->3
set_orb(t, 2 * start)
set_orb(t, 2 * start - 1)
! could encode matrix element here, but do it later for more efficiency
! call encode_part_sign(t,Root2 * (-1.0)**nMax, 1)
! also dont need to store it here, also later
! tempExcits(:,1) = t
! since only the 0 branch is taken, where there are now changes in
! the stepvector and the matrix element is just a sign i can directly
! go to the lowering semi stop
! i could write it specifically here, and i should for efficiency
! have to calc. weights here, which are only the normel single
! excitation weights
weights = init_singleWeight(csf_i, ende)
minusWeight = weights%proc%minus(negSwitches(semi), csf_i%B_real(semi), weights%dat)
plusWeight = weights%proc%plus(posSwitches(semi), csf_i%B_real(semi), weights%dat)
ASSERT(.not. isZero(ilut, semi))
ASSERT(minusWeight + plusWeight > 0.0_dp)
call encode_matrix_element(t, 0.0_dp, 2)
select case (csf_i%stepvector(semi))
case (3)
! have to check a few things with 3 semi stop
! did something wrong with matrix elements before here.. the 3
! semi-stop has an aditional sign...but all matrix elements
! the same..
call encode_matrix_element(t, (-1.0_dp)**(nOpen + 1.0_dp), 1)
! bullshit...
if (minusWeight > 0.0_dp .and. plusWeight > 0.0_dp) then
! both +1 and -1 excitations possible
! do 3->1 first -1 branch first
clr_orb(t, 2 * semi)
call setDeltaB(-1, t)
tempExcits(:, 1) = t
! then do 3->2: +1 branch(change from already set 1
clr_orb(t, 2 * semi - 1)
set_orb(t, 2 * semi)
call setDeltaB(1, t)
tempExcits(:, 2) = t
nExcits = 2
else if (near_zero(plusWeight)) then
! only -1 branch possible
! do 3->1 first -1 branch first
clr_orb(t, 2 * semi)
call setDeltaB(-1, t)
tempExcits(:, 1) = t
nExcits = 1
else if (near_zero(minusWeight)) then
! only +1 branch possible
! then do 3->2: +1 branch
clr_orb(t, 2 * semi - 1)
call getDoubleMatrixElement(2, 3, 0, gen_type%R, gen_type%R, bVal, 1.0_dp, tempWeight)
call encode_matrix_element(t, Root2 * tempWeight * ((-1.0_dp)**nOpen), 1)
call setDeltaB(1, t)
tempExcits(:, 1) = t
nExcits = 1
else
! something went wrong, or i have to cancel exciation.. todo how
call stop_all(this_routine, "something went wrong. shouldnt be here!")
end if
case (1)
! only one excitation possible, which also has to have
! non-zero weight or otherwise i wouldnt even be here
! and reuse already provided t
! 1 -> 0
clr_orb(t, 2 * semi - 1)
call getDoubleMatrixElement(0, 1, 0, gen_type%R, gen_type%R, bVal, 1.0_dp, tempWeight)
call setDeltaB(1, t)
! encode fullstart contribution and pseudo overlap region here
! too in one go. one additional -1 due to semistop
! no Root2 here, since it cancels with t from matEle,
call encode_matrix_element(t, Root2 * tempWeight * ((-1.0_dp)**(nOpen)), 1)
tempExcits(:, 1) = t
nExcits = 1
case (2)
! here we have
! 2 -> 0
clr_orb(t, 2 * semi)
call getDoubleMatrixElement(0, 2, 0, gen_type%R, gen_type%R, bVal, 1.0_dp, tempWeight)
call setDeltaB(-1, t)
! encode fullstart contribution and pseudo overlap region here
! too in one go. one additional -1 due to semistop
! no Root2 here, since it cancels with t from matEle,
call encode_matrix_element(t, Root2 * tempWeight * ((-1.0_dp)**(nOpen)), 1)
tempExcits(:, 1) = t
nExcits = 1
end select
! and then we have to do just a regular single excitation
excitInfo%currentGen = excitInfo%lastGen
do iOrb = semi + 1, ende - 1
call singleUpdate(ilut, csf_i, iOrb, excitInfo, posSwitches, negSwitches, &
weights, tempExcits, nExcits)
end do
! and do a regular end step
call singleEnd(ilut, csf_i, excitInfo, tempExcits, nExcits, excitations)
! that should be it...
end subroutine calcFullStartRaising