diff --git a/src/read_dump.cpp b/src/read_dump.cpp index 2327f0ba13..4125d5512a 100644 --- a/src/read_dump.cpp +++ b/src/read_dump.cpp @@ -90,8 +90,12 @@ void ReadDump::command(int narg, char **arg) store_files(1,&arg[0]); bigint nstep = ATOBIGINT(arg[1]); - fields_and_keywords(narg-2,&arg[2]); - setup_reader(); + + int nremain = narg - 2; + if (nremain) nremain = fields_and_keywords(narg-nremain,&arg[narg-nremain]); + else nremain = fields_and_keywords(0,NULL); + if (nremain) setup_reader(narg-nremain,&arg[narg-nremain]); + else setup_reader(0,NULL); // find the snapshot and read/bcast/process header info @@ -116,7 +120,7 @@ void ReadDump::command(int narg, char **arg) bigint natoms_prev = atom->natoms; atoms(); - if (me == 0) close(); + if (me == 0) reader->close_file(); // print out stats @@ -172,7 +176,7 @@ void ReadDump::store_files(int nstr, char **str) /* ---------------------------------------------------------------------- */ -void ReadDump::setup_reader() +void ReadDump::setup_reader(int narg, char **arg) { // allocate snapshot field buffer @@ -192,6 +196,8 @@ void ReadDump::setup_reader() // unrecognized style else error->all(FLERR,"Invalid dump reader style"); + + if (narg > 0) reader->settings(narg,arg); } /* ---------------------------------------------------------------------- @@ -213,8 +219,7 @@ bigint ReadDump::seek(bigint nrequest, int exact) for (ifile = 0; ifile < nfile; ifile++) { ntimestep = -1; - open(files[ifile]); - reader->file(fp); + reader->open_file(files[ifile]); while (1) { eofflag = reader->read_time(ntimestep); if (eofflag) break; @@ -222,13 +227,13 @@ bigint ReadDump::seek(bigint nrequest, int exact) reader->skip(); } if (ntimestep >= nrequest) break; - close(); + reader->close_file(); } currentfile = ifile; if (ntimestep < nrequest) ntimestep = -1; if (exact && ntimestep != nrequest) ntimestep = -1; - if (ntimestep < 0) close(); + if (ntimestep < 0) reader->close_file(); } MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world); @@ -258,8 +263,7 @@ bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip) for (ifile = currentfile; ifile < nfile; ifile++) { ntimestep = -1; - if (ifile != currentfile) open(files[ifile]); - reader->file(fp); + if (ifile != currentfile) reader->open_file(files[ifile]); while (1) { eofflag = reader->read_time(ntimestep); if (iskip == nskip) iskip = 0; @@ -271,7 +275,7 @@ bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip) else if (iskip < nskip) reader->skip(); else break; } - if (eofflag) close(); + if (eofflag) reader->close_file(); else break; } @@ -279,7 +283,7 @@ 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) close(); + if (ntimestep < 0) reader->close_file(); } MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world); @@ -334,10 +338,15 @@ void ReadDump::header(int fieldinfo) 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 && !triclinic) || - (!triclinic_snap && triclinic)) - error->one(FLERR,"Read_dump triclinic status does not match simulation"); + if (triclinic_snap < 0 && boxflag > 0) + error->all(FLERR,"No box information in dump. You have to 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 field and scaling info @@ -358,7 +367,7 @@ void ReadDump::header(int fieldinfo) // 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) { + if (scaled == SCALED && triclinic == 1) { int flag = 0; if (xflag != scaled) flag = 1; if (yflag != scaled) flag = 1; @@ -511,9 +520,9 @@ void ReadDump::atoms() process arg list for dump file fields and optional keywords ------------------------------------------------------------------------- */ -void ReadDump::fields_and_keywords(int narg, char **arg) +int ReadDump::fields_and_keywords(int narg, char **arg) { - // per-field vectors, leave space for ID + TYPE + // per-field vectors, leave space for ID and TYPE fieldtype = new int[narg+2]; fieldlabel = new char*[narg+2]; @@ -628,11 +637,14 @@ void ReadDump::fields_and_keywords(int narg, char **arg) 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"); + + return narg-iarg; } /* ---------------------------------------------------------------------- @@ -727,9 +739,8 @@ void ReadDump::process_atoms(int n) MPI_Allreduce(ucflag,ucflag_all,n,MPI_INT,MPI_SUM,world); int nlocal_previous = atom->nlocal; - double lamda[3],one[3]; - double *coord; - + double one[3]; + for (i = 0; i < n; i++) { if (ucflag_all[i]) continue; @@ -866,43 +877,3 @@ double ReadDump::zfield(int i, int j) return fields[i][j]*zprd + zlo; } -/* ---------------------------------------------------------------------- - proc 0 opens dump file - test if gzipped -------------------------------------------------------------------------- */ - -void ReadDump::open(char *file) -{ - compressed = 0; - char *suffix = file + strlen(file) - 3; - if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1; - if (!compressed) fp = fopen(file,"r"); - else { -#ifdef LAMMPS_GZIP - char gunzip[128]; - sprintf(gunzip,"gunzip -c %s",file); - fp = popen(gunzip,"r"); -#else - error->one(FLERR,"Cannot open gzipped file"); -#endif - } - - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open file %s",file); - error->one(FLERR,str); - } -} - -/* ---------------------------------------------------------------------- - close current dump file - only called by proc 0 -------------------------------------------------------------------------- */ - -void ReadDump::close() -{ - if (fp == NULL) return; - if (compressed) pclose(fp); - else fclose(fp); - fp = NULL; -} diff --git a/src/read_dump.h b/src/read_dump.h index 5cdfc0f932..fd3a0c798c 100644 --- a/src/read_dump.h +++ b/src/read_dump.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -34,12 +34,12 @@ class ReadDump : protected Pointers { void command(int, char **); void store_files(int, char **); - void setup_reader(); + void setup_reader(int, char **); bigint seek(bigint, int); void header(int); bigint next(bigint, bigint, int, int); void atoms(); - void fields_and_keywords(int, char **); + int fields_and_keywords(int, char **); private: int me,nprocs; @@ -57,7 +57,6 @@ private: int trimflag,purgeflag; int scaledflag; // user setting for coordinate scaling int scaled; // actual setting for coordinate scaling - int compressed; // flag for dump file compression char *readerstyle; // style of dump files to read int nfield; // # of fields to extract from dump file @@ -86,9 +85,6 @@ private: double xfield(int, int); double yfield(int, int); double zfield(int, int); - - void open(char *); - void close(); }; } diff --git a/src/reader.cpp b/src/reader.cpp index b8947716a5..79dbb1c046 100644 --- a/src/reader.cpp +++ b/src/reader.cpp @@ -12,20 +12,58 @@ ------------------------------------------------------------------------- */ #include "stdio.h" +#include "string.h" #include "reader.h" +#include "error.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -Reader::Reader(LAMMPS *lmp) : Pointers(lmp) {} +Reader::Reader(LAMMPS *lmp) : Pointers(lmp) +{ + fp = NULL; +} /* ---------------------------------------------------------------------- - set file ptr - caller opens/closes dump files + try to open given file + generic version for ASCII files that may be compressed ------------------------------------------------------------------------- */ -void Reader::file(FILE *fpcaller) +void Reader::open_file(const char *file) { - fp = fpcaller; + if (fp != NULL) close_file(); + + compressed = 0; + const char *suffix = file + strlen(file) - 3; + if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1; + if (!compressed) fp = fopen(file,"r"); + else { +#ifdef LAMMPS_GZIP + char gunzip[1024]; + sprintf(gunzip,"gunzip -c %s",file); + fp = popen(gunzip,"r"); +#else + error->one(FLERR,"Cannot open gzipped file"); +#endif + } + + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open file %s",file); + error->one(FLERR,str); + } +} + +/* ---------------------------------------------------------------------- + close current file if open + generic version for ASCII files that may be compressed +------------------------------------------------------------------------- */ + +void Reader::close_file() +{ + if (fp == NULL) return; + if (compressed) pclose(fp); + else fclose(fp); + fp = NULL; } diff --git a/src/reader.h b/src/reader.h index 28c6c35709..1c53562cfd 100644 --- a/src/reader.h +++ b/src/reader.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -16,7 +16,6 @@ #ifndef LMP_READER_H #define LMP_READER_H -#include "stdio.h" #include "pointers.h" namespace LAMMPS_NS { @@ -26,16 +25,20 @@ class Reader : protected Pointers { Reader(class LAMMPS *); virtual ~Reader() {} + virtual void settings(int, char**) {}; + 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; virtual void read_atoms(int, int, double **) = 0; - void file(FILE *); + virtual void open_file(const char *); + virtual void close_file(); protected: - FILE *fp; // pointer to file opened by caller + FILE *fp; // pointer to opened file or pipe + int compressed; // flag for dump file compression }; } diff --git a/src/reader_native.cpp b/src/reader_native.cpp index 7c17588fea..85dc19dfe9 100644 --- a/src/reader_native.cpp +++ b/src/reader_native.cpp @@ -101,10 +101,10 @@ void ReaderNative::skip() ------------------------------------------------------------------------- */ bigint ReaderNative::read_header(double box[3][3], int &triclinic, - int fieldinfo, int nfield, - int *fieldtype, char **fieldlabel, - int scaledflag, int &fieldflag, - int &xflag, int &yflag, int &zflag) + int fieldinfo, int nfield, + int *fieldtype, char **fieldlabel, + int scaledflag, int &fieldflag, + int &xflag, int &yflag, int &zflag) { bigint natoms; read_lines(2); @@ -254,11 +254,11 @@ bigint ReaderNative::read_header(double box[3][3], int &triclinic, delete [] labels; - // set fieldflag = 1 if any unfound fields + // set fieldflag = -1 if any unfound fields fieldflag = 0; for (int i = 0; i < nfield; i++) - if (fieldindex[i] < 0) fieldflag = 1; + if (fieldindex[i] < 0) fieldflag = -1; // create internal vector of word ptrs for future parsing of per-atom lines @@ -279,7 +279,7 @@ void ReaderNative::read_atoms(int n, int nfield, double **fields) int i,m; char *eof; - for (int i = 0; i < n; i++) { + for (i = 0; i < n; i++) { eof = fgets(line,MAXLINE,fp); if (eof == NULL) error->one(FLERR,"Unexpected end of dump file"); diff --git a/src/reader_xyz.cpp b/src/reader_xyz.cpp new file mode 100644 index 0000000000..c2ba6dd155 --- /dev/null +++ b/src/reader_xyz.cpp @@ -0,0 +1,215 @@ +/* ---------------------------------------------------------------------- + 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: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "stdlib.h" +#include "reader_xyz.h" +#include "atom.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define MAXLINE 1024 // max line length in dump file + +enum{ID,TYPE,X,Y,Z}; + +/* ---------------------------------------------------------------------- */ + +ReaderXYZ::ReaderXYZ(LAMMPS *lmp) : Reader(lmp) +{ + line = new char[MAXLINE]; + fieldindex = NULL; + nstep = 0; +} + +/* ---------------------------------------------------------------------- */ + +ReaderXYZ::~ReaderXYZ() +{ + delete [] line; + memory->destroy(fieldindex); +} + +/* ---------------------------------------------------------------------- + read and return time stamp from dump file + if first read reaches end-of-file, return 1 so caller can open next file + only called by proc 0 +------------------------------------------------------------------------- */ + +int ReaderXYZ::read_time(bigint &ntimestep) +{ + char *eof = fgets(line,MAXLINE,fp); + if (eof == NULL) return 1; + + // first line has to have the number of atoms + + natoms = ATOBIGINT(line); + if (natoms < 1) + error->one(FLERR,"Dump file is incorrectly formatted"); + + // skip over comment/title line + + read_lines(1); + + // fake time step numbers + + ntimestep = nstep; + + // count this frame + + ++nstep; + return 0; +} + +/* ---------------------------------------------------------------------- + skip snapshot from timestamp onward + only called by proc 0 +------------------------------------------------------------------------- */ + +void ReaderXYZ::skip() +{ + // invoke read_lines() in chunks no larger than MAXSMALLINT + + int nchunk; + bigint nremain = natoms; + while (nremain) { + nchunk = MIN(nremain,MAXSMALLINT); + read_lines(nchunk); + nremain -= nchunk; + } +} + +/* ---------------------------------------------------------------------- + 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 + 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 + only called by proc 0 +------------------------------------------------------------------------- */ + +bigint ReaderXYZ::read_header(double box[3][3], int &triclinic, + int fieldinfo, int nfield, + int *fieldtype, char **fieldlabel, + int scaledflag, int &fieldflag, + int &xflag, int &yflag, int &zflag) +{ + + nid = 0; + + // signal that we have no box info at all + + triclinic = -1; + + // if no field info requested, just return + + if (!fieldinfo) return natoms; + + 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. + + xflag = scaledflag; + yflag = scaledflag; + zflag = scaledflag; + + // copy fieldtype list for supported fields + + fieldflag = 0; + for (int i = 0; i < nfield; i++) { + if ( (fieldtype[i] == X) || + (fieldtype[i] == Y) || + (fieldtype[i] == Z) || + (fieldtype[i] == ID) || + (fieldtype[i] == TYPE) ) { + fieldindex[i] = fieldtype[i]; + } else { + fieldflag = 1; + } + } + + return natoms; +} + +/* ---------------------------------------------------------------------- + read N atom lines from dump file + stores appropriate values in fields array + return 0 if success, 1 if error + only called by proc 0 +------------------------------------------------------------------------- */ + +void ReaderXYZ::read_atoms(int n, int nfield, double **fields) +{ + int i,m; + char *eof; + int mytype; + double myx, myy, myz; + + for (i = 0; i < n; i++) { + eof = fgets(line,MAXLINE,fp); + if (eof == NULL) error->one(FLERR,"Unexpected end of dump file"); + + ++nid; + sscanf(line,"%*s%lg%lg%lg", &myx, &myy, &myz); + + // XXX: we could insert an element2type translation here + // XXX: for now we flag unrecognized types as type 0, + // XXX: which should trigger an error, if LAMMPS uses it. + mytype = atoi(line); + + for (m = 0; m < nfield; m++) { + switch (fieldindex[m]) { + case X: + fields[i][m] = myx; + break; + case Y: + fields[i][m] = myy; + break; + case Z: + fields[i][m] = myz; + break; + case ID: + fields[i][m] = nid; + break; + case TYPE: + fields[i][m] = mytype; + break; + } + } + } +} + +/* ---------------------------------------------------------------------- + read N lines from dump file + only last one is saved in line + return NULL if end-of-file error, else non-NULL + only called by proc 0 +------------------------------------------------------------------------- */ + +void ReaderXYZ::read_lines(int n) +{ + char *eof; + for (int i = 0; i < n; i++) eof = fgets(line,MAXLINE,fp); + if (eof == NULL) error->all(FLERR,"Unexpected end of dump file"); +} diff --git a/src/reader_xyz.h b/src/reader_xyz.h new file mode 100644 index 0000000000..cdd760b5cc --- /dev/null +++ b/src/reader_xyz.h @@ -0,0 +1,54 @@ +/* -*- 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 Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#ifdef READER_CLASS + +ReaderStyle(xyz,ReaderXYZ) + +#else + +#ifndef LMP_READER_XYZ_H +#define LMP_READER_XYZ_H + +#include "reader.h" + +namespace LAMMPS_NS { + +class ReaderXYZ : public Reader { + public: + ReaderXYZ(class LAMMPS *); + ~ReaderXYZ(); + + int read_time(bigint &); + void skip(); + bigint read_header(double [3][3], int &, int, int, int *, char **, + int, int &, int &, int &, int &); + void read_atoms(int, int, double **); + +private: + char *line; // line read from dump file + bigint nstep; // current (time) step number + bigint natoms; // current number of atoms + bigint nid; // current atom id. + + int *fieldindex; // mapping of input fields to dump + + void read_lines(int); +}; + +} + +#endif +#endif diff --git a/src/rerun.cpp b/src/rerun.cpp index ef0acdb81e..80baf123b2 100644 --- a/src/rerun.cpp +++ b/src/rerun.cpp @@ -105,18 +105,23 @@ void Rerun::command(int narg, char **arg) } else error->all(FLERR,"Illegal rerun command"); } - if (strcmp(arg[iarg],"dump") != 0) 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 + // pass list of filenames to ReadDump + // along with post-"dump" args and post-"format" args ReadDump *rd = new ReadDump(lmp); rd->store_files(nfile,arg); - rd->fields_and_keywords(narg-iarg-1,&arg[iarg+1]); - rd->setup_reader(); + if (nremain) + nremain = rd->fields_and_keywords(narg-nremain,&arg[narg-nremain]); + else nremain = rd->fields_and_keywords(0,NULL); + if (nremain) rd->setup_reader(narg-nremain,&arg[narg-nremain]); + else rd->setup_reader(0,NULL); // perform the psuedo run // invoke lmp->init() only once