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

Merge pull request #3475 from akohlmey/refactor-testing

Browse Source 
Refactor computes and fixes to use array of structs instead of many index arrays
This commit is contained in:
Axel Kohlmeyer
2022-10-27 17:07:52 -04:00
committed by GitHub
parent 850a303f34 d96bfcdf06
commit 1de74508c3
25 changed files with 2312 additions and 2481 deletions
Download Patch File Download Diff File Expand all files Collapse all files

537
src/EXTRA-FIX/fix_ave_correlate_long.cpp
View File

@ -26,9 +26,11 @@
#include "arg_info.h"
#include "citeme.h"
#include "comm.h"
#include "compute.h"
#include "error.h"
#include "input.h"
#include "math_special.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
@ -39,63 +41,72 @@
using namespace LAMMPS_NS;
using namespace FixConst;
using MathSpecial::powint;
enum{AUTO,UPPER,LOWER,AUTOUPPER,AUTOLOWER,FULL};
enum { AUTO, UPPER, LOWER, AUTOUPPER, AUTOLOWER, FULL };
static const char cite_fix_ave_correlate_long[] =
"fix ave/correlate/long command: doi:10.1063/1.3491098\n\n"
"@Article{Ramirez10,\n"
" author = {Jorge Rami{\'}rez and Sathish K. Sukumaran and Bart Vorselaars and Alexei E. Likhtman},\n"
" title = {Efficient on the Fly Calculation of Time Correlation Functions in Computer Simulations},"
" journal = {J.~Chem.\\ Phys.},\n"
" year = 2010,\n"
" volume = 133,\n"
" number = 15,\n"
" pages = {154103}\n"
"}\n\n";
"fix ave/correlate/long command: doi:10.1063/1.3491098\n\n"
"@Article{Ramirez10,\n"
" author = {Jorge Rami{\'}rez and Sathish K. Sukumaran and Bart Vorselaars and Alexei E. "
"Likhtman},\n"
" title = {Efficient on the Fly Calculation of Time Correlation Functions in Computer "
"Simulations},"
" journal = {J.~Chem.\\ Phys.},\n"
" year = 2010,\n"
" volume = 133,\n"
" number = 15,\n"
" pages = {154103}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS * lmp, int narg, char **arg):
Fix (lmp, narg, arg)
FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), shift(nullptr), shift2(nullptr), correlation(nullptr),
accumulator(nullptr), accumulator2(nullptr), ncorrelation(nullptr), naccumulator(nullptr),
insertindex(nullptr), fp(nullptr), cvalues(nullptr)
{
if (lmp->citeme) lmp->citeme->add(cite_fix_ave_correlate_long);
// At least nevery nfrez and one value are needed
if (narg < 6) error->all(FLERR,"Illegal fix ave/correlate/long command");
// At least nevery nfreq and one value are needed
if (narg < 6) utils::missing_cmd_args(FLERR, "fix ave/correlate/long", error);
MPI_Comm_rank(world,&me);
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
nfreq = utils::inumeric(FLERR,arg[4],false,lmp);
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
nfreq = utils::inumeric(FLERR, arg[4], false, lmp);
restart_global = 1;
global_freq = nfreq;
time_depend = 1;
// parse values until one isn't recognized
// expand args if any have wildcard character "*"
which = new int[narg-5];
argindex = new int[narg-5];
ids = new char*[narg-5];
value2index = new int[narg-5];
nvalues = 0;
int expand = 0;
char **earg;
int nargnew = utils::expand_args(FLERR, narg - 5, &arg[5], 0, earg, lmp);
int iarg = 5;
while (iarg < narg) {
if (earg != &arg[5]) expand = 1;
arg = earg;
// parse values
int iarg = 0;
while (iarg < nargnew) {
ArgInfo argi(arg[iarg]);
value_t val;
if (argi.get_type() == ArgInfo::NONE) break;
if ((argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_dim() > 1))
error->all(FLERR,"Illegal fix ave/correlate/long command");
error->all(FLERR, "Unknown fix ave/correlate/long data type: {}", arg[iarg]);
which[nvalues] = argi.get_type();
argindex[nvalues] = argi.get_index1();
ids[nvalues] = argi.copy_name();
val.which = argi.get_type();
val.argindex = argi.get_index1();
val.id = argi.get_name();
val.val.c = nullptr;
nvalues++;
values.push_back(val);
iarg++;
}
nvalues = values.size();
// optional args
@ -103,128 +114,136 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS * lmp, int narg, char **arg):
startstep = 0;
fp = nullptr;
overwrite = 0;
numcorrelators=20;
numcorrelators = 20;
p = 16;
m = 2;
char *title1 = nullptr;
char *title2 = nullptr;
while (iarg < narg) {
if (strcmp(arg[iarg],"type") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/correlate/long command");
if (strcmp(arg[iarg+1],"auto") == 0) type = AUTO;
else if (strcmp(arg[iarg+1],"upper") == 0) type = UPPER;
else if (strcmp(arg[iarg+1],"lower") == 0) type = LOWER;
else if (strcmp(arg[iarg+1],"auto/upper") == 0) type = AUTOUPPER;
else if (strcmp(arg[iarg+1],"auto/lower") == 0) type = AUTOLOWER;
else if (strcmp(arg[iarg+1],"full") == 0) type = FULL;
else error->all(FLERR,"Illegal fix ave/correlate/long command");
while (iarg < nargnew) {
if (strcmp(arg[iarg], "type") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long type", error);
if (strcmp(arg[iarg + 1], "auto") == 0)
type = AUTO;
else if (strcmp(arg[iarg + 1], "upper") == 0)
type = UPPER;
else if (strcmp(arg[iarg + 1], "lower") == 0)
type = LOWER;
else if (strcmp(arg[iarg + 1], "auto/upper") == 0)
type = AUTOUPPER;
else if (strcmp(arg[iarg + 1], "auto/lower") == 0)
type = AUTOLOWER;
else if (strcmp(arg[iarg + 1], "full") == 0)
type = FULL;
else
error->all(FLERR, "Unknown fix ave/correlate/long type: {}");
iarg += 2;
} else if (strcmp(arg[iarg],"start") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/correlate/long command");
startstep = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "start") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long start", error);
startstep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"ncorr") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/correlate/long command");
numcorrelators = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "ncorr") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long ncorr", error);
numcorrelators = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"nlen") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/correlate/long command");
p = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "nlen") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long nlen", error);
p = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"ncount") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/correlate/long command");
m = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
} else if (strcmp(arg[iarg], "ncount") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long ncount", error);
m = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"file") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/correlate/long command");
if (me == 0) {
fp = fopen(arg[iarg+1],"w");
} else if (strcmp(arg[iarg], "file") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long file", error);
if (comm->me == 0) {
fp = fopen(arg[iarg + 1], "w");
if (fp == nullptr)
error->one(FLERR,"Cannot open fix ave/correlate/long file {}: {}",
arg[iarg+1],utils::getsyserror());
error->one(FLERR, "Cannot open fix ave/correlate/long file {}: {}", arg[iarg + 1],
utils::getsyserror());
}
iarg += 2;
} else if (strcmp(arg[iarg],"overwrite") == 0) {
} 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/correlate/long command");
} else if (strcmp(arg[iarg], "title1") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long title1", error);
delete[] title1;
title1 = utils::strdup(arg[iarg+1]);
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/correlate/long command");
} else if (strcmp(arg[iarg], "title2") == 0) {
if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate/long title2", error);
delete[] title2;
title2 = utils::strdup(arg[iarg+1]);
title2 = utils::strdup(arg[iarg + 1]);
iarg += 2;
} else error->all(FLERR,"Illegal fix ave/correlate/long command");
} else
error->all(FLERR, "Unknown fix ave/correlate/long keyword: {}", arg[iarg]);
}
if (p % m != 0) error->all(FLERR,"fix_correlator: p mod m must be 0");
dmin = p/m;
length = numcorrelators*p;
if (p % m != 0) error->all(FLERR, "Fix ave/correlate/long: nlen must be divisible by ncount");
dmin = p / m;
length = numcorrelators * p;
npcorr = 0;
kmax = 0;
// setup and error check
// for fix inputs, check that fix frequency is acceptable
if (nevery <= 0 || nfreq <= 0)
error->all(FLERR,"Illegal fix ave/correlate/long command");
if (nfreq % nevery)
error->all(FLERR,"Illegal fix ave/correlate/long command");
if (nevery <= 0) error->all(FLERR, "Illegal fix ave/correlate/long nevery value: {}", nevery);
if (nfreq <= 0) error->all(FLERR, "Illegal fix ave/correlate/long nfreq value: {}", nfreq);
if (nfreq % nevery) error->all(FLERR, "Inconsistent fix ave/correlate/long nevery/nfreq values");
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/correlate/long does not exist");
if (argindex[i] == 0 && modify->compute[icompute]->scalar_flag == 0)
error->all(FLERR,"Fix ave/correlate/long compute does not calculate a scalar");
if (argindex[i] && modify->compute[icompute]->vector_flag == 0)
error->all(FLERR,"Fix ave/correlate/long compute does not calculate a vector");
if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector)
error->all(FLERR,"Fix ave/correlate/long compute vector is accessed out-of-range");
for (auto &val : values) {
} else if (which[i] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/correlate/long does not exist");
if (argindex[i] == 0 && modify->fix[ifix]->scalar_flag == 0)
error->all(FLERR,"Fix ave/correlate/long fix does not calculate a scalar");
if (argindex[i] && modify->fix[ifix]->vector_flag == 0)
error->all(FLERR,"Fix ave/correlate/long fix does not calculate a vector");
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector)
error->all(FLERR,"Fix ave/correlate/long fix vector is accessed out-of-range");
if (nevery % modify->fix[ifix]->global_freq)
error->all(FLERR,"Fix for fix ave/correlate/long not computed at compatible time");
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 fix ave/correlate/long does not exist", val.id);
if (val.argindex == 0 && val.val.c->scalar_flag == 0)
error->all(FLERR, "Fix ave/correlate/long compute {} does not calculate a scalar", val.id);
if (val.argindex && val.val.c->vector_flag == 0)
error->all(FLERR, "Fix ave/correlate/long compute {} does not calculate a vector", val.id);
if (val.argindex && val.argindex > val.val.c->size_vector)
error->all(FLERR, "Fix ave/correlate/long compute {} vector is accessed out-of-range",
val.id);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/correlate/long does not exist");
if (input->variable->equalstyle(ivariable) == 0)
error->all(FLERR,"Fix ave/correlate/long variable is not equal-style variable");
} 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 fix ave/correlate/long does not exist", val.id);
if (val.argindex == 0 && val.val.f->scalar_flag == 0)
error->all(FLERR, "Fix ave/correlate/long fix {} does not calculate a scalar", val.id);
if (val.argindex && val.val.f->vector_flag == 0)
error->all(FLERR, "Fix ave/correlate/long fix {} does not calculate a vector", val.id);
if (val.argindex && val.argindex > val.val.f->size_vector)
error->all(FLERR, "Fix ave/correlate/long fix {} vector is accessed out-of-range", val.id);
if (nevery % val.val.f->global_freq)
error->all(FLERR, "Fix {} for fix ave/correlate/long not computed at compatible time",
val.id);
} 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 fix ave/correlate/long does not exist", val.id);
if (val.argindex == 0 && input->variable->equalstyle(val.val.v) == 0)
error->all(FLERR, "Fix ave/correlate/long variable {} is not equal-style variable", val.id);
if (val.argindex && input->variable->vectorstyle(val.val.v) == 0)
error->all(FLERR, "Fix ave/correlate/long variable {} is not vector-style variable",
val.id);
}
}
// npair = # of correlation pairs to calculate
if (type == AUTO) npair = nvalues;
if (type == UPPER || type == LOWER) npair = nvalues*(nvalues-1)/2;
if (type == AUTOUPPER || type == AUTOLOWER) npair = nvalues*(nvalues+1)/2;
if (type == FULL) npair = nvalues*nvalues;
if (type == UPPER || type == LOWER) npair = nvalues * (nvalues - 1) / 2;
if (type == AUTOUPPER || type == AUTOLOWER) npair = nvalues * (nvalues + 1) / 2;
if (type == FULL) npair = nvalues * nvalues;
// 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,"# Time-correlated data for fix %s\n",id);
if (title2) fprintf(fp,"%s\n",title2);
@ -232,38 +251,49 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS * lmp, int narg, char **arg):
fprintf(fp,"# Time");
if (type == AUTO)
for (int i = 0; i < nvalues; i++)
fprintf(fp," %s*%s",arg[5+i],arg[5+i]);
fprintf(fp," %s*%s",earg[i],earg[i]);
else if (type == UPPER)
for (int i = 0; i < nvalues; i++)
for (int j = i+1; j < nvalues; j++)
fprintf(fp," %s*%s",arg[5+i],arg[5+j]);
fprintf(fp," %s*%s",earg[i],earg[j]);
else if (type == LOWER)
for (int i = 0; i < nvalues; i++)
for (int j = 0; j < i-1; j++)
fprintf(fp," %s*%s",arg[5+i],arg[5+j]);
fprintf(fp," %s*%s",earg[i],earg[j]);
else if (type == AUTOUPPER)
for (int i = 0; i < nvalues; i++)
for (int j = i; j < nvalues; j++)
fprintf(fp," %s*%s",arg[5+i],arg[5+j]);
fprintf(fp," %s*%s",earg[i],earg[j]);
else if (type == AUTOLOWER)
for (int i = 0; i < nvalues; i++)
for (int j = 0; j < i; j++)
fprintf(fp," %s*%s",arg[5+i],arg[5+j]);
fprintf(fp," %s*%s",earg[i],earg[j]);
else if (type == FULL)
for (int i = 0; i < nvalues; i++)
for (int j = 0; j < nvalues; j++)
fprintf(fp," %s*%s",arg[5+i],arg[5+j]);
fprintf(fp," %s*%s",earg[i],earg[j]);
fprintf(fp,"\n");
}
if (ferror(fp))
error->one(FLERR,"Error writing ave/correlate/long header: {}", utils::getsyserror());
filepos = platform::ftell(fp);
}
delete[] title1;
delete[] title2;
// if wildcard expansion occurred, free earg memory from expand_args()
// wait to do this until after file comment lines are printed
if (expand) {
for (int i = 0; i < nargnew; i++) delete[] earg[i];
memory->sfree(earg);
}
// allocate and initialize memory for calculated values and correlators
memory->create(values,nvalues,"correlator:values");
memory->create(cvalues,nvalues,"correlator:values");
memory->create(shift,npair,numcorrelators,p,"correlator:shift");
memory->create(shift2,npair,numcorrelators,p,"correlator:shift2"); //NOT OPTMAL
memory->create(correlation,npair,numcorrelators,p,"correlator:correlation");
@ -276,10 +306,10 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS * lmp, int narg, char **arg):
memory->create(t,length,"correlator:t");
memory->create(f,npair,length,"correlator:f");
for (int i=0;i<npair;i++)
for (int j=0;j<numcorrelators;j++) {
for (unsigned int k=0;k<p;k++) {
shift[i][j][k]=-2E10;
for (int i=0; i < npair; i++)
for (int j=0; j < numcorrelators; j++) {
for (unsigned int k=0; k < p; k++) {
shift[i][j][k]=-2e10;
shift2[i][j][k]=0.0;
correlation[i][j][k]=0.0;
}
@ -287,16 +317,15 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS * lmp, int narg, char **arg):
accumulator2[i][j]=0.0;
}
for (int i=0;i<numcorrelators;i++) {
for (unsigned int j=0;j<p;j++) ncorrelation[i][j]=0;
for (int i=0; i < numcorrelators; i++) {
for (unsigned int j=0; j < p; j++) ncorrelation[i][j]=0;
naccumulator[i]=0;
insertindex[i]=0;
}
for (int i=0;i<length;i++) t[i]=0.0;
for (int i=0;i<npair;i++)
for (int j=0;j<length;j++) f[i][j]=0.0;
for (int i=0; i < length; i++) t[i]=0.0;
for (int i=0; i < npair; i++)
for (int j=0; j < length; j++) f[i][j]=0.0;
// nvalid = next step on which end_of_step does something
// add nvalid to all computes that store invocation times
@ -314,13 +343,7 @@ FixAveCorrelateLong::FixAveCorrelateLong(LAMMPS * lmp, int narg, char **arg):
FixAveCorrelateLong::~FixAveCorrelateLong()
{
delete[] which;
delete[] argindex;
delete[] value2index;
for (int i = 0; i < nvalues; i++) delete[] ids[i];
delete[] ids;
memory->destroy(values);
memory->destroy(cvalues);
memory->destroy(shift);
memory->destroy(shift2);
memory->destroy(correlation);
@ -332,7 +355,7 @@ FixAveCorrelateLong::~FixAveCorrelateLong()
memory->destroy(t);
memory->destroy(f);
if (fp && me == 0) fclose(fp);
if (fp && comm->me == 0) fclose(fp);
}
/* ---------------------------------------------------------------------- */
@ -350,24 +373,22 @@ void FixAveCorrelateLong::init()
{
// set current indices for all computes,fixes,variables
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/correlate/long does not exist");
value2index[i] = icompute;
for (auto &val : values) {
} else if (which[i] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/correlate/long does not exist");
value2index[i] = ifix;
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 fix ave/correlate/long does not exist", val.id);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/correlate/long does not exist");
value2index[i] = ivariable;
} 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 fix ave/correlate/long does not exist", val.id);
} 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 fix ave/correlate/long does not exist", val.id);
}
}
@ -392,9 +413,6 @@ void FixAveCorrelateLong::setup(int /*vflag*/)
void FixAveCorrelateLong::end_of_step()
{
int i,m;
double scalar;
// skip if not step which requires doing something
bigint ntimestep = update->ntimestep;
@ -406,43 +424,51 @@ void FixAveCorrelateLong::end_of_step()
modify->clearstep_compute();
for (i = 0; i < nvalues; i++) {
m = value2index[i];
scalar = 0.0;
int i = 0;
for (auto &val : values) {
double scalar = 0.0;
// invoke compute if not previously invoked
if (which[i] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[m];
if (val.which == ArgInfo::COMPUTE) {
if (argindex[i] == 0) {
if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) {
compute->compute_scalar();
compute->invoked_flag |= Compute::INVOKED_SCALAR;
if (val.argindex == 0) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_SCALAR)) {
val.val.c->compute_scalar();
val.val.c->invoked_flag |= Compute::INVOKED_SCALAR;
}
scalar = compute->scalar;
scalar = 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;
}
scalar = compute->vector[argindex[i]-1];
scalar = val.val.c->vector[val.argindex-1];
}
// access fix fields, guaranteed to be ready
// access fix fields, guaranteed to be ready
} else if (which[i] == ArgInfo::FIX) {
if (argindex[i] == 0)
scalar = modify->fix[m]->compute_scalar();
} else if (val.which == ArgInfo::FIX) {
if (val.argindex == 0)
scalar = val.val.f->compute_scalar();
else
scalar = modify->fix[m]->compute_vector(argindex[i]-1);
scalar = val.val.f->compute_vector(val.argindex-1);
// evaluate equal-style variable
// evaluate equal-style or vector-style variable
} else if (which[i] == ArgInfo::VARIABLE)
scalar = input->variable->compute_equal(m);
} else if (val.which == ArgInfo::VARIABLE) {
if (val.argindex == 0)
scalar = input->variable->compute_equal(val.val.v);
else {
double *varvec;
int nvec = input->variable->compute_vector(val.val.v,&varvec);
int index = val.argindex;
if (nvec < index) scalar = 0.0;
else scalar = varvec[index-1];
}
}
values[i] = scalar;
cvalues[i++] = scalar;
}
// fistindex = index in values ring of earliest time sample
@ -456,19 +482,25 @@ void FixAveCorrelateLong::end_of_step()
if (ntimestep % nfreq) return;
// output result to file
evaluate();
if (fp && me == 0) {
if (fp && comm->me == 0) {
clearerr(fp);
if (overwrite) platform::fseek(fp,filepos);
fmt::print(fp,"# Timestep: {}\n", ntimestep);
for (unsigned int i=0;i<npcorr;++i) {
for (unsigned int i=0; i < npcorr; ++i) {
fprintf(fp, "%lg ", t[i]*update->dt*nevery);
for (int j=0;j<npair;++j) {
for (int j=0; j < npair; ++j) {
fprintf(fp, "%lg ", f[j][i]);
}
fprintf(fp, "\n");
}
if (ferror(fp))
error->one(FLERR,"Error writing out fix ave/correlate/long data: {}", utils::getsyserror());
fflush(fp);
if (overwrite) {
bigint fileend = platform::ftell(fp);
if ((fileend > 0) && (platform::ftruncate(fp,fileend)))
@ -481,20 +513,20 @@ void FixAveCorrelateLong::evaluate() {
unsigned int jm=0;
// First correlator
for (unsigned int j=0;j<p;++j) {
for (unsigned int j=0; j < p; ++j) {
if (ncorrelation[0][j] > 0) {
t[jm] = j;
for (int i=0;i<npair;++i)
for (int i=0;i < npair; ++i)
f[i][jm] = correlation[i][0][j]/ncorrelation[0][j];
++jm;
}
}
// Subsequent correlators
for (int k=1;k<kmax;++k) {
for (unsigned int j=dmin;j<p;++j) {
for (int k=1; k < kmax; ++k) {
for (unsigned int j=dmin; j < p; ++j) {
if (ncorrelation[k][j]>0) {
t[jm] = j * pow((double)m, k);
t[jm] = j * powint((double)m, k);
for (int i=0;i<npair;++i)
f[i][jm] = correlation[i][k][j] / ncorrelation[k][j];
++jm;
@ -517,35 +549,35 @@ void FixAveCorrelateLong::accumulate()
if (update->ntimestep <= last_accumulated_step) return;
if (type == AUTO) {
for (i=0; i<nvalues;i++) add(i,values[i]);
for (i=0; i < nvalues; i++) add(i,cvalues[i]);
} else if (type == UPPER) {
ipair = 0;
for (i=0;i<nvalues;i++)
for (j=i+1;j<nvalues;j++) add(ipair++,values[i],values[j]);
for (i=0; i < nvalues; i++)
for (j=i+1; j < nvalues; j++) add(ipair++,cvalues[i],cvalues[j]);
} else if (type == LOWER) {
ipair = 0;
for (i=0;i<nvalues;i++)
for (j=0;j<i;j++) add(ipair++,values[i],values[j]);
for (i=0; i < nvalues; i++)
for (j=0; j < i; j++) add(ipair++,cvalues[i],cvalues[j]);
} else if (type == AUTOUPPER) {
ipair = 0;
for (i=0;i<nvalues;i++)
for (j=i;j<nvalues;j++) {
if (i==j) add(ipair++,values[i]);
else add(ipair++,values[i],values[j]);
for (i=0; i < nvalues; i++)
for (j=i; j < nvalues; j++) {
if (i == j) add(ipair++,cvalues[i]);
else add(ipair++,cvalues[i],cvalues[j]);
}
} else if (type == AUTOLOWER) {
ipair = 0;
for (i=0;i<nvalues;i++)
for (j=0;j<=i;j++) {
if (i==j) add(ipair++,values[i]);
else add(ipair++,values[i],values[j]);
for (i=0; i < nvalues; i++)
for (j=0; j <=i; j++) {
if (i==j) add(ipair++,cvalues[i]);
else add(ipair++,cvalues[i],cvalues[j]);
}
} else if (type == FULL) {
ipair = 0;
for (i=0;i<nvalues;i++)
for (j=0;j<nvalues;j++) {
if (i==j) add(ipair++,values[i]);
else add(ipair++,values[i],values[j]);
for (i=0; i < nvalues; i++)
for (j=0; j < nvalues; j++) {
if (i == j) add(ipair++,cvalues[i]);
else add(ipair++,cvalues[i],cvalues[j]);
}
}
last_accumulated_step = update->ntimestep;
@ -565,38 +597,38 @@ void FixAveCorrelateLong::add(const int i, const double w, const int k) {
// Add to accumulator and, if needed, add to next correlator
accumulator[i][k] += w;
if (i==0) ++naccumulator[k];
if (i == 0) ++naccumulator[k];
if (naccumulator[k]==m) {
add(i,accumulator[i][k]/m, k+1);
accumulator[i][k]=0;
if (i==npair-1) naccumulator[k]=0;
if (i == npair-1) naccumulator[k]=0;
}
// Calculate correlation function
unsigned int ind1=insertindex[k];
if (k==0) { // First correlator is different
if (k == 0) { // First correlator is different
int ind2=ind1;
for (unsigned int j=0;j<p;++j) {
for (unsigned int j=0; j < p; ++j) {
if (shift[i][k][ind2] > -1e10) {
correlation[i][k][j]+= shift[i][k][ind1]*shift[i][k][ind2];
if (i==0) ++ncorrelation[k][j];
}
--ind2;
if (ind2<0) ind2+=p;
if (ind2 < 0) ind2+=p;
}
} else {
int ind2=ind1-dmin;
for (unsigned int j=dmin;j<p;++j) {
if (ind2<0) ind2+=p;
if (ind2 < 0) ind2 += p;
if (shift[i][k][ind2] > -1e10) {
correlation[i][k][j]+= shift[i][k][ind1]*shift[i][k][ind2];
if (i==0) ++ncorrelation[k][j];
if (i == 0) ++ncorrelation[k][j];
}
--ind2;
}
}
if (i==npair-1) {
if (i == npair-1) {
++insertindex[k];
if (insertindex[k]==p) insertindex[k]=0;
}
@ -606,8 +638,7 @@ void FixAveCorrelateLong::add(const int i, const double w, const int k) {
/* ----------------------------------------------------------------------
Add 2 scalar values to the cross-correlator k of pair i
------------------------------------------------------------------------- */
void FixAveCorrelateLong::add(const int i, const double wA, const double wB,
const int k) {
void FixAveCorrelateLong::add(const int i, const double wA, const double wB, const int k) {
if (k == numcorrelators) return;
if (k > kmax) kmax=k;
@ -616,21 +647,21 @@ void FixAveCorrelateLong::add(const int i, const double wA, const double wB,
accumulator[i][k] += wA;
accumulator2[i][k] += wB;
if (i==0) ++naccumulator[k];
if (naccumulator[k]==m) {
if (i == 0) ++naccumulator[k];
if (naccumulator[k] == m) {
add(i,accumulator[i][k]/m, accumulator2[i][k]/m,k+1);
accumulator[i][k]=0;
accumulator2[i][k]=0;
if (i==npair-1) naccumulator[k]=0;
if (i == npair-1) naccumulator[k]=0;
}
unsigned int ind1=insertindex[k];
if (k==0) {
if (k == 0) {
int ind2=ind1;
for (unsigned int j=0;j<p;++j) {
for (unsigned int j=0; j < p; ++j) {
if (shift[i][k][ind2] > -1e10) {
correlation[i][k][j]+= shift[i][k][ind1]*shift2[i][k][ind2];
if (i==0) ++ncorrelation[k][j];
if (i == 0) ++ncorrelation[k][j];
}
--ind2;
if (ind2<0) ind2+=p;
@ -638,19 +669,19 @@ void FixAveCorrelateLong::add(const int i, const double wA, const double wB,
}
else {
int ind2=ind1-dmin;
for (unsigned int j=dmin;j<p;++j) {
if (ind2<0) ind2+=p;
for (unsigned int j=dmin; j < p; ++j) {
if (ind2 < 0) ind2+=p;
if (shift[i][k][ind2] > -1e10) {
correlation[i][k][j]+= shift[i][k][ind1]*shift2[i][k][ind2];
if (i==0) ++ncorrelation[k][j];
if (i == 0) ++ncorrelation[k][j];
}
--ind2;
}
}
if (i==npair-1) {
if (i == npair-1) {
++insertindex[k];
if (insertindex[k]==p) insertindex[k]=0;
if (insertindex[k] == p) insertindex[k]=0;
}
}
@ -696,20 +727,20 @@ double FixAveCorrelateLong::memory_usage() {
------------------------------------------------------------------------- */
// Save everything except t and f
void FixAveCorrelateLong::write_restart(FILE *fp) {
if (me == 0) {
if (comm->me == 0) {
int nsize = 3*npair*numcorrelators*p + 2*npair*numcorrelators
+ numcorrelators*p + 2*numcorrelators + 6;
int n=0;
double *list;
memory->create(list,nsize,"correlator:list");
list[n++]=npair;
list[n++]=numcorrelators;
list[n++]=p;
list[n++]=m;
list[n++] = npair;
list[n++] = numcorrelators;
list[n++] = p;
list[n++] = m;
list[n++] = last_accumulated_step;
for (int i=0;i<npair;i++)
for (int j=0;j<numcorrelators;j++) {
for (unsigned int k=0;k<p;k++) {
for (int i=0; i < npair; i++)
for (int j=0; j < numcorrelators; j++) {
for (unsigned int k=0; k < p; k++) {
list[n++]=shift[i][j][k];
list[n++]=shift2[i][j][k];
list[n++]=correlation[i][j][k];
@ -717,8 +748,8 @@ void FixAveCorrelateLong::write_restart(FILE *fp) {
list[n++]=accumulator[i][j];
list[n++]=accumulator2[i][j];
}
for (int i=0;i<numcorrelators;i++) {
for (unsigned int j=0;j<p;j++) list[n++]=ncorrelation[i][j];
for (int i=0; i<numcorrelators; i++) {
for (unsigned int j=0; j < p; j++) list[n++]=ncorrelation[i][j];
list[n++]=naccumulator[i];
list[n++]=insertindex[i];
}
@ -730,7 +761,6 @@ void FixAveCorrelateLong::write_restart(FILE *fp) {
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
@ -738,18 +768,17 @@ void FixAveCorrelateLong::restart(char *buf)
{
int n = 0;
auto list = (double *) buf;
int npairin = static_cast<int> (list[n++]);
int npairin = static_cast<int>(list[n++]);
int numcorrelatorsin = static_cast<int> (list[n++]);
int pin = static_cast<int> (list[n++]);
int min = static_cast<int> (list[n++]);
last_accumulated_step = static_cast<int> (list[n++]);
int pin = static_cast<int>(list[n++]);
int min = static_cast<int>(list[n++]);
last_accumulated_step = static_cast<int>(list[n++]);
if ((npairin!=npair) || (numcorrelatorsin!=numcorrelators)
|| (pin!=(int)p) || (min!=(int)m))
error->all(FLERR,"Fix ave/correlate/long: restart and input data are different");
if ((npairin!=npair) || (numcorrelatorsin!=numcorrelators) || (pin!=(int)p) || (min!=(int)m))
error->all(FLERR, "Fix ave/correlate/long: restart and input data are different");
for (int i=0;i<npair;i++)
for (int j=0;j<numcorrelators;j++) {
for (int i=0; i < npair; i++)
for (int j=0; j < numcorrelators; j++) {
for (unsigned int k=0;k<p;k++) {
shift[i][j][k] = list[n++];
shift2[i][j][k] = list[n++];
@ -758,10 +787,10 @@ void FixAveCorrelateLong::restart(char *buf)
accumulator[i][j] = list[n++];
accumulator2[i][j] = list[n++];
}
for (int i=0;i<numcorrelators;i++) {
for (unsigned int j=0;j<p;j++)
for (int i=0; i < numcorrelators; i++) {
for (unsigned int j=0; j < p; j++)
ncorrelation[i][j] = static_cast<unsigned long int>(list[n++]);
naccumulator[i] = static_cast<unsigned int> (list[n++]);
insertindex[i] = static_cast<unsigned int> (list[n++]);
naccumulator[i] = static_cast<unsigned int>(list[n++]);
insertindex[i] = static_cast<unsigned int>(list[n++]);
}
}

19
src/EXTRA-FIX/fix_ave_correlate_long.h
View File

@ -58,18 +58,27 @@ class FixAveCorrelateLong : public Fix {
int length; // Length of result arrays
int kmax; // Maximum correlator attained during simulation
int me, nvalues;
int nfreq;
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, nfreq;
bigint nvalid, nvalid_last, last_accumulated_step;
int *which, *argindex, *value2index;
char **ids;
FILE *fp;
int type, startstep, overwrite;
bigint filepos;
int npair; // number of correlation pairs to calculate
double *values;
double *cvalues;
void accumulate();
void evaluate();

5
src/compute_chunk_spread_atom.cpp
View File

@ -63,8 +63,7 @@ ComputeChunkSpreadAtom(LAMMPS *lmp, int narg, char **arg) :
val.id = argi.get_name();
val.val.c = nullptr;
if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE)
|| (argi.get_dim() > 1))
if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) || (argi.get_dim() > 1))
error->all(FLERR,"Illegal compute chunk/spread/atom argument: {}", arg[iarg]);
values.push_back(val);
@ -153,7 +152,7 @@ void ComputeChunkSpreadAtom::init()
{
init_chunk();
// set indices of all computes,fixes,variables
// store references of all computes and fixes
for (auto &val : values) {
if (val.which == ArgInfo::COMPUTE) {

406
src/compute_global_atom.cpp
View File

@ -26,19 +26,12 @@
using namespace LAMMPS_NS;
enum{VECTOR,ARRAY};
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
ComputeGlobalAtom::ComputeGlobalAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
idref(nullptr), which(nullptr), argindex(nullptr), value2index(nullptr), ids(nullptr),
indices(nullptr), varatom(nullptr), vecglobal(nullptr)
Compute(lmp, narg, arg), indices(nullptr), varatom(nullptr), vecglobal(nullptr)
{
if (narg < 5) error->all(FLERR,"Illegal compute global/atom command");
if (narg < 5) utils::missing_cmd_args(FLERR,"compute global/atom", error);
peratom_flag = 1;
@ -47,13 +40,13 @@ ComputeGlobalAtom::ComputeGlobalAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = 3;
ArgInfo argi(arg[iarg]);
whichref = argi.get_type();
indexref = argi.get_index1();
idref = argi.copy_name();
reference.which = argi.get_type();
reference.argindex = argi.get_index1();
reference.id = argi.get_name();
if ((whichref == ArgInfo::UNKNOWN) || (whichref == ArgInfo::NONE)
if ((reference.which == ArgInfo::UNKNOWN) || (reference.which == ArgInfo::NONE)
|| (argi.get_dim() > 1))
error->all(FLERR,"Illegal compute global/atom command");
error->all(FLERR,"Illegal compute global/atom index property: {}", arg[iarg]);
iarg++;
@ -68,24 +61,21 @@ ComputeGlobalAtom::ComputeGlobalAtom(LAMMPS *lmp, int narg, char **arg) :
// parse values until one isn't recognized
which = new int[nargnew];
argindex = new int[nargnew];
ids = new char*[nargnew];
value2index = new int[nargnew];
nvalues = 0;
values.clear();
for (iarg = 0; iarg < nargnew; iarg++) {
ArgInfo argi2(arg[iarg]);
which[nvalues] = argi2.get_type();
argindex[nvalues] = argi2.get_index1();
ids[nvalues] = argi2.copy_name();
value_t val;
val.which = argi2.get_type();
val.argindex = argi2.get_index1();
val.id = argi2.get_name();
if ((which[nvalues] == ArgInfo::UNKNOWN) || (which[nvalues] == ArgInfo::NONE)
if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE)
|| (argi2.get_dim() > 1))
error->all(FLERR,"Illegal compute slice command");
error->all(FLERR,"Illegal compute global/atom global property: {}", arg[iarg]);
nvalues++;
values.push_back(val);
}
// if wildcard expansion occurred, free earg memory from expand_args()
@ -95,94 +85,97 @@ ComputeGlobalAtom::ComputeGlobalAtom(LAMMPS *lmp, int narg, char **arg) :
memory->sfree(earg);
}
// setup and error check both index arg and values
// setup and error check for both, index arg and values
if (whichref == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(idref);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute global/atom does not exist");
if (reference.which == ArgInfo::COMPUTE) {
reference.val.c = modify->get_compute_by_id(reference.id);
if (!reference.val.c)
error->all(FLERR,"Compute ID {} for compute global/atom index", reference.id);
if (!modify->compute[icompute]->peratom_flag)
error->all(FLERR,"Compute global/atom compute does not "
"calculate a per-atom vector or array");
if (indexref == 0 &&
modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,"Compute global/atom compute does not "
"calculate a per-atom vector");
if (indexref && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,"Compute global/atom compute does not "
"calculate a per-atom array");
if (indexref && indexref > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR,
"Compute global/atom compute array is accessed out-of-range");
if (!reference.val.c->peratom_flag)
error->all(FLERR,"Compute global/atom compute {} does not calculate a per-atom "
"vector or array", reference.id);
if ((reference.argindex == 0) && (reference.val.c->size_peratom_cols != 0))
error->all(FLERR,"Compute global/atom compute {} does not calculate a per-atom "
"vector", reference.id);
if (reference.argindex && (reference.val.c->size_peratom_cols == 0))
error->all(FLERR,"Compute global/atom compute does not calculate a per-atom "
"array", reference.id);
if (reference.argindex && (reference.argindex > reference.val.c->size_peratom_cols))
error->all(FLERR, "Compute global/atom compute array {} is accessed out-of-range",
reference.id);
} else if (whichref == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(idref);
if (!ifix)
error->all(FLERR,"Fix ID {} for compute global/atom does not exist", idref);
if (!ifix->peratom_flag)
error->all(FLERR,"Compute global/atom fix {} does not calculate a per-atom vector or array", idref);
if (indexref == 0 && (ifix->size_peratom_cols != 0))
error->all(FLERR,"Compute global/atom fix {} does not calculate a per-atom vector", idref);
if (indexref && (ifix->size_peratom_cols == 0))
error->all(FLERR,"Compute global/atom fix {} does not calculate a per-atom array", idref);
if (indexref && indexref > ifix->size_peratom_cols)
error->all(FLERR, "Compute global/atom fix {} array is accessed out-of-range", idref);
} else if (reference.which == ArgInfo::FIX) {
reference.val.f =modify->get_fix_by_id(reference.id);
if (!reference.val.f)
error->all(FLERR,"Fix ID {} for compute global/atom does not exist", reference.id);
if (!reference.val.f->peratom_flag)
error->all(FLERR,"Compute global/atom fix {} does not calculate a per-atom vector "
"or array", reference.id);
if (reference.argindex == 0 && (reference.val.f->size_peratom_cols != 0))
error->all(FLERR,"Compute global/atom fix {} does not calculate a per-atom vector",
reference.id);
if (reference.argindex && (reference.val.f->size_peratom_cols == 0))
error->all(FLERR,"Compute global/atom fix {} does not calculate a per-atom array",
reference.id);
if (reference.argindex && (reference.argindex > reference.val.f->size_peratom_cols))
error->all(FLERR, "Compute global/atom fix {} array is accessed out-of-range", reference.id);
} else if (whichref == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(idref);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute global/atom does not exist");
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Compute global/atom variable is not "
"atom-style variable");
} else if (reference.which == ArgInfo::VARIABLE) {
reference.val.v = input->variable->find(reference.id.c_str());
if (reference.val.v < 0)
error->all(FLERR,"Variable name {} for compute global/atom index does not exist",
reference.id);
if (input->variable->atomstyle(reference.val.v) == 0)
error->all(FLERR,"Compute global/atom index variable {} is not atom-style variable",
reference.id);
}
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 compute global/atom does not exist");
if (argindex[i] == 0) {
if (!modify->compute[icompute]->vector_flag)
error->all(FLERR,"Compute global/atom compute does not "
"calculate a global vector");
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 compute global/atom does not exist", val.id);
if (val.argindex == 0) {
if (!val.val.c->vector_flag)
error->all(FLERR,"Compute ID {} for global/atom compute does not calculate "
"a global vector", val.id);
} else {
if (!modify->compute[icompute]->array_flag)
error->all(FLERR,"Compute global/atom compute does not "
"calculate a global array");
if (argindex[i] > modify->compute[icompute]->size_array_cols)
error->all(FLERR,"Compute global/atom compute array is "
"accessed out-of-range");
if (!val.val.c->array_flag)
error->all(FLERR,"Compute ID {} for global/atom compute does not calculate "
"a global array", val.id);
if (val.argindex > val.val.c->size_array_cols)
error->all(FLERR,"Compute global/atom compute {} array is accessed out-of-range", val.id);
}
} else if (which[i] == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(ids[i]);
if (!ifix)
error->all(FLERR,"Fix ID for compute global/atom does not exist");
if (argindex[i] == 0) {
if (!ifix->vector_flag)
error->all(FLERR,"Compute global/atom fix {} does not calculate a global vector", ids[i]);
} 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 compute global/atom does not exist", val.id);
if (val.argindex == 0) {
if (!val.val.f->vector_flag)
error->all(FLERR,"Fix ID {} for compute global/atom compute does not calculate "
"a global vector", val.id);
} else {
if (!ifix->array_flag)
error->all(FLERR,"Compute global/atom fix {} does not calculate a global array", ids[i]);
if (argindex[i] > ifix->size_array_cols)
error->all(FLERR,"Compute global/atom fix {} array is accessed out-of-range", ids[i]);
if (!val.val.f->array_flag)
error->all(FLERR,"Fix ID {} for compute global/atom compute does not calculate "
"a global array", val.id);
if (val.argindex > val.val.f->size_array_cols)
error->all(FLERR,"Compute global/atom fix {} array is accessed out-of-range", val.id);
}
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute global/atom does not exist");
if (input->variable->vectorstyle(ivariable) == 0)
error->all(FLERR,"Compute global/atom variable is not vector-style 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 compute global/atom does not exist", val.id);
if (input->variable->vectorstyle(val.val.v) == 0)
error->all(FLERR,"Compute global/atom variable {} is not vector-style variable", val.id);
}
}
// this compute produces either a peratom vector or array
if (nvalues == 1) size_peratom_cols = 0;
else size_peratom_cols = nvalues;
if (values.size() == 1) size_peratom_cols = 0;
else size_peratom_cols = values.size();
nmax = maxvector = 0;
vector_atom = nullptr;
@ -193,14 +186,6 @@ ComputeGlobalAtom::ComputeGlobalAtom(LAMMPS *lmp, int narg, char **arg) :
ComputeGlobalAtom::~ComputeGlobalAtom()
{
delete [] idref;
delete [] which;
delete [] argindex;
for (int m = 0; m < nvalues; m++) delete [] ids[m];
delete [] ids;
delete [] value2index;
memory->destroy(indices);
memory->destroy(varatom);
memory->destroy(vecglobal);
@ -212,44 +197,38 @@ ComputeGlobalAtom::~ComputeGlobalAtom()
void ComputeGlobalAtom::init()
{
// set indices of all computes,fixes,variables
// store references of all computes, fixes, or variables
if (whichref == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(idref);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute global/atom does not exist");
ref2index = icompute;
} else if (whichref == ArgInfo::FIX) {
int ifix = modify->find_fix(idref);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute global/atom does not exist");
ref2index = ifix;
} else if (whichref == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(idref);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute global/atom does not exist");
ref2index = ivariable;
if (reference.which == ArgInfo::COMPUTE) {
reference.val.c = modify->get_compute_by_id(reference.id);
if (!reference.val.c)
error->all(FLERR,"Compute ID {} for compute global/atom index does not exist", reference.id);
} else if (reference.which == ArgInfo::FIX) {
reference.val.f = modify->get_fix_by_id(reference.id);
if (reference.val.f)
error->all(FLERR,"Fix ID {} for compute global/atom index does not exist", reference.id);
} else if (reference.which == ArgInfo::VARIABLE) {
reference.val.v = input->variable->find(reference.id.c_str());
if (reference.val.v < 0)
error->all(FLERR,"Variable name {} for compute global/atom index does not exist",
reference.id);
}
for (int m = 0; m < nvalues; m++) {
if (which[m] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute global/atom does not exist");
value2index[m] = icompute;
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 compute global/atom does not exist", val.id);
} else if (which[m] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute global/atom does not exist");
value2index[m] = ifix;
} 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 compute global/atom does not exist", val.id);
} else if (which[m] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute global/atom "
"does not exist");
value2index[m] = ivariable;
} 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 compute global/atom does not exist", val.id);
}
}
}
@ -258,8 +237,6 @@ void ComputeGlobalAtom::init()
void ComputeGlobalAtom::compute_peratom()
{
int i,j;
invoked_peratom = update->ntimestep;
// grow indices and output vector or array if necessary
@ -268,16 +245,16 @@ void ComputeGlobalAtom::compute_peratom()
nmax = atom->nmax;
memory->destroy(indices);
memory->create(indices,nmax,"global/atom:indices");
if (whichref == ArgInfo::VARIABLE) {
if (reference.which == ArgInfo::VARIABLE) {
memory->destroy(varatom);
memory->create(varatom,nmax,"global/atom:varatom");
}
if (nvalues == 1) {
if (values.size() == 1) {
memory->destroy(vector_atom);
memory->create(vector_atom,nmax,"global/atom:vector_atom");
} else {
memory->destroy(array_atom);
memory->create(array_atom,nmax,nvalues,"global/atom:array_atom");
memory->create(array_atom,nmax,values.size(),"global/atom:array_atom");
}
}
@ -286,81 +263,79 @@ void ComputeGlobalAtom::compute_peratom()
// indices are decremented from 1 to N -> 0 to N-1
int *mask = atom->mask;
int nlocal = atom->nlocal;
const int nlocal = atom->nlocal;
if (whichref == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[ref2index];
if (reference.which == ArgInfo::COMPUTE) {
if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM;
if (!(reference.val.c->invoked_flag & Compute::INVOKED_PERATOM)) {
reference.val.c->compute_peratom();
reference.val.c->invoked_flag |= Compute::INVOKED_PERATOM;
}
if (indexref == 0) {
double *compute_vector = compute->vector_atom;
for (i = 0; i < nlocal; i++)
if (reference.argindex == 0) {
double *compute_vector = reference.val.c->vector_atom;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
indices[i] = static_cast<int> (compute_vector[i]) - 1;
} else {
double **compute_array = compute->array_atom;
int im1 = indexref - 1;
for (i = 0; i < nlocal; i++)
double **compute_array = reference.val.c->array_atom;
int im1 = reference.argindex - 1;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
indices[i] = static_cast<int> (compute_array[i][im1]) - 1;
}
} else if (whichref == ArgInfo::FIX) {
if (update->ntimestep % modify->fix[ref2index]->peratom_freq)
error->all(FLERR,"Fix used in compute global/atom not "
"computed at compatible time");
Fix *fix = modify->fix[ref2index];
} else if (reference.which == ArgInfo::FIX) {
if (update->ntimestep % reference.val.f->peratom_freq)
error->all(FLERR,"Fix {} used in compute global/atom not computed at compatible time",
reference.id);
if (indexref == 0) {
double *fix_vector = fix->vector_atom;
for (i = 0; i < nlocal; i++)
if (reference.argindex == 0) {
double *fix_vector = reference.val.f->vector_atom;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
indices[i] = static_cast<int> (fix_vector[i]) - 1;
} else {
double **fix_array = fix->array_atom;
int im1 = indexref - 1;
for (i = 0; i < nlocal; i++)
double **fix_array = reference.val.f->array_atom;
int im1 = reference.argindex - 1;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
indices[i] = static_cast<int> (fix_array[i][im1]) - 1;
}
} else if (whichref == ArgInfo::VARIABLE) {
input->variable->compute_atom(ref2index,igroup,varatom,1,0);
for (i = 0; i < nlocal; i++)
} else if (reference.which == ArgInfo::VARIABLE) {
input->variable->compute_atom(reference.val.v, igroup, varatom, 1, 0);
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
indices[i] = static_cast<int> (varatom[i]) - 1;
}
// loop over values to fill output vector or array
for (int m = 0; m < nvalues; m++) {
int m = 0;
for (auto &val : values) {
// output = vector
if (argindex[m] == 0) {
if (val.argindex == 0) {
int vmax;
double *source;
if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[value2index[m]];
if (val.which == ArgInfo::COMPUTE) {
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;
}
source = compute->vector;
vmax = compute->size_vector;
source = val.val.c->vector;
vmax = val.val.c->size_vector;
} else if (which[m] == ArgInfo::FIX) {
if (update->ntimestep % modify->fix[value2index[m]]->peratom_freq)
error->all(FLERR,"Fix used in compute global/atom not computed at compatible time");
Fix *fix = modify->fix[value2index[m]];
vmax = fix->size_vector;
} else if (val.which == ArgInfo::FIX) {
if (update->ntimestep % val.val.f->peratom_freq)
error->all(FLERR,"Fix {} used in compute global/atom not computed at compatible time",
val.id);
vmax = reference.val.f->size_vector;
if (vmax > maxvector) {
maxvector = vmax;
@ -368,28 +343,28 @@ void ComputeGlobalAtom::compute_peratom()
memory->create(vecglobal,maxvector,"global/atom:vecglobal");
}
for (i = 0; i < vmax; i++)
vecglobal[i] = fix->compute_vector(i);
for (int i = 0; i < vmax; i++)
vecglobal[i] = val.val.f->compute_vector(i);
source = vecglobal;
} else if (which[m] == ArgInfo::VARIABLE) {
vmax = input->variable->compute_vector(value2index[m],&source);
} else if (val.which == ArgInfo::VARIABLE) {
vmax = input->variable->compute_vector(val.val.v, &source);
}
if (nvalues == 1) {
for (i = 0; i < nlocal; i++) {
if (values.size() == 1) {
for (int i = 0; i < nlocal; i++) {
vector_atom[i] = 0.0;
if (mask[i] & groupbit) {
j = indices[i];
int j = indices[i];
if (j >= 0 && j < vmax) vector_atom[i] = source[j];
}
}
} else {
for (i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) {
array_atom[i][m] = 0.0;
if (mask[i] & groupbit) {
j = indices[i];
int j = indices[i];
if (j >= 0 && j < vmax) array_atom[i][m] = source[j];
}
}
@ -400,18 +375,17 @@ void ComputeGlobalAtom::compute_peratom()
} else {
int vmax;
double *source;
int col = argindex[m] - 1;
int col = val.argindex - 1;
if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[value2index[m]];
if (val.which == ArgInfo::COMPUTE) {
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;
}
double **compute_array = compute->array;
vmax = compute->size_array_rows;
double **compute_array = val.val.c->array;
vmax = val.val.c->size_array_rows;
if (vmax > maxvector) {
maxvector = vmax;
@ -419,17 +393,16 @@ void ComputeGlobalAtom::compute_peratom()
memory->create(vecglobal,maxvector,"global/atom:vecglobal");
}
for (i = 0; i < vmax; i++)
for (int i = 0; i < vmax; i++)
vecglobal[i] = compute_array[i][col];
source = vecglobal;
} else if (which[m] == ArgInfo::FIX) {
if (update->ntimestep % modify->fix[value2index[m]]->peratom_freq)
error->all(FLERR,"Fix used in compute global/atom not "
"computed at compatible time");
Fix *fix = modify->fix[value2index[m]];
vmax = fix->size_array_rows;
} else if (val.which == ArgInfo::FIX) {
if (update->ntimestep % val.val.f->peratom_freq)
error->all(FLERR,"Fix {} used in compute global/atom not computed at compatible time",
val.id);
vmax = val.val.f->size_array_rows;
if (vmax > maxvector) {
maxvector = vmax;
@ -437,34 +410,35 @@ void ComputeGlobalAtom::compute_peratom()
memory->create(vecglobal,maxvector,"global/atom:vecglobal");
}
for (i = 0; i < vmax; i++)
vecglobal[i] = fix->compute_array(i,col);
for (int i = 0; i < vmax; i++)
vecglobal[i] = val.val.f->compute_array(i,col);
source = vecglobal;
} else if (which[m] == ArgInfo::VARIABLE) {
vmax = input->variable->compute_vector(value2index[m],&source);
} else if (val.which == ArgInfo::VARIABLE) {
vmax = input->variable->compute_vector(val.val.v, &source);
}
if (nvalues == 1) {
for (i = 0; i < nlocal; i++) {
if (values.size() == 1) {
for (int i = 0; i < nlocal; i++) {
vector_atom[i] = 0.0;
if (mask[i] & groupbit) {
j = indices[i];
int j = indices[i];
if (j >= 0 && j < vmax) vector_atom[i] = source[j];
}
}
} else {
for (i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) {
array_atom[i][m] = 0.0;
if (mask[i] & groupbit) {
j = indices[i];
int j = indices[i];
if (j >= 0 && j < vmax) array_atom[i][m] = source[j];
}
}
}
}
}
++m;
}
/* ----------------------------------------------------------------------
@ -473,7 +447,7 @@ void ComputeGlobalAtom::compute_peratom()
double ComputeGlobalAtom::memory_usage()
{
double bytes = (double)nmax*nvalues * sizeof(double);
double bytes = (double)nmax*values.size() * sizeof(double);
bytes += (double)nmax * sizeof(int); // indices
if (varatom) bytes += (double)nmax * sizeof(double); // varatom
bytes += (double)maxvector * sizeof(double); // vecglobal

18
src/compute_global_atom.h
View File

@ -33,12 +33,18 @@ class ComputeGlobalAtom : public Compute {
double memory_usage() override;
protected:
int whichref, indexref, ref2index;
char *idref;
int nvalues;
int *which, *argindex, *value2index;
char **ids;
struct value_t {
int which;
int argindex;
std::string id;
union {
class Compute *c;
class Fix *f;
int v;
} val;
};
std::vector<value_t> values;
value_t reference;
int nmax, maxvector;
int *indices;

396
src/compute_reduce.cpp
View File

@ -15,6 +15,7 @@
#include "arg_info.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
@ -34,16 +35,15 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), nvalues(0), which(nullptr), argindex(nullptr), flavor(nullptr),
value2index(nullptr), ids(nullptr), onevec(nullptr), replace(nullptr), indices(nullptr),
Compute(lmp, narg, arg), nvalues(0), onevec(nullptr), replace(nullptr), indices(nullptr),
owner(nullptr), idregion(nullptr), region(nullptr), varatom(nullptr)
{
int iarg = 0;
if (strcmp(style, "reduce") == 0) {
if (narg < 5) error->all(FLERR, "Illegal compute reduce command");
if (narg < 5) utils::missing_cmd_args(FLERR, "compute reduce", error);
iarg = 3;
} else if (strcmp(style, "reduce/region") == 0) {
if (narg < 6) error->all(FLERR, "Illegal compute reduce/region command");
if (narg < 6) utils::missing_cmd_args(FLERR, "compute reduce/region", error);
if (!domain->get_region_by_id(arg[3]))
error->all(FLERR, "Region {} for compute reduce/region does not exist", arg[3]);
idregion = utils::strdup(arg[3]);
@ -67,11 +67,9 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg], "aveabs") == 0)
mode = AVEABS;
else
error->all(FLERR, "Illegal compute {} operation {}", style, arg[iarg]);
error->all(FLERR, "Unknown compute {} mode: {}", style, arg[iarg]);
iarg++;
MPI_Comm_rank(world, &me);
// expand args if any have wildcard character "*"
int expand = 0;
@ -81,95 +79,92 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
if (earg != &arg[iarg]) expand = 1;
arg = earg;
// parse values until one isn't recognized
// parse values
which = new int[nargnew];
argindex = new int[nargnew];
flavor = new int[nargnew];
ids = new char *[nargnew];
value2index = new int[nargnew];
for (int i = 0; i < nargnew; ++i) {
which[i] = argindex[i] = flavor[i] = value2index[i] = ArgInfo::UNKNOWN;
ids[i] = nullptr;
}
values.clear();
nvalues = 0;
for (int iarg = 0; iarg < nargnew; ++iarg) {
value_t val;
iarg = 0;
while (iarg < nargnew) {
ids[nvalues] = nullptr;
val.id = "";
val.flavor = 0;
val.val.c = nullptr;
if (strcmp(arg[iarg], "x") == 0) {
which[nvalues] = ArgInfo::X;
argindex[nvalues++] = 0;
val.which = ArgInfo::X;
val.argindex = 0;
} else if (strcmp(arg[iarg], "y") == 0) {
which[nvalues] = ArgInfo::X;
argindex[nvalues++] = 1;
val.which = ArgInfo::X;
val.argindex = 1;
} else if (strcmp(arg[iarg], "z") == 0) {
which[nvalues] = ArgInfo::X;
argindex[nvalues++] = 2;
val.which = ArgInfo::X;
val.argindex = 2;
} else if (strcmp(arg[iarg], "vx") == 0) {
which[nvalues] = ArgInfo::V;
argindex[nvalues++] = 0;
val.which = ArgInfo::V;
val.argindex = 0;
} else if (strcmp(arg[iarg], "vy") == 0) {
which[nvalues] = ArgInfo::V;
argindex[nvalues++] = 1;
val.which = ArgInfo::V;
val.argindex = 1;
} else if (strcmp(arg[iarg], "vz") == 0) {
which[nvalues] = ArgInfo::V;
argindex[nvalues++] = 2;
val.which = ArgInfo::V;
val.argindex = 2;
} else if (strcmp(arg[iarg], "fx") == 0) {
which[nvalues] = ArgInfo::F;
argindex[nvalues++] = 0;
val.which = ArgInfo::F;
val.argindex = 0;
} else if (strcmp(arg[iarg], "fy") == 0) {
which[nvalues] = ArgInfo::F;
argindex[nvalues++] = 1;
val.which = ArgInfo::F;
val.argindex = 1;
} else if (strcmp(arg[iarg], "fz") == 0) {
which[nvalues] = ArgInfo::F;
argindex[nvalues++] = 2;
val.which = ArgInfo::F;
val.argindex = 2;
} else {
ArgInfo argi(arg[iarg]);
which[nvalues] = argi.get_type();
argindex[nvalues] = argi.get_index1();
ids[nvalues] = argi.copy_name();
val.which = argi.get_type();
val.argindex = argi.get_index1();
val.id = argi.get_name();
if ((which[nvalues] == ArgInfo::UNKNOWN) || (argi.get_dim() > 1))
error->all(FLERR, "Illegal compute reduce command");
if ((val.which == ArgInfo::UNKNOWN) || (argi.get_dim() > 1))
error->all(FLERR, "Illegal compute {} argument: {}", style, arg[iarg]);
if (which[nvalues] == ArgInfo::NONE) break;
nvalues++;
if (val.which == ArgInfo::NONE) break;
}
iarg++;
values.push_back(val);
}
// optional args
nvalues = values.size();
replace = new int[nvalues];
for (int i = 0; i < nvalues; i++) replace[i] = -1;
for (int i = 0; i < nvalues; ++i) replace[i] = -1;
std::string mycmd = "compute ";
mycmd += style;
while (iarg < nargnew) {
for (int iarg = nvalues; iarg < nargnew; iarg++) {
if (strcmp(arg[iarg], "replace") == 0) {
if (iarg + 3 > narg) error->all(FLERR, "Illegal compute reduce command");
if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, mycmd + " replace", error);
if (mode != MINN && mode != MAXX)
error->all(FLERR, "Compute reduce replace requires min or max mode");
error->all(FLERR, "Compute {} replace requires min or max mode", style);
int col1 = utils::inumeric(FLERR, arg[iarg + 1], false, lmp) - 1;
int col2 = utils::inumeric(FLERR, arg[iarg + 2], false, lmp) - 1;
if (col1 < 0 || col1 >= nvalues || col2 < 0 || col2 >= nvalues)
error->all(FLERR, "Illegal compute reduce command");
if (col1 == col2) error->all(FLERR, "Illegal compute reduce command");
if (replace[col1] >= 0 || replace[col2] >= 0)
error->all(FLERR, "Invalid replace values in compute reduce");
if ((col1 < 0) || (col1 >= nvalues))
error->all(FLERR, "Invalid compute {} replace first column index {}", style, col1);
if ((col2 < 0) || (col2 >= nvalues))
error->all(FLERR, "Invalid compute {} replace second column index {}", style, col2);
if (col1 == col2) error->all(FLERR, "Compute {} replace columns must be different");
if ((replace[col1] >= 0) || (replace[col2] >= 0))
error->all(FLERR, "Compute {} replace column already used for another replacement");
replace[col1] = col2;
iarg += 3;
iarg += 2;
} else
error->all(FLERR, "Illegal compute reduce command");
error->all(FLERR, "Unknown compute {} keyword: {}", style, arg[iarg]);
}
// delete replace if not set
// delete replace list if not set
int flag = 0;
for (int i = 0; i < nvalues; i++)
@ -188,68 +183,67 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
// setup and error check
for (int i = 0; i < nvalues; i++) {
if (which[i] == ArgInfo::X || which[i] == ArgInfo::V || which[i] == ArgInfo::F)
flavor[i] = PERATOM;
for (auto &val : values) {
if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F)
val.flavor = PERATOM;
else if (which[i] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[i]);
if (icompute < 0) error->all(FLERR, "Compute ID for compute reduce does not exist");
if (modify->compute[icompute]->peratom_flag) {
flavor[i] = PERATOM;
if (argindex[i] == 0 && modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,
"Compute reduce compute does not "
"calculate a per-atom vector");
if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,
"Compute reduce compute does not "
"calculate a per-atom array");
if (argindex[i] && argindex[i] > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR, "Compute reduce compute array is accessed out-of-range");
} else if (modify->compute[icompute]->local_flag) {
flavor[i] = LOCAL;
if (argindex[i] == 0 && modify->compute[icompute]->size_local_cols != 0)
error->all(FLERR,
"Compute reduce compute does not "
"calculate a local vector");
if (argindex[i] && modify->compute[icompute]->size_local_cols == 0)
error->all(FLERR,
"Compute reduce compute does not "
"calculate a local array");
if (argindex[i] && argindex[i] > modify->compute[icompute]->size_local_cols)
error->all(FLERR, "Compute reduce compute array is accessed out-of-range");
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 compute {} does not exist", val.id, style);
if (val.val.c->peratom_flag) {
val.flavor = PERATOM;
if (val.argindex == 0 && val.val.c->size_peratom_cols != 0)
error->all(FLERR, "Compute {} compute {} does not calculate a per-atom vector", style,
val.id);
if (val.argindex && val.val.c->size_peratom_cols == 0)
error->all(FLERR, "Compute {} compute {} does not calculate a per-atom array", style,
val.id);
if (val.argindex && val.argindex > val.val.c->size_peratom_cols)
error->all(FLERR, "Compute {} compute {} array is accessed out-of-range", style, val.id);
} else if (val.val.c->local_flag) {
val.flavor = LOCAL;
if (val.argindex == 0 && val.val.c->size_local_cols != 0)
error->all(FLERR, "Compute {} compute {} does not calculate a local vector", style,
val.id);
if (val.argindex && val.val.c->size_local_cols == 0)
error->all(FLERR, "Compute {} compute {} does not calculate a local array", style,
val.id);
if (val.argindex && val.argindex > val.val.c->size_local_cols)
error->all(FLERR, "Compute {} compute {} array is accessed out-of-range", style, val.id);
} else
error->all(FLERR, "Compute reduce compute calculates global values");
error->all(FLERR, "Compute {} compute {} calculates global values", style, val.id);
} else if (which[i] == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(ids[i]);
if (!ifix) error->all(FLERR, "Fix ID {} for compute reduce does not exist", ids[i]);
if (ifix->peratom_flag) {
flavor[i] = PERATOM;
if (argindex[i] == 0 && (ifix->size_peratom_cols != 0))
error->all(FLERR, "Compute reduce fix {} does not calculate a per-atom vector", ids[i]);
if (argindex[i] && (ifix->size_peratom_cols == 0))
error->all(FLERR, "Compute reduce fix {} does not calculate a per-atom array", ids[i]);
if (argindex[i] && (argindex[i] > ifix->size_peratom_cols))
error->all(FLERR, "Compute reduce fix {} array is accessed out-of-range", ids[i]);
} else if (ifix->local_flag) {
flavor[i] = LOCAL;
if (argindex[i] == 0 && (ifix->size_local_cols != 0))
error->all(FLERR, "Compute reduce fix {} does not calculate a local vector", ids[i]);
if (argindex[i] && (ifix->size_local_cols == 0))
error->all(FLERR, "Compute reduce fix {} does not calculate a local array", ids[i]);
if (argindex[i] && (argindex[i] > ifix->size_local_cols))
error->all(FLERR, "Compute reduce fix {} array is accessed out-of-range", ids[i]);
} 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 compute {} does not exist", val.id, style);
if (val.val.f->peratom_flag) {
val.flavor = PERATOM;
if (val.argindex == 0 && (val.val.f->size_peratom_cols != 0))
error->all(FLERR, "Compute {} fix {} does not calculate a per-atom vector", style,
val.id);
if (val.argindex && (val.val.f->size_peratom_cols == 0))
error->all(FLERR, "Compute {} fix {} does not calculate a per-atom array", style, val.id);
if (val.argindex && (val.argindex > val.val.f->size_peratom_cols))
error->all(FLERR, "Compute {} fix {} array is accessed out-of-range", style, val.id);
} else if (val.val.f->local_flag) {
val.flavor = LOCAL;
if (val.argindex == 0 && (val.val.f->size_local_cols != 0))
error->all(FLERR, "Compute {} fix {} does not calculate a local vector", style, val.id);
if (val.argindex && (val.val.f->size_local_cols == 0))
error->all(FLERR, "Compute {} fix {} does not calculate a local array", style, val.id);
if (val.argindex && (val.argindex > val.val.f->size_local_cols))
error->all(FLERR, "Compute {} fix {} array is accessed out-of-range", style, val.id);
} else
error->all(FLERR, "Compute reduce fix {} calculates global values", ids[i]);
error->all(FLERR, "Compute {} fix {} calculates global values", style, val.id);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0) error->all(FLERR, "Variable name for compute reduce does not exist");
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR, "Compute reduce variable is not atom-style variable");
flavor[i] = PERATOM;
} 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 compute {} does not exist", val.id, style);
if (input->variable->atomstyle(val.val.v) == 0)
error->all(FLERR, "Compute {} variable {} is not atom-style variable", style, val.id);
val.flavor = PERATOM;
}
}
@ -284,12 +278,6 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
ComputeReduce::~ComputeReduce()
{
delete[] which;
delete[] argindex;
delete[] flavor;
for (int m = 0; m < nvalues; m++) delete[] ids[m];
delete[] ids;
delete[] value2index;
delete[] replace;
delete[] idregion;
@ -307,24 +295,21 @@ void ComputeReduce::init()
{
// set indices of all computes,fixes,variables
for (int m = 0; m < nvalues; m++) {
if (which[m] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0) error->all(FLERR, "Compute ID for compute reduce does not exist");
value2index[m] = icompute;
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 compute {} does not exist", val.id, style);
} else if (which[m] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0) error->all(FLERR, "Fix ID for compute reduce does not exist");
value2index[m] = ifix;
} 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 compute {} does not exist", val.id, style);
} else if (which[m] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0) error->all(FLERR, "Variable name for compute reduce does not exist");
value2index[m] = ivariable;
} else
value2index[m] = ArgInfo::UNKNOWN;
} 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 compute {} does not exist", val.id, style);
}
}
// set index and check validity of region
@ -383,14 +368,14 @@ void ComputeReduce::compute_vector()
for (int m = 0; m < nvalues; m++)
if (replace[m] < 0) {
pairme.value = onevec[m];
pairme.proc = me;
pairme.proc = comm->me;
MPI_Allreduce(&pairme, &pairall, 1, MPI_DOUBLE_INT, MPI_MINLOC, world);
vector[m] = pairall.value;
owner[m] = pairall.proc;
}
for (int m = 0; m < nvalues; m++)
if (replace[m] >= 0) {
if (me == owner[replace[m]]) vector[m] = compute_one(m, indices[replace[m]]);
if (comm->me == owner[replace[m]]) vector[m] = compute_one(m, indices[replace[m]]);
MPI_Bcast(&vector[m], 1, MPI_DOUBLE, owner[replace[m]], world);
}
}
@ -404,14 +389,14 @@ void ComputeReduce::compute_vector()
for (int m = 0; m < nvalues; m++)
if (replace[m] < 0) {
pairme.value = onevec[m];
pairme.proc = me;
pairme.proc = comm->me;
MPI_Allreduce(&pairme, &pairall, 1, MPI_DOUBLE_INT, MPI_MAXLOC, world);
vector[m] = pairall.value;
owner[m] = pairall.proc;
}
for (int m = 0; m < nvalues; m++)
if (replace[m] >= 0) {
if (me == owner[replace[m]]) vector[m] = compute_one(m, indices[replace[m]]);
if (comm->me == owner[replace[m]]) vector[m] = compute_one(m, indices[replace[m]]);
MPI_Bcast(&vector[m], 1, MPI_DOUBLE, owner[replace[m]], world);
}
}
@ -436,24 +421,21 @@ void ComputeReduce::compute_vector()
double ComputeReduce::compute_one(int m, int flag)
{
int i;
// invoke the appropriate attribute,compute,fix,variable
// for flag = -1, compute scalar quantity by scanning over atom properties
// only include atoms in group for atom properties and per-atom quantities
index = -1;
int vidx = value2index[m];
auto &val = values[m];
// initialization in case it has not yet been run, e.g. when
// the compute was invoked right after it has been created
if (vidx == ArgInfo::UNKNOWN) {
init();
vidx = value2index[m];
if ((val.which == ArgInfo::COMPUTE) || (val.which == ArgInfo::FIX)) {
if (val.val.c == nullptr) init();
}
int aidx = argindex[m];
int aidx = val.argindex;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -461,77 +443,76 @@ double ComputeReduce::compute_one(int m, int flag)
if (mode == MINN) one = BIG;
if (mode == MAXX) one = -BIG;
if (which[m] == ArgInfo::X) {
if (val.which == ArgInfo::X) {
double **x = atom->x;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one, x[i][aidx], i);
} else
one = x[flag][aidx];
} else if (which[m] == ArgInfo::V) {
} else if (val.which == ArgInfo::V) {
double **v = atom->v;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one, v[i][aidx], i);
} else
one = v[flag][aidx];
} else if (which[m] == ArgInfo::F) {
} else if (val.which == ArgInfo::F) {
double **f = atom->f;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one, f[i][aidx], i);
} else
one = f[flag][aidx];
// invoke compute if not previously invoked
} else if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[vidx];
} else if (val.which == ArgInfo::COMPUTE) {
if (flavor[m] == PERATOM) {
if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM;
if (val.flavor == PERATOM) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) {
val.val.c->compute_peratom();
val.val.c->invoked_flag |= Compute::INVOKED_PERATOM;
}
if (aidx == 0) {
double *comp_vec = compute->vector_atom;
double *comp_vec = val.val.c->vector_atom;
int n = nlocal;
if (flag < 0) {
for (i = 0; i < n; i++)
for (int i = 0; i < n; i++)
if (mask[i] & groupbit) combine(one, comp_vec[i], i);
} else
one = comp_vec[flag];
} else {
double **carray_atom = compute->array_atom;
double **carray_atom = val.val.c->array_atom;
int n = nlocal;
int aidxm1 = aidx - 1;
if (flag < 0) {
for (i = 0; i < n; i++)
for (int i = 0; i < n; i++)
if (mask[i] & groupbit) combine(one, carray_atom[i][aidxm1], i);
} else
one = carray_atom[flag][aidxm1];
}
} else if (flavor[m] == LOCAL) {
if (!(compute->invoked_flag & Compute::INVOKED_LOCAL)) {
compute->compute_local();
compute->invoked_flag |= Compute::INVOKED_LOCAL;
} else if (val.flavor == LOCAL) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_LOCAL)) {
val.val.c->compute_local();
val.val.c->invoked_flag |= Compute::INVOKED_LOCAL;
}
if (aidx == 0) {
double *comp_vec = compute->vector_local;
int n = compute->size_local_rows;
double *comp_vec = val.val.c->vector_local;
int n = val.val.c->size_local_rows;
if (flag < 0)
for (i = 0; i < n; i++) combine(one, comp_vec[i], i);
for (int i = 0; i < n; i++) combine(one, comp_vec[i], i);
else
one = comp_vec[flag];
} else {
double **carray_local = compute->array_local;
int n = compute->size_local_rows;
double **carray_local = val.val.c->array_local;
int n = val.val.c->size_local_rows;
int aidxm1 = aidx - 1;
if (flag < 0)
for (i = 0; i < n; i++) combine(one, carray_local[i][aidxm1], i);
for (int i = 0; i < n; i++) combine(one, carray_local[i][aidxm1], i);
else
one = carray_local[flag][aidxm1];
}
@ -539,46 +520,42 @@ double ComputeReduce::compute_one(int m, int flag)
// access fix fields, check if fix frequency is a match
} else if (which[m] == ArgInfo::FIX) {
if (update->ntimestep % modify->fix[vidx]->peratom_freq)
error->all(FLERR,
"Fix used in compute reduce not "
"computed at compatible time");
Fix *fix = modify->fix[vidx];
} else if (val.which == ArgInfo::FIX) {
if (update->ntimestep % val.val.f->peratom_freq)
error->all(FLERR, "Fix {} used in compute {} not computed at compatible time", val.id, style);
if (flavor[m] == PERATOM) {
if (val.flavor == PERATOM) {
if (aidx == 0) {
double *fix_vector = fix->vector_atom;
int n = nlocal;
double *fix_vector = val.val.f->vector_atom;
if (flag < 0) {
for (i = 0; i < n; i++)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one, fix_vector[i], i);
} else
one = fix_vector[flag];
} else {
double **fix_array = fix->array_atom;
double **fix_array = val.val.f->array_atom;
int aidxm1 = aidx - 1;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one, fix_array[i][aidxm1], i);
} else
one = fix_array[flag][aidxm1];
}
} else if (flavor[m] == LOCAL) {
} else if (val.flavor == LOCAL) {
if (aidx == 0) {
double *fix_vector = fix->vector_local;
int n = fix->size_local_rows;
double *fix_vector = val.val.f->vector_local;
int n = val.val.f->size_local_rows;
if (flag < 0)
for (i = 0; i < n; i++) combine(one, fix_vector[i], i);
for (int i = 0; i < n; i++) combine(one, fix_vector[i], i);
else
one = fix_vector[flag];
} else {
double **fix_array = fix->array_local;
int n = fix->size_local_rows;
double **fix_array = val.val.f->array_local;
int n = val.val.f->size_local_rows;
int aidxm1 = aidx - 1;
if (flag < 0)
for (i = 0; i < n; i++) combine(one, fix_array[i][aidxm1], i);
for (int i = 0; i < n; i++) combine(one, fix_array[i][aidxm1], i);
else
one = fix_array[flag][aidxm1];
}
@ -586,16 +563,16 @@ double ComputeReduce::compute_one(int m, int flag)
// evaluate atom-style variable
} else if (which[m] == ArgInfo::VARIABLE) {
} else if (val.which == ArgInfo::VARIABLE) {
if (atom->nmax > maxatom) {
maxatom = atom->nmax;
memory->destroy(varatom);
memory->create(varatom, maxatom, "reduce:varatom");
}
input->variable->compute_atom(vidx, igroup, varatom, 1, 0);
input->variable->compute_atom(val.val.v, igroup, varatom, 1, 0);
if (flag < 0) {
for (i = 0; i < nlocal; i++)
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one, varatom[i], i);
} else
one = varatom[flag];
@ -608,31 +585,28 @@ double ComputeReduce::compute_one(int m, int flag)
bigint ComputeReduce::count(int m)
{
int vidx = value2index[m];
if (which[m] == ArgInfo::X || which[m] == ArgInfo::V || which[m] == ArgInfo::F)
auto &val = values[m];
if ((val.which == ArgInfo::X) || (val.which == ArgInfo::V) || (val.which == ArgInfo::F))
return group->count(igroup);
else if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[vidx];
if (flavor[m] == PERATOM) {
else if (val.which == ArgInfo::COMPUTE) {
if (val.flavor == PERATOM) {
return group->count(igroup);
} else if (flavor[m] == LOCAL) {
bigint ncount = compute->size_local_rows;
} else if (val.flavor == LOCAL) {
bigint ncount = val.val.c->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount, &ncountall, 1, MPI_LMP_BIGINT, MPI_SUM, world);
return ncountall;
}
} else if (which[m] == ArgInfo::FIX) {
Fix *fix = modify->fix[vidx];
if (flavor[m] == PERATOM) {
} else if (val.which == ArgInfo::FIX) {
if (val.flavor == PERATOM) {
return group->count(igroup);
} else if (flavor[m] == LOCAL) {
bigint ncount = fix->size_local_rows;
} else if (val.flavor == LOCAL) {
bigint ncount = val.val.f->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount, &ncountall, 1, MPI_LMP_BIGINT, MPI_SUM, world);
return ncountall;
}
} else if (which[m] == ArgInfo::VARIABLE)
} else if (val.which == ArgInfo::VARIABLE)
return group->count(igroup);
bigint dummy = 0;

22
src/compute_reduce.h
View File

@ -37,30 +37,38 @@ class ComputeReduce : public Compute {
double memory_usage() override;
protected:
int me;
int mode, nvalues;
int *which, *argindex, *flavor, *value2index;
char **ids;
struct value_t {
int which;
int argindex;
std::string id;
int flavor;
union {
class Compute *c;
class Fix *f;
int v;
} val;
};
std::vector<value_t> values;
double *onevec;
int *replace, *indices, *owner;
int index;
char *idregion;
class Region *region;
int maxatom;
double *varatom;
struct Pair {
struct valpair {
double value;
int proc;
};
Pair pairme, pairall;
valpair pairme, pairall;
virtual double compute_one(int, int);
virtual bigint count(int);
void combine(double &, double, int);
};
} // namespace LAMMPS_NS
#endif
#endif

231
src/compute_reduce_chunk.cpp
View File

@ -30,19 +30,17 @@
using namespace LAMMPS_NS;
enum{SUM,MINN,MAXX};
enum{ SUM, MINN, MAXX };
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
which(nullptr), argindex(nullptr), value2index(nullptr), idchunk(nullptr), ids(nullptr),
vlocal(nullptr), vglobal(nullptr), alocal(nullptr), aglobal(nullptr), varatom(nullptr)
Compute(lmp, narg, arg), idchunk(nullptr), vlocal(nullptr), vglobal(nullptr),
alocal(nullptr), aglobal(nullptr), varatom(nullptr), cchunk(nullptr), ichunk(nullptr)
{
if (narg < 6) error->all(FLERR,"Illegal compute reduce/chunk command");
if (narg < 6) utils::missing_cmd_args(FLERR,"compute reduce/chunk", error);
// ID of compute chunk/atom
@ -54,7 +52,7 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[4],"sum") == 0) mode = SUM;
else if (strcmp(arg[4],"min") == 0) mode = MINN;
else if (strcmp(arg[4],"max") == 0) mode = MAXX;
else error->all(FLERR,"Illegal compute reduce/chunk command");
else error->all(FLERR,"Unknown compute reduce/chunk mode: {}", arg[4]);
int iarg = 5;
@ -67,30 +65,22 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) :
if (earg != &arg[iarg]) expand = 1;
arg = earg;
// parse values until one isn't recognized
which = new int[nargnew];
argindex = new int[nargnew];
ids = new char*[nargnew];
value2index = new int[nargnew];
for (int i=0; i < nargnew; ++i) {
which[i] = argindex[i] = value2index[i] = ArgInfo::UNKNOWN;
ids[i] = nullptr;
}
nvalues = 0;
// parse values
values.clear();
for (iarg = 0; iarg < nargnew; iarg++) {
ArgInfo argi(arg[iarg]);
which[nvalues] = argi.get_type();
argindex[nvalues] = argi.get_index1();
ids[nvalues] = argi.copy_name();
value_t val;
val.which = argi.get_type();
val.argindex = argi.get_index1();
val.id = argi.get_name();
val.val.c = nullptr;
if ((which[nvalues] == ArgInfo::UNKNOWN) || (which[nvalues] == ArgInfo::NONE)
|| (argi.get_dim() > 1))
error->all(FLERR,"Illegal compute reduce/chunk command");
if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) || (argi.get_dim() > 1))
error->all(FLERR,"Illegal compute reduce/chunk argument: {}", arg[iarg]);
nvalues++;
values.push_back(val);
}
// if wildcard expansion occurred, free earg memory from expand_args()
@ -102,58 +92,55 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) :
// error check
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 compute reduce/chunk does not exist");
if (!modify->compute[icompute]->peratom_flag)
error->all(FLERR,"Compute reduce/chunk compute does not "
"calculate per-atom values");
if (argindex[i] == 0 &&
modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,"Compute reduce/chunk compute does not "
"calculate a per-atom vector");
if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,"Compute reduce/chunk compute does not "
"calculate a per-atom array");
if (argindex[i] &&
argindex[i] > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR,
"Compute reduce/chunk compute array is accessed out-of-range");
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 compute reduce/chunk does not exist", val.id);
if (!val.val.c->peratom_flag)
error->all(FLERR,"Compute reduce/chunk compute {} does not calculate per-atom values",
val.id);
if ((val.argindex == 0) && (val.val.c->size_peratom_cols != 0))
error->all(FLERR,"Compute reduce/chunk compute {} does not calculate a per-atom vector",
val.id);
if (val.argindex && (val.val.c->size_peratom_cols == 0))
error->all(FLERR,"Compute reduce/chunk compute {} does not calculate a per-atom array",
val.id);
if (val.argindex && (val.argindex > val.val.c->size_peratom_cols))
error->all(FLERR, "Compute reduce/chunk compute array {} is accessed out-of-range", val.id);
} else if (which[i] == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(ids[i]);
if (!ifix)
error->all(FLERR,"Fix ID {} for compute reduce/chunk does not exist", ids[i]);
if (!ifix->peratom_flag)
error->all(FLERR,"Compute reduce/chunk fix {} does not calculate per-atom values", ids[i]);
if ((argindex[i] == 0) && (ifix->size_peratom_cols != 0))
error->all(FLERR,"Compute reduce/chunk fix {} does not calculate a per-atom vector", ids[i]);
if (argindex[i] && (ifix->size_peratom_cols == 0))
error->all(FLERR,"Compute reduce/chunk fix {} does not calculate a per-atom array", ids[i]);
if (argindex[i] && (argindex[i] > ifix->size_peratom_cols))
error->all(FLERR,"Compute reduce/chunk fix {} array is accessed out-of-range", ids[i]);
} 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 compute reduce/chunk does not exist", val.id);
if (!val.val.f->peratom_flag)
error->all(FLERR,"Compute reduce/chunk fix {} does not calculate per-atom values", val.id);
if ((val.argindex == 0) && (val.val.f->size_peratom_cols != 0))
error->all(FLERR,"Compute reduce/chunk fix {} does not calculate a per-atom vector", val.id);
if (val.argindex && (val.val.f->size_peratom_cols == 0))
error->all(FLERR,"Compute reduce/chunk fix {} does not calculate a per-atom array", val.id);
if (val.argindex && (val.argindex > val.val.f->size_peratom_cols))
error->all(FLERR,"Compute reduce/chunk fix {} array is accessed out-of-range", val.id);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute reduce/chunk does not exist");
if (input->variable->atomstyle(ivariable) == 0)
} 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 compute reduce/chunk does not exist", val.id);
if (input->variable->atomstyle(val.val.v) == 0)
error->all(FLERR,"Compute reduce/chunk variable is not atom-style variable");
}
}
// this compute produces either a vector or array
if (nvalues == 1) {
if (values.size() == 1) {
vector_flag = 1;
size_vector_variable = 1;
extvector = 0;
} else {
array_flag = 1;
size_array_rows_variable = 1;
size_array_cols = nvalues;
size_array_cols = values.size();
extarray = 0;
}
@ -175,13 +162,7 @@ ComputeReduceChunk::ComputeReduceChunk(LAMMPS *lmp, int narg, char **arg) :
ComputeReduceChunk::~ComputeReduceChunk()
{
delete [] idchunk;
delete [] which;
delete [] argindex;
for (int m = 0; m < nvalues; m++) delete [] ids[m];
delete [] ids;
delete [] value2index;
delete[] idchunk;
memory->destroy(vlocal);
memory->destroy(vglobal);
@ -199,24 +180,21 @@ void ComputeReduceChunk::init()
// set indices of all computes,fixes,variables
for (int m = 0; m < nvalues; m++) {
if (which[m] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute reduce/chunk does not exist");
value2index[m] = icompute;
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 compute reduce/chunk does not exist", val.id);
} else if (which[m] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute reduce/chunk does not exist");
value2index[m] = ifix;
} 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 compute reduce/chunk does not exist", val.id);
} else if (which[m] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute reduce/chunk does not exist");
value2index[m] = ivariable;
} 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 compute reduce/chunk does not exist", val.id);
}
}
}
@ -225,13 +203,10 @@ void ComputeReduceChunk::init()
void ComputeReduceChunk::init_chunk()
{
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for "
"compute reduce/chunk");
cchunk = dynamic_cast<ComputeChunkAtom *>(modify->compute[icompute]);
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Compute reduce/chunk does not use chunk/atom compute");
cchunk = dynamic_cast<ComputeChunkAtom *>(modify->get_compute_by_id(idchunk));
if (!cchunk)
error->all(FLERR,"Compute chunk/atom {} does not exist or is incorrect style for "
"compute reduce/chunk", idchunk);
}
/* ---------------------------------------------------------------------- */
@ -295,26 +270,23 @@ void ComputeReduceChunk::compute_array()
memory->destroy(alocal);
memory->destroy(aglobal);
maxchunk = nchunk;
memory->create(alocal,maxchunk,nvalues,"reduce/chunk:alocal");
memory->create(aglobal,maxchunk,nvalues,"reduce/chunk:aglobal");
memory->create(alocal,maxchunk,values.size(),"reduce/chunk:alocal");
memory->create(aglobal,maxchunk,values.size(),"reduce/chunk:aglobal");
array = aglobal;
}
// perform local reduction of all peratom values
for (int m = 0; m < nvalues; m++) compute_one(m,&alocal[0][m],nvalues);
for (std::size_t m = 0; m < values.size(); m++) compute_one(m,&alocal[0][m],values.size());
// reduce the per-chunk values across all procs
if (mode == SUM)
MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*nvalues,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*values.size(),MPI_DOUBLE,MPI_SUM,world);
else if (mode == MINN)
MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*nvalues,
MPI_DOUBLE,MPI_MIN,world);
MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*values.size(),MPI_DOUBLE,MPI_MIN,world);
else if (mode == MAXX)
MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*nvalues,
MPI_DOUBLE,MPI_MAX,world);
MPI_Allreduce(&alocal[0][0],&aglobal[0][0],nchunk*values.size(),MPI_DOUBLE,MPI_MAX,world);
}
/* ---------------------------------------------------------------------- */
@ -323,7 +295,7 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride)
{
// initialize per-chunk values in accumulation vector
for (int i = 0; i < nvalues*nchunk; i += nstride) vchunk[i] = initvalue;
for (std::size_t i = 0; i < values.size()*nchunk; i += nstride) vchunk[i] = initvalue;
// loop over my atoms
// use peratom input and chunk ID of each atom to update vector
@ -331,27 +303,22 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride)
int *mask = atom->mask;
int nlocal = atom->nlocal;
auto &val = values[m];
int index = -1;
int vidx = value2index[m];
// initialization in case it has not yet been run, e.g. when
// the compute was invoked right after it has been created
if (vidx == ArgInfo::UNKNOWN) {
init();
vidx = value2index[m];
}
if (val.val.c == nullptr) init();
if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[vidx];
if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM;
if (val.which == ArgInfo::COMPUTE) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) {
val.val.c->compute_peratom();
val.val.c->invoked_flag |= Compute::INVOKED_PERATOM;
}
if (argindex[m] == 0) {
double *vcompute = compute->vector_atom;
if (val.argindex == 0) {
double *vcompute = val.val.c->vector_atom;
for (int i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1;
@ -359,8 +326,8 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride)
combine(vchunk[index*nstride],vcompute[i]);
}
} else {
double **acompute = compute->array_atom;
int argindexm1 = argindex[m] - 1;
double **acompute = val.val.c->array_atom;
int argindexm1 = val.argindex - 1;
for (int i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1;
@ -371,13 +338,12 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride)
// access fix fields, check if fix frequency is a match
} else if (which[m] == ArgInfo::FIX) {
Fix *fix = modify->fix[vidx];
if (update->ntimestep % fix->peratom_freq)
} else if (val.which == ArgInfo::FIX) {
if (update->ntimestep % val.val.f->peratom_freq)
error->all(FLERR,"Fix used in compute reduce/chunk not computed at compatible time");
if (argindex[m] == 0) {
double *vfix = fix->vector_atom;
if (val.argindex == 0) {
double *vfix = val.val.f->vector_atom;
for (int i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1;
@ -385,8 +351,8 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride)
combine(vchunk[index*nstride],vfix[i]);
}
} else {
double **afix = fix->array_atom;
int argindexm1 = argindex[m] - 1;
double **afix = val.val.f->array_atom;
int argindexm1 = val.argindex - 1;
for (int i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1;
@ -397,14 +363,14 @@ void ComputeReduceChunk::compute_one(int m, double *vchunk, int nstride)
// evaluate atom-style variable
} else if (which[m] == ArgInfo::VARIABLE) {
} else if (val.which == ArgInfo::VARIABLE) {
if (atom->nmax > maxatom) {
memory->destroy(varatom);
maxatom = atom->nmax;
memory->create(varatom,maxatom,"reduce/chunk:varatom");
}
input->variable->compute_atom(vidx,igroup,varatom,1,0);
input->variable->compute_atom(val.val.v,igroup,varatom,1,0);
for (int i = 0; i < nlocal; i++) {
if (!(mask[i] & groupbit)) continue;
index = ichunk[i]-1;
@ -449,11 +415,8 @@ void ComputeReduceChunk::lock_enable()
void ComputeReduceChunk::lock_disable()
{
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = dynamic_cast<ComputeChunkAtom *>(modify->compute[icompute]);
cchunk->lockcount--;
}
cchunk = dynamic_cast<ComputeChunkAtom *>(modify->get_compute_by_id(idchunk));
if (cchunk) cchunk->lockcount--;
}
/* ----------------------------------------------------------------------
@ -491,7 +454,7 @@ void ComputeReduceChunk::unlock(Fix *fixptr)
double ComputeReduceChunk::memory_usage()
{
double bytes = (bigint) maxatom * sizeof(double);
if (nvalues == 1) bytes += (double) maxchunk * 2 * sizeof(double);
else bytes += (double) maxchunk * nvalues * 2 * sizeof(double);
if (values.size() == 1) bytes += (double) maxchunk * 2 * sizeof(double);
else bytes += (double) maxchunk * values.size() * 2 * sizeof(double);
return bytes;
}

21
src/compute_reduce_chunk.h
View File

@ -41,12 +41,21 @@ class ComputeReduceChunk : public Compute {
double memory_usage() override;
private:
int mode, nvalues;
int *which, *argindex, *value2index;
char *idchunk;
char **ids;
struct value_t {
int which;
int argindex;
std::string id;
union {
class Compute *c;
class Fix *f;
int v;
} val;
};
std::vector<value_t> values;
int nchunk;
char *idchunk;
int mode, nchunk;
int maxchunk, maxatom;
double initvalue;
double *vlocal, *vglobal;
@ -60,8 +69,6 @@ class ComputeReduceChunk : public Compute {
void compute_one(int, double *, int);
void combine(double &, double);
};
} // namespace LAMMPS_NS
#endif
#endif

151
src/compute_reduce_region.cpp
View File

@ -55,62 +55,58 @@ double ComputeReduceRegion::compute_one(int m, int flag)
// only include atoms in group
index = -1;
auto &val = values[m];
// initialization in case it has not yet been run, e.g. when
// the compute was invoked right after it has been created
if ((val.which == ArgInfo::COMPUTE) || (val.which == ArgInfo::FIX)) {
if (val.val.c == nullptr) init();
}
int aidx = val.argindex;
double **x = atom->x;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int n = value2index[m];
// initialization in case it has not yet been run,
// e.g. when invoked
if (n == ArgInfo::UNKNOWN) {
init();
n = value2index[m];
}
int j = argindex[m];
double one = 0.0;
if (mode == MINN) one = BIG;
if (mode == MAXX) one = -BIG;
if (which[m] == ArgInfo::X) {
if (val.which == ArgInfo::X) {
if (flag < 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
combine(one, x[i][j], i);
combine(one, x[i][aidx], i);
} else
one = x[flag][j];
} else if (which[m] == ArgInfo::V) {
one = x[flag][aidx];
} else if (val.which == ArgInfo::V) {
double **v = atom->v;
if (flag < 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
combine(one, v[i][j], i);
combine(one, v[i][aidx], i);
} else
one = v[flag][j];
} else if (which[m] == ArgInfo::F) {
one = v[flag][aidx];
} else if (val.which == ArgInfo::F) {
double **f = atom->f;
if (flag < 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
combine(one, f[i][j], i);
combine(one, f[i][aidx], i);
} else
one = f[flag][j];
one = f[flag][aidx];
// invoke compute if not previously invoked
} else if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[n];
if (flavor[m] == PERATOM) {
if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM;
} else if (val.which == ArgInfo::COMPUTE) {
if (val.flavor == 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) {
double *compute_vector = compute->vector_atom;
if (aidx == 0) {
double *compute_vector = val.val.c->vector_atom;
if (flag < 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
@ -118,48 +114,48 @@ double ComputeReduceRegion::compute_one(int m, int flag)
} else
one = compute_vector[flag];
} else {
double **compute_array = compute->array_atom;
int jm1 = j - 1;
double **compute_array = val.val.c->array_atom;
int aidxm1 = aidx - 1;
if (flag < 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
combine(one, compute_array[i][jm1], i);
combine(one, compute_array[i][aidxm1], i);
} else
one = compute_array[flag][jm1];
one = compute_array[flag][aidxm1];
}
} else if (flavor[m] == LOCAL) {
if (!(compute->invoked_flag & Compute::INVOKED_LOCAL)) {
compute->compute_local();
compute->invoked_flag |= Compute::INVOKED_LOCAL;
} else if (val.flavor == 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) {
double *compute_vector = compute->vector_local;
if (aidx == 0) {
double *compute_vector = val.val.c->vector_local;
if (flag < 0)
for (int i = 0; i < compute->size_local_rows; i++) combine(one, compute_vector[i], i);
for (int i = 0; i < val.val.c->size_local_rows; i++) combine(one, compute_vector[i], i);
else
one = compute_vector[flag];
} else {
double **compute_array = compute->array_local;
int jm1 = j - 1;
double **compute_array = val.val.c->array_local;
int aidxm1 = aidx - 1;
if (flag < 0)
for (int i = 0; i < compute->size_local_rows; i++) combine(one, compute_array[i][jm1], i);
for (int i = 0; i < val.val.c->size_local_rows; i++)
combine(one, compute_array[i][aidxm1], i);
else
one = compute_array[flag][jm1];
one = compute_array[flag][aidxm1];
}
}
// check if fix frequency is a match
} else if (which[m] == ArgInfo::FIX) {
if (update->ntimestep % modify->fix[n]->peratom_freq)
error->all(FLERR, "Fix used in compute reduce not computed at compatible time");
Fix *fix = modify->fix[n];
} else if (val.which == ArgInfo::FIX) {
if (update->ntimestep % val.val.f->peratom_freq)
error->all(FLERR, "Fix {} used in compute {} not computed at compatible time", val.id, style);
if (flavor[m] == PERATOM) {
if (j == 0) {
double *fix_vector = fix->vector_atom;
if (val.flavor == PERATOM) {
if (aidx == 0) {
double *fix_vector = val.val.f->vector_atom;
if (flag < 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
@ -167,43 +163,44 @@ double ComputeReduceRegion::compute_one(int m, int flag)
} else
one = fix_vector[flag];
} else {
double **fix_array = fix->array_atom;
int jm1 = j - 1;
double **fix_array = val.val.f->array_atom;
int aidxm1 = aidx - 1;
if (flag < 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
combine(one, fix_array[i][jm1], i);
combine(one, fix_array[i][aidxm1], i);
} else
one = fix_array[flag][jm1];
one = fix_array[flag][aidxm1];
}
} else if (flavor[m] == LOCAL) {
if (j == 0) {
double *fix_vector = fix->vector_local;
} else if (val.flavor == LOCAL) {
if (aidx == 0) {
double *fix_vector = val.val.f->vector_local;
if (flag < 0)
for (int i = 0; i < fix->size_local_rows; i++) combine(one, fix_vector[i], i);
for (int i = 0; i < val.val.f->size_local_rows; i++) combine(one, fix_vector[i], i);
else
one = fix_vector[flag];
} else {
double **fix_array = fix->array_local;
int jm1 = j - 1;
double **fix_array = val.val.f->array_local;
int aidxm1 = aidx - 1;
if (flag < 0)
for (int i = 0; i < fix->size_local_rows; i++) combine(one, fix_array[i][jm1], i);
for (int i = 0; i < val.val.f->size_local_rows; i++)
combine(one, fix_array[i][aidxm1], i);
else
one = fix_array[flag][jm1];
one = fix_array[flag][aidxm1];
}
}
// evaluate atom-style variable
} else if (which[m] == ArgInfo::VARIABLE) {
} else if (val.which == ArgInfo::VARIABLE) {
if (atom->nmax > maxatom) {
maxatom = atom->nmax;
memory->destroy(varatom);
memory->create(varatom, maxatom, "reduce/region:varatom");
}
input->variable->compute_atom(n, igroup, varatom, 1, 0);
input->variable->compute_atom(val.val.v, igroup, varatom, 1, 0);
if (flag < 0) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2]))
@ -219,31 +216,29 @@ double ComputeReduceRegion::compute_one(int m, int flag)
bigint ComputeReduceRegion::count(int m)
{
int n = value2index[m];
auto &val = values[m];
if (which[m] == ArgInfo::X || which[m] == ArgInfo::V || which[m] == ArgInfo::F)
if (val.which == ArgInfo::X || val.which == ArgInfo::V || val.which == ArgInfo::F)
return group->count(igroup, region);
else if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[n];
if (flavor[m] == PERATOM) {
else if (val.which == ArgInfo::COMPUTE) {
if (val.flavor == PERATOM) {
return group->count(igroup, region);
} else if (flavor[m] == LOCAL) {
bigint ncount = compute->size_local_rows;
} else if (val.flavor == LOCAL) {
bigint ncount = val.val.c->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount, &ncountall, 1, MPI_DOUBLE, MPI_SUM, world);
return ncountall;
}
} else if (which[m] == ArgInfo::FIX) {
Fix *fix = modify->fix[n];
if (flavor[m] == PERATOM) {
} else if (val.which == ArgInfo::FIX) {
if (val.flavor == PERATOM) {
return group->count(igroup, region);
} else if (flavor[m] == LOCAL) {
bigint ncount = fix->size_local_rows;
} else if (val.flavor == LOCAL) {
bigint ncount = val.val.f->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount, &ncountall, 1, MPI_DOUBLE, MPI_SUM, world);
return ncountall;
}
} else if (which[m] == ArgInfo::VARIABLE)
} else if (val.which == ArgInfo::VARIABLE)
return group->count(igroup, region);
bigint dummy = 0;

336
src/compute_slice.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -27,96 +26,84 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
nvalues(0), which(nullptr), argindex(nullptr), value2index(nullptr), ids(nullptr)
ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal compute slice command");
if (narg < 7) utils::missing_cmd_args(FLERR, "compute slice", error);
MPI_Comm_rank(world,&me);
nstart = utils::inumeric(FLERR, arg[3], false, lmp);
nstop = utils::inumeric(FLERR, arg[4], false, lmp);
nskip = utils::inumeric(FLERR, arg[5], false, lmp);
nstart = utils::inumeric(FLERR,arg[3],false,lmp);
nstop = utils::inumeric(FLERR,arg[4],false,lmp);
nskip = utils::inumeric(FLERR,arg[5],false,lmp);
if (nstart < 1) error->all(FLERR, "Invalid compute slice nstart value {} < 1", nstart);
if (nstop < nstart) error->all(FLERR, "Invalid compute slice nstop value {} < {}", nstop, nstart);
if (nskip < 1) error->all(FLERR, "Invalid compute slice nskip value < 1: {}", nskip);
if (nstart < 1 || nstop < nstart || nskip < 1)
error->all(FLERR,"Illegal compute slice command");
// parse remaining values until one isn't recognized
which = new int[narg-6];
argindex = new int[narg-6];
ids = new char*[narg-6];
value2index = new int[narg-6];
nvalues = 0;
// parse values
values.clear();
for (int iarg = 6; iarg < narg; iarg++) {
ArgInfo argi(arg[iarg]);
which[nvalues] = argi.get_type();
argindex[nvalues] = argi.get_index1();
ids[nvalues] = argi.copy_name();
value_t val;
val.which = argi.get_type();
val.argindex = argi.get_index1();
val.id = argi.get_name();
val.val.c = nullptr;
if ((which[nvalues] == ArgInfo::UNKNOWN) || (which[nvalues] == ArgInfo::NONE)
|| (argi.get_dim() > 1))
error->all(FLERR,"Illegal compute slice command");
if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) || (argi.get_dim() > 1))
error->all(FLERR, "Illegal compute slice argument: {}", arg[iarg]);
nvalues++;
values.push_back(val);
}
// setup and error check
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 compute slice does not exist");
if (modify->compute[icompute]->vector_flag) {
if (argindex[i])
error->all(FLERR,"Compute slice compute does not "
"calculate a global array");
if (nstop > modify->compute[icompute]->size_vector)
error->all(FLERR,"Compute slice compute vector is "
"accessed out-of-range");
} else if (modify->compute[icompute]->array_flag) {
if (argindex[i] == 0)
error->all(FLERR,"Compute slice compute does not "
"calculate a global vector");
if (argindex[i] > modify->compute[icompute]->size_array_cols)
error->all(FLERR,"Compute slice compute array is "
"accessed out-of-range");
if (nstop > modify->compute[icompute]->size_array_rows)
error->all(FLERR,"Compute slice compute array is "
"accessed out-of-range");
} else error->all(FLERR,"Compute slice compute does not calculate "
"global vector or array");
} else if (which[i] == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(ids[i]);
if (!ifix)
error->all(FLERR,"Fix ID {} for compute slice does not exist", ids[i]);
if (ifix->vector_flag) {
if (argindex[i])
error->all(FLERR,"Compute slice fix {} does not calculate a global array", ids[i]);
if (nstop > ifix->size_vector)
error->all(FLERR,"Compute slice fix {} vector is accessed out-of-range", ids[i]);
} else if (ifix->array_flag) {
if (argindex[i] == 0)
error->all(FLERR,"Compute slice fix {} does not calculate a global vector", ids[i]);
if (argindex[i] > ifix->size_array_cols)
error->all(FLERR,"Compute slice fix {} array is accessed out-of-range", ids[i]);
if (nstop > ifix->size_array_rows)
error->all(FLERR,"Compute slice fix {} array is accessed out-of-range", ids[i]);
} else error->all(FLERR,"Compute slice fix {} does not calculate global vector or array", ids[i]);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute slice does not exist");
if (argindex[i] == 0 && input->variable->vectorstyle(ivariable) == 0)
error->all(FLERR,"Compute slice variable is not vector-style variable");
if (argindex[i])
error->all(FLERR,"Compute slice vector variable cannot be indexed");
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 compute slice does not exist", val.id);
if (val.val.c->vector_flag) {
if (val.argindex)
error->all(FLERR, "Compute slice compute {} does not calculate a global array", val.id);
if (nstop > val.val.c->size_vector)
error->all(FLERR, "Compute slice compute {} vector is accessed out-of-range", val.id);
} else if (val.val.c->array_flag) {
if (val.argindex == 0)
error->all(FLERR, "Compute slice compute {} does not calculate a global vector", val.id);
if (val.argindex > val.val.c->size_array_cols)
error->all(FLERR, "Compute slice compute {} array is accessed out-of-range", val.id);
if (nstop > val.val.c->size_array_rows)
error->all(FLERR, "Compute slice compute {} array is accessed out-of-range", val.id);
} else {
error->all(FLERR, "Compute slice compute {} does not calculate global vector or array",
val.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 compute slice does not exist", val.id);
if (val.val.f->vector_flag) {
if (val.argindex)
error->all(FLERR, "Compute slice fix {} does not calculate a global array", val.id);
if (nstop > val.val.f->size_vector)
error->all(FLERR, "Compute slice fix {} vector is accessed out-of-range", val.id);
} else if (val.val.f->array_flag) {
if (val.argindex == 0)
error->all(FLERR, "Compute slice fix {} does not calculate a global vector", val.id);
if (val.argindex > val.val.f->size_array_cols)
error->all(FLERR, "Compute slice fix {} array is accessed out-of-range", val.id);
if (nstop > val.val.f->size_array_rows)
error->all(FLERR, "Compute slice fix {} array is accessed out-of-range", val.id);
} else {
error->all(FLERR, "Compute slice fix {} does not calculate global vector or array", val.id);
}
} 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 compute slice does not exist", val.id);
if (val.argindex == 0 && input->variable->vectorstyle(val.val.v) == 0)
error->all(FLERR, "Compute slice variable {} is not vector-style variable", val.id);
if (val.argindex)
error->all(FLERR, "Compute slice vector variable {} cannot be indexed", val.id);
}
}
@ -124,68 +111,65 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) :
// for vector, set intensive/extensive to mirror input values
// for array, set intensive if all input values are intensive, else extensive
if (nvalues == 1) {
if (values.size() == 1) {
auto &val = values[0];
vector_flag = 1;
size_vector = (nstop-nstart) / nskip;
memory->create(vector,size_vector,"slice:vector");
size_vector = (nstop - nstart) / nskip;
memory->create(vector, size_vector, "slice:vector");
if (which[0] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[0]);
if (argindex[0] == 0) {
extvector = modify->compute[icompute]->extvector;
if (modify->compute[icompute]->extvector == -1) {
if (val.which == ArgInfo::COMPUTE) {
if (val.argindex == 0) {
extvector = val.val.c->extvector;
if (val.val.c->extvector == -1) {
extlist = new int[size_vector];
int j = 0;
for (int i = nstart; i < nstop; i += nskip)
extlist[j++] = modify->compute[icompute]->extlist[i-1];
for (int i = nstart; i < nstop; i += nskip) extlist[j++] = val.val.c->extlist[i - 1];
}
} else extvector = modify->compute[icompute]->extarray;
} else if (which[0] == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(ids[0]);
if (argindex[0] == 0) {
extvector = ifix->extvector;
if (ifix->extvector == -1) {
} else
extvector = val.val.c->extarray;
} else if (val.which == ArgInfo::FIX) {
if (val.argindex == 0) {
extvector = val.val.f->extvector;
if (val.val.f->extvector == -1) {
extlist = new int[size_vector];
int j = 0;
for (int i = nstart; i < nstop; i += nskip)
extlist[j++] = ifix->extlist[i-1];
for (int i = nstart; i < nstop; i += nskip) extlist[j++] = val.val.f->extlist[i - 1];
}
} else extvector = ifix->extarray;
} else if (which[0] == ArgInfo::VARIABLE) {
} else
extvector = val.val.f->extarray;
} else if (val.which == ArgInfo::VARIABLE) {
extvector = 0;
}
} else {
array_flag = 1;
size_array_rows = (nstop-nstart) / nskip;
size_array_cols = nvalues;
memory->create(array,size_array_rows,size_array_cols,"slice:array");
size_array_rows = (nstop - nstart) / nskip;
size_array_cols = values.size();
memory->create(array, size_array_rows, size_array_cols, "slice:array");
extarray = 0;
for (int i = 0; i < nvalues; i++) {
if (which[i] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[i]);
if (argindex[i] == 0) {
if (modify->compute[icompute]->extvector == 1) extarray = 1;
if (modify->compute[icompute]->extvector == -1) {
for (int j = 0; j < modify->compute[icompute]->size_vector; j++)
if (modify->compute[icompute]->extlist[j]) extarray = 1;
for (auto &val : values) {
if (val.which == ArgInfo::COMPUTE) {
if (val.argindex == 0) {
if (val.val.c->extvector == 1) extarray = 1;
if (val.val.c->extvector == -1) {
for (int j = 0; j < val.val.c->size_vector; j++)
if (val.val.c->extlist[j]) extarray = 1;
}
} else {
if (modify->compute[icompute]->extarray) extarray = 1;
if (val.val.c->extarray) extarray = 1;
}
} else if (which[i] == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(ids[i]);
if (argindex[i] == 0) {
if (ifix->extvector == 1) extarray = 1;
if (ifix->extvector == -1) {
for (int j = 0; j < ifix->size_vector; j++)
if (ifix->extlist[j]) extarray = 1;
} else if (val.which == ArgInfo::FIX) {
if (val.argindex == 0) {
if (val.val.f->extvector == 1) extarray = 1;
if (val.val.f->extvector == -1) {
for (int j = 0; j < val.val.f->size_vector; j++)
if (val.val.f->extlist[j]) extarray = 1;
}
} else {
if (ifix->extarray) extarray = 1;
if (val.val.f->extarray) extarray = 1;
}
} else if (which[i] == ArgInfo::VARIABLE) {
} else if (val.which == ArgInfo::VARIABLE) {
// variable is always intensive, does not change extarray
}
}
@ -196,12 +180,7 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) :
ComputeSlice::~ComputeSlice()
{
delete [] which;
delete [] argindex;
for (int m = 0; m < nvalues; m++) delete [] ids[m];
delete [] ids;
delete [] value2index;
delete [] extlist;
delete[] extlist;
memory->destroy(vector);
memory->destroy(array);
@ -213,22 +192,17 @@ void ComputeSlice::init()
{
// set indices and check validity of all computes,fixes
for (int m = 0; m < nvalues; m++) {
if (which[m] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute slice does not exist");
value2index[m] = icompute;
} else if (which[m] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute slice does not exist");
value2index[m] = ifix;
} else if (which[m] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute slice does not exist");
value2index[m] = ivariable;
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 compute slice does not exist", val.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 compute slice does not exist", val.id);
} 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 compute slice does not exist", val.id);
}
}
}
@ -239,7 +213,7 @@ void ComputeSlice::compute_vector()
{
invoked_vector = update->ntimestep;
extract_one(0,vector,1);
extract_one(0, vector, 1);
}
/* ---------------------------------------------------------------------- */
@ -248,8 +222,7 @@ void ComputeSlice::compute_array()
{
invoked_array = update->ntimestep;
for (int m = 0; m < nvalues; m++)
extract_one(0,&array[m][0],nvalues);
for (std::size_t m = 0; m < values.size(); m++) extract_one(0, &array[m][0], values.size());
}
/* ----------------------------------------------------------------------
@ -259,72 +232,67 @@ void ComputeSlice::compute_array()
void ComputeSlice::extract_one(int m, double *vec, int stride)
{
int i,j;
auto &val = values[m];
// invoke the appropriate compute if needed
if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[value2index[m]];
if (argindex[m] == 0) {
if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
compute->compute_vector();
compute->invoked_flag |= Compute::INVOKED_VECTOR;
if (val.which == ArgInfo::COMPUTE) {
if (val.argindex == 0) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_VECTOR)) {
val.val.c->compute_vector();
val.val.c->invoked_flag |= Compute::INVOKED_VECTOR;
}
double *cvector = compute->vector;
j = 0;
for (i = nstart; i < nstop; i += nskip) {
vec[j] = cvector[i-1];
double *cvector = val.val.c->vector;
int j = 0;
for (int i = nstart; i < nstop; i += nskip) {
vec[j] = cvector[i - 1];
j += 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;
}
double **carray = compute->array;
int icol = argindex[m]-1;
j = 0;
for (i = nstart; i < nstop; i += nskip) {
vec[j] = carray[i-1][icol];
double **carray = val.val.c->array;
int icol = val.argindex - 1;
int j = 0;
for (int i = nstart; i < nstop; i += nskip) {
vec[j] = carray[i - 1][icol];
j += stride;
}
}
// access fix fields, check if fix frequency is a match
// access fix fields, check if fix frequency is a match
} else if (which[m] == ArgInfo::FIX) {
if (update->ntimestep % modify->fix[value2index[m]]->global_freq)
error->all(FLERR,"Fix used in compute slice not "
"computed at compatible time");
Fix *fix = modify->fix[value2index[m]];
} else if (val.which == ArgInfo::FIX) {
if (update->ntimestep % val.val.f->global_freq)
error->all(FLERR, "Fix {} used in compute slice not computed at compatible time", val.id);
if (argindex[m] == 0) {
j = 0;
for (i = nstart; i < nstop; i += nskip) {
vec[j] = fix->compute_vector(i-1);
if (val.argindex == 0) {
int j = 0;
for (int i = nstart; i < nstop; i += nskip) {
vec[j] = val.val.f->compute_vector(i - 1);
j += stride;
}
} else {
int icol = argindex[m]-1;
j = 0;
for (i = nstart; i < nstop; i += nskip) {
vec[j] = fix->compute_array(i-1,icol);
int icol = val.argindex - 1;
int j = 0;
for (int i = nstart; i < nstop; i += nskip) {
vec[j] = val.val.f->compute_array(i - 1, icol);
j += stride;
}
}
// invoke vector-style variable
} else if (which[m] == ArgInfo::VARIABLE) {
} else if (val.which == ArgInfo::VARIABLE) {
double *varvec;
int nvec = input->variable->compute_vector(value2index[m],&varvec);
if (nvec < nstop)
error->all(FLERR,"Compute slice variable is not long enough");
j = 0;
for (i = nstart; i < nstop; i += nskip) {
vec[j] = varvec[i-1];
int nvec = input->variable->compute_vector(val.val.v, &varvec);
if (nvec < nstop) error->all(FLERR, "Compute slice variable {} is not long enough", val.id);
int j = 0;
for (int i = nstart; i < nstop; i += nskip) {
vec[j] = varvec[i - 1];
j += stride;
}
}

16
src/compute_slice.h
View File

@ -33,10 +33,18 @@ class ComputeSlice : public Compute {
void compute_array() override;
private:
int me;
int nstart, nstop, nskip, nvalues;
int *which, *argindex, *value2index;
char **ids;
struct value_t {
int which;
int argindex;
std::string id;
union {
class Compute *c;
class Fix *f;
int v;
} val;
};
std::vector<value_t> values;
int nstart, nstop, nskip;
void extract_one(int, double *, int);
};

279
src/fix_ave_atom.cpp
View File

@ -31,149 +31,135 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
nvalues(0), which(nullptr), argindex(nullptr), value2index(nullptr),
ids(nullptr), array(nullptr)
FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), array(nullptr)
{
if (narg < 7) error->all(FLERR,"Illegal fix ave/atom command");
if (narg < 7) utils::missing_cmd_args(FLERR, "fix ave/atom", error);
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
nrepeat = utils::inumeric(FLERR,arg[4],false,lmp);
peratom_freq = utils::inumeric(FLERR,arg[5],false,lmp);
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
nrepeat = utils::inumeric(FLERR, arg[4], false, lmp);
peratom_freq = utils::inumeric(FLERR, arg[5], false, lmp);
time_depend = 1;
nvalues = narg - 6;
// expand args if any have wildcard character "*"
// this can reset nvalues
int expand = 0;
char **earg;
nvalues = utils::expand_args(FLERR,nvalues,&arg[6],1,earg,lmp);
int nvalues = utils::expand_args(FLERR, narg - 6, &arg[6], 1, earg, lmp);
if (earg != &arg[6]) expand = 1;
arg = earg;
// parse values
which = new int[nvalues];
argindex = new int[nvalues];
ids = new char*[nvalues];
value2index = new int[nvalues];
values.clear();
for (int i = 0; i < nvalues; i++) {
ids[i] = nullptr;
if (strcmp(arg[i],"x") == 0) {
which[i] = ArgInfo::X;
argindex[i] = 0;
} else if (strcmp(arg[i],"y") == 0) {
which[i] = ArgInfo::X;
argindex[i] = 1;
} else if (strcmp(arg[i],"z") == 0) {
which[i] = ArgInfo::X;
argindex[i] = 2;
value_t val;
val.id = "";
val.val.c = nullptr;
} else if (strcmp(arg[i],"vx") == 0) {
which[i] = ArgInfo::V;
argindex[i] = 0;
} else if (strcmp(arg[i],"vy") == 0) {
which[i] = ArgInfo::V;
argindex[i] = 1;
} else if (strcmp(arg[i],"vz") == 0) {
which[i] = ArgInfo::V;
argindex[i] = 2;
if (strcmp(arg[i], "x") == 0) {
val.which = ArgInfo::X;
val.argindex = 0;
} else if (strcmp(arg[i], "y") == 0) {
val.which = ArgInfo::X;
val.argindex = 1;
} else if (strcmp(arg[i], "z") == 0) {
val.which = ArgInfo::X;
val.argindex = 2;
} else if (strcmp(arg[i],"fx") == 0) {
which[i] = ArgInfo::F;
argindex[i] = 0;
} else if (strcmp(arg[i],"fy") == 0) {
which[i] = ArgInfo::F;
argindex[i] = 1;
} else if (strcmp(arg[i],"fz") == 0) {
which[i] = ArgInfo::F;
argindex[i] = 2;
} else if (strcmp(arg[i], "vx") == 0) {
val.which = ArgInfo::V;
val.argindex = 0;
} else if (strcmp(arg[i], "vy") == 0) {
val.which = ArgInfo::V;
val.argindex = 1;
} else if (strcmp(arg[i], "vz") == 0) {
val.which = ArgInfo::V;
val.argindex = 2;
} else if (strcmp(arg[i], "fx") == 0) {
val.which = ArgInfo::F;
val.argindex = 0;
} else if (strcmp(arg[i], "fy") == 0) {
val.which = ArgInfo::F;
val.argindex = 1;
} else if (strcmp(arg[i], "fz") == 0) {
val.which = ArgInfo::F;
val.argindex = 2;
} else {
ArgInfo argi(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();
if ((which[i] == ArgInfo::UNKNOWN) || (which[i] == ArgInfo::NONE)
|| (argi.get_dim() > 1))
error->all(FLERR,"Illegal fix ave/atom command");
if ((val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) || (argi.get_dim() > 1))
error->all(FLERR, "Invalid fix ave/atom argument: {}", arg[i]);
}
values.push_back(val);
}
// if wildcard expansion occurred, free earg memory from exapnd_args()
if (expand) {
for (int i = 0; i < nvalues; i++) delete [] earg[i];
for (int i = 0; i < nvalues; i++) delete[] earg[i];
memory->sfree(earg);
}
// setup and error check
// for fix inputs, check that fix frequency is acceptable
if (nevery <= 0 || nrepeat <= 0 || peratom_freq <= 0)
error->all(FLERR,"Illegal fix ave/atom command");
if (nevery <= 0) error->all(FLERR,"Illegal fix ave/atom nevery value: {}", nevery);
if (nrepeat <= 0) error->all(FLERR,"Illegal fix ave/atom nrepeat value: {}", nrepeat);
if (peratom_freq <= 0) error->all(FLERR,"Illegal fix ave/atom nfreq value: {}", peratom_freq);
if (peratom_freq % nevery || nrepeat*nevery > peratom_freq)
error->all(FLERR,"Illegal fix ave/atom command");
error->all(FLERR,"Inconsistent fix ave/atom nevery/nrepeat/nfreq values");
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/atom does not exist");
if (modify->compute[icompute]->peratom_flag == 0)
error->all(FLERR,
"Fix ave/atom compute does not calculate per-atom values");
if (argindex[i] == 0 &&
modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,"Fix ave/atom compute does not "
"calculate a per-atom vector");
if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,"Fix ave/atom compute does not "
"calculate a per-atom array");
if (argindex[i] &&
argindex[i] > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR,"Fix ave/atom compute array is accessed out-of-range");
for (auto &val : values) {
} else if (which[i] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/atom does not exist");
if (modify->fix[ifix]->peratom_flag == 0)
error->all(FLERR,"Fix ave/atom fix does not calculate per-atom values");
if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
error->all(FLERR,
"Fix ave/atom fix does not calculate a per-atom vector");
if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
error->all(FLERR,
"Fix ave/atom fix does not calculate a per-atom array");
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
error->all(FLERR,"Fix ave/atom fix array is accessed out-of-range");
if (nevery % modify->fix[ifix]->peratom_freq)
error->all(FLERR,
"Fix for fix ave/atom not computed at compatible time");
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 fix ave/atom does not exist", val.id);
if (val.val.c->peratom_flag == 0)
error->all(FLERR, "Fix ave/atom compute {} does not calculate per-atom values", val.id);
if (val.argindex == 0 && val.val.c->size_peratom_cols != 0)
error->all(FLERR,"Fix ave/atom compute {} does not calculate a per-atom vector", val.id);
if (val.argindex && val.val.c->size_peratom_cols == 0)
error->all(FLERR,"Fix ave/atom compute {} does not calculate a per-atom array", val.id);
if (val.argindex && val.argindex > val.val.c->size_peratom_cols)
error->all(FLERR,"Fix ave/atom compute {} array is accessed out-of-range", val.id);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/atom does not exist");
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Fix ave/atom variable is not atom-style variable");
} 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 fix ave/atom does not exist", val.id);
if (val.val.f->peratom_flag == 0)
error->all(FLERR, "Fix ave/atom fix {} does not calculate per-atom values", val.id);
if (val.argindex == 0 && val.val.f->size_peratom_cols != 0)
error->all(FLERR, "Fix ave/atom fix {} does not calculate a per-atom vector", val.id);
if (val.argindex && val.val.f->size_peratom_cols == 0)
error->all(FLERR, "Fix ave/atom fix {} does not calculate a per-atom array", val.id);
if (val.argindex && val.argindex > val.val.f->size_peratom_cols)
error->all(FLERR,"Fix ave/atom fix {} array is accessed out-of-range", val.id);
if (nevery % val.val.f->peratom_freq)
error->all(FLERR, "Fix {} for fix ave/atom not computed at compatible time", val.id);
} 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 fix ave/atom does not exist", val.id);
if (input->variable->atomstyle(val.val.v) == 0)
error->all(FLERR,"Fix ave/atom variable {} is not atom-style variable", val.id);
}
}
// this fix produces either a per-atom vector or array
peratom_flag = 1;
if (nvalues == 1) size_peratom_cols = 0;
else size_peratom_cols = nvalues;
if (values.size() == 1) size_peratom_cols = 0;
else size_peratom_cols = values.size();
// perform initial allocation of atom-based array
// register with Atom class
@ -186,7 +172,7 @@ FixAveAtom::FixAveAtom(LAMMPS *lmp, int narg, char **arg) :
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
for (int m = 0; m < nvalues; m++)
for (std::size_t m = 0; m < values.size(); m++)
array[i][m] = 0.0;
// nvalid = next step on which end_of_step does something
@ -207,13 +193,6 @@ FixAveAtom::~FixAveAtom()
// unregister callback to this fix from Atom class
atom->delete_callback(id,Atom::GROW);
delete [] which;
delete [] argindex;
for (int m = 0; m < nvalues; m++) delete [] ids[m];
delete [] ids;
delete [] value2index;
memory->destroy(array);
}
@ -232,26 +211,20 @@ void FixAveAtom::init()
{
// set indices and check validity of all computes,fixes,variables
for (int m = 0; m < nvalues; m++) {
if (which[m] == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID for fix ave/atom does not exist");
value2index[m] = icompute;
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 fix ave/atom does not exist", val.id);
} else if (which[m] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/atom does not exist");
value2index[m] = ifix;
} 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 fix ave/atom does not exist", val.id);
} else if (which[m] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/atom does not exist");
value2index[m] = ivariable;
} else value2index[m] = -1;
} 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 fix ave/atom does not exist", val.id);
}
}
// need to reset nvalid if nvalid < ntimestep b/c minimize was performed
@ -276,8 +249,6 @@ void FixAveAtom::setup(int /*vflag*/)
void FixAveAtom::end_of_step()
{
int i,j,m,n;
// skip if not step which requires doing something
bigint ntimestep = update->ntimestep;
@ -289,8 +260,8 @@ void FixAveAtom::end_of_step()
int nlocal = atom->nlocal;
if (irepeat == 0)
for (i = 0; i < nlocal; i++)
for (m = 0; m < nvalues; m++)
for (int i = 0; i < nlocal; i++)
for (std::size_t m = 0; m < values.size(); m++)
array[i][m] = 0.0;
// accumulate results of attributes,computes,fixes,variables to local copy
@ -300,55 +271,54 @@ void FixAveAtom::end_of_step()
int *mask = atom->mask;
for (m = 0; m < nvalues; m++) {
n = value2index[m];
j = argindex[m];
int i, j, m = 0;
for (auto &val : values) {
j = val.argindex;
if (which[m] == ArgInfo::X) {
if (val.which == ArgInfo::X) {
double **x = atom->x;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) array[i][m] += x[i][j];
} else if (which[m] == ArgInfo::V) {
} else if (val.which == ArgInfo::V) {
double **v = atom->v;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) array[i][m] += v[i][j];
} else if (which[m] == ArgInfo::F) {
} else if (val.which == ArgInfo::F) {
double **f = atom->f;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) array[i][m] += f[i][j];
// invoke compute if not previously invoked
} else if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[n];
if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM;
} else if (val.which == ArgInfo::COMPUTE) {
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) {
double *compute_vector = compute->vector_atom;
double *compute_vector = val.val.c->vector_atom;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) array[i][m] += compute_vector[i];
} else {
int jm1 = j - 1;
double **compute_array = compute->array_atom;
double **compute_array = val.val.c->array_atom;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) array[i][m] += compute_array[i][jm1];
}
// access fix fields, guaranteed to be ready
} else if (which[m] == ArgInfo::FIX) {
} else if (val.which == ArgInfo::FIX) {
if (j == 0) {
double *fix_vector = modify->fix[n]->vector_atom;
double *fix_vector = val.val.f->vector_atom;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) array[i][m] += fix_vector[i];
} else {
int jm1 = j - 1;
double **fix_array = modify->fix[n]->array_atom;
double **fix_array = val.val.f->array_atom;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) array[i][m] += fix_array[i][jm1];
}
@ -356,10 +326,11 @@ void FixAveAtom::end_of_step()
// evaluate atom-style variable
// final argument = 1 sums result to array
} else if (which[m] == ArgInfo::VARIABLE) {
if (array) input->variable->compute_atom(n,igroup,&array[0][m],nvalues,1);
else input->variable->compute_atom(n,igroup,nullptr,nvalues,1);
} else if (val.which == ArgInfo::VARIABLE) {
if (array) input->variable->compute_atom(val.val.v,igroup,&array[0][m],values.size(),1);
else input->variable->compute_atom(val.val.v,igroup,nullptr,values.size(),1);
}
++m;
}
// done if irepeat < nrepeat
@ -382,7 +353,7 @@ void FixAveAtom::end_of_step()
double repeat = nrepeat;
for (i = 0; i < nlocal; i++)
for (m = 0; m < nvalues; m++)
for (m = 0; m < (int)values.size(); m++)
array[i][m] /= repeat;
}
@ -393,7 +364,7 @@ void FixAveAtom::end_of_step()
double FixAveAtom::memory_usage()
{
double bytes;
bytes = (double)atom->nmax*nvalues * sizeof(double);
bytes = (double)atom->nmax*values.size() * sizeof(double);
return bytes;
}
@ -403,7 +374,7 @@ double FixAveAtom::memory_usage()
void FixAveAtom::grow_arrays(int nmax)
{
memory->grow(array,nmax,nvalues,"fix_ave/atom:array");
memory->grow(array,nmax,values.size(),"fix_ave/atom:array");
array_atom = array;
if (array) vector_atom = array[0];
else vector_atom = nullptr;
@ -415,7 +386,7 @@ void FixAveAtom::grow_arrays(int nmax)
void FixAveAtom::copy_arrays(int i, int j, int /*delflag*/)
{
for (int m = 0; m < nvalues; m++)
for (std::size_t m = 0; m < values.size(); m++)
array[j][m] = array[i][m];
}
@ -425,8 +396,8 @@ void FixAveAtom::copy_arrays(int i, int j, int /*delflag*/)
int FixAveAtom::pack_exchange(int i, double *buf)
{
for (int m = 0; m < nvalues; m++) buf[m] = array[i][m];
return nvalues;
for (std::size_t m = 0; m < values.size(); m++) buf[m] = array[i][m];
return values.size();
}
/* ----------------------------------------------------------------------
@ -435,8 +406,8 @@ int FixAveAtom::pack_exchange(int i, double *buf)
int FixAveAtom::unpack_exchange(int nlocal, double *buf)
{
for (int m = 0; m < nvalues; m++) array[nlocal][m] = buf[m];
return nvalues;
for (std::size_t m = 0; m < values.size(); m++) array[nlocal][m] = buf[m];
return values.size();
}
/* ----------------------------------------------------------------------

15
src/fix_ave_atom.h
View File

@ -40,11 +40,20 @@ class FixAveAtom : public Fix {
int unpack_exchange(int, double *) override;
private:
int nvalues;
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 nrepeat, irepeat;
bigint nvalid, nvalid_last;
int *which, *argindex, *value2index;
char **ids;
double **array;
bigint nextvalid();

333
src/fix_ave_chunk.cpp
View File

@ -33,29 +33,24 @@
using namespace LAMMPS_NS;
using namespace FixConst;
enum{SCALAR,VECTOR};
enum{SAMPLE,ALL};
enum{NOSCALE,ATOM};
enum{ONE,RUNNING,WINDOW};
enum { SCALAR, VECTOR };
enum { SAMPLE, ALL };
enum { NOSCALE, ATOM };
enum { ONE, RUNNING, WINDOW };
/* ---------------------------------------------------------------------- */
FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
nvalues(0), nrepeat(0),
which(nullptr), argindex(nullptr), value2index(nullptr), ids(nullptr),
fp(nullptr), idchunk(nullptr), varatom(nullptr),
count_one(nullptr), count_many(nullptr), count_sum(nullptr),
values_one(nullptr), values_many(nullptr), values_sum(nullptr),
count_total(nullptr), count_list(nullptr),
values_total(nullptr), values_list(nullptr)
Fix(lmp, narg, arg), nvalues(0), nrepeat(0), fp(nullptr), idchunk(nullptr), varatom(nullptr),
count_one(nullptr), count_many(nullptr), count_sum(nullptr), values_one(nullptr),
values_many(nullptr), values_sum(nullptr), count_total(nullptr), count_list(nullptr),
values_total(nullptr), values_list(nullptr)
{
if (narg < 7) error->all(FLERR,"Illegal fix ave/chunk command");
if (narg < 7) utils::missing_cmd_args(FLERR, "fix ave/chunk", error);
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);
idchunk = utils::strdup(arg[6]);
@ -63,82 +58,81 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
no_change_box = 1;
time_depend = 1;
char * group = arg[1];
char *group = arg[1];
// expand args if any have wildcard character "*"
int expand = 0;
char **earg;
int nargnew = utils::expand_args(FLERR,narg-7,&arg[7],1,earg,lmp);
int nargnew = utils::expand_args(FLERR, narg - 7, &arg[7], 1, earg, lmp);
if (earg != &arg[7]) expand = 1;
arg = earg;
// parse values until one isn't recognized
which = new int[nargnew];
argindex = new int[nargnew];
ids = new char*[nargnew];
value2index = new int[nargnew];
densityflag = 0;
int iarg = 0;
values.clear();
while (iarg < nargnew) {
ids[nvalues] = nullptr;
value_t val;
val.id = "";
val.val.c = nullptr;
if (strcmp(arg[iarg],"vx") == 0) {
which[nvalues] = ArgInfo::V;
argindex[nvalues++] = 0;
val.which = ArgInfo::V;
val.which = 0;
} else if (strcmp(arg[iarg],"vy") == 0) {
which[nvalues] = ArgInfo::V;
argindex[nvalues++] = 1;
val.which = ArgInfo::V;
val.which = 1;
} else if (strcmp(arg[iarg],"vz") == 0) {
which[nvalues] = ArgInfo::V;
argindex[nvalues++] = 2;
val.which = ArgInfo::V;
val.which = 2;
} else if (strcmp(arg[iarg],"fx") == 0) {
which[nvalues] = ArgInfo::F;
argindex[nvalues++] = 0;
val.which = ArgInfo::F;
val.which = 0;
} else if (strcmp(arg[iarg],"fy") == 0) {
which[nvalues] = ArgInfo::F;
argindex[nvalues++] = 1;
val.which = ArgInfo::F;
val.which = 1;
} else if (strcmp(arg[iarg],"fz") == 0) {
which[nvalues] = ArgInfo::F;
argindex[nvalues++] = 2;
val.which = ArgInfo::F;
val.which = 2;
} else if (strcmp(arg[iarg],"mass") == 0) {
which[nvalues] = ArgInfo::MASS;
argindex[nvalues++] = 0;
val.which = ArgInfo::MASS;
val.which = 0;
} else if (strcmp(arg[iarg],"density/number") == 0) {
densityflag = 1;
which[nvalues] = ArgInfo::DENSITY_NUMBER;
argindex[nvalues++] = 0;
val.which = ArgInfo::DENSITY_NUMBER;
val.which = 0;
} else if (strcmp(arg[iarg],"density/mass") == 0) {
densityflag = 1;
which[nvalues] = ArgInfo::DENSITY_MASS;
argindex[nvalues++] = 0;
val.which = ArgInfo::DENSITY_MASS;
val.which = 0;
} else if (strcmp(arg[iarg],"temp") == 0) {
which[nvalues] = ArgInfo::TEMPERATURE;
argindex[nvalues++] = 0;
val.which = ArgInfo::TEMPERATURE;
val.which = 0;
} else {
ArgInfo argi(arg[iarg]);
if (argi.get_type() == ArgInfo::NONE) break;
if ((argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_dim() > 1))
error->all(FLERR,"Invalid fix ave/chunk command");
error->all(FLERR,"Unknown fix ave/chunk data value: {}", arg[iarg]);
which[nvalues] = argi.get_type();
argindex[nvalues] = argi.get_index1();
ids[nvalues] = argi.copy_name();
nvalues++;
val.which = argi.get_type();
val.argindex = argi.get_index1();
val.id = argi.get_name();
}
values.push_back(val);
iarg++;
}
if (nvalues == 0) error->all(FLERR,"No values in fix ave/chunk command");
nvalues = values.size();
if (nvalues == 0) error->all(FLERR, "No values in fix ave/chunk command");
// optional args
@ -159,7 +153,7 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
while (iarg < nargnew) {
if (strcmp(arg[iarg],"norm") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk norm", error);
if (strcmp(arg[iarg+1],"all") == 0) {
normflag = ALL;
scaleflag = ATOM;
@ -169,127 +163,124 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg+1],"none") == 0) {
normflag = SAMPLE;
scaleflag = NOSCALE;
} else error->all(FLERR,"Illegal fix ave/chunk command");
} else error->all(FLERR,"Unknown fix ave/chunk norm mode: {}", arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"ave") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk 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/chunk command");
else error->all(FLERR,"Unknown fix ave/chunk ave mode: {}", arg[iarg+1]);
if (ave == WINDOW) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk ave window", error);
nwindow = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/chunk command");
if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/chunk number of windows: {}", nwindow);
}
iarg += 2;
if (ave == WINDOW) iarg++;
} else if (strcmp(arg[iarg],"bias") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk bias", error);
biasflag = 1;
id_bias = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"adof") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk adof", error);
adof = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"cdof") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk cdof", error);
cdof = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"file") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk file", error);
if (comm->me == 0) {
fp = fopen(arg[iarg+1],"w");
if (fp == nullptr)
error->one(FLERR,"Cannot open fix ave/chunk file {}: {}",
arg[iarg+1], utils::getsyserror());
error->one(FLERR, "Cannot open fix ave/chunk file {}: {}",
arg[iarg+1], utils::getsyserror());
}
iarg += 2;
} else if (strcmp(arg[iarg],"overwrite") == 0) {
overwrite = 1;
iarg += 1;
} else if (strcmp(arg[iarg],"format") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk format", error);
delete[] format_user;
format_user = utils::strdup(arg[iarg+1]);
format = format_user;
iarg += 2;
} else if (strcmp(arg[iarg],"title1") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk 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/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk 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/chunk command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/chunk title3", error);
delete[] title3;
title3 = utils::strdup(arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal fix ave/chunk command");
} else error->all(FLERR,"Unknown fix ave/chunk keyword: {}", arg[iarg]);
}
// setup and error check
if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
error->all(FLERR,"Illegal fix ave/chunk command");
if (nevery <= 0) error->all(FLERR,"Illegal fix ave/chunk nevery value: {}", nevery);
if (nrepeat <= 0) error->all(FLERR,"Illegal fix ave/chunk nrepeat value: {}", nrepeat);
if (nfreq <= 0) error->all(FLERR,"Illegal fix ave/chunk nfreq value: {}", nfreq);
if (nfreq % nevery || nrepeat*nevery > nfreq)
error->all(FLERR,"Illegal fix ave/chunk command");
error->all(FLERR,"Inconsistent fix ave/chunk nevery/nrepeat/nfreq values");
if (ave != RUNNING && overwrite)
error->all(FLERR,"Illegal fix ave/chunk command");
error->all(FLERR,"Fix ave/chunk overwrite keyword requires ave running setting");
if (biasflag) {
int i = modify->find_compute(id_bias);
if (i < 0)
error->all(FLERR,"Could not find compute ID for temperature bias");
tbias = modify->compute[i];
tbias = modify->get_compute_by_id(id_bias);
if (!tbias) error->all(FLERR,"Could not find compute ID {} for temperature bias", id_bias);
if (tbias->tempflag == 0)
error->all(FLERR,"Bias compute does not calculate temperature");
error->all(FLERR,"Bias compute {} does not calculate temperature", id_bias);
if (tbias->tempbias == 0)
error->all(FLERR,"Bias compute does not calculate a velocity bias");
error->all(FLERR,"Bias compute {} does not calculate a velocity bias", id_bias);
}
for (int i = 0; i < nvalues; i++) {
if (which[i] == ArgInfo::COMPUTE) {
auto icompute = modify->get_compute_by_id(ids[i]);
if (!icompute)
error->all(FLERR,"Compute ID {} for fix ave/chunk does not exist",ids[i]);
if (icompute->peratom_flag == 0)
error->all(FLERR,"Fix ave/chunk compute {} does not calculate per-atom values",ids[i]);
if (argindex[i] == 0 && (icompute->size_peratom_cols != 0))
error->all(FLERR,"Fix ave/chunk compute {} does not calculate a per-atom vector",ids[i]);
if (argindex[i] && (icompute->size_peratom_cols == 0))
error->all(FLERR,"Fix ave/chunk compute {} does not calculate a per-atom array",ids[i]);
if (argindex[i] && (argindex[i] > icompute->size_peratom_cols))
error->all(FLERR,"Fix ave/chunk compute {} vector is accessed out-of-range",ids[i]);
for (auto &val : values) {
} else if (which[i] == ArgInfo::FIX) {
auto ifix = modify->get_fix_by_id(ids[i]);
if (!ifix)
error->all(FLERR, "Fix ID {} for fix ave/chunk does not exist",ids[i]);
if (ifix->peratom_flag == 0)
error->all(FLERR, "Fix ave/chunk fix {} does not calculate per-atom values",ids[i]);
if (argindex[i] == 0 && (ifix->size_peratom_cols != 0))
error->all(FLERR, "Fix ave/chunk fix {} does not calculate a per-atom vector",ids[i]);
if (argindex[i] && (ifix->size_peratom_cols == 0))
error->all(FLERR, "Fix ave/chunk fix {} does not calculate a per-atom array",ids[i]);
if (argindex[i] && argindex[i] > ifix->size_peratom_cols)
error->all(FLERR,"Fix ave/chunk fix {} vector is accessed out-of-range",ids[i]);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name {} for fix ave/chunk does not exist",ids[i]);
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Fix ave/chunk variable {} is not atom-style variable",ids[i]);
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 fix ave/chunk does not exist",val.id);
if (val.val.c->peratom_flag == 0)
error->all(FLERR,"Fix ave/chunk compute {} does not calculate per-atom values",val.id);
if (val.argindex == 0 && (val.val.c->size_peratom_cols != 0))
error->all(FLERR,"Fix ave/chunk compute {} does not calculate a per-atom vector",val.id);
if (val.argindex && (val.val.c->size_peratom_cols == 0))
error->all(FLERR,"Fix ave/chunk compute {} does not calculate a per-atom array",val.id);
if (val.argindex && (val.argindex > val.val.c->size_peratom_cols))
error->all(FLERR,"Fix ave/chunk compute {} vector is accessed out-of-range",val.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 fix ave/chunk does not exist",val.id);
if (val.val.f->peratom_flag == 0)
error->all(FLERR, "Fix ave/chunk fix {} does not calculate per-atom values",val.id);
if (val.argindex == 0 && (val.val.f->size_peratom_cols != 0))
error->all(FLERR, "Fix ave/chunk fix {} does not calculate a per-atom vector",val.id);
if (val.argindex && (val.val.f->size_peratom_cols == 0))
error->all(FLERR, "Fix ave/chunk fix {} does not calculate a per-atom array",val.id);
if (val.argindex && val.argindex > val.val.f->size_peratom_cols)
error->all(FLERR,"Fix ave/chunk fix {} vector is accessed out-of-range",val.id);
} 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 fix ave/chunk does not exist",val.id);
if (input->variable->atomstyle(val.val.v) == 0)
error->all(FLERR,"Fix ave/chunk variable {} is not atom-style variable",val.id);
}
}
@ -396,12 +387,6 @@ FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
FixAveChunk::~FixAveChunk()
{
delete[] which;
delete[] argindex;
for (int i = 0; i < nvalues; i++) delete[] ids[i];
delete[] ids;
delete[] value2index;
if (fp && comm->me == 0) fclose(fp);
memory->destroy(varatom);
@ -427,10 +412,6 @@ FixAveChunk::~FixAveChunk()
}
delete[] idchunk;
which = nullptr;
argindex = nullptr;
ids = nullptr;
value2index = nullptr;
fp = nullptr;
varatom = nullptr;
count_one = nullptr;
@ -463,42 +444,34 @@ void FixAveChunk::init()
// set indices and check validity of all computes,fixes,variables
// check that fix frequency is acceptable
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for fix ave/chunk");
cchunk = dynamic_cast<ComputeChunkAtom *>(modify->compute[icompute]);
cchunk = dynamic_cast<ComputeChunkAtom *>(modify->get_compute_by_id(idchunk));
if (!cchunk)
error->all(FLERR,"Chunk/atom compute {} does not exist or is "
"incorrect style for fix ave/chunk",idchunk);
if (biasflag) {
int i = modify->find_compute(id_bias);
if (i < 0)
error->all(FLERR,"Could not find compute ID for temperature bias");
tbias = modify->compute[i];
tbias = modify->get_compute_by_id(id_bias);
if (!tbias)
error->all(FLERR,"Could not find compute ID {} for temperature bias", id_bias);
}
for (int m = 0; m < nvalues; m++) {
if (which[m] == ArgInfo::COMPUTE) {
icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID for fix ave/chunk does not exist");
value2index[m] = icompute;
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 fix ave/chunk does not exist", val.id);
} else if (which[m] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/chunk does not exist");
value2index[m] = ifix;
} 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 fix ave/chunk does not exist", val.id);
if (nevery % modify->fix[ifix]->peratom_freq)
error->all(FLERR,
"Fix for fix ave/chunk not computed at compatible time");
if (nevery % val.val.f->peratom_freq)
error->all(FLERR, "Fix {} for fix ave/chunk not computed at compatible time", val.id);
} else if (which[m] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/chunk does not exist");
value2index[m] = ivariable;
} else value2index[m] = -1;
} 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 fix ave/chunk does not exist", val.id);
}
}
// need to reset nvalid if nvalid < ntimestep b/c minimize was performed
@ -527,7 +500,7 @@ void FixAveChunk::setup(int /*vflag*/)
void FixAveChunk::end_of_step()
{
int i,j,m,n,index;
int i,j,m,index;
// skip if not step which requires doing something
@ -609,15 +582,15 @@ void FixAveChunk::end_of_step()
modify->clearstep_compute();
for (m = 0; m < nvalues; m++) {
n = value2index[m];
j = argindex[m];
m = 0;
for (auto &val : values) {
j = val.argindex;
// V,F adds velocities,forces to values
if (which[m] == ArgInfo::V || which[m] == ArgInfo::F) {
if (val.which == ArgInfo::V || val.which == ArgInfo::F) {
double **attribute;
if (which[m] == ArgInfo::V) attribute = atom->v;
if (val.which == ArgInfo::V) attribute = atom->v;
else attribute = atom->f;
for (i = 0; i < nlocal; i++)
@ -628,7 +601,7 @@ void FixAveChunk::end_of_step()
// DENSITY_NUMBER adds 1 to values
} else if (which[m] == ArgInfo::DENSITY_NUMBER) {
} else if (val.which == ArgInfo::DENSITY_NUMBER) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
@ -638,8 +611,7 @@ void FixAveChunk::end_of_step()
// DENSITY_MASS or MASS adds mass to values
} else if ((which[m] == ArgInfo::DENSITY_MASS)
|| (which[m] == ArgInfo::MASS)) {
} else if ((val.which == ArgInfo::DENSITY_MASS) || (val.which == ArgInfo::MASS)) {
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
@ -661,7 +633,7 @@ void FixAveChunk::end_of_step()
// TEMPERATURE adds KE to values
// subtract and restore velocity bias if requested
} else if (which[m] == ArgInfo::TEMPERATURE) {
} else if (val.which == ArgInfo::TEMPERATURE) {
if (biasflag) {
if (tbias->invoked_scalar != ntimestep) tbias->compute_scalar();
@ -695,14 +667,13 @@ void FixAveChunk::end_of_step()
// COMPUTE adds its scalar or vector component to values
// invoke compute if not previously invoked
} else if (which[m] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[n];
if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM;
} else if (val.which == ArgInfo::COMPUTE) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_PERATOM)) {
val.val.c->compute_peratom();
val.val.c->invoked_flag |= Compute::INVOKED_PERATOM;
}
double *vector = compute->vector_atom;
double **array = compute->array_atom;
double *vector = val.val.c->vector_atom;
double **array = val.val.c->array_atom;
int jm1 = j - 1;
for (i = 0; i < nlocal; i++)
@ -715,9 +686,9 @@ void FixAveChunk::end_of_step()
// FIX adds its scalar or vector component to values
// access fix fields, guaranteed to be ready
} else if (which[m] == ArgInfo::FIX) {
double *vector = modify->fix[n]->vector_atom;
double **array = modify->fix[n]->array_atom;
} else if (val.which == ArgInfo::FIX) {
double *vector = val.val.f->vector_atom;
double **array = val.val.f->array_atom;
int jm1 = j - 1;
for (i = 0; i < nlocal; i++)
@ -730,14 +701,14 @@ void FixAveChunk::end_of_step()
// VARIABLE adds its per-atom quantities to values
// evaluate atom-style variable
} else if (which[m] == ArgInfo::VARIABLE) {
} else if (val.which == ArgInfo::VARIABLE) {
if (atom->nmax > maxvar) {
maxvar = atom->nmax;
memory->destroy(varatom);
memory->create(varatom,maxvar,"ave/chunk:varatom");
}
input->variable->compute_atom(n,igroup,varatom,1,0);
input->variable->compute_atom(val.val.v,igroup,varatom,1,0);
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
@ -745,6 +716,7 @@ void FixAveChunk::end_of_step()
values_one[index][m] += varatom[i];
}
}
++m;
}
// process the current sample
@ -784,14 +756,14 @@ void FixAveChunk::end_of_step()
for (m = 0; m < nchunk; m++) {
if (count_many[m] > 0.0)
for (j = 0; j < nvalues; j++) {
if (which[j] == ArgInfo::TEMPERATURE) {
if (values[j].which == ArgInfo::TEMPERATURE) {
values_many[m][j] += mvv2e*values_one[m][j] /
((cdof + adof*count_many[m]) * boltz);
} else if (which[j] == ArgInfo::DENSITY_NUMBER) {
} else if (values[j].which == ArgInfo::DENSITY_NUMBER) {
if (volflag == SCALAR) values_one[m][j] /= chunk_volume_scalar;
else values_one[m][j] /= chunk_volume_vec[m];
values_many[m][j] += values_one[m][j];
} else if (which[j] == ArgInfo::DENSITY_MASS) {
} else if (values[j].which == ArgInfo::DENSITY_MASS) {
if (volflag == SCALAR) values_one[m][j] /= chunk_volume_scalar;
else values_one[m][j] /= chunk_volume_vec[m];
values_many[m][j] += mv2d*values_one[m][j];
@ -854,13 +826,13 @@ void FixAveChunk::end_of_step()
for (m = 0; m < nchunk; m++) {
if (count_sum[m] > 0.0)
for (j = 0; j < nvalues; j++) {
if (which[j] == ArgInfo::TEMPERATURE) {
if (values[j].which == ArgInfo::TEMPERATURE) {
values_sum[m][j] *= mvv2e/((repeat*cdof + adof*count_sum[m])*boltz);
} else if (which[j] == ArgInfo::DENSITY_NUMBER) {
} else if (values[j].which == ArgInfo::DENSITY_NUMBER) {
if (volflag == SCALAR) values_sum[m][j] /= chunk_volume_scalar;
else values_sum[m][j] /= chunk_volume_vec[m];
values_sum[m][j] /= repeat;
} else if (which[j] == ArgInfo::DENSITY_MASS) {
} else if (values[j].which == ArgInfo::DENSITY_MASS) {
if (volflag == SCALAR) values_sum[m][j] /= chunk_volume_scalar;
else values_sum[m][j] /= chunk_volume_vec[m];
values_sum[m][j] *= mv2d/repeat;
@ -1045,8 +1017,7 @@ void FixAveChunk::allocate()
if (ave == WINDOW) {
memory->create(count_list,nwindow,nchunk,"ave/chunk:count_list");
memory->create(values_list,nwindow,nchunk,nvalues,
"ave/chunk:values_list");
memory->create(values_list,nwindow,nchunk,nvalues,"ave/chunk:values_list");
}
// reinitialize regrown count/values total since they accumulate

17
src/fix_ave_chunk.h
View File

@ -36,15 +36,24 @@ class FixAveChunk : public Fix {
double memory_usage() override;
private:
int 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;
int normflag, scaleflag, overwrite, biasflag, colextra;
bigint nvalid, nvalid_last;
double adof, cdof;
char *format, *format_user;
char *tstring, *sstring, *id_bias;
int *which, *argindex, *value2index;
char **ids;
class Compute *tbias; // ptr to additional bias compute
FILE *fp;

261
src/fix_ave_correlate.cpp
View File

@ -21,6 +21,7 @@
#include "fix_ave_correlate.h"
#include "arg_info.h"
#include "comm.h"
#include "compute.h"
#include "error.h"
#include "input.h"
@ -34,23 +35,20 @@
using namespace LAMMPS_NS;
using namespace FixConst;
enum{ONE,RUNNING};
enum{AUTO,UPPER,LOWER,AUTOUPPER,AUTOLOWER,FULL};
enum { ONE, RUNNING };
enum { AUTO, UPPER, LOWER, AUTOUPPER, AUTOLOWER, FULL };
/* ---------------------------------------------------------------------- */
FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
Fix (lmp, narg, arg),
nvalues(0), which(nullptr), argindex(nullptr), value2index(nullptr), ids(nullptr), fp(nullptr),
count(nullptr), values(nullptr), corr(nullptr), save_count(nullptr), save_corr(nullptr)
FixAveCorrelate::FixAveCorrelate(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), fp(nullptr), count(nullptr), cvalues(nullptr), corr(nullptr),
save_count(nullptr), save_corr(nullptr)
{
if (narg < 7) error->all(FLERR,"Illegal fix ave/correlate command");
if (narg < 7) utils::missing_cmd_args(FLERR, "fix ave/correlate", 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);
time_depend = 1;
global_freq = nfreq;
@ -59,34 +57,31 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
int expand = 0;
char **earg;
int nargnew = utils::expand_args(FLERR,narg-6,&arg[6],0,earg,lmp);
int nargnew = utils::expand_args(FLERR, narg - 6, &arg[6], 0, earg, lmp);
if (earg != &arg[6]) expand = 1;
arg = earg;
// parse values until one isn't recognized
which = new int[nargnew];
argindex = new int[nargnew];
ids = new char*[nargnew];
value2index = new int[nargnew];
nvalues = 0;
// parse values
int iarg = 0;
while (iarg < nargnew) {
ArgInfo argi(arg[iarg]);
value_t val;
if (argi.get_type() == ArgInfo::NONE) break;
if ((argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_dim() > 1))
error->all(FLERR,"Invalid fix ave/correlate command");
error->all(FLERR, "Unknown fix ave/correlate data type: {}", arg[iarg]);
which[nvalues] = argi.get_type();
argindex[nvalues] = argi.get_index1();
ids[nvalues] = argi.copy_name();
val.which = argi.get_type();
val.argindex = argi.get_index1();
val.id = argi.get_name();
val.val.c = nullptr;
nvalues++;
values.push_back(val);
iarg++;
}
nvalues = values.size();
// optional args
@ -102,32 +97,32 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
while (iarg < nargnew) {
if (strcmp(arg[iarg],"type") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/correlate command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate type", error);
if (strcmp(arg[iarg+1],"auto") == 0) type = AUTO;
else if (strcmp(arg[iarg+1],"upper") == 0) type = UPPER;
else if (strcmp(arg[iarg+1],"lower") == 0) type = LOWER;
else if (strcmp(arg[iarg+1],"auto/upper") == 0) type = AUTOUPPER;
else if (strcmp(arg[iarg+1],"auto/lower") == 0) type = AUTOLOWER;
else if (strcmp(arg[iarg+1],"full") == 0) type = FULL;
else error->all(FLERR,"Illegal fix ave/correlate command");
else error->all(FLERR,"Unknown fix ave/correlate type: {}");
iarg += 2;
} else if (strcmp(arg[iarg],"ave") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/correlate command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate ave", error);
if (strcmp(arg[iarg+1],"one") == 0) ave = ONE;
else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING;
else error->all(FLERR,"Illegal fix ave/correlate command");
else error->all(FLERR,"Unknown fix ave/correlate ave mode: {}", arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"start") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/correlate command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate start", error);
startstep = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"prefactor") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/correlate command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate prefactor", error);
prefactor = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"file") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/correlate command");
if (me == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate file", error);
if (comm->me == 0) {
fp = fopen(arg[iarg+1],"w");
if (fp == nullptr)
error->one(FLERR,"Cannot open fix ave/correlate file {}:"" {}",
@ -138,75 +133,67 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
overwrite = 1;
iarg += 1;
} else if (strcmp(arg[iarg],"title1") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/correlate command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate 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/correlate command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate 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/correlate command");
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix ave/correlate title3", error);
delete[] title3;
title3 = utils::strdup(arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal fix ave/correlate command");
} else error->all(FLERR,"Unkown fix ave/correlate keyword: {}", arg[iarg]);
}
// setup and error check
// for fix inputs, check that fix frequency is acceptable
if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
error->all(FLERR,"Illegal fix ave/correlate command");
if (nfreq % nevery)
error->all(FLERR,"Illegal fix ave/correlate command");
if (ave == ONE && nfreq < (nrepeat-1)*nevery)
error->all(FLERR,"Illegal fix ave/correlate command");
if (nevery <= 0) error->all(FLERR,"Illegal fix ave/correlate nevery value: {}", nevery);
if (nrepeat <= 0) error->all(FLERR,"Illegal fix ave/correlate nrepeat value: {}", nrepeat);
if (nfreq <= 0) error->all(FLERR,"Illegal fix ave/correlate nfreq value: {}", nfreq);
if (nfreq % nevery || nrepeat*nevery > nfreq)
error->all(FLERR,"Inconsistent fix ave/correlate nevery/nrepeat/nfreq values");
if (ave != RUNNING && overwrite)
error->all(FLERR,"Illegal fix ave/correlate command");
error->all(FLERR,"Fix ave/correlate overwrite keyword requires ave running setting");
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/correlate does not exist");
if (argindex[i] == 0 && modify->compute[icompute]->scalar_flag == 0)
error->all(FLERR,
"Fix ave/correlate compute does not calculate a scalar");
if (argindex[i] && modify->compute[icompute]->vector_flag == 0)
error->all(FLERR,
"Fix ave/correlate compute does not calculate a vector");
if (argindex[i] && argindex[i] > modify->compute[icompute]->size_vector)
error->all(FLERR,"Fix ave/correlate compute vector "
"is accessed out-of-range");
for (auto &val : values) {
} else if (which[i] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/correlate does not exist");
if (argindex[i] == 0 && modify->fix[ifix]->scalar_flag == 0)
error->all(FLERR,"Fix ave/correlate fix does not calculate a scalar");
if (argindex[i] && modify->fix[ifix]->vector_flag == 0)
error->all(FLERR,"Fix ave/correlate fix does not calculate a vector");
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_vector)
error->all(FLERR,
"Fix ave/correlate fix vector is accessed out-of-range");
if (nevery % modify->fix[ifix]->global_freq)
error->all(FLERR,"Fix for fix ave/correlate "
"not computed at compatible time");
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 fix ave/correlate does not exist", val.id);
if (val.argindex == 0 && val.val.c->scalar_flag == 0)
error->all(FLERR, "Fix ave/correlate compute {} does not calculate a scalar", val.id);
if (val.argindex && val.val.c->vector_flag == 0)
error->all(FLERR, "Fix ave/correlate compute {} does not calculate a vector", val.id);
if (val.argindex && val.argindex > val.val.c->size_vector)
error->all(FLERR, "Fix ave/correlate compute {} vector is accessed out-of-range", val.id);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/correlate does not exist");
if (argindex[i] == 0 && input->variable->equalstyle(ivariable) == 0)
error->all(FLERR,
"Fix ave/correlate variable is not equal-style variable");
if (argindex[i] && input->variable->vectorstyle(ivariable) == 0)
error->all(FLERR,
"Fix ave/correlate variable is not vector-style variable");
} 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 fix ave/correlate does not exist", val.id);
if (val.argindex == 0 && val.val.f->scalar_flag == 0)
error->all(FLERR, "Fix ave/correlate fix {} does not calculate a scalar", val.id);
if (val.argindex && val.val.f->vector_flag == 0)
error->all(FLERR, "Fix ave/correlate fix {} does not calculate a vector", val.id);
if (val.argindex && val.argindex > val.val.f->size_vector)
error->all(FLERR, "Fix ave/correlate fix {} vector is accessed out-of-range", val.id);
if (nevery % val.val.f->global_freq)
error->all(FLERR, "Fix {} for fix ave/correlate not computed at compatible time", val.id);
} 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 fix ave/correlate does not exist", val.id);
if (val.argindex == 0 && input->variable->equalstyle(val.val.v) == 0)
error->all(FLERR, "Fix ave/correlate variable {} is not equal-style variable", val.id);
if (val.argindex && input->variable->vectorstyle(val.val.v) == 0)
error->all(FLERR, "Fix ave/correlate variable {} is not vector-style variable", val.id);
}
}
@ -219,7 +206,7 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
// 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,"# Time-correlated data for fix %s\n",id);
@ -254,7 +241,7 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
fprintf(fp,"\n");
}
if (ferror(fp))
error->one(FLERR,"Error writing file header");
error->one(FLERR,"Error writing ave/correlate header: {}", utils::getsyserror());
filepos = platform::ftell(fp);
}
@ -275,7 +262,7 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
// set count and corr to zero since they accumulate
// also set save versions to zero in case accessed via compute_array()
memory->create(values,nrepeat,nvalues,"ave/correlate:values");
memory->create(cvalues,nrepeat,nvalues,"ave/correlate:cvalues");
memory->create(count,nrepeat,"ave/correlate:count");
memory->create(save_count,nrepeat,"ave/correlate:save_count");
memory->create(corr,nrepeat,npair,"ave/correlate:corr");
@ -312,19 +299,13 @@ FixAveCorrelate::FixAveCorrelate(LAMMPS * lmp, int narg, char **arg):
FixAveCorrelate::~FixAveCorrelate()
{
delete[] which;
delete[] argindex;
delete[] value2index;
for (int i = 0; i < nvalues; i++) delete[] ids[i];
delete[] ids;
memory->destroy(values);
memory->destroy(cvalues);
memory->destroy(count);
memory->destroy(save_count);
memory->destroy(corr);
memory->destroy(save_corr);
if (fp && me == 0) fclose(fp);
if (fp && comm->me == 0) fclose(fp);
}
/* ---------------------------------------------------------------------- */
@ -342,24 +323,21 @@ void FixAveCorrelate::init()
{
// set current indices for all computes,fixes,variables
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/correlate does not exist");
value2index[i] = icompute;
for (auto &val : values) {
} else if (which[i] == ArgInfo::FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/correlate does not exist");
value2index[i] = ifix;
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 fix ave/correlate does not exist", val.id);
} else if (which[i] == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/correlate does not exist");
value2index[i] = ivariable;
} 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 fix ave/correlate does not exist", val.id);
} 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 fix ave/correlate does not exist", val.id);
}
}
@ -387,8 +365,7 @@ void FixAveCorrelate::setup(int /*vflag*/)
void FixAveCorrelate::end_of_step()
{
int i,j,m;
double scalar;
int i,j;
// skip if not step which requires doing something
@ -406,51 +383,51 @@ void FixAveCorrelate::end_of_step()
lastindex++;
if (lastindex == nrepeat) lastindex = 0;
for (i = 0; i < nvalues; i++) {
m = value2index[i];
i = 0;
for (auto &val : values) {
double scalar = 0.0;
// invoke compute if not previously invoked
if (which[i] == ArgInfo::COMPUTE) {
Compute *compute = modify->compute[m];
if (val.which == ArgInfo::COMPUTE) {
if (argindex[i] == 0) {
if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) {
compute->compute_scalar();
compute->invoked_flag |= Compute::INVOKED_SCALAR;
if (val.argindex == 0) {
if (!(val.val.c->invoked_flag & Compute::INVOKED_SCALAR)) {
val.val.c->compute_scalar();
val.val.c->invoked_flag |= Compute::INVOKED_SCALAR;
}
scalar = compute->scalar;
scalar = 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;
}
scalar = compute->vector[argindex[i]-1];
scalar = val.val.c->vector[val.argindex-1];
}
// access fix fields, guaranteed to be ready
} else if (which[i] == ArgInfo::FIX) {
if (argindex[i] == 0)
scalar = modify->fix[m]->compute_scalar();
} else if (val.which == ArgInfo::FIX) {
if (val.argindex == 0)
scalar = val.val.f->compute_scalar();
else
scalar = modify->fix[m]->compute_vector(argindex[i]-1);
scalar = val.val.f->compute_vector(val.argindex-1);
// evaluate equal-style or vector-style variable
} else if (which[i] == ArgInfo::VARIABLE) {
if (argindex[i] == 0)
scalar = input->variable->compute_equal(m);
} else if (val.which == ArgInfo::VARIABLE) {
if (val.argindex == 0)
scalar = input->variable->compute_equal(val.val.v);
else {
double *varvec;
int nvec = input->variable->compute_vector(m,&varvec);
int index = argindex[i];
int nvec = input->variable->compute_vector(val.val.v,&varvec);
int index = val.argindex;
if (nvec < index) scalar = 0.0;
else scalar = varvec[index-1];
}
}
values[lastindex][i] = scalar;
cvalues[lastindex][i++] = scalar;
}
// fistindex = index in values ring of earliest time sample
@ -484,7 +461,7 @@ void FixAveCorrelate::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,nrepeat);
@ -499,7 +476,7 @@ void FixAveCorrelate::end_of_step()
fprintf(fp,"\n");
}
if (ferror(fp))
error->one(FLERR,"Error writing out correlation data");
error->one(FLERR,"Error writing out fix ave/correlate data: {}", utils::getsyserror());
fflush(fp);
@ -539,7 +516,7 @@ void FixAveCorrelate::accumulate()
for (k = 0; k < nsample; k++) {
ipair = 0;
for (i = 0; i < nvalues; i++) {
corr[k][ipair++] += values[m][i]*values[n][i];
corr[k][ipair++] += cvalues[m][i]*cvalues[n][i];
}
m--;
if (m < 0) m = nrepeat-1;
@ -550,7 +527,7 @@ void FixAveCorrelate::accumulate()
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = i+1; j < nvalues; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
corr[k][ipair++] += cvalues[m][i]*cvalues[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
@ -560,7 +537,7 @@ void FixAveCorrelate::accumulate()
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = 0; j < i; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
corr[k][ipair++] += cvalues[m][i]*cvalues[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
@ -570,7 +547,7 @@ void FixAveCorrelate::accumulate()
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = i; j < nvalues; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
corr[k][ipair++] += cvalues[m][i]*cvalues[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
@ -580,7 +557,7 @@ void FixAveCorrelate::accumulate()
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = 0; j <= i; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
corr[k][ipair++] += cvalues[m][i]*cvalues[n][j];
m--;
if (m < 0) m = nrepeat-1;
}
@ -590,7 +567,7 @@ void FixAveCorrelate::accumulate()
ipair = 0;
for (i = 0; i < nvalues; i++)
for (j = 0; j < nvalues; j++)
corr[k][ipair++] += values[m][i]*values[n][j];
corr[k][ipair++] += cvalues[m][i]*cvalues[n][j];
m--;
if (m < 0) m = nrepeat-1;
}

19
src/fix_ave_correlate.h
View File

@ -35,11 +35,20 @@ class FixAveCorrelate : public Fix {
double compute_array(int, int) override;
private:
int me, nvalues;
int nrepeat, nfreq;
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;
bigint nvalid, nvalid_last;
int *which, *argindex, *value2index;
char **ids;
FILE *fp;
int type, ave, startstep, overwrite;
@ -52,7 +61,7 @@ class FixAveCorrelate : public Fix {
int npair; // number of correlation pairs to calculate
int *count;
double **values, **corr;
double **cvalues, **corr;
int *save_count; // saved values at Nfreq for output via compute_array()
double **save_corr;

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

15
src/fix_ave_time.cpp
View File

@ -92,14 +92,12 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
value_t val;
val.keyword = arg[i];
val.which = argi.get_type();
key2col[arg[i]] = i;
if ((argi.get_type() == ArgInfo::NONE)
|| (argi.get_type() == ArgInfo::UNKNOWN)
|| (argi.get_dim() > 1))
if ((val.which == ArgInfo::NONE) || (val.which == ArgInfo::UNKNOWN) || (argi.get_dim() > 1))
error->all(FLERR,"Invalid fix ave/time argument: {}", arg[i]);
val.which = argi.get_type();
val.argindex = argi.get_index1();
val.varlen = 0;
val.offcol = 0;
@ -123,12 +121,9 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
// for fix inputs, check that fix frequency is acceptable
// set variable_length if any compute is variable length
if (nevery <= 0)
error->all(FLERR,"Illegal fix ave/time nevery value: {}", nevery);
if (nrepeat <= 0)
error->all(FLERR,"Illegal fix ave/time nrepeat value: {}", nrepeat);
if (nfreq <= 0)
error->all(FLERR,"Illegal fix ave/time nfreq value: {}", nfreq);
if (nevery <= 0) error->all(FLERR,"Illegal fix ave/time nevery value: {}", nevery);
if (nrepeat <= 0) error->all(FLERR,"Illegal fix ave/time nrepeat value: {}", nrepeat);
if (nfreq <= 0) error->all(FLERR,"Illegal fix ave/time nfreq value: {}", nfreq);
if (nfreq % nevery || nrepeat*nevery > nfreq)
error->all(FLERR,"Inconsistent fix ave/time nevery/nrepeat/nfreq values");
if (ave != RUNNING && overwrite)

623
src/fix_store_state.cpp
View File

File diff suppressed because it is too large Load Diff

25
src/fix_store_state.h
View File

@ -44,10 +44,22 @@ class FixStoreState : public Fix {
int maxsize_restart() override;
private:
int nvalues;
int *which, *argindex, *value2index;
char **ids;
double **values; // archived atom properties
struct value_t {
int which;
int argindex;
std::string id;
union {
class Compute *c;
class Fix *f;
int v;
int d;
int i;
} val;
void (FixStoreState::* pack_choice)(int); // ptr to pack function
};
std::vector<value_t> values;
double **avalues; // archived atom properties
double *vbuf; // 1d ptr to values
int comflag;
@ -56,9 +68,6 @@ class FixStoreState : public Fix {
int kflag, cfv_flag, firstflag;
int cfv_any; // 1 if any compute/fix/variable specified
typedef void (FixStoreState::*FnPtrPack)(int);
FnPtrPack *pack_choice; // ptrs to pack functions
void pack_id(int);
void pack_molecule(int);
void pack_type(int);
@ -113,8 +122,6 @@ class FixStoreState : public Fix {
void pack_tqy(int);
void pack_tqz(int);
};
} // namespace LAMMPS_NS
#endif
#endif

62
unittest/commands/test_compute_chunk.cpp
View File

@ -246,8 +246,8 @@ TEST_F(ComputeChunkTest, ChunkComputes)
command("compute trq all torque/chunk mols");
command("compute vcm all vcm/chunk mols");
command("fix hist1 all ave/time 1 1 1 c_ang[*] c_com[*] c_dip[*] c_gyr c_mom[*] c_omg[*] "
"c_trq[*] c_vcm[*] mode vector file all_chunk.dat format %15.8f");
command("fix hist2 all ave/time 1 1 1 c_tmp mode vector file vec_chunk.dat format %15.8f");
"c_trq[*] c_vcm[*] mode vector");
command("fix hist2 all ave/time 1 1 1 c_tmp mode vector");
command("run 0 post no");
END_HIDE_OUTPUT();
auto cang = get_array("ang");
@ -314,6 +314,64 @@ TEST_F(ComputeChunkTest, ChunkComputes)
EXPECT_NEAR(ctmp[4], -0.06062938, EPSILON);
EXPECT_NEAR(ctmp[5], -0.09219489, EPSILON);
}
TEST_F(ComputeChunkTest, ChunkSpreadGlobal)
{
if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP();
BEGIN_HIDE_OUTPUT();
command("pair_style lj/cut/coul/cut 10.0");
command("pair_coeff * * 0.01 3.0");
command("bond_style harmonic");
command("bond_coeff * 100.0 1.5");
command("compute gyr all gyration/chunk mols");
command("compute spr all chunk/spread/atom mols c_gyr");
command("compute glb all global/atom c_mols c_gyr");
command("variable odd atom ((id+1)%2)+1");
command("compute odd all global/atom v_odd c_gyr");
command("fix ave all ave/atom 1 1 1 c_spr c_glb c_odd");
command("run 0 post no");
END_HIDE_OUTPUT();
const int natoms = lammps_get_natoms(lmp);
auto cgyr = get_vector("gyr");
auto cspr = get_peratom("spr");
auto cglb = get_peratom("glb");
auto codd = get_peratom("odd");
auto ctag = get_peratom("tags");
for (int i = 0; i < natoms; ++i) {
EXPECT_EQ(cspr[i], cgyr[chunkmol[(int)ctag[i]] - 1]);
EXPECT_EQ(cglb[i], cgyr[chunkmol[(int)ctag[i]] - 1]);
EXPECT_EQ(codd[i], cgyr[(((int)ctag[i] + 1) % 2)]);
}
}
TEST_F(ComputeChunkTest, ChunkReduce)
{
if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP();
BEGIN_HIDE_OUTPUT();
command("pair_style lj/cut/coul/cut 10.0");
command("pair_coeff * * 0.01 3.0");
command("bond_style harmonic");
command("bond_coeff * 100.0 1.5");
command("compute prp all property/chunk mols count");
command("variable one atom 1");
command("compute red all reduce/chunk mols sum v_one");
command("fix ave all ave/time 1 1 1 c_prp c_red mode vector");
command("run 0 post no");
END_HIDE_OUTPUT();
const int nchunks = get_scalar("mols");
auto cprp = get_vector("prp");
auto cred = get_vector("red");
for (int i = 0; i < nchunks; ++i)
EXPECT_EQ(cprp[i], cred[i]);
}
} // namespace LAMMPS_NS
int main(int argc, char **argv)