Temperature inconsistency between L-J and reaxff potential

Hello everyone,

My system contains fluid and metal slab so that I want to simulate the chemical reaction between them.

My procedure is to first equilibrate the system at 300K with L-J potential to prevent the possible reaction (because it happens very fast if I equilibrate it with reaxff potential), then write out the data file, which is the initial structure for the second reaction process with reaxff potential.

However, when I check the log, I find that there is a huge temperature inconsistency at the beginning of second reaction process. Here is the log file.

Step Temp
0 308.49983
1000 1362.5331
2000 751.99976
3000 528.48456
4000 450.15296
5000 410.03117
6000 388.0228
7000 363.29414
8000 358.55172
9000 342.41522
10000 324.99101

And the reaction actually happens even faster in the beginning this time because from the point of reaxff the initial temperature is far more higher than 300K. This huge temperature inconsistency is really unacceptable here. (I use the exactly same berendsen thermostat for both potentials)

Is there any possible way to avoid it?
Thank you.


One possible way to prevent reactions at the beginning of the simulation is to add extra harmonic bonds between atoms from the same molecules. Once the system is equilibrated and the temperature is reasonable, the extra bonds can be removed.


This is tricky business, since adding the bonds will change the neighborlist unless you also change the special_bonds setting. With the default setting, all non-bonded pair interactions that are directly bonded (1-2 pair), or connected with one (1-3 pair) or via two bonds (1-4 pair) are excluded. This makes a lot of sense for molecular force fields like Amber, CHARMM, etc. However ReaxFF needs to have those pairs included or else the interactions are computed wrong since ReaxFF does determine the network of bonded interactions internally and dynamically from the non-bonded pairs in every step.

The difference you observe is due to a significant difference in potential energy which means that your choice of LJ potential parameters and/or force field is quite far from those that you have with ReaxFF.
The way to avoid this is to find/derive “better” parameters. Otherwise your “equilibration” is actually counterproductive and will move your system away from equilibrium.

So the “temperature inconsistency” is actually not there, but rather you have a potential energy inconsistency (which is much worse).

Thanks for your idea. Is it compatible to add two bonds (one from extra and one from reaxff itself)?
Besides, what if the molecular structure is constrained by additional angels and dihedrals?

I get it. This orginates from how these two potentials define the molecules energy.
I am just thinking, how about equilibrating both fluid and metal slab with reaxff but in two seperate simulation box? After I get the postion and velocity, I can combine them together as I want. Because they are both acquired under reaxff potential, I assume there will be no temperature inconsisitency? Sounds reasonable?

Well, following @akohlmey’s comment, this is may be more tricky that I initially though. But the idea anyway is not to generate a physically correct molecular structure, but rather to prevent the chemical reactions from occurring during the very first steps of the simulation by forcing the atoms to remain bonded together. Only after these additional bonds are removed, the system should relax toward the correct molecular structure.

I’ve never tried it myself, but I remember that people have (successfully) used similar trick.


No. The problem here is that each system is equilibrated with a free surface, but then will suddenly be put in contact with something else. You will have to leave a bit of a gap so that there are no overlaps and then you will have to do another equilibration to equilibrate the interface and likely will have to content with a shockwaves caused by the high potential energy from joining of the system or from imploding the extra vacuum you leave there.

The way to go for your specific system is to define two groups, one for the metal, one for the fluid. Then run equilibration with a time integration and thermostatting to one of the groups only. Then equilibrate the other group only for a bit and alternate a couple of times. That should allow each part of the system to equilibrate while the other is frozen but because of the alternating, even the interface area can adjust to the other material being present without mixing.

OK. I will try to equilibrate two groups separately in the same simulation box.
Just asking, if I lower down the fix qeq/reax frequency during equilibration, say from every 1 timestep to 1000 timestep, will the reaction be slowed down so that I can get the initial structure and temperature as I want?
Thanks you sir.

Yes. This is a more tricky work than I thought.

My instinct says no. Your results show that the ReaxFF force field, including partial charges, is not quite compatible with the LJ force field.

What you suggest would freeze the partial charges per thousand steps, at one less-compatible distribution, instead of letting them vary between those less-compatible distributions per time step. Such a change doesn’t seem like it would change the incompatibility you’re seeing.

You can certainly try it and see!

At that point you could just do a single point optimization of the charges (= run 0) and then unfix the fix and remove the check for fix qeq in the reaxff pair style settings. Not sure how much that would influence any reactions, since the charge equilibration is only part of the ReaxFF force field.

Now I am getting very curious to learn what kind of reaction would be that fast that it will significantly taint your results already during equilbration. Very few explosives have reactions that are that fast. I would be very concerned that there is something else very wrong.

If ReaxFF Sodium Is Dropped In A Bathtub Of ReaxFF Water, Will LAMMPS Explode?


I will set the initial charge to be zero after L-J equlibration for further test.

Setting charges to zero is before running ReaxFF is a bad thing to do:

I am simulating the catalytic decompsition of hydrogen peroxide on silver, nickel and gold etc. metals.
The reason why I think it happens so fast at the beginning is that I first equilibrate the system at 300K for 60ps, following by an annealing process which heats the system up to 500K and cools it down to 300K in 10ps, then continue the reaction process for 500ps. The initial number of H2O2 molecules is 200. But when I check the number after equilibration, it drops to only about 70.
If I skip the equilibration, just run the reaction simulation, although the initial H2O2 number will be 200 as I want, but the temperature will go up to about 1000K in the beginning which will may accelerate the reaction more.
Considering the above results, I thought that the equlibration process is still needed, but must find a way to prevent the reaction happening as much as posssible.

Oh there is always something bad waitting for me. :upside_down_face:
I am wondering now that is it wrong to set the inital charge to be zero in the equilibration?

LAMMPS can be a complex machine, and like any other complex machine, there are often many ways to break things and only a few ways to get things right.

A useful exercise is to open your script, assume that every line contains some mistake (or is some mistake) that will kill your simulation, and either convince yourself that the line is correct or delete it.

Yes, I am trying to get used to this painful but interesting conquer journey. :smile: