git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8199 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/atom.cpp

@@ -402,7 +402,8 @@ void Atom::modify_params(int narg, char **arg)
       else if (strcmp(arg[iarg+1],"hash") == 0) map_style = 2;
       else error->all(FLERR,"Illegal atom_modify command");
       if (domain->box_exist) 
-	error->all(FLERR,"Atom_modify map command after simulation box is defined");
+	error->all(FLERR,
+		   "Atom_modify map command after simulation box is defined");
       iarg += 2;
     } else if (strcmp(arg[iarg],"first") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
@@ -670,7 +671,7 @@ int Atom::map_find_hash(int global)
 /* ----------------------------------------------------------------------
    add unique tags to any atoms with tag = 0
    new tags are grouped by proc and start after max current tag
-   called after creating new atoms 
+   called after creating new atoms
 ------------------------------------------------------------------------- */
 
 void Atom::tag_extend()
@@ -1051,7 +1052,8 @@ void Atom::data_dihedrals(int n, char *buf)
 	atom4 <= 0 || atom4 > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Dihedrals section of data file");
     if (itype <= 0 || itype > ndihedraltypes)
-      error->one(FLERR,"Invalid dihedral type in Dihedrals section of data file");
+      error->one(FLERR,
+		 "Invalid dihedral type in Dihedrals section of data file");
     if ((m = map(atom2)) >= 0) {
       dihedral_type[m][num_dihedral[m]] = itype;
       dihedral_atom1[m][num_dihedral[m]] = atom1;
@@ -1110,7 +1112,8 @@ void Atom::data_impropers(int n, char *buf)
 	atom4 <= 0 || atom4 > map_tag_max)
       error->one(FLERR,"Invalid atom ID in Impropers section of data file");
     if (itype <= 0 || itype > nimpropertypes)
-      error->one(FLERR,"Invalid improper type in Impropers section of data file");
+      error->one(FLERR,
+		 "Invalid improper type in Impropers section of data file");
     if ((m = map(atom2)) >= 0) {
       improper_type[m][num_improper[m]] = itype;
       improper_atom1[m][num_improper[m]] = atom1;
@@ -1177,7 +1180,8 @@ void Atom::set_mass(const char *str)
   int n = sscanf(str,"%d %lg",&itype,&mass_one);
   if (n != 2) error->all(FLERR,"Invalid mass line in data file");
 
-  if (itype < 1 || itype > ntypes) error->all(FLERR,"Invalid type for mass set");
+  if (itype < 1 || itype > ntypes) 
+    error->all(FLERR,"Invalid type for mass set");
 
   mass[itype] = mass_one;
   mass_setflag[itype] = 1;
@@ -1193,7 +1197,8 @@ void Atom::set_mass(const char *str)
 void Atom::set_mass(int itype, double value)
 {
   if (mass == NULL) error->all(FLERR,"Cannot set mass for this atom style");
-  if (itype < 1 || itype > ntypes) error->all(FLERR,"Invalid type for mass set");
+  if (itype < 1 || itype > ntypes) 
+    error->all(FLERR,"Invalid type for mass set");
 
   mass[itype] = value;
   mass_setflag[itype] = 1;

src/create_atoms.cpp

@@ -28,8 +28,8 @@
 
 using namespace LAMMPS_NS;
 
-#define BIG      1.0e30
+#define BIG 1.0e30
-#define EPSILON  1.0e-6
+#define EPSILON 1.0e-6
 
 enum{BOX,REGION,SINGLE,RANDOM};
 
@@ -192,7 +192,6 @@ void CreateAtoms::command(int narg, char **arg)
   // reset simulation now that more atoms are defined
   // add tags for newly created atoms if possible
   // if global map exists, reset it
-  // if a molecular system, set nspecial to 0 for new atoms
 
   if (atom->natoms > MAXTAGINT) atom->tag_enable = 0;
   if (atom->natoms <= MAXTAGINT) atom->tag_extend();
@@ -202,14 +201,18 @@ void CreateAtoms::command(int narg, char **arg)
     atom->map_init();
     atom->map_set();
   }
-  if (atom->molecular) {
+
-    int **nspecial = atom->nspecial;
+  // if a molecular system, set nspecial to 0 for new atoms
-    for (int i = nlocal_previous; i < atom->nlocal; i++) {
+  // NOTE: 31May12, don't think this is needed, avec->create_atom() does it
-      nspecial[i][0] = 0;
+
-      nspecial[i][1] = 0;
+  //if (atom->molecular) {
-      nspecial[i][2] = 0;
+  //  int **nspecial = atom->nspecial;
-    }
+  //  for (int i = nlocal_previous; i < atom->nlocal; i++) {
-  }
+  //    nspecial[i][0] = 0;
+  //    nspecial[i][1] = 0;
+  //    nspecial[i][2] = 0;
+  //  }
+  //}
 }
 
 /* ----------------------------------------------------------------------

src/read_dump.cpp

@@ -0,0 +1,672 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Timothy Sirk (U Vermont)
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "string.h"
+#include "stdlib.h"
+#include "read_dump.h"
+#include "read_dump_native.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "update.h"
+#include "domain.h"
+#include "comm.h"
+#include "irregular.h"
+#include "error.h"
+#include "memory.h"
+
+using namespace LAMMPS_NS;
+
+#define CHUNK 1024
+#define EPSILON 1.0e-6
+
+enum{ID,TYPE,X,Y,Z,VX,VY,VZ,IX,IY,IZ};
+enum{UNSET,UNSCALED,SCALED};
+enum{NATIVE};
+
+/* ---------------------------------------------------------------------- */
+
+ReadDump::ReadDump(LAMMPS *lmp) : Pointers(lmp)
+{
+  MPI_Comm_rank(world,&me);
+  MPI_Comm_size(world,&nprocs);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ReadDump::command(int narg, char **arg)
+{
+  if (narg < 2) error->all(FLERR,"Illegal read_dump command");
+
+  nstep = ATOBIGINT(arg[1]);
+
+  // per-field vectors
+
+  int firstfield = 2;
+  fieldtype = new int[narg];
+  fieldlabel = new char*[narg];
+
+  // scan ahead to see if "add yes" keyword/value is used
+  // requires extra "type" field from from dump file
+  // add id and type fields as needed
+
+  int iarg;
+  for (iarg = firstfield; iarg < narg; iarg++)
+    if (strcmp(arg[iarg],"add") == 0)
+      if (iarg < narg-1 && strcmp(arg[iarg+1],"yes") == 0) break;
+
+  nfield = 0;
+  fieldtype[nfield++] = ID;
+  if (iarg < narg) fieldtype[nfield++] = TYPE;
+
+  // parse fields
+
+  iarg = firstfield;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"x") == 0) fieldtype[nfield++] = X;
+    else if (strcmp(arg[iarg],"y") == 0) fieldtype[nfield++] = Y;
+    else if (strcmp(arg[iarg],"z") == 0) fieldtype[nfield++] = Z;
+    else if (strcmp(arg[iarg],"vx") == 0) fieldtype[nfield++] = VX;
+    else if (strcmp(arg[iarg],"vy") == 0) fieldtype[nfield++] = VY;
+    else if (strcmp(arg[iarg],"vz") == 0) fieldtype[nfield++] = VZ;
+    else if (strcmp(arg[iarg],"ix") == 0) fieldtype[nfield++] = IX;
+    else if (strcmp(arg[iarg],"iy") == 0) fieldtype[nfield++] = IY;
+    else if (strcmp(arg[iarg],"iz") == 0) fieldtype[nfield++] = IZ;
+    else break;
+    iarg++;
+  }
+
+  dimension = domain->dimension;
+  triclinic = domain->triclinic;
+
+  if (fieldtype[nfield-1] == ID || fieldtype[nfield-1] == TYPE)
+    error->all(FLERR,"Illegal read_dump command");
+
+  if (dimension == 2) {
+    for (int i = 0; i < nfield; i++)
+      if (fieldtype[i] == Z || fieldtype[i] == VZ || fieldtype[i] == IZ)
+	error->all(FLERR,"Illegal read_dump command");
+  }
+
+  for (int i = 0; i < nfield; i++)
+    for (int j = i+1; j < nfield; j++)
+      if (fieldtype[i] == fieldtype[j])
+	error->all(FLERR,"Duplicate fields in read_dump command");
+
+  // parse optional args
+
+  boxflag = 1;
+  replaceflag = 1;
+  purgeflag = 0;
+  trimflag = 0;
+  addflag = 0;
+  for (int i = 0; i < nfield; i++) fieldlabel[i] = NULL;
+  scaledflag = UNSCALED;
+  format = NATIVE;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"box") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
+      if (strcmp(arg[iarg+1],"yes") == 0) boxflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) boxflag = 0;
+      else error->all(FLERR,"Illegal read_dump command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"replace") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
+      if (strcmp(arg[iarg+1],"yes") == 0) replaceflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) replaceflag = 0;
+      else error->all(FLERR,"Illegal read_dump command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"purge") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
+      if (strcmp(arg[iarg+1],"yes") == 0) purgeflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) purgeflag = 0;
+      else error->all(FLERR,"Illegal read_dump command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"trim") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
+      if (strcmp(arg[iarg+1],"yes") == 0) trimflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) trimflag = 0;
+      else error->all(FLERR,"Illegal read_dump command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"add") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
+      if (strcmp(arg[iarg+1],"yes") == 0) addflag = 1;
+      else if (strcmp(arg[iarg+1],"no") == 0) addflag = 0;
+      else error->all(FLERR,"Illegal read_dump command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"label") == 0) {
+      if (iarg+3 > narg) error->all(FLERR,"Illegal read_dump command");
+      int i;
+      for (i = 0; i < nfield; i++)
+	if (strcmp(arg[firstfield+i],arg[iarg+1]) == 0) break;
+      if (i == nfield) error->all(FLERR,"Illegal read_dump command");
+      fieldlabel[i] = arg[iarg+2];
+      iarg += 3;
+    } else if (strcmp(arg[iarg],"scaled") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
+      if (strcmp(arg[iarg+1],"yes") == 0) scaledflag = SCALED;
+      else if (strcmp(arg[iarg+1],"no") == 0) scaledflag = UNSCALED;
+      else error->all(FLERR,"Illegal read_dump command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"format") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
+      if (strcmp(arg[iarg+1],"native") == 0) format = NATIVE;
+      else error->all(FLERR,"Illegal read_dump command");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal read_dump command");
+  }
+
+  if (purgeflag && (replaceflag || trimflag))
+    error->all(FLERR,"If read_dump purges it cannot replace or trim");
+
+  // allocate snapshot field buffer
+
+  memory->create(fields,CHUNK,nfield,"read_dump:fields");
+
+  // create reader class
+  // could make this a parent class and customize with other readers
+
+  if (format == NATIVE) reader = new ReadDumpNative(lmp);
+
+  // proc 0 opens dump file and scans to correct snapshot
+  // after scan these values are set, so Bcast them:
+  // nsnapatoms, box[3][3], scaled
+  // NOTE: fieldlabel is just ptrs to input args in read_dump command
+  //       will not persist if want to use labels in rerun() command
+
+  if (me == 0) {
+    if (screen) fprintf(screen,"Scanning dump file ...\n");
+    open(arg[0]);
+    reader->init(fp);
+    reader->scan(nstep,nfield,fieldtype,fieldlabel,scaledflag,
+		 nsnapatoms,box,scaled);
+  }
+
+  MPI_Bcast(&nsnapatoms,1,MPI_LMP_BIGINT,0,world);
+  MPI_Bcast(&box[0][0],9,MPI_DOUBLE,0,world);
+  MPI_Bcast(&scaled,1,MPI_INT,0,world);
+
+  // for scaled coords and triclinic box:
+  // yindex,zindex = index of Y and Z fields
+  // already known to exist because checked in scan()
+  // needed for unscaling to absolute coords in xfield(), yfield(), zfield()
+
+  if (scaled == SCALED && triclinic) {
+    for (int i = 0; i < nfield; i++) {
+      if (fieldtype[i] == Y) yindex = i;
+      if (fieldtype[i] == Z) zindex = i;
+    }
+  }
+
+  // make local copy of snapshot box params
+
+  xlo = box[0][0];
+  xhi = box[0][1];
+  ylo = box[1][0];
+  yhi = box[1][1];
+  zlo = box[2][0];
+  zhi = box[2][1];
+  xprd = xhi - xlo;
+  yprd = yhi - ylo;
+  zprd = zhi - zlo;
+  if (triclinic) {
+    xy = box[0][2];
+    xz = box[1][2];
+    yz = box[2][2];
+  }
+
+  // reset timestep to nstep
+
+  char *tstr[1];
+  char str[32];
+  sprintf(str,BIGINT_FORMAT,nstep);
+  tstr[0] = str;
+  update->reset_timestep(1,tstr);
+
+  // reset simulation box from snapshot box parameters if requested
+  // do it now, so if adding atoms, procs will have correct sub-domains
+  // call domain->reset_box() later,
+  //   since can't shrink wrap until atom coords change and atoms are added
+
+  if (boxflag) {
+    domain->boxlo[0] = xlo;
+    domain->boxhi[0] = xhi;
+    domain->boxlo[1] = ylo;
+    domain->boxhi[1] = yhi;
+    if (dimension == 3) {
+      domain->boxlo[2] = zlo;
+      domain->boxhi[2] = zhi;
+    }
+    if (triclinic) {
+      domain->xy = xy;
+      if (dimension == 3) {
+	domain->xz = xz;
+	domain->yz = yz;
+      }
+    }
+
+    domain->set_initial_box();
+    domain->set_global_box();
+    comm->set_proc_grid();
+    domain->set_local_box();
+  }
+
+  // read in the snapshot
+
+  if (me == 0)
+    if (screen) fprintf(screen,"Reading snapshot from dump file ...\n");
+
+  // counters
+
+  bigint natoms_prev = atom->natoms;
+  npurge = nreplace = ntrim = nadd = 0;
+
+  // if purgeflag set, delete all current atoms
+
+  if (purgeflag) {
+    if (atom->map_style) atom->map_clear();
+    npurge = atom->nlocal;
+    atom->nlocal = atom->nghost = 0;
+    atom->natoms = 0;
+  }
+
+  // to match existing atoms to dump atoms:
+  // must build map if not a molecular system
+
+  int mapflag = 0;
+  if (atom->map_style == 0) {
+    mapflag = 1;
+    atom->map_style = 1;
+    atom->map_init();
+    atom->map_set();
+  }
+
+  // uflag[i] = 1 for each owned atom appearing in dump
+  // ucflag = similar flag for each chunk atom, used in process_atoms()
+
+  int nlocal = atom->nlocal;
+  memory->create(uflag,nlocal,"read_dump:uflag");
+  for (int i = 0; i < nlocal; i++) uflag[i] = 0;
+  memory->create(ucflag,CHUNK,"read_dump:ucflag");
+  memory->create(ucflag_all,CHUNK,"read_dump:ucflag");
+  
+  // read, broadcast, and process atoms from snapshot in chunks
+
+  addproc = -1;
+
+  int nchunk;
+  bigint nread = 0;
+  while (nread < nsnapatoms) {
+    nchunk = MIN(nsnapatoms-nread,CHUNK);
+    if (me == 0) reader->read(nchunk,fields);
+    MPI_Bcast(&fields[0][0],nchunk*nfield,MPI_DOUBLE,0,world);
+    process_atoms(nchunk);
+    nread += nchunk;
+  }
+
+  // if addflag set, add tags to new atoms if possible
+
+  if (addflag) {
+    bigint nblocal = atom->nlocal;
+    MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
+    if (atom->natoms < 0 || atom->natoms > MAXBIGINT)
+      error->all(FLERR,"Too many total atoms");
+    if (atom->natoms > MAXTAGINT) atom->tag_enable = 0;
+    if (atom->natoms <= MAXTAGINT) atom->tag_extend();
+  }
+
+  // if trimflag set, delete atoms not replaced by snapshot atoms
+
+  if (trimflag) {
+    delete_atoms(uflag);
+    bigint nblocal = atom->nlocal;
+    MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
+  }
+
+  // delete atom map if created it above
+  // else reinitialize map for current atoms
+  // do this before migrating atoms to new procs via Irregular
+
+  if (mapflag) {
+    atom->map_delete();
+    atom->map_style = 0;
+  } else {
+    atom->nghost = 0;
+    atom->map_init();
+    atom->map_set();
+  }
+
+  // close dump file
+
+  if (me == 0) {
+    if (compressed) pclose(fp);
+    else fclose(fp);
+  }
+
+  // move atoms back inside simulation box and to new processors
+  // use remap() instead of pbc() in case atoms moved a long distance
+  // adjust image flags of all atoms (old and new) based on current box
+  // use irregular() in case atoms moved a long distance
+
+  double **x = atom->x;
+  int *image = atom->image;
+  nlocal = atom->nlocal;
+  for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
+
+  if (triclinic) domain->x2lamda(atom->nlocal);
+  domain->reset_box();
+  Irregular *irregular = new Irregular(lmp);
+  irregular->migrate_atoms();
+  delete irregular;
+  if (triclinic) domain->lamda2x(atom->nlocal);
+
+  domain->print_box("  ");
+
+  // clean up
+
+  delete reader;
+  delete [] fieldtype;
+  delete [] fieldlabel;
+  memory->destroy(fields);
+  memory->destroy(uflag);
+  memory->destroy(ucflag);
+  memory->destroy(ucflag_all);
+
+  // print out stats
+
+  bigint npurge_all,nreplace_all,ntrim_all,nadd_all;
+
+  bigint tmp;
+  tmp = npurge;
+  MPI_Allreduce(&tmp,&npurge_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
+  tmp = nreplace;
+  MPI_Allreduce(&tmp,&nreplace_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
+  tmp = ntrim;
+  MPI_Allreduce(&tmp,&ntrim_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
+  tmp = nadd;
+  MPI_Allreduce(&tmp,&nadd_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
+
+  if (me == 0) {
+    if (screen) {
+      fprintf(screen,"  " BIGINT_FORMAT " atoms before read\n",natoms_prev);
+      fprintf(screen,"  " BIGINT_FORMAT " atoms in snapshot\n",nsnapatoms);
+      fprintf(screen,"  " BIGINT_FORMAT " atoms purged\n",npurge_all);
+      fprintf(screen,"  " BIGINT_FORMAT " atoms replaced\n",nreplace_all);
+      fprintf(screen,"  " BIGINT_FORMAT " atoms trimmed\n",ntrim_all);
+      fprintf(screen,"  " BIGINT_FORMAT " atoms added\n",nadd_all);
+      fprintf(screen,"  " BIGINT_FORMAT " atoms after read\n",atom->natoms);
+    }
+    if (logfile) {
+      fprintf(logfile,"  " BIGINT_FORMAT " atoms before read\n",natoms_prev);
+      fprintf(logfile,"  " BIGINT_FORMAT " atoms in snapshot\n",nsnapatoms);
+      fprintf(logfile,"  " BIGINT_FORMAT " atoms purged\n",npurge_all);
+      fprintf(logfile,"  " BIGINT_FORMAT " atoms replaced\n",nreplace_all);
+      fprintf(logfile,"  " BIGINT_FORMAT " atoms trimmed\n",ntrim_all);
+      fprintf(logfile,"  " BIGINT_FORMAT " atoms added\n",nadd_all);
+      fprintf(logfile,"  " BIGINT_FORMAT " atoms after read\n",atom->natoms);
+    }
+  }
+}
+
+
+/* ----------------------------------------------------------------------
+   process each of N atoms in chunk read from dump file
+   if in replace mode and atom ID matches current atom,
+     overwrite atom info with fields from dump file
+   if in add mode and atom ID does not match any current atom,
+     create new atom with dump file field values,
+     and assign to a proc in round-robin manner
+   use round-robin method, b/c atom coords may not be inside simulation box
+------------------------------------------------------------------------- */
+
+void ReadDump::process_atoms(int n)
+{
+  int i,m,ifield,itype;
+  int xbox,ybox,zbox;
+
+  double **x = atom->x;
+  double **v = atom->v;
+  int *image = atom->image;
+  int nlocal = atom->nlocal;
+
+  for (i = 0; i < n; i++) {
+    ucflag[i] = 0;
+
+    // map() call is invalid if purged all atoms
+    // setting m = -1 forces new atom not to match
+
+    if (!purgeflag) m = atom->map(static_cast<int> (fields[i][0]));
+    else m = -1;
+    if (m < 0 || m >= nlocal) continue;
+
+    ucflag[i] = 1;
+    uflag[m] = 1;
+
+    if (replaceflag) {
+      nreplace++;
+
+      // current image flags
+      
+      xbox = (image[m] & 1023) - 512;
+      ybox = (image[m] >> 10 & 1023) - 512;
+      zbox = (image[m] >> 20) - 512;
+    
+      // overwrite atom attributes with field info
+      // start from field 1 since 0 = id, 1 will be skipped if type
+      
+      for (ifield = 1; ifield < nfield; ifield++) {
+	switch (fieldtype[ifield]) {
+	case X:
+	  x[m][0] = xfield(i,ifield);
+	  break;
+	case Y:
+	  x[m][1] = yfield(i,ifield);
+	  break;
+	case Z:
+	  x[m][2] = zfield(i,ifield);
+	  break;
+	case VX:
+	  v[m][0] = fields[i][ifield];
+	  break;
+	case VY:
+	  v[m][1] = fields[i][ifield];
+	  break;
+	case VZ:
+	  v[m][2] = fields[i][ifield];
+	  break;
+	case IX:
+	  xbox = static_cast<int> (fields[i][ifield]);
+	  break;
+	case IY:
+	  ybox = static_cast<int> (fields[i][ifield]);
+	  break;
+	case IZ:
+	  zbox = static_cast<int> (fields[i][ifield]);
+	  break;
+	}
+      }
+      
+      // replace image flag in case changed by ix,iy,iz fields
+      
+      image[m] = (xbox << 20) | (ybox << 10) | zbox;
+    }
+  }
+
+  // create any atoms in chunk that no processor owned
+  // add atoms in round-robin sequence on processors
+  // cannot do it geometrically b/c dump coords may not be in simulation box
+
+  if (!addflag) return;
+
+  MPI_Allreduce(ucflag,ucflag_all,n,MPI_INT,MPI_SUM,world);
+
+  double lamda[3],one[3];
+  double *coord;
+
+  for (i = 0; i < n; i++) {
+    if (ucflag_all[i]) continue;
+
+    // each processor adds every Pth atom
+
+    addproc++;
+    if (addproc == nprocs) addproc = 0;
+    if (addproc != me) continue;
+
+    // create type and coord fields from dump file
+    // coord = 0.0 unless corresponding dump file field was specified
+
+    one[0] = one[1] = one[2] = 0.0;
+    for (ifield = 1; ifield < nfield; ifield++) {
+      switch (fieldtype[ifield]) {
+      case TYPE:
+	itype = static_cast<int> (fields[i][ifield]);
+	break;
+      case X:
+	one[0] = xfield(i,ifield);
+	break;
+      case Y:
+	one[1] = yfield(i,ifield);
+	break;
+      case Z:
+	one[2] = zfield(i,ifield);
+	break;
+      }
+    }
+
+    // create the atom on proc that owns it
+    // reset v,image ptrs in case they are reallocated
+
+    atom->avec->create_atom(itype,one);
+    nadd++;
+    
+    v = atom->v;
+    image = atom->image;
+    m = atom->nlocal;
+
+    // set atom attributes from other dump file fields
+    // xyzbox = 512 is default value set by create_atom()
+
+    xbox = ybox = zbox = 512;
+
+    for (ifield = 1; ifield < nfield; ifield++) {
+      switch (fieldtype[ifield]) {
+      case VX:
+	v[m][0] = fields[i][ifield];
+	break;
+      case VY:
+	v[m][1] = fields[i][ifield];
+	break;
+      case VZ:
+	v[m][2] = fields[i][ifield];
+	break;
+      case IX:
+	xbox = static_cast<int> (fields[i][ifield]);
+	break;
+      case IY:
+	ybox = static_cast<int> (fields[i][ifield]);
+	break;
+      case IZ:
+	zbox = static_cast<int> (fields[i][ifield]);
+	break;
+      }
+
+      // replace image flag in case changed by ix,iy,iz fields
+      
+      image[m] = (xbox << 20) | (ybox << 10) | zbox;
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   delete atoms not flagged as replaced by dump atoms
+------------------------------------------------------------------------- */
+
+void ReadDump::delete_atoms(int *uflag)
+{
+  AtomVec *avec = atom->avec;
+  int nlocal = atom->nlocal;
+
+  int i = 0;
+  while (i < nlocal) {
+    if (uflag[i] == 0) {
+      avec->copy(nlocal-1,i,1);
+      uflag[i] = uflag[nlocal-1];
+      nlocal--;
+      ntrim++;
+    } else i++;
+  }
+
+  atom->nlocal = nlocal;
+}
+
+/* ----------------------------------------------------------------------
+   convert XYZ fields in dump file into absolute, unscaled coordinates
+   depends on scaled vs unscaled and triclinic vs orthogonal
+   does not depend on wrapped vs unwrapped
+------------------------------------------------------------------------- */
+
+double ReadDump::xfield(int i, int j)
+{
+  if (scaled == UNSCALED) return fields[i][j];
+  else if (!triclinic) return fields[i][j]*xprd + xlo;
+  else if (dimension == 2) 
+    return xprd*fields[i][j] + xy*fields[i][yindex] + xlo;
+  return xprd*fields[i][j] + xy*fields[i][yindex] + xz*fields[i][zindex] + xlo;
+}
+
+double ReadDump::yfield(int i, int j)
+{
+  if (scaled == UNSCALED) return fields[i][j];
+  else if (!triclinic) return fields[i][j]*yprd + ylo;
+  else if (dimension == 2) return yprd*fields[i][j] + ylo;
+  return yprd*fields[i][j] + yz*fields[i][zindex] + ylo;
+}
+
+double ReadDump::zfield(int i, int j)
+{
+  if (scaled == UNSCALED) return fields[i][j];
+  return fields[i][j]*zprd + zlo;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 opens dump file
+   test if gzipped
+------------------------------------------------------------------------- */
+
+void ReadDump::open(char *file)
+{
+  compressed = 0;
+  char *suffix = file + strlen(file) - 3;
+  if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1;
+  if (!compressed) fp = fopen(file,"r");
+  else {
+#ifdef LAMMPS_GZIP
+    char gunzip[128];
+    sprintf(gunzip,"gunzip -c %s",file);
+    fp = popen(gunzip,"r");
+#else
+    error->one(FLERR,"Cannot open gzipped file");
+#endif
+  }
+
+  if (fp == NULL) {
+    char str[128];
+    sprintf(str,"Cannot open file %s",file);
+    error->one(FLERR,str);
+  }
+}

src/read_dump.h

@@ -0,0 +1,84 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://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.
+
+   Contributed by Timothy Sirk
+------------------------------------------------------------------------- */
+
+#ifdef COMMAND_CLASS
+
+CommandStyle(read_dump,ReadDump)
+
+#else
+
+#ifndef LMP_READ_DUMP_H
+#define LMP_READ_DUMP_H
+
+#include "stdio.h"
+#include "pointers.h"
+
+namespace LAMMPS_NS {
+
+class ReadDump : protected Pointers {
+ public:
+  ReadDump(class LAMMPS *);
+  ~ReadDump() {}
+  void command(int, char **);
+
+private:
+  int me,nprocs;
+  FILE *fp;
+  int dimension;
+  int triclinic;
+
+  bigint nstep;            // timestep to find in dump file
+  int boxflag;             // use dump file box params
+  int replaceflag,addflag; // flags for processing dump snapshot atoms
+  int trimflag,purgeflag;
+  int scaledflag;          // user setting for coordinate scaling
+  int scaled;              // actual setting for coordinate scaling
+  int format;              // style of dump file
+  int compressed;          // flag for dump file compression
+
+  int nfield;              // # of fields to extract from dump file
+  int *fieldtype;          // type of each field = X,VY,IZ,etc
+  char **fieldlabel;       // user specified label for field
+  double **fields;         // per-atom field values
+
+  double box[3][3];         // dump file box parameters
+  double xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz;
+  double xprd,yprd,zprd;
+
+  bigint nsnapatoms;        // # of atoms in dump file shapshot
+
+  int npurge,nreplace,ntrim,nadd;     // stats on processed atoms
+  int addproc;                        // proc that should add next atom
+  int yindex,zindex;                  // field index for Y,Z coords
+
+  int *uflag;               // set to 1 if snapshot atom matches owned atom
+  int *ucflag,*ucflag_all;  // set to 1 if snapshot chunk atom was processed
+
+  class ReadDumpNative *reader;     // class that reads native dump file
+
+  void process_atoms(int);
+  void delete_atoms(int *);
+
+  double xfield(int, int);
+  double yfield(int, int);
+  double zfield(int, int);
+
+  void open(char *);
+};
+
+}
+
+#endif
+#endif

src/read_dump_native.cpp

@@ -0,0 +1,341 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://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 "string.h"
+#include "stdlib.h"
+#include "read_dump_native.h"
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MAXLINE 1024        // max line length in dump file
+
+enum{ID,TYPE,X,Y,Z,VX,VY,VZ,IX,IY,IZ};
+enum{UNSET,UNSCALED,SCALED};
+
+/* ---------------------------------------------------------------------- */
+
+ReadDumpNative::ReadDumpNative(LAMMPS *lmp) : Pointers(lmp)
+{
+  dimension = domain->dimension;
+  triclinic = domain->triclinic;
+
+  line = new char[MAXLINE];
+  words = NULL;
+  fieldindex = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ReadDumpNative::~ReadDumpNative()
+{
+  delete [] line;
+  delete [] words;
+  delete [] fieldindex;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ReadDumpNative::init(FILE *fpcaller)
+{
+  fp = fpcaller;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 scans dump file until reaching snapshot with timestamp = Nstep
+   extract natoms and box bounds from snapshot
+   set fieldindex for specified fields and overall scaled setting
+   error check on current vs new box and fields
+   NOTE: error checking should maybe be moved to parent
+------------------------------------------------------------------------- */
+
+void ReadDumpNative::scan(bigint nstep, 
+			  int nfield_caller, int *fieldtype, 
+			  char **fieldlabel, int scaledflag,
+			  bigint &natoms, double box[3][3], int &scaled)
+{
+  int nchunk,triclinic_snap,s_index,u_index,su_index;
+  bigint ntimestep,nremain;
+  char *bc,*names;
+
+  nfield = nfield_caller;
+
+  read_lines(1);
+
+  while (1) {
+    if (strstr(line,"ITEM: TIMESTEP") != line)
+      error->one(FLERR,"Incorrectly formatted dump file");
+    read_lines(1);
+    sscanf(line,BIGINT_FORMAT,&ntimestep);
+
+    if (ntimestep > nstep)
+      error->one(FLERR,"Dump file does not contain requested snapshot");
+
+    read_lines(2);
+    sscanf(line,BIGINT_FORMAT,&natoms);
+
+    // skip snapshot
+    // invoke read_lines() in chunks no larger than MAXSMALLINT
+
+    if (ntimestep < nstep) {
+      read_lines(5);
+
+      nremain = natoms;
+      while (nremain) {
+	nchunk = MIN(nremain,MAXSMALLINT);
+	read_lines(nchunk);
+	nremain -= nchunk;
+      }
+
+      read_lines(1);
+    } else break;
+  }
+
+  // found correct snapshot
+  // read box size and boundary conditions
+
+  triclinic_snap = 0;
+  box[0][2] = box[1][2] = box[2][2] = 0.0;
+  read_lines(1);
+  bc = &line[strlen("ITEM: BOX BOUNDS ")];
+  if (bc[0] == 'x') {
+    triclinic_snap = 1;
+    bc = &bc[9];
+  }
+  
+  char boundstr[9];
+  domain->boundary_string(boundstr);
+  if (strstr(bc,boundstr) != bc) 
+    error->warning(FLERR,"Read_dump boundary flags do not match simulation");
+  
+  read_lines(1);
+  if (!triclinic_snap) sscanf(line,"%lg %lg",&box[0][0],&box[0][1]);
+  else sscanf(line,"%lg %lg %lg",&box[0][0],&box[0][1],&box[0][2]);
+  read_lines(1);
+  if (!triclinic_snap) sscanf(line,"%lg %lg",&box[1][0],&box[1][1]);
+  else sscanf(line,"%lg %lg %lg",&box[1][0],&box[1][1],&box[1][2]);
+  read_lines(1);
+  if (!triclinic_snap) sscanf(line,"%lg %lg",&box[2][0],&box[2][1]);
+  else sscanf(line,"%lg %lg %lg",&box[2][0],&box[2][1],&box[2][2]);
+  
+  // read ITEM: ATOMS line
+  // labels = column labels
+  
+  read_lines(1);
+  names = &line[strlen("ITEM: ATOMS ")];
+  
+  nwords = atom->count_words(names);
+  char **labels = new char*[nwords];
+  labels[0] = strtok(names," \t\n\r\f");
+  if (labels[0] == NULL) 
+    error->one(FLERR,"Incorrect atom format in dump file");
+  for (int m = 1; m < nwords; m++) {
+    labels[m] = strtok(NULL," \t\n\r\f");
+    if (labels[m] == NULL) 
+      error->one(FLERR,"Incorrect atom format in dump file");
+  }
+  
+  // match each field with column
+  
+  fieldindex = new int[nfield];
+  int xflag = UNSET;
+  int yflag = UNSET;
+  int zflag = UNSET;
+  
+  for (int i = 0; i < nfield; i++) {
+    if (fieldlabel[i]) {
+      fieldindex[i] = find_label(fieldlabel[i],nwords,labels);
+      if (fieldtype[i] == X) xflag = scaledflag;
+      else if (fieldtype[i] == Y) yflag = scaledflag;
+      else if (fieldtype[i] == Z) zflag = scaledflag;
+    }
+    
+    else if (fieldtype[i] == ID)
+      fieldindex[i] = find_label("id",nwords,labels);
+    else if (fieldtype[i] == TYPE)
+      fieldindex[i] = find_label("type",nwords,labels);
+    
+    else if (fieldtype[i] == X) {
+      fieldindex[i] = find_label("x",nwords,labels);
+      xflag = UNSCALED;
+      if (fieldindex[i] < 0) {
+	fieldindex[i] = nwords;
+	s_index = find_label("xs",nwords,labels);
+	u_index = find_label("xu",nwords,labels);
+	su_index = find_label("xsu",nwords,labels);
+	if (s_index >= 0 && s_index < fieldindex[i]) {
+	  fieldindex[i] = s_index;
+	  xflag = SCALED;
+	}
+	if (u_index >= 0 && u_index < fieldindex[i]) {
+	  fieldindex[i] = u_index;
+	  xflag = UNSCALED;
+	}
+	if (su_index >= 0 && su_index < fieldindex[i]) {
+	  fieldindex[i] = su_index;
+	  xflag = SCALED;
+	}
+      }
+      if (fieldindex[i] == nwords) fieldindex[i] = -1;
+      
+    } else if (fieldtype[i] == Y) {
+      fieldindex[i] = find_label("y",nwords,labels);
+      yflag = UNSCALED;
+      if (fieldindex[i] < 0) {
+	fieldindex[i] = nwords;
+	s_index = find_label("ys",nwords,labels);
+	u_index = find_label("yu",nwords,labels);
+	su_index = find_label("ysu",nwords,labels);
+	if (s_index >= 0 && s_index < fieldindex[i]) {
+	  fieldindex[i] = s_index;
+	  yflag = SCALED;
+	}
+	if (u_index >= 0 && u_index < fieldindex[i]) {
+	  fieldindex[i] = u_index;
+	  yflag = UNSCALED;
+	}
+	if (su_index >= 0 && su_index < fieldindex[i]) {
+	  fieldindex[i] = su_index;
+	  yflag = SCALED;
+	}
+      }
+      if (fieldindex[i] == nwords) fieldindex[i] = -1;
+      
+    } else if (fieldtype[i] == Z) {
+      fieldindex[i] = find_label("z",nwords,labels);
+      zflag = UNSCALED;
+      if (fieldindex[i] < 0) {
+	fieldindex[i] = nwords;
+	s_index = find_label("zs",nwords,labels);
+	u_index = find_label("zu",nwords,labels);
+	su_index = find_label("zsu",nwords,labels);
+	if (s_index >= 0 && s_index < fieldindex[i]) {
+	  fieldindex[i] = s_index;
+	  zflag = SCALED;
+	}
+	if (u_index >= 0 && u_index < fieldindex[i]) {
+	  fieldindex[i] = u_index;
+	  zflag = UNSCALED;
+	}
+	if (su_index >= 0 && su_index < fieldindex[i]) {
+	  fieldindex[i] = su_index;
+	  zflag = SCALED;
+	}
+      }
+      if (fieldindex[i] == nwords) fieldindex[i] = -1;
+      
+    } else if (fieldtype[i] == VX)
+      fieldindex[i] = find_label("vx",nwords,labels);
+    else if (fieldtype[i] == VY)
+      fieldindex[i] = find_label("vy",nwords,labels);
+    else if (fieldtype[i] == VZ)
+      fieldindex[i] = find_label("vz",nwords,labels);
+    else if (fieldtype[i] == IX)
+      fieldindex[i] = find_label("ix",nwords,labels);
+    else if (fieldtype[i] == IY)
+      fieldindex[i] = find_label("iy",nwords,labels);
+    else if (fieldtype[i] == IZ)
+      fieldindex[i] = find_label("iz",nwords,labels);
+  }
+  
+  // error checks
+  
+  if ((triclinic_snap && !triclinic) || 
+      (!triclinic_snap && triclinic))
+    error->one(FLERR,"Read_dump triclinic setting does not match simulation");
+  
+  for (int i = 0; i < nfield; i++)
+    if (fieldindex[i] < 0) 
+      error->one(FLERR,"Read_dump field not found in dump file");
+  
+  // set overall scaling of coordinates
+  // error if x,y,z scaling is not the same
+  
+  scaled = MAX(xflag,yflag);
+  scaled = MAX(zflag,scaled);
+  if ((xflag != UNSET && xflag != scaled) ||
+      (yflag != UNSET && yflag != scaled) ||
+      (zflag != UNSET && zflag != scaled))
+    error->one(FLERR,"Read_dump x,y,z fields do not have consistent scaling");
+  
+  // scaled, triclinic coords require all 3 x,y,z fields to perform unscaling
+  
+  if (scaled == SCALED && triclinic) {
+    int flag = 0;
+    if (xflag != scaled) flag = 1;
+    if (yflag != scaled) flag = 1;
+    if (dimension == 3 && zflag != scaled) flag = 1;
+    if (flag)
+      error->one(FLERR,"All read_dump x,y,z fields must be specified for "
+		 "scaled, triclinic coords");
+  }
+  
+  delete [] labels;
+  
+  // create vector of word ptrs for future parsing of per-atom lines
+  
+  words = new char*[nwords];
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads N atom lines from dump file
+   stores appropriate values in fields array
+------------------------------------------------------------------------- */
+
+void ReadDumpNative::read(int n, double **fields)
+{
+  int i,m;
+  char *eof;
+
+  for (int i = 0; i < n; i++) {
+    eof = fgets(line,MAXLINE,fp);
+    if (eof == NULL) error->one(FLERR,"Unexpected end of dump file");
+
+    // tokenize the line
+
+    words[0] = strtok(line," \t\n\r\f");
+    for (m = 1; m < nwords; m++)
+      words[m] = strtok(NULL," \t\n\r\f");
+
+    // convert selected fields to floats
+
+    for (m = 0; m < nfield; m++)
+      fields[i][m] = atof(words[fieldindex[m]]);
+  }
+}
+
+/* ----------------------------------------------------------------------
+   match label to any of N labels
+   return index of match or -1 if no match
+------------------------------------------------------------------------- */
+
+int ReadDumpNative::find_label(const char *label, int n, char **labels)
+{
+  for (int i = 0; i < n; i++)
+    if (strcmp(label,labels[i]) == 0) return i;
+  return -1;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads N lines from file
+   only last one is saved in line
+------------------------------------------------------------------------- */
+
+void ReadDumpNative::read_lines(int n)
+{
+  char *eof;
+  for (int i = 0; i < n; i++) eof = fgets(line,MAXLINE,fp);
+  if (eof == NULL) error->one(FLERR,"Unexpected end of dump file");
+}

src/read_dump_native.h

@@ -0,0 +1,51 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://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.
+
+   Contributed by Timothy Sirk
+------------------------------------------------------------------------- */
+
+#ifndef LMP_READ_DUMP_NATIVE_H
+#define LMP_READ_DUMP_NATIVE_H
+
+#include "stdio.h"
+#include "pointers.h"
+
+namespace LAMMPS_NS {
+
+class ReadDumpNative : protected Pointers {
+ public:
+  ReadDumpNative(class LAMMPS *);
+  ~ReadDumpNative();
+
+  void init(FILE *);
+  void scan(bigint, int, int *, char **, int, bigint &, double [3][3], int &);
+  void read(int, double **);
+
+private:
+  FILE *fp;                // pointer to file opened by caller
+  char *line;              // line read from dump file
+  int dimension;
+  int triclinic;
+
+  int nwords;              // # of per-atom columns in dump file
+  char **words;            // ptrs to words in parsed per-atom line
+
+  int nfield;              // # of fields to extract for each atom
+  int *fieldindex;         // index into words for each field
+
+  int find_label(const char *, int, char **);
+  void read_lines(int);
+};
+
+}
+
+#endif