real(dp) function get_j_opt(nI, corr_J)
! routine to evaluate Hongjuns J-optimization formulas
integer, intent(in) :: nI(nel)
real(dp), intent(in) :: corr_J
integer :: i, j, a, spin_p, spin_q, b
integer(n_int) :: ilut(0:niftot)
type(symmetry) :: p_sym, q_sym, a_sym, b_sym, k_sym
integer :: src(2), tgt(2)
real(dp) :: sgn
real(dp) :: two, rpa, exchange, sum_3, sum_hel
call EncodeBitDet(nI, ilut)
get_j_opt = 0.0_dp
two = 0.0_dp
rpa = 0.0_dp
exchange = 0.0_dp
sum_3 = 0.0_dp
sum_hel = 0.0_dp
if (.not. t_symmetric) then
do i = 1, nel! -1
do j = 1, nel
! i only have a contribution if the spins of nI and nJ
! are not the same!
if (.not. same_spin(nI(i), nI(j))) then
! and then I need to loop over the holes, but due to
! momentum conservation, only once!
do a = 1, nBasis
! if a is empty
if (IsNotOcc(ilut, a)) then
b = get_orb_from_kpoints(nI(i), nI(j), a)
if (IsNotOcc(ilut, b) .and. .not. same_spin(a, b)) then
p_sym = G1(nI(i))%sym
q_sym = G1(nI(j))%sym
spin_p = get_spin_pn(nI(i))
spin_q = get_spin_pn(nI(j))
! and now i have to think how to correly
! choose the momenta involved
if (same_spin(nI(i), a)) then
a_sym = G1(a)%sym
b_sym = G1(b)%sym
else
a_sym = G1(b)%sym
b_sym = G1(a)%sym
end if
k_sym = SymTable(p_sym%s, SymConjTab(a_sym%s))
! since i loop over all possible i,j i do not need
! the sum like below i think
get_j_opt = get_j_opt + &
two_body_contrib(corr_J, p_sym, a_sym) + &
three_body_rpa_contrib(corr_J, p_sym, a_sym, spin_p) + &
three_body_exchange_contrib(nI, corr_J, p_sym, q_sym, a_sym, spin_q)
two = two + two_body_contrib(corr_J, p_sym, a_sym)
rpa = rpa + three_body_rpa_contrib(corr_J, p_sym, a_sym, spin_p)
exchange = exchange + three_body_exchange_contrib(nI, corr_J, p_sym, q_sym, a_sym, spin_p)
end if
end if
end do
end if
end do
end do
else
do i = 1, nel! - 1
do j = 1, nel
! i only have a contribution if the spins of nI and nJ
! are not the same!
if (.not. same_spin(nI(i), nI(j))) then
! and then I need to loop over the holes, but due to
! momentum conservation, only once!
do a = 1, nBasis
! if a is empty
if (IsNotOcc(ilut, a)) then
b = get_orb_from_kpoints(nI(i), nI(j), a)
if (IsNotOcc(ilut, b) .and. .not. same_spin(a, b)) then! .and. a < b) then
src = [min(nI(i), nI(j)), max(nI(i), nI(j))]
tgt = [min(a, b), max(a, b)]
if (is_beta(src(1))) then
p_sym = G1(src(1))%sym
q_sym = G1(src(2))%sym
else
p_sym = G1(src(2))%sym
q_sym = G1(src(1))%sym
end if
spin_p = get_spin_pn(src(1))
spin_q = get_spin_pn(src(2))
if (same_spin(src(1), tgt(1))) then
a_sym = G1(tgt(1))%sym
b_sym = G1(tgt(2))%sym
sgn = 1.0_dp
else
a_sym = G1(tgt(2))%sym
b_sym = G1(tgt(1))%sym
sgn = -1.0_dp
end if
k_sym = SymTable(p_sym%s, SymConjTab(a_sym%s))
! since i loop over all possible i,j i do not need
! the sum like below i think
get_j_opt = get_j_opt + ( &
two_body_contrib(corr_J, p_sym, a_sym) + &
two_body_contrib(corr_J, q_sym, b_sym) + &
three_body_rpa_contrib(corr_J, p_sym, a_sym, spin_p) + &
three_body_rpa_contrib(corr_J, q_sym, b_sym, spin_q) + &
three_body_exchange_contrib(nI, corr_J, p_sym, q_sym, a_sym, spin_p) + &
three_body_exchange_contrib(nI, corr_J, q_sym, p_sym, b_sym, spin_q))
end if
end if
end do
end if
end do
end do
end if
end function get_j_opt