How to perform non-reactive MD simulation at a high temperature by ReaxFF

[Lance_Feng]

Dear LAMMPS Users,

Basically, I want to perform a non-reactive MD simulation by ReaxFF at a very high temperature (i.e. 2000K or 3000K) to equilibrate the system. More specifically, it is about hydrocarbon oxidation (say methane oxidation) ReaxFF simulations. However, it is obviously that the system will react at such a high temperature. I searched the LAMMPS mail list and found several similar questions asked by other users but I am still very confused and don’t know how to do it.

I read several relevant papers doing non-reactive MD simulation to equilibrate the hydrocarbon oxidation system and got below descriptions:

1. “The initial system was firstly equilibrated at 298 K for 50 ps to get the reasonable initial configuration. After the equilibrium, the system temperature was further raised to 2000 K from 298 K within 50 ps at uniform rate. The NVT canonical ensemble was employed for these two stages to control the volume and temperature of the system. During the first two stages, no reaction between molecules was observed. After these two stages, microcanonical ensemble (NVE) was used to explore the detailed reaction process between CH4 and O2.”

2. “All the systems were minimized by molecular dynamics at low temperature (5 K) using the NVE ensemble. Then, they were equilibrated with the NVT ensemble for 10 ps at 1000 K using a time step of 0.1 fs.”

3. “The simulation procedure began with minimization of each system using low-temperature MD. Next, the system was equilibrated at 2500 K for 100 ps using a NVT MD simulation with a MD time step of 0.1 fs. During the equilibration simulations, the C-O and H-O bond parameters were switched off to prevent reactions from occurring during the equilibration of the system.”

4. “Reactions between O-C and O-H were prevented during the 100-ps equilibration period by eliminating the bond parameters describing these interactions from the force field.”

Therefore, based on the above descriptions, I believe that there is a certain way to perform a non-reactive ReaxFF simulation especially by switching off or eliminating some parameters. I am wondering that anyone know how to do it? It seems that I need to modify the force field file a little bit. My another question is that why point 3 & 4 only turn off the C-O and H-O reactions but not all bond reactions?

Sorry for the long question. I would appreciate it very much if anyone can help me with it. Many thanks in advance!

Kind Regards,

Lance

[akohlmey]

Dear LAMMPS Users,

Basically, I want to perform a non-reactive MD simulation by ReaxFF at a
very high temperature (i.e. 2000K or 3000K) to equilibrate the system. More
specifically, it is about hydrocarbon oxidation (say methane oxidation)
ReaxFF simulations. However, it is obviously that the system will react at
such a high temperature. I searched the LAMMPS mail list and found several
similar questions asked by other users but I am still very confused and
don't know how to do it.

I read several relevant papers doing non-reactive MD simulation to
equilibrate the hydrocarbon oxidation system and got below descriptions:

1. "The initial system was firstly equilibrated at 298 K for 50 ps to get
the reasonable initial configuration. After the equilibrium, the system
temperature was further raised to 2000 K from 298 K within 50 ps at uniform
rate. The NVT canonical ensemble was employed for these two stages to
control the volume and temperature of the system. During the first two
stages, no reaction between molecules was observed. After these two stages,
microcanonical ensemble (NVE) was used to explore the detailed reaction
process between CH4 and O2."

2. "All the systems were minimized by molecular dynamics at low
temperature (5 K) using the NVE ensemble. Then, they were equilibrated with
the NVT ensemble for 10 ps at 1000 K using a time step of 0.1 fs."

3. "The simulation procedure began with minimization of each system using
low-temperature MD. Next, the system was equilibrated at 2500 K for 100 ps
using a NVT MD simulation with a MD time step of 0.1 fs. During the
equilibration simulations, the C-O and H-O bond parameters were switched
off to prevent reactions from occurring during the equilibration of the
system."

4. "Reactions between O-C and O-H were prevented during the 100-ps
equilibration period by eliminating the bond parameters describing these
interactions from the force field."

Therefore, based on the above descriptions, I believe that there is a
certain way to perform a non-reactive ReaxFF simulation especially by
switching off or eliminating some parameters. I am wondering that anyone
know how to do it? It seems that I need to modify the force field file a
little bit. My another question is that why point 3 & 4 only turn off the
C-O and H-O reactions but not all bond reactions?

Sorry for the long question. I would appreciate it very much if anyone can
help me with it. Many thanks in advance!

​an alternate approach would be to define harmonic restraints via defining
explicit bonds (and angles/dihedrals if needed), but doing the
calculations​ without exclusions, i.e. using:

special_bonds lj/coul 1.0 1.0 1.0

​if the bond lengths for the explicit bonds are chosen well, there should
be little difference for the system during equilibration and you can
smoothly turn the restrains on or off through modifying the force constant.
with the all-inclusive special_bonds settings, you are still evaluating the
full ReaxFF interactions, and just have the bond restraints added on top of
it.

playing around with the force field parameters may work, too, but requires
a sufficiently deep understanding of ReaxFF and the meaning of individual
parameters. there are a few folks with ReaxFF experience on this list, so
they may want to comment further.
​
​axel.​

[Lance_Feng]

Dear Axel,

Thank you very much for your reply and suggestions! I will have a try.

Kind Regards,

Lance

[srdjan_pusara]

Hi Alex , I have one question regarding your advice to define explicit bonds , and then use special_bonds lj/coul 1.0 1.0 1.0 .

How can I define explicit bonds , in reaxff simulations atom style is charge , I could not google how to define explicite bonds. Is it by create bonds command ? But create bonds cannot be used with atom style charge ?

Also on a page for special bonds command it says that LAMMPS will ignore the special bonds settings when many body potentials are used , so I dont understand how this apeoach works :

“These weighting factors are NOT used by pair styles that compute many-body interactions, since the “bonds” that result from such interactions are not permanent, but are created and broken dynamically as atom conformations change. Examples of pair styles in this category are EAM, MEAM, Stillinger-Weber, Tersoff, COMB, AIREBO, and ReaxFF. In fact, it generally makes no sense to define permanent bonds between atoms that interact via these potentials, though such bonds may exist elsewhere in your system, e.g. when using the pair_style hybrid command. Thus LAMMPS ignores special_bonds settings when manybody potentials are calculated.”

I have one more question, if during equilibration I choose to turn off bond parameters in ReaxFF in order to prevent reaction to happen , if I have relatively simple system , and only reaction can be bond forming between C and O atom , that means that in ReaxFF parameters file I should set C-O bond parameters to 0 ? Can someome help me how I exactly need to modify bond parameter in order to prevent reaction from happening ?

Thank you .

[akohlmey]

  Hi Alex , I have one question regarding your advice to define
explicit bonds , and then use special_bonds lj/coul 1.0 1.0 1.0 .

    How can I define explicit bonds , in reaxff simulations atom style is
charge , I could not google how to define explicite bonds. Is it by create
bonds command ? But create bonds cannot be used with atom style charge ?

​in the case you are referring to, where the goal is to protect selected
bonds with an additional restraining force, the move convenient way to add
such bonds is via a data file.

yes, you cannnot use atom style charge, but you can use atom style full,
which is a superset of atom style charge *and* allows to store bond info.​

    Also on a page for special bonds command it says that LAMMPS will
ignore the special bonds settings when many body potentials are used , so I
dont understand how this apeoach works :

  *"These weighting factors are NOT used by pair styles
<http://lammps.sandia.gov/doc/pair_style.html> that compute many-body
interactions, since the “bonds” that result from such interactions are not
permanent, but are created and broken dynamically as atom conformations
change. Examples of pair styles in this category are EAM, MEAM,
Stillinger-Weber, Tersoff, COMB, AIREBO, and ReaxFF. In fact, it generally
makes no sense to define permanent bonds between atoms that interact via
these potentials, though such bonds may exist elsewhere in your system,
e.g. when using the pair_style hybrid
<http://lammps.sandia.gov/doc/pair_hybrid.html> command. Thus LAMMPS
ignores special_bonds settings when manybody potentials are calculated."*

​this paragraph is missing one important detail. it is correct to say,

that many-body styles are ignoring the scaling factors from the
special_bonds command.
however​, that is not the whole story. when you have a system with explicit
bonds, the presence of these bonds will affect the construction of the
neighbor list. specifically, for special_bond scaling factors of 0.0 (the
default), bonded pairs will be removed from the neighbor list, which will
render the many-body pair force computation invalid, since those force
kernels not only look at individual pairs, but also use the neighbor list
to look up neighbors of neighbors and beyond. hence the need for
special_bonds lj/cut 1.0 1.0 1.0

I have one more question, if during equilibration I choose to turn off
bond parameters in ReaxFF in order to prevent reaction to happen , if I
have relatively simple system , and only reaction can be bond forming
between C and O atom , that means that in ReaxFF parameters file I should
set C-O bond parameters to 0 ? Can someome help me how I exactly need to
modify bond parameter in order to prevent reaction from happening ?

reaxff is a highly complex force field. i advise against editing the
parameter files.

axel.