sigma bond order calculation in ReaxFF reaxc_bond_orders.cpp

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS reax/c users,

While inspecting the code for bond order calculation in reaxc_bond_orders.cpp and comparing to the supplemental info from Chenoweth2008 paper I found

BO_s = (1.0 + bo_cut) * exp( C12 );

in the “int BOp” function, but from the paper sup. info. (eq. 2) I believe it should be

BO_s = exp( C12 );

of course bo_cut is often very small, but still doesn’t seem compatible with specs to me, and could make results go astray after a large number of interactions w/ somewhat high bo_cut. It seems to be there since from the first version of this file.

Does anyone has a clue as to why include bo_cut in sigma BO calculation?

Thanks,
Fabio Campolim
Ph.D. student, UFABC, Brazil

2

I will take a look, but please note that bo_cut is a parameter from
the force field file and is not something users can change. It is a
fixed value, and neither control settings nor input commands can
change it to something else.

Ray

3

Hi Ray, thanks for your reply.

ReaxFF is central to my Ph.D. thesis so I’m taking a close look at the code in order to fully understand what each parameter and equation does. Establishing correspondence between the code and van Duin’s papers is not straightforward sometimes. Also, we do want to optimize force field parameters at some point using LAMMPS and DFT calculations/experimental data.

I see now that all cutoffs (sbp, twbp and thbp) are being subtracted at some point from bond orders, like in valence and torsion angles cpp:

BOA_ij = bo_ij->BO - control->thb_cut;

So I imagine they are correct, but I can’t explain them, and they’re not documented.

There are a few other differences between Chenoweth2008 supplemental info and LAMMPS code I have noted, some of the equations are quite different. I’ll be documenting them and perhaps we can discuss if you like.

Thanks,
Fabio