Potential energy surfaces (PESs) for use in dynamics calculations of few-atom reactive systems are commonly modeled as functional forms fitting or interpolating a set of ab initio energies computed at many nuclear configurations. An automated procedure is here proposed for optimal configuration-space sampling in generating this set of energies as part of the grid-empowered molecular simulator GEMS (Laganà et al., J. Grid Comput. 2010, 8, 571-586). The scheme is based on a space-reduced formulation of the so-called bond-order variables allowing for a balanced representation of the attractive and repulsive regions of a diatom configuration space. Uniform grids based on space-reduced bond-order variables are proven to outperform those defined on the more conventional bond-length variables in converging the fitted/interpolated PES to the computed ab initio one with increasing number of grid points. Benchmarks are performed on the one- and three-dimensional prototype systems H2 and H3 using both a local-interpolation (modified Shepard) and a global-fitting (Aguado-Paniagua) scheme.