[lammps-users] Help: pair_style reax/c with ReaxFF-lg force field cannot reproduce the published benchmark result

Dear Lammps Users,

I’m trying to employ the ReaxFF-lg force field [1] by using the “pair_style reax/c” command in LAMMPS-3Mar20. A benchmark case of carbon dioxide crystal at 0K was tested. But the equilibrium lattice constant cannot match with the original publication [1]. The input and log files are attached below.

I prepared the 3*3*3 cells of carbon dioxide lattice (with space group Pa-3) in the file “data.CO2_ICSD_31077”, in which the initial lattice length is 6 Angstrom. Then this system is configured using the ReaxFF-lg force field and relaxed to equilibrium state at 0K. The lattice length shrinks to 5.84 Angstrom, which is not consistent with the 5.7 Angstrom reported in the table 2 in reference [1].

[1] Liu, L., et al. (2011). "ReaxFF- l g: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials." 115(40): 11016-11022.

I have excluded two possible reasons for this discrepancy:
1. The “ffield.reax.lg” found in the “lammps-3Mar20\potentials”folder is consistent with the Supplement Information file in reference [1], except that several irrelevant elements were removed in “ffield.reax.lg”.
2. The default control parameters seem not influence the result. I tried “NULL” and “lmp_control” to replace the cfile in the command “pair_style reax/c cfile lgvdw yes”. The lmp_control file is found in the “lammps-3Mar20\examples\reax\RDX” folder. The result does not change much.

This discrepancy may indicate that current usage of the ReaxFF-lg force field with the “pair_style reax/c” command is not correct. Have anyone here faced similar problem, or point out to a potential reason for this? Any comments are appreciated.


LIANG Tengfei

data.CO2_ICSD_31077 (11.5 KB)

ffield.reax.lg (11.6 KB)

in_Crystal_0K.CO2 (580 Bytes)

log.lammps (160 KB)

why don’t you use fix box/relax and minimize instead of trying to do an MD at 0K?
that would be the method of choice for a “relaxation” at 0K.
with thermostatting the system to 0K using fix langevin you are freezing out the potential for movement. you must have kinetic energy in your system for atoms to move when doing MD.
also, are you certain, that the box changes isotropically?

Dear Ray,

Thanks for your response!

According to the reference [1], the CO2 crystal at 0K is a training case of the ReaxFF-lg potenital field.

I think no matter it’s suitable or not, the CO2 crystal was simulated using the ReaxFF-lg potenital in reference [1]. But now I cannot reproduce the same result using the pair_style reax/c command in LAMMPS.

In addition, I have excluded the influence of the statistical error. The temperature near 0K, there is little thermal fluctuation.

[1] Liu, L., et al. (2011). “ReaxFF- l g: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials.” 115(40): 11016-11022.



Is the ReaxFF-lg potential suitable for CO2?


Dear Alex,

Thanks for your suggestion! I have run a new simulation using the fix box/relax and minimize. The dimensions in three axis are uncoupled. The input and log files are attached above. It’s much faster than my original MD relaxation procedure! But the resulted lattice constant is still around 17.5/3=5.83 Angstrom, which is larger than the 5.7 Angstrom reported in the original publication. For this simple system, such a discrepancy is considerable, which indicates that the potenital field we implement is different from the original publication!



log.lammps (3.58 KB)

in_Crystal_0K.CO2_fix_box_relax (529 Bytes)

Dear Tengfei,

I agree that the difference between your calculation of the lattice constant and that of Liu et al. strongly suggests that there are one or more differences between the energy model you have specified in LAMMPS and what Liu et al. created. I commend you for recognizing that. However, since this ReaxFF energy model is not one of the 8 painstakingly validated examples that we provide with LAMMPS, we are unable to help you further. Instead, I recommend that you contact the authors and figure out which of the many, many possible causes of discrepancy are causing these differences. I know from personal experience that without a code-to-code comparison, it is impossible to isolate the reason(s) for differences in complex energy models such as ReaxFF.


Dear Aidan,

Thanks for your suggestions. I will try to contact the authors to figure out the reasons. O(∩_∩)O