diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 2d0e438f5b..2357f0d317 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -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"); } + diff --git a/src/GRANULAR/fix_pour.h b/src/GRANULAR/fix_pour.h index 96fd44bcde..2aabc299c8 100644 --- a/src/GRANULAR/fix_pour.h +++ b/src/GRANULAR/fix_pour.h @@ -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 **); }; } diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index ec1db7f0f9..0cf4867a56 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 0bdb686e9b..5681219a50 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -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;