diff --git a/doc/src/create_atoms.rst b/doc/src/create_atoms.rst
index 2addd89fcb..2037c15518 100644
--- a/doc/src/create_atoms.rst
+++ b/doc/src/create_atoms.rst
@@ -28,7 +28,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 *ratio* or *subset* or *remap* or *var* or *set* or *radscale* or *meshmode* or *rotate* or *overlap* or *maxtry* or *units*
+* keyword = *mol* or *basis* or *ratio* or *subset* or *group* or *remap* or *var* or *set* or *radscale* or *meshmode* or *rotate* or *overlap* or *maxtry* or *units*
 
   .. parsed-literal::
 
@@ -44,6 +44,7 @@ Syntax
        *subset* values = Nsubset seed
          Nsubset = # of lattice sites to populate randomly
          seed = random # seed (positive integer)
+       *group* value = group name
        *remap* value = *yes* or *no*
        *var* value = name = variable name to evaluate for test of atom creation
        *set* values = dim name
@@ -83,7 +84,7 @@ Examples
 
    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 3 single 0 0 5 group newatom
    create_atoms 1 box var v set x xpos set y ypos
    create_atoms 2 random 50 12345 NULL overlap 2.0 maxtry 50
    create_atoms 1 mesh open_box.stl meshmode qrand 0.1 units box
@@ -395,6 +396,12 @@ correct number of particles are inserted, in a perfectly random
 fashion.  Which lattice sites are selected will change with the number
 of processors used.
 
+.. versionadded:: TBD
+
+The *group* keyword adds the newly created atoms to the named
+:doc:`group <group>`.  If the group does not yet exist it will be
+created.  All created atoms are always added to the group "all".
+
 The *remap* keyword only applies to the *single* style.  If it is set
 to *yes*, then if the specified position is outside the simulation
 box, it will mapped back into the box, assuming the relevant
diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp
index bac3a556eb..93f4519d3c 100644
--- a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@@ -23,6 +23,7 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
+#include "group.h"
 #include "input.h"
 #include "irregular.h"
 #include "lattice.h"
@@ -64,6 +65,16 @@ enum { NONE, RATIO, SUBSET };
 enum { BISECTION, QUASIRANDOM };
 
 static constexpr const char *mesh_name[] = {"recursive bisection", "quasi-random"};
+
+/* ---------------------------------------------------------------------- */
+
+CreateAtoms::CreateAtoms(LAMMPS *lmp) :
+    Command(lmp), basistype(nullptr), xmol(nullptr), vstr(nullptr), xstr(nullptr), ystr(nullptr),
+    zstr(nullptr), groupname(nullptr), region(nullptr), onemol(nullptr), ranmol(nullptr),
+    ranlatt(nullptr)
+{
+}
+
 /* ---------------------------------------------------------------------- */
 
 CreateAtoms::~CreateAtoms()
@@ -80,7 +91,6 @@ CreateAtoms::~CreateAtoms()
   delete ranmol;
   delete ranlatt;
 }
-CreateAtoms::CreateAtoms(LAMMPS *lmp) : Command(lmp), basistype(nullptr) {}
 
 /* ---------------------------------------------------------------------- */
 
@@ -179,8 +189,10 @@ void CreateAtoms::command(int narg, char **arg)
   mesh_density = 1.0;
 
   nbasis = domain->lattice->nbasis;
-  basistype = new int[nbasis];
-  for (int i = 0; i < nbasis; i++) basistype[i] = ntype;
+  if (nbasis > 0) {
+    basistype = new int[nbasis];
+    for (int i = 0; i < nbasis; i++) basistype[i] = ntype;
+  }
 
   while (iarg < narg) {
     if (strcmp(arg[iarg], "basis") == 0) {
@@ -273,6 +285,15 @@ void CreateAtoms::command(int narg, char **arg)
       if ((nsubset <= 0) || (subsetseed <= 0))
         error->all(FLERR, Error::NOPOINTER, "Illegal create_atoms subset settings");
       iarg += 3;
+    } else if (strcmp(arg[iarg], "group") == 0) {
+      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "create_atoms group", error);
+      if (strcmp(arg[iarg + 1], "none") == 0) {
+        delete[] groupname;
+        groupname = nullptr;
+      } else {
+        groupname = utils::strdup(arg[iarg + 1]);
+      }
+      iarg += 2;
     } else if (strcmp(arg[iarg], "overlap") == 0) {
       if (style != RANDOM)
         error->all(FLERR, iarg, "Create_atoms overlap can only be used with random style");
@@ -675,18 +696,6 @@ void CreateAtoms::command(int narg, char **arg)
     if (domain->triclinic) domain->lamda2x(atom->nlocal);
   }
 
-  // clean up
-
-  delete ranmol;
-  delete ranlatt;
-
-  delete[] basistype;
-  delete[] vstr;
-  delete[] xstr;
-  delete[] ystr;
-  delete[] zstr;
-  if (mode == MOLECULE) memory->destroy(xmol);
-
   // for MOLECULE mode:
   // create special bond lists for molecular systems,
   //   but not for atom style template
@@ -699,6 +708,13 @@ void CreateAtoms::command(int narg, char **arg)
     }
   }
 
+  // add atoms to group
+
+  if (groupname) {
+    int groupbit = group->bitmask[group->find_or_create(groupname)];
+    for (int i = nlocal_previous; i < atom->nlocal; ++i) atom->mask[i] |= groupbit;
+  }
+
   // print status
 
   MPI_Barrier(world);
diff --git a/src/create_atoms.h b/src/create_atoms.h
index ca2a8779a9..a1881f38b9 100644
--- a/src/create_atoms.h
+++ b/src/create_atoms.h
@@ -46,6 +46,7 @@ class CreateAtoms : public Command {
 
   int varflag, vvar, xvar, yvar, zvar;
   char *vstr, *xstr, *ystr, *zstr;
+  char *groupname;
 
   int ilo, ihi, jlo, jhi, klo, khi;