  We provide a formulation of potential non-relativistic quantum chromodynamics
(pNRQCD) suitable for calculating binding energies and matrix elements of
generic hadron and multi-hadron states made of heavy quarks in $SU(N_c)$ gauge
theory using quantum Monte Carlo techniques. We compute masses of quarkonium
and triply-heavy baryons in order to study the perturbative convergence of
pNRQCD and validate our numerical methods. Further, we study $SU(N_c)$ models
of composite dark matter and provide simple power series fits to our pNRQCD
results that can be used to relate dark meson and baryon masses to the
fundamental parameters of these models. For many systems comprised entirely of
heavy quarks, the quantum Monte Carlo methods employed here are less
computationally demanding than lattice field theory methods, although they
introduce additional perturbative approximations. The formalism presented here
may therefore be particularly useful for predicting composite dark matter
properties for a wide range of $N_c$ and heavy fermion masses.
