init_two_body_trancorr_fac_matrix Subroutine

    subroutine init_two_body_trancorr_fac_matrix()
        integer :: i, j
        type(symmetry) :: sym_i, sym_j

        ! for more efficiency, precompute the two-body factor for all possible
        ! symmetry symbols
        if (allocated(two_body_transcorr_factor_matrix)) deallocate(two_body_transcorr_factor_matrix)

        allocate(two_body_transcorr_factor_matrix(nBasis / 2, nBasis / 2), source=0.0_dp)

        ! loop over spatial orbitals
        do i = 1, nBasis / 2
            sym_i = G1(2 * i)%sym
            do j = 1, nBasis / 2
                sym_j = G1(2 * j)%Sym
                two_body_transcorr_factor_matrix(sym_j%s, sym_i%s) = &
                    bhub / omega &
                    * ((exp(trans_corr_param_2body) - 1.0_dp) * epsilon_kvec(sym_i) &
                        + (exp(-trans_corr_param_2body) - 1.0_dp) * epsilon_kvec(sym_j))
            end do
        end do

    end subroutine init_two_body_trancorr_fac_matrix