Python Interface to Lammps scatter_atoms wont work as expected

Dear Lammps Users,
As a follow up to my previous post titled “ Creating Multiple PyLammps instances … “
I thought that the issue stems from creating multiple instances however that might not be case.
Here I create a lammps instance with two different force fields for test purposes:
  #This is the opls part:

    lmp_opls.command("units real")
    lmp_opls.command("atom_style full")
    lmp_opls.command("pair_style lj/cut/coul/debye 10.0 8.0")
    lmp_opls.command("bond_style harmonic")
    lmp_opls.command("angle_style harmonic")
    lmp_opls.command("dihedral_style opls")
    lmp_opls.command("special_bonds lj/coul 0.0 0.0 0.5")
    lmp_opls.command("neighbor 2.0 bin")
    lmp_opls.command("neigh_modify every 1 delay 0 check yes")
    lmp_opls.command("timestep 1.0")
    lmp_opls.command("thermo 100”)
    lmp_opls.command(“read_data DATAFILE.lmp”)

  #This is for the reaxff part
    lmp.command("units real")
    lmp.command("atom_style charge")
    lmp.command("pair_style reax/c NULL")
    lmp.command("pair_coeff * * %s s" (ffield, typeorder))
    lmp.command("thermo 100")
    symbols = typeorder.split()
    for i,symbol in enumerate(symbols):
        lmp.command(" mass %d f" ( i+1, symbol_to_mass(symbol)))
    lmp.command("fix qeq all qeq/reax 1 0.0 10.0 1e-6 reax/c” )
#run test simulation with reaxff
    lmp.command('fix run_sim all %s temp %f f 100\.0 '(simtype, tempi, tempf)
#run test simulation with opls
    lmp_opls.command('fix run_sim all %s temp %f f 100\.0 '(simtype, tempi, tempf)

up to here everything works fine.
for both opls and reaxff instances reads same system.
What I want to do here is to update positions of atoms in the “reaxff" simulation with the positions of atoms from the "opls" simulation

To test this first I gather “x” from opls:
data = lmp_opls.gather_atoms("x",1,3)
data[0] = 1.2 # modify first atom’s first position (x coordinate) and scatter back to “opls”
lmp_opls.scatter_atoms("x",1,3,data)
#to check if it is done gather back:
data = lmp_opls.gather_atoms("x",1,3)
print '---> ','me ' ,me, data[0:3] # print
output is:
---> me 0 [13.877459877839492, 16.682006077036945, 14.031749182727161]
---> me 0 [1.2, 16.682006077036945, 14.031749182727161]
It looks like it works.
Now I try same thing with “reaxff” instance
        datar = lmp_reax.gather_atoms("x",1,3)
        print'---->','me Reax' ,me, datar[0:3]
        datar[0] = 1.3
        lmp_reax.scatter_atoms("x",1,3,datar)
        datar = lmp_reax.gather_atoms("x",1,3)
        print'---->','me Reax' ,me, datar[0:3]
But now it gives this error:

WARNING: Library error in lammps_scatter_atoms (../library.cpp:1275)

here relevant line from the library.cpp:

    if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0)
      flag = 1;
    if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
    if (lmp->atom->map_style == 0) flag = 1;
    if (flag) {
      if (lmp->comm->me == 0)
        lmp->error->warning(FLERR,"Library error in lammps_scatter_atoms");
      return;
    }

Probably the problem is with map_style here but I cant come up with a solution.
Helps greatly appreciated.
Regards.

Dundar Yilmaz, Ph. D.
Associate Research Professor
Department of Mechanical and Nuclear Engineering
138 Research East Building
Penn State University
University Park, PA 16802

The atom map is created by default only for systems with permanent bonds. You need to request it explicitly for ReaxFF simulation through atom_modify command.

Link: https://lammps.sandia.gov/doc/atom_modify.html

Michał

Thank you Michal
But I never use atom_modify command with reaxff simulations on LAMMPS. Why Lammps() instance checks this to scatter atoms?

By default, the atoms on a given thread are identified by their local ID and this is enough for ReaxFF. However, scatter_atoms() (and some other functions/commands) requires the global IDs of atoms. An atom map links local IDs with global ones.

Michał