/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   https://www.lammps.org/, 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 "fix_pour.h"

#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix_gravity.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "modify.h"
#include "molecule.h"
#include "random_park.h"
#include "region.h"
#include "region_block.h"
#include "region_cylinder.h"
#include "update.h"

#include <cmath>
#include <cstring>

using namespace LAMMPS_NS;
using namespace FixConst;
using MathConst::MY_2PI;
using MathConst::MY_4PI3;
using MathConst::MY_PI;

enum { ATOM, MOLECULE };
enum { ONE, RANGE, POLY };

static constexpr double EPSILON = 0.001;
static constexpr double SMALL = 1.0e-10;

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

FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
    Fix(lmp, narg, arg), radius_poly(nullptr), frac_poly(nullptr), idrigid(nullptr),
    idshake(nullptr), idregion(nullptr), region(nullptr), onemols(nullptr), molfrac(nullptr),
    coords(nullptr), imageflags(nullptr), fixrigid(nullptr), fixshake(nullptr), recvcounts(nullptr),
    displs(nullptr), random(nullptr), random2(nullptr)
{
  if (narg < 6) error->all(FLERR, "Illegal fix pour command");

  if (lmp->kokkos) error->all(FLERR, "Cannot yet use fix pour with the KOKKOS package");

  time_depend = 1;

  if (!atom->radius_flag || !atom->rmass_flag)
    error->all(FLERR, "Fix pour requires atom attributes radius, rmass");

  // required args

  ninsert = utils::inumeric(FLERR, arg[3], false, lmp);
  ntype = utils::inumeric(FLERR, arg[4], false, lmp);
  seed = utils::inumeric(FLERR, arg[5], false, lmp);

  if (seed <= 0) error->all(FLERR, "Illegal fix pour command");

  // read options from end of input line

  options(narg - 6, &arg[6]);

  // error check on type

  if (mode == ATOM && (ntype <= 0 || ntype > atom->ntypes))
    error->all(FLERR, "Invalid atom type in fix pour command");

  // error checks on region and its extent being inside simulation box

  if (!region) error->all(FLERR, "Must specify a region in fix pour");
  if (region->bboxflag == 0)
    error->all(FLERR, "Fix pour region {} does not support a bounding box", idregion);
  if (region->dynamic_check()) error->all(FLERR, "Fix pour region {} cannot be dynamic", idregion);

  if (strcmp(region->style, "block") == 0) {
    auto block = dynamic_cast<RegBlock *>(region);
    region_style = 1;
    xlo = block->xlo;
    xhi = block->xhi;
    ylo = block->ylo;
    yhi = block->yhi;
    zlo = block->zlo;
    zhi = block->zhi;
    if (xlo < domain->boxlo[0] || xhi > domain->boxhi[0] || ylo < domain->boxlo[1] ||
        yhi > domain->boxhi[1] || zlo < domain->boxlo[2] || zhi > domain->boxhi[2])
      error->all(FLERR, "Insertion region extends outside simulation box");
  } else if (strcmp(region->style, "cylinder") == 0) {
    auto cylinder = dynamic_cast<RegCylinder *>(region);
    region_style = 2;
    char axis = cylinder->axis;
    xc = cylinder->c1;
    yc = cylinder->c2;
    rc = cylinder->radius;
    zlo = cylinder->lo;
    zhi = cylinder->hi;
    if (axis != 'z') error->all(FLERR, "Must use a z-axis cylinder region with fix pour");
    if (xc - rc < domain->boxlo[0] || xc + rc > domain->boxhi[0] || yc - rc < domain->boxlo[1] ||
        yc + rc > domain->boxhi[1] || zlo < domain->boxlo[2] || zhi > domain->boxhi[2])
      error->all(FLERR, "Insertion region extends outside simulation box");
  } else
    error->all(FLERR, "Must use a block or cylinder region with fix pour");

  if (region_style == 2 && domain->dimension == 2)
    error->all(FLERR, "Must use a block region with fix pour for 2d simulations");

  // error check and further setup for mode = MOLECULE

  if (atom->tag_enable == 0) error->all(FLERR, "Cannot use fix_pour unless atoms have IDs");

  if (mode == MOLECULE) {
    for (int i = 0; i < nmol; i++) {
      if (onemols[i]->xflag == 0) error->all(FLERR, "Fix pour molecule must have coordinates");
      if (onemols[i]->typeflag == 0) error->all(FLERR, "Fix pour molecule must have atom types");
      if (ntype + onemols[i]->ntypes <= 0 || ntype + onemols[i]->ntypes > atom->ntypes)
        error->all(FLERR, "Invalid atom type in fix pour mol command");

      if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols)
        error->all(FLERR, "Fix pour molecule template ID must be same as atom style template ID");
      onemols[i]->check_attributes(0);

      // fix pour uses geoemetric center of molecule for insertion

      onemols[i]->compute_center();
    }
  }

  if (rigidflag && mode == ATOM) error->all(FLERR, "Cannot use fix pour rigid and not molecule");
  if (shakeflag && mode == ATOM) error->all(FLERR, "Cannot use fix pour shake and not molecule");
  if (rigidflag && shakeflag) error->all(FLERR, "Cannot use fix pour rigid and shake");

  // setup of coords and imageflags array

  if (mode == ATOM)
    natom_max = 1;
  else {
    natom_max = 0;
    for (int i = 0; i < nmol; i++) natom_max = MAX(natom_max, onemols[i]->natoms);
  }
  memory->create(coords, natom_max, 4, "pour:coords");
  memory->create(imageflags, natom_max, "pour:imageflags");

  // find max atom and molecule IDs just once

  if (idnext) find_maxid();

  // random number generator, same for all procs
  // warm up the generator 30x to avoid correlations in first-particle
  // positions if runs are repeated with consecutive seeds

  random = new RanPark(lmp, seed);
  for (int ii = 0; ii < 30; ii++) random->uniform();

  // allgather arrays

  MPI_Comm_rank(world, &me);
  MPI_Comm_size(world, &nprocs);
  recvcounts = new int[nprocs];
  displs = new int[nprocs];

  // grav = gravity in distance/time^2 units
  // assume grav = -magnitude at this point, enforce in init()

  auto fixlist = modify->get_fix_by_style("^gravity");
  if (fixlist.size() != 1)
    error->all(FLERR, "There must be exactly one fix gravity defined for fix pour");
  auto fixgrav = dynamic_cast<FixGravity *>(fixlist.front());

  grav = -fixgrav->magnitude * force->ftm2v;

  // nfreq = timesteps between insertions
  // should be time for a particle to fall from top of insertion region
  //   to bottom, taking into account that the region may be moving
  // set these 2 eqs equal to each other, solve for smallest positive t
  //   x = zhi + vz*t + 1/2 grav t^2
  //   x = zlo + rate*t
  //   gives t = [-(vz-rate) - sqrt((vz-rate)^2 - 2*grav*(zhi-zlo))] / grav
  //   where zhi-zlo > 0, grav < 0, and vz & rate can be either > 0 or < 0

  double v_relative, delta;
  if (domain->dimension == 3) {
    v_relative = vz - rate;
    delta = zhi - zlo;
  } else {
    v_relative = vy - rate;
    delta = yhi - ylo;
  }
  double t = (-v_relative - sqrt(v_relative * v_relative - 2.0 * grav * delta)) / grav;
  nfreq = static_cast<int>(t / update->dt + 0.5);

  // 1st insertion on next timestep

  force_reneighbor = 1;
  next_reneighbor = update->ntimestep + 1;
  nfirst = next_reneighbor;
  ninserted = 0;

  // nper = # to insert each time
  // depends on specified volume fraction
  // volume = volume of insertion region
  // volume_one = volume of inserted particle (with max possible radius)
  // in 3d, insure dy >= 1, for quasi-2d simulations

  double volume, volume_one = 1.0;

  molradius_max = 0.0;
  if (mode == MOLECULE) {
    for (int i = 0; i < nmol; i++) molradius_max = MAX(molradius_max, onemols[i]->molradius);
  }

  if (domain->dimension == 3) {
    if (region_style == 1) {
      double dy = yhi - ylo;
      if (dy < 1.0) dy = 1.0;
      volume = (xhi - xlo) * dy * (zhi - zlo);
    } else
      volume = MY_PI * rc * rc * (zhi - zlo);
    if (mode == MOLECULE) {
      volume_one = MY_4PI3 * molradius_max * molradius_max * molradius_max;
    } else if (dstyle == ONE || dstyle == RANGE) {
      volume_one = MY_4PI3 * radius_max * radius_max * radius_max;
    } else if (dstyle == POLY) {
      volume_one = 0.0;
      for (int i = 0; i < npoly; i++)
        volume_one += (MY_4PI3 * radius_poly[i] * radius_poly[i] * radius_poly[i]) * frac_poly[i];
    }
  } else {
    volume = (xhi - xlo) * (yhi - ylo);
    if (mode == MOLECULE) {
      volume_one = MY_PI * molradius_max * molradius_max;
    } else if (dstyle == ONE || dstyle == RANGE) {
      volume_one = MY_PI * radius_max * radius_max;
    } else if (dstyle == POLY) {
      volume_one = 0.0;
      for (int i = 0; i < npoly; i++)
        volume_one += (MY_PI * radius_poly[i] * radius_poly[i]) * frac_poly[i];
    }
  }

  nper = static_cast<int>(volfrac * volume / volume_one);
  if (nper == 0) error->all(FLERR, "Fix pour insertion count per timestep is 0");
  int nfinal = update->ntimestep + 1 + ((bigint) ninsert - 1) / nper * nfreq;

  // print stats

  if (me == 0)
    utils::logmesg(lmp, "Particle insertion: {} every {} steps, {} by step {}\n", nper, nfreq,
                   ninsert, nfinal);
}

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

FixPour::~FixPour()
{
  delete random;
  delete[] molfrac;
  delete[] idrigid;
  delete[] idshake;
  delete[] radius_poly;
  delete[] frac_poly;
  memory->destroy(coords);
  memory->destroy(imageflags);
  delete[] recvcounts;
  delete[] displs;
}

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

int FixPour::setmask()
{
  int mask = 0;
  mask |= PRE_EXCHANGE;
  return mask;
}

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

void FixPour::init()
{
  if (domain->triclinic) error->all(FLERR, "Cannot use fix pour with triclinic box");

  region = domain->get_region_by_id(idregion);
  if (!region) error->all(FLERR, "Fix pour region {} does not exist", idregion);

  // insure gravity fix (still) exists
  // for 3d must point in -z, for 2d must point in -y
  // else insertion cannot work

  auto fixlist = modify->get_fix_by_style("^gravity");
  if (fixlist.size() != 1)
    error->all(FLERR, "There must be exactly one fix gravity defined for fix pour");
  auto fixgrav = dynamic_cast<FixGravity *>(fixlist.front());
  if (fixgrav->varflag != FixGravity::CONSTANT)
    error->all(FLERR, "Fix gravity for fix pour must be constant");

  double xgrav = fixgrav->xgrav;
  double ygrav = fixgrav->ygrav;
  double zgrav = fixgrav->zgrav;

  if (domain->dimension == 3) {
    if (fabs(xgrav) > EPSILON || fabs(ygrav) > EPSILON || fabs(zgrav + 1.0) > EPSILON)
      error->all(FLERR, "Gravity must point in -z to use with fix pour in 3d");
  } else {
    if (fabs(xgrav) > EPSILON || fabs(ygrav + 1.0) > EPSILON || fabs(zgrav) > EPSILON)
      error->all(FLERR, "Gravity must point in -y to use with fix pour in 2d");
  }

  double gnew = -fixgrav->magnitude * force->ftm2v;
  if (gnew != grav) error->all(FLERR, "Gravity changed since fix pour was created");

  // if rigidflag defined, check for rigid/small fix
  // its molecule template must be same as this one

  if (rigidflag) {
    fixrigid = modify->get_fix_by_id(idrigid);
    if (!fixrigid) error->all(FLERR, "Fix pour rigid fix {} does not exist", idrigid);
    int tmp;
    if (onemols != (Molecule **) fixrigid->extract("onemol", tmp))
      error->all(FLERR, "Fix pour and fix rigid/small not using same molecule template ID");
  }

  // if shakeflag defined, check for SHAKE fix
  // its molecule template must be same as this one

  if (shakeflag) {
    fixshake = modify->get_fix_by_id(idshake);
    if (!fixshake) error->all(FLERR, "Fix pour shake fix {} does not exist", idshake);
    int tmp;
    if (onemols != (Molecule **) fixshake->extract("onemol", tmp))
      error->all(FLERR, "Fix pour and fix shake not using same molecule template ID");
  }
}

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

void FixPour::setup_pre_exchange()
{
  if (ninserted < ninsert)
    next_reneighbor = update->ntimestep + 1;
  else
    next_reneighbor = 0;
}

/* ----------------------------------------------------------------------
   perform particle insertion
------------------------------------------------------------------------- */

void FixPour::pre_exchange()
{
  int i, m, flag, nlocalprev, imol, natom;
  double r[3], rotmat[3][3], quat[4], vnew[3];
  double *newcoord;

  // just return if should not be called on this timestep

  if (next_reneighbor != update->ntimestep) return;

  // clear ghost count (and atom map) and any ghost bonus data
  //   internal to AtomVec
  // same logic as beginning of Comm::exchange()
  // do it now b/c inserting atoms will overwrite ghost atoms

  if (atom->map_style != Atom::MAP_NONE) atom->map_clear();
  atom->nghost = 0;
  atom->avec->clear_bonus();

  // find current max atom and molecule IDs on every insertion step

  if (!idnext) find_maxid();

  // nnew = # of particles (atoms or molecules) to insert this timestep

  int nnew = nper;
  if (ninserted + nnew > ninsert) nnew = ninsert - ninserted;

  // lo/hi current = z (or y) bounds of insertion region this timestep

  int dimension = domain->dimension;
  if (dimension == 3) {
    lo_current = zlo + (update->ntimestep - nfirst) * update->dt * rate;
    hi_current = zhi + (update->ntimestep - nfirst) * update->dt * rate;
  } else {
    lo_current = ylo + (update->ntimestep - nfirst) * update->dt * rate;
    hi_current = yhi + (update->ntimestep - nfirst) * update->dt * rate;
  }

  // ncount = # of my atoms that overlap the insertion region
  // nprevious = total of ncount across all procs

  int ncount = 0;
  for (i = 0; i < atom->nlocal; i++)
    if (overlap(i)) ncount++;

  int nprevious;
  MPI_Allreduce(&ncount, &nprevious, 1, MPI_INT, MPI_SUM, world);

  // xmine is for my atoms
  // xnear is for atoms from all procs + atoms to be inserted

  double **xmine, **xnear;
  memory->create(xmine, ncount, 4, "fix_pour:xmine");
  memory->create(xnear, nprevious + nnew * natom_max, 4, "fix_pour:xnear");
  int nnear = nprevious;

  // setup for allgatherv

  int n = 4 * ncount;
  MPI_Allgather(&n, 1, MPI_INT, recvcounts, 1, MPI_INT, world);

  displs[0] = 0;
  for (int iproc = 1; iproc < nprocs; iproc++)
    displs[iproc] = displs[iproc - 1] + recvcounts[iproc - 1];

  // load up xmine array

  double **x = atom->x;
  double *radius = atom->radius;

  ncount = 0;
  for (i = 0; i < atom->nlocal; i++)
    if (overlap(i)) {
      xmine[ncount][0] = x[i][0];
      xmine[ncount][1] = x[i][1];
      xmine[ncount][2] = x[i][2];
      xmine[ncount][3] = radius[i];
      ncount++;
    }

  // perform allgatherv to acquire list of nearby particles on all procs

  double *ptr = nullptr;
  if (ncount) ptr = xmine[0];
  MPI_Allgatherv(ptr, 4 * ncount, MPI_DOUBLE, xnear[0], recvcounts, displs, MPI_DOUBLE, world);

  // insert new particles into xnear list, one by one
  // check against all nearby atoms and previously inserted ones
  // if there is an overlap then try again at same z (3d) or y (2d) coord
  // else insert by adding to xnear list
  // max = maximum # of insertion attempts for all particles
  // h = height, biased to give uniform distribution in time of insertion
  // for MOLECULE mode:
  //   coords = coords of all atoms in particle
  //   perform random rotation around center pt
  //   apply PBC so final coords are inside box
  //   store image flag modified due to PBC

  int success;
  double radtmp, delx, dely, delz, rsq, radsum, rn, h;
  double coord[3];

  double denstmp;
  double *sublo = domain->sublo;
  double *subhi = domain->subhi;

  int nsuccess = 0;
  int attempt = 0;
  int maxiter = nnew * maxattempt;

  while (nsuccess < nnew) {
    rn = random->uniform();
    h = hi_current - rn * rn * (hi_current - lo_current);
    if (mode == ATOM) radtmp = radius_sample();

    success = 0;
    while (attempt < maxiter) {
      attempt++;
      xyz_random(h, coord);

      if (mode == ATOM) {
        natom = 1;
        coords[0][0] = coord[0];
        coords[0][1] = coord[1];
        coords[0][2] = coord[2];
        coords[0][3] = radtmp;
        imageflags[0] = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX;
      } else {
        double rng = random->uniform();
        imol = 0;
        while (rng > molfrac[imol]) imol++;
        natom = onemols[imol]->natoms;
        if (dimension == 3) {
          r[0] = random->uniform() - 0.5;
          r[1] = random->uniform() - 0.5;
          r[2] = random->uniform() - 0.5;
        } else {
          r[0] = r[1] = 0.0;
          r[2] = 1.0;
        }
        double theta = random->uniform() * MY_2PI;
        MathExtra::norm3(r);
        MathExtra::axisangle_to_quat(r, theta, quat);
        MathExtra::quat_to_mat(quat, rotmat);
        for (i = 0; i < natom; i++) {
          MathExtra::matvec(rotmat, onemols[imol]->dx[i], coords[i]);
          coords[i][0] += coord[0];
          coords[i][1] += coord[1];
          coords[i][2] += coord[2];

          // coords[3] = particle radius
          // default to 0.5, if radii not defined in Molecule
          //   same as atom->avec->create_atom(), invoked below

          if (onemols[imol]->radiusflag)
            coords[i][3] = onemols[imol]->radius[i];
          else
            coords[i][3] = 0.5;

          imageflags[i] = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX;
          domain->remap(coords[i], imageflags[i]);
        }
      }

      // if any pair of atoms overlap, try again
      // use minimum_image() to account for PBC

      for (m = 0; m < natom; m++) {
        for (i = 0; i < nnear; i++) {
          delx = coords[m][0] - xnear[i][0];
          dely = coords[m][1] - xnear[i][1];
          delz = coords[m][2] - xnear[i][2];
          domain->minimum_image(delx, dely, delz);
          rsq = delx * delx + dely * dely + delz * delz;
          radsum = coords[m][3] + xnear[i][3];
          if (rsq <= radsum * radsum) break;
        }
        if (i < nnear) break;
      }
      if (m == natom) {
        success = 1;
        break;
      }
    }

    if (!success) break;

    // proceed with insertion

    nsuccess++;
    nlocalprev = atom->nlocal;

    // add all atoms in particle to xnear

    for (m = 0; m < natom; m++) {
      xnear[nnear][0] = coords[m][0];
      xnear[nnear][1] = coords[m][1];
      xnear[nnear][2] = coords[m][2];
      xnear[nnear][3] = coords[m][3];
      nnear++;
    }

    // choose random velocity for new particle
    // used for every atom in molecule
    // z velocity set to what velocity would be if particle
    //   had fallen from top of insertion region
    //   this gives continuous stream of atoms
    //   solution for v from these 2 eqs, after eliminate t:
    //     v = vz + grav*t
    //     coord[2] = hi_current + vz*t + 1/2 grav t^2

    if (dimension == 3) {
      vnew[0] = vxlo + random->uniform() * (vxhi - vxlo);
      vnew[1] = vylo + random->uniform() * (vyhi - vylo);
      vnew[2] = -sqrt(vz * vz + 2.0 * grav * (coord[2] - hi_current));
    } else {
      vnew[0] = vxlo + random->uniform() * (vxhi - vxlo);
      vnew[1] = -sqrt(vy * vy + 2.0 * grav * (coord[1] - hi_current));
      vnew[2] = 0.0;
    }

    // check if new atoms are in my sub-box or above it if I am highest proc
    // if so, add atom to my list via create_atom()
    // initialize additional info about the atoms
    // set group mask to "all" plus fix group

    for (m = 0; m < natom; m++) {
      if (mode == ATOM) denstmp = density_lo + random->uniform() * (density_hi - density_lo);
      newcoord = coords[m];

      flag = 0;
      if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && newcoord[1] >= sublo[1] &&
          newcoord[1] < subhi[1] && newcoord[2] >= sublo[2] && newcoord[2] < subhi[2])
        flag = 1;
      else if (dimension == 3 && newcoord[2] >= domain->boxhi[2]) {
        if (comm->layout != Comm::LAYOUT_TILED) {
          if (comm->myloc[2] == comm->procgrid[2] - 1 && newcoord[0] >= sublo[0] &&
              newcoord[0] < subhi[0] && newcoord[1] >= sublo[1] && newcoord[1] < subhi[1])
            flag = 1;
        } else {
          if (comm->mysplit[2][1] == 1.0 && newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
              newcoord[1] >= sublo[1] && newcoord[1] < subhi[1])
            flag = 1;
        }
      } else if (dimension == 2 && newcoord[1] >= domain->boxhi[1]) {
        if (comm->layout != Comm::LAYOUT_TILED) {
          if (comm->myloc[1] == comm->procgrid[1] - 1 && newcoord[0] >= sublo[0] &&
              newcoord[0] < subhi[0])
            flag = 1;
        } else {
          if (comm->mysplit[1][1] == 1.0 && newcoord[0] >= sublo[0] && newcoord[0] < subhi[0])
            flag = 1;
        }
      }

      if (flag) {
        if (mode == ATOM)
          atom->avec->create_atom(ntype, coords[m]);
        else
          atom->avec->create_atom(ntype + onemols[imol]->type[m], coords[m]);
        int n = atom->nlocal - 1;
        atom->tag[n] = maxtag_all + m + 1;
        if (mode == MOLECULE) {
          if (atom->molecule_flag) {
            if (onemols[imol]->moleculeflag) {
              atom->molecule[n] = maxmol_all + onemols[imol]->molecule[m];
            } else {
              atom->molecule[n] = maxmol_all + 1;
            }
          }
          if (atom->molecular == Atom::TEMPLATE) {
            atom->molindex[n] = 0;
            atom->molatom[n] = m;
          }
        }
        atom->mask[n] = 1 | groupbit;
        atom->image[n] = imageflags[m];
        atom->v[n][0] = vnew[0];
        atom->v[n][1] = vnew[1];
        atom->v[n][2] = vnew[2];
        if (mode == ATOM) {
          radtmp = newcoord[3];
          atom->radius[n] = radtmp;
          atom->rmass[n] = 4.0 * MY_PI / 3.0 * radtmp * radtmp * radtmp * denstmp;
        } else {
          onemols[imol]->quat_external = quat;
          atom->add_molecule_atom(onemols[imol], m, n, maxtag_all);
        }

        modify->create_attribute(n);
      }
    }

    // FixRigidSmall::set_molecule stores rigid body attributes
    //   coord is new position of geometric center of mol, not COM
    // FixShake::set_molecule stores shake info for molecule

    if (rigidflag)
      fixrigid->set_molecule(nlocalprev, maxtag_all, imol, coord, vnew, quat);
    else if (shakeflag)
      fixshake->set_molecule(nlocalprev, maxtag_all, imol, coord, vnew, quat);

    maxtag_all += natom;
    if (mode == MOLECULE && atom->molecule_flag) {
      if (onemols[imol]->moleculeflag) {
        maxmol_all += onemols[imol]->nmolecules;
      } else {
        maxmol_all++;
      }
    }
  }

  // warn if not successful with all insertions b/c too many attempts

  int ninserted_atoms = nnear - nprevious;
  int ninserted_mols = ninserted_atoms / natom;
  ninserted += ninserted_mols;
  if (ninserted_mols < nnew && me == 0) error->warning(FLERR, "Less insertions than requested");

  // reset global natoms,nbonds,etc
  // increment maxtag_all and maxmol_all if necessary
  // if global map exists, reset it now instead of waiting for comm
  //   since other pre-exchange fixes may use it
  //   invoke map_init() b/c atom count has grown

  if (ninserted_atoms) {
    atom->natoms += ninserted_atoms;
    if (atom->natoms < 0) error->all(FLERR, "Too many total atoms");
    if (mode == MOLECULE) {
      atom->nbonds += (bigint) onemols[imol]->nbonds * ninserted_mols;
      atom->nangles += (bigint) onemols[imol]->nangles * ninserted_mols;
      atom->ndihedrals += (bigint) onemols[imol]->ndihedrals * ninserted_mols;
      atom->nimpropers += (bigint) onemols[imol]->nimpropers * ninserted_mols;
    }
    if (maxtag_all >= MAXTAGINT) error->all(FLERR, "New atom IDs exceed maximum allowed ID");
  }

  // rebuild atom map

  if (atom->map_style != Atom::MAP_NONE) {
    if (success) atom->map_init();
    atom->map_set();
  }

  // free local memory

  memory->destroy(xmine);
  memory->destroy(xnear);

  // next timestep to insert

  if (ninserted < ninsert)
    next_reneighbor += nfreq;
  else
    next_reneighbor = 0;
}

/* ----------------------------------------------------------------------
   maxtag_all = current max atom ID for all atoms
   maxmol_all = current max molecule ID for all atoms
------------------------------------------------------------------------- */

void FixPour::find_maxid()
{
  tagint *tag = atom->tag;
  tagint *molecule = atom->molecule;
  int nlocal = atom->nlocal;

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

  if (mode == MOLECULE && molecule) {
    max = 0;
    for (int i = 0; i < nlocal; i++) max = MAX(max, molecule[i]);
    MPI_Allreduce(&max, &maxmol_all, 1, MPI_LMP_TAGINT, MPI_MAX, world);
  }
}

/* ----------------------------------------------------------------------
   check if particle i could overlap with a particle inserted into region
   return 1 if yes, 0 if no
   for ATOM mode, use delta with maximum size for inserted atoms
   for MOLECULE mode, use delta with max radius of inserted molecules
   if ignore line/tri set, ignore line or tri particles
   account for PBC in overlap decision via outside() and minimum_image()
------------------------------------------------------------------------- */

int FixPour::overlap(int i)
{
  double delta;

  if (ignoreflag) {
    if (ignoreline && atom->line[i] >= 0) return 0;
    if (ignoretri && atom->tri[i] >= 0) return 0;
  }

  if (mode == ATOM)
    delta = atom->radius[i] + radius_max;
  else
    delta = atom->radius[i] + molradius_max;

  double *x = atom->x[i];

  if (domain->dimension == 3) {
    if (region_style == 1) {
      if (outside(0, x[0], xlo - delta, xhi + delta)) return 0;
      if (outside(1, x[1], ylo - delta, yhi + delta)) return 0;
      if (outside(2, x[2], lo_current - delta, hi_current + delta)) return 0;
    } else {
      double delx = x[0] - xc;
      double dely = x[1] - yc;
      double delz = 0.0;
      domain->minimum_image(delx, dely, delz);
      double rsq = delx * delx + dely * dely;
      double r = rc + delta;
      if (rsq > r * r) return 0;
      if (outside(2, x[2], lo_current - delta, hi_current + delta)) return 0;
    }
  } else {
    if (outside(0, x[0], xlo - delta, xhi + delta)) return 0;
    if (outside(1, x[1], lo_current - delta, hi_current + delta)) return 0;
  }

  return 1;
}

/* ----------------------------------------------------------------------
   check if value is inside/outside lo/hi bounds in dimension
   account for PBC if needed
   return 1 if value is outside, 0 if inside
------------------------------------------------------------------------- */

bool FixPour::outside(int dim, double value, double lo, double hi)
{
  double boxlo = domain->boxlo[dim];
  double boxhi = domain->boxhi[dim];

  // check for value inside/outside range, ignoring periodicity
  // if inside or dim is non-periodic, only this test is needed

  bool outside_range = (value < lo || value > hi);
  if (!outside_range || !domain->periodicity[dim]) return outside_range;

  // for periodic dimension:
  // must perform additional tests if range wraps around the periodic box

  bool outside_pbc_range = true;

  if ((lo < boxlo && hi > boxhi) || (hi - lo) > domain->prd[dim]) {
    // value is always inside
    outside_pbc_range = false;
  } else if (lo < boxlo) {
    // lower boundary crosses periodic boundary
    outside_pbc_range = (value > hi && value < lo + domain->prd[dim]);
  } else if (hi > boxhi) {
    // upper boundary crosses periodic boundary
    outside_pbc_range = (value < lo && value > hi - domain->prd[dim]);
  }

  return outside_pbc_range;
}

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

void FixPour::xyz_random(double h, double *coord)
{
  if (domain->dimension == 3) {
    if (region_style == 1) {
      coord[0] = xlo + random->uniform() * (xhi - xlo);
      coord[1] = ylo + random->uniform() * (yhi - ylo);
      coord[2] = h;
    } else {
      double r1, r2;
      while (true) {
        r1 = random->uniform() - 0.5;
        r2 = random->uniform() - 0.5;
        if (r1 * r1 + r2 * r2 < 0.25) break;
      }
      coord[0] = xc + 2.0 * r1 * rc;
      coord[1] = yc + 2.0 * r2 * rc;
      coord[2] = h;
    }
  } else {
    coord[0] = xlo + random->uniform() * (xhi - xlo);
    coord[1] = h;
    coord[2] = 0.0;
  }
}

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

double FixPour::radius_sample()
{
  if (dstyle == ONE) return radius_one;
  if (dstyle == RANGE) return radius_lo + random->uniform() * (radius_hi - radius_lo);

  double value = random->uniform();

  int i = 0;
  double sum = 0.0;
  while (sum < value) {
    sum += frac_poly[i];
    i++;
  }
  return radius_poly[i - 1];
}

/* ----------------------------------------------------------------------
   parse optional parameters at end of input line
------------------------------------------------------------------------- */

void FixPour::options(int narg, char **arg)
{
  // defaults

  mode = ATOM;
  rigidflag = 0;
  shakeflag = 0;
  idnext = 0;
  ignoreflag = ignoreline = ignoretri = 0;
  dstyle = ONE;
  radius_max = radius_one = 0.5;
  density_lo = density_hi = 1.0;
  volfrac = 0.25;
  maxattempt = 50;
  rate = 0.0;
  vxlo = vxhi = vylo = vyhi = vy = vz = 0.0;

  int iarg = 0;
  while (iarg < narg) {
    if (strcmp(arg[iarg], "region") == 0) {
      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command");
      region = domain->get_region_by_id(arg[iarg + 1]);
      if (!region) error->all(FLERR, "Fix pour region {} does not exist", arg[iarg + 1]);
      idregion = utils::strdup(arg[iarg + 1]);
      iarg += 2;
    } else if (strcmp(arg[iarg], "mol") == 0) {
      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command");
      int imol = atom->find_molecule(arg[iarg + 1]);
      if (imol == -1) error->all(FLERR, "Molecule template ID for fix pour does not exist");
      mode = MOLECULE;
      onemols = &atom->molecules[imol];
      nmol = onemols[0]->nset;
      delete[] molfrac;
      molfrac = new double[nmol];
      molfrac[0] = 1.0 / nmol;
      for (int i = 1; i < nmol - 1; i++) molfrac[i] = molfrac[i - 1] + 1.0 / nmol;
      molfrac[nmol - 1] = 1.0;
      iarg += 2;
    } else if (strcmp(arg[iarg], "molfrac") == 0) {
      if (mode != MOLECULE) error->all(FLERR, "Illegal fix pour command");
      if (iarg + nmol + 1 > narg) error->all(FLERR, "Illegal fix pour command");
      molfrac[0] = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
      for (int i = 1; i < nmol; i++)
        molfrac[i] = molfrac[i - 1] + utils::numeric(FLERR, arg[iarg + i + 1], false, lmp);
      if (molfrac[nmol - 1] < 1.0 - EPSILON || molfrac[nmol - 1] > 1.0 + EPSILON)
        error->all(FLERR, "Illegal fix pour command");
      molfrac[nmol - 1] = 1.0;
      iarg += nmol + 1;

    } else if (strcmp(arg[iarg], "rigid") == 0) {
      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command");
      delete[] idrigid;
      idrigid = utils::strdup(arg[iarg + 1]);
      rigidflag = 1;
      iarg += 2;
    } else if (strcmp(arg[iarg], "shake") == 0) {
      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command");
      delete[] idshake;
      idshake = utils::strdup(arg[iarg + 1]);
      shakeflag = 1;
      iarg += 2;

    } else if (strcmp(arg[iarg], "id") == 0) {
      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command");
      if (strcmp(arg[iarg + 1], "max") == 0)
        idnext = 0;
      else if (strcmp(arg[iarg + 1], "next") == 0)
        idnext = 1;
      else
        error->all(FLERR, "Illegal fix pour command");
      iarg += 2;

    } else if (strcmp(arg[iarg], "ignore") == 0) {
      if (atom->line_flag) ignoreline = 1;
      if (atom->tri_flag) ignoretri = 1;
      if (ignoreline || ignoretri) ignoreflag = 1;
      iarg += 1;

    } else if (strcmp(arg[iarg], "diam") == 0) {
      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command");
      if (strcmp(arg[iarg + 1], "one") == 0) {
        if (iarg + 3 > narg) error->all(FLERR, "Illegal fix pour command");
        dstyle = ONE;
        radius_one = 0.5 * utils::numeric(FLERR, arg[iarg + 2], false, lmp);
        radius_max = radius_one;
        iarg += 3;
      } else if (strcmp(arg[iarg + 1], "range") == 0) {
        if (iarg + 4 > narg) error->all(FLERR, "Illegal fix pour command");
        dstyle = RANGE;
        radius_lo = 0.5 * utils::numeric(FLERR, arg[iarg + 2], false, lmp);
        radius_hi = 0.5 * utils::numeric(FLERR, arg[iarg + 3], false, lmp);
        if (radius_lo > radius_hi) error->all(FLERR, "Illegal fix pour command");
        radius_max = radius_hi;
        iarg += 4;
      } else if (strcmp(arg[iarg + 1], "poly") == 0) {
        if (iarg + 3 > narg) error->all(FLERR, "Illegal fix pour command");
        dstyle = POLY;
        npoly = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
        if (npoly <= 0) error->all(FLERR, "Illegal fix pour command");
        if (iarg + 3 + 2 * npoly > narg) error->all(FLERR, "Illegal fix pour command");
        radius_poly = new double[npoly];
        frac_poly = new double[npoly];
        iarg += 3;
        radius_max = 0.0;
        for (int i = 0; i < npoly; i++) {
          radius_poly[i] = 0.5 * utils::numeric(FLERR, arg[iarg++], false, lmp);
          frac_poly[i] = utils::numeric(FLERR, arg[iarg++], false, lmp);
          if (radius_poly[i] <= 0.0 || frac_poly[i] < 0.0)
            error->all(FLERR, "Illegal fix pour command");
          radius_max = MAX(radius_max, radius_poly[i]);
        }
        double sum = 0.0;
        for (int i = 0; i < npoly; i++) sum += frac_poly[i];
        if (fabs(sum - 1.0) > SMALL)
          error->all(FLERR, "Fix pour polydisperse fractions do not sum to 1.0");
      } else
        error->all(FLERR, "Illegal fix pour command");

    } else if (strcmp(arg[iarg], "dens") == 0) {
      if (iarg + 3 > narg) error->all(FLERR, "Illegal fix pour command");
      density_lo = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
      density_hi = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
      if (density_lo > density_hi) error->all(FLERR, "Illegal fix pour command");
      iarg += 3;
    } else if (strcmp(arg[iarg], "vol") == 0) {
      if (iarg + 3 > narg) error->all(FLERR, "Illegal fix pour command");
      volfrac = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
      maxattempt = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
      iarg += 3;
    } else if (strcmp(arg[iarg], "rate") == 0) {
      if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command");
      rate = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
      iarg += 2;
    } else if (strcmp(arg[iarg], "vel") == 0) {
      if (domain->dimension == 3) {
        if (iarg + 6 > narg) error->all(FLERR, "Illegal fix pour command");
        vxlo = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
        vxhi = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
        vylo = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
        vyhi = utils::numeric(FLERR, arg[iarg + 4], false, lmp);
        if (vxlo > vxhi || vylo > vyhi) error->all(FLERR, "Illegal fix pour command");
        vz = utils::numeric(FLERR, arg[iarg + 5], false, lmp);
        iarg += 6;
      } else {
        if (iarg + 4 > narg) error->all(FLERR, "Illegal fix pour command");
        vxlo = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
        vxhi = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
        vy = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
        vz = 0.0;
        if (vxlo > vxhi) error->all(FLERR, "Illegal fix pour command");
        iarg += 4;
      }
    } else
      error->all(FLERR, "Illegal fix pour command");
  }
}

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

void FixPour::reset_dt()
{
  error->all(FLERR, "Cannot change timestep with fix pour");
}

/* ----------------------------------------------------------------------
   extract particle radius for atom type = itype
------------------------------------------------------------------------- */

void *FixPour::extract(const char *str, int &itype)
{
  if (strcmp(str, "radius") == 0) {
    if (mode == ATOM) {
      if (itype == ntype)
        oneradius = radius_max;
      else
        oneradius = 0.0;

    } else {

      // loop over onemols molecules
      // skip a molecule with no atoms as large as itype

      oneradius = 0.0;
      for (int i = 0; i < nmol; i++) {
        if (itype > ntype + onemols[i]->ntypes) continue;
        double *radius = onemols[i]->radius;
        int *type = onemols[i]->type;
        int natoms = onemols[i]->natoms;

        // check radii of atoms in Molecule with matching types
        // default to 0.5, if radii not defined in Molecule
        //   same as atom->avec->create_atom(), invoked in pre_exchange()

        for (int i = 0; i < natoms; i++)
          if (type[i] + ntype == itype) {
            if (radius)
              oneradius = MAX(oneradius, radius[i]);
            else
              oneradius = MAX(oneradius, 0.5);
          }
      }
    }
    itype = 0;
    return &oneradius;
  }
  return nullptr;
}