Recently I was using reaxff to simulate fluids around a metal slab. The problem is that it was just stuck in the middle of running. I found that there were nan in the thermo output and it continued to run with no more output until I manually checked. The last lines looks like this,
Step Temp Press KinEng PotEng TotEng
168840 612.86829 243.2728 1825.0171 -93146.435 -91321.418
168841 610.58581 215.77556 1818.2203 -93139.95 -91321.73
168842 608.39073 120.71219 1811.6837 -93133.474 -91321.791
168843 nan nan nan -93127.986 nan
You will find it interesting that nan seems to be related to velocities because there is no problem with potential energy.
I searched online and though it was the problem of my initial structure. Then I used simulated annealing to get a new structure and it still happened again. I have tested under 300K, 600K and 1000K. The problem only happened at 600K and 1000K. I tried with larger Tdamp for berendsen thermostat, but it didn’twork for me.
How should I fix this?
Thanks a lot!
Not very unusual. This means that some force was bad and had a NaN.
When the issue happens after over 150000 steps it cannot be the initial structure.
This hints at using too large a timestep. For higher temperatures, you have higher velocities and that means, that the discretization of the differential equations in the time integration may become too large and then atoms “can bump too far into each other” so that the time integration diverges.
A Berendsen thermostat is rarely a good choice. If you need to quickly remove, dissipate or instill kinetic energy and dissipative thermostat like one of the many langevin variants should be more effective and if you want to run in equilibrium, a Nose-Hoover thermostat is usually very much preferred.
Thanks for the detailed explanation. I will try with smaller timestep and N-H thermostat.
Have a nice day!
Just to update. It seems that all problems are not related to my timestep and thermo settings because all my following attemps failed.
After switched to another force field potential, things went smoothly well.
Neither of the papers you refer to have any mention of water. So it is not likely that the metal-water interaction is properly represented. In fact, I would worry about water in general.
Have you run the metal slab by itself at the given temperatures to see if it remains stable and reproduces expected properties?
Same for the water. Mind you at high temperature you are close to or beyond the super critical range where water will behave rather different from normal.
Thanks for your reply professor. I can’t be more happy to communicate with master like you.
Yes, I run the metal slab with both force field without water. Although it didn’t remain the slab structure (to be like a nano cluster), it was still held together rather than spread in the box.
So, is it an issue of force field? But I checked their reaxff force filed files, elements like C H O are all in the file. How should I understand this as a beginner of reaxff?
Thanks in advance.
I am neither a professor not a “master”.
It cannot be spread around without a solvent, but the fact that it melts is not a good sign if your simulation is maintaining a temperature significantly below the melting point.
Difficult to say. To rule out any issues from other simulation settings and confirming that you are using the parameters correctly, you need to reproduce the reference data from the paper describing the parameter set.
Unfortunately that doesn’t say anything about their suitability. There is a habit of ReaxFF parameter developers to just make a copy of an existing file created for a different purpose and then just modify/update the entries relevant to the problem at hand. Thus you can have elements in the file that may be correct on their own, but are not confirmed for the cross interactions.
Using ReaxFF is tricky business sometimes. The only way to be certain is to make many tests of specific test systems (e.g. confirm behavior of the bulk metal against some other type of potential, say EAM, them compare behavior of those for a slab, same for water as solvent and then a combined system where you employ pair style hybrid instead of ReaxFF).
At the risk of saying something very general (and therefore possibly wrong in your situation) and mundane (and therefore something you possibly already know), if you are not trying to study a chemical reaction with force field molecular dynamics, you should not use a chemically reactive force field like ReaxFF.
If I were trying to study a chemical reaction with force field molecular dynamics, I would always stop and ask myself if I would be better off using density functional theory instead. But that’s just me.
In this particular case, the paper you cite says that “The curve for the [pure Ag] simple cubic structure is poorly matched but is very high in energy”. So perhaps it is not surprising that the metal slab in your simulation is not behaving realistically.
It would be good for you to look up papers on molecular dynamics for water-metal systems and see what other people are doing. There are almost certainly already good force fields out there which won’t even require pair hybrid.
All the best!
Thanks for detailed explanation. I am trying with L-J and EAM force field to first equilibrate the system and then switching to reaxff for chemical reaction.
Thanks for reply.
DFT may be a better choice, but sadly now I am not capable of that.
Yeah, I think they developed the potential parameters for their condition ONLY. That’s why it didn’t work properly for my case.
My priority is chemical reaction so that I have to use reaxff. Anyway, I will find more papers to find a general idea here in order to do my simulation. (The trickiest but also the most interesting part )
Have a nice day!