1 // clang-format off
2 /* ----------------------------------------------------------------------
3    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4    https://www.lammps.org/, Sandia National Laboratories
5    Steve Plimpton, sjplimp@sandia.gov
6 
7    Copyright (2003) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 
15 #include "fix_ave_chunk.h"
16 
17 #include "arg_info.h"
18 #include "atom.h"
19 #include "comm.h"
20 #include "compute.h"
21 #include "compute_chunk_atom.h"
22 #include "domain.h"
23 #include "error.h"
24 #include "force.h"
25 #include "input.h"
26 #include "memory.h"
27 #include "modify.h"
28 #include "update.h"
29 #include "variable.h"
30 
31 #include <cstring>
32 #include <unistd.h>
33 
34 using namespace LAMMPS_NS;
35 using namespace FixConst;
36 
37 enum{SCALAR,VECTOR};
38 enum{SAMPLE,ALL};
39 enum{NOSCALE,ATOM};
40 enum{ONE,RUNNING,WINDOW};
41 
42 
43 /* ---------------------------------------------------------------------- */
44 
FixAveChunk(LAMMPS * lmp,int narg,char ** arg)45 FixAveChunk::FixAveChunk(LAMMPS *lmp, int narg, char **arg) :
46   Fix(lmp, narg, arg),
47   nvalues(0), nrepeat(0),
48   which(nullptr), argindex(nullptr), value2index(nullptr), ids(nullptr),
49   fp(nullptr), idchunk(nullptr), varatom(nullptr),
50   count_one(nullptr), count_many(nullptr), count_sum(nullptr),
51   values_one(nullptr), values_many(nullptr), values_sum(nullptr),
52   count_total(nullptr), count_list(nullptr),
53   values_total(nullptr), values_list(nullptr)
54 {
55   if (narg < 7) error->all(FLERR,"Illegal fix ave/chunk command");
56 
57   nevery = utils::inumeric(FLERR,arg[3],false,lmp);
58   nrepeat = utils::inumeric(FLERR,arg[4],false,lmp);
59   nfreq = utils::inumeric(FLERR,arg[5],false,lmp);
60 
61   idchunk = utils::strdup(arg[6]);
62 
63   global_freq = nfreq;
64   no_change_box = 1;
65 
66   char * group = arg[1];
67 
68   // expand args if any have wildcard character "*"
69 
70   int expand = 0;
71   char **earg;
72   int nargnew = utils::expand_args(FLERR,narg-7,&arg[7],1,earg,lmp);
73 
74   if (earg != &arg[7]) expand = 1;
75   arg = earg;
76 
77   // parse values until one isn't recognized
78 
79   which = new int[nargnew];
80   argindex = new int[nargnew];
81   ids = new char*[nargnew];
82   value2index = new int[nargnew];
83   densityflag = 0;
84 
85   int iarg = 0;
86   while (iarg < nargnew) {
87 
88     ids[nvalues] = nullptr;
89 
90     if (strcmp(arg[iarg],"vx") == 0) {
91       which[nvalues] = ArgInfo::V;
92       argindex[nvalues++] = 0;
93     } else if (strcmp(arg[iarg],"vy") == 0) {
94       which[nvalues] = ArgInfo::V;
95       argindex[nvalues++] = 1;
96     } else if (strcmp(arg[iarg],"vz") == 0) {
97       which[nvalues] = ArgInfo::V;
98       argindex[nvalues++] = 2;
99 
100     } else if (strcmp(arg[iarg],"fx") == 0) {
101       which[nvalues] = ArgInfo::F;
102       argindex[nvalues++] = 0;
103     } else if (strcmp(arg[iarg],"fy") == 0) {
104       which[nvalues] = ArgInfo::F;
105       argindex[nvalues++] = 1;
106     } else if (strcmp(arg[iarg],"fz") == 0) {
107       which[nvalues] = ArgInfo::F;
108       argindex[nvalues++] = 2;
109 
110     } else if (strcmp(arg[iarg],"density/number") == 0) {
111       densityflag = 1;
112       which[nvalues] = ArgInfo::DENSITY_NUMBER;
113       argindex[nvalues++] = 0;
114     } else if (strcmp(arg[iarg],"density/mass") == 0) {
115       densityflag = 1;
116       which[nvalues] = ArgInfo::DENSITY_MASS;
117       argindex[nvalues++] = 0;
118     } else if (strcmp(arg[iarg],"mass") == 0) {
119       which[nvalues] = ArgInfo::MASS;
120       argindex[nvalues++] = 0;
121     } else if (strcmp(arg[iarg],"temp") == 0) {
122       which[nvalues] = ArgInfo::TEMPERATURE;
123       argindex[nvalues++] = 0;
124 
125     } else {
126       ArgInfo argi(arg[iarg]);
127 
128       if (argi.get_type() == ArgInfo::NONE) break;
129       if ((argi.get_type() == ArgInfo::UNKNOWN) || (argi.get_dim() > 1))
130         error->all(FLERR,"Invalid fix ave/chunk command");
131 
132       which[nvalues] = argi.get_type();
133       argindex[nvalues] = argi.get_index1();
134       ids[nvalues] = argi.copy_name();
135 
136       nvalues++;
137     }
138     iarg++;
139   }
140 
141   if (nvalues == 0) error->all(FLERR,"No values in fix ave/chunk command");
142 
143   // optional args
144 
145   normflag = ALL;
146   scaleflag = ATOM;
147   ave = ONE;
148   nwindow = 0;
149   biasflag = 0;
150   id_bias = nullptr;
151   adof = domain->dimension;
152   cdof = 0.0;
153   overwrite = 0;
154   format_user = nullptr;
155   format = (char *) " %g";
156   char *title1 = nullptr;
157   char *title2 = nullptr;
158   char *title3 = nullptr;
159 
160   while (iarg < nargnew) {
161     if (strcmp(arg[iarg],"norm") == 0) {
162       if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
163       if (strcmp(arg[iarg+1],"all") == 0) {
164         normflag = ALL;
165         scaleflag = ATOM;
166       } else if (strcmp(arg[iarg+1],"sample") == 0) {
167         normflag = SAMPLE;
168         scaleflag = ATOM;
169       } else if (strcmp(arg[iarg+1],"none") == 0) {
170         normflag = SAMPLE;
171         scaleflag = NOSCALE;
172       } else error->all(FLERR,"Illegal fix ave/chunk command");
173       iarg += 2;
174     } else if (strcmp(arg[iarg],"ave") == 0) {
175       if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
176       if (strcmp(arg[iarg+1],"one") == 0) ave = ONE;
177       else if (strcmp(arg[iarg+1],"running") == 0) ave = RUNNING;
178       else if (strcmp(arg[iarg+1],"window") == 0) ave = WINDOW;
179       else error->all(FLERR,"Illegal fix ave/chunk command");
180       if (ave == WINDOW) {
181         if (iarg+3 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
182         nwindow = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
183         if (nwindow <= 0) error->all(FLERR,"Illegal fix ave/chunk command");
184       }
185       iarg += 2;
186       if (ave == WINDOW) iarg++;
187 
188     } else if (strcmp(arg[iarg],"bias") == 0) {
189       if (iarg+2 > narg)
190         error->all(FLERR,"Illegal fix ave/chunk command");
191       biasflag = 1;
192       id_bias = utils::strdup(arg[iarg+1]);
193       iarg += 2;
194     } else if (strcmp(arg[iarg],"adof") == 0) {
195       if (iarg+2 > narg)
196         error->all(FLERR,"Illegal fix ave/chunk command");
197       adof = utils::numeric(FLERR,arg[iarg+1],false,lmp);
198       iarg += 2;
199     } else if (strcmp(arg[iarg],"cdof") == 0) {
200       if (iarg+2 > narg)
201         error->all(FLERR,"Illegal fix ave/chunk command");
202       cdof = utils::numeric(FLERR,arg[iarg+1],false,lmp);
203       iarg += 2;
204 
205     } else if (strcmp(arg[iarg],"file") == 0) {
206       if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
207       if (comm->me == 0) {
208         fp = fopen(arg[iarg+1],"w");
209         if (fp == nullptr)
210           error->one(FLERR,"Cannot open fix ave/chunk file {}: {}",
211                                        arg[iarg+1], utils::getsyserror());
212       }
213       iarg += 2;
214     } else if (strcmp(arg[iarg],"overwrite") == 0) {
215       overwrite = 1;
216       iarg += 1;
217     } else if (strcmp(arg[iarg],"format") == 0) {
218       if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
219       delete [] format_user;
220       format_user = utils::strdup(arg[iarg+1]);
221       format = format_user;
222       iarg += 2;
223     } else if (strcmp(arg[iarg],"title1") == 0) {
224       if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
225       delete [] title1;
226       title1 = utils::strdup(arg[iarg+1]);
227       iarg += 2;
228     } else if (strcmp(arg[iarg],"title2") == 0) {
229       if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
230       delete [] title2;
231       title2 = utils::strdup(arg[iarg+1]);
232       iarg += 2;
233     } else if (strcmp(arg[iarg],"title3") == 0) {
234       if (iarg+2 > narg) error->all(FLERR,"Illegal fix ave/chunk command");
235       delete [] title3;
236       title3 = utils::strdup(arg[iarg+1]);
237       iarg += 2;
238     } else error->all(FLERR,"Illegal fix ave/chunk command");
239   }
240 
241   // setup and error check
242 
243   if (nevery <= 0 || nrepeat <= 0 || nfreq <= 0)
244     error->all(FLERR,"Illegal fix ave/chunk command");
245   if (nfreq % nevery || nrepeat*nevery > nfreq)
246     error->all(FLERR,"Illegal fix ave/chunk command");
247   if (ave != RUNNING && overwrite)
248     error->all(FLERR,"Illegal fix ave/chunk command");
249 
250   if (biasflag) {
251     int i = modify->find_compute(id_bias);
252     if (i < 0)
253       error->all(FLERR,"Could not find compute ID for temperature bias");
254     tbias = modify->compute[i];
255     if (tbias->tempflag == 0)
256       error->all(FLERR,"Bias compute does not calculate temperature");
257     if (tbias->tempbias == 0)
258       error->all(FLERR,"Bias compute does not calculate a velocity bias");
259   }
260 
261   for (int i = 0; i < nvalues; i++) {
262     if (which[i] == ArgInfo::COMPUTE) {
263       int icompute = modify->find_compute(ids[i]);
264       if (icompute < 0)
265         error->all(FLERR,"Compute ID for fix ave/chunk does not exist");
266       if (modify->compute[icompute]->peratom_flag == 0)
267         error->all(FLERR,"Fix ave/chunk compute does not "
268                    "calculate per-atom values");
269       if (argindex[i] == 0 &&
270           modify->compute[icompute]->size_peratom_cols != 0)
271         error->all(FLERR,"Fix ave/chunk compute does not "
272                    "calculate a per-atom vector");
273       if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
274         error->all(FLERR,"Fix ave/chunk compute does not "
275                    "calculate a per-atom array");
276       if (argindex[i] &&
277           argindex[i] > modify->compute[icompute]->size_peratom_cols)
278         error->all(FLERR,
279                    "Fix ave/chunk compute vector is accessed out-of-range");
280 
281     } else if (which[i] == ArgInfo::FIX) {
282       int ifix = modify->find_fix(ids[i]);
283       if (ifix < 0)
284         error->all(FLERR,"Fix ID for fix ave/chunk does not exist");
285       if (modify->fix[ifix]->peratom_flag == 0)
286         error->all(FLERR,
287                    "Fix ave/chunk fix does not calculate per-atom values");
288       if (argindex[i] == 0 && modify->fix[ifix]->size_peratom_cols != 0)
289         error->all(FLERR,
290                    "Fix ave/chunk fix does not calculate a per-atom vector");
291       if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
292         error->all(FLERR,
293                    "Fix ave/chunk fix does not calculate a per-atom array");
294       if (argindex[i] && argindex[i] > modify->fix[ifix]->size_peratom_cols)
295         error->all(FLERR,"Fix ave/chunk fix vector is accessed out-of-range");
296     } else if (which[i] == ArgInfo::VARIABLE) {
297       int ivariable = input->variable->find(ids[i]);
298       if (ivariable < 0)
299         error->all(FLERR,"Variable name for fix ave/chunk does not exist");
300       if (input->variable->atomstyle(ivariable) == 0)
301         error->all(FLERR,"Fix ave/chunk variable is not atom-style variable");
302     }
303   }
304 
305   // increment lock counter in compute chunk/atom
306   // only if nrepeat > 1 or ave = RUNNING/WINDOW,
307   //   so that locking spans multiple timesteps
308 
309   int icompute = modify->find_compute(idchunk);
310   if (icompute < 0)
311     error->all(FLERR,"Chunk/atom compute does not exist for fix ave/chunk");
312   cchunk = (ComputeChunkAtom *) modify->compute[icompute];
313   if (strcmp(cchunk->style,"chunk/atom") != 0)
314     error->all(FLERR,"Fix ave/chunk does not use chunk/atom compute");
315 
316   if (nrepeat > 1 || ave == RUNNING || ave == WINDOW) cchunk->lockcount++;
317   lockforever = 0;
318 
319   // print file comment lines
320 
321   if (fp && comm->me == 0) {
322     clearerr(fp);
323     if (title1) fprintf(fp,"%s\n",title1);
324     else fprintf(fp,"# Chunk-averaged data for fix %s and group %s\n",
325                  id, group);
326     if (title2) fprintf(fp,"%s\n",title2);
327     else fprintf(fp,"# Timestep Number-of-chunks Total-count\n");
328     if (title3) fprintf(fp,"%s\n",title3);
329     else {
330       int compress = cchunk->compress;
331       int ncoord = cchunk->ncoord;
332       if (!compress) {
333         if (ncoord == 0) fprintf(fp,"# Chunk Ncount");
334         else if (ncoord == 1) fprintf(fp,"# Chunk Coord1 Ncount");
335         else if (ncoord == 2) fprintf(fp,"# Chunk Coord1 Coord2 Ncount");
336         else if (ncoord == 3)
337           fprintf(fp,"# Chunk Coord1 Coord2 Coord3 Ncount");
338       } else {
339         if (ncoord == 0) fprintf(fp,"# Chunk OrigID Ncount");
340         else if (ncoord == 1) fprintf(fp,"# Chunk OrigID Coord1 Ncount");
341         else if (ncoord == 2) fprintf(fp,"# Chunk OrigID Coord1 Coord2 Ncount");
342         else if (ncoord == 3)
343           fprintf(fp,"# Chunk OrigID Coord1 Coord2 Coord3 Ncount");
344       }
345       for (int i = 0; i < nvalues; i++) fprintf(fp," %s",earg[i]);
346       fprintf(fp,"\n");
347     }
348     if (ferror(fp))
349       error->one(FLERR,"Error writing file header");
350 
351     filepos = ftell(fp);
352   }
353 
354   delete [] title1;
355   delete [] title2;
356   delete [] title3;
357 
358   // if wildcard expansion occurred, free earg memory from expand_args()
359   // wait to do this until after file comment lines are printed
360 
361   if (expand) {
362     for (int i = 0; i < nargnew; i++) delete [] earg[i];
363     memory->sfree(earg);
364   }
365 
366   // this fix produces a global array
367   // size_array_rows is variable and set by allocate()
368 
369   int compress = cchunk->compress;
370   int ncoord = cchunk->ncoord;
371   colextra = compress + ncoord;
372 
373   array_flag = 1;
374   size_array_cols = colextra + 1 + nvalues;
375   size_array_rows_variable = 1;
376   extarray = 0;
377 
378   // initializations
379 
380   irepeat = 0;
381   iwindow = window_limit = 0;
382   normcount = 0;
383 
384   maxvar = 0;
385   varatom = nullptr;
386 
387   count_one = count_many = count_sum = count_total = nullptr;
388   count_list = nullptr;
389   values_one = values_many = values_sum = values_total = nullptr;
390   values_list = nullptr;
391 
392   maxchunk = 0;
393   nchunk = 1;
394   allocate();
395 
396   // nvalid = next step on which end_of_step does something
397   // add nvalid to all computes that store invocation times
398   // since don't know a priori which are invoked by this fix
399   // once in end_of_step() can set timestep for ones actually invoked
400 
401   nvalid_last = -1;
402   nvalid = nextvalid();
403   modify->addstep_compute_all(nvalid);
404 }
405 
406 /* ---------------------------------------------------------------------- */
407 
~FixAveChunk()408 FixAveChunk::~FixAveChunk()
409 {
410   delete [] which;
411   delete [] argindex;
412   for (int i = 0; i < nvalues; i++) delete [] ids[i];
413   delete [] ids;
414   delete [] value2index;
415 
416   if (fp && comm->me == 0) fclose(fp);
417 
418   memory->destroy(varatom);
419   memory->destroy(count_one);
420   memory->destroy(count_many);
421   memory->destroy(count_sum);
422   memory->destroy(count_total);
423   memory->destroy(count_list);
424   memory->destroy(values_one);
425   memory->destroy(values_many);
426   memory->destroy(values_sum);
427   memory->destroy(values_total);
428   memory->destroy(values_list);
429 
430   // decrement lock counter in compute chunk/atom, it if still exists
431 
432   if (nrepeat > 1 || ave == RUNNING || ave == WINDOW) {
433     int icompute = modify->find_compute(idchunk);
434     if (icompute >= 0) {
435       cchunk = (ComputeChunkAtom *) modify->compute[icompute];
436       if (ave == RUNNING || ave == WINDOW) cchunk->unlock(this);
437       cchunk->lockcount--;
438     }
439   }
440 
441   delete [] idchunk;
442   which = nullptr;
443   argindex = nullptr;
444   ids = nullptr;
445   value2index = nullptr;
446   fp = nullptr;
447   varatom = nullptr;
448   count_one = nullptr;
449   count_many = nullptr;
450   count_sum = nullptr;
451   count_total = nullptr;
452   count_list = nullptr;
453   values_one = nullptr;
454   values_many = nullptr;
455   values_sum = nullptr;
456   values_total = nullptr;
457   values_list = nullptr;
458   idchunk = nullptr;
459   cchunk = nullptr;
460 }
461 
462 /* ---------------------------------------------------------------------- */
463 
setmask()464 int FixAveChunk::setmask()
465 {
466   int mask = 0;
467   mask |= END_OF_STEP;
468   return mask;
469 }
470 
471 /* ---------------------------------------------------------------------- */
472 
init()473 void FixAveChunk::init()
474 {
475   // set indices and check validity of all computes,fixes,variables
476   // check that fix frequency is acceptable
477 
478   int icompute = modify->find_compute(idchunk);
479   if (icompute < 0)
480     error->all(FLERR,"Chunk/atom compute does not exist for fix ave/chunk");
481   cchunk = (ComputeChunkAtom *) modify->compute[icompute];
482 
483   if (biasflag) {
484     int i = modify->find_compute(id_bias);
485     if (i < 0)
486       error->all(FLERR,"Could not find compute ID for temperature bias");
487     tbias = modify->compute[i];
488   }
489 
490   for (int m = 0; m < nvalues; m++) {
491     if (which[m] == ArgInfo::COMPUTE) {
492       icompute = modify->find_compute(ids[m]);
493       if (icompute < 0)
494         error->all(FLERR,"Compute ID for fix ave/chunk does not exist");
495       value2index[m] = icompute;
496 
497     } else if (which[m] == ArgInfo::FIX) {
498       int ifix = modify->find_fix(ids[m]);
499       if (ifix < 0)
500         error->all(FLERR,"Fix ID for fix ave/chunk does not exist");
501       value2index[m] = ifix;
502 
503       if (nevery % modify->fix[ifix]->peratom_freq)
504         error->all(FLERR,
505                    "Fix for fix ave/chunk not computed at compatible time");
506 
507     } else if (which[m] == ArgInfo::VARIABLE) {
508       int ivariable = input->variable->find(ids[m]);
509       if (ivariable < 0)
510         error->all(FLERR,"Variable name for fix ave/chunk does not exist");
511       value2index[m] = ivariable;
512 
513     } else value2index[m] = -1;
514   }
515 
516   // need to reset nvalid if nvalid < ntimestep b/c minimize was performed
517 
518   if (nvalid < update->ntimestep) {
519     irepeat = 0;
520     nvalid = nextvalid();
521     modify->addstep_compute_all(nvalid);
522   }
523 }
524 
525 /* ----------------------------------------------------------------------
526    only does averaging if nvalid = current timestep
527    do not call setup_chunks(), even though fix ave/spatial called setup_bins()
528    b/c could cause nchunk to change if Nfreq epoch crosses 2 runs
529    does mean that if change_box is used between runs to change box size,
530      that nchunk may not track it
531 ------------------------------------------------------------------------- */
532 
setup(int)533 void FixAveChunk::setup(int /*vflag*/)
534 {
535   end_of_step();
536 }
537 
538 /* ---------------------------------------------------------------------- */
539 
end_of_step()540 void FixAveChunk::end_of_step()
541 {
542   int i,j,m,n,index;
543 
544   // skip if not step which requires doing something
545   // error check if timestep was reset in an invalid manner
546 
547   bigint ntimestep = update->ntimestep;
548   if (ntimestep < nvalid_last || ntimestep > nvalid)
549     error->all(FLERR,"Invalid timestep reset for fix ave/chunk");
550   if (ntimestep != nvalid) return;
551   nvalid_last = nvalid;
552 
553   // first sample within single Nfreq epoch
554   // zero out arrays that accumulate over many samples, but not across epochs
555   // invoke setup_chunks() to determine current nchunk
556   //   re-allocate per-chunk arrays if needed
557   // invoke lock() in two cases:
558   //   if nrepeat > 1: so nchunk cannot change until Nfreq epoch is over,
559   //     will be unlocked on last repeat of this Nfreq
560   //   if ave = RUNNING/WINDOW and not yet locked:
561   //     set forever, will be unlocked in fix destructor
562   // wrap setup_chunks in clearstep/addstep b/c it may invoke computes
563   //   both nevery and nfreq are future steps,
564   //   since call below to cchunk->ichunk()
565   //     does not re-invoke internal cchunk compute on this same step
566 
567   if (irepeat == 0) {
568     if (cchunk->computeflag) modify->clearstep_compute();
569     nchunk = cchunk->setup_chunks();
570     if (cchunk->computeflag) {
571       modify->addstep_compute(ntimestep+nevery);
572       modify->addstep_compute(ntimestep+nfreq);
573     }
574     allocate();
575     if (nrepeat > 1 && ave == ONE)
576       cchunk->lock(this,ntimestep,ntimestep+((bigint)nrepeat-1)*nevery);
577     else if ((ave == RUNNING || ave == WINDOW) && !lockforever) {
578       cchunk->lock(this,update->ntimestep,-1);
579       lockforever = 1;
580     }
581     for (m = 0; m < nchunk; m++) {
582       count_many[m] = count_sum[m] = 0.0;
583       for (i = 0; i < nvalues; i++) values_many[m][i] = 0.0;
584     }
585 
586   // if any DENSITY requested, invoke setup_chunks() on each sampling step
587   // nchunk will not change but bin volumes might, e.g. for NPT simulation
588 
589   } else if (densityflag) {
590     cchunk->setup_chunks();
591   }
592 
593   // zero out arrays for one sample
594 
595   for (m = 0; m < nchunk; m++) {
596     count_one[m] = 0.0;
597     for (i = 0; i < nvalues; i++) values_one[m][i] = 0.0;
598   }
599 
600   // compute chunk/atom assigns atoms to chunk IDs
601   // extract ichunk index vector from compute
602   // ichunk = 1 to Nchunk for included atoms, 0 for excluded atoms
603   // wrap compute_ichunk in clearstep/addstep b/c it may invoke computes
604 
605   if (cchunk->computeflag) modify->clearstep_compute();
606 
607   cchunk->compute_ichunk();
608   int *ichunk = cchunk->ichunk;
609 
610   if (cchunk->computeflag) modify->addstep_compute(ntimestep+nevery);
611 
612   // perform the computation for one sample
613   // count # of atoms in each bin
614   // accumulate results of attributes,computes,fixes,variables to local copy
615   // sum within each chunk, only include atoms in fix group
616   // compute/fix/variable may invoke computes so wrap with clear/add
617 
618   int *mask = atom->mask;
619   int nlocal = atom->nlocal;
620 
621   for (i = 0; i < nlocal; i++)
622     if (mask[i] & groupbit && ichunk[i] > 0)
623       count_one[ichunk[i]-1]++;
624 
625   modify->clearstep_compute();
626 
627   for (m = 0; m < nvalues; m++) {
628     n = value2index[m];
629     j = argindex[m];
630 
631     // V,F adds velocities,forces to values
632 
633     if (which[m] == ArgInfo::V || which[m] == ArgInfo::F) {
634       double **attribute;
635       if (which[m] == ArgInfo::V) attribute = atom->v;
636       else attribute = atom->f;
637 
638       for (i = 0; i < nlocal; i++)
639         if (mask[i] & groupbit && ichunk[i] > 0) {
640           index = ichunk[i]-1;
641           values_one[index][m] += attribute[i][j];
642         }
643 
644     // DENSITY_NUMBER adds 1 to values
645 
646     } else if (which[m] == ArgInfo::DENSITY_NUMBER) {
647 
648       for (i = 0; i < nlocal; i++)
649         if (mask[i] & groupbit && ichunk[i] > 0) {
650           index = ichunk[i]-1;
651           values_one[index][m] += 1.0;
652         }
653 
654     // DENSITY_MASS or MASS adds mass to values
655 
656     } else if ((which[m] == ArgInfo::DENSITY_MASS)
657                || (which[m] == ArgInfo::MASS)) {
658       int *type = atom->type;
659       double *mass = atom->mass;
660       double *rmass = atom->rmass;
661 
662       if (rmass) {
663         for (i = 0; i < nlocal; i++)
664           if (mask[i] & groupbit && ichunk[i] > 0) {
665             index = ichunk[i]-1;
666             values_one[index][m] += rmass[i];
667           }
668       } else {
669         for (i = 0; i < nlocal; i++)
670           if (mask[i] & groupbit && ichunk[i] > 0) {
671             index = ichunk[i]-1;
672             values_one[index][m] += mass[type[i]];
673           }
674       }
675 
676     // TEMPERATURE adds KE to values
677     // subtract and restore velocity bias if requested
678 
679     } else if (which[m] == ArgInfo::TEMPERATURE) {
680 
681       if (biasflag) {
682         if (tbias->invoked_scalar != ntimestep) tbias->compute_scalar();
683         tbias->remove_bias_all();
684       }
685 
686       double **v = atom->v;
687       int *type = atom->type;
688       double *mass = atom->mass;
689       double *rmass = atom->rmass;
690 
691       if (rmass) {
692         for (i = 0; i < nlocal; i++)
693           if (mask[i] & groupbit && ichunk[i] > 0) {
694             index = ichunk[i]-1;
695             values_one[index][m] +=
696               (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * rmass[i];
697           }
698       } else {
699         for (i = 0; i < nlocal; i++)
700           if (mask[i] & groupbit && ichunk[i] > 0) {
701             index = ichunk[i]-1;
702             values_one[index][m] +=
703               (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
704               mass[type[i]];
705           }
706       }
707 
708       if (biasflag) tbias->restore_bias_all();
709 
710     // COMPUTE adds its scalar or vector component to values
711     // invoke compute if not previously invoked
712 
713     } else if (which[m] == ArgInfo::COMPUTE) {
714       Compute *compute = modify->compute[n];
715       if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
716         compute->compute_peratom();
717         compute->invoked_flag |= Compute::INVOKED_PERATOM;
718       }
719       double *vector = compute->vector_atom;
720       double **array = compute->array_atom;
721       int jm1 = j - 1;
722 
723       for (i = 0; i < nlocal; i++)
724         if (mask[i] & groupbit && ichunk[i] > 0) {
725           index = ichunk[i]-1;
726           if (j == 0) values_one[index][m] += vector[i];
727           else values_one[index][m] += array[i][jm1];
728         }
729 
730     // FIX adds its scalar or vector component to values
731     // access fix fields, guaranteed to be ready
732 
733     } else if (which[m] == ArgInfo::FIX) {
734       double *vector = modify->fix[n]->vector_atom;
735       double **array = modify->fix[n]->array_atom;
736       int jm1 = j - 1;
737 
738       for (i = 0; i < nlocal; i++)
739         if (mask[i] & groupbit && ichunk[i] > 0) {
740           index = ichunk[i]-1;
741           if (j == 0) values_one[index][m] += vector[i];
742           else values_one[index][m] += array[i][jm1];
743         }
744 
745     // VARIABLE adds its per-atom quantities to values
746     // evaluate atom-style variable
747 
748     } else if (which[m] == ArgInfo::VARIABLE) {
749       if (atom->nmax > maxvar) {
750         maxvar = atom->nmax;
751         memory->destroy(varatom);
752         memory->create(varatom,maxvar,"ave/chunk:varatom");
753       }
754 
755       input->variable->compute_atom(n,igroup,varatom,1,0);
756 
757       for (i = 0; i < nlocal; i++)
758         if (mask[i] & groupbit && ichunk[i] > 0) {
759           index = ichunk[i]-1;
760           values_one[index][m] += varatom[i];
761         }
762     }
763   }
764 
765   // process the current sample
766   // if normflag = ALL, accumulate values,count separately to many
767   // if normflag = SAMPLE, one = value/count, accumulate one to many
768   //   count is MPI summed here, value is MPI summed below across samples
769   //   exception is TEMPERATURE: normalize by DOF
770   //   exception is DENSITY_NUMBER:
771   //     normalize by bin volume, not by atom count
772   //   exception is DENSITY_MASS:
773   //     scale by mv2d, normalize by bin volume, not by atom count
774   //   exception is scaleflag = NOSCALE (norm = NONE):
775   //     no normalize by atom count
776   //     check last so other options can take precedence
777 
778   double mvv2e = force->mvv2e;
779   double mv2d = force->mv2d;
780   double boltz = force->boltz;
781 
782   if (normflag == ALL) {
783     for (m = 0; m < nchunk; m++) {
784       count_many[m] += count_one[m];
785       for (j = 0; j < nvalues; j++)
786         values_many[m][j] += values_one[m][j];
787     }
788   } else if (normflag == SAMPLE) {
789     MPI_Allreduce(count_one,count_many,nchunk,MPI_DOUBLE,MPI_SUM,world);
790 
791     if (cchunk->chunk_volume_vec) {
792       volflag = VECTOR;
793       chunk_volume_vec = cchunk->chunk_volume_vec;
794     } else {
795       volflag = SCALAR;
796       chunk_volume_scalar = cchunk->chunk_volume_scalar;
797     }
798 
799     for (m = 0; m < nchunk; m++) {
800       if (count_many[m] > 0.0)
801         for (j = 0; j < nvalues; j++) {
802           if (which[j] == ArgInfo::TEMPERATURE) {
803             values_many[m][j] += mvv2e*values_one[m][j] /
804               ((cdof + adof*count_many[m]) * boltz);
805           } else if (which[j] == ArgInfo::DENSITY_NUMBER) {
806             if (volflag == SCALAR) values_one[m][j] /= chunk_volume_scalar;
807             else values_one[m][j] /= chunk_volume_vec[m];
808             values_many[m][j] += values_one[m][j];
809           } else if (which[j] == ArgInfo::DENSITY_MASS) {
810             if (volflag == SCALAR) values_one[m][j] /= chunk_volume_scalar;
811             else values_one[m][j] /= chunk_volume_vec[m];
812             values_many[m][j] += mv2d*values_one[m][j];
813           } else if (scaleflag == NOSCALE) {
814             values_many[m][j] += values_one[m][j];
815           } else {
816             values_many[m][j] += values_one[m][j]/count_many[m];
817           }
818         }
819       count_sum[m] += count_many[m];
820     }
821   }
822 
823   // done if irepeat < nrepeat
824   // else reset irepeat and nvalid
825 
826   irepeat++;
827   if (irepeat < nrepeat) {
828     nvalid += nevery;
829     modify->addstep_compute(nvalid);
830     return;
831   }
832 
833   irepeat = 0;
834   nvalid = ntimestep+nfreq - ((bigint)nrepeat-1)*nevery;
835   modify->addstep_compute(nvalid);
836 
837   // unlock compute chunk/atom at end of Nfreq epoch
838   // do not unlock if ave = RUNNING or WINDOW
839 
840   if (nrepeat > 1 && ave == ONE) cchunk->unlock(this);
841 
842   // time average across samples
843   // if normflag = ALL, final is total value / total count
844   //   exception is TEMPERATURE: normalize by DOF for total count
845   //   exception is DENSITY_NUMBER:
846   //     normalize by final bin_volume and repeat, not by total count
847   //   exception is DENSITY_MASS:
848   //     scale by mv2d, normalize by bin volume and repeat, not by total count
849   //   exception is scaleflag == NOSCALE:
850   //     normalize by repeat, not by total count
851   //     check last so other options can take precedence
852   // if normflag = SAMPLE, final is sum of ave / repeat
853 
854   double repeat = nrepeat;
855 
856   if (normflag == ALL) {
857     MPI_Allreduce(count_many,count_sum,nchunk,MPI_DOUBLE,MPI_SUM,world);
858     MPI_Allreduce(&values_many[0][0],&values_sum[0][0],nchunk*nvalues,
859                   MPI_DOUBLE,MPI_SUM,world);
860 
861     if (cchunk->chunk_volume_vec) {
862       volflag = VECTOR;
863       chunk_volume_vec = cchunk->chunk_volume_vec;
864     } else {
865       volflag = SCALAR;
866       chunk_volume_scalar = cchunk->chunk_volume_scalar;
867     }
868 
869     for (m = 0; m < nchunk; m++) {
870       if (count_sum[m] > 0.0)
871         for (j = 0; j < nvalues; j++) {
872           if (which[j] == ArgInfo::TEMPERATURE) {
873             values_sum[m][j] *= mvv2e/((repeat*cdof + adof*count_sum[m])*boltz);
874           } else if (which[j] == ArgInfo::DENSITY_NUMBER) {
875             if (volflag == SCALAR) values_sum[m][j] /= chunk_volume_scalar;
876             else values_sum[m][j] /= chunk_volume_vec[m];
877             values_sum[m][j] /= repeat;
878           } else if (which[j] == ArgInfo::DENSITY_MASS) {
879             if (volflag == SCALAR) values_sum[m][j] /= chunk_volume_scalar;
880             else values_sum[m][j] /= chunk_volume_vec[m];
881             values_sum[m][j] *= mv2d/repeat;
882           } else if (scaleflag == NOSCALE) {
883             values_sum[m][j] /= repeat;
884           } else {
885             values_sum[m][j] /= count_sum[m];
886           }
887         }
888       count_sum[m] /= repeat;
889     }
890   } else if (normflag == SAMPLE) {
891     MPI_Allreduce(&values_many[0][0],&values_sum[0][0],nchunk*nvalues,
892                   MPI_DOUBLE,MPI_SUM,world);
893     for (m = 0; m < nchunk; m++) {
894       for (j = 0; j < nvalues; j++) values_sum[m][j] /= repeat;
895       count_sum[m] /= repeat;
896     }
897   }
898 
899   // if ave = ONE, only single Nfreq timestep value is needed
900   // if ave = RUNNING, combine with all previous Nfreq timestep values
901   // if ave = WINDOW, comine with nwindow most recent Nfreq timestep values
902 
903   if (ave == ONE) {
904     for (m = 0; m < nchunk; m++) {
905       for (i = 0; i < nvalues; i++)
906         values_total[m][i] = values_sum[m][i];
907       count_total[m] = count_sum[m];
908     }
909     normcount = 1;
910 
911   } else if (ave == RUNNING) {
912     for (m = 0; m < nchunk; m++) {
913       for (i = 0; i < nvalues; i++)
914         values_total[m][i] += values_sum[m][i];
915       count_total[m] += count_sum[m];
916     }
917     normcount++;
918 
919   } else if (ave == WINDOW) {
920     for (m = 0; m < nchunk; m++) {
921       for (i = 0; i < nvalues; i++) {
922         values_total[m][i] += values_sum[m][i];
923         if (window_limit) values_total[m][i] -= values_list[iwindow][m][i];
924         values_list[iwindow][m][i] = values_sum[m][i];
925       }
926       count_total[m] += count_sum[m];
927       if (window_limit) count_total[m] -= count_list[iwindow][m];
928       count_list[iwindow][m] = count_sum[m];
929     }
930 
931     iwindow++;
932     if (iwindow == nwindow) {
933       iwindow = 0;
934       window_limit = 1;
935     }
936     if (window_limit) normcount = nwindow;
937     else normcount = iwindow;
938   }
939 
940   // output result to file
941 
942   if (fp && comm->me == 0) {
943     clearerr(fp);
944     if (overwrite) fseek(fp,filepos,SEEK_SET);
945     double count = 0.0;
946     for (m = 0; m < nchunk; m++) count += count_total[m];
947     fprintf(fp,BIGINT_FORMAT " %d %g\n",ntimestep,nchunk,count);
948 
949     int compress = cchunk->compress;
950     int *chunkID = cchunk->chunkID;
951     int ncoord = cchunk->ncoord;
952     double **coord = cchunk->coord;
953 
954     if (!compress) {
955       if (ncoord == 0) {
956         for (m = 0; m < nchunk; m++) {
957           fprintf(fp,"  %d %g",m+1,count_total[m]/normcount);
958           for (i = 0; i < nvalues; i++)
959             fprintf(fp,format,values_total[m][i]/normcount);
960           fprintf(fp,"\n");
961         }
962       } else if (ncoord == 1) {
963         for (m = 0; m < nchunk; m++) {
964           fprintf(fp,"  %d %g %g",m+1,coord[m][0],
965                   count_total[m]/normcount);
966           for (i = 0; i < nvalues; i++)
967             fprintf(fp,format,values_total[m][i]/normcount);
968           fprintf(fp,"\n");
969         }
970       } else if (ncoord == 2) {
971         for (m = 0; m < nchunk; m++) {
972           fprintf(fp,"  %d %g %g %g",m+1,coord[m][0],coord[m][1],
973                   count_total[m]/normcount);
974           for (i = 0; i < nvalues; i++)
975             fprintf(fp,format,values_total[m][i]/normcount);
976           fprintf(fp,"\n");
977         }
978       } else if (ncoord == 3) {
979         for (m = 0; m < nchunk; m++) {
980           fprintf(fp,"  %d %g %g %g %g",m+1,
981                   coord[m][0],coord[m][1],coord[m][2],count_total[m]/normcount);
982           for (i = 0; i < nvalues; i++)
983             fprintf(fp,format,values_total[m][i]/normcount);
984           fprintf(fp,"\n");
985         }
986       }
987     } else {
988       if (ncoord == 0) {
989         for (m = 0; m < nchunk; m++) {
990           fprintf(fp,"  %d %d %g",m+1,chunkID[m],count_total[m]/normcount);
991           for (i = 0; i < nvalues; i++)
992             fprintf(fp,format,values_total[m][i]/normcount);
993           fprintf(fp,"\n");
994         }
995       } else if (ncoord == 1) {
996         for (m = 0; m < nchunk; m++) {
997           j = chunkID[m];
998           fprintf(fp,"  %d %d %g %g",m+1,j,coord[j-1][0],
999                   count_total[m]/normcount);
1000           for (i = 0; i < nvalues; i++)
1001             fprintf(fp,format,values_total[m][i]/normcount);
1002           fprintf(fp,"\n");
1003         }
1004       } else if (ncoord == 2) {
1005         for (m = 0; m < nchunk; m++) {
1006           j = chunkID[m];
1007           fprintf(fp,"  %d %d %g %g %g",m+1,j,coord[j-1][0],coord[j-1][1],
1008                   count_total[m]/normcount);
1009           for (i = 0; i < nvalues; i++)
1010             fprintf(fp,format,values_total[m][i]/normcount);
1011           fprintf(fp,"\n");
1012         }
1013       } else if (ncoord == 3) {
1014         for (m = 0; m < nchunk; m++) {
1015           j = chunkID[m];
1016           fprintf(fp,"  %d %d %g %g %g %g",m+1,j,coord[j-1][0],
1017                   coord[j-1][1],coord[j-1][2],count_total[m]/normcount);
1018           for (i = 0; i < nvalues; i++)
1019             fprintf(fp,format,values_total[m][i]/normcount);
1020           fprintf(fp,"\n");
1021         }
1022       }
1023     }
1024     if (ferror(fp))
1025       error->one(FLERR,"Error writing averaged chunk data");
1026 
1027     fflush(fp);
1028 
1029     if (overwrite) {
1030       long fileend = ftell(fp);
1031       if ((fileend > 0) && (ftruncate(fileno(fp),fileend)))
1032         perror("Error while tuncating output");
1033     }
1034   }
1035 }
1036 
1037 /* ----------------------------------------------------------------------
1038    allocate all per-chunk vectors
1039 ------------------------------------------------------------------------- */
1040 
allocate()1041 void FixAveChunk::allocate()
1042 {
1043   size_array_rows = nchunk;
1044 
1045   // reallocate chunk arrays if needed
1046 
1047   if (nchunk > maxchunk) {
1048     maxchunk = nchunk;
1049     memory->grow(count_one,nchunk,"ave/chunk:count_one");
1050     memory->grow(count_many,nchunk,"ave/chunk:count_many");
1051     memory->grow(count_sum,nchunk,"ave/chunk:count_sum");
1052     memory->grow(count_total,nchunk,"ave/chunk:count_total");
1053 
1054     memory->grow(values_one,nchunk,nvalues,"ave/chunk:values_one");
1055     memory->grow(values_many,nchunk,nvalues,"ave/chunk:values_many");
1056     memory->grow(values_sum,nchunk,nvalues,"ave/chunk:values_sum");
1057     memory->grow(values_total,nchunk,nvalues,"ave/chunk:values_total");
1058 
1059     // only allocate count and values list for ave = WINDOW
1060 
1061     if (ave == WINDOW) {
1062       memory->create(count_list,nwindow,nchunk,"ave/chunk:count_list");
1063       memory->create(values_list,nwindow,nchunk,nvalues,
1064                      "ave/chunk:values_list");
1065     }
1066 
1067     // reinitialize regrown count/values total since they accumulate
1068 
1069     int i,m;
1070     for (m = 0; m < nchunk; m++) {
1071       for (i = 0; i < nvalues; i++) values_total[m][i] = 0.0;
1072       count_total[m] = 0.0;
1073     }
1074   }
1075 }
1076 
1077 /* ----------------------------------------------------------------------
1078    return I,J array value
1079    if I exceeds current nchunks, return 0.0 instead of generating an error
1080    columns 1 to colextra = chunkID + ncoord
1081    next column = count, remaining columns = Nvalues
1082 ------------------------------------------------------------------------- */
1083 
compute_array(int i,int j)1084 double FixAveChunk::compute_array(int i, int j)
1085 {
1086   if (values_total == nullptr) return 0.0;
1087   if (i >= nchunk) return 0.0;
1088   if (j < colextra) {
1089     if (cchunk->compress) {
1090       if (j == 0) return (double) cchunk->chunkID[i];
1091       return cchunk->coord[i][j-1];
1092     } else return cchunk->coord[i][j];
1093   }
1094   j -= colextra + 1;
1095   if (!normcount) return 0.0;
1096   if (j < 0) return count_total[i]/normcount;
1097   return values_total[i][j]/normcount;
1098 }
1099 
1100 /* ----------------------------------------------------------------------
1101    calculate nvalid = next step on which end_of_step does something
1102    can be this timestep if multiple of nfreq and nrepeat = 1
1103    else backup from next multiple of nfreq
1104 ------------------------------------------------------------------------- */
1105 
nextvalid()1106 bigint FixAveChunk::nextvalid()
1107 {
1108   bigint nvalid = (update->ntimestep/nfreq)*nfreq + nfreq;
1109   if (nvalid-nfreq == update->ntimestep && nrepeat == 1)
1110     nvalid = update->ntimestep;
1111   else
1112     nvalid -= ((bigint)nrepeat-1)*nevery;
1113   if (nvalid < update->ntimestep) nvalid += nfreq;
1114   return nvalid;
1115 }
1116 
1117 /* ----------------------------------------------------------------------
1118    memory usage of varatom and bins
1119 ------------------------------------------------------------------------- */
1120 
memory_usage()1121 double FixAveChunk::memory_usage()
1122 {
1123   double bytes = (double)maxvar * sizeof(double);         // varatom
1124   bytes += (double)4*maxchunk * sizeof(double);           // count one,many,sum,total
1125   bytes += (double)nvalues*maxchunk * sizeof(double);     // values one,many,sum,total
1126   bytes += (double)nwindow*maxchunk * sizeof(double);          // count_list
1127   bytes += (double)nwindow*maxchunk*nvalues * sizeof(double);  // values_list
1128   return bytes;
1129 }
1130