About ReaxFF implementation

LAMMPS LAMMPS Development
1

Dear LAMMPS users/dev,

I have some questions regarding REAXFF implementation.

We are trying in my group to implement the correction described in this paper [https://pubs.acs.org/doi/full/10.1021/acs.jpclett.9b02810]. It aims to correct the energy drift observed in the NVE ensemble by using tapered functions on the bond orders. We believe such modifications would be useful for the LAMMPS community, and for my research project because I observe some unphysical rare events that may be related to this energy drift (I can develop if you care).

I think we’re on the right track, but some things are eluding us, especially with regards to the shifts applied to the bond orders:

  • BO_s = (1.0 + bo_cut) * exp( C12 ); in reaxff_bond_orders.cpp
  • BOA_ij = bo_ij->BO - control->thb_cut; in reaxff_torsion_angles.cpp
  • BOA_jk = bo_jk->BO - control->thb_cut; in reaxff_torsion_angles.cpp

I think we need to understand these details before we can properly apply the functions. Does anybody here could provide some additional explanations regarding these shifts? The question has already been raised here : sigma bond order calculation in ReaxFF reaxc_bond_orders.cpp

Best regards,

Thibaut Flottat

2

Thibaut,

this is a very specific question about the specific implementation of ReaxFF in LAMMPS. I think your best chance to get a competent answer is to contact the original author of the package: ‪H. Metin Aktulga‬ - ‪Google Scholar‬

Axel.

3

This is something that exists in the original ReaxFF Fortran code, so the best source of information would be Adri van Duin himself.

Also, I would be careful with implementing these changes. If you compare potential energy between old and new implementation (Fig.1a and Fig.2a) you will notice that the minimum is significantly changed. If you still aim to incorporate this into LAMMPS, please do it through a keyword, so old behaviour can be maintained and the tapering range should be adjustable as the values chosen by the authors clearly changed bond energies.

2 Likes
4

Axel,

Thank you for your suggestion, I have already contacted him (it is the first thing I did), he did not answer yet. So I guessed some dev here could have an answer.

Actually, I think I understood why BO_s is shifted by bo_cut: to have BO_s continuous at the cutoff. The derivative is not, but becomes so with the correction by Furman. My question now is: is it a good idea to apply the tapered function to the shifted BO_s ? It looks like we add a correction over a correction.

For the shifted BOA_ij (and other valence/torsion terms), I am not sure it is for the same purpose.

I contacted Furman, and get the modified code. The implementation in PuReMd looks older/different than in the code Furman sent me. I have to dig further in the code.

–

Mkanski,

Yes, you are right, I will contact him.
I have to test if it corrects unphysical behaviours I observe in my simulations. If yes, and if I propose an incorporation, it will be through keyword.