reax and reaxc not matching example output

Dear LAMMPS users,

I would like to use the ReaxFF force field to study carbon nanotubes.
However I a problem when I try to reproduce the output that is given in the examples.

I have am using the latest SVN version (rev. 8549) in serial.
I am running on a Debian GNU/Linux 64-bit system, this my compilation line directly after checkout
cd lib/reax; echo “reax_SYSLIB = -lgfortran” > Makefile.lammps; make -f Makefile.gfortran; cd …/…/src/STUBS; make; cd …; make yes-USER-REAXC; make yes-reax; make -j4 serial
To run the example :
cd …/examples/reax; …/…/src/lmp_serial < in.reaxc.rdx

I find similar problems with both USER-REAXC and REAX, but even among them the output is not the same.

For “USER-REAXC” :
The reference file is
http://petveturas.com/tmp/lammps_reax/log.reaxc.rdx.9Jan12.linux.1
Using the serial svn version I find
http://petveturas.com/tmp/lammps_reax/log.reaxc.rdx
Since there seems to be a different definition of the potential between the logfile I modified the input file to match the reference output.
http://petveturas.com/tmp/lammps_reax/log.reaxc.rdx.modified

And for “REAX” :
The reference file is
http://petveturas.com/tmp/lammps_reax/log.reax.rdx.9Jan12.linux.1
I find
http://petveturas.com/tmp/lammps_reax/log.reax.rdx
Again, changed the input
http://petveturas.com/tmp/lammps_reax/log.reax.rdx.modified

All my outputs are different but none are like the reference file.
Could somebody please help see my mistake?
At least with the steps above it should be easily reproducible.

The reference files are from 02.jan.2012, has some default setting in LAMMPS changed since?

Many thanks!

Best,
Jaap

Dear LAMMPS users,

I would like to use the ReaxFF force field to study carbon nanotubes.
However I a problem when I try to reproduce the output that is given in the
examples.

no surprise there. there have been several tweaks and bugfixes
to the reax implementations in LAMMPS since those example
runs were done and the outputs have not been updated. you
can see some mentioning of the changes here:
http://lammps.sandia.gov/bug.html

I have am using the latest SVN version (rev. 8549) in serial.
I am running on a Debian GNU/Linux 64-bit system, this my compilation line
directly after checkout
cd lib/reax; echo "reax_SYSLIB = -lgfortran" > Makefile.lammps; make -f
Makefile.gfortran; cd ../../src/STUBS; make; cd ..; make yes-USER-REAXC;
make yes-reax; make -j4 serial
To run the example :
cd ../examples/reax; ../../src/lmp_serial < in.reaxc.rdx

I find similar problems with both USER-REAXC and REAX, but even among them
the output is not the same.

For "USER-REAXC" :
The reference file is
http://petveturas.com/tmp/lammps_reax/log.reaxc.rdx.9Jan12.linux.1
Using the serial svn version I find
http://petveturas.com/tmp/lammps_reax/log.reaxc.rdx
Since there seems to be a different definition of the potential between the
logfile I modified the input file to match the reference output.
http://petveturas.com/tmp/lammps_reax/log.reaxc.rdx.modified

And for "REAX" :
The reference file is
http://petveturas.com/tmp/lammps_reax/log.reax.rdx.9Jan12.linux.1
I find
http://petveturas.com/tmp/lammps_reax/log.reax.rdx
Again, changed the input
http://petveturas.com/tmp/lammps_reax/log.reax.rdx.modified

All my outputs are different but none are like the reference file.
Could somebody please help see my mistake?
At least with the steps above it should be easily reproducible.

The reference files are from 02.jan.2012, has some default setting in LAMMPS
changed since?

yes.

axel.

Dear Alex Kohlmeyer,

Ok thanks your reply.
Only this doesn’t explain to me why the two implementations should differ.
I found that the differences are in columns 3,4,5 and 10, i.e. E_pair, TotEng, Press and ev.
The temperatures are also different but is probably a result of the different potentials.

The results first differ from the reference output upon going from revision 7975 to 7976.
This is simply because the example changes.

$ svn diff --summarize -r 7975:7976
M examples/reax/ffield.reax
M examples/reax/in.reax.tatb
M examples/reax/README
M examples/reax/in.reax.rdx
M examples/reax/in.reaxc.tatb
M examples/reax/in.reaxc.rdx
M examples/reax/control.reax_c.tatb
A examples/reax/control.reax_c.rdx

After this the two implementations reax and USER-REAXC still give the same results
However, the two implementations start to differ from revision 8149.

$ svn diff --summarize -r 8148:8149
M src/USER-REAXC/pair_reax_c.cpp
M src/USER-REAXC/reaxc_valence_angles.cpp
M src/USER-REAXC/reaxc_control.cpp
M src/USER-REAXC/reaxc_types.h

If I understand the changes (below) correctly this is a bugfix because the 3-body interaction cutoff was not squared before comparison.

$ diff lmp_814{8,9}/src/USER-REAXC/pair_reax_c.cpp
191a192

control->thb_cutsq = 0.00001;

$ diff lmp_814{8,9}/src/USER-REAXC/reaxc_valence_angles.cpp
248c248
< (bo_ij->BO * bo_jk->BO > 0.001) ) {

(bo_ij->BO * bo_jk->BO > control->thb_cutsq) ) {

$ diff lmp_814{8,9}/src/USER-REAXC/reaxc_control.cpp
75a76

control->thb_cutsq = 0.00001;
211a213,216
}
else if( strcmp(tmp[0], “thb_cutoff_sq”) == 0 ) {
val = atof(tmp[1]);
control->thb_cutsq = val;

$ diff lmp_814{8,9}/src/USER-REAXC/reaxc_types.h
529a530

real thb_cutsq;

Does mean that the ‘reax’ implementation contains a bug?
Which implementation should I trust/use?

Thanks.

Best,
Jaap

Does mean that the 'reax' implementation contains a bug?

possibly.

Which implementation should I trust/use?

i am not an expert on this, but reax/c is generally
recommended. it is more actively maintained
and has the better performance, too.

hopefully aidan, who has been the caretaker of
the fortran version, can comment on the reax status
vs. reax/c.

axel.

Dear Alex Kohlmeyer,

I think I have at least some more evidence that the USER-REAXC code is more reliable.
The original reason I started to do more extensive testing was that I found a problem in energy minimization.
Now I think this is related to a discontinuity in the energy. This discontinuity does not exist in the USER-REAXC code in my specific example.

I attached a working example, with input and processing / plot script as well as the resulting figure.
The system is a carbon nanotube with one chemisorbed oxygen atom.
The attached figure shows the energy as a function of the x-coordinate of the oxygen atom around the ‘optimized’ geometry.
The energies are slightly different but more importantly the ‘reax’ code has discontinuity near the minimum, which would explain the difficulty I find when optimizing the geometry.

Best,
Jaap

discont_test.tgz (57.5 KB)

reax_vs_reaxc.pdf (27.5 KB)

Dear Alex Kohlmeyer,

Sorry to keep bothering you but apparently my previous test was not at all representative.
I did the same test starting from the ‘reax/c’ minimum and now both show some discontinuity.

So I suppose the most likely thing is that there is simply something wrong with my input.

I specify reax as

pair_style reax
pair_coeff * * ffield.reax 1 3

and reax/c as

pair_style reax/c NULL
pair_coeff * * ffield.reax 1 3
fix 2 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c

I use “atom_style charge”, but specify a neutral system (zeros in the third column of the data file)

created from sample.xyz
361 atoms
2 atom types

0.0 +38.7908357334 xlo xhi
0.0 100.0000000000000000 ylo yhi
0.0 100.0000000000000000 zlo zhi

Atoms

1 1 0. 0.546583 0.000442 4.011575
2 1 0. 1.980005 -0.003243 4.011961
3 1 0. 2.702379 1.244105 3.817130
4 1 0. -0.171568 1.250286 3.813417

Is there any obvious reason why this should fail?

Thanks again.

Best,
Jaap

p.s. below is the complete input file I used to compute the total energy with reax

Here is also the attachment with the test.

discont_test2.tgz (47.9 KB)

ex.pdf (21 KB)

Hi Jaap,

The discontinuity resulted from the ffield file that was used, rather
than from the implementation of reax/c. Attached pdf contains two
curves obtained with your script; the red curve with
ffield.reax.nielson (the one you used which was developed by Nielson
and is not in the standard LAMMPS distribution) and the blue with
Chenoweth's CHO. Obviously ffield.reax.cho does not have the
discontinuity problem, and there may be something wrong with the
Nielson ffield.

Best,
Ray

CNT.pdf (18 KB)

Hi Jaap,

The reasons for results not matching the examples with reax/c are due
to the following:

1. The default force field "ffield.reax" in the examples dir was
changed from "ffield.reax.mattsson" to "ffield.reax.rdx" sometime
between Feb and May. (These two ffields are located in the
/potentials directory.) Changing from ffield.reax.rdx to
ffield.reax.mattson gives you very similar results, but not quite
exact.

2. The definition of thb_cutoff_sq in the source code was changed
from "equal to thb_cutoff" to "equal to thb_cutoff * thb_cutoff" in
May to correspond to Prof. Adri van Duin's serial ReaxFF code.
Setting "thb_cutoff_sq" to 0.001 in the control file will exactly
reproduce the RDX examples with Reax/C.

Best,
Ray

Dear Ray Shan,

Thank you for your responses.

However I still find a similar discontinuity, also with the ‘ffield.reax.cho’ force field (ex.pdf).
I again attached a similar script (discont_test3.tgz) but for another geometry.
I find these discontinuities whenever the geometry optimization gets “stuck”.

After separating the 14 different terms (contributions.pdf) in the ‘reax’ forcefield it seems that my problem comes from the conjugation term. I attached a script to reproduce this graph (reax_sep.tgz).

Thanks for any suggestions !

Best,
Jaap

discont_test3.tgz (42.6 KB)

ex.pdf (21 KB)

contributions.pdf (67.2 KB)

reax_sep.tgz (115 KB)

ReaxFF, like many potentials, contains small discontinuities. Mostly, these are associated with omission of bonds less than a certain bond order. This one is not tiny, but I suspect it is a normal part of the potential.

Hi Jaap,

I am not sure why that is. But you are probing the
energy-displacement curve within +/- 0.01 Angstroms, which may just be
too small. The ReaxFF force field or the conjugation term may not be
perfectly smooth within such a small displacement.

Ray