modified versions of creating atoms on subset of lattice, ditto for set type/fraction

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
+The *ratio* and *subset* keywords can be used in conjunction with the
-*region* styles to limit the total number of particles inserted. The
+*box* or *region* styles to limit the total number of particles
-specified number of particles N are placed randomly on the available
+inserted.  The lattice defines a set of *Nlatt* eligible sites for
-lattice points.
+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

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 <variable>` to be used as the specified value(s).  If the
+keywords allow an :doc:`atom-style or atomfile-style variable
-value is a variable, it should be specified as v\_name, where name is
+<variable>` to be used as the specified value(s).  If the value is a
-the variable name.  In this case, the variable will be evaluated, and
+variable, it should be specified as v\_name, where name is the
-its resulting per-atom value used to determine the value assigned to
+variable name.  In this case, the variable will be evaluated, and its
-each selected atom.  Note that the per-atom value from the variable
+resulting per-atom value used to determine the value assigned to each
-will be ignored for atoms that are not selected via the *style* and
+selected atom.  Note that the per-atom value from the variable will be
-*ID* settings explained above.  A simple way to use per-atom values
+ignored for atoms that are not selected via the *style* and *ID*
-from the variable to reset a property for all atoms is to use style
+settings explained above.  A simple way to use per-atom values from
-*atom* with *ID* = "\*"; this selects all atom IDs.
+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 <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
+to be exactly the specified fraction (0 <= *fraction* <= 1), but
-close.  Random numbers are used in such a way that a particular atom
+should be statistically close.  Random numbers are used in such a way
-is changed or not changed, regardless of how many processors are being
+that a particular atom is changed or not changed, regardless of how
-used.  This keyword does not allow use of an atom-style variable.
+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 <atom_style>` being used must support the use of molecule
+Keywords *type/ratio* and *type/subset" also set the atom type for a
-IDs.
+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
+Keyword *mol* sets the molecule ID for all selected atoms.  The
-all selected atoms.  For *charge*\ , the :doc:`atom style <atom_style>`
+:doc:`atom style <atom_style>` being used must support the use of
-being used must support the use of atomic charge. Keywords *vx*\ , *vy*\ ,
+molecule IDs.
-and *vz* set the velocities of all selected atoms.
+
+Keywords *x*\ , *y*\ , *z*\ , and *charge* set the coordinates or
+charge of all selected atoms.  For *charge*\ , the :doc:`atom style
+<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
+selected atoms.  The magnitude of the dipole moment is set by the
-by the length of this orientation vector.
+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

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

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) {}

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;
+  subsetflag = NONE;
-  flag = NULL;
+  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;
+      if (nsubset <= 0 || subsetseed <= 0) 
-      else {
+        error->all(FLERR,"Illegal create_atoms command");
-        if (me == 0) error->warning(FLERR,"Specifying an 'subset' value of "
-                                 "'0' is equivalent to no 'subset' keyword");
-        subsetflag = 0;
-      }
       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 = NULL;
-    ranlatt = new RanMars(lmp,subsetseed+me);
+  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<int> (xmin) - 1;
   jlo = static_cast<int> (ymin) - 1;
   klo = static_cast<int> (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
+  // rough estimate of total time used for create atoms
-  // 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.
   // 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
+  // maxestimate = time in hours
-  testimate *= static_cast<double>(khi-klo+1);
+  
-  testimate *= static_cast<double>(jhi-jlo+1);
+  double estimate = 2.5e-8/3600.0;
-  testimate *= static_cast<double>(ihi-ilo+1);
+  estimate *= static_cast<double> (khi-klo+1);
-  testimate *= static_cast<double>(nbasis);
+  estimate *= static_cast<double> (jhi-jlo+1);
+  estimate *= static_cast<double> (ihi-ilo+1);
+  estimate *= static_cast<double> (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,16 +818,73 @@ void CreateAtoms::add_lattice()
                          "approx. %.2f hours to complete\n",maxestimate);
   }
 
+  // count lattice sites on each proc
+
+  nlatt_overflow = 0;
+  loop_lattice(COUNT);
+
+  // nadd = # of atoms each proc will insert (estimated if subsetflag)
+
+  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");
+
+  bigint nadd;
+
+  if (subsetflag == NONE) {
+    if (nprocs == 1) nadd = nlatt;
+    else nadd = static_cast<bigint> (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<bigint> (subsetfrac * bnlattall);
+    if (nsubset > bnlattall)
+      error->all(FLERR,"Create_atoms subset size > # of lattice sites");
+    if (nprocs == 1) nadd = nsubset;
+    else nadd = static_cast<bigint> (LB_FACTOR * nsubset/bnlattall * nlatt);
+  }
+  
+  // allocate atom arrays to size N, rounded up by AtomVec->DELTA
+
+  bigint nbig = atom->avec->roundup(nadd + atom->nlocal);
+  int n = static_cast<int> (nbig);
+  atom->avec->grow(n);
+
+  // add atoms or molecules
+  // if no subset: add to all lattice sites
+  // if subset: count lattice sites, select random subset, then add
+
+  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);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   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::loop_lattice(int action)
+{
   int i,j,k,m;
 
-  // one pass for default mode, two passes for subset mode:
+  const double * const * const basis = domain->lattice->basis;
-  // first pass: count how many particles will be inserted
+
-  // second pass: filter to N number of particles (and insert)
+  nlatt = 0;
-  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++) {
@@ -861,135 +920,30 @@ void CreateAtoms::add_lattice()
               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
+          // this proc owns the lattice site
-            if (subsetflag && pass == 0) nlatt++;
+          // perform action: add, just count, add if flagged
-            else {
+          // add = add an atom or entire molecule to my list of atoms
-              if (!subsetflag || flag[iflag++] == 1)
+
+          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);
           }
+
+          nlatt++;
         }
       }
     }
   }
 }
-}
-
-/* ----------------------------------------------------------------------
-   define a random mask to insert only N particles on lattice
-------------------------------------------------------------------------- */
-
-void CreateAtoms::get_subset()
-{
-  enum{ATOMS,HOLES};
-  int i,j,temp,irand,offlag,npicks,pickmode,mynpicks,mysubset;
-  double myrand;
-  int *allnlatts;
-  int *allsubsets;
-  int *localpicks;
-  double *proc_sects;
-
-  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");
-
-  MPI_Allgather(&nlatt, 1, MPI_INT, &allnlatts[0], 1, MPI_INT, world);
-
-  int ntotal = 0;
-  for (i = 0; i < nprocs; i++)
-    ntotal += allnlatts[i];
-
-  if (nsubset > ntotal)
-     error->all(FLERR,"Attempting to insert more particles than available lattice points");
-
-  // 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 (nsubset > ntotal/2) {
-    pickmode = HOLES;
-    npicks = ntotal - nsubset;
-  } else {
-    pickmode = ATOMS;
-    npicks = nsubset;
-  }
-
-  mynpicks = npicks/nprocs;
-  if (me == 0) mynpicks = npicks - (nprocs-1)*(mynpicks);
-
-  for (i = 0; i < nprocs; i++)
-    localpicks[i] = 0;
-
-  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--;
-        }
-      }
-    }
-    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

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 ilo,ihi,jlo,jhi,klo,khi;
-  int *flag;                  // used to insert N number of particles on lattice
+
+  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
 };

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);
+  }
+}

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;

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<int> (nbig);
       atom->avec->grow(n);
 
       domain->boxlo[0] = boxlo[0]; domain->boxhi[0] = boxhi[0];

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<int> (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]);
   }

src/replicate.cpp

@@ -223,6 +223,12 @@ void Replicate::command(int narg, char **arg)
   else n = static_cast<int> (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<int> (nbig);
   atom->avec->grow(n);
   n = atom->nmax;
 

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]);
+        ranpark->reset(seed,x[i]);
-        if (random->uniform() > fraction) continue;
+        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<bigint> (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]);
+          ranpark->reset(seed,x[i]);
-          mu[i][0] = random->uniform() - 0.5;
+          mu[i][0] = ranpark->uniform() - 0.5;
-          mu[i][1] = random->uniform() - 0.5;
+          mu[i][1] = ranpark->uniform() - 0.5;
-          mu[i][2] = random->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]);
+          ranpark->reset(seed,x[i]);
-          mu[i][0] = random->uniform() - 0.5;
+          mu[i][0] = ranpark->uniform() - 0.5;
-          mu[i][1] = random->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]);
+          ranpark->reset(seed,x[i]);
-          sp[i][0] = random->uniform() - 0.5;
+          sp[i][0] = ranpark->uniform() - 0.5;
-          sp[i][1] = random->uniform() - 0.5;
+          sp[i][1] = ranpark->uniform() - 0.5;
-          sp[i][2] = random->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]);
+          ranpark->reset(seed,x[i]);
-          sp[i][0] = random->uniform() - 0.5;
+          sp[i][0] = ranpark->uniform() - 0.5;
-          sp[i][1] = random->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]);
+          ranpark->reset(seed,x[i]);
-          s = random->uniform();
+          s = ranpark->uniform();
           t1 = sqrt(1.0-s);
           t2 = sqrt(s);
-          theta1 = 2.0*MY_PI*random->uniform();
+          theta1 = 2.0*MY_PI*ranpark->uniform();
-          theta2 = 2.0*MY_PI*random->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]);
+          ranpark->reset(seed,x[i]);
-          theta2 = MY_PI*random->uniform();
+          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]);
+        ranpark->reset(seed,x[i]);
-        avec_line->bonus[atom->line[i]].theta = MY_2PI*random->uniform();
+        avec_line->bonus[atom->line[i]].theta = MY_2PI*ranpark->uniform();
         count++;
       }
     }
   }
 
-  delete random;
+  delete ranpark;
+  delete ranmars;
 }
 
 /* ---------------------------------------------------------------------- */

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;