diff --git a/src/fix_store.cpp b/src/fix_store.cpp index 9da5f6870d..45eeaaa502 100644 --- a/src/fix_store.cpp +++ b/src/fix_store.cpp @@ -29,6 +29,7 @@ FixStore::FixStore(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (narg != 5) error->all(FLERR,"Illegal fix store command"); // syntax: id group style 0/1 nvalue + // 0/1 flag = not-store or store values in restart file restart_peratom = force->inumeric(FLERR,arg[3]); nvalues = force->inumeric(FLERR,arg[4]); @@ -46,7 +47,7 @@ FixStore::FixStore(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) if (restart_peratom) atom->add_callback(1); // zero the storage - // since may be exchanged before filled by caller, e.g. fix store/force + // since may be exchanged before filled by caller int nlocal = atom->nlocal; if (vecflag) diff --git a/src/lammps.cpp b/src/lammps.cpp index 8a43b4f9aa..5d138617b0 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -537,6 +537,9 @@ void LAMMPS::destroy() delete atom; // atom must come after modify, neighbor // since fixes delete callbacks in atom delete timer; + + modify = NULL; // necessary since input->variable->varreader + // will be destructed later } /* ---------------------------------------------------------------------- diff --git a/src/read_data.cpp b/src/read_data.cpp index a45d560e83..9d96ccaea1 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -941,7 +941,7 @@ void ReadData::mass() next = strchr(buf,'\n'); *next = '\0'; atom->set_mass(buf); - buf += strlen(buf) + 1; + buf = next + 1; } delete [] original; } @@ -950,7 +950,6 @@ void ReadData::mass() void ReadData::paircoeffs() { - int i,m; char *next; char *buf = new char[atom->ntypes*MAXLINE]; @@ -958,13 +957,12 @@ void ReadData::paircoeffs() if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (i = 0; i < atom->ntypes; i++) { + for (int i = 0; i < atom->ntypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - m = strlen(buf) + 1; parse_coeffs(buf,NULL,1); force->pair->coeff(narg,arg); - buf += m; + buf = next + 1; } delete [] original; } @@ -973,7 +971,7 @@ void ReadData::paircoeffs() void ReadData::pairIJcoeffs() { - int i,j,m; + int i,j; char *next; int nsq = atom->ntypes* (atom->ntypes+1) / 2; @@ -987,10 +985,9 @@ void ReadData::pairIJcoeffs() for (j = i; j < atom->ntypes; j++) { next = strchr(buf,'\n'); *next = '\0'; - m = strlen(buf) + 1; parse_coeffs(buf,NULL,0); force->pair->coeff(narg,arg); - buf += m; + buf = next + 1; } delete [] original; } @@ -999,7 +996,6 @@ void ReadData::pairIJcoeffs() void ReadData::bondcoeffs() { - int i,m; char *next; char *buf = new char[atom->nbondtypes*MAXLINE]; @@ -1007,13 +1003,12 @@ void ReadData::bondcoeffs() if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (i = 0; i < atom->nbondtypes; i++) { + for (int i = 0; i < atom->nbondtypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - m = strlen(buf) + 1; parse_coeffs(buf,NULL,0); force->bond->coeff(narg,arg); - buf += m; + buf = next + 1; } delete [] original; } @@ -1022,7 +1017,6 @@ void ReadData::bondcoeffs() void ReadData::anglecoeffs(int which) { - int i,m; char *next; char *buf = new char[atom->nangletypes*MAXLINE]; @@ -1030,15 +1024,14 @@ void ReadData::anglecoeffs(int which) if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (i = 0; i < atom->nangletypes; i++) { + for (int i = 0; i < atom->nangletypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - m = strlen(buf) + 1; if (which == 0) parse_coeffs(buf,NULL,0); else if (which == 1) parse_coeffs(buf,"bb",0); else if (which == 2) parse_coeffs(buf,"ba",0); force->angle->coeff(narg,arg); - buf += m; + buf = next + 1; } delete [] original; } @@ -1047,7 +1040,6 @@ void ReadData::anglecoeffs(int which) void ReadData::dihedralcoeffs(int which) { - int i,m; char *next; char *buf = new char[atom->ndihedraltypes*MAXLINE]; @@ -1055,10 +1047,9 @@ void ReadData::dihedralcoeffs(int which) if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (i = 0; i < atom->ndihedraltypes; i++) { + for (int i = 0; i < atom->ndihedraltypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - m = strlen(buf) + 1; if (which == 0) parse_coeffs(buf,NULL,0); else if (which == 1) parse_coeffs(buf,"mbt",0); else if (which == 2) parse_coeffs(buf,"ebt",0); @@ -1066,7 +1057,7 @@ void ReadData::dihedralcoeffs(int which) else if (which == 4) parse_coeffs(buf,"aat",0); else if (which == 5) parse_coeffs(buf,"bb13",0); force->dihedral->coeff(narg,arg); - buf += m; + buf = next + 1; } delete [] original; } @@ -1075,7 +1066,6 @@ void ReadData::dihedralcoeffs(int which) void ReadData::impropercoeffs(int which) { - int i,m; char *next; char *buf = new char[atom->nimpropertypes*MAXLINE]; @@ -1083,14 +1073,13 @@ void ReadData::impropercoeffs(int which) if (eof) error->all(FLERR,"Unexpected end of data file"); char *original = buf; - for (i = 0; i < atom->nimpropertypes; i++) { + for (int i = 0; i < atom->nimpropertypes; i++) { next = strchr(buf,'\n'); *next = '\0'; - m = strlen(buf) + 1; if (which == 0) parse_coeffs(buf,NULL,0); else if (which == 1) parse_coeffs(buf,"aa",0); force->improper->coeff(narg,arg); - buf += m; + buf = next + 1; } delete [] original; } @@ -1102,7 +1091,7 @@ void ReadData::impropercoeffs(int which) void ReadData::fix(int ifix, char *line, bigint nlines) { - int i,m,nchunk,eof; + int nchunk,eof; bigint nread = 0; while (nread < nlines) { diff --git a/src/variable.cpp b/src/variable.cpp index 406da2474a..2d94d08d27 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -22,10 +22,12 @@ #include "update.h" #include "group.h" #include "domain.h" +#include "comm.h" #include "region.h" #include "modify.h" #include "compute.h" #include "fix.h" +#include "fix_store.h" #include "output.h" #include "thermo.h" #include "random_mars.h" @@ -41,14 +43,16 @@ using namespace MathConst; #define VARDELTA 4 #define MAXLEVEL 4 #define MAXLINE 256 +#define CHUNK 1024 #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a) -enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,FILEVAR,AFILEVAR,EQUAL,ATOM}; +enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,SCALARFILE,ATOMFILE,EQUAL,ATOM}; enum{ARG,OP}; // customize by adding a function -// if add before OR, also set precedence level and precedence length in *.h +// if add before OR, +// also set precedence level in constructor and precedence length in *.h enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY, NOT,EQ,NE,LT,LE,GT,GE,AND,OR, @@ -276,39 +280,41 @@ void Variable::set(int narg, char **arg) data[nvar] = new char*[num[nvar]]; copy(1,&arg[2],data[nvar]); - // FILEVAR for strings or numbers + // SCALARFILE for strings or numbers // which = 1st value + // data = 1 value, string to eval } else if (strcmp(arg[1],"file") == 0) { if (narg != 3) error->all(FLERR,"Illegal variable command"); if (find(arg[0]) >= 0) return; if (nvar == maxvar) grow(); - style[nvar] = FILEVAR; + style[nvar] = SCALARFILE; num[nvar] = 1; which[nvar] = 0; pad[nvar] = 0; data[nvar] = new char*[num[nvar]]; data[nvar][0] = new char[MAXLINE]; - reader[nvar] = new VarReader(lmp,arg[2],FILEVAR); - int flag = reader[nvar]->read(data[nvar][0]); + reader[nvar] = new VarReader(lmp,arg[0],arg[2],SCALARFILE); + int flag = reader[nvar]->read_scalar(data[nvar][0]); if (flag) error->all(FLERR,"File variable could not read value"); - // AFILEVAR for numbers + // ATOMFILE for numbers // which = 1st value + // data = NULL - } else if (strcmp(arg[1],"afile") == 0) { + } else if (strcmp(arg[1],"atomfile") == 0) { if (narg != 3) error->all(FLERR,"Illegal variable command"); if (find(arg[0]) >= 0) return; if (nvar == maxvar) grow(); - style[nvar] = AFILEVAR; + style[nvar] = ATOMFILE; num[nvar] = 1; which[nvar] = 0; pad[nvar] = 0; - //data[nvar] = new char*[num[nvar]]; - //data[nvar][0] = new char[MAXLINE]; - reader[nvar] = new VarReader(lmp,arg[2],AFILEVAR); - int flag = reader[nvar]->read(data[nvar][0]); - if (flag) error->all(FLERR,"Afile variable could not read values"); + data[nvar] = new char*[num[nvar]]; + data[nvar][0] = NULL; + reader[nvar] = new VarReader(lmp,arg[0],arg[2],ATOMFILE); + int flag = reader[nvar]->read_peratom(); + if (flag) error->all(FLERR,"Atomfile variable could not read values"); // EQUAL // remove pre-existing var if also style EQUAL (allows it to be reset) @@ -428,11 +434,22 @@ int Variable::next(int narg, char **arg) } } - } else if (istyle == FILEVAR) { + } else if (istyle == SCALARFILE) { for (int iarg = 0; iarg < narg; iarg++) { ivar = find(arg[iarg]); - int done = reader[ivar]->read(data[ivar][0]); + int done = reader[ivar]->read_scalar(data[ivar][0]); + if (done) { + flag = 1; + remove(ivar); + } + } + + } else if (istyle == ATOMFILE) { + + for (int iarg = 0; iarg < narg; iarg++) { + ivar = find(arg[iarg]); + int done = reader[ivar]->read_peratom(); if (done) { flag = 1; remove(ivar); @@ -484,11 +501,13 @@ int Variable::next(int narg, char **arg) /* ---------------------------------------------------------------------- return ptr to the data text associated with a variable - if INDEX or WORLD or UNIVERSE or STRING var, return ptr to stored string + if INDEX or WORLD or UNIVERSE or STRING or SCALARFILE var, + return ptr to stored string if LOOP or ULOOP var, write int to data[0] and return ptr to string if EQUAL var, evaluate variable and put result in str - if ATOM var, return NULL - return NULL if no variable or which is bad, caller must respond + if ATOM or ATOMFILE var, return NULL + return NULL if no variable with name or which value is bad, + caller must respond ------------------------------------------------------------------------- */ char *Variable::retrieve(char *name) @@ -500,7 +519,7 @@ char *Variable::retrieve(char *name) char *str; if (style[ivar] == INDEX || style[ivar] == WORLD || style[ivar] == UNIVERSE || style[ivar] == STRING || - style[ivar] == FILEVAR) { + style[ivar] == SCALARFILE) { str = data[ivar][which[ivar]]; } else if (style[ivar] == LOOP || style[ivar] == ULOOP) { char result[16]; @@ -524,7 +543,7 @@ char *Variable::retrieve(char *name) data[ivar][1] = new char[n]; strcpy(data[ivar][1],result); str = data[ivar][1]; - } else if (style[ivar] == ATOM) return NULL; + } else if (style[ivar] == ATOM || style[ivar] == ATOMFILE) return NULL; return str; } @@ -555,7 +574,7 @@ double Variable::compute_equal(char *str) } /* ---------------------------------------------------------------------- - compute result of atom-style variable evaluation + compute result of atom-style and atomfile-atyle variable evaluation only computed for atoms in igroup, else result is 0.0 answers are placed every stride locations into result if sumflag, add variable values to existing result @@ -565,30 +584,53 @@ void Variable::compute_atom(int ivar, int igroup, double *result, int stride, int sumflag) { Tree *tree; - double tmp = evaluate(data[ivar][0],&tree); - tmp = collapse_tree(tree); + double *vstore; + + if (style[ivar] == ATOM) { + double tmp = evaluate(data[ivar][0],&tree); + tmp = collapse_tree(tree); + } else vstore = reader[ivar]->fix->vstore; int groupbit = group->bitmask[igroup]; int *mask = atom->mask; int nlocal = atom->nlocal; - if (sumflag == 0) { - int m = 0; - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) result[m] = eval_tree(tree,i); - else result[m] = 0.0; - m += stride; + if (style[ivar] == ATOM) { + if (sumflag == 0) { + int m = 0; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) result[m] = eval_tree(tree,i); + else result[m] = 0.0; + m += stride; + } + + } else { + int m = 0; + for (int i = 0; i < nlocal; i++) { + if (mask[i] && groupbit) result[m] += eval_tree(tree,i); + m += stride; + } } } else { - int m = 0; - for (int i = 0; i < nlocal; i++) { - if (mask[i] && groupbit) result[m] += eval_tree(tree,i); - m += stride; + if (sumflag == 0) { + int m = 0; + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) result[m] = vstore[i]; + else result[m] = 0.0; + m += stride; + } + + } else { + int m = 0; + for (int i = 0; i < nlocal; i++) { + if (mask[i] && groupbit) result[m] += vstore[i]; + m += stride; + } } } - free_tree(tree); + if (style[ivar] == ATOM) free_tree(tree); } /* ---------------------------------------------------------------------- @@ -614,17 +656,18 @@ int Variable::equalstyle(int ivar) } /* ---------------------------------------------------------------------- - return 1 if variable is ATOM style, 0 if not + return 1 if variable is ATOM or ATOMFILE style, 0 if not ------------------------------------------------------------------------- */ int Variable::atomstyle(int ivar) { - if (style[ivar] == ATOM) return 1; + if (style[ivar] == ATOM || style[ivar] == ATOMFILE) return 1; return 0; } /* ---------------------------------------------------------------------- remove Nth variable from list and compact list + delete reader explicitly if it exists ------------------------------------------------------------------------- */ void Variable::remove(int n) @@ -633,6 +676,7 @@ void Variable::remove(int n) if (style[n] == LOOP || style[n] == ULOOP) delete [] data[n][0]; else for (int i = 0; i < num[n]; i++) delete [] data[n][i]; delete [] data[n]; + delete reader[n]; for (int i = n+1; i < nvar; i++) { names[i-1] = names[i]; @@ -983,6 +1027,7 @@ double Variable::evaluate(char *str, Tree **tree) newtree->type = ATOMARRAY; newtree->array = compute->vector_atom; newtree->nstride = 1; + newtree->selfalloc = 0; newtree->left = newtree->middle = newtree->right = NULL; treestack[ntreestack++] = newtree; @@ -1013,6 +1058,7 @@ double Variable::evaluate(char *str, Tree **tree) else newtree->array = NULL; newtree->nstride = compute->size_peratom_cols; + newtree->selfalloc = 0; newtree->left = newtree->middle = newtree->right = NULL; treestack[ntreestack++] = newtree; @@ -1158,6 +1204,7 @@ double Variable::evaluate(char *str, Tree **tree) newtree->type = ATOMARRAY; newtree->array = fix->vector_atom; newtree->nstride = 1; + newtree->selfalloc = 0; newtree->left = newtree->middle = newtree->right = NULL; treestack[ntreestack++] = newtree; @@ -1182,6 +1229,7 @@ double Variable::evaluate(char *str, Tree **tree) else newtree->array = NULL; newtree->nstride = fix->size_peratom_cols; + newtree->selfalloc = 0; newtree->left = newtree->middle = newtree->right = NULL; treestack[ntreestack++] = newtree; @@ -1216,9 +1264,9 @@ double Variable::evaluate(char *str, Tree **tree) i = ptr-str+1; } - // v_name = scalar from non atom-style global scalar + // v_name = scalar from non atom/atomfile variable - if (nbracket == 0 && style[ivar] != ATOM) { + if (nbracket == 0 && style[ivar] != ATOM && style[ivar] != ATOMFILE) { char *var = retrieve(id); if (var == NULL) @@ -1231,7 +1279,8 @@ double Variable::evaluate(char *str, Tree **tree) treestack[ntreestack++] = newtree; } else argstack[nargstack++] = atof(var); - // v_name = vector from atom-style per-atom vector + // v_name = per-atom vector from atom-style variable + // evaluate the atom-style variable as newtree } else if (nbracket == 0 && style[ivar] == ATOM) { @@ -1242,7 +1291,22 @@ double Variable::evaluate(char *str, Tree **tree) evaluate(data[ivar][0],&newtree); treestack[ntreestack++] = newtree; - // v_name[N] = scalar from atom-style per-atom vector + // v_name = per-atom vector from atomfile-style variable + + } else if (nbracket == 0 && style[ivar] == ATOMFILE) { + + if (tree == NULL) + error->all(FLERR,"Atomfile-style variable in " + "equal-style variable formula"); + Tree *newtree = new Tree(); + newtree->type = ATOMARRAY; + newtree->array = reader[ivar]->fix->vstore; + newtree->nstride = 1; + newtree->selfalloc = 0; + newtree->left = newtree->middle = newtree->right = NULL; + treestack[ntreestack++] = newtree; + + // v_name[N] = scalar from atom-style variable // compute the per-atom variable in result // use peratom2global to extract single value from result @@ -1255,6 +1319,13 @@ double Variable::evaluate(char *str, Tree **tree) tree,treestack,ntreestack,argstack,nargstack); memory->destroy(result); + // v_name[N] = scalar from atomfile-style variable + + } else if (nbracket && style[ivar] == ATOMFILE) { + + peratom2global(1,NULL,reader[ivar]->fix->vstore,1,index, + tree,treestack,ntreestack,argstack,nargstack); + } else error->all(FLERR,"Mismatched variable in variable formula"); delete [] id; @@ -2234,6 +2305,10 @@ void Variable::free_tree(Tree *tree) if (tree->left) free_tree(tree->left); if (tree->middle) free_tree(tree->middle); if (tree->right) free_tree(tree->right); + + if (tree->type == ATOMARRAY && tree->selfalloc) + memory->destroy(tree->array); + delete tree; } @@ -2925,7 +3000,7 @@ int Variable::region_function(char *id) contents = str between parentheses with one,two,three args return 0 if not a match, 1 if successfully processed customize by adding a special function: - sum(x),min(x),max(x),ave(x),trap(x),gmask(x),rmask(x),grmask(x,y) + sum(x),min(x),max(x),ave(x),trap(x),gmask(x),rmask(x),grmask(x,y),next(x) ------------------------------------------------------------------------- */ int Variable::special_function(char *word, char *contents, Tree **tree, @@ -3159,7 +3234,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, newtree->left = newtree->middle = newtree->right = NULL; treestack[ntreestack++] = newtree; - // file variable special function + // special function for file-style or atomfile-stlye variables } else if (strcmp(word,"next") == 0) { if (narg != 1) @@ -3168,21 +3243,49 @@ int Variable::special_function(char *word, char *contents, Tree **tree, int ivar = find(arg1); if (ivar == -1) error->all(FLERR,"Variable ID in variable formula does not exist"); - if (style[ivar] != FILEVAR) - error->all(FLERR,"Invalid variable in special function next"); - - double value = atof(data[ivar][0]); - reader[ivar]->read(data[ivar][0]); + // SCALARFILE has single current value, read next one // save value in tree or on argstack - if (tree) { + if (style[ivar] == SCALARFILE) { + double value = atof(data[ivar][0]); + int done = reader[ivar]->read_scalar(data[ivar][0]); + if (done) remove(ivar); + + if (tree) { + Tree *newtree = new Tree(); + newtree->type = VALUE; + newtree->value = value; + newtree->left = newtree->middle = newtree->right = NULL; + treestack[ntreestack++] = newtree; + } else argstack[nargstack++] = value; + + // ATOMFILE has one value per atom, only valid for + // save values in tree + + } else if (style[ivar] == ATOMFILE) { + if (tree == NULL) + error->all(FLERR,"Atomfile variable in equal-style variable formula"); + + // copy current per-atom values into result so can read next ones + // set selfalloc = 1 so will be deleted by free_tree() after eval + + double *result; + memory->create(result,atom->nlocal,"variable:result"); + memcpy(result,reader[ivar]->fix->vstore,atom->nlocal*sizeof(double)); + + int done = reader[ivar]->read_peratom(); + if (done) remove(ivar); + Tree *newtree = new Tree(); - newtree->type = VALUE; - newtree->value = value; + newtree->type = ATOMARRAY; + newtree->array = result; + newtree->nstride = 1; + newtree->selfalloc = 1; newtree->left = newtree->middle = newtree->right = NULL; treestack[ntreestack++] = newtree; - } else argstack[nargstack++] = value; + + } else error->all(FLERR,"Invalid variable style in special function next"); } delete [] arg1; @@ -3290,6 +3393,7 @@ void Variable::atom_vector(char *word, Tree **tree, Tree *newtree = new Tree(); newtree->type = ATOMARRAY; newtree->nstride = 3; + newtree->selfalloc = 0; newtree->left = newtree->middle = newtree->right = NULL; treestack[ntreestack++] = newtree; @@ -3664,12 +3768,15 @@ unsigned int Variable::data_mask(char *str) } /* ---------------------------------------------------------------------- - class to read variable values from a file, line by line + class to read variable values from a file + for flag = SCALARFILE, reads one value per line + for flag = ATOMFILE, reads set of one value per atom ------------------------------------------------------------------------- */ -VarReader::VarReader(LAMMPS *lmp, char *file, int flag) : Pointers(lmp) +VarReader::VarReader(LAMMPS *lmp, char *name, char *file, int flag) : + Pointers(lmp) { - MPI_Comm_rank(world,&me); + me = comm->me; style = flag; if (me == 0) { @@ -3679,6 +3786,37 @@ VarReader::VarReader(LAMMPS *lmp, char *file, int flag) : Pointers(lmp) sprintf(str,"Cannot open file variable file %s",file); error->one(FLERR,str); } + } else fp = NULL; + + // if atomfile-style variable, must store per-atom values read from file + // allocate a new fix STORE, so they persist + // id = variable-ID + VARIABLE_STORE, fix group = all + + fix = NULL; + id_fix = NULL; + buffer = NULL; + + if (style == ATOMFILE) { + if (atom->map_style == 0) + error->all(FLERR, + "Cannot use atomfile-style variable unless atom map exists"); + + int n = strlen(name) + strlen("_VARIABLE_STORE") + 1; + id_fix = new char[n]; + strcpy(id_fix,name); + strcat(id_fix,"_VARIABLE_STORE"); + + char **newarg = new char*[5]; + newarg[0] = id_fix; + newarg[1] = (char *) "all"; + newarg[2] = (char *) "STORE"; + newarg[3] = (char *) "0"; + newarg[4] = (char *) "1"; + modify->add_fix(5,newarg); + fix = (FixStore *) modify->fix[modify->nfix-1]; + delete [] newarg; + + buffer = new char[CHUNK*MAXLINE]; } } @@ -3687,16 +3825,24 @@ VarReader::VarReader(LAMMPS *lmp, char *file, int flag) : Pointers(lmp) VarReader::~VarReader() { if (me == 0) fclose(fp); + + // check modify in case all fixes have already been deleted + + if (fix) { + if (modify) modify->delete_fix(id_fix); + delete [] id_fix; + delete [] buffer; + } } /* ---------------------------------------------------------------------- - read for FILEVAR style - read next value from file into str + read for SCALARFILE style + read next value from file into str for file-style variable strip comments, skip blank lines return 0 if successful, 1 if end-of-file ------------------------------------------------------------------------- */ -int VarReader::read(char *str) +int VarReader::read_scalar(char *str) { int n; char *ptr; @@ -3721,10 +3867,49 @@ int VarReader::read(char *str) } /* ---------------------------------------------------------------------- - read for AFILEVAR style + read snapshot of per-atom values from file + into str for atomfile-style variable + return 0 if successful, 1 if end-of-file ------------------------------------------------------------------------- */ -int VarReader::read(double *) +int VarReader::read_peratom() { + int i,m,tagdata,nchunk,eof; + char *next; + double value; + + // set all per-atom values to 0.0 + // values that appear in file will overwrite this + + double *vstore = fix->vstore; + + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) vstore[i] = 0.0; + + int map_tag_max = atom->map_tag_max; + + // NOTE: read this value from file + bigint nlines = atom->natoms; + + bigint nread = 0; + while (nread < nlines) { + nchunk = MIN(nlines-nread,CHUNK); + eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer); + if (eof) return 1; + + char *buf = buffer; + for (i = 0; i < nchunk; i++) { + next = strchr(buf,'\n'); + *next = '\0'; + sscanf(buf,"%d %lg",&tagdata,&value); + if (tagdata <= 0 || tagdata > map_tag_max) + error->one(FLERR,"Invalid atom ID in variable file"); + if ((m = atom->map(tagdata)) >= 0) vstore[m] = value; + buf = next + 1; + } + + nread += nchunk; + } + return 0; } diff --git a/src/variable.h b/src/variable.h index 9a5071490c..0d70232233 100644 --- a/src/variable.h +++ b/src/variable.h @@ -41,7 +41,7 @@ class Variable : protected Pointers { private: int nvar; // # of defined variables - int maxvar; // max # of variables arrays can hold + int maxvar; // max # of variables following lists can hold char **names; // name of each variable int *style; // style of each variable int *num; // # of values for each variable @@ -60,13 +60,14 @@ class Variable : protected Pointers { int me; struct Tree { // parse tree for atom-style variables - double value; - double *array; - int *iarray; - int type; - int nstride; - int ivalue1,ivalue2; - Tree *left,*middle,*right; + double value; // single scalar + double *array; // per-atom or per-type list of doubles + int *iarray; // per-atom list of ints + int type; // operation, see enum{} in variable.cpp + int nstride; // stride between atoms if array is a 2d array + int selfalloc; // 1 if array is allocated here, else 0 + int ivalue1,ivalue2; // extra values for needed for gmask,rmask,grmask + Tree *left,*middle,*right; // ptrs further down tree }; void remove(int); @@ -96,14 +97,18 @@ class Variable : protected Pointers { class VarReader : protected Pointers { public: - VarReader(class LAMMPS *, char *, int); + class FixStore *fix; + char *id_fix; + + VarReader(class LAMMPS *, char *, char *, int); ~VarReader(); - int read(char *); - int read(double *); + int read_scalar(char *); + int read_peratom(); private: int me,style; FILE *fp; + char *buffer; }; }