ReaxFF pyrolysis at 2500 K: occasional temperature spikes and using fix nve/limit (how to judge reliability / reduce spikes)?

[ghuanhuan9]

Hi all,

I am performing ReaxFF MD simulations to study the thermal decomposition of a new energetic material (DAP-4) at 2500 K.

I have observed occasional instantaneous temperature spikes. While I understand that bond breaking in energetic materials releases significant kinetic energy, I am concerned about the reliability of the thermostatting strategy and potential artifacts from fix nve/limit.

LAMMPS: [LAMMPS (27 Jun 2024)]

What I see

1.With Berendsen-type thermostat + fix nve/limit 0.1 (earlier tests), I sometimes get large spikes up to ~5000–6000 K (Fig1) .

2.With the current setup (Langevin + fix nve/limit 0.1), spikes are smaller (about~ 3000-3600K) (Fig2).

3.In this 50-ps run, fix nve/limit 0.1 is active for ~5000 steps out of 500000 (~1%) (Fig3).

4.Impact on Reaction Initiation: I also compared the time evolution of the total number of species between the two thermostats (Fig4).

Main questions (only two)

1. For ReaxFF pyrolysis at 2500 K, are brief temperature spikes (even if rare) considered acceptable, or do they indicate numerical instability that can bias products/pathways?
2. To avoid temperature fluctuations, is fix langevin + fix nve/limit 0.05 a reasonable production approach? With ~2% of steps being limited, is it likely to bias kinetics/chemistry? If not recommended, what is the first thing should change to reduce spikes (timestep, QEq settings, neighbor rebuild, thermostat damping, etc.)?
   
   

input.in (2.1 KB)

Thank you for your time and advice.

Bowen

[akohlmey]

@ghuanhuan9

are brief temperature spikes (even if rare) considered acceptable, or do they indicate numerical instability that can bias products/pathways?

That is a question for somebody who has experience in running ReaxFF simulations under those conditions. There is obviously nobody here participating in this forum that has this experience and is willing to share it with you. There are many “lurkers” in public forums like this one, so it is no big surprise to me. As a non-expert, I am always suspicious and since you already using a rather short timestep, I would suspect that there may be a problem with the force field parameters that you are using. They may not be suitable for your simulation conditions or just not suitable at all.

@ghuanhuan9

is fix langevin + fix nve/limit 0.05 a reasonable production approach?

No. Fix nve/limit should only be used during equilibration.

@ghuanhuan9

If not recommended, what is the first thing should change to reduce spikes (timestep, QEq settings, neighbor rebuild, thermostat damping, etc.)?

As mentioned above, you have to check with what experts in that field of research have been doing. From my perspective, you are leaving out the most suspicious component and that would be the force field parameters. Since all your other choices already seem conservative, the only other option would be to reduce the time step from 0.1fs to an even smaller number. But that is quite suspicious.