init_three_body_const_mat Subroutine

    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