ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

convert fix ave/histo and fix ave/histo/weight

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-10-14 17:49:56 -04:00
parent 7d9076de4d
commit b3557b81c1
3 changed files with 442 additions and 546 deletions
Download Patch File Download Diff File Expand all files Collapse all files

644
src/fix_ave_histo.cpp
View File

@ -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]);
}
}

18
src/fix_ave_histo.h
View File

@ -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<value_t> 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;

326
src/fix_ave_histo_weight.cpp
View File

@ -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))