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

convert fix ave/correlate and fix ave/correlate/long

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-10-16 11:47:11 -04:00
parent 26468eccca
commit ec121ab812
4 changed files with 431 additions and 407 deletions
Download Patch File Download Diff File Expand all files Collapse all files

539
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) {
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::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;
} 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 (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::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;
@ -405,44 +423,52 @@ void FixAveCorrelateLong::end_of_step()
// compute/fix/variable may invoke computes so wrap with clear/add
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();

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) {
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::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;
} 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 (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::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;