jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11142 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2014-01-07 00:13:57 +00:00
parent bbc2a764c3
commit e1ebfdda1c
4 changed files with 593 additions and 229 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

693
src/GRANULAR/fix_pour.cpp
View File

@ -20,13 +20,16 @@
#include "force.h"
#include "update.h"
#include "comm.h"
#include "molecule.h"
#include "modify.h"
#include "fix_gravity.h"
#include "fix_rigid_small.h"
#include "domain.h"
#include "region.h"
#include "region_block.h"
#include "region_cylinder.h"
#include "random_park.h"
#include "math_extra.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
@ -37,6 +40,7 @@ using namespace MathConst;
#define EPSILON 0.001
enum{ATOM,MOLECULE};
enum{ONE,RANGE,POLY};
/* ---------------------------------------------------------------------- */
@ -59,102 +63,9 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
if (seed <= 0) error->all(FLERR,"Illegal fix pour command");
// option defaults
// read options from end of input line
int iregion = -1;
dstyle = ONE;
radius_max = radius_one = 0.5;
radius_poly = frac_poly = NULL;
density_lo = density_hi = 1.0;
volfrac = 0.25;
maxattempt = 50;
rate = 0.0;
vxlo = vxhi = vylo = vyhi = vy = vz = 0.0;
// optional args
int iarg = 6;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1) error->all(FLERR,"Fix pour region ID does not exist");
iarg += 2;
} else if (strcmp(arg[iarg],"diam") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
if (strcmp(arg[iarg+1],"one") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = ONE;
radius_one = 0.5 * force->numeric(FLERR,arg[iarg+2]);
radius_max = radius_one;
iarg += 3;
} else if (strcmp(arg[iarg+1],"range") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = RANGE;
radius_lo = 0.5 * force->numeric(FLERR,arg[iarg+2]);
radius_hi = 0.5 * force->numeric(FLERR,arg[iarg+3]);
radius_max = radius_hi;
iarg += 4;
} else if (strcmp(arg[iarg+1],"poly") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = POLY;
npoly = force->inumeric(FLERR,arg[iarg+2]);
if (npoly <= 0) error->all(FLERR,"Illegal fix pour command");
if (iarg+3 + 2*npoly > narg)
error->all(FLERR,"Illegal fix pour command");
radius_poly = new double[npoly];
frac_poly = new double[npoly];
iarg += 3;
radius_max = 0.0;
for (int i = 0; i < npoly; i++) {
radius_poly[i] = 0.5 * force->numeric(FLERR,arg[iarg++]);
frac_poly[i] = force->numeric(FLERR,arg[iarg++]);
if (radius_poly[i] <= 0.0 || frac_poly[i] < 0.0)
error->all(FLERR,"Illegal fix pour command");
radius_max = MAX(radius_max,radius_poly[i]);
}
double sum = 0.0;
for (int i = 0; i < npoly; i++) sum += frac_poly[i];
if (sum != 1.0) {
printf("SUM %g\n",sum);
error->all(FLERR,"Fix pour polydisperse fractions do not sum to 1.0");
}
} else error->all(FLERR,"Illegal fix pour command");
} else if (strcmp(arg[iarg],"dens") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
density_lo = force->numeric(FLERR,arg[iarg+1]);
density_hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"vol") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
volfrac = force->numeric(FLERR,arg[iarg+1]);
maxattempt = force->inumeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"rate") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
rate = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"vel") == 0) {
if (domain->dimension == 3) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix pour command");
vxlo = force->numeric(FLERR,arg[iarg+1]);
vxhi = force->numeric(FLERR,arg[iarg+2]);
vylo = force->numeric(FLERR,arg[iarg+3]);
vyhi = force->numeric(FLERR,arg[iarg+4]);
vz = force->numeric(FLERR,arg[iarg+5]);
iarg += 6;
} else {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix pour command");
vxlo = force->numeric(FLERR,arg[iarg+1]);
vxhi = force->numeric(FLERR,arg[iarg+2]);
vy = force->numeric(FLERR,arg[iarg+3]);
vz = 0.0;
iarg += 4;
}
} else error->all(FLERR,"Illegal fix pour command");
}
options(narg-6,&arg[6]);
// error checks on region and its extent being inside simulation box
@ -196,6 +107,35 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
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) {
if (atom->molecule_flag == 0)
error->all(FLERR,"Fix pour requires atom attribute molecule");
if (onemol->xflag == 0)
error->all(FLERR,"Fix pour molecule must have coordinates");
if (onemol->typeflag == 0)
error->all(FLERR,"Fix pour molecule must have atom types");
// fix pour uses geoemetric center of molecule for insertion
onemol->compute_center();
}
// setup of coords and imageflags array
if (mode == ATOM) natom = 1;
else natom = onemol->natoms;
memory->create(coords,natom,4,"pour:coords");
memory->create(imageflags,natom,"pour:imageflags");
// 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);
@ -260,9 +200,12 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
if (dy < 1.0) dy = 1.0;
volume = (xhi-xlo) * dy * (zhi-zlo);
} else volume = MY_PI*rc*rc * (zhi-zlo);
if (dstyle == ONE || dstyle == RANGE)
if (mode == MOLECULE) {
double molradius = onemol->molradius;
volume_one = 4.0/3.0 * MY_PI * molradius*molradius*molradius;
} else if (dstyle == ONE || dstyle == RANGE) {
volume_one = 4.0/3.0 * MY_PI * radius_max*radius_max*radius_max;
else if (dstyle == POLY) {
} else if (dstyle == POLY) {
volume_one = 0.0;
for (int i = 0; i < npoly; i++)
volume_one += (4.0/3.0 * MY_PI *
@ -270,9 +213,12 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
}
} else {
volume = (xhi-xlo) * (yhi-ylo);
if (dstyle == ONE || dstyle == RANGE)
if (mode == MOLECULE) {
double molradius = onemol->molradius;
volume_one = MY_PI * molradius*molradius;
} else if (dstyle == ONE || dstyle == RANGE) {
volume_one = MY_PI * radius_max*radius_max;
else if (dstyle == POLY) {
} 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];
@ -301,8 +247,11 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
FixPour::~FixPour()
{
delete random;
delete [] idrigid;
delete [] radius_poly;
delete [] frac_poly;
memory->destroy(coords);
memory->destroy(imageflags);
delete [] recvcounts;
delete [] displs;
}
@ -352,6 +301,21 @@ void FixPour::init()
double gnew = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v;
if (gnew != grav)
error->all(FLERR,"Gravity changed since fix pour was created");
// if rigidflag defined, check for rigid/small fix
// its molecule template must be same as mine
fixrigid = NULL;
if (rigidflag) {
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"rigid/small") == 0)
fixrigid = (FixRigidSmall *) modify->fix[i];
if (!fixrigid)
error->all(FLERR,"Fix pour rigid requires fix rigid/small definition");
if (onemol != fixrigid->onemol)
error->all(FLERR,
"Fix pour and fix rigid/small not using same molecule ID");
}
}
/* ----------------------------------------------------------------------
@ -360,20 +324,27 @@ void FixPour::init()
void FixPour::pre_exchange()
{
int i;
int i,j,m,flag,nlocalprev;
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;
// nnew = # to insert this timestep
// find current max atom and molecule IDs if necessary
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
if (domain->dimension == 3) {
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 {
@ -396,7 +367,7 @@ void FixPour::pre_exchange()
double **xmine,**xnear;
memory->create(xmine,ncount,4,"fix_pour:xmine");
memory->create(xnear,nprevious+nnew,4,"fix_pour:xnear");
memory->create(xnear,nprevious+nnew*natom,4,"fix_pour:xnear");
int nnear = nprevious;
// setup for allgatherv
@ -430,135 +401,221 @@ void FixPour::pre_exchange()
MPI_Allgatherv(ptr,4*ncount,MPI_DOUBLE,
xnear[0],recvcounts,displs,MPI_DOUBLE,world);
// insert new atoms into xnear list, one by one
// 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 coord[3],radtmp,delx,dely,delz,rsq,radsum,rn,h;
int attempt = 0;
int max = nnew * maxattempt;
int ntotal = nprevious+nnew;
while (nnear < ntotal) {
rn = random->uniform();
h = hi_current - rn*rn * (hi_current-lo_current);
radtmp = radius_sample();
success = 0;
while (attempt < max) {
attempt++;
xyz_random(h,coord);
for (i = 0; i < nnear; i++) {
delx = coord[0] - xnear[i][0];
dely = coord[1] - xnear[i][1];
delz = coord[2] - xnear[i][2];
rsq = delx*delx + dely*dely + delz*delz;
radsum = radtmp + xnear[i][3];
if (rsq <= radsum*radsum) break;
}
if (i == nnear) {
success = 1;
break;
}
}
if (success) {
xnear[nnear][0] = coord[0];
xnear[nnear][1] = coord[1];
xnear[nnear][2] = coord[2];
xnear[nnear][3] = radtmp;
nnear++;
} else break;
}
// warn if not all insertions were performed
ninserted += nnear-nprevious;
if (nnear - nprevious < nnew && me == 0)
error->warning(FLERR,"Less insertions than requested",0);
// check if new atom is in my sub-box or above it if I'm highest proc
// if so, add to my list via create_atom()
// initialize info about the atom
// type, diameter, density set from fix parameters
// group mask set to "all" plus fix group
// 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
// set npartner for new atom to 0 (assume not touching any others)
AtomVec *avec = atom->avec;
int j,m,flag;
double denstmp,vxtmp,vytmp,vztmp;
double *sublo = domain->sublo;
double *subhi = domain->subhi;
double radtmp,delx,dely,delz,rsq,radsum,rn,h;
double coord[3],xcm[3];
int nfix = modify->nfix;
Fix **fix = modify->fix;
for (i = nprevious; i < nnear; i++) {
coord[0] = xnear[i][0];
coord[1] = xnear[i][1];
coord[2] = xnear[i][2];
radtmp = xnear[i][3];
denstmp = density_lo + random->uniform() * (density_hi-density_lo);
if (domain->dimension == 3) {
vxtmp = vxlo + random->uniform() * (vxhi-vxlo);
vytmp = vylo + random->uniform() * (vyhi-vylo);
vztmp = -sqrt(vz*vz + 2.0*grav*(coord[2]-hi_current));
AtomVec *avec = atom->avec;
double denstmp,vxtmp,vytmp,vztmp;
double *sublo = domain->sublo;
double *subhi = domain->subhi;
int attempt = 0;
int maxiter = nnew * maxattempt;
int ntotal = nprevious + nnew*natom;
while (nnear < ntotal) {
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) {
coords[0][0] = coord[0];
coords[0][1] = coord[1];
coords[0][2] = coord[2];
coords[0][3] = radtmp;
imageflags[0] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
} else {
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,onemol->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 not defined in Molecule
// same as atom->avec->create_atom() when invoked below
if (onemol->radiusflag) coords[i][3] = onemol->radius[i];
else coords[i][3] = 0.5;
imageflags[i] = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
domain->remap(coords[i],imageflags[i]);
}
}
// if any pair of atoms overlap, try again
// use minimum_image() to account for PBC
for (m = 0; m < natom; m++) {
for (i = 0; i < nnear; i++) {
delx = coords[m][0] - xnear[i][0];
dely = coords[m][1] - xnear[i][1];
delz = coords[m][2] - xnear[i][2];
domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
radsum = coords[m][3] + xnear[i][3];
if (rsq <= radsum*radsum) break;
}
if (i < nnear) break;
}
if (m == natom) {
success = 1;
break;
}
}
if (!success) break;
// proceed with insertion
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 {
vxtmp = vxlo + random->uniform() * (vxhi-vxlo);
vytmp = -sqrt(vy*vy + 2.0*grav*(coord[1]-hi_current));
vztmp = 0.0;
vnew[0] = vxlo + random->uniform() * (vxhi-vxlo);
vnew[1] = -sqrt(vy*vy + 2.0*grav*(coord[1]-hi_current));
vnew[2] = 0.0;
}
flag = 0;
if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) flag = 1;
else if (domain->dimension == 3 && coord[2] >= domain->boxhi[2] &&
comm->myloc[2] == comm->procgrid[2]-1 &&
coord[0] >= sublo[0] && coord[0] < subhi[0] &&
coord[1] >= sublo[1] && coord[1] < subhi[1]) flag = 1;
else if (domain->dimension == 2 && coord[1] >= domain->boxhi[1] &&
comm->myloc[1] == comm->procgrid[1]-1 &&
coord[0] >= sublo[0] && coord[0] < subhi[0]) flag = 1;
// 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
if (flag) {
avec->create_atom(ntype,coord);
m = atom->nlocal - 1;
atom->type[m] = ntype;
atom->radius[m] = radtmp;
atom->rmass[m] = 4.0*MY_PI/3.0 * radtmp*radtmp*radtmp * denstmp;
atom->mask[m] = 1 | groupbit;
atom->v[m][0] = vxtmp;
atom->v[m][1] = vytmp;
atom->v[m][2] = vztmp;
for (j = 0; j < nfix; j++)
if (fix[j]->create_attribute) fix[j]->set_arrays(m);
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] &&
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 (dimension == 2 && newcoord[1] >= domain->boxhi[1] &&
comm->myloc[1] == comm->procgrid[1]-1 &&
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(onemol->type[m],coords[m]);
int n = atom->nlocal - 1;
atom->tag[n] = maxtag_all + m+1;
if (mode == MOLECULE) atom->molecule[n] = maxmol_all;
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 atom->add_molecule_atom(onemol,m,n,maxtag_all);
for (j = 0; j < nfix; j++)
if (fix[j]->create_attribute) fix[j]->set_arrays(n);
}
}
// FixRigidSmall::set_molecule stores rigid body attributes
// coord is new position of geometric center of mol, not COM
if (fixrigid)
fixrigid->set_molecule(nlocalprev,maxtag_all,coord,vnew,quat);
maxtag_all += natom;
if (mode == MOLECULE) maxmol_all++;
}
// reset global natoms
// set tag # of new particles beyond all previous atoms
// if global map exists, reset it now instead of waiting for comm
// since deleting atoms messes up ghosts
// warn if not successful with all insertions b/c too many attempts
if (nnear - nprevious > 0) {
atom->natoms += nnear - nprevious;
if (atom->tag_enable) {
atom->tag_extend();
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
int ninserted_atoms = nnear - nprevious;
int ninserted_mols = ninserted_atoms / natom;
ninserted += ninserted_mols;
if (ninserted_mols < nnew && me == 0)
error->warning(FLERR,"Less insertions than requested",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 adding atoms messes up ghosts
if (ninserted_atoms) {
atom->natoms += ninserted_atoms;
if (mode == MOLECULE) {
atom->nbonds += onemol->nbonds * ninserted_mols;
atom->nangles += onemol->nangles * ninserted_mols;
atom->ndihedrals += onemol->ndihedrals * ninserted_mols;
atom->nimpropers += onemol->nimpropers * ninserted_mols;
}
//if (idnext) {
// maxtag_all += ninserted_atoms;
// if (mode == MOLECULE) maxmol_all += ninserted_mols;
// }
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
}
@ -573,38 +630,99 @@ void FixPour::pre_exchange()
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()
{
int *tag = atom->tag;
int *molecule = atom->molecule;
int nlocal = atom->nlocal;
int max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
MPI_Allreduce(&max,&maxtag_all,1,MPI_INT,MPI_MAX,world);
if (mode == MOLECULE) {
max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]);
MPI_Allreduce(&max,&maxmol_all,1,MPI_INT,MPI_MAX,world);
}
}
/* ----------------------------------------------------------------------
check if particle i could overlap with a particle inserted into region
return 1 if yes, 0 if no
use radius_hi = maximum size for inserted particle
for ATOM mode, use delta with maximum size for inserted atoms
for MOLECULE mode, use delta with radius of inserted molecules
account for PBC in overlap decision via outside() and minimum_image()
------------------------------------------------------------------------- */
int FixPour::overlap(int i)
{
double delta = radius_max + atom->radius[i];
double **x = atom->x;
double delta;
if (mode == ATOM) delta = atom->radius[i] + radius_max;
else delta = atom->radius[i] + onemol->molradius;
double *boxlo = domain->boxlo;
double *boxhi = domain->boxhi;
double *prd = domain->prd;
int *periodicity = domain->periodicity;
double *x = atom->x[i];
if (domain->dimension == 3) {
if (region_style == 1) {
if (x[i][0] < xlo-delta || x[i][0] > xhi+delta ||
x[i][1] < ylo-delta || x[i][1] > yhi+delta ||
x[i][2] < lo_current-delta || x[i][2] > hi_current+delta) return 0;
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 {
if (x[i][2] < lo_current-delta || x[i][2] > hi_current+delta) return 0;
double delx = x[i][0] - xc;
double dely = x[i][1] - yc;
double delx = x[0] - xc;
double dely = x[1] - yc;
double delz = 0.0;
domain->minimum_image(delx,dely,delz);
double rsq = delx*delx + dely*dely;
double r = rc + delta;
if (rsq > r*r) return 0;
if (outside(2,x[2],lo_current-delta,hi_current+delta)) return 0;
}
} else {
if (x[i][0] < xlo-delta || x[i][0] > xhi+delta ||
x[i][1] < lo_current-delta || x[i][1] > hi_current+delta) return 0;
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
------------------------------------------------------------------------- */
int FixPour::outside(int dim, double value, double lo, double hi)
{
double boxlo = domain->boxlo[dim];
double boxhi = domain->boxhi[dim];
if (domain->periodicity[dim]) {
if (lo < boxlo && hi > boxhi) {
return 0;
} else if (lo < boxlo) {
if (value > hi && value < lo + domain->prd[dim]) return 1;
} else if (hi > boxhi) {
if (value > hi - domain->prd[dim] && value < lo) return 1;
} else {
if (value < lo || value > hi) return 1;
}
}
if (value < lo || value > hi) return 1;
return 0;
}
/* ---------------------------------------------------------------------- */
void FixPour::xyz_random(double h, double *coord)
@ -651,9 +769,136 @@ double FixPour::radius_sample()
return radius_poly[i-1];
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
------------------------------------------------------------------------- */
void FixPour::options(int narg, char **arg)
{
// defaults
iregion = -1;
mode = ATOM;
rigidflag = 0;
idrigid = NULL;
idnext = 0;
dstyle = ONE;
radius_max = radius_one = 0.5;
radius_poly = frac_poly = NULL;
density_lo = density_hi = 1.0;
volfrac = 0.25;
maxattempt = 50;
rate = 0.0;
vxlo = vxhi = vylo = vyhi = vy = vz = 0.0;
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1) error->all(FLERR,"Fix pour region ID does not exist");
iarg += 2;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
int imol = atom->find_molecule(arg[iarg+1]);
if (imol == -1)
error->all(FLERR,"Molecule ID for fix pour does not exist");
mode = MOLECULE;
onemol = atom->molecules[imol];
iarg += 2;
} else if (strcmp(arg[iarg],"rigid") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour 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],"id") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
if (strcmp(arg[iarg+1],"max") == 0) idnext = 0;
else if (strcmp(arg[iarg+1],"next") == 0) idnext = 1;
else error->all(FLERR,"Illegal fix pour command");
iarg += 2;
} else if (strcmp(arg[iarg],"diam") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
if (strcmp(arg[iarg+1],"one") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = ONE;
radius_one = 0.5 * force->numeric(FLERR,arg[iarg+2]);
radius_max = radius_one;
iarg += 3;
} else if (strcmp(arg[iarg+1],"range") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = RANGE;
radius_lo = 0.5 * force->numeric(FLERR,arg[iarg+2]);
radius_hi = 0.5 * force->numeric(FLERR,arg[iarg+3]);
radius_max = radius_hi;
iarg += 4;
} else if (strcmp(arg[iarg+1],"poly") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
dstyle = POLY;
npoly = force->inumeric(FLERR,arg[iarg+2]);
if (npoly <= 0) error->all(FLERR,"Illegal fix pour command");
if (iarg+3 + 2*npoly > narg)
error->all(FLERR,"Illegal fix pour command");
radius_poly = new double[npoly];
frac_poly = new double[npoly];
iarg += 3;
radius_max = 0.0;
for (int i = 0; i < npoly; i++) {
radius_poly[i] = 0.5 * force->numeric(FLERR,arg[iarg++]);
frac_poly[i] = force->numeric(FLERR,arg[iarg++]);
if (radius_poly[i] <= 0.0 || frac_poly[i] < 0.0)
error->all(FLERR,"Illegal fix pour command");
radius_max = MAX(radius_max,radius_poly[i]);
}
double sum = 0.0;
for (int i = 0; i < npoly; i++) sum += frac_poly[i];
if (sum != 1.0)
error->all(FLERR,"Fix pour polydisperse fractions do not sum to 1.0");
} else error->all(FLERR,"Illegal fix pour command");
} else if (strcmp(arg[iarg],"dens") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
density_lo = force->numeric(FLERR,arg[iarg+1]);
density_hi = force->numeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"vol") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command");
volfrac = force->numeric(FLERR,arg[iarg+1]);
maxattempt = force->inumeric(FLERR,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"rate") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
rate = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"vel") == 0) {
if (domain->dimension == 3) {
if (iarg+6 > narg) error->all(FLERR,"Illegal fix pour command");
vxlo = force->numeric(FLERR,arg[iarg+1]);
vxhi = force->numeric(FLERR,arg[iarg+2]);
vylo = force->numeric(FLERR,arg[iarg+3]);
vyhi = force->numeric(FLERR,arg[iarg+4]);
vz = force->numeric(FLERR,arg[iarg+5]);
iarg += 6;
} else {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix pour command");
vxlo = force->numeric(FLERR,arg[iarg+1]);
vxhi = force->numeric(FLERR,arg[iarg+2]);
vy = force->numeric(FLERR,arg[iarg+3]);
vz = 0.0;
iarg += 4;
}
} else error->all(FLERR,"Illegal fix pour command");
}
}
/* ---------------------------------------------------------------------- */
void FixPour::reset_dt()
{
error->all(FLERR,"Cannot change timestep with fix pour");
}

15
src/GRANULAR/fix_pour.h
View File

@ -44,7 +44,7 @@ class FixPour : public Fix {
private:
int ninsert,ntype,seed;
int dstyle,npoly;
int iregion,mode,idnext,dstyle,npoly,rigidflag;
double radius_one,radius_max;
double radius_lo,radius_hi;
double *radius_poly,*frac_poly;
@ -57,16 +57,27 @@ class FixPour : public Fix {
double xlo,xhi,ylo,yhi,zlo,zhi;
double xc,yc,rc;
double grav;
char *idrigid;
class Molecule *onemol;
int natom;
double **coords;
int *imageflags;
class FixRigidSmall *fixrigid;
int me,nprocs;
int *recvcounts,*displs;
int nfreq,nfirst,ninserted,nper;
double lo_current,hi_current;
class RanPark *random;
int maxtag_all,maxmol_all;
class RanPark *random,*random2;
void find_maxid();
int overlap(int);
int outside(int, double, double, double);
void xyz_random(double, double *);
double radius_sample();
void options(int, char **);
};
}

108
src/RIGID/fix_rigid_small.cpp
View File

@ -21,6 +21,7 @@
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "molecule.h"
#include "domain.h"
#include "update.h"
#include "respa.h"
@ -122,6 +123,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
int seed;
langflag = 0;
infile = NULL;
onemol = NULL;
int iarg = 4;
while (iarg < narg) {
@ -138,7 +140,6 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Fix rigid/small langevin period must be > 0.0");
if (seed <= 0) error->all(FLERR,"Illegal fix rigid/small command");
iarg += 5;
} else if (strcmp(arg[iarg],"infile") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
delete [] infile;
@ -147,10 +148,32 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
strcpy(infile,arg[iarg+1]);
restart_file = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command");
int imol = atom->find_molecule(arg[iarg+1]);
if (imol == -1)
error->all(FLERR,"Molecule ID for fix rigid/small does not exist");
onemol = atom->molecules[imol];
iarg += 2;
} else error->all(FLERR,"Illegal fix rigid/small command");
}
// error check and further setup for Molecule template
if (onemol) {
if (onemol->xflag == 0)
error->all(FLERR,"Fix rigid/small molecule must have coordinates");
if (onemol->typeflag == 0)
error->all(FLERR,"Fix rigid/small molecule must have atom types");
// fix rigid/small uses center, masstotal, COM, inertia of molecule
onemol->compute_center();
onemol->compute_mass();
onemol->compute_com();
onemol->compute_inertia();
}
// create rigid bodies based on molecule ID
// sets bodytag for owned atoms
// body attributes are computed later by setup_bodies()
@ -222,7 +245,6 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
MPI_Allreduce(&one,&nbody,1,MPI_INT,MPI_SUM,world);
bigint atomall;
MPI_Allreduce(&atomone,&atomall,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (nbody == 0) error->all(FLERR,"No rigid bodies defined");
if (me == 0) {
if (screen) {
@ -1381,9 +1403,11 @@ void FixRigidSmall::create_bodies()
comm->ring(m,sizeof(double),buf,3,ring_farthest,NULL);
// find maxextent of rsqfar across all procs
// if defined, include molecule->maxextent
MPI_Allreduce(&rsqfar,&maxextent,1,MPI_DOUBLE,MPI_MAX,world);
maxextent = sqrt(maxextent);
if (onemol) maxextent = MAX(maxextent,onemol->maxextent);
// clean up
@ -1547,6 +1571,10 @@ void FixRigidSmall::setup_bodies_static()
MPI_Allreduce(&flag,&extended,1,MPI_INT,MPI_MAX,world);
}
// extended = 1 if using molecule template with finite-size particles
if (onemol && onemol->radiusflag) extended = 1;
// grow extended arrays and set extended flags for each particle
// orientflag = 4 if any particle stores ellipsoid or tri orientation
// orientflag = 1 if any particle stores line orientation
@ -2386,6 +2414,80 @@ void FixRigidSmall::copy_arrays(int i, int j, int delflag)
bodyown[j] = bodyown[i];
}
/* ----------------------------------------------------------------------
initialize a molecule inserted by another fix, e.g. deposit or pour
nlocalprev = # of atoms on this proc before molecule inserted
tagprev = atom ID previous to new atoms in the molecule
xgeom = geometric center of new molecule
vcm = COM velocity of new molecule
quat = rotation of new molecule (around geometric center)
relative to template in Molecule class
------------------------------------------------------------------------- */
void FixRigidSmall::set_molecule(int nlocalprev, int tagprev,
double *xgeom, double *vcm, double *quat)
{
double ctr2com[3],ctr2com_rotate[3];
double rotmat[3][3];
int nlocal = atom->nlocal;
if (nlocalprev == nlocal) return;
double **x = atom->x;
int *tag = atom->tag;
for (int i = nlocalprev; i < nlocal; i++) {
bodytag[i] = tagprev + onemol->comatom;
if (tag[i]-tagprev == onemol->comatom) bodyown[i] = nlocal_body;
int m = tag[i]-tagprev-1;
displace[i][0] = onemol->dxbody[m][0];
displace[i][1] = onemol->dxbody[m][1];
displace[i][2] = onemol->dxbody[m][2];
eflags[i] = 0;
if (onemol->radiusflag) {
eflags[i] |= SPHERE;
eflags[i] |= OMEGA;
eflags[i] |= TORQUE;
}
if (bodyown[i] >= 0) {
if (nlocal_body == nmax_body) grow_body();
Body *b = &body[nlocal_body];
b->mass = onemol->masstotal;
// new COM = Q (onemol->xcm - onemol->center) + xgeom
// Q = rotation matrix associated with quat
MathExtra::quat_to_mat(quat,rotmat);
MathExtra::sub3(onemol->com,onemol->center,ctr2com);
MathExtra::matvec(rotmat,ctr2com,ctr2com_rotate);
MathExtra::add3(ctr2com_rotate,xgeom,b->xcm);
b->vcm[0] = vcm[0];
b->vcm[1] = vcm[1];
b->vcm[2] = vcm[2];
b->inertia[0] = onemol->inertia[0];
b->inertia[1] = onemol->inertia[1];
b->inertia[2] = onemol->inertia[2];
// final quat is product of insertion quat and original quat
// true even if insertion rotation was not around COM
MathExtra::quatquat(quat,onemol->quat,b->quat);
MathExtra::q_to_exyz(b->quat,b->ex_space,b->ey_space,b->ez_space);
b->angmom[0] = b->angmom[1] = b->angmom[2] = 0.0;
b->omega[0] = b->omega[1] = b->omega[2] = 0.0;
b->image = ((tagint) IMGMAX << IMG2BITS) |
((tagint) IMGMAX << IMGBITS) | IMGMAX;
b->ilocal = i;
nlocal_body++;
}
}
}
/* ----------------------------------------------------------------------
initialize one atom's array values, called when atom is created
------------------------------------------------------------------------- */

6
src/RIGID/fix_rigid_small.h
View File

@ -33,6 +33,10 @@ class FixRigidSmall : public Fix {
static FixRigidSmall *frsptr;
// Molecules being added on-the-fly as rigid bodies
class Molecule *onemol;
FixRigidSmall(class LAMMPS *, int, char **);
virtual ~FixRigidSmall();
virtual int setmask();
@ -68,6 +72,8 @@ class FixRigidSmall : public Fix {
double compute_scalar();
double memory_usage();
void set_molecule(int, int, double *, double *, double *);
protected:
int me,nprocs;
double dtv,dtf,dtq;