Proteins are found in solution as ensembles of conformations in dynamic equilibrium. Exploration of functional motions occurring on μs-ms time scales by Molecular Dynamics (MD) simulations still remains computationally challenging. Alternatively, normal mode (NM) analysis is a well suited method to characterize intrinsic slow collective motions, often associated to protein function, but the absence of anharmonic effects preclude a proper characterization of conformational distributions in a multidimensional NM space. Using both methods jointly appears as an attractive approach that allows an extended sampling of the conformational space. In line with this view the MDeNM (Molecular Dynamics with excited Normal Modes) method presented here consists in multiple-replicas short MD simulations in which motions described by a given subset of low frequency NMs are kinetically excited. This is achieved by adding additional atomic velocities along several randomly determined linear combinations of NM vectors, thus allowing an efficient coupling between slow and fast motions. The relatively high-energy conformations generated with MDeNM are further relaxed with standard MD simulations enabling to determine free energy landscapes. Two widely studied proteins were selected as examples: hen egg lysozyme and HIV-1 protease. In both cases, MDeNM provides a larger extent of sampling in a few nanoseconds, outperforming long standard MD simulations. It was noticed a high degree of correlation with motions inferred from experimental sources (X-ray, EPR and NMR), and with free energy estimations obtained by metadynamics. Finally, the large sets of conformations obtained with MDeNM can be used in a better characterization of relevant dynamical populations, allowing a better interpretation of experimental data such as SAXS and NMR spectra.