definition of bonds and angles when using reax/c

Hi all,

I would like to model organic reactions in various environments (notably water and organic solvents such as DMSO) with ReaxFF.

Due to the fact that in general with ReaxFF each atom is defined by only one atom type, my feeling is that it might be relatively difficult to find parameters that provide at the same time an accurate description for the reacting species, and that reflects well solute-solvent interactions. In addition I want to avoid reactions between molecules of solvent, or between the solvent and the solute (I have noticed that this problem often occurs leading to the decomposition of many solvent molecules during simulations).
To solve this, my first idea was to combine two force fields: a not-reacting one for the environment, and ReaxFF for the reacting species. Unfortunately, this is not straightforward as ReaxFF cannot be combined with another pair_style (due to the hybrid/reax incompatibility) which prevents replacing several ReaxFF terms (coulomb, vdw, and h-bond) with another pair_style.

So, I tried to define bonds (bond_style harmonic) and valence angles (angle_style harmonic) for all solvent molecules, and I turned off the corresponding ReaxFF terms by setting the value of several selected parameters to 0 in the ffield file (e.g. pval1 for Eval). Actually, new atom types must be added in the ffield file to describe the solvent molecules (I have attached the ffied file that I use, see the atom types 12 and 13 which corresponds to H and O respectively).

When I try to apply this procedure in simulations several problems that I’m not able to explain appear. First, I have noticed a segmentation fault that seems linked to the evaluation of the hbond energy. There is a kind of discontinuity in the evolution of this term: after reaching a given threshold, the hbond term becomes immediately very large (-nan). The error can be reproduced with the first input files attached (in.FF2, and data.FF2). In this simulation, there are two water molecules, one is fully described with ReaxFF terms, while the other is described with two bonds, one angle and the h-bond, Coulomb, and vdw terms of ReaxFF . The second problem concerns larger systems and can be reproduced with the files in.waterbox and data.waterbox (simply a box containing 500 water molecules defined with harmonic bonds, harmonic valence angles, and Coulomb, vdw and hbond ReaxFF terms). It is linked with convergence problem of the qeq fix. In fact, it seems that the problem is due to the coulomb term which rises very fast. It is worth to notice that this problem is not due to the starting atomic configuration since the qeq fix converged when no bonds and angles were defined.

So far, I did not found any satisfying explanation for any of these errors (are they linked together or not? is it even possible to define bonds and angles when using reax/c? ...)

Thanks in adavance for all the comments!

Pierre

ffield.reax (32.5 KB)

data.FF2 (495 Bytes)

data.waterbox (76.6 KB)

in.waterbox (1.79 KB)

in.FF2 (1.85 KB)

So, I tried to define bonds (bond_style harmonic) and valence angles (angle_style harmonic) for all solvent molecules, and I >turned off the corresponding ReaxFF terms by setting the value of several selected parameters to 0 in the ffield file (e.g. pval1 >for Eval). Actually, new atom types must be added in the ffield file to describe the solvent molecules (I have attached the ffied >file that I use, see the atom types 12 and 13 which corresponds to H and O respectively).

I can imagine all kinds of pitfalls in this approach, both
conceptually and mistakes
you might make in tweaking the ReaxFF file. I'll let Aidan comment on
whether this is possible or a good idea.

Steve

You have attempted to create a new parameterization of ReaxFF by manually
editing the ffield file. You are asking us why the resultant LAMMPS
simulations are exhibiting various kinds of numerical instabilities. Since
we don't see these problems with existing ffield files, I have to assume
this is due to a bad ReaxFF parameters. If you can demonstrate any of
these problems with a ffield file that is known to be "well-posed", I'll
take a look.

Aidan

Hi all,

I would like to model organic reactions in various environments (notably water and organic solvents such as DMSO) with ReaxFF.

Due to the fact that in general with ReaxFF each atom is defined by only one atom type, my feeling is that it might be relatively difficult to find parameters that provide at the same time an accurate description for the reacting species, and that reflects well solute-solvent interactions. In addition I want to avoid reactions between molecules of solvent, or between the solvent and the solute (I have noticed that this problem often occurs leading to the decomposition of many solvent molecules during simulations).
To solve this, my first idea was to combine two force fields: a not-reacting one for the environment, and ReaxFF for the reacting species. Unfortunately, this is not straightforward as ReaxFF cannot be combined with another pair_style (due to the hybrid/reax incompatibility) which prevents replacing several ReaxFF terms (coulomb, vdw, and h-bond) with another pair_style.

i don't think that even if there was no restriction on reax with
hybrid pair styles
that this would work correctly for you. pair style hybrid computes a simple
mechanical coupling, however you also need a proper electrostatic
embedding - same as with QM/MM couplings - to have your (polar!) solvent
properly interact with the reax part of your system.

So, I tried to define bonds (bond_style harmonic) and valence angles (angle_style harmonic) for all solvent molecules, and I turned off the corresponding ReaxFF terms by setting the value of several selected parameters to 0 in the ffield file (e.g. pval1 for Eval). Actually, new atom types must be added in the ffield file to describe the solvent molecules (I have attached the ffied file that I use, see the atom types 12 and 13 which corresponds to H and O respectively).

When I try to apply this procedure in simulations several problems that I’m not able to explain appear. First, I have noticed a segmentation fault that seems linked to the evaluation of the hbond energy. There is a kind of discontinuity in the evolution of this term: after reaching a given threshold, the hbond term becomes immediately very large (-nan). The error can be reproduced with the first input files attached (in.FF2, and data.FF2). In this simulation, there are two water molecules, one is fully described with ReaxFF terms, while the other is described with two bonds, one angle and the h-bond, Coulomb, and vdw terms of ReaxFF . The second problem concerns larger systems and can be reproduced with the files in.waterbox and data.waterbox (simply a box containing 500 water molecules defined with harmonic bonds, harmonic valence angles, and Coulomb, vdw and hbond ReaxFF terms). It is linked with convergence problem of the qeq fix. In fact, it seems that the problem is due to the coulomb term which rises very fast. It is worth to notice that this problem is not due to the starting atomic configuration since the qeq fix converged when no bonds and angles were defined.

So far, I did not found any satisfying explanation for any of these errors (are they linked together or not? is it even possible to define bonds and angles when using reax/c? ...)

i am surprised that you even managed to get it to work this far.

axel.

Thank you for those comments!
First, I understand that turning off several terms of ReaxFF and replacing them with others (of a different nature) as I try to do, can lead to conceptual difficulties. Notably because all terms are consistent together due to the parametrization. But at the moment, my main interrogation is knowing if it is even technically feasible or not.

Following Aidan's suggestion, I have tested whether the qeq convergence problem was still present when using one of the ffield file provided with LAMMPS. I took the C/H/O combustion force field of Chenoweth et al. (I did not modify it at all). The simulation of the box containting 500 water molecules still fail (at the very begining) if bonds and angles are defined (my input files are in attachements, as well as the corresponding log file).

In fact I've noticed that when bonds are defined, LAMMPS seems to automatically turn off ReaxFF terms in which bond orders take part, that is all terms except the vdw and the coulomb ones.
You can see for example in the attached log file that the value for the ReaxFF bonding term (eb) at step 0 is abnormally low (-89.8 kcal/mol) for a system containing 500 water molecules. I believe that this value is only due to interactions between atoms of distinct molecules that are close enough so that significant bond orders are computed, and for which harmonic bonds are not defined.
Is this suppression of ReaxFF "bonding" terms an expected behavior? And if yes, is there a way to circumvent it and allow LAMMPS treating at the same time a bond_style and all ReaxFF terms?

A last point, when I decrease the cut-off value of the qeq fix to 4.5 or 4.0 A (I know it's very low), the rapid increase of the Coulomb term still occurs, but a bit later in the simulation.

Pierre

log.waterbox (3.71 KB)

data.waterbox (106 KB)

in.waterbox (1.79 KB)

You appear to be calculating double contributions to bond and angle
energies: once from ReaxFF and one from the harmonic potentials that you
have specified. In addition, this ReaxFF potential is for energetic
materials. I do not expect it to perform well for liquid water.
Nonetheless, you should first trying running just ReaxFF, without any
extra bond and angle terms.

Sorry Aidan, I forgot to mention that I use this waterbox simulation as a test, I do not try to extract any meaningful physical value from it at this point.
Of course, I have tried the same simulation without defining bonds and angles, and it works well (once more I do not know if it reproduces well the behavior of liquid water but the simulation runs normally).
I would like to model non-reacting solvent molecules surrounding reactive species. One of the major problem is that solvent molecules tend to overreact (with each other or with the solute), and this holds true with all the potential files that I have tested so far. I think that defining fixed bonds and angles for solvents molecules could make them non-reactive. That is why I am interested in using ReaxFF in combination with a bond_style and an angle_style, which seems to be disallowed by LAMMPS as it can be seen from the log file of my previous message. So, I wonder if this is expected or not (there is no error or warning message but ReaxFF bonding terms are obviously set to 0 when bonds are defined)? And if yes, if there is a way to circumvent it?

Thank you again!

Pierre

I see. You will have to work a little harder to demonstrate this issue.
Use a small box, don't run dynamics, and carefully separate out the values
of the energies that ReaxFF is computing from of the fixed bonds and
angles. Only then can you say for sure what is happening. I do not expect
the ReaxFF calculation to change in response to the presence of fixed
bonds and angles.

That is why I am interested in using ReaxFF in combination with a bond_style and an angle_style, which seems to be disallowed > by LAMMPS as it can be seen from the log file of my previous message.

I do not think it is disallowed to use a bond or angle style with
ReaxFF. It is just that ReaxFF knows
nothing about what the bond/angle styles compute, so you have problems
if you are computing both ReaxFF
and explicit bonds for some atoms.

The right way to do your problem is probably to only have ReaxFF
operation on a subset of the atoms,
which would be easy using pair hybrid. E.g. assign ReaxFF to some
atoms, and no pair style (or a simple LJ pair style) to others,
and then also define bonds and angles for non-ReaxFF atoms.

All other pair styles allow that, but ReaxFF does not yet.
Unfortunately both the Fortran and C++ ReaxFF were written
without pair hybrid in mind.

Steve

Hi,

Following Aidan's comment, I have checked the influence of defining bonds and angles when running minimizations on one water molecule with pair_style reax/c.
I have attahced my input and log files. Files named *.reaxb correspond to minimizations in which bonds and angles were defined, on the contrary no bond or angle is defined in the *.reax files.
The separation of all contributions to the potential energy is included in the log files. Comparing log.reax with log.reaxb, it clearly appears that all terms of the pair_style reax/c are affected by the definition of bonds. At first I would say that LAMMPS does not compute any of the reax/c terms between pairs of atoms that are linked together by the definition of a bond or an angle. Actually the ea contribution (Eunder + Eover) is the only one that does not go to 0. Let me also point out that when bonds are defined, bonded atoms are considered by LAMMPS as special neighbors (as it can be seen from the last section of the log.reaxb file).

Knowing this, is there a way to force LAMMPS computing reax/c terms despite the definition of bonds/angles?

Pierre

data.reax (260 Bytes)

log.reaxb (7.07 KB)

data.reaxb (358 Bytes)

log.reax (6.21 KB)

in.reaxb (1.74 KB)

in.reax (1.24 KB)

Pierre,

Thanks for making this clear. It looks like you are right. The
special-neighbors factor is being used within reax/c. We'll take a look at
this more closely and get back to you.

Aidan

Hi,

Following Aidan's comment, I have checked the influence of defining bonds and angles when running minimizations on one water molecule with pair_style reax/c.
I have attahced my input and log files. Files named *.reaxb correspond to minimizations in which bonds and angles were defined, on the contrary no bond or angle is defined in the *.reax files.
The separation of all contributions to the potential energy is included in the log files. Comparing log.reax with log.reaxb, it clearly appears that all terms of the pair_style reax/c are affected by the definition of bonds. At first I would say that LAMMPS does not compute any of the reax/c terms between pairs of atoms that are linked together by the definition of a bond or an angle. Actually the ea contribution (Eunder + Eover) is the only one that does not go to 0. Let me also point out that when bonds are defined, bonded atoms are considered by LAMMPS as special neighbors (as it can be seen from the last section of the log.reaxb file).

Knowing this, is there a way to force LAMMPS computing reax/c terms despite the definition of bonds/angles?

have you tried this setting?:

special_bonds lj/coul 1.0 1.0 1.0

axel.

Sorry Aidan, I forgot to mention that I use this waterbox simulation as a test, I do not try to extract any meaningful physical value from it at this point.
Of course, I have tried the same simulation without defining bonds and angles, and it works well (once more I do not know if it reproduces well the behavior of liquid water but the simulation runs normally).
I would like to model non-reacting solvent molecules surrounding reactive species. One of the major problem is that solvent molecules tend to overreact (with each other or with the solute), and this holds true with all the potential files that I have tested so far. I think that defining fixed bonds and angles for solvents molecules could make them non-reactive. That is why I am interested in using ReaxFF in combination with a bond_style and an angle_style, which seems to be disallowed by LAMMPS as it can be seen from the log file of my previous message. So, I wonder if this is expected or not (there is no error or warning message but ReaxFF bonding terms are obviously set to 0 when bonds are defined)? And if yes, if there is a way to circumvent it?

when the exclusion factor is 0.0 for a given interaction,
then the pair is completely removed from the neighbor list.

the default for special_bonds is 0.0 0.0 1.0

axel.

That makes senses. Settings the 1st and 2nd special factors to any value
other than zero will probably eliminate the issue.

Indeed, I just tried Axel's suggestion and it seems to work.
Thank you all!

Pierre