1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 #include "atom.h"
16 #include "atom_vec.h"
17 #include "style_atom.h" // IWYU pragma: keep
18
19 #include "comm.h"
20 #include "compute.h"
21 #include "domain.h"
22 #include "error.h"
23 #include "fix.h"
24 #include "force.h"
25 #include "group.h"
26 #include "input.h"
27 #include "math_const.h"
28 #include "memory.h"
29 #include "modify.h"
30 #include "molecule.h"
31 #include "neighbor.h"
32 #include "update.h"
33 #include "variable.h"
34
35 #include "library.h"
36
37 #include <algorithm>
38 #include <cstring>
39
40 #ifdef LMP_INTEL
41 #include "neigh_request.h"
42 #endif
43
44 #ifdef LMP_GPU
45 #include "fix_gpu.h"
46 #include <cmath>
47 #endif
48
49 using namespace LAMMPS_NS;
50 using namespace MathConst;
51
52 #define DELTA 1
53 #define DELTA_PERATOM 64
54 #define EPSILON 1.0e-6
55
56 /* ---------------------------------------------------------------------- */
57
58 /** \class LAMMPS_NS::Atom
59 * \brief Class to provide access to atom data
60
61 \verbatim embed:rst
62 The Atom class provides access to atom style related global settings and
63 per-atom data that is stored with atoms and migrates with them from
64 sub-domain to sub-domain as atoms move around. This includes topology
65 data, which is stored with either one specific atom or all atoms involved
66 depending on the settings of the :doc:`newton command <newton>`.
67
68 The actual per-atom data is allocated and managed by one of the various
69 classes derived from the AtomVec class as determined by
70 the :doc:`atom_style command <atom_style>`. The pointers in the Atom class
71 are updated by the AtomVec class as needed.
72 \endverbatim
73 */
74
75 /** Atom class constructor
76 *
77 * This resets and initialized all kinds of settings,
78 * parameters, and pointer variables for per-atom arrays.
79 * This also initializes the factory for creating
80 * instances of classes derived from the AtomVec base
81 * class, which correspond to the selected atom style.
82 *
83 * \param lmp pointer to the base LAMMPS class */
84
Atom(LAMMPS * lmp)85 Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
86 {
87 natoms = 0;
88 nlocal = nghost = nmax = 0;
89 ntypes = 0;
90 nellipsoids = nlines = ntris = nbodies = 0;
91 nbondtypes = nangletypes = ndihedraltypes = nimpropertypes = 0;
92 nbonds = nangles = ndihedrals = nimpropers = 0;
93
94 firstgroupname = nullptr;
95 sortfreq = 1000;
96 nextsort = 0;
97 userbinsize = 0.0;
98 maxbin = maxnext = 0;
99 binhead = nullptr;
100 next = permute = nullptr;
101
102 // data structure with info on per-atom vectors/arrays
103
104 nperatom = maxperatom = 0;
105 peratom = nullptr;
106
107 // --------------------------------------------------------------------
108 // 1st customization section: customize by adding new per-atom variables
109
110 tag = nullptr;
111 type = mask = nullptr;
112 image = nullptr;
113 x = v = f = nullptr;
114
115 // charged and dipolar particles
116
117 q = nullptr;
118 mu = nullptr;
119
120 // finite-size particles
121
122 omega = angmom = torque = nullptr;
123 radius = rmass = nullptr;
124 ellipsoid = line = tri = body = nullptr;
125
126 // molecular systems
127
128 molecule = nullptr;
129 molindex = molatom = nullptr;
130
131 bond_per_atom = extra_bond_per_atom = 0;
132 num_bond = nullptr;
133 bond_type = nullptr;
134 bond_atom = nullptr;
135
136 angle_per_atom = extra_angle_per_atom = 0;
137 num_angle = nullptr;
138 angle_type = nullptr;
139 angle_atom1 = angle_atom2 = angle_atom3 = nullptr;
140
141 dihedral_per_atom = extra_dihedral_per_atom = 0;
142 num_dihedral = nullptr;
143 dihedral_type = nullptr;
144 dihedral_atom1 = dihedral_atom2 = dihedral_atom3 = dihedral_atom4 = nullptr;
145
146 improper_per_atom = extra_improper_per_atom = 0;
147 num_improper = nullptr;
148 improper_type = nullptr;
149 improper_atom1 = improper_atom2 = improper_atom3 = improper_atom4 = nullptr;
150
151 maxspecial = 1;
152 nspecial = nullptr;
153 special = nullptr;
154
155 // PERI package
156
157 vfrac = s0 = nullptr;
158 x0 = nullptr;
159
160 // SPIN package
161
162 sp = fm = fm_long = nullptr;
163
164 // EFF and AWPMD packages
165
166 spin = nullptr;
167 eradius = ervel = erforce = nullptr;
168 ervelforce = nullptr;
169 cs = csforce = vforce = nullptr;
170 etag = nullptr;
171
172 // CG-DNA package
173
174 id5p = nullptr;
175
176 // DPD-REACT package
177
178 uCond = uMech = uChem = uCG = uCGnew = nullptr;
179 duChem = dpdTheta = nullptr;
180
181 // MESO package
182
183 cc = cc_flux = nullptr;
184 edpd_temp = edpd_flux = edpd_cv = nullptr;
185
186 // MESONT package
187
188 length = nullptr;
189 buckling = nullptr;
190 bond_nt = nullptr;
191
192 // MACHDYN package
193
194 contact_radius = nullptr;
195 smd_data_9 = nullptr;
196 smd_stress = nullptr;
197 eff_plastic_strain = nullptr;
198 eff_plastic_strain_rate = nullptr;
199 damage = nullptr;
200
201 // SPH package
202
203 rho = drho = esph = desph = cv = nullptr;
204 vest = nullptr;
205
206 // DIELECTRIC package
207
208 area = ed = em = epsilon = curvature = q_unscaled = nullptr;
209
210 // end of customization section
211 // --------------------------------------------------------------------
212
213 // user-defined molecules
214
215 nmolecule = 0;
216 molecules = nullptr;
217
218 // custom atom arrays
219
220 nivector = ndvector = niarray = ndarray = 0;
221 ivector = nullptr;
222 dvector = nullptr;
223 iarray = nullptr;
224 darray = nullptr;
225 icols = dcols = nullptr;
226 ivname = dvname = ianame = daname = nullptr;
227
228 // initialize atom style and array existence flags
229
230 set_atomflag_defaults();
231
232 // initialize peratom data structure
233
234 peratom_create();
235
236 // ntype-length arrays
237
238 mass = nullptr;
239 mass_setflag = nullptr;
240
241 // callback lists & extra restart info
242
243 nextra_grow = nextra_restart = nextra_border = 0;
244 extra_grow = extra_restart = extra_border = nullptr;
245 nextra_grow_max = nextra_restart_max = nextra_border_max = 0;
246 nextra_store = 0;
247 extra = nullptr;
248
249 // default atom ID and mapping values
250
251 tag_enable = 1;
252 map_style = map_user = MAP_NONE;
253 map_tag_max = -1;
254 map_maxarray = map_nhash = map_nbucket = -1;
255
256 max_same = 0;
257 sametag = nullptr;
258 map_array = nullptr;
259 map_bucket = nullptr;
260 map_hash = nullptr;
261
262 unique_tags = nullptr;
263
264 atom_style = nullptr;
265 avec = nullptr;
266
267 avec_map = new AtomVecCreatorMap();
268
269 #define ATOM_CLASS
270 #define AtomStyle(key,Class) \
271 (*avec_map)[#key] = &avec_creator<Class>;
272 #include "style_atom.h" // IWYU pragma: keep
273 #undef AtomStyle
274 #undef ATOM_CLASS
275 }
276
277 /* ---------------------------------------------------------------------- */
278
~Atom()279 Atom::~Atom()
280 {
281 delete [] atom_style;
282 delete avec;
283 delete avec_map;
284
285 delete [] firstgroupname;
286 memory->destroy(binhead);
287 memory->destroy(next);
288 memory->destroy(permute);
289
290 memory->destroy(tag);
291 memory->destroy(type);
292 memory->destroy(mask);
293 memory->destroy(image);
294 memory->destroy(x);
295 memory->destroy(v);
296 memory->destroy(f);
297
298 // delete peratom data struct
299
300 for (int i = 0; i < nperatom; i++)
301 delete [] peratom[i].name;
302 memory->sfree(peratom);
303
304 // delete custom atom arrays
305
306 for (int i = 0; i < nivector; i++) {
307 delete [] ivname[i];
308 memory->destroy(ivector[i]);
309 }
310 for (int i = 0; i < ndvector; i++) {
311 delete [] dvname[i];
312 if (dvector) // (needed for Kokkos)
313 memory->destroy(dvector[i]);
314 }
315 for (int i = 0; i < niarray; i++) {
316 delete [] ianame[i];
317 memory->destroy(iarray[i]);
318 }
319 for (int i = 0; i < ndarray; i++) {
320 delete [] daname[i];
321 memory->destroy(darray[i]);
322 }
323
324 memory->sfree(ivname);
325 memory->sfree(dvname);
326 memory->sfree(ianame);
327 memory->sfree(daname);
328 memory->sfree(ivector);
329 memory->sfree(dvector);
330 memory->sfree(iarray);
331 memory->sfree(darray);
332 memory->sfree(icols);
333 memory->sfree(dcols);
334
335 // delete user-defined molecules
336
337 for (int i = 0; i < nmolecule; i++) delete molecules[i];
338 memory->sfree(molecules);
339
340 // delete per-type arrays
341
342 delete [] mass;
343 delete [] mass_setflag;
344
345 // delete extra arrays
346
347 memory->destroy(extra_grow);
348 memory->destroy(extra_restart);
349 memory->destroy(extra_border);
350 memory->destroy(extra);
351
352 // delete mapping data structures
353
354 Atom::map_delete();
355
356 delete unique_tags;
357 }
358
359 /* ----------------------------------------------------------------------
360 copy modify settings from old Atom class to current Atom class
361 ------------------------------------------------------------------------- */
362
settings(Atom * old)363 void Atom::settings(Atom *old)
364 {
365 tag_enable = old->tag_enable;
366 map_user = old->map_user;
367 map_style = old->map_style;
368 sortfreq = old->sortfreq;
369 userbinsize = old->userbinsize;
370 if (old->firstgroupname)
371 firstgroupname = utils::strdup(old->firstgroupname);
372 }
373
374 /* ----------------------------------------------------------------------
375 one-time creation of peratom data structure
376 ------------------------------------------------------------------------- */
377
peratom_create()378 void Atom::peratom_create()
379 {
380 for (int i = 0; i < nperatom; i++)
381 delete [] peratom[i].name;
382 memory->sfree(peratom);
383
384 peratom = nullptr;
385 nperatom = maxperatom = 0;
386
387 // --------------------------------------------------------------------
388 // 2nd customization section: add peratom variables here, order does not matter
389 // register tagint & imageint variables as INT or BIGINT
390
391 int tagintsize = INT;
392 if (sizeof(tagint) == 8) tagintsize = BIGINT;
393 int imageintsize = INT;
394 if (sizeof(imageint) == 8) imageintsize = BIGINT;
395
396 add_peratom("id",&tag,tagintsize,0);
397 add_peratom("type",&type,INT,0);
398 add_peratom("mask",&mask,INT,0);
399 add_peratom("image",&image,imageintsize,0);
400
401 add_peratom("x",&x,DOUBLE,3);
402 add_peratom("v",&v,DOUBLE,3);
403 add_peratom("f",&f,DOUBLE,3,1); // set per-thread flag
404
405 add_peratom("rmass",&rmass,DOUBLE,0);
406 add_peratom("q",&q,DOUBLE,0);
407 add_peratom("mu",&mu,DOUBLE,4);
408 add_peratom("mu3",&mu,DOUBLE,3); // just first 3 values of mu[4]
409
410 // finite size particles
411
412 add_peratom("radius",&radius,DOUBLE,0);
413 add_peratom("omega",&omega,DOUBLE,3);
414 add_peratom("torque",&torque,DOUBLE,3,1); // set per-thread flag
415 add_peratom("angmom",&angmom,DOUBLE,3);
416
417 add_peratom("ellipsoid",&ellipsoid,INT,0);
418 add_peratom("line",&line,INT,0);
419 add_peratom("tri",&tri,INT,0);
420 add_peratom("body",&body,INT,0);
421
422 // MOLECULE package
423
424 add_peratom("molecule",&molecule,tagintsize,0);
425 add_peratom("molindex",&molindex,INT,0);
426 add_peratom("molatom",&molatom,INT,0);
427
428 add_peratom("nspecial",&nspecial,INT,3);
429 add_peratom_vary("special",&special,tagintsize,&maxspecial,&nspecial,3);
430
431 add_peratom("num_bond",&num_bond,INT,0);
432 add_peratom_vary("bond_type",&bond_type,INT,&bond_per_atom,&num_bond);
433 add_peratom_vary("bond_atom",&bond_atom,tagintsize,&bond_per_atom,&num_bond);
434
435 add_peratom("num_angle",&num_angle,INT,0);
436 add_peratom_vary("angle_type",&angle_type,INT,&angle_per_atom,&num_angle);
437 add_peratom_vary("angle_atom1",&angle_atom1,tagintsize,
438 &angle_per_atom,&num_angle);
439 add_peratom_vary("angle_atom2",&angle_atom2,tagintsize,
440 &angle_per_atom,&num_angle);
441 add_peratom_vary("angle_atom3",&angle_atom3,tagintsize,
442 &angle_per_atom,&num_angle);
443
444 add_peratom("num_dihedral",&num_dihedral,INT,0);
445 add_peratom_vary("dihedral_type",&dihedral_type,INT,
446 &dihedral_per_atom,&num_dihedral);
447 add_peratom_vary("dihedral_atom1",&dihedral_atom1,tagintsize,
448 &dihedral_per_atom,&num_dihedral);
449 add_peratom_vary("dihedral_atom2",&dihedral_atom2,tagintsize,
450 &dihedral_per_atom,&num_dihedral);
451 add_peratom_vary("dihedral_atom3",&dihedral_atom3,tagintsize,
452 &dihedral_per_atom,&num_dihedral);
453 add_peratom_vary("dihedral_atom4",&dihedral_atom4,tagintsize,
454 &dihedral_per_atom,&num_dihedral);
455
456 add_peratom("num_improper",&num_improper,INT,0);
457 add_peratom_vary("improper_type",&improper_type,INT,
458 &improper_per_atom,&num_improper);
459 add_peratom_vary("improper_atom1",&improper_atom1,tagintsize,
460 &improper_per_atom,&num_improper);
461 add_peratom_vary("improper_atom2",&improper_atom2,tagintsize,
462 &improper_per_atom,&num_improper);
463 add_peratom_vary("improper_atom3",&improper_atom3,tagintsize,
464 &improper_per_atom,&num_improper);
465 add_peratom_vary("improper_atom4",&improper_atom4,tagintsize,
466 &improper_per_atom,&num_improper);
467
468 // PERI package
469
470 add_peratom("vfrac",&vfrac,DOUBLE,0);
471 add_peratom("s0",&s0,DOUBLE,0);
472 add_peratom("x0",&x0,DOUBLE,3);
473
474 // SPIN package
475
476 add_peratom("sp",&sp,DOUBLE,4);
477 add_peratom("fm",&fm,DOUBLE,3,1);
478 add_peratom("fm_long",&fm_long,DOUBLE,3,1);
479
480 // EFF package
481
482 add_peratom("spin",&spin,INT,0);
483 add_peratom("eradius",&eradius,DOUBLE,0);
484 add_peratom("ervel",&ervel,DOUBLE,0);
485 add_peratom("erforce",&erforce,DOUBLE,0,1); // set per-thread flag
486
487 // AWPMD package
488
489 add_peratom("cs",&cs,DOUBLE,2);
490 add_peratom("csforce",&csforce,DOUBLE,2);
491 add_peratom("vforce",&vforce,DOUBLE,3);
492 add_peratom("ervelforce",&ervelforce,DOUBLE,0);
493 add_peratom("etag",&etag,INT,0);
494
495 // CG-DNA package
496
497 add_peratom("id5p",&id5p,tagintsize,0);
498
499 // DPD-REACT package
500
501 add_peratom("dpdTheta",&dpdTheta,DOUBLE,0);
502 add_peratom("uCond",&uCond,DOUBLE,0);
503 add_peratom("uMech",&uMech,DOUBLE,0);
504 add_peratom("uChem",&uChem,DOUBLE,0);
505 add_peratom("uCG",&uCG,DOUBLE,0);
506 add_peratom("uCGnew",&uCGnew,DOUBLE,0);
507 add_peratom("duChem",&duChem,DOUBLE,0);
508
509 // MESO package
510
511 add_peratom("edpd_cv",&edpd_cv,DOUBLE,0);
512 add_peratom("edpd_temp",&edpd_temp,DOUBLE,0);
513 add_peratom("vest_temp",&vest_temp,DOUBLE,0);
514 add_peratom("edpd_flux",&edpd_flux,DOUBLE,0,1); // set per-thread flag
515 add_peratom("cc",&cc,DOUBLE,1);
516 add_peratom("cc_flux",&cc_flux,DOUBLE,1,1); // set per-thread flag
517
518 // MESONT package
519
520 add_peratom("length",&length,DOUBLE,0);
521 add_peratom("buckling",&buckling,INT,0);
522 add_peratom("bond_nt",&bond_nt,tagintsize,2);
523
524 // SPH package
525
526 add_peratom("rho",&rho,DOUBLE,0);
527 add_peratom("drho",&drho,DOUBLE,0,1); // set per-thread flag
528 add_peratom("esph",&esph,DOUBLE,0);
529 add_peratom("desph",&desph,DOUBLE,0,1); // set per-thread flag
530 add_peratom("vest",&vest,DOUBLE,3);
531 add_peratom("cv",&cv,DOUBLE,0);
532
533 // MACHDYN package
534
535 add_peratom("contact_radius",&contact_radius,DOUBLE,0);
536 add_peratom("smd_data_9",&smd_data_9,DOUBLE,1);
537 add_peratom("smd_stress",&smd_stress,DOUBLE,1);
538 add_peratom("eff_plastic_strain",&eff_plastic_strain,DOUBLE,0);
539 add_peratom("eff_plastic_strain_rate",&eff_plastic_strain_rate,DOUBLE,0);
540 add_peratom("damage",&damage,DOUBLE,0);
541
542 // DIELECTRIC package
543
544 add_peratom("area",&area,DOUBLE,0);
545 add_peratom("ed",&ed,DOUBLE,0);
546 add_peratom("em",&em,DOUBLE,0);
547 add_peratom("epsilon",&epsilon,DOUBLE,0);
548 add_peratom("curvature",&curvature,DOUBLE,0);
549 add_peratom("q_unscaled",&q_unscaled,DOUBLE,0);
550
551 // end of customization section
552 // --------------------------------------------------------------------
553 }
554
555 /* ----------------------------------------------------------------------
556 add info for a single per-atom vector/array to PerAtom data struct
557 cols = 0: per-atom vector
558 cols = N: static per-atom array with N columns
559 use add_peratom_vary() when column count varies per atom
560 ------------------------------------------------------------------------- */
561
add_peratom(const char * name,void * address,int datatype,int cols,int threadflag)562 void Atom::add_peratom(const char *name, void *address,
563 int datatype, int cols, int threadflag)
564 {
565 if (nperatom == maxperatom) {
566 maxperatom += DELTA_PERATOM;
567 peratom = (PerAtom *)
568 memory->srealloc(peratom,maxperatom*sizeof(PerAtom),"atom:peratom");
569 }
570
571 peratom[nperatom].name = utils::strdup(name);
572 peratom[nperatom].address = address;
573 peratom[nperatom].datatype = datatype;
574 peratom[nperatom].cols = cols;
575 peratom[nperatom].threadflag = threadflag;
576 peratom[nperatom].address_length = nullptr;
577
578 nperatom++;
579 }
580
581 /* ----------------------------------------------------------------------
582 change the column count of an existing peratom array entry
583 allows atom_style to specify column count as an argument
584 see atom_style tdpd as an example
585 ------------------------------------------------------------------------- */
586
add_peratom_change_columns(const char * name,int cols)587 void Atom::add_peratom_change_columns(const char *name, int cols)
588 {
589 for (int i = 0; i < nperatom; i++) {
590 if (strcmp(name,peratom[i].name) == 0) {
591 peratom[i].cols = cols;
592 return;
593 }
594 }
595 error->all(FLERR,"Could not find name of peratom array for column change");
596 }
597
598 /* ----------------------------------------------------------------------
599 add info for a single per-atom array to PerAtom data struct
600 cols = address of int variable with max columns per atom
601 for collength = 0:
602 length = address of peratom vector with column count per atom
603 e.g. num_bond
604 for collength = N:
605 length = address of peratom array with column count per atom
606 collength = index of column (1 to N) in peratom array with count
607 e.g. nspecial
608 ------------------------------------------------------------------------- */
609
add_peratom_vary(const char * name,void * address,int datatype,int * cols,void * length,int collength)610 void Atom::add_peratom_vary(const char *name, void *address,
611 int datatype, int *cols, void *length, int collength)
612 {
613 if (nperatom == maxperatom) {
614 maxperatom += DELTA_PERATOM;
615 peratom = (PerAtom *)
616 memory->srealloc(peratom,maxperatom*sizeof(PerAtom),"atom:peratom");
617 }
618
619 peratom[nperatom].name = utils::strdup(name);
620 peratom[nperatom].address = address;
621 peratom[nperatom].datatype = datatype;
622 peratom[nperatom].cols = -1;
623 peratom[nperatom].threadflag = 0;
624 peratom[nperatom].address_maxcols = cols;
625 peratom[nperatom].address_length = length;
626 peratom[nperatom].collength = collength;
627
628 nperatom++;
629 }
630
631 /* ----------------------------------------------------------------------
632 add info for a single per-atom array to PerAtom data struct
633 ------------------------------------------------------------------------- */
634
set_atomflag_defaults()635 void Atom::set_atomflag_defaults()
636 {
637 // --------------------------------------------------------------------
638 // 3rd customization section: customize by adding new flag
639 // identical list as 2nd customization in atom.h
640
641 sphere_flag = ellipsoid_flag = line_flag = tri_flag = body_flag = 0;
642 peri_flag = electron_flag = 0;
643 wavepacket_flag = sph_flag = 0;
644 molecule_flag = molindex_flag = molatom_flag = 0;
645 q_flag = mu_flag = 0;
646 rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0;
647 vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0;
648 cs_flag = csforce_flag = vforce_flag = ervelforce_flag = etag_flag = 0;
649 rho_flag = esph_flag = cv_flag = vest_flag = 0;
650 dpd_flag = edpd_flag = tdpd_flag = 0;
651 sp_flag = 0;
652 x0_flag = 0;
653 smd_flag = damage_flag = 0;
654 mesont_flag = 0;
655 contact_radius_flag = smd_data_9_flag = smd_stress_flag = 0;
656 eff_plastic_strain_flag = eff_plastic_strain_rate_flag = 0;
657
658 pdscale = 1.0;
659 }
660
661 /* ----------------------------------------------------------------------
662 create an AtomVec style
663 called from lammps.cpp, input script, restart file, replicate
664 ------------------------------------------------------------------------- */
665
create_avec(const std::string & style,int narg,char ** arg,int trysuffix)666 void Atom::create_avec(const std::string &style, int narg, char **arg, int trysuffix)
667 {
668 delete[] atom_style;
669 if (avec) delete avec;
670 atom_style = nullptr;
671 avec = nullptr;
672
673 // unset atom style and array existence flags
674 // may have been set by old avec
675
676 set_atomflag_defaults();
677
678 // create instance of AtomVec
679 // use grow() to initialize atom-based arrays to length 1
680 // so that x[0][0] can always be referenced even if proc has no atoms
681
682 int sflag;
683 avec = new_avec(style,trysuffix,sflag);
684 avec->store_args(narg,arg);
685 avec->process_args(narg,arg);
686 avec->grow(1);
687
688 if (sflag) {
689 std::string estyle = style + "/";
690 if (sflag == 1) estyle += lmp->suffix;
691 else estyle += lmp->suffix2;
692 atom_style = utils::strdup(estyle);
693 } else {
694 atom_style = utils::strdup(style);
695 }
696
697 // if molecular system:
698 // atom IDs must be defined
699 // force atom map to be created
700 // map style will be reset to array vs hash to by map_init()
701
702 molecular = avec->molecular;
703 if ((molecular != Atom::ATOMIC) && (tag_enable == 0))
704 error->all(FLERR,"Atom IDs must be used for molecular systems");
705 if (molecular != Atom::ATOMIC) map_style = MAP_YES;
706 }
707
708 /* ----------------------------------------------------------------------
709 generate an AtomVec class, first with suffix appended
710 ------------------------------------------------------------------------- */
711
new_avec(const std::string & style,int trysuffix,int & sflag)712 AtomVec *Atom::new_avec(const std::string &style, int trysuffix, int &sflag)
713 {
714 if (trysuffix && lmp->suffix_enable) {
715 if (lmp->suffix) {
716 sflag = 1;
717 std::string estyle = style + "/" + lmp->suffix;
718 if (avec_map->find(estyle) != avec_map->end()) {
719 AtomVecCreator &avec_creator = (*avec_map)[estyle];
720 return avec_creator(lmp);
721 }
722 }
723
724 if (lmp->suffix2) {
725 sflag = 2;
726 std::string estyle = style + "/" + lmp->suffix2;
727 if (avec_map->find(estyle) != avec_map->end()) {
728 AtomVecCreator &avec_creator = (*avec_map)[estyle];
729 return avec_creator(lmp);
730 }
731 }
732 }
733
734 sflag = 0;
735 if (avec_map->find(style) != avec_map->end()) {
736 AtomVecCreator &avec_creator = (*avec_map)[style];
737 return avec_creator(lmp);
738 }
739
740 error->all(FLERR,utils::check_packages_for_style("atom",style,lmp));
741 return nullptr;
742 }
743
744 /* ----------------------------------------------------------------------
745 one instance per AtomVec style in style_atom.h
746 ------------------------------------------------------------------------- */
747
748 template <typename T>
avec_creator(LAMMPS * lmp)749 AtomVec *Atom::avec_creator(LAMMPS *lmp)
750 {
751 return new T(lmp);
752 }
753
754 /* ---------------------------------------------------------------------- */
755
init()756 void Atom::init()
757 {
758 // delete extra array since it doesn't persist past first run
759
760 if (nextra_store) {
761 memory->destroy(extra);
762 extra = nullptr;
763 nextra_store = 0;
764 }
765
766 // check arrays that are atom type in length
767
768 check_mass(FLERR);
769
770 // setup of firstgroup
771
772 if (firstgroupname) {
773 firstgroup = group->find(firstgroupname);
774 if (firstgroup < 0)
775 error->all(FLERR,"Could not find atom_modify first group ID");
776 } else firstgroup = -1;
777
778 // init AtomVec
779
780 avec->init();
781 }
782
783 /* ---------------------------------------------------------------------- */
784
setup()785 void Atom::setup()
786 {
787 // setup bins for sorting
788 // cannot do this in init() because uses neighbor cutoff
789
790 if (sortfreq > 0) setup_sort_bins();
791 }
792
793 /* ----------------------------------------------------------------------
794 return ptr to AtomVec class if matches style or to matching hybrid sub-class
795 return nullptr if no match
796 ------------------------------------------------------------------------- */
797
style_match(const char * style)798 AtomVec *Atom::style_match(const char *style)
799 {
800 if (strcmp(atom_style,style) == 0) return avec;
801 else if (strcmp(atom_style,"hybrid") == 0) {
802 auto avec_hybrid = (AtomVecHybrid *) avec;
803 for (int i = 0; i < avec_hybrid->nstyles; i++)
804 if (strcmp(avec_hybrid->keywords[i],style) == 0)
805 return avec_hybrid->styles[i];
806 }
807 return nullptr;
808 }
809
810 /* ----------------------------------------------------------------------
811 modify parameters of the atom style
812 some options can only be invoked before simulation box is defined
813 first and sort options cannot be used together
814 ------------------------------------------------------------------------- */
815
modify_params(int narg,char ** arg)816 void Atom::modify_params(int narg, char **arg)
817 {
818 if (narg == 0) error->all(FLERR,"Illegal atom_modify command");
819
820 int iarg = 0;
821 while (iarg < narg) {
822 if (strcmp(arg[iarg],"id") == 0) {
823 if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
824 if (domain->box_exist)
825 error->all(FLERR,
826 "Atom_modify id command after simulation box is defined");
827 if (strcmp(arg[iarg+1],"yes") == 0) tag_enable = 1;
828 else if (strcmp(arg[iarg+1],"no") == 0) tag_enable = 0;
829 else error->all(FLERR,"Illegal atom_modify command");
830 iarg += 2;
831 } else if (strcmp(arg[iarg],"map") == 0) {
832 if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
833 if (domain->box_exist)
834 error->all(FLERR,
835 "Atom_modify map command after simulation box is defined");
836 if (strcmp(arg[iarg+1],"array") == 0) map_user = 1;
837 else if (strcmp(arg[iarg+1],"hash") == 0) map_user = 2;
838 else if (strcmp(arg[iarg+1],"yes") == 0) map_user = 3;
839 else error->all(FLERR,"Illegal atom_modify command");
840 map_style = map_user;
841 iarg += 2;
842 } else if (strcmp(arg[iarg],"first") == 0) {
843 if (iarg+2 > narg) error->all(FLERR,"Illegal atom_modify command");
844 if (strcmp(arg[iarg+1],"all") == 0) {
845 delete [] firstgroupname;
846 firstgroupname = nullptr;
847 } else {
848 firstgroupname = utils::strdup(arg[iarg+1]);
849 sortfreq = 0;
850 }
851 iarg += 2;
852 } else if (strcmp(arg[iarg],"sort") == 0) {
853 if (iarg+3 > narg) error->all(FLERR,"Illegal atom_modify command");
854 sortfreq = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
855 userbinsize = utils::numeric(FLERR,arg[iarg+2],false,lmp);
856 if (sortfreq < 0 || userbinsize < 0.0)
857 error->all(FLERR,"Illegal atom_modify command");
858 if (sortfreq >= 0 && firstgroupname)
859 error->all(FLERR,"Atom_modify sort and first options "
860 "cannot be used together");
861 iarg += 3;
862 } else error->all(FLERR,"Illegal atom_modify command");
863 }
864 }
865
866 /* ----------------------------------------------------------------------
867 check that atom IDs are valid
868 error if any atom ID < 0 or atom ID = MAXTAGINT
869 if any atom ID > 0, error if any atom ID == 0
870 if any atom ID > 0, error if tag_enable = 0
871 if all atom IDs = 0, tag_enable must be 0
872 if max atom IDs < natoms, must be duplicates
873 OK if max atom IDs > natoms
874 NOTE: not fully checking that atom IDs are unique
875 ------------------------------------------------------------------------- */
876
tag_check()877 void Atom::tag_check()
878 {
879 tagint min = MAXTAGINT;
880 tagint max = 0;
881
882 for (int i = 0; i < nlocal; i++) {
883 min = MIN(min,tag[i]);
884 max = MAX(max,tag[i]);
885 }
886
887 tagint minall,maxall;
888 MPI_Allreduce(&min,&minall,1,MPI_LMP_TAGINT,MPI_MIN,world);
889 MPI_Allreduce(&max,&maxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
890
891 if (minall < 0) error->all(FLERR,"One or more Atom IDs is negative");
892 if (maxall >= MAXTAGINT) error->all(FLERR,"One or more atom IDs is too big");
893 if (maxall > 0 && minall == 0)
894 error->all(FLERR,"One or more atom IDs is zero");
895 if (maxall > 0 && tag_enable == 0)
896 error->all(FLERR,"Non-zero atom IDs with atom_modify id = no");
897 if (maxall == 0 && natoms && tag_enable)
898 error->all(FLERR,"All atom IDs = 0 but atom_modify id = yes");
899 if (tag_enable && maxall < natoms)
900 error->all(FLERR,"Duplicate atom IDs exist");
901 }
902
903 /* ----------------------------------------------------------------------
904 add unique tags to any atoms with tag = 0
905 new tags are grouped by proc and start after max current tag
906 called after creating new atoms
907 error if new tags will exceed MAXTAGINT
908 ------------------------------------------------------------------------- */
909
tag_extend()910 void Atom::tag_extend()
911 {
912 // maxtag_all = max tag for all atoms
913
914 tagint maxtag = 0;
915 for (int i = 0; i < nlocal; i++) maxtag = MAX(maxtag,tag[i]);
916 tagint maxtag_all;
917 MPI_Allreduce(&maxtag,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
918
919 // DEBUG: useful for generating 64-bit IDs even for small systems
920 // use only when LAMMPS is compiled with BIGBIG
921
922 //maxtag_all += 1000000000000;
923
924 // notag = # of atoms I own with no tag (tag = 0)
925 // notag_sum = # of total atoms on procs <= me with no tag
926
927 bigint notag = 0;
928 for (int i = 0; i < nlocal; i++) if (tag[i] == 0) notag++;
929
930 bigint notag_total;
931 MPI_Allreduce(¬ag,¬ag_total,1,MPI_LMP_BIGINT,MPI_SUM,world);
932 if (notag_total >= MAXTAGINT)
933 error->all(FLERR,"New atom IDs exceed maximum allowed ID");
934
935 bigint notag_sum;
936 MPI_Scan(¬ag,¬ag_sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
937
938 // itag = 1st new tag that my untagged atoms should use
939
940 tagint itag = maxtag_all + notag_sum - notag + 1;
941 for (int i = 0; i < nlocal; i++) if (tag[i] == 0) tag[i] = itag++;
942 }
943
944 /* ----------------------------------------------------------------------
945 check that atom IDs span range from 1 to Natoms inclusive
946 return 0 if mintag != 1 or maxtag != Natoms
947 return 1 if OK
948 doesn't actually check if all tag values are used
949 ------------------------------------------------------------------------- */
950
tag_consecutive()951 int Atom::tag_consecutive()
952 {
953 tagint idmin = MAXTAGINT;
954 tagint idmax = 0;
955
956 for (int i = 0; i < nlocal; i++) {
957 idmin = MIN(idmin,tag[i]);
958 idmax = MAX(idmax,tag[i]);
959 }
960 tagint idminall,idmaxall;
961 MPI_Allreduce(&idmin,&idminall,1,MPI_LMP_TAGINT,MPI_MIN,world);
962 MPI_Allreduce(&idmax,&idmaxall,1,MPI_LMP_TAGINT,MPI_MAX,world);
963
964 if (idminall != 1 || idmaxall != natoms) return 0;
965 return 1;
966 }
967
968 /* ----------------------------------------------------------------------
969 check that bonus data settings are valid
970 error if number of atoms with ellipsoid/line/tri/body flags
971 are consistent with global setting.
972 ------------------------------------------------------------------------- */
973
bonus_check()974 void Atom::bonus_check()
975 {
976 bigint local_ellipsoids = 0, local_lines = 0, local_tris = 0;
977 bigint local_bodies = 0, num_global;
978
979 for (int i = 0; i < nlocal; ++i) {
980 if (ellipsoid && (ellipsoid[i] >=0)) ++local_ellipsoids;
981 if (line && (line[i] >=0)) ++local_lines;
982 if (tri && (tri[i] >=0)) ++local_tris;
983 if (body && (body[i] >=0)) ++local_bodies;
984 }
985
986 MPI_Allreduce(&local_ellipsoids,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
987 if (nellipsoids != num_global)
988 error->all(FLERR,"Inconsistent 'ellipsoids' header value and number of "
989 "atoms with enabled ellipsoid flags");
990
991 MPI_Allreduce(&local_lines,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
992 if (nlines != num_global)
993 error->all(FLERR,"Inconsistent 'lines' header value and number of "
994 "atoms with enabled line flags");
995
996 MPI_Allreduce(&local_tris,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
997 if (ntris != num_global)
998 error->all(FLERR,"Inconsistent 'tris' header value and number of "
999 "atoms with enabled tri flags");
1000
1001 MPI_Allreduce(&local_bodies,&num_global,1,MPI_LMP_BIGINT,MPI_SUM,world);
1002 if (nbodies != num_global)
1003 error->all(FLERR,"Inconsistent 'bodies' header value and number of "
1004 "atoms with enabled body flags");
1005 }
1006
1007 /* ----------------------------------------------------------------------
1008 deallocate molecular topology arrays
1009 done before realloc with (possibly) new 2nd dimension set to
1010 correctly initialized per-atom values, e.g. bond_per_atom
1011 needs to be called whenever 2nd dimensions are changed
1012 and these arrays are already pre-allocated,
1013 e.g. due to grow(1) in create_avec()
1014 ------------------------------------------------------------------------- */
1015
deallocate_topology()1016 void Atom::deallocate_topology()
1017 {
1018 memory->destroy(atom->bond_type);
1019 memory->destroy(atom->bond_atom);
1020 atom->bond_type = nullptr;
1021 atom->bond_atom = nullptr;
1022
1023 memory->destroy(atom->angle_type);
1024 memory->destroy(atom->angle_atom1);
1025 memory->destroy(atom->angle_atom2);
1026 memory->destroy(atom->angle_atom3);
1027 atom->angle_type = nullptr;
1028 atom->angle_atom1 = atom->angle_atom2 = atom->angle_atom3 = nullptr;
1029
1030 memory->destroy(atom->dihedral_type);
1031 memory->destroy(atom->dihedral_atom1);
1032 memory->destroy(atom->dihedral_atom2);
1033 memory->destroy(atom->dihedral_atom3);
1034 memory->destroy(atom->dihedral_atom4);
1035 atom->dihedral_type = nullptr;
1036 atom->dihedral_atom1 = atom->dihedral_atom2 =
1037 atom->dihedral_atom3 = atom->dihedral_atom4 = nullptr;
1038
1039 memory->destroy(atom->improper_type);
1040 memory->destroy(atom->improper_atom1);
1041 memory->destroy(atom->improper_atom2);
1042 memory->destroy(atom->improper_atom3);
1043 memory->destroy(atom->improper_atom4);
1044 atom->improper_type = nullptr;
1045 atom->improper_atom1 = atom->improper_atom2 =
1046 atom->improper_atom3 = atom->improper_atom4 = nullptr;
1047 }
1048
1049 /* ----------------------------------------------------------------------
1050 unpack N lines from Atom section of data file
1051 call style-specific routine to parse line
1052 ------------------------------------------------------------------------- */
1053
data_atoms(int n,char * buf,tagint id_offset,tagint mol_offset,int type_offset,int shiftflag,double * shift)1054 void Atom::data_atoms(int n, char *buf, tagint id_offset, tagint mol_offset,
1055 int type_offset, int shiftflag, double *shift)
1056 {
1057 int m,xptr,iptr;
1058 imageint imagedata;
1059 double xdata[3],lamda[3];
1060 double *coord;
1061 char *next;
1062
1063 next = strchr(buf,'\n');
1064 *next = '\0';
1065 int nwords = utils::trim_and_count_words(buf);
1066 *next = '\n';
1067
1068 if (nwords != avec->size_data_atom && nwords != avec->size_data_atom + 3)
1069 error->all(FLERR,"Incorrect atom format in data file");
1070
1071 char **values = new char*[nwords];
1072
1073 // set bounds for my proc
1074 // if periodic and I am lo/hi proc, adjust bounds by EPSILON
1075 // insures all data atoms will be owned even with round-off
1076
1077 int triclinic = domain->triclinic;
1078
1079 double epsilon[3];
1080 if (triclinic) epsilon[0] = epsilon[1] = epsilon[2] = EPSILON;
1081 else {
1082 epsilon[0] = domain->prd[0] * EPSILON;
1083 epsilon[1] = domain->prd[1] * EPSILON;
1084 epsilon[2] = domain->prd[2] * EPSILON;
1085 }
1086
1087 double sublo[3],subhi[3];
1088 if (triclinic == 0) {
1089 sublo[0] = domain->sublo[0]; subhi[0] = domain->subhi[0];
1090 sublo[1] = domain->sublo[1]; subhi[1] = domain->subhi[1];
1091 sublo[2] = domain->sublo[2]; subhi[2] = domain->subhi[2];
1092 } else {
1093 sublo[0] = domain->sublo_lamda[0]; subhi[0] = domain->subhi_lamda[0];
1094 sublo[1] = domain->sublo_lamda[1]; subhi[1] = domain->subhi_lamda[1];
1095 sublo[2] = domain->sublo_lamda[2]; subhi[2] = domain->subhi_lamda[2];
1096 }
1097
1098 if (comm->layout != Comm::LAYOUT_TILED) {
1099 if (domain->xperiodic) {
1100 if (comm->myloc[0] == 0) sublo[0] -= epsilon[0];
1101 if (comm->myloc[0] == comm->procgrid[0]-1) subhi[0] += epsilon[0];
1102 }
1103 if (domain->yperiodic) {
1104 if (comm->myloc[1] == 0) sublo[1] -= epsilon[1];
1105 if (comm->myloc[1] == comm->procgrid[1]-1) subhi[1] += epsilon[1];
1106 }
1107 if (domain->zperiodic) {
1108 if (comm->myloc[2] == 0) sublo[2] -= epsilon[2];
1109 if (comm->myloc[2] == comm->procgrid[2]-1) subhi[2] += epsilon[2];
1110 }
1111
1112 } else {
1113 if (domain->xperiodic) {
1114 if (comm->mysplit[0][0] == 0.0) sublo[0] -= epsilon[0];
1115 if (comm->mysplit[0][1] == 1.0) subhi[0] += epsilon[0];
1116 }
1117 if (domain->yperiodic) {
1118 if (comm->mysplit[1][0] == 0.0) sublo[1] -= epsilon[1];
1119 if (comm->mysplit[1][1] == 1.0) subhi[1] += epsilon[1];
1120 }
1121 if (domain->zperiodic) {
1122 if (comm->mysplit[2][0] == 0.0) sublo[2] -= epsilon[2];
1123 if (comm->mysplit[2][1] == 1.0) subhi[2] += epsilon[2];
1124 }
1125 }
1126
1127 // xptr = which word in line starts xyz coords
1128 // iptr = which word in line starts ix,iy,iz image flags
1129
1130 xptr = avec->xcol_data - 1;
1131 int imageflag = 0;
1132 if (nwords > avec->size_data_atom) imageflag = 1;
1133 if (imageflag) iptr = nwords - 3;
1134
1135 // loop over lines of atom data
1136 // tokenize the line into values
1137 // extract xyz coords and image flags
1138 // remap atom into simulation box
1139 // if atom is in my sub-domain, unpack its values
1140
1141 int flagx = 0, flagy = 0, flagz = 0;
1142 for (int i = 0; i < n; i++) {
1143 next = strchr(buf,'\n');
1144
1145 for (m = 0; m < nwords; m++) {
1146 buf += strspn(buf," \t\n\r\f");
1147 buf[strcspn(buf," \t\n\r\f")] = '\0';
1148 if (strlen(buf) == 0)
1149 error->all(FLERR,"Incorrect atom format in data file");
1150 values[m] = buf;
1151 buf += strlen(buf)+1;
1152 }
1153
1154 int imx = 0, imy = 0, imz = 0;
1155 if (imageflag) {
1156 imx = utils::inumeric(FLERR,values[iptr],false,lmp);
1157 imy = utils::inumeric(FLERR,values[iptr+1],false,lmp);
1158 imz = utils::inumeric(FLERR,values[iptr+2],false,lmp);
1159 if ((domain->dimension == 2) && (imz != 0))
1160 error->all(FLERR,"Z-direction image flag must be 0 for 2d-systems");
1161 if ((!domain->xperiodic) && (imx != 0)) { flagx = 1; imx = 0; }
1162 if ((!domain->yperiodic) && (imy != 0)) { flagy = 1; imy = 0; }
1163 if ((!domain->zperiodic) && (imz != 0)) { flagz = 1; imz = 0; }
1164 }
1165 imagedata = ((imageint) (imx + IMGMAX) & IMGMASK) |
1166 (((imageint) (imy + IMGMAX) & IMGMASK) << IMGBITS) |
1167 (((imageint) (imz + IMGMAX) & IMGMASK) << IMG2BITS);
1168
1169 xdata[0] = utils::numeric(FLERR,values[xptr],false,lmp);
1170 xdata[1] = utils::numeric(FLERR,values[xptr+1],false,lmp);
1171 xdata[2] = utils::numeric(FLERR,values[xptr+2],false,lmp);
1172 if (shiftflag) {
1173 xdata[0] += shift[0];
1174 xdata[1] += shift[1];
1175 xdata[2] += shift[2];
1176 }
1177
1178 domain->remap(xdata,imagedata);
1179 if (triclinic) {
1180 domain->x2lamda(xdata,lamda);
1181 coord = lamda;
1182 } else coord = xdata;
1183
1184 if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
1185 coord[1] >= sublo[1] && coord[1] < subhi[1] &&
1186 coord[2] >= sublo[2] && coord[2] < subhi[2]) {
1187 avec->data_atom(xdata,imagedata,values);
1188 if (id_offset) tag[nlocal-1] += id_offset;
1189 if (mol_offset) molecule[nlocal-1] += mol_offset;
1190 if (type_offset) {
1191 type[nlocal-1] += type_offset;
1192 if (type[nlocal-1] > ntypes)
1193 error->one(FLERR,"Invalid atom type in Atoms section of data file");
1194 }
1195 }
1196
1197 buf = next + 1;
1198 }
1199
1200 // warn if reading data with non-zero image flags for non-periodic boundaries.
1201 // we may want to turn this into an error at some point, since this essentially
1202 // creates invalid position information that works by accident most of the time.
1203
1204 if (comm->me == 0) {
1205 if (flagx)
1206 error->warning(FLERR,"Non-zero imageflag(s) in x direction for "
1207 "non-periodic boundary reset to zero");
1208 if (flagy)
1209 error->warning(FLERR,"Non-zero imageflag(s) in y direction for "
1210 "non-periodic boundary reset to zero");
1211 if (flagz)
1212 error->warning(FLERR,"Non-zero imageflag(s) in z direction for "
1213 "non-periodic boundary reset to zero");
1214 }
1215
1216 delete [] values;
1217 }
1218
1219 /* ----------------------------------------------------------------------
1220 unpack N lines from Velocity section of data file
1221 check that atom IDs are > 0 and <= map_tag_max
1222 call style-specific routine to parse line
1223 ------------------------------------------------------------------------- */
1224
data_vels(int n,char * buf,tagint id_offset)1225 void Atom::data_vels(int n, char *buf, tagint id_offset)
1226 {
1227 int j,m;
1228 tagint tagdata;
1229 char *next;
1230
1231 next = strchr(buf,'\n');
1232 *next = '\0';
1233 int nwords = utils::trim_and_count_words(buf);
1234 *next = '\n';
1235
1236 if (nwords != avec->size_data_vel)
1237 error->all(FLERR,"Incorrect velocity format in data file");
1238
1239 char **values = new char*[nwords];
1240
1241 // loop over lines of atom velocities
1242 // tokenize the line into values
1243 // if I own atom tag, unpack its values
1244
1245 for (int i = 0; i < n; i++) {
1246 next = strchr(buf,'\n');
1247
1248 for (j = 0; j < nwords; j++) {
1249 buf += strspn(buf," \t\n\r\f");
1250 buf[strcspn(buf," \t\n\r\f")] = '\0';
1251 values[j] = buf;
1252 buf += strlen(buf)+1;
1253 }
1254
1255 tagdata = ATOTAGINT(values[0]) + id_offset;
1256 if (tagdata <= 0 || tagdata > map_tag_max)
1257 error->one(FLERR,"Invalid atom ID in Velocities section of data file");
1258 if ((m = map(tagdata)) >= 0) avec->data_vel(m,&values[1]);
1259
1260 buf = next + 1;
1261 }
1262
1263 delete [] values;
1264 }
1265
1266 /* ----------------------------------------------------------------------
1267 process N bonds read into buf from data files
1268 if count is non-nullptr, just count bonds per atom
1269 else store them with atoms
1270 check that atom IDs are > 0 and <= map_tag_max
1271 ------------------------------------------------------------------------- */
1272
data_bonds(int n,char * buf,int * count,tagint id_offset,int type_offset)1273 void Atom::data_bonds(int n, char *buf, int *count, tagint id_offset,
1274 int type_offset)
1275 {
1276 int m,tmp,itype,rv;
1277 tagint atom1,atom2;
1278 char *next;
1279 int newton_bond = force->newton_bond;
1280
1281 for (int i = 0; i < n; i++) {
1282 next = strchr(buf,'\n');
1283 *next = '\0';
1284 rv = sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT,
1285 &tmp,&itype,&atom1,&atom2);
1286 if (rv != 4)
1287 error->one(FLERR,"Incorrect format of Bonds section in data file");
1288 if (id_offset) {
1289 atom1 += id_offset;
1290 atom2 += id_offset;
1291 }
1292 itype += type_offset;
1293
1294 if ((atom1 <= 0) || (atom1 > map_tag_max) ||
1295 (atom2 <= 0) || (atom2 > map_tag_max) || (atom1 == atom2))
1296 error->one(FLERR,"Invalid atom ID in Bonds section of data file");
1297 if (itype <= 0 || itype > nbondtypes)
1298 error->one(FLERR,"Invalid bond type in Bonds section of data file");
1299 if ((m = map(atom1)) >= 0) {
1300 if (count) count[m]++;
1301 else {
1302 bond_type[m][num_bond[m]] = itype;
1303 bond_atom[m][num_bond[m]] = atom2;
1304 num_bond[m]++;
1305 avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset);
1306 }
1307 }
1308 if (newton_bond == 0) {
1309 if ((m = map(atom2)) >= 0) {
1310 if (count) count[m]++;
1311 else {
1312 bond_type[m][num_bond[m]] = itype;
1313 bond_atom[m][num_bond[m]] = atom1;
1314 num_bond[m]++;
1315 avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset);
1316 }
1317 }
1318 }
1319 buf = next + 1;
1320 }
1321 }
1322
1323 /* ----------------------------------------------------------------------
1324 process N angles read into buf from data files
1325 if count is non-nullptr, just count angles per atom
1326 else store them with atoms
1327 check that atom IDs are > 0 and <= map_tag_max
1328 ------------------------------------------------------------------------- */
1329
data_angles(int n,char * buf,int * count,tagint id_offset,int type_offset)1330 void Atom::data_angles(int n, char *buf, int *count, tagint id_offset,
1331 int type_offset)
1332 {
1333 int m,tmp,itype,rv;
1334 tagint atom1,atom2,atom3;
1335 char *next;
1336 int newton_bond = force->newton_bond;
1337
1338 for (int i = 0; i < n; i++) {
1339 next = strchr(buf,'\n');
1340 *next = '\0';
1341 rv = sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
1342 &tmp,&itype,&atom1,&atom2,&atom3);
1343 if (rv != 5)
1344 error->one(FLERR,"Incorrect format of Angles section in data file");
1345 if (id_offset) {
1346 atom1 += id_offset;
1347 atom2 += id_offset;
1348 atom3 += id_offset;
1349 }
1350 itype += type_offset;
1351
1352 if ((atom1 <= 0) || (atom1 > map_tag_max) ||
1353 (atom2 <= 0) || (atom2 > map_tag_max) ||
1354 (atom3 <= 0) || (atom3 > map_tag_max) ||
1355 (atom1 == atom2) || (atom1 == atom3) || (atom2 == atom3))
1356 error->one(FLERR,"Invalid atom ID in Angles section of data file");
1357 if (itype <= 0 || itype > nangletypes)
1358 error->one(FLERR,"Invalid angle type in Angles section of data file");
1359 if ((m = map(atom2)) >= 0) {
1360 if (count) count[m]++;
1361 else {
1362 angle_type[m][num_angle[m]] = itype;
1363 angle_atom1[m][num_angle[m]] = atom1;
1364 angle_atom2[m][num_angle[m]] = atom2;
1365 angle_atom3[m][num_angle[m]] = atom3;
1366 num_angle[m]++;
1367 }
1368 }
1369 if (newton_bond == 0) {
1370 if ((m = map(atom1)) >= 0) {
1371 if (count) count[m]++;
1372 else {
1373 angle_type[m][num_angle[m]] = itype;
1374 angle_atom1[m][num_angle[m]] = atom1;
1375 angle_atom2[m][num_angle[m]] = atom2;
1376 angle_atom3[m][num_angle[m]] = atom3;
1377 num_angle[m]++;
1378 }
1379 }
1380 if ((m = map(atom3)) >= 0) {
1381 if (count) count[m]++;
1382 else {
1383 angle_type[m][num_angle[m]] = itype;
1384 angle_atom1[m][num_angle[m]] = atom1;
1385 angle_atom2[m][num_angle[m]] = atom2;
1386 angle_atom3[m][num_angle[m]] = atom3;
1387 num_angle[m]++;
1388 }
1389 }
1390 }
1391 buf = next + 1;
1392 }
1393 }
1394
1395 /* ----------------------------------------------------------------------
1396 process N dihedrals read into buf from data files
1397 if count is non-nullptr, just count diihedrals per atom
1398 else store them with atoms
1399 check that atom IDs are > 0 and <= map_tag_max
1400 ------------------------------------------------------------------------- */
1401
data_dihedrals(int n,char * buf,int * count,tagint id_offset,int type_offset)1402 void Atom::data_dihedrals(int n, char *buf, int *count, tagint id_offset,
1403 int type_offset)
1404 {
1405 int m,tmp,itype,rv;
1406 tagint atom1,atom2,atom3,atom4;
1407 char *next;
1408 int newton_bond = force->newton_bond;
1409
1410 for (int i = 0; i < n; i++) {
1411 next = strchr(buf,'\n');
1412 *next = '\0';
1413 rv = sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT
1414 " " TAGINT_FORMAT " " TAGINT_FORMAT,
1415 &tmp,&itype,&atom1,&atom2,&atom3,&atom4);
1416 if (rv != 6)
1417 error->one(FLERR,"Incorrect format of Dihedrals section in data file");
1418 if (id_offset) {
1419 atom1 += id_offset;
1420 atom2 += id_offset;
1421 atom3 += id_offset;
1422 atom4 += id_offset;
1423 }
1424 itype += type_offset;
1425
1426 if ((atom1 <= 0) || (atom1 > map_tag_max) ||
1427 (atom2 <= 0) || (atom2 > map_tag_max) ||
1428 (atom3 <= 0) || (atom3 > map_tag_max) ||
1429 (atom4 <= 0) || (atom4 > map_tag_max) ||
1430 (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
1431 (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
1432 error->one(FLERR,"Invalid atom ID in Dihedrals section of data file");
1433 if (itype <= 0 || itype > ndihedraltypes)
1434 error->one(FLERR,
1435 "Invalid dihedral type in Dihedrals section of data file");
1436 if ((m = map(atom2)) >= 0) {
1437 if (count) count[m]++;
1438 else {
1439 dihedral_type[m][num_dihedral[m]] = itype;
1440 dihedral_atom1[m][num_dihedral[m]] = atom1;
1441 dihedral_atom2[m][num_dihedral[m]] = atom2;
1442 dihedral_atom3[m][num_dihedral[m]] = atom3;
1443 dihedral_atom4[m][num_dihedral[m]] = atom4;
1444 num_dihedral[m]++;
1445 }
1446 }
1447 if (newton_bond == 0) {
1448 if ((m = map(atom1)) >= 0) {
1449 if (count) count[m]++;
1450 else {
1451 dihedral_type[m][num_dihedral[m]] = itype;
1452 dihedral_atom1[m][num_dihedral[m]] = atom1;
1453 dihedral_atom2[m][num_dihedral[m]] = atom2;
1454 dihedral_atom3[m][num_dihedral[m]] = atom3;
1455 dihedral_atom4[m][num_dihedral[m]] = atom4;
1456 num_dihedral[m]++;
1457 }
1458 }
1459 if ((m = map(atom3)) >= 0) {
1460 if (count) count[m]++;
1461 else {
1462 dihedral_type[m][num_dihedral[m]] = itype;
1463 dihedral_atom1[m][num_dihedral[m]] = atom1;
1464 dihedral_atom2[m][num_dihedral[m]] = atom2;
1465 dihedral_atom3[m][num_dihedral[m]] = atom3;
1466 dihedral_atom4[m][num_dihedral[m]] = atom4;
1467 num_dihedral[m]++;
1468 }
1469 }
1470 if ((m = map(atom4)) >= 0) {
1471 if (count) count[m]++;
1472 else {
1473 dihedral_type[m][num_dihedral[m]] = itype;
1474 dihedral_atom1[m][num_dihedral[m]] = atom1;
1475 dihedral_atom2[m][num_dihedral[m]] = atom2;
1476 dihedral_atom3[m][num_dihedral[m]] = atom3;
1477 dihedral_atom4[m][num_dihedral[m]] = atom4;
1478 num_dihedral[m]++;
1479 }
1480 }
1481 }
1482 buf = next + 1;
1483 }
1484 }
1485
1486 /* ----------------------------------------------------------------------
1487 process N impropers read into buf from data files
1488 if count is non-nullptr, just count impropers per atom
1489 else store them with atoms
1490 check that atom IDs are > 0 and <= map_tag_max
1491 ------------------------------------------------------------------------- */
1492
data_impropers(int n,char * buf,int * count,tagint id_offset,int type_offset)1493 void Atom::data_impropers(int n, char *buf, int *count, tagint id_offset,
1494 int type_offset)
1495 {
1496 int m,tmp,itype,rv;
1497 tagint atom1,atom2,atom3,atom4;
1498 char *next;
1499 int newton_bond = force->newton_bond;
1500
1501 for (int i = 0; i < n; i++) {
1502 next = strchr(buf,'\n');
1503 *next = '\0';
1504 rv = sscanf(buf,"%d %d "
1505 TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT,
1506 &tmp,&itype,&atom1,&atom2,&atom3,&atom4);
1507 if (rv != 6)
1508 error->one(FLERR,"Incorrect format of Impropers section in data file");
1509 if (id_offset) {
1510 atom1 += id_offset;
1511 atom2 += id_offset;
1512 atom3 += id_offset;
1513 atom4 += id_offset;
1514 }
1515 itype += type_offset;
1516
1517 if ((atom1 <= 0) || (atom1 > map_tag_max) ||
1518 (atom2 <= 0) || (atom2 > map_tag_max) ||
1519 (atom3 <= 0) || (atom3 > map_tag_max) ||
1520 (atom4 <= 0) || (atom4 > map_tag_max) ||
1521 (atom1 == atom2) || (atom1 == atom3) || (atom1 == atom4) ||
1522 (atom2 == atom3) || (atom2 == atom4) || (atom3 == atom4))
1523 error->one(FLERR,"Invalid atom ID in Impropers section of data file");
1524 if (itype <= 0 || itype > nimpropertypes)
1525 error->one(FLERR,
1526 "Invalid improper type in Impropers section of data file");
1527 if ((m = map(atom2)) >= 0) {
1528 if (count) count[m]++;
1529 else {
1530 improper_type[m][num_improper[m]] = itype;
1531 improper_atom1[m][num_improper[m]] = atom1;
1532 improper_atom2[m][num_improper[m]] = atom2;
1533 improper_atom3[m][num_improper[m]] = atom3;
1534 improper_atom4[m][num_improper[m]] = atom4;
1535 num_improper[m]++;
1536 }
1537 }
1538 if (newton_bond == 0) {
1539 if ((m = map(atom1)) >= 0) {
1540 if (count) count[m]++;
1541 else {
1542 improper_type[m][num_improper[m]] = itype;
1543 improper_atom1[m][num_improper[m]] = atom1;
1544 improper_atom2[m][num_improper[m]] = atom2;
1545 improper_atom3[m][num_improper[m]] = atom3;
1546 improper_atom4[m][num_improper[m]] = atom4;
1547 num_improper[m]++;
1548 }
1549 }
1550 if ((m = map(atom3)) >= 0) {
1551 if (count) count[m]++;
1552 else {
1553 improper_type[m][num_improper[m]] = itype;
1554 improper_atom1[m][num_improper[m]] = atom1;
1555 improper_atom2[m][num_improper[m]] = atom2;
1556 improper_atom3[m][num_improper[m]] = atom3;
1557 improper_atom4[m][num_improper[m]] = atom4;
1558 num_improper[m]++;
1559 }
1560 }
1561 if ((m = map(atom4)) >= 0) {
1562 if (count) count[m]++;
1563 else {
1564 improper_type[m][num_improper[m]] = itype;
1565 improper_atom1[m][num_improper[m]] = atom1;
1566 improper_atom2[m][num_improper[m]] = atom2;
1567 improper_atom3[m][num_improper[m]] = atom3;
1568 improper_atom4[m][num_improper[m]] = atom4;
1569 num_improper[m]++;
1570 }
1571 }
1572 }
1573 buf = next + 1;
1574 }
1575 }
1576
1577 /* ----------------------------------------------------------------------
1578 unpack N lines from atom-style specific bonus section of data file
1579 check that atom IDs are > 0 and <= map_tag_max
1580 call style-specific routine to parse line
1581 ------------------------------------------------------------------------- */
1582
data_bonus(int n,char * buf,AtomVec * avec_bonus,tagint id_offset)1583 void Atom::data_bonus(int n, char *buf, AtomVec *avec_bonus, tagint id_offset)
1584 {
1585 int j,m,tagdata;
1586 char *next;
1587
1588 next = strchr(buf,'\n');
1589 *next = '\0';
1590 int nwords = utils::trim_and_count_words(buf);
1591 *next = '\n';
1592
1593 if (nwords != avec_bonus->size_data_bonus)
1594 error->all(FLERR,"Incorrect bonus data format in data file");
1595
1596 char **values = new char*[nwords];
1597
1598 // loop over lines of bonus atom data
1599 // tokenize the line into values
1600 // if I own atom tag, unpack its values
1601
1602 for (int i = 0; i < n; i++) {
1603 next = strchr(buf,'\n');
1604
1605 for (j = 0; j < nwords; j++) {
1606 buf += strspn(buf," \t\n\r\f");
1607 buf[strcspn(buf," \t\n\r\f")] = '\0';
1608 values[j] = buf;
1609 buf += strlen(buf)+1;
1610 }
1611
1612 tagdata = ATOTAGINT(values[0]) + id_offset;
1613 if (tagdata <= 0 || tagdata > map_tag_max)
1614 error->one(FLERR,"Invalid atom ID in Bonus section of data file");
1615
1616 // ok to call child's data_atom_bonus() method thru parent avec_bonus,
1617 // since data_bonus() was called with child ptr, and method is virtual
1618
1619 if ((m = map(tagdata)) >= 0) avec_bonus->data_atom_bonus(m,&values[1]);
1620
1621 buf = next + 1;
1622 }
1623
1624 delete [] values;
1625 }
1626
1627 /* ----------------------------------------------------------------------
1628 unpack N bodies from Bodies section of data file
1629 each body spans multiple lines
1630 check that atom IDs are > 0 and <= map_tag_max
1631 call style-specific routine to parse line
1632 ------------------------------------------------------------------------- */
1633
data_bodies(int n,char * buf,AtomVec * avec_body,tagint id_offset)1634 void Atom::data_bodies(int n, char *buf, AtomVec *avec_body, tagint id_offset)
1635 {
1636 int j,m,nvalues,tagdata,ninteger,ndouble;
1637
1638 int maxint = 0;
1639 int maxdouble = 0;
1640 int *ivalues = nullptr;
1641 double *dvalues = nullptr;
1642
1643 if (!unique_tags) unique_tags = new std::set<tagint>;
1644
1645 // loop over lines of body data
1646 // if I own atom tag, tokenize lines into ivalues/dvalues, call data_body()
1647 // else skip values
1648
1649 for (int i = 0; i < n; i++) {
1650 buf += strspn(buf," \t\n\r\f");
1651 buf[strcspn(buf," \t\n\r\f")] = '\0';
1652 tagdata = utils::tnumeric(FLERR,buf,false,lmp) + id_offset;
1653 buf += strlen(buf)+1;
1654
1655 if (tagdata <= 0 || tagdata > map_tag_max)
1656 error->one(FLERR,"Invalid atom ID in Bodies section of data file");
1657
1658 if (unique_tags->find(tagdata) == unique_tags->end())
1659 unique_tags->insert(tagdata);
1660 else
1661 error->one(FLERR,"Duplicate atom ID in Bodies section of data file");
1662
1663 buf += strspn(buf," \t\n\r\f");
1664 buf[strcspn(buf," \t\n\r\f")] = '\0';
1665 ninteger = utils::inumeric(FLERR,buf,false,lmp);
1666 buf += strlen(buf)+1;
1667
1668 buf += strspn(buf," \t\n\r\f");
1669 buf[strcspn(buf," \t\n\r\f")] = '\0';
1670 ndouble = utils::inumeric(FLERR,buf,false,lmp);
1671 buf += strlen(buf)+1;
1672
1673 if ((m = map(tagdata)) >= 0) {
1674 if (ninteger > maxint) {
1675 delete [] ivalues;
1676 maxint = ninteger;
1677 ivalues = new int[maxint];
1678 }
1679 if (ndouble > maxdouble) {
1680 delete [] dvalues;
1681 maxdouble = ndouble;
1682 dvalues = new double[maxdouble];
1683 }
1684
1685 for (j = 0; j < ninteger; j++) {
1686 buf += strspn(buf," \t\n\r\f");
1687 buf[strcspn(buf," \t\n\r\f")] = '\0';
1688 ivalues[j] = utils::inumeric(FLERR,buf,false,lmp);
1689 buf += strlen(buf)+1;
1690 }
1691
1692 for (j = 0; j < ndouble; j++) {
1693 buf += strspn(buf," \t\n\r\f");
1694 buf[strcspn(buf," \t\n\r\f")] = '\0';
1695 dvalues[j] = utils::numeric(FLERR,buf,false,lmp);
1696 buf += strlen(buf)+1;
1697 }
1698
1699 avec_body->data_body(m,ninteger,ndouble,ivalues,dvalues);
1700
1701 } else {
1702 nvalues = ninteger + ndouble; // number of values to skip
1703 for (j = 0; j < nvalues; j++) {
1704 buf += strspn(buf," \t\n\r\f");
1705 buf[strcspn(buf," \t\n\r\f")] = '\0';
1706 buf += strlen(buf)+1;
1707 }
1708 }
1709 }
1710
1711 delete [] ivalues;
1712 delete [] dvalues;
1713 }
1714
1715 /* ----------------------------------------------------------------------
1716 init per-atom fix/compute/variable values for newly created atoms
1717 called from create_atoms, read_data, read_dump,
1718 lib::lammps_create_atoms()
1719 fixes, computes, variables may or may not exist when called
1720 ------------------------------------------------------------------------- */
1721
data_fix_compute_variable(int nprev,int nnew)1722 void Atom::data_fix_compute_variable(int nprev, int nnew)
1723 {
1724 for (int m = 0; m < modify->nfix; m++) {
1725 Fix *fix = modify->fix[m];
1726 if (fix->create_attribute)
1727 for (int i = nprev; i < nnew; i++)
1728 fix->set_arrays(i);
1729 }
1730
1731 for (int m = 0; m < modify->ncompute; m++) {
1732 Compute *compute = modify->compute[m];
1733 if (compute->create_attribute)
1734 for (int i = nprev; i < nnew; i++)
1735 compute->set_arrays(i);
1736 }
1737
1738 for (int i = nprev; i < nnew; i++)
1739 input->variable->set_arrays(i);
1740 }
1741
1742 /* ----------------------------------------------------------------------
1743 allocate arrays of length ntypes
1744 only done after ntypes is set
1745 ------------------------------------------------------------------------- */
1746
allocate_type_arrays()1747 void Atom::allocate_type_arrays()
1748 {
1749 if (avec->mass_type == AtomVec::PER_TYPE) {
1750 mass = new double[ntypes+1];
1751 mass_setflag = new int[ntypes+1];
1752 for (int itype = 1; itype <= ntypes; itype++) mass_setflag[itype] = 0;
1753 }
1754 }
1755
1756 /* ----------------------------------------------------------------------
1757 set a mass and flag it as set
1758 called from reading of data file
1759 type_offset may be used when reading multiple data files
1760 ------------------------------------------------------------------------- */
1761
set_mass(const char * file,int line,const char * str,int type_offset)1762 void Atom::set_mass(const char *file, int line, const char *str, int type_offset)
1763 {
1764 if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style");
1765
1766 int itype;
1767 double mass_one;
1768 int n = sscanf(str,"%d %lg",&itype,&mass_one);
1769 if (n != 2) error->all(file,line,"Invalid mass line in data file");
1770 itype += type_offset;
1771
1772 if (itype < 1 || itype > ntypes)
1773 error->all(file,line,"Invalid type for mass set");
1774
1775 mass[itype] = mass_one;
1776 mass_setflag[itype] = 1;
1777
1778 if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value");
1779 }
1780
1781 /* ----------------------------------------------------------------------
1782 set a mass and flag it as set
1783 called from EAM pair routine
1784 ------------------------------------------------------------------------- */
1785
set_mass(const char * file,int line,int itype,double value)1786 void Atom::set_mass(const char *file, int line, int itype, double value)
1787 {
1788 if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style");
1789 if (itype < 1 || itype > ntypes)
1790 error->all(file,line,"Invalid type for mass set");
1791
1792 mass[itype] = value;
1793 mass_setflag[itype] = 1;
1794
1795 if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value");
1796 }
1797
1798 /* ----------------------------------------------------------------------
1799 set one or more masses and flag them as set
1800 called from reading of input script
1801 ------------------------------------------------------------------------- */
1802
set_mass(const char * file,int line,int,char ** arg)1803 void Atom::set_mass(const char *file, int line, int /*narg*/, char **arg)
1804 {
1805 if (mass == nullptr) error->all(file,line,"Cannot set mass for this atom style");
1806
1807 int lo,hi;
1808 utils::bounds(file,line,arg[0],1,ntypes,lo,hi,error);
1809 if (lo < 1 || hi > ntypes) error->all(file,line,"Invalid type for mass set");
1810
1811 for (int itype = lo; itype <= hi; itype++) {
1812 mass[itype] = utils::numeric(FLERR,arg[1],false,lmp);
1813 mass_setflag[itype] = 1;
1814
1815 if (mass[itype] <= 0.0) error->all(file,line,"Invalid mass value");
1816 }
1817 }
1818
1819 /* ----------------------------------------------------------------------
1820 set all masses
1821 called from reading of restart file, also from ServerMD
1822 ------------------------------------------------------------------------- */
1823
set_mass(double * values)1824 void Atom::set_mass(double *values)
1825 {
1826 for (int itype = 1; itype <= ntypes; itype++) {
1827 mass[itype] = values[itype];
1828 mass_setflag[itype] = 1;
1829 }
1830 }
1831
1832 /* ----------------------------------------------------------------------
1833 check that all per-atom-type masses have been set
1834 ------------------------------------------------------------------------- */
1835
check_mass(const char * file,int line)1836 void Atom::check_mass(const char *file, int line)
1837 {
1838 if (mass == nullptr) return;
1839 if (rmass_flag) return;
1840 for (int itype = 1; itype <= ntypes; itype++)
1841 if (mass_setflag[itype] == 0)
1842 error->all(file,line,"Not all per-type masses are set");
1843 }
1844
1845 /* ----------------------------------------------------------------------
1846 check that radii of all particles of itype are the same
1847 return 1 if true, else return 0
1848 also return the radius value for that type
1849 ------------------------------------------------------------------------- */
1850
radius_consistency(int itype,double & rad)1851 int Atom::radius_consistency(int itype, double &rad)
1852 {
1853 double value = -1.0;
1854 int flag = 0;
1855 for (int i = 0; i < nlocal; i++) {
1856 if (type[i] != itype) continue;
1857 if (value < 0.0) value = radius[i];
1858 else if (value != radius[i]) flag = 1;
1859 }
1860
1861 int flagall;
1862 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
1863 if (flagall) return 0;
1864
1865 MPI_Allreduce(&value,&rad,1,MPI_DOUBLE,MPI_MAX,world);
1866 return 1;
1867 }
1868
1869 /* ----------------------------------------------------------------------
1870 check that shape of all particles of itype are the same
1871 return 1 if true, else return 0
1872 also return the 3 shape params for itype
1873 ------------------------------------------------------------------------- */
1874
shape_consistency(int itype,double & shapex,double & shapey,double & shapez)1875 int Atom::shape_consistency(int itype,
1876 double &shapex, double &shapey, double &shapez)
1877 {
1878 double zero[3] = {0.0, 0.0, 0.0};
1879 double one[3] = {-1.0, -1.0, -1.0};
1880 double *shape;
1881
1882 auto avec_ellipsoid = (AtomVecEllipsoid *) style_match("ellipsoid");
1883 auto bonus = avec_ellipsoid->bonus;
1884
1885 int flag = 0;
1886 for (int i = 0; i < nlocal; i++) {
1887 if (type[i] != itype) continue;
1888 if (ellipsoid[i] < 0) shape = zero;
1889 else shape = bonus[ellipsoid[i]].shape;
1890
1891 if (one[0] < 0.0) {
1892 one[0] = shape[0];
1893 one[1] = shape[1];
1894 one[2] = shape[2];
1895 } else if (one[0] != shape[0] || one[1] != shape[1] || one[2] != shape[2])
1896 flag = 1;
1897 }
1898
1899 int flagall;
1900 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
1901 if (flagall) return 0;
1902
1903 double oneall[3];
1904 MPI_Allreduce(one,oneall,3,MPI_DOUBLE,MPI_MAX,world);
1905 shapex = oneall[0];
1906 shapey = oneall[1];
1907 shapez = oneall[2];
1908 return 1;
1909 }
1910
1911 /* ----------------------------------------------------------------------
1912 add a new molecule template = set of molecules
1913 ------------------------------------------------------------------------- */
1914
add_molecule(int narg,char ** arg)1915 void Atom::add_molecule(int narg, char **arg)
1916 {
1917 if (narg < 1) error->all(FLERR,"Illegal molecule command");
1918
1919 if (find_molecule(arg[0]) >= 0)
1920 error->all(FLERR,"Reuse of molecule template ID");
1921
1922 // 1st molecule in set stores nset = # of mols, others store nset = 0
1923 // ifile = count of molecules in set
1924 // index = argument index where next molecule starts, updated by constructor
1925
1926 int ifile = 1;
1927 int index = 1;
1928 while (1) {
1929 molecules = (Molecule **)
1930 memory->srealloc(molecules,(nmolecule+1)*sizeof(Molecule *),
1931 "atom::molecules");
1932 molecules[nmolecule] = new Molecule(lmp,narg,arg,index);
1933 molecules[nmolecule]->nset = 0;
1934 molecules[nmolecule-ifile+1]->nset++;
1935 nmolecule++;
1936 if (molecules[nmolecule-1]->last) break;
1937 ifile++;
1938 }
1939 }
1940
1941 /* ----------------------------------------------------------------------
1942 find first molecule in set with template ID
1943 return -1 if does not exist
1944 ------------------------------------------------------------------------- */
1945
find_molecule(char * id)1946 int Atom::find_molecule(char *id)
1947 {
1948 if (id == nullptr) return -1;
1949 int imol;
1950 for (imol = 0; imol < nmolecule; imol++)
1951 if (strcmp(id,molecules[imol]->id) == 0) return imol;
1952 return -1;
1953 }
1954
1955 /* ----------------------------------------------------------------------
1956 add info to current atom ilocal from molecule template onemol and its iatom
1957 offset = atom ID preceding IDs of atoms in this molecule
1958 called by fixes and commands that add molecules
1959 ------------------------------------------------------------------------- */
1960
add_molecule_atom(Molecule * onemol,int iatom,int ilocal,tagint offset)1961 void Atom::add_molecule_atom(Molecule *onemol, int iatom,
1962 int ilocal, tagint offset)
1963 {
1964 if (onemol->qflag && q_flag) q[ilocal] = onemol->q[iatom];
1965 if (onemol->radiusflag && radius_flag) radius[ilocal] = onemol->radius[iatom];
1966 if (onemol->rmassflag && rmass_flag) rmass[ilocal] = onemol->rmass[iatom];
1967 else if (rmass_flag)
1968 rmass[ilocal] = 4.0*MY_PI/3.0 *
1969 radius[ilocal]*radius[ilocal]*radius[ilocal];
1970 if (onemol->bodyflag) {
1971 body[ilocal] = 0; // as if a body read from data file
1972 onemol->avec_body->data_body(ilocal,onemol->nibody,onemol->ndbody,
1973 onemol->ibodyparams,onemol->dbodyparams);
1974 onemol->avec_body->set_quat(ilocal,onemol->quat_external);
1975 }
1976
1977 if (molecular != Atom::MOLECULAR) return;
1978
1979 // add bond topology info
1980 // for molecular atom styles, but not atom style template
1981
1982 if (avec->bonds_allow) {
1983 num_bond[ilocal] = onemol->num_bond[iatom];
1984 for (int i = 0; i < num_bond[ilocal]; i++) {
1985 bond_type[ilocal][i] = onemol->bond_type[iatom][i];
1986 bond_atom[ilocal][i] = onemol->bond_atom[iatom][i] + offset;
1987 }
1988 }
1989
1990 if (avec->angles_allow) {
1991 num_angle[ilocal] = onemol->num_angle[iatom];
1992 for (int i = 0; i < num_angle[ilocal]; i++) {
1993 angle_type[ilocal][i] = onemol->angle_type[iatom][i];
1994 angle_atom1[ilocal][i] = onemol->angle_atom1[iatom][i] + offset;
1995 angle_atom2[ilocal][i] = onemol->angle_atom2[iatom][i] + offset;
1996 angle_atom3[ilocal][i] = onemol->angle_atom3[iatom][i] + offset;
1997 }
1998 }
1999
2000 if (avec->dihedrals_allow) {
2001 num_dihedral[ilocal] = onemol->num_dihedral[iatom];
2002 for (int i = 0; i < num_dihedral[ilocal]; i++) {
2003 dihedral_type[ilocal][i] = onemol->dihedral_type[iatom][i];
2004 dihedral_atom1[ilocal][i] = onemol->dihedral_atom1[iatom][i] + offset;
2005 dihedral_atom2[ilocal][i] = onemol->dihedral_atom2[iatom][i] + offset;
2006 dihedral_atom3[ilocal][i] = onemol->dihedral_atom3[iatom][i] + offset;
2007 dihedral_atom4[ilocal][i] = onemol->dihedral_atom4[iatom][i] + offset;
2008 }
2009 }
2010
2011 if (avec->impropers_allow) {
2012 num_improper[ilocal] = onemol->num_improper[iatom];
2013 for (int i = 0; i < num_improper[ilocal]; i++) {
2014 improper_type[ilocal][i] = onemol->improper_type[iatom][i];
2015 improper_atom1[ilocal][i] = onemol->improper_atom1[iatom][i] + offset;
2016 improper_atom2[ilocal][i] = onemol->improper_atom2[iatom][i] + offset;
2017 improper_atom3[ilocal][i] = onemol->improper_atom3[iatom][i] + offset;
2018 improper_atom4[ilocal][i] = onemol->improper_atom4[iatom][i] + offset;
2019 }
2020 }
2021
2022 if (onemol->specialflag) {
2023 nspecial[ilocal][0] = onemol->nspecial[iatom][0];
2024 nspecial[ilocal][1] = onemol->nspecial[iatom][1];
2025 int n = nspecial[ilocal][2] = onemol->nspecial[iatom][2];
2026 for (int i = 0; i < n; i++)
2027 special[ilocal][i] = onemol->special[iatom][i] + offset;
2028 }
2029 }
2030
2031 /* ----------------------------------------------------------------------
2032 reorder owned atoms so those in firstgroup appear first
2033 called by comm->exchange() if atom_modify first group is set
2034 only owned atoms exist at this point, no ghost atoms
2035 ------------------------------------------------------------------------- */
2036
first_reorder()2037 void Atom::first_reorder()
2038 {
2039 // insure there is one extra atom location at end of arrays for swaps
2040
2041 if (nlocal == nmax) avec->grow(0);
2042
2043 // loop over owned atoms
2044 // nfirst = index of first atom not in firstgroup
2045 // when find firstgroup atom out of place, swap it with atom nfirst
2046
2047 int bitmask = group->bitmask[firstgroup];
2048 nfirst = 0;
2049 while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
2050
2051 for (int i = 0; i < nlocal; i++) {
2052 if (mask[i] & bitmask && i > nfirst) {
2053 avec->copy(i,nlocal,0);
2054 avec->copy(nfirst,i,0);
2055 avec->copy(nlocal,nfirst,0);
2056 while (nfirst < nlocal && mask[nfirst] & bitmask) nfirst++;
2057 }
2058 }
2059 }
2060
2061 /* ----------------------------------------------------------------------
2062 perform spatial sort of atoms within my sub-domain
2063 always called between comm->exchange() and comm->borders()
2064 don't have to worry about clearing/setting atom->map since done in comm
2065 ------------------------------------------------------------------------- */
2066
sort()2067 void Atom::sort()
2068 {
2069 int i,m,n,ix,iy,iz,ibin,empty;
2070
2071 // set next timestep for sorting to take place
2072
2073 nextsort = (update->ntimestep/sortfreq)*sortfreq + sortfreq;
2074
2075 // re-setup sort bins if needed
2076
2077 if (domain->box_change) setup_sort_bins();
2078 if (nbins == 1) return;
2079
2080 // reallocate per-atom vectors if needed
2081
2082 if (nlocal > maxnext) {
2083 memory->destroy(next);
2084 memory->destroy(permute);
2085 maxnext = atom->nmax;
2086 memory->create(next,maxnext,"atom:next");
2087 memory->create(permute,maxnext,"atom:permute");
2088 }
2089
2090 // insure there is one extra atom location at end of arrays for swaps
2091
2092 if (nlocal == nmax) avec->grow(0);
2093
2094 // bin atoms in reverse order so linked list will be in forward order
2095
2096 for (i = 0; i < nbins; i++) binhead[i] = -1;
2097
2098 for (i = nlocal-1; i >= 0; i--) {
2099 ix = static_cast<int> ((x[i][0]-bboxlo[0])*bininvx);
2100 iy = static_cast<int> ((x[i][1]-bboxlo[1])*bininvy);
2101 iz = static_cast<int> ((x[i][2]-bboxlo[2])*bininvz);
2102 ix = MAX(ix,0);
2103 iy = MAX(iy,0);
2104 iz = MAX(iz,0);
2105 ix = MIN(ix,nbinx-1);
2106 iy = MIN(iy,nbiny-1);
2107 iz = MIN(iz,nbinz-1);
2108 ibin = iz*nbiny*nbinx + iy*nbinx + ix;
2109 next[i] = binhead[ibin];
2110 binhead[ibin] = i;
2111 }
2112
2113 // permute = desired permutation of atoms
2114 // permute[I] = J means Ith new atom will be Jth old atom
2115
2116 n = 0;
2117 for (m = 0; m < nbins; m++) {
2118 i = binhead[m];
2119 while (i >= 0) {
2120 permute[n++] = i;
2121 i = next[i];
2122 }
2123 }
2124
2125 // current = current permutation, just reuse next vector
2126 // current[I] = J means Ith current atom is Jth old atom
2127
2128 int *current = next;
2129 for (i = 0; i < nlocal; i++) current[i] = i;
2130
2131 // reorder local atom list, when done, current = permute
2132 // perform "in place" using copy() to extra atom location at end of list
2133 // inner while loop processes one cycle of the permutation
2134 // copy before inner-loop moves an atom to end of atom list
2135 // copy after inner-loop moves atom at end of list back into list
2136 // empty = location in atom list that is currently empty
2137
2138 for (i = 0; i < nlocal; i++) {
2139 if (current[i] == permute[i]) continue;
2140 avec->copy(i,nlocal,0);
2141 empty = i;
2142 while (permute[empty] != i) {
2143 avec->copy(permute[empty],empty,0);
2144 empty = current[empty] = permute[empty];
2145 }
2146 avec->copy(nlocal,empty,0);
2147 current[empty] = permute[empty];
2148 }
2149
2150 // sanity check that current = permute
2151
2152 //int flag = 0;
2153 //for (i = 0; i < nlocal; i++)
2154 // if (current[i] != permute[i]) flag = 1;
2155 //int flagall;
2156 //MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
2157 //if (flagall) error->all(FLERR,"Atom sort did not operate correctly");
2158 }
2159
2160 /* ----------------------------------------------------------------------
2161 setup bins for spatial sorting of atoms
2162 ------------------------------------------------------------------------- */
2163
setup_sort_bins()2164 void Atom::setup_sort_bins()
2165 {
2166 // binsize:
2167 // user setting if explicitly set
2168 // default = 1/2 of neighbor cutoff
2169 // check if neighbor cutoff = 0.0
2170 // and in that case, disable sorting
2171
2172 double binsize = 0.0;
2173 if (userbinsize > 0.0) binsize = userbinsize;
2174 else if (neighbor->cutneighmax > 0.0) binsize = 0.5 * neighbor->cutneighmax;
2175
2176 if ((binsize == 0.0) && (sortfreq > 0)) {
2177 sortfreq = 0;
2178 if (comm->me == 0)
2179 error->warning(FLERR,"No pairwise cutoff or binsize set. "
2180 "Atom sorting therefore disabled.");
2181 return;
2182 }
2183
2184 double bininv = 1.0/binsize;
2185
2186 // nbin xyz = local bins
2187 // bbox lo/hi = bounding box of my sub-domain
2188
2189 if (domain->triclinic)
2190 domain->bbox(domain->sublo_lamda,domain->subhi_lamda,bboxlo,bboxhi);
2191 else {
2192 bboxlo[0] = domain->sublo[0];
2193 bboxlo[1] = domain->sublo[1];
2194 bboxlo[2] = domain->sublo[2];
2195 bboxhi[0] = domain->subhi[0];
2196 bboxhi[1] = domain->subhi[1];
2197 bboxhi[2] = domain->subhi[2];
2198 }
2199
2200 nbinx = static_cast<int> ((bboxhi[0]-bboxlo[0]) * bininv);
2201 nbiny = static_cast<int> ((bboxhi[1]-bboxlo[1]) * bininv);
2202 nbinz = static_cast<int> ((bboxhi[2]-bboxlo[2]) * bininv);
2203 if (domain->dimension == 2) nbinz = 1;
2204
2205 if (nbinx == 0) nbinx = 1;
2206 if (nbiny == 0) nbiny = 1;
2207 if (nbinz == 0) nbinz = 1;
2208
2209 bininvx = nbinx / (bboxhi[0]-bboxlo[0]);
2210 bininvy = nbiny / (bboxhi[1]-bboxlo[1]);
2211 bininvz = nbinz / (bboxhi[2]-bboxlo[2]);
2212
2213 #ifdef LMP_INTEL
2214 int intel_neigh = 0;
2215 if (neighbor->nrequest) {
2216 if (neighbor->requests[0]->intel) intel_neigh = 1;
2217 } else if (neighbor->old_nrequest)
2218 if (neighbor->old_requests[0]->intel) intel_neigh = 1;
2219 if (intel_neigh && userbinsize == 0.0) {
2220 if (neighbor->binsizeflag) bininv = 1.0/neighbor->binsize_user;
2221
2222 double nx_low = neighbor->bboxlo[0];
2223 double ny_low = neighbor->bboxlo[1];
2224 double nz_low = neighbor->bboxlo[2];
2225 double nxbbox = neighbor->bboxhi[0] - nx_low;
2226 double nybbox = neighbor->bboxhi[1] - ny_low;
2227 double nzbbox = neighbor->bboxhi[2] - nz_low;
2228 int nnbinx = static_cast<int> (nxbbox * bininv);
2229 int nnbiny = static_cast<int> (nybbox * bininv);
2230 int nnbinz = static_cast<int> (nzbbox * bininv);
2231 if (domain->dimension == 2) nnbinz = 1;
2232
2233 if (nnbinx == 0) nnbinx = 1;
2234 if (nnbiny == 0) nnbiny = 1;
2235 if (nnbinz == 0) nnbinz = 1;
2236
2237 double binsizex = nxbbox/nnbinx;
2238 double binsizey = nybbox/nnbiny;
2239 double binsizez = nzbbox/nnbinz;
2240
2241 bininvx = 1.0 / binsizex;
2242 bininvy = 1.0 / binsizey;
2243 bininvz = 1.0 / binsizez;
2244
2245 int lxo = (bboxlo[0] - nx_low) * bininvx;
2246 int lyo = (bboxlo[1] - ny_low) * bininvy;
2247 int lzo = (bboxlo[2] - nz_low) * bininvz;
2248 bboxlo[0] = nx_low + static_cast<double>(lxo) / bininvx;
2249 bboxlo[1] = ny_low + static_cast<double>(lyo) / bininvy;
2250 bboxlo[2] = nz_low + static_cast<double>(lzo) / bininvz;
2251 nbinx = static_cast<int>((bboxhi[0] - bboxlo[0]) * bininvx) + 1;
2252 nbiny = static_cast<int>((bboxhi[1] - bboxlo[1]) * bininvy) + 1;
2253 nbinz = static_cast<int>((bboxhi[2] - bboxlo[2]) * bininvz) + 1;
2254 bboxhi[0] = bboxlo[0] + static_cast<double>(nbinx) / bininvx;
2255 bboxhi[1] = bboxlo[1] + static_cast<double>(nbiny) / bininvy;
2256 bboxhi[2] = bboxlo[2] + static_cast<double>(nbinz) / bininvz;
2257 }
2258 #endif
2259
2260 #ifdef LMP_GPU
2261 if (userbinsize == 0.0) {
2262 int ifix = modify->find_fix("package_gpu");
2263 if (ifix >= 0) {
2264 const double subx = domain->subhi[0] - domain->sublo[0];
2265 const double suby = domain->subhi[1] - domain->sublo[1];
2266 const double subz = domain->subhi[2] - domain->sublo[2];
2267
2268 FixGPU *fix = static_cast<FixGPU *>(modify->fix[ifix]);
2269 binsize = fix->binsize(subx, suby, subz, atom->nlocal,
2270 neighbor->cutneighmax);
2271 bininv = 1.0 / binsize;
2272
2273 nbinx = static_cast<int> (ceil(subx * bininv));
2274 nbiny = static_cast<int> (ceil(suby * bininv));
2275 nbinz = static_cast<int> (ceil(subz * bininv));
2276 if (domain->dimension == 2) nbinz = 1;
2277
2278 if (nbinx == 0) nbinx = 1;
2279 if (nbiny == 0) nbiny = 1;
2280 if (nbinz == 0) nbinz = 1;
2281
2282 bininvx = bininv;
2283 bininvy = bininv;
2284 bininvz = bininv;
2285 }
2286 }
2287 #endif
2288
2289 if (1.0*nbinx*nbiny*nbinz > INT_MAX)
2290 error->one(FLERR,"Too many atom sorting bins");
2291
2292 nbins = nbinx*nbiny*nbinz;
2293
2294 // reallocate per-bin memory if needed
2295
2296 if (nbins > maxbin) {
2297 memory->destroy(binhead);
2298 maxbin = nbins;
2299 memory->create(binhead,maxbin,"atom:binhead");
2300 }
2301 }
2302
2303 /* ----------------------------------------------------------------------
2304 register a callback to a fix so it can manage atom-based arrays
2305 happens when fix is created
2306 flag = 0 for grow, 1 for restart, 2 for border comm
2307 ------------------------------------------------------------------------- */
2308
add_callback(int flag)2309 void Atom::add_callback(int flag)
2310 {
2311 int ifix;
2312
2313 // find the fix
2314 // if find null pointer:
2315 // it's this one, since it is being replaced and has just been deleted
2316 // at this point in re-creation
2317 // if don't find null pointer:
2318 // i is set to nfix = new one currently being added at end of list
2319
2320 for (ifix = 0; ifix < modify->nfix; ifix++)
2321 if (modify->fix[ifix] == nullptr) break;
2322
2323 // add callback to lists and sort, reallocating if necessary
2324 // sorting is required in cases where fixes were replaced as it ensures atom
2325 // data is read/written/transfered in the same order that fixes are called
2326
2327 if (flag == GROW) {
2328 if (nextra_grow == nextra_grow_max) {
2329 nextra_grow_max += DELTA;
2330 memory->grow(extra_grow,nextra_grow_max,"atom:extra_grow");
2331 }
2332 extra_grow[nextra_grow] = ifix;
2333 nextra_grow++;
2334 std::sort(extra_grow, extra_grow + nextra_grow);
2335 } else if (flag == RESTART) {
2336 if (nextra_restart == nextra_restart_max) {
2337 nextra_restart_max += DELTA;
2338 memory->grow(extra_restart,nextra_restart_max,"atom:extra_restart");
2339 }
2340 extra_restart[nextra_restart] = ifix;
2341 nextra_restart++;
2342 std::sort(extra_restart, extra_restart + nextra_restart);
2343 } else if (flag == BORDER) {
2344 if (nextra_border == nextra_border_max) {
2345 nextra_border_max += DELTA;
2346 memory->grow(extra_border,nextra_border_max,"atom:extra_border");
2347 }
2348 extra_border[nextra_border] = ifix;
2349 nextra_border++;
2350 std::sort(extra_border, extra_border + nextra_border);
2351 }
2352 }
2353
2354 /* ----------------------------------------------------------------------
2355 unregister a callback to a fix
2356 happens when fix is deleted, called by its destructor
2357 flag = 0 for grow, 1 for restart
2358 ------------------------------------------------------------------------- */
2359
delete_callback(const char * id,int flag)2360 void Atom::delete_callback(const char *id, int flag)
2361 {
2362 if (id == nullptr) return;
2363
2364 int ifix = modify->find_fix(id);
2365
2366 // compact the list of callbacks
2367
2368 if (flag == GROW) {
2369 int match;
2370 for (match = 0; match < nextra_grow; match++)
2371 if (extra_grow[match] == ifix) break;
2372 if ((nextra_grow == 0) || (match == nextra_grow))
2373 error->all(FLERR,"Trying to delete non-existent Atom::grow() callback");
2374 for (int i = match; i < nextra_grow-1; i++)
2375 extra_grow[i] = extra_grow[i+1];
2376 nextra_grow--;
2377
2378 } else if (flag == RESTART) {
2379 int match;
2380 for (match = 0; match < nextra_restart; match++)
2381 if (extra_restart[match] == ifix) break;
2382 if ((nextra_restart == 0) || (match == nextra_restart))
2383 error->all(FLERR,"Trying to delete non-existent Atom::restart() callback");
2384 for (int i = match; i < nextra_restart-1; i++)
2385 extra_restart[i] = extra_restart[i+1];
2386 nextra_restart--;
2387
2388 } else if (flag == BORDER) {
2389 int match;
2390 for (match = 0; match < nextra_border; match++)
2391 if (extra_border[match] == ifix) break;
2392 if ((nextra_border == 0) || (match == nextra_border))
2393 error->all(FLERR,"Trying to delete non-existent Atom::border() callback");
2394 for (int i = match; i < nextra_border-1; i++)
2395 extra_border[i] = extra_border[i+1];
2396 nextra_border--;
2397 }
2398 }
2399
2400 /* ----------------------------------------------------------------------
2401 decrement ptrs in callback lists to fixes beyond the deleted ifix
2402 happens after fix is deleted
2403 ------------------------------------------------------------------------- */
2404
update_callback(int ifix)2405 void Atom::update_callback(int ifix)
2406 {
2407 for (int i = 0; i < nextra_grow; i++)
2408 if (extra_grow[i] > ifix) extra_grow[i]--;
2409 for (int i = 0; i < nextra_restart; i++)
2410 if (extra_restart[i] > ifix) extra_restart[i]--;
2411 for (int i = 0; i < nextra_border; i++)
2412 if (extra_border[i] > ifix) extra_border[i]--;
2413 }
2414
2415 /* ----------------------------------------------------------------------
2416 find custom per-atom vector with name
2417 return index if found, -1 if not found
2418 lists of names can have NULL entries if previously removed
2419 return flag = 0/1 for int/double
2420 return cols = 0/N for vector/array where N = # of columns
2421 ------------------------------------------------------------------------- */
2422
find_custom(const char * name,int & flag,int & cols)2423 int Atom::find_custom(const char *name, int &flag, int &cols)
2424 {
2425 if (name == nullptr) return -1;
2426
2427 for (int i = 0; i < nivector; i++)
2428 if (ivname[i] && strcmp(ivname[i],name) == 0) {
2429 flag = 0;
2430 cols = 0;
2431 return i;
2432 }
2433
2434 for (int i = 0; i < ndvector; i++)
2435 if (dvname[i] && strcmp(dvname[i],name) == 0) {
2436 flag = 1;
2437 cols = 0;
2438 return i;
2439 }
2440
2441 for (int i = 0; i < niarray; i++)
2442 if (ianame[i] && strcmp(ianame[i],name) == 0) {
2443 flag = 0;
2444 cols = icols[i];
2445 return i;
2446 }
2447
2448 for (int i = 0; i < ndarray; i++)
2449 if (daname[i] && strcmp(daname[i],name) == 0) {
2450 flag = 1;
2451 cols = dcols[i];
2452 return i;
2453 }
2454
2455 return -1;
2456 }
2457
2458 /** \brief Add a custom per-atom property with the given name and type and size
2459 \verbatim embed:rst
2460
2461 This function will add a custom per-atom property with one or more values
2462 with the name "name" to the list of custom properties.
2463 This function is called, e.g. from :doc:`fix property/atom <fix_property_atom>`.
2464 \endverbatim
2465 * \param name Name of the property (w/o a "i_" or "d_" or "i2_" or "d2_" prefix)
2466 * \param flag Data type of property: 0 for int, 1 for double
2467 * \param cols Number of values: 0 for a single value, 1 or more for a vector of values
2468 * \return index of property in the respective list of properties
2469 */
add_custom(const char * name,int flag,int cols)2470 int Atom::add_custom(const char *name, int flag, int cols)
2471 {
2472 int index = -1;
2473
2474 if ((flag == 0) && (cols == 0)) {
2475 index = nivector;
2476 nivector++;
2477 ivname = (char **) memory->srealloc(ivname,nivector*sizeof(char *),"atom:ivname");
2478 ivname[index] = utils::strdup(name);
2479 ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *),"atom:ivector");
2480 memory->create(ivector[index],nmax,"atom:ivector");
2481
2482 } else if ((flag == 1) && (cols == 0)) {
2483 index = ndvector;
2484 ndvector++;
2485 dvname = (char **) memory->srealloc(dvname,ndvector*sizeof(char *),"atom:dvname");
2486 dvname[index] = utils::strdup(name);
2487 dvector = (double **) memory->srealloc(dvector,ndvector*sizeof(double *),"atom:dvector");
2488 memory->create(dvector[index],nmax,"atom:dvector");
2489
2490 } else if ((flag == 0) && (cols > 0)) {
2491 index = niarray;
2492 niarray++;
2493 ianame = (char **) memory->srealloc(ianame,niarray*sizeof(char *),"atom:ianame");
2494 ianame[index] = utils::strdup(name);
2495 iarray = (int ***) memory->srealloc(iarray,niarray*sizeof(int **),"atom:iarray");
2496 memory->create(iarray[index],nmax,cols,"atom:iarray");
2497
2498 icols = (int *) memory->srealloc(icols,niarray*sizeof(int),"atom:icols");
2499 icols[index] = cols;
2500
2501 } else if ((flag == 1) && (cols > 0)) {
2502 index = ndarray;
2503 ndarray++;
2504 daname = (char **) memory->srealloc(daname,ndarray*sizeof(char *),"atom:daname");
2505 daname[index] = utils::strdup(name);
2506 darray = (double ***) memory->srealloc(darray,ndarray*sizeof(double **),"atom:darray");
2507 memory->create(darray[index],nmax,cols,"atom:darray");
2508
2509 dcols = (int *) memory->srealloc(dcols,ndarray*sizeof(int),"atom:dcols");
2510 dcols[index] = cols;
2511 }
2512 if (index < 0)
2513 error->all(FLERR,"Invalid call to Atom::add_custom()");
2514 return index;
2515 }
2516
2517 /*! \brief Remove a custom per-atom property of a given type and size
2518 *
2519 \verbatim embed:rst
2520 This will remove a property that was requested, e.g. by the
2521 :doc:`fix property/atom <fix_property_atom>` command. It frees the
2522 allocated memory and sets the pointer to ``nullptr`` for the entry in
2523 the list so it can be reused. The lists of these pointers are never
2524 compacted or shrink, so that indices to name mappings remain valid.
2525 \endverbatim
2526 * \param index Index of property in the respective list of properties
2527 * \param flag Data type of property: 0 for int, 1 for double
2528 * \param cols Number of values: 0 for a single value, 1 or more for a vector of values
2529 */
remove_custom(int index,int flag,int cols)2530 void Atom::remove_custom(int index, int flag, int cols)
2531 {
2532 if (flag == 0 && cols == 0) {
2533 memory->destroy(ivector[index]);
2534 ivector[index] = nullptr;
2535 delete [] ivname[index];
2536 ivname[index] = nullptr;
2537
2538 } else if (flag == 1 && cols == 0) {
2539 memory->destroy(dvector[index]);
2540 dvector[index] = nullptr;
2541 delete [] dvname[index];
2542 dvname[index] = nullptr;
2543
2544 } else if (flag == 0 && cols) {
2545 memory->destroy(iarray[index]);
2546 iarray[index] = nullptr;
2547 delete [] ianame[index];
2548 ianame[index] = nullptr;
2549
2550 } else if (flag == 1 && cols) {
2551 memory->destroy(darray[index]);
2552 darray[index] = nullptr;
2553 delete [] daname[index];
2554 daname[index] = nullptr;
2555 }
2556 }
2557
2558 /** Provide access to internal data of the Atom class by keyword
2559 *
2560 \verbatim embed:rst
2561
2562 This function is a way to access internal per-atom data. This data is
2563 distributed across MPI ranks and thus only the data for "local" atoms
2564 are expected to be available. Whether also data for "ghost" atoms is
2565 stored and up-to-date depends on various simulation settings.
2566
2567 This table lists a large part of the supported names, their data types,
2568 length of the data area, and a short description.
2569
2570 .. list-table::
2571 :header-rows: 1
2572 :widths: auto
2573
2574 * - Name
2575 - Type
2576 - Items per atom
2577 - Description
2578 * - mass
2579 - double
2580 - 1
2581 - per-type mass. This array is **NOT** a per-atom array
2582 but of length ``ntypes+1``, element 0 is ignored.
2583 * - id
2584 - tagint
2585 - 1
2586 - atom ID of the particles
2587 * - type
2588 - int
2589 - 1
2590 - atom type of the particles
2591 * - mask
2592 - int
2593 - 1
2594 - bitmask for mapping to groups. Individual bits are set
2595 to 0 or 1 for each group.
2596 * - image
2597 - imageint
2598 - 1
2599 - 3 image flags encoded into a single integer.
2600 See :cpp:func:`lammps_encode_image_flags`.
2601 * - x
2602 - double
2603 - 3
2604 - x-, y-, and z-coordinate of the particles
2605 * - v
2606 - double
2607 - 3
2608 - x-, y-, and z-component of the velocity of the particles
2609 * - f
2610 - double
2611 - 3
2612 - x-, y-, and z-component of the force on the particles
2613 * - molecule
2614 - int
2615 - 1
2616 - molecule ID of the particles
2617 * - q
2618 - double
2619 - 1
2620 - charge of the particles
2621 * - mu
2622 - double
2623 - 3
2624 - dipole moment of the particles
2625 * - omega
2626 - double
2627 - 3
2628 - x-, y-, and z-component of rotational velocity of the particles
2629 * - angmom
2630 - double
2631 - 3
2632 - x-, y-, and z-component of angular momentum of the particles
2633 * - torque
2634 - double
2635 - 3
2636 - x-, y-, and z-component of the torque on the particles
2637 * - radius
2638 - double
2639 - 1
2640 - radius of the (extended) particles
2641 * - rmass
2642 - double
2643 - 1
2644 - per-atom mass of the particles. ``nullptr`` if per-type masses are
2645 used. See the :cpp:func:`rmass_flag<lammps_extract_setting>` setting.
2646 * - ellipsoid
2647 - int
2648 - 1
2649 - 1 if the particle is an ellipsoidal particle, 0 if not
2650 * - line
2651 - int
2652 - 1
2653 - 1 if the particle is a line particle, 0 if not
2654 * - tri
2655 - int
2656 - 1
2657 - 1 if the particle is a triangulated particle, 0 if not
2658 * - body
2659 - int
2660 - 1
2661 - 1 if the particle is a body particle, 0 if not
2662 * - i_name
2663 - int
2664 - 1
2665 - single integer value defined by fix property/atom vector name
2666 * - d_name
2667 - double
2668 - 1
2669 - single double value defined by fix property/atom vector name
2670 * - i2_name
2671 - int
2672 - n
2673 - N integer values defined by fix property/atom array name
2674 * - d2_name
2675 - double
2676 - n
2677 - N double values defined by fix property/atom array name
2678
2679 *See also*
2680 :cpp:func:`lammps_extract_atom`
2681
2682 \endverbatim
2683 *
2684 * \sa extract_datatype
2685 *
2686 * \param name string with the keyword of the desired property.
2687 Typically the name of the pointer variable returned
2688 * \return pointer to the requested data cast to ``void *`` or ``nullptr`` */
2689
extract(const char * name)2690 void *Atom::extract(const char *name)
2691 {
2692 // --------------------------------------------------------------------
2693 // 4th customization section: customize by adding new variable name
2694 // if new variable is from a package, add package comment
2695
2696 if (strcmp(name,"mass") == 0) return (void *) mass;
2697
2698 if (strcmp(name,"id") == 0) return (void *) tag;
2699 if (strcmp(name,"type") == 0) return (void *) type;
2700 if (strcmp(name,"mask") == 0) return (void *) mask;
2701 if (strcmp(name,"image") == 0) return (void *) image;
2702 if (strcmp(name,"x") == 0) return (void *) x;
2703 if (strcmp(name,"v") == 0) return (void *) v;
2704 if (strcmp(name,"f") == 0) return (void *) f;
2705 if (strcmp(name,"molecule") == 0) return (void *) molecule;
2706 if (strcmp(name,"q") == 0) return (void *) q;
2707 if (strcmp(name,"mu") == 0) return (void *) mu;
2708 if (strcmp(name,"omega") == 0) return (void *) omega;
2709 if (strcmp(name,"angmom") == 0) return (void *) angmom;
2710 if (strcmp(name,"torque") == 0) return (void *) torque;
2711 if (strcmp(name,"radius") == 0) return (void *) radius;
2712 if (strcmp(name,"rmass") == 0) return (void *) rmass;
2713 if (strcmp(name,"ellipsoid") == 0) return (void *) ellipsoid;
2714 if (strcmp(name,"line") == 0) return (void *) line;
2715 if (strcmp(name,"tri") == 0) return (void *) tri;
2716 if (strcmp(name,"body") == 0) return (void *) body;
2717
2718 if (strcmp(name,"vfrac") == 0) return (void *) vfrac;
2719 if (strcmp(name,"s0") == 0) return (void *) s0;
2720 if (strcmp(name,"x0") == 0) return (void *) x0;
2721
2722 if (strcmp(name,"spin") == 0) return (void *) spin;
2723 if (strcmp(name,"eradius") == 0) return (void *) eradius;
2724 if (strcmp(name,"ervel") == 0) return (void *) ervel;
2725 if (strcmp(name,"erforce") == 0) return (void *) erforce;
2726 if (strcmp(name,"ervelforce") == 0) return (void *) ervelforce;
2727 if (strcmp(name,"cs") == 0) return (void *) cs;
2728 if (strcmp(name,"csforce") == 0) return (void *) csforce;
2729 if (strcmp(name,"vforce") == 0) return (void *) vforce;
2730 if (strcmp(name,"etag") == 0) return (void *) etag;
2731
2732 if (strcmp(name,"rho") == 0) return (void *) rho;
2733 if (strcmp(name,"drho") == 0) return (void *) drho;
2734 if (strcmp(name,"esph") == 0) return (void *) esph;
2735 if (strcmp(name,"desph") == 0) return (void *) desph;
2736 if (strcmp(name,"cv") == 0) return (void *) cv;
2737 if (strcmp(name,"vest") == 0) return (void *) vest;
2738
2739 // MESONT package
2740
2741 if (strcmp(name,"length") == 0) return (void *) length;
2742 if (strcmp(name,"buckling") == 0) return (void *) buckling;
2743 if (strcmp(name,"bond_nt") == 0) return (void *) bond_nt;
2744
2745 // MACHDYN package
2746
2747 if (strcmp(name, "contact_radius") == 0) return (void *) contact_radius;
2748 if (strcmp(name, "smd_data_9") == 0) return (void *) smd_data_9;
2749 if (strcmp(name, "smd_stress") == 0) return (void *) smd_stress;
2750 if (strcmp(name, "eff_plastic_strain") == 0)
2751 return (void *) eff_plastic_strain;
2752 if (strcmp(name, "eff_plastic_strain_rate") == 0)
2753 return (void *) eff_plastic_strain_rate;
2754 if (strcmp(name, "damage") == 0) return (void *) damage;
2755
2756 // DPD-REACT pakage
2757
2758 if (strcmp(name,"dpdTheta") == 0) return (void *) dpdTheta;
2759
2760 // DPD-MESO package
2761
2762 if (strcmp(name,"edpd_temp") == 0) return (void *) edpd_temp;
2763
2764 // DIELECTRIC package
2765
2766 if (strcmp(name,"area") == 0) return (void *) area;
2767 if (strcmp(name,"ed") == 0) return (void *) ed;
2768 if (strcmp(name,"em") == 0) return (void *) em;
2769 if (strcmp(name,"epsilon") == 0) return (void *) epsilon;
2770 if (strcmp(name,"curvature") == 0) return (void *) curvature;
2771 if (strcmp(name,"q_unscaled") == 0) return (void *) q_unscaled;
2772
2773 // end of customization section
2774 // --------------------------------------------------------------------
2775
2776 // custom vectors and arrays
2777
2778 if (utils::strmatch(name,"^[id]2?_")) {
2779 int which = 0, array = 0;
2780 if (name[0] == 'd') which = 1;
2781 if (name[1] == '2') array = 1;
2782
2783 int index,flag,cols;
2784 if (!array) index = find_custom(&name[2],flag,cols);
2785 else index = find_custom(&name[3],flag,cols);
2786
2787 if (index < 0) return nullptr;
2788 if (which != flag) return nullptr;
2789 if ((!array && cols) || (array && !cols)) return nullptr;
2790
2791 if (!which && !array) return (void *) ivector[index];
2792 if (which && !array) return (void *) dvector[index];
2793 if (!which && array) return (void *) iarray[index];
2794 if (which && array) return (void *) darray[index];
2795 }
2796
2797 return nullptr;
2798 }
2799
2800 /** Provide data type info about internal data of the Atom class
2801 *
2802 \verbatim embed:rst
2803
2804 .. versionadded:: 18Sep2020
2805
2806 \endverbatim
2807 *
2808 * \sa extract
2809 *
2810 * \param name string with the keyword of the desired property.
2811 * \return data type constant for desired property or -1 */
2812
extract_datatype(const char * name)2813 int Atom::extract_datatype(const char *name)
2814 {
2815 // --------------------------------------------------------------------
2816 // 5th customization section: customize by adding new variable name
2817
2818 if (strcmp(name,"mass") == 0) return LAMMPS_DOUBLE;
2819
2820 if (strcmp(name,"id") == 0) return LAMMPS_TAGINT;
2821 if (strcmp(name,"type") == 0) return LAMMPS_INT;
2822 if (strcmp(name,"mask") == 0) return LAMMPS_INT;
2823 if (strcmp(name,"image") == 0) return LAMMPS_TAGINT;
2824 if (strcmp(name,"x") == 0) return LAMMPS_DOUBLE_2D;
2825 if (strcmp(name,"v") == 0) return LAMMPS_DOUBLE_2D;
2826 if (strcmp(name,"f") == 0) return LAMMPS_DOUBLE_2D;
2827 if (strcmp(name,"molecule") == 0) return LAMMPS_TAGINT;
2828 if (strcmp(name,"q") == 0) return LAMMPS_DOUBLE;
2829 if (strcmp(name,"mu") == 0) return LAMMPS_DOUBLE_2D;
2830 if (strcmp(name,"omega") == 0) return LAMMPS_DOUBLE_2D;
2831 if (strcmp(name,"angmom") == 0) return LAMMPS_DOUBLE_2D;
2832 if (strcmp(name,"torque") == 0) return LAMMPS_DOUBLE_2D;
2833 if (strcmp(name,"radius") == 0) return LAMMPS_DOUBLE;
2834 if (strcmp(name,"rmass") == 0) return LAMMPS_DOUBLE;
2835 if (strcmp(name,"ellipsoid") == 0) return LAMMPS_INT;
2836 if (strcmp(name,"line") == 0) return LAMMPS_INT;
2837 if (strcmp(name,"tri") == 0) return LAMMPS_INT;
2838 if (strcmp(name,"body") == 0) return LAMMPS_INT;
2839
2840 if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE;
2841 if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE;
2842 if (strcmp(name,"x0") == 0) return LAMMPS_DOUBLE_2D;
2843
2844 if (strcmp(name,"spin") == 0) return LAMMPS_INT;
2845 if (strcmp(name,"eradius") == 0) return LAMMPS_DOUBLE;
2846 if (strcmp(name,"ervel") == 0) return LAMMPS_DOUBLE;
2847 if (strcmp(name,"erforce") == 0) return LAMMPS_DOUBLE;
2848 if (strcmp(name,"ervelforce") == 0) return LAMMPS_DOUBLE;
2849 if (strcmp(name,"cs") == 0) return LAMMPS_DOUBLE_2D;
2850 if (strcmp(name,"csforce") == 0) return LAMMPS_DOUBLE_2D;
2851 if (strcmp(name,"vforce") == 0) return LAMMPS_DOUBLE_2D;
2852 if (strcmp(name,"etag") == 0) return LAMMPS_INT;
2853
2854 if (strcmp(name,"rho") == 0) return LAMMPS_DOUBLE;
2855 if (strcmp(name,"drho") == 0) return LAMMPS_DOUBLE;
2856 if (strcmp(name,"esph") == 0) return LAMMPS_DOUBLE;
2857 if (strcmp(name,"desph") == 0) return LAMMPS_DOUBLE;
2858 if (strcmp(name,"cv") == 0) return LAMMPS_DOUBLE;
2859 if (strcmp(name,"vest") == 0) return LAMMPS_DOUBLE_2D;
2860
2861 // MESONT package
2862
2863 if (strcmp(name,"length") == 0) return LAMMPS_DOUBLE;
2864 if (strcmp(name,"buckling") == 0) return LAMMPS_INT;
2865 if (strcmp(name,"bond_nt") == 0) return LAMMPS_TAGINT_2D;
2866
2867 // MACHDYN package
2868
2869 if (strcmp(name, "contact_radius") == 0) return LAMMPS_DOUBLE;
2870 if (strcmp(name, "smd_data_9") == 0) return LAMMPS_DOUBLE_2D;
2871 if (strcmp(name, "smd_stress") == 0) return LAMMPS_DOUBLE_2D;
2872 if (strcmp(name, "eff_plastic_strain") == 0) return LAMMPS_DOUBLE;
2873 if (strcmp(name, "eff_plastic_strain_rate") == 0) return LAMMPS_DOUBLE;
2874 if (strcmp(name, "damage") == 0) return LAMMPS_DOUBLE;
2875
2876 // DPD-REACT package
2877
2878 if (strcmp(name,"dpdTheta") == 0) return LAMMPS_DOUBLE;
2879
2880 // DPD-MESO package
2881
2882 if (strcmp(name,"edpd_temp") == 0) return LAMMPS_DOUBLE;
2883
2884 // DIELECTRIC package
2885
2886 if (strcmp(name,"area") == 0) return LAMMPS_DOUBLE;
2887 if (strcmp(name,"ed") == 0) return LAMMPS_DOUBLE;
2888 if (strcmp(name,"em") == 0) return LAMMPS_DOUBLE;
2889 if (strcmp(name,"epsilon") == 0) return LAMMPS_DOUBLE;
2890 if (strcmp(name,"curvature") == 0) return LAMMPS_DOUBLE;
2891 if (strcmp(name,"q_unscaled") == 0) return LAMMPS_DOUBLE;
2892
2893 // end of customization section
2894 // --------------------------------------------------------------------
2895
2896 // custom vectors and arrays
2897
2898 if (utils::strmatch(name,"^[id]2?_")) {
2899 int which = 0, array = 0;
2900 if (name[0] == 'd') which = 1;
2901 if (name[1] == '2') array = 1;
2902
2903 int index,flag,cols;
2904 if (!array) index = find_custom(&name[2],flag,cols);
2905 else index = find_custom(&name[3],flag,cols);
2906
2907 if (index < 0) return -1;
2908 if (which != flag) return -1;
2909 if ((!array && cols) || (array && !cols)) return -1;
2910
2911 if (which == 0) return LAMMPS_INT;
2912 else return LAMMPS_DOUBLE;
2913 }
2914
2915 return -1;
2916 }
2917
2918 /* ----------------------------------------------------------------------
2919 return # of bytes of allocated memory
2920 call to avec tallies per-atom vectors
2921 add in global to local mapping storage
2922 ------------------------------------------------------------------------- */
2923
memory_usage()2924 double Atom::memory_usage()
2925 {
2926 double bytes = avec->memory_usage();
2927
2928 bytes += (double)max_same*sizeof(int);
2929 if (map_style == MAP_ARRAY)
2930 bytes += memory->usage(map_array,map_maxarray);
2931 else if (map_style == MAP_HASH) {
2932 bytes += (double)map_nbucket*sizeof(int);
2933 bytes += (double)map_nhash*sizeof(HashElem);
2934 }
2935 if (maxnext) {
2936 bytes += memory->usage(next,maxnext);
2937 bytes += memory->usage(permute,maxnext);
2938 }
2939
2940 return bytes;
2941 }
2942
2943