about reax and reaxc bond energy discrepancy

LAMMPS LAMMPS Mailing List Mirror
1

Dear all,

I made recently a comparison between original Prof. van Duin’s reax library and reaxc library implemented in lammps and found a difference in bond energy term. For instance, for my Ni/ZrO2 interface consisting of 34 atoms bond energy calculated with reax is -7587.6595 kcal/mol, while with reax/c it is -7754.7632 kcal/mol. Cell size was 3.64A x 5.15A x 43A.

Having checked two source codes I found that the difference might come from the way the bond energy is calculated in Bonds() function in reaxc_bonds.cpp file. Specifically, in the double loop running over all pairs of bonded atoms (say atoms i and j) when “system->my_atoms[i].orig_id == system->my_atoms[j].orig_id” (sorry for this) energy of this i-j bond must be probably divided by 2. After this small modification I’ve obtained with reaxc bond energy -7587.6610 kcal/mol, which is identical to the original reax value.

When I increased simulation cell size to 40A x 41A x 43A I got same energies with unmodified reaxc and with my modifications in it. So probably the discrepancy is caused by the bonds between an atom and its images across periodic boundaries.

Could you please let me know which approach (reax or reaxc) is correct in the case of small simulation cell?

2

Aidan and Ray can comment.

Steve

3

Hi Albert,

This does indeed sound like an incorrect calculation of bond energy when an atom is bonded to its periodic image. Can you send us a very simple example? Ideally, it would include a “replicate n n n” command and Ebond/n^3will vary when n is increased from 1 to 2.

Aidan

4

Dear Aidan,

thank you for your comments. I’ve prepared minimum setup demonstrating the issue based on a pure fcc aluminium structure. Size of the simulation box is 4.05A along all three axes. The simulation box contains only 4 atoms. Below I list bond energies obtained with different n values using two different reaxff libraries.

n Bond energy (reax), kcal/mol Bond energy (reax/c), kcal/mol
1x1x1 -551.36073 -551.96218
2x2x2 -4410.8858 -4410.8858

The situation is similar to what you predicted. With reax library bond energy scales perfectly with n as E_2 = 2^3 * E_1, while with reax/c for n=1 bond energy is lower than what we could predict from the n=2 case (E_2 > 2^3 * E_1). With larger values of n bond energies are same, calculated with reax and reaxc.

I attach one file with atomic configuration (struc.lmp), one reaxff potential file (ffield.reax), which is actually provided in lammps, and two input scripts - in.Al_reax and in.Al_reaxc, which use reax and reaxc libraries, respectively. The latter two read the former ones. I hope you can reproduce this behaviour on your side with these files.

struc.lmp (366 Bytes)

ffield.reax.Al (21.7 KB)

in.Al_reaxc (994 Bytes)

in.Al_reax (993 Bytes)

5

Thanks, Albert. I will take a look and let you know.

Thanks,
Ray

6

Albert,

The attached source file eliminates the double-counting of bonds for self-bonded atoms. I have checked this into the repo, and something similar for USER-REAXC/reaxc_torsion_angles.cpp. The same issue probably occurs in reaxc_hydrogen_bonds.cpp, but without an actual example, I don’t want to attempt to fix it.

Aidan

reaxc_bonds.cpp (5.03 KB)

7

Thank you very much for the modification!

The problem doesn’t appear any more.