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_vec_body.h"
16 #include "style_body.h" // IWYU pragma: keep
17
18 #include "atom.h"
19 #include "body.h"
20 #include "error.h"
21 #include "fix.h"
22 #include "memory.h"
23 #include "modify.h"
24 #include "my_pool_chunk.h"
25
26 #include <cstring>
27
28 using namespace LAMMPS_NS;
29
30 /* ---------------------------------------------------------------------- */
31
AtomVecBody(LAMMPS * lmp)32 AtomVecBody::AtomVecBody(LAMMPS *lmp) : AtomVec(lmp)
33 {
34 molecular = Atom::ATOMIC;
35 bonus_flag = 1;
36
37 // first 3 sizes do not include values from body itself
38 // 1st,2nd body counts are added in process_args() via body style
39 // 3rd body count is added in size_restart_bonus()
40 // size_data_bonus is not used by Atom class for body style
41
42 size_forward_bonus = 4;
43 size_border_bonus = 10;
44 size_restart_bonus_one = 10;
45 size_data_bonus = 0;
46
47 atom->body_flag = 1;
48 atom->rmass_flag = 1;
49 atom->angmom_flag = atom->torque_flag = 1;
50 atom->radius_flag = 1;
51
52 nlocal_bonus = nghost_bonus = nmax_bonus = 0;
53 bonus = nullptr;
54
55 bptr = nullptr;
56
57 if (sizeof(double) == sizeof(int)) intdoubleratio = 1;
58 else if (sizeof(double) == 2*sizeof(int)) intdoubleratio = 2;
59 else error->all(FLERR,"Internal error in atom_style body");
60
61 // strings with peratom variables to include in each AtomVec method
62 // strings cannot contain fields in corresponding AtomVec default strings
63 // order of fields in a string does not matter
64 // except: fields_data_atom & fields_data_vel must match data file
65
66 fields_grow = (char *) "radius rmass angmom torque body";
67 fields_copy = (char *) "radius rmass angmom";
68 fields_comm = (char *) "";
69 fields_comm_vel = (char *) "angmom";
70 fields_reverse = (char *) "torque";
71 fields_border = (char *) "radius rmass";
72 fields_border_vel = (char *) "radius rmass angmom";
73 fields_exchange = (char *) "radius rmass angmom";
74 fields_restart = (char *) "radius rmass angmom";
75 fields_create = (char *) "radius rmass angmom body";
76 fields_data_atom = (char *) "id type body rmass x";
77 fields_data_vel = (char *) "id v angmom";
78 }
79
80 /* ---------------------------------------------------------------------- */
81
~AtomVecBody()82 AtomVecBody::~AtomVecBody()
83 {
84 int nall = nlocal_bonus + nghost_bonus;
85 for (int i = 0; i < nall; i++) {
86 icp->put(bonus[i].iindex);
87 dcp->put(bonus[i].dindex);
88 }
89 memory->sfree(bonus);
90
91 delete bptr;
92 }
93
94 /* ----------------------------------------------------------------------
95 process additional args
96 instantiate Body class
97 set size_forward and size_border to max sizes
98 ------------------------------------------------------------------------- */
99
process_args(int narg,char ** arg)100 void AtomVecBody::process_args(int narg, char **arg)
101 {
102 // suppress unused parameter warning dependent on style_body.h
103
104 (void)(arg);
105
106 if (narg < 1) error->all(FLERR,"Invalid atom_style body command");
107
108 if (0) {
109 bptr = nullptr;
110
111 #define BODY_CLASS
112 #define BodyStyle(key,Class) \
113 } else if (strcmp(arg[0],#key) == 0) { \
114 bptr = new Class(lmp,narg,arg);
115 #include "style_body.h" // IWYU pragma: keep
116 #undef BodyStyle
117 #undef BODY_CLASS
118
119 } else error->all(FLERR,utils::
120 check_packages_for_style("body",arg[0],lmp).c_str());
121
122 bptr->avec = this;
123 icp = bptr->icp;
124 dcp = bptr->dcp;
125
126 // max size of forward/border and exchange comm
127 // bptr values = max number of additional ivalues/dvalues from Body class
128
129 size_forward_bonus += bptr->size_forward;
130 size_border_bonus += bptr->size_border;
131 maxexchange = bptr->maxexchange;
132
133 setup_fields();
134 }
135
136 /* ----------------------------------------------------------------------
137 set local copies of all grow ptrs used by this class, except defaults
138 needed in replicate when 2 atom classes exist and it calls pack_restart()
139 ------------------------------------------------------------------------- */
140
grow_pointers()141 void AtomVecBody::grow_pointers()
142 {
143 body = atom->body;
144 rmass = atom->rmass;
145 radius = atom->radius;
146 angmom = atom->angmom;
147 }
148
149 /* ----------------------------------------------------------------------
150 grow bonus data structure
151 ------------------------------------------------------------------------- */
152
grow_bonus()153 void AtomVecBody::grow_bonus()
154 {
155 nmax_bonus = grow_nmax_bonus(nmax_bonus);
156 if (nmax_bonus < 0)
157 error->one(FLERR,"Per-processor system is too big");
158
159 bonus = (Bonus *) memory->srealloc(bonus,nmax_bonus*sizeof(Bonus),
160 "atom:bonus");
161 }
162
163 /* ----------------------------------------------------------------------
164 copy atom I info to atom J
165 if delflag and atom J has bonus data, then delete it
166 ------------------------------------------------------------------------- */
167
copy_bonus(int i,int j,int delflag)168 void AtomVecBody::copy_bonus(int i, int j, int delflag)
169 {
170 // if deleting atom J via delflag and J has bonus data, then delete it
171
172 if (delflag && body[j] >= 0) {
173 int k = body[j];
174 icp->put(bonus[k].iindex);
175 dcp->put(bonus[k].dindex);
176 copy_bonus_all(nlocal_bonus-1,k);
177 nlocal_bonus--;
178 }
179
180 // if atom I has bonus data, reset I's bonus.ilocal to loc J
181 // do NOT do this if self-copy (I=J) since I's bonus data is already deleted
182
183 if (body[i] >= 0 && i != j) bonus[body[i]].ilocal = j;
184 body[j] = body[i];
185 }
186
187 /* ----------------------------------------------------------------------
188 copy bonus data from I to J, effectively deleting the J entry
189 also reset body that points to I to now point to J
190 ------------------------------------------------------------------------- */
191
copy_bonus_all(int i,int j)192 void AtomVecBody::copy_bonus_all(int i, int j)
193 {
194 body[bonus[i].ilocal] = j;
195 memcpy(&bonus[j],&bonus[i],sizeof(Bonus));
196 }
197
198 /* ----------------------------------------------------------------------
199 clear ghost info in bonus data
200 called before ghosts are recommunicated in comm and irregular
201 ------------------------------------------------------------------------- */
202
clear_bonus()203 void AtomVecBody::clear_bonus()
204 {
205 int nall = nlocal_bonus + nghost_bonus;
206 for (int i = nlocal_bonus; i < nall; i++) {
207 icp->put(bonus[i].iindex);
208 dcp->put(bonus[i].dindex);
209 }
210 nghost_bonus = 0;
211
212 if (atom->nextra_grow)
213 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
214 modify->fix[atom->extra_grow[iextra]]->clear_bonus();
215 }
216
217 /* ---------------------------------------------------------------------- */
218
pack_comm_bonus(int n,int * list,double * buf)219 int AtomVecBody::pack_comm_bonus(int n, int *list, double *buf)
220 {
221 int i,j,m;
222 double *quat;
223
224 m = 0;
225 for (i = 0; i < n; i++) {
226 j = list[i];
227 if (body[j] >= 0) {
228 quat = bonus[body[j]].quat;
229 buf[m++] = quat[0];
230 buf[m++] = quat[1];
231 buf[m++] = quat[2];
232 buf[m++] = quat[3];
233 m += bptr->pack_comm_body(&bonus[body[j]],&buf[m]);
234 }
235 }
236
237 return m;
238 }
239
240 /* ---------------------------------------------------------------------- */
241
unpack_comm_bonus(int n,int first,double * buf)242 void AtomVecBody::unpack_comm_bonus(int n, int first, double *buf)
243 {
244 int i,m,last;
245 double *quat;
246
247 m = 0;
248 last = first + n;
249 for (i = first; i < last; i++) {
250 if (body[i] >= 0) {
251 quat = bonus[body[i]].quat;
252 quat[0] = buf[m++];
253 quat[1] = buf[m++];
254 quat[2] = buf[m++];
255 quat[3] = buf[m++];
256 m += bptr->unpack_comm_body(&bonus[body[i]],&buf[m]);
257 }
258 }
259 }
260
261 /* ---------------------------------------------------------------------- */
262
pack_border_bonus(int n,int * list,double * buf)263 int AtomVecBody::pack_border_bonus(int n, int *list, double *buf)
264 {
265 int i,j,m;
266 double *quat,*inertia;
267
268 m = 0;
269 for (i = 0; i < n; i++) {
270 j = list[i];
271 if (body[j] < 0) buf[m++] = ubuf(0).d;
272 else {
273 buf[m++] = ubuf(1).d;
274 quat = bonus[body[j]].quat;
275 inertia = bonus[body[j]].inertia;
276 buf[m++] = quat[0];
277 buf[m++] = quat[1];
278 buf[m++] = quat[2];
279 buf[m++] = quat[3];
280 buf[m++] = inertia[0];
281 buf[m++] = inertia[1];
282 buf[m++] = inertia[2];
283 buf[m++] = ubuf(bonus[body[j]].ninteger).d;
284 buf[m++] = ubuf(bonus[body[j]].ndouble).d;
285 m += bptr->pack_border_body(&bonus[body[j]],&buf[m]);
286 }
287 }
288
289 return m;
290 }
291
292 /* ---------------------------------------------------------------------- */
293
unpack_border_bonus(int n,int first,double * buf)294 int AtomVecBody::unpack_border_bonus(int n, int first, double *buf)
295 {
296 int i,j,m,last;
297 double *quat,*inertia;
298
299 m = 0;
300 last = first + n;
301 for (i = first; i < last; i++) {
302 body[i] = (int) ubuf(buf[m++]).i;
303 if (body[i] == 0) body[i] = -1;
304 else {
305 j = nlocal_bonus + nghost_bonus;
306 if (j == nmax_bonus) grow_bonus();
307 quat = bonus[j].quat;
308 inertia = bonus[j].inertia;
309 quat[0] = buf[m++];
310 quat[1] = buf[m++];
311 quat[2] = buf[m++];
312 quat[3] = buf[m++];
313 inertia[0] = buf[m++];
314 inertia[1] = buf[m++];
315 inertia[2] = buf[m++];
316 bonus[j].ninteger = (int) ubuf(buf[m++]).i;
317 bonus[j].ndouble = (int) ubuf(buf[m++]).i;
318 // corresponding put() calls are in clear_bonus()
319 bonus[j].ivalue = icp->get(bonus[j].ninteger,bonus[j].iindex);
320 bonus[j].dvalue = dcp->get(bonus[j].ndouble,bonus[j].dindex);
321 m += bptr->unpack_border_body(&bonus[j],&buf[m]);
322 bonus[j].ilocal = i;
323 body[i] = j;
324 nghost_bonus++;
325 }
326 }
327
328 return m;
329 }
330
331 /* ----------------------------------------------------------------------
332 pack data for atom I for sending to another proc
333 xyz must be 1st 3 values, so comm::exchange() can test on them
334 ------------------------------------------------------------------------- */
335
pack_exchange_bonus(int i,double * buf)336 int AtomVecBody::pack_exchange_bonus(int i, double *buf)
337 {
338 int m = 0;
339
340 if (body[i] < 0) buf[m++] = ubuf(0).d;
341 else {
342 buf[m++] = ubuf(1).d;
343 int j = body[i];
344 double *quat = bonus[j].quat;
345 double *inertia = bonus[j].inertia;
346 buf[m++] = quat[0];
347 buf[m++] = quat[1];
348 buf[m++] = quat[2];
349 buf[m++] = quat[3];
350 buf[m++] = inertia[0];
351 buf[m++] = inertia[1];
352 buf[m++] = inertia[2];
353 buf[m++] = ubuf(bonus[j].ninteger).d;
354 buf[m++] = ubuf(bonus[j].ndouble).d;
355 memcpy(&buf[m],bonus[j].ivalue,bonus[j].ninteger*sizeof(int));
356 if (intdoubleratio == 1) m += bonus[j].ninteger;
357 else m += (bonus[j].ninteger+1)/2;
358 memcpy(&buf[m],bonus[j].dvalue,bonus[j].ndouble*sizeof(double));
359 m += bonus[j].ndouble;
360 }
361
362 return m;
363 }
364
365 /* ---------------------------------------------------------------------- */
366
unpack_exchange_bonus(int ilocal,double * buf)367 int AtomVecBody::unpack_exchange_bonus(int ilocal, double *buf)
368 {
369 int m = 0;
370
371 body[ilocal] = (int) ubuf(buf[m++]).i;
372 if (body[ilocal] == 0) body[ilocal] = -1;
373 else {
374 if (nlocal_bonus == nmax_bonus) grow_bonus();
375 double *quat = bonus[nlocal_bonus].quat;
376 double *inertia = bonus[nlocal_bonus].inertia;
377 quat[0] = buf[m++];
378 quat[1] = buf[m++];
379 quat[2] = buf[m++];
380 quat[3] = buf[m++];
381 inertia[0] = buf[m++];
382 inertia[1] = buf[m++];
383 inertia[2] = buf[m++];
384 bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i;
385 bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i;
386 // corresponding put() calls are in copy()
387 bonus[nlocal_bonus].ivalue = icp->get(bonus[nlocal_bonus].ninteger,
388 bonus[nlocal_bonus].iindex);
389 bonus[nlocal_bonus].dvalue = dcp->get(bonus[nlocal_bonus].ndouble,
390 bonus[nlocal_bonus].dindex);
391 memcpy(bonus[nlocal_bonus].ivalue,&buf[m],
392 bonus[nlocal_bonus].ninteger*sizeof(int));
393 if (intdoubleratio == 1) m += bonus[nlocal_bonus].ninteger;
394 else m += (bonus[nlocal_bonus].ninteger+1)/2;
395 memcpy(bonus[nlocal_bonus].dvalue,&buf[m],
396 bonus[nlocal_bonus].ndouble*sizeof(double));
397 m += bonus[nlocal_bonus].ndouble;
398
399 bonus[nlocal_bonus].ilocal = ilocal;
400 body[ilocal] = nlocal_bonus++;
401 }
402
403 return m;
404 }
405
406 /* ----------------------------------------------------------------------
407 unpack data for one atom from restart file including bonus data
408 ------------------------------------------------------------------------- */
409
size_restart_bonus()410 int AtomVecBody::size_restart_bonus()
411 {
412 int i;
413
414 int n = 0;
415 int nlocal = atom->nlocal;
416 for (i = 0; i < nlocal; i++) {
417 if (body[i] >= 0) {
418 n += size_restart_bonus_one;
419 if (intdoubleratio == 1) n += bonus[body[i]].ninteger;
420 else n += (bonus[body[i]].ninteger+1)/2;
421 n += bonus[body[i]].ndouble;
422 } else n++;
423 }
424
425 return n;
426 }
427
428 /* ----------------------------------------------------------------------
429 pack atom I's data for restart file including extra quantities
430 xyz must be 1st 3 values, so that read_restart can test on them
431 molecular types may be negative, but write as positive
432 ------------------------------------------------------------------------- */
433
pack_restart_bonus(int i,double * buf)434 int AtomVecBody::pack_restart_bonus(int i, double *buf)
435 {
436 int m = 0;
437
438 if (body[i] < 0) buf[m++] = ubuf(0).d;
439 else {
440 buf[m++] = ubuf(1).d;
441 int j = body[i];
442 double *quat = bonus[j].quat;
443 double *inertia = bonus[j].inertia;
444 buf[m++] = quat[0];
445 buf[m++] = quat[1];
446 buf[m++] = quat[2];
447 buf[m++] = quat[3];
448 buf[m++] = inertia[0];
449 buf[m++] = inertia[1];
450 buf[m++] = inertia[2];
451 buf[m++] = ubuf(bonus[j].ninteger).d;
452 buf[m++] = ubuf(bonus[j].ndouble).d;
453 memcpy(&buf[m],bonus[j].ivalue,bonus[j].ninteger*sizeof(int));
454 if (intdoubleratio == 1) m += bonus[j].ninteger;
455 else m += (bonus[j].ninteger+1)/2;
456 memcpy(&buf[m],bonus[j].dvalue,bonus[j].ndouble*sizeof(double));
457 m += bonus[j].ndouble;
458 }
459
460 return m;
461 }
462
463 /* ----------------------------------------------------------------------
464 unpack data for one atom from restart file including bonus data
465 ------------------------------------------------------------------------- */
466
unpack_restart_bonus(int ilocal,double * buf)467 int AtomVecBody::unpack_restart_bonus(int ilocal, double *buf)
468 {
469 int m = 0;
470
471 body[ilocal] = (int) ubuf(buf[m++]).i;
472 if (body[ilocal] == 0) body[ilocal] = -1;
473 else {
474 if (nlocal_bonus == nmax_bonus) grow_bonus();
475 double *quat = bonus[nlocal_bonus].quat;
476 double *inertia = bonus[nlocal_bonus].inertia;
477 quat[0] = buf[m++];
478 quat[1] = buf[m++];
479 quat[2] = buf[m++];
480 quat[3] = buf[m++];
481 inertia[0] = buf[m++];
482 inertia[1] = buf[m++];
483 inertia[2] = buf[m++];
484 bonus[nlocal_bonus].ninteger = (int) ubuf(buf[m++]).i;
485 bonus[nlocal_bonus].ndouble = (int) ubuf(buf[m++]).i;
486 bonus[nlocal_bonus].ivalue = icp->get(bonus[nlocal_bonus].ninteger,
487 bonus[nlocal_bonus].iindex);
488 bonus[nlocal_bonus].dvalue = dcp->get(bonus[nlocal_bonus].ndouble,
489 bonus[nlocal_bonus].dindex);
490 memcpy(bonus[nlocal_bonus].ivalue,&buf[m],
491 bonus[nlocal_bonus].ninteger*sizeof(int));
492 if (intdoubleratio == 1) m += bonus[nlocal_bonus].ninteger;
493 else m += (bonus[nlocal_bonus].ninteger+1)/2;
494 memcpy(bonus[nlocal_bonus].dvalue,&buf[m],
495 bonus[nlocal_bonus].ndouble*sizeof(double));
496 m += bonus[nlocal_bonus].ndouble;
497 bonus[nlocal_bonus].ilocal = ilocal;
498 body[ilocal] = nlocal_bonus++;
499 }
500
501 return m;
502 }
503
504 /* ----------------------------------------------------------------------
505 create one atom of itype at coord
506 set other values to defaults
507 ------------------------------------------------------------------------- */
508
create_atom_post(int ilocal)509 void AtomVecBody::create_atom_post(int ilocal)
510 {
511 radius[ilocal] = 0.5;
512 rmass[ilocal] = 1.0;
513 body[ilocal] = -1;
514 }
515
516 /* ----------------------------------------------------------------------
517 modify what AtomVec::data_atom() just unpacked
518 or initialize other atom quantities
519 ------------------------------------------------------------------------- */
520
data_atom_post(int ilocal)521 void AtomVecBody::data_atom_post(int ilocal)
522 {
523 body_flag = body[ilocal];
524 if (body_flag == 0) body_flag = -1;
525 else if (body_flag == 1) body_flag = 0;
526 else error->one(FLERR,"Invalid body flag in Atoms section of data file");
527 body[ilocal] = body_flag;
528
529 if (rmass[ilocal] <= 0.0)
530 error->one(FLERR,"Invalid density in Atoms section of data file");
531
532 radius[ilocal] = 0.5;
533 angmom[ilocal][0] = 0.0;
534 angmom[ilocal][1] = 0.0;
535 angmom[ilocal][2] = 0.0;
536 }
537
538 /* ----------------------------------------------------------------------
539 unpack one body from Bodies section of data file
540 ------------------------------------------------------------------------- */
541
data_body(int m,int ninteger,int ndouble,int * ivalues,double * dvalues)542 void AtomVecBody::data_body(int m, int ninteger, int ndouble,
543 int *ivalues, double *dvalues)
544 {
545 if (body[m])
546 error->one(FLERR,"Assigning body parameters to non-body atom");
547 if (nlocal_bonus == nmax_bonus) grow_bonus();
548 bonus[nlocal_bonus].ilocal = m;
549 bptr->data_body(nlocal_bonus,ninteger,ndouble,ivalues,dvalues);
550 body[m] = nlocal_bonus++;
551 }
552
553 /* ----------------------------------------------------------------------
554 return # of bytes of allocated memory
555 ------------------------------------------------------------------------- */
556
memory_usage_bonus()557 double AtomVecBody::memory_usage_bonus()
558 {
559 double bytes = 0;
560 bytes += (double)nmax_bonus*sizeof(Bonus);
561 bytes += icp->size() + dcp->size();
562
563 int nall = nlocal_bonus + nghost_bonus;
564 for (int i = 0; i < nall; i++) {
565 if (body[i] >= 0) {
566 bytes += (double)bonus[body[i]].ninteger * sizeof(int);
567 bytes += (double)bonus[body[i]].ndouble * sizeof(double);
568 }
569 }
570
571 return bytes;
572 }
573
574 /* ----------------------------------------------------------------------
575 modify values for AtomVec::pack_data() to pack
576 ------------------------------------------------------------------------- */
577
pack_data_pre(int ilocal)578 void AtomVecBody::pack_data_pre(int ilocal)
579 {
580 body_flag = body[ilocal];
581
582 if (body_flag < 0) body[ilocal] = 0;
583 else body[ilocal] = 1;
584 }
585
586 /* ----------------------------------------------------------------------
587 pack bonus body info for writing to data file
588 if buf is nullptr, just return buffer size
589 ------------------------------------------------------------------------- */
590
pack_data_bonus(double * buf,int)591 int AtomVecBody::pack_data_bonus(double *buf, int /*flag*/)
592 {
593 int i;
594
595 tagint *tag = atom->tag;
596 int nlocal = atom->nlocal;
597
598 int m = 0;
599 for (i = 0; i < nlocal; i++) {
600 if (body[i] < 0) continue;
601 int n = bptr->pack_data_body(tag[i],body[i],buf);
602 m += n;
603 if (buf) buf += n;
604 }
605
606 return m;
607 }
608
609 /* ----------------------------------------------------------------------
610 write bonus body info to data file
611 ------------------------------------------------------------------------- */
612
write_data_bonus(FILE * fp,int n,double * buf,int)613 void AtomVecBody::write_data_bonus(FILE *fp, int n, double *buf, int /*flag*/)
614 {
615 int i = 0;
616 while (i < n) {
617 i += bptr->write_data_body(fp,&buf[i]);
618 }
619 }
620
621 /* ----------------------------------------------------------------------
622 unmodify values packed by AtomVec::pack_data()
623 ------------------------------------------------------------------------- */
624
pack_data_post(int ilocal)625 void AtomVecBody::pack_data_post(int ilocal)
626 {
627 body[ilocal] = body_flag;
628 }
629
630 /* ----------------------------------------------------------------------
631 body computes its size based on ivalues/dvalues and returns it
632 ------------------------------------------------------------------------- */
633
radius_body(int ninteger,int ndouble,int * ivalues,double * dvalues)634 double AtomVecBody::radius_body(int ninteger, int ndouble,
635 int *ivalues, double *dvalues)
636 {
637 return bptr->radius_body(ninteger,ndouble,ivalues,dvalues);
638 }
639
640 /* ----------------------------------------------------------------------
641 reset quat orientation for atom M to quat_external
642 called by Atom:add_molecule_atom()
643 ------------------------------------------------------------------------- */
644
set_quat(int m,double * quat_external)645 void AtomVecBody::set_quat(int m, double *quat_external)
646 {
647 if (body[m] < 0) error->one(FLERR,"Assigning quat to non-body atom");
648 double *quat = bonus[body[m]].quat;
649 quat[0] = quat_external[0]; quat[1] = quat_external[1];
650 quat[2] = quat_external[2]; quat[3] = quat_external[3];
651 }
652
653 /* ----------------------------------------------------------------------
654 debug method for sanity checking of own/bonus data pointers
655 ------------------------------------------------------------------------- */
656
657 /*
658 void AtomVecBody::check(int flag)
659 {
660 for (int i = 0; i < atom->nlocal; i++) {
661 if (body[i] >= 0 && body[i] >= nlocal_bonus) {
662 printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
663 errorx->one(FLERR,"BAD AAA");
664 }
665 }
666 for (int i = atom->nlocal; i < atom->nlocal+atom->nghost; i++) {
667 if (body[i] >= 0 &&
668 (body[i] < nlocal_bonus ||
669 body[i] >= nlocal_bonus+nghost_bonus)) {
670 printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
671 errorx->one(FLERR,"BAD BBB");
672 }
673 }
674 for (int i = 0; i < nlocal_bonus; i++) {
675 if (bonus[i].ilocal < 0 || bonus[i].ilocal >= atom->nlocal) {
676 printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
677 errorx->one(FLERR,"BAD CCC");
678 }
679 }
680 for (int i = 0; i < nlocal_bonus; i++) {
681 if (body[bonus[i].ilocal] != i) {
682 printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
683 errorx->one(FLERR,"BAD DDD");
684 }
685 }
686 for (int i = nlocal_bonus; i < nlocal_bonus+nghost_bonus; i++) {
687 if (bonus[i].ilocal < atom->nlocal ||
688 bonus[i].ilocal >= atom->nlocal+atom->nghost) {
689 printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
690 errorx->one(FLERR,"BAD EEE");
691 }
692 }
693 for (int i = nlocal_bonus; i < nlocal_bonus+nghost_bonus; i++) {
694 if (body[bonus[i].ilocal] != i) {
695 printf("Proc %d, step %ld, flag %d\n",comm->me,update->ntimestep,flag);
696 errorx->one(FLERR,"BAD FFF");
697 }
698 }
699 }
700 */
701