hybrid of gcmc and nvt and atom indices

when combining fix gcmc and fix nvt into a hybrid run,
I always encounter "missing bond atoms" errors if the
number of molecules inserted by gcmc is very low.

A pure gcmc run under the same conditions *will run fine*,
combination with nvt in the style of

  ...
  fix md all nvt temp $t $t 100
  fix mc water gcmc $k $k $k 0 29491 $t $u 0.15 mol SPC region Central maxangle 75 full_energy
  ...
(t = 300, k = 100, "Central"= region full box except a waver at 0 < z < 10)

will certainly fail somewhere within a few thousand time steps
with an error like:

Step Atoms Temp C1 C2 PotEng KinEng TotEng Press Density
        0 3400 70.520741 300 0 -701295.15 714.50017 -700580.65 704.92812 0.24425807
      100 3415 427.32775 1221.8534 -66.232905 -705336.2 4348.6947 -700987.5 1154.2998 0.24462181
...
     1400 3442 291.62362 276.70332 -90.296453 -708282.3 2991.1739 -705291.13 -1514.6458 0.24527655
     1500 3439 296.18342 283.63146 -205.49677 -708277.97 3035.2951 -705242.67 -1849.1854 0.2452038
     1600 3457 293.31906 281.86667 -137.91897 -708138.09 3021.679 -705116.41 -1385.0266 0.24564029
ERROR on proc 0: Bond atoms 3977 3978 missing on proc 0 at step 1631 (../neigh_bond.cpp:65)

As can be seen, the numbers of the missing atoms
are much higher than the current number of atoms in the
system (which is: a water vapor above a silicate wafer).

This does, oddly enough, never happen if the density of the
added molecule (in fact, this is flexible spc water without
shake constraints) is above some threshold.

Is it possible that the integrator, when looking up intramolecular
terms, is confused by "holes" in the non-consecutive sequence of atom numbers? But why does this happen in very thin systems only?

(Tested with the latests stable LAMMPS and an actual git version,
I can provide input files if necessary.)

Thanks & Regards

M.

Please send me a small example that exhibits the behavior. I will take a look.