subroutine fill_sings_2rdm_guga(spawn, ilutI, ilutJ, sign_i, sign_j, mat_ele, rdm_ind)
! this is the routine where I test the filling of 2-RDM based on
! single excitations to mimic the workflow in the stochastic
! RDM sampling
type(rdm_spawn_t), intent(inout) :: spawn
integer(n_int), intent(in) :: ilutI(0:GugaBits%len_tot), ilutJ(0:GugaBits%len_tot)
real(dp), intent(in) :: sign_i(:), sign_j(:), mat_ele
integer(int_rdm), intent(in) :: rdm_ind
integer :: i, a, st, en, step_i(nSpatOrbs), step_j(nSpatOrbs), &
delta_b(nSpatOrbs), gen, d_i, d_j, delta, &
b_i(nSpatOrbs), b_j(nSpatOrbs)
real(dp) :: occ_i(nSpatOrbs), occ_j(nSpatOrbs)
real(dp) :: botCont, topCont, tempWeight, prod, StartCont, EndCont
integer :: iO, jO, step
! here I have to fill W + R/L, single overlap RR/LL
! and full-start/stop RL with no change in the double overlap region
! i also have to correct the coupling coefficient here effifiently
! for this it would be best to have both CSFs I and J involved.
! and i have to figure out all correct index combinations here
! for all the entries the singles contribute to.
! ok this is the final routine i need to consider i guess..
! or hope.. I need the matrix elements and indices, to which a
! certain type of single excitations contributes to..
! I need to effectively recalculate the correct matrix element
! and store it in the correct RDM place.
! essentially I need to loop over the contracted index and
! correctly take the coupling coefficient into account.
call extract_1_rdm_ind(rdm_ind, i, a)
st = min(i, a)
en = max(i, a)
! do i have access to current_stepvector here? I think so..
! no I dont! or just to be save
! this is essentially a mimic of the function
! calc_integral_contribution_single in guga_excitations
! but wait a minute..
! in the stochastic excitation generation, I calculate the
! full Hamiltonian matrix element..
! I do not have access to the coupling coefficient for
! i, a anymore.. damn..
! and in general, I always calculate the full matrix element..
! i have to take this into account!
! or change the stochastic excitation generation, to also yield
! this information..
! or "just" recalculate everything..
! i could use the x1 element storage for this quantity..
! this would simplify things!
! for simplicity it is best to have these quantities:
call calc_csf_i(ilutI, step_i, b_i, occ_i)
call calc_csf_i(ilutJ, step_j, b_j, occ_j)
delta_b = b_i - b_j
! calculate the bottom contribution depending on the excited stepvalue
select case (step_i(st))
case (0)
! this implicates a raising st:
if (isOne(ilutJ, st)) then
call getDoubleMatrixElement(1, 0, 0, gen_type%L, gen_type%R, real(b_i(st), dp), &
1.0_dp, x1_element=botCont)
else
call getDoubleMatrixElement(2, 0, 0, gen_type%L, gen_type%R, real(b_i(st), dp), &
1.0_dp, x1_element=botCont)
end if
StartCont = 0.0_dp
gen = gen_type%R
case (3)
! implies lowering st
if (isOne(ilutJ, st)) then
! need tA(0,2)
botCont = funA_0_2overR2(real(b_i(st), dp))
else
! need -tA(2,0)
botCont = minFunA_2_0_overR2(real(b_i(st), dp))
end if
StartCont = 1.0_dp
gen = gen_type%L
case (1)
botCont = funA_m1_1_overR2(real(b_i(st), dp))
! check which generator
if (isZero(ilutJ, st)) then
botCont = -botCont
StartCont = 0.0_dp
gen = gen_type%L
else
StartCont = 1.0_dp
gen = gen_type%R
end if
case (2)
botCont = funA_3_1_overR2(real(b_i(st), dp))
if (isThree(ilutJ, st)) then
botCont = -botCont
StartCont = 1.0_dp
gen = gen_type%R
else
StartCont = 0.0_dp
gen = gen_type%L
end if
end select
! do top contribution also already
select case (step_i(en))
case (0)
if (isOne(ilutJ, en)) then
topCont = funA_2_0_overR2(real(b_i(en), dp))
else
topCont = minFunA_0_2_overR2(real(b_i(en), dp))
end if
EndCont = 0.0_dp
case (3)
if (isOne(ilutJ, en)) then
topCont = minFunA_2_0_overR2(real(b_i(en), dp))
else
topCont = funA_0_2overR2(real(b_i(en), dp))
end if
EndCont = 1.0_dp
case (1)
topCont = funA_2_0_overR2(real(b_i(en), dp))
if (isThree(ilutJ, en)) then
topCont = -topCont
EndCont = 1.0_dp
else
EndCont = 0.0_dp
end if
case (2)
topCont = funA_0_2overR2(real(b_i(en), dp))
if (isZero(ilutJ, en)) then
topCont = -topCont
EndCont = 0.0_dp
else
EndCont = 1.0_dp
end if
end select
! depending on i and j calulate the corresponding single and double
! integral weights and check if they are non-zero...
! gets quite involved... :( need to loop over all orbitals
! have to reset prod inside the loop each time!
do iO = 1, st - 1
! no contribution if not occupied.
if (step_i(iO) == 0) cycle
! else it gets a contrbution weighted with orbital occupation
! first easy part:
! i think also singles I should store symmetrically! as
! the conjugated element will of course never be stored!
! W + R/L contribution
call add_to_rdm_spawn_t(spawn, iO, iO, i, a, &
occ_i(iO) * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, i, a, iO, iO, &
occ_i(iO) * sign_i * sign_j * mat_ele, .true.)
! exchange contribution:
! wtf is this b_i == 0 check? that does not make any sense..
if (step_i(iO) == 3 .or. ((occ_i(iO) .isclose.1.0_dp) .and. b_i(iO) == 0)) then
! then it is easy:
! just figure out correct indices
call add_to_rdm_spawn_t(spawn, i, iO, iO, a, &
-occ_i(iO) / 2.0 * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, iO, a, i, iO, &
-occ_i(iO) / 2.0 * sign_i * sign_j * mat_ele, .true.)
else
! otherwise i have to recalc the x1 element
step = step_i(iO)
call getDoubleMatrixElement(step, step, -1, gen_type%L, gen_type%R, real(b_i(iO), dp), &
1.0_dp, x1_element=prod)
! and then do the remaining:
do jO = iO + 1, st - 1
! need the stepvalue entries to correctly access the mixed
! generator matrix elements
step = step_i(jO)
call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, &
real(b_i(jO), dp), 1.0_dp, x1_element=tempWeight)
prod = prod * tempWeight
end do
prod = prod * botCont
call add_to_rdm_spawn_t(spawn, i, iO, iO, a, &
(prod - occ_i(iO) / 2.0) * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, iO, a, i, iO, &
(prod - occ_i(iO) / 2.0) * sign_i * sign_j * mat_ele, .true.)
end if
end do
! start segment: only W + R/L
! but this depends on the type of excitation here.. or?
call add_to_rdm_spawn_t(spawn, st, st, i, a, &
StartCont * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, i, a, st, st, &
StartCont * sign_i * sign_j * mat_ele, .true.)
! loop over excitation range:
do iO = st + 1, en - 1
! W + R/L
call add_to_rdm_spawn_t(spawn, iO, iO, i, a, &
occ_i(iO) * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, i, a, iO, iO, &
occ_i(iO) * sign_i * sign_j * mat_ele, .true.)
! exchange type:
! oh damn I need the delta-B value here..
d_i = step_i(iO)
d_j = step_j(iO)
delta = delta_b(iO - 1)
prod = getDoubleContribution(d_j, d_i, delta, gen, real(b_i(iO), dp))
call add_to_rdm_spawn_t(spawn, i, iO, iO, a, &
prod * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, iO, a, i, iO, &
prod * sign_i * sign_j * mat_ele, .true.)
end do
! end contribution
call add_to_rdm_spawn_t(spawn, en, en, i, a, &
EndCont * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, i, a, en, en, &
EndCont * sign_i * sign_j * mat_ele, .true.)
! loop above:
do iO = en + 1, nSpatOrbs
if (step_i(iO) == 0) cycle
! W + R/L contribution
call add_to_rdm_spawn_t(spawn, iO, iO, i, a, &
occ_i(iO) * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, i, a, iO, iO, &
occ_i(iO) * sign_i * sign_j * mat_ele, .true.)
! exchange
if (step_i(iO) == 3 .or. (b_i(iO) == 1 .and. step_i(iO) == 1)) then
! only x0 contribution
! then it is easy:
! just figure out correct indices
call add_to_rdm_spawn_t(spawn, i, iO, iO, a, &
-occ_i(iO) / 2.0 * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, iO, a, i, iO, &
-occ_i(iO) / 2.0 * sign_i * sign_j * mat_ele, .true.)
else
prod = 1.0_dp
do jO = en + 1, iO - 1
step = step_i(jO)
call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, real(b_i(jO), dp), &
1.0_dp, x1_element=tempWeight)
prod = prod * tempWeight
end do
step = step_i(iO)
call getMixedFullStop(step, step, 0, real(b_i(iO), dp), &
x1_element=tempWeight)
prod = prod * tempWeight
prod = prod * topCont
call add_to_rdm_spawn_t(spawn, i, iO, iO, a, &
(prod - occ_i(iO) / 2.0) * sign_i * sign_j * mat_ele, .true.)
call add_to_rdm_spawn_t(spawn, iO, a, i, iO, &
(prod - occ_i(iO) / 2.0) * sign_i * sign_j * mat_ele, .true.)
end if
end do
end subroutine fill_sings_2rdm_guga