Merge branch 'merge_oct23_2015' into lammps-icms

# Resolved Conflicts:
#	src/molecule.h
#	src/respa.cpp

src/USER-INTEL/fix_intel.cpp

@@ -664,7 +664,7 @@ int FixIntel::set_host_affinity(const int nomp)
   FILE *p;
   char cmd[512];
   char readbuf[INTEL_MAX_HOST_CORE_COUNT*5];
-  sprintf(cmd, "lscpu -p=cpu,core,socket | grep -v '#' |"
+  sprintf(cmd, "lscpu -p | grep -v '#' |"
 	  "sort -t, -k 3,3n -k 2,2n | awk -F, '{print $1}'");
   p = popen(cmd, "r");
   if (p == NULL) return -1;

src/atom.cpp

@@ -950,7 +950,8 @@ void Atom::data_vels(int n, char *buf, tagint id_offset)
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
-void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset)
+void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset,
+                      int type_offset) 
 {
   int m,tmp,itype;
   tagint atom1,atom2;
@@ -966,6 +967,7 @@ void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset)
       atom1 += id_offset;
       atom2 += id_offset;
     }
+    itype += type_offset;
 
     if (atom1 <= 0 || atom1 > map_tag_max ||
         atom2 <= 0 || atom2 > map_tag_max)
@@ -1001,7 +1003,8 @@ void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset)
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
-void Atom::data_angles(int n, char *buf, int *count, tagint id_offset)
+void Atom::data_angles(int n, char *buf, int *count, tagint id_offset,
+                       int type_offset) 
 {
   int m,tmp,itype;
   tagint atom1,atom2,atom3;
@@ -1018,6 +1021,7 @@ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset)
       atom2 += id_offset;
       atom3 += id_offset;
     }
+    itype += type_offset;
 
     if (atom1 <= 0 || atom1 > map_tag_max ||
         atom2 <= 0 || atom2 > map_tag_max ||
@@ -1068,7 +1072,8 @@ void Atom::data_angles(int n, char *buf, int *count, tagint id_offset)
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
-void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset)
+void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset,
+                          int type_offset)
 {
   int m,tmp,itype;
   tagint atom1,atom2,atom3,atom4;
@@ -1087,6 +1092,7 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset)
       atom3 += id_offset;
       atom4 += id_offset;
     }
+    itype += type_offset;
 
     if (atom1 <= 0 || atom1 > map_tag_max ||
         atom2 <= 0 || atom2 > map_tag_max ||
@@ -1153,7 +1159,8 @@ void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset)
    check that atom IDs are > 0 and <= map_tag_max
 ------------------------------------------------------------------------- */
 
-void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset)
+void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset,
+                          int type_offset)
 {
   int m,tmp,itype;
   tagint atom1,atom2,atom3,atom4;
@@ -1172,6 +1179,7 @@ void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset)
       atom3 += id_offset;
       atom4 += id_offset;
     }
+    itype += type_offset;
 
     if (atom1 <= 0 || atom1 > map_tag_max ||
         atom2 <= 0 || atom2 > map_tag_max ||
@@ -1496,17 +1504,17 @@ void Atom::add_molecule(int narg, char **arg)
   if (find_molecule(arg[0]) >= 0) 
     error->all(FLERR,"Reuse of molecule template ID");
 
-  // may over-allocate if not all args are mol files, but OK for srealloc
-
-  molecules = (Molecule **)
-    memory->srealloc(molecules,(nmolecule+narg-1)*sizeof(Molecule *),
-                     "atom::molecules");
-
   // 1st molecule in set stores nset = # of mols, others store nset = 0
+  // ifile = count of molecules in set
+  // index = argument index where next molecule starts, updated by constructor
 
   int ifile = 1;
+  int index = 1;
   while (1) {
-    molecules[nmolecule] = new Molecule(lmp,narg,arg,ifile);
+    molecules = (Molecule **)
+      memory->srealloc(molecules,(nmolecule+1)*sizeof(Molecule *),
+                       "atom::molecules");
+    molecules[nmolecule] = new Molecule(lmp,narg,arg,index);
     molecules[nmolecule]->nset = 0;
     molecules[nmolecule-ifile+1]->nset++;
     nmolecule++;

src/atom.h

@@ -213,10 +213,10 @@ class Atom : protected Pointers {
   void data_atoms(int, char *, tagint, int, int, double *);
   void data_vels(int, char *, tagint);
 
-  void data_bonds(int, char *, int *, tagint);
+  void data_bonds(int, char *, int *, tagint, int);
-  void data_angles(int, char *, int *, tagint);
+  void data_angles(int, char *, int *, tagint, int);
-  void data_dihedrals(int, char *, int *, tagint);
+  void data_dihedrals(int, char *, int *, tagint, int);
-  void data_impropers(int, char *, int *, tagint);
+  void data_impropers(int, char *, int *, tagint, int);
 
   void data_bonus(int, char *, class AtomVec *, tagint);
   void data_bodies(int, char *, class AtomVecBody *, tagint);

src/molecule.cpp

@@ -35,11 +35,12 @@ using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
-Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int ifile) : Pointers(lmp)
+Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) : 
+  Pointers(lmp)
 {
   me = comm->me;
 
-  if (ifile >= narg) error->all(FLERR,"Illegal molecule command");
+  if (index >= narg) error->all(FLERR,"Illegal molecule command");
 
   int n = strlen(arg[0]) + 1;
   id = new char[n];
@@ -50,21 +51,14 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int ifile) : Pointers(lmp)
       error->all(FLERR,"Molecule template ID must be "
                  "alphanumeric or underscore characters");
 
-  // scan args past ifile to reach optional args
+  // parse args until reach unknown arg (next file)
-  // set last = 1 if no more files in list
-
-  last = 0;
-  int iarg = ifile+1;
-  while (iarg < narg) {
-    if (strcmp(arg[iarg],"offset") == 0) break;
-    iarg++;
-  }
-  if (iarg == ifile+1) last = 1;
-
-  // parse optional args
 
   toffset = 0;
   boffset = aoffset = doffset = ioffset = 0;
+  sizescale = 1.0;
+
+  int ifile = index;
+  int iarg = ifile+1;
 
   while (iarg < narg) {
     if (strcmp(arg[iarg],"offset") == 0) {
@@ -78,9 +72,45 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int ifile) : Pointers(lmp)
           doffset < 0 || ioffset < 0) 
         error->all(FLERR,"Illegal molecule command");
       iarg += 6;
-    } else error->all(FLERR,"Illegal molecule command");
+    } else if (strcmp(arg[iarg],"toff") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+      toffset = force->inumeric(FLERR,arg[iarg+1]);
+      if (toffset < 0) error->all(FLERR,"Illegal molecule command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"boff") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+      boffset = force->inumeric(FLERR,arg[iarg+1]);
+      if (boffset < 0) error->all(FLERR,"Illegal molecule command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"aoff") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+      aoffset = force->inumeric(FLERR,arg[iarg+1]);
+      if (aoffset < 0) error->all(FLERR,"Illegal molecule command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"doff") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+      doffset = force->inumeric(FLERR,arg[iarg+1]);
+      if (doffset < 0) error->all(FLERR,"Illegal molecule command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"ioff") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+      ioffset = force->inumeric(FLERR,arg[iarg+1]);
+      if (ioffset < 0) error->all(FLERR,"Illegal molecule command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"scale") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal molecule command");
+      sizescale = force->numeric(FLERR,arg[iarg+1]);
+      if (sizescale <= 0.0) error->all(FLERR,"Illegal molecule command");
+      iarg += 2;
+    } else break;
   }
 
+  index = iarg;
+
+  // last molecule if have scanned all args
+
+  if (iarg == narg) last = 1;
+
   // initialize all fields to empty
 
   initialize();
@@ -392,10 +422,14 @@ void Molecule::read(int flag)
     else if (strstr(line,"mass")) {
       massflag = 1;
       sscanf(line,"%lg",&masstotal);
+      masstotal *= sizescale*sizescale*sizescale;
     }
     else if (strstr(line,"com")) {
       comflag = 1;
       sscanf(line,"%lg %lg %lg",&com[0],&com[1],&com[2]);
+      com[0] *= sizescale;
+      com[1] *= sizescale;
+      com[2] *= sizescale;
       if (domain->dimension == 2 && com[2] != 0.0)
         error->all(FLERR,"Molecule file z center-of-mass must be 0.0 for 2d");
     }
@@ -404,6 +438,12 @@ void Molecule::read(int flag)
       sscanf(line,"%lg %lg %lg %lg %lg %lg",
              &itensor[0],&itensor[1],&itensor[2],
              &itensor[3],&itensor[4],&itensor[5]);
+      itensor[0] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+      itensor[1] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+      itensor[2] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+      itensor[3] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+      itensor[4] *= sizescale*sizescale*sizescale*sizescale*sizescale;
+      itensor[5] *= sizescale*sizescale*sizescale*sizescale*sizescale;
     }
 
     else break;
@@ -414,8 +454,10 @@ void Molecule::read(int flag)
   if (natoms < 1) error->all(FLERR,"No 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) error->all(FLERR,"Invalid dihedral count in molecule file");
+  if (ndihedrals < 0) 
-  if (nimpropers < 0) error->all(FLERR,"Invalid improper count in molecule file");
+    error->all(FLERR,"Invalid dihedral count in molecule file");
+  if (nimpropers < 0) 
+    error->all(FLERR,"Invalid improper count in molecule file");
 
   // count = vector for tallying bonds,angles,etc per atom
 
@@ -544,6 +586,9 @@ void Molecule::coords(char *line)
         error->all(FLERR,"Invalid Coords section in molecule file");
     }
     sscanf(line,"%d %lg %lg %lg",&tmp,&x[i][0],&x[i][1],&x[i][2]);
+    x[i][0] *= sizescale;
+    x[i][1] *= sizescale;
+    x[i][2] *= sizescale;
   }
 
   if (domain->dimension == 2) {
@@ -614,6 +659,7 @@ void Molecule::diameters(char *line)
         error->all(FLERR,"Invalid Diameters section in molecule file");
     }
     sscanf(line,"%d %lg",&tmp,&radius[i]);
+    radius[i] *= sizescale;
     radius[i] *= 0.5;
     maxradius = MAX(maxradius,radius[i]);
   }
@@ -638,6 +684,7 @@ void Molecule::masses(char *line)
         error->all(FLERR,"Invalid Masses section in molecule file");
     }
     sscanf(line,"%d %lg",&tmp,&rmass[i]);
+    rmass[i] *= sizescale*sizescale*sizescale;
   }
 
   for (int i = 0; i < natoms; i++)

src/molecule.h

@@ -102,7 +102,7 @@ class Molecule : protected Pointers {
   double **dxbody;     // displacement of each atom relative to COM
                        // in body frame (diagonalized interia tensor)
 
-  Molecule(class LAMMPS *, int, char **, int);
+  Molecule(class LAMMPS *, int, char **, int &);
   ~Molecule();
   void compute_center();
   void compute_mass();
@@ -115,6 +115,8 @@ class Molecule : protected Pointers {
   FILE *fp;
   int *count;
   int toffset,boffset,aoffset,doffset,ioffset;
+  int autospecial;
+  double sizescale;
 
   void read(int);
   void coords(char *);

src/read_data.cpp

@@ -279,7 +279,7 @@ void ReadData::command(int narg, char **arg)
 
   // -----------------------------------------------------------------
 
-  // perform 1-pass read if no molecular topoogy in file
+  // perform 1-pass read if no molecular topology in file
   // perform 2-pass read if molecular topology,
   //   first pass calculates max topology/atom
 
@@ -301,7 +301,7 @@ void ReadData::command(int narg, char **arg)
   triclinic = 0;
   keyword[0] = '\0';
 
-  int nlocal_previous = atom->nlocal;
+  nlocal_previous = atom->nlocal;
   int firstpass = 1;
 
   while (1) {
@@ -315,7 +315,7 @@ void ReadData::command(int narg, char **arg)
     
     // read header info
     
-    header();
+    header(firstpass);
 
     // problem setup using info from header
     // only done once, if firstpass and first data file
@@ -648,7 +648,7 @@ void ReadData::command(int narg, char **arg)
     // will also observe extra settings even if bond/etc topology not in file
     // leaves other atom arrays unchanged, since already nmax in length
 
-    atom->deallocate_topology();
+    if (addflag == NONE) atom->deallocate_topology();
     atom->avec->grow(atom->nmax);
   }
 
@@ -788,7 +788,7 @@ void ReadData::command(int narg, char **arg)
    some logic differs if adding atoms
 ------------------------------------------------------------------------- */
 
-void ReadData::header()
+void ReadData::header(int firstpass)
 {
   int n;
   char *ptr;
@@ -856,7 +856,7 @@ void ReadData::header()
     if (strstr(line,"atoms")) {
       sscanf(line,BIGINT_FORMAT,&natoms);
       if (addflag == NONE) atom->natoms = natoms;
-      else atom->natoms += natoms;
+      else if (firstpass) atom->natoms += natoms;
 
     // check for these first
     // otherwise "triangles" will be matched as "angles"
@@ -1093,6 +1093,8 @@ void ReadData::velocities()
 
 void ReadData::bonds(int firstpass)
 {
+  int nchunk,eof;
+
   if (me == 0) {
     if (firstpass) {
       if (screen) fprintf(screen,"  scanning bonds ...\n");
@@ -1114,32 +1116,38 @@ void ReadData::bonds(int firstpass)
 
   // read and process bonds
 
-  int nchunk,eof;
   bigint nread = 0;
-  bigint nbonds = atom->nbonds;
 
   while (nread < nbonds) {
     nchunk = MIN(nbonds-nread,CHUNK);
     eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
     if (eof) error->all(FLERR,"Unexpected end of data file");
-    atom->data_bonds(nchunk,buffer,count,id_offset);
+    atom->data_bonds(nchunk,buffer,count,id_offset,boffset);
     nread += nchunk;
   }
 
   // if firstpass: tally max bond/atom and return
+  // if addflag = NONE, store max bond/atom with extra
+  // else just check actual max does not exceed existing max
 
   if (firstpass) {
     int max = 0;
-    for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
+    for (int i = nlocal_previous; i < nlocal; i++) max = MAX(max,count[i]);
     int maxall;
     MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
-    maxall += atom->extra_bond_per_atom;
+    if (addflag == NONE) maxall += atom->extra_bond_per_atom;
 
     if (me == 0) {
       if (screen) fprintf(screen,"  %d = max bonds/atom\n",maxall);
       if (logfile) fprintf(logfile,"  %d = max bonds/atom\n",maxall);
     }
-    atom->bond_per_atom = maxall;
+
+    if (addflag != NONE) {
+      if (maxall > atom->bond_per_atom) 
+        error->all(FLERR,"Subsequent read data induced "
+                   "too many bonds per atom");
+    } else atom->bond_per_atom = maxall;
+
     memory->destroy(count);
     return;
   }
@@ -1147,7 +1155,7 @@ void ReadData::bonds(int firstpass)
   // if 2nd pass: check that bonds were assigned correctly
 
   bigint n = 0;
-  for (int i = 0; i < nlocal; i++) n += atom->num_bond[i];
+  for (int i = nlocal_previous; i < nlocal; i++) n += atom->num_bond[i];
   bigint sum;
   MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
   int factor = 1;
@@ -1158,7 +1166,7 @@ void ReadData::bonds(int firstpass)
     if (logfile) fprintf(logfile,"  " BIGINT_FORMAT " bonds\n",sum/factor);
   }
 
-  if (sum != factor*atom->nbonds)
+  if (sum != factor*nbonds)
     error->all(FLERR,"Bonds assigned incorrectly");
 }
 
@@ -1168,6 +1176,8 @@ void ReadData::bonds(int firstpass)
 
 void ReadData::angles(int firstpass)
 {
+  int nchunk,eof;
+
   if (me == 0) {
     if (firstpass) {
       if (screen) fprintf(screen,"  scanning angles ...\n");
@@ -1189,32 +1199,38 @@ void ReadData::angles(int firstpass)
 
   // read and process angles
 
-  int nchunk,eof;
   bigint nread = 0;
-  bigint nangles = atom->nangles;
 
   while (nread < nangles) {
     nchunk = MIN(nangles-nread,CHUNK);
     eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
     if (eof) error->all(FLERR,"Unexpected end of data file");
-    atom->data_angles(nchunk,buffer,count,id_offset);
+    atom->data_angles(nchunk,buffer,count,id_offset,aoffset);
     nread += nchunk;
   }
 
   // if firstpass: tally max angle/atom and return
+  // if addflag = NONE, store max angle/atom with extra
+  // else just check actual max does not exceed existing max
 
   if (firstpass) {
     int max = 0;
-    for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
+    for (int i = nlocal_previous; i < nlocal; i++) max = MAX(max,count[i]);
     int maxall;
     MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
-    maxall += atom->extra_angle_per_atom;
+    if (addflag == NONE) maxall += atom->extra_angle_per_atom;
 
     if (me == 0) {
       if (screen) fprintf(screen,"  %d = max angles/atom\n",maxall);
       if (logfile) fprintf(logfile,"  %d = max angles/atom\n",maxall);
     }
-    atom->angle_per_atom = maxall;
+
+    if (addflag != NONE) {
+      if (maxall > atom->angle_per_atom) 
+        error->all(FLERR,"Subsequent read data induced "
+                   "too many angles per atom");
+    } else atom->angle_per_atom = maxall;
+
     memory->destroy(count);
     return;
   }
@@ -1222,7 +1238,7 @@ void ReadData::angles(int firstpass)
   // if 2nd pass: check that angles were assigned correctly
 
   bigint n = 0;
-  for (int i = 0; i < nlocal; i++) n += atom->num_angle[i];
+  for (int i = nlocal_previous; i < nlocal; i++) n += atom->num_angle[i];
   bigint sum;
   MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
   int factor = 1;
@@ -1233,7 +1249,7 @@ void ReadData::angles(int firstpass)
     if (logfile) fprintf(logfile,"  " BIGINT_FORMAT " angles\n",sum/factor);
   }
 
-  if (sum != factor*atom->nangles)
+  if (sum != factor*nangles)
     error->all(FLERR,"Angles assigned incorrectly");
 }
 
@@ -1243,6 +1259,8 @@ void ReadData::angles(int firstpass)
 
 void ReadData::dihedrals(int firstpass)
 {
+  int nchunk,eof;
+
   if (me == 0) {
     if (firstpass) {
       if (screen) fprintf(screen,"  scanning dihedrals ...\n");
@@ -1264,31 +1282,38 @@ void ReadData::dihedrals(int firstpass)
 
   // read and process dihedrals
 
-  int nchunk,eof;
   bigint nread = 0;
-  bigint ndihedrals = atom->ndihedrals;
 
   while (nread < ndihedrals) {
     nchunk = MIN(ndihedrals-nread,CHUNK);
     eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
     if (eof) error->all(FLERR,"Unexpected end of data file");
-    atom->data_dihedrals(nchunk,buffer,count,id_offset);
+    atom->data_dihedrals(nchunk,buffer,count,id_offset,doffset);
     nread += nchunk;
   }
 
   // if firstpass: tally max dihedral/atom and return
+  // if addflag = NONE, store max dihedral/atom with extra
+  // else just check actual max does not exceed existing max
 
   if (firstpass) {
     int max = 0;
     for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
     int maxall;
     MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
-    maxall += atom->extra_dihedral_per_atom;
+    if (addflag == NONE) maxall += atom->extra_dihedral_per_atom;
 
     if (me == 0) {
       if (screen) fprintf(screen,"  %d = max dihedrals/atom\n",maxall);
       if (logfile) fprintf(logfile,"  %d = max dihedrals/atom\n",maxall);
     }
+
+    if (addflag != NONE) {
+      if (maxall > atom->dihedral_per_atom) 
+        error->all(FLERR,"Subsequent read data induced "
+                   "too many dihedrals per atom");
+    } else atom->dihedral_per_atom = maxall;
+
     atom->dihedral_per_atom = maxall;
     memory->destroy(count);
     return;
@@ -1297,7 +1322,7 @@ void ReadData::dihedrals(int firstpass)
   // if 2nd pass: check that dihedrals were assigned correctly
 
   bigint n = 0;
-  for (int i = 0; i < nlocal; i++) n += atom->num_dihedral[i];
+  for (int i = nlocal_previous; i < nlocal; i++) n += atom->num_dihedral[i];
   bigint sum;
   MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
   int factor = 1;
@@ -1308,7 +1333,7 @@ void ReadData::dihedrals(int firstpass)
     if (logfile) fprintf(logfile,"  " BIGINT_FORMAT " dihedrals\n",sum/factor);
   }
 
-  if (sum != factor*atom->ndihedrals)
+  if (sum != factor*ndihedrals)
     error->all(FLERR,"Dihedrals assigned incorrectly");
 }
 
@@ -1318,6 +1343,8 @@ void ReadData::dihedrals(int firstpass)
 
 void ReadData::impropers(int firstpass)
 {
+  int nchunk,eof;
+
   if (me == 0) {
     if (firstpass) {
       if (screen) fprintf(screen,"  scanning impropers ...\n");
@@ -1339,32 +1366,38 @@ void ReadData::impropers(int firstpass)
 
   // read and process impropers
 
-  int nchunk,eof;
   bigint nread = 0;
-  bigint nimpropers = atom->nimpropers;
 
   while (nread < nimpropers) {
     nchunk = MIN(nimpropers-nread,CHUNK);
     eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
     if (eof) error->all(FLERR,"Unexpected end of data file");
-    atom->data_impropers(nchunk,buffer,count,id_offset);
+    atom->data_impropers(nchunk,buffer,count,id_offset,ioffset);
     nread += nchunk;
   }
 
   // if firstpass: tally max improper/atom and return
+  // if addflag = NONE, store max improper/atom
+  // else just check it does not exceed existing max
 
   if (firstpass) {
     int max = 0;
-    for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
+    for (int i = nlocal_previous; i < nlocal; i++) max = MAX(max,count[i]);
     int maxall;
     MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
-    maxall += atom->extra_improper_per_atom;
+    if (addflag == NONE) maxall += atom->extra_improper_per_atom;
 
     if (me == 0) {
       if (screen) fprintf(screen,"  %d = max impropers/atom\n",maxall);
       if (logfile) fprintf(logfile,"  %d = max impropers/atom\n",maxall);
     }
-    atom->improper_per_atom = maxall;
+
+    if (addflag != NONE) {
+      if (maxall > atom->improper_per_atom) 
+        error->all(FLERR,"Subsequent read data induced "
+                   "too many impropers per atom");
+    } else atom->improper_per_atom = maxall;
+
     memory->destroy(count);
     return;
   }
@@ -1372,7 +1405,7 @@ void ReadData::impropers(int firstpass)
   // if 2nd pass: check that impropers were assigned correctly
 
   bigint n = 0;
-  for (int i = 0; i < nlocal; i++) n += atom->num_improper[i];
+  for (int i = nlocal_previous; i < nlocal; i++) n += atom->num_improper[i];
   bigint sum;
   MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
   int factor = 1;
@@ -1383,7 +1416,7 @@ void ReadData::impropers(int firstpass)
     if (logfile) fprintf(logfile,"  " BIGINT_FORMAT " impropers\n",sum/factor);
   }
 
-  if (sum != factor*atom->nimpropers)
+  if (sum != factor*nimpropers)
     error->all(FLERR,"Impropers assigned incorrectly");
 }
 

src/read_data.h

@@ -41,6 +41,7 @@ class ReadData : protected Pointers {
 
   bigint id_offset;
 
+  int nlocal_previous;
   bigint natoms;
   bigint nbonds,nangles,ndihedrals,nimpropers;
   int ntypes;
@@ -81,7 +82,7 @@ class ReadData : protected Pointers {
   void open(char *);
   void scan(int &, int &, int &, int &);
   int reallocate(int **, int, int);
-  void header();
+  void header(int);
   void parse_keyword(int);
   void skip_lines(bigint);
   void parse_coeffs(char *, const char *, int, int, int);

src/respa.cpp

@@ -733,6 +733,7 @@ void Respa::recurse(int ilevel)
       modify->pre_reverse(eflag,vflag);
       timer->stamp(Timer::MODIFY);
     }
+
     if (newton[ilevel]) {
       comm->reverse_comm();
       timer->stamp(Timer::COMM);

src/version.h

@@ -1 +1 @@
-#define LAMMPS_VERSION "22 Oct 2015"
+#define LAMMPS_VERSION "23 Oct 2015"