Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/src/reader_native.cpp
f16f9393ef09lammm-devel

Last Change
Raw File

reader_native.cpp
View Options

/* ----------------------------------------------------------------------
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 <string.h>
#include <stdlib.h>
#include "reader_native.h"
#include "atom.h"
#include "memory.h"
#include "error.h"
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,Q,IX,IY,IZ,FX,FY,FZ};
enum{UNSET,NOSCALE_NOWRAP,NOSCALE_WRAP,SCALE_NOWRAP,SCALE_WRAP};
/* ---------------------------------------------------------------------- */
ReaderNative::ReaderNative(LAMMPS *lmp) : Reader(lmp)
{
line = new char[MAXLINE];
words = NULL;
fieldindex = NULL;
}
/* ---------------------------------------------------------------------- */
ReaderNative::~ReaderNative()
{
delete [] line;
delete [] words;
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 ReaderNative::read_time(bigint &ntimestep)
{
char *eof = fgets(line,MAXLINE,fp);
if (eof == NULL) return 1;
if (strstr(line,"ITEM: TIMESTEP") != line)
error->one(FLERR,"Dump file is incorrectly formatted");
read_lines(1);
sscanf(line,BIGINT_FORMAT,&ntimestep);
return 0;
}
/* ----------------------------------------------------------------------
skip snapshot from timestamp onward
only called by proc 0
------------------------------------------------------------------------- */
void ReaderNative::skip()
{
read_lines(2);
bigint natoms;
sscanf(line,BIGINT_FORMAT,&natoms);
read_lines(5);
// 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 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 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 scaleflag, int wrapflag, int &fieldflag,
int &xflag, int &yflag, int &zflag)
{
bigint natoms;
read_lines(2);
sscanf(line,BIGINT_FORMAT,&natoms);
triclinic = 0;
box[0][2] = box[1][2] = box[2][2] = 0.0;
read_lines(1);
if (line[strlen("ITEM: BOX BOUNDS ")] == 'x') triclinic = 1;
read_lines(1);
if (!triclinic) sscanf(line,"%lg %lg",&box[0][0],&box[0][1]);
else sscanf(line,"%lg %lg %lg",&box[0][0],&box[0][1],&box[0][2]);
read_lines(1);
if (!triclinic) sscanf(line,"%lg %lg",&box[1][0],&box[1][1]);
else sscanf(line,"%lg %lg %lg",&box[1][0],&box[1][1],&box[1][2]);
read_lines(1);
if (!triclinic) sscanf(line,"%lg %lg",&box[2][0],&box[2][1]);
else sscanf(line,"%lg %lg %lg",&box[2][0],&box[2][1],&box[2][2]);
read_lines(1);
// if no field info requested, just return
if (!fieldinfo) return natoms;
// exatract column labels and match to requested fields
char *labelline = &line[strlen("ITEM: ATOMS ")];
nwords = atom->count_words(labelline);
char **labels = new char*[nwords];
labels[0] = strtok(labelline," \t\n\r\f");
if (labels[0] == NULL) {
delete[] labels;
return 1;
}
for (int m = 1; m < nwords; m++) {
labels[m] = strtok(NULL," \t\n\r\f");
if (labels[m] == NULL) {
delete[] labels;
return 1;
}
}
// 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 scaleflag + wrapflag (if fieldlabel set) or column label
memory->create(fieldindex,nfield,"read_dump:fieldindex");
int s_index,u_index,su_index;
xflag = UNSET;
yflag = UNSET;
zflag = UNSET;
for (int i = 0; i < nfield; i++) {
if (fieldlabel[i]) {
fieldindex[i] = find_label(fieldlabel[i],nwords,labels);
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)
fieldindex[i] = find_label("id",nwords,labels);
else if (fieldtype[i] == TYPE)
fieldindex[i] = find_label("type",nwords,labels);
else if (fieldtype[i] == X) {
fieldindex[i] = find_label("x",nwords,labels);
xflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("xs",nwords,labels);
u_index = find_label("xu",nwords,labels);
su_index = find_label("xsu",nwords,labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
xflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
xflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
xflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords) fieldindex[i] = -1;
} else if (fieldtype[i] == Y) {
fieldindex[i] = find_label("y",nwords,labels);
yflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("ys",nwords,labels);
u_index = find_label("yu",nwords,labels);
su_index = find_label("ysu",nwords,labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
yflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
yflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
yflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords) fieldindex[i] = -1;
} else if (fieldtype[i] == Z) {
fieldindex[i] = find_label("z",nwords,labels);
zflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("zs",nwords,labels);
u_index = find_label("zu",nwords,labels);
su_index = find_label("zsu",nwords,labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
zflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
zflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
zflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords) fieldindex[i] = -1;
} else if (fieldtype[i] == VX)
fieldindex[i] = find_label("vx",nwords,labels);
else if (fieldtype[i] == VY)
fieldindex[i] = find_label("vy",nwords,labels);
else if (fieldtype[i] == VZ)
fieldindex[i] = find_label("vz",nwords,labels);
else if (fieldtype[i] == FX)
fieldindex[i] = find_label("fx",nwords,labels);
else if (fieldtype[i] == FY)
fieldindex[i] = find_label("fy",nwords,labels);
else if (fieldtype[i] == FZ)
fieldindex[i] = find_label("fz",nwords,labels);
else if (fieldtype[i] == Q)
fieldindex[i] = find_label("q",nwords,labels);
else if (fieldtype[i] == IX)
fieldindex[i] = find_label("ix",nwords,labels);
else if (fieldtype[i] == IY)
fieldindex[i] = find_label("iy",nwords,labels);
else if (fieldtype[i] == IZ)
fieldindex[i] = find_label("iz",nwords,labels);
}
delete [] labels;
// set fieldflag = -1 if any unfound fields
fieldflag = 0;
for (int i = 0; i < nfield; i++)
if (fieldindex[i] < 0) fieldflag = -1;
// create internal vector of word ptrs for future parsing of per-atom lines
words = new char*[nwords];
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 ReaderNative::read_atoms(int n, int nfield, double **fields)
{
int i,m;
char *eof;
for (i = 0; i < n; i++) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of dump file");
// tokenize the line
words[0] = strtok(line," \t\n\r\f");
for (m = 1; m < nwords; m++)
words[m] = strtok(NULL," \t\n\r\f");
// convert selected fields to floats
for (m = 0; m < nfield; m++)
fields[i][m] = atof(words[fieldindex[m]]);
}
}
/* ----------------------------------------------------------------------
match label to any of N labels
return index of match or -1 if no match
------------------------------------------------------------------------- */
int ReaderNative::find_label(const char *label, int n, char **labels)
{
for (int i = 0; i < n; i++)
if (strcmp(label,labels[i]) == 0) return i;
return -1;
}
/* ----------------------------------------------------------------------
read N lines from dump file
only last one is saved in line
only called by proc 0
------------------------------------------------------------------------- */
void ReaderNative::read_lines(int n)
{
char *eof = NULL;
if (n <= 0) return;
for (int i = 0; i < n; i++) eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of dump file");
}
c4science · Help