complete moving packages and styles around
This commit is contained in:
Axel Kohlmeyer
907
src/fix_deposit.cpp
Normal file
907
src/fix_deposit.cpp
Normal file
@ -0,0 +1,907 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
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_deposit.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "atom_vec.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fix.h"
|
||||
#include "lattice.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 "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
using namespace MathConst;
|
||||
|
||||
enum{ATOM,MOLECULE};
|
||||
enum{DIST_UNIFORM,DIST_GAUSSIAN};
|
||||
|
||||
#define EPSILON 1.0e6
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg), idregion(nullptr), idrigid(nullptr),
|
||||
idshake(nullptr), onemols(nullptr), molfrac(nullptr), coords(nullptr), imageflags(nullptr),
|
||||
fixrigid(nullptr), fixshake(nullptr), random(nullptr)
|
||||
{
|
||||
if (narg < 7) error->all(FLERR,"Illegal fix deposit command");
|
||||
|
||||
restart_global = 1;
|
||||
time_depend = 1;
|
||||
|
||||
// required args
|
||||
|
||||
ninsert = utils::inumeric(FLERR,arg[3],false,lmp);
|
||||
ntype = utils::inumeric(FLERR,arg[4],false,lmp);
|
||||
nfreq = utils::inumeric(FLERR,arg[5],false,lmp);
|
||||
seed = utils::inumeric(FLERR,arg[6],false,lmp);
|
||||
|
||||
if (seed <= 0) error->all(FLERR,"Illegal fix deposit command");
|
||||
|
||||
// read options from end of input line
|
||||
|
||||
options(narg-7,&arg[7]);
|
||||
|
||||
// error check on type
|
||||
|
||||
if (mode == ATOM && (ntype <= 0 || ntype > atom->ntypes))
|
||||
error->all(FLERR,"Invalid atom type in fix deposit command");
|
||||
|
||||
// error checks on region and its extent being inside simulation box
|
||||
|
||||
if (iregion == -1) error->all(FLERR,"Must specify a region in fix deposit");
|
||||
if (domain->regions[iregion]->bboxflag == 0)
|
||||
error->all(FLERR,"Fix deposit region does not support a bounding box");
|
||||
if (domain->regions[iregion]->dynamic_check())
|
||||
error->all(FLERR,"Fix deposit region cannot be dynamic");
|
||||
|
||||
xlo = domain->regions[iregion]->extent_xlo;
|
||||
xhi = domain->regions[iregion]->extent_xhi;
|
||||
ylo = domain->regions[iregion]->extent_ylo;
|
||||
yhi = domain->regions[iregion]->extent_yhi;
|
||||
zlo = domain->regions[iregion]->extent_zlo;
|
||||
zhi = domain->regions[iregion]->extent_zhi;
|
||||
|
||||
if (domain->triclinic == 0) {
|
||||
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,"Deposition region extends outside simulation box");
|
||||
} else {
|
||||
if (xlo < domain->boxlo_bound[0] || xhi > domain->boxhi_bound[0] ||
|
||||
ylo < domain->boxlo_bound[1] || yhi > domain->boxhi_bound[1] ||
|
||||
zlo < domain->boxlo_bound[2] || zhi > domain->boxhi_bound[2])
|
||||
error->all(FLERR,"Deposition region extends outside simulation box");
|
||||
}
|
||||
|
||||
// error check and further setup for mode = MOLECULE
|
||||
|
||||
if (atom->tag_enable == 0)
|
||||
error->all(FLERR,"Cannot use fix_deposit unless atoms have IDs");
|
||||
|
||||
if (mode == MOLECULE) {
|
||||
for (int i = 0; i < nmol; i++) {
|
||||
if (onemols[i]->xflag == 0)
|
||||
error->all(FLERR,"Fix deposit molecule must have coordinates");
|
||||
if (onemols[i]->typeflag == 0)
|
||||
error->all(FLERR,"Fix deposit 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 deposit mol command");
|
||||
|
||||
if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols)
|
||||
error->all(FLERR,"Fix deposit molecule template ID must be same "
|
||||
"as atom_style template ID");
|
||||
onemols[i]->check_attributes(0);
|
||||
|
||||
// fix deposit uses geoemetric center of molecule for insertion
|
||||
|
||||
onemols[i]->compute_center();
|
||||
}
|
||||
}
|
||||
|
||||
if (rigidflag && mode == ATOM)
|
||||
error->all(FLERR,"Cannot use fix deposit rigid and not molecule");
|
||||
if (shakeflag && mode == ATOM)
|
||||
error->all(FLERR,"Cannot use fix deposit shake and not molecule");
|
||||
if (rigidflag && shakeflag)
|
||||
error->all(FLERR,"Cannot use fix deposit 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,3,"deposit:coords");
|
||||
memory->create(imageflags,natom_max,"deposit:imageflags");
|
||||
|
||||
// setup scaling
|
||||
|
||||
double xscale,yscale,zscale;
|
||||
if (scaleflag) {
|
||||
xscale = domain->lattice->xlattice;
|
||||
yscale = domain->lattice->ylattice;
|
||||
zscale = domain->lattice->zlattice;
|
||||
}
|
||||
else xscale = yscale = zscale = 1.0;
|
||||
|
||||
// apply scaling to all input parameters with dist/vel units
|
||||
|
||||
if (domain->dimension == 2) {
|
||||
lo *= yscale;
|
||||
hi *= yscale;
|
||||
rate *= yscale;
|
||||
} else {
|
||||
lo *= zscale;
|
||||
hi *= zscale;
|
||||
rate *= zscale;
|
||||
}
|
||||
deltasq *= xscale*xscale;
|
||||
nearsq *= xscale*xscale;
|
||||
vxlo *= xscale;
|
||||
vxhi *= xscale;
|
||||
vylo *= yscale;
|
||||
vyhi *= yscale;
|
||||
vzlo *= zscale;
|
||||
vzhi *= zscale;
|
||||
xmid *= xscale;
|
||||
ymid *= yscale;
|
||||
zmid *= zscale;
|
||||
sigma *= xscale; // same as in region sphere
|
||||
tx *= xscale;
|
||||
ty *= yscale;
|
||||
tz *= zscale;
|
||||
|
||||
// find current max atom and molecule IDs if necessary
|
||||
|
||||
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();
|
||||
|
||||
// set up reneighboring
|
||||
|
||||
force_reneighbor = 1;
|
||||
next_reneighbor = update->ntimestep + 1;
|
||||
nfirst = next_reneighbor;
|
||||
ninserted = 0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixDeposit::~FixDeposit()
|
||||
{
|
||||
delete random;
|
||||
delete [] molfrac;
|
||||
delete [] idrigid;
|
||||
delete [] idshake;
|
||||
delete [] idregion;
|
||||
memory->destroy(coords);
|
||||
memory->destroy(imageflags);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixDeposit::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= PRE_EXCHANGE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixDeposit::init()
|
||||
{
|
||||
// set index and check validity of region
|
||||
|
||||
iregion = domain->find_region(idregion);
|
||||
if (iregion == -1)
|
||||
error->all(FLERR,"Region ID for fix deposit does not exist");
|
||||
|
||||
// if rigidflag defined, check for rigid/small fix
|
||||
// its molecule template must be same as this one
|
||||
|
||||
fixrigid = nullptr;
|
||||
if (rigidflag) {
|
||||
int ifix = modify->find_fix(idrigid);
|
||||
if (ifix < 0) error->all(FLERR,"Fix deposit rigid fix does not exist");
|
||||
fixrigid = modify->fix[ifix];
|
||||
int tmp;
|
||||
if (onemols != (Molecule **) fixrigid->extract("onemol",tmp))
|
||||
error->all(FLERR,
|
||||
"Fix deposit 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
|
||||
|
||||
fixshake = nullptr;
|
||||
if (shakeflag) {
|
||||
int ifix = modify->find_fix(idshake);
|
||||
if (ifix < 0) error->all(FLERR,"Fix deposit shake fix does not exist");
|
||||
fixshake = modify->fix[ifix];
|
||||
int tmp;
|
||||
if (onemols != (Molecule **) fixshake->extract("onemol",tmp))
|
||||
error->all(FLERR,"Fix deposit and fix shake not using "
|
||||
"same molecule template ID");
|
||||
}
|
||||
|
||||
// for finite size spherical particles:
|
||||
// warn if near < 2 * maxrad of existing and inserted particles
|
||||
// since may lead to overlaps
|
||||
// if inserted molecule does not define diameters,
|
||||
// use AtomVecSphere::create_atom() default radius = 0.5
|
||||
|
||||
if (atom->radius_flag) {
|
||||
double *radius = atom->radius;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
double maxrad = 0.0;
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
maxrad = MAX(maxrad,radius[i]);
|
||||
|
||||
double maxradall;
|
||||
MPI_Allreduce(&maxrad,&maxradall,1,MPI_DOUBLE,MPI_MAX,world);
|
||||
|
||||
double maxradinsert = 0.0;
|
||||
if (mode == MOLECULE) {
|
||||
for (int i = 0; i < nmol; i++) {
|
||||
if (onemols[i]->radiusflag)
|
||||
maxradinsert = MAX(maxradinsert,onemols[i]->maxradius);
|
||||
else maxradinsert = MAX(maxradinsert,0.5);
|
||||
}
|
||||
} else maxradinsert = 0.5;
|
||||
|
||||
double separation = MAX(2.0*maxradinsert,maxradall+maxradinsert);
|
||||
if (sqrt(nearsq) < separation && comm->me == 0)
|
||||
error->warning(FLERR,"Fix deposit near setting < possible "
|
||||
"overlap separation {}",separation);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixDeposit::setup_pre_exchange()
|
||||
{
|
||||
if (ninserted < ninsert) next_reneighbor = update->ntimestep+1;
|
||||
else next_reneighbor = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
perform particle insertion
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixDeposit::pre_exchange()
|
||||
{
|
||||
int i,m,n,nlocalprev,imol,natom,flag,flagall;
|
||||
double coord[3],lamda[3],delx,dely,delz,rsq;
|
||||
double r[3],vnew[3],rotmat[3][3],quat[4];
|
||||
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();
|
||||
|
||||
// compute current offset = bottom of insertion volume
|
||||
|
||||
double offset = 0.0;
|
||||
if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate;
|
||||
|
||||
double *sublo,*subhi;
|
||||
if (domain->triclinic == 0) {
|
||||
sublo = domain->sublo;
|
||||
subhi = domain->subhi;
|
||||
} else {
|
||||
sublo = domain->sublo_lamda;
|
||||
subhi = domain->subhi_lamda;
|
||||
}
|
||||
|
||||
// find current max atom and molecule IDs if necessary
|
||||
|
||||
if (!idnext) find_maxid();
|
||||
|
||||
// attempt an insertion until successful
|
||||
|
||||
int dimension = domain->dimension;
|
||||
|
||||
int success = 0;
|
||||
int attempt = 0;
|
||||
while (attempt < maxattempt) {
|
||||
attempt++;
|
||||
|
||||
// choose random position for new particle within region
|
||||
if (distflag == DIST_UNIFORM) {
|
||||
do {
|
||||
coord[0] = xlo + random->uniform() * (xhi-xlo);
|
||||
coord[1] = ylo + random->uniform() * (yhi-ylo);
|
||||
coord[2] = zlo + random->uniform() * (zhi-zlo);
|
||||
} while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0);
|
||||
} else if (distflag == DIST_GAUSSIAN) {
|
||||
do {
|
||||
coord[0] = xmid + random->gaussian() * sigma;
|
||||
coord[1] = ymid + random->gaussian() * sigma;
|
||||
coord[2] = zmid + random->gaussian() * sigma;
|
||||
} while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0);
|
||||
} else error->all(FLERR,"Unknown particle distribution in fix deposit");
|
||||
|
||||
// adjust vertical coord by offset
|
||||
|
||||
if (dimension == 2) coord[1] += offset;
|
||||
else coord[2] += offset;
|
||||
|
||||
// if global, reset vertical coord to be lo-hi above highest atom
|
||||
// if local, reset vertical coord to be lo-hi above highest "nearby" atom
|
||||
// local computation computes lateral distance between 2 particles w/ PBC
|
||||
// when done, have final coord of atom or center pt of molecule
|
||||
|
||||
if (globalflag || localflag) {
|
||||
int dim;
|
||||
double max,maxall,delx,dely,delz,rsq;
|
||||
|
||||
if (dimension == 2) {
|
||||
dim = 1;
|
||||
max = domain->boxlo[1];
|
||||
} else {
|
||||
dim = 2;
|
||||
max = domain->boxlo[2];
|
||||
}
|
||||
|
||||
double **x = atom->x;
|
||||
int nlocal = atom->nlocal;
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
if (localflag) {
|
||||
delx = coord[0] - x[i][0];
|
||||
dely = coord[1] - x[i][1];
|
||||
delz = 0.0;
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
if (dimension == 2) rsq = delx*delx;
|
||||
else rsq = delx*delx + dely*dely;
|
||||
if (rsq > deltasq) continue;
|
||||
}
|
||||
if (x[i][dim] > max) max = x[i][dim];
|
||||
}
|
||||
|
||||
MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world);
|
||||
if (dimension == 2)
|
||||
coord[1] = maxall + lo + random->uniform()*(hi-lo);
|
||||
else
|
||||
coord[2] = maxall + lo + random->uniform()*(hi-lo);
|
||||
}
|
||||
|
||||
// coords = coords of all atoms
|
||||
// for molecule, perform random rotation around center pt
|
||||
// apply PBC so final coords are inside box
|
||||
// also modify image flags due to PBC
|
||||
|
||||
if (mode == ATOM) {
|
||||
natom = 1;
|
||||
coords[0][0] = coord[0];
|
||||
coords[0][1] = coord[1];
|
||||
coords[0][2] = coord[2];
|
||||
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) {
|
||||
if (orientflag) {
|
||||
r[0] = rx;
|
||||
r[1] = ry;
|
||||
r[2] = rz;
|
||||
} else {
|
||||
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];
|
||||
|
||||
imageflags[i] = ((imageint) IMGMAX << IMG2BITS) |
|
||||
((imageint) IMGMAX << IMGBITS) | IMGMAX;
|
||||
domain->remap(coords[i],imageflags[i]);
|
||||
}
|
||||
}
|
||||
|
||||
// check distance between any existing atom and any inserted atom
|
||||
// if less than near, try again
|
||||
// use minimum_image() to account for PBC
|
||||
|
||||
double **x = atom->x;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
flag = 0;
|
||||
for (m = 0; m < natom; m++) {
|
||||
for (i = 0; i < nlocal; i++) {
|
||||
delx = coords[m][0] - x[i][0];
|
||||
dely = coords[m][1] - x[i][1];
|
||||
delz = coords[m][2] - x[i][2];
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
if (rsq < nearsq) flag = 1;
|
||||
}
|
||||
}
|
||||
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
|
||||
if (flagall) continue;
|
||||
|
||||
// proceed with insertion
|
||||
|
||||
nlocalprev = atom->nlocal;
|
||||
|
||||
// choose random velocity for new particle
|
||||
// used for every atom in molecule
|
||||
|
||||
vnew[0] = vxlo + random->uniform() * (vxhi-vxlo);
|
||||
vnew[1] = vylo + random->uniform() * (vyhi-vylo);
|
||||
vnew[2] = vzlo + random->uniform() * (vzhi-vzlo);
|
||||
|
||||
// if target specified, change velocity vector accordingly
|
||||
|
||||
if (targetflag) {
|
||||
double vel = sqrt(vnew[0]*vnew[0] + vnew[1]*vnew[1] + vnew[2]*vnew[2]);
|
||||
delx = tx - coord[0];
|
||||
dely = ty - coord[1];
|
||||
delz = tz - coord[2];
|
||||
double rsq = delx*delx + dely*dely + delz*delz;
|
||||
if (rsq > 0.0) {
|
||||
double rinv = sqrt(1.0/rsq);
|
||||
vnew[0] = delx*rinv*vel;
|
||||
vnew[1] = dely*rinv*vel;
|
||||
vnew[2] = delz*rinv*vel;
|
||||
}
|
||||
}
|
||||
|
||||
// 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 (domain->triclinic) {
|
||||
domain->x2lamda(coords[m],lamda);
|
||||
newcoord = lamda;
|
||||
} else 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]);
|
||||
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 == MOLECULE) {
|
||||
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);
|
||||
|
||||
success = 1;
|
||||
break;
|
||||
}
|
||||
|
||||
// warn if not successful b/c too many attempts
|
||||
|
||||
if (!success && comm->me == 0)
|
||||
error->warning(FLERR,"Particle deposition was unsuccessful");
|
||||
|
||||
// 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 (success) {
|
||||
atom->natoms += natom;
|
||||
if (atom->natoms < 0)
|
||||
error->all(FLERR,"Too many total atoms");
|
||||
if (mode == MOLECULE) {
|
||||
atom->nbonds += onemols[imol]->nbonds;
|
||||
atom->nangles += onemols[imol]->nangles;
|
||||
atom->ndihedrals += onemols[imol]->ndihedrals;
|
||||
atom->nimpropers += onemols[imol]->nimpropers;
|
||||
}
|
||||
maxtag_all += natom;
|
||||
if (maxtag_all >= MAXTAGINT)
|
||||
error->all(FLERR,"New atom IDs exceed maximum allowed ID");
|
||||
if (mode == MOLECULE && atom->molecule_flag) {
|
||||
if (onemols[imol]->moleculeflag) {
|
||||
maxmol_all += onemols[imol]->nmolecules;
|
||||
} else {
|
||||
maxmol_all++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// rebuild atom map
|
||||
|
||||
if (atom->map_style != Atom::MAP_NONE) {
|
||||
if (success) atom->map_init();
|
||||
atom->map_set();
|
||||
}
|
||||
|
||||
// next timestep to insert
|
||||
// next_reneighbor = 0 if done
|
||||
|
||||
if (success) ninserted++;
|
||||
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 FixDeposit::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);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
parse optional parameters at end of input line
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixDeposit::options(int narg, char **arg)
|
||||
{
|
||||
// defaults
|
||||
|
||||
iregion = -1;
|
||||
idregion = nullptr;
|
||||
mode = ATOM;
|
||||
molfrac = nullptr;
|
||||
rigidflag = 0;
|
||||
idrigid = nullptr;
|
||||
shakeflag = 0;
|
||||
idshake = nullptr;
|
||||
idnext = 0;
|
||||
globalflag = localflag = 0;
|
||||
lo = hi = deltasq = 0.0;
|
||||
nearsq = 0.0;
|
||||
maxattempt = 10;
|
||||
rateflag = 0;
|
||||
vxlo = vxhi = vylo = vyhi = vzlo = vzhi = 0.0;
|
||||
distflag = DIST_UNIFORM;
|
||||
sigma = 1.0;
|
||||
xmid = ymid = zmid = 0.0;
|
||||
scaleflag = 1;
|
||||
targetflag = 0;
|
||||
orientflag = 0;
|
||||
rx = 0.0;
|
||||
ry = 0.0;
|
||||
rz = 0.0;
|
||||
|
||||
int iarg = 0;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"region") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
iregion = domain->find_region(arg[iarg+1]);
|
||||
if (iregion == -1)
|
||||
error->all(FLERR,"Region ID for fix deposit does not exist");
|
||||
idregion = utils::strdup(arg[iarg+1]);
|
||||
iarg += 2;
|
||||
|
||||
} else if (strcmp(arg[iarg],"mol") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
int imol = atom->find_molecule(arg[iarg+1]);
|
||||
if (imol == -1)
|
||||
error->all(FLERR,"Molecule template ID for fix deposit 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 deposit command");
|
||||
if (iarg+nmol+1 > narg) error->all(FLERR,"Illegal fix deposit 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 deposit 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 deposit 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 deposit 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 deposit 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 deposit command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"global") == 0) {
|
||||
if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
globalflag = 1;
|
||||
localflag = 0;
|
||||
lo = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
hi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg],"local") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
localflag = 1;
|
||||
globalflag = 0;
|
||||
lo = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
hi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
deltasq = utils::numeric(FLERR,arg[iarg+3],false,lmp) *
|
||||
utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
iarg += 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"near") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
nearsq = utils::numeric(FLERR,arg[iarg+1],false,lmp) *
|
||||
utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"attempt") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
maxattempt = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"rate") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
rateflag = 1;
|
||||
rate = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"vx") == 0) {
|
||||
if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
vxlo = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
vxhi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg],"vy") == 0) {
|
||||
if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
vylo = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
vyhi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg],"vz") == 0) {
|
||||
if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
vzlo = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
vzhi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg],"orient") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
orientflag = 1;
|
||||
rx = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
ry = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
rz = utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
if (domain->dimension == 2 && (rx != 0.0 || ry != 0.0))
|
||||
error->all(FLERR,"Illegal fix deposit orient settings");
|
||||
if (rx == 0.0 && ry == 0.0 && rz == 0.0)
|
||||
error->all(FLERR,"Illegal fix deposit orient settings");
|
||||
iarg += 4;
|
||||
} else if (strcmp(arg[iarg],"units") == 0) {
|
||||
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
|
||||
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
|
||||
else error->all(FLERR,"Illegal fix deposit command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"gaussian") == 0) {
|
||||
if (iarg+5 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
xmid = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
ymid = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
zmid = utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
sigma = utils::numeric(FLERR,arg[iarg+4],false,lmp);
|
||||
distflag = DIST_GAUSSIAN;
|
||||
iarg += 5;
|
||||
} else if (strcmp(arg[iarg],"target") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
|
||||
tx = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
ty = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
tz = utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
targetflag = 1;
|
||||
iarg += 4;
|
||||
} else error->all(FLERR,"Illegal fix deposit command");
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
pack entire state of Fix into one write
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixDeposit::write_restart(FILE *fp)
|
||||
{
|
||||
int n = 0;
|
||||
double list[5];
|
||||
list[n++] = random->state();
|
||||
list[n++] = ninserted;
|
||||
list[n++] = ubuf(nfirst).d;
|
||||
list[n++] = ubuf(next_reneighbor).d;
|
||||
list[n++] = ubuf(update->ntimestep).d;
|
||||
|
||||
if (comm->me == 0) {
|
||||
int size = n * sizeof(double);
|
||||
fwrite(&size,sizeof(int),1,fp);
|
||||
fwrite(list,sizeof(double),n,fp);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
use state info from restart file to restart the Fix
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixDeposit::restart(char *buf)
|
||||
{
|
||||
int n = 0;
|
||||
double *list = (double *) buf;
|
||||
|
||||
seed = static_cast<int>(list[n++]);
|
||||
ninserted = static_cast<int>(list[n++]);
|
||||
nfirst = static_cast<bigint>(ubuf(list[n++]).i);
|
||||
next_reneighbor = static_cast<bigint>(ubuf(list[n++]).i);
|
||||
|
||||
bigint ntimestep_restart = static_cast<bigint>(ubuf(list[n++]).i);
|
||||
if (ntimestep_restart != update->ntimestep)
|
||||
error->all(FLERR,"Must not reset timestep when restarting this fix");
|
||||
|
||||
random->reset(seed);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
extract particle radius for atom type = itype
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void *FixDeposit::extract(const char *str, int &itype)
|
||||
{
|
||||
if (strcmp(str,"radius") == 0) {
|
||||
if (mode == ATOM) {
|
||||
if (itype == ntype) oneradius = 0.5;
|
||||
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
Block a user