Pair Style Temperature instability

Hi all,

As many of you know we are trying to create a new pair style, which in many ways has been surprisingly challenging. The current problem comes as a result of a solution to a variety of other smaller challenges. In our pair style bonds are no longer completely independent of one another. The minimum depth, distance, and force constant for each bond depends on the remaining valence of the atom it is attached to. So if you have an ion that wants three bonds worth of bonds and has two bonds worth of bonds other than the bond in question, the bond in question bond would be a single bond and would have the ideal length depth and distance set accordingly. While this sounds really good in principle, has great thermochemical accuracy, it turns out the forces are extremely sensitive to the small changes in the other bonds. In effect a small deviation in one of the other bonds makes the current bond think its a little further out from equilibrium than it actually is. Then in the next step something else will be a little further out than it really is and so forth like a round robin. Minimization works perfectly but when you give it a little push its like riding around on greased rails. Even though the temperature thermosetting is set to every time step the temperature still reaches 1000-2000K for a room temperature simulation.

We think the solution is a kind of damping factor where we can average over the normal vibrational oscillations of the bonding. This way each bond would see the equilibrium position rather than the instantaneous position of the other bonds. I am guessing a factor on the order of 10-100 would be sufficient assuming a time step in the range of 0.1 to 0.5 fs.

Two questions:

  1. Does this make sense as a solution to the question of the excess temperature and disorder that is occurring?
  2. How would I code such a thing?

thanks much

Before you modify your force field, I suggest you run dynamics

with it (not minimization). If you cannot run NVE with a suitable

timestep and conserve energy for reasonable timescales, then

it indicates the forces/energies are not consistent. If that

is the case, then minimization will often have hard-to-debug



Hi Matthew,

This sounds to me like an energy conservation problem where the forces might have a wrong sign or a slightly incorrect derivation. Could you run NVE with force terms turned on one by one? That would probably help in isolating the responsible force term. If all derivatives are correct, could you do a time step sensitivity test? It may be the case that you just need a smaller time step considering the complexity of the functional form.

No comment for your proposed solution, but it may be difficult to code it up. If I am understanding you correctly, then you need information from 10-100 steps earlier to solve the current time step. This is not something LAMMPS is set up to do. Maybe you could try fix store/state? You will need to read the doc page to see if it will work for you.