A Lagrangian generalization of time-reversible Born- Oppenheimer molecular dynamics [Niklasson et al., Phys. Rev. Lett. vol. 97, 123001 (2006)] is presented. The Lagrangian includes extended electronic degrees of freedom as auxiliary dynamical variables in addition to the nuclear coordinates and momenta. While the nuclear degrees of freedom propagate on the Born- Oppenheimer potential energy surface, the extended auxiliary electronic degrees of freedom evolve as a harmonic oscillator centered around the adiabatic propagation of the self- consistent ground state. The formulation enables the application of higher-order symplectic or geometric integration schemes that are stable and energy conserving even under incomplete self-consistency convergence. I will demonstrate how the extended Born-Oppenheimer molecular dynamics improves the accuracy by over an order of magnitude compared to previous formulations at the same level of computational cost.