Let \(f: (M,g)\to (N,h)\) be a smooth map between two Riemannian manifolds and \(G:N\to\mathbb{R}\) be a given function. The authors study the following energy functional \(E_G(f)={1\over 2}\int[|df|^2- 2G(f)]dv_g\), and call \(f\) the harmonic map with potential \(G\) if \(f\) satisfies the Euler-Lagrange equation \(\tau(f)+\nabla G(f)=0\). Some variational properties and some existence results for the functional \(E_G(f)\) are established. It is proved that if \(G\in C^0(S^2)\) is a nonconstant function, then \(E_G(f)\) does not attain its minimum in any non-trivial homotopy class of maps \(f:S^2\to S^2\). There exists a smooth map \(f: T^2\to S^2\), of degree 1, and a potential \(G: S^2\to\mathbb{R}\) of class \(C^{1,\beta}\) for all \(0\leq\beta<1\), such that \(f\) is harmonic with potential \(G\). Moreover, \(f\) can be chosen to be a constant map on an open set of \(T^2\).