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