A semi-automatic sampling and fitting procedure for generating sum-of-product (Born–Oppenheimer) potential energy surfaces based on a high-dimensional model representation is presented. The adaptive sampling procedure and subsequent fitting rely on energies only and can be used for re-fitting existing analytic potential energy surfaces in the sum-of-product form or for direct fits from ab initio computations. The method is tested by fitting ground electronic state potential energy surfaces for small to medium sized semi-rigid molecules, i.e., HFCO, HONO, and HCOOH, based on ab initio computations at the coupled-cluster single double and perturbative triples-F12/cc-pVTZ-F12 or MP2/aug-cc-pVTZ levels of theory. Vibrational eigenstates are computed using block improved relaxation in the Heidelberg multi-configurational time dependent Hartree package and compared to available experimental and theoretical data. The new potential energy surfaces are compared to the best ones currently available for these molecules in terms of accuracy, including resulting vibrational states, required number of sampling points, and number of fitting parameters. The present procedure leads to compact expansions and scales well with the number of dimensions for simple potentials such as single or double wells.