create_atoms subset: make 'fully parallel', fix bugs and typos

src/create_atoms.cpp

@@ -839,65 +839,61 @@ void CreateAtoms::add_lattice()
 void CreateAtoms::get_subset()
 {
   enum{ATOMS,HOLES};
-  int i,j,temp,irand,mysubset,npicks,pickmode;
+  int i,j,temp,irand,npicks,mynpicks,pickmode;
   double myrand;
-  if (nprocs > 1) {
+  int allnlatts[nprocs];
-    int allnlatts[nprocs];
+  int allsubsets[nprocs];
-    int allsubsets[nprocs];
+  int localpicks[nprocs];
-    double proc_sects[nprocs];
+  double proc_sects[nprocs];
 
-    MPI_Allgather(&nlatt, 1, MPI_INT, &allnlatts, 1, MPI_INT, world);
+  MPI_Allgather(&nlatt, 1, MPI_INT, &allnlatts, 1, MPI_INT, world);
 
-    if (me == 0) {
+  int ntotal = 0;
-      int ntotal = 0;
+  for (i = 0; i < nprocs; i++)
-      for (i = 0; i < nprocs; i++)
+    ntotal += allnlatts[i];
-        ntotal += allnlatts[i];
 
-      if (nsubset > ntotal)
+  if (nsubset > ntotal)
-         error->one(FLERR,"Attempting to insert more particles than available lattice points");
+     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
+  // define regions of unity based on a proc's fraction of total lattice points
-      proc_sects[0] = allnlatts[0]/ntotal;
+  proc_sects[0] = (double) allnlatts[0] / (double) ntotal;
-      for (i = 1; i < nprocs; i++)
+  for (i = 1; i < nprocs; i++)
-        proc_sects[i] = proc_sects[i-1] + (double) allnlatts[i] / (double) ntotal;
+    proc_sects[i] = proc_sects[i-1] + (double) allnlatts[i] / (double) ntotal;
 
-      if (nsubset > ntotal/2) {
+  if (nsubset > ntotal/2) {
-        pickmode = HOLES;
+    pickmode = HOLES;
-        npicks = ntotal - nsubset;
+    npicks = ntotal - nsubset;
-      } else {
-        pickmode = ATOMS;
-        npicks = nsubset;
-      }
-
-      for (i = 0; i < nprocs; i++)
-        allsubsets[i] = 0;
-
-      for (i = 0; i < npicks; i++) {
-        myrand = ranlatt->uniform();
-        for (j = 0; j < nprocs; j++)
-          if (myrand < proc_sects[j]) {
-            allsubsets[j]++;
-            break;
-          }
-      }
-
-      if (pickmode == HOLES)
-        for (i = 0; i < nprocs; i++)
-          allsubsets[i] = allnlatts[i] - allsubsets[i];
-    }
-
-    MPI_Scatter(&allsubsets, 1, MPI_INT, &mysubset, 1, MPI_INT, 0, world);
   } else {
-    if (nsubset > nlatt)
+    pickmode = ATOMS;
-       error->one(FLERR,"Attempting to insert more particles than available lattice points");
+    npicks = nsubset;
-    mysubset = 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];
+
   // each processor chooses its random lattice points
   memory->create(flag,nlatt,"create_atoms:flag");
 
   for (i = 0; i < nlatt; i++)
-    if (i < mysubset)
+    if (i < allsubsets[me])
       flag[i] = 1;
     else
       flag[i] = 0;

src/create_atoms.h

@@ -162,7 +162,7 @@ E: Attempting to insert more particles than available lattice points
 
 Self-explanatory.
 
-W: Specifying an 'randnpos' value of '0' is equivalent to no 'randnpos' keyword
+W: Specifying an 'subset' value of '0' is equivalent to no 'subset' keyword
 
 Self-explanatory.