From ec121ab812253e358cc1561ce574cdbefe7df783 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 16 Oct 2022 11:47:11 -0400 Subject: [PATCH] convert fix ave/correlate and fix ave/correlate/long --- src/EXTRA-FIX/fix_ave_correlate_long.cpp | 539 ++++++++++++----------- src/EXTRA-FIX/fix_ave_correlate_long.h | 19 +- src/fix_ave_correlate.cpp | 261 +++++------ src/fix_ave_correlate.h | 19 +- 4 files changed, 431 insertions(+), 407 deletions(-) diff --git a/src/EXTRA-FIX/fix_ave_correlate_long.cpp b/src/EXTRA-FIX/fix_ave_correlate_long.cpp index 9ed8a1e466..d81fce1673 100644 --- a/src/EXTRA-FIX/fix_ave_correlate_long.cpp +++ b/src/EXTRA-FIX/fix_ave_correlate_long.cpp @@ -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;idestroy(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;idt*nevery); - for (int j=0;jone(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 0) { t[jm] = j; - for (int i=0;i0) { - t[jm] = j * pow((double)m, k); + t[jm] = j * powint((double)m, k); for (int i=0;intimestep <= last_accumulated_step) return; if (type == AUTO) { - for (i=0; intimestep; @@ -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 -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 -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 -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 -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 (list[n++]); + int npairin = static_cast(list[n++]); int numcorrelatorsin = static_cast (list[n++]); - int pin = static_cast (list[n++]); - int min = static_cast (list[n++]); - last_accumulated_step = static_cast (list[n++]); + int pin = static_cast(list[n++]); + int min = static_cast(list[n++]); + last_accumulated_step = static_cast(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(list[n++]); - naccumulator[i] = static_cast (list[n++]); - insertindex[i] = static_cast (list[n++]); + naccumulator[i] = static_cast(list[n++]); + insertindex[i] = static_cast(list[n++]); } } diff --git a/src/EXTRA-FIX/fix_ave_correlate_long.h b/src/EXTRA-FIX/fix_ave_correlate_long.h index 0c540aa727..0495664934 100644 --- a/src/EXTRA-FIX/fix_ave_correlate_long.h +++ b/src/EXTRA-FIX/fix_ave_correlate_long.h @@ -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 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(); diff --git a/src/fix_ave_correlate.cpp b/src/fix_ave_correlate.cpp index 17182c8667..30d3e10e67 100644 --- a/src/fix_ave_correlate.cpp +++ b/src/fix_ave_correlate.cpp @@ -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; } diff --git a/src/fix_ave_correlate.h b/src/fix_ave_correlate.h index 517740a121..5dc49ed343 100644 --- a/src/fix_ave_correlate.h +++ b/src/fix_ave_correlate.h @@ -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 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;