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.
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)
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.
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.
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.
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.
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.
