Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80132634
fix_ave_chunk.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Aug 30, 01:49
Size
35 KB
Mime Type
text/x-c
Expires
Sun, Sep 1, 01:49 (2 d)
Engine
blob
Format
Raw Data
Handle
20301216
Attached To
rLAMMPS lammps
fix_ave_chunk.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include "fix_ave_chunk.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "modify.h"
#include "compute.h"
#include "compute_chunk_atom.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
enum{V,F,DENSITY_NUMBER,DENSITY_MASS,MASS,TEMPERATURE,COMPUTE,FIX,VARIABLE};
enum{SAMPLE,ALL};
enum{NOSCALE,ATOM};
enum{ONE,RUNNING,WINDOW};
#define INVOKED_PERATOM 8
/* ---------------------------------------------------------------------- */
FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
nvalues(0), nrepeat(0),
which(NULL), argindex(NULL), value2index(NULL), ids(NULL),
fp(NULL), idchunk(NULL), varatom(NULL),
count_one(NULL), count_many(NULL), count_sum(NULL),
values_one(NULL), values_many(NULL), values_sum(NULL),
count_total(NULL), count_list(NULL),
values_total(NULL), values_list(NULL)
{
if (narg < 7) error->all(FLERR,"Illegal fix ave/chunk command");
MPI_Comm_rank(world,&me);
nevery = force->inumeric(FLERR,arg[3]);
nrepeat = force->inumeric(FLERR,arg[4]);
nfreq = force->inumeric(FLERR,arg[5]);
int n = strlen(arg[6]) + 1;
idchunk = new char[n];
strcpy(idchunk,arg[6]);
global_freq = nfreq;
no_change_box = 1;
// expand args if any have wildcard character "*"
int expand = 0;
char **earg;
int nargnew = input->expand_args(narg-7,&arg[7],1,earg);
if (earg != &arg[7]) expand = 1;
arg = earg;
// parse values until one isn't recognized
which = new int[nargnew];
argindex = new int[nargnew];
ids = new char*[nargnew];
value2index = new int[nargnew];
int iarg = 0;
while (iarg < nargnew) {
ids[nvalues] = NULL;
if (strcmp(arg[iarg],"vx") == 0) {
which[nvalues] = V;
argindex[nvalues++] = 0;
} else if (strcmp(arg[iarg],"vy") == 0) {
which[nvalues] = V;
argindex[nvalues++] = 1;
} else if (strcmp(arg[iarg],"vz") == 0) {
which[nvalues] = V;
argindex[nvalues++] = 2;
} else if (strcmp(arg[iarg],"fx") == 0) {
which[nvalues] = F;
argindex[nvalues++] = 0;
} else if (strcmp(arg[iarg],"fy") == 0) {
which[nvalues] = F;
argindex[nvalues++] = 1;
} else if (strcmp(arg[iarg],"fz") == 0) {
which[nvalues] = F;
argindex[nvalues++] = 2;
} else if (strcmp(arg[iarg],"density/number") == 0) {
which[nvalues] = DENSITY_NUMBER;
argindex[nvalues++] = 0;
} else if (strcmp(arg[iarg],"density/mass") == 0) {
which[nvalues] = DENSITY_MASS;
argindex[nvalues++] = 0;
} else if (strcmp(arg[iarg],"mass") == 0) {
which[nvalues] = MASS;
argindex[nvalues++] = 0;
} else if (strcmp(arg[iarg],"temp") == 0) {
which[nvalues] = TEMPERATURE;
argindex[nvalues++] = 0;
} else if (strncmp(arg[iarg],"c_",2) == 0 ||
strncmp(arg[iarg],"f_",2) == 0 ||
strncmp(arg[iarg],"v_",2) == 0) {
if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
int n = strlen(arg[iarg]);
char *suffix = new char[n];
strcpy(suffix,&arg[iarg][2]);
char *ptr = strchr(suffix,'[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Illegal fix ave/chunk command");
argindex[nvalues] = atoi(ptr+1);
*ptr = '\0';
} else argindex[nvalues] = 0;
n = strlen(suffix) + 1;
ids[nvalues] = new char[n];
strcpy(ids[nvalues],suffix);
nvalues++;
delete [] suffix;
} else break;
iarg++;
}
if (nvalues == 0) error->all(FLERR,"No values in fix ave/chunk command");
// optional args
normflag = ALL;
scaleflag = ATOM;
ave = ONE;
nwindow = 0;
biasflag = 0;
id_bias = NULL;
adof = domain->dimension;
cdof = 0.0;
overwrite = 0;
format_user = NULL;
format = (char *) " %g";
char *title1 = NULL;
char *title2 = NULL;
char *title3 = NULL;
while (iarg < nargnew) {
if (strcmp(arg[iarg],"norm") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (strcmp(arg[iarg+1],"all") == 0) {
normflag = ALL;
scaleflag = ATOM;
} else if (strcmp(arg[iarg+1],"sample") == 0) {
normflag = SAMPLE;
scaleflag = ATOM;
} else if (strcmp(arg[iarg+1],"none") == 0) {
normflag = SAMPLE;
scaleflag = NOSCALE;
} else error->all(FLERR,"Illegal fix ave/chunk command");
iarg += 2;
} else if (strcmp(arg[iarg],"ave") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (strcmp(arg[iarg+1],"one") == 0) ave = ONE;
else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING;
else if (strcmp(arg[iarg+1],"window") == 0) ave = WINDOW;
else error->all(FLERR,"Illegal fix ave/chunk command");
if (ave == WINDOW) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
nwindow = force->inumeric(FLERR,arg[iarg+2]);
if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/chunk command");
}
iarg += 2;
if (ave == WINDOW) iarg++;
} else if (strcmp(arg[iarg],"bias") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/chunk command");
biasflag = 1;
int n = strlen(arg[iarg+1]) + 1;
id_bias = new char[n];
strcpy(id_bias,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"adof") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/chunk command");
adof = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"cdof") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal fix ave/chunk command");
cdof = force->numeric(FLERR,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"file") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
if (me == 0) {
fp = fopen(arg[iarg+1],"w");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open fix ave/chunk file %s",arg[iarg+1]);
error->one(FLERR,str);
}
}
iarg += 2;
} else if (strcmp(arg[iarg],"overwrite") == 0) {
overwrite = 1;
iarg += 1;
} else if (strcmp(arg[iarg],"format") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
delete [] format_user;
int n = strlen(arg[iarg+1]) + 2;
format_user = new char[n];
sprintf(format_user," %s",arg[iarg+1]);
format = format_user;
iarg += 2;
} else if (strcmp(arg[iarg],"title1") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
delete [] title1;
int n = strlen(arg[iarg+1]) + 1;
title1 = new char[n];
strcpy(title1,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"title2") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
delete [] title2;
int n = strlen(arg[iarg+1]) + 1;
title2 = new char[n];
strcpy(title2,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"title3") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
delete [] title3;
int n = strlen(arg[iarg+1]) + 1;
title3 = new char[n];
strcpy(title3,arg[iarg+1]);
iarg += 2;
} else error->all(FLERR,"Illegal fix ave/chunk command");
}
// setup and error check
if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
error->all(FLERR,"Illegal fix ave/chunk command");
if (nfreq % nevery || nrepeat*nevery > nfreq)
error->all(FLERR,"Illegal fix ave/chunk command");
if (ave != RUNNING && overwrite)
error->all(FLERR,"Illegal fix ave/chunk command");
if (biasflag) {
int i = modify->find_compute(id_bias);
if (i < 0)
error->all(FLERR,"Could not find compute ID for temperature bias");
tbias = modify->compute[i];
if (tbias->tempflag == 0)
error->all(FLERR,"Bias compute does not calculate temperature");
if (tbias->tempbias == 0)
error->all(FLERR,"Bias compute does not calculate a velocity bias");
}
for (int i = 0; i < nvalues; i++) {
if (which[i] == COMPUTE) {
int icompute = modify->find_compute(ids[i]);
if (icompute < 0)
error->all(FLERR,"Compute ID for fix ave/chunk does not exist");
if (modify->compute[icompute]->peratom_flag == 0)
error->all(FLERR,"Fix ave/chunk compute does not "
"calculate per-atom values");
if (argindex[i] == 0 &&
modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,"Fix ave/chunk compute does not "
"calculate a per-atom vector");
if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,"Fix ave/chunk compute does not "
"calculate a per-atom array");
if (argindex[i] &&
argindex[i] > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR,
"Fix ave/chunk compute vector is accessed out-of-range");
} else if (which[i] == FIX) {
int ifix = modify->find_fix(ids[i]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/chunk does not exist");
if (modify->fix[ifix]->peratom_flag == 0)
error->all(FLERR,
"Fix ave/chunk fix does not calculate per-atom values");
if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
error->all(FLERR,
"Fix ave/chunk fix does not calculate a per-atom vector");
if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
error->all(FLERR,
"Fix ave/chunk fix does not calculate a per-atom array");
if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
error->all(FLERR,"Fix ave/chunk fix vector is accessed out-of-range");
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/chunk does not exist");
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Fix ave/chunk variable is not atom-style variable");
}
}
// increment lock counter in compute chunk/atom
// only if nrepeat > 1 or ave = RUNNING/WINDOW,
// so that locking spans multiple timesteps
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for fix ave/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style,"chunk/atom") != 0)
error->all(FLERR,"Fix ave/chunk does not use chunk/atom compute");
if (nrepeat > 1 || ave == RUNNING || ave == WINDOW) cchunk->lockcount++;
lockforever = 0;
// print file comment lines
if (fp && me == 0) {
clearerr(fp);
if (title1) fprintf(fp,"%s\n",title1);
else fprintf(fp,"# Chunk-averaged data for fix %s and group %s\n",
id,arg[1]);
if (title2) fprintf(fp,"%s\n",title2);
else fprintf(fp,"# Timestep Number-of-chunks Total-count\n");
if (title3) fprintf(fp,"%s\n",title3);
else {
int compress = cchunk->compress;
int ncoord = cchunk->ncoord;
if (!compress) {
if (ncoord == 0) fprintf(fp,"# Chunk Ncount");
else if (ncoord == 1) fprintf(fp,"# Chunk Coord1 Ncount");
else if (ncoord == 2) fprintf(fp,"# Chunk Coord1 Coord2 Ncount");
else if (ncoord == 3)
fprintf(fp,"# Chunk Coord1 Coord2 Coord3 Ncount");
} else {
if (ncoord == 0) fprintf(fp,"# Chunk OrigID Ncount");
else if (ncoord == 1) fprintf(fp,"# Chunk OrigID Coord1 Ncount");
else if (ncoord == 2) fprintf(fp,"# Chunk OrigID Coord1 Coord2 Ncount");
else if (ncoord == 3)
fprintf(fp,"# Chunk OrigID Coord1 Coord2 Coord3 Ncount");
}
for (int i = 0; i < nvalues; i++) fprintf(fp," %s",earg[i]);
fprintf(fp,"\n");
}
if (ferror(fp))
error->one(FLERR,"Error writing file header");
filepos = ftell(fp);
}
delete [] title1;
delete [] title2;
delete [] title3;
// 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);
}
// this fix produces a global array
// size_array_rows is variable and set by allocate()
int compress = cchunk->compress;
int ncoord = cchunk->ncoord;
colextra = compress + ncoord;
array_flag = 1;
size_array_cols = colextra + 1 + nvalues;
size_array_rows_variable = 1;
extarray = 0;
// initializations
irepeat = 0;
iwindow = window_limit = 0;
normcount = 0;
maxvar = 0;
varatom = NULL;
count_one = count_many = count_sum = count_total = NULL;
count_list = NULL;
values_one = values_many = values_sum = values_total = NULL;
values_list = NULL;
maxchunk = 0;
nchunk = 1;
allocate();
// nvalid = next step on which end_of_step does something
// add nvalid to all computes that store invocation times
// since don't know a priori which are invoked by this fix
// once in end_of_step() can set timestep for ones actually invoked
nvalid_last = -1;
nvalid = nextvalid();
modify->addstep_compute_all(nvalid);
}
/* ---------------------------------------------------------------------- */
FixAveChunk::~FixAveChunk()
{
delete [] which;
delete [] argindex;
for (int i = 0; i < nvalues; i++) delete [] ids[i];
delete [] ids;
delete [] value2index;
if (fp && me == 0) fclose(fp);
memory->destroy(varatom);
memory->destroy(count_one);
memory->destroy(count_many);
memory->destroy(count_sum);
memory->destroy(count_total);
memory->destroy(count_list);
memory->destroy(values_one);
memory->destroy(values_many);
memory->destroy(values_sum);
memory->destroy(values_total);
memory->destroy(values_list);
// decrement lock counter in compute chunk/atom, it if still exists
if (nrepeat > 1 || ave == RUNNING || ave == WINDOW) {
int icompute = modify->find_compute(idchunk);
if (icompute >= 0) {
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (ave == RUNNING || ave == WINDOW) cchunk->unlock(this);
cchunk->lockcount--;
}
}
delete [] idchunk;
which = NULL;
argindex = NULL;
ids = NULL;
value2index = NULL;
fp = NULL;
varatom = NULL;
count_one = NULL;
count_many = NULL;
count_sum = NULL;
count_total = NULL;
count_list = NULL;
values_one = NULL;
values_many = NULL;
values_sum = NULL;
values_total = NULL;
values_list = NULL;
idchunk = NULL;
cchunk = NULL;
}
/* ---------------------------------------------------------------------- */
int FixAveChunk::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixAveChunk::init()
{
// set indices and check validity of all computes,fixes,variables
// check that fix frequency is acceptable
int icompute = modify->find_compute(idchunk);
if (icompute < 0)
error->all(FLERR,"Chunk/atom compute does not exist for fix ave/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (biasflag) {
int i = modify->find_compute(id_bias);
if (i < 0)
error->all(FLERR,"Could not find compute ID for temperature bias");
tbias = modify->compute[i];
}
for (int m = 0; m < nvalues; m++) {
if (which[m] == COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID for fix ave/chunk does not exist");
value2index[m] = icompute;
} else if (which[m] == FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID for fix ave/chunk does not exist");
value2index[m] = ifix;
if (nevery % modify->fix[ifix]->peratom_freq)
error->all(FLERR,
"Fix for fix ave/chunk not computed at compatible time");
} else if (which[m] == VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
error->all(FLERR,"Variable name for fix ave/chunk does not exist");
value2index[m] = ivariable;
} else value2index[m] = -1;
}
// need to reset nvalid if nvalid < ntimestep b/c minimize was performed
if (nvalid < update->ntimestep) {
irepeat = 0;
nvalid = nextvalid();
modify->addstep_compute_all(nvalid);
}
}
/* ----------------------------------------------------------------------
only does averaging if nvalid = current timestep
do not call setup_chunks(), even though fix ave/spatial called setup_bins()
b/c could cause nchunk to change if Nfreq epoch crosses 2 runs
does mean that if change_box is used between runs to change box size,
that nchunk may not track it
------------------------------------------------------------------------- */
void FixAveChunk::setup(int vflag)
{
end_of_step();
}
/* ---------------------------------------------------------------------- */
void FixAveChunk::end_of_step()
{
int i,j,m,n,index;
// skip if not step which requires doing something
// error check if timestep was reset in an invalid manner
bigint ntimestep = update->ntimestep;
if (ntimestep < nvalid_last || ntimestep > nvalid)
error->all(FLERR,"Invalid timestep reset for fix ave/chunk");
if (ntimestep != nvalid) return;
nvalid_last = nvalid;
// first sample within single Nfreq epoch
// zero out arrays that accumulate over many samples, but not across epochs
// invoke setup_chunks() to determine current nchunk
// re-allocate per-chunk arrays if needed
// invoke lock() in two cases:
// if nrepeat > 1: so nchunk cannot change until Nfreq epoch is over,
// will be unlocked on last repeat of this Nfreq
// if ave = RUNNING/WINDOW and not yet locked:
// set forever, will be unlocked in fix destructor
// wrap setup_chunks in clearstep/addstep b/c it may invoke computes
// both nevery and nfreq are future steps,
// since call below to cchunk->ichunk()
// does not re-invoke internal cchunk compute on this same step
if (irepeat == 0) {
if (cchunk->computeflag) modify->clearstep_compute();
nchunk = cchunk->setup_chunks();
if (cchunk->computeflag) {
modify->addstep_compute(ntimestep+nevery);
modify->addstep_compute(ntimestep+nfreq);
}
allocate();
if (nrepeat > 1 && ave == ONE)
cchunk->lock(this,ntimestep,ntimestep+(nrepeat-1)*nevery);
else if ((ave == RUNNING || ave == WINDOW) && !lockforever) {
cchunk->lock(this,update->ntimestep,-1);
lockforever = 1;
}
for (m = 0; m < nchunk; m++) {
count_many[m] = count_sum[m] = 0.0;
for (i = 0; i < nvalues; i++) values_many[m][i] = 0.0;
}
}
// zero out arrays for one sample
for (m = 0; m < nchunk; m++) {
count_one[m] = 0.0;
for (i = 0; i < nvalues; i++) values_one[m][i] = 0.0;
}
// compute chunk/atom assigns atoms to chunk IDs
// extract ichunk index vector from compute
// ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
// wrap compute_ichunk in clearstep/addstep b/c it may invoke computes
if (cchunk->computeflag) modify->clearstep_compute();
cchunk->compute_ichunk();
int *ichunk = cchunk->ichunk;
if (cchunk->computeflag) modify->addstep_compute(ntimestep+nevery);
// perform the computation for one sample
// count # of atoms in each bin
// accumulate results of attributes,computes,fixes,variables to local copy
// sum within each chunk, only include atoms in fix group
// compute/fix/variable may invoke computes so wrap with clear/add
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0)
count_one[ichunk[i]-1]++;
modify->clearstep_compute();
for (m = 0; m < nvalues; m++) {
n = value2index[m];
j = argindex[m];
// V,F adds velocities,forces to values
if (which[m] == V || which[m] == F) {
double **attribute;
if (which[m] == V) attribute = atom->v;
else attribute = atom->f;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
values_one[index][m] += attribute[i][j];
}
// DENSITY_NUMBER adds 1 to values
} else if (which[m] == DENSITY_NUMBER) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
values_one[index][m] += 1.0;
}
// DENSITY_MASS or MASS adds mass to values
} else if (which[m] == DENSITY_MASS || which[m] == MASS) {
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
values_one[index][m] += rmass[i];
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
values_one[index][m] += mass[type[i]];
}
}
// TEMPERATURE adds KE to values
// subtract and restore velocity bias if requested
} else if (which[m] == TEMPERATURE) {
if (biasflag) {
if (tbias->invoked_scalar != ntimestep) tbias->compute_scalar();
tbias->remove_bias_all();
}
double **v = atom->v;
int *type = atom->type;
double *mass = atom->mass;
double *rmass = atom->rmass;
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
values_one[index][m] +=
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
values_one[index][m] +=
(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[type[i]];
}
}
if (biasflag) tbias->restore_bias_all();
// COMPUTE adds its scalar or vector component to values
// invoke compute if not previously invoked
} else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[n];
if (!(compute->invoked_flag & INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= INVOKED_PERATOM;
}
double *vector = compute->vector_atom;
double **array = compute->array_atom;
int jm1 = j - 1;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
if (j == 0) values_one[index][m] += vector[i];
else values_one[index][m] += array[i][jm1];
}
// FIX adds its scalar or vector component to values
// access fix fields, guaranteed to be ready
} else if (which[m] == FIX) {
double *vector = modify->fix[n]->vector_atom;
double **array = modify->fix[n]->array_atom;
int jm1 = j - 1;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
if (j == 0) values_one[index][m] += vector[i];
else values_one[index][m] += array[i][jm1];
}
// VARIABLE adds its per-atom quantities to values
// evaluate atom-style variable
} else if (which[m] == VARIABLE) {
if (atom->nmax > maxvar) {
maxvar = atom->nmax;
memory->destroy(varatom);
memory->create(varatom,maxvar,"ave/chunk:varatom");
}
input->variable->compute_atom(n,igroup,varatom,1,0);
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit && ichunk[i] > 0) {
index = ichunk[i]-1;
values_one[index][m] += varatom[i];
}
}
}
// process the current sample
// if normflag = ALL, accumulate values,count separately to many
// if normflag = SAMPLE, one = value/count, accumulate one to many
// count is MPI summed here, value is MPI summed below across samples
// exception is TEMPERATURE: normalize by DOF
// exception is DENSITYs: no normalize by atom count
// exception is scaleflag = NOSCALE : no normalize by atom count
// check last so other options can take precedence
double mvv2e = force->mvv2e;
double boltz = force->boltz;
if (normflag == ALL) {
for (m = 0; m < nchunk; m++) {
count_many[m] += count_one[m];
for (j = 0; j < nvalues; j++)
values_many[m][j] += values_one[m][j];
}
} else if (normflag == SAMPLE) {
MPI_Allreduce(count_one,count_many,nchunk,MPI_DOUBLE,MPI_SUM,world);
for (m = 0; m < nchunk; m++) {
if (count_many[m] > 0.0)
for (j = 0; j < nvalues; j++) {
if (which[j] == TEMPERATURE)
values_many[m][j] += mvv2e*values_one[m][j] /
((cdof + adof*count_many[m]) * boltz);
else if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS ||
scaleflag == NOSCALE)
values_many[m][j] += values_one[m][j];
else
values_many[m][j] += values_one[m][j]/count_many[m];
}
count_sum[m] += count_many[m];
}
}
// done if irepeat < nrepeat
// else reset irepeat and nvalid
irepeat++;
if (irepeat < nrepeat) {
nvalid += nevery;
modify->addstep_compute(nvalid);
return;
}
irepeat = 0;
nvalid = ntimestep+nfreq - (nrepeat-1)*nevery;
modify->addstep_compute(nvalid);
// unlock compute chunk/atom at end of Nfreq epoch
// do not unlock if ave = RUNNING or WINDOW
if (nrepeat > 1 && ave == ONE) cchunk->unlock(this);
// time average across samples
// if normflag = ALL, final is total value / total count
// exception is TEMPERATURE: normalize by DOF for total count
// exception is DENSITYs: normalize by repeat, not total count
// exception is scaleflag == NOSCALE: normalize by repeat, not total count
// check last so other options can take precedence
// if normflag = SAMPLE, final is sum of ave / repeat
double repeat = nrepeat;
double mv2d = force->mv2d;
if (normflag == ALL) {
MPI_Allreduce(count_many,count_sum,nchunk,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&values_many[0][0],&values_sum[0][0],nchunk*nvalues,
MPI_DOUBLE,MPI_SUM,world);
for (m = 0; m < nchunk; m++) {
if (count_sum[m] > 0.0)
for (j = 0; j < nvalues; j++) {
if (which[j] == TEMPERATURE)
values_sum[m][j] *= mvv2e / ((cdof + adof*count_sum[m]) * boltz);
else if (which[j] == DENSITY_MASS)
values_sum[m][j] *= mv2d/repeat;
else if (which[j] == DENSITY_NUMBER || scaleflag == NOSCALE)
values_sum[m][j] /= repeat;
else values_sum[m][j] /= count_sum[m];
}
count_sum[m] /= repeat;
}
} else if (normflag == SAMPLE) {
MPI_Allreduce(&values_many[0][0],&values_sum[0][0],nchunk*nvalues,
MPI_DOUBLE,MPI_SUM,world);
for (m = 0; m < nchunk; m++) {
for (j = 0; j < nvalues; j++) values_sum[m][j] /= repeat;
count_sum[m] /= repeat;
}
}
// DENSITYs are additionally normalized by chunk volume
// use scalar or vector values for volume(s)
// if chunks are not spatial bins, chunk_volume_scalar = 1.0
for (j = 0; j < nvalues; j++)
if (which[j] == DENSITY_NUMBER || which[j] == DENSITY_MASS) {
if (cchunk->chunk_volume_vec) {
double *chunk_volume_vec = cchunk->chunk_volume_vec;
for (m = 0; m < nchunk; m++)
values_sum[m][j] /= chunk_volume_vec[m];
} else {
double chunk_volume_scalar = cchunk->chunk_volume_scalar;
for (m = 0; m < nchunk; m++)
values_sum[m][j] /= chunk_volume_scalar;
}
}
// if ave = ONE, only single Nfreq timestep value is needed
// if ave = RUNNING, combine with all previous Nfreq timestep values
// if ave = WINDOW, comine with nwindow most recent Nfreq timestep values
if (ave == ONE) {
for (m = 0; m < nchunk; m++) {
for (i = 0; i < nvalues; i++)
values_total[m][i] = values_sum[m][i];
count_total[m] = count_sum[m];
}
normcount = 1;
} else if (ave == RUNNING) {
for (m = 0; m < nchunk; m++) {
for (i = 0; i < nvalues; i++)
values_total[m][i] += values_sum[m][i];
count_total[m] += count_sum[m];
}
normcount++;
} else if (ave == WINDOW) {
for (m = 0; m < nchunk; m++) {
for (i = 0; i < nvalues; i++) {
values_total[m][i] += values_sum[m][i];
if (window_limit) values_total[m][i] -= values_list[iwindow][m][i];
values_list[iwindow][m][i] = values_sum[m][i];
}
count_total[m] += count_sum[m];
if (window_limit) count_total[m] -= count_list[iwindow][m];
count_list[iwindow][m] = count_sum[m];
}
iwindow++;
if (iwindow == nwindow) {
iwindow = 0;
window_limit = 1;
}
if (window_limit) normcount = nwindow;
else normcount = iwindow;
}
// output result to file
if (fp && me == 0) {
clearerr(fp);
if (overwrite) fseek(fp,filepos,SEEK_SET);
double count = 0.0;
for (m = 0; m < nchunk; m++) count += count_total[m];
fprintf(fp,BIGINT_FORMAT " %d %g\n",ntimestep,nchunk,count);
int compress = cchunk->compress;
int *chunkID = cchunk->chunkID;
int ncoord = cchunk->ncoord;
double **coord = cchunk->coord;
if (!compress) {
if (ncoord == 0) {
for (m = 0; m < nchunk; m++) {
fprintf(fp," %d %g",m+1,count_total[m]/normcount);
for (i = 0; i < nvalues; i++)
fprintf(fp,format,values_total[m][i]/normcount);
fprintf(fp,"\n");
}
} else if (ncoord == 1) {
for (m = 0; m < nchunk; m++) {
fprintf(fp," %d %g %g",m+1,coord[m][0],
count_total[m]/normcount);
for (i = 0; i < nvalues; i++)
fprintf(fp,format,values_total[m][i]/normcount);
fprintf(fp,"\n");
}
} else if (ncoord == 2) {
for (m = 0; m < nchunk; m++) {
fprintf(fp," %d %g %g %g",m+1,coord[m][0],coord[m][1],
count_total[m]/normcount);
for (i = 0; i < nvalues; i++)
fprintf(fp,format,values_total[m][i]/normcount);
fprintf(fp,"\n");
}
} else if (ncoord == 3) {
for (m = 0; m < nchunk; m++) {
fprintf(fp," %d %g %g %g %g",m+1,
coord[m][0],coord[m][1],coord[m][2],count_total[m]/normcount);
for (i = 0; i < nvalues; i++)
fprintf(fp,format,values_total[m][i]/normcount);
fprintf(fp,"\n");
}
}
} else {
int j;
if (ncoord == 0) {
for (m = 0; m < nchunk; m++) {
fprintf(fp," %d %d %g",m+1,chunkID[m],count_total[m]/normcount);
for (i = 0; i < nvalues; i++)
fprintf(fp,format,values_total[m][i]/normcount);
fprintf(fp,"\n");
}
} else if (ncoord == 1) {
for (m = 0; m < nchunk; m++) {
j = chunkID[m];
fprintf(fp," %d %d %g %g",m+1,j,coord[j-1][0],
count_total[m]/normcount);
for (i = 0; i < nvalues; i++)
fprintf(fp,format,values_total[m][i]/normcount);
fprintf(fp,"\n");
}
} else if (ncoord == 2) {
for (m = 0; m < nchunk; m++) {
j = chunkID[m];
fprintf(fp," %d %d %g %g %g",m+1,j,coord[j-1][0],coord[j-1][1],
count_total[m]/normcount);
for (i = 0; i < nvalues; i++)
fprintf(fp,format,values_total[m][i]/normcount);
fprintf(fp,"\n");
}
} else if (ncoord == 3) {
for (m = 0; m < nchunk; m++) {
j = chunkID[m];
fprintf(fp," %d %d %g %g %g %g",m+1,j,coord[j-1][0],
coord[j-1][1],coord[j-1][2],count_total[m]/normcount);
for (i = 0; i < nvalues; i++)
fprintf(fp,format,values_total[m][i]/normcount);
fprintf(fp,"\n");
}
}
}
if (ferror(fp))
error->one(FLERR,"Error writing averaged chunk data");
fflush(fp);
if (overwrite) {
long fileend = ftell(fp);
if (fileend > 0) ftruncate(fileno(fp),fileend);
}
}
}
/* ----------------------------------------------------------------------
allocate all per-chunk vectors
------------------------------------------------------------------------- */
void FixAveChunk::allocate()
{
size_array_rows = nchunk;
// reallocate chunk arrays if needed
if (nchunk > maxchunk) {
maxchunk = nchunk;
memory->grow(count_one,nchunk,"ave/chunk:count_one");
memory->grow(count_many,nchunk,"ave/chunk:count_many");
memory->grow(count_sum,nchunk,"ave/chunk:count_sum");
memory->grow(count_total,nchunk,"ave/chunk:count_total");
memory->grow(values_one,nchunk,nvalues,"ave/chunk:values_one");
memory->grow(values_many,nchunk,nvalues,"ave/chunk:values_many");
memory->grow(values_sum,nchunk,nvalues,"ave/chunk:values_sum");
memory->grow(values_total,nchunk,nvalues,"ave/chunk:values_total");
// only allocate count and values list for ave = WINDOW
if (ave == WINDOW) {
memory->create(count_list,nwindow,nchunk,"ave/chunk:count_list");
memory->create(values_list,nwindow,nchunk,nvalues,
"ave/chunk:values_list");
}
// reinitialize regrown count/values total since they accumulate
int i,m;
for (m = 0; m < nchunk; m++) {
for (i = 0; i < nvalues; i++) values_total[m][i] = 0.0;
count_total[m] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
return I,J array value
if I exceeds current nchunks, return 0.0 instead of generating an error
columns 1 to colextra = chunkID + ncoord
next column = count, remaining columns = Nvalues
------------------------------------------------------------------------- */
double FixAveChunk::compute_array(int i, int j)
{
if (values_total == NULL) return 0.0;
if (i >= nchunk) return 0.0;
if (j < colextra) {
if (cchunk->compress) {
if (j == 0) return (double) cchunk->chunkID[i];
return cchunk->coord[i][j-1];
} else return cchunk->coord[i][j];
}
j -= colextra + 1;
if (!normcount) return 0.0;
if (j < 0) return count_total[i]/normcount;
return values_total[i][j]/normcount;
}
/* ----------------------------------------------------------------------
calculate nvalid = next step on which end_of_step does something
can be this timestep if multiple of nfreq and nrepeat = 1
else backup from next multiple of nfreq
------------------------------------------------------------------------- */
bigint FixAveChunk::nextvalid()
{
bigint nvalid = (update->ntimestep/nfreq)*nfreq + nfreq;
if (nvalid-nfreq == update->ntimestep && nrepeat == 1)
nvalid = update->ntimestep;
else
nvalid -= (nrepeat-1)*nevery;
if (nvalid < update->ntimestep) nvalid += nfreq;
return nvalid;
}
/* ----------------------------------------------------------------------
memory usage of varatom and bins
------------------------------------------------------------------------- */
double FixAveChunk::memory_usage()
{
double bytes = maxvar * sizeof(double); // varatom
bytes += 4*maxchunk * sizeof(double); // count one,many,sum,total
bytes += nvalues*maxchunk * sizeof(double); // values one,many,sum,total
bytes += nwindow*maxchunk * sizeof(double); // count_list
bytes += nwindow*maxchunk*nvalues * sizeof(double); // values_list
return bytes;
}
Event Timeline
Log In to Comment