From b3557b81c14cddac04524f96726332dda3ec7b52 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 14 Oct 2022 17:49:56 -0400 Subject: [PATCH] convert fix ave/histo and fix ave/histo/weight --- src/fix_ave_histo.cpp | 644 +++++++++++++++-------------------- src/fix_ave_histo.h | 18 +- src/fix_ave_histo_weight.cpp | 326 +++++++++--------- 3 files changed, 442 insertions(+), 546 deletions(-) diff --git a/src/fix_ave_histo.cpp b/src/fix_ave_histo.cpp index c2acd652ad..46c26acad7 100644 --- a/src/fix_ave_histo.cpp +++ b/src/fix_ave_histo.cpp @@ -16,6 +16,7 @@ #include "arg_info.h" #include "atom.h" +#include "comm.h" #include "compute.h" #include "error.h" #include "input.h" @@ -29,29 +30,24 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{ONE,RUNNING}; -enum{SCALAR,VECTOR,WINDOW}; -enum{DEFAULT,GLOBAL,PERATOM,LOCAL}; -enum{IGNORE,END,EXTRA}; - +enum { ONE, RUNNING }; +enum { SCALAR, VECTOR, WINDOW }; +enum { DEFAULT, GLOBAL, PERATOM, LOCAL }; +enum { IGNORE, END, EXTRA }; #define BIG 1.0e20 /* ---------------------------------------------------------------------- */ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - nvalues(0), which(nullptr), argindex(nullptr), value2index(nullptr), - ids(nullptr), fp(nullptr), stats_list(nullptr), - bin(nullptr), bin_total(nullptr), bin_all(nullptr), bin_list(nullptr), - coord(nullptr), vector(nullptr) + Fix(lmp, narg, arg), nvalues(0), fp(nullptr), stats_list(nullptr), bin(nullptr), + bin_total(nullptr), bin_all(nullptr), bin_list(nullptr), coord(nullptr), vector(nullptr) { - if (narg < 10) error->all(FLERR,"Illegal fix ave/histo command"); + auto mycmd = fmt::format("fix {}", style); + if (narg < 10) utils::missing_cmd_args(FLERR, mycmd, error); - MPI_Comm_rank(world,&me); - - nevery = utils::inumeric(FLERR,arg[3],false,lmp); - nrepeat = utils::inumeric(FLERR,arg[4],false,lmp); - nfreq = utils::inumeric(FLERR,arg[5],false,lmp); + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + nrepeat = utils::inumeric(FLERR, arg[4], false, lmp); + nfreq = utils::inumeric(FLERR, arg[5], false, lmp); global_freq = nfreq; vector_flag = 1; @@ -63,9 +59,9 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : dynamic_group_allow = 1; time_depend = 1; - lo = utils::numeric(FLERR,arg[6],false,lmp); - hi = utils::numeric(FLERR,arg[7],false,lmp); - nbins = utils::inumeric(FLERR,arg[8],false,lmp); + lo = utils::numeric(FLERR, arg[6], false, lmp); + hi = utils::numeric(FLERR, arg[7], false, lmp); + nbins = utils::inumeric(FLERR, arg[8], false, lmp); // scan values to count them // then read options so know mode = SCALAR/VECTOR before re-reading values @@ -91,7 +87,7 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : } else break; } - if (nvalues == 0) error->all(FLERR,"No values in fix ave/histo command"); + if (nvalues == 0) error->all(FLERR,"No values in {} command", mycmd); options(iarg,narg,arg); @@ -100,70 +96,65 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : int expand = 0; char **earg; - nvalues = utils::expand_args(FLERR,nvalues,&arg[9],mode,earg,lmp); + nvalues = utils::expand_args(FLERR, nvalues, &arg[9], mode, earg, lmp); if (earg != &arg[9]) expand = 1; arg = earg; // parse values - which = new int[nvalues]; - argindex = new int[nvalues]; - value2index = new int[nvalues]; - ids = new char*[nvalues]; + values.clear(); for (int i = 0; i < nvalues; i++) { + value_t val; + val.id = ""; + val.val.c = nullptr; + if (strcmp(arg[i],"x") == 0) { - which[i] = ArgInfo::X; - argindex[i] = 0; - ids[i] = nullptr; + val.which = ArgInfo::X; + val.argindex = 0; } else if (strcmp(arg[i],"y") == 0) { - which[i] = ArgInfo::X; - argindex[i] = 1; - ids[i] = nullptr; + val.which = ArgInfo::X; + val.argindex = 1; } else if (strcmp(arg[i],"z") == 0) { - which[i] = ArgInfo::X; - argindex[i] = 2; - ids[i] = nullptr; + val.which = ArgInfo::X; + val.argindex = 2; } else if (strcmp(arg[i],"vx") == 0) { - which[i] = ArgInfo::V; - argindex[i] = 0; - ids[i] = nullptr; + val.which = ArgInfo::V; + val.argindex = 0; } else if (strcmp(arg[i],"vy") == 0) { - which[i] = ArgInfo::V; - argindex[i] = 1; - ids[i] = nullptr; + val.which = ArgInfo::V; + val.argindex = 1; } else if (strcmp(arg[i],"vz") == 0) { - which[i] = ArgInfo::V; - argindex[i] = 2; - ids[i] = nullptr; + val.which = ArgInfo::V; + val.argindex = 2; } else if (strcmp(arg[i],"fx") == 0) { - which[i] = ArgInfo::F; - argindex[i] = 0; - ids[i] = nullptr; + val.which = ArgInfo::F; + val.argindex = 0; } else if (strcmp(arg[i],"fy") == 0) { - which[i] = ArgInfo::F; - argindex[i] = 1; - ids[i] = nullptr; + val.which = ArgInfo::F; + val.argindex = 1; } else if (strcmp(arg[i],"fz") == 0) { - which[i] = ArgInfo::F; - argindex[i] = 2; - ids[i] = nullptr; + val.which = ArgInfo::F; + val.argindex = 2; } else { ArgInfo argi(arg[i]); if (argi.get_type() == ArgInfo::NONE) break; if ((argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_dim() > 1)) - error->all(FLERR,"Invalid fix ave/histo command"); + error->all(FLERR,"Invalid {} argument: {}", mycmd, arg[i]); - which[i] = argi.get_type(); - argindex[i] = argi.get_index1(); - ids[i] = argi.copy_name(); + val.which = argi.get_type(); + val.argindex = argi.get_index1(); + val.id = argi.get_name(); } + values.push_back(val); } + if (nvalues != (int)values.size()) + error->all(FLERR, "Could not parse value data consistently for {}", mycmd); // if wildcard expansion occurred, free earg memory from expand_args() @@ -175,69 +166,67 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : // check input args for kind consistency // all inputs must all be global, per-atom, or local - if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0) - error->all(FLERR,"Illegal fix ave/histo command"); + if (nevery <= 0) + error->all(FLERR,"Illegal {} nevery value: {}", mycmd, nevery); + if (nrepeat <= 0) + error->all(FLERR,"Illegal {} nrepeat value: {}", mycmd, nrepeat); + if (nfreq <= 0) + error->all(FLERR,"Illegal {} nfreq value: {}", mycmd, nfreq); if (nfreq % nevery || nrepeat*nevery > nfreq) - error->all(FLERR,"Illegal fix ave/histo command"); - if (lo >= hi) error->all(FLERR,"Illegal fix ave/histo command"); - if (nbins <= 0) error->all(FLERR,"Illegal fix ave/histo command"); + error->all(FLERR,"Inconsistent {} nevery/nrepeat/nfreq values", mycmd); if (ave != RUNNING && overwrite) - error->all(FLERR,"Illegal fix ave/histo command"); + error->all(FLERR,"{} overwrite keyword requires ave running setting", mycmd); int kindglobal,kindperatom,kindlocal; - - for (int i = 0; i < nvalues; i++) { + for (auto &val : values) { kindglobal = kindperatom = kindlocal = 0; - if ((which[i] == ArgInfo::X) || (which[i] == ArgInfo::V) - || (which[i] == ArgInfo::F)) { + if ((val.which == ArgInfo::X) || (val.which == ArgInfo::V) || (val.which == ArgInfo::F)) { kindperatom = 1; - } else if (which[i] == ArgInfo::COMPUTE) { - int c_id = modify->find_compute(ids[i]); - if (c_id < 0) error->all(FLERR,"Fix ave/histo input is invalid compute"); - Compute *compute = modify->compute[c_id]; + } else if (val.which == ArgInfo::COMPUTE) { + val.val.c = modify->get_compute_by_id(val.id); + if (!val.val.c) error->all(FLERR,"Compute ID {} for {} does not exist", val.id, mycmd); // computes can produce multiple kinds of output - if (compute->scalar_flag || compute->vector_flag || compute->array_flag) + if (val.val.c->scalar_flag || val.val.c->vector_flag || val.val.c->array_flag) kindglobal = 1; - if (compute->peratom_flag) kindperatom = 1; - if (compute->local_flag) kindlocal = 1; + if (val.val.c->peratom_flag) kindperatom = 1; + if (val.val.c->local_flag) kindlocal = 1; - } else if (which[i] == ArgInfo::FIX) { - int f_id = modify->find_fix(ids[i]); - if (f_id < 0) error->all(FLERR,"Fix ave/histo input is invalid fix"); - Fix *fix = modify->fix[f_id]; + } else if (val.which == ArgInfo::FIX) { + val.val.f = modify->get_fix_by_id(val.id); + if (!val.val.f) error->all(FLERR,"Fix ID {} for {} does not exist", val.id, mycmd); // fixes can produce multiple kinds of output - if (fix->scalar_flag || fix->vector_flag || fix->array_flag) + if (val.val.f->scalar_flag || val.val.f->vector_flag || val.val.f->array_flag) kindglobal = 1; - if (fix->peratom_flag) kindperatom = 1; - if (fix->local_flag) kindlocal = 1; + if (val.val.f->peratom_flag) kindperatom = 1; + if (val.val.f->local_flag) kindlocal = 1; - } else if (which[i] == ArgInfo::VARIABLE) { - int ivariable = input->variable->find(ids[i]); - if (ivariable < 0) - error->all(FLERR,"Fix ave/histo input is invalid variable"); + } else if (val.which == ArgInfo::VARIABLE) { + val.val.v = input->variable->find(val.id.c_str()); + if (val.val.v < 0) + error->all(FLERR,"Variable name {} for {} does not exist", val.id, mycmd); // variables only produce one kind of output - if (input->variable->equalstyle(ivariable)) kindglobal = 1; - else if (input->variable->atomstyle(ivariable)) kindperatom = 1; - else error->all(FLERR,"Fix ave/histo input is invalid kind of variable"); + if (input->variable->equalstyle(val.val.v)) kindglobal = 1; + else if (input->variable->atomstyle(val.val.v)) kindperatom = 1; + else error->all(FLERR,"{} variable {} is incompatible style", mycmd, val.id); } if (kind == DEFAULT) { if (kindglobal + kindperatom + kindlocal > 1) - error->all(FLERR,"Fix ave/histo input kind is ambiguous"); + error->all(FLERR,"{} input kind is ambiguous", mycmd); if (kindglobal) kind = GLOBAL; if (kindperatom) kind = PERATOM; if (kindlocal) kind = LOCAL; } else if (kind == GLOBAL) { if (!kindglobal) - error->all(FLERR,"Fix ave/histo input kind is invalid"); + error->all(FLERR,"{} input kind is not global", mycmd); } else if (kind == PERATOM) { if (!kindperatom) - error->all(FLERR,"Fix ave/histo input kind is invalid"); + error->all(FLERR,"{} input kind is not peratom", mycmd); } else if (kind == LOCAL) { if (!kindlocal) - error->all(FLERR,"Fix ave/histo input kind is invalid"); + error->all(FLERR,"{} input kind is not local", mycmd); } } @@ -245,183 +234,112 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : // for fix inputs, check that fix frequency is acceptable if (kind == PERATOM && mode == SCALAR) - error->all(FLERR, - "Fix ave/histo cannot input per-atom values in scalar mode"); + error->all(FLERR, "{} cannot process per-atom values in scalar mode", mycmd); if (kind == LOCAL && mode == SCALAR) - error->all(FLERR,"Fix ave/histo cannot input local values in scalar mode"); + error->all(FLERR,"{} cannot process local values in scalar mode", mycmd); - for (int i = 0; i < nvalues; i++) { - if (which[i] == ArgInfo::COMPUTE && kind == GLOBAL && mode == SCALAR) { - int icompute = modify->find_compute(ids[i]); - if (icompute < 0) - error->all(FLERR,"Compute ID for fix ave/histo does not exist"); - if (argindex[i] == 0 && modify->compute[icompute]->scalar_flag == 0) - error->all(FLERR, - "Fix ave/histo compute does not calculate a global scalar"); - if (argindex[i] && modify->compute[icompute]->vector_flag == 0) - error->all(FLERR, - "Fix ave/histo compute does not calculate a global vector"); - if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector) - error->all(FLERR, - "Fix ave/histo compute vector is accessed out-of-range"); + for (auto &val : values) { + if (val.which == ArgInfo::COMPUTE && kind == GLOBAL && mode == SCALAR) { + if (val.argindex == 0 && val.val.c->scalar_flag == 0) + error->all(FLERR, "{} compute {} does not calculate a global scalar", mycmd, val.id); + if (val.argindex && val.val.c->vector_flag == 0) + error->all(FLERR, "{} compute {} does not calculate a global vector", mycmd, val.id); + if (val.argindex && val.argindex > val.val.c->size_vector) + error->all(FLERR, "{} compute {} vector is accessed out-of-range", mycmd, val.id); - } else if (which[i] == ArgInfo::COMPUTE && kind == GLOBAL && mode == VECTOR) { - int icompute = modify->find_compute(ids[i]); - if (icompute < 0) - error->all(FLERR,"Compute ID for fix ave/histo does not exist"); - if (argindex[i] == 0 && modify->compute[icompute]->vector_flag == 0) - error->all(FLERR, - "Fix ave/histo compute does not calculate a global vector"); - if (argindex[i] && modify->compute[icompute]->array_flag == 0) - error->all(FLERR, - "Fix ave/histo compute does not calculate a global array"); - if (argindex[i] && - argindex[i] > modify->compute[icompute]->size_array_cols) - error->all(FLERR, - "Fix ave/histo compute array is accessed out-of-range"); + } else if (val.which == ArgInfo::COMPUTE && kind == GLOBAL && mode == VECTOR) { + if (val.argindex == 0 && val.val.c->vector_flag == 0) + error->all(FLERR, "{} compute {} does not calculate a global vector", mycmd, val.id); + if (val.argindex && val.val.c->array_flag == 0) + error->all(FLERR, "{} compute {} does not calculate a global array", mycmd, val.id); + if (val.argindex && val.argindex > val.val.c->size_array_cols) + error->all(FLERR, "{} compute {} array is accessed out-of-range", mycmd, val.id); - } else if (which[i] == ArgInfo::COMPUTE && kind == PERATOM) { - int icompute = modify->find_compute(ids[i]); - if (icompute < 0) - error->all(FLERR,"Compute ID for fix ave/histo does not exist"); - if (modify->compute[icompute]->peratom_flag == 0) - error->all(FLERR, - "Fix ave/histo compute does not calculate per-atom values"); - if (argindex[i] == 0 && - modify->compute[icompute]->size_peratom_cols != 0) - error->all(FLERR,"Fix ave/histo compute does not " - "calculate a per-atom vector"); - if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0) - error->all(FLERR,"Fix ave/histo compute does not " - "calculate a per-atom array"); - if (argindex[i] && - argindex[i] > modify->compute[icompute]->size_peratom_cols) - error->all(FLERR, - "Fix ave/histo compute array is accessed out-of-range"); + } else if (val.which == ArgInfo::COMPUTE && kind == PERATOM) { + if (val.val.c->peratom_flag == 0) + error->all(FLERR, "{} compute {} does not calculate per-atom values", mycmd, val.id); + if (val.argindex == 0 && val.val.c->size_peratom_cols != 0) + error->all(FLERR, "{} compute {} does not calculate a per-atom vector", mycmd, val.id); + if (val.argindex && val.val.c->size_peratom_cols == 0) + error->all(FLERR, "{} compute {} does not calculate a per-atom array", mycmd, val.id); + if (val.argindex && val.argindex > val.val.c->size_peratom_cols) + error->all(FLERR, "{} compute {} array is accessed out-of-range", mycmd, val.id); - } else if (which[i] == ArgInfo::COMPUTE && kind == LOCAL) { - int icompute = modify->find_compute(ids[i]); - if (icompute < 0) - error->all(FLERR,"Compute ID for fix ave/histo does not exist"); - if (modify->compute[icompute]->local_flag == 0) - error->all(FLERR, - "Fix ave/histo compute does not calculate local values"); - if (argindex[i] == 0 && - modify->compute[icompute]->size_local_cols != 0) - error->all(FLERR,"Fix ave/histo compute does not " - "calculate a local vector"); - if (argindex[i] && modify->compute[icompute]->size_local_cols == 0) - error->all(FLERR,"Fix ave/histo compute does not " - "calculate a local array"); - if (argindex[i] && - argindex[i] > modify->compute[icompute]->size_local_cols) - error->all(FLERR, - "Fix ave/histo compute array is accessed out-of-range"); + } else if (val.which == ArgInfo::COMPUTE && kind == LOCAL) { + if (val.val.c->local_flag == 0) + error->all(FLERR, "{} compute {} does not calculate local values", mycmd, val.id); + if (val.argindex == 0 && val.val.c->size_local_cols != 0) + error->all(FLERR, "{} compute {} does not calculate a local vector", mycmd, val.id); + if (val.argindex && val.val.c->size_local_cols == 0) + error->all(FLERR, "{} compute {} does not calculate a local array", mycmd, val.id); + if (val.argindex && val.argindex > val.val.c->size_local_cols) + error->all(FLERR, "{} compute {} array is accessed out-of-range", mycmd, val.id); - } else if (which[i] == ArgInfo::FIX && kind == GLOBAL && mode == SCALAR) { - int ifix = modify->find_fix(ids[i]); - if (ifix < 0) - error->all(FLERR,"Fix ID for fix ave/histo does not exist"); - if (argindex[i] == 0 && modify->fix[ifix]->scalar_flag == 0) - error->all(FLERR, - "Fix ave/histo fix does not calculate a global scalar"); - if (argindex[i] && modify->fix[ifix]->vector_flag == 0) - error->all(FLERR, - "Fix ave/histo fix does not calculate a global vector"); - if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector) - error->all(FLERR,"Fix ave/histo fix vector is accessed out-of-range"); - if (nevery % modify->fix[ifix]->global_freq) - error->all(FLERR, - "Fix for fix ave/histo not computed at compatible time"); + } else if (val.which == ArgInfo::FIX && kind == GLOBAL && mode == SCALAR) { + if (val.argindex == 0 && val.val.f->scalar_flag == 0) + error->all(FLERR, "{} fix {} does not calculate a global scalar", mycmd, val.id); + if (val.argindex && val.val.f->vector_flag == 0) + error->all(FLERR, "{} fix {} does not calculate a global vector", mycmd, val.id); + if (val.argindex && val.argindex > val.val.f->size_vector) + error->all(FLERR, "{} fix {} vector is accessed out-of-range", mycmd, val.id); + if (nevery % val.val.f->global_freq) + error->all(FLERR, "Fix {} for {} not computed at compatible time", val.id, mycmd); - } else if (which[i] == ArgInfo::FIX && kind == GLOBAL && mode == VECTOR) { - int ifix = modify->find_fix(ids[i]); - if (ifix < 0) - error->all(FLERR,"Fix ID for fix ave/histo does not exist"); - if (argindex[i] == 0 && modify->fix[ifix]->vector_flag == 0) - error->all(FLERR, - "Fix ave/histo fix does not calculate a global vector"); - if (argindex[i] && modify->fix[ifix]->array_flag == 0) - error->all(FLERR,"Fix ave/histo fix does not calculate a global array"); - if (argindex[i] && argindex[i] > modify->fix[ifix]->size_array_cols) - error->all(FLERR,"Fix ave/histo fix array is accessed out-of-range"); - if (nevery % modify->fix[ifix]->global_freq) - error->all(FLERR, - "Fix for fix ave/histo not computed at compatible time"); + } else if (val.which == ArgInfo::FIX && kind == GLOBAL && mode == VECTOR) { + if (val.argindex == 0 && val.val.f->vector_flag == 0) + error->all(FLERR, "{} fix {} does not calculate a global vector", mycmd, val.id); + if (val.argindex && val.val.f->array_flag == 0) + error->all(FLERR, "{} fix {} does not calculate a global array", mycmd, val.id); + if (val.argindex && val.argindex > val.val.f->size_array_cols) + error->all(FLERR, "{} fix {} array is accessed out-of-range", mycmd, val.id); + if (nevery % val.val.f->global_freq) + error->all(FLERR, "Fix {} for {} not computed at compatible time", val.id, mycmd); - } else if (which[i] == ArgInfo::FIX && kind == PERATOM) { - int ifix = modify->find_fix(ids[i]); - if (ifix < 0) - error->all(FLERR,"Fix ID for fix ave/histo does not exist"); - if (modify->fix[ifix]->peratom_flag == 0) - error->all(FLERR, - "Fix ave/histo fix does not calculate per-atom values"); - if (argindex[i] == 0 && - modify->fix[ifix]->size_peratom_cols != 0) - error->all(FLERR,"Fix ave/histo fix does not " - "calculate a per-atom vector"); - if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0) - error->all(FLERR,"Fix ave/histo fix does not " - "calculate a per-atom array"); - if (argindex[i] && - argindex[i] > modify->fix[ifix]->size_peratom_cols) - error->all(FLERR,"Fix ave/histo fix array is accessed out-of-range"); - if (nevery % modify->fix[ifix]->global_freq) - error->all(FLERR, - "Fix for fix ave/histo not computed at compatible time"); + } else if (val.which == ArgInfo::FIX && kind == PERATOM) { + if (val.val.f->peratom_flag == 0) + error->all(FLERR, "{} fix {} does not calculate per-atom values", mycmd, val.id); + if (val.argindex == 0 && val.val.f->size_peratom_cols != 0) + error->all(FLERR," {} fix {} does not calculate a per-atom vector", mycmd, val.id); + if (val.argindex && val.val.f->size_peratom_cols == 0) + error->all(FLERR, "{} fix {} does not ""calculate a per-atom array", mycmd, val.id); + if (val.argindex && val.argindex > val.val.f->size_peratom_cols) + error->all(FLERR, "{} fix {} array is accessed out-of-range", mycmd, val.id); + if (nevery % val.val.f->global_freq) + error->all(FLERR, "Fix {} for {} not computed at compatible time", val.id, mycmd); - } else if (which[i] == ArgInfo::FIX && kind == LOCAL) { - int ifix = modify->find_fix(ids[i]); - if (ifix < 0) - error->all(FLERR,"Fix ID for fix ave/histo does not exist"); - if (modify->fix[ifix]->local_flag == 0) - error->all(FLERR,"Fix ave/histo fix does not calculate local values"); - if (argindex[i] == 0 && - modify->fix[ifix]->size_local_cols != 0) - error->all(FLERR,"Fix ave/histo fix does not " - "calculate a local vector"); - if (argindex[i] && modify->fix[ifix]->size_local_cols == 0) - error->all(FLERR,"Fix ave/histo fix does not " - "calculate a local array"); - if (argindex[i] && - argindex[i] > modify->fix[ifix]->size_local_cols) - error->all(FLERR,"Fix ave/histo fix array is accessed out-of-range"); - if (nevery % modify->fix[ifix]->global_freq) - error->all(FLERR, - "Fix for fix ave/histo not computed at compatible time"); + } else if (val.which == ArgInfo::FIX && kind == LOCAL) { + if (val.val.f->local_flag == 0) + error->all(FLERR, "{} fix {} does not calculate local values", mycmd, val.id); + if (val.argindex == 0 && val.val.f->size_local_cols != 0) + error->all(FLERR, "{} fix {} does not calculate a local vector", mycmd, val.id); + if (val.argindex && val.val.f->size_local_cols == 0) + error->all(FLERR, "{} fix does not calculate a local array", mycmd, val.id); + if (val.argindex && val.argindex > val.val.f->size_local_cols) + error->all(FLERR, "{} fix {} array is accessed out-of-range", mycmd, val.id); + if (nevery % val.val.f->global_freq) + error->all(FLERR, "Fix {} for {} not computed at compatible time", val.id, mycmd); - } else if (which[i] == ArgInfo::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 (val.which == ArgInfo::VARIABLE && kind == GLOBAL && mode == SCALAR) { + if (val.argindex == 0 && input->variable->equalstyle(val.val.v) == 0) + error->all(FLERR,"{} variable {} is not equal-style variable", mycmd, val.id); + if (val.argindex && input->variable->vectorstyle(val.val.v) == 0) + error->all(FLERR,"{} variable {} is not vector-style variable" , mycmd, val.id); - } else if (which[i] == ArgInfo::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 (val.which == ArgInfo::VARIABLE && kind == GLOBAL && mode == VECTOR) { + if (val.argindex == 0 && input->variable->vectorstyle(val.val.v) == 0) + error->all(FLERR,"{} variable {} is not vector-style variable", mycmd, val.id); + if (val.argindex) error->all(FLERR,"{} variable {} cannot be indexed", mycmd, val.id); - } else if (which[i] == ArgInfo::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"); + } else if (val.which == ArgInfo::VARIABLE && kind == PERATOM) { + if (val.argindex == 0 && input->variable->atomstyle(val.val.v) == 0) + error->all(FLERR,"{} variable {} is not atom-style variable", mycmd, val.id); + if (val.argindex) error->all(FLERR,"{} variable {} cannot be indexed", mycmd, val.id); } } // print file comment lines - if (fp && me == 0) { + if (fp && comm->me == 0) { clearerr(fp); if (title1) fprintf(fp,"%s\n",title1); else fprintf(fp,"# Histogrammed data for fix %s\n",id); @@ -431,9 +349,7 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : if (title3) fprintf(fp,"%s\n",title3); else fprintf(fp,"# Bin Coord Count Count/Total\n"); - if (ferror(fp)) - error->one(FLERR,"Error writing file header"); - + if (ferror(fp)) error->one(FLERR,"Error writing file header: {}", utils::getsyserror()); filepos = platform::ftell(fp); } @@ -502,13 +418,7 @@ FixAveHisto::FixAveHisto(LAMMPS *lmp, int narg, char **arg) : FixAveHisto::~FixAveHisto() { - delete[] which; - delete[] argindex; - delete[] value2index; - for (int i = 0; i < nvalues; i++) delete[] ids[i]; - delete[] ids; - - if (fp && me == 0) fclose(fp); + if (fp && comm->me == 0) fclose(fp); delete[] bin; delete[] bin_total; @@ -532,26 +442,20 @@ int FixAveHisto::setmask() void FixAveHisto::init() { - // set current indices for all computes,fixes,variables + auto mycmd = fmt::format("fix {}", style); - for (int i = 0; i < nvalues; i++) { - if (which[i] == ArgInfo::COMPUTE) { - int icompute = modify->find_compute(ids[i]); - if (icompute < 0) - error->all(FLERR,"Compute ID for fix ave/histo does not exist"); - value2index[i] = icompute; + // update indices/pointers for all computes,fixes,variables - } else if (which[i] == ArgInfo::FIX) { - int ifix = modify->find_fix(ids[i]); - if (ifix < 0) - error->all(FLERR,"Fix ID for fix ave/histo does not exist"); - value2index[i] = ifix; - - } else if (which[i] == ArgInfo::VARIABLE) { - int ivariable = input->variable->find(ids[i]); - if (ivariable < 0) - error->all(FLERR,"Variable name for fix ave/histo does not exist"); - value2index[i] = ivariable; + for (auto &val : values) { + if (val.which == ArgInfo::COMPUTE) { + val.val.c = modify->get_compute_by_id(val.id); + if (!val.val.c) error->all(FLERR,"Compute ID {} for {} does not exist", val.id, mycmd); + } else if (val.which == ArgInfo::FIX) { + val.val.f = modify->get_fix_by_id(val.id); + if (!val.val.f) error->all(FLERR,"Fix ID {} for {} does not exist", val.id, mycmd); + } else if (val.which == ArgInfo::VARIABLE) { + val.val.v = input->variable->find(val.id.c_str()); + if (val.val.v < 0) error->all(FLERR,"Variable name {} for {} does not exist", val.id, mycmd); } } @@ -577,8 +481,6 @@ void FixAveHisto::setup(int /*vflag*/) void FixAveHisto::end_of_step() { - int i,j,m; - // skip if not step which requires doing something bigint ntimestep = update->ntimestep; @@ -591,7 +493,7 @@ void FixAveHisto::end_of_step() stats[0] = stats[1] = 0.0; stats[2] = BIG; stats[3] = -BIG; - for (i = 0; i < nbins; i++) bin[i] = 0.0; + for (int i = 0; i < nbins; i++) bin[i] = 0.0; } // accumulate results of computes,fixes,variables to local copy @@ -599,132 +501,128 @@ void FixAveHisto::end_of_step() modify->clearstep_compute(); - for (i = 0; i < nvalues; i++) { - m = value2index[i]; - j = argindex[i]; + for (auto &val : values) { + int j = val.argindex; // atom attributes - if (which[i] == ArgInfo::X) + if (val.which == ArgInfo::X) bin_atoms(&atom->x[0][j],3); - else if (which[i] == ArgInfo::V) + else if (val.which == ArgInfo::V) bin_atoms(&atom->v[0][j],3); - else if (which[i] == ArgInfo::F) + else if (val.which == ArgInfo::F) bin_atoms(&atom->f[0][j],3); // invoke compute if not previously invoked - if (which[i] == ArgInfo::COMPUTE) { - Compute *compute = modify->compute[m]; + if (val.which == ArgInfo::COMPUTE) { if (kind == GLOBAL && mode == SCALAR) { if (j == 0) { - if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) { - compute->compute_scalar(); - compute->invoked_flag |= Compute::INVOKED_SCALAR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_SCALAR)) { + val.val.c->compute_scalar(); + val.val.c->invoked_flag |= Compute::INVOKED_SCALAR; } - bin_one(compute->scalar); + bin_one(val.val.c->scalar); } else { - if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { - compute->compute_vector(); - compute->invoked_flag |= Compute::INVOKED_VECTOR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) { + val.val.c->compute_vector(); + val.val.c->invoked_flag |= Compute::INVOKED_VECTOR; } - bin_one(compute->vector[j-1]); + bin_one(val.val.c->vector[j-1]); } } else if (kind == GLOBAL && mode == VECTOR) { if (j == 0) { - if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { - compute->compute_vector(); - compute->invoked_flag |= Compute::INVOKED_VECTOR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) { + val.val.c->compute_vector(); + val.val.c->invoked_flag |= Compute::INVOKED_VECTOR; } - bin_vector(compute->size_vector,compute->vector,1); + bin_vector(val.val.c->size_vector,val.val.c->vector,1); } else { - if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { - compute->compute_array(); - compute->invoked_flag |= Compute::INVOKED_ARRAY; + if (!(val.val.c->invoked_flag & Compute::INVOKED_ARRAY)) { + val.val.c->compute_array(); + val.val.c->invoked_flag |= Compute::INVOKED_ARRAY; } - if (compute->array) - bin_vector(compute->size_array_rows,&compute->array[0][j-1], - compute->size_array_cols); + if (val.val.c->array) + bin_vector(val.val.c->size_array_rows,&val.val.c->array[0][j-1], + val.val.c->size_array_cols); } } else if (kind == PERATOM) { - if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { - compute->compute_peratom(); - compute->invoked_flag |= Compute::INVOKED_PERATOM; + if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) { + val.val.c->compute_peratom(); + val.val.c->invoked_flag |= Compute::INVOKED_PERATOM; } if (j == 0) - bin_atoms(compute->vector_atom,1); - else if (compute->array_atom) - bin_atoms(&compute->array_atom[0][j-1],compute->size_peratom_cols); + bin_atoms(val.val.c->vector_atom,1); + else if (val.val.c->array_atom) + bin_atoms(&val.val.c->array_atom[0][j-1],val.val.c->size_peratom_cols); } else if (kind == LOCAL) { - if (!(compute->invoked_flag & Compute::INVOKED_LOCAL)) { - compute->compute_local(); - compute->invoked_flag |= Compute::INVOKED_LOCAL; + if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) { + val.val.c->compute_local(); + val.val.c->invoked_flag |= Compute::INVOKED_LOCAL; } if (j == 0) - bin_vector(compute->size_local_rows,compute->vector_local,1); - else if (compute->array_local) - bin_vector(compute->size_local_rows,&compute->array_local[0][j-1], - compute->size_local_cols); + bin_vector(val.val.c->size_local_rows,val.val.c->vector_local,1); + else if (val.val.c->array_local) + bin_vector(val.val.c->size_local_rows,&val.val.c->array_local[0][j-1], + val.val.c->size_local_cols); } // access fix fields, guaranteed to be ready - } else if (which[i] == ArgInfo::FIX) { - - Fix *fix = modify->fix[m]; + } else if (val.which == ArgInfo::FIX) { if (kind == GLOBAL && mode == SCALAR) { - if (j == 0) bin_one(fix->compute_scalar()); - else bin_one(fix->compute_vector(j-1)); + if (j == 0) bin_one(val.val.f->compute_scalar()); + else bin_one(val.val.f->compute_vector(j-1)); } else if (kind == GLOBAL && mode == VECTOR) { if (j == 0) { - int n = fix->size_vector; - for (i = 0; i < n; i++) bin_one(fix->compute_vector(i)); + int n = val.val.f->size_vector; + for (int i = 0; i < n; i++) bin_one(val.val.f->compute_vector(i)); } else { - int n = fix->size_vector; - for (i = 0; i < n; i++) bin_one(fix->compute_array(i,j-1)); + int n = val.val.f->size_vector; + for (int i = 0; i < n; i++) bin_one(val.val.f->compute_array(i,j-1)); } } else if (kind == PERATOM) { - if (j == 0) bin_atoms(fix->vector_atom,1); - else if (fix->array_atom) - bin_atoms(fix->array_atom[j-1],fix->size_peratom_cols); + if (j == 0) bin_atoms(val.val.f->vector_atom,1); + else if (val.val.f->array_atom) + bin_atoms(val.val.f->array_atom[j-1],val.val.f->size_peratom_cols); } else if (kind == LOCAL) { - if (j == 0) bin_vector(fix->size_local_rows,fix->vector_local,1); - else if (fix->array_local) - bin_vector(fix->size_local_rows,&fix->array_local[0][j-1], - fix->size_local_cols); + if (j == 0) bin_vector(val.val.f->size_local_rows,val.val.f->vector_local,1); + else if (val.val.f->array_local) + bin_vector(val.val.f->size_local_rows,&val.val.f->array_local[0][j-1], + val.val.f->size_local_cols); } // evaluate equal-style or vector-style or atom-style variable - } else if (which[i] == ArgInfo::VARIABLE) { + } else if (val.which == ArgInfo::VARIABLE) { if (kind == GLOBAL && mode == SCALAR) { - if (j == 0) bin_one(input->variable->compute_equal(m)); + if (j == 0) bin_one(input->variable->compute_equal(val.val.v)); else { double *varvec; - int nvec = input->variable->compute_vector(m,&varvec); + int nvec = input->variable->compute_vector(val.val.v,&varvec); if (nvec < j) bin_one(0.0); else bin_one(varvec[j-1]); } } else if (kind == GLOBAL && mode == VECTOR) { double *varvec; - int nvec = input->variable->compute_vector(m,&varvec); + int nvec = input->variable->compute_vector(val.val.v,&varvec); bin_vector(nvec,varvec,1); - } else if (which[i] == ArgInfo::VARIABLE && kind == PERATOM) { + } else if (val.which == ArgInfo::VARIABLE && kind == PERATOM) { if (atom->nmax > maxatom) { memory->destroy(vector); maxatom = atom->nmax; memory->create(vector,maxatom,"ave/histo:vector"); } - input->variable->compute_atom(m,igroup,vector,1,0); + input->variable->compute_atom(val.val.v,igroup,vector,1,0); bin_atoms(vector,1); } } @@ -756,7 +654,7 @@ void FixAveHisto::end_of_step() stats[1] = stats_all[1]; stats[2] = stats_all[2]; stats[3] = stats_all[3]; - for (i = 0; i < nbins; i++) bin[i] = bin_all[i]; + for (int i = 0; i < nbins; i++) bin[i] = bin_all[i]; } // if ave = ONE, only single Nfreq timestep value is needed @@ -768,16 +666,17 @@ void FixAveHisto::end_of_step() stats_total[1] = stats[1]; stats_total[2] = stats[2]; stats_total[3] = stats[3]; - for (i = 0; i < nbins; i++) bin_total[i] = bin[i]; + for (int i = 0; i < nbins; i++) bin_total[i] = bin[i]; } else if (ave == RUNNING) { stats_total[0] += stats[0]; stats_total[1] += stats[1]; stats_total[2] = MIN(stats_total[2],stats[2]); stats_total[3] = MAX(stats_total[3],stats[3]); - for (i = 0; i < nbins; i++) bin_total[i] += bin[i]; + for (int i = 0; i < nbins; i++) bin_total[i] += bin[i]; } else if (ave == WINDOW) { + int m; stats_total[0] += stats[0]; if (window_limit) stats_total[0] -= stats_list[iwindow][0]; stats_list[iwindow][0] = stats[0]; @@ -790,14 +689,14 @@ void FixAveHisto::end_of_step() stats_list[iwindow][2] = stats[2]; stats_total[2] = stats_list[0][2]; - for (i = 1; i < m; i++) + for (int i = 1; i < m; i++) stats_total[2] = MIN(stats_total[2],stats_list[i][2]); stats_list[iwindow][3] = stats[3]; stats_total[3] = stats_list[0][3]; - for (i = 1; i < m; i++) + for (int i = 1; i < m; i++) stats_total[3] = MAX(stats_total[3],stats_list[i][3]); - for (i = 0; i < nbins; i++) { + for (int i = 0; i < nbins; i++) { bin_total[i] += bin[i]; if (window_limit) bin_total[i] -= bin_list[iwindow][i]; bin_list[iwindow][i] = bin[i]; @@ -812,17 +711,17 @@ void FixAveHisto::end_of_step() // output result to file - if (fp && me == 0) { + if (fp && comm->me == 0) { clearerr(fp); if (overwrite) platform::fseek(fp,filepos); fmt::print(fp,"{} {} {} {} {} {}\n",ntimestep,nbins, stats_total[0],stats_total[1],stats_total[2],stats_total[3]); if (stats_total[0] != 0.0) - for (i = 0; i < nbins; i++) + for (int i = 0; i < nbins; i++) fprintf(fp,"%d %g %g %g\n", i+1,coord[i],bin_total[i],bin_total[i]/stats_total[0]); else - for (i = 0; i < nbins; i++) + for (int i = 0; i < nbins; i++) fprintf(fp,"%d %g %g %g\n",i+1,coord[i],0.0,0.0); if (ferror(fp)) @@ -937,73 +836,74 @@ void FixAveHisto::options(int iarg, int narg, char **arg) title3 = nullptr; // optional args - + auto mycmd = fmt::format("fix {}", style); + while (iarg < narg) { if (strcmp(arg[iarg],"file") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); - if (me == 0) { + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " file", error); + if (comm->me == 0) { fp = fopen(arg[iarg+1],"w"); if (fp == nullptr) - error->one(FLERR,"Cannot open fix ave/histo file {}: {}", - arg[iarg+1], utils::getsyserror()); + error->one(FLERR, "Cannot open fix ave/histo file {}: {}", + arg[iarg+1], utils::getsyserror()); } iarg += 2; } else if (strcmp(arg[iarg],"kind") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " kind", error); if (strcmp(arg[iarg+1],"global") == 0) kind = GLOBAL; else if (strcmp(arg[iarg+1],"peratom") == 0) kind = PERATOM; else if (strcmp(arg[iarg+1],"local") == 0) kind = LOCAL; - else error->all(FLERR,"Illegal fix ave/histo command"); + else error->all(FLERR,"Unknown fix ave/histo kind option: {}", arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"ave") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " ave", error); if (strcmp(arg[iarg+1],"one") == 0) ave = ONE; else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING; else if (strcmp(arg[iarg+1],"window") == 0) ave = WINDOW; - else error->all(FLERR,"Illegal fix ave/histo command"); + else error->all(FLERR,"Unknown fix ave/histo ave option: {}", arg[iarg+1]); if (ave == WINDOW) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+3 > narg) utils::missing_cmd_args(FLERR, mycmd + " ave window", error); nwindow = utils::inumeric(FLERR,arg[iarg+2],false,lmp); - if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/histo command"); + if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/histo ave window size: {}", nwindow); } iarg += 2; if (ave == WINDOW) iarg++; } else if (strcmp(arg[iarg],"start") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " start", error); startstep = utils::inumeric(FLERR,arg[iarg+1],false,lmp); iarg += 2; } else if (strcmp(arg[iarg],"mode") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " mode", error); if (strcmp(arg[iarg+1],"scalar") == 0) mode = SCALAR; else if (strcmp(arg[iarg+1],"vector") == 0) mode = VECTOR; - else error->all(FLERR,"Illegal fix ave/histo command"); + else error->all(FLERR,"Unknown fix ave/histo mode option: {}", arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"beyond") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " beyond", error); if (strcmp(arg[iarg+1],"ignore") == 0) beyond = IGNORE; else if (strcmp(arg[iarg+1],"end") == 0) beyond = END; else if (strcmp(arg[iarg+1],"extra") == 0) beyond = EXTRA; - else error->all(FLERR,"Illegal fix ave/histo command"); + else error->all(FLERR,"Unknown fix ave/histo beyond option: {}", arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"overwrite") == 0) { overwrite = 1; iarg += 1; } else if (strcmp(arg[iarg],"title1") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " title1", error); delete[] title1; title1 = utils::strdup(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"title2") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " title2", error); delete[] title2; title2 = utils::strdup(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"title3") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/histo command"); + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, mycmd + " title3", error); delete[] title3; title3 = utils::strdup(arg[iarg+1]); iarg += 2; - } else error->all(FLERR,"Illegal fix ave/histo command"); + } else error->all(FLERR,"Unknown {} option: {}", mycmd, arg[iarg]); } } diff --git a/src/fix_ave_histo.h b/src/fix_ave_histo.h index ed64c8562d..620fb3b321 100644 --- a/src/fix_ave_histo.h +++ b/src/fix_ave_histo.h @@ -36,11 +36,21 @@ class FixAveHisto : public Fix { double compute_array(int, int) override; protected: - int me, nvalues; - int nrepeat, nfreq, irepeat; + struct value_t { + int which; // type of data: COMPUTE, FIX, VARIABLE + int argindex; // 1-based index if data is vector, else 0 + std::string id; // compute/fix/variable ID + union { + class Compute *c; + class Fix *f; + int v; + } val; + }; + std::vector values; + + int nvalues, nrepeat, nfreq, irepeat; bigint nvalid, nvalid_last; - int *which, *argindex, *value2index; - char **ids; + FILE *fp; double lo, hi, binsize, bininv; int kind, beyond, overwrite; diff --git a/src/fix_ave_histo_weight.cpp b/src/fix_ave_histo_weight.cpp index f90551b52f..8051d93b7c 100644 --- a/src/fix_ave_histo_weight.cpp +++ b/src/fix_ave_histo_weight.cpp @@ -19,6 +19,7 @@ #include "arg_info.h" #include "atom.h" +#include "comm.h" #include "compute.h" #include "error.h" #include "fix.h" @@ -31,68 +32,61 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{ONE,RUNNING}; -enum{SCALAR,VECTOR,WINDOW}; -enum{DEFAULT,GLOBAL,PERATOM,LOCAL}; -enum{IGNORE,END,EXTRA}; -enum{SINGLE,VALUE}; +enum { ONE, RUNNING }; +enum { SCALAR, VECTOR, WINDOW }; +enum { DEFAULT, GLOBAL, PERATOM, LOCAL }; +enum { IGNORE, END, EXTRA }; +enum { SINGLE, VALUE }; #define BIG 1.0e20 /* ---------------------------------------------------------------------- */ FixAveHistoWeight::FixAveHistoWeight(LAMMPS *lmp, int narg, char **arg) : - FixAveHisto(lmp, narg, arg) + FixAveHisto(lmp, narg, arg) { // nvalues = 2 required for histo/weight - if (nvalues != 2) error->all(FLERR,"Illegal fix ave/histo/weight command"); + if (nvalues != 2) + error->all(FLERR, "Illegal fix ave/histo/weight command: must have two data arguments"); // check that length of 2 values is the same - int size[2] = {0,0}; + int size[2] = {0, 0}; for (int i = 0; i < nvalues; i++) { - if (which[i] == ArgInfo::X || which[i] == ArgInfo::V || which[i] == ArgInfo::F) { + auto &val = values[i]; + if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F) { size[i] = atom->nlocal; - } else if (which[i] == ArgInfo::COMPUTE && kind == GLOBAL && mode == SCALAR) { - int icompute = modify->find_compute(ids[i]); - size[i] = modify->compute[icompute]->size_vector; - } else if (which[i] == ArgInfo::COMPUTE && kind == GLOBAL && mode == VECTOR) { - int icompute = modify->find_compute(ids[i]); - size[i] = modify->compute[icompute]->size_array_rows; - } else if (which[i] == ArgInfo::COMPUTE && kind == PERATOM) { + } else if (val.which == ArgInfo::COMPUTE && kind == GLOBAL && mode == SCALAR) { + size[i] = val.val.c->size_vector; + } else if (val.which == ArgInfo::COMPUTE && kind == GLOBAL && mode == VECTOR) { + size[i] = val.val.c->size_array_rows; + } else if (val.which == ArgInfo::COMPUTE && kind == PERATOM) { size[i] = atom->nlocal; - } else if (which[i] == ArgInfo::COMPUTE && kind == LOCAL) { - int icompute = modify->find_compute(ids[i]); - size[i] = modify->compute[icompute]->size_local_rows; - } else if (which[i] == ArgInfo::FIX && kind == GLOBAL && mode == SCALAR) { - int ifix = modify->find_fix(ids[i]); - size[i] = modify->fix[ifix]->size_vector; - } else if (which[i] == ArgInfo::FIX && kind == GLOBAL && mode == VECTOR) { - int ifix = modify->find_fix(ids[i]); - size[i]= modify->fix[ifix]->size_array_rows; - } else if (which[i] == ArgInfo::FIX && kind == PERATOM) { + } else if (val.which == ArgInfo::COMPUTE && kind == LOCAL) { + size[i] = val.val.c->size_local_rows; + } else if (val.which == ArgInfo::FIX && kind == GLOBAL && mode == SCALAR) { + size[i] = val.val.f->size_vector; + } else if (val.which == ArgInfo::FIX && kind == GLOBAL && mode == VECTOR) { + size[i]= val.val.f->size_array_rows; + } else if (val.which == ArgInfo::FIX && kind == PERATOM) { size[i] = atom->nlocal; - } else if (which[i] == ArgInfo::FIX && kind == LOCAL) { - int ifix = modify->find_fix(ids[i]); - size[i] = modify->fix[ifix]->size_local_rows; - } else if (which[i] == ArgInfo::VARIABLE && kind == PERATOM) { + } else if (val.which == ArgInfo::FIX && kind == LOCAL) { + size[i] = val.val.f->size_local_rows; + } else if (val.which == ArgInfo::VARIABLE && kind == PERATOM) { size[i] = atom->nlocal; } } if (size[0] != size[1]) - error->all(FLERR,"Fix ave/histo/weight value and weight vector " - "lengths do not match"); + error->all(FLERR,"Fix ave/histo/weight value and weight vector lengths do not match"); } /* ---------------------------------------------------------------------- */ void FixAveHistoWeight::end_of_step() { - int i,j,m; - // skip if not step which requires doing something bigint ntimestep = update->ntimestep; @@ -105,7 +99,7 @@ void FixAveHistoWeight::end_of_step() stats[0] = stats[1] = 0.0; stats[2] = BIG; stats[3] = -BIG; - for (i = 0; i < nbins; i++) bin[i] = 0.0; + for (int i = 0; i < nbins; i++) bin[i] = 0.0; } // first calculate weight factors, then bin single value @@ -119,265 +113,256 @@ void FixAveHistoWeight::end_of_step() double weight = 0.0; double *weights = nullptr; int stride = 0; - i = 1; - - m = value2index[i]; - j = argindex[i]; + auto &val = values[1]; + int j = val.argindex; // atom attributes - if (which[i] == ArgInfo::X) { + if (val.which == ArgInfo::X) { weights = &atom->x[0][j]; stride = 3; - } else if (which[i] == ArgInfo::V) { + } else if (val.which == ArgInfo::V) { weights = &atom->v[0][j]; stride = 3; bin_atoms(&atom->v[0][j],3); - } else if (which[i] == ArgInfo::F) { + } else if (val.which == ArgInfo::F) { weights = &atom->f[0][j]; stride = 3; } // invoke compute if not previously invoked - if (which[i] == ArgInfo::COMPUTE) { - - Compute *compute = modify->compute[m]; + if (val.which == ArgInfo::COMPUTE) { if (kind == GLOBAL && mode == SCALAR) { if (j == 0) { - if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) { - compute->compute_scalar(); - compute->invoked_flag |= Compute::INVOKED_SCALAR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_SCALAR)) { + val.val.c->compute_scalar(); + val.val.c->invoked_flag |= Compute::INVOKED_SCALAR; } - weight = compute->scalar; + weight = val.val.c->scalar; } else { - if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { - compute->compute_vector(); - compute->invoked_flag |= Compute::INVOKED_VECTOR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) { + val.val.c->compute_vector(); + val.val.c->invoked_flag |= Compute::INVOKED_VECTOR; } - weight = compute->vector[j-1]; + weight = val.val.c->vector[j-1]; } } else if (kind == GLOBAL && mode == VECTOR) { if (j == 0) { - if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { - compute->compute_vector(); - compute->invoked_flag |= Compute::INVOKED_VECTOR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) { + val.val.c->compute_vector(); + val.val.c->invoked_flag |= Compute::INVOKED_VECTOR; } - weights = compute->vector; + weights = val.val.c->vector; stride = 1; } else { - if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { - compute->compute_array(); - compute->invoked_flag |= Compute::INVOKED_ARRAY; + if (!(val.val.c->invoked_flag & Compute::INVOKED_ARRAY)) { + val.val.c->compute_array(); + val.val.c->invoked_flag |= Compute::INVOKED_ARRAY; } - if (compute->array) weights = &compute->array[0][j-1]; - stride = compute->size_array_cols; + if (val.val.c->array) weights = &val.val.c->array[0][j-1]; + stride = val.val.c->size_array_cols; } } else if (kind == PERATOM) { - if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { - compute->compute_peratom(); - compute->invoked_flag |= Compute::INVOKED_PERATOM; + if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) { + val.val.c->compute_peratom(); + val.val.c->invoked_flag |= Compute::INVOKED_PERATOM; } if (j == 0) { - weights = compute->vector_atom; + weights = val.val.c->vector_atom; stride = 1; - } else if (compute->array_atom) { - weights = &compute->array_atom[0][j-1]; - stride = compute->size_peratom_cols; + } else if (val.val.c->array_atom) { + weights = &val.val.c->array_atom[0][j-1]; + stride = val.val.c->size_peratom_cols; } } else if (kind == LOCAL) { - if (!(compute->invoked_flag & Compute::INVOKED_LOCAL)) { - compute->compute_local(); - compute->invoked_flag |= Compute::INVOKED_LOCAL; + if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) { + val.val.c->compute_local(); + val.val.c->invoked_flag |= Compute::INVOKED_LOCAL; } if (j == 0) { - weights = compute->vector_local; + weights = val.val.c->vector_local; stride = 1; - } else if (compute->array_local) { - weights = &compute->array_local[0][j-1]; - stride = compute->size_local_cols; + } else if (val.val.c->array_local) { + weights = &val.val.c->array_local[0][j-1]; + stride = val.val.c->size_local_cols; } } // access fix fields, guaranteed to be ready - } else if (which[i] == ArgInfo::FIX) { - - Fix *fix = modify->fix[m]; + } else if (val.which == ArgInfo::FIX) { if (kind == GLOBAL && mode == SCALAR) { - if (j == 0) weight = fix->compute_scalar(); - else weight = fix->compute_vector(j-1); + if (j == 0) weight = val.val.f->compute_scalar(); + else weight = val.val.f->compute_vector(j-1); } else if (kind == GLOBAL && mode == VECTOR) { error->all(FLERR,"Fix ave/histo/weight option not yet supported"); // NOTE: need to allocate local storage if (j == 0) { - int n = fix->size_vector; - for (i = 0; i < n; i++) weights[n] = fix->compute_vector(i); + int n = val.val.f->size_vector; + for (int i = 0; i < n; i++) weights[n] = val.val.f->compute_vector(i); } else { - int n = fix->size_vector; - for (i = 0; i < n; i++) weights[n] = fix->compute_array(i,j-1); + int n = val.val.f->size_vector; + for (int i = 0; i < n; i++) weights[n] = val.val.f->compute_array(i,j-1); } } else if (kind == PERATOM) { if (j == 0) { - weights = fix->vector_atom; + weights = val.val.f->vector_atom; stride = 1; - } else if (fix->array_atom) { - weights = fix->array_atom[j-1]; - stride = fix->size_peratom_cols; + } else if (val.val.f->array_atom) { + weights = val.val.f->array_atom[j-1]; + stride = val.val.f->size_peratom_cols; } } else if (kind == LOCAL) { if (j == 0) { - weights = fix->vector_local; + weights = val.val.f->vector_local; stride = 1; - } else if (fix->array_local) { - weights = &fix->array_local[0][j-1]; - stride = fix->size_local_cols; + } else if (val.val.f->array_local) { + weights = &val.val.f->array_local[0][j-1]; + stride = val.val.f->size_local_cols; } } // evaluate equal-style variable - } else if (which[i] == ArgInfo::VARIABLE && kind == GLOBAL) { - weight = input->variable->compute_equal(m); + } else if (val.which == ArgInfo::VARIABLE && kind == GLOBAL) { + weight = input->variable->compute_equal(val.val.v); - } else if (which[i] == ArgInfo::VARIABLE && kind == PERATOM) { + } else if (val.which == ArgInfo::VARIABLE && kind == PERATOM) { if (atom->nmax > maxatom) { memory->destroy(vector); maxatom = atom->nmax; memory->create(vector,maxatom,"ave/histo/weight:vector"); } - input->variable->compute_atom(m,igroup,vector,1,0); + input->variable->compute_atom(val.val.v,igroup,vector,1,0); weights = vector; stride = 1; } // bin values using weights, values are 1st value (i = 0) - i = 0; - m = value2index[i]; - j = argindex[i]; + val = values[0]; + j = val.argindex; // atom attributes - if (which[i] == ArgInfo::X && weights != nullptr) + if (val.which == ArgInfo::X && weights != nullptr) bin_atoms_weights(&atom->x[0][j],3,weights,stride); - else if (which[i] == ArgInfo::V && weights != nullptr) + else if (val.which == ArgInfo::V && weights != nullptr) bin_atoms_weights(&atom->v[0][j],3,weights,stride); - else if (which[i] == ArgInfo::F && weights != nullptr) + else if (val.which == ArgInfo::F && weights != nullptr) bin_atoms_weights(&atom->f[0][j],3,weights,stride); // invoke compute if not previously invoked - if (which[i] == ArgInfo::COMPUTE) { - Compute *compute = modify->compute[m]; + if (val.which == ArgInfo::COMPUTE) { + if (kind == GLOBAL && mode == SCALAR) { if (j == 0) { - if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) { - compute->compute_scalar(); - compute->invoked_flag |= Compute::INVOKED_SCALAR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_SCALAR)) { + val.val.c->compute_scalar(); + val.val.c->invoked_flag |= Compute::INVOKED_SCALAR; } - bin_one_weights(compute->scalar,weight); + bin_one_weights(val.val.c->scalar,weight); } else { - if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { - compute->compute_vector(); - compute->invoked_flag |= Compute::INVOKED_VECTOR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) { + val.val.c->compute_vector(); + val.val.c->invoked_flag |= Compute::INVOKED_VECTOR; } - bin_one_weights(compute->vector[j-1],weight); + bin_one_weights(val.val.c->vector[j-1],weight); } } else if (kind == GLOBAL && mode == VECTOR) { if (j == 0) { - if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { - compute->compute_vector(); - compute->invoked_flag |= Compute::INVOKED_VECTOR; + if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) { + val.val.c->compute_vector(); + val.val.c->invoked_flag |= Compute::INVOKED_VECTOR; } - bin_vector_weights(compute->size_vector,compute->vector,1, + bin_vector_weights(val.val.c->size_vector,val.val.c->vector,1, weights,stride); } else { - if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { - compute->compute_array(); - compute->invoked_flag |= Compute::INVOKED_ARRAY; + if (!(val.val.c->invoked_flag & Compute::INVOKED_ARRAY)) { + val.val.c->compute_array(); + val.val.c->invoked_flag |= Compute::INVOKED_ARRAY; } - if (compute->array) - bin_vector_weights(compute->size_array_rows,&compute->array[0][j-1], - compute->size_array_cols,weights,stride); + if (val.val.c->array) + bin_vector_weights(val.val.c->size_array_rows,&val.val.c->array[0][j-1], + val.val.c->size_array_cols,weights,stride); } } else if (kind == PERATOM) { - if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { - compute->compute_peratom(); - compute->invoked_flag |= Compute::INVOKED_PERATOM; + if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) { + val.val.c->compute_peratom(); + val.val.c->invoked_flag |= Compute::INVOKED_PERATOM; } if (j == 0) - bin_atoms_weights(compute->vector_atom,1,weights, stride); - else if (compute->array_atom) - bin_atoms_weights(&compute->array_atom[0][j-1], - compute->size_peratom_cols,weights,stride); + bin_atoms_weights(val.val.c->vector_atom,1,weights, stride); + else if (val.val.c->array_atom) + bin_atoms_weights(&val.val.c->array_atom[0][j-1], + val.val.c->size_peratom_cols,weights,stride); } else if (kind == LOCAL) { - if (!(compute->invoked_flag & Compute::INVOKED_LOCAL)) { - compute->compute_local(); - compute->invoked_flag |= Compute::INVOKED_LOCAL; + if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) { + val.val.c->compute_local(); + val.val.c->invoked_flag |= Compute::INVOKED_LOCAL; } if (j == 0) - bin_vector_weights(compute->size_local_rows, - compute->vector_local,1,weights,stride); - else if (compute->array_local) - bin_vector_weights(compute->size_local_rows, - &compute->array_local[0][j-1], - compute->size_local_cols,weights,stride); + bin_vector_weights(val.val.c->size_local_rows, + val.val.c->vector_local,1,weights,stride); + else if (val.val.c->array_local) + bin_vector_weights(val.val.c->size_local_rows, + &val.val.c->array_local[0][j-1], + val.val.c->size_local_cols,weights,stride); } // access fix fields, guaranteed to be ready - } else if (which[i] == ArgInfo::FIX) { - - Fix *fix = modify->fix[m]; + } else if (val.which == ArgInfo::FIX) { if (kind == GLOBAL && mode == SCALAR) { - if (j == 0) bin_one_weights(fix->compute_scalar(),weight); - else bin_one_weights(fix->compute_vector(j-1),weight); + if (j == 0) bin_one_weights(val.val.f->compute_scalar(),weight); + else bin_one_weights(val.val.f->compute_vector(j-1),weight); } else if (kind == GLOBAL && mode == VECTOR) { if (j == 0) { - int n = fix->size_vector; - for (i = 0; i < n; i++) - bin_one_weights(fix->compute_vector(i),weights[i*stride]); + int n = val.val.f->size_vector; + for (int i = 0; i < n; i++) + bin_one_weights(val.val.f->compute_vector(i),weights[i*stride]); } else { - int n = fix->size_vector; - for (i = 0; i < n; i++) - bin_one_weights(fix->compute_array(i,j-1),weights[i*stride]); + int n = val.val.f->size_vector; + for (int i = 0; i < n; i++) + bin_one_weights(val.val.f->compute_array(i,j-1),weights[i*stride]); } } else if (kind == PERATOM) { if (j == 0) - bin_atoms_weights(fix->vector_atom,1,weights,stride); - else if (fix->array_atom) - bin_atoms_weights(fix->array_atom[j-1],fix->size_peratom_cols, + bin_atoms_weights(val.val.f->vector_atom,1,weights,stride); + else if (val.val.f->array_atom) + bin_atoms_weights(val.val.f->array_atom[j-1],val.val.f->size_peratom_cols, weights,stride); } else if (kind == LOCAL) { - if (j == 0) bin_vector_weights(fix->size_local_rows,fix->vector_local,1, + if (j == 0) bin_vector_weights(val.val.f->size_local_rows,val.val.f->vector_local,1, weights,stride); - else if (fix->array_local) - bin_vector_weights(fix->size_local_rows,&fix->array_local[0][j-1], - fix->size_local_cols,weights,stride); + else if (val.val.f->array_local) + bin_vector_weights(val.val.f->size_local_rows,&val.val.f->array_local[0][j-1], + val.val.f->size_local_cols,weights,stride); } // evaluate equal-style variable - } else if (which[i] == ArgInfo::VARIABLE && kind == GLOBAL) { - bin_one_weights(input->variable->compute_equal(m),weight); + } else if (val.which == ArgInfo::VARIABLE && kind == GLOBAL) { + bin_one_weights(input->variable->compute_equal(val.val.v),weight); - } else if (which[i] == ArgInfo::VARIABLE && kind == PERATOM) { + } else if (val.which == ArgInfo::VARIABLE && kind == PERATOM) { if (atom->nmax > maxatom) { memory->destroy(vector); maxatom = atom->nmax; memory->create(vector,maxatom,"ave/histo/weight:vector"); } - input->variable->compute_atom(m,igroup,vector,1,0); + input->variable->compute_atom(val.val.v,igroup,vector,1,0); bin_atoms_weights(vector,1,weights,stride); } @@ -409,7 +394,7 @@ void FixAveHistoWeight::end_of_step() stats[1] = stats_all[1]; stats[2] = stats_all[2]; stats[3] = stats_all[3]; - for (i = 0; i < nbins; i++) bin[i] = bin_all[i]; + for (int i = 0; i < nbins; i++) bin[i] = bin_all[i]; } // if ave = ONE, only single Nfreq timestep value is needed @@ -421,14 +406,14 @@ void FixAveHistoWeight::end_of_step() stats_total[1] = stats[1]; stats_total[2] = stats[2]; stats_total[3] = stats[3]; - for (i = 0; i < nbins; i++) bin_total[i] = bin[i]; + for (int i = 0; i < nbins; i++) bin_total[i] = bin[i]; } else if (ave == RUNNING) { stats_total[0] += stats[0]; stats_total[1] += stats[1]; stats_total[2] = MIN(stats_total[2],stats[2]); stats_total[3] = MAX(stats_total[3],stats[3]); - for (i = 0; i < nbins; i++) bin_total[i] += bin[i]; + for (int i = 0; i < nbins; i++) bin_total[i] += bin[i]; } else if (ave == WINDOW) { stats_total[0] += stats[0]; @@ -438,19 +423,20 @@ void FixAveHistoWeight::end_of_step() if (window_limit) stats_total[1] -= stats_list[iwindow][1]; stats_list[iwindow][1] = stats[1]; + int m; if (window_limit) m = nwindow; else m = iwindow+1; stats_list[iwindow][2] = stats[2]; stats_total[2] = stats_list[0][2]; - for (i = 1; i < m; i++) + for (int i = 1; i < m; i++) stats_total[2] = MIN(stats_total[2],stats_list[i][2]); stats_list[iwindow][3] = stats[3]; stats_total[3] = stats_list[0][3]; - for (i = 1; i < m; i++) + for (int i = 1; i < m; i++) stats_total[3] = MAX(stats_total[3],stats_list[i][3]); - for (i = 0; i < nbins; i++) { + for (int i = 0; i < nbins; i++) { bin_total[i] += bin[i]; if (window_limit) bin_total[i] -= bin_list[iwindow][i]; bin_list[iwindow][i] = bin[i]; @@ -465,17 +451,17 @@ void FixAveHistoWeight::end_of_step() // output result to file - if (fp && me == 0) { + if (fp && comm->me == 0) { clearerr(fp); if (overwrite) platform::fseek(fp,filepos); fmt::print(fp,"{} {} {} {} {} {}\n",ntimestep,nbins, stats_total[0],stats_total[1],stats_total[2],stats_total[3]); if (stats_total[0] != 0.0) - for (i = 0; i < nbins; i++) + for (int i = 0; i < nbins; i++) fprintf(fp,"%d %g %g %g\n", i+1,coord[i],bin_total[i],bin_total[i]/stats_total[0]); else - for (i = 0; i < nbins; i++) + for (int i = 0; i < nbins; i++) fprintf(fp,"%d %g %g %g\n",i+1,coord[i],0.0,0.0); if (ferror(fp))