pure function calcDiagExchangeGUGA_nI(iOrb, jOrb, nI) result(exchange)
! calculates the exchange contribution to the diagonal matrix elements
! this is the implementation if only nI is provided
integer, intent(in) :: iOrb, jOrb, nI(nEl)
real(dp) :: exchange
real(dp) :: bVector(nEl)
integer :: i
! the b-vector is also needed for these calculations:
bVector = calcB_vector_nI(nI)
! probably could use current b vector.. or reference b vector even...
! yes definitly. no not really since this is also used for general
! diagonal matrix elements not only the current determinant in the
! fciqmc loop
exchange = 1.0_dp
! then i need the exchange product term between orbitals s and p
do i = iOrb + 1, jOrb - 1
if (.not. isDouble(nI, i)) then
if (is_beta(nI(i))) then
exchange = exchange * functionA(bVector(i), 2.0_dp, 0.0_dp) &
* functionA(bVector(i), -1.0_dp, 1.0_dp)
else
exchange = exchange * functionA(bVector(i), 0.0_dp, 2.0_dp) * &
functionA(bVector(i), 3.0_dp, 1.0_dp)
end if
end if
end do
if (is_beta(nI(iOrb))) then
exchange = exchange * functionA(bVector(iOrb), 2.0_dp, 0.0_dp)
if (is_beta(nI(jOrb))) then
exchange = exchange * functionA(bVector(jOrb), -1.0_dp, 1.0_dp)
else
exchange = -exchange * functionA(bVector(jOrb), 3.0_dp, 1.0_dp)
end if
else
exchange = exchange * functionA(bVector(iOrb), 0.0_dp, 2.0_dp)
if (is_beta(nI(jOrb))) then
exchange = -exchange * functionA(bVector(jOrb), -1.0_dp, 1.0_dp)
else
exchange = exchange * functionA(bVector(jOrb), 3.0_dp, 1.0_dp)
end if
end if
end function calcDiagExchangeGUGA_nI