fix_deposit.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, Jun 5, 12:05

fix_deposit.cpp
View Options

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

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

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

	#include <math.h>
	#include <stdlib.h>
	#include <string.h>
	#include "fix_deposit.h"
	#include "atom.h"
	#include "atom_vec.h"
	#include "molecule.h"
	#include "force.h"
	#include "update.h"
	#include "modify.h"
	#include "fix.h"
	#include "comm.h"
	#include "domain.h"
	#include "lattice.h"
	#include "region.h"
	#include "random_park.h"
	#include "math_extra.h"
	#include "math_const.h"
	#include "memory.h"
	#include "error.h"

	using namespace LAMMPS_NS;
	using namespace FixConst;
	using namespace MathConst;

	enum{ATOM,MOLECULE};
	enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files
	enum{DIST_UNIFORM,DIST_GAUSSIAN};

	#define EPSILON 1.0e6

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

	FixDeposit::FixDeposit(LAMMPS lmp, int narg, char *arg) :
	Fix(lmp, narg, arg), idregion(NULL), idrigid(NULL),
	idshake(NULL), onemols(NULL), molfrac(NULL), coords(NULL), imageflags(NULL),
	fixrigid(NULL), fixshake(NULL), random(NULL)
	{
	if (narg < 7) error->all(FLERR,"Illegal fix deposit command");

	restart_global = 1;
	time_depend = 1;

	// required args

	ninsert = force->inumeric(FLERR,arg[3]);
	ntype = force->inumeric(FLERR,arg[4]);
	nfreq = force->inumeric(FLERR,arg[5]);
	seed = force->inumeric(FLERR,arg[6]);

	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 == 2 && 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 = xscalexscale;
	nearsq = xscalexscale;
	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

	random = new RanPark(lmp,seed);

	// 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 = NULL;
	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 = NULL;
	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) {
	char str[128];
	sprintf(str,"Fix deposit near setting < possible overlap separation %g",
	separation);
	error->warning(FLERR,str);
	}
	}
	}

	/* ----------------------------------------------------------------------
	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 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

	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 = delxdelx + delydely;
	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) {
	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 = delxdelx + delydely + 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 = delxdelx + delydely + delz*delz;
	if (rsq > 0.0) {
	double rinv = sqrt(1.0/rsq);
	vnew[0] = delxrinvvel;
	vnew[1] = delyrinvvel;
	vnew[2] = delzrinvvel;
	}
	}

	// 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 != 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 != 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) atom->molecule[n] = maxmol_all+1;
	if (atom->molecular == 2) {
	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);

	// old code: unsuccessful if no proc performed insertion of an atom
	// don't think that check is necessary
	// if get this far, should always be succesful
	// would be hard to undo partial insertion for a molecule
	// better to check how many atoms could be inserted (w/out inserting)
	// then sum to insure all are inserted, before doing actual insertion
	// MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world);

	success = 1;
	break;
	}

	// warn if not successful b/c too many attempts

	if (!success && comm->me == 0)
	error->warning(FLERR,"Particle deposition was unsuccessful",0);

	// 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) maxmol_all++;
	if (atom->map_style) {
	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 = NULL;
	mode = ATOM;
	molfrac = NULL;
	rigidflag = 0;
	idrigid = NULL;
	shakeflag = 0;
	idshake = NULL;
	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;

	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");
	int n = strlen(arg[iarg+1]) + 1;
	idregion = new char[n];
	strcpy(idregion,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] = force->numeric(FLERR,arg[iarg+1]);
	for (int i = 1; i < nmol; i++)
	molfrac[i] = molfrac[i-1] + force->numeric(FLERR,arg[iarg+i+1]);
	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");
	int n = strlen(arg[iarg+1]) + 1;
	delete [] idrigid;
	idrigid = new char[n];
	strcpy(idrigid,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");
	int n = strlen(arg[iarg+1]) + 1;
	delete [] idshake;
	idshake = new char[n];
	strcpy(idshake,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 = force->numeric(FLERR,arg[iarg+1]);
	hi = force->numeric(FLERR,arg[iarg+2]);
	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 = force->numeric(FLERR,arg[iarg+1]);
	hi = force->numeric(FLERR,arg[iarg+2]);
	deltasq = force->numeric(FLERR,arg[iarg+3]) *
	force->numeric(FLERR,arg[iarg+3]);
	iarg += 4;

	} else if (strcmp(arg[iarg],"near") == 0) {
	if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
	nearsq = force->numeric(FLERR,arg[iarg+1]) *
	force->numeric(FLERR,arg[iarg+1]);
	iarg += 2;
	} else if (strcmp(arg[iarg],"attempt") == 0) {
	if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
	maxattempt = force->inumeric(FLERR,arg[iarg+1]);
	iarg += 2;
	} else if (strcmp(arg[iarg],"rate") == 0) {
	if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
	rateflag = 1;
	rate = force->numeric(FLERR,arg[iarg+1]);
	iarg += 2;
	} else if (strcmp(arg[iarg],"vx") == 0) {
	if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
	vxlo = force->numeric(FLERR,arg[iarg+1]);
	vxhi = force->numeric(FLERR,arg[iarg+2]);
	iarg += 3;
	} else if (strcmp(arg[iarg],"vy") == 0) {
	if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
	vylo = force->numeric(FLERR,arg[iarg+1]);
	vyhi = force->numeric(FLERR,arg[iarg+2]);
	iarg += 3;
	} else if (strcmp(arg[iarg],"vz") == 0) {
	if (iarg+3 > narg) error->all(FLERR,"Illegal fix deposit command");
	vzlo = force->numeric(FLERR,arg[iarg+1]);
	vzhi = force->numeric(FLERR,arg[iarg+2]);
	iarg += 3;
	} 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 = force->numeric(FLERR,arg[iarg+1]);
	ymid = force->numeric(FLERR,arg[iarg+2]);
	zmid = force->numeric(FLERR,arg[iarg+3]);
	sigma = force->numeric(FLERR,arg[iarg+4]);
	distflag = DIST_GAUSSIAN;
	iarg += 5;
	} else if (strcmp(arg[iarg],"target") == 0) {
	if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
	tx = force->numeric(FLERR,arg[iarg+1]);
	ty = force->numeric(FLERR,arg[iarg+2]);
	tz = force->numeric(FLERR,arg[iarg+3]);
	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[4];
	list[n++] = random->state();
	list[n++] = ninserted;
	list[n++] = nfirst;
	list[n++] = next_reneighbor;

	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<int> (list[n++]);
	next_reneighbor = static_cast<int> (list[n++]);

	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 NULL;
	}

fix_deposit.cppNo OneTemporaryActions

File Metadata

fix_deposit.cppView Options

Event Timeline

fix_deposit.cpp
No OneTemporary
Actions

fix_deposit.cpp
View Options