faster, simpler, 'more completely random' way to do this

src/create_atoms.cpp

@@ -56,7 +56,7 @@ void CreateAtoms::command(int narg, char **arg)
 {
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
-  
+
   if (domain->box_exist == 0)
     error->all(FLERR,"Create_atoms command before simulation box is defined");
   if (modify->nfix_restart_peratom)
@@ -110,13 +110,13 @@ void CreateAtoms::command(int narg, char **arg)
   remapflag = 0;
   mode = ATOM;
   int molseed;
-  int nboxseed;
+  int insertseed;
   varflag = 0;
   vstr = xstr = ystr = zstr = NULL;
   quatone[0] = quatone[1] = quatone[2] = 0.0;
-  nlattpts = created_Nmask = nbox = nboxflag = 0;
+  nlattpts = ninsert = insertflag = 0;
   Nmask = NULL;
-  
+
   nbasis = domain->lattice->nbasis;
   basistype = new int[nbasis];
   for (int i = 0; i < nbasis; i++) basistype[i] = ntype;
@@ -200,9 +200,9 @@ void CreateAtoms::command(int narg, char **arg)
       iarg += 5;
     } else if (strcmp(arg[iarg],"insert") == 0) {
       if (iarg+3 > narg) error->all(FLERR,"Illegal create_atoms command");
-      nbox = force->inumeric(FLERR,arg[iarg+1]);
+      ninsert = force->inumeric(FLERR,arg[iarg+1]);
-      nboxseed = force->inumeric(FLERR,arg[iarg+2]);
+      insertseed = force->inumeric(FLERR,arg[iarg+2]);
-      nboxflag = 1;
+      insertflag = 1;
       iarg += 3;
     } else error->all(FLERR,"Illegal create_atoms command");
   }
@@ -241,10 +241,10 @@ void CreateAtoms::command(int narg, char **arg)
     ranmol = new RanMars(lmp,molseed+comm->me);
   }
 
-  if (nboxflag) {
+  if (insertflag) {
-    ranbox = new RanMars(lmp,nboxseed+comm->me);
+    ranbox = new RanMars(lmp,insertseed+comm->me);
   }
-  
+
   // error check and further setup for variable test
 
   if (!vstr && (xstr || ystr || zstr))
@@ -381,7 +381,7 @@ void CreateAtoms::command(int narg, char **arg)
   else add_lattice();
 
   // init per-atom fix/compute/variable values for created atoms
-  
+
   atom->data_fix_compute_variable(nlocal_previous,atom->nlocal);
 
   // set new total # of atoms and error check
@@ -542,7 +542,7 @@ void CreateAtoms::command(int narg, char **arg)
   delete [] ystr;
   delete [] zstr;
   delete [] Nmask;
-  
+
   // print status
 
   if (comm->me == 0) {
@@ -764,7 +764,7 @@ void CreateAtoms::add_lattice()
   double *coord;
 
   int i,j,k,m;
-  
+
   // first pass: count how many particles will be inserted
   // second pass: filter to N number of particles (and insert)
   int maskcntr = 0;
@@ -774,35 +774,35 @@ void CreateAtoms::add_lattice()
       for (j = jlo; j <= jhi; j++)
         for (i = ilo; i <= ihi; i++)
           for (m = 0; m < nbasis; m++) {
-    
+
             x[0] = i + basis[m][0];
             x[1] = j + basis[m][1];
             x[2] = k + basis[m][2];
-    
+
             // convert from lattice coords to box coords
-    
+
             domain->lattice->lattice2box(x[0],x[1],x[2]);
-    
+
             // if a region was specified, test if atom is in it
-    
+
             if (style == REGION)
               if (!domain->regions[nregion]->match(x[0],x[1],x[2])) continue;
-    
+
             // 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;
-    
+
             // add the atom or entire molecule to my list of atoms
             if (pass == 0) nlattpts++;
             else {
@@ -821,72 +821,79 @@ void CreateAtoms::add_lattice()
 void CreateAtoms::lattice_mask()
 {
 
-  // Nmask will be used to choose which lattice points to insert
+  if (ninsert == 0 && insertflag && me == 0)
-  Nmask = new int[nlattpts];
+    error->warning(FLERR,"Specifying an 'insert' value of '0' is equivalent to no 'insert' keyword");
-  for (int i = 0; i < nlattpts; i++) 
+
-    Nmask[i] = 1;
+  if (ninsert > 0) {
-  
+    int nboxme;
-  if (nbox == 0 && nboxflag && me == 0) error->warning(FLERR,"Specifying an 'insert' value of '0' is equivalent to no 'insert' keyword");
-  
-  int nboxme = 0;
-  if (nbox > 0) {
     if (nprocs > 1) {
       int *allnlattpts = new int[nprocs]();
-      
+      int allnboxmes[nprocs];
+      int cumsum_lattpts[nprocs];
+      for (size_t i = 0; i < nprocs; i++)
+        allnboxmes[i] = 0;
+
       MPI_Allgather(&nlattpts, 1, MPI_INT, allnlattpts, 1, MPI_INT, world);
-      
+
-      int total_lattpts = 0;
+      if (me == 0) {
-      for (int i = 0; i < nprocs; i++)
+        int total_lattpts = 0;
-        total_lattpts += allnlattpts[i];
+        for (int i = 0; i < nprocs; i++) {
-      
+          total_lattpts += allnlattpts[i];
-      if (nbox > total_lattpts) error->all(FLERR,"Attempting to insert more particles than available lattice points");
+          cumsum_lattpts[i] = total_lattpts;
-      if (nbox < 0) error->all(FLERR,"Cannot insert a negative number of particles (in this universe)");
-      
-      // adjust nboxme based on personal number of eligible lattice points
-      
-      nboxme = round(nbox*nlattpts/total_lattpts);
-      
-      // limit nboxme to available lattice points on this processor
-      
-      if (nboxme > nlattpts) nboxme = nlattpts;
-      
-      // readjust so that all nboxme's add to exactly nbox
-      
-      int sum_nboxme;
-      MPI_Allreduce(&nboxme,&sum_nboxme,1,MPI_INT,MPI_SUM,world);
-      
-      if (nbox != sum_nboxme) {
-      
-        int diff = nbox - sum_nboxme;
-        int diff_sign = 1;
-        if (diff < 0) {
-          diff_sign = -1;
-          diff = -diff;
         }
-        
+
-        // go around adding or subtracting one from all processors (if there's space) until equal
+
-        int which_proc = 0;
+        if (ninsert > total_lattpts)
-        int success_me = 0, success_all = 0;
+          error->all(FLERR,"Attempting to insert more particles than available lattice points");
-        int failed_me = 0, failed_all = 0;
+        if (ninsert < 0)
-        while (success_all < diff && failed_all < nprocs) {
+          error->all(FLERR,"Cannot insert a negative number of particles (in this universe)");
-          if (me == which_proc++ && !failed_me)
+
-            if (nboxme < nlattpts) {
+        // using proc 0, let's insert N particles onto available lattice points by instead
-              nboxme += diff_sign;
+        // poking [nlattpts - N] holes into the lattice, randomly
-              success_me++;
+        int allNmask[total_lattpts];
-            } else {
+        for (int i = 0; i < total_lattpts; i++)
-              failed_me == 1;
+          allNmask[i] = 1;
+        int nholes = total_lattpts - ninsert;
+        int noptions = total_lattpts;
+        for (int i = 0; i < nholes; i++) {
+          int hindex = ceil(ranbox->uniform()*noptions);
+          int optcount = 0;
+          for (int j = 0; j < total_lattpts; j++) {
+            if (allNmask[j] == 1) {
+              optcount++;
+              if (optcount == hindex) {
+                allNmask[j] = 0;
+                noptions--;
+                break;
+              }
             }
-          MPI_Allreduce(&success_me,&success_all,1,MPI_INT,MPI_SUM,world);
+          }
-          MPI_Allreduce(&failed_me,&failed_all,1,MPI_INT,MPI_SUM,world);
+        }
-          if (which_proc > nprocs)
+
-            which_proc = 0;
+        // tell individual procs their nboxme (personal # to insert)
+        int iproc = 0;
+        while (allnlattpts[iproc] == 0 && iproc < nprocs)
+          iproc++;
+        for (int i = 0; i < total_lattpts; i++) {
+          if (i == cumsum_lattpts[iproc]) {
+            iproc++;
+            while (allnlattpts[iproc] == 0 && iproc < nprocs)
+              iproc++;
+          }
+          allnboxmes[iproc] += allNmask[i];
         }
-        // should have guaranteed success at this point
-        if (failed_all == nprocs) error->warning(FLERR,"This is an uncaught error. Ask developer");
       }
+      MPI_Scatter(&allnboxmes, 1, MPI_INT, &nboxme, 1, MPI_INT, 0, world);
       delete [] allnlattpts;
-    } else nboxme = nbox;
+    } else nboxme = ninsert;
-    
+
+    // probably faster to have individual processors 're-choose' their random points
+    // Nmask will be used to indicate which lattice points to insert
+    if (nlattpts < nboxme) printf("WHOAA\n"); //debug
+    Nmask = new int[nlattpts];
+    for (int i = 0; i < nlattpts; i++)
+      Nmask[i] = 1;
+
     // let's insert N particles onto available lattice points by instead
     // poking [nlattpts - N] holes into the lattice, randomly
     int nholes = nlattpts - nboxme;
@@ -901,14 +908,12 @@ void CreateAtoms::lattice_mask()
             Nmask[j] = 0;
             noptions--;
             break;
-          }  
+          }
         }
       }
     }
   }
-  
+
-  created_Nmask = 1;
-  
 }
 
 /* ----------------------------------------------------------------------

src/create_atoms.h

@@ -43,7 +43,7 @@ class CreateAtoms : protected Pointers {
   class Molecule *onemol;
   class RanMars *ranmol;
   class RanMars *ranbox;
-  
+
   int triclinic;
   double sublo[3],subhi[3];   // epsilon-extended proc sub-box for adding atoms
 
@@ -54,7 +54,7 @@ class CreateAtoms : protected Pointers {
   void add_molecule(double *, double * = NULL);
   int nlattpts;                // number of eligible lattice points
   int *Nmask;                  // used to insert N number of particles on lattice
-  int created_Nmask,nbox,nboxflag; 
+  int ninsert,insertflag;
   int vartest(double *);       // evaluate a variable with new atom position
 };