diff --git a/src/compute_slice.cpp b/src/compute_slice.cpp index 700b2545c9..d2c3c18de4 100644 --- a/src/compute_slice.cpp +++ b/src/compute_slice.cpp @@ -18,13 +18,15 @@ #include "modify.h" #include "fix.h" #include "group.h" +#include "input.h" +#include "variable.h" #include "memory.h" #include "error.h" #include "force.h" using namespace LAMMPS_NS; -enum{COMPUTE,FIX}; +enum{COMPUTE,FIX,VARIABLE}; #define INVOKED_VECTOR 2 #define INVOKED_ARRAY 4 @@ -55,9 +57,11 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) : for (int iarg = 6; iarg < narg; iarg++) { if (strncmp(arg[iarg],"c_",2) == 0 || - strncmp(arg[iarg],"f_",2) == 0) { + strncmp(arg[iarg],"f_",2) == 0 || + strncmp(arg[iarg],"v_",2) == 0) { if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE; else if (arg[iarg][0] == 'f') which[nvalues] = FIX; + else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE; int n = strlen(arg[iarg]); char *suffix = new char[n]; @@ -89,36 +93,53 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Compute ID for compute slice does not exist"); if (modify->compute[icompute]->vector_flag) { if (argindex[i]) - error->all(FLERR,"Compute slice compute does not calculate a global array"); + error->all(FLERR,"Compute slice compute does not " + "calculate a global array"); if (nstop > modify->compute[icompute]->size_vector) - error->all(FLERR,"Compute slice compute vector is accessed out-of-range"); + error->all(FLERR,"Compute slice compute vector is " + "accessed out-of-range"); } else if (modify->compute[icompute]->array_flag) { if (argindex[i] == 0) - error->all(FLERR,"Compute slice compute does not calculate a global vector"); + error->all(FLERR,"Compute slice compute does not " + "calculate a global vector"); if (argindex[i] > modify->compute[icompute]->size_array_cols) - error->all(FLERR,"Compute slice compute array is accessed out-of-range"); + error->all(FLERR,"Compute slice compute array is " + "accessed out-of-range"); if (nstop > modify->compute[icompute]->size_array_rows) - error->all(FLERR,"Compute slice compute array is accessed out-of-range"); + error->all(FLERR,"Compute slice compute array is " + "accessed out-of-range"); } else error->all(FLERR,"Compute slice compute does not calculate " "global vector or array"); + } else if (which[i] == FIX) { int ifix = modify->find_fix(ids[i]); if (ifix < 0) error->all(FLERR,"Fix ID for compute slice does not exist"); if (modify->fix[ifix]->vector_flag) { if (argindex[i]) - error->all(FLERR,"Compute slice fix does not calculate a global array"); + error->all(FLERR,"Compute slice fix does not " + "calculate a global array"); if (nstop > modify->fix[ifix]->size_vector) error->all(FLERR,"Compute slice fix vector is accessed out-of-range"); } else if (modify->fix[ifix]->array_flag) { if (argindex[i] == 0) - error->all(FLERR,"Compute slice fix does not calculate a global vector"); + error->all(FLERR,"Compute slice fix does not " + "calculate a global vector"); if (argindex[i] > modify->fix[ifix]->size_array_cols) error->all(FLERR,"Compute slice fix array is accessed out-of-range"); if (nstop > modify->fix[ifix]->size_array_rows) error->all(FLERR,"Compute slice fix array is accessed out-of-range"); } else error->all(FLERR,"Compute slice fix does not calculate " "global vector or array"); + + } else if (which[i] == VARIABLE) { + int ivariable = input->variable->find(ids[i]); + if (ivariable < 0) + error->all(FLERR,"Variable name for compute slice does not exist"); + if (argindex[i] == 0 && input->variable->vectorstyle(ivariable) == 0) + error->all(FLERR,"Compute slice variable is not vector-style variable"); + if (argindex[i]) + error->all(FLERR,"Compute slice vector variable cannot be indexed"); } } @@ -157,6 +178,8 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) : extlist[j++] = modify->fix[ifix]->extlist[i-1]; } } else extvector = modify->fix[ifix]->extarray; + } else if (which[0] == VARIABLE) { + extvector = 0; } } else { @@ -189,6 +212,8 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) : } else { if (modify->fix[ifix]->extarray) extarray = 1; } + } else if (which[i] == VARIABLE) { + // variable is always intensive, does not change extarray } } } @@ -225,6 +250,11 @@ void ComputeSlice::init() if (ifix < 0) error->all(FLERR,"Fix ID for compute slice does not exist"); value2index[m] = ifix; + } else if (which[m] == VARIABLE) { + int ivariable = input->variable->find(ids[m]); + if (ivariable < 0) + error->all(FLERR,"Variable name for compute slice does not exist"); + value2index[m] = ivariable; } } } @@ -292,7 +322,8 @@ void ComputeSlice::extract_one(int m, double *vec, int stride) } else if (which[m] == FIX) { if (update->ntimestep % modify->fix[value2index[m]]->global_freq) - error->all(FLERR,"Fix used in compute slice not computed at compatible time"); + error->all(FLERR,"Fix used in compute slice not " + "computed at compatible time"); Fix *fix = modify->fix[value2index[m]]; if (argindex[m] == 0) { @@ -309,5 +340,18 @@ void ComputeSlice::extract_one(int m, double *vec, int stride) j += stride; } } + + // invoke vector-style variable + + } else if (which[m] == VARIABLE) { + double *varvec; + int nvec = input->variable->compute_vector(value2index[m],&varvec); + if (nvec < nstop) + error->all(FLERR,"Compute slice variable is not long enough"); + j = 0; + for (i = nstart; i < nstop; i += nskip) { + vec[j] = varvec[i-1]; + j += stride; + } } } diff --git a/src/domain.cpp b/src/domain.cpp index 217f30c4ed..a4a6a1feed 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -498,6 +498,17 @@ void Domain::pbc() int *mask = atom->mask; imageint *image = atom->image; + // verify owned atoms all have valid numerical coords + // may not if computed pairwise force between 2 atoms at same location + + int flag = 0; + for (i = 0; i < nlocal; i++) + if (!ISFINITE(x[i][0]) || !ISFINITE(x[i][1]) || !ISFINITE(x[i][2])) + flag = 1; + if (flag) error->one(FLERR,"Non-numeric atom coords - simulation unstable"); + + // setup for PBC checks + if (triclinic == 0) { lo = boxlo; hi = boxhi; @@ -508,6 +519,8 @@ void Domain::pbc() period = prd_lamda; } + // apply PBC to each owned atom + for (i = 0; i < nlocal; i++) { if (xperiodic) { if (x[i][0] < lo[0]) { diff --git a/src/fix_ave_correlate.cpp b/src/fix_ave_correlate.cpp index 336723943c..f2dac0ac4b 100644 --- a/src/fix_ave_correlate.cpp +++ b/src/fix_ave_correlate.cpp @@ -216,9 +216,12 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg): int ivariable = input->variable->find(ids[i]); if (ivariable < 0) error->all(FLERR,"Variable name for fix ave/correlate does not exist"); - if (input->variable->equalstyle(ivariable) == 0) + if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0) error->all(FLERR, "Fix ave/correlate variable is not equal-style variable"); + if (argindex[i] && input->variable->vectorstyle(ivariable) == 0) + error->all(FLERR, + "Fix ave/correlate variable is not vector-style variable"); } } @@ -443,10 +446,19 @@ void FixAveCorrelate::end_of_step() else scalar = modify->fix[m]->compute_vector(argindex[i]-1); - // evaluate equal-style variable + // evaluate equal-style or vector-style variable - } else if (which[i] == VARIABLE) - scalar = input->variable->compute_equal(m); + } else if (which[i] == VARIABLE) { + if (argindex[i] == 0) + scalar = input->variable->compute_equal(m); + else { + double *varvec; + int nvec = input->variable->compute_vector(m,&varvec); + int index = argindex[i]; + if (nvec < index) scalar = 0.0; + else scalar = varvec[index-1]; + } + } values[lastindex][i] = scalar; } diff --git a/src/fix_ave_histo.cpp b/src/fix_ave_histo.cpp index 190a68596a..902934be0f 100644 --- a/src/fix_ave_histo.cpp +++ b/src/fix_ave_histo.cpp @@ -422,15 +422,32 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Fix for fix ave/histo not computed at compatible time"); - } else if (which[i] == VARIABLE && kind == GLOBAL) { + } else if (which[i] == VARIABLE && kind == GLOBAL && mode == SCALAR) { int ivariable = input->variable->find(ids[i]); if (ivariable < 0) error->all(FLERR,"Variable name for fix ave/histo does not exist"); + if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0) + error->all(FLERR,"Fix ave/histo variable is not equal-style variable"); + if (argindex[i] && input->variable->vectorstyle(ivariable) == 0) + error->all(FLERR,"Fix ave/histo variable is not vector-style variable"); + + } else if (which[i] == VARIABLE && kind == GLOBAL && mode == VECTOR) { + int ivariable = input->variable->find(ids[i]); + if (ivariable < 0) + error->all(FLERR,"Variable name for fix ave/histo does not exist"); + if (argindex[i] == 0 && input->variable->vectorstyle(ivariable) == 0) + error->all(FLERR,"Fix ave/histo variable is not vector-style variable"); + if (argindex[i]) + error->all(FLERR,"Fix ave/histo variable cannot be indexed"); } else if (which[i] == VARIABLE && kind == PERATOM) { int ivariable = input->variable->find(ids[i]); if (ivariable < 0) error->all(FLERR,"Variable name for fix ave/histo does not exist"); + if (argindex[i] == 0 && input->variable->atomstyle(ivariable) == 0) + error->all(FLERR,"Fix ave/histo variable is not atom-style variable"); + if (argindex[i]) + error->all(FLERR,"Fix ave/histo variable cannot be indexed"); } } @@ -719,19 +736,32 @@ void FixAveHisto::end_of_step() fix->size_local_cols); } - // evaluate equal-style or atom-style variable + // evaluate equal-style or vector-style or atom-style variable - } else if (which[i] == VARIABLE && kind == GLOBAL) { - bin_one(input->variable->compute_equal(m)); + } else if (which[i] == VARIABLE) { + if (kind == GLOBAL && mode == SCALAR) { + if (j == 0) bin_one(input->variable->compute_equal(m)); + else { + double *varvec; + int nvec = input->variable->compute_vector(m,&varvec); + if (nvec < j) bin_one(0.0); + else bin_one(varvec[j-1]); + } - } else if (which[i] == VARIABLE && kind == PERATOM) { - if (atom->nlocal > maxatom) { - memory->destroy(vector); - maxatom = atom->nmax; - memory->create(vector,maxatom,"ave/histo:vector"); + } else if (kind == GLOBAL && mode == VECTOR) { + double *varvec; + int nvec = input->variable->compute_vector(m,&varvec); + bin_vector(nvec,varvec,1); + + } else if (which[i] == VARIABLE && kind == PERATOM) { + if (atom->nlocal > maxatom) { + memory->destroy(vector); + maxatom = atom->nmax; + memory->create(vector,maxatom,"ave/histo:vector"); + } + input->variable->compute_atom(m,igroup,vector,1,0); + bin_atoms(vector,1); } - input->variable->compute_atom(m,igroup,vector,1,0); - bin_atoms(vector,1); } } diff --git a/src/fix_ave_time.cpp b/src/fix_ave_time.cpp index 5fe07faa2f..3537e36d54 100644 --- a/src/fix_ave_time.cpp +++ b/src/fix_ave_time.cpp @@ -243,14 +243,24 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR, "Fix for fix ave/time not computed at compatible time"); - } else if (which[i] == VARIABLE) { + } else if (which[i] == VARIABLE && mode == SCALAR) { int ivariable = input->variable->find(ids[i]); if (ivariable < 0) error->all(FLERR,"Variable name for fix ave/time does not exist"); - if (input->variable->equalstyle(ivariable) == 0) + if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0) error->all(FLERR,"Fix ave/time variable is not equal-style variable"); - if (mode == VECTOR) - error->all(FLERR,"Fix ave/time cannot use variable with vector mode"); + if (argindex[i] && input->variable->vectorstyle(ivariable) == 0) + error->all(FLERR,"Fix ave/time variable is not vector-style variable"); + + } else if (which[i] == VARIABLE && mode == VECTOR) { + int ivariable = input->variable->find(ids[i]); + if (ivariable < 0) + error->all(FLERR,"Variable name for fix ave/time does not exist"); + if (argindex[i] == 0 && input->variable->vectorstyle(ivariable) == 0) + error->all(FLERR,"Fix ave/time variable is not vector-style variable"); + if (argindex[i]) + error->all(FLERR,"Fix ave/time mode vector variable cannot be indexed"); + varlen[i] = 1; } } @@ -282,7 +292,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : if (any_variable_length && (nrepeat > 1 || ave == RUNNING || ave == WINDOW)) { for (int i = 0; i < nvalues; i++) - if (varlen[i]) { + if (varlen[i] && which[i] == COMPUTE) { int icompute = modify->find_compute(ids[i]); modify->compute[icompute]->lock_enable(); } @@ -358,8 +368,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : if (argindex[0] == 0) extscalar = fix->extscalar; else if (fix->extvector >= 0) extscalar = fix->extvector; else extscalar = fix->extlist[argindex[0]-1]; - } else if (which[0] == VARIABLE) + } else if (which[0] == VARIABLE) { extscalar = 0; + } } else { vector_flag = 1; @@ -377,8 +388,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : if (argindex[i] == 0) extlist[i] = fix->extscalar; else if (fix->extvector >= 0) extlist[i] = fix->extvector; else extlist[i] = fix->extlist[argindex[i]-1]; - } else if (which[i] == VARIABLE) + } else if (which[i] == VARIABLE) { extlist[i] = 0; + } } } @@ -405,6 +417,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : for (int i = 0; i < nrows; i++) extlist[i] = fix->extlist[i]; } } else extvector = fix->extarray; + } else if (which[0] == VARIABLE) { + extlist = new int[nrows]; + for (int i = 0; i < nrows; i++) extlist[i] = 0; } } else { @@ -422,6 +437,8 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) : Fix *fix = modify->fix[modify->find_fix(ids[i])]; if (argindex[i] == 0) value = fix->extvector; else value = fix->extarray; + } else if (which[i] == VARIABLE) { + value = 0; } if (value == -1) error->all(FLERR,"Fix ave/time cannot set output array " @@ -518,13 +535,11 @@ void FixAveTime::init() if (icompute < 0) error->all(FLERR,"Compute ID for fix ave/time does not exist"); value2index[i] = icompute; - } else if (which[i] == FIX) { int ifix = modify->find_fix(ids[i]); if (ifix < 0) error->all(FLERR,"Fix ID for fix ave/time does not exist"); value2index[i] = ifix; - } else if (which[i] == VARIABLE) { int ivariable = input->variable->find(ids[i]); if (ivariable < 0) @@ -576,8 +591,7 @@ void FixAveTime::invoke_scalar(bigint ntimestep) double scalar; // zero if first sample within single Nfreq epoch - // NOTE: doc this - // are not checking for returned length, just initialize it + // if any input is variable length, initialize current length // check for exceeding length is done below if (irepeat == 0) { @@ -599,6 +613,7 @@ void FixAveTime::invoke_scalar(bigint ntimestep) m = value2index[i]; // invoke compute if not previously invoked + // insure no out-of-range access to variable-length compute vector if (which[i] == COMPUTE) { Compute *compute = modify->compute[m]; @@ -614,9 +629,6 @@ void FixAveTime::invoke_scalar(bigint ntimestep) compute->compute_vector(); compute->invoked_flag |= INVOKED_VECTOR; } - - // insure no out-of-range access to variable-length compute vector - if (varlen[i] && compute->size_vector < argindex[i]) scalar = 0.0; else scalar = compute->vector[argindex[i]-1]; } @@ -629,10 +641,19 @@ void FixAveTime::invoke_scalar(bigint ntimestep) else scalar = modify->fix[m]->compute_vector(argindex[i]-1); - // evaluate equal-style variable + // evaluate equal-style or vector-style variable + // insure no out-of-range access to vector-style variable - } else if (which[i] == VARIABLE) - scalar = input->variable->compute_equal(m); + } else if (which[i] == VARIABLE) { + if (argindex[i] == 0) + scalar = input->variable->compute_equal(m); + else { + double *varvec; + int nvec = input->variable->compute_vector(m,&varvec); + if (nvec < argindex[i]) scalar = 0.0; + else scalar = varvec[argindex[i]-1]; + } + } // add value to vector or just set directly if offcol is set @@ -750,7 +771,7 @@ void FixAveTime::invoke_vector(bigint ntimestep) bigint ntimestep = update->ntimestep; int lockforever_flag = 0; for (i = 0; i < nvalues; i++) { - if (!varlen[i]) continue; + if (!varlen[i] || which[i] != COMPUTE) continue; if (nrepeat > 1 && ave == ONE) { Compute *compute = modify->compute[value2index[i]]; compute->lock(this,ntimestep,ntimestep+(nrepeat-1)*nevery); @@ -812,6 +833,18 @@ void FixAveTime::invoke_vector(bigint ntimestep) for (i = 0; i < nrows; i++) column[i] = fix->compute_array(i,icol); } + + // evaluate vector-style variable + // insure nvec = nrows, else error + // could be different on this timestep than when column_length(1) set nrows + + } else if (which[j] == VARIABLE) { + double *varvec; + int nvec = input->variable->compute_vector(m,&varvec); + if (nvec != nrows) + error->all(FLERR,"Fix ave/time vector-style variable changed length"); + for (i = 0; i < nrows; i++) + column[i] = varvec[i]; } // add columns of values to array or just set directly if offcol is set @@ -935,7 +968,9 @@ int FixAveTime::column_length(int dynamic) int ifix = modify->find_fix(ids[i]); if (argindex[i] == 0) lengthone = modify->fix[ifix]->size_vector; else lengthone = modify->fix[ifix]->size_array_rows; - } + } else if (which[i] == VARIABLE) { + // variables are always varlen = 1, so dynamic + } if (length == 0) length = lengthone; else if (lengthone != length) error->all(FLERR,"Fix ave/time columns are inconsistent lengths"); @@ -945,15 +980,20 @@ int FixAveTime::column_length(int dynamic) // determine new nrows for dynamic values // either all must be the same // or must match other static values - // don't need to check if not MODE = VECTOR, just inovke lock_length() + // don't need to check if not MODE = VECTOR, just invoke lock_length() if (dynamic) { length = 0; for (int i = 0; i < nvalues; i++) { if (varlen[i] == 0) continue; m = value2index[i]; - Compute *compute = modify->compute[m]; - lengthone = compute->lock_length(); + if (which[i] == COMPUTE) { + Compute *compute = modify->compute[m]; + lengthone = compute->lock_length(); + } else if (which[i] == VARIABLE) { + double *varvec; + lengthone = input->variable->compute_vector(m,&varvec); + } if (mode == SCALAR) continue; if (all_variable_length) { if (length == 0) length = lengthone; diff --git a/src/fix_vector.cpp b/src/fix_vector.cpp index 7abb2f6454..d4bd3087d6 100644 --- a/src/fix_vector.cpp +++ b/src/fix_vector.cpp @@ -111,8 +111,10 @@ FixVector::FixVector(LAMMPS *lmp, int narg, char **arg) : int ivariable = input->variable->find(ids[i]); if (ivariable < 0) error->all(FLERR,"Variable name for fix vector does not exist"); - if (input->variable->equalstyle(ivariable) == 0) + if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0) error->all(FLERR,"Fix vector variable is not equal-style variable"); + if (argindex[i] && input->variable->vectorstyle(ivariable) == 0) + error->all(FLERR,"Fix vector variable is not vector-style variable"); } } @@ -290,10 +292,18 @@ void FixVector::end_of_step() else result[i] = modify->fix[m]->compute_vector(argindex[i]-1); - // evaluate equal-style variable + // evaluate equal-style or vector-style variable } else if (which[i] == VARIABLE) - result[i] = input->variable->compute_equal(m); + if (argindex[i] == 0) + result[i] = input->variable->compute_equal(m); + else { + double *varvec; + int nvec = input->variable->compute_vector(m,&varvec); + int index = argindex[i]; + if (nvec < index) result[i] = 0.0; + else result[i] = varvec[index-1]; + } } // trigger computes on next needed step diff --git a/src/thermo.cpp b/src/thermo.cpp index f95f8c9bbd..be83ef8e94 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -911,11 +911,14 @@ void Thermo::parse_fields(char *str) n = input->variable->find(id); if (n < 0) error->all(FLERR,"Could not find thermo custom variable name"); - if (input->variable->equalstyle(n) == 0) + if (argindex1[nfield] == 0 && input->variable->equalstyle(n) == 0) error->all(FLERR, "Thermo custom variable is not equal-style variable"); - if (argindex1[nfield]) - error->all(FLERR,"Thermo custom variable cannot be indexed"); + if (argindex1[nfield] && input->variable->vectorstyle(n) == 0) + error->all(FLERR, + "Thermo custom variable is not vector-style variable"); + if (argindex2[nfield]) + error->all(FLERR,"Thermo custom variable cannot have two indices"); field2index[nfield] = add_variable(id); addfield(copy,&Thermo::compute_variable,FLOAT); @@ -1474,7 +1477,17 @@ void Thermo::compute_fix() void Thermo::compute_variable() { - dvalue = input->variable->compute_equal(variables[field2index[ifield]]); + int iarg = argindex1[ifield]; + + if (iarg == 0) + dvalue = input->variable->compute_equal(variables[field2index[ifield]]); + else { + double *varvec; + int nvec = + input->variable->compute_vector(variables[field2index[ifield]],&varvec); + if (nvec < iarg) dvalue = 0.0; + else dvalue = varvec[iarg-1]; + } } /* ---------------------------------------------------------------------- diff --git a/src/variable.cpp b/src/variable.cpp index d024223beb..cd41735eda 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -52,7 +52,7 @@ using namespace MathConst; #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a) enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,GETENV, - SCALARFILE,ATOMFILE,FORMAT,EQUAL,ATOM,PYTHON}; + SCALARFILE,ATOMFILE,FORMAT,EQUAL,ATOM,VECTOR,PYTHON}; enum{ARG,OP}; // customize by adding a function @@ -65,7 +65,7 @@ enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY, RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,LOGFREQ2, STRIDE,STRIDE2,VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,GRMASK, IS_ACTIVE,IS_DEFINED,IS_AVAILABLE, - VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY}; + VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY,VECTORARRAY}; // customize by adding a special function @@ -92,6 +92,7 @@ Variable::Variable(LAMMPS *lmp) : Pointers(lmp) pad = NULL; reader = NULL; data = NULL; + vecs = NULL; eval_in_progress = NULL; @@ -125,6 +126,7 @@ Variable::~Variable() if (style[i] == LOOP || style[i] == ULOOP) delete [] data[i][0]; else for (int j = 0; j < num[i]; j++) delete [] data[i][j]; delete [] data[i]; + if (style[i] == VECTOR) memory->destroy(vecs[i].values); } memory->sfree(names); memory->destroy(style); @@ -133,6 +135,7 @@ Variable::~Variable() memory->destroy(pad); memory->sfree(reader); memory->sfree(data); + memory->sfree(vecs); memory->destroy(eval_in_progress); @@ -428,6 +431,30 @@ void Variable::set(int narg, char **arg) copy(1,&arg[2],data[nvar]); } + // VECTOR + // replace pre-existing var if also style VECTOR (allows it to be reset) + // num = 1, which = 1st value + // data = 1 value, string to eval + + } else if (strcmp(arg[1],"vector") == 0) { + if (narg != 3) error->all(FLERR,"Illegal variable command"); + int ivar = find(arg[0]); + if (ivar >= 0) { + if (style[ivar] != VECTOR) + error->all(FLERR,"Cannot redefine variable as a different style"); + delete [] data[ivar][0]; + copy(1,&arg[2],data[ivar]); + replaceflag = 1; + } else { + if (nvar == maxvar) grow(); + style[nvar] = VECTOR; + num[nvar] = 1; + which[nvar] = 0; + pad[nvar] = 0; + data[nvar] = new char*[num[nvar]]; + copy(1,&arg[2],data[nvar]); + } + // PYTHON // replace pre-existing var if also style PYTHON (allows it to be reset) // num = 2, which = 1st value @@ -530,12 +557,12 @@ int Variable::next(int narg, char **arg) error->all(FLERR,"All variables in next command must be same style"); } - // invalid styles: STRING, EQUAL, WORLD, ATOM, GETENV, FORMAT, PYTHON + // invalid styles: STRING, EQUAL, WORLD, ATOM, VECTOR, GETENV, FORMAT, PYTHON int istyle = style[find(arg[0])]; if (istyle == STRING || istyle == EQUAL || istyle == WORLD || - istyle == GETENV || istyle == ATOM || istyle == FORMAT || - istyle == PYTHON) + istyle == GETENV || istyle == ATOM || istyle == VECTOR || + istyle == FORMAT || istyle == PYTHON) error->all(FLERR,"Invalid variable style with next command"); // if istyle = UNIVERSE or ULOOP, insure all such variables are incremented @@ -712,6 +739,17 @@ int Variable::atomstyle(int ivar) return 0; } +/* ---------------------------------------------------------------------- + return 1 if variable is VECTOR style, 0 if not + this is checked before call to compute_vector() to return a vector of doubles +------------------------------------------------------------------------- */ + +int Variable::vectorstyle(int ivar) +{ + if (style[ivar] == VECTOR) return 1; + return 0; +} + /* ---------------------------------------------------------------------- check if variable with name is PYTHON and matches funcname called by Python class before it invokes a Python function @@ -737,7 +775,7 @@ char *Variable::pythonstyle(char *name, char *funcname) if FORMAT, evaluate its variable and put formatted result in str if GETENV, query environment and put result in str if PYTHON, evaluate Python function, it will put result in str - if ATOM or ATOMFILE, return NULL + if ATOM or ATOMFILE or VECTOR, return NULL return NULL if no variable with name, or which value is bad, caller must respond ------------------------------------------------------------------------- */ @@ -797,8 +835,9 @@ char *Variable::retrieve(char *name) error->all(FLERR,"Python variable does not match Python function"); python->invoke_function(ifunc,data[ivar][1]); str = data[ivar][1]; - } else if (style[ivar] == ATOM || style[ivar] == ATOMFILE) return NULL; - + } else if (style[ivar] == ATOM || style[ivar] == ATOMFILE || + style[ivar] == VECTOR) return NULL; + eval_in_progress[ivar] = 0; return str; @@ -859,6 +898,7 @@ void Variable::compute_atom(int ivar, int igroup, eval_in_progress[ivar] = 1; if (style[ivar] == ATOM) { + treetype = ATOM; evaluate(data[ivar][0],&tree); collapse_tree(tree); } else vstore = reader[ivar]->fixstore->vstore; @@ -908,10 +948,58 @@ void Variable::compute_atom(int ivar, int igroup, } if (style[ivar] == ATOM) free_tree(tree); - eval_in_progress[ivar] = 0; } +/* ---------------------------------------------------------------------- + compute result of vector-style variable evaluation + return length of vector and result pointer to vector values + if length == 0 or -1 (mismatch), generate an error + if variable already computed on this timestep, just return + else evaulate the formula and its length, store results in VecVar entry +------------------------------------------------------------------------- */ + +int Variable::compute_vector(int ivar, double **result) +{ + Tree *tree; + if (vecs[ivar].currentstep == update->ntimestep) { + *result = vecs[ivar].values; + return vecs[ivar].n; + } + + if (eval_in_progress[ivar]) + error->all(FLERR,"Variable has circular dependency"); + eval_in_progress[ivar] = 1; + + treetype = VECTOR; + evaluate(data[ivar][0],&tree); + collapse_tree(tree); + int nlen = size_tree_vector(tree); + if (nlen == 0) error->all(FLERR,"Vector-style variable has zero length"); + if (nlen < 0) error->all(FLERR, + "Inconsistent lengths in vector-style variable"); + + // (re)allocate space for results if necessary + + if (nlen > vecs[ivar].nmax) { + memory->destroy(vecs[ivar].values); + vecs[ivar].nmax = nlen; + memory->create(vecs[ivar].values,vecs[ivar].nmax,"variable:values"); + } + + vecs[ivar].n = nlen; + vecs[ivar].currentstep = update->ntimestep; + double *vec = vecs[ivar].values; + for (int i = 0; i < nlen; i++) + vec[i] = eval_tree(tree,i); + + free_tree(tree); + eval_in_progress[ivar] = 0; + + *result = vec; + return nlen; +} + /* ---------------------------------------------------------------------- save copy of EQUAL style ivar formula in copy allocate copy here, later equal_restore() call will free it @@ -997,6 +1085,13 @@ void Variable::grow() data = (char ***) memory->srealloc(data,maxvar*sizeof(char **),"var:data"); + vecs = (VecVar *) memory->srealloc(vecs,maxvar*sizeof(VecVar),"var:vecvar"); + for (int i = old; i < maxvar; i++) { + vecs[i].nmax = 0; + vecs[i].currentstep = -1; + vecs[i].values = NULL; + } + memory->grow(eval_in_progress,maxvar,"var:eval_in_progress"); for (int i = 0; i < maxvar; i++) eval_in_progress[i] = 0; } @@ -1017,7 +1112,8 @@ void Variable::copy(int narg, char **from, char **to) /* ---------------------------------------------------------------------- recursive evaluation of a string str - str is an equal-style or atom-style formula containing one or more items: + str is an equal-style or atom-style or vector-style formula + containing one or more items: number = 0.0, -5.45, 2.8e-4, ... constant = PI, version, yes, no, on, off thermo keyword = ke, vol, atoms, ... @@ -1034,7 +1130,7 @@ void Variable::copy(int narg, char **from, char **to) variable = v_name, v_name[i] equal-style variables passes in tree = NULL: evaluate the formula, return result as a double - atom-style variable passes in tree = non-NULL: + atom-style and vector-style variables pass in tree = non-NULL: parse the formula but do not evaluate it create a parse tree and return it ------------------------------------------------------------------------- */ @@ -1145,11 +1241,17 @@ double Variable::evaluate(char *str, Tree **tree) // compute // ---------------- - if (strncmp(word,"c_",2) == 0) { + if (strncmp(word,"c_",2) == 0 || strncmp(word,"C_",2) == 0) { if (domain->box_exist == 0) error->all(FLERR, "Variable evaluation before simulation box is defined"); + // uppercase used to force access of + // global vector vs global scalar, and global array vs global vector + + int lowercase = 1; + if (word[0] == 'C') lowercase = 0; + int icompute = modify->find_compute(word+2); if (icompute < 0) error->all(FLERR,"Invalid compute ID in variable formula"); @@ -1176,9 +1278,9 @@ double Variable::evaluate(char *str, Tree **tree) } } - // c_ID = scalar from global scalar + // c_ID = scalar from global scalar, must be lowercase - if (nbracket == 0 && compute->scalar_flag) { + if (nbracket == 0 && compute->scalar_flag && lowercase) { if (update->whichflag == 0) { if (compute->invoked_scalar != update->ntimestep) @@ -1199,9 +1301,9 @@ double Variable::evaluate(char *str, Tree **tree) treestack[ntreestack++] = newtree; } else argstack[nargstack++] = value1; - // c_ID[i] = scalar from global vector + // c_ID[i] = scalar from global vector, must be lowercase - } else if (nbracket == 1 && compute->vector_flag) { + } else if (nbracket == 1 && compute->vector_flag && lowercase) { if (index1 > compute->size_vector && compute->size_vector_variable == 0) @@ -1228,9 +1330,9 @@ double Variable::evaluate(char *str, Tree **tree) treestack[ntreestack++] = newtree; } else argstack[nargstack++] = value1; - // c_ID[i][j] = scalar from global array + // c_ID[i][j] = scalar from global array, must be lowercase - } else if (nbracket == 2 && compute->array_flag) { + } else if (nbracket == 2 && compute->array_flag && lowercase) { if (index1 > compute->size_array_rows && compute->size_array_rows_variable == 0) @@ -1260,6 +1362,68 @@ double Variable::evaluate(char *str, Tree **tree) treestack[ntreestack++] = newtree; } else argstack[nargstack++] = value1; + // c_ID = vector from global vector, lowercase or uppercase + + } else if (nbracket == 0 && compute->vector_flag) { + + if (tree == NULL) + error->all(FLERR, + "Compute global vector in equal-style variable formula"); + if (treetype == ATOM) + error->all(FLERR, + "Compute global vector in atom-style variable formula"); + if (compute->size_vector == 0) + error->all(FLERR,"Variable formula compute vector is zero length"); + if (update->whichflag == 0) { + if (compute->invoked_vector != update->ntimestep) + error->all(FLERR,"Compute used in variable between runs " + "is not current"); + } else if (!(compute->invoked_flag & INVOKED_VECTOR)) { + compute->compute_vector(); + compute->invoked_flag |= INVOKED_VECTOR; + } + + Tree *newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = compute->vector; + newtree->nvector = compute->size_vector; + newtree->nstride = 1; + newtree->selfalloc = 0; + newtree->first = newtree->second = NULL; + newtree->nextra = 0; + treestack[ntreestack++] = newtree; + + // c_ID[i] = vector from global array, lowercase or uppercase + + } else if (nbracket == 1 && compute->array_flag) { + + if (tree == NULL) + error->all(FLERR, + "Compute global vector in equal-style variable formula"); + if (treetype == ATOM) + error->all(FLERR, + "Compute global vector in atom-style variable formula"); + if (compute->size_array_rows == 0) + error->all(FLERR,"Variable formula compute array is zero length"); + if (update->whichflag == 0) { + if (compute->invoked_array != update->ntimestep) + error->all(FLERR,"Compute used in variable between runs " + "is not current"); + } else if (!(compute->invoked_flag & INVOKED_ARRAY)) { + compute->compute_array(); + compute->invoked_flag |= INVOKED_ARRAY; + } + + Tree *newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = &compute->array[0][index1-1]; + newtree->nvector = compute->size_array_rows; + newtree->nstride = compute->size_array_cols; + newtree->selfalloc = 0; + newtree->first = newtree->second = NULL; + newtree->nextra = 0; + treestack[ntreestack++] = newtree; + // c_ID[i] = scalar from per-atom vector } else if (nbracket == 1 && compute->peratom_flag && @@ -1311,6 +1475,9 @@ double Variable::evaluate(char *str, Tree **tree) if (tree == NULL) error->all(FLERR, "Per-atom compute in equal-style variable formula"); + if (treetype == VECTOR) + error->all(FLERR, + "Per-atom compute in vector-style variable formula"); if (update->whichflag == 0) { if (compute->invoked_peratom != update->ntimestep) error->all(FLERR,"Compute used in variable between runs " @@ -1337,6 +1504,9 @@ double Variable::evaluate(char *str, Tree **tree) if (tree == NULL) error->all(FLERR, "Per-atom compute in equal-style variable formula"); + if (treetype == VECTOR) + error->all(FLERR, + "Per-atom compute in vector-style variable formula"); if (index1 > compute->size_peratom_cols) error->all(FLERR,"Variable formula compute array " "is accessed out-of-range"); @@ -1367,11 +1537,17 @@ double Variable::evaluate(char *str, Tree **tree) // fix // ---------------- - } else if (strncmp(word,"f_",2) == 0) { + } else if (strncmp(word,"f_",2) == 0 || strncmp(word,"F_",2) == 0) { if (domain->box_exist == 0) error->all(FLERR, "Variable evaluation before simulation box is defined"); + // uppercase used to force access of + // global vector vs global scalar, and global array vs global vector + + int lowercase = 1; + if (word[0] == 'F') lowercase = 0; + int ifix = modify->find_fix(word+2); if (ifix < 0) error->all(FLERR,"Invalid fix ID in variable formula"); Fix *fix = modify->fix[ifix]; @@ -1397,9 +1573,9 @@ double Variable::evaluate(char *str, Tree **tree) } } - // f_ID = scalar from global scalar + // f_ID = scalar from global scalar, must be lowercase - if (nbracket == 0 && fix->scalar_flag) { + if (nbracket == 0 && fix->scalar_flag && lowercase) { if (update->whichflag > 0 && update->ntimestep % fix->global_freq) error->all(FLERR,"Fix in variable not computed at compatible time"); @@ -1414,9 +1590,9 @@ double Variable::evaluate(char *str, Tree **tree) treestack[ntreestack++] = newtree; } else argstack[nargstack++] = value1; - // f_ID[i] = scalar from global vector + // f_ID[i] = scalar from global vector, must be lowercase - } else if (nbracket == 1 && fix->vector_flag) { + } else if (nbracket == 1 && fix->vector_flag && lowercase) { if (index1 > fix->size_vector && fix->size_vector_variable == 0) @@ -1435,9 +1611,9 @@ double Variable::evaluate(char *str, Tree **tree) treestack[ntreestack++] = newtree; } else argstack[nargstack++] = value1; - // f_ID[i][j] = scalar from global array + // f_ID[i][j] = scalar from global array, must be lowercase - } else if (nbracket == 2 && fix->array_flag) { + } else if (nbracket == 2 && fix->array_flag && lowercase) { if (index1 > fix->size_array_rows && fix->size_array_rows_variable == 0) @@ -1459,6 +1635,68 @@ double Variable::evaluate(char *str, Tree **tree) treestack[ntreestack++] = newtree; } else argstack[nargstack++] = value1; + // f_ID = vector from global vector, lowercase or uppercase + + } else if (nbracket == 0 && fix->vector_flag) { + + if (update->whichflag > 0 && update->ntimestep % fix->global_freq) + error->all(FLERR,"Fix in variable not computed at compatible time"); + if (tree == NULL) + error->all(FLERR,"Fix global vector in " + "equal-style variable formula"); + if (treetype == ATOM) + error->all(FLERR,"Fix global vector in " + "atom-style variable formula"); + if (fix->size_vector == 0) + error->all(FLERR,"Variable formula fix vector is zero length"); + + int nvec = fix->size_vector; + double *vec; + memory->create(vec,nvec,"variable:values"); + for (int m = 0; m < nvec; m++) + vec[m] = fix->compute_vector(m); + + Tree *newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = vec; + newtree->nvector = nvec; + newtree->nstride = 1; + newtree->selfalloc = 1; + newtree->first = newtree->second = NULL; + newtree->nextra = 0; + treestack[ntreestack++] = newtree; + + // f_ID[i] = vector from global array, lowercase or uppercase + + } else if (nbracket == 1 && fix->array_flag) { + + if (update->whichflag > 0 && update->ntimestep % fix->global_freq) + error->all(FLERR,"Fix in variable not computed at compatible time"); + if (tree == NULL) + error->all(FLERR,"Fix global vector in " + "equal-style variable formula"); + if (treetype == ATOM) + error->all(FLERR,"Fix global vector in " + "atom-style variable formula"); + if (fix->size_array_rows == 0) + error->all(FLERR,"Variable formula fix array is zero length"); + + int nvec = fix->size_array_rows; + double *vec; + memory->create(vec,nvec,"variable:values"); + for (int m = 0; m < nvec; m++) + vec[m] = fix->compute_array(m,index1-1); + + Tree *newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = vec; + newtree->nvector = nvec; + newtree->nstride = 1; + newtree->selfalloc = 1; + newtree->first = newtree->second = NULL; + newtree->nextra = 0; + treestack[ntreestack++] = newtree; + // f_ID[i] = scalar from per-atom vector } else if (nbracket == 1 && fix->peratom_flag && @@ -1568,9 +1806,10 @@ double Variable::evaluate(char *str, Tree **tree) i = ptr-str+1; } - // v_name = scalar from non atom/atomfile variable + // v_name = scalar from non atom/atomfile and non vector-style variable - if (nbracket == 0 && style[ivar] != ATOM && style[ivar] != ATOMFILE) { + if (nbracket == 0 && style[ivar] != ATOM && style[ivar] != ATOMFILE && + style[ivar] != VECTOR) { char *var = retrieve(word+2); if (var == NULL) @@ -1592,6 +1831,10 @@ double Variable::evaluate(char *str, Tree **tree) if (tree == NULL) error->all(FLERR, "Atom-style variable in equal-style variable formula"); + if (treetype == VECTOR) + error->all(FLERR, + "Atom-style variable in vector-style variable formula"); + Tree *newtree; evaluate(data[ivar][0],&newtree); treestack[ntreestack++] = newtree; @@ -1603,6 +1846,10 @@ double Variable::evaluate(char *str, Tree **tree) if (tree == NULL) error->all(FLERR,"Atomfile-style variable in " "equal-style variable formula"); + if (treetype == VECTOR) + error->all(FLERR, + "Atomfile-style variable in vector-style variable formula"); + Tree *newtree = new Tree(); newtree->type = ATOMARRAY; newtree->array = reader[ivar]->fixstore->vstore; @@ -1612,6 +1859,31 @@ double Variable::evaluate(char *str, Tree **tree) newtree->nextra = 0; treestack[ntreestack++] = newtree; + // v_name = vector from vector-style variable + // evaluate the vector-style variable, put result in newtree + + } else if (nbracket == 0 && style[ivar] == VECTOR) { + + if (tree == NULL) + error->all(FLERR, + "Vector-style variable in equal-style variable formula"); + if (treetype == ATOM) + error->all(FLERR, + "Vector-style variable in atom-style variable formula"); + + double *vec; + int nvec = compute_vector(ivar,&vec); + + Tree *newtree = new Tree(); + newtree->type = VECTORARRAY; + newtree->array = vec; + newtree->nvector = nvec; + newtree->nstride = 1; + newtree->selfalloc = 0; + newtree->first = newtree->second = NULL; + newtree->nextra = 0; + 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 @@ -1632,6 +1904,26 @@ double Variable::evaluate(char *str, Tree **tree) peratom2global(1,NULL,reader[ivar]->fixstore->vstore,1,index, tree,treestack,ntreestack,argstack,nargstack); + // v_name[N] = scalar from vector-style variable + // compute the vector-style variable, extract single value + + } else if (nbracket && style[ivar] == VECTOR) { + + double *vec; + int nvec = compute_vector(ivar,&vec); + if (index <= 0 || index > nvec) + error->all(FLERR,"Invalid index into vector-style variable"); + int m = index; // convert from tagint to int + + if (tree) { + Tree *newtree = new Tree(); + newtree->type = VALUE; + newtree->value = vec[m-1]; + newtree->first = newtree->second = NULL; + newtree->nextra = 0; + treestack[ntreestack++] = newtree; + } else argstack[nargstack++] = vec[m-1]; + } else error->all(FLERR,"Mismatched variable in variable formula"); // ---------------- @@ -1892,7 +2184,7 @@ double Variable::evaluate(char *str, Tree **tree) one-time collapse of an atom-style variable parse tree tree was created by one-time parsing of formula string via evaluate() only keep tree nodes that depend on - ATOMARRAY, TYPEARRAY, INTARRAY, BIGINTARRAY + ATOMARRAY, TYPEARRAY, INTARRAY, BIGINTARRAY, VECTOR remainder is converted to single VALUE this enables optimal eval_tree loop over atoms customize by adding a function: @@ -1912,6 +2204,7 @@ double Variable::collapse_tree(Tree *tree) if (tree->type == TYPEARRAY) return 0.0; if (tree->type == INTARRAY) return 0.0; if (tree->type == BIGINTARRAY) return 0.0; + if (tree->type == VECTORARRAY) return 0.0; if (tree->type == ADD) { arg1 = collapse_tree(tree->first); @@ -2403,7 +2696,8 @@ double Variable::collapse_tree(Tree *tree) } /* ---------------------------------------------------------------------- - evaluate an atom-style variable parse tree for atom I + evaluate an atom-style or vector-style variable parse tree + index I = atom I or vector index I tree was created by one-time parsing of formula string via evaulate() customize by adding a function: sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan(), @@ -2422,6 +2716,7 @@ double Variable::eval_tree(Tree *tree, int i) if (tree->type == TYPEARRAY) return tree->array[atom->type[i]]; if (tree->type == INTARRAY) return (double) tree->iarray[i*tree->nstride]; if (tree->type == BIGINTARRAY) return (double) tree->barray[i*tree->nstride]; + if (tree->type == VECTORARRAY) return tree->array[i*tree->nstride]; if (tree->type == ADD) return eval_tree(tree->first,i) + eval_tree(tree->second,i); @@ -2728,6 +3023,40 @@ double Variable::eval_tree(Tree *tree, int i) return 0.0; } +/* ---------------------------------------------------------------------- + scan entire tree, find size of vectors for vector-style variable + return N for consistent vector size + return 0 for no vector size, caller flags as error + return -1 for inconsistent vector size, caller flags as error +------------------------------------------------------------------------- */ + +int Variable::size_tree_vector(Tree *tree) +{ + int nsize = 0; + if (tree->type == VECTORARRAY) nsize = tree->nvector; + if (tree->first) nsize = compare_tree_vector(nsize,size_tree_vector(tree->first)); + if (tree->second) nsize = compare_tree_vector(nsize,size_tree_vector(tree->second)); + if (tree->nextra) { + for (int i = 0; i < tree->nextra; i++) + nsize = compare_tree_vector(nsize,size_tree_vector(tree->extra[i])); + } + return nsize; +} + +/* ---------------------------------------------------------------------- + compare size of two vectors for vector-style variable + return positive size if same or one has no size 0 + return -1 error if one is already error or not same positive size +------------------------------------------------------------------------- */ + +int Variable::compare_tree_vector(int i, int j) +{ + if (i < 0 || j < 0) return -1; + if (i == 0 || j == 0) return MAX(i,j); + if (i != j) return -1; + return i; +} + /* ---------------------------------------------------------------------- */ void Variable::free_tree(Tree *tree) @@ -2739,9 +3068,7 @@ void Variable::free_tree(Tree *tree) delete [] tree->extra; } - if (tree->type == ATOMARRAY && tree->selfalloc) - memory->destroy(tree->array); - + if (tree->selfalloc) memory->destroy(tree->array); delete tree; } @@ -3532,9 +3859,12 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Compute *compute = NULL; Fix *fix = NULL; + int ivar = -1; int index,nvec,nstride; char *ptr1,*ptr2; + // argument is compute + if (strstr(args[0],"c_") == args[0]) { ptr1 = strchr(args[0],'['); if (ptr1) { @@ -3574,6 +3904,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree, nstride = compute->size_array_cols; } else error->all(FLERR,"Mismatched compute in variable formula"); + // argument is fix + } else if (strstr(args[0],"f_") == args[0]) { ptr1 = strchr(args[0],'['); if (ptr1) { @@ -3600,6 +3932,28 @@ int Variable::special_function(char *word, char *contents, Tree **tree, nstride = fix->size_array_cols; } else error->all(FLERR,"Mismatched fix in variable formula"); + // argument is vector-style variable + + } else if (strstr(args[0],"v_") == args[0]) { + ptr1 = strchr(args[0],'['); + if (ptr1) { + ptr2 = ptr1; + index = (int) int_between_brackets(ptr2,0); + *ptr1 = '\0'; + } else index = 0; + + if (index) + error->all(FLERR,"Invalid special function in variable formula"); + ivar = find(&args[0][2]); + if (ivar < 0) error->all(FLERR,"Invalid special function in variable formula"); + if (style[ivar] != VECTOR) + error->all(FLERR,"Mis-matched special function variable in variable formula"); + if (eval_in_progress[ivar]) error->all(FLERR,"Variable has circular dependency"); + + double *vec; + nvec = compute_vector(ivar,&vec); + nstride = 1; + } else error->all(FLERR,"Invalid special function in variable formula"); value = 0.0; @@ -3661,6 +4015,28 @@ int Variable::special_function(char *word, char *contents, Tree **tree, } } + if (ivar >= 0) { + double one; + double *vec = vecs[ivar].values; + for (int i = 0; i < nvec; i++) { + one = vec[i]; + if (method == SUM) value += one; + else if (method == XMIN) value = MIN(value,one); + else if (method == XMAX) value = MAX(value,one); + else if (method == AVE) value += one; + else if (method == TRAP) value += one; + else if (method == SLOPE) { + if (nvec > 1) xvalue = (double) i / (nvec-1); + else xvalue = 0.0; + sx += xvalue; + sy += one; + sxx += xvalue*xvalue; + sxy += xvalue*one; + } + } + if (method == TRAP) value -= 0.5*vec[0] + 0.5*vec[nvec-1]; + } + if (method == AVE) value /= nvec; if (method == SLOPE) { diff --git a/src/variable.h b/src/variable.h index 05ea3f57e9..eecf905015 100644 --- a/src/variable.h +++ b/src/variable.h @@ -35,12 +35,14 @@ class Variable : protected Pointers { int equalstyle(int); int atomstyle(int); + int vectorstyle(int); char *pythonstyle(char *, char *); char *retrieve(char *); double compute_equal(int); double compute_equal(char *); void compute_atom(int, int, double *, int, int); + int compute_vector(int, double **); tagint int_between_brackets(char *&, int); double evaluate_boolean(char *); @@ -63,7 +65,15 @@ class Variable : protected Pointers { class VarReader **reader; // variable that reads from file char ***data; // str value of each variable's values - int *eval_in_progress; // flag if evaluation of variable is in progress + struct VecVar { + int n,nmax; + bigint currentstep; + double *values; + }; + VecVar *vecs; + + int *eval_in_progress; // flag if evaluation of variable is in progress + int treetype; // ATOM or VECTOR flag for formula evaluation class RanMars *randomequal; // random number generator for equal-style vars class RanMars *randomatom; // random number generator for atom-style vars @@ -73,12 +83,13 @@ class Variable : protected Pointers { class Python *python; // ptr to embedded Python interpreter - struct Tree { // parse tree for atom-style variables + struct Tree { // parse tree for atom-style or vector-style variables double value; // single scalar double *array; // per-atom or per-type list of doubles int *iarray; // per-atom list of ints bigint *barray; // per-atom list of bigints int type; // operation, see enum{} in variable.cpp + int nvector; // length of array for vector-style variable 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 @@ -94,6 +105,8 @@ class Variable : protected Pointers { double evaluate(char *, Tree **); double collapse_tree(Tree *); double eval_tree(Tree *, int); + int size_tree_vector(Tree *); + int compare_tree_vector(int, int); void free_tree(Tree *); int find_matching_paren(char *, int, char *&); int math_function(char *, char *, Tree **, Tree **, int &, double *, int &);