diff --git a/doc/src/read_dump.txt b/doc/src/read_dump.txt index fbcc89a453..df357a5018 100644 --- a/doc/src/read_dump.txt +++ b/doc/src/read_dump.txt @@ -22,7 +22,8 @@ field = {x} or {y} or {z} or {vx} or {vy} or {vz} or {q} or {ix} or {iy} or {iz} {ix},{iy},{iz} = image flags in each dimension {fx},{fy},{fz} = force components :pre zero or more keyword/value pairs may be appended :l -keyword = {box} or {replace} or {purge} or {trim} or {add} or {label} or {scaled} or {wrapped} or {format} :l +keyword = {nfile} or {box} or {replace} or {purge} or {trim} or {add} or {label} or {scaled} or {wrapped} or {format} :l + {nfile} value = Nfiles = how many parallel dump files exist {box} value = {yes} or {no} = replace simulation box with dump box {replace} value = {yes} or {no} = overwrite atoms with dump atoms {purge} value = {yes} or {no} = delete all atoms before adding dump atoms @@ -85,9 +86,18 @@ command, after the dump snapshot is read. If the dump filename specified as {file} ends with ".gz", the dump file is read in gzipped format. You cannot (yet) read a dump file -that was written in binary format with a ".bin" suffix, or to multiple -files via the "%" option in the dump file name. See the -"dump"_dump.html command for details. +that was written in binary format with a ".bin" suffix. + +You can read dump files that were written (in parallel) to multiple +files via the "%" wild-card character in the dump file name. If any +specified dump file name contains a "%", they must all contain it. +See the "dump"_dump.html command for details. +The "%" wild-card character is only supported by the {native} format +for dump files, described next. + +If reading parallel dump files, you must also use the {nfile} keyword +to tell LAMMPS how many parallel files exist, via its specified +{Nfiles} value. The format of the dump file is selected through the {format} keyword. If specified, it must be the last keyword used, since all remaining diff --git a/doc/src/rerun.txt b/doc/src/rerun.txt index 9a92a7374f..08a5698615 100644 --- a/doc/src/rerun.txt +++ b/doc/src/rerun.txt @@ -89,12 +89,23 @@ this auxiliary information should be defined in the usual way, e.g. in a data file read in by a "read_data"_read_data.html command, before using the rerun command. +Also note that the frequency of thermodynamic or dump output from the +rerun simulation will depend on settings made in the rerun script, the +same as for output from any LAMMPS simulation. See further info below +as to what that means if the timesteps for snapshots read from dump +files do not match the specified output frequency. + :line If more than one dump file is specified, the dump files are read one after the other. It is assumed that snapshot timesteps will be in ascending order. If a snapshot is encountered that is not in -ascending order, it will cause the rerun command to complete. +ascending order, it will skip the snapshot until it reads one that is. +This allows skipping of a duplicate snapshot (same timestep), +e.g. that appeared at the end of one file and beginning of the next. +However if you specify a series of dump files in an incorrect order +(with respect to the timesteps they contain), you may skip large +numbers of snapshots The {first}, {last}, {every}, {skip} keywords determine which snapshots are read from the dump file(s). Snapshots are skipped until @@ -177,12 +188,12 @@ a timestep it expects to be, LAMMPS will flag an error. The various forms of LAMMPS output, as defined by the "thermo_style"_thermo_style.html, "thermo"_thermo.html, -"dump"_dump.html, and "restart"_restart.html commands occur on -specific timesteps. If successive dump snapshots skip those -timesteps, then no output will be produced. E.g. if you request -thermodynamic output every 100 steps, but the dump file snapshots are -every 1000 steps, then you will only see thermodynamic output every -1000 steps. +"dump"_dump.html, and "restart"_restart.html commands occur with +specified frequency, e.g. every N steps. If the timestep for a dump +snapshot is not a multiple of N, then it will be read and processed, +but no output will be produced. If you want output for every dump +snapshot, you can simply use N=1 for an output frequency, e.g. for +thermodyanmic output or new dump file output. :line diff --git a/src/read_dump.cpp b/src/read_dump.cpp index 326369c17e..551cc63c9e 100644 --- a/src/read_dump.cpp +++ b/src/read_dump.cpp @@ -43,8 +43,7 @@ using namespace LAMMPS_NS; -#define CHUNK 1024 -#define EPSILON 1.0e-6 +#define CHUNK 16384 // also in reader_native.cpp @@ -65,17 +64,21 @@ ReadDump::ReadDump(LAMMPS *lmp) : Pointers(lmp) nfile = 0; files = NULL; + nnew = maxnew = 0; nfield = 0; fieldtype = NULL; fieldlabel = NULL; fields = NULL; + buf = NULL; int n = strlen("native") + 1; readerstyle = new char[n]; strcpy(readerstyle,"native"); - reader = NULL; - fp = NULL; + nreader = 0; + readers = NULL; + nsnapatoms = NULL; + clustercomm = MPI_COMM_NULL; } /* ---------------------------------------------------------------------- */ @@ -90,7 +93,13 @@ ReadDump::~ReadDump() delete [] readerstyle; memory->destroy(fields); - delete reader; + memory->destroy(buf); + + for (int i = 0; i < nreader; i++) delete readers[i]; + delete [] readers; + delete [] nsnapatoms; + + MPI_Comm_free(&clustercomm); } /* ---------------------------------------------------------------------- */ @@ -134,13 +143,20 @@ void ReadDump::command(int narg, char **arg) bigint natoms_prev = atom->natoms; atoms(); - if (me == 0) reader->close_file(); + if (filereader) + for (int i = 0; i < nreader; i++) + readers[i]->close_file(); // print out stats - bigint npurge_all,nreplace_all,ntrim_all,nadd_all; + bigint nsnap_all,npurge_all,nreplace_all,ntrim_all,nadd_all; + + bigint tmp = 0; + if (filereader) + for (int i = 0; i < nreader; i++) + tmp += nsnapatoms[i]; + MPI_Allreduce(&tmp,&nsnap_all,1,MPI_LMP_BIGINT,MPI_SUM,world); - bigint tmp; tmp = npurge; MPI_Allreduce(&tmp,&npurge_all,1,MPI_LMP_BIGINT,MPI_SUM,world); tmp = nreplace; @@ -155,7 +171,7 @@ void ReadDump::command(int narg, char **arg) 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 in snapshot\n",nsnap_all); 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); @@ -164,7 +180,7 @@ void ReadDump::command(int narg, char **arg) } 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 in snapshot\n",nsnap_all); 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); @@ -181,10 +197,22 @@ void ReadDump::store_files(int nstr, char **str) nfile = nstr; files = new char*[nfile]; + // either all or none of files must have '%' wild-card + for (int i = 0; i < nfile; i++) { int n = strlen(str[i]) + 1; files[i] = new char[n]; strcpy(files[i],str[i]); + + if (i == 0) { + if (strchr(files[i],'%')) multiproc = 1; + else multiproc = 0; + } else { + if (multiproc && !strchr(files[i],'%')) + error->all(FLERR,"All read_dump files must be serial or parallel"); + if (!multiproc && strchr(files[i],'%')) + error->all(FLERR,"All read_dump files must be serial or parallel"); + } } } @@ -192,18 +220,47 @@ void ReadDump::store_files(int nstr, char **str) void ReadDump::setup_reader(int narg, char **arg) { - // allocate snapshot field buffer + // setup serial or parallel file reading + // multiproc = 0: only one file to read from, only proc 0 is a reader + // multiproc_nfile >= nprocs: every proc reads one or more files + // multiproc_nfile < nprocs: multiproc_nfile readers, create clusters + // see read_dump.h for explanation of these variables - memory->create(fields,CHUNK,nfield,"read_dump:fields"); + if (multiproc == 0) { + nreader = 1; + firstfile = -1; + MPI_Comm_dup(world,&clustercomm); + } else if (multiproc_nfile >= nprocs) { + firstfile = static_cast ((bigint) me * multiproc_nfile/nprocs); + int lastfile = static_cast ((bigint) (me+1) * multiproc_nfile/nprocs); + nreader = lastfile - firstfile; + MPI_Comm_split(world,me,0,&clustercomm); + } else if (multiproc_nfile < nprocs) { + nreader = 1; + int icluster = static_cast ((bigint) me * multiproc_nfile/nprocs); + firstfile = icluster; + MPI_Comm_split(world,icluster,0,&clustercomm); + } - // create reader class + MPI_Comm_rank(clustercomm,&me_cluster); + MPI_Comm_size(clustercomm,&nprocs_cluster); + if (me_cluster == 0) filereader = 1; + else filereader = 0; + + readers = new Reader*[nreader]; + nsnapatoms = new bigint[nreader]; + + // create Nreader reader classes per reader // match readerstyle to options in style_reader.h if (0) return; // dummy line to enable else-if macro expansion #define READER_CLASS #define ReaderStyle(key,Class) \ - else if (strcmp(readerstyle,#key) == 0) reader = new Class(lmp); + else if (strcmp(readerstyle,#key) == 0) { \ + for (int i = 0; i < nreader; i++) \ + readers[i] = new Class(lmp); \ + } #include "style_reader.h" #undef READER_CLASS @@ -211,9 +268,11 @@ void ReadDump::setup_reader(int narg, char **arg) else error->all(FLERR,"Unknown dump reader style"); - // pass any arguments to reader + // pass any arguments to readers - if (narg > 0) reader->settings(narg,arg); + if (narg > 0 && filereader) + for (int i = 0; i < nreader; i++) + readers[i]->settings(narg,arg); } /* ---------------------------------------------------------------------- @@ -228,6 +287,8 @@ bigint ReadDump::seek(bigint nrequest, int exact) int ifile,eofflag; bigint ntimestep; + // proc 0 finds the timestep in its first reader + if (me == 0) { // exit file loop when dump timestep >= nrequest @@ -235,24 +296,75 @@ bigint ReadDump::seek(bigint nrequest, int exact) for (ifile = 0; ifile < nfile; ifile++) { ntimestep = -1; - reader->open_file(files[ifile]); + if (multiproc) { + char *ptr = strchr(files[ifile],'%'); + char *multiname = new char[strlen(files[ifile]) + 16]; + *ptr = '\0'; + sprintf(multiname,"%s%d%s",files[ifile],0,ptr+1); + *ptr = '%'; + readers[0]->open_file(multiname); + delete [] multiname; + } else readers[0]->open_file(files[ifile]); + while (1) { - eofflag = reader->read_time(ntimestep); + eofflag = readers[0]->read_time(ntimestep); if (eofflag) break; if (ntimestep >= nrequest) break; - reader->skip(); + readers[0]->skip(); } + if (ntimestep >= nrequest) break; - reader->close_file(); + readers[0]->close_file(); } currentfile = ifile; if (ntimestep < nrequest) ntimestep = -1; if (exact && ntimestep != nrequest) ntimestep = -1; - if (ntimestep < 0) reader->close_file(); } + // proc 0 broadcasts timestep and currentfile to all procs + MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world); + MPI_Bcast(¤tfile,1,MPI_INT,0,world); + + // if ntimestep < 0: + // all filereader procs close all their files and return + + if (ntimestep < 0) { + if (filereader) + for (int i = 0; i < nreader; i++) + readers[i]->close_file(); + return ntimestep; + } + + // for multiproc mode: + // all filereader procs search for same ntimestep in currentfile + + if (multiproc && filereader) { + for (int i = 0; i < nreader; i++) { + if (me == 0 && i == 0) continue; // proc 0, reader 0 already found it + char *ptr = strchr(files[currentfile],'%'); + char *multiname = new char[strlen(files[currentfile]) + 16]; + *ptr = '\0'; + sprintf(multiname,"%s%d%s",files[currentfile],firstfile+i,ptr+1); + *ptr = '%'; + readers[i]->open_file(multiname); + delete [] multiname; + + bigint step; + while (1) { + eofflag = readers[i]->read_time(step); + if (eofflag) break; + if (step == ntimestep) break; + readers[i]->skip(); + } + + if (eofflag) + error->one(FLERR,"Read dump parallel files " + "do not all have same timestep"); + } + } + return ntimestep; } @@ -270,6 +382,8 @@ bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip) int ifile,eofflag; bigint ntimestep; + // proc 0 finds the timestep in its first reader + if (me == 0) { // exit file loop when dump timestep matches all criteria @@ -279,19 +393,34 @@ bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip) for (ifile = currentfile; ifile < nfile; ifile++) { ntimestep = -1; - if (ifile != currentfile) reader->open_file(files[ifile]); + if (ifile != currentfile) { + if (multiproc) { + char *ptr = strchr(files[ifile],'%'); + char *multiname = new char[strlen(files[ifile]) + 16]; + *ptr = '\0'; + sprintf(multiname,"%s%d%s",files[ifile],0,ptr+1); + *ptr = '%'; + readers[0]->open_file(multiname); + delete [] multiname; + } else readers[0]->open_file(files[ifile]); + } + while (1) { - eofflag = reader->read_time(ntimestep); + eofflag = readers[0]->read_time(ntimestep); + if (eofflag) break; + if (ntimestep > nlast) break; + if (ntimestep <= ncurrent) { + readers[0]->skip(); + continue; + } if (iskip == nskip) iskip = 0; iskip++; - if (eofflag) break; - if (ntimestep <= ncurrent) break; - if (ntimestep > nlast) break; - if (nevery && ntimestep % nevery) reader->skip(); - else if (iskip < nskip) reader->skip(); + if (nevery && ntimestep % nevery) readers[0]->skip(); + else if (iskip < nskip) readers[0]->skip(); else break; } - if (eofflag) reader->close_file(); + + if (eofflag) readers[0]->close_file(); else break; } @@ -299,10 +428,50 @@ bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip) if (eofflag) ntimestep = -1; if (ntimestep <= ncurrent) ntimestep = -1; if (ntimestep > nlast) ntimestep = -1; - if (ntimestep < 0) reader->close_file(); } + // proc 0 broadcasts timestep and currentfile to all procs + MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world); + MPI_Bcast(¤tfile,1,MPI_INT,0,world); + + // if ntimestep < 0: + // all filereader procs close all their files and return + + if (ntimestep < 0) { + for (int i = 0; i < nreader; i++) + readers[i]->close_file(); + return ntimestep; + } + + // for multiproc mode: + // all filereader procs search for same ntimestep in currentfile + + if (multiproc && filereader) { + for (int i = 0; i < nreader; i++) { + if (me == 0 && i == 0) continue; + char *ptr = strchr(files[currentfile],'%'); + char *multiname = new char[strlen(files[currentfile]) + 16]; + *ptr = '\0'; + sprintf(multiname,"%s%d%s",files[currentfile],firstfile+i,ptr+1); + *ptr = '%'; + readers[i]->open_file(multiname); + delete [] multiname; + + bigint step; + while (1) { + eofflag = readers[i]->read_time(step); + if (eofflag) break; + if (step == ntimestep) break; + readers[i]->skip(); + } + + if (eofflag) + error->one(FLERR,"Read dump parallel files " + "do not all have same timestep"); + } + } + return ntimestep; } @@ -316,15 +485,17 @@ void ReadDump::header(int fieldinfo) int triclinic_snap; int fieldflag,xflag,yflag,zflag; - if (me == 0) - nsnapatoms = reader->read_header(box,triclinic_snap, - fieldinfo,nfield,fieldtype,fieldlabel, - scaleflag,wrapflag,fieldflag, - xflag,yflag,zflag); + if (filereader) { + for (int i = 0; i < nreader; i++) + nsnapatoms[i] = readers[i]->read_header(box,triclinic_snap,fieldinfo, + nfield,fieldtype,fieldlabel, + scaleflag,wrapflag,fieldflag, + xflag,yflag,zflag); + } - MPI_Bcast(&nsnapatoms,1,MPI_LMP_BIGINT,0,world); - MPI_Bcast(&triclinic_snap,1,MPI_INT,0,world); - MPI_Bcast(&box[0][0],9,MPI_DOUBLE,0,world); + MPI_Bcast(nsnapatoms,nreader,MPI_LMP_BIGINT,0,clustercomm); + MPI_Bcast(&triclinic_snap,1,MPI_INT,0,clustercomm); + MPI_Bcast(&box[0][0],9,MPI_DOUBLE,0,clustercomm); // local copy of snapshot box parameters // used in xfield,yfield,zfield when converting dump atom to absolute coords @@ -358,10 +529,10 @@ void ReadDump::header(int fieldinfo) if (!fieldinfo) return; - MPI_Bcast(&fieldflag,1,MPI_INT,0,world); - MPI_Bcast(&xflag,1,MPI_INT,0,world); - MPI_Bcast(&yflag,1,MPI_INT,0,world); - MPI_Bcast(&zflag,1,MPI_INT,0,world); + MPI_Bcast(&fieldflag,1,MPI_INT,0,clustercomm); + MPI_Bcast(&xflag,1,MPI_INT,0,clustercomm); + MPI_Bcast(&yflag,1,MPI_INT,0,clustercomm); + MPI_Bcast(&zflag,1,MPI_INT,0,clustercomm); // error check on current vs new box and fields // triclinic_snap < 0 means no box info in file @@ -430,7 +601,9 @@ void ReadDump::header(int fieldinfo) } } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + read and process one snapshot of atoms +------------------------------------------------------------------------- */ void ReadDump::atoms() { @@ -447,8 +620,23 @@ void ReadDump::atoms() atom->natoms = 0; } - // to match existing atoms to dump atoms: + // read all the snapshot atoms into fields + // each proc will own an arbitrary subset of atoms + + read_atoms(); + + // migrate old owned atoms to new procs based on atom IDs + // not necessary if purged all old atoms or if only 1 proc + + if (!purgeflag && nprocs > 1) migrate_old_atoms(); + + // migrate new snapshot atoms to same new procs based on atom IDs + // not necessary if purged all old atoms or if only 1 proc + + if (!purgeflag && nprocs > 1) migrate_new_atoms(); + // must build map if not a molecular system + // this will be needed to match new atoms to old atoms int mapflag = 0; if (atom->map_style == 0) { @@ -457,57 +645,14 @@ void ReadDump::atoms() atom->map_set(); } - // uflag[i] = 1 for each owned atom appearing in dump - // ucflag = similar flag for each chunk atom, used in process_atoms() + // each proc now owns both old and new info for same subset of atoms + // update each local atom with new info - 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"); + process_atoms(); - // read, broadcast, and process atoms from snapshot in chunks + // check that atom IDs are valid - addproc = -1; - - int nchunk; - bigint nread = 0; - while (nread < nsnapatoms) { - nchunk = MIN(nsnapatoms-nread,CHUNK); - if (me == 0) reader->read_atoms(nchunk,nfield,fields); - MPI_Bcast(&fields[0][0],nchunk*nfield,MPI_DOUBLE,0,world); - process_atoms(nchunk); - nread += nchunk; - } - - // if addflag = YESADD or KEEPADD, update total atom count - - if (addflag == YESADD || addflag == KEEPADD) { - bigint nblocal = atom->nlocal; - MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); - } - - // if addflag = YESADD, assign IDs to new snapshot atoms - - if (addflag == YESADD) { - if (atom->natoms < 0 || atom->natoms >= MAXBIGINT) - error->all(FLERR,"Too many total atoms"); - if (atom->tag_enable) atom->tag_extend(); - } - - // if trimflag set, delete atoms not replaced by snapshot atoms - - if (trimflag) { - delete_atoms(); - bigint nblocal = atom->nlocal; - MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); - } - - // can now delete uflag arrays - - memory->destroy(uflag); - memory->destroy(ucflag); - memory->destroy(ucflag_all); + atom->tag_check(); // delete atom map if created it above // else reinitialize map for current atoms @@ -548,209 +693,152 @@ void ReadDump::atoms() domain->set_local_box(); } - // 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 + // migrate atoms to their new owing proc, based on atom coords - double **x = atom->x; - imageint *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(1); - delete irregular; - if (triclinic) domain->lamda2x(atom->nlocal); - - // check that atom IDs are valid - - atom->tag_check(); + migrate_atoms_by_coords(); } /* ---------------------------------------------------------------------- - process arg list for dump file fields and optional keywords + read all the snapshot atoms into fields + done in different ways for multiproc no/yes and # of procs < or >= nprocs + nnew = # of snapshot atoms this proc stores ------------------------------------------------------------------------- */ -int ReadDump::fields_and_keywords(int narg, char **arg) +void ReadDump::read_atoms() { - // per-field vectors, leave space for ID and TYPE + int count,nread,nsend,nrecv,otherproc; + bigint nsnap,ntotal,ofirst,olast,rfirst,rlast,lo,hi; + MPI_Request request; + MPI_Status status; - fieldtype = new int[narg+2]; - fieldlabel = new char*[narg+2]; + // one reader per cluster of procs + // each reading proc reads one file and splits data across cluster + // cluster can be all procs or a subset - // add id and type fields as needed - // scan ahead to see if "add yes/keep" keyword/value is used - // requires extra "type" field from from dump file + if (!multiproc || multiproc_nfile < nprocs) { + nsnap = nsnapatoms[0]; - int iarg; - for (iarg = 0; iarg < narg; iarg++) - if (strcmp(arg[iarg],"add") == 0) - if (iarg < narg-1 && (strcmp(arg[iarg+1],"yes") == 0 || - strcmp(arg[iarg+1],"keep") == 0)) break; + if (filereader) { + if (!buf) memory->create(buf,CHUNK,nfield,"read_dump:buf"); - nfield = 0; - fieldtype[nfield++] = ID; - if (iarg < narg) fieldtype[nfield++] = TYPE; + otherproc = 0; + ofirst = (bigint) otherproc * nsnap/nprocs_cluster; + olast = (bigint) (otherproc+1) * nsnap/nprocs_cluster; + if (olast-ofirst > MAXSMALLINT) + error->one(FLERR,"Read dump snapshot is too large for a proc"); + nnew = static_cast (olast - ofirst); - // parse fields + if (nnew > maxnew || maxnew == 0) { + memory->destroy(fields); + maxnew = MAX(nnew,1); // avoid NULL ptr + memory->create(fields,maxnew,nfield,"read_dump:fields"); + } - iarg = 0; - while (iarg < narg) { - int type = whichtype(arg[iarg]); - if (type < 0) break; - if (type == Q && !atom->q_flag) - error->all(FLERR,"Read dump of atom property that isn't allocated"); - fieldtype[nfield++] = type; - iarg++; - } + ntotal = 0; + while (ntotal < nsnap) { + nread = MIN(CHUNK,nsnap-ntotal); + readers[0]->read_atoms(nread,nfield,buf); + rfirst = ntotal; + rlast = ntotal + nread; - // check for no fields + nsend = 0; + while (nsend < nread) { + lo = MAX(ofirst,rfirst); + hi = MIN(olast,rlast); + if (otherproc) // send to otherproc or copy to self + MPI_Send(&buf[nsend][0],(hi-lo)*nfield,MPI_DOUBLE, + otherproc,0,clustercomm); + else + memcpy(&fields[rfirst][0],&buf[nsend][0], + (hi-lo)*nfield*sizeof(double)); + nsend += hi-lo; + if (hi == olast) { + otherproc++; + ofirst = (bigint) otherproc * nsnap/nprocs_cluster; + olast = (bigint) (otherproc+1) * nsnap/nprocs_cluster; + } + } - if (fieldtype[nfield-1] == ID || fieldtype[nfield-1] == TYPE) - error->all(FLERR,"Illegal read_dump command"); + ntotal += nread; + } - if (dimension == 2) { - for (int i = 0; i < nfield; i++) - if (fieldtype[i] == Z || fieldtype[i] == VZ || - fieldtype[i] == IZ || fieldtype[i] == FZ) - error->all(FLERR,"Illegal read_dump command"); - } + } else { + ofirst = (bigint) me_cluster * nsnap/nprocs_cluster; + olast = (bigint) (me_cluster+1) * nsnap/nprocs_cluster; + if (olast-ofirst > MAXSMALLINT) + error->one(FLERR,"Read dump snapshot is too large for a proc"); + nnew = static_cast (olast - ofirst); + if (nnew > maxnew || maxnew == 0) { + memory->destroy(fields); + maxnew = MAX(nnew,1); // avoid NULL ptr + memory->create(fields,maxnew,nfield,"read_dump:fields"); + } - 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"); + nrecv = 0; + while (nrecv < nnew) { + MPI_Irecv(&fields[nrecv][0],(nnew-nrecv)*nfield,MPI_DOUBLE,0,0, + clustercomm,&request); + MPI_Wait(&request,&status); + MPI_Get_count(&status,MPI_DOUBLE,&count); + nrecv += count/nfield; + } + } - // parse optional args + // every proc is a filereader, reads one or more files + // each proc keeps all data it reads, no communication required - boxflag = 1; - replaceflag = 1; - purgeflag = 0; - trimflag = 0; - addflag = NOADD; - for (int i = 0; i < nfield; i++) fieldlabel[i] = NULL; - scaleflag = 0; - wrapflag = 1; + } else if (multiproc_nfile >= nprocs) { + bigint sum = 0; + for (int i = 0; i < nreader; i++) + sum += nsnapatoms[i]; + if (sum > MAXSMALLINT) + error->one(FLERR,"Read dump snapshot is too large for a proc"); + nnew = static_cast (sum); + if (nnew > maxnew || maxnew == 0) { + memory->destroy(fields); + maxnew = MAX(nnew,1); // avoid NULL ptr + memory->create(fields,maxnew,nfield,"read_dump:fields"); + } - 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 = YESADD; - else if (strcmp(arg[iarg+1],"no") == 0) addflag = NOADD; - else if (strcmp(arg[iarg+1],"keep") == 0) addflag = KEEPADD; - 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 type = whichtype(arg[iarg+1]); - int i; - for (i = 0; i < nfield; i++) - if (type == fieldtype[i]) break; - if (i == nfield) error->all(FLERR,"Illegal read_dump command"); - int n = strlen(arg[iarg+2]) + 1; - fieldlabel[i] = new char[n]; - strcpy(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) scaleflag = 1; - else if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0; - else error->all(FLERR,"Illegal read_dump command"); - iarg += 2; - } else if (strcmp(arg[iarg],"wrapped") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command"); - if (strcmp(arg[iarg+1],"yes") == 0) wrapflag = 1; - else if (strcmp(arg[iarg+1],"no") == 0) wrapflag = 0; - 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"); - delete [] readerstyle; - int n = strlen(arg[iarg+1]) + 1; - readerstyle = new char[n]; - strcpy(readerstyle,arg[iarg+1]); - iarg += 2; - break; - } else error->all(FLERR,"Illegal read_dump command"); - } - - if (purgeflag && (replaceflag || trimflag)) - error->all(FLERR,"If read_dump purges it cannot replace or trim"); - if (addflag == KEEPADD && atom->tag_enable == 0) - error->all(FLERR,"Read_dump cannot use 'add keep' without atom IDs"); - - return narg-iarg; + nnew = 0; + for (int i = 0; i < nreader; i++) { + nsnap = nsnapatoms[i]; + ntotal = 0; + while (ntotal < nsnap) { + nread = MIN(CHUNK,nsnap-ntotal); + readers[i]->read_atoms(nread,nfield,&fields[nnew+ntotal]); + ntotal += nread; + } + nnew += nsnap; + } + } } /* ---------------------------------------------------------------------- - check if str is a field argument - if yes, return index of which - if not, return -1 -------------------------------------------------------------------------- */ - -int ReadDump::whichtype(char *str) -{ - int type = -1; - if (strcmp(str,"id") == 0) type = ID; - else if (strcmp(str,"type") == 0) type = TYPE; - else if (strcmp(str,"x") == 0) type = X; - else if (strcmp(str,"y") == 0) type = Y; - else if (strcmp(str,"z") == 0) type = Z; - else if (strcmp(str,"vx") == 0) type = VX; - else if (strcmp(str,"vy") == 0) type = VY; - else if (strcmp(str,"vz") == 0) type = VZ; - else if (strcmp(str,"q") == 0) type = Q; - else if (strcmp(str,"ix") == 0) type = IX; - else if (strcmp(str,"iy") == 0) type = IY; - else if (strcmp(str,"iz") == 0) type = IZ; - else if (strcmp(str,"fx") == 0) type = FX; - else if (strcmp(str,"fy") == 0) type = FY; - else if (strcmp(str,"fz") == 0) type = FZ; - return type; -} - -/* ---------------------------------------------------------------------- - process each of N atoms in chunk read from dump file + update info for each old atom I own based on snapshot info 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 + if in add mode and atom ID does not match any old atom, + create new atom with dump file field values ------------------------------------------------------------------------- */ -void ReadDump::process_atoms(int n) +void ReadDump::process_atoms() { int i,m,ifield,itype; int xbox,ybox,zbox; tagint mtag; + int *updateflag,*newflag; + + // updateflag[i] = flag for old atoms, 1 if updated, else 0 + // newflag[i] = flag for new atoms, 0 if used to update old atom, else 1 + + int nlocal = atom->nlocal; + memory->create(updateflag,nlocal,"read_dump:updateflag"); + for (int i = 0; i < nlocal; i++) updateflag[i] = 0; + memory->create(newflag,nnew,"read_dump:newflag"); + for (int i = 0; i < nnew; i++) newflag[i] = 1; + + // loop over new atoms double **x = atom->x; double **v = atom->v; @@ -758,11 +846,9 @@ void ReadDump::process_atoms(int n) double **f = atom->f; tagint *tag = atom->tag; imageint *image = atom->image; - int nlocal = atom->nlocal; tagint map_tag_max = atom->map_tag_max; - for (i = 0; i < n; i++) { - ucflag[i] = 0; + for (i = 0; i < nnew; i++) { // check if new atom matches one I own // setting m = -1 forces new atom not to match @@ -774,8 +860,8 @@ void ReadDump::process_atoms(int n) else m = -1; if (m < 0 || m >= nlocal) continue; - ucflag[i] = 1; - uflag[m] = 1; + updateflag[m] = 1; + newflag[i] = 0; if (replaceflag) { nreplace++; @@ -843,12 +929,39 @@ void ReadDump::process_atoms(int n) } } - // 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 - // check that dump file snapshot has atom type field + // if trimflag set, delete atoms not updated by snapshot atoms - if (addflag == NOADD) return; + if (trimflag) { + AtomVec *avec = atom->avec; + + int i = 0; + while (i < nlocal) { + if (!updateflag[i]) { + avec->copy(nlocal-1,i,1); + updateflag[i] = updateflag[nlocal-1]; + nlocal--; + ntrim++; + } else i++; + } + + atom->nlocal = nlocal; + bigint nblocal = atom->nlocal; + MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); + } + + // done if cannot add new atoms + + if (addflag == NOADD) { + memory->destroy(updateflag); + memory->destroy(newflag); + return; + } + + // ---------------------------------------------------- + // create new atoms for dump file atoms with ID that matches no old atom + // ---------------------------------------------------- + + // first check that dump file snapshot has atom type field int tflag = 0; for (ifield = 0; ifield < nfield; ifield++) @@ -856,19 +969,11 @@ void ReadDump::process_atoms(int n) if (!tflag) error->all(FLERR,"Cannot add atoms if dump file does not store atom type"); - MPI_Allreduce(ucflag,ucflag_all,n,MPI_INT,MPI_SUM,world); - int nlocal_previous = atom->nlocal; double one[3]; - 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; + for (i = 0; i < nnew; i++) { + if (!newflag[i]) continue; // create type and coord fields from dump file // coord = 0.0 unless corresponding dump file field was specified @@ -937,7 +1042,7 @@ void ReadDump::process_atoms(int n) break; } - // replace image flag in case changed by ix,iy,iz fields + // reset image flag in case changed by ix,iy,iz fields image[m] = ((imageint) (xbox + IMGMAX) & IMGMASK) | (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | @@ -945,31 +1050,287 @@ void ReadDump::process_atoms(int n) } } + // if addflag = YESADD or KEEPADD, update total atom count + + if (addflag == YESADD || addflag == KEEPADD) { + bigint nblocal = atom->nlocal; + MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); + } + + // if addflag = YESADD, + // assign consistent IDs to new snapshot atoms across all procs + + if (addflag == YESADD) { + if (atom->natoms < 0 || atom->natoms >= MAXBIGINT) + error->all(FLERR,"Too many total atoms"); + if (atom->tag_enable) atom->tag_extend(); + } + // init per-atom fix/compute/variable values for created atoms atom->data_fix_compute_variable(nlocal_previous,atom->nlocal); + + // free allocated vectors + + memory->destroy(updateflag); + memory->destroy(newflag); } /* ---------------------------------------------------------------------- - delete atoms not flagged as replaced by dump atoms + migrate old atoms to new procs based on atom IDs + use migrate_atoms() with explicit processor assignments ------------------------------------------------------------------------- */ -void ReadDump::delete_atoms() +void ReadDump::migrate_old_atoms() { - AtomVec *avec = atom->avec; + tagint *tag = atom->tag; 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++; + int *procassign; + memory->create(procassign,nlocal,"read_dump:procassign"); + for (int i = 0; i < nlocal; i++) + procassign[i] = tag[i] % nprocs; + + Irregular *irregular = new Irregular(lmp); + irregular->migrate_atoms(1,1,procassign); + delete irregular; + + memory->destroy(procassign); +} + +/* ---------------------------------------------------------------------- + migrate new atoms to same new procs based on atom IDs +------------------------------------------------------------------------- */ + +void ReadDump::migrate_new_atoms() +{ + tagint mtag; + int *procassign; + double **newfields; + + memory->create(procassign,nnew,"read_dump:procassign"); + for (int i = 0; i < nnew; i++) { + mtag = static_cast (fields[i][0]); + procassign[i] = mtag % nprocs; } - atom->nlocal = nlocal; + Irregular *irregular = new Irregular(lmp); + int nrecv = irregular->create_data(nnew,procassign,1); + int newmaxnew = MAX(nrecv,maxnew); + newmaxnew = MAX(newmaxnew,1); // avoid NULL ptr + memory->create(newfields,newmaxnew,nfield,"read_dump:newfields"); + irregular->exchange_data((char *) &fields[0][0],nfield*sizeof(double), + (char *) &newfields[0][0]); + irregular->destroy_data(); + delete irregular; + + memory->destroy(fields); + memory->destroy(procassign); + + // point fields at newfields + + fields = newfields; + maxnew = newmaxnew; + nnew = nrecv; +} + +/* ---------------------------------------------------------------------- + migrate final atoms to new procs based on atom coords + use migrate_atoms() with implicit processor assignments based on atom coords + 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 +------------------------------------------------------------------------- */ + +void ReadDump::migrate_atoms_by_coords() +{ + double **x = atom->x; + imageint *image = atom->image; + int 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(1); + delete irregular; + if (triclinic) domain->lamda2x(atom->nlocal); +} + +/* ---------------------------------------------------------------------- + process arg list for dump file fields and optional keywords +------------------------------------------------------------------------- */ + +int ReadDump::fields_and_keywords(int narg, char **arg) +{ + // per-field vectors, leave space for ID and TYPE + + fieldtype = new int[narg+2]; + fieldlabel = new char*[narg+2]; + + // add id and type fields as needed + // scan ahead to see if "add yes/keep" keyword/value is used + // requires extra "type" field from from dump file + + int iarg; + for (iarg = 0; iarg < narg; iarg++) + if (strcmp(arg[iarg],"add") == 0) + if (iarg < narg-1 && (strcmp(arg[iarg+1],"yes") == 0 || + strcmp(arg[iarg+1],"keep") == 0)) break; + + nfield = 0; + fieldtype[nfield++] = ID; + if (iarg < narg) fieldtype[nfield++] = TYPE; + + // parse fields + + iarg = 0; + while (iarg < narg) { + int type = whichtype(arg[iarg]); + if (type < 0) break; + if (type == Q && !atom->q_flag) + error->all(FLERR,"Read dump of atom property that isn't allocated"); + fieldtype[nfield++] = type; + iarg++; + } + + // check for no fields + + 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 || fieldtype[i] == FZ) + 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 + + multiproc_nfile = 0; + boxflag = 1; + replaceflag = 1; + purgeflag = 0; + trimflag = 0; + addflag = NOADD; + for (int i = 0; i < nfield; i++) fieldlabel[i] = NULL; + scaleflag = 0; + wrapflag = 1; + + while (iarg < narg) { + if (strcmp(arg[iarg],"nfile") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command"); + multiproc_nfile = force->inumeric(FLERR,arg[iarg+1]); + iarg += 2; + } else 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 = YESADD; + else if (strcmp(arg[iarg+1],"no") == 0) addflag = NOADD; + else if (strcmp(arg[iarg+1],"keep") == 0) addflag = KEEPADD; + 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 type = whichtype(arg[iarg+1]); + int i; + for (i = 0; i < nfield; i++) + if (type == fieldtype[i]) break; + if (i == nfield) error->all(FLERR,"Illegal read_dump command"); + int n = strlen(arg[iarg+2]) + 1; + fieldlabel[i] = new char[n]; + strcpy(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) scaleflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0; + else error->all(FLERR,"Illegal read_dump command"); + iarg += 2; + } else if (strcmp(arg[iarg],"wrapped") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command"); + if (strcmp(arg[iarg+1],"yes") == 0) wrapflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) wrapflag = 0; + 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"); + delete [] readerstyle; + int n = strlen(arg[iarg+1]) + 1; + readerstyle = new char[n]; + strcpy(readerstyle,arg[iarg+1]); + iarg += 2; + break; + } else error->all(FLERR,"Illegal read_dump command"); + } + + if (multiproc == 0 && multiproc_nfile) + error->all(FLERR,"Dump file is not a multi-proc file"); + if (multiproc && multiproc_nfile == 0) + error->all(FLERR,"Dump file is a multi-proc file"); + + if (purgeflag && (replaceflag || trimflag)) + error->all(FLERR,"If read_dump purges it cannot replace or trim"); + if (addflag == KEEPADD && atom->tag_enable == 0) + error->all(FLERR,"Read_dump cannot use 'add keep' without atom IDs"); + + return narg-iarg; +} + +/* ---------------------------------------------------------------------- + check if str is a field argument + if yes, return index of which + if not, return -1 +------------------------------------------------------------------------- */ + +int ReadDump::whichtype(char *str) +{ + int type = -1; + if (strcmp(str,"id") == 0) type = ID; + else if (strcmp(str,"type") == 0) type = TYPE; + else if (strcmp(str,"x") == 0) type = X; + else if (strcmp(str,"y") == 0) type = Y; + else if (strcmp(str,"z") == 0) type = Z; + else if (strcmp(str,"vx") == 0) type = VX; + else if (strcmp(str,"vy") == 0) type = VY; + else if (strcmp(str,"vz") == 0) type = VZ; + else if (strcmp(str,"q") == 0) type = Q; + else if (strcmp(str,"ix") == 0) type = IX; + else if (strcmp(str,"iy") == 0) type = IY; + else if (strcmp(str,"iz") == 0) type = IZ; + else if (strcmp(str,"fx") == 0) type = FX; + else if (strcmp(str,"fy") == 0) type = FY; + else if (strcmp(str,"fz") == 0) type = FZ; + return type; } /* ---------------------------------------------------------------------- diff --git a/src/read_dump.h b/src/read_dump.h index eb70b99808..73cb0a8995 100644 --- a/src/read_dump.h +++ b/src/read_dump.h @@ -22,6 +22,7 @@ CommandStyle(read_dump,ReadDump) #ifndef LMP_READ_DUMP_H #define LMP_READ_DUMP_H +#include #include #include "pointers.h" @@ -43,15 +44,26 @@ class ReadDump : protected Pointers { private: int me,nprocs; - FILE *fp; - int dimension; + char **files; // list of input dump files to process + int nfile; // # of dump files to process (each may be parallel) + int currentfile; // current open file (0 to nfile-1) + + MPI_Comm clustercomm; // comm for proc cluster that reads/shares a file + int me_cluster,nprocs_cluster; // proc ID and count for my read cluster + + int multiproc; // 0 = each dump file is a single file + // 1 = each dump file is parallel (multiple files) + int multiproc_nfile; // number of parallel files in one dump file + + int nreader; // # of parallel dump files read by my cluster + int firstfile; // index of 1st dump file my cluster reads + // (0 to multiproc_nfile-1) + int filereader; // 1 if this proc reads from a dump file(s) + + int dimension; // same as in Domain int triclinic; - int nfile; // # of dump files to process - char **files; // list of file names - int currentfile; // currently open file - int boxflag; // overwrite simulation with dump file box params int replaceflag,addflag; // flags for processing dump snapshot atoms int trimflag,purgeflag; @@ -59,10 +71,13 @@ private: int wrapflag; // user 0/1 if dump file coords are unwrapped/wrapped char *readerstyle; // style of dump files to read + int nnew; // # of dump file atoms this proc owns 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 + int maxnew; // allocation size of fields array + double **buf; // read buffer int scaled; // 0/1 if dump file coords are unscaled/scaled int wrapped; // 0/1 if dump file coords are unwrapped/wrapped @@ -71,20 +86,26 @@ private: double xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz; // dump snapshot box params double xprd,yprd,zprd; - bigint nsnapatoms; // # of atoms in dump file shapshot + bigint *nsnapatoms; // # of atoms in one snapshot from + // one (parallel) dump file + // nreader-length vector b/c a reader proc + // may read from multiple parallel dump files 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 Reader **readers; // class that reads a dump file + // nreader-length list of readers if proc reads + // from multiple parallel dump files - class Reader *reader; // class that reads dump file + void read_atoms(); + void process_atoms(); + void migrate_old_atoms(); + void migrate_new_atoms(); + void migrate_atoms_by_coords(); + void setup_multiproc(); int whichtype(char *); - void process_atoms(int); - void delete_atoms(); double xfield(int, int); double yfield(int, int); diff --git a/src/read_dump2.cpp b/src/read_dump2.cpp new file mode 100644 index 0000000000..37bc3e3589 --- /dev/null +++ b/src/read_dump2.cpp @@ -0,0 +1,1011 @@ +/* ---------------------------------------------------------------------- + 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 (ARL) +------------------------------------------------------------------------- */ + +// lmptype.h must be first b/c this file uses MAXBIGINT and includes mpi.h +// due to OpenMPI bug which sets INT64_MAX via its mpi.h +// before lmptype.h can set flags to insure it is done correctly + +#include "lmptype.h" +#include +#include +#include +#include "read_dump2.h" +#include "reader.h" +#include "style_reader.h" +#include "atom.h" +#include "atom_vec.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "compute.h" +#include "domain.h" +#include "comm.h" +#include "force.h" +#include "irregular.h" +#include "input.h" +#include "variable.h" +#include "error.h" +#include "memory.h" + +using namespace LAMMPS_NS; + +#define CHUNK 16384 + +// also in reader_native.cpp + +enum{ID,TYPE,X,Y,Z,VX,VY,VZ,Q,IX,IY,IZ,FX,FY,FZ}; +enum{UNSET,NOSCALE_NOWRAP,NOSCALE_WRAP,SCALE_NOWRAP,SCALE_WRAP}; +enum{NOADD,YESADD,KEEPADD}; + +/* ---------------------------------------------------------------------- */ + +ReadDump2::ReadDump2(LAMMPS *lmp) : Pointers(lmp) +{ + MPI_Comm_rank(world,&me); + MPI_Comm_size(world,&nprocs); + + dimension = domain->dimension; + triclinic = domain->triclinic; + + nfile = 0; + files = NULL; + + nfield = 0; + fieldtype = NULL; + fieldlabel = NULL; + fields = NULL; + + int n = strlen("native") + 1; + readerstyle = new char[n]; + strcpy(readerstyle,"native"); + + reader = NULL; + fp = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ReadDump2::~ReadDump2() +{ + for (int i = 0; i < nfile; i++) delete [] files[i]; + delete [] files; + for (int i = 0; i < nfield; i++) delete [] fieldlabel[i]; + delete [] fieldlabel; + delete [] fieldtype; + delete [] readerstyle; + + memory->destroy(fields); + delete reader; +} + +/* ---------------------------------------------------------------------- */ + +void ReadDump2::command(int narg, char **arg) +{ + if (domain->box_exist == 0) + error->all(FLERR,"Read_dump command before simulation box is defined"); + + if (narg < 2) error->all(FLERR,"Illegal read_dump command"); + + store_files(1,&arg[0]); + bigint nstep = force->bnumeric(FLERR,arg[1]); + + int nremain = narg - 2; + if (nremain) nremain = fields_and_keywords(nremain,&arg[narg-nremain]); + else nremain = fields_and_keywords(0,NULL); + if (nremain) setup_reader(nremain,&arg[narg-nremain]); + else setup_reader(0,NULL); + + // find the snapshot and read/bcast/process header info + + if (me == 0 && screen) fprintf(screen,"Scanning dump file ...\n"); + + bigint ntimestep = seek(nstep,1); + if (ntimestep < 0) + error->all(FLERR,"Dump file does not contain requested snapshot"); + header(1); + + // reset timestep to nstep + + update->reset_timestep(nstep); + + // counters + + // read in the snapshot and reset system + + if (me == 0 && screen) + fprintf(screen,"Reading snapshot from dump file ...\n"); + + bigint natoms_prev = atom->natoms; + atoms(); + + if (me == 0) reader->close_file(); + + // 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); + + domain->print_box(" "); + + 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); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void ReadDump2::store_files(int nstr, char **str) +{ + nfile = nstr; + files = new char*[nfile]; + + for (int i = 0; i < nfile; i++) { + int n = strlen(str[i]) + 1; + files[i] = new char[n]; + strcpy(files[i],str[i]); + } +} + +/* ---------------------------------------------------------------------- */ + +void ReadDump2::setup_reader(int narg, char **arg) +{ + // allocate snapshot field buffer + + memory->create(fields,CHUNK,nfield,"read_dump:fields"); + + // create reader class + // match readerstyle to options in style_reader.h + + if (0) return; // dummy line to enable else-if macro expansion + +#define READER_CLASS +#define ReaderStyle(key,Class) \ + else if (strcmp(readerstyle,#key) == 0) reader = new Class(lmp); +#include "style_reader.h" +#undef READER_CLASS + + // unrecognized style + + else error->all(FLERR,"Unknown dump reader style"); + + // pass any arguments to reader + + if (narg > 0) reader->settings(narg,arg); +} + +/* ---------------------------------------------------------------------- + seek Nrequest timestep in one or more dump files + if exact = 1, must find exactly Nrequest + if exact = 0, find first step >= Nrequest + return matching ntimestep or -1 if did not find a match +------------------------------------------------------------------------- */ + +bigint ReadDump2::seek(bigint nrequest, int exact) +{ + int ifile,eofflag; + bigint ntimestep; + + if (me == 0) { + + // exit file loop when dump timestep >= nrequest + // or files exhausted + + for (ifile = 0; ifile < nfile; ifile++) { + ntimestep = -1; + reader->open_file(files[ifile]); + while (1) { + eofflag = reader->read_time(ntimestep); + if (eofflag) break; + if (ntimestep >= nrequest) break; + reader->skip(); + } + if (ntimestep >= nrequest) break; + reader->close_file(); + } + + currentfile = ifile; + if (ntimestep < nrequest) ntimestep = -1; + if (exact && ntimestep != nrequest) ntimestep = -1; + if (ntimestep < 0) reader->close_file(); + } + + MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world); + return ntimestep; +} + +/* ---------------------------------------------------------------------- + find next matching snapshot in one or more dump files + Ncurrent = current timestep from last snapshot + Nlast = match no timestep bigger than Nlast + Nevery = only match timesteps that are a multiple of Nevery + Nskip = skip every this many timesteps + return matching ntimestep or -1 if did not find a match +------------------------------------------------------------------------- */ + +bigint ReadDump2::next(bigint ncurrent, bigint nlast, int nevery, int nskip) +{ + int ifile,eofflag; + bigint ntimestep; + + if (me == 0) { + + // exit file loop when dump timestep matches all criteria + // or files exhausted + + int iskip = 0; + + for (ifile = currentfile; ifile < nfile; ifile++) { + ntimestep = -1; + if (ifile != currentfile) reader->open_file(files[ifile]); + while (1) { + eofflag = reader->read_time(ntimestep); + + // new code logic to match new parallel read_dump + if (eofflag) break; + if (ntimestep > nlast) break; + if (ntimestep <= ncurrent) { + reader->skip(); + continue; + } + if (iskip == nskip) iskip = 0; + iskip++; + if (nevery && ntimestep % nevery) reader->skip(); + else if (iskip < nskip) reader->skip(); + else break; + + // old code logic + //if (iskip == nskip) iskip = 0; + //iskip++; + //if (eofflag) break; + //if (ntimestep <= ncurrent) break; + //if (ntimestep > nlast) break; + //if (nevery && ntimestep % nevery) reader->skip(); + //else if (iskip < nskip) reader->skip(); + //else break; + } + if (eofflag) reader->close_file(); + else break; + } + + currentfile = ifile; + if (eofflag) ntimestep = -1; + if (ntimestep <= ncurrent) ntimestep = -1; + if (ntimestep > nlast) ntimestep = -1; + if (ntimestep < 0) reader->close_file(); + } + + MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world); + return ntimestep; +} + +/* ---------------------------------------------------------------------- + read and broadcast and store snapshot header info + set nsnapatoms = # of atoms in snapshot +------------------------------------------------------------------------- */ + +void ReadDump2::header(int fieldinfo) +{ + int triclinic_snap; + int fieldflag,xflag,yflag,zflag; + + if (me == 0) + nsnapatoms = reader->read_header(box,triclinic_snap, + fieldinfo,nfield,fieldtype,fieldlabel, + scaleflag,wrapflag,fieldflag, + xflag,yflag,zflag); + + MPI_Bcast(&nsnapatoms,1,MPI_LMP_BIGINT,0,world); + MPI_Bcast(&triclinic_snap,1,MPI_INT,0,world); + MPI_Bcast(&box[0][0],9,MPI_DOUBLE,0,world); + + // local copy of snapshot box parameters + // used in xfield,yfield,zfield when converting dump atom to absolute coords + + xlo = box[0][0]; + xhi = box[0][1]; + ylo = box[1][0]; + yhi = box[1][1]; + zlo = box[2][0]; + zhi = box[2][1]; + if (triclinic_snap) { + xy = box[0][2]; + xz = box[1][2]; + yz = box[2][2]; + double xdelta = MIN(0.0,xy); + xdelta = MIN(xdelta,xz); + xdelta = MIN(xdelta,xy+xz); + xlo = xlo - xdelta; + xdelta = MAX(0.0,xy); + xdelta = MAX(xdelta,xz); + xdelta = MAX(xdelta,xy+xz); + xhi = xhi - xdelta; + ylo = ylo - MIN(0.0,yz); + yhi = yhi - MAX(0.0,yz); + } + xprd = xhi - xlo; + yprd = yhi - ylo; + zprd = zhi - zlo; + + // done if not checking fields + + if (!fieldinfo) return; + + MPI_Bcast(&fieldflag,1,MPI_INT,0,world); + MPI_Bcast(&xflag,1,MPI_INT,0,world); + MPI_Bcast(&yflag,1,MPI_INT,0,world); + MPI_Bcast(&zflag,1,MPI_INT,0,world); + + // error check on current vs new box and fields + // triclinic_snap < 0 means no box info in file + + if (triclinic_snap < 0 && boxflag > 0) + error->all(FLERR,"No box information in dump, must use 'box no'"); + if (triclinic_snap >= 0) { + if ((triclinic_snap && !triclinic) || + (!triclinic_snap && triclinic)) + error->one(FLERR,"Read_dump triclinic status does not match simulation"); + } + + // error check on requested fields exisiting in dump file + + if (fieldflag < 0) + error->one(FLERR,"Read_dump field not found in dump file"); + + // all explicitly requested x,y,z must have consistent scaling & wrapping + + int value = MAX(xflag,yflag); + value = MAX(zflag,value); + if ((xflag != UNSET && xflag != value) || + (yflag != UNSET && yflag != value) || + (zflag != UNSET && zflag != value)) + error->one(FLERR, + "Read_dump xyz fields do not have consistent scaling/wrapping"); + + // set scaled/wrapped based on xyz flags + + value = UNSET; + if (xflag != UNSET) value = xflag; + if (yflag != UNSET) value = yflag; + if (zflag != UNSET) value = zflag; + + if (value == UNSET) { + scaled = wrapped = 0; + } else if (value == NOSCALE_NOWRAP) { + scaled = wrapped = 0; + } else if (value == NOSCALE_WRAP) { + scaled = 0; + wrapped = 1; + } else if (value == SCALE_NOWRAP) { + scaled = 1; + wrapped = 0; + } else if (value == SCALE_WRAP) { + scaled = wrapped = 1; + } + + // scaled, triclinic coords require all 3 x,y,z fields, to perform unscaling + // set yindex,zindex = column index of Y and Z fields in fields array + // needed for unscaling to absolute coords in xfield(), yfield(), zfield() + + if (scaled && triclinic == 1) { + int flag = 0; + if (xflag == UNSET) flag = 1; + if (yflag == UNSET) flag = 1; + if (dimension == 3 && zflag == UNSET) flag = 1; + if (flag) + error->one(FLERR,"All read_dump x,y,z fields must be specified for " + "scaled, triclinic coords"); + + for (int i = 0; i < nfield; i++) { + if (fieldtype[i] == Y) yindex = i; + if (fieldtype[i] == Z) zindex = i; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void ReadDump2::atoms() +{ + // initialize counters + + 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_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_atoms(nchunk,nfield,fields); + MPI_Bcast(&fields[0][0],nchunk*nfield,MPI_DOUBLE,0,world); + process_atoms(nchunk); + nread += nchunk; + } + + // if addflag = YESADD, assign IDs to new snapshot atoms + + if (addflag == YESADD) { + 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->tag_enable) atom->tag_extend(); + } + + // if trimflag set, delete atoms not replaced by snapshot atoms + + if (trimflag) { + delete_atoms(); + bigint nblocal = atom->nlocal; + MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); + } + + // can now delete uflag arrays + + memory->destroy(uflag); + memory->destroy(ucflag); + memory->destroy(ucflag_all); + + // 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(); + } + + // overwrite simulation box with dump snapshot box if requested + // reallocate processors to box + + 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(0); + domain->set_local_box(); + } + + // 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; + imageint *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(1); + delete irregular; + if (triclinic) domain->lamda2x(atom->nlocal); + + // check that atom IDs are valid + + atom->tag_check(); +} + +/* ---------------------------------------------------------------------- + process arg list for dump file fields and optional keywords +------------------------------------------------------------------------- */ + +int ReadDump2::fields_and_keywords(int narg, char **arg) +{ + // per-field vectors, leave space for ID and TYPE + + fieldtype = new int[narg+2]; + fieldlabel = new char*[narg+2]; + + // add id and type fields as needed + // scan ahead to see if "add yes/keep" keyword/value is used + // requires extra "type" field from from dump file + + int iarg; + for (iarg = 0; iarg < narg; iarg++) + if (strcmp(arg[iarg],"add") == 0) + if (iarg < narg-1 && (strcmp(arg[iarg+1],"yes") == 0 || + strcmp(arg[iarg+1],"keep") == 0)) break; + + nfield = 0; + fieldtype[nfield++] = ID; + if (iarg < narg) fieldtype[nfield++] = TYPE; + + // parse fields + + iarg = 0; + while (iarg < narg) { + int type = whichtype(arg[iarg]); + if (type < 0) break; + if (type == Q && !atom->q_flag) + error->all(FLERR,"Read dump of atom property that isn't allocated"); + fieldtype[nfield++] = type; + iarg++; + } + + // check for no fields + + 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 || fieldtype[i] == FZ) + 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 = NOADD; + for (int i = 0; i < nfield; i++) fieldlabel[i] = NULL; + scaleflag = 0; + wrapflag = 1; + + 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 = YESADD; + else if (strcmp(arg[iarg+1],"no") == 0) addflag = NOADD; + else if (strcmp(arg[iarg+1],"keep") == 0) addflag = KEEPADD; + 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 type = whichtype(arg[iarg+1]); + int i; + for (i = 0; i < nfield; i++) + if (type == fieldtype[i]) break; + if (i == nfield) error->all(FLERR,"Illegal read_dump command"); + int n = strlen(arg[iarg+2]) + 1; + fieldlabel[i] = new char[n]; + strcpy(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) scaleflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0; + else error->all(FLERR,"Illegal read_dump command"); + iarg += 2; + } else if (strcmp(arg[iarg],"wrapped") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command"); + if (strcmp(arg[iarg+1],"yes") == 0) wrapflag = 1; + else if (strcmp(arg[iarg+1],"no") == 0) wrapflag = 0; + 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"); + delete [] readerstyle; + int n = strlen(arg[iarg+1]) + 1; + readerstyle = new char[n]; + strcpy(readerstyle,arg[iarg+1]); + iarg += 2; + break; + } else error->all(FLERR,"Illegal read_dump command"); + } + + if (purgeflag && (replaceflag || trimflag)) + error->all(FLERR,"If read_dump purges it cannot replace or trim"); + if (addflag == KEEPADD && atom->tag_enable == 0) + error->all(FLERR,"Read_dump cannot use 'add keep' without atom IDs"); + + return narg-iarg; +} + +/* ---------------------------------------------------------------------- + check if str is a field argument + if yes, return index of which + if not, return -1 +------------------------------------------------------------------------- */ + +int ReadDump2::whichtype(char *str) +{ + int type = -1; + if (strcmp(str,"id") == 0) type = ID; + else if (strcmp(str,"type") == 0) type = TYPE; + else if (strcmp(str,"x") == 0) type = X; + else if (strcmp(str,"y") == 0) type = Y; + else if (strcmp(str,"z") == 0) type = Z; + else if (strcmp(str,"vx") == 0) type = VX; + else if (strcmp(str,"vy") == 0) type = VY; + else if (strcmp(str,"vz") == 0) type = VZ; + else if (strcmp(str,"q") == 0) type = Q; + else if (strcmp(str,"ix") == 0) type = IX; + else if (strcmp(str,"iy") == 0) type = IY; + else if (strcmp(str,"iz") == 0) type = IZ; + else if (strcmp(str,"fx") == 0) type = FX; + else if (strcmp(str,"fy") == 0) type = FY; + else if (strcmp(str,"fz") == 0) type = FZ; + return type; +} + +/* ---------------------------------------------------------------------- + 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 ReadDump2::process_atoms(int n) +{ + int i,m,ifield,itype; + int xbox,ybox,zbox; + tagint mtag; + + double **x = atom->x; + double **v = atom->v; + double *q = atom->q; + double **f = atom->f; + tagint *tag = atom->tag; + imageint *image = atom->image; + int nlocal = atom->nlocal; + tagint map_tag_max = atom->map_tag_max; + + for (i = 0; i < n; i++) { + ucflag[i] = 0; + + // check if new atom matches one I own + // setting m = -1 forces new atom not to match + // NOTE: atom ID in fields is stored as double, not as ubuf + // so can only cast it to tagint, thus cannot be full 64-bit ID + + mtag = static_cast (fields[i][0]); + if (mtag <= map_tag_max) m = atom->map(mtag); + else m = -1; + if (m < 0 || m >= nlocal) continue; + + ucflag[i] = 1; + uflag[m] = 1; + + if (replaceflag) { + nreplace++; + + // current image flags + + xbox = (image[m] & IMGMASK) - IMGMAX; + ybox = (image[m] >> IMGBITS & IMGMASK) - IMGMAX; + zbox = (image[m] >> IMG2BITS) - IMGMAX; + + // 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 Q: + q[m] = 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 (fields[i][ifield]); + break; + case IY: + ybox = static_cast (fields[i][ifield]); + break; + case IZ: + zbox = static_cast (fields[i][ifield]); + break; + case FX: + f[m][0] = fields[i][ifield]; + break; + case FY: + f[m][1] = fields[i][ifield]; + break; + case FZ: + f[m][2] = fields[i][ifield]; + break; + } + } + + // replace image flag in case changed by ix,iy,iz fields or unwrapping + + if (!wrapped) xbox = ybox = zbox = 0; + + image[m] = ((imageint) (xbox + IMGMAX) & IMGMASK) | + (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | + (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS); + } + } + + // 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 + // check that dump file snapshot has atom type field + + if (addflag == NOADD) return; + + int tflag = 0; + for (ifield = 0; ifield < nfield; ifield++) + if (fieldtype[ifield] == TYPE) tflag = 1; + if (!tflag) + error->all(FLERR,"Cannot add atoms if dump file does not store atom type"); + + MPI_Allreduce(ucflag,ucflag_all,n,MPI_INT,MPI_SUM,world); + + int nlocal_previous = atom->nlocal; + double one[3]; + + 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 + + itype = 0; + one[0] = one[1] = one[2] = 0.0; + for (ifield = 1; ifield < nfield; ifield++) { + switch (fieldtype[ifield]) { + case TYPE: + itype = static_cast (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 + + m = atom->nlocal; + atom->avec->create_atom(itype,one); + nadd++; + + tag = atom->tag; + v = atom->v; + q = atom->q; + image = atom->image; + + // set atom attributes from other dump file fields + + xbox = ybox = zbox = 0; + + for (ifield = 0; ifield < nfield; ifield++) { + switch (fieldtype[ifield]) { + case ID: + if (addflag == KEEPADD) + tag[m] = static_cast (fields[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 Q: + q[m] = fields[i][ifield]; + break; + case IX: + xbox = static_cast (fields[i][ifield]); + break; + case IY: + ybox = static_cast (fields[i][ifield]); + break; + case IZ: + zbox = static_cast (fields[i][ifield]); + break; + } + + // replace image flag in case changed by ix,iy,iz fields + + image[m] = ((imageint) (xbox + IMGMAX) & IMGMASK) | + (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | + (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS); + } + } + + // init per-atom fix/compute/variable values for created atoms + + atom->data_fix_compute_variable(nlocal_previous,atom->nlocal); +} + +/* ---------------------------------------------------------------------- + delete atoms not flagged as replaced by dump atoms +------------------------------------------------------------------------- */ + +void ReadDump2::delete_atoms() +{ + 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 ReadDump2::xfield(int i, int j) +{ + if (!scaled) 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 ReadDump2::yfield(int i, int j) +{ + if (!scaled) 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 ReadDump2::zfield(int i, int j) +{ + if (!scaled) return fields[i][j]; + return fields[i][j]*zprd + zlo; +} diff --git a/src/read_dump2.h b/src/read_dump2.h new file mode 100644 index 0000000000..27105860f5 --- /dev/null +++ b/src/read_dump2.h @@ -0,0 +1,172 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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_dump2,ReadDump2) + +#else + +#ifndef LMP_READ_DUMP2_H +#define LMP_READ_DUMP2_H + +#include +#include "pointers.h" + +namespace LAMMPS_NS { + +class ReadDump2 : protected Pointers { + public: + ReadDump2(class LAMMPS *); + ~ReadDump2(); + void command(int, char **); + + void store_files(int, char **); + void setup_reader(int, char **); + bigint seek(bigint, int); + void header(int); + bigint next(bigint, bigint, int, int); + void atoms(); + int fields_and_keywords(int, char **); + +private: + int me,nprocs; + FILE *fp; + + int dimension; + int triclinic; + + int nfile; // # of dump files to process + char **files; // list of file names + int currentfile; // currently open file + + int boxflag; // overwrite simulation with dump file box params + int replaceflag,addflag; // flags for processing dump snapshot atoms + int trimflag,purgeflag; + int scaleflag; // user 0/1 if dump file coords are unscaled/scaled + int wrapflag; // user 0/1 if dump file coords are unwrapped/wrapped + char *readerstyle; // style of dump files to read + + 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 + + int scaled; // 0/1 if dump file coords are unscaled/scaled + int wrapped; // 0/1 if dump file coords are unwrapped/wrapped + + double box[3][3]; // dump file box parameters + double xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz; // dump snapshot box params + 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 Reader *reader; // class that reads dump file + + int whichtype(char *); + void process_atoms(int); + void delete_atoms(); + + double xfield(int, int); + double yfield(int, int); + double zfield(int, int); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Read_dump command before simulation box is defined + +The read_dump command cannot be used before a read_data, read_restart, +or create_box command. + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Dump file does not contain requested snapshot + +Self-explanatory. + +E: Unknown dump reader style + +The choice of dump reader style via the format keyword is unknown. + +E: No box information in dump, must use 'box no' + +UNDOCUMENTED + +E: Read_dump triclinic status does not match simulation + +Both the dump snapshot and the current LAMMPS simulation must +be using either an orthogonal or triclinic box. + +E: Read_dump field not found in dump file + +Self-explanatory. + +E: Read_dump xyz fields do not have consistent scaling/wrapping + +Self-explanatory. + +E: All read_dump x,y,z fields must be specified for scaled, triclinic coords + +For triclinic boxes and scaled coordinates you must specify all 3 of +the x,y,z fields, else LAMMPS cannot reconstruct the unscaled +coordinates. + +E: Too many total atoms + +See the setting for bigint in the src/lmptype.h file. + +E: Read dump of atom property that isn't allocated + +Self-explanatory. + +E: Duplicate fields in read_dump command + +Self-explanatory. + +E: If read_dump purges it cannot replace or trim + +These operations are not compatible. See the read_dump doc +page for details. + +E: Read_dump cannot use 'add keep' without atom IDs + +UNDOCUMENTED + +E: Cannot add atoms if dump file does not store atom type + +UNDOCUMENTED + +U: No box information in dump. You have to use 'box no' + +Self-explanatory. + +*/ diff --git a/src/read_restart.h b/src/read_restart.h index 4cee62a1cf..23d6ec3fba 100644 --- a/src/read_restart.h +++ b/src/read_restart.h @@ -31,11 +31,13 @@ class ReadRestart : protected Pointers { void command(int, char **); private: - int me,nprocs,nprocs_file,multiproc_file; + int me,nprocs; FILE *fp; - int multiproc; // 0 = proc 0 writes for all - // else # of procs writing files + int multiproc; // 0 = restart file is a single file + // 1 = restart file is parallel (multiple files) + int multiproc_file; // # of parallel files in restart + int nprocs_file; // total # of procs that wrote restart file // MPI-IO values diff --git a/src/reader.cpp b/src/reader.cpp index 22a21812e6..cf344b37b3 100644 --- a/src/reader.cpp +++ b/src/reader.cpp @@ -18,6 +18,8 @@ using namespace LAMMPS_NS; +// only proc 0 calls methods of this class, except for constructor/destructor + /* ---------------------------------------------------------------------- */ Reader::Reader(LAMMPS *lmp) : Pointers(lmp) diff --git a/src/rerun.cpp b/src/rerun.cpp index f8a37b5946..063ea882c2 100644 --- a/src/rerun.cpp +++ b/src/rerun.cpp @@ -57,7 +57,7 @@ void Rerun::command(int narg, char **arg) if (nfile == 0 || nfile == narg) error->all(FLERR,"Illegal rerun command"); // parse optional args up until "dump" - // user MAXBIGINT -1 so Output can add 1 to it and still be a big int + // use MAXBIGINT -1 so Output can add 1 to it and still be a big int bigint first = 0; bigint last = MAXBIGINT - 1; diff --git a/src/rerun2.cpp b/src/rerun2.cpp new file mode 100644 index 0000000000..1747f63b97 --- /dev/null +++ b/src/rerun2.cpp @@ -0,0 +1,193 @@ +/* ---------------------------------------------------------------------- + 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 +#include +#include "rerun2.h" +#include "read_dump2.h" +#include "domain.h" +#include "update.h" +#include "integrate.h" +#include "modify.h" +#include "output.h" +#include "finish.h" +#include "timer.h" +#include "error.h" +#include "force.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +Rerun2::Rerun2(LAMMPS *lmp) : Pointers(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void Rerun2::command(int narg, char **arg) +{ + if (domain->box_exist == 0) + error->all(FLERR,"Rerun command before simulation box is defined"); + + if (narg < 2) error->all(FLERR,"Illegal rerun command"); + + // list of dump files = args until a keyword + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"first") == 0) break; + if (strcmp(arg[iarg],"last") == 0) break; + if (strcmp(arg[iarg],"every") == 0) break; + if (strcmp(arg[iarg],"skip") == 0) break; + if (strcmp(arg[iarg],"start") == 0) break; + if (strcmp(arg[iarg],"stop") == 0) break; + if (strcmp(arg[iarg],"dump") == 0) break; + iarg++; + } + int nfile = iarg; + if (nfile == 0 || nfile == narg) error->all(FLERR,"Illegal rerun command"); + + // parse optional args up until "dump" + // use MAXBIGINT -1 so Output can add 1 to it and still be a big int + + bigint first = 0; + bigint last = MAXBIGINT - 1; + int nevery = 0; + int nskip = 1; + int startflag = 0; + int stopflag = 0; + bigint start = -1; + bigint stop = -1; + + while (iarg < narg) { + if (strcmp(arg[iarg],"first") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal rerun command"); + first = force->bnumeric(FLERR,arg[iarg+1]); + if (first < 0) error->all(FLERR,"Illegal rerun command"); + iarg += 2; + } else if (strcmp(arg[iarg],"last") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal rerun command"); + last = force->bnumeric(FLERR,arg[iarg+1]); + if (last < 0) error->all(FLERR,"Illegal rerun command"); + iarg += 2; + } else if (strcmp(arg[iarg],"every") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal rerun command"); + nevery = force->inumeric(FLERR,arg[iarg+1]); + if (nevery < 0) error->all(FLERR,"Illegal rerun command"); + iarg += 2; + } else if (strcmp(arg[iarg],"skip") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal rerun command"); + nskip = force->inumeric(FLERR,arg[iarg+1]); + if (nskip <= 0) error->all(FLERR,"Illegal rerun command"); + iarg += 2; + } else if (strcmp(arg[iarg],"start") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal rerun command"); + startflag = 1; + start = force->bnumeric(FLERR,arg[iarg+1]); + if (start < 0) error->all(FLERR,"Illegal rerun command"); + iarg += 2; + } else if (strcmp(arg[iarg],"stop") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal rerun command"); + stopflag = 1; + stop = force->bnumeric(FLERR,arg[iarg+1]); + if (stop < 0) error->all(FLERR,"Illegal rerun command"); + iarg += 2; + } else if (strcmp(arg[iarg],"dump") == 0) { + break; + } else error->all(FLERR,"Illegal rerun command"); + } + + int nremain = narg - iarg - 1; + if (nremain <= 0) error->all(FLERR,"Illegal rerun command"); + if (first > last) error->all(FLERR,"Illegal rerun command"); + if (startflag && stopflag && start > stop) + error->all(FLERR,"Illegal rerun command"); + + // pass list of filenames to ReadDump + // along with post-"dump" args and post-"format" args + + ReadDump2 *rd = new ReadDump2(lmp); + + rd->store_files(nfile,arg); + if (nremain) + nremain = rd->fields_and_keywords(nremain,&arg[narg-nremain]); + else nremain = rd->fields_and_keywords(0,NULL); + if (nremain) rd->setup_reader(nremain,&arg[narg-nremain]); + else rd->setup_reader(0,NULL); + + // perform the pseudo run + // invoke lmp->init() only once + // read all relevant snapshots + // use setup_minimal() since atoms are already owned by correct procs + // addstep_compute_all() insures energy/virial computed on every snapshot + + update->whichflag = 1; + + if (startflag) update->beginstep = update->firststep = start; + else update->beginstep = update->firststep = first; + if (stopflag) update->endstep = update->laststep = stop; + else update->endstep = update->laststep = last; + + int firstflag = 1; + int ndump = 0; + + lmp->init(); + + timer->init(); + timer->barrier_start(); + + bigint ntimestep = rd->seek(first,0); + if (ntimestep < 0) + error->all(FLERR,"Rerun dump file does not contain requested snapshot"); + + while (1) { + ndump++; + rd->header(firstflag); + update->reset_timestep(ntimestep); + rd->atoms(); + modify->init(); + update->integrate->setup_minimal(1); + modify->end_of_step(); + if (firstflag) output->setup(); + else if (output->next) output->write(ntimestep); + + firstflag = 0; + ntimestep = rd->next(ntimestep,last,nevery,nskip); + if (stopflag && ntimestep > stop) + error->all(FLERR,"Read rerun dump file timestep > specified stop"); + if (ntimestep < 0) break; + } + + // insure thermo output on last dump timestep + + output->next_thermo = update->ntimestep; + output->write(update->ntimestep); + + timer->barrier_stop(); + + update->integrate->cleanup(); + + // set update->nsteps to ndump for Finish stats to print + + update->nsteps = ndump; + + Finish finish(lmp); + finish.end(1); + + update->whichflag = 0; + update->firststep = update->laststep = 0; + update->beginstep = update->endstep = 0; + + // clean-up + + delete rd; +} diff --git a/src/rerun2.h b/src/rerun2.h new file mode 100644 index 0000000000..fa4e2c49a6 --- /dev/null +++ b/src/rerun2.h @@ -0,0 +1,59 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef COMMAND_CLASS + +CommandStyle(rerun2,Rerun2) + +#else + +#ifndef LMP_RERUN2_H +#define LMP_RERUN2_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class Rerun2 : protected Pointers { + public: + Rerun2(class LAMMPS *); + void command(int, char **); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Rerun command before simulation box is defined + +The rerun command cannot be used before a read_data, read_restart, or +create_box command. + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Rerun dump file does not contain requested snapshot + +Self-explanatory. + +E: Read rerun dump file timestep > specified stop + +Self-explanatory. + +*/