prepare for multiple data files

not yet tested

src/atom.cpp

@@ -313,6 +313,7 @@ Atom::~Atom()
   memory->sfree(molecules);
 
   // delete label map
+
   delete lmap;
 
   // delete per-type arrays

src/label_map.cpp

@@ -15,6 +15,7 @@
 
 #include "force.h"
 #include "memory.h"
+#include "error.h"
 
 #include <cstring>
 
@@ -39,6 +40,12 @@ LabelMap::~LabelMap()
   atypelabel.clear();
   dtypelabel.clear();
   itypelabel.clear();
+
+  delete [] lmap2lmap.atom;
+  delete [] lmap2lmap.bond;
+  delete [] lmap2lmap.angle;
+  delete [] lmap2lmap.dihedral;
+  delete [] lmap2lmap.improper;
 }
 
 /* ----------------------------------------------------------------------
@@ -76,54 +83,127 @@ void LabelMap::allocate_type_labels()
     for (int i = 0; i < nimpropertypes; i++)
       itypelabel[i] = fmt::format("{}",i);
   }
+
+  // allocate space for lmap2lmap
+
+  lmap2lmap.atom = new int[natomtypes];
+  lmap2lmap.bond = new int[nbondtypes];
+  lmap2lmap.angle = new int[nangletypes];
+  lmap2lmap.dihedral = new int[ndihedraltypes];
+  lmap2lmap.improper = new int[nimpropertypes];
 }
 
 /* ----------------------------------------------------------------------
-   copy lmap1 to lmap2
+   copy another map (lmap2) into this one
+   if label already exists, leave in place
+   else, put new label in next available slot
 ------------------------------------------------------------------------- */
 
-void LabelMap::copy_lmap(LabelMap *lmap1, LabelMap *lmap2)
+void LabelMap::merge_lmap(class LabelMap *lmap2)
 {
-  int ncopy;
+  for (int i = 0; i < lmap2->natomtypes; i++)
-
+    find_or_create(lmap2->typelabel[i],typelabel,natomtypes);
-  ncopy = min(lmap1->natomtypes, lmap2->natomtypes);
-  for (int i = 0; i < ncopy; i++)
-    lmap2->typelabel[i] = lmap1->typelabel[i];
 
   if (force->bond) {
-    ncopy = min(lmap1->nbondtypes, lmap2->nbondtypes);
+    for (int i = 0; i < lmap2->nbondtypes; i++)
-    for (int i = 0; i < ncopy; i++)
+      find_or_create(lmap2->btypelabel[i],btypelabel,nbondtypes);
-      lmap1->btypelabel[i] = lmap2->btypelabel[i];
   }
 
   if (force->angle) {
-    ncopy = min(lmap1->nangletypes, lmap2->nangletypes);
+    for (int i = 0; i < lmap2->nangletypes; i++)
-    for (int i = 0; i < ncopy; i++)
+      find_or_create(lmap2->atypelabel[i],atypelabel,nangletypes);
-      lmap1->atypelabel[i] = lmap2->atypelabel[i];
   }
 
   if (force->dihedral) {
-    ncopy = min(lmap1->ndihedraltypes, lmap2->ndihedraltypes);
+    for (int i = 0; i < lmap2->ndihedraltypes; i++)
-    for (int i = 0; i < ncopy; i++)
+      find_or_create(lmap2->dtypelabel[i],dtypelabel,ndihedraltypes);
-      lmap1->dtypelabel[i] = lmap2->dtypelabel[i];
   }
 
   if (force->improper) {
-    ncopy = min(lmap1->nimpropertypes, lmap2->nimpropertypes);
+    for (int i = 0; i < lmap2->nimpropertypes; i++)
-    for (int i = 0; i < ncopy; i++)
+      find_or_create(lmap2->itypelabel[i],itypelabel,nimpropertypes);
-      lmap1->itypelabel[i] = lmap2->itypelabel[i];
   }
 }
 
+/* ----------------------------------------------------------------------
+   get mapping between this label map and another (lmap2)
+   values of lmap2lmap point to equivalent indices in lmap2
+------------------------------------------------------------------------- */
+
+void LabelMap::create_lmap2lmap(class LabelMap *lmap2)
+{
+  int type;
+
+  for (int i = 0; i < natomtypes; i++) {
+    type = find(typelabel[i],lmap2->typelabel,lmap2->natomtypes);
+    lmap2lmap.atom[i] = type - 1;
+  }
+
+  if (force->bond) {
+    lmap2lmap.bond = new int[nbondtypes];
+    for (int i = 0; i < nbondtypes; i++) {
+      type = find(btypelabel[i],lmap2->btypelabel,lmap2->nbondtypes);
+      lmap2lmap.bond[i] = type - 1;
+    }
+  }
+
+  if (force->angle) {
+    for (int i = 0; i < nangletypes; i++) {
+      type = find(atypelabel[i],lmap2->atypelabel,lmap2->nangletypes);
+      lmap2lmap.angle[i] = type - 1;
+    }
+  }
+
+  if (force->dihedral) {
+    for (int i = 0; i < ndihedraltypes; i++) {
+      type = find(dtypelabel[i],lmap2->dtypelabel,lmap2->ndihedraltypes);
+      lmap2lmap.dihedral[i] = type - 1;
+    }
+  }
+
+  if (force->improper) {
+    for (int i = 0; i < nimpropertypes; i++) {
+      type = find(itypelabel[i],lmap2->itypelabel,lmap2->nimpropertypes);
+      lmap2lmap.improper[i] = type - 1;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   find type label with name or create type if it doesn't exist
+   return numeric type
+------------------------------------------------------------------------- */
+
+int LabelMap::find_or_create(std::string mylabel, std::vector<std::string> labels, int ntypes)
+{
+  for (int i = 0; i < ntypes; i++)
+    if (labels[i] == mylabel) return i+1;
+
+  // if no match found, create new label at next available index
+  // label map assumed to be intialized with numeric index
+  // user labels are assumed to be alphanumeric (not a number)
+
+  for (int i = 0; i < ntypes; i++) {
+    if (stoi(labels[i]) != i+1) {
+      labels[i] = mylabel;
+      return i+1;
+    }
+  }
+
+  // if label cannot be found or created, need more space reserved
+
+  error->all(FLERR,"Topology type exceeds system topology type");
+}
+
 /* ----------------------------------------------------------------------
    find integer type given a type label
    return -1 if type not yet defined
 ------------------------------------------------------------------------- */
 
-int LabelMap::find_type(char *mytypelabel, char **typelabelarray, int num_types)
+int LabelMap::find(std::string mylabel, std::vector<std::string> labels, int ntypes)
 {
-  for (int i = 0; i < num_types; i++) {
+  for (int i = 0; i < ntypes; i++)
-    if (typelabelarray[i] && strcmp(mytypelabel,typelabelarray[i]) == 0) return i+1;
+    if (labels[i] == mylabel) return i+1;
-  }
+
   return -1;
 }

src/label_map.h

@@ -25,12 +25,26 @@ class LabelMap : protected Pointers {
    std::vector<std::string> typelabel,btypelabel,atypelabel;
    std::vector<std::string> dtypelabel,itypelabel;
 
+   // per-type data struct mapping this label map to another
+
+   struct Lmap2Lmap {
+     int *atom;
+     int *bond;
+     int *angle;
+     int *dihedral;
+     int *improper;
+   };
+
+   Lmap2Lmap lmap2lmap;
+
    LabelMap(LAMMPS *lmp);
    ~LabelMap();
 
    void allocate_type_labels();
-   void copy_lmap(class LabelMap *, class LabelMap *);
+   void merge_lmap(class LabelMap *); // copy another lmap into this one
-   int find_type(char *, char **, int);
+   void create_lmap2lmap(class LabelMap *); // index mapping between two lmaps
+   int find_or_create(std::string, std::vector<std::string>, int); // look up type or create new type
+   int find(std::string, std::vector<std::string>, int); // look up type index
 
  protected:
 };
@@ -41,5 +55,9 @@ class LabelMap : protected Pointers {
 
 /* ERROR/WARNING messages:
 
+E: Topology type exceeds system topology type
+
+The number of bond, angle, etc types exceeds the system setting. See
+the create_box or read_data command for how to specify these values.
 
 */

src/read_data.cpp

@@ -1981,13 +1981,13 @@ void ReadData::typelabels(std::vector<std::string> &mytypelabel, int myntypes)
   }
   delete [] typelabel;
 
-  // if first data file, assign this label map to atom class
+  // if first data file, assign this read_data label map to atom class
   // else, determine mapping to let labels override numeric types
 
   if (addflag == NONE) {
-    lmap->copy_lmap(lmap,atom->lmap);
+    atom->lmap->merge_lmap(lmap);
   } else {
-    ; // get lmap2lmap mapping. ...in progress...
+    lmap->create_lmap2lmap(atom->lmap);
   }
 }