Accurate and efficient modeling of flexible molecules in realistic environments is key for understanding most molecular mechanisms in chemistry, biology, and materials science. Existing methods for simulations of molecules consisting of more than a few dozen atoms involve a tradeoff between accuracy and efficiency. With the advent of high-level quantum chemical methods and machine learning, it becomes possible to develop novel simulation techniques that are free of this tradeoff. The goal of this project is to enable studying flexible molecules containing up to 100-200 atoms with accuracy close to the “gold standard” provided by high-level quantum chemical approximations to the solution of the Schrödinger equation. In practice, we will enable essentially exact molecular dynamics simulations where both electronic and nuclear degrees of freedom will be treated fully quantum-mechanically. This ambitious task requires building highly accurate and data-efficient machine learning models by constructing data-driven molecular descriptors, learning their metric, and creating optimal reference datasets as well as the development of appropriate enhanced sampling techniques. Our methods will be applied and tested on a number of molecules of practical interest, such as tripeptides, molecular switches, and drugs. The developed formalism will be implemented in an open source code making it available for other researchers and industry.