diff --git a/src/molecule.cpp b/src/molecule.cpp
index 5d71db68e6..644e3ff79c 100644
--- a/src/molecule.cpp
+++ b/src/molecule.cpp
@@ -55,14 +55,10 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
 
   if (index >= narg) error->all(FLERR,"Illegal molecule command");
 
-  int n = strlen(arg[0]) + 1;
-  id = new char[n];
-  strcpy(id,arg[0]);
-
-  for (int i = 0; i < n-1; i++)
-    if (!isalnum(id[i]) && id[i] != '_')
-      error->all(FLERR,"Molecule template ID must be "
-                 "alphanumeric or underscore characters");
+  id = utils::strdup(arg[0]);
+  if (!utils::is_id(id))
+    error->all(FLERR,"Molecule template ID must have only "
+               "alphanumeric or underscore characters");
 
   // parse args until reach unknown arg (next file)
 
@@ -129,16 +125,21 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
 
   initialize();
 
-  // scan file for sizes of all fields and allocate them
+  // scan file for sizes of all fields and allocate storage for them
 
-  if (me == 0) open(arg[ifile]);
+  if (me == 0) {
+    fp = fopen(arg[ifile],"r");
+    if (fp == nullptr)
+      error->one(FLERR,fmt::format("Cannot open molecule file {}: {}",
+                                   arg[ifile], utils::getsyserror()));
+  }
   read(0);
   if (me == 0) fclose(fp);
   allocate();
 
   // read file again to populate all fields
 
-  if (me == 0) open(arg[ifile]);
+  if (me == 0) fp = fopen(arg[ifile],"r");
   read(1);
   if (me == 0) fclose(fp);
 
@@ -397,8 +398,8 @@ void Molecule::compute_inertia()
 
 void Molecule::read(int flag)
 {
-  char line[MAXLINE],keyword[MAXLINE];
-  char *eof,*ptr;
+  char line[MAXLINE];
+  char *eof;
 
   // skip 1st line of file
 
@@ -415,15 +416,14 @@ void Molecule::read(int flag)
 
     readline(line);
 
-    // trim anything from '#' onward
-    // if line is blank, continue
+    // trim comments. if line is blank, continue
 
-    if ((ptr = strchr(line,'#'))) *ptr = '\0';
-    if (strspn(line," \t\n\r") == strlen(line)) continue;
+    auto text = utils::trim(utils::trim_comment(line));
+    if (text.empty()) continue;
 
     // search line for header keywords and set corresponding variable
     try {
-      ValueTokenizer values(line);
+      ValueTokenizer values(text);
 
       int nmatch = values.count();
       int nwant = 0;
@@ -485,10 +485,16 @@ void Molecule::read(int flag)
         nibody = values.next_int();
         ndbody = values.next_int();
         nwant = 3;
-      } else break;
-
+      } else {
+        // unknown header keyword
+        if (utils::strmatch(text,"^\\d+\\s+\\S+")) {
+          values.next_int();
+          auto keyword = values.next_string();
+          error->one(FLERR,fmt::format("Invalid header keyword: {}",keyword));
+        } else break;
+      }
       if (nmatch != nwant)
-        error->one(FLERR,"Invalid header in molecule file");
+        error->one(FLERR,"Invalid header line format in molecule file");
     } catch (TokenizerException &e) {
       error->one(FLERR, fmt::format("Invalid header in molecule file\n"
                                     "{}", e.what()));
@@ -498,7 +504,7 @@ void Molecule::read(int flag)
   // error checks
 
   if (natoms < 1)
-    error->all(FLERR,"No count or invalid atom count in molecule file");
+    error->all(FLERR,"No atoms or invalid atom count in molecule file");
   if (nbonds < 0) error->all(FLERR,"Invalid bond count in molecule file");
   if (nangles < 0) error->all(FLERR,"Invalid angle count in molecule file");
   if (ndihedrals < 0)
@@ -512,108 +518,113 @@ void Molecule::read(int flag)
 
   // grab keyword and skip next line
 
-  parse_keyword(0,line,keyword);
+  std::string keyword = parse_keyword(0,line);
   readline(line);
 
   // loop over sections of molecule file
 
-  while (strlen(keyword) > 0) {
-    if (strcmp(keyword,"Coords") == 0) {
+  while (!keyword.empty()) {
+    if (keyword == "Coords") {
       xflag = 1;
       if (flag) coords(line);
-      else skip_lines(natoms,line);
-    } else if (strcmp(keyword,"Types") == 0) {
+      else skip_lines(natoms,line,keyword);
+    } else if (keyword == "Types") {
       typeflag = 1;
       if (flag) types(line);
-      else skip_lines(natoms,line);
-    } else if (strcmp(keyword,"Molecules") == 0) {
+      else skip_lines(natoms,line,keyword);
+    } else if (keyword == "Molecules") {
       moleculeflag = 1;
       if (flag) molecules(line);
-      else skip_lines(natoms,line);
-    } else if (strcmp(keyword,"Fragments") == 0) {
+      else skip_lines(natoms,line,keyword);
+    } else if (keyword == "Fragments") {
       if (nfragments == 0)
-        error->all(FLERR,"Molecule file has fragments but no nfragments setting");
+        error->all(FLERR,"Found Fragments section but no nfragments setting in header");
       fragmentflag = 1;
       if (flag) fragments(line);
-      else skip_lines(nfragments,line);
-    } else if (strcmp(keyword,"Charges") == 0) {
+      else skip_lines(nfragments,line,keyword);
+    } else if (keyword == "Charges") {
       qflag = 1;
       if (flag) charges(line);
-      else skip_lines(natoms,line);
-    } else if (strcmp(keyword,"Diameters") == 0) {
+      else skip_lines(natoms,line,keyword);
+    } else if (keyword == "Diameters") {
       radiusflag = 1;
       if (flag) diameters(line);
-      else skip_lines(natoms,line);
-    } else if (strcmp(keyword,"Masses") == 0) {
+      else skip_lines(natoms,line,keyword);
+    } else if (keyword == "Masses") {
       rmassflag = 1;
       if (flag) masses(line);
-      else skip_lines(natoms,line);
+      else skip_lines(natoms,line,keyword);
 
-    } else if (strcmp(keyword,"Bonds") == 0) {
+    } else if (keyword == "Bonds") {
       if (nbonds == 0)
-        error->all(FLERR,"Molecule file has bonds but no nbonds setting");
+        error->all(FLERR,"Found Bonds section but no nbonds setting in header");
       bondflag = tag_require = 1;
       bonds(flag,line);
-    } else if (strcmp(keyword,"Angles") == 0) {
+    } else if (keyword == "Angles") {
       if (nangles == 0)
-        error->all(FLERR,"Molecule file has angles but no nangles setting");
+        error->all(FLERR,"Found Angles section but no nangles setting in header");
       angleflag = tag_require = 1;
       angles(flag,line);
-    } else if (strcmp(keyword,"Dihedrals") == 0) {
-      if (ndihedrals == 0) error->all(FLERR,"Molecule file has dihedrals "
-                                      "but no ndihedrals setting");
+    } else if (keyword == "Dihedrals") {
+      if (ndihedrals == 0) error->all(FLERR,"Found Dihedrals section"
+                                      "but no ndihedrals setting in header");
       dihedralflag = tag_require = 1;
       dihedrals(flag,line);
-    } else if (strcmp(keyword,"Impropers") == 0) {
-      if (nimpropers == 0) error->all(FLERR,"Molecule file has impropers "
-                                      "but no nimpropers setting");
+    } else if (keyword == "Impropers") {
+      if (nimpropers == 0) error->all(FLERR,"Found Impropers section"
+                                      "but no nimpropers setting in header");
       improperflag = tag_require = 1;
       impropers(flag,line);
 
-    } else if (strcmp(keyword,"Special Bond Counts") == 0) {
+    } else if (keyword == "Special Bond Counts") {
       nspecialflag = 1;
       nspecial_read(flag,line);
-    } else if (strcmp(keyword,"Special Bonds") == 0) {
+    } else if (keyword == "Special Bonds") {
       specialflag = tag_require = 1;
       if (flag) special_read(line);
-      else skip_lines(natoms,line);
+      else skip_lines(natoms,line,keyword);
 
-    } else if (strcmp(keyword,"Shake Flags") == 0) {
+    } else if (keyword == "Shake Flags") {
       shakeflagflag = 1;
       if (flag) shakeflag_read(line);
-      else skip_lines(natoms,line);
-    } else if (strcmp(keyword,"Shake Atoms") == 0) {
+      else skip_lines(natoms,line,keyword);
+    } else if (keyword == "Shake Atoms") {
       shakeatomflag = tag_require = 1;
       if (shaketypeflag) shakeflag = 1;
       if (!shakeflagflag)
-        error->all(FLERR,"Molecule file shake flags not before shake atoms");
+        error->all(FLERR,"Shake Flags section must come before Shake Atoms section");
       if (flag) shakeatom_read(line);
-      else skip_lines(natoms,line);
-    } else if (strcmp(keyword,"Shake Bond Types") == 0) {
+      else skip_lines(natoms,line,keyword);
+    } else if (keyword == "Shake Bond Types") {
       shaketypeflag = 1;
       if (shakeatomflag) shakeflag = 1;
       if (!shakeflagflag)
-        error->all(FLERR,"Molecule file shake flags not before shake bonds");
+        error->all(FLERR,"Shake Flags section must come before Shake Bonds section");
       if (flag) shaketype_read(line);
-      else skip_lines(natoms,line);
+      else skip_lines(natoms,line,keyword);
 
-    } else if (strcmp(keyword,"Body Integers") == 0) {
+    } else if (keyword == "Body Integers") {
       if (bodyflag == 0 || nibody == 0)
-        error->all(FLERR,"Molecule file has body params "
-                   "but no setting for them");
+        error->all(FLERR,"Found Body Integers section but no setting in header");
       ibodyflag = 1;
       body(flag,0,line);
-    } else if (strcmp(keyword,"Body Doubles") == 0) {
+    } else if (keyword == "Body Doubles") {
       if (bodyflag == 0 || ndbody == 0)
-        error->all(FLERR,"Molecule file has body params "
-                   "but no setting for them");
+        error->all(FLERR,"Found Body Doubles section but no setting in header");
       dbodyflag = 1;
       body(flag,1,line);
+    } else {
 
-    } else error->one(FLERR,fmt::format("Unknown section '{}' in molecule "
-                                        "file", keyword));
+      // Error: Either a too long/short section or a typo in the keyword
 
-    parse_keyword(1,line,keyword);
+      if (utils::strmatch(keyword,"^[A-Za-z ]+$"))
+        error->one(FLERR,fmt::format("Unknown section '{}' in molecule "
+                                     "file\n",keyword));
+      else error->one(FLERR,fmt::format("Unexpected line in molecule file "
+                                        "while looking for the next "
+                                        "section:\n{}",line));
+    }
+    keyword = parse_keyword(1,line);
   }
 
   // error check
@@ -678,11 +689,14 @@ void Molecule::coords(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
-      if (values.count() != 4) error->one(FLERR,"Invalid Coords section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 4)
+        error->all(FLERR,fmt::format("Invalid line in Coords section of "
+                                     "molecule file: {}",line));
 
       int iatom = values.next_int() - 1;
-      if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Coords section in molecule file");
+      if (iatom < 0 || iatom >= natoms)
+        error->all(FLERR,"Invalid atom index in Coords section of molecule file");
       count[iatom]++;
       x[iatom][0] = values.next_double();
       x[iatom][1] = values.next_double();
@@ -693,17 +707,19 @@ void Molecule::coords(char *line)
       x[iatom][2] *= sizescale;
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid Coords section in molecule file\n"
-                                  "{}", e.what()));
+    error->all(FLERR,fmt::format("Invalid line in Coords section of "
+                                 "molecule file: {}\n{}",e.what(),line));
   }
 
   for (int i = 0; i < natoms; i++)
-    if (count[i] == 0) error->all(FLERR,"Invalid Coords section in molecule file");
+    if (count[i] == 0) error->all(FLERR,fmt::format("Atom {} missing in Coords "
+                                                    "section of molecule file",i+1));
 
   if (domain->dimension == 2) {
     for (int i = 0; i < natoms; i++)
       if (x[i][2] != 0.0)
-        error->all(FLERR,"Molecule file z coord must be 0.0 for 2d");
+        error->all(FLERR,fmt::format("Z coord in molecule file for atom {} "
+                                     "must be 0.0 for 2d-simulation.",i+1));
   }
 }
 
@@ -719,29 +735,33 @@ void Molecule::types(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
-      if (values.count() != 2) error->one(FLERR,"Invalid Types section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 2)
+        error->all(FLERR,fmt::format("Invalid line in Types section of "
+                                     "molecule file: {}",line));
 
       int iatom = values.next_int() - 1;
-      if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Types section in molecule file");
+      if (iatom < 0 || iatom >= natoms)
+        error->all(FLERR,"Invalid atom index in Types section of molecule file");
       count[iatom]++;
       type[iatom] = values.next_int();
       type[iatom] += toffset;
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid Types section in molecule file\n"
-                                  "{}", e.what()));
+    error->all(FLERR, fmt::format("Invalid line in Types section of "
+                                  "molecule file: {}\n{}", e.what(),line));
   }
 
-  for (int i = 0; i < natoms; i++)
-    if (count[i] == 0) error->all(FLERR,"Invalid Types section in molecule file");
+  for (int i = 0; i < natoms; i++) {
+    if (count[i] == 0) error->all(FLERR,fmt::format("Atom {} missing in Types "
+                                                    "section of molecule file",i+1));
 
-  for (int i = 0; i < natoms; i++)
     if ((type[i] <= 0) || (domain->box_exist && (type[i] > atom->ntypes)))
-      error->all(FLERR,"Invalid atom type in molecule file");
+      error->all(FLERR,fmt::format("Invalid atom type {} for atom {} "
+                                   "in molecule file",type[i],i+1));
 
-  for (int i = 0; i < natoms; i++)
     ntypes = MAX(ntypes,type[i]);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -755,26 +775,31 @@ void Molecule::molecules(char *line)
   try {
     for (int i = 0; i < natoms; i++) {
       readline(line);
-      ValueTokenizer values(line);
-      if (values.count() != 2) error->one(FLERR,"Invalid Molecules section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 2)
+        error->all(FLERR,fmt::format("Invalid line in Molecules section of "
+                                     "molecule file: {}",line));
 
       int iatom = values.next_int() - 1;
-      if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Molecules section in molecule file");
+      if (iatom < 0 || iatom >= natoms)
+        error->all(FLERR,"Invalid atom index in Molecules section of molecule file");
       count[iatom]++;
       molecule[iatom] = values.next_tagint();
       // molecule[iatom] += moffset; // placeholder for possible molecule offset
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid Molecules section in molecule file\n"
-                                  "{}", e.what()));
+    error->all(FLERR, fmt::format("Invalid line in Molecules section of "
+                                  "molecule file: {}\n{}",e.what(),line));
   }
 
   for (int i = 0; i < natoms; i++)
-    if (count[i] == 0) error->all(FLERR,"Invalid Molecules section in molecule file");
+    if (count[i] == 0) error->all(FLERR,fmt::format("Atom {} missing in Molecules "
+                                                    "section of molecule file",i+1));
 
   for (int i = 0; i < natoms; i++)
-    if (molecule[i] <= 0)
-      error->all(FLERR,"Invalid molecule ID in molecule file");
+    if (molecule[i] < 0)
+      error->all(FLERR,fmt::format("Invalid molecule ID {} for atom {} "
+                                   "in molecule file",molecule[i],i+1));
 
   for (int i = 0; i < natoms; i++)
     nmolecules = MAX(nmolecules,molecule[i]);
@@ -790,23 +815,26 @@ void Molecule::fragments(char *line)
     for (int i = 0; i < nfragments; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
+      ValueTokenizer values(utils::trim_comment(line));
 
       if ((int)values.count() > natoms+1)
-        error->one(FLERR,"Invalid atom ID in Fragments section of molecule file");
+        error->all(FLERR,"Too many atoms per fragment in Fragments "
+                   "section of molecule file");
 
       fragmentnames[i] = values.next_string();
 
       while (values.has_next()) {
-        int atomID = values.next_int();
-        if (atomID <= 0 || atomID > natoms)
-          error->one(FLERR,"Invalid atom ID in Fragments section of molecule file");
-        fragmentmask[i][atomID-1] = 1;
+        int iatom = values.next_int()-1;
+        if (iatom < 0 || iatom >= natoms)
+          error->all(FLERR,fmt::format("Invalid atom ID {} for fragment {} in "
+                                       "Fragments section of molecule file",
+                                       iatom+1, fragmentnames[i]));
+        fragmentmask[i][iatom] = 1;
       }
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid atom ID in Fragments section of molecule file\n"
-                                  "{}", e.what()));
+    error->all(FLERR, fmt::format("Invalid atom ID in Fragments section of "
+                                  "molecule file: {}\n{}", e.what(),line));
   }
 }
 
@@ -821,21 +849,26 @@ void Molecule::charges(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
-      if ((int)values.count() != 2) error->one(FLERR,"Invalid Charges section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if ((int)values.count() != 2)
+        error->all(FLERR,fmt::format("Invalid line in Charges section of "
+                                     "molecule file: {}",line));
 
       int iatom = values.next_int() - 1;
-      if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Charges section in molecule file");
+      if (iatom < 0 || iatom >= natoms)
+        error->all(FLERR,"Invalid atom index in Charges section of molecule file");
+
       count[iatom]++;
       q[iatom] = values.next_double();
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid Charges section in molecule file\n"
-                                  "{}", e.what()));
+    error->all(FLERR, fmt::format("Invalid line in Charges section of "
+                                  "molecule file: {}.\n{}",e.what(),line));
   }
 
   for (int i = 0; i < natoms; i++)
-    if (count[i] == 0) error->all(FLERR,"Invalid Charges section in molecule file");
+    if (count[i] == 0)  error->all(FLERR,fmt::format("Atom {} missing in Charges "
+                                                    "section of molecule file",i+1));
 }
 
 /* ----------------------------------------------------------------------
@@ -850,11 +883,13 @@ void Molecule::diameters(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
-      if (values.count() != 2) error->one(FLERR,"Invalid Diameters section in molecule file");
-
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 2)
+        error->all(FLERR,fmt::format("Invalid line in Diameters section of "
+                                     "molecule file: {}",line));
       int iatom = values.next_int() - 1;
-      if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Diameters section in molecule file");
+      if (iatom < 0 || iatom >= natoms)
+        error->all(FLERR,"Invalid atom index in Diameters section of molecule file");
       count[iatom]++;
       radius[iatom] = values.next_double();
       radius[iatom] *= sizescale;
@@ -862,16 +897,17 @@ void Molecule::diameters(char *line)
       maxradius = MAX(maxradius,radius[iatom]);
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid Diameters section in molecule file\n"
-                                  "{}", e.what()));
+    error->all(FLERR, fmt::format("Invalid line in Diameters section of "
+                                  "molecule file: {}\n{}",e.what(),line));
   }
 
-  for (int i = 0; i < natoms; i++)
-    if (count[i] == 0) error->all(FLERR,"Invalid Diameters section in molecule file");
-
-  for (int i = 0; i < natoms; i++)
+  for (int i = 0; i < natoms; i++) {
+    if (count[i] == 0)  error->all(FLERR,fmt::format("Atom {} missing in Diameters "
+                                                    "section of molecule file",i+1));
     if (radius[i] < 0.0)
-      error->all(FLERR,"Invalid atom diameter in molecule file");
+      error->all(FLERR,fmt::format("Invalid atom diameter {} for atom {} "
+                                   "in molecule file", radius[i], i+1));
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -885,25 +921,30 @@ void Molecule::masses(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
-      if (values.count() != 2) error->one(FLERR,"Invalid Masses section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 2)
+        error->all(FLERR,fmt::format("Invalid line in Masses section of "
+                                     "molecule file: {}",line));
 
       int iatom = values.next_int() - 1;
-      if (iatom < 0 || iatom >= natoms) error->one(FLERR,"Invalid Masses section in molecule file");
+      if (iatom < 0 || iatom >= natoms)
+        error->all(FLERR,"Invalid atom index in Masses section of molecule file");
       count[iatom]++;
       rmass[iatom] = values.next_double();
       rmass[iatom] *= sizescale*sizescale*sizescale;
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid Masses section in molecule file\n"
-                                  "{}", e.what()));
+    error->all(FLERR, fmt::format("Invalid line in Masses section of "
+                                  "molecule file: {}\n{}",e.what(),line));
   }
 
-  for (int i = 0; i < natoms; i++)
-    if (count[i] == 0) error->all(FLERR,"Invalid Masses section in molecule file");
-
-  for (int i = 0; i < natoms; i++)
-    if (rmass[i] <= 0.0) error->all(FLERR,"Invalid atom mass in molecule file");
+  for (int i = 0; i < natoms; i++) {
+    if (count[i] == 0)   error->all(FLERR,fmt::format("Atom {} missing in Masses "
+                                                    "section of molecule file",i+1));
+    if (rmass[i] <= 0.0)
+      error->all(FLERR,fmt::format("Invalid atom mass {} for atom {} "
+                                   "in molecule file", radius[i], i+1));
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -929,24 +970,26 @@ void Molecule::bonds(int flag, char *line)
     readline(line);
 
     try {
-      ValueTokenizer values(line);
-      if (values.count() != 4) error->one(FLERR,"Invalid Bonds section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 4)
+        error->all(FLERR,fmt::format("Invalid line in Bonds section of "
+                                     "molecule file: {}",line));
       values.next_int();
       itype = values.next_int();
       atom1 = values.next_tagint();
       atom2 = values.next_tagint();
     } catch (TokenizerException &e) {
-      error->one(FLERR, fmt::format("Invalid Bonds section in molecule file\n"
-                                    "{}", e.what()));
+      error->all(FLERR, fmt::format("Invalid line in Bonds section of "
+                                  "molecule file: {}\n{}",e.what(),line));
     }
 
     itype += boffset;
 
     if ((atom1 <= 0) || (atom1 > natoms) ||
         (atom2 <= 0) || (atom2 > natoms) || (atom1 == atom2))
-      error->one(FLERR,"Invalid atom ID in Bonds section of molecule file");
+      error->all(FLERR,"Invalid atom ID in Bonds section of molecule file");
     if ((itype <= 0) || (domain->box_exist && (itype > atom->nbondtypes)))
-      error->one(FLERR,"Invalid bond type in Bonds section of molecule file");
+      error->all(FLERR,"Invalid bond type in Bonds section of molecule file");
 
     if (flag) {
       m = atom1-1;
@@ -997,16 +1040,18 @@ void Molecule::angles(int flag, char *line)
     readline(line);
 
     try {
-      ValueTokenizer values(line);
-      if (values.count() != 5) error->one(FLERR,"Invalid Angles section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 5)
+        error->all(FLERR,fmt::format("Invalid line in Angles section of "
+                                     "molecule file: {}",line));
       values.next_int();
       itype = values.next_int();
       atom1 = values.next_tagint();
       atom2 = values.next_tagint();
       atom3 = values.next_tagint();
     } catch (TokenizerException &e) {
-      error->one(FLERR, fmt::format("Invalid Angles section in molecule file\n"
-                                    "{}", e.what()));
+      error->all(FLERR, fmt::format("Invalid line in Angles section of "
+                                    "molecule file: {}\n{}",e.what(),line));
     }
 
     itype += aoffset;
@@ -1015,9 +1060,9 @@ void Molecule::angles(int flag, char *line)
         (atom2 <= 0) || (atom2 > natoms) ||
         (atom3 <= 0) || (atom3 > natoms) ||
         (atom1 == atom2) || (atom1 == atom3) || (atom2 == atom3))
-      error->one(FLERR,"Invalid atom ID in Angles section of molecule file");
+      error->all(FLERR,"Invalid atom ID in Angles section of molecule file");
     if ((itype <= 0) || (domain->box_exist && (itype > atom->nangletypes)))
-      error->one(FLERR,"Invalid angle type in Angles section of molecule file");
+      error->all(FLERR,"Invalid angle type in Angles section of molecule file");
 
     if (flag) {
       m = atom2-1;
@@ -1081,8 +1126,11 @@ void Molecule::dihedrals(int flag, char *line)
     readline(line);
 
     try {
-      ValueTokenizer values(line);
-      if (values.count() != 6) error->one(FLERR,"Invalid Dihedrals section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 6)
+        error->all(FLERR,fmt::format("Invalid line in Dihedrals section of "
+                                     "molecule file: {}",line));
+
       values.next_int();
       itype = values.next_int();
       atom1 = values.next_tagint();
@@ -1090,8 +1138,8 @@ void Molecule::dihedrals(int flag, char *line)
       atom3 = values.next_tagint();
       atom4 = values.next_tagint();
     } catch (TokenizerException &e) {
-      error->one(FLERR, fmt::format("Invalid Dihedrals section in molecule file\n"
-                                    "{}", e.what()));
+      error->all(FLERR, fmt::format("Invalid line in Dihedrals section of "
+                                    "molecule file: {}\n{}",e.what(),line));
     }
 
     itype += doffset;
@@ -1102,10 +1150,10 @@ void Molecule::dihedrals(int flag, char *line)
         (atom4 <= 0) || (atom4 > natoms) ||
         (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
         (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
-      error->one(FLERR,
+      error->all(FLERR,
                  "Invalid atom ID in dihedrals section of molecule file");
     if ((itype <= 0) || (domain->box_exist && (itype > atom->ndihedraltypes)))
-      error->one(FLERR,"Invalid dihedral type in Dihedrals section of molecule file");
+      error->all(FLERR,"Invalid dihedral type in Dihedrals section of molecule file");
 
     if (flag) {
       m = atom2-1;
@@ -1180,8 +1228,10 @@ void Molecule::impropers(int flag, char *line)
     readline(line);
 
     try {
-      ValueTokenizer values(line);
-      if (values.count() != 6) error->one(FLERR,"Invalid Impropers section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 6)
+                error->all(FLERR,fmt::format("Invalid line in Impropers section of "
+                                     "molecule file: {}",line));
       values.next_int();
       itype = values.next_int();
       atom1 = values.next_tagint();
@@ -1189,8 +1239,8 @@ void Molecule::impropers(int flag, char *line)
       atom3 = values.next_tagint();
       atom4 = values.next_tagint();
     } catch (TokenizerException &e) {
-      error->one(FLERR, fmt::format("Invalid Impropers section in molecule file\n"
-                                    "{}", e.what()));
+      error->all(FLERR, fmt::format("Invalid line in Impropers section of "
+                                    "molecule file: {}\n{}",e.what(),line));
     }
 
     itype += ioffset;
@@ -1201,10 +1251,10 @@ void Molecule::impropers(int flag, char *line)
         (atom4 <= 0) || (atom4 > natoms) ||
         (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
         (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
-      error->one(FLERR,
+      error->all(FLERR,
                  "Invalid atom ID in impropers section of molecule file");
     if ((itype <= 0) || (domain->box_exist && (itype > atom->nimpropertypes)))
-      error->one(FLERR,"Invalid improper type in Impropers section of molecule file");
+      error->all(FLERR,"Invalid improper type in Impropers section of molecule file");
 
     if (flag) {
       m = atom2-1;
@@ -1273,15 +1323,17 @@ void Molecule::nspecial_read(int flag, char *line)
     int c1, c2, c3;
 
     try {
-      ValueTokenizer values(line);
-      if (values.count() != 4) error->one(FLERR,"Invalid Special Bond Counts section in molecule file");
+      ValueTokenizer values(utils::trim_comment(line));
+      if (values.count() != 4)
+        error->all(FLERR,fmt::format("Invalid line in Special Bond Counts section of "
+                                     "molecule file: {}",line));
       values.next_int();
       c1 = values.next_tagint();
       c2 = values.next_tagint();
       c3 = values.next_tagint();
     } catch (TokenizerException &e) {
-      error->one(FLERR, fmt::format("Invalid Special Bond Counts section in molecule file\n"
-                                    "{}", e.what()));
+      error->all(FLERR, fmt::format("Invalid line in Special Bond Counts section of "
+                                    "molecule file: {}\n{}",e.what(),line));
     }
 
     if (flag) {
@@ -1302,11 +1354,11 @@ void Molecule::special_read(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
+      ValueTokenizer values(utils::trim_comment(line));
       int nwords = values.count();
 
       if (nwords != nspecial[i][2]+1)
-        error->one(FLERR,"Molecule file special list "
+        error->all(FLERR,"Molecule file special list "
                         "does not match special count");
 
       values.next_int(); // ignore
@@ -1315,12 +1367,13 @@ void Molecule::special_read(char *line)
         special[i][m-1] = values.next_tagint();
         if (special[i][m-1] <= 0 || special[i][m-1] > natoms ||
             special[i][m-1] == i+1)
-          error->one(FLERR,"Invalid special atom index in molecule file");
+          error->all(FLERR,"Invalid atom index in Special Bonds "
+                     "section of molecule file");
       }
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid Molecule file special list\n"
-                                  "{}", e.what()));
+    error->all(FLERR, fmt::format("Invalid line in Special Bonds section of "
+                                  "molecule file: {}\n{}",e.what(),line));
   }
 }
 
@@ -1351,7 +1404,7 @@ void Molecule::special_generate()
         nspecial[i][0]++;
         nspecial[atom2][0]++;
         if (count[i] >= atom->maxspecial || count[atom2] >= atom->maxspecial)
-          error->one(FLERR,"Molecule auto special bond generation overflow");
+          error->all(FLERR,"Molecule auto special bond generation overflow");
         tmpspecial[i][count[i]++] = atom2 + 1;
         tmpspecial[atom2][count[atom2]++] = i + 1;
       }
@@ -1363,7 +1416,7 @@ void Molecule::special_generate()
         atom1 = i;
         atom2 = bond_atom[i][j];
         if (count[atom1] >= atom->maxspecial)
-          error->one(FLERR,"Molecule auto special bond generation overflow");
+          error->all(FLERR,"Molecule auto special bond generation overflow");
         tmpspecial[i][count[atom1]++] = atom2;
       }
     }
@@ -1386,7 +1439,7 @@ void Molecule::special_generate()
         }
         if (!dedup) {
           if (count[i] >= atom->maxspecial)
-            error->one(FLERR,"Molecule auto special bond generation overflow");
+            error->all(FLERR,"Molecule auto special bond generation overflow");
           tmpspecial[i][count[i]++] = tmpspecial[tmpspecial[i][m]-1][j];
           nspecial[i][1]++;
         }
@@ -1410,7 +1463,7 @@ void Molecule::special_generate()
         }
         if (!dedup) {
           if (count[i] >= atom->maxspecial)
-            error->one(FLERR,"Molecule auto special bond generation overflow");
+            error->all(FLERR,"Molecule auto special bond generation overflow");
           tmpspecial[i][count[i]++] = tmpspecial[tmpspecial[i][m]-1][j];
           nspecial[i][2]++;
         }
@@ -1440,22 +1493,22 @@ void Molecule::shakeflag_read(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
+      ValueTokenizer values(utils::trim_comment(line));
 
       if (values.count() != 2)
-        error->one(FLERR,"Invalid Shake Flags section in molecule file");
+        error->all(FLERR,"Invalid Shake Flags section in molecule file");
 
       values.next_int();
       shake_flag[i] = values.next_int();
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid Shake Flags section in molecule file\n"
+    error->all(FLERR, fmt::format("Invalid Shake Flags section in molecule file\n"
                                   "{}", e.what()));
   }
 
   for (int i = 0; i < natoms; i++)
     if (shake_flag[i] < 0 || shake_flag[i] > 4)
-      error->one(FLERR,"Invalid shake flag in molecule file");
+      error->all(FLERR,"Invalid shake flag in molecule file");
 }
 
 /* ----------------------------------------------------------------------
@@ -1469,7 +1522,7 @@ void Molecule::shakeatom_read(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
+      ValueTokenizer values(utils::trim_comment(line));
       nmatch = values.count();
 
       switch (shake_flag[i]) {
@@ -1511,15 +1564,15 @@ void Molecule::shakeatom_read(char *line)
           break;
 
         default:
-          error->one(FLERR,"Invalid shake atom in molecule file");
+          error->all(FLERR,"Invalid shake atom in molecule file");
       }
 
       if (nmatch != nwant)
-        error->one(FLERR,"Invalid shake atom in molecule file");
+        error->all(FLERR,"Invalid shake atom in molecule file");
     }
 
   } catch (TokenizerException &e) {
-    error->one(FLERR,fmt::format("Invalid shake atom in molecule file\n"
+    error->all(FLERR,fmt::format("Invalid shake atom in molecule file\n"
                                  "{}", e.what()));
   }
 
@@ -1528,7 +1581,7 @@ void Molecule::shakeatom_read(char *line)
     if (m == 1) m = 3;
     for (int j = 0; j < m; j++)
       if (shake_atom[i][j] <= 0 || shake_atom[i][j] > natoms)
-        error->one(FLERR,"Invalid shake atom in molecule file");
+        error->all(FLERR,"Invalid shake atom in molecule file");
   }
 }
 
@@ -1543,7 +1596,7 @@ void Molecule::shaketype_read(char *line)
     for (int i = 0; i < natoms; i++) {
       readline(line);
 
-      ValueTokenizer values(line);
+      ValueTokenizer values(utils::trim_comment(line));
       nmatch = values.count();
 
       switch (shake_flag[i]) {
@@ -1582,14 +1635,14 @@ void Molecule::shaketype_read(char *line)
           break;
 
         default:
-          error->one(FLERR,"Invalid shake type data in molecule file");
+          error->all(FLERR,"Invalid shake type data in molecule file");
       }
 
       if (nmatch != nwant)
-        error->one(FLERR,"Invalid shake type data in molecule file");
+        error->all(FLERR,"Invalid shake type data in molecule file");
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid shake type data in molecule file\n",
+    error->all(FLERR, fmt::format("Invalid shake type data in molecule file\n",
                                   "{}", e.what()));
   }
 
@@ -1598,10 +1651,10 @@ void Molecule::shaketype_read(char *line)
     if (m == 1) m = 3;
     for (int j = 0; j < m-1; j++)
       if (shake_type[i][j] <= 0)
-        error->one(FLERR,"Invalid shake bond type in molecule file");
+        error->all(FLERR,"Invalid shake bond type in molecule file");
     if (shake_flag[i] == 1)
       if (shake_type[i][2] <= 0)
-        error->one(FLERR,"Invalid shake angle type in molecule file");
+        error->all(FLERR,"Invalid shake angle type in molecule file");
   }
 }
 
@@ -1621,13 +1674,13 @@ void Molecule::body(int flag, int pflag, char *line)
     while (nword < nparam) {
       readline(line);
 
-      ValueTokenizer values(line);
+      ValueTokenizer values(utils::trim_comment(line));
       int ncount = values.count();
 
       if (ncount == 0)
-        error->one(FLERR,"Too few values in body section of molecule file");
+        error->all(FLERR,"Too few values in body section of molecule file");
       if (nword+ncount > nparam)
-        error->one(FLERR,"Too many values in body section of molecule file");
+        error->all(FLERR,"Too many values in body section of molecule file");
 
       if (flag) {
         if (pflag == 0) {
@@ -1642,7 +1695,7 @@ void Molecule::body(int flag, int pflag, char *line)
       } else nword += ncount;
     }
   } catch (TokenizerException &e) {
-    error->one(FLERR, fmt::format("Invalid body params in molecule file\n",
+    error->all(FLERR, fmt::format("Invalid body params in molecule file\n",
                                   "{}", e.what()));
   }
 }
@@ -1939,18 +1992,6 @@ void Molecule::deallocate()
   memory->destroy(dbodyparams);
 }
 
-/* ----------------------------------------------------------------------
-   open molecule file
-------------------------------------------------------------------------- */
-
-void Molecule::open(char *file)
-{
-  fp = fopen(file,"r");
-  if (fp == nullptr)
-    error->one(FLERR,fmt::format("Cannot open molecule file {}: {}",
-                                 file, utils::getsyserror()));
-}
-
 /* ----------------------------------------------------------------------
    read and bcast a line
 ------------------------------------------------------------------------- */
@@ -1973,8 +2014,9 @@ void Molecule::readline(char *line)
    flag = 1, line has already been read
 ------------------------------------------------------------------------- */
 
-void Molecule::parse_keyword(int flag, char *line, char *keyword)
+std::string Molecule::parse_keyword(int flag, char *line)
 {
+  char line2[MAXLINE];
   if (flag) {
 
     // read upto non-blank line plus 1 following line
@@ -1986,42 +2028,38 @@ void Molecule::parse_keyword(int flag, char *line, char *keyword)
       while (eof == 0 && strspn(line," \t\n\r") == strlen(line)) {
         if (fgets(line,MAXLINE,fp) == nullptr) eof = 1;
       }
-      if (fgets(keyword,MAXLINE,fp) == nullptr) eof = 1;
+      if (fgets(line2,MAXLINE,fp) == nullptr) eof = 1;
     }
 
     // if eof, set keyword empty and return
 
     MPI_Bcast(&eof,1,MPI_INT,0,world);
     if (eof) {
-      keyword[0] = '\0';
-      return;
+      return std::string("");
     }
 
     // bcast keyword line to all procs
 
-    int n;
-    if (me == 0) n = strlen(line) + 1;
-    MPI_Bcast(&n,1,MPI_INT,0,world);
-    MPI_Bcast(line,n,MPI_CHAR,0,world);
+    MPI_Bcast(line,MAXLINE,MPI_CHAR,0,world);
   }
 
-  // copy non-whitespace and non-comment portion of line into keyword
+  // return non-whitespace and non-comment portion of line
 
-  int start = strspn(line," \t\n\r");
-  int stop = strcspn(line,"#") - 1;
-  while (line[stop] == ' ' || line[stop] == '\t'
-         || line[stop] == '\n' || line[stop] == '\r') stop--;
-  line[stop+1] = '\0';
-  strcpy(keyword,&line[start]);
+  return utils::trim(utils::trim_comment(line));
 }
 
 /* ----------------------------------------------------------------------
-   skip N lines of file
+   skip N lines of file. Check if non-numeric content (e.g. keyword).
 ------------------------------------------------------------------------- */
 
-void Molecule::skip_lines(int n, char *line)
+void Molecule::skip_lines(int n, char *line, const std::string &section)
 {
-  for (int i = 0; i < n; i++) readline(line);
+  for (int i = 0; i < n; i++) {
+    readline(line);
+    if (utils::strmatch(utils::trim(utils::trim_comment(line)),"^[A-Za-z ]+$"))
+      error->one(FLERR,fmt::format("Unexpected line in molecule file while "
+                                   "skipping {} section:\n{}",section,line));
+  }
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/molecule.h b/src/molecule.h
index 3316cd3f1f..00b2b123a0 100644
--- a/src/molecule.h
+++ b/src/molecule.h
@@ -161,10 +161,9 @@ class Molecule : protected Pointers {
   void allocate();
   void deallocate();
 
-  void open(char *);
   void readline(char *);
-  void parse_keyword(int, char *, char *);
-  void skip_lines(int, char *);
+  std::string parse_keyword(int, char *);
+  void skip_lines(int, char *, const std::string &);
 
   // void print();
 };
diff --git a/src/utils.cpp b/src/utils.cpp
index a762446150..4b9cc8a5f6 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -854,7 +854,7 @@ std::vector<std::string> utils::split_words(const std::string &text)
 ------------------------------------------------------------------------- */
 
 bool utils::is_integer(const std::string &str) {
-  if (str.size() == 0) {
+  if (str.empty()) {
     return false;
   }
 
@@ -870,7 +870,7 @@ bool utils::is_integer(const std::string &str) {
 ------------------------------------------------------------------------- */
 
 bool utils::is_double(const std::string &str) {
-  if (str.size() == 0) {
+  if (str.empty()) {
     return false;
   }
 
@@ -883,6 +883,22 @@ bool utils::is_double(const std::string &str) {
   return true;
 }
 
+/* ----------------------------------------------------------------------
+   Return whether string is a valid ID string
+------------------------------------------------------------------------- */
+
+bool utils::is_id(const std::string &str) {
+  if (str.empty()) {
+    return false;
+  }
+
+  for (auto c : str) {
+    if (isalnum(c) || (c == '_')) continue;
+    return false;
+  }
+  return true;
+}
+
 /* ----------------------------------------------------------------------
    strip off leading part of path, return just the filename
 ------------------------------------------------------------------------- */
diff --git a/src/utils.h b/src/utils.h
index 70a4c08cd3..eab81f1343 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -335,6 +335,14 @@ namespace LAMMPS_NS {
 
     bool is_double(const std::string &str);
 
+    /** Check if string is a valid ID
+     * ID strings may contain only letters, numbers, and underscores.
+     *
+     * \param str string that should be checked
+     * \return true, if string contains valid id, false otherwise */
+
+    bool is_id(const std::string &str);
+
     /** Try to detect pathname from FILE pointer.
      *
      * Currently only supported on Linux, otherwise will report "(unknown)".
diff --git a/unittest/formats/CMakeLists.txt b/unittest/formats/CMakeLists.txt
index 3c1ed66c06..7b96cf2d32 100644
--- a/unittest/formats/CMakeLists.txt
+++ b/unittest/formats/CMakeLists.txt
@@ -7,6 +7,10 @@ add_executable(test_image_flags test_image_flags.cpp)
 target_link_libraries(test_image_flags PRIVATE lammps GTest::GMock GTest::GTest)
 add_test(NAME ImageFlags COMMAND test_image_flags WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
 
+add_executable(test_molecule_file test_molecule_file.cpp)
+target_link_libraries(test_molecule_file PRIVATE lammps GTest::GMock GTest::GTest)
+add_test(NAME MoleculeFile COMMAND test_molecule_file WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+
 add_executable(test_pair_unit_convert test_pair_unit_convert.cpp)
 target_link_libraries(test_pair_unit_convert PRIVATE lammps GTest::GMock GTest::GTest)
 add_test(NAME PairUnitConvert COMMAND test_pair_unit_convert WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
diff --git a/unittest/formats/test_molecule_file.cpp b/unittest/formats/test_molecule_file.cpp
new file mode 100644
index 0000000000..3bfc0dc4fa
--- /dev/null
+++ b/unittest/formats/test_molecule_file.cpp
@@ -0,0 +1,313 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://lammps.sandia.gov/, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "atom.h"
+#include "info.h"
+#include "input.h"
+#include "lammps.h"
+#include "molecule.h"
+#include "utils.h"
+#include "gmock/gmock.h"
+#include "gtest/gtest.h"
+
+#include <cstdio>
+#include <mpi.h>
+#include <string>
+
+using namespace LAMMPS_NS;
+
+using testing::MatchesRegex;
+using testing::StrEq;
+
+using utils::split_words;
+
+#define test_name test_info_->name()
+
+#if defined(OMPI_MAJOR_VERSION)
+const bool have_openmpi = true;
+#else
+const bool have_openmpi = false;
+#endif
+
+#define TEST_FAILURE(errmsg, ...)                                 \
+    if (Info::has_exceptions()) {                                 \
+        ::testing::internal::CaptureStdout();                     \
+        ASSERT_ANY_THROW({__VA_ARGS__});                          \
+        auto mesg = ::testing::internal::GetCapturedStdout();     \
+        ASSERT_THAT(mesg, MatchesRegex(errmsg));                  \
+    } else {                                                      \
+        if (!have_openmpi) {                                      \
+            ::testing::internal::CaptureStdout();                 \
+            ASSERT_DEATH({__VA_ARGS__}, "");                      \
+            auto mesg = ::testing::internal::GetCapturedStdout(); \
+            ASSERT_THAT(mesg, MatchesRegex(errmsg));              \
+        }                                                         \
+    }
+
+static void create_molecule_files()
+{
+    // create molecule files
+    const char h2o_file[] = "# Water molecule. SPC/E model.\n\n3 atoms\n2 bonds\n1 angles\n\n"
+                            "Coords\n\n1 1.12456 0.09298 1.27452\n"
+                            "2 1.53683 0.75606 1.89928\n3 0.49482 0.56390 0.65678\n\n"
+                            "Types\n\n1 1\n2 2\n3 2\n\n"
+                            "Charges\n\n1 -0.8472\n2 0.4236\n3 0.4236\n\n"
+                            "Bonds\n\n1 1 1 2\n2 1 1 3\n\n"
+                            "Angles\n\n1 1 2 1 3\n\n"
+                            "Shake Flags\n\n1 1\n2 1\n3 1\n\n"
+                            "Shake Atoms\n\n1 1 2 3\n2 1 2 3\n3 1 2 3\n\n"
+                            "Shake Bond Types\n\n1 1 1 1\n2 1 1 1\n3 1 1 1\n\n"
+                            "Special Bond Counts\n\n1 2 0 0\n2 1 1 0\n3 1 1 0\n\n"
+                            "Special Bonds\n\n1 2 3\n2 1 3\n3 1 2\n\n";
+    const char co2_file[] = "# CO2 molecule file. TraPPE model.\n\n"
+                            "3 atoms\n2 bonds\n1 angles\n\n"
+                            "Coords\n\n1 0.0 0.0 0.0\n2 -1.16 0.0 0.0\n3 1.16 0.0 0.0\n\n"
+                            "Types\n\n1 1\n2 2\n3 2\n\n"
+                            "Charges\n\n1 0.7\n2 -0.35\n3 -0.35\n\n"
+                            "Bonds\n\n1 1 1 2\n2 1 1 3\n\n"
+                            "Angles\n\n1 1 2 1 3\n\n"
+                            "Special Bond Counts\n\n1 2 0 0\n2 1 1 0\n3 1 1 0\n\n"
+                            "Special Bonds\n\n1 2 3\n2 1 3\n3 1 2\n\n";
+
+    FILE *fp = fopen("tmp.h2o.mol", "w");
+    if (fp) {
+        fputs(h2o_file, fp);
+        fclose(fp);
+    }
+    rename("tmp.h2o.mol", "h2o.mol");
+    fp = fopen("tmp.co2.mol", "w");
+    if (fp) {
+        fputs(co2_file, fp);
+        fclose(fp);
+    }
+    rename("tmp.co2.mol", "co2.mol");
+}
+
+// whether to print verbose output (i.e. not capturing LAMMPS screen output).
+bool verbose = false;
+
+class MoleculeFileTest : public ::testing::Test {
+protected:
+    LAMMPS *lmp;
+
+    void SetUp() override
+    {
+        const char *args[] = {"MoleculeFileTest", "-log", "none", "-echo", "screen", "-nocite"};
+        char **argv        = (char **)args;
+        int argc           = sizeof(args) / sizeof(char *);
+        if (!verbose) ::testing::internal::CaptureStdout();
+        lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
+        create_molecule_files();
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        ASSERT_NE(lmp, nullptr);
+    }
+
+    void TearDown() override
+    {
+        if (!verbose) ::testing::internal::CaptureStdout();
+        delete lmp;
+        if (!verbose) ::testing::internal::GetCapturedStdout();
+        remove("h2o.mol");
+        remove("co2.mol");
+    }
+
+    void run_mol_cmd(const std::string &name, const std::string &args, const std::string &content)
+    {
+        std::string file = name + ".mol";
+        FILE *fp         = fopen(file.c_str(), "w");
+        fputs(content.c_str(), fp);
+        fclose(fp);
+
+        lmp->input->one(fmt::format("molecule {} {} {}", name, file, args));
+        remove(file.c_str());
+    }
+};
+
+TEST_F(MoleculeFileTest, nofile)
+{
+    TEST_FAILURE(".*Cannot open molecule file nofile.mol.*",
+                 lmp->input->one("molecule 1 nofile.mol"););
+}
+
+TEST_F(MoleculeFileTest, badid)
+{
+    TEST_FAILURE(".*Molecule template ID must have only "
+                 "alphanumeric or underscore characters.*",
+                 lmp->input->one("molecule @mol nofile.mol"););
+}
+
+TEST_F(MoleculeFileTest, badargs)
+{
+    TEST_FAILURE(".*Illegal molecule command.*",
+                 run_mol_cmd(test_name, "offset 1 2 3 4",
+                             "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"););
+    TEST_FAILURE(
+        ".*Illegal molecule command.*",
+        run_mol_cmd(test_name, "toff", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"););
+    TEST_FAILURE(
+        ".*Illegal molecule command.*",
+        run_mol_cmd(test_name, "boff", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"););
+    TEST_FAILURE(
+        ".*Illegal molecule command.*",
+        run_mol_cmd(test_name, "aoff", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"););
+    TEST_FAILURE(
+        ".*Illegal molecule command.*",
+        run_mol_cmd(test_name, "doff", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"););
+    TEST_FAILURE(
+        ".*Illegal molecule command.*",
+        run_mol_cmd(test_name, "ioff", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"););
+    TEST_FAILURE(
+        ".*Illegal molecule command.*",
+        run_mol_cmd(test_name, "scale", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n"););
+    remove("badargs.mol");
+}
+
+TEST_F(MoleculeFileTest, noatom)
+{
+    TEST_FAILURE(".*ERROR: No atoms or invalid atom count in molecule file.*",
+                 run_mol_cmd(test_name, "",
+                             "Comment\n0 atoms\n1 bonds\n\n"
+                             " Coords\n\nBonds\n\n 1 1 2\n"););
+    remove("noatom.mol");
+}
+
+TEST_F(MoleculeFileTest, empty)
+{
+    TEST_FAILURE(".*ERROR: Unexpected end of molecule file.*",
+                 run_mol_cmd(test_name, "", "Comment\n\n"););
+    remove("empty.mol");
+}
+
+TEST_F(MoleculeFileTest, nospecial)
+{
+    TEST_FAILURE(".*ERROR: Cannot auto-generate special bonds before simulation box is defined.*",
+                 run_mol_cmd(test_name, "",
+                             "Comment\n3 atoms\n\n2 bonds\n\n"
+                             " Coords\n\n 1 1.0 1.0 1.0\n 2 1.0 1.0 0.0\n 3 1.0 0.0 1.0\n"
+                             " Bonds\n\n 1 1 1 2\n 2 1 1 3\n"););
+    remove("nospecial.mol");
+}
+
+TEST_F(MoleculeFileTest, minimal)
+{
+    ::testing::internal::CaptureStdout();
+    run_mol_cmd(test_name, "", "Comment\n1 atoms\n\n Coords\n\n 1 0.0 0.0 0.0\n");
+    auto output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    ASSERT_THAT(output, MatchesRegex(".*Read molecule template.*1 molecules.*1 atoms.*0 bonds.*"));
+}
+
+TEST_F(MoleculeFileTest, twomols)
+{
+    ::testing::internal::CaptureStdout();
+    run_mol_cmd(test_name, "",
+                "Comment\n2 atoms\n\n"
+                " Coords\n\n 1 0.0 0.0 0.0\n 2 0.0 0.0 1.0\n"
+                " Molecules\n\n 1 1\n 2 2\n\n Types\n\n 1 1\n 2 2\n\n");
+    auto output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    ASSERT_THAT(output, MatchesRegex(".*Read molecule template.*2 molecules.*2 atoms "
+                                     "with max type 2.*0 bonds.*"));
+}
+
+TEST_F(MoleculeFileTest, twofiles)
+{
+    ::testing::internal::CaptureStdout();
+    lmp->input->one("molecule twomols h2o.mol co2.mol offset 2 1 1 0 0");
+    auto output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    ASSERT_THAT(output, MatchesRegex(".*Read molecule template twomols:.*1 molecules.*3 atoms "
+                                     "with max type 2.*2 bonds with max type 1.*"
+                                     "1 angles with max type 1.*0 dihedrals.*"
+                                     ".*Read molecule template twomols:.*1 molecules.*3 atoms "
+                                     "with max type 4.*2 bonds with max type 2.*"
+                                     "1 angles with max type 2.*0 dihedrals.*"));
+}
+
+TEST_F(MoleculeFileTest, bonds)
+{
+    ::testing::internal::CaptureStdout();
+    lmp->input->one("atom_style bond");
+    lmp->input->one("region box block 0 1 0 1 0 1");
+    lmp->input->one("create_box 2 box bond/types 2 extra/bond/per/atom 2 "
+                    "extra/special/per/atom 4");
+    run_mol_cmd(test_name, "",
+                "Comment\n"
+                "4 atoms\n"
+                "2 bonds\n\n"
+                " Coords\n\n"
+                " 1 1.0 1.0 1.0\n"
+                " 2 1.0 1.0 0.0\n"
+                " 3 1.0 0.0 1.0\n"
+                " 4 1.0 0.0 0.0\n"
+                " Types\n\n"
+                " 1 1\n"
+                " 2 1\n"
+                " 3 2\n"
+                " 4 2\n\n"
+                " Bonds\n\n"
+                " 1 1 1 2\n"
+                " 2 2 1 3\n\n");
+    auto output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    ASSERT_THAT(output, MatchesRegex(".*Read molecule template.*1 molecules.*4 atoms.*type.*2.*"
+                                     "2 bonds.*type.*2.*0 angles.*"));
+
+    ::testing::internal::CaptureStdout();
+    lmp->input->one("mass * 2.0");
+    lmp->input->one("create_atoms 0 single 0.5 0.5 0.5 mol bonds 67235");
+    output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+    ASSERT_THAT(output, MatchesRegex(".*Created 4 atoms.*"));
+
+    ::testing::internal::CaptureStdout();
+    Molecule *mol = lmp->atom->molecules[0];
+    ASSERT_EQ(mol->natoms, 4);
+    ASSERT_EQ(lmp->atom->natoms, 4);
+    mol->compute_mass();
+    mol->compute_com();
+    ASSERT_DOUBLE_EQ(mol->masstotal, 8.0);
+    EXPECT_DOUBLE_EQ(mol->com[0], 1.0);
+    EXPECT_DOUBLE_EQ(mol->com[1], 0.5);
+    EXPECT_DOUBLE_EQ(mol->com[2], 0.5);
+    EXPECT_DOUBLE_EQ(mol->maxextent, sqrt(2));
+    EXPECT_EQ(mol->comatom, 1);
+    output = ::testing::internal::GetCapturedStdout();
+    if (verbose) std::cout << output;
+}
+
+int main(int argc, char **argv)
+{
+    MPI_Init(&argc, &argv);
+    ::testing::InitGoogleMock(&argc, argv);
+
+    if (have_openmpi && !LAMMPS_NS::Info::has_exceptions())
+        std::cout << "Warning: using OpenMPI without exceptions. "
+                     "Death tests will be skipped\n";
+
+    // handle arguments passed via environment variable
+    if (const char *var = getenv("TEST_ARGS")) {
+        std::vector<std::string> env = split_words(var);
+        for (auto arg : env) {
+            if (arg == "-v") {
+                verbose = true;
+            }
+        }
+    }
+
+    if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true;
+
+    int rv = RUN_ALL_TESTS();
+    MPI_Finalize();
+    return rv;
+}
diff --git a/unittest/utils/test_utils.cpp b/unittest/utils/test_utils.cpp
index a0e4022d2c..5b509318c9 100644
--- a/unittest/utils/test_utils.cpp
+++ b/unittest/utils/test_utils.cpp
@@ -287,6 +287,61 @@ TEST(Utils, signed_double_and_d_exponential)
     ASSERT_FALSE(utils::is_double("-10D-22"));
 }
 
+TEST(Utils, valid_id1)
+{
+    ASSERT_TRUE(utils::is_id("abc"));
+}
+
+TEST(Utils, valid_id2)
+{
+    ASSERT_TRUE(utils::is_id("123"));
+}
+
+TEST(Utils, valid_id3)
+{
+    ASSERT_TRUE(utils::is_id("abc123"));
+}
+
+TEST(Utils, valid_id4)
+{
+    ASSERT_TRUE(utils::is_id("abc_123"));
+}
+
+TEST(Utils, valid_id5)
+{
+    ASSERT_TRUE(utils::is_id("123_abc"));
+}
+
+TEST(Utils, valid_id6)
+{
+    ASSERT_TRUE(utils::is_id("_123"));
+}
+
+TEST(Utils, valid_id7)
+{
+    ASSERT_TRUE(utils::is_id("___"));
+}
+
+TEST(Utils, invalid_id1)
+{
+    ASSERT_FALSE(utils::is_id("+abc"));
+}
+
+TEST(Utils, invalid_id2)
+{
+    ASSERT_FALSE(utils::is_id("a[1]"));
+}
+
+TEST(Utils, invalid_id3)
+{
+    ASSERT_FALSE(utils::is_id("b(c)"));
+}
+
+TEST(Utils, invalid_id4)
+{
+    ASSERT_FALSE(utils::is_id("a$12"));
+}
+
 TEST(Utils, strmatch_beg)
 {
     ASSERT_TRUE(utils::strmatch("rigid/small/omp", "^rigid"));