A formulation of a direct, iterative method for obtaining the lowest eigenvalues and eigenvectors of a Hessian matrix is presented. Similar to the iterative schemes in electronic structure configuration interaction calculations (methods due to Lanczos, Davidson, and others), the mass-weighted Hessian matrix K is not constructed explicitly; instead, its operation on a basis vector (a direction in the 3N Cartesian configuration space of the atoms) is computed based on the principles of dynamical equations of motion. By noting that the time derivative of the gradient vector in the harmonic force field is related to the particles’ momenta via dg∕dt=Kp, a Hessian-vector product can be computed on the fly by finite differentiation of the gradient along the direction specified by the p vector. Thus, only two evaluations of the gradient are required per Davidson-like iteration per root, which leads to a linear scaling behavior of the computational effort with the number of atoms (provided that the evaluation of the gradient scales linearly). Preliminary results are presented for a 27 000-atom He4 nanodroplet.