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