A dynamical generalization of the instantaneous normal mode (INM) theory of liquid state dynamics is presented. Due to anharmonicities the eigenvalues and eigenvectors of the Hessian matrix change with time. Therefore, regular INM theory gives a description of molecular dynamics valid only for short times. Starting out from the classical equations of motion the velocity correlation function is expressed in terms of a series of propagation matrices. These are calculated by diagonalizing the Hessian matrix at configurations equidistant in time along a short piece of trajectory. Correlation functions calculated by this normal mode propagation (NMP) for a representative selection of atomic systems agree quantitatively with results, from molecular-dynamics simulation.