The recent progress of linear-scaling or O(N) methods in density functional theory (DFT) is remarkable. Given this, we might expect that first-principles molecular dynamics (FPMD) simulations based on DFT could treat more realistic and complex systems using the O(N) technique. However, very few examples of O(N) FPMD simulations exist to date, and information on the accuracy and reliability of the simulations is very limited. In this paper, we show that efficient and robust O(N) FPMD simulations are now possible by the combination of the extended Lagrangian Born-Oppenheimer molecular dynamics method, which was recently proposed by Niklasson ( Phys. Rev. Lett. 2008, 100, 123004), and the density matrix method as an O(N) technique. Using our linear-scaling DFT code Conquest, we investigate the reliable calculation conditions for accurate O(N) FPMD and demonstrate that we are now able to do practical, reliable self-consistent FPMD simulations of a very large system containing 32768 atoms.