Hi all,
I was going to simulate a polymer in a solvent (water). To make sure
that the time step that I had chosen for rRESPA run_style was OK, I ran
some nve simulations. I found out that even when I use Verlet
integrator, the energy does not remain constant. The time step chosen
for the Verlet run_style was 0.5 fs, which is supposed to conserve the
energy. Can anybody help me with this one? What do you think is the
problem?
Thanks
Mehdi Hamaneh
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
involving hydrogens.
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.
Paul
Mehdi,
A few more questions and things to consider:
Are you using an LJ cutoff or a shift? What is the
cutoff?
Are you using charges? If so, are you using PPPM?
If none of these suggestions seem to solve the
problem, you might send us your input script and data
file (if not too large).
Paul