Lost atom with LAMMPS calculator

Hi,

I tried to use the LAMMPS calculator in ASE to compute the energy of twisted bilayer graphene. I used the input script below.

import ase
from ase.io import lammpsdata
from ase.calculators.lammpslib import LAMMPSlib
from mpi4py import MPI
data = lammpsdata.read_lammps_data('original.data',style="atomic")
data.set_pbc([1,1,0])
cmds = ["pair_style hybrid/overlay rebo kolmogorov/crespi/z 16",
        "pair_coeff * * rebo CH.rebo C C",
        "pair_coeff 1 2 kolmogorov/crespi/z CC.KC C C"]
head = ['units metal',
        'box tilt large',
        'atom_style atomic',
        'atom_modify     map yes',
        'processors * * 1']
lammps =  LAMMPSlib(lmpcmds=cmds,atom_types={'H':1,'He':2},atom_type_masses={'H':12.011,'He':12.011},log_file='lammps.log',keep_alive=True,lammps_header=head)
data.calc = lammps
print("Energy = {:23.17e}".format(data.get_potential_energy()))

It works fine if I run the script with 6 cores. But the error (please see below) occurs when I run the script with more than 8 cores.

LAMMPS (23 Jun 2022 - Update 1)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
  using 1 OpenMP thread(s) per MPI task
units metal
box tilt large
atom_style atomic
atom_modify     map yes
processors * * 1
boundary p p f
box tilt large
region cell prism    0 61.780847 0 53.503783 0 20.0     30.890424 0.0 0.0     units box
create_box 2 cell
Created triclinic box = (0 0 0) to (61.780847 53.503783 20) with tilt (30.890424 0 0)
WARNING: Triclinic box skew is large (src/domain.cpp:224)
  2 by 4 by 1 MPI processor grid
echo none
Created 2524 atoms
  using lattice units in triclinic box = (0 0 0) to (61.780847 53.503783 20) with tilt (30.890424 0 0)
  create_atoms CPU = 0.004 seconds
Setting atom values ...
  1 settings made for type
Setting atom values ...
  1 settings made for type
Setting atom values ...
  1 settings made for type
...
pair_style hybrid/overlay rebo kolmogorov/crespi/z 16
pair_coeff * * rebo CH.rebo C C
Reading rebo potential file CH.rebo with DATE: 2018-7-3
pair_coeff 1 2 kolmogorov/crespi/z CC.KC C C
mass 1 12.010999996910237896940998325590
mass 2 12.010999996910237896940998325590
variable pxx equal pxx
variable pyy equal pyy
variable pzz equal pzz
variable pxy equal pxy
variable pxz equal pxz
variable pyz equal pyz
thermo_style custom pe pxx emol ecoul
variable fx atom fx
variable fy atom fy
variable fz atom fz
variable pe equal pe
neigh_modify delay 0 every 1 check yes
echo none
run 0
WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
Neighbor list info ...
  update every 1 steps, delay 0 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 18
  ghost atom cutoff = 18
  binsize = 9, bins = 11 6 3
  3 neighbor lists, perpetual/occasional/extra = 3 0 0
  (1) pair rebo, perpetual
      attributes: full, newton on, ghost
      pair build: full/bin/ghost
      stencil: full/ghost/bin/3d
      bin: standard
  (2) pair kolmogorov/crespi/z, perpetual, skip from (3)
      attributes: half, newton on
      pair build: skip
      stencil: none
      bin: none
  (3) neighbor class addition, perpetual
      attributes: half, newton on
      pair build: half/bin/newton/tri
      stencil: half/bin/3d/tri
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.63 | 5.635 | 5.65 Mbytes
    PotEng          Pxx           E_mol          E_coul
ERROR: Lost atoms: original 2524 current 1909 (src/thermo.cpp:481)
Last command: run 0

I already tested the input script with LAMMPS and it works well. It seems that there is an issue with an interface between LAMMPSlib and ASE. If anyone has a comment or thought, please let me know.

Thank you,