Therefore, nonadiabatic dynamics simulations can be performed using only energies and gradients. This simplicity makes it attractive for excited-state methods, where NACs are not available or too expensive to compute. It is worth mentioning that the Landau-Zener hopping probability reflects the inverse proportionality of the NACs on the energy gap of the PE surfaces. Here we demonstrate that combining the interpolation capability of Py-Surf with the Landau-Zener algorithm opens up another attractive and cheap doorway for surface-hopping simulations. Having global PE surfaces at hand from the interpolation, gradients can be calculated from the fitted energies either by a finite-difference scheme or computing the analytic gradient of the fitting algorithm itself. The latter is common in machine learning libraries that support automated gradient generators. As the prediction of energies by interpolation is computationally cheap, it is also inexpensive to predict the gradients of the PE surfaces based on the interpolation. This approach, hereafter called energy-only simulations, allows one to perform nonadiabatic molecular dynamics simulation using only energies as input from ab-initio calculations. Nuclear gradients needed for the propagation are computed numerically or analytically derived from the fitting algorithm itself. This scheme is of particular interest for electronic structure methods for which the analytic gradients are not available. The general workflow has the following steps; i) Filling the database using typical sampling approaches, e.g. Wigner sampling, or performing nonadiabatic dynamics simulations as it is done in this work. ii) After the conformational space is sampled nonadiabatic dynamics are performed using only the fitting algorithm until the trajectory leaves the conformational space, where the fit is accurate. If this occurs within the simulation time, it is necessary to extend the database by performing new ab-initio calculations. Additionally, a so-called adaptive sampling scheme can be applied. Hereby, the fit is used, whenever it is accurate, else ab-initio calculations are performed and are used to extend the database. In the ideal case, the database contains only the datapoints needed to accurately fit the relevant conformational space, and should therefore not show an exponential increase with the number of degrees of freedom. However, this depends on the sampling approach and on the system at hand. The proposed workflow can be in principle applied in all cases as long as Landau-Zener surface hopping is valid, i.e. it is exact for the two state case and can be successfully used for transitions that mainly involve two states, which are reasonably well separated from other electronic states, whereas it conceptually fails when more electronic states are within a small energy range. Additionally, the selection of the training data for an accurate fit can be considered as another limiting factor, and so far no general solution to this problem exists.