Ways to start ab initio MD from classical MD

Question

I am running molecular dynamics simulations of water for testing purposes. The box is quite small, if you ask a guy running classical MD, and relatively large, if you ask a DFT guy: I have 58 water molecules in periodic boundary conditions.

To save CPU time, I am optimizing my cell with a classical force field before running the ab initio MD. I equilibrate the system classically at 300K for 1 ns, then take the last snapshot and use it as input for ab initio MD. My ab initio MD is regular DFT-based Born-Oppenheimer MD with a plane wave basis set and PAW (pseudo)potentials (VASP is the code). In both classical and ab initio simulations I am keeping the temperature constant at 300K using a velocity-rescaling thermostat.

I am surveying two different ways to make the transition between classical and ab initio:

1. Take the initial velocities and positions from the classical trajectory and import them as initial configuration for the ab initio simulation
2. Freeze the system to zero temperature keeping the classical positions, import that to the DFT code, then quickly (I'm doing it in 0.5 ps at the moment) heat up to 300K

I was hoping that both strategies would lead to the same average energy after a short (say 10 ps) equilibration period, especially considering that the starting configuration is exactly the same (same initial positions) except for the mentioned temperature trick (initial velocities differ). This is not the case. The figure below shows that the simulation where the system is frozen and then quickly heated up finds an energy region about 1 eV lower in energy than the other one, where the velocities where imported from classical MD.

My questions are:

1. whether this is to be expected;
2. are there known successful strategies to optimize the transition from classical to ab initio MD;
3. and could you point me towards pertinent literature on the matter?

Edit:

I have been running some more tests and -with the limited data I have at the moment- it seems that this might be a system-specific problem. A test with methanol instead of water in a box of the same size showed that the two different initial velocity schemes quickly converge to the same average energy. However, the classical configuration was very close to the quantum one in the case of methanol, that is, the energy at t = 0 was very close to the average energy after convergence. Water is a notoriously difficult system, so perhaps this problem is more or less water-specific. If no answers are added I will try and post one based on my results once I'm done with all the testing.

Answer 1

after a short (say 10 ps) equilibration period

1. You say yourself that the re-equilibration period is 'short'. Have you tried waiting longer to see if the two re-initialized systems converge, and if so at what rate?
2. Velocity-rescaling is a notoriously naive thermostat. Maybe you could replace it with something more realistic? (Berendsen, Nose-Hoover, etc.)
3. If you're at all concerned about sampling a 'large' system, you can make the ergodic principle and parallel computing work for you: Sample multiple re-initialized realizations from the tail of your classical trajectory (taking snapshots of the classically equilibrated system). Then you can run multiple instances of your quantum code, each using a different re-initialization phase-space configuration, and average the outputs. Since the separate realizations evolve independently of each-other, the simulations are embarrassingly parallel.

Sorry I can't offer any authoritative sources...