Due to their high computational cost, the AIMD simulations with hybrid functionals were performed with a smaller cubic box of dimensions 9.939 3 Å 3 containing 32 water molecules. A multiple time step (MTS) scheme was used to further accelerate the simulations, with an inner time step of δt = 15 a.u. and an outer time step of ∆t = n • δt, where the time step ratio n is chosen to maintain sufficient energy conservation. At inner time steps, fast force components are given by a delta-ML model that predicts PBE0 forces based on the LDA (F inner = F LDA + ∆F PBE0-LDA ML ), while total forces are corrected at the outer time step with their slow components (F outer = F Minnesota -F inner ) to fully recover the higher-level Minnesota forces . In this approach, ML serves only as a low-level surrogate operating on shorter timescales without impacting the target DFT level. Note that the inner PBE0 level does not need to match the outer Minnesota level entirely, but should be close enough so that their difference slowly varies in time and dynamically decouples from fast force components. Ultimately, the Minnesota level is recovered at larger physical time steps by construction, ensuring that the structural and dynamical properties are not affected , unlike in ML-potential MD.