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 8d2547fb10..350d2787f7 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/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. + +*/