simtk.org

I was troubleshooting the fact that our system ends up with a very large Ep and Ek after running local energy minimization, and I came across a very strange observation.

I tried advancing the integrator step-by-step after running local energy minimization. I'm using Variable Langevin Integrator. With a system I'm currently simulating, I'm observing the following timesteps and energies:

At the end of energy minimization, Ek = 1.5 kT/particle; Ep = 0.02 kT/particle
After one timestep, the time advances by 815 fs. Kinetic and potential energy shoot up: Ek=94.92 kT/particle; Ep=12.68
After the second timestep, the time advances by 17 fs, energies are: Ek=0.97 Ep=12.83

We expect the timestep to be at most 50-100fs given the potentials we have. So it looks like the first step of the integrator was very large, which leads to large integration error and hence to very large energies. This is with error_tol of 0.001 - issue persists with error_tol of 0.0001 as well, and we really need 0.01 for this system.

It also looks like if I decrease tolerance of energy minimization(0.001 vs 0.3), then the first step of variable langevin integrator ends up even larger (1500 fs for example), and the energy shoots up to 1000+kT/particle. Not sure why exactly this happens - maybe variableLangevin measures the forces to be too small and adjusts the timestep to be too large? On the contrary, if I don't do the energy minimization at all, in this particular system everything goes just fine: the system starts and continues with 40fs timestep, extra 3kT/particle of potential energy we've been trying to minimize is pumped into kinetic energy and gets absorbed by the thermostat. Yet we often need to dissipate much more Ep at the start of a simulation, and then this becomes an issue.

I tried using setStepSize of variableLangevinIntegrator, but it does nothing, and the description warns you that this might be the case.

I can try to make a self-contained example, but it will take time to make it because I'm running this using our polymer simulation library.

I can see why this would happen. The variable step size integrator assumes there's nothing "special" about the current coordinates. It assumes the forces at the start of the step are representative of the forces in the local region of configuration space. Or more technically, it uses an error estimate that's only accurate to first order.

When you've just done an energy minimization so you're exactly at the bottom of an energy minimum, that assumption isn't correct. The force is tiny, and to first order it doesn't grow as you move away from the minimum. The lowest non-vanishing term in the energy is second order in the distance from a local minimum. So that leads the integrator to think it can safely take a larger step than it can.

And of course, the more precisely you minimize the energy (so the closer you are to the exact minimum), the worse the problem gets.

That makes sense! And that does affect the start of nearly every simulation we've done in the past, but for most of them it won't be a problem.

I will work around this by running a fixed timestep integrator for this project, where precise alignment at the starting conformation is very relevant.

It would be generally useful for us to have an easy workaround that I can just plug into the library, rather than forcing everyone to do a few steps with a different integrator and switch afterwards. Would it be possible to add a seed timestep argument for the Variable Langevin integrator? It looks like many problems happen in the first timestep only, so seeding it with a different timestep may help it a lot. Or maybe maximum timestep?

A maximum step size option sounds like a really good idea. I just opened a feature request for it: https://github.com/openmm/openmm/issues/2383.

Thank you!

I thought of another workaround: randomize particle positions by a small value after minimization and before starting a simulation. It should be sufficient for now.

The maximum step size feature is now implemented: https://github.com/openmm/openmm/pull/2384

Thanks for a quick fix!