1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 This file is from LAMMPS
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 #include <stdlib.h>
47 #include "atom_vec_bond.h"
48 #include "atom.h"
49 #include "comm.h"
50 #include "domain.h"
51 #include "modify.h"
52 #include "fix.h"
53 #include "memory.h"
54 #include "error.h"
55
56 using namespace LAMMPS_NS;
57
58 #define DELTA 10000
59
60 /* ---------------------------------------------------------------------- */
61
AtomVecBond(LAMMPS * lmp)62 AtomVecBond::AtomVecBond(LAMMPS *lmp) : AtomVec(lmp)
63 {
64 molecular = 1;
65 bonds_allow = 1;
66 mass_type = 1;
67
68 comm_x_only = comm_f_only = 1;
69 size_forward = 3;
70 size_reverse = 3;
71 size_border = 7;
72 size_velocity = 3;
73 size_data_atom = 6;
74 size_data_vel = 4;
75 xcol_data = 4;
76
77 atom->molecule_flag = 1;
78 }
79
80 /* ----------------------------------------------------------------------
81 grow atom arrays
82 n = 0 grows arrays by DELTA
83 n > 0 allocates arrays to size n
84 ------------------------------------------------------------------------- */
85
grow(int n)86 void AtomVecBond::grow(int n)
87 {
88 if (n == 0) nmax += DELTA;
89 else nmax = n;
90 atom->nmax = nmax;
91 if (nmax < 0 || nmax > MAXSMALLINT)
92 error->one(FLERR,"Per-processor system is too big");
93
94 tag = memory->grow(atom->tag,nmax,"atom:tag");
95 type = memory->grow(atom->type,nmax,"atom:type");
96 mask = memory->grow(atom->mask,nmax,"atom:mask");
97 image = memory->grow(atom->image,nmax,"atom:image");
98 x = memory->grow(atom->x,nmax,3,"atom:x");
99 v = memory->grow(atom->v,nmax,3,"atom:v");
100 f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
101
102 molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
103
104 nspecial = memory->grow(atom->nspecial,nmax,3,"atom:nspecial");
105 special = memory->grow(atom->special,nmax,atom->maxspecial,"atom:special");
106
107 num_bond = memory->grow(atom->num_bond,nmax,"atom:num_bond");
108 bond_type = memory->grow(atom->bond_type,nmax,atom->bond_per_atom,
109 "atom:bond_type");
110 bond_atom = memory->grow(atom->bond_atom,nmax,atom->bond_per_atom,
111 "atom:bond_atom");
112
113 if (atom->nextra_grow)
114 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
115 modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
116 }
117
118 /* ----------------------------------------------------------------------
119 reset local array ptrs
120 ------------------------------------------------------------------------- */
121
grow_reset()122 void AtomVecBond::grow_reset()
123 {
124 tag = atom->tag; type = atom->type;
125 mask = atom->mask; image = atom->image;
126 x = atom->x; v = atom->v; f = atom->f;
127 molecule = atom->molecule;
128 nspecial = atom->nspecial; special = atom->special;
129 num_bond = atom->num_bond; bond_type = atom->bond_type;
130 bond_atom = atom->bond_atom;
131 }
132
133 /* ----------------------------------------------------------------------
134 copy atom I info to atom J
135 ------------------------------------------------------------------------- */
136
copy(int i,int j,int delflag)137 void AtomVecBond::copy(int i, int j, int delflag)
138 {
139 int k;
140
141 tag[j] = tag[i];
142 type[j] = type[i];
143 mask[j] = mask[i];
144 image[j] = image[i];
145 x[j][0] = x[i][0];
146 x[j][1] = x[i][1];
147 x[j][2] = x[i][2];
148 v[j][0] = v[i][0];
149 v[j][1] = v[i][1];
150 v[j][2] = v[i][2];
151
152 molecule[j] = molecule[i];
153
154 num_bond[j] = num_bond[i];
155 for (k = 0; k < num_bond[j]; k++) {
156 bond_type[j][k] = bond_type[i][k];
157 bond_atom[j][k] = bond_atom[i][k];
158 }
159
160 nspecial[j][0] = nspecial[i][0];
161 nspecial[j][1] = nspecial[i][1];
162 nspecial[j][2] = nspecial[i][2];
163 for (k = 0; k < nspecial[j][2]; k++) special[j][k] = special[i][k];
164
165 if (atom->nextra_grow)
166 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
167 modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j,delflag);
168 }
169
170 /* ---------------------------------------------------------------------- */
171
pack_comm(int n,int * list,double * buf,int pbc_flag,int * pbc)172 int AtomVecBond::pack_comm(int n, int *list, double *buf,
173 int pbc_flag, int *pbc)
174 {
175 int i,j,m;
176 double dx,dy,dz;
177
178 m = 0;
179 if (pbc_flag == 0) {
180 for (i = 0; i < n; i++) {
181 j = list[i];
182 buf[m++] = x[j][0];
183 buf[m++] = x[j][1];
184 buf[m++] = x[j][2];
185 }
186 } else {
187 if (domain->triclinic == 0) {
188 dx = pbc[0]*domain->xprd;
189 dy = pbc[1]*domain->yprd;
190 dz = pbc[2]*domain->zprd;
191 } else {
192 dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
193 dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
194 dz = pbc[2]*domain->zprd;
195 }
196 for (i = 0; i < n; i++) {
197 j = list[i];
198 buf[m++] = x[j][0] + dx;
199 buf[m++] = x[j][1] + dy;
200 buf[m++] = x[j][2] + dz;
201 }
202 }
203 return m;
204 }
205
206 /* ---------------------------------------------------------------------- */
207
pack_comm_vel(int n,int * list,double * buf,int pbc_flag,int * pbc)208 int AtomVecBond::pack_comm_vel(int n, int *list, double *buf,
209 int pbc_flag, int *pbc)
210 {
211 int i,j,m;
212 double dx,dy,dz,dvx,dvy,dvz;
213
214 m = 0;
215 if (pbc_flag == 0) {
216 for (i = 0; i < n; i++) {
217 j = list[i];
218 buf[m++] = x[j][0];
219 buf[m++] = x[j][1];
220 buf[m++] = x[j][2];
221 buf[m++] = v[j][0];
222 buf[m++] = v[j][1];
223 buf[m++] = v[j][2];
224 }
225 } else {
226 if (domain->triclinic == 0) {
227 dx = pbc[0]*domain->xprd;
228 dy = pbc[1]*domain->yprd;
229 dz = pbc[2]*domain->zprd;
230 } else {
231 dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
232 dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
233 dz = pbc[2]*domain->zprd;
234 }
235 if (!deform_vremap) {
236 for (i = 0; i < n; i++) {
237 j = list[i];
238 buf[m++] = x[j][0] + dx;
239 buf[m++] = x[j][1] + dy;
240 buf[m++] = x[j][2] + dz;
241 buf[m++] = v[j][0];
242 buf[m++] = v[j][1];
243 buf[m++] = v[j][2];
244 }
245 } else {
246 dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
247 dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
248 dvz = pbc[2]*h_rate[2];
249 for (i = 0; i < n; i++) {
250 j = list[i];
251 buf[m++] = x[j][0] + dx;
252 buf[m++] = x[j][1] + dy;
253 buf[m++] = x[j][2] + dz;
254 if (mask[i] & deform_groupbit) {
255 buf[m++] = v[j][0] + dvx;
256 buf[m++] = v[j][1] + dvy;
257 buf[m++] = v[j][2] + dvz;
258 } else {
259 buf[m++] = v[j][0];
260 buf[m++] = v[j][1];
261 buf[m++] = v[j][2];
262 }
263 }
264 }
265 }
266 return m;
267 }
268
269 /* ---------------------------------------------------------------------- */
270
unpack_comm(int n,int first,double * buf)271 void AtomVecBond::unpack_comm(int n, int first, double *buf)
272 {
273 int i,m,last;
274
275 m = 0;
276 last = first + n;
277 for (i = first; i < last; i++) {
278 x[i][0] = buf[m++];
279 x[i][1] = buf[m++];
280 x[i][2] = buf[m++];
281 }
282 }
283
284 /* ---------------------------------------------------------------------- */
285
unpack_comm_vel(int n,int first,double * buf)286 void AtomVecBond::unpack_comm_vel(int n, int first, double *buf)
287 {
288 int i,m,last;
289
290 m = 0;
291 last = first + n;
292 for (i = first; i < last; i++) {
293 x[i][0] = buf[m++];
294 x[i][1] = buf[m++];
295 x[i][2] = buf[m++];
296 v[i][0] = buf[m++];
297 v[i][1] = buf[m++];
298 v[i][2] = buf[m++];
299 }
300 }
301
302 /* ---------------------------------------------------------------------- */
303
pack_reverse(int n,int first,double * buf)304 int AtomVecBond::pack_reverse(int n, int first, double *buf)
305 {
306 int i,m,last;
307
308 m = 0;
309 last = first + n;
310 for (i = first; i < last; i++) {
311 buf[m++] = f[i][0];
312 buf[m++] = f[i][1];
313 buf[m++] = f[i][2];
314 }
315 return m;
316 }
317
318 /* ---------------------------------------------------------------------- */
319
unpack_reverse(int n,int * list,double * buf)320 void AtomVecBond::unpack_reverse(int n, int *list, double *buf)
321 {
322 int i,j,m;
323
324 m = 0;
325 for (i = 0; i < n; i++) {
326 j = list[i];
327 f[j][0] += buf[m++];
328 f[j][1] += buf[m++];
329 f[j][2] += buf[m++];
330 }
331 }
332
333 /* ---------------------------------------------------------------------- */
334
pack_border(int n,int * list,double * buf,int pbc_flag,int * pbc)335 int AtomVecBond::pack_border(int n, int *list, double *buf,
336 int pbc_flag, int *pbc)
337 {
338 int i,j,m;
339 double dx,dy,dz;
340
341 m = 0;
342 if (pbc_flag == 0) {
343 for (i = 0; i < n; i++) {
344 j = list[i];
345 buf[m++] = x[j][0];
346 buf[m++] = x[j][1];
347 buf[m++] = x[j][2];
348 buf[m++] = ubuf(tag[j]).d;
349 buf[m++] = ubuf(type[j]).d;
350 buf[m++] = ubuf(mask[j]).d;
351 buf[m++] = ubuf(molecule[j]).d;
352 }
353 } else {
354 if (domain->triclinic == 0) {
355 dx = pbc[0]*domain->xprd;
356 dy = pbc[1]*domain->yprd;
357 dz = pbc[2]*domain->zprd;
358 } else {
359 dx = pbc[0];
360 dy = pbc[1];
361 dz = pbc[2];
362 }
363 for (i = 0; i < n; i++) {
364 j = list[i];
365 buf[m++] = x[j][0] + dx;
366 buf[m++] = x[j][1] + dy;
367 buf[m++] = x[j][2] + dz;
368 buf[m++] = ubuf(tag[j]).d;
369 buf[m++] = ubuf(type[j]).d;
370 buf[m++] = ubuf(mask[j]).d;
371 buf[m++] = ubuf(molecule[j]).d;
372 }
373 }
374
375 if (atom->nextra_border)
376 for (int iextra = 0; iextra < atom->nextra_border; iextra++)
377 m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
378
379 return m;
380 }
381
382 /* ---------------------------------------------------------------------- */
383
pack_border_vel(int n,int * list,double * buf,int pbc_flag,int * pbc)384 int AtomVecBond::pack_border_vel(int n, int *list, double *buf,
385 int pbc_flag, int *pbc)
386 {
387 int i,j,m;
388 double dx,dy,dz,dvx,dvy,dvz;
389
390 m = 0;
391 if (pbc_flag == 0) {
392 for (i = 0; i < n; i++) {
393 j = list[i];
394 buf[m++] = x[j][0];
395 buf[m++] = x[j][1];
396 buf[m++] = x[j][2];
397 buf[m++] = ubuf(tag[j]).d;
398 buf[m++] = ubuf(type[j]).d;
399 buf[m++] = ubuf(mask[j]).d;
400 buf[m++] = ubuf(molecule[j]).d;
401 buf[m++] = v[j][0];
402 buf[m++] = v[j][1];
403 buf[m++] = v[j][2];
404 }
405 } else {
406 if (domain->triclinic == 0) {
407 dx = pbc[0]*domain->xprd;
408 dy = pbc[1]*domain->yprd;
409 dz = pbc[2]*domain->zprd;
410 } else {
411 dx = pbc[0];
412 dy = pbc[1];
413 dz = pbc[2];
414 }
415 if (!deform_vremap) {
416 for (i = 0; i < n; i++) {
417 j = list[i];
418 buf[m++] = x[j][0] + dx;
419 buf[m++] = x[j][1] + dy;
420 buf[m++] = x[j][2] + dz;
421 buf[m++] = ubuf(tag[j]).d;
422 buf[m++] = ubuf(type[j]).d;
423 buf[m++] = ubuf(mask[j]).d;
424 buf[m++] = ubuf(molecule[j]).d;
425 buf[m++] = v[j][0];
426 buf[m++] = v[j][1];
427 buf[m++] = v[j][2];
428 }
429 } else {
430 dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
431 dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
432 dvz = pbc[2]*h_rate[2];
433 for (i = 0; i < n; i++) {
434 j = list[i];
435 buf[m++] = x[j][0] + dx;
436 buf[m++] = x[j][1] + dy;
437 buf[m++] = x[j][2] + dz;
438 buf[m++] = ubuf(tag[j]).d;
439 buf[m++] = ubuf(type[j]).d;
440 buf[m++] = ubuf(mask[j]).d;
441 buf[m++] = ubuf(molecule[j]).d;
442 if (mask[i] & deform_groupbit) {
443 buf[m++] = v[j][0] + dvx;
444 buf[m++] = v[j][1] + dvy;
445 buf[m++] = v[j][2] + dvz;
446 } else {
447 buf[m++] = v[j][0];
448 buf[m++] = v[j][1];
449 buf[m++] = v[j][2];
450 }
451 }
452 }
453 }
454
455 if (atom->nextra_border)
456 for (int iextra = 0; iextra < atom->nextra_border; iextra++)
457 m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
458
459 return m;
460 }
461
462 /* ---------------------------------------------------------------------- */
463
pack_border_hybrid(int n,int * list,double * buf)464 int AtomVecBond::pack_border_hybrid(int n, int *list, double *buf)
465 {
466 int i,j,m;
467
468 m = 0;
469 for (i = 0; i < n; i++) {
470 j = list[i];
471 buf[m++] = molecule[j];
472 }
473 return m;
474 }
475
476 /* ---------------------------------------------------------------------- */
477
unpack_border(int n,int first,double * buf)478 void AtomVecBond::unpack_border(int n, int first, double *buf)
479 {
480 int i,m,last;
481
482 m = 0;
483 last = first + n;
484 for (i = first; i < last; i++) {
485 if (i == nmax) grow(0);
486 x[i][0] = buf[m++];
487 x[i][1] = buf[m++];
488 x[i][2] = buf[m++];
489 tag[i] = (int) ubuf(buf[m++]).i;
490 type[i] = (int) ubuf(buf[m++]).i;
491 mask[i] = (int) ubuf(buf[m++]).i;
492 molecule[i] = (int) ubuf(buf[m++]).i;
493 }
494
495 if (atom->nextra_border)
496 for (int iextra = 0; iextra < atom->nextra_border; iextra++)
497 m += modify->fix[atom->extra_border[iextra]]->
498 unpack_border(n,first,&buf[m]);
499 }
500
501 /* ---------------------------------------------------------------------- */
502
unpack_border_vel(int n,int first,double * buf)503 void AtomVecBond::unpack_border_vel(int n, int first, double *buf)
504 {
505 int i,m,last;
506
507 m = 0;
508 last = first + n;
509 for (i = first; i < last; i++) {
510 if (i == nmax) grow(0);
511 x[i][0] = buf[m++];
512 x[i][1] = buf[m++];
513 x[i][2] = buf[m++];
514 tag[i] = (int) ubuf(buf[m++]).i;
515 type[i] = (int) ubuf(buf[m++]).i;
516 mask[i] = (int) ubuf(buf[m++]).i;
517 molecule[i] = (int) ubuf(buf[m++]).i;
518 v[i][0] = buf[m++];
519 v[i][1] = buf[m++];
520 v[i][2] = buf[m++];
521 }
522
523 if (atom->nextra_border)
524 for (int iextra = 0; iextra < atom->nextra_border; iextra++)
525 m += modify->fix[atom->extra_border[iextra]]->
526 unpack_border(n,first,&buf[m]);
527 }
528
529 /* ---------------------------------------------------------------------- */
530
unpack_border_hybrid(int n,int first,double * buf)531 int AtomVecBond::unpack_border_hybrid(int n, int first, double *buf)
532 {
533 int i,m,last;
534
535 m = 0;
536 last = first + n;
537 for (i = first; i < last; i++)
538 molecule[i] = (int) ubuf(buf[m++]).i;
539 return m;
540 }
541
542 /* ----------------------------------------------------------------------
543 pack data for atom I for sending to another proc
544 xyz must be 1st 3 values, so comm::exchange() can test on them
545 ------------------------------------------------------------------------- */
546
pack_exchange(int i,double * buf)547 int AtomVecBond::pack_exchange(int i, double *buf)
548 {
549 int k;
550
551 int m = 1;
552 buf[m++] = x[i][0];
553 buf[m++] = x[i][1];
554 buf[m++] = x[i][2];
555 buf[m++] = v[i][0];
556 buf[m++] = v[i][1];
557 buf[m++] = v[i][2];
558 buf[m++] = ubuf(tag[i]).d;
559 buf[m++] = ubuf(type[i]).d;
560 buf[m++] = ubuf(mask[i]).d;
561 buf[m++] = ubuf(image[i]).d;
562
563 buf[m++] = ubuf(molecule[i]).d;
564
565 buf[m++] = ubuf(num_bond[i]).d;
566 for (k = 0; k < num_bond[i]; k++) {
567 buf[m++] = ubuf(bond_type[i][k]).d;
568 buf[m++] = ubuf(bond_atom[i][k]).d;
569 }
570
571 buf[m++] = ubuf(nspecial[i][0]).d;
572 buf[m++] = ubuf(nspecial[i][1]).d;
573 buf[m++] = ubuf(nspecial[i][2]).d;
574 for (k = 0; k < nspecial[i][2]; k++) buf[m++] = ubuf(special[i][k]).d;
575
576 if (atom->nextra_grow)
577 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
578 m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
579
580 buf[0] = m;
581 return m;
582 }
583
584 /* ---------------------------------------------------------------------- */
585
unpack_exchange(double * buf)586 int AtomVecBond::unpack_exchange(double *buf)
587 {
588 int k;
589
590 int nlocal = atom->nlocal;
591 if (nlocal == nmax) grow(0);
592
593 int m = 1;
594 x[nlocal][0] = buf[m++];
595 x[nlocal][1] = buf[m++];
596 x[nlocal][2] = buf[m++];
597 v[nlocal][0] = buf[m++];
598 v[nlocal][1] = buf[m++];
599 v[nlocal][2] = buf[m++];
600 tag[nlocal] = (int) ubuf(buf[m++]).i;
601 type[nlocal] = (int) ubuf(buf[m++]).i;
602 mask[nlocal] = (int) ubuf(buf[m++]).i;
603 image[nlocal] = (tagint) ubuf(buf[m++]).i;
604
605 molecule[nlocal] = (int) ubuf(buf[m++]).i;
606
607 num_bond[nlocal] = (int) ubuf(buf[m++]).i;
608 for (k = 0; k < num_bond[nlocal]; k++) {
609 bond_type[nlocal][k] = (int) ubuf(buf[m++]).i;
610 bond_atom[nlocal][k] = (int) ubuf(buf[m++]).i;
611 }
612
613 nspecial[nlocal][0] = (int) ubuf(buf[m++]).i;
614 nspecial[nlocal][1] = (int) ubuf(buf[m++]).i;
615 nspecial[nlocal][2] = (int) ubuf(buf[m++]).i;
616 for (k = 0; k < nspecial[nlocal][2]; k++)
617 special[nlocal][k] = (int) ubuf(buf[m++]).i;
618
619 if (atom->nextra_grow)
620 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
621 m += modify->fix[atom->extra_grow[iextra]]->
622 unpack_exchange(nlocal,&buf[m]);
623
624 atom->nlocal++;
625 return m;
626 }
627
628 /* ----------------------------------------------------------------------
629 size of restart data for all atoms owned by this proc
630 include extra data stored by fixes
631 ------------------------------------------------------------------------- */
632
size_restart()633 int AtomVecBond::size_restart()
634 {
635 int i;
636
637 int nlocal = atom->nlocal;
638 int n = 0;
639 for (i = 0; i < nlocal; i++)
640 n += 13 + 2*num_bond[i];
641
642 if (atom->nextra_restart)
643 for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
644 for (i = 0; i < nlocal; i++)
645 n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
646
647 return n;
648 }
649
650 /* ----------------------------------------------------------------------
651 pack atom I's data for restart file including extra quantities
652 xyz must be 1st 3 values, so that read_restart can test on them
653 molecular types may be negative, but write as positive
654 ------------------------------------------------------------------------- */
655
pack_restart(int i,double * buf)656 int AtomVecBond::pack_restart(int i, double *buf)
657 {
658 int k;
659
660 int m = 1;
661 buf[m++] = x[i][0];
662 buf[m++] = x[i][1];
663 buf[m++] = x[i][2];
664 buf[m++] = ubuf(tag[i]).d;
665 buf[m++] = ubuf(type[i]).d;
666 buf[m++] = ubuf(mask[i]).d;
667 buf[m++] = ubuf(image[i]).d;
668 buf[m++] = v[i][0];
669 buf[m++] = v[i][1];
670 buf[m++] = v[i][2];
671
672 buf[m++] = ubuf(molecule[i]).d;
673
674 buf[m++] = ubuf(num_bond[i]).d;
675 for (k = 0; k < num_bond[i]; k++) {
676 buf[m++] = ubuf(MAX(bond_type[i][k],-bond_type[i][k])).d;
677 buf[m++] = ubuf(bond_atom[i][k]).d;
678 }
679
680 if (atom->nextra_restart)
681 for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
682 m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
683
684 buf[0] = m;
685 return m;
686 }
687
688 /* ----------------------------------------------------------------------
689 unpack data for one atom from restart file including extra quantities
690 ------------------------------------------------------------------------- */
691
unpack_restart(double * buf)692 int AtomVecBond::unpack_restart(double *buf)
693 {
694 int k;
695
696 int nlocal = atom->nlocal;
697 if (nlocal == nmax) {
698 grow(0);
699 if (atom->nextra_store)
700 memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra");
701 }
702
703 int m = 1;
704 x[nlocal][0] = buf[m++];
705 x[nlocal][1] = buf[m++];
706 x[nlocal][2] = buf[m++];
707 tag[nlocal] = (int) ubuf(buf[m++]).i;
708 type[nlocal] = (int) ubuf(buf[m++]).i;
709 mask[nlocal] = (int) ubuf(buf[m++]).i;
710 image[nlocal] = (tagint) ubuf(buf[m++]).i;
711 v[nlocal][0] = buf[m++];
712 v[nlocal][1] = buf[m++];
713 v[nlocal][2] = buf[m++];
714
715 molecule[nlocal] = (int) ubuf(buf[m++]).i;
716
717 num_bond[nlocal] = (int) ubuf(buf[m++]).i;
718 for (k = 0; k < num_bond[nlocal]; k++) {
719 bond_type[nlocal][k] = (int) ubuf(buf[m++]).i;
720 bond_atom[nlocal][k] = (int) ubuf(buf[m++]).i;
721 }
722
723 nspecial[nlocal][0] = nspecial[nlocal][1] = nspecial[nlocal][2] = 0;
724
725 double **extra = atom->extra;
726 if (atom->nextra_store) {
727 int size = static_cast<int> (buf[0]) - m;
728 for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
729 }
730
731 atom->nlocal++;
732 return m;
733 }
734
735 /* ----------------------------------------------------------------------
736 create one atom of itype at coord
737 set other values to defaults
738 ------------------------------------------------------------------------- */
739
create_atom(int itype,double * coord)740 void AtomVecBond::create_atom(int itype, double *coord)
741 {
742 int nlocal = atom->nlocal;
743 if (nlocal == nmax) grow(0);
744
745 tag[nlocal] = 0;
746 type[nlocal] = itype;
747 x[nlocal][0] = coord[0];
748 x[nlocal][1] = coord[1];
749 x[nlocal][2] = coord[2];
750 mask[nlocal] = 1;
751 image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
752 ((tagint) IMGMAX << IMGBITS) | IMGMAX;
753 v[nlocal][0] = 0.0;
754 v[nlocal][1] = 0.0;
755 v[nlocal][2] = 0.0;
756
757 molecule[nlocal] = 0;
758 num_bond[nlocal] = 0;
759 nspecial[nlocal][0] = nspecial[nlocal][1] = nspecial[nlocal][2] = 0;
760
761 atom->nlocal++;
762 }
763
764 /* ----------------------------------------------------------------------
765 unpack one line from Atoms section of data file
766 initialize other atom quantities
767 ------------------------------------------------------------------------- */
768
data_atom(double * coord,tagint imagetmp,char ** values)769 void AtomVecBond::data_atom(double *coord, tagint imagetmp, char **values)
770 {
771 int nlocal = atom->nlocal;
772 if (nlocal == nmax) grow(0);
773
774 tag[nlocal] = atoi(values[0]);
775 if (tag[nlocal] <= 0)
776 error->one(FLERR,"Invalid atom ID in Atoms section of data file");
777
778 molecule[nlocal] = atoi(values[1]);
779
780 type[nlocal] = atoi(values[2]);
781 if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
782 error->one(FLERR,"Invalid atom type in Atoms section of data file");
783
784 x[nlocal][0] = coord[0];
785 x[nlocal][1] = coord[1];
786 x[nlocal][2] = coord[2];
787
788 image[nlocal] = imagetmp;
789
790 mask[nlocal] = 1;
791 v[nlocal][0] = 0.0;
792 v[nlocal][1] = 0.0;
793 v[nlocal][2] = 0.0;
794 num_bond[nlocal] = 0;
795
796 atom->nlocal++;
797 }
798
799 /* ----------------------------------------------------------------------
800 unpack hybrid quantities from one line in Atoms section of data file
801 initialize other atom quantities for this sub-style
802 ------------------------------------------------------------------------- */
803
data_atom_hybrid(int nlocal,char ** values)804 int AtomVecBond::data_atom_hybrid(int nlocal, char **values)
805 {
806 molecule[nlocal] = atoi(values[0]);
807
808 num_bond[nlocal] = 0;
809
810 return 1;
811 }
812
813 /* ----------------------------------------------------------------------
814 pack atom info for data file including 3 image flags
815 ------------------------------------------------------------------------- */
816
pack_data(double ** buf)817 void AtomVecBond::pack_data(double **buf)
818 {
819 int nlocal = atom->nlocal;
820 for (int i = 0; i < nlocal; i++) {
821 buf[i][0] = ubuf(tag[i]).d;
822 buf[i][1] = ubuf(molecule[i]).d;
823 buf[i][2] = ubuf(type[i]).d;
824 buf[i][3] = x[i][0];
825 buf[i][4] = x[i][1];
826 buf[i][5] = x[i][2];
827 buf[i][6] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
828 buf[i][7] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
829 buf[i][8] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
830 }
831 }
832
833 /* ----------------------------------------------------------------------
834 pack hybrid atom info for data file
835 ------------------------------------------------------------------------- */
836
pack_data_hybrid(int i,double * buf)837 int AtomVecBond::pack_data_hybrid(int i, double *buf)
838 {
839 buf[0] = ubuf(molecule[i]).d;
840 return 1;
841 }
842
843 /* ----------------------------------------------------------------------
844 write atom info to data file including 3 image flags
845 ------------------------------------------------------------------------- */
846
write_data(FILE * fp,int n,double ** buf)847 void AtomVecBond::write_data(FILE *fp, int n, double **buf)
848 {
849 for (int i = 0; i < n; i++)
850 fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %d %d %d\n",
851 (int) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
852 (int) ubuf(buf[i][2]).i,
853 buf[i][3],buf[i][4],buf[i][5],
854 (int) ubuf(buf[i][6]).i,(int) ubuf(buf[i][7]).i,
855 (int) ubuf(buf[i][8]).i);
856 }
857
858 /* ----------------------------------------------------------------------
859 write hybrid atom info to data file
860 ------------------------------------------------------------------------- */
861
write_data_hybrid(FILE * fp,double * buf)862 int AtomVecBond::write_data_hybrid(FILE *fp, double *buf)
863 {
864 fprintf(fp," %d",(int) ubuf(buf[0]).i);
865 return 1;
866 }
867
868 /* ----------------------------------------------------------------------
869 return # of bytes of allocated memory
870 ------------------------------------------------------------------------- */
871
memory_usage()872 bigint AtomVecBond::memory_usage()
873 {
874 bigint bytes = 0;
875
876 if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax);
877 if (atom->memcheck("type")) bytes += memory->usage(type,nmax);
878 if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax);
879 if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
880 if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
881 if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
882 if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
883
884 if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
885 if (atom->memcheck("nspecial")) bytes += memory->usage(nspecial,nmax,3);
886 if (atom->memcheck("special"))
887 bytes += memory->usage(special,nmax,atom->maxspecial);
888
889 if (atom->memcheck("num_bond")) bytes += memory->usage(num_bond,nmax);
890 if (atom->memcheck("bond_type"))
891 bytes += memory->usage(bond_type,nmax,atom->bond_per_atom);
892 if (atom->memcheck("bond_atom"))
893 bytes += memory->usage(bond_atom,nmax,atom->bond_per_atom);
894
895 return bytes;
896 }
897