Gradient errors with qeq/reax

Hello,

I am trying to use the reax/c pair_type with the original ReaxFF CHO combustion parameters to benchmark a code we are developing for saddle point optimization on molecular systems. (Aside: I’m not concerned with getting accurate energetics or even geometries, I just need inexpensive molecular PESs that have the same general features as more expensive PESs from electronic structure theory calculations). I am calling LAMMPS using the Python package ASE, which uses LAMMPS’s official Python interface.

I am running into an issue where the gradients calculated by LAMMPS seem inconsistent with numerically-calculated gradients, but only when I use the qeq/reax fix.

Here is an example script that illustrates the problem:

#!/usr/bin/env python

from ase.io import read

from ase.calculators.lammpslib import LAMMPSlib

atoms = read(‘test.xyz’)

atoms.center(vacuum=0.)

calc = LAMMPSlib(lammps_header=[‘units real’,

‘atom_style charge’,

‘atom_modify map array sort 0 0’,

‘boundary f f f’],

lmpcmds=[‘pair_style reax/c NULL checkqeq no’,

‘pair_coeff * * ffield.reax.cho O C H’,

‘fix 1 all qeq/reax 1 0.0 10.0 1e-6 reax/c’,

],

log_file=‘test.log’,

keep_alive=True)

atoms.set_calculator(calc)

f_analytical = atoms.get_forces()

f_numerical = calc.calculate_numerical_forces(atoms, 1e-3)

print(‘Ratio:\n’, f_analytical / f_numerical)

I’m using the ffield.reax.cho file present in lammps/examples/reax/CHO for these tests.

Here is the geometry file, test.xyz:

12

Properties=species:S:1:pos:R:3:Z:I:1 pbc=“F F F”

O 1.26710000 -0.43400000 0.52740000 8

O -1.78860000 -0.84880000 1.60600000 8

C -2.07510000 0.25800000 -0.45370000 6

C 2.44520000 0.03660000 -0.05010000 6

C -1.25040000 -0.35960000 0.62620000 6

H -1.41130000 0.64530000 -1.25440000 1

H -2.67080000 1.09740000 -0.03780000 1

H -2.75920000 -0.50240000 -0.88510000 1

H 2.83330000 0.88930000 0.54400000 1

H 3.19960000 -0.77660000 -0.06620000 1

H 2.24030000 0.37160000 -1.08760000 1

H -0.03010000 -0.37680000 0.53130000 1

And here is the output:

Ratio:

[[0.99971936 1.00013575 1.00011809]

[1.0000494 1.00002897 0.99998253]

[0.99964268 0.99946184 0.99944392]

[0.99994863 1.00010895 0.99999966]

[1.00005802 1.00004018 1.00004601]

[0.99920468 0.96715718 0.9917202 ]

[1.01810746 0.99643867 0.999678 ]

[0.86276394 0.99884058 0.95186258]

[1.00010121 0.99955586 1.0000197 ]

[1.00015064 0.99993186 1.00469049]

[0.99972746 0.99969436 0.99883464]

[1.01393212 0.99842161 0.99967939]]

Note that some of the components of the gradient disagree by up to 14%. Removing the qeq/reax fix seems to mostly eliminate this issue. Changing the finite displacement step size doesn’t help (in this script it is 1e-3 Angstrom), and changing the qeq/reax tolerance (in this script 1e-6) doesn’t help. This is actually a relatively modest example; I see much larger errors on other systems, but those systems have many more atoms, and I wanted to avoid filling this email with too many poorly-formatted numbers.

These errors are substantial enough to cause serious convergence issues with my saddle point optimizing code, which relies on accurate gradients to calculate finite difference Hessian-vector products.

Thanks,

Eric Hermes

CCing Ray Shan to see if he has ideas on the fix qeq/reax contribution.

Steve

not sure, if this is relevant, but we've recently added changes to
LAMMPS, which allow to use fix qeq/shielded instead of fix qeq/reax
with pair style reax/c in order to work around cases, where fix
qeq/reax would crash but fix qeq/shielded would work.

fix qeq/shielded is a re-implementation and and part of the
generalized QEQ package and thus may be "cleaner" than qeq/reax
in my tests, there were only minimal differences due to floating point
math precision limitations, but it may be worth checking, whether
there is a different behavior in this particular case.

axel.

Axel,

I tried running the test with `qeq/shielded`, and I got pretty much the same errors in the gradients. In fact, with `qeq/shielded`, I was able to improve convergence by lowering tolerance to 1e-12 (which I realize is absurdly small) and increasing the maximum number of iterations to 1000, which let me decrease the finite difference step size to 1e-5 Angstrom. This still results in the exact same errors.

Here is my new input:
from ase.io import read
from ase.calculators.lammpslib import LAMMPSlib

atoms = read('test.xyz')
atoms.center(vacuum=0.)

calc = LAMMPSlib(lammps_header=['units real',
                                'atom_style charge',
                                'atom_modify map array sort 0 0',
                                'boundary f f f'],
                 lmpcmds=['pair_style reax/c NULL checkqeq no',
                          'pair_coeff * * ffield.reax.cho O C H',
# 'fix 1 all qeq/reax 1 0.0 10.0 1e-12 reax/c'
                          'fix 1 all qeq/shielded 1 10 1e-12 1000 reax/c'
                          ],
                 log_file='test.log',
                 keep_alive=True)
atoms.set_calculator(calc)

f_analytical = atoms.get_forces()
f_numerical = calc.calculate_numerical_forces(atoms, 1e-5)
print('Ratio:\n', f_analytical / f_numerical)

And the new output:
Ratio:
[[0.99977128 1.00014625 1.00013576]
[1.00007552 1.00004205 1.00004367]
[0.99965457 0.99950924 0.99950043]
[1.00005591 1.00006447 1.00005361]
[1.00011362 1.00003254 1.00003975]
[0.99918412 0.97405552 0.99171004]
[1.01947762 0.99642427 0.99975086]
[0.86082277 0.99877919 0.95351369]
[1.00011399 1.00024284 1.0000426 ]
[1.00023009 1.00006462 1.0022227 ]
[0.99991101 0.99991704 0.99913163]
[1.01345352 0.99854881 0.99968067]]

Incidentally, I also tried using the old Fortran-based `pair_style reax` with the latest stable release of LAMMPS, and I still saw the same errors in the gradients. I’m starting to think this might be a fundamental ReaxFF issue, rather than a LAMMPS issue.

Eric