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 // NOTE: allow for bin center to be variables for sphere/cylinder
16 
17 #include "compute_chunk_atom.h"
18 
19 #include "arg_info.h"
20 #include "atom.h"
21 #include "comm.h"
22 #include "domain.h"
23 #include "error.h"
24 #include "fix.h"
25 #include "fix_store.h"
26 #include "group.h"
27 #include "input.h"
28 #include "lattice.h"
29 #include "math_const.h"
30 #include "memory.h"
31 #include "modify.h"
32 #include "region.h"
33 #include "update.h"
34 #include "variable.h"
35 
36 #include <cmath>
37 #include <cstring>
38 #include <map>
39 #include <utility>
40 
41 using namespace LAMMPS_NS;
42 using namespace MathConst;
43 
44 enum{LOWER,CENTER,UPPER,COORD};
45 enum{BOX,LATTICE,REDUCED};
46 enum{NODISCARD,MIXED,YESDISCARD};
47 enum{ONCE,NFREQ,EVERY};              // used in several files
48 enum{LIMITMAX,LIMITEXACT};
49 
50 #define IDMAX 1024*1024
51 
52 /* ---------------------------------------------------------------------- */
53 
ComputeChunkAtom(LAMMPS * lmp,int narg,char ** arg)54 ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
55   Compute(lmp, narg, arg),
56   chunk_volume_vec(nullptr), coord(nullptr), ichunk(nullptr), chunkID(nullptr),
57   cfvid(nullptr), idregion(nullptr), region(nullptr), cchunk(nullptr), fchunk(nullptr),
58   varatom(nullptr), id_fix(nullptr), fixstore(nullptr), lockfix(nullptr), chunk(nullptr),
59   exclude(nullptr), hash(nullptr)
60 {
61   if (narg < 4) error->all(FLERR,"Illegal compute chunk/atom command");
62 
63   peratom_flag = 1;
64   scalar_flag = 1;
65   extscalar = 0;
66   size_peratom_cols = 0;
67   create_attribute = 1;
68 
69   // chunk style and its args
70 
71   int iarg = 0;
72 
73   binflag = 0;
74   ncoord = 0;
75   cfvid = nullptr;
76 
77   if (strcmp(arg[3],"bin/1d") == 0) {
78     binflag = 1;
79     which = ArgInfo::BIN1D;
80     ncoord = 1;
81     iarg = 4;
82     readdim(narg,arg,iarg,0);
83     iarg += 3;
84   } else if (strcmp(arg[3],"bin/2d") == 0) {
85     binflag = 1;
86     which = ArgInfo::BIN2D;
87     ncoord = 2;
88     iarg = 4;
89     readdim(narg,arg,iarg,0);
90     readdim(narg,arg,iarg+3,1);
91     iarg += 6;
92   } else if (strcmp(arg[3],"bin/3d") == 0) {
93     binflag = 1;
94     which = ArgInfo::BIN3D;
95     ncoord = 3;
96     iarg = 4;
97     readdim(narg,arg,iarg,0);
98     readdim(narg,arg,iarg+3,1);
99     readdim(narg,arg,iarg+6,2);
100     iarg += 9;
101 
102   } else if (strcmp(arg[3],"bin/sphere") == 0) {
103     binflag = 1;
104     which = ArgInfo::BINSPHERE;
105     ncoord = 1;
106     iarg = 4;
107     if (iarg+6 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
108     sorigin_user[0] = utils::numeric(FLERR,arg[iarg],false,lmp);
109     sorigin_user[1] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
110     sorigin_user[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
111     sradmin_user = utils::numeric(FLERR,arg[iarg+3],false,lmp);
112     sradmax_user = utils::numeric(FLERR,arg[iarg+4],false,lmp);
113     nsbin = utils::inumeric(FLERR,arg[iarg+5],false,lmp);
114     iarg += 6;
115   } else if (strcmp(arg[3],"bin/cylinder") == 0) {
116     binflag = 1;
117     which = ArgInfo::BINCYLINDER;
118     ncoord = 2;
119     iarg = 4;
120     readdim(narg,arg,iarg,0);
121     iarg += 3;
122     if (dim[0] == 0) {
123       cdim1 = 1;
124       cdim2 = 2;
125     } else if (dim[0] == 1) {
126       cdim1 = 0;
127       cdim2 = 2;
128     } else {
129       cdim1 = 0;
130       cdim2 = 1;
131     }
132     if (iarg+5 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
133     corigin_user[dim[0]] = 0.0;
134     corigin_user[cdim1] = utils::numeric(FLERR,arg[iarg],false,lmp);
135     corigin_user[cdim2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
136     cradmin_user = utils::numeric(FLERR,arg[iarg+2],false,lmp);
137     cradmax_user = utils::numeric(FLERR,arg[iarg+3],false,lmp);
138 
139     ncbin = utils::inumeric(FLERR,arg[iarg+4],false,lmp);
140     iarg += 5;
141 
142   } else if (strcmp(arg[3],"type") == 0) {
143     which = ArgInfo::TYPE;
144     iarg = 4;
145   } else if (strcmp(arg[3],"molecule") == 0) {
146     which = ArgInfo::MOLECULE;
147     iarg = 4;
148 
149   } else {
150     ArgInfo argi(arg[3]);
151 
152     which = argi.get_type();
153     argindex = argi.get_index1();
154     cfvid = argi.copy_name();
155 
156     if ((which == ArgInfo::UNKNOWN) || (which == ArgInfo::NONE)
157         || (argi.get_dim() > 1))
158       error->all(FLERR,"Illegal compute chunk/atom command");
159     iarg = 4;
160   }
161 
162   // optional args
163 
164   regionflag = 0;
165   idregion = nullptr;
166   nchunksetflag = 0;
167   nchunkflag = EVERY;
168   limit = 0;
169   limitstyle = LIMITMAX;
170   limitfirst = 0;
171   idsflag = EVERY;
172   compress = 0;
173   int discardsetflag = 0;
174   discard = MIXED;
175   minflag[0] = LOWER;
176   minflag[1] = LOWER;
177   minflag[2] = LOWER;
178   maxflag[0] = UPPER;
179   maxflag[1] = UPPER;
180   maxflag[2] = UPPER;
181   scaleflag = LATTICE;
182   pbcflag = 0;
183 
184   while (iarg < narg) {
185     if (strcmp(arg[iarg],"region") == 0) {
186       if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
187       int iregion = domain->find_region(arg[iarg+1]);
188       if (iregion == -1)
189         error->all(FLERR,"Region ID for compute chunk/atom does not exist");
190       idregion = utils::strdup(arg[iarg+1]);
191       regionflag = 1;
192       iarg += 2;
193     } else if (strcmp(arg[iarg],"nchunk") == 0) {
194       if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
195       if (strcmp(arg[iarg+1],"once") == 0) nchunkflag = ONCE;
196       else if (strcmp(arg[iarg+1],"every") == 0) nchunkflag = EVERY;
197       else error->all(FLERR,"Illegal compute chunk/atom command");
198       nchunksetflag = 1;
199       iarg += 2;
200     } else if (strcmp(arg[iarg],"limit") == 0) {
201       if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
202       limit = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
203       if (limit < 0) error->all(FLERR,"Illegal compute chunk/atom command");
204       if (limit && !compress) limitfirst = 1;
205       iarg += 2;
206       if (limit) {
207         if (iarg > narg)
208           error->all(FLERR,"Illegal compute chunk/atom command");
209         if (strcmp(arg[iarg],"max") == 0) limitstyle = LIMITMAX;
210         else if (strcmp(arg[iarg],"exact") == 0) limitstyle = LIMITEXACT;
211         else error->all(FLERR,"Illegal compute chunk/atom command");
212         iarg++;
213       }
214     } else if (strcmp(arg[iarg],"ids") == 0) {
215       if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
216       if (strcmp(arg[iarg+1],"once") == 0) idsflag = ONCE;
217       else if (strcmp(arg[iarg+1],"nfreq") == 0) idsflag = NFREQ;
218       else if (strcmp(arg[iarg+1],"every") == 0) idsflag = EVERY;
219       else error->all(FLERR,"Illegal compute chunk/atom command");
220       iarg += 2;
221     } else if (strcmp(arg[iarg],"compress") == 0) {
222       if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
223       else if (strcmp(arg[iarg+1],"no") == 0) compress = 0;
224       else if (strcmp(arg[iarg+1],"yes") == 0) compress = 1;
225       else error->all(FLERR,"Illegal compute chunk/atom command");
226       iarg += 2;
227     } else if (strcmp(arg[iarg],"discard") == 0) {
228       if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
229       if (strcmp(arg[iarg+1],"mixed") == 0) discard = MIXED;
230       else if (strcmp(arg[iarg+1],"no") == 0) discard = NODISCARD;
231       else if (strcmp(arg[iarg+1],"yes") == 0) discard = YESDISCARD;
232       else error->all(FLERR,"Illegal compute chunk/atom command");
233       discardsetflag = 1;
234       iarg += 2;
235     } else if (strcmp(arg[iarg],"bound") == 0) {
236       if (iarg+4 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
237       int idim = 0;
238       if (strcmp(arg[iarg+1],"x") == 0) idim = 0;
239       else if (strcmp(arg[iarg+1],"y") == 0) idim = 1;
240       else if (strcmp(arg[iarg+1],"z") == 0) idim = 2;
241       else error->all(FLERR,"Illegal compute chunk/atom command");
242       minflag[idim] = COORD;
243       if (strcmp(arg[iarg+2],"lower") == 0) minflag[idim] = LOWER;
244       else minvalue[idim] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
245       maxflag[idim] = COORD;
246       if (strcmp(arg[iarg+3],"upper") == 0) maxflag[idim] = UPPER;
247       else maxvalue[idim] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
248       iarg += 4;
249     } else if (strcmp(arg[iarg],"units") == 0) {
250       if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
251       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = BOX;
252       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = LATTICE;
253       else if (strcmp(arg[iarg+1],"reduced") == 0) scaleflag = REDUCED;
254       else error->all(FLERR,"Illegal compute chunk/atom command");
255       iarg += 2;
256     } else if (strcmp(arg[iarg],"pbc") == 0) {
257       if (iarg+2 > narg) error->all(FLERR,"Illegal compute chunk/atom command");
258       if (strcmp(arg[iarg+1],"no") == 0) pbcflag = 0;
259       else if (strcmp(arg[iarg+1],"yes") == 0) pbcflag = 1;
260       else error->all(FLERR,"Illegal compute chunk/atom command");
261       iarg += 2;
262     } else error->all(FLERR,"Illegal compute chunk/atom command");
263   }
264 
265   // set nchunkflag and discard to default values if not explicitly set
266   // for binning style, also check in init() if simulation box is static,
267   //   which sets nchunkflag = ONCE
268 
269   if (!nchunksetflag) {
270     if (binflag) {
271       if (scaleflag == REDUCED) nchunkflag = ONCE;
272       else nchunkflag = EVERY;
273     }
274     if (which == ArgInfo::TYPE) nchunkflag = ONCE;
275     if (which == ArgInfo::MOLECULE) {
276       if (regionflag) nchunkflag = EVERY;
277       else nchunkflag = ONCE;
278     }
279     if (compress) nchunkflag = EVERY;
280   }
281 
282   if (!discardsetflag) {
283     if (binflag) discard = MIXED;
284     else discard = YESDISCARD;
285   }
286 
287   // error checks
288 
289   if (which == ArgInfo::MOLECULE && !atom->molecule_flag)
290     error->all(FLERR,"Compute chunk/atom molecule for non-molecular system");
291 
292   if (!binflag && discard == MIXED)
293     error->all(FLERR,"Compute chunk/atom without bins "
294                "cannot use discard mixed");
295   if (which == ArgInfo::BIN1D && delta[0] <= 0.0)
296     error->all(FLERR,"Illegal compute chunk/atom command");
297   if (which == ArgInfo::BIN2D && (delta[0] <= 0.0 || delta[1] <= 0.0))
298     error->all(FLERR,"Illegal compute chunk/atom command");
299   if (which == ArgInfo::BIN2D && (dim[0] == dim[1]))
300     error->all(FLERR,"Illegal compute chunk/atom command");
301   if (which == ArgInfo::BIN3D &&
302       (delta[0] <= 0.0 || delta[1] <= 0.0 || delta[2] <= 0.0))
303       error->all(FLERR,"Illegal compute chunk/atom command");
304   if (which == ArgInfo::BIN3D &&
305       (dim[0] == dim[1] || dim[1] == dim[2] || dim[0] == dim[2]))
306       error->all(FLERR,"Illegal compute chunk/atom command");
307   if (which == ArgInfo::BINSPHERE) {
308     if (domain->dimension == 2 && sorigin_user[2] != 0.0)
309       error->all(FLERR,"Compute chunk/atom sphere z origin must be 0.0 for 2d");
310     if (sradmin_user < 0.0 || sradmin_user >= sradmax_user || nsbin < 1)
311       error->all(FLERR,"Illegal compute chunk/atom command");
312   }
313   if (which == ArgInfo::BINCYLINDER) {
314     if (delta[0] <= 0.0)
315       error->all(FLERR,"Illegal compute chunk/atom command");
316     if (domain->dimension == 2 && dim[0] != 2)
317       error->all(FLERR,"Compute chunk/atom cylinder axis must be z for 2d");
318     if (cradmin_user < 0.0 || cradmin_user >= cradmax_user || ncbin < 1)
319       error->all(FLERR,"Illegal compute chunk/atom command");
320   }
321 
322   if (which == ArgInfo::COMPUTE) {
323     int icompute = modify->find_compute(cfvid);
324     if (icompute < 0)
325       error->all(FLERR,"Compute ID for compute chunk /atom does not exist");
326     if (modify->compute[icompute]->peratom_flag == 0)
327       error->all(FLERR,
328                  "Compute chunk/atom compute does not calculate "
329                  "per-atom values");
330     if (argindex == 0 &&
331         modify->compute[icompute]->size_peratom_cols != 0)
332       error->all(FLERR,"Compute chunk/atom compute does not "
333                  "calculate a per-atom vector");
334     if (argindex && modify->compute[icompute]->size_peratom_cols == 0)
335       error->all(FLERR,"Compute chunk/atom compute does not "
336                  "calculate a per-atom array");
337     if (argindex &&
338         argindex > modify->compute[icompute]->size_peratom_cols)
339       error->all(FLERR,"Compute chunk/atom compute array is "
340                  "accessed out-of-range");
341   }
342 
343   if (which == ArgInfo::FIX) {
344     int ifix = modify->find_fix(cfvid);
345     if (ifix < 0)
346       error->all(FLERR,"Fix ID for compute chunk/atom does not exist");
347     if (modify->fix[ifix]->peratom_flag == 0)
348       error->all(FLERR,"Compute chunk/atom fix does not calculate "
349                  "per-atom values");
350     if (argindex == 0 && modify->fix[ifix]->size_peratom_cols != 0)
351       error->all(FLERR,
352                  "Compute chunk/atom fix does not calculate a per-atom vector");
353     if (argindex && modify->fix[ifix]->size_peratom_cols == 0)
354       error->all(FLERR,
355                  "Compute chunk/atom fix does not calculate a per-atom array");
356     if (argindex && argindex > modify->fix[ifix]->size_peratom_cols)
357       error->all(FLERR,"Compute chunk/atom fix array is accessed out-of-range");
358   }
359 
360   if (which == ArgInfo::VARIABLE) {
361     int ivariable = input->variable->find(cfvid);
362     if (ivariable < 0)
363       error->all(FLERR,"Variable name for compute chunk/atom does not exist");
364     if (input->variable->atomstyle(ivariable) == 0)
365       error->all(FLERR,"Compute chunk/atom variable is not "
366                  "atom-style variable");
367   }
368 
369   // setup scaling
370 
371   if (binflag) {
372     if (domain->triclinic == 1 && scaleflag != REDUCED)
373       error->all(FLERR,"Compute chunk/atom for triclinic boxes "
374                  "requires units reduced");
375   }
376 
377   if (scaleflag == LATTICE) {
378     xscale = domain->lattice->xlattice;
379     yscale = domain->lattice->ylattice;
380     zscale = domain->lattice->zlattice;
381   } else xscale = yscale = zscale = 1.0;
382 
383   // apply scaling factors and cylinder dims orthogonal to axis
384 
385   if (binflag) {
386     double scale = 1.0;
387     if (which == ArgInfo::BIN1D || which == ArgInfo::BIN2D
388         || which == ArgInfo::BIN3D || which == ArgInfo::BINCYLINDER) {
389       if (which == ArgInfo::BIN1D || which == ArgInfo::BINCYLINDER) ndim = 1;
390       if (which == ArgInfo::BIN2D) ndim = 2;
391       if (which == ArgInfo::BIN3D) ndim = 3;
392       for (int idim = 0; idim < ndim; idim++) {
393         if (dim[idim] == 0) scale = xscale;
394         else if (dim[idim] == 1) scale = yscale;
395         else if (dim[idim] == 2) scale = zscale;
396         delta[idim] *= scale;
397         invdelta[idim] = 1.0/delta[idim];
398         if (originflag[idim] == COORD) origin[idim] *= scale;
399         if (minflag[idim] == COORD) minvalue[idim] *= scale;
400         if (maxflag[idim] == COORD) maxvalue[idim] *= scale;
401       }
402     } else if (which == ArgInfo::BINSPHERE) {
403       sorigin_user[0] *= xscale;
404       sorigin_user[1] *= yscale;
405       sorigin_user[2] *= zscale;
406       sradmin_user *= xscale;     // radii are scaled by xscale
407       sradmax_user *= xscale;
408     } else if (which == ArgInfo::BINCYLINDER) {
409       if (dim[0] == 0) {
410         corigin_user[cdim1] *= yscale;
411         corigin_user[cdim2] *= zscale;
412         cradmin_user *= yscale;     // radii are scaled by first non-axis dim
413         cradmax_user *= yscale;
414       } else if (dim[0] == 1) {
415         corigin_user[cdim1] *= xscale;
416         corigin_user[cdim2] *= zscale;
417         cradmin_user *= xscale;
418         cradmax_user *= xscale;
419       } else {
420         corigin_user[cdim1] *= xscale;
421         corigin_user[cdim2] *= yscale;
422         cradmin_user *= xscale;
423         cradmax_user *= xscale;
424       }
425     }
426   }
427 
428   // initialize chunk vector and per-chunk info
429 
430   nmax = 0;
431   chunk = nullptr;
432   nmaxint = -1;
433   ichunk = nullptr;
434   exclude = nullptr;
435 
436   nchunk = 0;
437   chunk_volume_scalar = 1.0;
438   chunk_volume_vec = nullptr;
439   coord = nullptr;
440   chunkID = nullptr;
441 
442   // computeflag = 1 if this compute might invoke another compute
443   // during assign_chunk_ids()
444 
445   if (which == ArgInfo::COMPUTE || which == ArgInfo::FIX || which == ArgInfo::VARIABLE) computeflag = 1;
446   else computeflag = 0;
447 
448   // other initializations
449 
450   invoked_setup = -1;
451   invoked_ichunk = -1;
452 
453   id_fix = nullptr;
454   fixstore = nullptr;
455 
456   if (compress) hash = new std::map<tagint,int>();
457   else hash = nullptr;
458 
459   maxvar = 0;
460   varatom = nullptr;
461 
462   lockcount = 0;
463   lockfix = nullptr;
464 
465   if (which == ArgInfo::MOLECULE) molcheck = 1;
466   else molcheck = 0;
467 }
468 
469 /* ---------------------------------------------------------------------- */
470 
~ComputeChunkAtom()471 ComputeChunkAtom::~ComputeChunkAtom()
472 {
473   // check nfix in case all fixes have already been deleted
474 
475   if (id_fix && modify->nfix) modify->delete_fix(id_fix);
476   delete [] id_fix;
477 
478   memory->destroy(chunk);
479   memory->destroy(ichunk);
480   memory->destroy(exclude);
481   memory->destroy(chunk_volume_vec);
482   memory->destroy(coord);
483   memory->destroy(chunkID);
484 
485   delete [] idregion;
486   delete [] cfvid;
487   delete hash;
488 
489   memory->destroy(varatom);
490 }
491 
492 /* ---------------------------------------------------------------------- */
493 
init()494 void ComputeChunkAtom::init()
495 {
496   // set and check validity of region
497 
498   if (regionflag) {
499     int iregion = domain->find_region(idregion);
500     if (iregion == -1)
501       error->all(FLERR,"Region ID for compute chunk/atom does not exist");
502     region = domain->regions[iregion];
503   }
504 
505   // set compute,fix,variable
506 
507   if (which == ArgInfo::COMPUTE) {
508     int icompute = modify->find_compute(cfvid);
509     if (icompute < 0)
510       error->all(FLERR,"Compute ID for compute chunk/atom does not exist");
511     cchunk = modify->compute[icompute];
512   } else if (which == ArgInfo::FIX) {
513     int ifix = modify->find_fix(cfvid);
514     if (ifix < 0)
515       error->all(FLERR,"Fix ID for compute chunk/atom does not exist");
516     fchunk = modify->fix[ifix];
517   } else if (which == ArgInfo::VARIABLE) {
518     int ivariable = input->variable->find(cfvid);
519     if (ivariable < 0)
520       error->all(FLERR,"Variable name for compute chunk/atom does not exist");
521     vchunk = ivariable;
522   }
523 
524   // for style MOLECULE, check that no mol IDs exceed MAXSMALLINT
525   // don't worry about group or optional region
526 
527   if (which == ArgInfo::MOLECULE) {
528     tagint *molecule = atom->molecule;
529     int nlocal = atom->nlocal;
530     tagint maxone = -1;
531     for (int i = 0; i < nlocal; i++)
532       if (molecule[i] > maxone) maxone = molecule[i];
533     tagint maxall;
534     MPI_Allreduce(&maxone,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
535     if (maxall > MAXSMALLINT)
536       error->all(FLERR,"Molecule IDs too large for compute chunk/atom");
537   }
538 
539   // for binning, if nchunkflag not already set, set it to ONCE or EVERY
540   // depends on whether simulation box size is static or dynamic
541   // reset invoked_setup if this is not first run and box just became static
542 
543   if (binflag && !nchunksetflag && !compress && scaleflag != REDUCED) {
544     if (domain->box_change_size == 0) {
545       if (nchunkflag == EVERY && invoked_setup >= 0) invoked_setup = -1;
546       nchunkflag = ONCE;
547     } else nchunkflag = EVERY;
548   }
549 
550   // require nchunkflag = ONCE if idsflag = ONCE
551   // b/c nchunk cannot change if chunk IDs are frozen
552   // can't check until now since nchunkflag may have been adjusted in init()
553 
554   if (idsflag == ONCE && nchunkflag != ONCE)
555     error->all(FLERR,"Compute chunk/atom ids once but nchunk is not once");
556 
557   // create/destroy fix STORE for persistent chunk IDs as needed
558   // need to do this if idsflag = ONCE or locks will be used by other commands
559   // need to wait until init() so that fix command(s) are in place
560   //   they increment lockcount if they lock this compute
561   // fixstore ID = compute-ID + COMPUTE_STORE, fix group = compute group
562   // fixstore initializes all values to 0.0
563 
564   if ((idsflag == ONCE || lockcount) && !fixstore) {
565     id_fix = utils::strdup(id + std::string("_COMPUTE_STORE"));
566     fixstore = (FixStore *) modify->add_fix(fmt::format("{} {} STORE peratom 1 1",
567                                                         id_fix, group->names[igroup]));
568   }
569 
570   if ((idsflag != ONCE && !lockcount) && fixstore) {
571     modify->delete_fix(id_fix);
572     fixstore = nullptr;
573   }
574 }
575 
576 /* ----------------------------------------------------------------------
577    invoke setup_chunks and/or compute_ichunk if only done ONCE
578    so that nchunks or chunk IDs are assigned when this compute was specified
579      as opposed to first times compute_peratom() or compute_ichunk() is called
580 ------------------------------------------------------------------------- */
581 
setup()582 void ComputeChunkAtom::setup()
583 {
584   if (nchunkflag == ONCE) setup_chunks();
585   if (idsflag == ONCE) compute_ichunk();
586   else invoked_ichunk = -1;
587 }
588 
589 /* ----------------------------------------------------------------------
590    only called by classes that use per-atom computes in standard way
591      dump, variable, thermo output, other computes, etc
592    not called by fix chunk or compute chunk commands
593      they invoke setup_chunks() and compute_ichunk() directly
594 ------------------------------------------------------------------------- */
595 
compute_peratom()596 void ComputeChunkAtom::compute_peratom()
597 {
598   invoked_peratom = update->ntimestep;
599 
600   // grow floating point chunk vector if necessary
601 
602   if (atom->nmax > nmax) {
603     memory->destroy(chunk);
604     nmax = atom->nmax;
605     memory->create(chunk,nmax,"chunk/atom:chunk");
606     vector_atom = chunk;
607   }
608 
609   setup_chunks();
610   compute_ichunk();
611 
612   // copy integer indices into floating-point chunk vector
613 
614   int nlocal = atom->nlocal;
615   for (int i = 0; i < nlocal; i++) chunk[i] = ichunk[i];
616 }
617 
618 
619 /* ----------------------------------------------------------------------
620    to return the number of chunks, we first need to make certain
621    that compute_peratom() has been called.
622 ------------------------------------------------------------------------- */
compute_scalar()623 double ComputeChunkAtom::compute_scalar()
624 {
625   if (invoked_peratom != update->ntimestep)
626     compute_peratom();
627   invoked_scalar = update->ntimestep;
628 
629   return (scalar = nchunk);
630 }
631 
632 /* ----------------------------------------------------------------------
633    set lock, so that nchunk will not change from startstep to stopstep
634    called by fix for duration of time it requires lock
635    OK if called by multiple fix commands
636      error if all callers do not have same duration
637      last caller holds the lock, so it can also unlock
638    stopstep can be positive for final step of finite-size time window
639    or can be -1 for infinite-size time window
640 ------------------------------------------------------------------------- */
641 
lock(Fix * fixptr,bigint startstep,bigint stopstep)642 void ComputeChunkAtom::lock(Fix *fixptr, bigint startstep, bigint stopstep)
643 {
644   if (lockfix == nullptr) {
645     lockfix = fixptr;
646     lockstart = startstep;
647     lockstop = stopstep;
648     return;
649   }
650 
651   if (startstep != lockstart || stopstep != lockstop)
652     error->all(FLERR,"Two fix commands using "
653                "same compute chunk/atom command in incompatible ways");
654 
655   // set lock to last calling Fix, since it will be last to unlock()
656 
657   lockfix = fixptr;
658 }
659 
660 /* ----------------------------------------------------------------------
661    unset lock
662    can only be done by fix command that holds the lock
663 ------------------------------------------------------------------------- */
664 
unlock(Fix * fixptr)665 void ComputeChunkAtom::unlock(Fix *fixptr)
666 {
667   if (fixptr != lockfix) return;
668   lockfix = nullptr;
669 }
670 
671 /* ----------------------------------------------------------------------
672    assign chunk IDs from 1 to Nchunk to every atom, or 0 if not in chunk
673 ------------------------------------------------------------------------- */
674 
compute_ichunk()675 void ComputeChunkAtom::compute_ichunk()
676 {
677   int i;
678 
679   // skip if already done on this step
680 
681   if (invoked_ichunk == update->ntimestep) return;
682 
683   // if old IDs persist via storage in fixstore, then just retrieve them
684   // yes if idsflag = ONCE, and already done once
685   //   or if idsflag = NFREQ and lock is in place and are on later timestep
686   // else proceed to recalculate per-atom chunk assignments
687 
688   const int nlocal = atom->nlocal;
689   int restore = 0;
690   if (idsflag == ONCE && invoked_ichunk >= 0) restore = 1;
691   if (idsflag == NFREQ && lockfix && update->ntimestep > lockstart) restore = 1;
692 
693   if (restore) {
694     invoked_ichunk = update->ntimestep;
695     double *vstore = fixstore->vstore;
696     for (i = 0; i < nlocal; i++) ichunk[i] = static_cast<int> (vstore[i]);
697     return;
698   }
699 
700   // assign chunk IDs to atoms
701   // will exclude atoms not in group or in optional region
702   // already invoked if this is same timestep as last setup_chunks()
703   // however, when between runs or using rerun, we need it again.
704 
705   if ((update->ntimestep > invoked_setup) || (invoked_ichunk < 0)) assign_chunk_ids();
706 
707   invoked_ichunk = update->ntimestep;
708 
709   // compress chunk IDs via hash of the original uncompressed IDs
710   // also apply discard rule except for binning styles which already did
711 
712   if (compress) {
713     if (binflag) {
714       for (i = 0; i < nlocal; i++) {
715         if (exclude[i]) continue;
716         if (hash->find(ichunk[i]) == hash->end()) exclude[i] = 1;
717         else ichunk[i] = hash->find(ichunk[i])->second;
718       }
719     } else if (discard == NODISCARD) {
720       for (i = 0; i < nlocal; i++) {
721         if (exclude[i]) continue;
722         if (hash->find(ichunk[i]) == hash->end()) ichunk[i] = nchunk;
723         else ichunk[i] = hash->find(ichunk[i])->second;
724       }
725     } else {
726       for (i = 0; i < nlocal; i++) {
727         if (exclude[i]) continue;
728         if (hash->find(ichunk[i]) == hash->end()) exclude[i] = 1;
729         else ichunk[i] = hash->find(ichunk[i])->second;
730       }
731     }
732 
733   // else if no compression apply discard rule by itself
734 
735   } else {
736     if (discard == NODISCARD) {
737       for (i = 0; i < nlocal; i++) {
738         if (exclude[i]) continue;
739         if (ichunk[i] < 1 || ichunk[i] > nchunk) ichunk[i] = nchunk;;
740       }
741     } else {
742       for (i = 0; i < nlocal; i++) {
743         if (exclude[i]) continue;
744         if (ichunk[i] < 1 || ichunk[i] > nchunk) exclude[i] = 1;
745       }
746     }
747   }
748 
749   // set ichunk = 0 for excluded atoms
750   // this should set any ichunk values which have not yet been set
751 
752   for (i = 0; i < nlocal; i++)
753     if (exclude[i]) ichunk[i] = 0;
754 
755   // if newly calculated IDs need to persist, store them in fixstore
756   // yes if idsflag = ONCE or idsflag = NFREQ and lock is in place
757 
758   if (idsflag == ONCE || (idsflag == NFREQ && lockfix)) {
759     double *vstore = fixstore->vstore;
760     for (i = 0; i < nlocal; i++) vstore[i] = ichunk[i];
761   }
762 
763   // one-time check if which = MOLECULE and
764   // any chunks do not contain all atoms in the molecule
765 
766   if (molcheck) {
767     check_molecules();
768     molcheck = 0;
769   }
770 }
771 
772 /* ----------------------------------------------------------------------
773    setup chunks
774    return nchunk = # of chunks
775      all atoms will be assigned a chunk ID from 1 to Nchunk, or 0
776    also setup any internal state needed to quickly assign atoms to chunks
777    called from compute_peratom() and also directly from
778      fix chunk and compute chunk commands
779 ------------------------------------------------------------------------- */
780 
setup_chunks()781 int ComputeChunkAtom::setup_chunks()
782 {
783   if (invoked_setup == update->ntimestep) return nchunk;
784 
785   // check if setup needs to be done
786   // no if lock is in place
787   // no if nchunkflag = ONCE, and already done once
788   // otherwise yes
789   // even if no, check if need to re-compute bin volumes
790   //   so that fix ave/chunk can do proper density normalization
791 
792   int flag = 0;
793   if (lockfix) flag = 1;
794   if (nchunkflag == ONCE && invoked_setup >= 0) flag = 1;
795 
796   if (flag) {
797     if (binflag && scaleflag == REDUCED && domain->box_change_size)
798       bin_volumes();
799     return nchunk;
800   }
801 
802   invoked_setup = update->ntimestep;
803 
804   // assign chunk IDs to atoms
805   // will exclude atoms not in group or in optional region
806   // for binning styles, need to setup bins and their volume first
807   //   else chunk_volume_scalar = entire box volume
808   // IDs are needed to scan for max ID and for compress()
809 
810   if (binflag) {
811     if (which == ArgInfo::BIN1D || which == ArgInfo::BIN2D
812         || which == ArgInfo::BIN3D)
813       nchunk = setup_xyz_bins();
814     else if (which == ArgInfo::BINSPHERE) nchunk = setup_sphere_bins();
815     else if (which == ArgInfo::BINCYLINDER) nchunk = setup_cylinder_bins();
816     bin_volumes();
817   } else {
818     chunk_volume_scalar = domain->xprd * domain->yprd;
819     if (domain->dimension == 3) chunk_volume_scalar *= domain->zprd;
820   }
821 
822   assign_chunk_ids();
823 
824   // set nchunk for chunk styles other than binning
825   // for styles other than TYPE, scan for max ID
826 
827   if (which == ArgInfo::TYPE) nchunk = atom->ntypes;
828   else if (!binflag) {
829 
830     int nlocal = atom->nlocal;
831     int hi = -1;
832     for (int i = 0; i < nlocal; i++) {
833       if (exclude[i]) continue;
834       if (ichunk[i] > hi) hi = ichunk[i];
835     }
836 
837     MPI_Allreduce(&hi,&nchunk,1,MPI_INT,MPI_MAX,world);
838     if (nchunk <= 0) nchunk = 1;
839   }
840 
841   // apply limit setting as well as compression of chunks with no atoms
842   // if limit is set, there are 3 cases:
843   //   no compression, limit specified before compression, or vice versa
844 
845   if (limit && !binflag) {
846     if (!compress) {
847       if (limitstyle == LIMITMAX) nchunk = MIN(nchunk,limit);
848       else if (limitstyle == LIMITEXACT) nchunk = limit;
849     } else if (limitfirst) {
850       nchunk = MIN(nchunk,limit);
851     }
852   }
853 
854   if (compress) compress_chunk_ids();
855 
856   if (limit && !binflag && compress) {
857     if (limitstyle == LIMITMAX) nchunk = MIN(nchunk,limit);
858     else if (limitstyle == LIMITEXACT) nchunk = limit;
859   }
860 
861   return nchunk;
862 }
863 
864 /* ----------------------------------------------------------------------
865    assign chunk IDs for all atoms, via ichunk vector
866    except excluded atoms, their chunk IDs are set to 0 later
867    also set exclude vector to 0/1 for all atoms
868      excluded atoms are those not in group or in optional region
869    called from compute_ichunk() and setup_chunks()
870 ------------------------------------------------------------------------- */
871 
assign_chunk_ids()872 void ComputeChunkAtom::assign_chunk_ids()
873 {
874   int i;
875 
876   // grow integer chunk index vector if necessary
877 
878   if (atom->nmax > nmaxint) {
879     memory->destroy(ichunk);
880     memory->destroy(exclude);
881     nmaxint = atom->nmax;
882     memory->create(ichunk,nmaxint,"chunk/atom:ichunk");
883     memory->create(exclude,nmaxint,"chunk/atom:exclude");
884   }
885 
886   // update region if necessary
887 
888   if (regionflag) region->prematch();
889 
890   // exclude = 1 if atom is not assigned to a chunk
891   // exclude atoms not in group or not in optional region
892 
893   double **x = atom->x;
894   int *mask = atom->mask;
895   const int nlocal = atom->nlocal;
896 
897   if (regionflag) {
898     for (i = 0; i < nlocal; i++) {
899       if (mask[i] & groupbit &&
900           region->match(x[i][0],x[i][1],x[i][2])) exclude[i] = 0;
901       else exclude[i] = 1;
902     }
903   } else {
904     for (i = 0; i < nlocal; i++) {
905       if (mask[i] & groupbit) exclude[i] = 0;
906       else exclude[i] = 1;
907     }
908   }
909 
910   // set ichunk to style value for included atoms
911   // binning styles apply discard rule, others do not yet
912 
913   if (binflag) {
914     if (which == ArgInfo::BIN1D) atom2bin1d();
915     else if (which == ArgInfo::BIN2D) atom2bin2d();
916     else if (which == ArgInfo::BIN3D) atom2bin3d();
917     else if (which == ArgInfo::BINSPHERE) atom2binsphere();
918     else if (which == ArgInfo::BINCYLINDER) atom2bincylinder();
919 
920   } else if (which == ArgInfo::TYPE) {
921     int *type = atom->type;
922     for (i = 0; i < nlocal; i++) {
923       if (exclude[i]) continue;
924       ichunk[i] = type[i];
925     }
926 
927   } else if (which == ArgInfo::MOLECULE) {
928     tagint *molecule = atom->molecule;
929     for (i = 0; i < nlocal; i++) {
930       if (exclude[i]) continue;
931       ichunk[i] = static_cast<int> (molecule[i]);
932     }
933 
934   } else if (which == ArgInfo::COMPUTE) {
935     if (!(cchunk->invoked_flag & Compute::INVOKED_PERATOM)) {
936       cchunk->compute_peratom();
937       cchunk->invoked_flag |= Compute::INVOKED_PERATOM;
938     }
939 
940     if (argindex == 0) {
941       double *vec = cchunk->vector_atom;
942       for (i = 0; i < nlocal; i++) {
943         if (exclude[i]) continue;
944         ichunk[i] = static_cast<int> (vec[i]);
945       }
946     } else {
947       double **array = cchunk->array_atom;
948       int argm1 = argindex - 1;
949       for (i = 0; i < nlocal; i++) {
950         if (exclude[i]) continue;
951         ichunk[i] = static_cast<int> (array[i][argm1]);
952       }
953     }
954 
955   } else if (which == ArgInfo::FIX) {
956     if (update->ntimestep % fchunk->peratom_freq)
957       error->all(FLERR,"Fix used in compute chunk/atom not "
958                  "computed at compatible time");
959 
960     if (argindex == 0) {
961       double *vec = fchunk->vector_atom;
962       for (i = 0; i < nlocal; i++) {
963         if (exclude[i]) continue;
964         ichunk[i] = static_cast<int> (vec[i]);
965       }
966     } else {
967       double **array = fchunk->array_atom;
968       int argm1 = argindex - 1;
969       for (i = 0; i < nlocal; i++) {
970         if (exclude[i]) continue;
971         ichunk[i] = static_cast<int> (array[i][argm1]);
972       }
973     }
974 
975   } else if (which == ArgInfo::VARIABLE) {
976     if (atom->nmax > maxvar) {
977       maxvar = atom->nmax;
978       memory->destroy(varatom);
979       memory->create(varatom,maxvar,"chunk/atom:varatom");
980     }
981 
982     input->variable->compute_atom(vchunk,igroup,varatom,1,0);
983     for (i = 0; i < nlocal; i++) {
984       if (exclude[i]) continue;
985       ichunk[i] = static_cast<int> (varatom[i]);
986     }
987   }
988 }
989 
990 /* ----------------------------------------------------------------------
991    compress chunk IDs currently assigned to atoms across all processors
992      by removing those with no atoms assigned
993    current assignment excludes atoms not in group or in optional region
994    current Nchunk = max ID
995    operation:
996      use hash to store list of populated IDs that I own
997      add new IDs to populated lists communicated from all other procs
998      final hash has global list of populated ideas
999    reset Nchunk = length of global list
1000    called by setup_chunks() when setting Nchunk
1001    remapping of chunk IDs to smaller Nchunk occurs later in compute_ichunk()
1002 ------------------------------------------------------------------------- */
1003 
compress_chunk_ids()1004 void ComputeChunkAtom::compress_chunk_ids()
1005 {
1006   hash->clear();
1007 
1008   // put my IDs into hash
1009 
1010   int nlocal = atom->nlocal;
1011   for (int i = 0; i < nlocal; i++) {
1012     if (exclude[i]) continue;
1013     if (hash->find(ichunk[i]) == hash->end()) (*hash)[ichunk[i]] = 0;
1014   }
1015 
1016   // n = # of my populated IDs
1017   // nall = n summed across all procs
1018 
1019   int n = hash->size();
1020   bigint nbone = n;
1021   bigint nball;
1022   MPI_Allreduce(&nbone,&nball,1,MPI_LMP_BIGINT,MPI_SUM,world);
1023 
1024   // create my list of populated IDs
1025 
1026   int *list = nullptr;
1027   memory->create(list,n,"chunk/atom:list");
1028 
1029   n = 0;
1030   std::map<tagint,int>::iterator pos;
1031   for (pos = hash->begin(); pos != hash->end(); ++pos)
1032     list[n++] = pos->first;
1033 
1034   // if nall < 1M, just allgather all ID lists on every proc
1035   // else perform ring comm
1036   // add IDs from all procs to my hash
1037 
1038   if (nball <= IDMAX) {
1039 
1040     // setup for allgatherv
1041 
1042     int nprocs = comm->nprocs;
1043     int nall = nball;
1044     int *recvcounts,*displs,*listall;
1045     memory->create(recvcounts,nprocs,"chunk/atom:recvcounts");
1046     memory->create(displs,nprocs,"chunk/atom:displs");
1047     memory->create(listall,nall,"chunk/atom:listall");
1048 
1049     MPI_Allgather(&n,1,MPI_INT,recvcounts,1,MPI_INT,world);
1050 
1051     displs[0] = 0;
1052     for (int iproc = 1; iproc < nprocs; iproc++)
1053       displs[iproc] = displs[iproc-1] + recvcounts[iproc-1];
1054 
1055     // allgatherv acquires list of populated IDs from all procs
1056 
1057     MPI_Allgatherv(list,n,MPI_INT,listall,recvcounts,displs,MPI_INT,world);
1058 
1059     // add all unique IDs in listall to my hash
1060 
1061     for (int i = 0; i < nall; i++)
1062       if (hash->find(listall[i]) == hash->end()) (*hash)[listall[i]] = 0;
1063 
1064     // clean up
1065 
1066     memory->destroy(recvcounts);
1067     memory->destroy(displs);
1068     memory->destroy(listall);
1069 
1070   } else {
1071     comm->ring(n,sizeof(int),list,1,idring,nullptr,(void *)this,0);
1072   }
1073 
1074   memory->destroy(list);
1075 
1076   // nchunk = length of hash containing populated IDs from all procs
1077 
1078   nchunk = hash->size();
1079 
1080   // reset hash value of each original chunk ID to ordered index
1081   //   ordered index = new compressed chunk ID (1 to Nchunk)
1082   //   leverages fact that map stores keys in ascending order
1083   // also allocate and set chunkID = list of original chunk IDs
1084   //   used by fix ave/chunk and compute property/chunk
1085 
1086   memory->destroy(chunkID);
1087   memory->create(chunkID,nchunk,"chunk/atom:chunkID");
1088 
1089   n = 0;
1090   for (pos = hash->begin(); pos != hash->end(); ++pos) {
1091     chunkID[n] = pos->first;
1092     (*hash)[pos->first] = ++n;
1093   }
1094 }
1095 
1096 /* ----------------------------------------------------------------------
1097    callback from comm->ring()
1098    cbuf = list of N chunk IDs from another proc
1099    loop over the list, add each to my hash
1100    hash ends up storing all unique IDs across all procs
1101 ------------------------------------------------------------------------- */
1102 
idring(int n,char * cbuf,void * ptr)1103 void ComputeChunkAtom::idring(int n, char *cbuf, void *ptr)
1104 {
1105   ComputeChunkAtom *cptr = (ComputeChunkAtom *)ptr;
1106   tagint *list = (tagint *) cbuf;
1107   std::map<tagint,int> *hash = cptr->hash;
1108   for (int i = 0; i < n; i++) (*hash)[list[i]] = 0;
1109 }
1110 
1111 /* ----------------------------------------------------------------------
1112    one-time check for which = MOLECULE to check
1113      if each chunk contains all atoms in the molecule
1114    issue warning if not
1115    note that this check is without regard to discard rule
1116    if discard == NODISCARD, there is no easy way to check that all
1117      atoms in an out-of-bounds molecule were added to a chunk,
1118      some could have been excluded by group or region, others not
1119 ------------------------------------------------------------------------- */
1120 
check_molecules()1121 void ComputeChunkAtom::check_molecules()
1122 {
1123   tagint *molecule = atom->molecule;
1124   int nlocal = atom->nlocal;
1125 
1126   int flag = 0;
1127 
1128   if (!compress) {
1129     for (int i = 0; i < nlocal; i++) {
1130       if (molecule[i] > 0 && molecule[i] <= nchunk &&
1131           ichunk[i] == 0) flag = 1;
1132     }
1133   } else {
1134     int molid;
1135     for (int i = 0; i < nlocal; i++) {
1136       molid = static_cast<int> (molecule[i]);
1137       if (hash->find(molid) != hash->end() && ichunk[i] == 0) flag = 1;
1138     }
1139   }
1140 
1141   int flagall;
1142   MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
1143   if (flagall && comm->me == 0)
1144     error->warning(FLERR,
1145                    "One or more chunks do not contain all atoms in molecule");
1146 }
1147 
1148 /* ----------------------------------------------------------------------
1149    setup xyz spatial bins and their extent and coordinates
1150    return nbins = # of bins, will become # of chunks
1151    called from setup_chunks()
1152 ------------------------------------------------------------------------- */
1153 
setup_xyz_bins()1154 int ComputeChunkAtom::setup_xyz_bins()
1155 {
1156   int i,j,k,m,n,idim;
1157   double lo,hi,coord1,coord2;
1158 
1159   // lo = bin boundary immediately below boxlo or minvalue
1160   // hi = bin boundary immediately above boxhi or maxvalue
1161   // allocate and initialize arrays based on new bin count
1162 
1163   double binlo[3],binhi[3];
1164   if (scaleflag == REDUCED) {
1165     binlo[0] = domain->boxlo_lamda[0];
1166     binlo[1] = domain->boxlo_lamda[1];
1167     binlo[2] = domain->boxlo_lamda[2];
1168     binhi[0] = domain->boxhi_lamda[0];
1169     binhi[1] = domain->boxhi_lamda[1];
1170     binhi[2] = domain->boxhi_lamda[2];
1171   } else {
1172     binlo[0] = domain->boxlo[0];
1173     binlo[1] = domain->boxlo[1];
1174     binlo[2] = domain->boxlo[2];
1175     binhi[0] = domain->boxhi[0];
1176     binhi[1] = domain->boxhi[1];
1177     binhi[2] = domain->boxhi[2];
1178   }
1179 
1180   if (minflag[0] == COORD) binlo[0] = minvalue[0];
1181   if (minflag[1] == COORD) binlo[1] = minvalue[1];
1182   if (minflag[2] == COORD) binlo[2] = minvalue[2];
1183   if (maxflag[0] == COORD) binhi[0] = maxvalue[0];
1184   if (maxflag[1] == COORD) binhi[1] = maxvalue[1];
1185   if (maxflag[2] == COORD) binhi[2] = maxvalue[2];
1186 
1187   int nbins = 1;
1188 
1189   for (m = 0; m < ndim; m++) {
1190     idim = dim[m];
1191     if (originflag[m] == LOWER) origin[m] = binlo[idim];
1192     else if (originflag[m] == UPPER) origin[m] = binhi[idim];
1193     else if (originflag[m] == CENTER)
1194       origin[m] = 0.5 * (binlo[idim] + binhi[idim]);
1195 
1196     if (origin[m] < binlo[idim]) {
1197       n = static_cast<int> ((binlo[idim] - origin[m]) * invdelta[m]);
1198       lo = origin[m] + n*delta[m];
1199     } else {
1200       n = static_cast<int> ((origin[m] - binlo[idim]) * invdelta[m]);
1201       lo = origin[m] - n*delta[m];
1202       if (lo > binlo[idim]) lo -= delta[m];
1203     }
1204     if (origin[m] < binhi[idim]) {
1205       n = static_cast<int> ((binhi[idim] - origin[m]) * invdelta[m]);
1206       hi = origin[m] + n*delta[m];
1207       if (hi < binhi[idim]) hi += delta[m];
1208     } else {
1209       n = static_cast<int> ((origin[m] - binhi[idim]) * invdelta[m]);
1210       hi = origin[m] - n*delta[m];
1211     }
1212 
1213     if (lo > hi) error->all(FLERR,"Invalid bin bounds in compute chunk/atom");
1214 
1215     offset[m] = lo;
1216     nlayers[m] = static_cast<int> ((hi-lo) * invdelta[m] + 0.5);
1217     nbins *= nlayers[m];
1218   }
1219 
1220   // allocate and set bin coordinates
1221 
1222   memory->destroy(coord);
1223   memory->create(coord,nbins,ndim,"chunk/atom:coord");
1224 
1225   if (ndim == 1) {
1226     for (i = 0; i < nlayers[0]; i++)
1227       coord[i][0] = offset[0] + (i+0.5)*delta[0];
1228   } else if (ndim == 2) {
1229     m = 0;
1230     for (i = 0; i < nlayers[0]; i++) {
1231       coord1 = offset[0] + (i+0.5)*delta[0];
1232       for (j = 0; j < nlayers[1]; j++) {
1233         coord[m][0] = coord1;
1234         coord[m][1] = offset[1] + (j+0.5)*delta[1];
1235         m++;
1236       }
1237     }
1238   } else if (ndim == 3) {
1239     m = 0;
1240     for (i = 0; i < nlayers[0]; i++) {
1241       coord1 = offset[0] + (i+0.5)*delta[0];
1242       for (j = 0; j < nlayers[1]; j++) {
1243         coord2 = offset[1] + (j+0.5)*delta[1];
1244         for (k = 0; k < nlayers[2]; k++) {
1245           coord[m][0] = coord1;
1246           coord[m][1] = coord2;
1247           coord[m][2] = offset[2] + (k+0.5)*delta[2];
1248           m++;
1249         }
1250       }
1251     }
1252   }
1253 
1254   return nbins;
1255 }
1256 
1257 /* ----------------------------------------------------------------------
1258    setup spherical spatial bins and their single coordinate
1259    return nsphere = # of bins, will become # of chunks
1260    called from setup_chunks()
1261 ------------------------------------------------------------------------- */
1262 
setup_sphere_bins()1263 int ComputeChunkAtom::setup_sphere_bins()
1264 {
1265   // convert sorigin_user to sorigin
1266   // sorigin,srad are always in box units, for orthogonal or triclinic domains
1267   // lamda2x works for either orthogonal or triclinic
1268 
1269   if (scaleflag == REDUCED) {
1270     domain->lamda2x(sorigin_user,sorigin);
1271     sradmin = sradmin_user * (domain->boxhi[0]-domain->boxlo[0]);
1272     sradmax = sradmax_user * (domain->boxhi[0]-domain->boxlo[0]);
1273   } else {
1274     sorigin[0] = sorigin_user[0];
1275     sorigin[1] = sorigin_user[1];
1276     sorigin[2] = sorigin_user[2];
1277     sradmin = sradmin_user;
1278     sradmax = sradmax_user;
1279   }
1280 
1281   // if pbcflag set, sradmax must be < 1/2 box in any periodic dim
1282   // treat orthogonal and triclinic the same
1283   // check every time bins are created
1284 
1285   if (pbcflag) {
1286     double *prd_half = domain->prd_half;
1287     int *periodicity = domain->periodicity;
1288     int flag = 0;
1289     if (periodicity[0] && sradmax > prd_half[0]) flag = 1;
1290     if (periodicity[1] && sradmax > prd_half[1]) flag = 1;
1291     if (domain->dimension == 3 &&
1292         periodicity[2] && sradmax > prd_half[2]) flag = 1;
1293     if (flag)
1294       error->all(FLERR,"Compute chunk/atom bin/sphere radius "
1295                  "is too large for periodic box");
1296   }
1297 
1298   sinvrad = nsbin / (sradmax-sradmin);
1299 
1300   // allocate and set bin coordinates
1301   // coord = midpt of radii for a spherical shell
1302 
1303   memory->destroy(coord);
1304   memory->create(coord,nsbin,1,"chunk/atom:coord");
1305 
1306   double rlo,rhi;
1307 
1308   for (int i = 0; i < nsbin; i++) {
1309     rlo = sradmin + i * (sradmax-sradmin) / nsbin;
1310     rhi = sradmin + (i+1) * (sradmax-sradmin) / nsbin;
1311     if (i == nsbin-1) rhi = sradmax;
1312     coord[i][0] = 0.5 * (rlo+rhi);
1313   }
1314 
1315   return nsbin;
1316 }
1317 
1318 /* ----------------------------------------------------------------------
1319    setup cylindrical spatial bins and their single coordinate
1320    return nsphere = # of bins, will become # of chunks
1321    called from setup_chunks()
1322 ------------------------------------------------------------------------- */
1323 
setup_cylinder_bins()1324 int ComputeChunkAtom::setup_cylinder_bins()
1325 {
1326   // setup bins along cylinder axis
1327   // ncplane = # of axis bins
1328 
1329   ncplane = setup_xyz_bins();
1330 
1331   // convert corigin_user to corigin
1332   // corigin is always in box units, for orthogonal or triclinic domains
1333   // lamda2x works for either orthogonal or triclinic
1334 
1335   if (scaleflag == REDUCED) {
1336     domain->lamda2x(corigin_user,corigin);
1337     cradmin = cradmin_user * (domain->boxhi[cdim1]-domain->boxlo[cdim1]);
1338     cradmax = cradmax_user * (domain->boxhi[cdim1]-domain->boxlo[cdim1]);
1339   } else {
1340     corigin[cdim1] = corigin_user[cdim1];
1341     corigin[cdim2] = corigin_user[cdim2];
1342     cradmin = cradmin_user;
1343     cradmax = cradmax_user;
1344   }
1345 
1346   // if pbcflag set, sradmax must be < 1/2 box in any periodic non-axis dim
1347   // treat orthogonal and triclinic the same
1348   // check every time bins are created
1349 
1350   if (pbcflag) {
1351     double *prd_half = domain->prd_half;
1352     int *periodicity = domain->periodicity;
1353     int flag = 0;
1354     if (periodicity[cdim1] && sradmax > prd_half[cdim1]) flag = 1;
1355     if (periodicity[cdim2] && sradmax > prd_half[cdim2]) flag = 1;
1356     if (flag)
1357       error->all(FLERR,"Compute chunk/atom bin/cylinder radius "
1358                  "is too large for periodic box");
1359   }
1360 
1361   cinvrad = ncbin / (cradmax-cradmin);
1362 
1363   // allocate and set radial bin coordinates
1364   // radial coord = midpt of radii for a cylindrical shell
1365   // axiscoord = saved bin coords along cylndrical axis
1366   // radcoord = saved bin coords in radial direction
1367 
1368   double **axiscoord = coord;
1369   memory->create(coord,ncbin,1,"chunk/atom:coord");
1370   double **radcoord = coord;
1371 
1372   double rlo,rhi;
1373 
1374   for (int i = 0; i < ncbin; i++) {
1375     rlo = cradmin + i * (cradmax-cradmin) / ncbin;
1376     rhi = cradmin + (i+1) * (cradmax-cradmin) / ncbin;
1377     if (i == ncbin-1) rhi = cradmax;
1378     coord[i][0] = 0.5 * (rlo+rhi);
1379   }
1380 
1381   // create array of combined coords for all bins
1382 
1383   memory->create(coord,ncbin*ncplane,2,"chunk/atom:coord");
1384   int m = 0;
1385   for (int i = 0; i < ncbin; i++)
1386     for (int j = 0; j < ncplane; j++) {
1387       coord[m][0] = radcoord[i][0];
1388       coord[m][1] = axiscoord[j][0];
1389       m++;
1390     }
1391   memory->destroy(axiscoord);
1392   memory->destroy(radcoord);
1393 
1394   return ncbin*ncplane;
1395 }
1396 
1397 /* ----------------------------------------------------------------------
1398    calculate chunk volumes = bin volumes
1399    scalar if all bins have same volume
1400    vector if per-bin volumes are different
1401 ------------------------------------------------------------------------- */
1402 
bin_volumes()1403 void ComputeChunkAtom::bin_volumes()
1404 {
1405   if (which == ArgInfo::BIN1D || which == ArgInfo::BIN2D
1406       || which == ArgInfo::BIN3D) {
1407     if (domain->dimension == 3)
1408       chunk_volume_scalar = domain->xprd * domain->yprd * domain->zprd;
1409     else chunk_volume_scalar = domain->xprd * domain->yprd;
1410     double *prd;
1411     if (scaleflag == REDUCED) prd = domain->prd_lamda;
1412     else prd = domain->prd;
1413     for (int m = 0; m < ndim; m++)
1414       chunk_volume_scalar *= delta[m]/prd[dim[m]];
1415 
1416   } else if (which == ArgInfo::BINSPHERE) {
1417     memory->destroy(chunk_volume_vec);
1418     memory->create(chunk_volume_vec,nchunk,"chunk/atom:chunk_volume_vec");
1419     double rlo,rhi,vollo,volhi;
1420     for (int i = 0; i < nchunk; i++) {
1421       rlo = sradmin + i * (sradmax-sradmin) / nsbin;
1422       rhi = sradmin + (i+1) * (sradmax-sradmin) / nsbin;
1423       if (i == nchunk-1) rhi = sradmax;
1424       vollo = 4.0/3.0 * MY_PI * rlo*rlo*rlo;
1425       volhi = 4.0/3.0 * MY_PI * rhi*rhi*rhi;
1426       chunk_volume_vec[i] = volhi - vollo;
1427     }
1428 
1429   } else if (which == ArgInfo::BINCYLINDER) {
1430     memory->destroy(chunk_volume_vec);
1431     memory->create(chunk_volume_vec,nchunk,"chunk/atom:chunk_volume_vec");
1432 
1433     // slabthick = delta of bins along cylinder axis
1434 
1435     double *prd;
1436     if (scaleflag == REDUCED) prd = domain->prd_lamda;
1437     else prd = domain->prd;
1438     double slabthick = domain->prd[dim[0]] * delta[0]/prd[dim[0]];
1439 
1440     // area lo/hi of concentric circles in radial direction
1441 
1442     int iradbin;
1443     double rlo,rhi,arealo,areahi;
1444     for (int i = 0; i < nchunk; i++) {
1445       iradbin = i / ncplane;
1446       rlo = cradmin + iradbin * (cradmax-cradmin) / ncbin;
1447       rhi = cradmin + (iradbin+1) * (cradmax-cradmin) / ncbin;
1448       if (iradbin == ncbin-1) rhi = cradmax;
1449       arealo = MY_PI * rlo*rlo;
1450       areahi = MY_PI * rhi*rhi;
1451       chunk_volume_vec[i] = (areahi-arealo) * slabthick;
1452     }
1453   }
1454 }
1455 
1456 /* ----------------------------------------------------------------------
1457    assign each atom to a 1d spatial bin (layer)
1458 ------------------------------------------------------------------------- */
1459 
atom2bin1d()1460 void ComputeChunkAtom::atom2bin1d()
1461 {
1462   int i,ibin;
1463   double *boxlo,*boxhi,*prd;
1464   double xremap;
1465 
1466   double **x = atom->x;
1467   int nlocal = atom->nlocal;
1468 
1469   int idim = dim[0];
1470   int nlayer1m1 = nlayers[0] - 1;
1471   int periodicity = domain->periodicity[idim];
1472 
1473   if (periodicity) {
1474     if (scaleflag == REDUCED) {
1475       boxlo = domain->boxlo_lamda;
1476       boxhi = domain->boxhi_lamda;
1477       prd = domain->prd_lamda;
1478     } else {
1479       boxlo = domain->boxlo;
1480       boxhi = domain->boxhi;
1481       prd = domain->prd;
1482     }
1483   }
1484 
1485   // remap each atom's relevant coord back into box via PBC if necessary
1486   // if scaleflag = REDUCED, box coords -> lamda coords
1487   // apply discard rule
1488 
1489   if (scaleflag == REDUCED) domain->x2lamda(nlocal);
1490 
1491   for (i = 0; i < nlocal; i++) {
1492     if (exclude[i]) continue;
1493 
1494     xremap = x[i][idim];
1495     if (periodicity) {
1496       if (xremap < boxlo[idim]) xremap += prd[idim];
1497       if (xremap >= boxhi[idim]) xremap -= prd[idim];
1498     }
1499 
1500     ibin = static_cast<int> ((xremap - offset[0]) * invdelta[0]);
1501     if (xremap < offset[0]) ibin--;
1502 
1503     if (discard == MIXED) {
1504       if (!minflag[idim]) ibin = MAX(ibin,0);
1505       else if (ibin < 0) {
1506         exclude[i] = 1;
1507         continue;
1508       }
1509       if (!maxflag[idim]) ibin = MIN(ibin,nlayer1m1);
1510       else if (ibin > nlayer1m1) {
1511         exclude[i] = 1;
1512         continue;
1513       }
1514     } else if (discard == NODISCARD) {
1515       ibin = MAX(ibin,0);
1516       ibin = MIN(ibin,nlayer1m1);
1517     } else if (ibin < 0 || ibin > nlayer1m1) {
1518       exclude[i] = 1;
1519       continue;
1520     }
1521 
1522     ichunk[i] = ibin+1;
1523   }
1524 
1525   if (scaleflag == REDUCED) domain->lamda2x(nlocal);
1526 }
1527 
1528 /* ----------------------------------------------------------------------
1529    assign each atom to a 2d spatial bin (pencil)
1530 ------------------------------------------------------------------------- */
1531 
atom2bin2d()1532 void ComputeChunkAtom::atom2bin2d()
1533 {
1534   int i,ibin,i1bin,i2bin;
1535   double *boxlo,*boxhi,*prd;
1536   double xremap,yremap;
1537 
1538   double **x = atom->x;
1539   int nlocal = atom->nlocal;
1540 
1541   int idim = dim[0];
1542   int jdim = dim[1];
1543   int nlayer1m1 = nlayers[0] - 1;
1544   int nlayer2m1 = nlayers[1] - 1;
1545   int *periodicity = domain->periodicity;
1546 
1547   if (periodicity[idim] || periodicity[jdim]) {
1548     if (scaleflag == REDUCED) {
1549       boxlo = domain->boxlo_lamda;
1550       boxhi = domain->boxhi_lamda;
1551       prd = domain->prd_lamda;
1552     } else {
1553       boxlo = domain->boxlo;
1554       boxhi = domain->boxhi;
1555       prd = domain->prd;
1556     }
1557   }
1558 
1559   // remap each atom's relevant coord back into box via PBC if necessary
1560   // if scaleflag = REDUCED, box coords -> lamda coords
1561   // apply discard rule
1562 
1563   if (scaleflag == REDUCED) domain->x2lamda(nlocal);
1564 
1565   for (i = 0; i < nlocal; i++) {
1566     if (exclude[i]) continue;
1567 
1568     xremap = x[i][idim];
1569     if (periodicity[idim]) {
1570       if (xremap < boxlo[idim]) xremap += prd[idim];
1571       if (xremap >= boxhi[idim]) xremap -= prd[idim];
1572     }
1573 
1574     i1bin = static_cast<int> ((xremap - offset[0]) * invdelta[0]);
1575     if (xremap < offset[0]) i1bin--;
1576 
1577     if (discard == MIXED) {
1578       if (!minflag[idim]) i1bin = MAX(i1bin,0);
1579       else if (i1bin < 0) {
1580         exclude[i] = 1;
1581         continue;
1582       }
1583       if (!maxflag[idim]) i1bin = MIN(i1bin,nlayer1m1);
1584       else if (i1bin > nlayer1m1) {
1585         exclude[i] = 1;
1586         continue;
1587       }
1588     } else if (discard == NODISCARD) {
1589       i1bin = MAX(i1bin,0);
1590       i1bin = MIN(i1bin,nlayer1m1);
1591     } else if (i1bin < 0 || i1bin > nlayer1m1) {
1592       exclude[i] = 1;
1593       continue;
1594     }
1595 
1596     yremap = x[i][jdim];
1597     if (periodicity[jdim]) {
1598       if (yremap < boxlo[jdim]) yremap += prd[jdim];
1599       if (yremap >= boxhi[jdim]) yremap -= prd[jdim];
1600     }
1601 
1602     i2bin = static_cast<int> ((yremap - offset[1]) * invdelta[1]);
1603     if (yremap < offset[1]) i2bin--;
1604 
1605     if (discard == MIXED) {
1606       if (!minflag[jdim]) i2bin = MAX(i2bin,0);
1607       else if (i2bin < 0) {
1608         exclude[i] = 1;
1609         continue;
1610       }
1611       if (!maxflag[jdim]) i2bin = MIN(i2bin,nlayer2m1);
1612       else if (i2bin > nlayer2m1) {
1613         exclude[i] = 1;
1614         continue;
1615       }
1616     } else if (discard == NODISCARD) {
1617       i2bin = MAX(i2bin,0);
1618       i2bin = MIN(i2bin,nlayer2m1);
1619     } else if (i2bin < 0 || i2bin > nlayer2m1) {
1620       exclude[i] = 1;
1621       continue;
1622     }
1623 
1624     ibin = i1bin*nlayers[1] + i2bin;
1625     ichunk[i] = ibin+1;
1626   }
1627 
1628   if (scaleflag == REDUCED) domain->lamda2x(nlocal);
1629 }
1630 
1631 /* ----------------------------------------------------------------------
1632    assign each atom to a 3d spatial bin (brick)
1633 ------------------------------------------------------------------------- */
1634 
atom2bin3d()1635 void ComputeChunkAtom::atom2bin3d()
1636 {
1637   int i,ibin,i1bin,i2bin,i3bin;
1638   double *boxlo,*boxhi,*prd;
1639   double xremap,yremap,zremap;
1640 
1641   double **x = atom->x;
1642   int nlocal = atom->nlocal;
1643 
1644   int idim = dim[0];
1645   int jdim = dim[1];
1646   int kdim = dim[2];
1647   int nlayer1m1 = nlayers[0] - 1;
1648   int nlayer2m1 = nlayers[1] - 1;
1649   int nlayer3m1 = nlayers[2] - 1;
1650   int *periodicity = domain->periodicity;
1651 
1652   if (periodicity[idim] || periodicity[jdim] || periodicity[kdim]) {
1653     if (scaleflag == REDUCED) {
1654       boxlo = domain->boxlo_lamda;
1655       boxhi = domain->boxhi_lamda;
1656       prd = domain->prd_lamda;
1657     } else {
1658       boxlo = domain->boxlo;
1659       boxhi = domain->boxhi;
1660       prd = domain->prd;
1661     }
1662   }
1663 
1664   // remap each atom's relevant coord back into box via PBC if necessary
1665   // if scaleflag = REDUCED, box coords -> lamda coords
1666   // apply discard rule
1667 
1668   if (scaleflag == REDUCED) domain->x2lamda(nlocal);
1669 
1670   for (i = 0; i < nlocal; i++) {
1671     if (exclude[i]) continue;
1672 
1673     xremap = x[i][idim];
1674     if (periodicity[idim]) {
1675       if (xremap < boxlo[idim]) xremap += prd[idim];
1676       if (xremap >= boxhi[idim]) xremap -= prd[idim];
1677     }
1678 
1679     i1bin = static_cast<int> ((xremap - offset[0]) * invdelta[0]);
1680     if (xremap < offset[0]) i1bin--;
1681 
1682     if (discard == MIXED) {
1683       if (!minflag[idim]) i1bin = MAX(i1bin,0);
1684       else if (i1bin < 0) {
1685         exclude[i] = 1;
1686         continue;
1687       }
1688       if (!maxflag[idim]) i1bin = MIN(i1bin,nlayer1m1);
1689       else if (i1bin > nlayer1m1) {
1690         exclude[i] = 1;
1691         continue;
1692       }
1693     } else if (discard == NODISCARD) {
1694       i1bin = MAX(i1bin,0);
1695       i1bin = MIN(i1bin,nlayer1m1);
1696     } else if (i1bin < 0 || i1bin > nlayer1m1) {
1697       exclude[i] = 1;
1698       continue;
1699     }
1700 
1701     yremap = x[i][jdim];
1702     if (periodicity[jdim]) {
1703       if (yremap < boxlo[jdim]) yremap += prd[jdim];
1704       if (yremap >= boxhi[jdim]) yremap -= prd[jdim];
1705     }
1706 
1707     i2bin = static_cast<int> ((yremap - offset[1]) * invdelta[1]);
1708     if (yremap < offset[1]) i2bin--;
1709 
1710     if (discard == MIXED) {
1711       if (!minflag[jdim]) i2bin = MAX(i2bin,0);
1712       else if (i2bin < 0) {
1713         exclude[i] = 1;
1714         continue;
1715       }
1716       if (!maxflag[jdim]) i2bin = MIN(i2bin,nlayer2m1);
1717       else if (i2bin > nlayer2m1) {
1718         exclude[i] = 1;
1719         continue;
1720       }
1721     } else if (discard == NODISCARD) {
1722       i2bin = MAX(i2bin,0);
1723       i2bin = MIN(i2bin,nlayer2m1);
1724     } else if (i2bin < 0 || i2bin > nlayer2m1) {
1725       exclude[i] = 1;
1726       continue;
1727     }
1728 
1729     zremap = x[i][kdim];
1730     if (periodicity[kdim]) {
1731       if (zremap < boxlo[kdim]) zremap += prd[kdim];
1732       if (zremap >= boxhi[kdim]) zremap -= prd[kdim];
1733     }
1734 
1735     i3bin = static_cast<int> ((zremap - offset[2]) * invdelta[2]);
1736     if (zremap < offset[2]) i3bin--;
1737 
1738     if (discard == MIXED) {
1739       if (!minflag[kdim]) i3bin = MAX(i3bin,0);
1740       else if (i3bin < 0) {
1741         exclude[i] = 1;
1742         continue;
1743       }
1744       if (!maxflag[kdim]) i3bin = MIN(i3bin,nlayer3m1);
1745       else if (i3bin > nlayer3m1) {
1746         exclude[i] = 1;
1747         continue;
1748       }
1749     } else if (discard == NODISCARD) {
1750       i3bin = MAX(i3bin,0);
1751       i3bin = MIN(i3bin,nlayer3m1);
1752     } else if (i3bin < 0 || i3bin > nlayer3m1) {
1753       exclude[i] = 1;
1754       continue;
1755     }
1756 
1757     ibin = i1bin*nlayers[1]*nlayers[2] + i2bin*nlayers[2] + i3bin;
1758     ichunk[i] = ibin+1;
1759   }
1760 
1761   if (scaleflag == REDUCED) domain->lamda2x(nlocal);
1762 }
1763 
1764 /* ----------------------------------------------------------------------
1765    assign each atom to a spherical bin
1766 ------------------------------------------------------------------------- */
1767 
atom2binsphere()1768 void ComputeChunkAtom::atom2binsphere()
1769 {
1770   int i,ibin;
1771   double dx,dy,dz,r;
1772   double xremap,yremap,zremap;
1773 
1774   double *boxlo = domain->boxlo;
1775   double *boxhi = domain->boxhi;
1776   double *prd = domain->prd;
1777   double *prd_half = domain->prd_half;
1778   int *periodicity = domain->periodicity;
1779 
1780   // remap each atom's relevant coords back into box via PBC if necessary
1781   // apply discard rule based on rmin and rmax
1782 
1783   double **x = atom->x;
1784   int nlocal = atom->nlocal;
1785 
1786   for (i = 0; i < nlocal; i++) {
1787     if (exclude[i]) continue;
1788 
1789     xremap = x[i][0];
1790     if (periodicity[0]) {
1791       while (xremap < boxlo[0]) {xremap += prd[0];}
1792       while (xremap >= boxhi[0]) {xremap -= prd[0];}
1793     }
1794     yremap = x[i][1];
1795     if (periodicity[1]) {
1796       while (yremap < boxlo[1]) {yremap += prd[1];}
1797       while (yremap >= boxhi[1]) {yremap -= prd[1];}
1798     }
1799     zremap = x[i][2];
1800     if (periodicity[2]) {
1801       while (zremap < boxlo[2]) {zremap += prd[2];}
1802       while (zremap >= boxhi[2]) {zremap -= prd[2];}
1803     }
1804 
1805     dx = xremap - sorigin[0];
1806     dy = yremap - sorigin[1];
1807     dz = zremap - sorigin[2];
1808 
1809     // if requested, apply PBC to distance from sphere center
1810     // treat orthogonal and triclinic the same
1811     //   with dx,dy,dz = lengths independent of each other
1812     // so do not use domain->minimum_image() which couples for triclinic
1813 
1814     if (pbcflag) {
1815       if (periodicity[0]) {
1816         while (fabs(dx) > prd_half[0]) {
1817           if (dx < 0.0) dx += prd[0];
1818           else dx -= prd[0];
1819         }
1820       }
1821       if (periodicity[1]) {
1822         while (fabs(dy) > prd_half[1]) {
1823           if (dy < 0.0) dy += prd[1];
1824           else dy -= prd[1];
1825         }
1826       }
1827       if (periodicity[2]) {
1828         while (fabs(dz) > prd_half[2]) {
1829           if (dz < 0.0) dz += prd[2];
1830           else dz -= prd[2];
1831         }
1832       }
1833     }
1834 
1835     r = sqrt(dx*dx + dy*dy + dz*dz);
1836 
1837     ibin = static_cast<int> ((r - sradmin) * sinvrad);
1838     if (r < sradmin) ibin--;
1839 
1840     if (discard == MIXED || discard == NODISCARD) {
1841       ibin = MAX(ibin,0);
1842       ibin = MIN(ibin,nchunk-1);
1843     } else if (ibin < 0 || ibin >= nchunk) {
1844       exclude[i] = 1;
1845       continue;
1846     }
1847 
1848     ichunk[i] = ibin+1;
1849   }
1850 }
1851 
1852 /* ----------------------------------------------------------------------
1853    assign each atom to a cylindrical bin
1854 ------------------------------------------------------------------------- */
1855 
atom2bincylinder()1856 void ComputeChunkAtom::atom2bincylinder()
1857 {
1858   int i,rbin,kbin;
1859   double d1,d2,r;
1860   double remap1,remap2;
1861 
1862   // first use atom2bin1d() to bin all atoms along cylinder axis
1863 
1864   atom2bin1d();
1865 
1866   // now bin in radial direction
1867   // kbin = bin along cylinder axis
1868   // rbin = bin in radial direction
1869 
1870   double *boxlo = domain->boxlo;
1871   double *boxhi = domain->boxhi;
1872   double *prd = domain->prd;
1873   double *prd_half = domain->prd_half;
1874   int *periodicity = domain->periodicity;
1875 
1876   // remap each atom's relevant coords back into box via PBC if necessary
1877   // apply discard rule based on rmin and rmax
1878 
1879   double **x = atom->x;
1880   int nlocal = atom->nlocal;
1881 
1882   for (i = 0; i < nlocal; i++) {
1883     if (exclude[i]) continue;
1884     kbin = ichunk[i] - 1;
1885 
1886     remap1 = x[i][cdim1];
1887     if (periodicity[cdim1]) {
1888       if (remap1 < boxlo[cdim1]) remap1 += prd[cdim1];
1889       if (remap1 >= boxhi[cdim1]) remap1 -= prd[cdim1];
1890     }
1891     remap2 = x[i][cdim2];
1892     if (periodicity[cdim2]) {
1893       if (remap2 < boxlo[cdim2]) remap2 += prd[cdim2];
1894       if (remap2 >= boxhi[cdim2]) remap2 -= prd[cdim2];
1895     }
1896 
1897     d1 = remap1 - corigin[cdim1];
1898     d2 = remap2 - corigin[cdim2];
1899 
1900     // if requested, apply PBC to distance from cylinder axis
1901     // treat orthogonal and triclinic the same
1902     //   with d1,d2 = lengths independent of each other
1903 
1904     if (pbcflag) {
1905       if (periodicity[cdim1]) {
1906         if (fabs(d1) > prd_half[cdim1]) {
1907           if (d1 < 0.0) d1 += prd[cdim1];
1908           else d1 -= prd[cdim1];
1909         }
1910       }
1911       if (periodicity[cdim2]) {
1912         if (fabs(d2) > prd_half[cdim2]) {
1913           if (d2 < 0.0) d2 += prd[cdim2];
1914           else d2 -= prd[cdim2];
1915         }
1916       }
1917     }
1918 
1919     r = sqrt(d1*d1 + d2*d2);
1920 
1921     rbin = static_cast<int> ((r - cradmin) * cinvrad);
1922     if (r < cradmin) rbin--;
1923 
1924     if (discard == MIXED || discard == NODISCARD) {
1925       rbin = MAX(rbin,0);
1926       rbin = MIN(rbin,ncbin-1);
1927     } else if (rbin < 0 || rbin >= ncbin) {
1928       exclude[i] = 1;
1929       continue;
1930     }
1931 
1932     // combine axis and radial bin indices to set ichunk
1933 
1934     ichunk[i] = rbin*ncplane + kbin + 1;
1935   }
1936 }
1937 
1938 /* ----------------------------------------------------------------------
1939    process args for one dimension of binning info
1940    set dim, originflag, origin, delta
1941 ------------------------------------------------------------------------- */
1942 
readdim(int narg,char ** arg,int iarg,int idim)1943 void ComputeChunkAtom::readdim(int narg, char **arg, int iarg, int idim)
1944 {
1945   if (narg < iarg+3) error->all(FLERR,"Illegal compute chunk/atom command");
1946   if (strcmp(arg[iarg],"x") == 0) dim[idim] = 0;
1947   else if (strcmp(arg[iarg],"y") == 0) dim[idim] = 1;
1948   else if (strcmp(arg[iarg],"z") == 0) dim[idim] = 2;
1949   else error->all(FLERR,"Illegal compute chunk/atom command");
1950 
1951   if (dim[idim] == 2 && domain->dimension == 2)
1952     error->all(FLERR,"Cannot use compute chunk/atom bin z for 2d model");
1953 
1954   if (strcmp(arg[iarg+1],"lower") == 0) originflag[idim] = LOWER;
1955   else if (strcmp(arg[iarg+1],"center") == 0) originflag[idim] = CENTER;
1956   else if (strcmp(arg[iarg+1],"upper") == 0) originflag[idim] = UPPER;
1957   else originflag[idim] = COORD;
1958   if (originflag[idim] == COORD)
1959     origin[idim] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
1960 
1961   delta[idim] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
1962 }
1963 
1964 /* ----------------------------------------------------------------------
1965    initialize one atom's storage values, called when atom is created
1966    just set chunkID to 0 for new atom
1967 ------------------------------------------------------------------------- */
1968 
set_arrays(int i)1969 void ComputeChunkAtom::set_arrays(int i)
1970 {
1971   if (!fixstore) return;
1972   double *vstore = fixstore->vstore;
1973   vstore[i] = 0.0;
1974 }
1975 
1976 /* ----------------------------------------------------------------------
1977    memory usage of local atom-based arrays and per-chunk arrays
1978    note: nchunk is actually 0 until first call
1979 ------------------------------------------------------------------------- */
1980 
memory_usage()1981 double ComputeChunkAtom::memory_usage()
1982 {
1983   double bytes = 2*MAX(nmaxint,0) * sizeof(int);   // ichunk,exclude
1984   bytes += (double)nmax * sizeof(double);                  // chunk
1985   bytes += (double)ncoord*nchunk * sizeof(double);         // coord
1986   if (compress) bytes += (double)nchunk * sizeof(int);     // chunkID
1987   return bytes;
1988 }
1989