The nonlinear optical response of molecules is calculated using classical equations of motion for the electronic density matrix coupled to nuclear coordinates. The equations apply to an arbitrary (gauge-invariant) parametrization of the electronic density matrix. The electronic problem is treated at the time-dependent Hartree-Fock level, but the extension to time-dependent density functional theory is straightforward. An expansion in the inverse nuclear masses provides a low-cost calculation of polarizabilities at intermediate frequencies thigh compared to nuclear vibrations and low compared to electronic transitions). This regime is particularly relevant for optical materials applications. Closed expressions are derived for the electronic and nuclear contributions to the lowest two polarizabilities alpha and beta.