An interpolation scheme for potential energy surfaces is presented. It employs a regular grid and finite element interpolation. The aim is the reduction of the computational expense for molecular dynamics simulation with a quantum chemical potential energy function. The methods used are described in detail. The feasibility is demonstrated and the efficiency and accuracy are evaluated for the photoisomerization of cis-stilbene in supercritical argon, using an ab initio configuration-interaction treatment for the first electronically excited state of the stilbene molecule and classical force fields for the solvent-solute interactions (quantum mechanical/molecular mechanical molecular dynamics). The number of required quantum chemical calculations of energy and gradients was substantially reduced compared to a simulation not using the interpolation scheme. On the other hand, the impact on the accuracy is insignificant.