We present a renormalon-subtracted formulation of potential nonrelativistic QCD (pNRQCD) for precision spectroscopy of heavy hadron systems, combining variational and Green’s function Monte Carlo (VMC/GFMC) methods with static two- and three-body potentials through NNLO in ${\alpha }_s$ (leading order in $1/m_Q$). Minimal renormalon subtraction (MRS) systematically sums leading factorially growing terms, thereby stabilizing perturbative convergence and reducing renormalization-scale dependence. We tune charm and bottom quark masses to spin-averaged $1S$ quarkonium states and predict ${\Omega }_{ccc}$, ${\Omega }_{ccb}$, ${\Omega }_{cbb}$, and ${\Omega }_{bbb}$ baryon masses, as well as QCD-stable baryons containing top quarks. NNLO MRS results undershoot lattice QCD by 125–175~MeV, a difference we attribute to contributions beyond the static perturbative potentials, including $1/m_Q$ corrections. Applying these methods to unequal-mass fully-heavy tetraquarks, we determine the critical heavy-to-light mass ratio for binding and compute binding energies across the mass-ratio landscape.