Reax/c and compute rdf command

Dear colleagues!

Recently I have been calculating the liquid carbon states with reax/c potential in LAMMPS.
And when calculating the RDF, I found a very strange bug: g(r) did not go to 1 for large r for the liquid carbon?
When I use others potentials and the same states (and same cell) for the liquid carbon, g(r) goes to 1 and all is OK.
I attached my script with source file and RDF image.
I am using LAMMPS-03Feb13. Could anyone comment on this,
or at least give some ideas what I'm doing wrong?

Thank you,

C_ReaxFF_in.dat (48.6 KB)

in.script (1.36 KB)

RDF.pdf (37.8 KB)

Hi Vitaly,

This is not a bug, rather this can mean two things:

1. The liquid carbon is not well equilibrated within the total of 2.5 ps, and
2. This particular force field description does not yield amorphous carbon.

I'd suggest you search the Literature if there is a ReaxFF force field
for amorphous carbon, and if there is one, equilibrate with a larger
timestep (0.01 is too small for a pure system).


Assuming you could generate a LARGE equilibrated system using any
other of the potentials that work, you could try using that as an
initial config into your reaxff simulation. The idea is to replace the
smoothing you get in the rdf from ave/time by the corresponding one
given to you by a larger system size. Hopefully if the system is large
enough you won't have to use the ave/time to get the rdf. Or at least
only very few steps will be required thus avoiding reaxff to "ruin"
the originally equilibrated liquid carbon.
Note: This is just an exercise to convince you there is no bug in the
rdf calculation when using reaxff.

Another thing you could try to verify there is no bug is to take an
equilibrated config from another potential
that works

Thank You very much for your assistance!

1) I have tried different methods of calculations. I used a cell with 512 atoms and with 3500 atoms,
and got the same result. I equilibrated my system within the total of 25-30 ps with different time steps
and got the same.
2) I had took the same source file and replaced the ReaxFF potential with Airebo, and then I got a normal RDF.
3) I calculated the RDF function in the package VMD with the coordinates and got a normal RDF! This result I absolutely do not understand!?! Why? This means, that the ReaxFF modelling gives
the correct coordinates in the simulation!
4) I will search the literature about ReaxFF and amorphous carbon, may be I can find something interesting...

Thanks again!

14.02.2013 21:02, Carlos Campana пишет:

Hi Vitaly and Carlos,

I was able to get a normal RDF by changing pair_style reax/c to a very
simple lj/cut:

pair_style lj/cut 10.0
pair_coeff 1 1 1.0 1.0 10.0

This is interesting since this definitely is a bug somewhere in the
code. I will take a look.


Dear Ray!

Did you find a bug or have not had the time?
I have run the calculations of RDF in the latest version of LAMMPS,
and again got the same bug :-(.


15.02.2013 4:32, Ray Shan пишет:


Sorry I totally forgot about this. I ran some more tests and both reax and reax/c yield the same PDF curve, which are different from a simple lj/cut. No other fixes were defined so nothing moved. I have no clue why this is the case. Also one peculiar thing is if I replicate the system, the PDF is closer to 1 at long distance. It implies something is wrong with periodic conditions, but the energy and stress are identical before and after replicate.