subroutine init_three_body_const_mat()
integer :: i, j
type(symmetry) :: sym_i, sym_j
if (allocated(three_body_const_mat)) deallocate(three_body_const_mat)
allocate(three_body_const_mat(nBasis / 2, nBasis / 2, -1:1), source=0.0_dp)
do i = 1, nBasis / 2
sym_i = G1(2 * i)%Sym
do j = 1, nBasis / 2
sym_j = G1(2 * j)%Sym
three_body_const_mat(sym_i%s, sym_j%s, -1) = &
-three_body_prefac &
* n_opp(-1) &
* real(epsilon_kvec(sym_i) + epsilon_kvec(sym_j), dp)
three_body_const_mat(sym_i%s, sym_j%s, 1) = &
-three_body_prefac &
* n_opp(1) &
* real(epsilon_kvec(sym_i) + epsilon_kvec(sym_j), dp)
end do
end do
end subroutine init_three_body_const_mat