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