From 90ff2eb6c97638be07f8bbddbbfa6daec081cc40 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 21 Jan 2020 13:31:47 -0700 Subject: [PATCH] modified versions of creating atoms on subset of lattice, ditto for set type/fraction --- doc/src/create_atoms.rst | 26 +++- doc/src/set.rst | 68 ++++++--- src/atom_vec.cpp | 15 ++ src/atom_vec.h | 1 + src/create_atoms.cpp | 302 +++++++++++++++++---------------------- src/create_atoms.h | 16 ++- src/random_mars.cpp | 132 ++++++++++++++++- src/random_mars.h | 1 + src/read_data.cpp | 6 + src/read_restart.cpp | 12 +- src/replicate.cpp | 6 + src/set.cpp | 136 ++++++++++++++---- src/set.h | 3 +- 13 files changed, 486 insertions(+), 238 deletions(-) diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst index fa4e6b7c59..595095e213 100644 --- a/doc/src/create_atoms.rst +++ b/doc/src/create_atoms.rst @@ -27,7 +27,7 @@ Syntax region-ID = create atoms within this region, use NULL for entire simulation box * zero or more keyword/value pairs may be appended -* keyword = *mol* or *basis* or *subset* or *remap* or *var* or *set* or *units* +* keyword = *mol* or *basis* or *ratio* or *subset* or *remap* or *var* or *set* or *rotate* or *units* .. parsed-literal:: @@ -37,7 +37,12 @@ Syntax *basis* values = M itype M = which basis atom itype = atom type (1-N) to assign to this basis atom - *subset* value = N = number of particles to add at lattice points + *ratio* values = frac seed + frac = fraction of lattice sites (0 to 1) to populate randomly + seed = random # seed (positive integer) + *subset* values = Nsubset seed + Nsubset = # of lattice sites to populate randomly + seed = random # seed (positive integer) *remap* value = *yes* or *no* *var* value = name = variable name to evaluate for test of atom creation *set* values = dim name @@ -60,6 +65,7 @@ Examples create_atoms 1 box create_atoms 3 region regsphere basis 2 3 + create_atoms 3 region regsphere basis 2 3 ratio 0.5 74637 create_atoms 3 single 0 0 5 create_atoms 1 box var v set x xpos set y ypos @@ -215,10 +221,18 @@ command for specifics on how basis atoms are defined for the unit cell of the lattice. By default, all created atoms are assigned the argument *type* as their atom type. -The *subset* keyword can be used in conjunction with the *box* or -*region* styles to limit the total number of particles inserted. The -specified number of particles N are placed randomly on the available -lattice points. +The *ratio* and *subset* keywords can be used in conjunction with the +*box* or *region* styles to limit the total number of particles +inserted. The lattice defines a set of *Nlatt* eligible sites for +inserting particles, which may be limited by the *region* style or the +*var* and *set* keywords. For the *ratio* keyword only the specified +fraction of them (0 <= *frac* <= 1) will be assigned particles. For +the *subset* keyword only the specified *Nsubset* of them will be +assigned particles. In both cases the assigned lattice sites are +chosen randomly. An iterative algorithm is used which insures the +correct number of particles are inserted, in a perfectly random +fashion. Which lattice sites are selected will change with the number +of processors used. The *remap* keyword only applies to the *single* style. If it is set to *yes*\ , then if the specified position is outside the simulation diff --git a/doc/src/set.rst b/doc/src/set.rst index 7a9bb4662c..087c3487e5 100644 --- a/doc/src/set.rst +++ b/doc/src/set.rst @@ -14,7 +14,7 @@ Syntax * style = *atom* or *type* or *mol* or *group* or *region* * ID = atom ID range or type range or mol ID range or group ID or region ID * one or more keyword/value pairs may be appended -* keyword = *type* or *type/fraction* or *mol* or *x* or *y* or *z* or *charge* or *dipole* or *dipole/random* or *quat* or *spin* or *spin/random* or *quat* or *quat/random* or *diameter* or *shape* or *length* or *tri* or *theta* or *theta/random* or *angmom* or *omega* or *mass* or *density* or *density/disc* or *volume* or *image* or *bond* or *angle* or *dihedral* or *improper* or *meso/e* or *meso/cv* or *meso/rho* or *smd/contact/radius* or *smd/mass/density* or *dpd/theta* or *edpd/temp* or *edpd/cv* or *cc* or *i\_name* or *d\_name* +* keyword = *type* or *type/fraction* or *type/ratio* or *type/subset* or *mol* or *x* or *y* or *z* or *charge* or *dipole* or *dipole/random* or *quat* or *spin* or *spin/random* or *quat* or *quat/random* or *diameter* or *shape* or *length* or *tri* or *theta* or *theta/random* or *angmom* or *omega* or *mass* or *density* or *density/disc* or *volume* or *image* or *bond* or *angle* or *dihedral* or *improper* or *meso/e* or *meso/cv* or *meso/rho* or *smd/contact/radius* or *smd/mass/density* or *dpd/theta* or *edpd/temp* or *edpd/cv* or *cc* or *i\_name* or *d\_name* .. parsed-literal:: @@ -22,7 +22,15 @@ Syntax value can be an atom-style variable (see below) *type/fraction* values = type fraction seed type = new atom type - fraction = fraction of selected atoms to set to new atom type + fraction = approximate fraction of selected atoms to set to new atom type + seed = random # seed (positive integer) + *type/ratio* values = type fraction seed + type = new atom type + fraction = exact fraction of selected atoms to set to new atom type + seed = random # seed (positive integer) + *type/subset* values = type Nsubset seed + type = new atom type + Nsubset = exact number of selected atoms to set to new atom type seed = random # seed (positive integer) *mol* value = molecule ID value can be an atom-style variable (see below) @@ -184,15 +192,16 @@ This section describes the keyword options for which properties to change, for the selected atoms. Note that except where explicitly prohibited below, all of the -keywords allow an :doc:`atom-style or atomfile-style variable ` to be used as the specified value(s). If the -value is a variable, it should be specified as v\_name, where name is -the variable name. In this case, the variable will be evaluated, and -its resulting per-atom value used to determine the value assigned to -each selected atom. Note that the per-atom value from the variable -will be ignored for atoms that are not selected via the *style* and -*ID* settings explained above. A simple way to use per-atom values -from the variable to reset a property for all atoms is to use style -*atom* with *ID* = "\*"; this selects all atom IDs. +keywords allow an :doc:`atom-style or atomfile-style variable +` to be used as the specified value(s). If the value is a +variable, it should be specified as v\_name, where name is the +variable name. In this case, the variable will be evaluated, and its +resulting per-atom value used to determine the value assigned to each +selected atom. Note that the per-atom value from the variable will be +ignored for atoms that are not selected via the *style* and *ID* +settings explained above. A simple way to use per-atom values from +the variable to reset a property for all atoms is to use style *atom* +with *ID* = "\*"; this selects all atom IDs. Atom-style variables can specify formulas with various mathematical functions, and include :doc:`thermo\_style ` command @@ -220,23 +229,36 @@ command. Keyword *type/fraction* sets the atom type for a fraction of the selected atoms. The actual number of atoms changed is not guaranteed -to be exactly the requested fraction, but should be statistically -close. Random numbers are used in such a way that a particular atom -is changed or not changed, regardless of how many processors are being -used. This keyword does not allow use of an atom-style variable. +to be exactly the specified fraction (0 <= *fraction* <= 1), but +should be statistically close. Random numbers are used in such a way +that a particular atom is changed or not changed, regardless of how +many processors are being used. This keyword does not allow use of an +atom-style variable. -Keyword *mol* sets the molecule ID for all selected atoms. The :doc:`atom style ` being used must support the use of molecule -IDs. +Keywords *type/ratio* and *type/subset" also set the atom type for a +fraction of the selected atoms. The actual number of atoms changed +will be exactly the requested number. For *type/ratio* the specified +fraction (0 <= *fraction* <= 1) determines the number. For +*type/subset*, the specified *Nsubset* is the number. An iterative +algorithm is used which insures the correct number of atoms are +selected, in a perfectly random fashion. Which atoms are selected +will change with the number of processors used. These keywords do not +allow use of an atom-style variable. -Keywords *x*\ , *y*\ , *z*\ , and *charge* set the coordinates or charge of -all selected atoms. For *charge*\ , the :doc:`atom style ` -being used must support the use of atomic charge. Keywords *vx*\ , *vy*\ , -and *vz* set the velocities of all selected atoms. +Keyword *mol* sets the molecule ID for all selected atoms. The +:doc:`atom style ` being used must support the use of +molecule IDs. + +Keywords *x*\ , *y*\ , *z*\ , and *charge* set the coordinates or +charge of all selected atoms. For *charge*\ , the :doc:`atom style +` being used must support the use of atomic +charge. Keywords *vx*\ , *vy*\ , and *vz* set the velocities of all +selected atoms. Keyword *dipole* uses the specified x,y,z values as components of a vector to set as the orientation of the dipole moment vectors of the -selected atoms. The magnitude of the dipole moment is set -by the length of this orientation vector. +selected atoms. The magnitude of the dipole moment is set by the +length of this orientation vector. Keyword *dipole/random* randomizes the orientation of the dipole moment vectors for the selected atoms and sets the magnitude of each diff --git a/src/atom_vec.cpp b/src/atom_vec.cpp index bc94e36e7a..17fa5b2fb7 100644 --- a/src/atom_vec.cpp +++ b/src/atom_vec.cpp @@ -19,6 +19,8 @@ #include "error.h" #include "utils.h" +#include "comm.h" + using namespace LAMMPS_NS; #define DELTA 16384 @@ -108,6 +110,19 @@ int AtomVec::grow_nmax_bonus(int nmax_bonus) return nmax_bonus; } +/* ---------------------------------------------------------------------- + roundup N so it is a multiple of DELTA + error if N exceeds 32-bit int, since will be used as arg to grow() +------------------------------------------------------------------------- */ + +bigint AtomVec::roundup(bigint n) +{ + if (n % DELTA) n = n/DELTA * DELTA + DELTA; + if (n > MAXSMALLINT) + error->one(FLERR,"Too many atoms created on one or more procs"); + return n; +} + /* ---------------------------------------------------------------------- unpack one line from Velocities section of data file ------------------------------------------------------------------------- */ diff --git a/src/atom_vec.h b/src/atom_vec.h index 2b57238c3b..db4139204a 100644 --- a/src/atom_vec.h +++ b/src/atom_vec.h @@ -57,6 +57,7 @@ class AtomVec : protected Pointers { virtual void grow(int) = 0; virtual void grow_reset() = 0; + bigint roundup(bigint); virtual void copy(int, int, int) = 0; virtual void clear_bonus() {} virtual void force_clear(int, size_t) {} diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 088512cb82..0461ec702a 100644 --- a/src/create_atoms.cpp +++ b/src/create_atoms.cpp @@ -39,9 +39,12 @@ using namespace MathConst; #define BIG 1.0e30 #define EPSILON 1.0e-6 +#define LB_FACTOR 1.1 enum{BOX,REGION,SINGLE,RANDOM}; enum{ATOM,MOLECULE}; +enum{COUNT,INSERT,INSERT_SELECTED}; +enum{NONE,RATIO,SUBSET}; /* ---------------------------------------------------------------------- */ @@ -108,12 +111,11 @@ void CreateAtoms::command(int narg, char **arg) remapflag = 0; mode = ATOM; int molseed; - int subsetseed; varflag = 0; vstr = xstr = ystr = zstr = NULL; quatone[0] = quatone[1] = quatone[2] = 0.0; - nlatt = nsubset = subsetflag = 0; - flag = NULL; + subsetflag = NONE; + int subsetseed; nbasis = domain->lattice->nbasis; basistype = new int[nbasis]; @@ -194,16 +196,21 @@ void CreateAtoms::command(int narg, char **arg) MathExtra::norm3(axisone); MathExtra::axisangle_to_quat(axisone,thetaone,quatone); iarg += 5; + } else if (strcmp(arg[iarg],"ratio") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command"); + subsetflag = RATIO; + subsetfrac = force->numeric(FLERR,arg[iarg+1]); + subsetseed = force->inumeric(FLERR,arg[iarg+2]); + if (subsetfrac <= 0.0 || subsetfrac > 1.0 || subsetseed <= 0) + error->all(FLERR,"Illegal create_atoms command"); + iarg += 3; } else if (strcmp(arg[iarg],"subset") == 0) { if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command"); - nsubset = force->inumeric(FLERR,arg[iarg+1]); + subsetflag = SUBSET; + nsubset = force->bnumeric(FLERR,arg[iarg+1]); subsetseed = force->inumeric(FLERR,arg[iarg+2]); - if (nsubset > 0) subsetflag = 1; - else { - if (me == 0) error->warning(FLERR,"Specifying an 'subset' value of " - "'0' is equivalent to no 'subset' keyword"); - subsetflag = 0; - } + if (nsubset <= 0 || subsetseed <= 0) + error->all(FLERR,"Illegal create_atoms command"); iarg += 3; } else error->all(FLERR,"Illegal create_atoms command"); } @@ -242,9 +249,8 @@ void CreateAtoms::command(int narg, char **arg) ranmol = new RanMars(lmp,molseed+me); } - if (subsetflag) { - ranlatt = new RanMars(lmp,subsetseed+me); - } + ranlatt = NULL; + if (subsetflag != NONE) ranlatt = new RanMars(lmp,subsetseed+me); // error check and further setup for variable test @@ -534,12 +540,13 @@ void CreateAtoms::command(int narg, char **arg) // clean up delete ranmol; + delete ranlatt; + if (domain->lattice) delete [] basistype; delete [] vstr; delete [] xstr; delete [] ystr; delete [] zstr; - memory->destroy(flag); // for MOLECULE mode: // create special bond lists for molecular systems, @@ -780,7 +787,6 @@ void CreateAtoms::add_lattice() // which can lead to missing atoms in rare cases // extra decrement of lo if min < 0, since static_cast(-1.5) = -1 - int ilo,ihi,jlo,jhi,klo,khi; ilo = static_cast (xmin) - 1; jlo = static_cast (ymin) - 1; klo = static_cast (zmin) - 1; @@ -792,22 +798,18 @@ void CreateAtoms::add_lattice() if (ymin < 0.0) jlo--; if (zmin < 0.0) klo--; - // iterate on 3d periodic lattice of unit cells using loop bounds - // iterate on nbasis atoms in each unit cell - // convert lattice coords to box coords - // add atom or molecule (on each basis point) if it meets all criteria - - const double * const * const basis = domain->lattice->basis; - - // rough estimate of total time used for create atoms. + // rough estimate of total time used for create atoms // one inner loop takes about 25ns on a typical desktop CPU core in 2019 - double testimate = 2.5e-8/3600.0; // convert seconds to hours - testimate *= static_cast(khi-klo+1); - testimate *= static_cast(jhi-jlo+1); - testimate *= static_cast(ihi-ilo+1); - testimate *= static_cast(nbasis); + // maxestimate = time in hours + + double estimate = 2.5e-8/3600.0; + estimate *= static_cast (khi-klo+1); + estimate *= static_cast (jhi-jlo+1); + estimate *= static_cast (ihi-ilo+1); + estimate *= static_cast (nbasis); + double maxestimate = 0.0; - MPI_Reduce(&testimate,&maxestimate,1,MPI_DOUBLE,MPI_MAX,0,world); + MPI_Reduce(&estimate,&maxestimate,1,MPI_DOUBLE,MPI_MAX,0,world); if ((comm->me == 0) && (maxestimate > 0.01)) { if (screen) fprintf(screen,"WARNING: create_atoms will take " @@ -816,181 +818,133 @@ void CreateAtoms::add_lattice() "approx. %.2f hours to complete\n",maxestimate); } - int i,j,k,m; + // count lattice sites on each proc - // one pass for default mode, two passes for subset mode: - // first pass: count how many particles will be inserted - // second pass: filter to N number of particles (and insert) - int iflag = 0; - int npass = 1; - if (subsetflag) npass = 2; - for (int pass = 0; pass < npass; pass++) { - if (pass == 1) get_subset(); - for (k = klo; k <= khi; k++) { - for (j = jlo; j <= jhi; j++) { - for (i = ilo; i <= ihi; i++) { - for (m = 0; m < nbasis; m++) { - double *coord; - double x[3],lamda[3]; + nlatt_overflow = 0; + loop_lattice(COUNT); - x[0] = i + basis[m][0]; - x[1] = j + basis[m][1]; - x[2] = k + basis[m][2]; + // nadd = # of atoms each proc will insert (estimated if subsetflag) - // convert from lattice coords to box coords + int overflow; + MPI_Allreduce(&nlatt_overflow,&overflow,1,MPI_INT,MPI_SUM,world); + if (overflow) + error->all(FLERR,"Create_atoms lattice size overflow on 1 or more procs"); - domain->lattice->lattice2box(x[0],x[1],x[2]); + bigint nadd; - // if a region was specified, test if atom is in it + if (subsetflag == NONE) { + if (nprocs == 1) nadd = nlatt; + else nadd = static_cast (LB_FACTOR * nlatt); + } else { + bigint bnlatt = nlatt; + bigint bnlattall; + MPI_Allreduce(&bnlatt,&bnlattall,1,MPI_LMP_BIGINT,MPI_SUM,world); + if (subsetflag == RATIO) + nsubset = static_cast (subsetfrac * bnlattall); + if (nsubset > bnlattall) + error->all(FLERR,"Create_atoms subset size > # of lattice sites"); + if (nprocs == 1) nadd = nsubset; + else nadd = static_cast (LB_FACTOR * nsubset/bnlattall * nlatt); + } + + // allocate atom arrays to size N, rounded up by AtomVec->DELTA - if (style == REGION) - if (!domain->regions[nregion]->match(x[0],x[1],x[2])) continue; + bigint nbig = atom->avec->roundup(nadd + atom->nlocal); + int n = static_cast (nbig); + atom->avec->grow(n); - // if variable test specified, eval variable + // add atoms or molecules + // if no subset: add to all lattice sites + // if subset: count lattice sites, select random subset, then add - if (varflag && vartest(x) == 0) continue; - - // test if atom/molecule position is in my subbox - - if (triclinic) { - domain->x2lamda(x,lamda); - coord = lamda; - } else coord = x; - - 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]) continue; - - // add the atom or entire molecule to my list of atoms - if (subsetflag && pass == 0) nlatt++; - else { - if (!subsetflag || flag[iflag++] == 1) - if (mode == ATOM) atom->avec->create_atom(basistype[m],x); - else if (quatone[0] == 0 && quatone[1] == 0 && quatone[2] == 0) - add_molecule(x); - else add_molecule(x,quatone); - } - } - } - } - } + if (subsetflag == NONE) loop_lattice(INSERT); + else { + memory->create(flag,nlatt,"create_atoms:flag"); + memory->create(next,nlatt,"create_atoms:next"); + ranlatt->select_subset(nsubset,nlatt,flag,next); + loop_lattice(INSERT_SELECTED); + memory->destroy(flag); + memory->destroy(next); } } /* ---------------------------------------------------------------------- - define a random mask to insert only N particles on lattice + iterate on 3d periodic lattice of unit cells using loop bounds + iterate on nbasis atoms in each unit cell + convert lattice coords to box coords + check if lattice point meets all criteria to be added + perform action on atom or molecule (on each basis point) if meets all criteria + actions = add, count, add if flagged ------------------------------------------------------------------------- */ -void CreateAtoms::get_subset() +void CreateAtoms::loop_lattice(int action) { - enum{ATOMS,HOLES}; - int i,j,temp,irand,offlag,npicks,pickmode,mynpicks,mysubset; - double myrand; - int *allnlatts; - int *allsubsets; - int *localpicks; - double *proc_sects; + int i,j,k,m; - memory->create(allnlatts,nprocs,"create_atoms:allnlatts"); - memory->create(allsubsets,nprocs,"create_atoms:allsubsets"); - memory->create(localpicks,nprocs,"create_atoms:localpicks"); - memory->create(proc_sects,nprocs,"create_atoms:proc_sects"); + const double * const * const basis = domain->lattice->basis; - MPI_Allgather(&nlatt, 1, MPI_INT, &allnlatts[0], 1, MPI_INT, world); + nlatt = 0; - int ntotal = 0; - for (i = 0; i < nprocs; i++) - ntotal += allnlatts[i]; + for (k = klo; k <= khi; k++) { + for (j = jlo; j <= jhi; j++) { + for (i = ilo; i <= ihi; i++) { + for (m = 0; m < nbasis; m++) { + double *coord; + double x[3],lamda[3]; + + x[0] = i + basis[m][0]; + x[1] = j + basis[m][1]; + x[2] = k + basis[m][2]; - if (nsubset > ntotal) - error->all(FLERR,"Attempting to insert more particles than available lattice points"); + // convert from lattice coords to box coords + + domain->lattice->lattice2box(x[0],x[1],x[2]); - // define regions of unity based on a proc's fraction of total lattice points - proc_sects[0] = (double) allnlatts[0] / (double) ntotal; - for (i = 1; i < nprocs; i++) - proc_sects[i] = proc_sects[i-1] + (double) allnlatts[i] / (double) ntotal; + // if a region was specified, test if atom is in it - if (nsubset > ntotal/2) { - pickmode = HOLES; - npicks = ntotal - nsubset; - } else { - pickmode = ATOMS; - npicks = nsubset; - } + if (style == REGION) + if (!domain->regions[nregion]->match(x[0],x[1],x[2])) continue; - mynpicks = npicks/nprocs; - if (me == 0) mynpicks = npicks - (nprocs-1)*(mynpicks); + // if variable test specified, eval variable + + if (varflag && vartest(x) == 0) continue; + + // test if atom/molecule position is in my subbox + + if (triclinic) { + domain->x2lamda(x,lamda); + coord = lamda; + } else coord = x; + + 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]) continue; + + // this proc owns the lattice site + // perform action: add, just count, add if flagged + // add = add an atom or entire molecule to my list of atoms - for (i = 0; i < nprocs; i++) - localpicks[i] = 0; + if (action == INSERT) { + if (mode == ATOM) atom->avec->create_atom(basistype[m],x); + else if (quatone[0] == 0 && quatone[1] == 0 && quatone[2] == 0) + add_molecule(x); + else add_molecule(x,quatone); + } else if (action == COUNT) { + if (nlatt == MAXSMALLINT) nlatt_overflow = 1; + } else if (action == INSERT_SELECTED && flag[nlatt]) { + if (mode == ATOM) atom->avec->create_atom(basistype[m],x); + else if (quatone[0] == 0 && quatone[1] == 0 && quatone[2] == 0) + add_molecule(x); + else add_molecule(x,quatone); + } - for (i = 0; i < mynpicks; i++) { - myrand = ranlatt->uniform(); - for (j = 0; j < nprocs; j++) - if (myrand < proc_sects[j]) { - localpicks[j]++; - break; - } - } - - MPI_Allreduce(&localpicks[0],&allsubsets[0],nprocs,MPI_INT,MPI_SUM,world); - - if (pickmode == HOLES) - for (i = 0; i < nprocs; i++) - allsubsets[i] = allnlatts[i] - allsubsets[i]; - - mysubset = allsubsets[me]; - - // it's possible, but statistically unlikely, that a proc was assigned too many - // proc 0 will fix this - offlag = 0; - npicks = 0; - for (i = 0; i < nprocs; i++) { - if (allsubsets[i] > allnlatts[i]) { - offlag = 1; - npicks += allsubsets[i] - allnlatts[i]; - } - } - - if (offlag == 1) { - if (me == 0) { - while (npicks > 0) { // while loop - myrand = ranlatt->uniform(); - for (j = 0; j < nprocs; j++) - if (myrand < proc_sects[j]) break; - if (allsubsets[j] < allnlatts[j]) { - allsubsets[j]++; - npicks--; + nlatt++; } } } - MPI_Scatter(&allsubsets[0], 1, MPI_INT, &mysubset, 1, MPI_INT, 0, world); } - - // each processor chooses its random lattice points - memory->create(flag,nlatt,"create_atoms:flag"); - - for (i = 0; i < nlatt; i++) - if (i < mysubset) - flag[i] = 1; - else - flag[i] = 0; - - // shuffle filled lattice points - for (i = nlatt-1; i > 0; --i) { - irand = floor(ranlatt->uniform()*(i+1)); - temp = flag[i]; - flag[i] = flag[irand]; - flag[irand] = temp; - } - - memory->destroy(allnlatts); - memory->destroy(allsubsets); - memory->destroy(localpicks); - memory->destroy(proc_sects); } - /* ---------------------------------------------------------------------- add a randomly rotated molecule with its center at center if quat_user set, perform requested rotation diff --git a/src/create_atoms.h b/src/create_atoms.h index 3ebb84db42..5375d8dc70 100644 --- a/src/create_atoms.h +++ b/src/create_atoms.h @@ -32,16 +32,24 @@ class CreateAtoms : protected Pointers { private: int me,nprocs; int ntype,style,mode,nregion,nbasis,nrandom,seed; + int remapflag; + int subsetflag; + bigint nsubset; + double subsetfrac; int *basistype; double xone[3],quatone[4]; - int remapflag; int varflag,vvar,xvar,yvar,zvar; char *vstr,*xstr,*ystr,*zstr; char *xstr_copy,*ystr_copy,*zstr_copy; - int nlatt,nsubset,subsetflag; - int *flag; // used to insert N number of particles on lattice + int ilo,ihi,jlo,jhi,klo,khi; + + int nlatt; // number of owned lattice sites + int nlatt_overflow; // 1 if local nlatt exceeds a 32-bit int + + int *flag; // flag subset of particles to insert on lattice + int *next; class Molecule *onemol; class RanMars *ranmol; @@ -53,7 +61,7 @@ class CreateAtoms : protected Pointers { void add_single(); void add_random(); void add_lattice(); - void get_subset(); + void loop_lattice(int); void add_molecule(double *, double * = NULL); int vartest(double *); // evaluate a variable with new atom position }; diff --git a/src/random_mars.cpp b/src/random_mars.cpp index a832b3640f..e43c054977 100644 --- a/src/random_mars.cpp +++ b/src/random_mars.cpp @@ -21,6 +21,8 @@ using namespace LAMMPS_NS; +enum{ADD,SUBTRACT}; + /* ---------------------------------------------------------------------- */ RanMars::RanMars(LAMMPS *lmp, int seed) : Pointers(lmp), @@ -149,7 +151,7 @@ double RanMars::besselexp(double theta, double alpha, double cp) { double first,v1,v2; - if (theta < 0.0 || alpha < 0.0 || alpha > 1.0) + if (theta < 0.0 || alpha < 0.0 || alpha < 1.0) error->all(FLERR,"Invalid Bessel exponential distribution parameters"); v1 = uniform(); @@ -166,3 +168,131 @@ double RanMars::besselexp(double theta, double alpha, double cp) return first; } + +/* ---------------------------------------------------------------------- + select random subset of size Ntarget out of Ntotal items + Ntotal = sum of Nmine across all procs + mark,next = vectors of length Nmine + return mark = 0 for unselected item, 1 for selected item + next = work vector used to store linked lists for 2 active sets of items + IMPORTANT: this method must be called simultaneously by all procs +------------------------------------------------------------------------- */ + +void RanMars::select_subset(bigint ntarget, int nmine, int *mark, int *next) +{ + int mode,index,oldindex,newvalue,nflip,which,niter; + int active[2],first[2],last[2]; + int newactive[2],newfirst[2],newlast[2]; + bigint nmark,nactive,nactiveall,nflipall,bnflip; + bigint activeall[2],bsum[3],bsumall[3]; + double thresh; + + active[0] = nmine; + active[1] = 0; + first[0] = 0; + first[1] = -1; + last[0] = nmine-1; + last[1] = -1; + + bigint bnmine = nmine; + bigint bnall; + MPI_Allreduce(&bnmine,&bnall,1,MPI_LMP_BIGINT,MPI_SUM,world); + activeall[0] = bnall; + + for (int i = 0; i < nmine; i++) mark[i] = 0; + for (int i = 0; i < nmine; i++) next[i] = i+1; + next[nmine-1] = -1; + + nmark = 0; + niter = 0; + + while (nmark != ntarget) { + + // choose to ADD or SUBTRACT from current nmark + // thresh = desired flips / size of active set + // nactive = size of current active set, only for debug output below + + if (ntarget-nmark > 0) { + mode = ADD; + // nactive = active[mode]; + thresh = 1.0 * (ntarget-nmark) / activeall[mode]; + } else { + mode = SUBTRACT; + // nactive = active[mode]; + thresh = 1.0 * (nmark-ntarget) / activeall[mode]; + } + + // bound thresh for RNG accuracy + + thresh = MAX(thresh,0.01); + thresh = MIN(thresh,0.99); + + // new empty active sets for next iteration + + newactive[0] = newactive[1] = 0; + newfirst[0] = newfirst[1] = -1; + newlast[0] = newlast[1] = -1; + + // index = first value in ADD or SUBTRACT set + + if (mode == ADD) newvalue = 1; + else if (mode == SUBTRACT) newvalue = 0; + index = first[mode]; + + // flip marks from 0 -> 1 (ADD) or 1 -> 0 (SUBTRACT) + // loop over active set via next vector = linked list + // flip each value based on RN < thresh + + nflip = 0; + while (index >= 0) { + if (uniform() < thresh) { + mark[index] = newvalue; + nflip++; + } + oldindex = index; + index = next[index]; + + // oldindex can now be appended to a new active set + // which = which of two new active sets to append to + + which = mark[oldindex]; + newactive[which]++; + if (newfirst[which] < 0) newfirst[which] = oldindex; + if (newlast[which] >= 0) next[newlast[which]] = oldindex; + newlast[which] = oldindex; + next[oldindex] = -1; + + // set active sets for next iteration to the new ones + // next vector is already updated + + active[0] = newactive[0]; + active[1] = newactive[1]; + first[0] = newfirst[0]; + first[1] = newfirst[1]; + last[0] = newlast[0]; + last[1] = newlast[1]; + } + + // update nmark and activeall + + bsum[0] = nflip; + bsum[1] = active[0]; + bsum[2] = active[1]; + bsum[3] = nactive; + MPI_Allreduce(&bsum,&bsumall,4,MPI_LMP_BIGINT,MPI_SUM,world); + nflipall = bsumall[0]; + activeall[0] = bsumall[1]; + activeall[1] = bsumall[2]; + nactiveall = bsumall[3]; + + if (mode == ADD) nmark += nflipall; + else if (mode == SUBTRACT) nmark -= nflipall; + + niter++; + + // DEBUG output of stats + + //if (comm->me == 0) printf("%d %ld %ld %g %ld\n", + // niter,nmark,nactiveall,thresh,nflipall); + } +} diff --git a/src/random_mars.h b/src/random_mars.h index af653e201b..03199c84e8 100644 --- a/src/random_mars.h +++ b/src/random_mars.h @@ -27,6 +27,7 @@ class RanMars : protected Pointers { double gaussian(double mu, double sigma); double rayleigh(double sigma); double besselexp(double theta, double alpha, double cp); + void select_subset(bigint, int, int *, int *); private: int save; diff --git a/src/read_data.cpp b/src/read_data.cpp index d558b87633..5d45e8fe05 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -440,6 +440,12 @@ void ReadData::command(int narg, char **arg) atom->allocate_type_arrays(); atom->deallocate_topology(); + + // allocate atom arrays to N, rounded up by increment of DELTA + + bigint nbig = n; + nbig = atom->avec->roundup(nbig); + n = static_cast (nbig); atom->avec->grow(n); domain->boxlo[0] = boxlo[0]; domain->boxhi[0] = boxhi[0]; diff --git a/src/read_restart.cpp b/src/read_restart.cpp index a31700990d..583343ef51 100644 --- a/src/read_restart.cpp +++ b/src/read_restart.cpp @@ -165,6 +165,12 @@ void ReadRestart::command(int narg, char **arg) atom->allocate_type_arrays(); atom->deallocate_topology(); + + // allocate atom arrays to size N, rounded up by AtomVec->DELTA + + bigint nbig = n; + nbig = atom->avec->roundup(nbig); + n = static_cast (nbig); atom->avec->grow(n); n = atom->nmax; @@ -211,13 +217,17 @@ void ReadRestart::command(int narg, char **arg) memory->create(buf,assignedChunkSize,"read_restart:buf"); mpiio->read((headerOffset+assignedChunkOffset),assignedChunkSize,buf); mpiio->close(); - if (!nextra) { // We can actually calculate number of atoms from assignedChunkSize + + // can calculate number of atoms from assignedChunkSize + + if (!nextra) { atom->nlocal = 1; // temporarily claim there is one atom... int perAtomSize = avec->size_restart(); // ...so we can get its size atom->nlocal = 0; // restore nlocal to zero atoms int atomCt = (int) (assignedChunkSize / perAtomSize); if (atomCt > atom->nmax) avec->grow(atomCt); } + m = 0; while (m < assignedChunkSize) m += avec->unpack_restart(&buf[m]); } diff --git a/src/replicate.cpp b/src/replicate.cpp index 1617ab0313..e6aefe49ad 100644 --- a/src/replicate.cpp +++ b/src/replicate.cpp @@ -223,6 +223,12 @@ void Replicate::command(int narg, char **arg) else n = static_cast (LB_FACTOR * atom->natoms / nprocs); atom->allocate_type_arrays(); + + // allocate atom arrays to size N, rounded up by AtomVec->DELTA + + bigint nbig = n; + nbig = atom->avec->roundup(nbig); + n = static_cast (nbig); atom->avec->grow(n); n = atom->nmax; diff --git a/src/set.cpp b/src/set.cpp index 3bf00063c5..9807b07ea4 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -30,6 +30,7 @@ #include "input.h" #include "variable.h" #include "random_park.h" +#include "random_mars.h" #include "math_extra.h" #include "math_const.h" #include "memory.h" @@ -41,7 +42,8 @@ using namespace MathConst; enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; -enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI, +enum{TYPE,TYPE_FRACTION,TYPE_RATIO,TYPE_SUBSET, + MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE,LENGTH,TRI, DIPOLE,DIPOLE_RANDOM,SPIN,SPIN_RANDOM,QUAT,QUAT_RANDOM, THETA,THETA_RANDOM,ANGMOM,OMEGA, DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER, @@ -109,6 +111,34 @@ void Set::command(int narg, char **arg) setrandom(TYPE_FRACTION); iarg += 4; + } else if (strcmp(arg[iarg],"type/ratio") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal set command"); + newtype = force->inumeric(FLERR,arg[iarg+1]); + fraction = force->numeric(FLERR,arg[iarg+2]); + ivalue = force->inumeric(FLERR,arg[iarg+3]); + if (newtype <= 0 || newtype > atom->ntypes) + error->all(FLERR,"Invalid value in set command"); + if (fraction < 0.0 || fraction > 1.0) + error->all(FLERR,"Invalid value in set command"); + if (ivalue <= 0) + error->all(FLERR,"Invalid random number seed in set command"); + setrandom(TYPE_RATIO); + iarg += 4; + + } else if (strcmp(arg[iarg],"type/subset") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal set command"); + newtype = force->inumeric(FLERR,arg[iarg+1]); + nsubset = force->bnumeric(FLERR,arg[iarg+2]); + ivalue = force->inumeric(FLERR,arg[iarg+3]); + if (newtype <= 0 || newtype > atom->ntypes) + error->all(FLERR,"Invalid value in set command"); + if (nsubset < 0) + error->all(FLERR,"Invalid value in set command"); + if (ivalue <= 0) + error->all(FLERR,"Invalid random number seed in set command"); + setrandom(TYPE_SUBSET); + iarg += 4; + } else if (strcmp(arg[iarg],"mol") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1); @@ -1001,23 +1031,72 @@ void Set::setrandom(int keyword) AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri"); AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body"); - RanPark *random = new RanPark(lmp,1); double **x = atom->x; int seed = ivalue; - // set fraction of atom types to newtype + RanPark *ranpark = new RanPark(lmp,1); + RanMars *ranmars = new RanMars(lmp,seed + comm->me); + + // set approx fraction of atom types to newtype if (keyword == TYPE_FRACTION) { int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) if (select[i]) { - random->reset(seed,x[i]); - if (random->uniform() > fraction) continue; + ranpark->reset(seed,x[i]); + if (ranpark->uniform() > fraction) continue; atom->type[i] = newtype; count++; } + // set exact count of atom types to newtype + // for TYPE_RATIO, exact = fraction out of total eligible + // for TYPE_SUBSET, exact = nsubset out of total eligible + + } else if (keyword == TYPE_RATIO || keyword == TYPE_SUBSET) { + int nlocal = atom->nlocal; + + // count = number of eligible atoms I own + + count = 0; + for (i = 0; i < nlocal; i++) + if (select[i]) count++; + + // convert specified fraction to nsubset + + bigint bcount = count; + bigint allcount; + MPI_Allreduce(&bcount,&allcount,1,MPI_LMP_BIGINT,MPI_SUM,world); + + if (keyword == TYPE_RATIO) { + nsubset = static_cast (fraction * allcount); + } else if (keyword == TYPE_SUBSET) { + if (nsubset > allcount) + error->all(FLERR,"Set type/subset value exceeds eligible atoms"); + } + + // make selection + + int *flag = memory->create(flag,count,"set:flag"); + int *work = memory->create(work,count,"set:work"); + + ranmars->select_subset(nsubset,count,flag,work); + + // change types of selected atoms + + count = 0; + for (i = 0; i < nlocal; i++) + if (select[i] && flag[i]) { + atom->type[i] = newtype; + count++; + } + + // clean up + + memory->destroy(flag); + memory->destroy(work); + // set dipole moments to random orientations in 3d or 2d // dipole length is determined by dipole type array @@ -1030,10 +1109,10 @@ void Set::setrandom(int keyword) if (domain->dimension == 3) { for (i = 0; i < nlocal; i++) if (select[i]) { - random->reset(seed,x[i]); - mu[i][0] = random->uniform() - 0.5; - mu[i][1] = random->uniform() - 0.5; - mu[i][2] = random->uniform() - 0.5; + ranpark->reset(seed,x[i]); + mu[i][0] = ranpark->uniform() - 0.5; + mu[i][1] = ranpark->uniform() - 0.5; + mu[i][2] = ranpark->uniform() - 0.5; msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]; scale = dvalue/sqrt(msq); mu[i][0] *= scale; @@ -1046,9 +1125,9 @@ void Set::setrandom(int keyword) } else { for (i = 0; i < nlocal; i++) if (select[i]) { - random->reset(seed,x[i]); - mu[i][0] = random->uniform() - 0.5; - mu[i][1] = random->uniform() - 0.5; + ranpark->reset(seed,x[i]); + mu[i][0] = ranpark->uniform() - 0.5; + mu[i][1] = ranpark->uniform() - 0.5; mu[i][2] = 0.0; msq = mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1]; scale = dvalue/sqrt(msq); @@ -1072,10 +1151,10 @@ void Set::setrandom(int keyword) if (domain->dimension == 3) { for (i = 0; i < nlocal; i++) if (select[i]) { - random->reset(seed,x[i]); - sp[i][0] = random->uniform() - 0.5; - sp[i][1] = random->uniform() - 0.5; - sp[i][2] = random->uniform() - 0.5; + ranpark->reset(seed,x[i]); + sp[i][0] = ranpark->uniform() - 0.5; + sp[i][1] = ranpark->uniform() - 0.5; + sp[i][2] = ranpark->uniform() - 0.5; sp_sq = sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1] + sp[i][2]*sp[i][2]; scale = 1.0/sqrt(sp_sq); sp[i][0] *= scale; @@ -1088,9 +1167,9 @@ void Set::setrandom(int keyword) } else { for (i = 0; i < nlocal; i++) if (select[i]) { - random->reset(seed,x[i]); - sp[i][0] = random->uniform() - 0.5; - sp[i][1] = random->uniform() - 0.5; + ranpark->reset(seed,x[i]); + sp[i][0] = ranpark->uniform() - 0.5; + sp[i][1] = ranpark->uniform() - 0.5; sp[i][2] = 0.0; sp_sq = sp[i][0]*sp[i][0] + sp[i][1]*sp[i][1]; scale = 1.0/sqrt(sp_sq); @@ -1120,12 +1199,12 @@ void Set::setrandom(int keyword) else error->one(FLERR,"Cannot set quaternion for atom that has none"); - random->reset(seed,x[i]); - s = random->uniform(); + ranpark->reset(seed,x[i]); + s = ranpark->uniform(); t1 = sqrt(1.0-s); t2 = sqrt(s); - theta1 = 2.0*MY_PI*random->uniform(); - theta2 = 2.0*MY_PI*random->uniform(); + theta1 = 2.0*MY_PI*ranpark->uniform(); + theta2 = 2.0*MY_PI*ranpark->uniform(); quat[0] = cos(theta2)*t2; quat[1] = sin(theta1)*t1; quat[2] = cos(theta1)*t1; @@ -1144,8 +1223,8 @@ void Set::setrandom(int keyword) else error->one(FLERR,"Cannot set quaternion for atom that has none"); - random->reset(seed,x[i]); - theta2 = MY_PI*random->uniform(); + ranpark->reset(seed,x[i]); + theta2 = MY_PI*ranpark->uniform(); quat[0] = cos(theta2); quat[1] = 0.0; quat[2] = 0.0; @@ -1162,14 +1241,15 @@ void Set::setrandom(int keyword) if (select[i]) { if (atom->line[i] < 0) error->one(FLERR,"Cannot set theta for atom that is not a line"); - random->reset(seed,x[i]); - avec_line->bonus[atom->line[i]].theta = MY_2PI*random->uniform(); + ranpark->reset(seed,x[i]); + avec_line->bonus[atom->line[i]].theta = MY_2PI*ranpark->uniform(); count++; } } } - delete random; + delete ranpark; + delete ranmars; } /* ---------------------------------------------------------------------- */ diff --git a/src/set.h b/src/set.h index e6ba2a3bce..6788849677 100644 --- a/src/set.h +++ b/src/set.h @@ -34,8 +34,9 @@ class Set : protected Pointers { int *select; int style,ivalue,newtype,count,index_custom; int ximage,yimage,zimage,ximageflag,yimageflag,zimageflag; - double dvalue,xvalue,yvalue,zvalue,wvalue,fraction; int cc_index; + bigint nsubset; + double dvalue,xvalue,yvalue,zvalue,wvalue,fraction; int varflag,varflag1,varflag2,varflag3,varflag4; int ivar1,ivar2,ivar3,ivar4;