diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp
index 7a8bf630df..00455096f8 100644
--- a/src/GRANULAR/fix_pour.cpp
+++ b/src/GRANULAR/fix_pour.cpp
@@ -131,7 +131,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
 
       if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols)
         error->all(FLERR, "Fix pour molecule template ID must be same as atom style template ID");
-      onemols[i]->check_attributes(0);
+      onemols[i]->check_attributes();
 
       // fix pour uses geoemetric center of molecule for insertion
 
diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp
index 835d049bc6..0cc6590d5a 100644
--- a/src/MC/fix_gcmc.cpp
+++ b/src/MC/fix_gcmc.cpp
@@ -177,7 +177,7 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) :
     if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols)
       error->all(FLERR,"Fix gcmc molecule template ID must be same "
                  "as atom_style template ID");
-    onemols[imol]->check_attributes(0);
+    onemols[imol]->check_attributes();
   }
 
   if (charge_flag && atom->q == nullptr)
diff --git a/src/MC/fix_widom.cpp b/src/MC/fix_widom.cpp
index 0a20e7adf3..d1bc5dfa58 100644
--- a/src/MC/fix_widom.cpp
+++ b/src/MC/fix_widom.cpp
@@ -153,7 +153,7 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
     if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols)
       error->all(FLERR,"Fix widom molecule template ID must be same "
                  "as atom_style template ID");
-    onemols[imol]->check_attributes(0);
+    onemols[imol]->check_attributes();
   }
 
   if (charge_flag && atom->q == nullptr)
diff --git a/src/REACTION/fix_bond_react.cpp b/src/REACTION/fix_bond_react.cpp
index 0383449bbd..036959b14a 100644
--- a/src/REACTION/fix_bond_react.cpp
+++ b/src/REACTION/fix_bond_react.cpp
@@ -471,8 +471,8 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
     open(files[i]);
     onemol = atom->molecules[unreacted_mol[i]];
     twomol = atom->molecules[reacted_mol[i]];
-    onemol->check_attributes(0);
-    twomol->check_attributes(0);
+    onemol->check_attributes();
+    twomol->check_attributes();
     get_molxspecials();
     read(i);
     fclose(fp);
diff --git a/src/atom.cpp b/src/atom.cpp
index 91b72841f0..3b9a47918d 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -1900,6 +1900,18 @@ int Atom::find_molecule(char *id)
   return -1;
 }
 
+/* ----------------------------------------------------------------------
+   return vector of molecules which match template ID
+------------------------------------------------------------------------- */
+
+std::vector<Molecule *>Atom::get_molecule_by_id(const std::string &id)
+{
+  std::vector<Molecule *> result;
+  for (int imol = 0; imol < nmolecule; ++imol)
+    if (id == molecules[imol]->id) result.push_back(molecules[imol]);
+  return result;
+}
+
 /* ----------------------------------------------------------------------
    add info to current atom ilocal from molecule template onemol and its iatom
    offset = atom ID preceding IDs of atoms in this molecule
@@ -1912,8 +1924,7 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom, int ilocal, tagint off
   if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom];
   if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom];
   else if (rmass_flag)
-    rmass[ilocal] = 4.0*MY_PI/3.0 *
-      radius[ilocal]*radius[ilocal]*radius[ilocal];
+    rmass[ilocal] = 4.0*MY_PI/3.0 * radius[ilocal]*radius[ilocal]*radius[ilocal];
   if (onemol->bodyflag) {
     body[ilocal] = 0;     // as if a body read from data file
     onemol->avec_body->data_body(ilocal,onemol->nibody,onemol->ndbody,
@@ -1923,10 +1934,8 @@ void Atom::add_molecule_atom(Molecule *onemol, int iatom, int ilocal, tagint off
 
   // initialize custom per-atom properties to zero if present
 
-  for (int i = 0; i < nivector; ++i)
-    ivector[i][ilocal] = 0;
-  for (int i = 0; i < ndvector; ++i)
-    dvector[i][ilocal] = 0.0;
+  for (int i = 0; i < nivector; ++i) ivector[i][ilocal] = 0;
+  for (int i = 0; i < ndvector; ++i) dvector[i][ilocal] = 0.0;
   for (int i = 0; i < niarray; ++i)
     for (int j = 0; j < icols[i]; ++j)
       iarray[i][ilocal][j] = 0;
diff --git a/src/atom.h b/src/atom.h
index df43898dbf..b5b52d7421 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -345,6 +345,7 @@ class Atom : protected Pointers {
 
   void add_molecule(int, char **);
   int find_molecule(char *);
+  std::vector<Molecule *>get_molecule_by_id(const std::string &);
   void add_molecule_atom(class Molecule *, int, int, tagint);
 
   void first_reorder();
diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp
index da352e3ae5..85e83d694f 100644
--- a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@@ -299,7 +299,7 @@ void CreateAtoms::command(int narg, char **arg)
     if (onemol->tag_require && !atom->tag_enable)
       error->all(FLERR, "Create_atoms molecule has atom IDs, but system does not");
 
-    onemol->check_attributes(0);
+    onemol->check_attributes();
 
     // use geometric center of molecule for insertion
     // molecule random number generator, different for each proc
diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp
index fd6eb3e36a..be30fbc97b 100644
--- a/src/fix_deposit.cpp
+++ b/src/fix_deposit.cpp
@@ -117,7 +117,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
       if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols)
         error->all(FLERR,"Fix deposit molecule template ID must be same "
                    "as atom_style template ID");
-      onemols[i]->check_attributes(0);
+      onemols[i]->check_attributes();
 
       // fix deposit uses geoemetric center of molecule for insertion
 
diff --git a/src/molecule.cpp b/src/molecule.cpp
index 816b2686ab..702bee45af 100644
--- a/src/molecule.cpp
+++ b/src/molecule.cpp
@@ -1671,64 +1671,48 @@ int Molecule::findfragment(const char *name)
 
 /* ----------------------------------------------------------------------
    error check molecule attributes and topology against system settings
-   flag = 0, just check this molecule
-   flag = 1, check all molecules in set, this is 1st molecule in set
 ------------------------------------------------------------------------- */
 
-void Molecule::check_attributes(int flag)
+void Molecule::check_attributes()
 {
-  int n = 1;
-  if (flag) n = nset;
-  int imol = atom->find_molecule(id);
+  // check per-atom attributes of molecule
+  // warn if not a match
 
-  for (int i = imol; i < imol+n; i++) {
-    Molecule *onemol = atom->molecules[imol];
+  int mismatch = 0;
+  if (qflag && !atom->q_flag) mismatch = 1;
+  if (radiusflag && !atom->radius_flag) mismatch = 1;
+  if (rmassflag && !atom->rmass_flag) mismatch = 1;
 
-    // check per-atom attributes of molecule
-    // warn if not a match
+  if (mismatch && me == 0)
+    error->warning(FLERR,"Molecule attributes do not match system attributes");
 
-    int mismatch = 0;
-    if (onemol->qflag && !atom->q_flag) mismatch = 1;
-    if (onemol->radiusflag && !atom->radius_flag) mismatch = 1;
-    if (onemol->rmassflag && !atom->rmass_flag) mismatch = 1;
+  // for all atom styles, check nbondtype,etc
 
-    if (mismatch && me == 0)
-      error->warning(FLERR,"Molecule attributes do not match system attributes");
+  mismatch = 0;
+  if (atom->nbondtypes < nbondtypes) mismatch = 1;
+  if (atom->nangletypes < nangletypes) mismatch = 1;
+  if (atom->ndihedraltypes < ndihedraltypes) mismatch = 1;
+  if (atom->nimpropertypes < nimpropertypes) mismatch = 1;
 
-    // for all atom styles, check nbondtype,etc
+  if (mismatch) error->all(FLERR,"Molecule topology type exceeds system topology type");
 
-    mismatch = 0;
-    if (atom->nbondtypes < onemol->nbondtypes) mismatch = 1;
-    if (atom->nangletypes < onemol->nangletypes) mismatch = 1;
-    if (atom->ndihedraltypes < onemol->ndihedraltypes) mismatch = 1;
-    if (atom->nimpropertypes < onemol->nimpropertypes) mismatch = 1;
+  // for molecular atom styles, check bond_per_atom,etc + maxspecial
+  // do not check for atom style template, since nothing stored per atom
 
-    if (mismatch)
-      error->all(FLERR,"Molecule topology type exceeds system topology type");
+  if (atom->molecular == Atom::MOLECULAR) {
+    if (atom->avec->bonds_allow && atom->bond_per_atom < bond_per_atom) mismatch = 1;
+    if (atom->avec->angles_allow && atom->angle_per_atom < angle_per_atom) mismatch = 1;
+    if (atom->avec->dihedrals_allow && atom->dihedral_per_atom < dihedral_per_atom) mismatch = 1;
+    if (atom->avec->impropers_allow && atom->improper_per_atom < improper_per_atom) mismatch = 1;
+    if (atom->maxspecial < maxspecial) mismatch = 1;
 
-    // for molecular atom styles, check bond_per_atom,etc + maxspecial
-    // do not check for atom style template, since nothing stored per atom
-
-    if (atom->molecular == Atom::MOLECULAR) {
-      if (atom->avec->bonds_allow &&
-          atom->bond_per_atom < onemol->bond_per_atom) mismatch = 1;
-      if (atom->avec->angles_allow &&
-          atom->angle_per_atom < onemol->angle_per_atom) mismatch = 1;
-      if (atom->avec->dihedrals_allow &&
-          atom->dihedral_per_atom < onemol->dihedral_per_atom) mismatch = 1;
-      if (atom->avec->impropers_allow &&
-          atom->improper_per_atom < onemol->improper_per_atom) mismatch = 1;
-      if (atom->maxspecial < onemol->maxspecial) mismatch = 1;
-
-      if (mismatch)
-        error->all(FLERR,"Molecule topology/atom exceeds system topology/atom");
-    }
-
-    // warn if molecule topology defined but no special settings
-
-    if (onemol->bondflag && !onemol->specialflag)
-      if (me == 0) error->warning(FLERR,"Molecule has bond topology but no special bond settings");
+    if (mismatch) error->all(FLERR,"Molecule topology/atom exceeds system topology/atom");
   }
+
+  // warn if molecule topology defined but no special settings
+
+  if (bondflag && !specialflag)
+    if (me == 0) error->warning(FLERR,"Molecule has bond topology but no special bond settings");
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/molecule.h b/src/molecule.h
index 69106e5a34..b4a8e642cc 100644
--- a/src/molecule.h
+++ b/src/molecule.h
@@ -125,7 +125,7 @@ class Molecule : protected Pointers {
   void compute_com();
   void compute_inertia();
   int findfragment(const char *);
-  void check_attributes(int);
+  void check_attributes();
 
  private:
   int me;
diff --git a/src/read_data.cpp b/src/read_data.cpp
index 1ec8210e5e..74bb7e73d3 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -866,7 +866,10 @@ void ReadData::command(int narg, char **arg)
   // insure nbondtypes,etc are still consistent with template molecules,
   //   in case data file re-defined them
 
-  if (atom->molecular == Atom::TEMPLATE) atom->avec->onemols[0]->check_attributes(1);
+  if (atom->molecular == Atom::TEMPLATE) {
+    int nset = MAX(1, atom->avec->onemols[0]->nset);
+    for (int i = 0; i < nset; ++i) atom->avec->onemols[i]->check_attributes();
+  }
 
   // if adding atoms, migrate atoms to new processors
   // use irregular() b/c box size could have changed dramaticaly
diff --git a/unittest/formats/test_molecule_file.cpp b/unittest/formats/test_molecule_file.cpp
index 17bd30a349..d84aa1c960 100644
--- a/unittest/formats/test_molecule_file.cpp
+++ b/unittest/formats/test_molecule_file.cpp
@@ -216,6 +216,9 @@ TEST_F(MoleculeFileTest, twomols)
     auto output = END_CAPTURE_OUTPUT();
     ASSERT_THAT(output, ContainsRegex(".*Read molecule template.*\n.*2 molecules.*\n"
                                       ".*0 fragments.*\n.*2 atoms with max type 2.*\n.*0 bonds.*"));
+    ASSERT_EQ(lmp->atom->nmolecule, 1);
+    auto mols = lmp->atom->get_molecule_by_id(test_name);
+    ASSERT_EQ(mols.size(), 1);
 }
 
 TEST_F(MoleculeFileTest, twofiles)
@@ -231,6 +234,17 @@ TEST_F(MoleculeFileTest, twofiles)
                       ".*Read molecule template twomols:.*\n.*1 molecules.*\n"
                       ".*0 fragments.*\n.*3 atoms with max type 4.*\n.*2 bonds with max type 2.*\n"
                       ".*1 angles with max type 2.*\n.*0 dihedrals.*"));
+    BEGIN_CAPTURE_OUTPUT();
+    command("molecule h2o moltest.h2o.mol");
+    command("molecule co2 moltest.co2.mol");
+    output = END_CAPTURE_OUTPUT();
+    ASSERT_EQ(lmp->atom->nmolecule, 4);
+    auto mols = lmp->atom->get_molecule_by_id("twomols");
+    ASSERT_EQ(mols.size(), 2);
+    mols = lmp->atom->get_molecule_by_id("h2o");
+    ASSERT_EQ(mols.size(), 1);
+    mols = lmp->atom->get_molecule_by_id("co2");
+    ASSERT_EQ(mols.size(), 1);
 }
 
 TEST_F(MoleculeFileTest, bonds)