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 /* ----------------------------------------------------------------------
16 Contributing author (triclinic and multi-neigh) : Pieter in 't Veld (SNL)
17 Contributing author (improved multi-neigh) : Joel Clemmer (SNL)
18 ------------------------------------------------------------------------- */
19
20 #include "neighbor.h"
21
22 #include "atom.h"
23 #include "atom_vec.h"
24 #include "citeme.h"
25 #include "comm.h"
26 #include "compute.h"
27 #include "domain.h"
28 #include "error.h"
29 #include "fix.h"
30 #include "force.h"
31 #include "group.h"
32 #include "memory.h"
33 #include "modify.h"
34 #include "nbin.h"
35 #include "neigh_list.h"
36 #include "neigh_request.h"
37 #include "npair.h"
38 #include "nstencil.h"
39 #include "ntopo.h"
40 #include "output.h"
41 #include "pair.h"
42 #include "pair_hybrid.h"
43 #include "respa.h"
44 #include "style_nbin.h" // IWYU pragma: keep
45 #include "style_npair.h" // IWYU pragma: keep
46 #include "style_nstencil.h" // IWYU pragma: keep
47 #include "style_ntopo.h" // IWYU pragma: keep
48 #include "tokenizer.h"
49 #include "update.h"
50
51 #include <cmath>
52 #include <cstring>
53
54 using namespace LAMMPS_NS;
55 using namespace NeighConst;
56
57 #define RQDELTA 1
58 #define EXDELTA 1
59 #define DELTA_PERATOM 64
60
61 #define BIG 1.0e20
62
63 enum{NONE,ALL,PARTIAL,TEMPLATE};
64
65 static const char cite_neigh_multi_old[] =
66 "neighbor multi/old command: doi:10.1016/j.cpc.2008.03.005\n\n"
67 "@Article{Intveld08,\n"
68 " author = {P.{\\,}J.~in{\\,}'t~Veld and S.{\\,}J.~Plimpton"
69 " and G.{\\,}S.~Grest},\n"
70 " title = {Accurate and Efficient Methods for Modeling Colloidal\n"
71 " Mixtures in an Explicit Solvent using Molecular Dynamics},\n"
72 " journal = {Comp.~Phys.~Comm.},\n"
73 " year = 2008,\n"
74 " volume = 179,\n"
75 " pages = {320--329}\n"
76 "}\n\n";
77
78 static const char cite_neigh_multi[] =
79 "neighbor multi command: doi:10.1016/j.cpc.2008.03.005, doi:10.1007/s40571-020-00361-2\n\n"
80 "@Article{Intveld08,\n"
81 " author = {P.{\\,}J.~in{\\,}'t~Veld and S.{\\,}J.~Plimpton"
82 " and G.{\\,}S.~Grest},\n"
83 " title = {Accurate and Efficient Methods for Modeling Colloidal\n"
84 " Mixtures in an Explicit Solvent using Molecular Dynamics},\n"
85 " journal = {Comp.~Phys.~Comm.},\n"
86 " year = 2008,\n"
87 " volume = 179,\n"
88 " pages = {320--329}\n"
89 "}\n\n"
90 "@article{Stratford2018,\n"
91 " author = {Stratford, Kevin and Shire, Tom and Hanley, Kevin},\n"
92 " title = {Implementation of multi-level contact detection in LAMMPS},\n"
93 " year = {2018}\n"
94 "}\n\n"
95 "@article{Shire2020,\n"
96 " author = {Shire, Tom and Hanley, Kevin J. and Stratford, Kevin},\n"
97 " title = {DEM simulations of polydisperse media: efficient contact\n"
98 " detection applied to investigate the quasi-static limit},\n"
99 " journal = {Computational Particle Mechanics},\n"
100 " year = {2020}\n"
101 "}\n\n";
102
103 //#define NEIGH_LIST_DEBUG 1
104
105 /* ---------------------------------------------------------------------- */
106
Neighbor(LAMMPS * lmp)107 Neighbor::Neighbor(LAMMPS *lmp) : Pointers(lmp),
108 pairclass(nullptr), pairnames(nullptr), pairmasks(nullptr)
109 {
110 MPI_Comm_rank(world,&me);
111 MPI_Comm_size(world,&nprocs);
112
113 firsttime = 1;
114
115 style = Neighbor::BIN;
116 every = 1;
117 delay = 10;
118 dist_check = 1;
119 pgsize = 100000;
120 oneatom = 2000;
121 binsizeflag = 0;
122 build_once = 0;
123 cluster_check = 0;
124 ago = -1;
125
126 cutneighmax = 0.0;
127 cutneighsq = nullptr;
128 cutneighghostsq = nullptr;
129 cuttype = nullptr;
130 cuttypesq = nullptr;
131 fixchecklist = nullptr;
132
133 // pairwise neighbor lists and associated data structs
134
135 nlist = 0;
136 lists = nullptr;
137
138 nbin = 0;
139 neigh_bin = nullptr;
140
141 nstencil = 0;
142 neigh_stencil = nullptr;
143
144 neigh_pair = nullptr;
145
146 nstencil_perpetual = 0;
147 slist = nullptr;
148
149 npair_perpetual = 0;
150 plist = nullptr;
151
152 nrequest = maxrequest = 0;
153 requests = nullptr;
154
155 old_nrequest = 0;
156 old_requests = nullptr;
157
158 old_style = style;
159 old_triclinic = 0;
160 old_pgsize = pgsize;
161 old_oneatom = oneatom;
162
163 binclass = nullptr;
164 binnames = nullptr;
165 binmasks = nullptr;
166 stencilclass = nullptr;
167 stencilnames = nullptr;
168 stencilmasks = nullptr;
169
170 // topology lists
171
172 bondwhich = anglewhich = dihedralwhich = improperwhich = NONE;
173
174 neigh_bond = nullptr;
175 neigh_angle = nullptr;
176 neigh_dihedral = nullptr;
177 neigh_improper = nullptr;
178
179 // coords at last neighboring
180
181 maxhold = 0;
182 xhold = nullptr;
183 lastcall = -1;
184 last_setup_bins = -1;
185
186 // pair exclusion list info
187
188 includegroup = 0;
189
190 nex_type = maxex_type = 0;
191 ex1_type = ex2_type = nullptr;
192 ex_type = nullptr;
193
194 nex_group = maxex_group = 0;
195 ex1_group = ex2_group = ex1_bit = ex2_bit = nullptr;
196
197 nex_mol = maxex_mol = 0;
198 ex_mol_group = ex_mol_bit = ex_mol_intra = nullptr;
199
200 // Multi data
201
202 type2collection = nullptr;
203 collection2cut = nullptr;
204 collection = nullptr;
205 cutcollectionsq = nullptr;
206 custom_collection_flag = 0;
207 interval_collection_flag = 0;
208 nmax_collection = 0;
209
210 // Kokkos setting
211
212 copymode = 0;
213 }
214
215 /* ---------------------------------------------------------------------- */
216
~Neighbor()217 Neighbor::~Neighbor()
218 {
219 if (copymode) return;
220
221 memory->destroy(cutneighsq);
222 memory->destroy(cutneighghostsq);
223 delete [] cuttype;
224 delete [] cuttypesq;
225 delete [] fixchecklist;
226
227 for (int i = 0; i < nlist; i++) delete lists[i];
228 for (int i = 0; i < nbin; i++) delete neigh_bin[i];
229 for (int i = 0; i < nstencil; i++) delete neigh_stencil[i];
230 for (int i = 0; i < nlist; i++) delete neigh_pair[i];
231 delete [] lists;
232 delete [] neigh_bin;
233 delete [] neigh_stencil;
234 delete [] neigh_pair;
235
236 delete [] slist;
237 delete [] plist;
238
239 for (int i = 0; i < nrequest; i++)
240 if (requests[i]) delete requests[i];
241 memory->sfree(requests);
242 for (int i = 0; i < old_nrequest; i++)
243 if (old_requests[i]) delete old_requests[i];
244 memory->sfree(old_requests);
245
246 delete [] binclass;
247 delete [] binnames;
248 delete [] binmasks;
249 delete [] stencilclass;
250 delete [] stencilnames;
251 delete [] stencilmasks;
252 delete [] pairclass;
253 delete [] pairnames;
254 delete [] pairmasks;
255
256 delete neigh_bond;
257 delete neigh_angle;
258 delete neigh_dihedral;
259 delete neigh_improper;
260
261 memory->destroy(xhold);
262
263 memory->destroy(ex1_type);
264 memory->destroy(ex2_type);
265 memory->destroy(ex_type);
266
267 memory->destroy(ex1_group);
268 memory->destroy(ex2_group);
269 delete [] ex1_bit;
270 delete [] ex2_bit;
271
272 memory->destroy(ex_mol_group);
273 delete [] ex_mol_bit;
274 memory->destroy(ex_mol_intra);
275
276 memory->destroy(type2collection);
277 memory->destroy(collection2cut);
278 memory->destroy(collection);
279 memory->destroy(cutcollectionsq);
280 }
281
282 /* ---------------------------------------------------------------------- */
283
init()284 void Neighbor::init()
285 {
286 int i,j,n;
287
288 ncalls = ndanger = 0;
289 dimension = domain->dimension;
290 triclinic = domain->triclinic;
291 newton_pair = force->newton_pair;
292
293 // error check
294
295 if (delay > 0 && (delay % every) != 0)
296 error->all(FLERR,"Neighbor delay must be 0 or multiple of every setting");
297
298 if (pgsize < 10*oneatom)
299 error->all(FLERR,"Neighbor page size must be >= 10x the one atom setting");
300
301 // ------------------------------------------------------------------
302 // settings
303
304 // bbox lo/hi ptrs = bounding box of entire domain, stored by Domain
305
306 if (triclinic == 0) {
307 bboxlo = domain->boxlo;
308 bboxhi = domain->boxhi;
309 } else {
310 bboxlo = domain->boxlo_bound;
311 bboxhi = domain->boxhi_bound;
312 }
313
314 // set neighbor cutoffs (force cutoff + skin)
315 // trigger determines when atoms migrate and neighbor lists are rebuilt
316 // needs to be non-zero for migration distance check
317 // even if pair = nullptr and no neighbor lists are used
318 // cutneigh = force cutoff + skin if cutforce > 0, else cutneigh = 0
319 // cutneighghost = pair cutghost if it requests it, else same as cutneigh
320
321 triggersq = 0.25*skin*skin;
322 boxcheck = 0;
323 if (domain->box_change && (domain->xperiodic || domain->yperiodic ||
324 (dimension == 3 && domain->zperiodic)))
325 boxcheck = 1;
326
327 n = atom->ntypes;
328 if (cutneighsq == nullptr) {
329 if (lmp->kokkos) init_cutneighsq_kokkos(n);
330 else memory->create(cutneighsq,n+1,n+1,"neigh:cutneighsq");
331 memory->create(cutneighghostsq,n+1,n+1,"neigh:cutneighghostsq");
332 cuttype = new double[n+1];
333 cuttypesq = new double[n+1];
334 }
335
336 double cutoff,delta,cut;
337 cutneighmin = BIG;
338 cutneighmax = 0.0;
339
340 for (i = 1; i <= n; i++) {
341 cuttype[i] = cuttypesq[i] = 0.0;
342 for (j = 1; j <= n; j++) {
343 if (force->pair) cutoff = sqrt(force->pair->cutsq[i][j]);
344 else cutoff = 0.0;
345 if (cutoff > 0.0) delta = skin;
346 else delta = 0.0;
347 cut = cutoff + delta;
348
349 cutneighsq[i][j] = cut*cut;
350 cuttype[i] = MAX(cuttype[i],cut);
351 cuttypesq[i] = MAX(cuttypesq[i],cut*cut);
352 cutneighmin = MIN(cutneighmin,cut);
353 cutneighmax = MAX(cutneighmax,cut);
354
355 if (force->pair && force->pair->ghostneigh) {
356 cut = force->pair->cutghost[i][j] + skin;
357 cutneighghostsq[i][j] = cut*cut;
358 } else cutneighghostsq[i][j] = cut*cut;
359 }
360 }
361 cutneighmaxsq = cutneighmax * cutneighmax;
362
363 // Define cutoffs for multi
364 if (style == Neighbor::MULTI) {
365 int icollection, jcollection;
366
367 // If collections not yet defined, create default map using types
368 if (!custom_collection_flag) {
369 ncollections = n;
370 interval_collection_flag = 0;
371 if (!type2collection)
372 memory->create(type2collection,n+1,"neigh:type2collection");
373 for (i = 1; i <= n; i++)
374 type2collection[i] = i-1;
375 }
376
377 memory->grow(cutcollectionsq, ncollections, ncollections, "neigh:cutcollectionsq");
378
379 // 3 possible ways of defining collections
380 // 1) Types are used to define collections
381 // Each collection loops through its owned types, and uses cutneighsq to calculate its cutoff
382 // 2) Collections are defined by intervals, point particles
383 // Types are first sorted into collections based on cutneighsq[i][i]
384 // Each collection loops through its owned types, and uses cutneighsq to calculate its cutoff
385 // 3) Collections are defined by intervals, finite particles
386 //
387
388 // Define collection cutoffs
389 for (i = 0; i < ncollections; i++)
390 for (j = 0; j < ncollections; j++)
391 cutcollectionsq[i][j] = 0.0;
392
393 if (!interval_collection_flag) {
394 finite_cut_flag = 0;
395 for (i = 1; i <= n; i++){
396 icollection = type2collection[i];
397 for (j = 1; j <= n; j++){
398 jcollection = type2collection[j];
399 if (cutneighsq[i][j] > cutcollectionsq[icollection][jcollection]) {
400 cutcollectionsq[icollection][jcollection] = cutneighsq[i][j];
401 cutcollectionsq[jcollection][icollection] = cutneighsq[i][j];
402 }
403 }
404 }
405 } else {
406 if (force->pair->finitecutflag) {
407 finite_cut_flag = 1;
408 // If cutoffs depend on finite atom sizes, use radii of intervals to find cutoffs
409 double ri, rj, tmp;
410 for (i = 0; i < ncollections; i++){
411 ri = collection2cut[i]*0.5;
412 for (j = 0; j < ncollections; j++){
413 rj = collection2cut[j]*0.5;
414 tmp = force->pair->radii2cut(ri, rj) + skin;
415 cutcollectionsq[i][j] = tmp*tmp;
416 }
417 }
418 } else {
419 finite_cut_flag = 0;
420
421 // Map types to collections
422 if (!type2collection)
423 memory->create(type2collection,n+1,"neigh:type2collection");
424
425 for (i = 1; i <= n; i++)
426 type2collection[i] = -1;
427
428 double cuttmp;
429 for (i = 1; i <= n; i++){
430 // Remove skin added to cutneighsq
431 cuttmp = sqrt(cutneighsq[i][i]) - skin;
432 for (icollection = 0; icollection < ncollections; icollection ++){
433 if (collection2cut[icollection] >= cuttmp) {
434 type2collection[i] = icollection;
435 break;
436 }
437 }
438
439 if (type2collection[i] == -1)
440 error->all(FLERR, "Pair cutoff exceeds interval cutoffs for multi");
441 }
442
443 // Define cutoffs
444 for (i = 1; i <= n; i++){
445 icollection = type2collection[i];
446 for (j = 1; j <= n; j++){
447 jcollection = type2collection[j];
448 if (cutneighsq[i][j] > cutcollectionsq[icollection][jcollection]) {
449 cutcollectionsq[icollection][jcollection] = cutneighsq[i][j];
450 cutcollectionsq[jcollection][icollection] = cutneighsq[i][j];
451 }
452 }
453 }
454 }
455 }
456 }
457
458 // rRESPA cutoffs
459
460 int respa = 0;
461 if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
462 if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
463 if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
464 }
465
466 if (respa) {
467 double *cut_respa = ((Respa *) update->integrate)->cutoff;
468 cut_inner_sq = (cut_respa[1] + skin) * (cut_respa[1] + skin);
469 cut_middle_sq = (cut_respa[3] + skin) * (cut_respa[3] + skin);
470 cut_middle_inside_sq = (cut_respa[0] - skin) * (cut_respa[0] - skin);
471 if (cut_respa[0]-skin < 0) cut_middle_inside_sq = 0.0;
472 }
473
474 // fixchecklist = other classes that can induce reneighboring in decide()
475
476 restart_check = 0;
477 if (output->restart_flag) restart_check = 1;
478
479 delete [] fixchecklist;
480 fixchecklist = nullptr;
481 fixchecklist = new int[modify->nfix];
482
483 fix_check = 0;
484 for (i = 0; i < modify->nfix; i++)
485 if (modify->fix[i]->force_reneighbor)
486 fixchecklist[fix_check++] = i;
487
488 must_check = 0;
489 if (restart_check || fix_check) must_check = 1;
490
491 // set special_flag for 1-2, 1-3, 1-4 neighbors
492 // flag[0] is not used, flag[1] = 1-2, flag[2] = 1-3, flag[3] = 1-4
493 // flag = 0 if both LJ/Coulomb special values are 0.0
494 // flag = 1 if both LJ/Coulomb special values are 1.0
495 // flag = 2 otherwise or if KSpace solver is enabled
496 // pairwise portion of KSpace solver uses all 1-2,1-3,1-4 neighbors
497 // or selected Coulomb-approixmation pair styles require it
498
499 if (force->special_lj[1] == 0.0 && force->special_coul[1] == 0.0)
500 special_flag[1] = 0;
501 else if (force->special_lj[1] == 1.0 && force->special_coul[1] == 1.0)
502 special_flag[1] = 1;
503 else special_flag[1] = 2;
504
505 if (force->special_lj[2] == 0.0 && force->special_coul[2] == 0.0)
506 special_flag[2] = 0;
507 else if (force->special_lj[2] == 1.0 && force->special_coul[2] == 1.0)
508 special_flag[2] = 1;
509 else special_flag[2] = 2;
510
511 if (force->special_lj[3] == 0.0 && force->special_coul[3] == 0.0)
512 special_flag[3] = 0;
513 else if (force->special_lj[3] == 1.0 && force->special_coul[3] == 1.0)
514 special_flag[3] = 1;
515 else special_flag[3] = 2;
516
517 // We cannot remove special neighbors with kspace or kspace-like pair styles
518 // as the exclusion needs to remove the full coulomb and not the damped interaction.
519 // Special treatment is required for hybrid pair styles since Force::pair_match()
520 // will only return a non-null pointer if there is only one substyle of the kind.
521
522 if (force->kspace) {
523 special_flag[1] = special_flag[2] = special_flag[3] = 2;
524 } else {
525 PairHybrid *ph = reinterpret_cast<PairHybrid *>(force->pair_match("^hybrid",0));
526 if (ph) {
527 int flag=0;
528 for (int isub=0; isub < ph->nstyles; ++isub) {
529 if (force->pair_match("coul/wolf",0,isub)
530 || force->pair_match("coul/dsf",0,isub)
531 || force->pair_match("coul/exclude",0)
532 || force->pair_match("thole",0,isub))
533 ++flag;
534 }
535 if (flag)
536 special_flag[1] = special_flag[2] = special_flag[3] = 2;
537 } else {
538 if (force->pair_match("coul/wolf",0)
539 || force->pair_match("coul/dsf",0)
540 || force->pair_match("coul/exclude",0)
541 || force->pair_match("thole",0))
542 special_flag[1] = special_flag[2] = special_flag[3] = 2;
543 }
544 }
545
546 // ------------------------------------------------------------------
547 // xhold array
548
549 // free if not needed for this run
550
551 if (dist_check == 0) {
552 memory->destroy(xhold);
553 maxhold = 0;
554 xhold = nullptr;
555 }
556
557 // first time allocation
558
559 if (dist_check) {
560 if (maxhold == 0) {
561 maxhold = atom->nmax;
562 memory->create(xhold,maxhold,3,"neigh:xhold");
563 }
564 }
565
566 // ------------------------------------------------------------------
567 // exclusion lists
568
569 // depend on type, group, molecule settings from neigh_modify
570 // warn if exclusions used with KSpace solver
571
572 n = atom->ntypes;
573
574 if (nex_type == 0 && nex_group == 0 && nex_mol == 0) exclude = 0;
575 else exclude = 1;
576
577 if (nex_type) {
578 if (lmp->kokkos)
579 init_ex_type_kokkos(n);
580 else {
581 memory->destroy(ex_type);
582 memory->create(ex_type,n+1,n+1,"neigh:ex_type");
583 }
584
585 for (i = 1; i <= n; i++)
586 for (j = 1; j <= n; j++)
587 ex_type[i][j] = 0;
588
589 for (i = 0; i < nex_type; i++) {
590 if (ex1_type[i] <= 0 || ex1_type[i] > n ||
591 ex2_type[i] <= 0 || ex2_type[i] > n)
592 error->all(FLERR,"Invalid atom type in neighbor exclusion list");
593 ex_type[ex1_type[i]][ex2_type[i]] = 1;
594 ex_type[ex2_type[i]][ex1_type[i]] = 1;
595 }
596 }
597
598 if (nex_group) {
599 if (lmp->kokkos)
600 init_ex_bit_kokkos();
601 else {
602 delete [] ex1_bit;
603 delete [] ex2_bit;
604 ex1_bit = new int[nex_group];
605 ex2_bit = new int[nex_group];
606 }
607
608 for (i = 0; i < nex_group; i++) {
609 ex1_bit[i] = group->bitmask[ex1_group[i]];
610 ex2_bit[i] = group->bitmask[ex2_group[i]];
611 }
612 }
613
614 if (nex_mol) {
615 if (lmp->kokkos)
616 init_ex_mol_bit_kokkos();
617 else {
618 delete [] ex_mol_bit;
619 ex_mol_bit = new int[nex_mol];
620 }
621
622 for (i = 0; i < nex_mol; i++)
623 ex_mol_bit[i] = group->bitmask[ex_mol_group[i]];
624 }
625
626 if (exclude && force->kspace && me == 0)
627 error->warning(FLERR,"Neighbor exclusions used with KSpace solver "
628 "may give inconsistent Coulombic energies");
629
630 if (lmp->kokkos)
631 set_binsize_kokkos();
632
633 // ------------------------------------------------------------------
634 // create pairwise lists
635 // one-time call to init_styles() to scan style files and setup
636 // init_pair() creates auxiliary classes: NBin, NStencil, NPair
637
638 if (firsttime) init_styles();
639 firsttime = 0;
640
641 int same = init_pair();
642
643 // invoke copy_neighbor_info() in Bin,Stencil,Pair classes
644 // copied once per run in case any cutoff, exclusion, special info changed
645
646 for (i = 0; i < nbin; i++) neigh_bin[i]->copy_neighbor_info();
647 for (i = 0; i < nstencil; i++) neigh_stencil[i]->copy_neighbor_info();
648 for (i = 0; i < nlist; i++)
649 if (neigh_pair[i]) neigh_pair[i]->copy_neighbor_info();
650
651 if (!same && comm->me == 0) print_pairwise_info();
652
653 // can now delete requests so next run can make new ones
654 // print_pairwise_info() made use of requests
655 // set of NeighLists now stores all needed info
656
657 for (i = 0; i < nrequest; i++) {
658 delete requests[i];
659 requests[i] = nullptr;
660 }
661 nrequest = 0;
662
663 // ------------------------------------------------------------------
664 // create topology lists
665 // instantiated topo styles can change from run to run
666
667 init_topology();
668 }
669
670 /* ----------------------------------------------------------------------
671 create and initialize lists of Nbin, Nstencil, NPair classes
672 lists have info on all classes in 3 style*.h files
673 cannot do this in constructor, b/c too early to instantiate classes
674 ------------------------------------------------------------------------- */
675
init_styles()676 void Neighbor::init_styles()
677 {
678 // extract info from NBin classes listed in style_nbin.h
679
680 nbclass = 0;
681
682 #define NBIN_CLASS
683 #define NBinStyle(key,Class,bitmasks) nbclass++;
684 #include "style_nbin.h" // IWYU pragma: keep
685 #undef NBinStyle
686 #undef NBIN_CLASS
687
688 binclass = new BinCreator[nbclass];
689 binnames = new char*[nbclass];
690 binmasks = new int[nbclass];
691 nbclass = 0;
692
693 #define NBIN_CLASS
694 #define NBinStyle(key,Class,bitmasks) \
695 binnames[nbclass] = (char *) #key; \
696 binclass[nbclass] = &bin_creator<Class>; \
697 binmasks[nbclass++] = bitmasks;
698 #include "style_nbin.h" // IWYU pragma: keep
699 #undef NBinStyle
700 #undef NBIN_CLASS
701
702 // extract info from NStencil classes listed in style_nstencil.h
703
704 nsclass = 0;
705
706 #define NSTENCIL_CLASS
707 #define NStencilStyle(key,Class,bitmasks) nsclass++;
708 #include "style_nstencil.h" // IWYU pragma: keep
709 #undef NStencilStyle
710 #undef NSTENCIL_CLASS
711
712 stencilclass = new StencilCreator[nsclass];
713 stencilnames = new char*[nsclass];
714 stencilmasks = new int[nsclass];
715 nsclass = 0;
716
717 #define NSTENCIL_CLASS
718 #define NStencilStyle(key,Class,bitmasks) \
719 stencilnames[nsclass] = (char *) #key; \
720 stencilclass[nsclass] = &stencil_creator<Class>; \
721 stencilmasks[nsclass++] = bitmasks;
722 #include "style_nstencil.h" // IWYU pragma: keep
723 #undef NStencilStyle
724 #undef NSTENCIL_CLASS
725
726 // extract info from NPair classes listed in style_npair.h
727
728 npclass = 0;
729
730 #define NPAIR_CLASS
731 #define NPairStyle(key,Class,bitmasks) npclass++;
732 #include "style_npair.h" // IWYU pragma: keep
733 #undef NPairStyle
734 #undef NPAIR_CLASS
735
736 pairclass = new PairCreator[npclass];
737 pairnames = new char*[npclass];
738 pairmasks = new int[npclass];
739 npclass = 0;
740
741 #define NPAIR_CLASS
742 #define NPairStyle(key,Class,bitmasks) \
743 pairnames[npclass] = (char *) #key; \
744 pairclass[npclass] = &pair_creator<Class>; \
745 pairmasks[npclass++] = bitmasks;
746 #include "style_npair.h" // IWYU pragma: keep
747 #undef NPairStyle
748 #undef NPAIR_CLASS
749 }
750
751 /* ----------------------------------------------------------------------
752 create and initialize NPair classes
753 ------------------------------------------------------------------------- */
754
init_pair()755 int Neighbor::init_pair()
756 {
757 int i,j,k,m;
758
759 // test if pairwise lists need to be re-created
760 // no need to re-create if:
761 // neigh style, triclinic, pgsize, oneatom have not changed
762 // current requests = old requests
763 // so just return:
764 // delete requests so next run can make new ones
765 // current set of NeighLists already stores all needed info
766 // requests are compared via identical() before:
767 // any requests are morphed using logic below
768 // any requests are added below, e.g. as parents of pair hybrid skip lists
769 // copy them via requests_new2old() BEFORE any changes made to requests
770 // necessary b/c morphs can change requestor settings (see comment below)
771
772 int same = 1;
773 if (style != old_style) same = 0;
774 if (triclinic != old_triclinic) same = 0;
775 if (pgsize != old_pgsize) same = 0;
776 if (oneatom != old_oneatom) same = 0;
777
778 if (nrequest != old_nrequest) same = 0;
779 else
780 for (i = 0; i < nrequest; i++)
781 if (requests[i]->identical(old_requests[i]) == 0) same = 0;
782
783 #ifdef NEIGH_LIST_DEBUG
784 if (comm->me == 0) printf("SAME flag %d\n",same);
785 #endif
786
787 if (same) return same;
788 requests_new2old();
789
790 // delete old lists since creating new ones
791
792 for (i = 0; i < nlist; i++) delete lists[i];
793 for (i = 0; i < nbin; i++) delete neigh_bin[i];
794 for (i = 0; i < nstencil; i++) delete neigh_stencil[i];
795 for (i = 0; i < nlist; i++) delete neigh_pair[i];
796 delete [] lists;
797 delete [] neigh_bin;
798 delete [] neigh_stencil;
799 delete [] neigh_pair;
800
801 // error check on requests
802 // do not allow occasional, ghost, bin list
803 // b/c it still uses variant of coord2bin() in NPair() method
804 // instead of atom2bin, this could cause error b/c stoms have
805 // moved out of proc domain by time occasional list is built
806 // solution would be to use a different NBin variant
807 // that used Npair::coord2bin(x,ix,iy,iz) (then delete it from NPair)
808 // and stored the ix,iy,iz values for all atoms (including ghosts)
809 // at time of binning when neighbor lists are rebuilt,
810 // similar to what vanilla Nbin::coord2atom() does now in atom2bin
811
812 if (style == Neighbor::BIN) {
813 for (i = 0; i < nrequest; i++)
814 if (requests[i]->occasional && requests[i]->ghost)
815 error->all(FLERR,"Cannot request an occasional binned neighbor list "
816 "with ghost info");
817 }
818
819 // morph requests in various ways
820 // purpose is to avoid duplicate or inefficient builds
821 // may add new requests if a needed request to derive from does not exist
822 // methods:
823 // (1) unique = create unique lists if cutoff is explicitly set
824 // (2) skip = create any new non-skip lists needed by pair hybrid skip lists
825 // (3) granular = adjust parent and skip lists for granular onesided usage
826 // (4) h/f = pair up any matching half/full lists
827 // (5) copy = convert as many lists as possible to copy lists
828 // order of morph methods matters:
829 // (3) after (2), b/c it adjusts lists created by (2)
830 // (4) after (2) and (3),
831 // b/c (2) may create new full lists, (3) may change them
832 // (5) last, after all lists are finalized, so all possible copies found
833
834 int nrequest_original = nrequest;
835
836 morph_unique();
837 morph_skip();
838 morph_granular(); // this method can change flags set by requestor
839 morph_halffull();
840 morph_copy();
841
842 // create new lists, one per request including added requests
843 // wait to allocate initial pages until copy lists are detected
844 // NOTE: can I allocate now, instead of down below?
845
846 nlist = nrequest;
847
848 lists = new NeighList*[nrequest];
849 neigh_bin = new NBin*[nrequest];
850 neigh_stencil = new NStencil*[nrequest];
851 neigh_pair = new NPair*[nrequest];
852
853 // allocate new lists
854 // pass list ptr back to requestor (except for Command class)
855 // only for original requests, not ones added by Neighbor class
856
857 for (i = 0; i < nrequest; i++) {
858 if (requests[i]->kokkos_host || requests[i]->kokkos_device)
859 create_kokkos_list(i);
860 else lists[i] = new NeighList(lmp);
861 lists[i]->index = i;
862 lists[i]->requestor = requests[i]->requestor;
863
864 if (requests[i]->pair) {
865 lists[i]->requestor_type = NeighList::PAIR;
866 } else if (requests[i]->fix) {
867 lists[i]->requestor_type = NeighList::FIX;
868 } else if (requests[i]->compute) {
869 lists[i]->requestor_type = NeighList::COMPUTE;
870 }
871
872 if (requests[i]->pair && i < nrequest_original) {
873 Pair *pair = (Pair *) requests[i]->requestor;
874 pair->init_list(requests[i]->id,lists[i]);
875 } else if (requests[i]->fix && i < nrequest_original) {
876 Fix *fix = (Fix *) requests[i]->requestor;
877 fix->init_list(requests[i]->id,lists[i]);
878 } else if (requests[i]->compute && i < nrequest_original) {
879 Compute *compute = (Compute *) requests[i]->requestor;
880 compute->init_list(requests[i]->id,lists[i]);
881 }
882 }
883
884 // invoke post_constructor() for all lists
885 // copies info from requests to lists, sets ptrs to related lists
886
887 for (i = 0; i < nrequest; i++)
888 lists[i]->post_constructor(requests[i]);
889
890 // assign Bin,Stencil,Pair style to each list
891
892 int flag;
893 for (i = 0; i < nrequest; i++) {
894 flag = choose_bin(requests[i]);
895 lists[i]->bin_method = flag;
896 if (flag < 0)
897 error->all(FLERR,"Requested neighbor bin option does not exist");
898
899 flag = choose_stencil(requests[i]);
900 lists[i]->stencil_method = flag;
901 if (flag < 0)
902 error->all(FLERR,"Requested neighbor stencil method does not exist");
903
904 flag = choose_pair(requests[i]);
905 lists[i]->pair_method = flag;
906 if (flag < 0)
907 error->all(FLERR,"Requested neighbor pair method does not exist");
908 }
909
910 // instantiate unique Bin,Stencil classes in neigh_bin & neigh_stencil vecs
911 // unique = only one of its style, or request unique flag set (custom cutoff)
912
913 nbin = 0;
914 for (i = 0; i < nrequest; i++) {
915 requests[i]->index_bin = -1;
916 flag = lists[i]->bin_method;
917 if (flag == 0) continue;
918 if (!requests[i]->unique) {
919 for (j = 0; j < nbin; j++)
920 if (neigh_bin[j]->istyle == flag &&
921 neigh_bin[j]->cutoff_custom == 0.0) break;
922 if (j < nbin) {
923 requests[i]->index_bin = j;
924 continue;
925 }
926 }
927
928 BinCreator &bin_creator = binclass[flag-1];
929 neigh_bin[nbin] = bin_creator(lmp);
930 neigh_bin[nbin]->post_constructor(requests[i]);
931 neigh_bin[nbin]->istyle = flag;
932
933 requests[i]->index_bin = nbin;
934 nbin++;
935 }
936
937 nstencil = 0;
938 for (i = 0; i < nrequest; i++) {
939 requests[i]->index_stencil = -1;
940 flag = lists[i]->stencil_method;
941 if (flag == 0) continue;
942 if (!requests[i]->unique) {
943 for (j = 0; j < nstencil; j++)
944 if (neigh_stencil[j]->istyle == flag &&
945 neigh_stencil[j]->cutoff_custom == 0.0) break;
946 if (j < nstencil) {
947 requests[i]->index_stencil = j;
948 continue;
949 }
950 }
951
952 StencilCreator &stencil_creator = stencilclass[flag-1];
953 neigh_stencil[nstencil] = stencil_creator(lmp);
954 neigh_stencil[nstencil]->post_constructor(requests[i]);
955 neigh_stencil[nstencil]->istyle = flag;
956
957 if (lists[i]->bin_method > 0) {
958 neigh_stencil[nstencil]->nb = neigh_bin[requests[i]->index_bin];
959 if (neigh_stencil[nstencil]->nb == nullptr)
960 error->all(FLERR,"Could not assign bin method to neighbor stencil");
961 }
962
963 requests[i]->index_stencil = nstencil;
964 nstencil++;
965 }
966
967 // instantiate one Pair class per list in neigh_pair vec
968
969 for (i = 0; i < nrequest; i++) {
970 requests[i]->index_pair = -1;
971 flag = lists[i]->pair_method;
972 if (flag == 0) {
973 neigh_pair[i] = nullptr;
974 continue;
975 }
976
977 PairCreator &pair_creator = pairclass[flag-1];
978 lists[i]->np = neigh_pair[i] = pair_creator(lmp);
979 neigh_pair[i]->post_constructor(requests[i]);
980 neigh_pair[i]->istyle = flag;
981
982 if (lists[i]->bin_method > 0) {
983 neigh_pair[i]->nb = neigh_bin[requests[i]->index_bin];
984 if (neigh_pair[i]->nb == nullptr)
985 error->all(FLERR,"Could not assign bin method to neighbor pair");
986 }
987 if (lists[i]->stencil_method > 0) {
988 neigh_pair[i]->ns = neigh_stencil[requests[i]->index_stencil];
989 if (neigh_pair[i]->ns == nullptr)
990 error->all(FLERR,"Could not assign stencil method to neighbor pair");
991 }
992
993 requests[i]->index_pair = i;
994 }
995
996 // allocate initial pages for each list, except if copy flag set
997
998 for (i = 0; i < nlist; i++) {
999 if (lists[i]->copy && !lists[i]->kk2cpu)
1000 continue;
1001 lists[i]->setup_pages(pgsize,oneatom);
1002 }
1003
1004 // first-time allocation of per-atom data for lists that are built and store
1005 // lists that do not store: copy
1006 // use atom->nmax for both grow() args
1007 // i.e. grow first time to expanded size to avoid future reallocs
1008 // also Kokkos list initialization
1009
1010 int maxatom = atom->nmax;
1011 for (i = 0; i < nlist; i++) {
1012 if (neigh_pair[i] && (!lists[i]->copy || lists[i]->kk2cpu))
1013 lists[i]->grow(maxatom,maxatom);
1014 }
1015
1016 // plist = indices of perpetual NPair classes
1017 // perpetual = non-occasional, re-built at every reneighboring
1018 // slist = indices of perpetual NStencil classes
1019 // perpetual = used by any perpetual NPair class
1020
1021 delete [] slist;
1022 delete [] plist;
1023 nstencil_perpetual = npair_perpetual = 0;
1024 slist = new int[nstencil];
1025 plist = new int[nlist];
1026
1027 for (i = 0; i < nlist; i++) {
1028 if (lists[i]->occasional == 0 && lists[i]->pair_method)
1029 plist[npair_perpetual++] = i;
1030 }
1031
1032 for (i = 0; i < nstencil; i++) {
1033 flag = 0;
1034 for (j = 0; j < npair_perpetual; j++)
1035 if (lists[plist[j]]->stencil_method == neigh_stencil[i]->istyle)
1036 flag = 1;
1037 if (flag) slist[nstencil_perpetual++] = i;
1038 }
1039
1040 // reorder plist vector if necessary
1041 // relevant for lists that are derived from a parent list:
1042 // half-full,copy,skip
1043 // the child index must appear in plist after the parent index
1044 // swap two indices within plist when dependency is mis-ordered
1045 // start double loop check again whenever a swap is made
1046 // done when entire double loop test results in no swaps
1047
1048 NeighList *ptr;
1049
1050 int done = 0;
1051 while (!done) {
1052 done = 1;
1053 for (i = 0; i < npair_perpetual; i++) {
1054 for (k = 0; k < 3; k++) {
1055 ptr = nullptr;
1056 if (k == 0) ptr = lists[plist[i]]->listcopy;
1057 if (k == 1) ptr = lists[plist[i]]->listskip;
1058 if (k == 2) ptr = lists[plist[i]]->listfull;
1059 if (ptr == nullptr) continue;
1060 for (m = 0; m < nrequest; m++)
1061 if (ptr == lists[m]) break;
1062 for (j = 0; j < npair_perpetual; j++)
1063 if (m == plist[j]) break;
1064 if (j < i) continue;
1065 int tmp = plist[i]; // swap I,J indices
1066 plist[i] = plist[j];
1067 plist[j] = tmp;
1068 done = 0;
1069 break;
1070 }
1071 if (!done) break;
1072 }
1073 }
1074
1075 // debug output
1076
1077 #ifdef NEIGH_LIST_DEBUG
1078 for (i = 0; i < nrequest; i++) lists[i]->print_attributes();
1079 #endif
1080
1081 return same;
1082 }
1083
1084 /* ----------------------------------------------------------------------
1085 scan NeighRequests to set additional flags
1086 only for custom cutoff lists
1087 ------------------------------------------------------------------------- */
1088
morph_unique()1089 void Neighbor::morph_unique()
1090 {
1091 NeighRequest *irq;
1092
1093 for (int i = 0; i < nrequest; i++) {
1094 irq = requests[i];
1095
1096 // if cut flag set by requestor, set unique flag
1097 // this forces Pair,Stencil,Bin styles to be instantiated separately
1098
1099 if (irq->cut) irq->unique = 1;
1100 }
1101 }
1102
1103 /* ----------------------------------------------------------------------
1104 scan NeighRequests to process all skip lists
1105 look for a matching non-skip list
1106 if one exists, point at it via skiplist
1107 else make new parent via copy_request() and point at it
1108 ------------------------------------------------------------------------- */
1109
morph_skip()1110 void Neighbor::morph_skip()
1111 {
1112 int i,j,inewton,jnewton;
1113 NeighRequest *irq,*jrq,*nrq;
1114
1115 for (i = 0; i < nrequest; i++) {
1116 irq = requests[i];
1117
1118 // only processing skip lists
1119
1120 if (!irq->skip) continue;
1121
1122 // these lists are created other ways, no need for skipping
1123 // halffull list and its full parent may both skip,
1124 // but are checked to insure matching skip info
1125
1126 if (irq->halffull) continue;
1127 if (irq->copy) continue;
1128
1129 // check all other lists
1130
1131 for (j = 0; j < nrequest; j++) {
1132 if (i == j) continue;
1133 jrq = requests[j];
1134
1135 // can only skip from a perpetual non-skip list
1136
1137 if (jrq->occasional) continue;
1138 if (jrq->skip) continue;
1139
1140 // both lists must be half, or both full
1141
1142 if (irq->half != jrq->half) continue;
1143 if (irq->full != jrq->full) continue;
1144
1145 // both lists must be newton on, or both newton off
1146 // IJ newton = 1 for newton on, 2 for newton off
1147
1148 inewton = irq->newton;
1149 if (inewton == 0) inewton = force->newton_pair ? 1 : 2;
1150 jnewton = jrq->newton;
1151 if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2;
1152 if (inewton != jnewton) continue;
1153
1154 // these flags must be same,
1155 // else 2 lists do not store same pairs
1156 // or their data structures are different
1157 // this includes custom cutoff set by requestor
1158 // NOTE: need check for 2 Kokkos flags?
1159
1160 if (irq->ghost != jrq->ghost) continue;
1161 if (irq->size != jrq->size) continue;
1162 if (irq->history != jrq->history) continue;
1163 if (irq->bond != jrq->bond) continue;
1164 if (irq->omp != jrq->omp) continue;
1165 if (irq->intel != jrq->intel) continue;
1166 if (irq->kokkos_host != jrq->kokkos_host) continue;
1167 if (irq->kokkos_device != jrq->kokkos_device) continue;
1168 if (irq->ssa != jrq->ssa) continue;
1169 if (irq->cut != jrq->cut) continue;
1170 if (irq->cutoff != jrq->cutoff) continue;
1171
1172 // 2 lists are a match
1173
1174 break;
1175 }
1176
1177 // if matching list exists, point to it
1178 // else create a new identical list except non-skip
1179 // for new list, set neigh = 1, skip = 0, no skip vec/array,
1180 // copy unique flag (since copy_request() will not do it)
1181 // note: parents of skip lists do not have associated history
1182 // b/c child skip lists have the associated history
1183
1184 if (j < nrequest) irq->skiplist = j;
1185 else {
1186 int newrequest = request(this,-1);
1187 irq->skiplist = newrequest;
1188
1189 nrq = requests[newrequest];
1190 nrq->copy_request(irq,0);
1191 nrq->pair = nrq->fix = nrq->compute = nrq->command = 0;
1192 nrq->neigh = 1;
1193 nrq->skip = 0;
1194 if (irq->unique) nrq->unique = 1;
1195 }
1196 }
1197 }
1198
1199 /* ----------------------------------------------------------------------
1200 scan NeighRequests just added by morph_skip for hybrid granular
1201 adjust newton/oneside parent settings if children require onesided skipping
1202 also set children off2on flag if parent becomes a newton off list
1203 this is needed because line/gran and tri/gran pair styles
1204 require onesided neigh lists and system newton on,
1205 but parent list must be newton off to enable the onesided skipping
1206 ------------------------------------------------------------------------- */
1207
morph_granular()1208 void Neighbor::morph_granular()
1209 {
1210 int i,j;
1211 NeighRequest *irq,*jrq;
1212
1213 for (i = 0; i < nrequest; i++) {
1214 irq = requests[i];
1215
1216 // only examine NeighRequests added by morph_skip()
1217 // only those with size attribute for granular systems
1218
1219 if (!irq->neigh) continue;
1220 if (!irq->size) continue;
1221
1222 // check children of this list
1223
1224 int onesided = -1;
1225 for (j = 0; j < nrequest; j++) {
1226 jrq = requests[j];
1227
1228 // only consider JRQ pair, size lists that skip from Irq list
1229
1230 if (!jrq->pair) continue;
1231 if (!jrq->size) continue;
1232 if (!jrq->skip || jrq->skiplist != i) continue;
1233
1234 // onesided = -1 if no children
1235 // onesided = 0/1 = child granonesided value if same for all children
1236 // onesided = 2 if children have different granonesided values
1237
1238 if (onesided < 0) onesided = jrq->granonesided;
1239 else if (onesided != jrq->granonesided) onesided = 2;
1240 if (onesided == 2) break;
1241 }
1242
1243 // if onesided = 2, parent has children with both granonesided = 0/1
1244 // force parent newton off (newton = 2) to enable onesided skip by child
1245 // set parent granonesided = 0, so it stores all neighs in usual manner
1246 // set off2on = 1 for all children, since they expect newton on lists
1247 // this is b/c granonesided only set by line/gran and tri/gran which
1248 // both require system newton on
1249
1250 if (onesided == 2) {
1251 irq->newton = 2;
1252 irq->granonesided = 0;
1253
1254 for (j = 0; j < nrequest; j++) {
1255 jrq = requests[j];
1256
1257 // only consider JRQ pair, size lists that skip from Irq list
1258
1259 if (!jrq->pair) continue;
1260 if (!jrq->size) continue;
1261 if (!jrq->skip || jrq->skiplist != i) continue;
1262
1263 jrq->off2on = 1;
1264 }
1265 }
1266 }
1267 }
1268
1269 /* ----------------------------------------------------------------------
1270 scan NeighRequests for possible half lists to derive from full lists
1271 if 2 requests match, set half list to derive from full list
1272 ------------------------------------------------------------------------- */
1273
morph_halffull()1274 void Neighbor::morph_halffull()
1275 {
1276 int i,j;
1277 NeighRequest *irq,*jrq;
1278
1279 for (i = 0; i < nrequest; i++) {
1280 irq = requests[i];
1281
1282 // only processing half lists
1283
1284 if (!irq->half) continue;
1285
1286 // these lists are created other ways, no need for halffull
1287 // do want to process skip lists
1288
1289 if (irq->copy) continue;
1290
1291 // check all other lists
1292
1293 for (j = 0; j < nrequest; j++) {
1294 if (i == j) continue;
1295 jrq = requests[j];
1296
1297 // can only derive from a perpetual full list
1298 // newton setting of derived list does not matter
1299
1300 if (jrq->occasional) continue;
1301 if (!jrq->full) continue;
1302
1303 // these flags must be same,
1304 // else 2 lists do not store same pairs
1305 // or their data structures are different
1306 // this includes custom cutoff set by requestor
1307
1308 if (irq->ghost != jrq->ghost) continue;
1309 if (irq->size != jrq->size) continue;
1310 if (irq->history != jrq->history) continue;
1311 if (irq->bond != jrq->bond) continue;
1312 if (irq->omp != jrq->omp) continue;
1313 if (irq->intel != jrq->intel) continue;
1314 if (irq->kokkos_host != jrq->kokkos_host) continue;
1315 if (irq->kokkos_device != jrq->kokkos_device) continue;
1316 if (irq->ssa != jrq->ssa) continue;
1317 if (irq->cut != jrq->cut) continue;
1318 if (irq->cutoff != jrq->cutoff) continue;
1319
1320 // skip flag must be same
1321 // if both are skip lists, skip info must match
1322
1323 if (irq->skip != jrq->skip) continue;
1324 if (irq->skip && irq->same_skip(jrq) == 0) continue;
1325
1326 // 2 lists are a match
1327
1328 break;
1329 }
1330
1331 // if matching list exists, point to it
1332
1333 if (j < nrequest) {
1334 irq->halffull = 1;
1335 irq->halffulllist = j;
1336 }
1337 }
1338 }
1339
1340 /* ----------------------------------------------------------------------
1341 scan NeighRequests for possible copies
1342 if 2 requests match, turn one into a copy of the other
1343 ------------------------------------------------------------------------- */
1344
morph_copy()1345 void Neighbor::morph_copy()
1346 {
1347 int i,j,inewton,jnewton;
1348 NeighRequest *irq,*jrq;
1349
1350 for (i = 0; i < nrequest; i++) {
1351 irq = requests[i];
1352
1353 // this list is already a copy list due to another morph method
1354
1355 if (irq->copy) continue;
1356
1357 // check all other lists
1358
1359 for (j = 0; j < nrequest; j++) {
1360 if (i == j) continue;
1361 jrq = requests[j];
1362
1363 // other list is already copied from this one
1364
1365 if (jrq->copy && jrq->copylist == i) continue;
1366
1367 // other list (jrq) to copy from must be perpetual
1368 // list that becomes a copy list (irq) can be perpetual or occasional
1369 // if both lists are perpetual, require j < i
1370 // to prevent circular dependence with 3 or more copies of a list
1371
1372 if (jrq->occasional) continue;
1373 if (!irq->occasional && j > i) continue;
1374
1375 // both lists must be half, or both full
1376
1377 if (irq->half != jrq->half) continue;
1378 if (irq->full != jrq->full) continue;
1379
1380 // both lists must be newton on, or both newton off
1381 // IJ newton = 1 for newton on, 2 for newton off
1382
1383 inewton = irq->newton;
1384 if (inewton == 0) inewton = force->newton_pair ? 1 : 2;
1385 jnewton = jrq->newton;
1386 if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2;
1387 if (inewton != jnewton) continue;
1388
1389 // ok for non-ghost list to copy from ghost list, but not vice versa
1390
1391 if (irq->ghost && !jrq->ghost) continue;
1392
1393 // do not copy from a list with respa middle/inner
1394 // b/c its outer list will not be complete
1395
1396 if (jrq->respamiddle) continue;
1397 if (jrq->respainner) continue;
1398
1399 // these flags must be same,
1400 // else 2 lists do not store same pairs
1401 // or their data structures are different
1402 // this includes custom cutoff set by requestor
1403 // no need to check omp b/c it stores same pairs
1404 // NOTE: need check for 2 Kokkos flags?
1405
1406 if (irq->size != jrq->size) continue;
1407 if (irq->history != jrq->history) continue;
1408 if (irq->bond != jrq->bond) continue;
1409 if (irq->intel != jrq->intel) continue;
1410 if (irq->kokkos_host && !jrq->kokkos_host) continue;
1411 if (irq->kokkos_device && !jrq->kokkos_device) continue;
1412 if (irq->ssa != jrq->ssa) continue;
1413 if (irq->cut != jrq->cut) continue;
1414 if (irq->cutoff != jrq->cutoff) continue;
1415
1416 // skip flag must be same
1417 // if both are skip lists, skip info must match
1418
1419 if (irq->skip != jrq->skip) continue;
1420 if (irq->skip && irq->same_skip(jrq) == 0) continue;
1421
1422 // 2 lists are a match
1423
1424 break;
1425 }
1426
1427 // turn list I into a copy of list J
1428 // do not copy a list from another copy list, but from its parent list
1429
1430 if (j < nrequest) {
1431 irq->copy = 1;
1432 if (jrq->copy) irq->copylist = jrq->copylist;
1433 else irq->copylist = j;
1434 }
1435 }
1436 }
1437
1438 /* ----------------------------------------------------------------------
1439 create and initialize NTopo classes
1440 ------------------------------------------------------------------------- */
1441
init_topology()1442 void Neighbor::init_topology()
1443 {
1444 int i,m;
1445
1446 if (atom->molecular == Atom::ATOMIC) return;
1447
1448 // set flags that determine which topology neighbor classes to use
1449 // these settings could change from run to run, depending on fixes defined
1450 // bonds,etc can only be broken for atom->molecular = Atom::MOLECULAR, not Atom::TEMPLATE
1451 // SHAKE sets bonds and angles negative
1452 // gcmc sets all bonds, angles, etc negative
1453 // bond_quartic sets bonds to 0
1454 // delete_bonds sets all interactions negative
1455
1456 int bond_off = 0;
1457 int angle_off = 0;
1458 for (i = 0; i < modify->nfix; i++)
1459 if (utils::strmatch(modify->fix[i]->style,"^shake")
1460 || utils::strmatch(modify->fix[i]->style,"^rattle"))
1461 bond_off = angle_off = 1;
1462 if (force->bond && force->bond_match("quartic")) bond_off = 1;
1463
1464 if (atom->avec->bonds_allow && atom->molecular == Atom::MOLECULAR) {
1465 for (i = 0; i < atom->nlocal; i++) {
1466 if (bond_off) break;
1467 for (m = 0; m < atom->num_bond[i]; m++)
1468 if (atom->bond_type[i][m] <= 0) bond_off = 1;
1469 }
1470 }
1471
1472 if (atom->avec->angles_allow && atom->molecular == Atom::MOLECULAR) {
1473 for (i = 0; i < atom->nlocal; i++) {
1474 if (angle_off) break;
1475 for (m = 0; m < atom->num_angle[i]; m++)
1476 if (atom->angle_type[i][m] <= 0) angle_off = 1;
1477 }
1478 }
1479
1480 int dihedral_off = 0;
1481 if (atom->avec->dihedrals_allow && atom->molecular == Atom::MOLECULAR) {
1482 for (i = 0; i < atom->nlocal; i++) {
1483 if (dihedral_off) break;
1484 for (m = 0; m < atom->num_dihedral[i]; m++)
1485 if (atom->dihedral_type[i][m] <= 0) dihedral_off = 1;
1486 }
1487 }
1488
1489 int improper_off = 0;
1490 if (atom->avec->impropers_allow && atom->molecular == Atom::MOLECULAR) {
1491 for (i = 0; i < atom->nlocal; i++) {
1492 if (improper_off) break;
1493 for (m = 0; m < atom->num_improper[i]; m++)
1494 if (atom->improper_type[i][m] <= 0) improper_off = 1;
1495 }
1496 }
1497
1498 for (i = 0; i < modify->nfix; i++)
1499 if ((strcmp(modify->fix[i]->style,"gcmc") == 0))
1500 bond_off = angle_off = dihedral_off = improper_off = 1;
1501
1502 // sync on/off settings across all procs
1503
1504 int onoff = bond_off;
1505 MPI_Allreduce(&onoff,&bond_off,1,MPI_INT,MPI_MAX,world);
1506 onoff = angle_off;
1507 MPI_Allreduce(&onoff,&angle_off,1,MPI_INT,MPI_MAX,world);
1508 onoff = dihedral_off;
1509 MPI_Allreduce(&onoff,&dihedral_off,1,MPI_INT,MPI_MAX,world);
1510 onoff = improper_off;
1511 MPI_Allreduce(&onoff,&improper_off,1,MPI_INT,MPI_MAX,world);
1512
1513 // instantiate NTopo classes
1514
1515 if (atom->avec->bonds_allow) {
1516 int old_bondwhich = bondwhich;
1517 if (atom->molecular == Atom::TEMPLATE) bondwhich = TEMPLATE;
1518 else if (bond_off) bondwhich = PARTIAL;
1519 else bondwhich = ALL;
1520 if (!neigh_bond || bondwhich != old_bondwhich) {
1521 delete neigh_bond;
1522 if (bondwhich == ALL)
1523 neigh_bond = new NTopoBondAll(lmp);
1524 else if (bondwhich == PARTIAL)
1525 neigh_bond = new NTopoBondPartial(lmp);
1526 else if (bondwhich == TEMPLATE)
1527 neigh_bond = new NTopoBondTemplate(lmp);
1528 }
1529 }
1530
1531 if (atom->avec->angles_allow) {
1532 int old_anglewhich = anglewhich;
1533 if (atom->molecular == Atom::TEMPLATE) anglewhich = TEMPLATE;
1534 else if (angle_off) anglewhich = PARTIAL;
1535 else anglewhich = ALL;
1536 if (!neigh_angle || anglewhich != old_anglewhich) {
1537 delete neigh_angle;
1538 if (anglewhich == ALL)
1539 neigh_angle = new NTopoAngleAll(lmp);
1540 else if (anglewhich == PARTIAL)
1541 neigh_angle = new NTopoAnglePartial(lmp);
1542 else if (anglewhich == TEMPLATE)
1543 neigh_angle = new NTopoAngleTemplate(lmp);
1544 }
1545 }
1546
1547 if (atom->avec->dihedrals_allow) {
1548 int old_dihedralwhich = dihedralwhich;
1549 if (atom->molecular == Atom::TEMPLATE) dihedralwhich = TEMPLATE;
1550 else if (dihedral_off) dihedralwhich = PARTIAL;
1551 else dihedralwhich = ALL;
1552 if (!neigh_dihedral || dihedralwhich != old_dihedralwhich) {
1553 delete neigh_dihedral;
1554 if (dihedralwhich == ALL)
1555 neigh_dihedral = new NTopoDihedralAll(lmp);
1556 else if (dihedralwhich == PARTIAL)
1557 neigh_dihedral = new NTopoDihedralPartial(lmp);
1558 else if (dihedralwhich == TEMPLATE)
1559 neigh_dihedral = new NTopoDihedralTemplate(lmp);
1560 }
1561 }
1562
1563 if (atom->avec->impropers_allow) {
1564 int old_improperwhich = improperwhich;
1565 if (atom->molecular == Atom::TEMPLATE) improperwhich = TEMPLATE;
1566 else if (improper_off) improperwhich = PARTIAL;
1567 else improperwhich = ALL;
1568 if (!neigh_improper || improperwhich != old_improperwhich) {
1569 delete neigh_improper;
1570 if (improperwhich == ALL)
1571 neigh_improper = new NTopoImproperAll(lmp);
1572 else if (improperwhich == PARTIAL)
1573 neigh_improper = new NTopoImproperPartial(lmp);
1574 else if (improperwhich == TEMPLATE)
1575 neigh_improper = new NTopoImproperTemplate(lmp);
1576 }
1577 }
1578 }
1579
1580 /* ----------------------------------------------------------------------
1581 output summary of pairwise neighbor list info
1582 only called by proc 0
1583 ------------------------------------------------------------------------- */
1584
print_pairwise_info()1585 void Neighbor::print_pairwise_info()
1586 {
1587 int i;
1588 NeighRequest *rq;
1589
1590 const double cutghost = MAX(cutneighmax,comm->cutghostuser);
1591
1592 double binsize, bbox[3];
1593 bbox[0] = bboxhi[0]-bboxlo[0];
1594 bbox[1] = bboxhi[1]-bboxlo[1];
1595 bbox[2] = bboxhi[2]-bboxlo[2];
1596 if (binsizeflag) binsize = binsize_user;
1597 else if (style == Neighbor::BIN) binsize = 0.5*cutneighmax;
1598 else binsize = 0.5*cutneighmin;
1599 if (binsize == 0.0) binsize = bbox[0];
1600
1601 int nperpetual = 0;
1602 int noccasional = 0;
1603 int nextra = 0;
1604 for (i = 0; i < nlist; i++) {
1605 if (lists[i]->pair_method == 0) nextra++;
1606 else if (lists[i]->occasional) noccasional++;
1607 else nperpetual++;
1608 }
1609
1610 std::string out = "Neighbor list info ...\n";
1611 out += fmt::format(" update every {} steps, delay {} steps, check {}\n",
1612 every,delay,dist_check ? "yes" : "no");
1613 out += fmt::format(" max neighbors/atom: {}, page size: {}\n",
1614 oneatom, pgsize);
1615 out += fmt::format(" master list distance cutoff = {:.8g}\n",cutneighmax);
1616 out += fmt::format(" ghost atom cutoff = {:.8g}\n",cutghost);
1617 if (style != Neighbor::NSQ)
1618 out += fmt::format(" binsize = {:.8g}, bins = {:g} {:g} {:g}\n",binsize,
1619 ceil(bbox[0]/binsize), ceil(bbox[1]/binsize),
1620 ceil(bbox[2]/binsize));
1621
1622 out += fmt::format(" {} neighbor lists, perpetual/occasional/extra = {} {} {}\n",
1623 nlist,nperpetual,noccasional,nextra);
1624
1625 for (i = 0; i < nlist; i++) {
1626 rq = requests[i];
1627 if (rq->pair) {
1628 char *pname = force->pair_match_ptr((Pair *) rq->requestor);
1629 if (pname) out += fmt::format(" ({}) pair {}",i+1,pname);
1630 else out += fmt::format(" ({}) pair (none)",i+1);
1631 } else if (rq->fix) {
1632 out += fmt::format(" ({}) fix {}",i+1,((Fix *) rq->requestor)->style);
1633 } else if (rq->compute) {
1634 out += fmt::format(" ({}) compute {}",i+1,((Compute *) rq->requestor)->style);
1635 } else if (rq->command) {
1636 out += fmt::format(" ({}) command {}",i+1,rq->command_style);
1637 } else if (rq->neigh) {
1638 out += fmt::format(" ({}) neighbor class addition",i+1);
1639 }
1640
1641 if (rq->occasional) out += ", occasional";
1642 else out += ", perpetual";
1643
1644 // order these to get single output of most relevant
1645
1646 if (rq->copy)
1647 out += fmt::format(", copy from ({})",rq->copylist+1);
1648 else if (rq->halffull)
1649 out += fmt::format(", half/full from ({})",rq->halffulllist+1);
1650 else if (rq->skip)
1651 out += fmt::format(", skip from ({})",rq->skiplist+1);
1652 out += "\n";
1653
1654 // list of neigh list attributes
1655
1656 out += " attributes: ";
1657 if (rq->half) out += "half";
1658 else if (rq->full) out += "full";
1659
1660 if (rq->newton == 0) {
1661 if (force->newton_pair) out += ", newton on";
1662 else out += ", newton off";
1663 } else if (rq->newton == 1) out += ", newton on";
1664 else if (rq->newton == 2) out += ", newton off";
1665
1666 if (rq->ghost) out += ", ghost";
1667 if (rq->size) out += ", size";
1668 if (rq->history) out += ", history";
1669 if (rq->granonesided) out += ", onesided";
1670 if (rq->respamiddle) out += ", respa outer/middle/inner";
1671 else if (rq->respainner) out += ", respa outer/inner";
1672 if (rq->bond) out += ", bond";
1673 if (rq->omp) out += ", omp";
1674 if (rq->intel) out += ", intel";
1675 if (rq->kokkos_device) out += ", kokkos_device";
1676 if (rq->kokkos_host) out += ", kokkos_host";
1677 if (rq->ssa) out += ", ssa";
1678 if (rq->cut) out += fmt::format(", cut {}",rq->cutoff);
1679 if (rq->off2on) out += ", off2on";
1680 out += "\n";
1681
1682 out += " ";
1683 if (lists[i]->pair_method == 0) out += "pair build: none\n";
1684 else out += fmt::format("pair build: {}\n",pairnames[lists[i]->pair_method-1]);
1685
1686 out += " ";
1687 if (lists[i]->stencil_method == 0) out += "stencil: none\n";
1688 else out += fmt::format("stencil: {}\n",stencilnames[lists[i]->stencil_method-1]);
1689
1690 out += " ";
1691 if (lists[i]->bin_method == 0) out += "bin: none\n";
1692 else out += fmt::format("bin: {}\n",binnames[lists[i]->bin_method-1]);
1693
1694 }
1695 utils::logmesg(lmp,out);
1696 }
1697
1698 /* ----------------------------------------------------------------------
1699 make copy of current requests and Neighbor params
1700 used to compare to when next run occurs
1701 ------------------------------------------------------------------------- */
1702
requests_new2old()1703 void Neighbor::requests_new2old()
1704 {
1705 for (int i = 0; i < old_nrequest; i++) delete old_requests[i];
1706 memory->sfree(old_requests);
1707
1708 old_nrequest = nrequest;
1709 old_requests = (NeighRequest **)
1710 memory->smalloc(old_nrequest*sizeof(NeighRequest *),
1711 "neighbor:old_requests");
1712
1713 for (int i = 0; i < old_nrequest; i++) {
1714 old_requests[i] = new NeighRequest(lmp);
1715 old_requests[i]->copy_request(requests[i],1);
1716 }
1717
1718 old_style = style;
1719 old_triclinic = triclinic;
1720 old_pgsize = pgsize;
1721 old_oneatom = oneatom;
1722 }
1723
1724 /* ----------------------------------------------------------------------
1725 find and return request made by classptr
1726 if not found or classpt = nullptr, return nullptr
1727 ------------------------------------------------------------------------- */
1728
find_request(void * classptr)1729 NeighRequest *Neighbor::find_request(void *classptr)
1730 {
1731 if (classptr == nullptr) return nullptr;
1732
1733 for (int i = 0; i < nrequest; i++)
1734 if (requests[i]->requestor == classptr) return requests[i];
1735
1736 return nullptr;
1737 }
1738
1739 /* ----------------------------------------------------------------------
1740 assign NBin class to a NeighList
1741 use neigh request settings to build mask
1742 match mask to list of masks of known Nbin classes
1743 return index+1 of match in list of masks
1744 return 0 for no binning
1745 return -1 if no match
1746 ------------------------------------------------------------------------- */
1747
choose_bin(NeighRequest * rq)1748 int Neighbor::choose_bin(NeighRequest *rq)
1749 {
1750 // no binning needed
1751
1752 if (style == Neighbor::NSQ) return 0;
1753 if (rq->skip || rq->copy || rq->halffull) return 0;
1754
1755 // use request settings to match exactly one NBin class mask
1756 // checks are bitwise using NeighConst bit masks
1757
1758 int mask;
1759
1760 for (int i = 0; i < nbclass; i++) {
1761 mask = binmasks[i];
1762
1763 // require match of these request flags and mask bits
1764 // (!A != !B) is effectively a logical xor
1765
1766 if (!rq->intel != !(mask & NB_INTEL)) continue;
1767 if (!rq->ssa != !(mask & NB_SSA)) continue;
1768 if (!rq->kokkos_device != !(mask & NB_KOKKOS_DEVICE)) continue;
1769 if (!rq->kokkos_host != !(mask & NB_KOKKOS_HOST)) continue;
1770
1771 // multi neighbor style require multi bin style
1772 if (style == Neighbor::MULTI) {
1773 if (!(mask & NB_MULTI)) continue;
1774 } else {
1775 if (!(mask & NB_STANDARD)) continue;
1776 }
1777
1778 return i+1;
1779 }
1780
1781 // error return if matched none
1782
1783 return -1;
1784 }
1785
1786 /* ----------------------------------------------------------------------
1787 assign NStencil class to a NeighList
1788 use neigh request settings to build mask
1789 match mask to list of masks of known NStencil classes
1790 return index+1 of match in list of masks
1791 return 0 for no binning
1792 return -1 if no match
1793 ------------------------------------------------------------------------- */
1794
choose_stencil(NeighRequest * rq)1795 int Neighbor::choose_stencil(NeighRequest *rq)
1796 {
1797 // no stencil creation needed
1798
1799 if (style == Neighbor::NSQ) return 0;
1800 if (rq->skip || rq->copy || rq->halffull) return 0;
1801
1802 // convert newton request to newtflag = on or off
1803
1804 int newtflag = 1;
1805 if (rq->newton == 0 && newton_pair) newtflag = 1;
1806 else if (rq->newton == 0 && !newton_pair) newtflag = 0;
1807 else if (rq->newton == 1) newtflag = 1;
1808 else if (rq->newton == 2) newtflag = 0;
1809
1810 // request a full stencil if building full neighbor list or newton is off
1811 int fullflag = 0;
1812 if (rq->full) fullflag = 1;
1813 if (!newtflag) fullflag = 1;
1814
1815 //printf("STENCIL RQ FLAGS: hff %d %d n %d g %d s %d newtflag %d fullflag %d\n",
1816 // rq->half,rq->full,rq->newton,rq->ghost,rq->ssa,
1817 // newtflag, fullflag);
1818
1819 // use request and system settings to match exactly one NStencil class mask
1820 // checks are bitwise using NeighConst bit masks
1821
1822 int mask;
1823
1824 for (int i = 0; i < nsclass; i++) {
1825 mask = stencilmasks[i];
1826
1827 //printf("III %d: half %d full %d ghost %d ssa %d\n",
1828 // i,mask & NS_HALF,mask & NS_FULL,mask & NS_GHOST,mask & NS_SSA);
1829
1830 // exactly one of half or full is set and must match
1831
1832 if (fullflag) {
1833 if (!(mask & NS_FULL)) continue;
1834 } else {
1835 if (!(mask & NS_HALF)) continue;
1836 }
1837
1838 // require match of these request flags and mask bits
1839 // (!A != !B) is effectively a logical xor
1840
1841 if (!rq->ghost != !(mask & NS_GHOST)) continue;
1842 if (!rq->ssa != !(mask & NS_SSA)) continue;
1843
1844 // neighbor style is one of BIN, MULTI_OLD, or MULTI and must match
1845
1846 if (style == Neighbor::BIN) {
1847 if (!(mask & NS_BIN)) continue;
1848 } else if (style == Neighbor::MULTI_OLD) {
1849 if (!(mask & NS_MULTI_OLD)) continue;
1850 } else if (style == Neighbor::MULTI) {
1851 if (!(mask & NS_MULTI)) continue;
1852 }
1853
1854 // dimension is 2 or 3 and must match
1855
1856 if (dimension == 2) {
1857 if (!(mask & NS_2D)) continue;
1858 } else if (dimension == 3) {
1859 if (!(mask & NS_3D)) continue;
1860 }
1861
1862 // domain triclinic flag is on or off and must match
1863
1864 if (triclinic) {
1865 if (!(mask & NS_TRI)) continue;
1866 } else if (!triclinic) {
1867 if (!(mask & NS_ORTHO)) continue;
1868 }
1869
1870 return i+1;
1871 }
1872
1873 // error return if matched none
1874
1875 return -1;
1876 }
1877
1878 /* ----------------------------------------------------------------------
1879 assign NPair class to a NeighList
1880 use neigh request settings to build mask
1881 match mask to list of masks of known NPair classes
1882 return index+1 of match in list of masks
1883 return 0 for no binning
1884 return -1 if no match
1885 ------------------------------------------------------------------------- */
1886
choose_pair(NeighRequest * rq)1887 int Neighbor::choose_pair(NeighRequest *rq)
1888 {
1889 // error check for includegroup with ghost neighbor request
1890
1891 if (includegroup && rq->ghost)
1892 error->all(FLERR,"Neighbor include group not allowed with ghost neighbors");
1893
1894 // convert newton request to newtflag = on or off
1895
1896 bool newtflag;
1897 if (rq->newton == 0 && newton_pair) newtflag = true;
1898 else if (rq->newton == 0 && !newton_pair) newtflag = false;
1899 else if (rq->newton == 1) newtflag = true;
1900 else if (rq->newton == 2) newtflag = false;
1901 else error->all(FLERR,"Illegal 'newton' flag in neighbor list request");
1902
1903 int molecular = atom->molecular;
1904
1905 //printf("PAIR RQ FLAGS: hf %d %d n %d g %d sz %d gos %d r %d b %d o %d i %d "
1906 // "kk %d %d ss %d dn %d sk %d cp %d hf %d oo %d\n",
1907 // rq->half,rq->full,rq->newton,rq->ghost,rq->size,
1908 // rq->granonesided,rq->respaouter,rq->bond,rq->omp,rq->intel,
1909 // rq->kokkos_host,rq->kokkos_device,rq->ssa,rq->dnum,
1910 // rq->skip,rq->copy,rq->halffull,rq->off2on);
1911
1912 // use request and system settings to match exactly one NPair class mask
1913 // checks are bitwise using NeighConst bit masks
1914
1915 int mask;
1916
1917 for (int i = 0; i < npclass; i++) {
1918 mask = pairmasks[i];
1919
1920 //printf(" PAIR NAMES i %d %d name %s mask %d\n",i,nrequest,
1921 // pairnames[i],pairmasks[i]);
1922
1923 // if copy request, no further checks needed, just return or continue
1924 // Kokkos device/host flags must also match in order to copy
1925
1926 if (rq->copy) {
1927 if (!(mask & NP_COPY)) continue;
1928 if (rq->kokkos_device || rq->kokkos_host) {
1929 if (!rq->kokkos_device != !(mask & NP_KOKKOS_DEVICE)) continue;
1930 if (!rq->kokkos_host != !(mask & NP_KOKKOS_HOST)) continue;
1931 }
1932 if (!requests[rq->copylist]->kokkos_device != !(mask & NP_KOKKOS_DEVICE)) continue;
1933 if (!requests[rq->copylist]->kokkos_host != !(mask & NP_KOKKOS_HOST)) continue;
1934 return i+1;
1935 }
1936
1937 // exactly one of half or full is set and must match
1938
1939 if (rq->half) {
1940 if (!(mask & NP_HALF)) continue;
1941 } else if (rq->full) {
1942 if (!(mask & NP_FULL)) continue;
1943 }
1944
1945 // newtflag is on or off and must match
1946
1947 if (newtflag) {
1948 if (!(mask & NP_NEWTON)) continue;
1949 } else if (!newtflag) {
1950 if (!(mask & NP_NEWTOFF)) continue;
1951 }
1952
1953 // if molecular on, do not match ATOMONLY (b/c a MOLONLY Npair exists)
1954 // if molecular off, do not match MOLONLY (b/c an ATOMONLY Npair exists)
1955
1956 if (molecular != Atom::ATOMIC) {
1957 if (mask & NP_ATOMONLY) continue;
1958 } else if (molecular == Atom::ATOMIC) {
1959 if (mask & NP_MOLONLY) continue;
1960 }
1961
1962 // require match of these request flags and mask bits
1963 // (!A != !B) is effectively a logical xor
1964
1965 if (!rq->ghost != !(mask & NP_GHOST)) continue;
1966 if (!rq->size != !(mask & NP_SIZE)) continue;
1967 if (!rq->respaouter != !(mask & NP_RESPA)) continue;
1968 if (!rq->granonesided != !(mask & NP_ONESIDE)) continue;
1969 if (!rq->respaouter != !(mask & NP_RESPA)) continue;
1970 if (!rq->bond != !(mask & NP_BOND)) continue;
1971 if (!rq->omp != !(mask & NP_OMP)) continue;
1972 if (!rq->intel != !(mask & NP_INTEL)) continue;
1973 if (!rq->kokkos_device != !(mask & NP_KOKKOS_DEVICE)) continue;
1974 if (!rq->kokkos_host != !(mask & NP_KOKKOS_HOST)) continue;
1975 if (!rq->ssa != !(mask & NP_SSA)) continue;
1976
1977 if (!rq->skip != !(mask & NP_SKIP)) continue;
1978
1979 if (!rq->halffull != !(mask & NP_HALF_FULL)) continue;
1980 if (!rq->off2on != !(mask & NP_OFF2ON)) continue;
1981
1982 // neighbor style is one of NSQ, BIN, MULTI_OLD, or MULTI and must match
1983
1984 if (style == Neighbor::NSQ) {
1985 if (!(mask & NP_NSQ)) continue;
1986 } else if (style == Neighbor::BIN) {
1987 if (!(mask & NP_BIN)) continue;
1988 } else if (style == Neighbor::MULTI_OLD) {
1989 if (!(mask & NP_MULTI_OLD)) continue;
1990 } else if (style == Neighbor::MULTI) {
1991 if (!(mask & NP_MULTI)) continue;
1992 }
1993
1994 // domain triclinic flag is on or off and must match
1995
1996 if (triclinic) {
1997 if (!(mask & NP_TRI)) continue;
1998 } else if (!triclinic) {
1999 if (!(mask & NP_ORTHO)) continue;
2000 }
2001
2002 return i+1;
2003 }
2004
2005 // error return if matched none
2006
2007 return -1;
2008 }
2009
2010 /* ----------------------------------------------------------------------
2011 called by other classes to request a pairwise neighbor list
2012 ------------------------------------------------------------------------- */
2013
request(void * requestor,int instance)2014 int Neighbor::request(void *requestor, int instance)
2015 {
2016 if (nrequest == maxrequest) {
2017 maxrequest += RQDELTA;
2018 requests = (NeighRequest **)
2019 memory->srealloc(requests,maxrequest*sizeof(NeighRequest *),
2020 "neighbor:requests");
2021 }
2022
2023 requests[nrequest] = new NeighRequest(lmp);
2024 requests[nrequest]->index = nrequest;
2025 requests[nrequest]->requestor = requestor;
2026 requests[nrequest]->requestor_instance = instance;
2027 nrequest++;
2028 return nrequest-1;
2029 }
2030
2031 /* ----------------------------------------------------------------------
2032 one instance per entry in style_neigh_bin.h
2033 ------------------------------------------------------------------------- */
2034
2035 template <typename T>
bin_creator(LAMMPS * lmp)2036 NBin *Neighbor::bin_creator(LAMMPS *lmp)
2037 {
2038 return new T(lmp);
2039 }
2040
2041 /* ----------------------------------------------------------------------
2042 one instance per entry in style_neigh_stencil.h
2043 ------------------------------------------------------------------------- */
2044
2045 template <typename T>
stencil_creator(LAMMPS * lmp)2046 NStencil *Neighbor::stencil_creator(LAMMPS *lmp)
2047 {
2048 return new T(lmp);
2049 }
2050
2051 /* ----------------------------------------------------------------------
2052 one instance per entry in style_neigh_pair.h
2053 ------------------------------------------------------------------------- */
2054
2055 template <typename T>
pair_creator(LAMMPS * lmp)2056 NPair *Neighbor::pair_creator(LAMMPS *lmp)
2057 {
2058 return new T(lmp);
2059 }
2060
2061 /* ----------------------------------------------------------------------
2062 setup neighbor binning and neighbor stencils
2063 called before run and every reneighbor if box size/shape changes
2064 only operates on perpetual lists
2065 build_one() operates on occasional lists
2066 ------------------------------------------------------------------------- */
2067
setup_bins()2068 void Neighbor::setup_bins()
2069 {
2070 // invoke setup_bins() for all NBin
2071 // actual binning is performed in build()
2072
2073 for (int i = 0; i < nbin; i++)
2074 neigh_bin[i]->setup_bins(style);
2075
2076 // invoke create_setup() and create() for all perpetual NStencil
2077 // same ops performed for occasional lists in build_one()
2078
2079 for (int i = 0; i < nstencil_perpetual; i++) {
2080 neigh_stencil[slist[i]]->create_setup();
2081 neigh_stencil[slist[i]]->create();
2082 }
2083
2084 last_setup_bins = update->ntimestep;
2085 }
2086
2087 /* ---------------------------------------------------------------------- */
2088
decide()2089 int Neighbor::decide()
2090 {
2091 if (must_check) {
2092 bigint n = update->ntimestep;
2093 if (restart_check && n == output->next_restart) return 1;
2094 for (int i = 0; i < fix_check; i++)
2095 if (n == modify->fix[fixchecklist[i]]->next_reneighbor) return 1;
2096 }
2097
2098 ago++;
2099 if (ago >= delay && ago % every == 0) {
2100 if (build_once) return 0;
2101 if (dist_check == 0) return 1;
2102 return check_distance();
2103 } else return 0;
2104 }
2105
2106 /* ----------------------------------------------------------------------
2107 if any atom moved trigger distance (half of neighbor skin) return 1
2108 shrink trigger distance if box size has changed
2109 conservative shrink procedure:
2110 compute distance each of 8 corners of box has moved since last reneighbor
2111 reduce skin distance by sum of 2 largest of the 8 values
2112 if reduced skin distance is negative, set to zero
2113 new trigger = 1/2 of reduced skin distance
2114 for orthogonal box, only need 2 lo/hi corners
2115 for triclinic, need all 8 corners since deformations can displace all 8
2116 ------------------------------------------------------------------------- */
2117
check_distance()2118 int Neighbor::check_distance()
2119 {
2120 double delx,dely,delz,rsq;
2121 double delta,deltasq,delta1,delta2;
2122
2123 if (boxcheck) {
2124 if (triclinic == 0) {
2125 delx = bboxlo[0] - boxlo_hold[0];
2126 dely = bboxlo[1] - boxlo_hold[1];
2127 delz = bboxlo[2] - boxlo_hold[2];
2128 delta1 = sqrt(delx*delx + dely*dely + delz*delz);
2129 delx = bboxhi[0] - boxhi_hold[0];
2130 dely = bboxhi[1] - boxhi_hold[1];
2131 delz = bboxhi[2] - boxhi_hold[2];
2132 delta2 = sqrt(delx*delx + dely*dely + delz*delz);
2133 delta = 0.5 * (skin - (delta1+delta2));
2134 if (delta < 0.0) delta = 0.0;
2135 deltasq = delta*delta;
2136 } else {
2137 domain->box_corners();
2138 delta1 = delta2 = 0.0;
2139 for (int i = 0; i < 8; i++) {
2140 delx = corners[i][0] - corners_hold[i][0];
2141 dely = corners[i][1] - corners_hold[i][1];
2142 delz = corners[i][2] - corners_hold[i][2];
2143 delta = sqrt(delx*delx + dely*dely + delz*delz);
2144 if (delta > delta1) delta1 = delta;
2145 else if (delta > delta2) delta2 = delta;
2146 }
2147 delta = 0.5 * (skin - (delta1+delta2));
2148 if (delta < 0.0) delta = 0.0;
2149 deltasq = delta*delta;
2150 }
2151 } else deltasq = triggersq;
2152
2153 double **x = atom->x;
2154 int nlocal = atom->nlocal;
2155 if (includegroup) nlocal = atom->nfirst;
2156
2157 int flag = 0;
2158 for (int i = 0; i < nlocal; i++) {
2159 delx = x[i][0] - xhold[i][0];
2160 dely = x[i][1] - xhold[i][1];
2161 delz = x[i][2] - xhold[i][2];
2162 rsq = delx*delx + dely*dely + delz*delz;
2163 if (rsq > deltasq) flag = 1;
2164 }
2165
2166 int flagall;
2167 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
2168 if (flagall && ago == MAX(every,delay)) ndanger++;
2169 return flagall;
2170 }
2171
2172 /* ----------------------------------------------------------------------
2173 build perpetual neighbor lists
2174 called at setup and every few timesteps during run or minimization
2175 topology lists also built if topoflag = 1 (Kokkos calls with topoflag=0)
2176 ------------------------------------------------------------------------- */
2177
build(int topoflag)2178 void Neighbor::build(int topoflag)
2179 {
2180 int i,m;
2181
2182 ago = 0;
2183 ncalls++;
2184 lastcall = update->ntimestep;
2185
2186 int nlocal = atom->nlocal;
2187 int nall = nlocal + atom->nghost;
2188 // rebuild collection array from scratch
2189 if (style == Neighbor::MULTI) build_collection(0);
2190
2191 // check that using special bond flags will not overflow neigh lists
2192
2193 if (nall > NEIGHMASK)
2194 error->one(FLERR,"Too many local+ghost atoms for neighbor list");
2195
2196 // store current atom positions and box size if needed
2197
2198 if (dist_check) {
2199 double **x = atom->x;
2200 if (includegroup) nlocal = atom->nfirst;
2201 if (atom->nmax > maxhold) {
2202 maxhold = atom->nmax;
2203 memory->destroy(xhold);
2204 memory->create(xhold,maxhold,3,"neigh:xhold");
2205 }
2206 for (i = 0; i < nlocal; i++) {
2207 xhold[i][0] = x[i][0];
2208 xhold[i][1] = x[i][1];
2209 xhold[i][2] = x[i][2];
2210 }
2211 if (boxcheck) {
2212 if (triclinic == 0) {
2213 boxlo_hold[0] = bboxlo[0];
2214 boxlo_hold[1] = bboxlo[1];
2215 boxlo_hold[2] = bboxlo[2];
2216 boxhi_hold[0] = bboxhi[0];
2217 boxhi_hold[1] = bboxhi[1];
2218 boxhi_hold[2] = bboxhi[2];
2219 } else {
2220 domain->box_corners();
2221 corners = domain->corners;
2222 for (i = 0; i < 8; i++) {
2223 corners_hold[i][0] = corners[i][0];
2224 corners_hold[i][1] = corners[i][1];
2225 corners_hold[i][2] = corners[i][2];
2226 }
2227 }
2228 }
2229 }
2230
2231 // bin atoms for all NBin instances
2232 // not just NBin associated with perpetual lists, also occasional lists
2233 // b/c cannot wait to bin occasional lists in build_one() call
2234 // if bin then, atoms may have moved outside of proc domain & bin extent,
2235 // leading to errors or even a crash
2236
2237 if (style != Neighbor::NSQ) {
2238 if (last_setup_bins < 0) setup_bins();
2239 for (i = 0; i < nbin; i++) {
2240 neigh_bin[i]->bin_atoms_setup(nall);
2241 neigh_bin[i]->bin_atoms();
2242 }
2243 }
2244
2245 // build pairwise lists for all perpetual NPair/NeighList
2246 // grow() with nlocal/nall args so that only realloc if have to
2247
2248 for (i = 0; i < npair_perpetual; i++) {
2249 m = plist[i];
2250 if (!lists[m]->copy || lists[m]->kk2cpu)
2251 lists[m]->grow(nlocal,nall);
2252 neigh_pair[m]->build_setup();
2253 neigh_pair[m]->build(lists[m]);
2254 }
2255
2256 // build topology lists for bonds/angles/etc
2257
2258 if ((atom->molecular != Atom::ATOMIC) && topoflag) build_topology();
2259 }
2260
2261 /* ----------------------------------------------------------------------
2262 build topology neighbor lists: bond, angle, dihedral, improper
2263 copy their list info back to Neighbor for access by bond/angle/etc classes
2264 ------------------------------------------------------------------------- */
2265
build_topology()2266 void Neighbor::build_topology()
2267 {
2268 if (force->bond) {
2269 neigh_bond->build();
2270 nbondlist = neigh_bond->nbondlist;
2271 bondlist = neigh_bond->bondlist;
2272 }
2273 if (force->angle) {
2274 neigh_angle->build();
2275 nanglelist = neigh_angle->nanglelist;
2276 anglelist = neigh_angle->anglelist;
2277 }
2278 if (force->dihedral) {
2279 neigh_dihedral->build();
2280 ndihedrallist = neigh_dihedral->ndihedrallist;
2281 dihedrallist = neigh_dihedral->dihedrallist;
2282 }
2283 if (force->improper) {
2284 neigh_improper->build();
2285 nimproperlist = neigh_improper->nimproperlist;
2286 improperlist = neigh_improper->improperlist;
2287 }
2288 }
2289
2290 /* ----------------------------------------------------------------------
2291 build a single occasional pairwise neighbor list indexed by I
2292 called by other classes
2293 ------------------------------------------------------------------------- */
2294
build_one(class NeighList * mylist,int preflag)2295 void Neighbor::build_one(class NeighList *mylist, int preflag)
2296 {
2297 // check if list structure is initialized
2298
2299 if (mylist == nullptr)
2300 error->all(FLERR,"Trying to build an occasional neighbor list "
2301 "before initialization completed");
2302
2303 // build_one() should never be invoked on a perpetual list
2304
2305 if (!mylist->occasional)
2306 error->all(FLERR,"Neighbor build one invoked on perpetual list");
2307
2308 // no need to build if already built since last re-neighbor
2309 // preflag is set by fix bond/create and fix bond/swap
2310 // b/c they invoke build_one() on same step neigh list is re-built,
2311 // but before re-build, so need to use ">" instead of ">="
2312
2313 NPair *np = neigh_pair[mylist->index];
2314
2315 if (preflag) {
2316 if (np->last_build > lastcall) return;
2317 } else {
2318 if (np->last_build >= lastcall) return;
2319 }
2320
2321 // if this is copy list and parent is occasional list,
2322 // or this is halffull and parent is occasional list,
2323 // or this is skip list and parent is occasional list,
2324 // insure parent is current
2325
2326 if (mylist->listcopy && mylist->listcopy->occasional)
2327 build_one(mylist->listcopy,preflag);
2328 if (mylist->listfull && mylist->listfull->occasional)
2329 build_one(mylist->listfull,preflag);
2330 if (mylist->listskip && mylist->listskip->occasional)
2331 build_one(mylist->listskip,preflag);
2332
2333 // create stencil if hasn't been created since last setup_bins() call
2334
2335 NStencil *ns = np->ns;
2336 if (ns && ns->last_stencil < last_setup_bins) {
2337 ns->create_setup();
2338 ns->create();
2339 }
2340
2341 // build the list
2342
2343 if (!mylist->copy || mylist->kk2cpu)
2344 mylist->grow(atom->nlocal,atom->nlocal+atom->nghost);
2345 np->build_setup();
2346 np->build(mylist);
2347 }
2348
2349 /* ----------------------------------------------------------------------
2350 set neighbor style and skin distance
2351 ------------------------------------------------------------------------- */
2352
set(int narg,char ** arg)2353 void Neighbor::set(int narg, char **arg)
2354 {
2355 if (narg != 2) error->all(FLERR,"Illegal neighbor command");
2356
2357 skin = utils::numeric(FLERR,arg[0],false,lmp);
2358 if (skin < 0.0) error->all(FLERR,"Illegal neighbor command");
2359
2360 if (strcmp(arg[1],"nsq") == 0) style = Neighbor::NSQ;
2361 else if (strcmp(arg[1],"bin") == 0) style = Neighbor::BIN;
2362 else if (strcmp(arg[1],"multi") == 0) {
2363 style = Neighbor::MULTI;
2364 ncollections = atom->ntypes;
2365 } else if (strcmp(arg[1],"multi/old") == 0) style = Neighbor::MULTI_OLD;
2366 else error->all(FLERR,"Illegal neighbor command");
2367
2368 if (style == Neighbor::MULTI_OLD && lmp->citeme) lmp->citeme->add(cite_neigh_multi_old);
2369 if (style == Neighbor::MULTI && lmp->citeme) lmp->citeme->add(cite_neigh_multi);
2370 }
2371
2372 /* ----------------------------------------------------------------------
2373 reset timestamps in all NeignBin, NStencil, NPair classes
2374 so that neighbor lists will rebuild properly with timestep change
2375 ditto for lastcall and last_setup_bins
2376 ------------------------------------------------------------------------- */
2377
reset_timestep(bigint)2378 void Neighbor::reset_timestep(bigint /*ntimestep*/)
2379 {
2380 for (int i = 0; i < nbin; i++)
2381 neigh_bin[i]->last_bin = -1;
2382 for (int i = 0; i < nstencil; i++)
2383 neigh_stencil[i]->last_stencil = -1;
2384 for (int i = 0; i < nlist; i++) {
2385 if (!neigh_pair[i]) continue;
2386 neigh_pair[i]->last_build = -1;
2387 }
2388
2389 lastcall = -1;
2390 last_setup_bins = -1;
2391 }
2392
2393 /* ----------------------------------------------------------------------
2394 modify parameters of the pair-wise neighbor build
2395 ------------------------------------------------------------------------- */
2396
modify_params(int narg,char ** arg)2397 void Neighbor::modify_params(int narg, char **arg)
2398 {
2399 int iarg = 0;
2400 while (iarg < narg) {
2401 if (strcmp(arg[iarg],"every") == 0) {
2402 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2403 every = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
2404 if (every <= 0) error->all(FLERR,"Illegal neigh_modify command");
2405 iarg += 2;
2406 } else if (strcmp(arg[iarg],"delay") == 0) {
2407 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2408 delay = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
2409 if (delay < 0) error->all(FLERR,"Illegal neigh_modify command");
2410 iarg += 2;
2411 } else if (strcmp(arg[iarg],"check") == 0) {
2412 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2413 if (strcmp(arg[iarg+1],"yes") == 0) dist_check = 1;
2414 else if (strcmp(arg[iarg+1],"no") == 0) dist_check = 0;
2415 else error->all(FLERR,"Illegal neigh_modify command");
2416 iarg += 2;
2417 } else if (strcmp(arg[iarg],"once") == 0) {
2418 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2419 if (strcmp(arg[iarg+1],"yes") == 0) build_once = 1;
2420 else if (strcmp(arg[iarg+1],"no") == 0) build_once = 0;
2421 else error->all(FLERR,"Illegal neigh_modify command");
2422 iarg += 2;
2423 } else if (strcmp(arg[iarg],"page") == 0) {
2424 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2425 old_pgsize = pgsize;
2426 pgsize = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
2427 iarg += 2;
2428 } else if (strcmp(arg[iarg],"one") == 0) {
2429 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2430 old_oneatom = oneatom;
2431 oneatom = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
2432 iarg += 2;
2433 } else if (strcmp(arg[iarg],"binsize") == 0) {
2434 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2435 binsize_user = utils::numeric(FLERR,arg[iarg+1],false,lmp);
2436 if (binsize_user <= 0.0) binsizeflag = 0;
2437 else binsizeflag = 1;
2438 iarg += 2;
2439 } else if (strcmp(arg[iarg],"cluster") == 0) {
2440 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2441 if (strcmp(arg[iarg+1],"yes") == 0) cluster_check = 1;
2442 else if (strcmp(arg[iarg+1],"no") == 0) cluster_check = 0;
2443 else error->all(FLERR,"Illegal neigh_modify command");
2444 iarg += 2;
2445
2446 } else if (strcmp(arg[iarg],"include") == 0) {
2447 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2448 includegroup = group->find(arg[iarg+1]);
2449 if (includegroup < 0)
2450 error->all(FLERR,"Invalid group ID in neigh_modify command");
2451 if (includegroup && (atom->firstgroupname == nullptr ||
2452 strcmp(arg[iarg+1],atom->firstgroupname) != 0))
2453 error->all(FLERR,
2454 "Neigh_modify include group != atom_modify first group");
2455 iarg += 2;
2456
2457 } else if (strcmp(arg[iarg],"exclude") == 0) {
2458 if (iarg+2 > narg) error->all(FLERR,"Illegal neigh_modify command");
2459
2460 if (strcmp(arg[iarg+1],"type") == 0) {
2461 if (iarg+4 > narg) error->all(FLERR,"Illegal neigh_modify command");
2462 if (nex_type == maxex_type) {
2463 maxex_type += EXDELTA;
2464 memory->grow(ex1_type,maxex_type,"neigh:ex1_type");
2465 memory->grow(ex2_type,maxex_type,"neigh:ex2_type");
2466 }
2467 ex1_type[nex_type] = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
2468 ex2_type[nex_type] = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
2469 nex_type++;
2470 iarg += 4;
2471
2472 } else if (strcmp(arg[iarg+1],"group") == 0) {
2473 if (iarg+4 > narg) error->all(FLERR,"Illegal neigh_modify command");
2474 if (nex_group == maxex_group) {
2475 maxex_group += EXDELTA;
2476 memory->grow(ex1_group,maxex_group,"neigh:ex1_group");
2477 memory->grow(ex2_group,maxex_group,"neigh:ex2_group");
2478 }
2479 ex1_group[nex_group] = group->find(arg[iarg+2]);
2480 ex2_group[nex_group] = group->find(arg[iarg+3]);
2481 if (ex1_group[nex_group] == -1 || ex2_group[nex_group] == -1)
2482 error->all(FLERR,"Invalid group ID in neigh_modify command");
2483 nex_group++;
2484 iarg += 4;
2485
2486 } else if (strcmp(arg[iarg+1],"molecule/inter") == 0 ||
2487 strcmp(arg[iarg+1],"molecule/intra") == 0) {
2488 if (iarg+3 > narg) error->all(FLERR,"Illegal neigh_modify command");
2489 if (atom->molecule_flag == 0)
2490 error->all(FLERR,"Neigh_modify exclude molecule "
2491 "requires atom attribute molecule");
2492 if (nex_mol == maxex_mol) {
2493 maxex_mol += EXDELTA;
2494 memory->grow(ex_mol_group,maxex_mol,"neigh:ex_mol_group");
2495 if (lmp->kokkos)
2496 grow_ex_mol_intra_kokkos();
2497 else
2498 memory->grow(ex_mol_intra,maxex_mol,"neigh:ex_mol_intra");
2499 }
2500 ex_mol_group[nex_mol] = group->find(arg[iarg+2]);
2501 if (ex_mol_group[nex_mol] == -1)
2502 error->all(FLERR,"Invalid group ID in neigh_modify command");
2503 if (strcmp(arg[iarg+1],"molecule/intra") == 0)
2504 ex_mol_intra[nex_mol] = 1;
2505 else
2506 ex_mol_intra[nex_mol] = 0;
2507 nex_mol++;
2508 iarg += 3;
2509
2510 } else if (strcmp(arg[iarg+1],"none") == 0) {
2511 nex_type = nex_group = nex_mol = 0;
2512 iarg += 2;
2513
2514 } else error->all(FLERR,"Illegal neigh_modify command");
2515 } else if (strcmp(arg[iarg],"collection/interval") == 0) {
2516 if (style != Neighbor::MULTI)
2517 error->all(FLERR,"Cannot use collection/interval command without multi setting");
2518
2519 if (iarg+2 > narg)
2520 error->all(FLERR,"Invalid collection/interval command");
2521 ncollections = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
2522 if (ncollections < 1)
2523 error->all(FLERR,"Invalid collection/interval command");
2524 if (iarg+1+ncollections > narg)
2525 error->all(FLERR,"Invalid collection/interval command");
2526
2527 int i;
2528
2529 // Invalidate old user cutoffs
2530
2531 comm->ncollections_cutoff = 0;
2532 interval_collection_flag = 1;
2533 custom_collection_flag = 1;
2534 memory->grow(collection2cut,ncollections,"neigh:collection2cut");
2535
2536 // Set upper cutoff for each collection
2537
2538 double cut_interval;
2539 for (i = 0; i < ncollections; i++){
2540 cut_interval = utils::numeric(FLERR,arg[iarg+2+i],false,lmp);
2541 collection2cut[i] = cut_interval;
2542
2543 if (i != 0)
2544 if (collection2cut[i-1] >= collection2cut[i])
2545 error->all(FLERR,"Nonsequential interval cutoffs in collection/interval setting");
2546 }
2547
2548 iarg += 2 + ncollections;
2549 } else if (strcmp(arg[iarg],"collection/type") == 0) {
2550 if (style != Neighbor::MULTI)
2551 error->all(FLERR,"Cannot use collection/type command without multi setting");
2552
2553 if (iarg+2 > narg)
2554 error->all(FLERR,"Invalid collection/type command");
2555 ncollections = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
2556 if (ncollections < 1)
2557 error->all(FLERR,"Invalid collection/interval command");
2558 if (iarg+1+ncollections > narg)
2559 error->all(FLERR,"Invalid collection/type command");
2560
2561 int ntypes = atom->ntypes;
2562 int nlo, nhi, i, k;
2563
2564 // Invalidate old user cutoffs
2565
2566 comm->ncollections_cutoff = 0;
2567 interval_collection_flag = 0;
2568 custom_collection_flag = 1;
2569 if (!type2collection)
2570 memory->create(type2collection,ntypes+1,"neigh:type2collection");
2571
2572 // Erase previous mapping
2573
2574 for (i = 1; i <= ntypes; i++)
2575 type2collection[i] = -1;
2576
2577 // For each custom range, define mapping for types in interval
2578
2579 for (i = 0; i < ncollections; i++){
2580 std::vector<std::string> words = Tokenizer(arg[iarg+2+i], ",").as_vector();
2581 for (const auto &word : words) {
2582 utils::bounds(FLERR,word,1,ntypes,nlo,nhi,error);
2583 for (k = nlo; k <= nhi; k++) {
2584 if (type2collection[k] != -1)
2585 error->all(FLERR,"Type specified more than once in collection/type commnd");
2586 type2collection[k] = i;
2587 }
2588 }
2589 }
2590
2591 // Check for undefined atom type
2592
2593 for (i = 1; i <= ntypes; i++){
2594 if (type2collection[i] == -1) {
2595 error->all(FLERR,"Type missing in collection/type commnd");
2596 }
2597 }
2598
2599 iarg += 2 + ncollections;
2600 } else error->all(FLERR,"Illegal neigh_modify command");
2601 }
2602 }
2603
2604 /* ----------------------------------------------------------------------
2605 convenience function to allow modifying parameters from a single string
2606 ------------------------------------------------------------------------- */
2607
modify_params(const std::string & modcmd)2608 void Neighbor::modify_params(const std::string &modcmd)
2609 {
2610 auto args = utils::split_words(modcmd);
2611 char **newarg = new char*[args.size()];
2612 int i=0;
2613 for (const auto &arg : args) {
2614 newarg[i++] = (char *)arg.c_str();
2615 }
2616 modify_params(args.size(),newarg);
2617 delete[] newarg;
2618 }
2619
2620 /* ----------------------------------------------------------------------
2621 remove the first group-group exclusion matching group1, group2
2622 ------------------------------------------------------------------------- */
2623
exclusion_group_group_delete(int group1,int group2)2624 void Neighbor::exclusion_group_group_delete(int group1, int group2)
2625 {
2626 int m, mlast;
2627 for (m = 0; m < nex_group; m++)
2628 if (ex1_group[m] == group1 && ex2_group[m] == group2 )
2629 break;
2630
2631 mlast = m;
2632 if (mlast == nex_group)
2633 error->all(FLERR,"Unable to find group-group exclusion");
2634
2635 for (m = mlast+1; m < nex_group; m++) {
2636 ex1_group[m-1] = ex1_group[m];
2637 ex2_group[m-1] = ex2_group[m];
2638 ex1_bit[m-1] = ex1_bit[m];
2639 ex2_bit[m-1] = ex2_bit[m];
2640 }
2641 nex_group--;
2642 }
2643
2644
2645 /* ----------------------------------------------------------------------
2646 return the value of exclude - used to check compatibility with GPU
2647 ------------------------------------------------------------------------- */
2648
exclude_setting()2649 int Neighbor::exclude_setting()
2650 {
2651 return exclude;
2652 }
2653
2654 /* ----------------------------------------------------------------------
2655 check if any of the old requested neighbor lists are full
2656 ------------------------------------------------------------------------- */
2657
any_full()2658 int Neighbor::any_full()
2659 {
2660 int any_full = 0;
2661 for (int i = 0; i < old_nrequest; i++) {
2662 if (old_requests[i]->full) any_full = 1;
2663 }
2664 return any_full;
2665 }
2666
2667 /* ----------------------------------------------------------------------
2668 populate collection array for multi starting at the index istart
2669 ------------------------------------------------------------------------- */
2670
build_collection(int istart)2671 void Neighbor::build_collection(int istart)
2672 {
2673 if (style != Neighbor::MULTI)
2674 error->all(FLERR, "Cannot define atom collections without neighbor style multi");
2675
2676 int nmax = atom->nlocal+atom->nghost;
2677 if (nmax > nmax_collection) {
2678 nmax_collection = nmax+DELTA_PERATOM;
2679 memory->grow(collection, nmax_collection, "neigh:collection");
2680 }
2681
2682 if (finite_cut_flag) {
2683 double cut;
2684 int icollection;
2685 for (int i = istart; i < nmax; i++){
2686 cut = force->pair->atom2cut(i);
2687 collection[i] = -1;
2688
2689 for (icollection = 0; icollection < ncollections; icollection++){
2690 if (collection2cut[icollection] >= cut) {
2691 collection[i] = icollection;
2692 break;
2693 }
2694 }
2695
2696 if (collection[i] == -1)
2697 error->one(FLERR, "Atom cutoff exceeds interval cutoffs for multi");
2698 }
2699 } else {
2700 int *type = atom->type;
2701 for (int i = istart; i < nmax; i++){
2702 collection[i] = type2collection[type[i]];
2703 }
2704 }
2705 }
2706
2707 /* ----------------------------------------------------------------------
2708 return # of bytes of allocated memory
2709 ------------------------------------------------------------------------- */
2710
memory_usage()2711 double Neighbor::memory_usage()
2712 {
2713 double bytes = 0;
2714 bytes += memory->usage(xhold,maxhold,3);
2715
2716 for (int i = 0; i < nlist; i++)
2717 if (lists[i]) bytes += lists[i]->memory_usage();
2718 for (int i = 0; i < nstencil; i++)
2719 bytes += neigh_stencil[i]->memory_usage();
2720 for (int i = 0; i < nbin; i++)
2721 bytes += neigh_bin[i]->memory_usage();
2722
2723 if (neigh_bond) bytes += neigh_bond->memory_usage();
2724 if (neigh_angle) bytes += neigh_angle->memory_usage();
2725 if (neigh_dihedral) bytes += neigh_dihedral->memory_usage();
2726 if (neigh_improper) bytes += neigh_improper->memory_usage();
2727
2728 return bytes;
2729 }
2730