"hbondchk failed" error while using hybrid pair style (reax/c & lj/cuc/coul/cut)

Greetings,

I am trying to run a simulation in which there are two solute molecules (modelled with reax/c), surrounded by solvent (modelled with OPLS). The simulation seems to work fine until about 1000-4000 iterations in (the exact number depends mostly on the time step), when LAMMPS crashes, producing the “hbondchk failed” error below. There seems to be nothing unusual in the dumped data prior to the error as far as I can see.
I’ve tried changing the safezone and mincap parameters, but this seems to have no effect

I’ve read that hbondchk errors can often be avoided by increasing skin distance; however, at best, this only seems to delay the crash

If I delete all the OPLS molecules (using delete_atoms), the error does not occur (at least, not within the number of iterations I’ve tried).

I have tried running this with both the latest stable version of LAMMPS (9/Dec/2014), and the 30/Jan/2015 version. The input file and data file are attached, and the output is below.

Might anyone be able to provide insight on what’s causing this and how to fix it?

The output:
LAMMPS (30 Jan 2015)
Reading data file …
orthogonal box = (10.3911 10.3911 10.3911) to (49.6089 49.6089 49.6089)
1 by 1 by 1 MPI processor grid
reading atoms …
126 atoms
scanning bonds …
4 = max bonds/atom
scanning angles …
6 = max angles/atom
scanning dihedrals …
12 = max dihedrals/atom
reading bonds …
70 bonds
reading angles …
124 angles
reading dihedrals …
144 dihedrals
Finding 1-2 1-3 1-4 neighbors …
4 = max # of 1-2 neighbors
7 = max # of 1-3 neighbors
14 = max # of 1-4 neighbors
16 = max # of special neighbors
56 atoms in group solute
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (…/pair.cpp:203)
Setting up run …
WARNING: Inconsistent image flags (…/domain.cpp:678)
Memory usage per processor = 20.5791 Mbytes
Step Temp PotEng TotEng Press Volume
0 363.15 -7024.806 -6889.4959 2674.505 60318
1 457.68833 -7065.7631 -6895.2279 2058.7727 60320.341
2 598.97361 -7123.8869 -6900.7085 1356.2881 60326.287
3 670.10842 -7146.2079 -6896.5247 1085.7373 60334.609
4 684.86777 -7153.1893 -6898.0067 782.01691 60344.835
5 672.92561 -7150.6515 -6899.9186 412.67248 60356.434
6 663.20755 -7146.3697 -6899.2578 410.78304 60368.76
7 667.55298 -7149.8988 -6901.1677 181.74156 60381.81
8 661.37186 -7148.1346 -6901.7066 -196.42309 60395.183

995 341.10028 -7278.7208 -7151.6264 256.41472 71430.528
996 357.86048 -7283.4336 -7150.0943 278.39294 71445.162
997 375.27844 -7290.1988 -7150.3696 67.43826 71460.471
998 376.89536 -7292.5313 -7152.0996 87.085242 71475.937
999 370.38435 -7289.9491 -7151.9434 57.758139 71491.608
1000 366.43567 -7288.5226 -7151.9882 -40.082169 71507.412
step1000-hbondchk failed: H=3 end(H)=130 str(H+1)=100
application called MPI_Abort(MPI_COMM_WORLD, -14) - process 0

ffield.reax (16.7 KB)

loop.data (11.2 KB)

loop.in (5.57 KB)

I don’t see anything obviously wrong with your script.

You don’t need hybrid/overlap - hybrid is fine.

A timestep of 1 fmsec is somewhat big. But reducing
it to 0.1 still gives the error, just after more steps.

The error is coming from within ReaxFF. Presumably
the water is altering the dynamics in a bad way.

Possibly Ray (CCd) can comment on why ReaxFF would
be throwing that error. You might also check that

two atoms (ReaxFF atom and water atom) are not
ever getting too close together.

Steve

A timestep of 1 fs is indeed too big for ReaxFF - I suggest no bigger than 0.2 fs.

If I change this line, “fix 1 all npt temp 363.15 363.15 100.0 iso 1.0 1.0 1000.0” to “fix 1 solute …” then no error occurs. So ReaxFF fails when it is building the hydrogen bond list for the NULL atom types. I will have to take a deeper look.

Ray

I don’t think ReaxFF should be “seeing” any of the atoms with NULL types. I.e. the neigh
list it is passed should only contain atoms that are the ReaxFF types. For both
the I and J indices. It should also never be doing a loop over all nlocal atoms, only
over what is in its passed neigh list, for both the I,J indices.

If you do fix 1 solute npt, then the ReaxFF atoms will not be time integrated, and thus
not move. So that is probably why the problem goes away - the config never changes.

Steve

I don't think ReaxFF should be "seeing" any of the atoms with NULL types.
I.e. the neigh
list it is passed should only contain atoms that are the ReaxFF types.
For both
the I and J indices. It should also never be doing a loop over all nlocal
atoms, only
over what is in its passed neigh list, for both the I,J indices.

Yes, that could be the problem in reax/c implementation, particularly in
reaxc_force.cpp that the hydrogen bond neigh list is looping over
nlocal+nghost.

If you do fix 1 solute npt, then the ReaxFF atoms will not be time
integrated, and thus
not move. So that is probably why the problem goes away - the config
never changes.

No, it is the opposite - ReaxFF atoms are in the group solute. So when
only ReaxFF atoms move there is no hbondchk error.

Ray

Hi,
Thanks to both of you for your response! It does appear the NULL atom types may have something to do with it, as the error does not occur if I put actual atom types there. Would you be able to mention any specific spots in the source where you think the issue may be happening? I may try to take a look at it.
Thanks,

Ryan

So the problem goes away if atom sorting is turned off (atom sort 0 0). This has to do with how reax/c checks the hydrogen bond list when it is being used with hybrid or hybrid/overlay. With sorting tuned off, the check does not throw false error. Please try if it works for you.

Cheers,
Ray