Let \((M,H,g_H;g)\) be a sub-Riemannian manifold and (N, h) be a Riemannian manifold. For a smooth map \(u: M \rightarrow N\) , we consider the energy functional \(E_G(u) = \frac{1}{2} \int _M[|\textrm{d}u_\text {H}|^2 - 2\,G(u)] \textrm{d}V_M\) , where \(\textrm{d}u_\text {H}\) is the horizontal differential of u, \(G:N\rightarrow \mathbb {R}\) is a smooth function on N. The critical maps of \(E_G(u)\) are referred to as subelliptic harmonic maps with potential G. In this paper, we investigate the existence problem for subelliptic harmonic maps with potentials by a subelliptic heat flow. Assuming that the target Riemannian manifold has nonpositive sectional curvature and the potential G satisfies various suitable conditions, we prove some Eells–Sampson-type existence results when the source manifold is either a step-2 sub-Riemannian manifold or a step-r sub-Riemannian manifold whose sub-Riemannian structure comes from a tense Riemannian foliation.