subroutine createSingleStart(ilut, csf_i, excitInfo, posSwitches, negSwitches, &
weightObj, tempExcits, nExcits)
! subroutine to create full single excitations starts for ilut
! allocates the necessary arrays and fills up first excitations,
! depending on stepvalue in ilut, the corresponding b value, and
! probabilitstic weight functions(to exclude excitations eventually
! leading to zero weight)
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(in) :: excitInfo
real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
type(WeightObj_t), intent(in) :: weightObj
integer(n_int), intent(out), allocatable :: tempExcits(:, :)
integer, intent(out) :: nExcits
character(*), parameter :: this_routine = "createSingleStart"
integer :: ierr, nmax, st, gen
integer(n_int) :: t(0:nifguga)
real(dp) :: minusWeight, plusWeight, tempWeight, bVal
! have already asserted that start and end values of stepvector and
! generator type are consistent to allow for an excitation.
! maybe still assert here in debug mode atleast
#ifdef DEBUG_
! also assert we are not calling it for a weight gen. accidently
ASSERT(excitInfo%currentGen /= 0)
ASSERT(isProperCSF_ilut(ilut))
! also check if calculated b vector really fits to ilut
if (excitInfo%currentGen == gen_type%R) then
ASSERT(.not. isThree(ilut, excitInfo%fullStart))
else if (excitInfo%currentGen == gen_type%L) then
ASSERT(.not. isZero(ilut, excitInfo%fullStart))
end if
! also have checked if atleast on branch way can lead to an excitaiton
#endif
! for more oversight
st = excitInfo%fullStart
gen = excitInfo%currentGen
bVal = csf_i%B_real(st)
! 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 = 4 + 4 * 2**count_open_orbs_ij(csf_i, st, excitInfo%fullEnd)
allocate(tempExcits(0:nifguga, nMax), stat=ierr)
! create start depending on stepvalue of ilut at start, b value,
! which is already calculated at the start of the excitation call and
! probabilistic weights
! for now just scratch it up in an inefficient way.. optimize later
! copy ilut into first row(or column?) of tempExcits
tempExcits(:, 1) = ilut
! i think determinants get initiated with zero matrix element, but not
! sure ... anyway set matrix element to 1
! maybe need temporary ilut storage
t = ilut
nExcits = 0
select case (csf_i%stepvector(st))
case (1)
! set corresponding orbital to 0 or 3 depending on generator type
if (gen == gen_type%R) then ! raising gen case
! set the alpha orbital also to 1 to make d=1 -> d'=3
set_orb(t, 2 * st)
! does it work like that:
! would have all necessary values and if statements to directly
! calculated matrix element here... maybe do it here after all
! to be more efficient...
tempWeight = getSingleMatrixElement(3, 1, +1, gen, bVal)
else ! lowering gen case
! clear beta orbital to make d=1 -> d'=0
clr_orb(t, 2 * st - 1)
tempWeight = getSingleMatrixElement(0, 1, +1, gen, bVal)
end if
! save number of already stored excitations
nExcits = nExcits + 1
! store matrix element
call encode_matrix_element(t, tempWeight, 1)
! set deltaB:
! for both cases deltaB = +1, set that as signed weight
! update: use previously unused flags to encode deltaB
call setDeltaB(1, t)
! and store it in list:
tempExcits(:, nExcits) = t
case (2)
if (gen == gen_type%R) then
set_orb(t, 2 * st - 1)
! matrix elements
tempWeight = getSingleMatrixElement(3, 2, -1, gen, bVal)
else
clr_orb(t, 2 * st)
! matrix elements
tempWeight = getSingleMatrixElement(0, 2, -1, gen, bVal)
end if
nExcits = nExcits + 1
call encode_matrix_element(t, tempWeight, 1)
call setDeltaB(-1, t)
tempExcits(:, nExcits) = t
case (0, 3)
! if the deltaB = -1 weights is > 0 this one always works
! write weight calculation function! is a overkill here,
! but makes it easier to use afterwards with stochastic excitaions
! b value restriction now implemented ind the weights...
! -> plus weight is 0 if bValue == 0
! use new weight objects here.
minusWeight = weightObj%proc%minus(negSwitches(st), bVal, weightObj%dat)
plusWeight = weightObj%proc%plus(posSwitches(st), bVal, weightObj%dat)
! still have to do some abort excitation routine if both weights
! are 0
ASSERT(.not. near_zero(minusWeight + plusWeight))
if (minusWeight > 0.0_dp) then
if (gen == gen_type%R) then
! set beta bit
set_orb(t, 2 * st - 1)
! matrix element
tempWeight = getSingleMatrixElement(1, 0, -1, gen, bVal)
else
! clear alpha bit
clr_orb(t, 2 * st)
! matrix element
tempWeight = getSingleMatrixElement(1, 3, -1, gen, bVal)
end if
nExcits = nExcits + 1
call encode_matrix_element(t, tempWeight, 1)
call setDeltaB(-1, t)
tempExcits(:, nExcits) = t
end if
! if plusWeight and the bValue allows it also a deltaB=+1 branch
! dont need bVal check anymore, since already implemented in
! weight calc.
if (plusWeight > 0.0_dp) then
! reset t
t = ilut
if (gen == gen_type%R) then
! set alpha bit
set_orb(t, 2 * st)
! matrix element
tempWeight = getSingleMatrixElement(2, 0, +1, gen, bVal)
else
! cleat beta bit
clr_orb(t, 2 * st - 1)
! matrix element
tempWeight = getSingleMatrixElement(2, 3, +1, gen, bVal)
end if
! update list and number
nExcits = nExcits + 1
call encode_matrix_element(t, tempWeight, 1)
call setDeltaB(+1, t)
tempExcits(:, nExcits) = t
end if
end select
end subroutine createSingleStart