A new methodology to perform hybrid empirical/ab-initio molecular dynamics is presented. The method combines the well-established hybrid ONIOM scheme with the recently developed ADMP (Atom-Centered Density Matrix Propagation) approach, where the one electron density matrix expanded in an atom-centered Gaussian basis set is propagated as electronic variables along with the classical nuclear degrees of freedom via an extended Lagrangian procedure. The unified and single-valued ONIOM expressions for the energy and energy derivatives allow for an implementation of conservative dynamics. It is found that atom-centered basis sets large enough to provide good chemical accuracy can be used even when electronic embedding is adopted in the QM/MM scheme, and this does not affect the well-behaved and conservative nature of the dynamics. The method contains very appealing features for the study of biological systems, including the ability to employ accurate density functionals, the freedom to choose a periodic or a cluster boundary condition for the system under study, asymptotic O(N) scaling through established techniques, and the ability to use reasonably large time-steps through the tensorial fictitious mass scheme. The general ADMP/ONIOM formalism is illustrated through a series of test calculations. A, simulated study of proton hopping inside the gramicidin A ion channel is also presented, to show the potential of the method in describing reactivity in large systems.