A highly accurate and efficient method for molecular global potential energy surface (PES) construction and fitting is demonstrated. An interpolating moving least-squares (IMLS) method using low-density ab initio potential, gradient, or Hessian values to compute PES parameters is shown to lead to an accurate and efficient PES representation. The method is automated and flexible so that a PES can be optimally generated for classical trajectories, spectroscopy, or other applications. Two main drivers for the fitting method have been developed thus far. The first is a PES generator designed primarily for spectroscopy applications. Using this method, the configuration space defined by a specified energy range is automatically fit to a predefined accuracy. A second approach is based on trajectory methods for computing reaction rates. In this approach, the configuration space that is dynamically accessible to a particular ensemble of trajectories is fit "on the fly." Results that are indicative of the accuracy, efficiency, and scalability will be presented.