ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
f7f5a3e3e4399f07338c390148ace5e85a0b7c0d
lammps/src/GRANULAR/fix_pour.cpp

1105 lines
38 KiB
C++
Raw Blame History

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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");
scalar_flag = 1;
extscalar = 0;
time_depend = 1;
initialize_flag = 0;
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 (onemols[i]->natoms <= 0) error->all(FLERR, "Fix pour molecule must have atoms");
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();
// 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");
// 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];
}
/* ---------------------------------------------------------------------- */
FixPour::~FixPour()
{
delete random;
delete[] molfrac;
delete[] idrigid;
delete[] idshake;
delete[] idregion;
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);
// Find gravity fix
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");
// Perform one time initialization operations
// outside of constructor in case fix definition precedes timestep definition/atom creation
if (!initialize_flag) {
initialize_flag = 1;
// find max atom and molecule IDs just once
if (idnext) find_maxid();
// grav = gravity in distance/time^2 units
// assume grav = -magnitude at this point, enforce in init()
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 = std::lround(t / update->dt);
// 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, ensure dy >= max diameter, 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 < (2 * radius_max)) dy = (2 * radius_max);
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);
}
// ensure gravity fix (still) exists
// for 3d must point in -z, for 2d must point in -y
// else insertion cannot work
double gnew = -fixgrav->magnitude * force->ftm2v;
if (gnew != grav) error->all(FLERR, "Gravity changed since fix pour was created");
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");
}
// 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(FLERR, 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, "Fewer 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;
// body particle molecule template must contain only one atom
atom->nbodies += (bigint) onemols[imol]->bodyflag * 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) {
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(FLERR, 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) utils::missing_cmd_args(FLERR, "pour region", error);
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) utils::missing_cmd_args(FLERR, "pour mol", error);
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, must specify mol keyword before molfrac");
if (iarg + nmol + 1 > narg) utils::missing_cmd_args(FLERR, "pour molfrac", error);
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: {}", molfrac[nmol - 1]);
molfrac[nmol - 1] = 1.0;
iarg += nmol + 1;
} else if (strcmp(arg[iarg], "rigid") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "pour rigid", error);
delete[] idrigid;
idrigid = utils::strdup(arg[iarg + 1]);
rigidflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg], "shake") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "pour shake", error);
delete[] idshake;
idshake = utils::strdup(arg[iarg + 1]);
shakeflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg], "id") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "pour id", error);
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 id argument: {}", arg[iarg + 1]);
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) utils::missing_cmd_args(FLERR, "pour diam", error);
if (strcmp(arg[iarg + 1], "one") == 0) {
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "pour diam one", error);
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) utils::missing_cmd_args(FLERR, "pour diam range", error);
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 radii: {} exceeds {}", radius_lo, radius_hi);
radius_max = radius_hi;
iarg += 4;
} else if (strcmp(arg[iarg + 1], "poly") == 0) {
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "pour diam poly", error);
dstyle = POLY;
npoly = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
if (npoly <= 0) error->all(FLERR, "Illegal fix pour poly number: {}", npoly);
if (iarg + 3 + 2 * npoly > narg) utils::missing_cmd_args(FLERR, "pour diam poly", error);
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)
error->all(FLERR, "Illegal fix pour poly radius: {}", radius_poly[i]);
if (frac_poly[i] < 0.0)
error->all(FLERR, "Illegal fix pour poly fraction: {}", frac_poly[i]);
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) utils::missing_cmd_args(FLERR, "pour dens", error);
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 densities: {} exceeds {}", density_lo, density_hi);
iarg += 3;
} else if (strcmp(arg[iarg], "vol") == 0) {
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "pour vol", error);
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) utils::missing_cmd_args(FLERR, "pour rate", error);
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) utils::missing_cmd_args(FLERR, "pour vel", error);
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) error->all(FLERR, "Illegal fix pour x velocities: {} exceeds {}", vxlo, vxhi);
if (vylo > vyhi) error->all(FLERR, "Illegal fix pour y velocities: {} exceeds {}", vylo, vyhi);
vz = utils::numeric(FLERR, arg[iarg + 5], false, lmp);
iarg += 6;
} else {
if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "pour vel", error);
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 velocities: {} exceeds {}", vxlo, vxhi);
iarg += 4;
}
} else
error->all(FLERR, "Illegal fix pour command argument: {}", arg[iarg]);
}
}
/* ----------------------------------------------------------------------
output number of successful insertions
------------------------------------------------------------------------- */
double FixPour::compute_scalar()
{
return ninserted;
}
/* ---------------------------------------------------------------------- */
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;
}
Reference in New Issue View Git Blame Copy Permalink