/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#include "atom.h"
#include "atom_vec.h"
#include "style_atom.h"  // IWYU pragma: keep

#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "group.h"
#include "input.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "molecule.h"
#include "neighbor.h"
#include "update.h"
#include "variable.h"

#include "library.h"

#include <algorithm>
#include <cstring>

#ifdef LMP_USER_INTEL
#include "neigh_request.h"
#endif

using namespace LAMMPS_NS;
using namespace MathConst;

#define DELTA 1
#define DELTA_PERATOM 64
#define EPSILON 1.0e-6

/* ---------------------------------------------------------------------- */

/** \class LAMMPS_NS::Atom
 *  \brief Class to provide access to atom data

\verbatim embed:rst
The Atom class provides access to atom style related global settings and
per-atom data that is stored with atoms and migrates with them from
sub-domain to sub-domain as atoms move around.  This includes topology
data, which is stored with either one specific atom or all atoms involved
depending on the settings of the :doc:`newton command <newton>`.

The actual per-atom data is allocated and managed by one of the various
classes derived from the AtomVec class as determined by
the :doc:`atom_style command <atom_style>`.  The pointers in the Atom class
are updated by the AtomVec class as needed.
\endverbatim
 */

/** Atom class constructor
 *
 * This resets and initialized all kinds of settings,
 * parameters, and pointer variables for per-atom arrays.
 * This also initializes the factory for creating
 * instances of classes derived from the AtomVec base
 * class, which correspond to the selected atom style.
 *
 * \param  lmp  pointer to the base LAMMPS class */

Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
{
  natoms = 0;
  nlocal = nghost = nmax = 0;
  ntypes = 0;
  nellipsoids = nlines = ntris = nbodies = 0;
  nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
  nbonds = nangles = ndihedrals = nimpropers = 0;

  firstgroupname = nullptr;
  sortfreq = 1000;
  nextsort = 0;
  userbinsize = 0.0;
  maxbin = maxnext = 0;
  binhead = nullptr;
  next = permute = nullptr;

  // data structure with info on per-atom vectors/arrays

  nperatom = maxperatom = 0;
  peratom = nullptr;

  // --------------------------------------------------------------------
  // 1st customization section: customize by adding new per-atom variables

  tag = nullptr;
  type = mask = nullptr;
  image = nullptr;
  x = v = f = nullptr;

  // charged and dipolar particles

  q = nullptr;
  mu = nullptr;

  // finite-size particles

  omega = angmom = torque = nullptr;
  radius = rmass = nullptr;
  ellipsoid = line = tri = body = nullptr;

  // molecular systems

  molecule = nullptr;
  molindex = molatom = nullptr;

  bond_per_atom =  extra_bond_per_atom = 0;
  num_bond = nullptr;
  bond_type = nullptr;
  bond_atom = nullptr;

  angle_per_atom = extra_angle_per_atom = 0;
  num_angle = nullptr;
  angle_type = nullptr;
  angle_atom1 = angle_atom2 = angle_atom3 = nullptr;

  dihedral_per_atom = extra_dihedral_per_atom = 0;
  num_dihedral = nullptr;
  dihedral_type = nullptr;
  dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = nullptr;

  improper_per_atom = extra_improper_per_atom = 0;
  num_improper = nullptr;
  improper_type = nullptr;
  improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = nullptr;

  maxspecial = 1;
  nspecial = nullptr;
  special = nullptr;

  // PERI package

  vfrac = s0 = nullptr;
  x0 = nullptr;

  // SPIN package

  sp = fm = fm_long = nullptr;

  // USER-EFF and USER-AWPMD packages

  spin = nullptr;
  eradius = ervel = erforce = nullptr;
  ervelforce = nullptr;
  cs = csforce = vforce = nullptr;
  etag = nullptr;

  // USER-DPD package

  uCond = uMech = uChem = uCG = uCGnew = nullptr;
  duChem = dpdTheta = nullptr;

  // USER-MESO package

  cc = cc_flux = nullptr;
  edpd_temp = edpd_flux = edpd_cv = nullptr;

  // USER-MESONT package

  length = nullptr;
  buckling = nullptr;
  bond_nt = nullptr;

  // USER-SMD package

  contact_radius = nullptr;
  smd_data_9 = nullptr;
  smd_stress = nullptr;
  eff_plastic_strain = nullptr;
  eff_plastic_strain_rate = nullptr;
  damage = nullptr;

  // USER-SPH package

  rho = drho = esph = desph = cv = nullptr;
  vest = nullptr;

  // end of customization section
  // --------------------------------------------------------------------

  // user-defined molecules

  nmolecule = 0;
  molecules = nullptr;

  // custom atom arrays

  nivector = ndvector = 0;
  ivector = nullptr;
  dvector = nullptr;
  iname = dname = nullptr;

  // initialize atom style and array existence flags

  set_atomflag_defaults();

  // initialize peratom data structure

  peratom_create();

  // ntype-length arrays

  mass = nullptr;
  mass_setflag = nullptr;

  // callback lists & extra restart info

  nextra_grow = nextra_restart = nextra_border = 0;
  extra_grow = extra_restart = extra_border = nullptr;
  nextra_grow_max = nextra_restart_max = nextra_border_max = 0;
  nextra_store = 0;
  extra = nullptr;

  // default atom ID and mapping values

  tag_enable = 1;
  map_style = map_user = MAP_NONE;
  map_tag_max = -1;
  map_maxarray = map_nhash = map_nbucket = -1;

  max_same = 0;
  sametag = nullptr;
  map_array = nullptr;
  map_bucket = nullptr;
  map_hash = nullptr;

  unique_tags = nullptr;

  atom_style = nullptr;
  avec = nullptr;

  avec_map = new AtomVecCreatorMap();

#define ATOM_CLASS
#define AtomStyle(key,Class) \
  (*avec_map)[#key] = &avec_creator<Class>;
#include "style_atom.h"
#undef AtomStyle
#undef ATOM_CLASS
}

/* ---------------------------------------------------------------------- */

Atom::~Atom()
{
  delete [] atom_style;
  delete avec;
  delete avec_map;

  delete [] firstgroupname;
  memory->destroy(binhead);
  memory->destroy(next);
  memory->destroy(permute);

  memory->destroy(tag);
  memory->destroy(type);
  memory->destroy(mask);
  memory->destroy(image);
  memory->destroy(x);
  memory->destroy(v);
  memory->destroy(f);

  // delete peratom data struct

  for (int i = 0; i < nperatom; i++)
    delete [] peratom[i].name;
  memory->sfree(peratom);

  // delete custom atom arrays

  for (int i = 0; i < nivector; i++) {
    delete [] iname[i];
    memory->destroy(ivector[i]);
  }
  if (dvector != nullptr) {
    for (int i = 0; i < ndvector; i++) {
      delete [] dname[i];
      memory->destroy(dvector[i]);
    }
  }

  memory->sfree(iname);
  memory->sfree(dname);
  memory->sfree(ivector);
  memory->sfree(dvector);

  // delete user-defined molecules

  for (int i = 0; i < nmolecule; i++) delete molecules[i];
  memory->sfree(molecules);

  // delete per-type arrays

  delete [] mass;
  delete [] mass_setflag;

  // delete extra arrays

  memory->destroy(extra_grow);
  memory->destroy(extra_restart);
  memory->destroy(extra_border);
  memory->destroy(extra);

  // delete mapping data structures

  map_delete();

  delete unique_tags;
}

/* ----------------------------------------------------------------------
   copy modify settings from old Atom class to current Atom class
------------------------------------------------------------------------- */

void Atom::settings(Atom *old)
{
  tag_enable = old->tag_enable;
  map_user = old->map_user;
  map_style = old->map_style;
  sortfreq = old->sortfreq;
  userbinsize = old->userbinsize;
  if (old->firstgroupname) {
    int n = strlen(old->firstgroupname) + 1;
    firstgroupname = new char[n];
    strcpy(firstgroupname,old->firstgroupname);
  }
}

/* ----------------------------------------------------------------------
   one-time creation of peratom data structure
------------------------------------------------------------------------- */

void Atom::peratom_create()
{
  for (int i = 0; i < nperatom; i++)
    delete [] peratom[i].name;
  memory->sfree(peratom);

  peratom = nullptr;
  nperatom = maxperatom = 0;

  // --------------------------------------------------------------------
  // 2nd customization section: add peratom variables here, order does not matter
  // register tagint & imageint variables as INT or BIGINT

  int tagintsize = INT;
  if (sizeof(tagint) == 8) tagintsize = BIGINT;
  int imageintsize = INT;
  if (sizeof(imageint) == 8) imageintsize = BIGINT;

  add_peratom("id",&tag,tagintsize,0);
  add_peratom("type",&type,INT,0);
  add_peratom("mask",&mask,INT,0);
  add_peratom("image",&image,imageintsize,0);

  add_peratom("x",&x,DOUBLE,3);
  add_peratom("v",&v,DOUBLE,3);
  add_peratom("f",&f,DOUBLE,3,1);      // set per-thread flag

  add_peratom("rmass",&rmass,DOUBLE,0);
  add_peratom("q",&q,DOUBLE,0);
  add_peratom("mu",&mu,DOUBLE,4);
  add_peratom("mu3",&mu,DOUBLE,3);     // just first 3 values of mu[4]

  // finite size particles

  add_peratom("radius",&radius,DOUBLE,0);
  add_peratom("omega",&omega,DOUBLE,3);
  add_peratom("torque",&torque,DOUBLE,3,1);    // set per-thread flag
  add_peratom("angmom",&angmom,DOUBLE,3);

  add_peratom("ellipsoid",&ellipsoid,INT,0);
  add_peratom("line",&line,INT,0);
  add_peratom("tri",&tri,INT,0);
  add_peratom("body",&body,INT,0);

  // MOLECULE package

  add_peratom("molecule",&molecule,tagintsize,0);
  add_peratom("molindex",&molindex,INT,0);
  add_peratom("molatom",&molatom,INT,0);

  add_peratom("nspecial",&nspecial,INT,3);
  add_peratom_vary("special",&special,tagintsize,&maxspecial,&nspecial,3);

  add_peratom("num_bond",&num_bond,INT,0);
  add_peratom_vary("bond_type",&bond_type,INT,&bond_per_atom,&num_bond);
  add_peratom_vary("bond_atom",&bond_atom,tagintsize,&bond_per_atom,&num_bond);

  add_peratom("num_angle",&num_angle,INT,0);
  add_peratom_vary("angle_type",&angle_type,INT,&angle_per_atom,&num_angle);
  add_peratom_vary("angle_atom1",&angle_atom1,tagintsize,
                   &angle_per_atom,&num_angle);
  add_peratom_vary("angle_atom2",&angle_atom2,tagintsize,
                   &angle_per_atom,&num_angle);
  add_peratom_vary("angle_atom3",&angle_atom3,tagintsize,
                   &angle_per_atom,&num_angle);

  add_peratom("num_dihedral",&num_dihedral,INT,0);
  add_peratom_vary("dihedral_type",&dihedral_type,INT,
                   &dihedral_per_atom,&num_dihedral);
  add_peratom_vary("dihedral_atom1",&dihedral_atom1,tagintsize,
                   &dihedral_per_atom,&num_dihedral);
  add_peratom_vary("dihedral_atom2",&dihedral_atom2,tagintsize,
                   &dihedral_per_atom,&num_dihedral);
  add_peratom_vary("dihedral_atom3",&dihedral_atom3,tagintsize,
                   &dihedral_per_atom,&num_dihedral);
  add_peratom_vary("dihedral_atom4",&dihedral_atom4,tagintsize,
                   &dihedral_per_atom,&num_dihedral);

  add_peratom("num_improper",&num_improper,INT,0);
  add_peratom_vary("improper_type",&improper_type,INT,
                   &improper_per_atom,&num_improper);
  add_peratom_vary("improper_atom1",&improper_atom1,tagintsize,
                   &improper_per_atom,&num_improper);
  add_peratom_vary("improper_atom2",&improper_atom2,tagintsize,
                   &improper_per_atom,&num_improper);
  add_peratom_vary("improper_atom3",&improper_atom3,tagintsize,
                   &improper_per_atom,&num_improper);
  add_peratom_vary("improper_atom4",&improper_atom4,tagintsize,
                   &improper_per_atom,&num_improper);

  // PERI package

  add_peratom("vfrac",&vfrac,DOUBLE,0);
  add_peratom("s0",&s0,DOUBLE,0);
  add_peratom("x0",&x0,DOUBLE,3);

  // SPIN package

  add_peratom("sp",&sp,DOUBLE,4);
  add_peratom("fm",&fm,DOUBLE,3,1);
  add_peratom("fm_long",&fm_long,DOUBLE,3,1);

  // USER-EFF package

  add_peratom("spin",&spin,INT,0);
  add_peratom("eradius",&eradius,DOUBLE,0);
  add_peratom("ervel",&ervel,DOUBLE,0);
  add_peratom("erforce",&erforce,DOUBLE,0,1);     // set per-thread flag

  // USER-AWPMD package

  add_peratom("cs",&cs,DOUBLE,2);
  add_peratom("csforce",&csforce,DOUBLE,2);
  add_peratom("vforce",&vforce,DOUBLE,3);
  add_peratom("ervelforce",&ervelforce,DOUBLE,0);
  add_peratom("etag",&etag,INT,0);

  // USER-DPD package

  add_peratom("dpdTheta",&dpdTheta,DOUBLE,0);
  add_peratom("uCond",&uCond,DOUBLE,0);
  add_peratom("uMech",&uMech,DOUBLE,0);
  add_peratom("uChem",&uChem,DOUBLE,0);
  add_peratom("uCG",&uCG,DOUBLE,0);
  add_peratom("uCGnew",&uCGnew,DOUBLE,0);
  add_peratom("duChem",&duChem,DOUBLE,0);

  // USER-MESO package

  add_peratom("edpd_cv",&edpd_cv,DOUBLE,0);
  add_peratom("edpd_temp",&edpd_temp,DOUBLE,0);
  add_peratom("vest_temp",&vest_temp,DOUBLE,0);
  add_peratom("edpd_flux",&edpd_flux,DOUBLE,0,1);     // set per-thread flag
  add_peratom("cc",&cc,DOUBLE,1);
  add_peratom("cc_flux",&cc_flux,DOUBLE,1,1);         // set per-thread flag

  // USER-MESONT package

  add_peratom("length",&length,DOUBLE,0);
  add_peratom("buckling",&buckling,INT,0);
  add_peratom("bond_nt",&bond_nt,tagintsize,2);

  // USER-SPH package

  add_peratom("rho",&rho,DOUBLE,0);
  add_peratom("drho",&drho,DOUBLE,0,1);               // set per-thread flag
  add_peratom("esph",&esph,DOUBLE,0);
  add_peratom("desph",&desph,DOUBLE,0,1);             // set per-thread flag
  add_peratom("vest",&vest,DOUBLE,3);
  add_peratom("cv",&cv,DOUBLE,0);

  // USER-SMD package

  add_peratom("contact_radius",&contact_radius,DOUBLE,0);
  add_peratom("smd_data_9",&smd_data_9,DOUBLE,1);
  add_peratom("smd_stress",&smd_stress,DOUBLE,1);
  add_peratom("eff_plastic_strain",&eff_plastic_strain,DOUBLE,0);
  add_peratom("eff_plastic_strain_rate",&eff_plastic_strain_rate,DOUBLE,0);
  add_peratom("damage",&damage,DOUBLE,0);

  // end of customization section
  // --------------------------------------------------------------------
}

/* ----------------------------------------------------------------------
   add info for a single per-atom vector/array to PerAtom data struct
   cols = 0: per-atom vector
   cols = N: static per-atom array with N columns
   use add_peratom_vary() when column count varies per atom
------------------------------------------------------------------------- */

void Atom::add_peratom(const char *name, void *address,
                       int datatype, int cols, int threadflag)
{
  if (nperatom == maxperatom) {
    maxperatom += DELTA_PERATOM;
    peratom = (PerAtom *)
      memory->srealloc(peratom,maxperatom*sizeof(PerAtom),"atom:peratom");
  }

  int n = strlen(name) + 1;
  peratom[nperatom].name = new char[n];
  strcpy(peratom[nperatom].name,name);
  peratom[nperatom].address = address;
  peratom[nperatom].datatype = datatype;
  peratom[nperatom].cols = cols;
  peratom[nperatom].threadflag = threadflag;
  peratom[nperatom].address_length = nullptr;

  nperatom++;
}

/* ----------------------------------------------------------------------
   change the column count of an existing peratom array entry
   allows atom_style to specify column count as an argument
   see atom_style tdpd as an example
------------------------------------------------------------------------- */

void Atom::add_peratom_change_columns(const char *name, int cols)
{
  for (int i = 0; i < nperatom; i++) {
    if (strcmp(name,peratom[i].name) == 0) {
            peratom[i].cols = cols;
            return;
    }
  }
  error->all(FLERR,"Could not find name of peratom array for column change");
}

/* ----------------------------------------------------------------------
   add info for a single per-atom array to PerAtom data struct
   cols = address of int variable with max columns per atom
   for collength = 0:
     length = address of peratom vector with column count per atom
     e.g. num_bond
   for collength = N:
     length = address of peratom array with column count per atom
     collength = index of column (1 to N) in peratom array with count
     e.g. nspecial
------------------------------------------------------------------------- */

void Atom::add_peratom_vary(const char *name, void *address,
                            int datatype, int *cols, void *length, int collength)
{
  if (nperatom == maxperatom) {
    maxperatom += DELTA_PERATOM;
    peratom = (PerAtom *)
      memory->srealloc(peratom,maxperatom*sizeof(PerAtom),"atom:peratom");
  }

  int n = strlen(name) + 1;
  peratom[nperatom].name = new char[n];
  strcpy(peratom[nperatom].name,name);
  peratom[nperatom].address = address;
  peratom[nperatom].datatype = datatype;
  peratom[nperatom].cols = -1;
  peratom[nperatom].threadflag = 0;
  peratom[nperatom].address_maxcols = cols;
  peratom[nperatom].address_length = length;
  peratom[nperatom].collength = collength;

  nperatom++;
}

/* ----------------------------------------------------------------------
   add info for a single per-atom array to PerAtom data struct
------------------------------------------------------------------------- */

void Atom::set_atomflag_defaults()
{
  // --------------------------------------------------------------------
  // 3rd customization section: customize by adding new flag
  // identical list as 2nd customization in atom.h

  sphere_flag = ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
  peri_flag = electron_flag = 0;
  wavepacket_flag = sph_flag = 0;
  molecule_flag = molindex_flag = molatom_flag = 0;
  q_flag = mu_flag = 0;
  rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
  vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
  cs_flag = csforce_flag = vforce_flag = ervelforce_flag = etag_flag = 0;
  rho_flag = esph_flag = cv_flag = vest_flag = 0;
  dpd_flag = edpd_flag = tdpd_flag = 0;
  sp_flag = 0;
  x0_flag = 0;
  smd_flag = damage_flag = 0;
  mesont_flag = 0;
  contact_radius_flag = smd_data_9_flag = smd_stress_flag = 0;
  eff_plastic_strain_flag = eff_plastic_strain_rate_flag = 0;

  pdscale = 1.0;
}

/* ----------------------------------------------------------------------
   create an AtomVec style
   called from lammps.cpp, input script, restart file, replicate
------------------------------------------------------------------------- */

void Atom::create_avec(const std::string &style, int narg, char **arg, int trysuffix)
{
  delete [] atom_style;
  if (avec) delete avec;
  atom_style = nullptr;
  avec = nullptr;

  // unset atom style and array existence flags
  // may have been set by old avec

  set_atomflag_defaults();

  // create instance of AtomVec
  // use grow() to initialize atom-based arrays to length 1
  //   so that x[0][0] can always be referenced even if proc has no atoms

  int sflag;
  avec = new_avec(style,trysuffix,sflag);
  avec->store_args(narg,arg);
  avec->process_args(narg,arg);
  avec->grow(1);

  if (sflag) {
    std::string estyle = style + "/";
    if (sflag == 1) estyle += lmp->suffix;
    else estyle += lmp->suffix2;
    atom_style = new char[estyle.size()+1];
    strcpy(atom_style,estyle.c_str());
  } else {
    atom_style = new char[style.size()+1];
    strcpy(atom_style,style.c_str());
  }

  // if molecular system:
  // atom IDs must be defined
  // force atom map to be created
  // map style will be reset to array vs hash to by map_init()

  molecular = avec->molecular;
  if (molecular && tag_enable == 0)
    error->all(FLERR,"Atom IDs must be used for molecular systems");
  if (molecular != Atom::ATOMIC) map_style = MAP_YES;
}

/* ----------------------------------------------------------------------
   generate an AtomVec class, first with suffix appended
------------------------------------------------------------------------- */

AtomVec *Atom::new_avec(const std::string &style, int trysuffix, int &sflag)
{
  if (trysuffix && lmp->suffix_enable) {
    if (lmp->suffix) {
      sflag = 1;
      std::string estyle = style + "/" + lmp->suffix;
      if (avec_map->find(estyle) != avec_map->end()) {
        AtomVecCreator &avec_creator = (*avec_map)[estyle];
        return avec_creator(lmp);
      }
    }

    if (lmp->suffix2) {
      sflag = 2;
      std::string estyle = style + "/" + lmp->suffix2;
      if (avec_map->find(estyle) != avec_map->end()) {
        AtomVecCreator &avec_creator = (*avec_map)[estyle];
        return avec_creator(lmp);
      }
    }
  }

  sflag = 0;
  if (avec_map->find(style) != avec_map->end()) {
    AtomVecCreator &avec_creator = (*avec_map)[style];
    return avec_creator(lmp);
  }

  error->all(FLERR,utils::check_packages_for_style("atom",style,lmp));
  return nullptr;
}

/* ----------------------------------------------------------------------
   one instance per AtomVec style in style_atom.h
------------------------------------------------------------------------- */

template <typename T>
AtomVec *Atom::avec_creator(LAMMPS *lmp)
{
  return new T(lmp);
}

/* ---------------------------------------------------------------------- */

void Atom::init()
{
  // delete extra array since it doesn't persist past first run

  if (nextra_store) {
    memory->destroy(extra);
    extra = nullptr;
    nextra_store = 0;
  }

  // check arrays that are atom type in length

  check_mass(FLERR);

  // setup of firstgroup

  if (firstgroupname) {
    firstgroup = group->find(firstgroupname);
    if (firstgroup < 0)
      error->all(FLERR,"Could not find atom_modify first group ID");
  } else firstgroup = -1;

  // init AtomVec

  avec->init();
}

/* ---------------------------------------------------------------------- */

void Atom::setup()
{
  // setup bins for sorting
  // cannot do this in init() because uses neighbor cutoff

  if (sortfreq > 0) setup_sort_bins();
}

/* ----------------------------------------------------------------------
   return ptr to AtomVec class if matches style or to matching hybrid sub-class
   return nullptr if no match
------------------------------------------------------------------------- */

AtomVec *Atom::style_match(const char *style)
{
  if (strcmp(atom_style,style) == 0) return avec;
  else if (strcmp(atom_style,"hybrid") == 0) {
    auto avec_hybrid = (AtomVecHybrid *) avec;
    for (int i = 0; i < avec_hybrid->nstyles; i++)
      if (strcmp(avec_hybrid->keywords[i],style) == 0)
        return avec_hybrid->styles[i];
  }
  return nullptr;
}

/* ----------------------------------------------------------------------
   modify parameters of the atom style
   some options can only be invoked before simulation box is defined
   first and sort options cannot be used together
------------------------------------------------------------------------- */

void Atom::modify_params(int narg, char **arg)
{
  if (narg == 0) error->all(FLERR,"Illegal atom_modify command");

  int iarg = 0;
  while (iarg < narg) {
    if (strcmp(arg[iarg],"id") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
      if (domain->box_exist)
        error->all(FLERR,
                   "Atom_modify id command after simulation box is defined");
      if (strcmp(arg[iarg+1],"yes") == 0) tag_enable = 1;
      else if (strcmp(arg[iarg+1],"no") == 0) tag_enable = 0;
      else error->all(FLERR,"Illegal atom_modify command");
      iarg += 2;
    } else if (strcmp(arg[iarg],"map") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
      if (domain->box_exist)
        error->all(FLERR,
                   "Atom_modify map command after simulation box is defined");
      if (strcmp(arg[iarg+1],"array") == 0) map_user = 1;
      else if (strcmp(arg[iarg+1],"hash") == 0) map_user = 2;
      else if (strcmp(arg[iarg+1],"yes") == 0) map_user = 3;
      else error->all(FLERR,"Illegal atom_modify command");
      map_style = map_user;
      iarg += 2;
    } else if (strcmp(arg[iarg],"first") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
      if (strcmp(arg[iarg+1],"all") == 0) {
        delete [] firstgroupname;
        firstgroupname = nullptr;
      } else {
        int n = strlen(arg[iarg+1]) + 1;
        firstgroupname = new char[n];
        strcpy(firstgroupname,arg[iarg+1]);
        sortfreq = 0;
      }
      iarg += 2;
    } else if (strcmp(arg[iarg],"sort") == 0) {
      if (iarg+3 > narg) error->all(FLERR,"Illegal atom_modify command");
      sortfreq = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
      userbinsize = utils::numeric(FLERR,arg[iarg+2],false,lmp);
      if (sortfreq < 0 || userbinsize < 0.0)
        error->all(FLERR,"Illegal atom_modify command");
      if (sortfreq >= 0 && firstgroupname)
        error->all(FLERR,"Atom_modify sort and first options "
                   "cannot be used together");
      iarg += 3;
    } else error->all(FLERR,"Illegal atom_modify command");
  }
}

/* ----------------------------------------------------------------------
   check that atom IDs are valid
   error if any atom ID < 0 or atom ID = MAXTAGINT
   if any atom ID > 0, error if any atom ID == 0
   if any atom ID > 0, error if tag_enable = 0
   if all atom IDs = 0, tag_enable must be 0
   if max atom IDs < natoms, must be duplicates
   OK if max atom IDs > natoms
   NOTE: not fully checking that atom IDs are unique
------------------------------------------------------------------------- */

void Atom::tag_check()
{
  tagint min = MAXTAGINT;
  tagint max = 0;

  for (int i = 0; i < nlocal; i++) {
    min = MIN(min,tag[i]);
    max = MAX(max,tag[i]);
  }

  tagint minall,maxall;
  MPI_Allreduce(&min,&minall,1,MPI_LMP_TAGINT,MPI_MIN,world);
  MPI_Allreduce(&max,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world);

  if (minall < 0) error->all(FLERR,"One or more Atom IDs is negative");
  if (maxall >= MAXTAGINT) error->all(FLERR,"One or more atom IDs is too big");
  if (maxall > 0 && minall == 0)
    error->all(FLERR,"One or more atom IDs is zero");
  if (maxall > 0 && tag_enable == 0)
    error->all(FLERR,"Non-zero atom IDs with atom_modify id = no");
  if (maxall == 0 && natoms && tag_enable)
    error->all(FLERR,"All atom IDs = 0 but atom_modify id = yes");
  if (tag_enable && maxall < natoms)
    error->all(FLERR,"Duplicate atom IDs exist");
}

/* ----------------------------------------------------------------------
   add unique tags to any atoms with tag = 0
   new tags are grouped by proc and start after max current tag
   called after creating new atoms
   error if new tags will exceed MAXTAGINT
------------------------------------------------------------------------- */

void Atom::tag_extend()
{
  // maxtag_all = max tag for all atoms

  tagint maxtag = 0;
  for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]);
  tagint maxtag_all;
  MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);

  // DEBUG: useful for generating 64-bit IDs even for small systems
  // use only when LAMMPS is compiled with BIGBIG

  //maxtag_all += 1000000000000;

  // notag = # of atoms I own with no tag (tag = 0)
  // notag_sum = # of total atoms on procs <= me with no tag

  bigint notag = 0;
  for (int i = 0; i < nlocal; i++) if (tag[i] == 0) notag++;

  bigint notag_total;
  MPI_Allreduce(&notag,&notag_total,1,MPI_LMP_BIGINT,MPI_SUM,world);
  if (notag_total >= MAXTAGINT)
    error->all(FLERR,"New atom IDs exceed maximum allowed ID");

  bigint notag_sum;
  MPI_Scan(&notag,&notag_sum,1,MPI_LMP_BIGINT,MPI_SUM,world);

  // itag = 1st new tag that my untagged atoms should use

  tagint itag = maxtag_all + notag_sum - notag + 1;
  for (int i = 0; i < nlocal; i++) if (tag[i] == 0) tag[i] = itag++;
}

/* ----------------------------------------------------------------------
   check that atom IDs span range from 1 to Natoms inclusive
   return 0 if mintag != 1 or maxtag != Natoms
   return 1 if OK
   doesn't actually check if all tag values are used
------------------------------------------------------------------------- */

int Atom::tag_consecutive()
{
  tagint idmin = MAXTAGINT;
  tagint idmax = 0;

  for (int i = 0; i < nlocal; i++) {
    idmin = MIN(idmin,tag[i]);
    idmax = MAX(idmax,tag[i]);
  }
  tagint idminall,idmaxall;
  MPI_Allreduce(&idmin,&idminall,1,MPI_LMP_TAGINT,MPI_MIN,world);
  MPI_Allreduce(&idmax,&idmaxall,1,MPI_LMP_TAGINT,MPI_MAX,world);

  if (idminall != 1 || idmaxall != natoms) return 0;
  return 1;
}

/* ----------------------------------------------------------------------
   check that bonus data settings are valid
   error if number of atoms with ellipsoid/line/tri/body flags
   are consistent with global setting.
------------------------------------------------------------------------- */

void Atom::bonus_check()
{
  bigint local_ellipsoids = 0, local_lines = 0, local_tris = 0;
  bigint local_bodies = 0, num_global;

  for (int i = 0; i < nlocal; ++i) {
    if (ellipsoid && (ellipsoid[i] >=0)) ++local_ellipsoids;
    if (line && (line[i] >=0)) ++local_lines;
    if (tri && (tri[i] >=0)) ++local_tris;
    if (body && (body[i] >=0)) ++local_bodies;
  }

  MPI_Allreduce(&local_ellipsoids,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
  if (nellipsoids != num_global)
    error->all(FLERR,"Inconsistent 'ellipsoids' header value and number of "
               "atoms with enabled ellipsoid flags");

  MPI_Allreduce(&local_lines,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
  if (nlines != num_global)
    error->all(FLERR,"Inconsistent 'lines' header value and number of "
               "atoms with enabled line flags");

  MPI_Allreduce(&local_tris,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
  if (ntris != num_global)
    error->all(FLERR,"Inconsistent 'tris' header value and number of "
               "atoms with enabled tri flags");

  MPI_Allreduce(&local_bodies,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
  if (nbodies != num_global)
    error->all(FLERR,"Inconsistent 'bodies' header value and number of "
               "atoms with enabled body flags");
}

/* ----------------------------------------------------------------------
   deallocate molecular topology arrays
   done before realloc with (possibly) new 2nd dimension set to
     correctly initialized per-atom values, e.g. bond_per_atom
   needs to be called whenever 2nd dimensions are changed
     and these arrays are already pre-allocated,
     e.g. due to grow(1) in create_avec()
------------------------------------------------------------------------- */

void Atom::deallocate_topology()
{
  memory->destroy(atom->bond_type);
  memory->destroy(atom->bond_atom);
  atom->bond_type = nullptr;
  atom->bond_atom = nullptr;

  memory->destroy(atom->angle_type);
  memory->destroy(atom->angle_atom1);
  memory->destroy(atom->angle_atom2);
  memory->destroy(atom->angle_atom3);
  atom->angle_type = nullptr;
  atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = nullptr;

  memory->destroy(atom->dihedral_type);
  memory->destroy(atom->dihedral_atom1);
  memory->destroy(atom->dihedral_atom2);
  memory->destroy(atom->dihedral_atom3);
  memory->destroy(atom->dihedral_atom4);
  atom->dihedral_type = nullptr;
  atom->dihedral_atom1 = atom->dihedral_atom2 =
    atom->dihedral_atom3 = atom->dihedral_atom4 = nullptr;

  memory->destroy(atom->improper_type);
  memory->destroy(atom->improper_atom1);
  memory->destroy(atom->improper_atom2);
  memory->destroy(atom->improper_atom3);
  memory->destroy(atom->improper_atom4);
  atom->improper_type = nullptr;
  atom->improper_atom1 = atom->improper_atom2 =
    atom->improper_atom3 = atom->improper_atom4 = nullptr;
}

/* ----------------------------------------------------------------------
   unpack N lines from Atom section of data file
   call style-specific routine to parse line
------------------------------------------------------------------------- */

void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
                      int type_offset, int shiftflag, double *shift)
{
  int m,xptr,iptr;
  imageint imagedata;
  double xdata[3],lamda[3];
  double *coord;
  char *next;

  next = strchr(buf,'\n');
  *next = '\0';
  int nwords = utils::trim_and_count_words(buf);
  *next = '\n';

  if (nwords != avec->size_data_atom && nwords != avec->size_data_atom + 3)
    error->all(FLERR,"Incorrect atom format in data file");

  char **values = new char*[nwords];

  // set bounds for my proc
  // if periodic and I am lo/hi proc, adjust bounds by EPSILON
  // insures all data atoms will be owned even with round-off

  int triclinic = domain->triclinic;

  double epsilon[3];
  if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON;
  else {
    epsilon[0] = domain->prd[0] * EPSILON;
    epsilon[1] = domain->prd[1] * EPSILON;
    epsilon[2] = domain->prd[2] * EPSILON;
  }

  double sublo[3],subhi[3];
  if (triclinic == 0) {
    sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
    sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
    sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
  } else {
    sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
    sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
    sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
  }

  if (comm->layout != Comm::LAYOUT_TILED) {
    if (domain->xperiodic) {
      if (comm->myloc[0] == 0) sublo[0] -= epsilon[0];
      if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += epsilon[0];
    }
    if (domain->yperiodic) {
      if (comm->myloc[1] == 0) sublo[1] -= epsilon[1];
      if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += epsilon[1];
    }
    if (domain->zperiodic) {
      if (comm->myloc[2] == 0) sublo[2] -= epsilon[2];
      if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += epsilon[2];
    }

  } else {
    if (domain->xperiodic) {
      if (comm->mysplit[0][0] == 0.0) sublo[0] -= epsilon[0];
      if (comm->mysplit[0][1] == 1.0) subhi[0] += epsilon[0];
    }
    if (domain->yperiodic) {
      if (comm->mysplit[1][0] == 0.0) sublo[1] -= epsilon[1];
      if (comm->mysplit[1][1] == 1.0) subhi[1] += epsilon[1];
    }
    if (domain->zperiodic) {
      if (comm->mysplit[2][0] == 0.0) sublo[2] -= epsilon[2];
      if (comm->mysplit[2][1] == 1.0) subhi[2] += epsilon[2];
    }
  }

  // xptr = which word in line starts xyz coords
  // iptr = which word in line starts ix,iy,iz image flags

  xptr = avec->xcol_data - 1;
  int imageflag = 0;
  if (nwords > avec->size_data_atom) imageflag = 1;
  if (imageflag) iptr = nwords - 3;

  // loop over lines of atom data
  // tokenize the line into values
  // extract xyz coords and image flags
  // remap atom into simulation box
  // if atom is in my sub-domain, unpack its values

  for (int i = 0; i < n; i++) {
    next = strchr(buf,'\n');

    values[0] = strtok(buf," \t\n\r\f");
    if (values[0] == nullptr)
      error->all(FLERR,"Incorrect atom format in data file");
    for (m = 1; m < nwords; m++) {
      values[m] = strtok(nullptr," \t\n\r\f");
      if (values[m] == nullptr)
        error->all(FLERR,"Incorrect atom format in data file");
    }

    int imx = 0;
    int imy = 0;
    int imz = 0;
    if (imageflag) {
      imx = utils::inumeric(FLERR,values[iptr],false,lmp);
      imy = utils::inumeric(FLERR,values[iptr+1],false,lmp);
      imz = utils::inumeric(FLERR,values[iptr+2],false,lmp);
      if ((domain->dimension == 2) && (imz != 0))
        error->all(FLERR,"Z-direction image flag must be 0 for 2d-systems");
    }
    imagedata = ((imageint) (imx + IMGMAX) & IMGMASK) |
        (((imageint) (imy + IMGMAX) & IMGMASK) << IMGBITS) |
        (((imageint) (imz + IMGMAX) & IMGMASK) << IMG2BITS);

    xdata[0] = utils::numeric(FLERR,values[xptr],false,lmp);
    xdata[1] = utils::numeric(FLERR,values[xptr+1],false,lmp);
    xdata[2] = utils::numeric(FLERR,values[xptr+2],false,lmp);
    if (shiftflag) {
      xdata[0] += shift[0];
      xdata[1] += shift[1];
      xdata[2] += shift[2];
    }

    domain->remap(xdata,imagedata);
    if (triclinic) {
      domain->x2lamda(xdata,lamda);
      coord = lamda;
    } else coord = xdata;

    if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
        coord[1] >= sublo[1] && coord[1] < subhi[1] &&
        coord[2] >= sublo[2] && coord[2] < subhi[2]) {
      avec->data_atom(xdata,imagedata,values);
      if (id_offset) tag[nlocal-1] += id_offset;
      if (mol_offset) molecule[nlocal-1] += mol_offset;
      if (type_offset) {
        type[nlocal-1] += type_offset;
        if (type[nlocal-1] > ntypes)
          error->one(FLERR,"Invalid atom type in Atoms section of data file");
      }
    }

    buf = next + 1;
  }

  delete [] values;
}

/* ----------------------------------------------------------------------
   unpack N lines from Velocity section of data file
   check that atom IDs are > 0 and <= map_tag_max
   call style-specific routine to parse line
------------------------------------------------------------------------- */

void Atom::data_vels(int n, char *buf, tagint id_offset)
{
  int j,m;
  tagint tagdata;
  char *next;

  next = strchr(buf,'\n');
  *next = '\0';
  int nwords = utils::trim_and_count_words(buf);
  *next = '\n';

  if (nwords != avec->size_data_vel)
    error->all(FLERR,"Incorrect velocity format in data file");

  char **values = new char*[nwords];

  // loop over lines of atom velocities
  // tokenize the line into values
  // if I own atom tag, unpack its values

  for (int i = 0; i < n; i++) {
    next = strchr(buf,'\n');

    values[0] = strtok(buf," \t\n\r\f");
    for (j = 1; j < nwords; j++)
      values[j] = strtok(nullptr," \t\n\r\f");

    tagdata = ATOTAGINT(values[0]) + id_offset;
    if (tagdata <= 0 || tagdata > map_tag_max)
      error->one(FLERR,"Invalid atom ID in Velocities section of data file");
    if ((m = map(tagdata)) >= 0) avec->data_vel(m,&values[1]);

    buf = next + 1;
  }

  delete [] values;
}

/* ----------------------------------------------------------------------
   process N bonds read into buf from data files
   if count is non-nullptr, just count bonds per atom
   else store them with atoms
   check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */

void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset,
                      int type_offset)
{
  int m,tmp,itype,rv;
  tagint atom1,atom2;
  char *next;
  int newton_bond = force->newton_bond;

  for (int i = 0; i < n; i++) {
    next = strchr(buf,'\n');
    *next = '\0';
    rv = sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT,
                &tmp,&itype,&atom1,&atom2);
    if (rv != 4)
      error->one(FLERR,"Incorrect format of Bonds section in data file");
    if (id_offset) {
      atom1 += id_offset;
      atom2 += id_offset;
    }
    itype += type_offset;

    if ((atom1 <= 0) || (atom1 > map_tag_max) ||
        (atom2 <= 0) || (atom2 > map_tag_max) || (atom1 == atom2))
      error->one(FLERR,"Invalid atom ID in Bonds section of data file");
    if (itype <= 0 || itype > nbondtypes)
      error->one(FLERR,"Invalid bond type in Bonds section of data file");
    if ((m = map(atom1)) >= 0) {
      if (count) count[m]++;
      else {
        bond_type[m][num_bond[m]] = itype;
        bond_atom[m][num_bond[m]] = atom2;
        num_bond[m]++;
      }
    }
    if (newton_bond == 0) {
      if ((m = map(atom2)) >= 0) {
        if (count) count[m]++;
        else {
          bond_type[m][num_bond[m]] = itype;
          bond_atom[m][num_bond[m]] = atom1;
          num_bond[m]++;
        }
      }
    }
    buf = next + 1;
  }
}

/* ----------------------------------------------------------------------
   process N angles read into buf from data files
   if count is non-nullptr, just count angles per atom
   else store them with atoms
   check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */

void Atom::data_angles(int n, char *buf, int *count, tagint id_offset,
                       int type_offset)
{
  int m,tmp,itype,rv;
  tagint atom1,atom2,atom3;
  char *next;
  int newton_bond = force->newton_bond;

  for (int i = 0; i < n; i++) {
    next = strchr(buf,'\n');
    *next = '\0';
    rv = sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
                &tmp,&itype,&atom1,&atom2,&atom3);
    if (rv != 5)
      error->one(FLERR,"Incorrect format of Angles section in data file");
    if (id_offset) {
      atom1 += id_offset;
      atom2 += id_offset;
      atom3 += id_offset;
    }
    itype += type_offset;

    if ((atom1 <= 0) || (atom1 > map_tag_max) ||
        (atom2 <= 0) || (atom2 > map_tag_max) ||
        (atom3 <= 0) || (atom3 > map_tag_max) ||
        (atom1 == atom2) || (atom1 == atom3) || (atom2 == atom3))
      error->one(FLERR,"Invalid atom ID in Angles section of data file");
    if (itype <= 0 || itype > nangletypes)
      error->one(FLERR,"Invalid angle type in Angles section of data file");
    if ((m = map(atom2)) >= 0) {
      if (count) count[m]++;
      else {
        angle_type[m][num_angle[m]] = itype;
        angle_atom1[m][num_angle[m]] = atom1;
        angle_atom2[m][num_angle[m]] = atom2;
        angle_atom3[m][num_angle[m]] = atom3;
        num_angle[m]++;
      }
    }
    if (newton_bond == 0) {
      if ((m = map(atom1)) >= 0) {
        if (count) count[m]++;
        else {
          angle_type[m][num_angle[m]] = itype;
          angle_atom1[m][num_angle[m]] = atom1;
          angle_atom2[m][num_angle[m]] = atom2;
          angle_atom3[m][num_angle[m]] = atom3;
          num_angle[m]++;
        }
      }
      if ((m = map(atom3)) >= 0) {
        if (count) count[m]++;
        else {
          angle_type[m][num_angle[m]] = itype;
          angle_atom1[m][num_angle[m]] = atom1;
          angle_atom2[m][num_angle[m]] = atom2;
          angle_atom3[m][num_angle[m]] = atom3;
          num_angle[m]++;
        }
      }
    }
    buf = next + 1;
  }
}

/* ----------------------------------------------------------------------
   process N dihedrals read into buf from data files
   if count is non-nullptr, just count diihedrals per atom
   else store them with atoms
   check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */

void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset,
                          int type_offset)
{
  int m,tmp,itype,rv;
  tagint atom1,atom2,atom3,atom4;
  char *next;
  int newton_bond = force->newton_bond;

  for (int i = 0; i < n; i++) {
    next = strchr(buf,'\n');
    *next = '\0';
    rv = sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT
                " " TAGINT_FORMAT " " TAGINT_FORMAT,
                &tmp,&itype,&atom1,&atom2,&atom3,&atom4);
    if (rv != 6)
      error->one(FLERR,"Incorrect format of Dihedrals section in data file");
    if (id_offset) {
      atom1 += id_offset;
      atom2 += id_offset;
      atom3 += id_offset;
      atom4 += id_offset;
    }
    itype += type_offset;

    if ((atom1 <= 0) || (atom1 > map_tag_max) ||
        (atom2 <= 0) || (atom2 > map_tag_max) ||
        (atom3 <= 0) || (atom3 > map_tag_max) ||
        (atom4 <= 0) || (atom4 > map_tag_max) ||
        (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
        (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
      error->one(FLERR,"Invalid atom ID in Dihedrals section of data file");
    if (itype <= 0 || itype > ndihedraltypes)
      error->one(FLERR,
                 "Invalid dihedral type in Dihedrals section of data file");
    if ((m = map(atom2)) >= 0) {
      if (count) count[m]++;
      else {
        dihedral_type[m][num_dihedral[m]] = itype;
        dihedral_atom1[m][num_dihedral[m]] = atom1;
        dihedral_atom2[m][num_dihedral[m]] = atom2;
        dihedral_atom3[m][num_dihedral[m]] = atom3;
        dihedral_atom4[m][num_dihedral[m]] = atom4;
        num_dihedral[m]++;
      }
    }
    if (newton_bond == 0) {
      if ((m = map(atom1)) >= 0) {
        if (count) count[m]++;
        else {
          dihedral_type[m][num_dihedral[m]] = itype;
          dihedral_atom1[m][num_dihedral[m]] = atom1;
          dihedral_atom2[m][num_dihedral[m]] = atom2;
          dihedral_atom3[m][num_dihedral[m]] = atom3;
          dihedral_atom4[m][num_dihedral[m]] = atom4;
          num_dihedral[m]++;
        }
      }
      if ((m = map(atom3)) >= 0) {
        if (count) count[m]++;
        else {
          dihedral_type[m][num_dihedral[m]] = itype;
          dihedral_atom1[m][num_dihedral[m]] = atom1;
          dihedral_atom2[m][num_dihedral[m]] = atom2;
          dihedral_atom3[m][num_dihedral[m]] = atom3;
          dihedral_atom4[m][num_dihedral[m]] = atom4;
          num_dihedral[m]++;
        }
      }
      if ((m = map(atom4)) >= 0) {
        if (count) count[m]++;
        else {
          dihedral_type[m][num_dihedral[m]] = itype;
          dihedral_atom1[m][num_dihedral[m]] = atom1;
          dihedral_atom2[m][num_dihedral[m]] = atom2;
          dihedral_atom3[m][num_dihedral[m]] = atom3;
          dihedral_atom4[m][num_dihedral[m]] = atom4;
          num_dihedral[m]++;
        }
      }
    }
    buf = next + 1;
  }
}

/* ----------------------------------------------------------------------
   process N impropers read into buf from data files
   if count is non-nullptr, just count impropers per atom
   else store them with atoms
   check that atom IDs are > 0 and <= map_tag_max
------------------------------------------------------------------------- */

void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset,
                          int type_offset)
{
  int m,tmp,itype,rv;
  tagint atom1,atom2,atom3,atom4;
  char *next;
  int newton_bond = force->newton_bond;

  for (int i = 0; i < n; i++) {
    next = strchr(buf,'\n');
    *next = '\0';
    rv = sscanf(buf,"%d %d "
                TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
                &tmp,&itype,&atom1,&atom2,&atom3,&atom4);
    if (rv != 6)
      error->one(FLERR,"Incorrect format of Impropers section in data file");
    if (id_offset) {
      atom1 += id_offset;
      atom2 += id_offset;
      atom3 += id_offset;
      atom4 += id_offset;
    }
    itype += type_offset;

    if ((atom1 <= 0) || (atom1 > map_tag_max) ||
        (atom2 <= 0) || (atom2 > map_tag_max) ||
        (atom3 <= 0) || (atom3 > map_tag_max) ||
        (atom4 <= 0) || (atom4 > map_tag_max) ||
        (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
        (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
      error->one(FLERR,"Invalid atom ID in Impropers section of data file");
    if (itype <= 0 || itype > nimpropertypes)
      error->one(FLERR,
                 "Invalid improper type in Impropers section of data file");
    if ((m = map(atom2)) >= 0) {
      if (count) count[m]++;
      else {
        improper_type[m][num_improper[m]] = itype;
        improper_atom1[m][num_improper[m]] = atom1;
        improper_atom2[m][num_improper[m]] = atom2;
        improper_atom3[m][num_improper[m]] = atom3;
        improper_atom4[m][num_improper[m]] = atom4;
        num_improper[m]++;
      }
    }
    if (newton_bond == 0) {
      if ((m = map(atom1)) >= 0) {
        if (count) count[m]++;
        else {
          improper_type[m][num_improper[m]] = itype;
          improper_atom1[m][num_improper[m]] = atom1;
          improper_atom2[m][num_improper[m]] = atom2;
          improper_atom3[m][num_improper[m]] = atom3;
          improper_atom4[m][num_improper[m]] = atom4;
          num_improper[m]++;
        }
      }
      if ((m = map(atom3)) >= 0) {
        if (count) count[m]++;
        else {
          improper_type[m][num_improper[m]] = itype;
          improper_atom1[m][num_improper[m]] = atom1;
          improper_atom2[m][num_improper[m]] = atom2;
          improper_atom3[m][num_improper[m]] = atom3;
          improper_atom4[m][num_improper[m]] = atom4;
          num_improper[m]++;
        }
      }
      if ((m = map(atom4)) >= 0) {
        if (count) count[m]++;
        else {
          improper_type[m][num_improper[m]] = itype;
          improper_atom1[m][num_improper[m]] = atom1;
          improper_atom2[m][num_improper[m]] = atom2;
          improper_atom3[m][num_improper[m]] = atom3;
          improper_atom4[m][num_improper[m]] = atom4;
          num_improper[m]++;
        }
      }
    }
    buf = next + 1;
  }
}

/* ----------------------------------------------------------------------
   unpack N lines from atom-style specific bonus section of data file
   check that atom IDs are > 0 and <= map_tag_max
   call style-specific routine to parse line
------------------------------------------------------------------------- */

void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset)
{
  int j,m,tagdata;
  char *next;

  next = strchr(buf,'\n');
  *next = '\0';
  int nwords = utils::trim_and_count_words(buf);
  *next = '\n';

  if (nwords != avec_bonus->size_data_bonus)
    error->all(FLERR,"Incorrect bonus data format in data file");

  char **values = new char*[nwords];

  // loop over lines of bonus atom data
  // tokenize the line into values
  // if I own atom tag, unpack its values

  for (int i = 0; i < n; i++) {
    next = strchr(buf,'\n');

    values[0] = strtok(buf," \t\n\r\f");
    for (j = 1; j < nwords; j++)
      values[j] = strtok(nullptr," \t\n\r\f");

    tagdata = ATOTAGINT(values[0]) + id_offset;
    if (tagdata <= 0 || tagdata > map_tag_max)
      error->one(FLERR,"Invalid atom ID in Bonus section of data file");

    // ok to call child's data_atom_bonus() method thru parent avec_bonus,
    // since data_bonus() was called with child ptr, and method is virtual

    if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,&values[1]);

    buf = next + 1;
  }

  delete [] values;
}

/* ----------------------------------------------------------------------
   unpack N bodies from Bodies section of data file
   each body spans multiple lines
   check that atom IDs are > 0 and <= map_tag_max
   call style-specific routine to parse line
------------------------------------------------------------------------- */

void Atom::data_bodies(int n, char *buf, AtomVec *avec_body, tagint id_offset)
{
  int j,m,nvalues,tagdata,ninteger,ndouble;

  int maxint = 0;
  int maxdouble = 0;
  int *ivalues = nullptr;
  double *dvalues = nullptr;

  if (!unique_tags) unique_tags = new std::set<tagint>;

  // loop over lines of body data
  // if I own atom tag, tokenize lines into ivalues/dvalues, call data_body()
  // else skip values

  for (int i = 0; i < n; i++) {
    if (i == 0) tagdata = ATOTAGINT(strtok(buf," \t\n\r\f")) + id_offset;
    else tagdata = ATOTAGINT(strtok(nullptr," \t\n\r\f")) + id_offset;

    if (tagdata <= 0 || tagdata > map_tag_max)
      error->one(FLERR,"Invalid atom ID in Bodies section of data file");

    if (unique_tags->find(tagdata) == unique_tags->end())
      unique_tags->insert(tagdata);
    else
      error->one(FLERR,"Duplicate atom ID in Bodies section of data file");

    ninteger = utils::inumeric(FLERR,strtok(nullptr," \t\n\r\f"),false,lmp);
    ndouble = utils::inumeric(FLERR,strtok(nullptr," \t\n\r\f"),false,lmp);

    if ((m = map(tagdata)) >= 0) {
      if (ninteger > maxint) {
        delete [] ivalues;
        maxint = ninteger;
        ivalues = new int[maxint];
      }
      if (ndouble > maxdouble) {
        delete [] dvalues;
        maxdouble = ndouble;
        dvalues = new double[maxdouble];
      }

      for (j = 0; j < ninteger; j++)
        ivalues[j] = utils::inumeric(FLERR,strtok(nullptr," \t\n\r\f"),false,lmp);
      for (j = 0; j < ndouble; j++)
        dvalues[j] = utils::numeric(FLERR,strtok(nullptr," \t\n\r\f"),false,lmp);

      avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues);

    } else {
      nvalues = ninteger + ndouble;    // number of values to skip
      for (j = 0; j < nvalues; j++)
        strtok(nullptr," \t\n\r\f");
    }
  }

  delete [] ivalues;
  delete [] dvalues;
}

/* ----------------------------------------------------------------------
   init per-atom fix/compute/variable values for newly created atoms
   called from create_atoms, read_data, read_dump,
     lib::lammps_create_atoms()
   fixes, computes, variables may or may not exist when called
------------------------------------------------------------------------- */

void Atom::data_fix_compute_variable(int nprev, int nnew)
{
  for (int m = 0; m < modify->nfix; m++) {
    Fix *fix = modify->fix[m];
    if (fix->create_attribute)
      for (int i = nprev; i < nnew; i++)
        fix->set_arrays(i);
  }

  for (int m = 0; m < modify->ncompute; m++) {
    Compute *compute = modify->compute[m];
    if (compute->create_attribute)
      for (int i = nprev; i < nnew; i++)
        compute->set_arrays(i);
  }

  for (int i = nprev; i < nnew; i++)
    input->variable->set_arrays(i);
}

/* ----------------------------------------------------------------------
   allocate arrays of length ntypes
   only done after ntypes is set
------------------------------------------------------------------------- */

void Atom::allocate_type_arrays()
{
  if (avec->mass_type == AtomVec::PER_TYPE) {
    mass = new double[ntypes+1];
    mass_setflag = new int[ntypes+1];
    for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0;
  }
}

/* ----------------------------------------------------------------------
   set a mass and flag it as set
   called from reading of data file
   type_offset may be used when reading multiple data files
------------------------------------------------------------------------- */

void Atom::set_mass(const char *file, int line, const char *str, int type_offset)
{
  if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style");

  int itype;
  double mass_one;
  int n = sscanf(str,"%d %lg",&itype,&mass_one);
  if (n != 2) error->all(file,line,"Invalid mass line in data file");
  itype += type_offset;

  if (itype < 1 || itype > ntypes)
    error->all(file,line,"Invalid type for mass set");

  mass[itype] = mass_one;
  mass_setflag[itype] = 1;

  if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value");
}

/* ----------------------------------------------------------------------
   set a mass and flag it as set
   called from EAM pair routine
------------------------------------------------------------------------- */

void Atom::set_mass(const char *file, int line, int itype, double value)
{
  if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style");
  if (itype < 1 || itype > ntypes)
    error->all(file,line,"Invalid type for mass set");

  mass[itype] = value;
  mass_setflag[itype] = 1;

  if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value");
}

/* ----------------------------------------------------------------------
   set one or more masses and flag them as set
   called from reading of input script
------------------------------------------------------------------------- */

void Atom::set_mass(const char *file, int line, int /*narg*/, char **arg)
{
  if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style");

  int lo,hi;
  utils::bounds(file,line,arg[0],1,ntypes,lo,hi,error);
  if (lo < 1 || hi > ntypes) error->all(file,line,"Invalid type for mass set");

  for (int itype = lo; itype <= hi; itype++) {
    mass[itype] = atof(arg[1]);
    mass_setflag[itype] = 1;

    if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value");
  }
}

/* ----------------------------------------------------------------------
   set all masses
   called from reading of restart file, also from ServerMD
------------------------------------------------------------------------- */

void Atom::set_mass(double *values)
{
  for (int itype = 1; itype <= ntypes; itype++) {
    mass[itype] = values[itype];
    mass_setflag[itype] = 1;
  }
}

/* ----------------------------------------------------------------------
   check that all per-atom-type masses have been set
------------------------------------------------------------------------- */

void Atom::check_mass(const char *file, int line)
{
  if (mass == nullptr) return;
  if (rmass_flag) return;
  for (int itype = 1; itype <= ntypes; itype++)
    if (mass_setflag[itype] == 0)
      error->all(file,line,"Not all per-type masses are set");
}

/* ----------------------------------------------------------------------
   check that radii of all particles of itype are the same
   return 1 if true, else return 0
   also return the radius value for that type
------------------------------------------------------------------------- */

int Atom::radius_consistency(int itype, double &rad)
{
  double value = -1.0;
  int flag = 0;
  for (int i = 0; i < nlocal; i++) {
    if (type[i] != itype) continue;
    if (value < 0.0) value = radius[i];
    else if (value != radius[i]) flag = 1;
  }

  int flagall;
  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
  if (flagall) return 0;

  MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world);
  return 1;
}

/* ----------------------------------------------------------------------
   check that shape of all particles of itype are the same
   return 1 if true, else return 0
   also return the 3 shape params for itype
------------------------------------------------------------------------- */

int Atom::shape_consistency(int itype,
                            double &shapex, double &shapey, double &shapez)
{
  double zero[3] = {0.0, 0.0, 0.0};
  double one[3] = {-1.0, -1.0, -1.0};
  double *shape;

  auto avec_ellipsoid = (AtomVecEllipsoid *) style_match("ellipsoid");
  auto bonus = avec_ellipsoid->bonus;

  int flag = 0;
  for (int i = 0; i < nlocal; i++) {
    if (type[i] != itype) continue;
    if (ellipsoid[i] < 0) shape = zero;
    else shape = bonus[ellipsoid[i]].shape;

    if (one[0] < 0.0) {
      one[0] = shape[0];
      one[1] = shape[1];
      one[2] = shape[2];
    } else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2])
      flag = 1;
  }

  int flagall;
  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
  if (flagall) return 0;

  double oneall[3];
  MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world);
  shapex = oneall[0];
  shapey = oneall[1];
  shapez = oneall[2];
  return 1;
}

/* ----------------------------------------------------------------------
   add a new molecule template = set of molecules
------------------------------------------------------------------------- */

void Atom::add_molecule(int narg, char **arg)
{
  if (narg < 1) error->all(FLERR,"Illegal molecule command");

  if (find_molecule(arg[0]) >= 0)
    error->all(FLERR,"Reuse of molecule template ID");

  // 1st molecule in set stores nset = # of mols, others store nset = 0
  // ifile = count of molecules in set
  // index = argument index where next molecule starts, updated by constructor

  int ifile = 1;
  int index = 1;
  while (1) {
    molecules = (Molecule **)
      memory->srealloc(molecules,(nmolecule+1)*sizeof(Molecule *),
                       "atom::molecules");
    molecules[nmolecule] = new Molecule(lmp,narg,arg,index);
    molecules[nmolecule]->nset = 0;
    molecules[nmolecule-ifile+1]->nset++;
    nmolecule++;
    if (molecules[nmolecule-1]->last) break;
    ifile++;
  }
}

/* ----------------------------------------------------------------------
   find first molecule in set with template ID
   return -1 if does not exist
------------------------------------------------------------------------- */

int Atom::find_molecule(char *id)
{
  if(id == nullptr) return -1;
  int imol;
  for (imol = 0; imol < nmolecule; imol++)
    if (strcmp(id,molecules[imol]->id) == 0) return imol;
  return -1;
}

/* ----------------------------------------------------------------------
   add info to current atom ilocal from molecule template onemol and its iatom
   offset = atom ID preceding IDs of atoms in this molecule
   called by fixes and commands that add molecules
------------------------------------------------------------------------- */

void Atom::add_molecule_atom(Molecule *onemol, int iatom,
                             int ilocal, tagint offset)
{
  if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom];
  if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom];
  if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom];
  else if (rmass_flag)
    rmass[ilocal] = 4.0*MY_PI/3.0 *
      radius[ilocal]*radius[ilocal]*radius[ilocal];
  if (onemol->bodyflag) {
    body[ilocal] = 0;     // as if a body read from data file
    onemol->avec_body->data_body(ilocal,onemol->nibody,onemol->ndbody,
                                 onemol->ibodyparams,onemol->dbodyparams);
    onemol->avec_body->set_quat(ilocal,onemol->quat_external);
  }

  if (molecular != 1) return;

  // add bond topology info
  // for molecular atom styles, but not atom style template

  if (avec->bonds_allow) {
    num_bond[ilocal] = onemol->num_bond[iatom];
    for (int i = 0; i < num_bond[ilocal]; i++) {
      bond_type[ilocal][i] = onemol->bond_type[iatom][i];
      bond_atom[ilocal][i] = onemol->bond_atom[iatom][i] + offset;
    }
  }

  if (avec->angles_allow) {
    num_angle[ilocal] = onemol->num_angle[iatom];
    for (int i = 0; i < num_angle[ilocal]; i++) {
      angle_type[ilocal][i] = onemol->angle_type[iatom][i];
      angle_atom1[ilocal][i] = onemol->angle_atom1[iatom][i] + offset;
      angle_atom2[ilocal][i] = onemol->angle_atom2[iatom][i] + offset;
      angle_atom3[ilocal][i] = onemol->angle_atom3[iatom][i] + offset;
    }
  }

  if (avec->dihedrals_allow) {
    num_dihedral[ilocal] = onemol->num_dihedral[iatom];
    for (int i = 0; i < num_dihedral[ilocal]; i++) {
      dihedral_type[ilocal][i] = onemol->dihedral_type[iatom][i];
      dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset;
      dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset;
      dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset;
      dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset;
    }
  }

  if (avec->impropers_allow) {
    num_improper[ilocal] = onemol->num_improper[iatom];
    for (int i = 0; i < num_improper[ilocal]; i++) {
      improper_type[ilocal][i] = onemol->improper_type[iatom][i];
      improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset;
      improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset;
      improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset;
      improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset;
    }
  }

  if (onemol->specialflag) {
    nspecial[ilocal][0] = onemol->nspecial[iatom][0];
    nspecial[ilocal][1] = onemol->nspecial[iatom][1];
    int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2];
    for (int i = 0; i < n; i++)
      special[ilocal][i] = onemol->special[iatom][i] + offset;
  }
}

/* ----------------------------------------------------------------------
   reorder owned atoms so those in firstgroup appear first
   called by comm->exchange() if atom_modify first group is set
   only owned atoms exist at this point, no ghost atoms
------------------------------------------------------------------------- */

void Atom::first_reorder()
{
  // insure there is one extra atom location at end of arrays for swaps

  if (nlocal == nmax) avec->grow(0);

  // loop over owned atoms
  // nfirst = index of first atom not in firstgroup
  // when find firstgroup atom out of place, swap it with atom nfirst

  int bitmask = group->bitmask[firstgroup];
  nfirst = 0;
  while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;

  for (int i = 0; i < nlocal; i++) {
    if (mask[i] & bitmask && i > nfirst) {
      avec->copy(i,nlocal,0);
      avec->copy(nfirst,i,0);
      avec->copy(nlocal,nfirst,0);
      while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
    }
  }
}

/* ----------------------------------------------------------------------
   perform spatial sort of atoms within my sub-domain
   always called between comm->exchange() and comm->borders()
   don't have to worry about clearing/setting atom->map since done in comm
------------------------------------------------------------------------- */

void Atom::sort()
{
  int i,m,n,ix,iy,iz,ibin,empty;

  // set next timestep for sorting to take place

  nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq;

  // re-setup sort bins if needed

  if (domain->box_change) setup_sort_bins();
  if (nbins == 1) return;

  // reallocate per-atom vectors if needed

  if (nlocal > maxnext) {
    memory->destroy(next);
    memory->destroy(permute);
    maxnext = atom->nmax;
    memory->create(next,maxnext,"atom:next");
    memory->create(permute,maxnext,"atom:permute");
  }

  // insure there is one extra atom location at end of arrays for swaps

  if (nlocal == nmax) avec->grow(0);

  // bin atoms in reverse order so linked list will be in forward order

  for (i = 0; i < nbins; i++) binhead[i] = -1;

  for (i = nlocal-1; i >= 0; i--) {
    ix = static_cast<int> ((x[i][0]-bboxlo[0])*bininvx);
    iy = static_cast<int> ((x[i][1]-bboxlo[1])*bininvy);
    iz = static_cast<int> ((x[i][2]-bboxlo[2])*bininvz);
    ix = MAX(ix,0);
    iy = MAX(iy,0);
    iz = MAX(iz,0);
    ix = MIN(ix,nbinx-1);
    iy = MIN(iy,nbiny-1);
    iz = MIN(iz,nbinz-1);
    ibin = iz*nbiny*nbinx + iy*nbinx + ix;
    next[i] = binhead[ibin];
    binhead[ibin] = i;
  }

  // permute = desired permutation of atoms
  // permute[I] = J means Ith new atom will be Jth old atom

  n = 0;
  for (m = 0; m < nbins; m++) {
    i = binhead[m];
    while (i >= 0) {
      permute[n++] = i;
      i = next[i];
    }
  }

  // current = current permutation, just reuse next vector
  // current[I] = J means Ith current atom is Jth old atom

  int *current = next;
  for (i = 0; i < nlocal; i++) current[i] = i;

  // reorder local atom list, when done, current = permute
  // perform "in place" using copy() to extra atom location at end of list
  // inner while loop processes one cycle of the permutation
  // copy before inner-loop moves an atom to end of atom list
  // copy after inner-loop moves atom at end of list back into list
  // empty = location in atom list that is currently empty

  for (i = 0; i < nlocal; i++) {
    if (current[i] == permute[i]) continue;
    avec->copy(i,nlocal,0);
    empty = i;
    while (permute[empty] != i) {
      avec->copy(permute[empty],empty,0);
      empty = current[empty] = permute[empty];
    }
    avec->copy(nlocal,empty,0);
    current[empty] = permute[empty];
  }

  // sanity check that current = permute

  //int flag = 0;
  //for (i = 0; i < nlocal; i++)
  //  if (current[i] != permute[i]) flag = 1;
  //int flagall;
  //MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
  //if (flagall) error->all(FLERR,"Atom sort did not operate correctly");
}

/* ----------------------------------------------------------------------
   setup bins for spatial sorting of atoms
------------------------------------------------------------------------- */

void Atom::setup_sort_bins()
{
  // binsize:
  // user setting if explicitly set
  // default = 1/2 of neighbor cutoff
  // check if neighbor cutoff = 0.0
  // and in that case, disable sorting

  double binsize = 0.0;
  if (userbinsize > 0.0) binsize = userbinsize;
  else if (neighbor->cutneighmax > 0.0) binsize = 0.5 * neighbor->cutneighmax;

  if ((binsize == 0.0) && (sortfreq > 0)) {
    sortfreq = 0;
    if (comm->me == 0)
          error->warning(FLERR,"No pairwise cutoff or binsize set. "
                         "Atom sorting therefore disabled.");
    return;
  }

  double bininv = 1.0/binsize;

  // nbin xyz = local bins
  // bbox lo/hi = bounding box of my sub-domain

  if (domain->triclinic)
    domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi);
  else {
    bboxlo[0] = domain->sublo[0];
    bboxlo[1] = domain->sublo[1];
    bboxlo[2] = domain->sublo[2];
    bboxhi[0] = domain->subhi[0];
    bboxhi[1] = domain->subhi[1];
    bboxhi[2] = domain->subhi[2];
  }

  nbinx = static_cast<int> ((bboxhi[0]-bboxlo[0]) * bininv);
  nbiny = static_cast<int> ((bboxhi[1]-bboxlo[1]) * bininv);
  nbinz = static_cast<int> ((bboxhi[2]-bboxlo[2]) * bininv);
  if (domain->dimension == 2) nbinz = 1;

  if (nbinx == 0) nbinx = 1;
  if (nbiny == 0) nbiny = 1;
  if (nbinz == 0) nbinz = 1;

  bininvx = nbinx / (bboxhi[0]-bboxlo[0]);
  bininvy = nbiny / (bboxhi[1]-bboxlo[1]);
  bininvz = nbinz / (bboxhi[2]-bboxlo[2]);

  #ifdef LMP_USER_INTEL
  int intel_neigh = 0;
  if (neighbor->nrequest) {
    if (neighbor->requests[0]->intel) intel_neigh = 1;
  } else if (neighbor->old_nrequest)
    if (neighbor->old_requests[0]->intel) intel_neigh = 1;
  if (intel_neigh && userbinsize == 0.0) {
    if (neighbor->binsizeflag) bininv = 1.0/neighbor->binsize_user;

    double nx_low = neighbor->bboxlo[0];
    double ny_low = neighbor->bboxlo[1];
    double nz_low = neighbor->bboxlo[2];
    double nxbbox = neighbor->bboxhi[0] - nx_low;
    double nybbox = neighbor->bboxhi[1] - ny_low;
    double nzbbox = neighbor->bboxhi[2] - nz_low;
    int nnbinx = static_cast<int> (nxbbox * bininv);
    int nnbiny = static_cast<int> (nybbox * bininv);
    int nnbinz = static_cast<int> (nzbbox * bininv);
    if (domain->dimension == 2) nnbinz = 1;

    if (nnbinx == 0) nnbinx = 1;
    if (nnbiny == 0) nnbiny = 1;
    if (nnbinz == 0) nnbinz = 1;

    double binsizex = nxbbox/nnbinx;
    double binsizey = nybbox/nnbiny;
    double binsizez = nzbbox/nnbinz;

    bininvx = 1.0 / binsizex;
    bininvy = 1.0 / binsizey;
    bininvz = 1.0 / binsizez;

    int lxo = (bboxlo[0] - nx_low) * bininvx;
    int lyo = (bboxlo[1] - ny_low) * bininvy;
    int lzo = (bboxlo[2] - nz_low) * bininvz;
    bboxlo[0] = nx_low + static_cast<double>(lxo) / bininvx;
    bboxlo[1] = ny_low + static_cast<double>(lyo) / bininvy;
    bboxlo[2] = nz_low + static_cast<double>(lzo) / bininvz;
    nbinx = static_cast<int>((bboxhi[0] - bboxlo[0]) * bininvx) + 1;
    nbiny = static_cast<int>((bboxhi[1] - bboxlo[1]) * bininvy) + 1;
    nbinz = static_cast<int>((bboxhi[2] - bboxlo[2]) * bininvz) + 1;
    bboxhi[0] = bboxlo[0] + static_cast<double>(nbinx) / bininvx;
    bboxhi[1] = bboxlo[1] + static_cast<double>(nbiny) / bininvy;
    bboxhi[2] = bboxlo[2] + static_cast<double>(nbinz) / bininvz;
  }
  #endif

  if (1.0*nbinx*nbiny*nbinz > INT_MAX)
    error->one(FLERR,"Too many atom sorting bins");

  nbins = nbinx*nbiny*nbinz;

  // reallocate per-bin memory if needed

  if (nbins > maxbin) {
    memory->destroy(binhead);
    maxbin = nbins;
    memory->create(binhead,maxbin,"atom:binhead");
  }
}

/* ----------------------------------------------------------------------
   register a callback to a fix so it can manage atom-based arrays
   happens when fix is created
   flag = 0 for grow, 1 for restart, 2 for border comm
------------------------------------------------------------------------- */

void Atom::add_callback(int flag)
{
  int ifix;

  // find the fix
  // if find null pointer:
  //   it's this one, since it is being replaced and has just been deleted
  //   at this point in re-creation
  // if don't find null pointer:
  //   i is set to nfix = new one currently being added at end of list

  for (ifix = 0; ifix < modify->nfix; ifix++)
    if (modify->fix[ifix] == nullptr) break;

  // add callback to lists and sort, reallocating if necessary
  // sorting is required in cases where fixes were replaced as it ensures atom
  // data is read/written/transfered in the same order that fixes are called

  if (flag == GROW) {
    if (nextra_grow == nextra_grow_max) {
      nextra_grow_max += DELTA;
      memory->grow(extra_grow,nextra_grow_max,"atom:extra_grow");
    }
    extra_grow[nextra_grow] = ifix;
    nextra_grow++;
    std::sort(extra_grow, extra_grow + nextra_grow);
  } else if (flag == RESTART) {
    if (nextra_restart == nextra_restart_max) {
      nextra_restart_max += DELTA;
      memory->grow(extra_restart,nextra_restart_max,"atom:extra_restart");
    }
    extra_restart[nextra_restart] = ifix;
    nextra_restart++;
    std::sort(extra_restart, extra_restart + nextra_restart);
  } else if (flag == BORDER) {
    if (nextra_border == nextra_border_max) {
      nextra_border_max += DELTA;
      memory->grow(extra_border,nextra_border_max,"atom:extra_border");
    }
    extra_border[nextra_border] = ifix;
    nextra_border++;
    std::sort(extra_border, extra_border + nextra_border);
  }
}

/* ----------------------------------------------------------------------
   unregister a callback to a fix
   happens when fix is deleted, called by its destructor
   flag = 0 for grow, 1 for restart
------------------------------------------------------------------------- */

void Atom::delete_callback(const char *id, int flag)
{
  if (id == nullptr) return;

  int ifix = modify->find_fix(id);

  // compact the list of callbacks

  if (flag == GROW) {
    int match;
    for (match = 0; match < nextra_grow; match++)
      if (extra_grow[match] == ifix) break;
    if ((nextra_grow == 0) || (match == nextra_grow))
      error->all(FLERR,"Trying to delete non-existent Atom::grow() callback");
    for (int i = match; i < nextra_grow-1; i++)
      extra_grow[i] = extra_grow[i+1];
    nextra_grow--;

  } else if (flag == RESTART) {
    int match;
    for (match = 0; match < nextra_restart; match++)
      if (extra_restart[match] == ifix) break;
    if ((nextra_restart == 0) || (match == nextra_restart))
      error->all(FLERR,"Trying to delete non-existent Atom::restart() callback");
    for (int i = match; i < nextra_restart-1; i++)
      extra_restart[i] = extra_restart[i+1];
    nextra_restart--;

  } else if (flag == BORDER) {
    int match;
    for (match = 0; match < nextra_border; match++)
      if (extra_border[match] == ifix) break;
    if ((nextra_border == 0) || (match == nextra_border))
      error->all(FLERR,"Trying to delete non-existent Atom::border() callback");
    for (int i = match; i < nextra_border-1; i++)
      extra_border[i] = extra_border[i+1];
    nextra_border--;
  }
}

/* ----------------------------------------------------------------------
   decrement ptrs in callback lists to fixes beyond the deleted ifix
   happens after fix is deleted
------------------------------------------------------------------------- */

void Atom::update_callback(int ifix)
{
  for (int i = 0; i < nextra_grow; i++)
    if (extra_grow[i] > ifix) extra_grow[i]--;
  for (int i = 0; i < nextra_restart; i++)
    if (extra_restart[i] > ifix) extra_restart[i]--;
  for (int i = 0; i < nextra_border; i++)
    if (extra_border[i] > ifix) extra_border[i]--;
}

/* ----------------------------------------------------------------------
   find custom per-atom vector with name
   return index if found, and flag = 0/1 for int/double
   return -1 if not found
------------------------------------------------------------------------- */

int Atom::find_custom(const char *name, int &flag)
{
  if(name == nullptr) return -1;

  for (int i = 0; i < nivector; i++)
    if (iname[i] && strcmp(iname[i],name) == 0) {
      flag = 0;
      return i;
    }

  for (int i = 0; i < ndvector; i++)
    if (dname[i] && strcmp(dname[i],name) == 0) {
      flag = 1;
      return i;
    }

  return -1;
}

/** \brief Add a custom per-atom property with the given name and type
\verbatim embed:rst

This function will add a custom per-atom property with the name "name"
as either list of int or double to the list of custom properties.  This
function is called, e.g. from :doc:`fix property/atom <fix_property_atom>`.
\endverbatim
 * \param name Name of the property (w/o a "d_" or "i_" prefix)
 * \param flag Data type of property: 0 for int, 1 for double
 * \return Index of property in the respective list of properties
 */
int Atom::add_custom(const char *name, int flag)
{
  int index;

  if (flag == 0) {
    index = nivector;
    nivector++;
    iname = (char **) memory->srealloc(iname,nivector*sizeof(char *),
                                       "atom:iname");
    int n = strlen(name) + 1;
    iname[index] = new char[n];
    strcpy(iname[index],name);
    ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *),
                                        "atom:ivector");
    memory->create(ivector[index],nmax,"atom:ivector");
  } else {
    index = ndvector;
    ndvector++;
    dname = (char **) memory->srealloc(dname,ndvector*sizeof(char *),
                                       "atom:dname");
    int n = strlen(name) + 1;
    dname[index] = new char[n];
    strcpy(dname[index],name);
    dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *),
                                           "atom:dvector");
    memory->create(dvector[index],nmax,"atom:dvector");
  }

  return index;
}

/*! \brief Remove a custom per-atom property of a given type
 *
\verbatim embed:rst
This will remove a property that was requested e.g. by the
:doc:`fix property/atom <fix_property_atom>` command.  It frees the
allocated memory and sets the pointer to ``nullptr`` to the entry in
the list can be reused. The lists of those pointers will never be
compacted or never shrink, so that index to name mappings remain valid.
\endverbatim
 *
 * \param flag whether the property is integer (=0) or double (=1)
 * \param index of that property in the respective list.
 */
void Atom::remove_custom(int flag, int index)
{
  if (flag == 0) {
    memory->destroy(ivector[index]);
    ivector[index] = nullptr;
    delete [] iname[index];
    iname[index] = nullptr;
  } else {
    memory->destroy(dvector[index]);
    dvector[index] = nullptr;
    delete [] dname[index];
    dname[index] = nullptr;
  }
}

/** Provide access to internal data of the Atom class by keyword
 *
\verbatim embed:rst

This function is a way to access internal per-atom data.  This data is
distributed across MPI ranks and thus only the data for "local" atoms
are expected to be available.  Whether also data for "ghost" atoms is
stored and up-to-date depends on various simulation settings.

This table lists a large part of the supported names, their data types,
length of the data area, and a short description.

.. list-table::
   :header-rows: 1
   :widths: auto

   * - Name
     - Type
     - Items per atom
     - Description
   * - mass
     - double
     - 1
     - per-type mass. This array is **NOT** a per-atom array
       but of length ``ntypes+1``, element 0 is ignored.
   * - id
     - tagint
     - 1
     - atom ID of the particles
   * - type
     - int
     - 1
     - atom type of the particles
   * - mask
     - int
     - 1
     - bitmask for mapping to groups. Individual bits are set
       to 0 or 1 for each group.
   * - image
     - imageint
     - 1
     - 3 image flags encoded into a single integer.
       See :cpp:func:`lammps_encode_image_flags`.
   * - x
     - double
     - 3
     - x-, y-, and z-coordinate of the particles
   * - v
     - double
     - 3
     - x-, y-, and z-component of the velocity of the particles
   * - f
     - double
     - 3
     - x-, y-, and z-component of the force on the particles
   * - molecule
     - int
     - 1
     - molecule ID of the particles
   * - q
     - double
     - 1
     - charge of the particles
   * - mu
     - double
     - 3
     - dipole moment of the particles
   * - omega
     - double
     - 3
     - x-, y-, and z-component of rotational velocity of the particles
   * - angmom
     - double
     - 3
     - x-, y-, and z-component of angular momentum of the particles
   * - torque
     - double
     - 3
     - x-, y-, and z-component of the torque on the particles
   * - radius
     - double
     - 1
     - radius of the (extended) particles
   * - rmass
     - double
     - 1
     - per-atom mass of the particles. ``nullptr`` if per-type masses are
       used. See the :cpp:func:`rmass_flag<lammps_extract_setting>` setting.
   * - ellipsoid
     - int
     - 1
     - 1 if the particle is an ellipsoidal particle, 0 if not
   * - line
     - int
     - 1
     - 1 if the particle is a line particle, 0 if not
   * - tri
     - int
     - 1
     - 1 if the particle is a triangulated particle, 0 if not
   * - body
     - int
     - 1
     - 1 if the particle is a body particle, 0 if not

\endverbatim
 *
 * \sa extract_datatype
 *
 * \param  name  string with the keyword of the desired property.
                 Typically the name of the pointer variable returned
 * \return       pointer to the requested data cast to ``void *`` or ``nullptr`` */

void *Atom::extract(const char *name)
{
  // --------------------------------------------------------------------
  // 4th customization section: customize by adding new variable name
  // please see the following function to set the type of the data
  // so that programs can detect it dynamically at run time.

  /* NOTE: this array is only of length ntypes+1 */
  if (strcmp(name,"mass") == 0) return (void *) mass;

  if (strcmp(name,"id") == 0) return (void *) tag;
  if (strcmp(name,"type") == 0) return (void *) type;
  if (strcmp(name,"mask") == 0) return (void *) mask;
  if (strcmp(name,"image") == 0) return (void *) image;
  if (strcmp(name,"x") == 0) return (void *) x;
  if (strcmp(name,"v") == 0) return (void *) v;
  if (strcmp(name,"f") == 0) return (void *) f;
  if (strcmp(name,"molecule") == 0) return (void *) molecule;
  if (strcmp(name,"q") == 0) return (void *) q;
  if (strcmp(name,"mu") == 0) return (void *) mu;
  if (strcmp(name,"omega") == 0) return (void *) omega;
  if (strcmp(name,"angmom") == 0) return (void *) angmom;
  if (strcmp(name,"torque") == 0) return (void *) torque;
  if (strcmp(name,"radius") == 0) return (void *) radius;
  if (strcmp(name,"rmass") == 0) return (void *) rmass;
  if (strcmp(name,"ellipsoid") == 0) return (void *) ellipsoid;
  if (strcmp(name,"line") == 0) return (void *) line;
  if (strcmp(name,"tri") == 0) return (void *) tri;
  if (strcmp(name,"body") == 0) return (void *) body;

  if (strcmp(name,"vfrac") == 0) return (void *) vfrac;
  if (strcmp(name,"s0") == 0) return (void *) s0;
  if (strcmp(name,"x0") == 0) return (void *) x0;

  if (strcmp(name,"spin") == 0) return (void *) spin;
  if (strcmp(name,"eradius") == 0) return (void *) eradius;
  if (strcmp(name,"ervel") == 0) return (void *) ervel;
  if (strcmp(name,"erforce") == 0) return (void *) erforce;
  if (strcmp(name,"ervelforce") == 0) return (void *) ervelforce;
  if (strcmp(name,"cs") == 0) return (void *) cs;
  if (strcmp(name,"csforce") == 0) return (void *) csforce;
  if (strcmp(name,"vforce") == 0) return (void *) vforce;
  if (strcmp(name,"etag") == 0) return (void *) etag;

  if (strcmp(name,"rho") == 0) return (void *) rho;
  if (strcmp(name,"drho") == 0) return (void *) drho;
  if (strcmp(name,"esph") == 0) return (void *) esph;
  if (strcmp(name,"desph") == 0) return (void *) desph;
  if (strcmp(name,"cv") == 0) return (void *) cv;
  if (strcmp(name,"vest") == 0) return (void *) vest;

  // USER-MESONT package
  if (strcmp(name,"length") == 0) return (void *) length;
  if (strcmp(name,"buckling") == 0) return (void *) buckling;
  if (strcmp(name,"bond_nt") == 0) return (void *) bond_nt;

  if (strcmp(name, "contact_radius") == 0) return (void *) contact_radius;
  if (strcmp(name, "smd_data_9") == 0) return (void *) smd_data_9;
  if (strcmp(name, "smd_stress") == 0) return (void *) smd_stress;
  if (strcmp(name, "eff_plastic_strain") == 0)
    return (void *) eff_plastic_strain;
  if (strcmp(name, "eff_plastic_strain_rate") == 0)
    return (void *) eff_plastic_strain_rate;
  if (strcmp(name, "damage") == 0) return (void *) damage;

  if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta;
  if (strcmp(name,"edpd_temp") == 0) return (void *) edpd_temp;

  // end of customization section
  // --------------------------------------------------------------------

  return nullptr;
}


/** Provide data type info about internal data of the Atom class
 *
\verbatim embed:rst

.. versionadded:: 18Sep2020

\endverbatim
 *
 * \sa extract
 *
 * \param  name  string with the keyword of the desired property.
 * \return       data type constant for desired property or -1 */

int Atom::extract_datatype(const char *name)
{
  // --------------------------------------------------------------------
  // 5th customization section: customize by adding new variable name

  if (strcmp(name,"mass") == 0) return LAMMPS_DOUBLE;

  if (strcmp(name,"id") == 0) return LAMMPS_TAGINT;
  if (strcmp(name,"type") == 0) return LAMMPS_INT;
  if (strcmp(name,"mask") == 0) return LAMMPS_INT;
  if (strcmp(name,"image") == 0) return LAMMPS_TAGINT;
  if (strcmp(name,"x") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"v") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"f") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"molecule") == 0) return LAMMPS_TAGINT;
  if (strcmp(name,"q") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"mu") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"omega") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"angmom") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"torque") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"radius") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"rmass") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"ellipsoid") == 0) return LAMMPS_INT;
  if (strcmp(name,"line") == 0) return LAMMPS_INT;
  if (strcmp(name,"tri") == 0) return LAMMPS_INT;
  if (strcmp(name,"body") == 0) return LAMMPS_INT;

  if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"x0") == 0) return LAMMPS_DOUBLE_2D;

  if (strcmp(name,"spin") == 0) return LAMMPS_INT;
  if (strcmp(name,"eradius") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"ervel") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"erforce") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"ervelforce") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"cs") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"csforce") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"vforce") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name,"etag") == 0) return LAMMPS_INT;

  if (strcmp(name,"rho") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"drho") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"esph") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"desph") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"cv") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"vest") == 0) return LAMMPS_DOUBLE_2D;

  // USER-MESONT package
  if (strcmp(name,"length") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"buckling") == 0) return LAMMPS_INT;
  if (strcmp(name,"bond_nt") == 0) return  LAMMPS_TAGINT_2D;

  if (strcmp(name, "contact_radius") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name, "smd_data_9") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name, "smd_stress") == 0) return LAMMPS_DOUBLE_2D;
  if (strcmp(name, "eff_plastic_strain") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name, "eff_plastic_strain_rate") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name, "damage") == 0) return LAMMPS_DOUBLE;

  if (strcmp(name,"dpdTheta") == 0) return LAMMPS_DOUBLE;
  if (strcmp(name,"edpd_temp") == 0) return LAMMPS_DOUBLE;

  // end of customization section
  // --------------------------------------------------------------------

  return -1;
}

/* ----------------------------------------------------------------------
   return # of bytes of allocated memory
   call to avec tallies per-atom vectors
   add in global to local mapping storage
------------------------------------------------------------------------- */

double Atom::memory_usage()
{
  double bytes = avec->memory_usage();

  bytes += max_same*sizeof(int);
  if (map_style == MAP_ARRAY)
    bytes += memory->usage(map_array,map_maxarray);
  else if (map_style == MAP_HASH) {
    bytes += map_nbucket*sizeof(int);
    bytes += map_nhash*sizeof(HashElem);
  }
  if (maxnext) {
    bytes += memory->usage(next,maxnext);
    bytes += memory->usage(permute,maxnext);
  }

  return bytes;
}