diff --git a/src/USER-MOLFILE/reader_molfile.cpp b/src/USER-MOLFILE/reader_molfile.cpp index a1156efc67..2e3aa89506 100644 --- a/src/USER-MOLFILE/reader_molfile.cpp +++ b/src/USER-MOLFILE/reader_molfile.cpp @@ -185,7 +185,7 @@ void ReaderMolfile::skip() read remaining header info: return natoms box bounds, triclinic (inferred), fieldflag (1 if any fields not found), - xyz flag = UNSET (not a requested field), SCALED, UNSCALED + xyz flags = from input scaleflag & wrapflag if fieldflag set: match Nfield fields to per-atom column labels allocate and set fieldindex = which column each field maps to @@ -198,7 +198,7 @@ void ReaderMolfile::skip() bigint ReaderMolfile::read_header(double box[3][3], int &triclinic, int fieldinfo, int nfield, int *fieldtype, char **fieldlabel, - int scaledflag, int &fieldflag, + int scaleflag, int wrapflag, int &fieldflag, int &xflag, int &yflag, int &zflag) { nid = 0; @@ -279,13 +279,15 @@ bigint ReaderMolfile::read_header(double box[3][3], int &triclinic, memory->create(fieldindex,nfield,"read_dump:fieldindex"); - // we know nothing about the scaling style of coordinates, - // so the caller has to set the proper flag. - xflag = scaledflag; - yflag = scaledflag; - zflag = scaledflag; + // we know nothing about the style of coordinates, + // so caller has to set the proper flags + + xflag = 2*scaleflag + wrapflag + 1; + yflag = 2*scaleflag + wrapflag + 1; + zflag = 2*scaleflag + wrapflag + 1; // copy fieldtype list for supported fields + fieldflag = 0; needvels = 0; for (int i = 0; i < nfield; i++) { diff --git a/src/USER-MOLFILE/reader_molfile.h b/src/USER-MOLFILE/reader_molfile.h index 31399d0e85..5450843d6d 100644 --- a/src/USER-MOLFILE/reader_molfile.h +++ b/src/USER-MOLFILE/reader_molfile.h @@ -36,7 +36,7 @@ class ReaderMolfile : public Reader { virtual int read_time(bigint &); virtual void skip(); virtual bigint read_header(double [3][3], int &, int, int, int *, char **, - int, int &, int &, int &, int &); + int, int, int &, int &, int &, int &); virtual void read_atoms(int, int, double **); virtual void open_file(const char *); diff --git a/src/read_dump.cpp b/src/read_dump.cpp index f7ddc232aa..c5bc964300 100644 --- a/src/read_dump.cpp +++ b/src/read_dump.cpp @@ -38,8 +38,10 @@ using namespace LAMMPS_NS; #define CHUNK 1024 #define EPSILON 1.0e-6 +// also in reader_native.cpp + enum{ID,TYPE,X,Y,Z,VX,VY,VZ,IX,IY,IZ}; -enum{UNSET,UNSCALED,SCALED}; +enum{UNSET,NOSCALE_NOWRAP,NOSCALE_WRAP,SCALE_NOWRAP,SCALE_WRAP}; /* ---------------------------------------------------------------------- */ @@ -298,7 +300,6 @@ bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip) /* ---------------------------------------------------------------------- read and broadcast and store snapshot header info set nsnapatoms = # of atoms in snapshot - set yindex,zindex ------------------------------------------------------------------------- */ void ReadDump::header(int fieldinfo) @@ -309,7 +310,8 @@ void ReadDump::header(int fieldinfo) if (me == 0) nsnapatoms = reader->read_header(box,triclinic_snap, fieldinfo,nfield,fieldtype,fieldlabel, - scaledflag,fieldflag,xflag,yflag,zflag); + scaleflag,wrapflag,fieldflag, + xflag,yflag,zflag); MPI_Bcast(&nsnapatoms,1,MPI_LMP_BIGINT,0,world); MPI_Bcast(&triclinic_snap,1,MPI_INT,0,world); @@ -363,30 +365,51 @@ void ReadDump::header(int fieldinfo) error->one(FLERR,"Read_dump triclinic status does not match simulation"); } - // error check field and scaling info + // error check on requested fields exisiting in dump file if (fieldflag < 0) error->one(FLERR,"Read_dump field not found in dump file"); - // set overall scaling of coordinates - // error if x,y,z scaling are not the same + // all explicitly requested x,y,z must have consistent scaling & wrapping - scaled = MAX(xflag,yflag); - scaled = MAX(zflag,scaled); - if ((xflag != UNSET && xflag != scaled) || - (yflag != UNSET && yflag != scaled) || - (zflag != UNSET && zflag != scaled)) - error->one(FLERR,"Read_dump x,y,z fields do not have consistent scaling"); + 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 == SCALED && triclinic == 1) { + if (scaled && triclinic == 1) { int flag = 0; - if (xflag != scaled) flag = 1; - if (yflag != scaled) flag = 1; - if (dimension == 3 && zflag != scaled) flag = 1; + 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"); @@ -597,7 +620,8 @@ int ReadDump::fields_and_keywords(int narg, char **arg) trimflag = 0; addflag = 0; for (int i = 0; i < nfield; i++) fieldlabel[i] = NULL; - scaledflag = UNSCALED; + scaleflag = 0; + wrapflag = 1; while (iarg < narg) { if (strcmp(arg[iarg],"box") == 0) { @@ -642,8 +666,14 @@ int ReadDump::fields_and_keywords(int narg, char **arg) iarg += 3; } else if (strcmp(arg[iarg],"scaled") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command"); - if (strcmp(arg[iarg+1],"yes") == 0) scaledflag = SCALED; - else if (strcmp(arg[iarg+1],"no") == 0) scaledflag = UNSCALED; + 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) { @@ -742,7 +772,9 @@ void ReadDump::process_atoms(int n) } } - // replace image flag in case changed by ix,iy,iz fields + // replace image flag in case changed by ix,iy,iz fields or unwrapping + + if (!wrapped) xbox = ybox = zbox = 0; image[m] = ((tagint) (xbox + IMGMAX) & IMGMASK) | (((tagint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) | @@ -877,7 +909,7 @@ void ReadDump::delete_atoms() double ReadDump::xfield(int i, int j) { - if (scaled == UNSCALED) return fields[i][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; @@ -886,7 +918,7 @@ double ReadDump::xfield(int i, int j) double ReadDump::yfield(int i, int j) { - if (scaled == UNSCALED) return fields[i][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; @@ -894,6 +926,6 @@ double ReadDump::yfield(int i, int j) double ReadDump::zfield(int i, int j) { - if (scaled == UNSCALED) return fields[i][j]; + if (!scaled) return fields[i][j]; return fields[i][j]*zprd + zlo; } diff --git a/src/read_dump.h b/src/read_dump.h index 66d78b93d3..3a779da1f9 100644 --- a/src/read_dump.h +++ b/src/read_dump.h @@ -55,8 +55,8 @@ private: int boxflag; // overwrite simulation with dump file box params int replaceflag,addflag; // flags for processing dump snapshot atoms int trimflag,purgeflag; - int scaledflag; // user setting for coordinate scaling - int scaled; // actual setting for coordinate scaling + int 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 @@ -64,6 +64,9 @@ private: 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; diff --git a/src/reader.h b/src/reader.h index 8801970d29..571016ba7f 100644 --- a/src/reader.h +++ b/src/reader.h @@ -30,7 +30,7 @@ class Reader : protected Pointers { virtual int read_time(bigint &) = 0; virtual void skip() = 0; virtual bigint read_header(double [3][3], int &, int, int, int *, char **, - int, int &, int &, int &, int &) = 0; + int, int, int &, int &, int &, int &) = 0; virtual void read_atoms(int, int, double **) = 0; virtual void open_file(const char *); diff --git a/src/reader_native.cpp b/src/reader_native.cpp index 85dc19dfe9..167be4ebe7 100644 --- a/src/reader_native.cpp +++ b/src/reader_native.cpp @@ -22,8 +22,10 @@ using namespace LAMMPS_NS; #define MAXLINE 1024 // max line length in dump file +// also in read_dump.cpp + enum{ID,TYPE,X,Y,Z,VX,VY,VZ,IX,IY,IZ}; -enum{UNSET,UNSCALED,SCALED}; +enum{UNSET,NOSCALE_NOWRAP,NOSCALE_WRAP,SCALE_NOWRAP,SCALE_WRAP}; /* ---------------------------------------------------------------------- */ @@ -90,20 +92,21 @@ void ReaderNative::skip() read remaining header info: return natoms box bounds, triclinic (inferred), fieldflag (1 if any fields not found), - xyz flag = UNSET (not a requested field), SCALED, UNSCALED + xyz flags = UNSET (not a requested field), SCALE/WRAP as in enum if fieldflag set: match Nfield fields to per-atom column labels allocate and set fieldindex = which column each field maps to fieldtype = X,VX,IZ etc fieldlabel = user-specified label or NULL if use fieldtype default - xyz flag = scaledflag if has fieldlabel name, else set by x,xs,xu,xsu + xyz flags = scaleflag+wrapflag if has fieldlabel name, + else set by x,xs,xu,xsu only called by proc 0 ------------------------------------------------------------------------- */ bigint ReaderNative::read_header(double box[3][3], int &triclinic, int fieldinfo, int nfield, int *fieldtype, char **fieldlabel, - int scaledflag, int &fieldflag, + int scaleflag, int wrapflag, int &fieldflag, int &xflag, int &yflag, int &zflag) { bigint natoms; @@ -147,7 +150,7 @@ bigint ReaderNative::read_header(double box[3][3], int &triclinic, // match each field with a column of per-atom data // if fieldlabel set, match with explicit column // else infer one or more column matches from fieldtype - // xyz flag set by scaledflag (if fieldlabel set) or column label + // xyz flag set by scaleflag + wrapflag (if fieldlabel set) or column label memory->create(fieldindex,nfield,"read_dump:fieldindex"); @@ -159,9 +162,9 @@ bigint ReaderNative::read_header(double box[3][3], int &triclinic, for (int i = 0; i < nfield; i++) { if (fieldlabel[i]) { fieldindex[i] = find_label(fieldlabel[i],nwords,labels); - if (fieldtype[i] == X) xflag = scaledflag; - else if (fieldtype[i] == Y) yflag = scaledflag; - else if (fieldtype[i] == Z) zflag = scaledflag; + if (fieldtype[i] == X) xflag = 2*scaleflag + wrapflag + 1; + else if (fieldtype[i] == Y) yflag = 2*scaleflag + wrapflag + 1; + else if (fieldtype[i] == Z) zflag = 2*scaleflag + wrapflag + 1; } else if (fieldtype[i] == ID) @@ -171,7 +174,7 @@ bigint ReaderNative::read_header(double box[3][3], int &triclinic, else if (fieldtype[i] == X) { fieldindex[i] = find_label("x",nwords,labels); - xflag = UNSCALED; + xflag = NOSCALE_WRAP; if (fieldindex[i] < 0) { fieldindex[i] = nwords; s_index = find_label("xs",nwords,labels); @@ -179,22 +182,22 @@ bigint ReaderNative::read_header(double box[3][3], int &triclinic, su_index = find_label("xsu",nwords,labels); if (s_index >= 0 && s_index < fieldindex[i]) { fieldindex[i] = s_index; - xflag = SCALED; + xflag = SCALE_WRAP; } if (u_index >= 0 && u_index < fieldindex[i]) { fieldindex[i] = u_index; - xflag = UNSCALED; + xflag = NOSCALE_NOWRAP; } if (su_index >= 0 && su_index < fieldindex[i]) { fieldindex[i] = su_index; - xflag = SCALED; + xflag = SCALE_NOWRAP; } } if (fieldindex[i] == nwords) fieldindex[i] = -1; } else if (fieldtype[i] == Y) { fieldindex[i] = find_label("y",nwords,labels); - yflag = UNSCALED; + yflag = NOSCALE_WRAP; if (fieldindex[i] < 0) { fieldindex[i] = nwords; s_index = find_label("ys",nwords,labels); @@ -202,22 +205,22 @@ bigint ReaderNative::read_header(double box[3][3], int &triclinic, su_index = find_label("ysu",nwords,labels); if (s_index >= 0 && s_index < fieldindex[i]) { fieldindex[i] = s_index; - yflag = SCALED; + yflag = SCALE_WRAP; } if (u_index >= 0 && u_index < fieldindex[i]) { fieldindex[i] = u_index; - yflag = UNSCALED; + yflag = NOSCALE_NOWRAP; } if (su_index >= 0 && su_index < fieldindex[i]) { fieldindex[i] = su_index; - yflag = SCALED; + yflag = SCALE_NOWRAP; } } if (fieldindex[i] == nwords) fieldindex[i] = -1; } else if (fieldtype[i] == Z) { fieldindex[i] = find_label("z",nwords,labels); - zflag = UNSCALED; + zflag = NOSCALE_WRAP; if (fieldindex[i] < 0) { fieldindex[i] = nwords; s_index = find_label("zs",nwords,labels); @@ -225,15 +228,15 @@ bigint ReaderNative::read_header(double box[3][3], int &triclinic, su_index = find_label("zsu",nwords,labels); if (s_index >= 0 && s_index < fieldindex[i]) { fieldindex[i] = s_index; - zflag = SCALED; + zflag = SCALE_WRAP; } if (u_index >= 0 && u_index < fieldindex[i]) { fieldindex[i] = u_index; - zflag = UNSCALED; + zflag = NOSCALE_NOWRAP; } if (su_index >= 0 && su_index < fieldindex[i]) { fieldindex[i] = su_index; - zflag = SCALED; + zflag = SCALE_NOWRAP; } } if (fieldindex[i] == nwords) fieldindex[i] = -1; diff --git a/src/reader_native.h b/src/reader_native.h index 8873dba4fc..066ae1c928 100644 --- a/src/reader_native.h +++ b/src/reader_native.h @@ -34,7 +34,7 @@ class ReaderNative : public Reader { int read_time(bigint &); void skip(); bigint read_header(double [3][3], int &, int, int, int *, char **, - int, int &, int &, int &, int &); + int, int, int &, int &, int &, int &); void read_atoms(int, int, double **); private: diff --git a/src/reader_xyz.cpp b/src/reader_xyz.cpp index e37bda0415..a8f425cf3a 100644 --- a/src/reader_xyz.cpp +++ b/src/reader_xyz.cpp @@ -98,7 +98,7 @@ void ReaderXYZ::skip() read remaining header info: return natoms box bounds, triclinic (inferred), fieldflag (1 if any fields not found), - xyz flag = UNSET (not a requested field), SCALED, UNSCALED + xyz flags = from input scaleflag & wrapflag if fieldflag set: match Nfield fields to per-atom column labels allocate and set fieldindex = which column each field maps to @@ -111,7 +111,7 @@ void ReaderXYZ::skip() bigint ReaderXYZ::read_header(double box[3][3], int &triclinic, int fieldinfo, int nfield, int *fieldtype, char **fieldlabel, - int scaledflag, int &fieldflag, + int scaleflag, int wrapflag, int &fieldflag, int &xflag, int &yflag, int &zflag) { @@ -128,11 +128,11 @@ bigint ReaderXYZ::read_header(double box[3][3], int &triclinic, memory->create(fieldindex,nfield,"read_dump:fieldindex"); // for xyz we know nothing about the style of coordinates, - // so the caller has to set the proper flag. + // so caller has to set the proper flags - xflag = scaledflag; - yflag = scaledflag; - zflag = scaledflag; + xflag = 2*scaleflag + wrapflag + 1; + yflag = 2*scaleflag + wrapflag + 1; + zflag = 2*scaleflag + wrapflag + 1; // copy fieldtype list for supported fields diff --git a/src/reader_xyz.h b/src/reader_xyz.h index ad4b56c396..698c6ef240 100644 --- a/src/reader_xyz.h +++ b/src/reader_xyz.h @@ -34,7 +34,7 @@ class ReaderXYZ : public Reader { int read_time(bigint &); void skip(); bigint read_header(double [3][3], int &, int, int, int *, char **, - int, int &, int &, int &, int &); + int, int, int &, int &, int &, int &); void read_atoms(int, int, double **); private: