Hi Mehdi. I have several thoughts on this that might
be useful to you:
1) 0.5 fs might not be small enough if you have high
frequency oscillations, i.e. bonds vibrations
2) You might consider using shake to constrain bonds
involving hydrogens. Typically, people use 2.0 fs
timesteps with shake. This is typically done without
respa. Doing it without respa is easier to set up and
ensure good energy conservation.
3) What criteria are you using for "good energy
conservation"? All criterion I've seen are somewhat
arbitrary. No matter what, there will be some
fluctuation in energy due to numerical drift.
4) If your initial configuration has not been well
minimized --- equilibrated, there will typically be
much more energy drift. An easier (better?) test would
be to run some equilibration, for say 10 ps, then
start your energy conservation test.
With these thoughts in mind, plus my experience with
similar systems, I'd recommend using shake on the
bonds involving hydrogrens (and water angles), use a 2
fs timestep (assuming you're near 300 K),
equilibrating before testing for energy conservation,
and coming up with a good metric for reasonable
conservation. The one I usually use is something like:
log(average(abs((E(t) - E_ave)/E_ave)) <= -2.5
The run should be at least 10 ps long. I don't
particularly like this criteria since it seems quite
arbitrary and biased in favor of larger systems, but
it seems to be what people use.
I hope this helps.