We elaborate a novel and efficient method to obtain multiphase Helmholtz free energies from molecular dynamics (MD) simulations over wide ranges of volume and temperature in materials that can be described by temperature-independent ion forces, with both higher accuracy and order-of-magnitude cost savings compared to direct thermodynamicintegration techniques. Our method leverages and significantly extends the technique of reversible scaling molecular dynamics (RSMD) proposed by de Koning et al. [Phys. Rev. Lett. 83, 3973 (1999)], which allows a free-energy difference in a given phase at constant volume to be calculated as a function of temperature from a single MD simulation. In mechanically stable solid phases, our approach carefully combines quasiharmonic lattice dynamics at low temperatures with an accurate and fully isolated RSMD simulation of the anharmonic vibrational free energy at high temperatures to produce a seamless free energy from zero temperature to above melt along constant-volume isochores. In the liquid, we combine a unique calculation of the free energy along a high-temperature reference isotherm with isochoric RSMD simulations from that temperature to below melt. In metastable solid phases that are mechanically unstable at low temperature, we use two-phase MD melt simulations together with the liquid free energy to obtain the solid free energy along the solidus melt line and then perform isochoric RSMD simulations to temperatures above and below that point. While our free-energy method is general, we have specifically adapted it here to the case of metals in which the ion forces are well described by model generalized pseudopotential theory (MGPT) multi-ion interatomic potentials, and additive electron-thermal free-energy contributions can be 2 included. Then using refined Ta6.8x MGPT potentials, we have converged total free energies and their components to very high and unprecedented sub−milli-Rydberg (mRy) numerical accuracy in the stable-bcc, liquid, and metastable-fcc phases of tantalum for volumes ranging from up to 26% expansion to nearly two-fold compression and for temperatures to 25,000 K. In turn, we have successfully used the free energies so obtained to calculate physically accurate thermodynamic properties and gain new insight into their behavior, including sensitive thermodynamic derivatives, bcc and fcc melt curves, and a multiphase equation of state for Ta over the same temperature range and for pressures as high as 600 GPa. We show that the anharmonic free-energy component in the bcc solid, although only 1-5 mRy in magnitude for Ta, can have a significant (15-20%) effect on thermal expansivity, the Grüneisen parameter, and melt temperatures. We further show that the electron-thermal free-energy component can similarly impact the specific heat and thermal expansivity in both the solid and the liquid, while only minimally affecting (to ≤ 3%) the bcc and fcc melt curves.