Possible Bug in REAX

Hi,
  I am trying to run some simple calculations with using REAX
potential. I ask for forces on atoms in two sightly different
configuration, and get very different answers. The two configurations
differ only in one coordinate of one atom in the 9th decimal place,
while the forces differ by as much as 200%. I am thinking this is not
the intended behavior from REAX, so it is perhaps a bug.
  The files needed to reproduce this are available temporarily at
http://dhruv.ccmr.cornell.edu/ashivni/
(You can get individual files or all of them zipped in the one zip
file). The contents of the files are fairly obvious:

1. data1 and data2 : Atomic configuration for the two cases; you can
check that they are the same except for a 9th place decimal digit
2. ffield.reax : reax force field (supplied with lammps
under the name ffield.reax.rdx)
3. in_lammps1 : input file1, that I run with lmp_serial <
in_lammps1 to get dump1
4. in_lammps2 : input file2, that I run with lmp_serial <
in_lammps2 to get dump2
5. lmpInfo.txt : output of lmp_serial -h on my machine
6. dump1 and dump2 : output on running the input files through lammps.
Clearly the output are vastly different.

Clearly, the dump1 and dump2 files give wildly different forces. I
have tried this on some older versions of lammps and on different
platforms and gotten the same results.

Any help will be appreciated.

Thanks.

PS: I am appending the individual files now, just for convenience.

<start data1>

(written by ASE)

20 atoms
1 atom types
0.0 6.644926445 xlo xhi
0.0 17.672900142 ylo yhi
0.0 10.000000000 zlo zhi

Atoms

     1 1 0 6.413575990 10.218395017 5.000000000
     2 1 0 5.013623267 9.931052352 5.000000000
     3 1 0 5.919965964 12.814578119 5.000000000
     4 1 0 5.452305541 7.454326900 5.000000000
     5 1 0 5.945858910 4.858469109 5.000000000
     6 1 0 5.004842877 6.033464424 5.000000000
     7 1 0 4.075302159 11.023442061 5.000000000
     8 1 0 4.518097454 12.431444959 5.000000000
     9 1 0 2.639284798 10.746849730 5.000000000
    10 1 0 2.127039183 9.358655346 5.000000000
    11 1 0 4.521375697 8.550176124 5.000000000
    12 1 0 3.093977637 8.314244487 5.000000000
    13 1 0 3.564349068 5.827347964 5.000000000
    14 1 0 2.581729778 6.926051352 5.000000000
    15 1 0 1.656848043 11.845577451 5.000000000
    16 1 0 0.216141669 11.639457957 5.000000000
    17 1 0 0.699638666 9.122720863 5.000000000
    18 1 0 0.702936745 5.241449706 5.000000000
    19 1 0 1.145710876 6.649455527 5.000000000
    20 1 0 0.207384918 7.741841980 5.000000000

<end data1>

Aidan or Ray can probably take a look at this.

Steve

Hi Ashivni,

Thanks for the detailed problem description and attaching all
necessary files. I was able to reproduce the inconsistency you saw
with reax/c. The difference in your original data file are the
coordinates of atom 15, so I grepped the info of atom 15 from the dump
outputs:

reax/c with ffield.reax.rdx

Dear Ray,
    It is possible that this is specific to the potential file
ffield.reax,rdx, however, just changing the cutoff parameter does not
remedy the situation. The following files have the cutoff parameter
changed, but produce similar error. data3 and data4 have coordinates
off by a minuscule amount, yet, when you use ffield.reax1 as the
parameter file, the forces are wildly different. I have set the cutoff
parameter to 0.00011 . You can also find the files at
http://dhruv.ccmr.cornell.edu/ashivni/

I am appending the files here for convenience.
Thanks!
Ashivni.

<start data3 >
data3 (written by ASE)

20 atoms
1 atom types
0.0 6.644770167 xlo xhi
0.0 17.672719689 ylo yhi
0.0 10.000000000 zlo zhi

Atoms

     1 1 0 6.413711328 10.218398917 5.000000000
     2 1 0 5.013745765 9.931105964 5.000000000
     3 1 0 5.919620999 12.814651197 5.000000000
     4 1 0 5.451960965 7.454343432 5.000000000
     5 1 0 5.946053070 4.858045828 5.000000000
     6 1 0 5.004860086 6.033134587 5.000000000
     7 1 0 4.075276776 11.023379343 5.000000000
     8 1 0 4.517761830 12.431528544 5.000000000
     9 1 0 2.639315492 10.746682369 5.000000000
    10 1 0 2.126964655 9.358552518 5.000000000
    11 1 0 4.521296442 8.550276041 5.000000000
    12 1 0 3.093921917 8.314166665 5.000000000
    13 1 0 3.564146426 5.827431703 5.000000000
    14 1 0 2.581576980 6.926034819 5.000000000
    15 1 0 1.656721073 11.845286845 5.000000000
    16 1 0 0.216030292 11.639580254 5.000000000
    17 1 0 0.699594017 9.122444395 5.000000000
    18 1 0 0.703127183 5.241189841 5.000000000
    19 1 0 1.145615333 6.649345181 5.000000000
    20 1 0 0.207148158 7.741618445 5.000000000

<end data3>

Hi Ashivni,

No, you misunderstood what I meant. I meant changing the "thb_cutoff"
in the control file, not changing parameters in the force field file.
Please see pair_style reax/c doc page on how to change the default
variable "thb_cutoff".

Note 1: Unless one is developing a new force field for a new type of
application, one should not change anything in any force field file.

Note 2: Changing control settings is only possible with reax/c, not with reax.

Cheers,
Ray

Ashivni,
You seem to be the kind of guy who likes to understand the "fine
print". That's a good thing. Why don't you ask Ari Van Duin for his
original reaxff fortran code and check if the same issues appear
there?
Carlos

I'd like to chime in and say that generally ReaxFF is parameterized
against energies of different states and bonding arrangements, and
bond stiffnesses are not usually included (at least not explicitly).
I've found that this often leads to sharp spikes in the force gradient
for carbon, and likely other elements.

As Ray said it's just a product of the force field.

-Ben

Benjamin Jensen
PhD Student
Mechanical Engineering - Engineering Mechanics Department
Michigan Technological University
Houghton, MI

Dear all,
   Thanks for the input and apologies for the silence on my end. I was
caught up in a lot of travel. Anyhow, it seems that playing with the
thb_cutoff parameter in the control file for reax/c does not fix the
problem. Attached is a case where I have set thb_cutoff = 0.002, and I
see the same behavior. Does this mean that reax will always behave
funny. Ben said in his last email that there will be sharp spikes in
the force field gradient; it seems that the gradient will have to be
about 10^9 (in force/distance for real units). Does that sound like it
should a reasonable behavior for any potential?

Here are the new files (ffield.reax is unchanged, they can be found
as usual on http://dhruv.ccmr.cornell.edu/ashivni/ )

<start data5>
data5 (written by ASE)

44 atoms
1 atom types
0.0 16.468921964 xlo xhi
0.0 17.518219581 ylo yhi
0.0 10.000000000 zlo zhi

Atoms

     1 1 0 16.271223711 8.334632314 5.000000000
     2 1 0 16.061572113 6.831026305 5.000000000
     3 1 0 15.134627435 9.183589276 5.000000000
     4 1 0 15.344275678 10.687195097 5.000000000
     5 1 0 13.795073065 8.726851345 5.000000000
     6 1 0 14.202420197 11.614799876 5.000000000
     7 1 0 14.698920948 6.272366983 5.000000000
     8 1 0 13.569611787 7.219684469 5.000000000
     9 1 0 12.654914418 9.561559773 5.000000000
    10 1 0 12.836613280 11.062568262 5.000000000
    11 1 0 11.307220081 9.086788948 5.000000000
    12 1 0 11.671085858 11.966342156 5.000000000
    13 1 0 12.192058947 6.691031843 5.000000000
    14 1 0 11.059054222 7.614840291 5.000000000
    15 1 0 10.140296339 9.893160509 5.000000000
    16 1 0 10.320715867 11.385416899 5.000000000
    17 1 0 8.808410539 9.358755568 5.000000000
    18 1 0 9.118721133 12.262300827 5.000000000
    19 1 0 9.603340608 5.560277664 5.000000000
    20 1 0 9.720855988 7.031963468 5.000000000
    21 1 0 8.324436923 4.870896262 5.000000000
    22 1 0 8.580173923 7.902922178 5.000000000
    23 1 0 7.654855303 10.219358167 5.000000000
    24 1 0 7.782967939 11.699685564 5.000000000
    25 1 0 6.356755325 9.615295431 5.000000000
    26 1 0 6.612497184 12.647323873 5.000000000
    27 1 0 7.153962367 5.818528856 5.000000000
    28 1 0 7.282074808 7.298856537 5.000000000
    29 1 0 5.818208032 5.255919748 5.000000000
    30 1 0 6.128521899 8.159462557 5.000000000
    31 1 0 5.216073755 10.486256230 5.000000000
    32 1 0 5.333594246 11.957941970 5.000000000
    33 1 0 3.877876254 9.903374933 5.000000000
    34 1 0 4.616214026 6.132806963 5.000000000
    35 1 0 4.796634242 7.625060374 5.000000000
    36 1 0 3.265843593 5.551880884 5.000000000
    37 1 0 3.629705972 8.431427370 5.000000000
    38 1 0 2.744871824 10.827183898 5.000000000
    39 1 0 1.367316541 10.298539699 5.000000000
    40 1 0 2.100314654 6.455649913 5.000000000
    41 1 0 2.282010833 7.956659347 5.000000000
    42 1 0 0.734505121 5.903421619 5.000000000
    43 1 0 1.141854610 8.791372183 5.000000000
    44 1 0 0.238005238 11.245854354 5.000000000

<end data5>

Ashivni,

Your estimate of 10^9 for the force gradient assumes the force is
continuous. Have you checked to see if the two geometries are straddling a
step change in the force? Most potentials have points of discontinuity in
the force, and sometimes also the energy. ReaxFF has both.

Aidan

Ashivni,

Have you discovered that when you change the “thb_cutoff” back to 0.001 for data5 and data6, the results match again? Before when you were using data1 and data2, you had to change thb_cutoff to anything else other than 0.001 to match the results. This supports Aidan’s suggestion that you are experiencing points of discontinuity in the force, which can be due to force field, control setting parameters, and maybe the source code. Per Carlos’ suggestion, perhaps the best way to verify this is to request a stand-alone ReaxFF version from Prof. van Duin, and use the exact same force field and control settings for comparison.

Ray

Allow me to add that although not generally mentioned in Reaxff papers, sometimes the value of thb_cutof used by Van Duin and colleagues to obtain certain parametrization is not the default of 0.001