We present a novel method for the calculation of the energy density of states D(E) for systems described by classical statistical mechanics. The method builds on an extension of a recently proposed strategy that allows the free energy profile of a canonical system to be recovered within a pre-assigned accuracy, [A. Laio and M. Parrinello, PNAS 2002]. The method allows a good control over the error on the recovered system entropy. This fact is exploited to obtain D(E) more efficiently by combining measurements at different temperatures. The accuracy and efficiency of the method are tested for the two-dimensional Ising model (up to size 50x50) by comparison with both exact results and previous studies. This method is a general one and should be applicable to more realistic model systems.It has long been recognized that all energy-related thermodynamic properties of a classical canonical system can be calculated once the energy density of states D(E) is known. In fact, starting from the partition function Z = dE D(E)e −βE , quantities like free energies, specific heats and phase transition temperatures can be computed in a straightforward manner. In principle, all canonical averages can be calculated through a multidimensional density of states. Due to this central role in equilibrium thermodynamics a variety of theoretical and computations studies have addressed the problem of how to obtain D(E) (or, equivalently, the entropy profile S(E) = ln(D(E)) in a reliable and efficient way.In principle, D(E) could be calculated from the histogram of the energies visited during a single "very long" constanttemperature simulation [1]. In practice, for any finite simulation only a limited energy windows is sampled so that the recovery of the system thermodynamics over a wide temperature range is unfeasible.Several alternative strategies have been developed to remedy this shortcoming. For example the multiple histogram reweighting technique relies on performing several simulations at different temperatures, so as to explore different (overlapping) energy intervals [2]. The various histograms are then optimally combined to obtain D(E) over the union of the energy intervals. Another successful family of techniques aims at obtaining D(E) by changing iteratively the probability with which the various energy levels are visited in stochastic dynamics until the recorded energy histogram is "flat". Such methods include entropic sampling [4], multicanonical and broad-histogram techniques techniques [3,5] and also the recent method of Wang and Landau [6].These and similar techniques have proved to be very valuable in a variety of contexts [7,8,9], but there is still ample scope for alternative approaches that could provide improved efficiency and better error control. Here we propose to modify and extend the metadynamics method recently introduced by two of us [10] for evaluating D(E). The algorithm we introduce allows a good error control on the explored energy range. Moreover, although within our approach D(E) could be reconstructed ...