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

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->open_file(files[ifile]);
-      reader->file(fp);
       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]);
+      if (ifile != currentfile) reader->open_file(files[ifile]);
-      reader->file(fp);
       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) ||
+  if (triclinic_snap < 0 && boxflag > 0)
-      (!triclinic_snap && triclinic))
+    error->all(FLERR,"No box information in dump. You have to use 'box no'");
-    error->one(FLERR,"Read_dump triclinic status does not match simulation");
+  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 one[3];
-  double *coord;
+
-  
   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;
-}

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();
 };
 
 }

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
+   try to open given file
-   caller opens/closes dump files
+   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;
 }

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
 };
 
 }

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 fieldinfo, int nfield,
-                                   int *fieldtype, char **fieldlabel,
+                                 int *fieldtype, char **fieldlabel,
-                                   int scaledflag, int &fieldflag,
+                                 int scaledflag, int &fieldflag,
-                                   int &xflag, int &yflag, int &zflag)
+                                 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");
 

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");
+}

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

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]);
+  if (nremain) 
-  rd->setup_reader();
+    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