Hello,
I have been testing the various energy contributions computed by the C and Fortran versions of reax that are implemented within LAMMPS (REAX and USER-REAXC packages). I have two observations to share:
- I see disagreement between the two versions in the VdW energy. It seems that the C version includes the VdW inner wall term, while the Fortran version does not. If we are following Chenoweth paper that is cited in the LAMMPS documentation, then the VdW inner wall term should not be included in the equations (i.e., the Fortran version is consistent, the C version is not). If we are following another reference for this term, could you point me to the correct reference?
Chenoweth, van Duin and Goddard, Journal of Physical Chemistry A, 112, 1040-1053 (2008)
- With regards to the hydrogen bond energy, it seems that we have again deviated from the Chenoweth paper where the hydrogen bond term uses a sin4 term, rather than a sin8 term used in Chenoweth. The C version uses the sin4 term. The Fortran has both the sin4 and sin**8 terms coded, where the functional form is controlled by the lhbnew flag in reax_poten.F as written below:
2227 if (lhbnew .eq. 0) then
2228 ehbh=(1.0-exphu1)dehb(ityhb)exphu2sin2sin2sin2sin2
2229 else
2230 ehbh=(1.0-exphu1)dehb(ityhb)exphu2sin2sin2
2231 endif
The lhbnew flag is set within the lammps/src/REAX/pair_reax.cpp file (ihbnew) and set to 1, indicating that the sin**4 functional form should be used. Here is the code from pair_reax.cpp:
48 cutmax = 0.0;
49 hbcut = 6.0;
50 ihbnew = 1;
51 iprune = 4;
52 ihb = 1;
53 chpot = 0;
However, when passing the variable back to the reax library, somehow it is reset back to lhbnew=0, thereby defaulting to the sin**8 functional form. Perhaps I have missed where this is set back to 0 in the LAMMPS source, but it seems that this is not the desired outcome. Could you please provide comments as to whether you observe similar behavior?
Thanks,
–Jim