Hello,
I've compared energies & forces on several peptides calculated using both
ReaxFF implementations in LAMMPS (reax and reax/c) with results from the
commercial SCM ReaxFF implementation and found significant differences in
valence angle energies and total forces on some atoms. Forces on most of the
atoms are equal in LAMMPS and SCM ReaxFF, but especially the hydrogen atoms in
-NH3 groups exhibit striking differences (not only in magnitude, but also in
direction of the force).
The smallest system I've been able to reproduce this on is CH3NH3 (protonated
methylamine; but the charge doesn't matter and I've failed to set total charge
in SCM ReaxFF so the test was done using a neutral molecule). I've put all the
files here: http://is.muni.cz/www/395845/CH3NH3.tar.gz
I've been able to reproduce the differences using both force fields I checked:
CHON_Gly.ff: doi:10.1021/jp108642r
HCONSB.ff: included with SCM ReaxFF, C/H/O: Chenoweth et al 2008; S/C/H: Castro
et al. 2011; N: B/N: based on Weismiller et al, 2010
Both implementations in LAMMPS provide the same results. I've discussed this
with the authors of SCM ReaxFF and they confirmed their software produces
results equal to original reax code they received from Adri van Duin in July
2009. The version.f file in the original code contains:
100 format ('ReaxFF version 2.0')
110 format(a20,'Fri Jun 26 10:42:10 EDT 2009')
Setting the p_val4 parameter for the H-H-H valence angle term in the force
field to zero makes the difference vanish. This more or less turns off the
dependence of the angle energy on the bond orders of atoms involved. I've
compared the bond orders of all bonds in the system, but found no difference
between LAMMPS and SCM ReaxFF. It seems something in the valence angle term
silently changed in the original reax version SCM ReaxFF is based on and
LAMMPS was not adapted accordingly. It can be just an implementation bug, too.
(I know that the valence angle term changed between the initial 2001
publication and the latest 2008 functional form - both LAMMPS and SCM ReaxFF
use the new version.).
Does anyone have any hints how to diagnose this further and get it fixed?
Best regards,
Tomáš Trnka
Laboratory of Computational Chemistry
National Centre for Biomolecular Research, Masaryk University
Brno, Czech Republic
Hi Tomas,
I have downloaded your tar ball and am taking a look. Do you have the
"fort.73" file from reaxff code that outputs itemized potential energies?
Thanks,
Ray
There is a bond-order cutoff value for valence angles that is hard-coded
in Adri's code, and also in pair reax and reax/c. The value in the LAMMPS
implementations is
1e-3, while Adri seems to prefer 1e-5. I am reluctant to change the LAMMPS
value, as it will break our regression tests, but if this is considered
truly wrong, I am okay with changing it.
Aidan
p.s. here is some more detailed information on this issue:
Here is the line of code in pair_reax that filters out weak angles:
lib/reax/reax_connect.F
if (bo(jindexbondlist)*bo(kindexbondlist).lt.0.001) goto 50
Here is the corresponding filter in reaxc:
lib/reaxc/reaxc_valence_angles.cpp:
if( (j < system->n) && (BOA_jk > 0.0) &&
(bo_ij->BO > control->thb_cut) &&
(bo_jk->BO > control->thb_cut) &&
(bo_ij->BO * bo_jk->BO > 0.001) ) {
Hi Tomas,
As Aidan mentioned, changing the hard coded value of "0.001" to "0.00001"
changes the angle energy from 20.099 to 20.138, which agrees with the
stand-alone reaxff valence angle energy of 20.14.
Ray
Dne Čt 17. května 2012 15:40:27, Shan, Tzu-Ray NMN napsal(a):
Hi Tomas,
I have downloaded your tar ball and am taking a look. Do you have the
"fort.73" file from reaxff code that outputs itemized potential energies?
Hello,
everything I have is in that tarball. Perhaps I need to enable writing the file
somehow, but I don't have a clue how. Do you know how to do it? Otherwise I'll
have to ask the SCM support.
2T
There is a bond-order cutoff value for valence angles that is hard-coded
in Adri's code, and also in pair reax and reax/c. The value in the LAMMPS
implementations is
1e-3, while Adri seems to prefer 1e-5. I am reluctant to change the LAMMPS
value, as it will break our regression tests, but if this is considered
truly wrong, I am okay with changing it.
Hello,
thanks for the explanation, changing the cutoff value to 1e-5 makes all the
force differences on all my test systems vanish. As far as I am concerned, the
value in LAMMPS can be kept as is (I'm running a patched version of reax/c
anyway so changing this is trivial); however, since it leads to really
significant differences in forces, it should probably be set at whatever value
the force fields were parameterised for. But possibly I've hit just a single
pathological case with the NH3 hydrogens and the differences are negligible on
99% of systems…
Thanks again for such a quick help!
Best regards,
Tomáš Trnka