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 /* ----------------------------------------------------------------------
47 Contributing author: Mike Brown (SNL)
48 ------------------------------------------------------------------------- */
49
50 #include <stdlib.h>
51 #include "atom_vec_ellipsoid.h"
52 #include "math_extra.h"
53 #include "atom.h"
54 #include "comm.h"
55 #include "force.h"
56 #include "domain.h"
57 #include "modify.h"
58 #include "fix.h"
59 #include "math_const.h"
60 #include "memory.h"
61 #include "error.h"
62
63 using namespace LAMMPS_NS;
64 using namespace MathConst;
65
66 #define DELTA 10000
67 #define DELTA_BONUS 10000
68
69 /* ---------------------------------------------------------------------- */
70
AtomVecEllipsoid(LAMMPS * lmp)71 AtomVecEllipsoid::AtomVecEllipsoid(LAMMPS *lmp) : AtomVec(lmp)
72 {
73 molecular = 0;
74
75 comm_x_only = comm_f_only = 0;
76 size_forward = 7;
77 size_reverse = 6;
78 size_border = 14;
79 size_velocity = 6;
80 size_data_atom = 7;
81 size_data_vel = 7;
82 size_data_bonus = 8;
83 xcol_data = 5;
84
85 atom->ellipsoid_flag = 1;
86 atom->rmass_flag = atom->angmom_flag = atom->torque_flag = 1;
87
88 nlocal_bonus = nghost_bonus = nmax_bonus = 0;
89 bonus = NULL;
90 }
91
92 /* ---------------------------------------------------------------------- */
93
~AtomVecEllipsoid()94 AtomVecEllipsoid::~AtomVecEllipsoid()
95 {
96 memory->sfree(bonus);
97 }
98
99 /* ----------------------------------------------------------------------
100 grow atom arrays
101 n = 0 grows arrays by DELTA
102 n > 0 allocates arrays to size n
103 ------------------------------------------------------------------------- */
104
grow(int n)105 void AtomVecEllipsoid::grow(int n)
106 {
107 if (n == 0) nmax += DELTA;
108 else nmax = n;
109 atom->nmax = nmax;
110 if (nmax < 0 || nmax > MAXSMALLINT)
111 error->one(FLERR,"Per-processor system is too big");
112
113 tag = memory->grow(atom->tag,nmax,"atom:tag");
114 type = memory->grow(atom->type,nmax,"atom:type");
115 mask = memory->grow(atom->mask,nmax,"atom:mask");
116 image = memory->grow(atom->image,nmax,"atom:image");
117 x = memory->grow(atom->x,nmax,3,"atom:x");
118 v = memory->grow(atom->v,nmax,3,"atom:v");
119 f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
120
121 rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
122 angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom");
123 torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
124 ellipsoid = memory->grow(atom->ellipsoid,nmax,"atom:ellipsoid");
125
126 if (atom->nextra_grow)
127 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
128 modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax);
129 }
130
131 /* ----------------------------------------------------------------------
132 reset local array ptrs
133 ------------------------------------------------------------------------- */
134
grow_reset()135 void AtomVecEllipsoid::grow_reset()
136 {
137 tag = atom->tag; type = atom->type;
138 mask = atom->mask; image = atom->image;
139 x = atom->x; v = atom->v; f = atom->f;
140 rmass = atom->rmass; angmom = atom->angmom; torque = atom->torque;
141 ellipsoid = atom->ellipsoid;
142 }
143
144 /* ----------------------------------------------------------------------
145 grow bonus data structure
146 ------------------------------------------------------------------------- */
147
grow_bonus()148 void AtomVecEllipsoid::grow_bonus()
149 {
150 nmax_bonus += DELTA_BONUS;
151 if (nmax_bonus < 0 || nmax_bonus > MAXSMALLINT)
152 error->one(FLERR,"Per-processor system is too big");
153
154 bonus = (Bonus *) memory->srealloc(bonus,nmax_bonus*sizeof(Bonus),
155 "atom:bonus");
156 }
157
158 /* ----------------------------------------------------------------------
159 copy atom I info to atom J
160 ------------------------------------------------------------------------- */
161
copy(int i,int j,int delflag)162 void AtomVecEllipsoid::copy(int i, int j, int delflag)
163 {
164 tag[j] = tag[i];
165 type[j] = type[i];
166 mask[j] = mask[i];
167 image[j] = image[i];
168 x[j][0] = x[i][0];
169 x[j][1] = x[i][1];
170 x[j][2] = x[i][2];
171 v[j][0] = v[i][0];
172 v[j][1] = v[i][1];
173 v[j][2] = v[i][2];
174
175 rmass[j] = rmass[i];
176 angmom[j][0] = angmom[i][0];
177 angmom[j][1] = angmom[i][1];
178 angmom[j][2] = angmom[i][2];
179
180 // if deleting atom J via delflag and J has bonus data, then delete it
181
182 if (delflag && ellipsoid[j] >= 0) {
183 copy_bonus(nlocal_bonus-1,ellipsoid[j]);
184 nlocal_bonus--;
185 }
186
187 // if atom I has bonus data, reset I's bonus.ilocal to loc J
188 // do NOT do this if self-copy (I=J) since I's bonus data is already deleted
189
190 if (ellipsoid[i] >= 0 && i != j) bonus[ellipsoid[i]].ilocal = j;
191 ellipsoid[j] = ellipsoid[i];
192
193 if (atom->nextra_grow)
194 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
195 modify->fix[atom->extra_grow[iextra]]->copy_arrays(i,j,delflag);
196 }
197
198 /* ----------------------------------------------------------------------
199 copy bonus data from I to J, effectively deleting the J entry
200 also reset ellipsoid that points to I to now point to J
201 ------------------------------------------------------------------------- */
202
copy_bonus(int i,int j)203 void AtomVecEllipsoid::copy_bonus(int i, int j)
204 {
205 ellipsoid[bonus[i].ilocal] = j;
206 memcpy(&bonus[j],&bonus[i],sizeof(Bonus));
207 }
208
209 /* ----------------------------------------------------------------------
210 clear ghost info in bonus data
211 called before ghosts are recommunicated in comm and irregular
212 ------------------------------------------------------------------------- */
213
clear_bonus()214 void AtomVecEllipsoid::clear_bonus()
215 {
216 nghost_bonus = 0;
217 }
218
219 /* ----------------------------------------------------------------------
220 set shape values in bonus data for particle I
221 oriented aligned with xyz axes
222 this may create or delete entry in bonus data
223 ------------------------------------------------------------------------- */
224
set_shape(int i,double shapex,double shapey,double shapez)225 void AtomVecEllipsoid::set_shape(int i,
226 double shapex, double shapey, double shapez)
227 {
228 if (ellipsoid[i] < 0) {
229 if (shapex == 0.0 && shapey == 0.0 && shapez == 0.0) return;
230 if (nlocal_bonus == nmax_bonus) grow_bonus();
231 double *shape = bonus[nlocal_bonus].shape;
232 double *quat = bonus[nlocal_bonus].quat;
233 shape[0] = shapex;
234 shape[1] = shapey;
235 shape[2] = shapez;
236 quat[0] = 1.0;
237 quat[1] = 0.0;
238 quat[2] = 0.0;
239 quat[3] = 0.0;
240 bonus[nlocal_bonus].ilocal = i;
241 ellipsoid[i] = nlocal_bonus++;
242 } else if (shapex == 0.0 && shapey == 0.0 && shapez == 0.0) {
243 copy_bonus(nlocal_bonus-1,ellipsoid[i]);
244 nlocal_bonus--;
245 ellipsoid[i] = -1;
246 } else {
247 double *shape = bonus[ellipsoid[i]].shape;
248 shape[0] = shapex;
249 shape[1] = shapey;
250 shape[2] = shapez;
251 }
252 }
253
254 /* ---------------------------------------------------------------------- */
255
pack_comm(int n,int * list,double * buf,int pbc_flag,int * pbc)256 int AtomVecEllipsoid::pack_comm(int n, int *list, double *buf,
257 int pbc_flag, int *pbc)
258 {
259 int i,j,m;
260 double dx,dy,dz;
261 double *quat;
262
263 m = 0;
264 if (pbc_flag == 0) {
265 for (i = 0; i < n; i++) {
266 j = list[i];
267 buf[m++] = x[j][0];
268 buf[m++] = x[j][1];
269 buf[m++] = x[j][2];
270 if (ellipsoid[j] >= 0) {
271 quat = bonus[ellipsoid[j]].quat;
272 buf[m++] = quat[0];
273 buf[m++] = quat[1];
274 buf[m++] = quat[2];
275 buf[m++] = quat[3];
276 }
277 }
278 } else {
279 if (domain->triclinic == 0) {
280 dx = pbc[0]*domain->xprd;
281 dy = pbc[1]*domain->yprd;
282 dz = pbc[2]*domain->zprd;
283 } else {
284 dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
285 dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
286 dz = pbc[2]*domain->zprd;
287 }
288 for (i = 0; i < n; i++) {
289 j = list[i];
290 buf[m++] = x[j][0] + dx;
291 buf[m++] = x[j][1] + dy;
292 buf[m++] = x[j][2] + dz;
293 if (ellipsoid[j] >= 0) {
294 quat = bonus[ellipsoid[j]].quat;
295 buf[m++] = quat[0];
296 buf[m++] = quat[1];
297 buf[m++] = quat[2];
298 buf[m++] = quat[3];
299 }
300 }
301 }
302 return m;
303 }
304
305 /* ---------------------------------------------------------------------- */
306
pack_comm_vel(int n,int * list,double * buf,int pbc_flag,int * pbc)307 int AtomVecEllipsoid::pack_comm_vel(int n, int *list, double *buf,
308 int pbc_flag, int *pbc)
309 {
310 int i,j,m;
311 double dx,dy,dz,dvx,dvy,dvz;
312 double *quat;
313
314 m = 0;
315 if (pbc_flag == 0) {
316 for (i = 0; i < n; i++) {
317 j = list[i];
318 buf[m++] = x[j][0];
319 buf[m++] = x[j][1];
320 buf[m++] = x[j][2];
321 if (ellipsoid[j] >= 0) {
322 quat = bonus[ellipsoid[j]].quat;
323 buf[m++] = quat[0];
324 buf[m++] = quat[1];
325 buf[m++] = quat[2];
326 buf[m++] = quat[3];
327 }
328 buf[m++] = v[j][0];
329 buf[m++] = v[j][1];
330 buf[m++] = v[j][2];
331 buf[m++] = angmom[j][0];
332 buf[m++] = angmom[j][1];
333 buf[m++] = angmom[j][2];
334 }
335 } else {
336 if (domain->triclinic == 0) {
337 dx = pbc[0]*domain->xprd;
338 dy = pbc[1]*domain->yprd;
339 dz = pbc[2]*domain->zprd;
340 } else {
341 dx = pbc[0]*domain->xprd + pbc[5]*domain->xy + pbc[4]*domain->xz;
342 dy = pbc[1]*domain->yprd + pbc[3]*domain->yz;
343 dz = pbc[2]*domain->zprd;
344 }
345 if (!deform_vremap) {
346 for (i = 0; i < n; i++) {
347 j = list[i];
348 buf[m++] = x[j][0] + dx;
349 buf[m++] = x[j][1] + dy;
350 buf[m++] = x[j][2] + dz;
351 if (ellipsoid[j] >= 0) {
352 quat = bonus[ellipsoid[j]].quat;
353 buf[m++] = quat[0];
354 buf[m++] = quat[1];
355 buf[m++] = quat[2];
356 buf[m++] = quat[3];
357 }
358 buf[m++] = v[j][0];
359 buf[m++] = v[j][1];
360 buf[m++] = v[j][2];
361 buf[m++] = angmom[j][0];
362 buf[m++] = angmom[j][1];
363 buf[m++] = angmom[j][2];
364 }
365 } else {
366 dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
367 dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
368 dvz = pbc[2]*h_rate[2];
369 for (i = 0; i < n; i++) {
370 j = list[i];
371 buf[m++] = x[j][0] + dx;
372 buf[m++] = x[j][1] + dy;
373 buf[m++] = x[j][2] + dz;
374 if (ellipsoid[j] >= 0) {
375 quat = bonus[ellipsoid[j]].quat;
376 buf[m++] = quat[0];
377 buf[m++] = quat[1];
378 buf[m++] = quat[2];
379 buf[m++] = quat[3];
380 }
381 if (mask[i] & deform_groupbit) {
382 buf[m++] = v[j][0] + dvx;
383 buf[m++] = v[j][1] + dvy;
384 buf[m++] = v[j][2] + dvz;
385 } else {
386 buf[m++] = v[j][0];
387 buf[m++] = v[j][1];
388 buf[m++] = v[j][2];
389 }
390 buf[m++] = angmom[j][0];
391 buf[m++] = angmom[j][1];
392 buf[m++] = angmom[j][2];
393 }
394 }
395 }
396 return m;
397 }
398
399 /* ---------------------------------------------------------------------- */
400
pack_comm_hybrid(int n,int * list,double * buf)401 int AtomVecEllipsoid::pack_comm_hybrid(int n, int *list, double *buf)
402 {
403 int i,j,m;
404 double *quat;
405
406 m = 0;
407 for (i = 0; i < n; i++) {
408 j = list[i];
409 if (ellipsoid[j] >= 0) {
410 quat = bonus[ellipsoid[j]].quat;
411 buf[m++] = quat[0];
412 buf[m++] = quat[1];
413 buf[m++] = quat[2];
414 buf[m++] = quat[3];
415 }
416 }
417 return m;
418 }
419
420 /* ---------------------------------------------------------------------- */
421
unpack_comm(int n,int first,double * buf)422 void AtomVecEllipsoid::unpack_comm(int n, int first, double *buf)
423 {
424 int i,m,last;
425 double *quat;
426
427 m = 0;
428 last = first + n;
429 for (i = first; i < last; i++) {
430 x[i][0] = buf[m++];
431 x[i][1] = buf[m++];
432 x[i][2] = buf[m++];
433 if (ellipsoid[i] >= 0) {
434 quat = bonus[ellipsoid[i]].quat;
435 quat[0] = buf[m++];
436 quat[1] = buf[m++];
437 quat[2] = buf[m++];
438 quat[3] = buf[m++];
439 }
440 }
441 }
442
443 /* ---------------------------------------------------------------------- */
444
unpack_comm_vel(int n,int first,double * buf)445 void AtomVecEllipsoid::unpack_comm_vel(int n, int first, double *buf)
446 {
447 int i,m,last;
448 double *quat;
449
450 m = 0;
451 last = first + n;
452 for (i = first; i < last; i++) {
453 x[i][0] = buf[m++];
454 x[i][1] = buf[m++];
455 x[i][2] = buf[m++];
456 if (ellipsoid[i] >= 0) {
457 quat = bonus[ellipsoid[i]].quat;
458 quat[0] = buf[m++];
459 quat[1] = buf[m++];
460 quat[2] = buf[m++];
461 quat[3] = buf[m++];
462 }
463 v[i][0] = buf[m++];
464 v[i][1] = buf[m++];
465 v[i][2] = buf[m++];
466 angmom[i][0] = buf[m++];
467 angmom[i][1] = buf[m++];
468 angmom[i][2] = buf[m++];
469 }
470 }
471
472 /* ---------------------------------------------------------------------- */
473
unpack_comm_hybrid(int n,int first,double * buf)474 int AtomVecEllipsoid::unpack_comm_hybrid(int n, int first, double *buf)
475 {
476 int i,m,last;
477 double *quat;
478
479 m = 0;
480 last = first + n;
481 for (i = first; i < last; i++) {
482 if (ellipsoid[i] >= 0) {
483 quat = bonus[ellipsoid[i]].quat;
484 quat[0] = buf[m++];
485 quat[1] = buf[m++];
486 quat[2] = buf[m++];
487 quat[3] = buf[m++];
488 }
489 }
490 return m;
491 }
492
493 /* ---------------------------------------------------------------------- */
494
pack_reverse(int n,int first,double * buf)495 int AtomVecEllipsoid::pack_reverse(int n, int first, double *buf)
496 {
497 int i,m,last;
498
499 m = 0;
500 last = first + n;
501 for (i = first; i < last; i++) {
502 buf[m++] = f[i][0];
503 buf[m++] = f[i][1];
504 buf[m++] = f[i][2];
505 buf[m++] = torque[i][0];
506 buf[m++] = torque[i][1];
507 buf[m++] = torque[i][2];
508 }
509 return m;
510 }
511
512 /* ---------------------------------------------------------------------- */
513
pack_reverse_hybrid(int n,int first,double * buf)514 int AtomVecEllipsoid::pack_reverse_hybrid(int n, int first, double *buf)
515 {
516 int i,m,last;
517
518 m = 0;
519 last = first + n;
520 for (i = first; i < last; i++) {
521 buf[m++] = torque[i][0];
522 buf[m++] = torque[i][1];
523 buf[m++] = torque[i][2];
524 }
525 return m;
526 }
527
528 /* ---------------------------------------------------------------------- */
529
unpack_reverse(int n,int * list,double * buf)530 void AtomVecEllipsoid::unpack_reverse(int n, int *list, double *buf)
531 {
532 int i,j,m;
533
534 m = 0;
535 for (i = 0; i < n; i++) {
536 j = list[i];
537 f[j][0] += buf[m++];
538 f[j][1] += buf[m++];
539 f[j][2] += buf[m++];
540 torque[j][0] += buf[m++];
541 torque[j][1] += buf[m++];
542 torque[j][2] += buf[m++];
543 }
544 }
545
546 /* ---------------------------------------------------------------------- */
547
unpack_reverse_hybrid(int n,int * list,double * buf)548 int AtomVecEllipsoid::unpack_reverse_hybrid(int n, int *list, double *buf)
549 {
550 int i,j,m;
551
552 m = 0;
553 for (i = 0; i < n; i++) {
554 j = list[i];
555 torque[j][0] += buf[m++];
556 torque[j][1] += buf[m++];
557 torque[j][2] += buf[m++];
558 }
559 return m;
560 }
561
562 /* ---------------------------------------------------------------------- */
563
pack_border(int n,int * list,double * buf,int pbc_flag,int * pbc)564 int AtomVecEllipsoid::pack_border(int n, int *list, double *buf,
565 int pbc_flag, int *pbc)
566 {
567 int i,j,m;
568 double dx,dy,dz;
569 double *shape,*quat;
570
571 m = 0;
572 if (pbc_flag == 0) {
573 for (i = 0; i < n; i++) {
574 j = list[i];
575 buf[m++] = x[j][0];
576 buf[m++] = x[j][1];
577 buf[m++] = x[j][2];
578 buf[m++] = ubuf(tag[j]).d;
579 buf[m++] = ubuf(type[j]).d;
580 buf[m++] = ubuf(mask[j]).d;
581 if (ellipsoid[j] < 0) buf[m++] = ubuf(0).d;
582 else {
583 buf[m++] = ubuf(1).d;
584 shape = bonus[ellipsoid[j]].shape;
585 quat = bonus[ellipsoid[j]].quat;
586 buf[m++] = shape[0];
587 buf[m++] = shape[1];
588 buf[m++] = shape[2];
589 buf[m++] = quat[0];
590 buf[m++] = quat[1];
591 buf[m++] = quat[2];
592 buf[m++] = quat[3];
593 }
594 }
595 } else {
596 if (domain->triclinic == 0) {
597 dx = pbc[0]*domain->xprd;
598 dy = pbc[1]*domain->yprd;
599 dz = pbc[2]*domain->zprd;
600 } else {
601 dx = pbc[0];
602 dy = pbc[1];
603 dz = pbc[2];
604 }
605 for (i = 0; i < n; i++) {
606 j = list[i];
607 buf[m++] = x[j][0] + dx;
608 buf[m++] = x[j][1] + dy;
609 buf[m++] = x[j][2] + dz;
610 buf[m++] = ubuf(tag[j]).d;
611 buf[m++] = ubuf(type[j]).d;
612 buf[m++] = ubuf(mask[j]).d;
613 if (ellipsoid[j] < 0) buf[m++] = ubuf(0).d;
614 else {
615 buf[m++] = ubuf(1).d;
616 shape = bonus[ellipsoid[j]].shape;
617 quat = bonus[ellipsoid[j]].quat;
618 buf[m++] = shape[0];
619 buf[m++] = shape[1];
620 buf[m++] = shape[2];
621 buf[m++] = quat[0];
622 buf[m++] = quat[1];
623 buf[m++] = quat[2];
624 buf[m++] = quat[3];
625 }
626 }
627 }
628
629 if (atom->nextra_border)
630 for (int iextra = 0; iextra < atom->nextra_border; iextra++)
631 m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
632
633 return m;
634 }
635
636 /* ---------------------------------------------------------------------- */
637
pack_border_vel(int n,int * list,double * buf,int pbc_flag,int * pbc)638 int AtomVecEllipsoid::pack_border_vel(int n, int *list, double *buf,
639 int pbc_flag, int *pbc)
640 {
641 int i,j,m;
642 double dx,dy,dz,dvx,dvy,dvz;
643 double *shape,*quat;
644
645 m = 0;
646 if (pbc_flag == 0) {
647 for (i = 0; i < n; i++) {
648 j = list[i];
649 buf[m++] = x[j][0];
650 buf[m++] = x[j][1];
651 buf[m++] = x[j][2];
652 buf[m++] = ubuf(tag[j]).d;
653 buf[m++] = ubuf(type[j]).d;
654 buf[m++] = ubuf(mask[j]).d;
655 if (ellipsoid[j] < 0) buf[m++] = ubuf(0).d;
656 else {
657 buf[m++] = ubuf(1).d;
658 shape = bonus[ellipsoid[j]].shape;
659 quat = bonus[ellipsoid[j]].quat;
660 buf[m++] = shape[0];
661 buf[m++] = shape[1];
662 buf[m++] = shape[2];
663 buf[m++] = quat[0];
664 buf[m++] = quat[1];
665 buf[m++] = quat[2];
666 buf[m++] = quat[3];
667 }
668 buf[m++] = v[j][0];
669 buf[m++] = v[j][1];
670 buf[m++] = v[j][2];
671 buf[m++] = angmom[j][0];
672 buf[m++] = angmom[j][1];
673 buf[m++] = angmom[j][2];
674 }
675 } else {
676 if (domain->triclinic == 0) {
677 dx = pbc[0]*domain->xprd;
678 dy = pbc[1]*domain->yprd;
679 dz = pbc[2]*domain->zprd;
680 } else {
681 dx = pbc[0];
682 dy = pbc[1];
683 dz = pbc[2];
684 }
685 if (!deform_vremap) {
686 for (i = 0; i < n; i++) {
687 j = list[i];
688 buf[m++] = x[j][0] + dx;
689 buf[m++] = x[j][1] + dy;
690 buf[m++] = x[j][2] + dz;
691 buf[m++] = ubuf(tag[j]).d;
692 buf[m++] = ubuf(type[j]).d;
693 buf[m++] = ubuf(mask[j]).d;
694 if (ellipsoid[j] < 0) buf[m++] = ubuf(0).d;
695 else {
696 buf[m++] = ubuf(1).d;
697 shape = bonus[ellipsoid[j]].shape;
698 quat = bonus[ellipsoid[j]].quat;
699 buf[m++] = shape[0];
700 buf[m++] = shape[1];
701 buf[m++] = shape[2];
702 buf[m++] = quat[0];
703 buf[m++] = quat[1];
704 buf[m++] = quat[2];
705 buf[m++] = quat[3];
706 }
707 buf[m++] = v[j][0];
708 buf[m++] = v[j][1];
709 buf[m++] = v[j][2];
710 buf[m++] = angmom[j][0];
711 buf[m++] = angmom[j][1];
712 buf[m++] = angmom[j][2];
713 }
714 } else {
715 dvx = pbc[0]*h_rate[0] + pbc[5]*h_rate[5] + pbc[4]*h_rate[4];
716 dvy = pbc[1]*h_rate[1] + pbc[3]*h_rate[3];
717 dvz = pbc[2]*h_rate[2];
718 for (i = 0; i < n; i++) {
719 j = list[i];
720 buf[m++] = x[j][0] + dx;
721 buf[m++] = x[j][1] + dy;
722 buf[m++] = x[j][2] + dz;
723 buf[m++] = ubuf(tag[j]).d;
724 buf[m++] = ubuf(type[j]).d;
725 buf[m++] = ubuf(mask[j]).d;
726 if (ellipsoid[j] < 0) buf[m++] = ubuf(0).d;
727 else {
728 buf[m++] = ubuf(1).d;
729 shape = bonus[ellipsoid[j]].shape;
730 quat = bonus[ellipsoid[j]].quat;
731 buf[m++] = shape[0];
732 buf[m++] = shape[1];
733 buf[m++] = shape[2];
734 buf[m++] = quat[0];
735 buf[m++] = quat[1];
736 buf[m++] = quat[2];
737 buf[m++] = quat[3];
738 }
739 if (mask[i] & deform_groupbit) {
740 buf[m++] = v[j][0] + dvx;
741 buf[m++] = v[j][1] + dvy;
742 buf[m++] = v[j][2] + dvz;
743 } else {
744 buf[m++] = v[j][0];
745 buf[m++] = v[j][1];
746 buf[m++] = v[j][2];
747 }
748 buf[m++] = angmom[j][0];
749 buf[m++] = angmom[j][1];
750 buf[m++] = angmom[j][2];
751 }
752 }
753 }
754
755 if (atom->nextra_border)
756 for (int iextra = 0; iextra < atom->nextra_border; iextra++)
757 m += modify->fix[atom->extra_border[iextra]]->pack_border(n,list,&buf[m]);
758
759 return m;
760 }
761
762 /* ---------------------------------------------------------------------- */
763
pack_border_hybrid(int n,int * list,double * buf)764 int AtomVecEllipsoid::pack_border_hybrid(int n, int *list, double *buf)
765 {
766 int i,j,m;
767 double *shape,*quat;
768
769 m = 0;
770 for (i = 0; i < n; i++) {
771 j = list[i];
772 if (ellipsoid[j] < 0) buf[m++] = ubuf(0).d;
773 else {
774 buf[m++] = ubuf(1).d;
775 shape = bonus[ellipsoid[j]].shape;
776 quat = bonus[ellipsoid[j]].quat;
777 buf[m++] = shape[0];
778 buf[m++] = shape[1];
779 buf[m++] = shape[2];
780 buf[m++] = quat[0];
781 buf[m++] = quat[1];
782 buf[m++] = quat[2];
783 buf[m++] = quat[3];
784 }
785 }
786 return m;
787 }
788
789 /* ---------------------------------------------------------------------- */
790
unpack_border(int n,int first,double * buf)791 void AtomVecEllipsoid::unpack_border(int n, int first, double *buf)
792 {
793 int i,j,m,last;
794 double *shape,*quat;
795
796 m = 0;
797 last = first + n;
798 for (i = first; i < last; i++) {
799 if (i == nmax) grow(0);
800 x[i][0] = buf[m++];
801 x[i][1] = buf[m++];
802 x[i][2] = buf[m++];
803 tag[i] = (int) ubuf(buf[m++]).i;
804 type[i] = (int) ubuf(buf[m++]).i;
805 mask[i] = (int) ubuf(buf[m++]).i;
806 ellipsoid[i] = (int) ubuf(buf[m++]).i;
807 if (ellipsoid[i] == 0) ellipsoid[i] = -1;
808 else {
809 j = nlocal_bonus + nghost_bonus;
810 if (j == nmax_bonus) grow_bonus();
811 shape = bonus[j].shape;
812 quat = bonus[j].quat;
813 shape[0] = buf[m++];
814 shape[1] = buf[m++];
815 shape[2] = buf[m++];
816 quat[0] = buf[m++];
817 quat[1] = buf[m++];
818 quat[2] = buf[m++];
819 quat[3] = buf[m++];
820 bonus[j].ilocal = i;
821 ellipsoid[i] = j;
822 nghost_bonus++;
823 }
824 }
825
826 if (atom->nextra_border)
827 for (int iextra = 0; iextra < atom->nextra_border; iextra++)
828 m += modify->fix[atom->extra_border[iextra]]->
829 unpack_border(n,first,&buf[m]);
830 }
831
832 /* ---------------------------------------------------------------------- */
833
unpack_border_vel(int n,int first,double * buf)834 void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf)
835 {
836 int i,j,m,last;
837 double *shape,*quat;
838
839 m = 0;
840 last = first + n;
841 for (i = first; i < last; i++) {
842 if (i == nmax) grow(0);
843 x[i][0] = buf[m++];
844 x[i][1] = buf[m++];
845 x[i][2] = buf[m++];
846 tag[i] = (int) ubuf(buf[m++]).i;
847 type[i] = (int) ubuf(buf[m++]).i;
848 mask[i] = (int) ubuf(buf[m++]).i;
849 ellipsoid[i] = (int) ubuf(buf[m++]).i;
850 if (ellipsoid[i] == 0) ellipsoid[i] = -1;
851 else {
852 j = nlocal_bonus + nghost_bonus;
853 if (j == nmax_bonus) grow_bonus();
854 shape = bonus[j].shape;
855 quat = bonus[j].quat;
856 shape[0] = buf[m++];
857 shape[1] = buf[m++];
858 shape[2] = buf[m++];
859 quat[0] = buf[m++];
860 quat[1] = buf[m++];
861 quat[2] = buf[m++];
862 quat[3] = buf[m++];
863 bonus[j].ilocal = i;
864 ellipsoid[i] = j;
865 nghost_bonus++;
866 }
867 v[i][0] = buf[m++];
868 v[i][1] = buf[m++];
869 v[i][2] = buf[m++];
870 angmom[i][0] = buf[m++];
871 angmom[i][1] = buf[m++];
872 angmom[i][2] = buf[m++];
873 }
874
875 if (atom->nextra_border)
876 for (int iextra = 0; iextra < atom->nextra_border; iextra++)
877 m += modify->fix[atom->extra_border[iextra]]->
878 unpack_border(n,first,&buf[m]);
879 }
880
881 /* ---------------------------------------------------------------------- */
882
unpack_border_hybrid(int n,int first,double * buf)883 int AtomVecEllipsoid::unpack_border_hybrid(int n, int first, double *buf)
884 {
885 int i,j,m,last;
886 double *shape,*quat;
887
888 m = 0;
889 last = first + n;
890 for (i = first; i < last; i++) {
891 ellipsoid[i] = (int) ubuf(buf[m++]).i;
892 if (ellipsoid[i] == 0) ellipsoid[i] = -1;
893 else {
894 j = nlocal_bonus + nghost_bonus;
895 if (j == nmax_bonus) grow_bonus();
896 shape = bonus[j].shape;
897 quat = bonus[j].quat;
898 shape[0] = buf[m++];
899 shape[1] = buf[m++];
900 shape[2] = buf[m++];
901 quat[0] = buf[m++];
902 quat[1] = buf[m++];
903 quat[2] = buf[m++];
904 quat[3] = buf[m++];
905 bonus[j].ilocal = i;
906 ellipsoid[i] = j;
907 nghost_bonus++;
908 }
909 }
910 return m;
911 }
912
913 /* ----------------------------------------------------------------------
914 pack data for atom I for sending to another proc
915 xyz must be 1st 3 values, so comm::exchange() can test on them
916 ------------------------------------------------------------------------- */
917
pack_exchange(int i,double * buf)918 int AtomVecEllipsoid::pack_exchange(int i, double *buf)
919 {
920 int m = 1;
921 buf[m++] = x[i][0];
922 buf[m++] = x[i][1];
923 buf[m++] = x[i][2];
924 buf[m++] = v[i][0];
925 buf[m++] = v[i][1];
926 buf[m++] = v[i][2];
927 buf[m++] = ubuf(tag[i]).d;
928 buf[m++] = ubuf(type[i]).d;
929 buf[m++] = ubuf(mask[i]).d;
930 buf[m++] = ubuf(image[i]).d;
931
932 buf[m++] = rmass[i];
933 buf[m++] = angmom[i][0];
934 buf[m++] = angmom[i][1];
935 buf[m++] = angmom[i][2];
936
937 if (ellipsoid[i] < 0) buf[m++] = ubuf(0).d;
938 else {
939 buf[m++] = ubuf(1).d;
940 int j = ellipsoid[i];
941 double *shape = bonus[j].shape;
942 double *quat = bonus[j].quat;
943 buf[m++] = shape[0];
944 buf[m++] = shape[1];
945 buf[m++] = shape[2];
946 buf[m++] = quat[0];
947 buf[m++] = quat[1];
948 buf[m++] = quat[2];
949 buf[m++] = quat[3];
950 }
951
952 if (atom->nextra_grow)
953 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
954 m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i,&buf[m]);
955
956 buf[0] = m;
957 return m;
958 }
959
960 /* ---------------------------------------------------------------------- */
961
unpack_exchange(double * buf)962 int AtomVecEllipsoid::unpack_exchange(double *buf)
963 {
964 int nlocal = atom->nlocal;
965 if (nlocal == nmax) grow(0);
966
967 int m = 1;
968 x[nlocal][0] = buf[m++];
969 x[nlocal][1] = buf[m++];
970 x[nlocal][2] = buf[m++];
971 v[nlocal][0] = buf[m++];
972 v[nlocal][1] = buf[m++];
973 v[nlocal][2] = buf[m++];
974 tag[nlocal] = (int) ubuf(buf[m++]).i;
975 type[nlocal] = (int) ubuf(buf[m++]).i;
976 mask[nlocal] = (int) ubuf(buf[m++]).i;
977 image[nlocal] = (tagint) ubuf(buf[m++]).i;
978
979 rmass[nlocal] = buf[m++];
980 angmom[nlocal][0] = buf[m++];
981 angmom[nlocal][1] = buf[m++];
982 angmom[nlocal][2] = buf[m++];
983
984 ellipsoid[nlocal] = (int) ubuf(buf[m++]).i;
985 if (ellipsoid[nlocal] == 0) ellipsoid[nlocal] = -1;
986 else {
987 if (nlocal_bonus == nmax_bonus) grow_bonus();
988 double *shape = bonus[nlocal_bonus].shape;
989 double *quat = bonus[nlocal_bonus].quat;
990 shape[0] = buf[m++];
991 shape[1] = buf[m++];
992 shape[2] = buf[m++];
993 quat[0] = buf[m++];
994 quat[1] = buf[m++];
995 quat[2] = buf[m++];
996 quat[3] = buf[m++];
997 bonus[nlocal_bonus].ilocal = nlocal;
998 ellipsoid[nlocal] = nlocal_bonus++;
999 }
1000
1001 if (atom->nextra_grow)
1002 for (int iextra = 0; iextra < atom->nextra_grow; iextra++)
1003 m += modify->fix[atom->extra_grow[iextra]]->
1004 unpack_exchange(nlocal,&buf[m]);
1005
1006 atom->nlocal++;
1007 return m;
1008 }
1009
1010 /* ----------------------------------------------------------------------
1011 size of restart data for all atoms owned by this proc
1012 include extra data stored by fixes
1013 ------------------------------------------------------------------------- */
1014
size_restart()1015 int AtomVecEllipsoid::size_restart()
1016 {
1017 int i;
1018
1019 int n = 0;
1020 int nlocal = atom->nlocal;
1021 for (i = 0; i < nlocal; i++)
1022 if (ellipsoid[i] >= 0) n += 23;
1023 else n += 16;
1024
1025 if (atom->nextra_restart)
1026 for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
1027 for (i = 0; i < nlocal; i++)
1028 n += modify->fix[atom->extra_restart[iextra]]->size_restart(i);
1029
1030 return n;
1031 }
1032
1033 /* ----------------------------------------------------------------------
1034 pack atom I's data for restart file including bonus data
1035 xyz must be 1st 3 values, so that read_restart can test on them
1036 molecular types may be negative, but write as positive
1037 ------------------------------------------------------------------------- */
1038
pack_restart(int i,double * buf)1039 int AtomVecEllipsoid::pack_restart(int i, double *buf)
1040 {
1041 int m = 1;
1042 buf[m++] = x[i][0];
1043 buf[m++] = x[i][1];
1044 buf[m++] = x[i][2];
1045 buf[m++] = ubuf(tag[i]).d;
1046 buf[m++] = ubuf(type[i]).d;
1047 buf[m++] = ubuf(mask[i]).d;
1048 buf[m++] = ubuf(image[i]).d;
1049 buf[m++] = v[i][0];
1050 buf[m++] = v[i][1];
1051 buf[m++] = v[i][2];
1052
1053 buf[m++] = rmass[i];
1054 buf[m++] = angmom[i][0];
1055 buf[m++] = angmom[i][1];
1056 buf[m++] = angmom[i][2];
1057
1058 if (ellipsoid[i] < 0) buf[m++] = ubuf(0).d;
1059 else {
1060 buf[m++] = ubuf(1).d;
1061 int j = ellipsoid[i];
1062 buf[m++] = bonus[j].shape[0];
1063 buf[m++] = bonus[j].shape[1];
1064 buf[m++] = bonus[j].shape[2];
1065 buf[m++] = bonus[j].quat[0];
1066 buf[m++] = bonus[j].quat[1];
1067 buf[m++] = bonus[j].quat[2];
1068 buf[m++] = bonus[j].quat[3];
1069 }
1070
1071 if (atom->nextra_restart)
1072 for (int iextra = 0; iextra < atom->nextra_restart; iextra++)
1073 m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i,&buf[m]);
1074
1075 buf[0] = m;
1076 return m;
1077 }
1078
1079 /* ----------------------------------------------------------------------
1080 unpack data for one atom from restart file including bonus data
1081 ------------------------------------------------------------------------- */
1082
unpack_restart(double * buf)1083 int AtomVecEllipsoid::unpack_restart(double *buf)
1084 {
1085 int nlocal = atom->nlocal;
1086 if (nlocal == nmax) {
1087 grow(0);
1088 if (atom->nextra_store)
1089 memory->grow(atom->extra,nmax,atom->nextra_store,"atom:extra");
1090 }
1091
1092 int m = 1;
1093 x[nlocal][0] = buf[m++];
1094 x[nlocal][1] = buf[m++];
1095 x[nlocal][2] = buf[m++];
1096 tag[nlocal] = (int) ubuf(buf[m++]).i;
1097 type[nlocal] = (int) ubuf(buf[m++]).i;
1098 mask[nlocal] = (int) ubuf(buf[m++]).i;
1099 image[nlocal] = (tagint) ubuf(buf[m++]).i;
1100 v[nlocal][0] = buf[m++];
1101 v[nlocal][1] = buf[m++];
1102 v[nlocal][2] = buf[m++];
1103
1104 rmass[nlocal] = buf[m++];
1105 angmom[nlocal][0] = buf[m++];
1106 angmom[nlocal][1] = buf[m++];
1107 angmom[nlocal][2] = buf[m++];
1108
1109 ellipsoid[nlocal] = (int) ubuf(buf[m++]).i;
1110 if (ellipsoid[nlocal] == 0) ellipsoid[nlocal] = -1;
1111 else {
1112 if (nlocal_bonus == nmax_bonus) grow_bonus();
1113 double *shape = bonus[nlocal_bonus].shape;
1114 double *quat = bonus[nlocal_bonus].quat;
1115 shape[0] = buf[m++];
1116 shape[1] = buf[m++];
1117 shape[2] = buf[m++];
1118 quat[0] = buf[m++];
1119 quat[1] = buf[m++];
1120 quat[2] = buf[m++];
1121 quat[3] = buf[m++];
1122 bonus[nlocal_bonus].ilocal = nlocal;
1123 ellipsoid[nlocal] = nlocal_bonus++;
1124 }
1125
1126 double **extra = atom->extra;
1127 if (atom->nextra_store) {
1128 int size = static_cast<int> (buf[0]) - m;
1129 for (int i = 0; i < size; i++) extra[nlocal][i] = buf[m++];
1130 }
1131
1132 atom->nlocal++;
1133 return m;
1134 }
1135
1136 /* ----------------------------------------------------------------------
1137 create one atom of itype at coord
1138 set other values to defaults
1139 ------------------------------------------------------------------------- */
1140
create_atom(int itype,double * coord)1141 void AtomVecEllipsoid::create_atom(int itype, double *coord)
1142 {
1143 int nlocal = atom->nlocal;
1144 if (nlocal == nmax) grow(0);
1145
1146 tag[nlocal] = 0;
1147 type[nlocal] = itype;
1148 x[nlocal][0] = coord[0];
1149 x[nlocal][1] = coord[1];
1150 x[nlocal][2] = coord[2];
1151 mask[nlocal] = 1;
1152 image[nlocal] = ((tagint) IMGMAX << IMG2BITS) |
1153 ((tagint) IMGMAX << IMGBITS) | IMGMAX;
1154 v[nlocal][0] = 0.0;
1155 v[nlocal][1] = 0.0;
1156 v[nlocal][2] = 0.0;
1157
1158 rmass[nlocal] = 1.0;
1159 angmom[nlocal][0] = 0.0;
1160 angmom[nlocal][1] = 0.0;
1161 angmom[nlocal][2] = 0.0;
1162 ellipsoid[nlocal] = -1;
1163
1164 atom->nlocal++;
1165 }
1166
1167 /* ----------------------------------------------------------------------
1168 unpack one line from Atoms section of data file
1169 initialize other atom quantities
1170 ------------------------------------------------------------------------- */
1171
data_atom(double * coord,tagint imagetmp,char ** values)1172 void AtomVecEllipsoid::data_atom(double *coord, tagint imagetmp, char **values)
1173 {
1174 int nlocal = atom->nlocal;
1175 if (nlocal == nmax) grow(0);
1176
1177 tag[nlocal] = atoi(values[0]);
1178 if (tag[nlocal] <= 0)
1179 error->one(FLERR,"Invalid atom ID in Atoms section of data file");
1180
1181 type[nlocal] = atoi(values[1]);
1182 if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
1183 error->one(FLERR,"Invalid atom type in Atoms section of data file");
1184
1185 ellipsoid[nlocal] = atoi(values[2]);
1186 if (ellipsoid[nlocal] == 0) ellipsoid[nlocal] = -1;
1187 else if (ellipsoid[nlocal] == 1) ellipsoid[nlocal] = 0;
1188 else error->one(FLERR,"Invalid atom type in Atoms section of data file");
1189
1190 rmass[nlocal] = atof(values[3]);
1191 if (rmass[nlocal] <= 0.0)
1192 error->one(FLERR,"Invalid density in Atoms section of data file");
1193
1194 x[nlocal][0] = coord[0];
1195 x[nlocal][1] = coord[1];
1196 x[nlocal][2] = coord[2];
1197
1198 image[nlocal] = imagetmp;
1199
1200 mask[nlocal] = 1;
1201 v[nlocal][0] = 0.0;
1202 v[nlocal][1] = 0.0;
1203 v[nlocal][2] = 0.0;
1204 angmom[nlocal][0] = 0.0;
1205 angmom[nlocal][1] = 0.0;
1206 angmom[nlocal][2] = 0.0;
1207
1208 atom->nlocal++;
1209 }
1210
1211 /* ----------------------------------------------------------------------
1212 unpack hybrid quantities from one line in Atoms section of data file
1213 initialize other atom quantities for this sub-style
1214 ------------------------------------------------------------------------- */
1215
data_atom_hybrid(int nlocal,char ** values)1216 int AtomVecEllipsoid::data_atom_hybrid(int nlocal, char **values)
1217 {
1218 ellipsoid[nlocal] = atoi(values[0]);
1219 if (ellipsoid[nlocal] == 0) ellipsoid[nlocal] = -1;
1220 else if (ellipsoid[nlocal] == 1) ellipsoid[nlocal] = 0;
1221 else error->one(FLERR,"Invalid atom type in Atoms section of data file");
1222
1223 rmass[nlocal] = atof(values[1]);
1224 if (rmass[nlocal] <= 0.0)
1225 error->one(FLERR,"Invalid density in Atoms section of data file");
1226
1227 return 2;
1228 }
1229
1230 /* ----------------------------------------------------------------------
1231 unpack one line from Ellipsoids section of data file
1232 ------------------------------------------------------------------------- */
1233
data_atom_bonus(int m,char ** values)1234 void AtomVecEllipsoid::data_atom_bonus(int m, char **values)
1235 {
1236 if (ellipsoid[m])
1237 error->one(FLERR,"Assigning ellipsoid parameters to non-ellipsoid atom");
1238
1239 if (nlocal_bonus == nmax_bonus) grow_bonus();
1240
1241 double *shape = bonus[nlocal_bonus].shape;
1242 shape[0] = 0.5 * atof(values[0]);
1243 shape[1] = 0.5 * atof(values[1]);
1244 shape[2] = 0.5 * atof(values[2]);
1245 if (shape[0] <= 0.0 || shape[1] <= 0.0 || shape[2] <= 0.0)
1246 error->one(FLERR,"Invalid shape in Ellipsoids section of data file");
1247
1248 double *quat = bonus[nlocal_bonus].quat;
1249 quat[0] = atof(values[3]);
1250 quat[1] = atof(values[4]);
1251 quat[2] = atof(values[5]);
1252 quat[3] = atof(values[6]);
1253 MathExtra::qnormalize(quat);
1254
1255 // reset ellipsoid mass
1256 // previously stored density in rmass
1257
1258 rmass[m] *= 4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2];
1259
1260 bonus[nlocal_bonus].ilocal = m;
1261 ellipsoid[m] = nlocal_bonus++;
1262 }
1263
1264 /* ----------------------------------------------------------------------
1265 unpack one line from Velocities section of data file
1266 ------------------------------------------------------------------------- */
1267
data_vel(int m,char ** values)1268 void AtomVecEllipsoid::data_vel(int m, char **values)
1269 {
1270 v[m][0] = atof(values[0]);
1271 v[m][1] = atof(values[1]);
1272 v[m][2] = atof(values[2]);
1273 angmom[m][0] = atof(values[3]);
1274 angmom[m][1] = atof(values[4]);
1275 angmom[m][2] = atof(values[5]);
1276 }
1277
1278 /* ----------------------------------------------------------------------
1279 unpack hybrid quantities from one line in Velocities section of data file
1280 ------------------------------------------------------------------------- */
1281
data_vel_hybrid(int m,char ** values)1282 int AtomVecEllipsoid::data_vel_hybrid(int m, char **values)
1283 {
1284 angmom[m][0] = atof(values[0]);
1285 angmom[m][1] = atof(values[1]);
1286 angmom[m][2] = atof(values[2]);
1287 return 3;
1288 }
1289
1290 /* ----------------------------------------------------------------------
1291 pack atom info for data file including 3 image flags
1292 ------------------------------------------------------------------------- */
1293
pack_data(double ** buf)1294 void AtomVecEllipsoid::pack_data(double **buf)
1295 {
1296 double *shape;
1297
1298 int nlocal = atom->nlocal;
1299 for (int i = 0; i < nlocal; i++) {
1300 buf[i][0] = ubuf(tag[i]).d;
1301 buf[i][1] = ubuf(type[i]).d;
1302 if (ellipsoid[i] < 0) buf[i][2] = ubuf(0).d;
1303 else buf[i][2] = ubuf(1).d;
1304 if (ellipsoid[i] < 0) buf[i][3] = rmass[i];
1305 else {
1306 shape = bonus[ellipsoid[i]].shape;
1307 buf[i][3] = rmass[i] / (4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2]);
1308 }
1309 buf[i][4] = x[i][0];
1310 buf[i][5] = x[i][1];
1311 buf[i][6] = x[i][2];
1312 buf[i][7] = ubuf((image[i] & IMGMASK) - IMGMAX).d;
1313 buf[i][8] = ubuf((image[i] >> IMGBITS & IMGMASK) - IMGMAX).d;
1314 buf[i][9] = ubuf((image[i] >> IMG2BITS) - IMGMAX).d;
1315 }
1316 }
1317
1318 /* ----------------------------------------------------------------------
1319 pack hybrid atom info for data file
1320 ------------------------------------------------------------------------- */
1321
pack_data_hybrid(int i,double * buf)1322 int AtomVecEllipsoid::pack_data_hybrid(int i, double *buf)
1323 {
1324 if (ellipsoid[i] < 0) buf[0] = ubuf(0).d;
1325 else buf[0] = ubuf(1).d;
1326 if (ellipsoid[i] < 0) buf[1] = rmass[i];
1327 else {
1328 double *shape = bonus[ellipsoid[i]].shape;
1329 buf[1] = rmass[i] / (4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2]);
1330 }
1331 return 2;
1332 }
1333
1334 /* ----------------------------------------------------------------------
1335 write atom info to data file including 3 image flags
1336 ------------------------------------------------------------------------- */
1337
write_data(FILE * fp,int n,double ** buf)1338 void AtomVecEllipsoid::write_data(FILE *fp, int n, double **buf)
1339 {
1340 for (int i = 0; i < n; i++)
1341 fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d\n",
1342 (int) ubuf(buf[i][0]).i,(int) ubuf(buf[i][1]).i,
1343 (int) ubuf(buf[i][2]).i,
1344 buf[i][3],buf[i][4],buf[i][5],buf[i][6],
1345 (int) ubuf(buf[i][7]).i,(int) ubuf(buf[i][8]).i,
1346 (int) ubuf(buf[i][9]).i);
1347 }
1348
1349 /* ----------------------------------------------------------------------
1350 write hybrid atom info to data file
1351 ------------------------------------------------------------------------- */
1352
write_data_hybrid(FILE * fp,double * buf)1353 int AtomVecEllipsoid::write_data_hybrid(FILE *fp, double *buf)
1354 {
1355 fprintf(fp," %d %-1.16e",(int) ubuf(buf[0]).i,buf[1]);
1356 return 2;
1357 }
1358
1359 /* ----------------------------------------------------------------------
1360 pack velocity info for data file
1361 ------------------------------------------------------------------------- */
1362
pack_vel(double ** buf)1363 void AtomVecEllipsoid::pack_vel(double **buf)
1364 {
1365 int nlocal = atom->nlocal;
1366 for (int i = 0; i < nlocal; i++) {
1367 buf[i][0] = ubuf(tag[i]).d;
1368 buf[i][1] = v[i][0];
1369 buf[i][2] = v[i][1];
1370 buf[i][3] = v[i][2];
1371 buf[i][4] = angmom[i][0];
1372 buf[i][5] = angmom[i][1];
1373 buf[i][6] = angmom[i][2];
1374 }
1375 }
1376
1377 /* ----------------------------------------------------------------------
1378 pack hybrid velocity info for data file
1379 ------------------------------------------------------------------------- */
1380
pack_vel_hybrid(int i,double * buf)1381 int AtomVecEllipsoid::pack_vel_hybrid(int i, double *buf)
1382 {
1383 buf[0] = angmom[i][0];
1384 buf[1] = angmom[i][1];
1385 buf[2] = angmom[i][2];
1386 return 3;
1387 }
1388
1389 /* ----------------------------------------------------------------------
1390 write velocity info to data file
1391 ------------------------------------------------------------------------- */
1392
write_vel(FILE * fp,int n,double ** buf)1393 void AtomVecEllipsoid::write_vel(FILE *fp, int n, double **buf)
1394 {
1395 for (int i = 0; i < n; i++)
1396 fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e\n",
1397 (int) ubuf(buf[i][0]).i,buf[i][1],buf[i][2],buf[i][3],
1398 buf[i][4],buf[i][5],buf[i][6]);
1399 }
1400
1401 /* ----------------------------------------------------------------------
1402 write hybrid velocity info to data file
1403 ------------------------------------------------------------------------- */
1404
write_vel_hybrid(FILE * fp,double * buf)1405 int AtomVecEllipsoid::write_vel_hybrid(FILE *fp, double *buf)
1406 {
1407 fprintf(fp," %-1.16e %-1.16e %-1.16e",buf[0],buf[1],buf[2]);
1408 return 3;
1409 }
1410
1411 /* ----------------------------------------------------------------------
1412 return # of bytes of allocated memory
1413 ------------------------------------------------------------------------- */
1414
memory_usage()1415 bigint AtomVecEllipsoid::memory_usage()
1416 {
1417 bigint bytes = 0;
1418
1419 if (atom->memcheck("tag")) bytes += memory->usage(tag,nmax);
1420 if (atom->memcheck("type")) bytes += memory->usage(type,nmax);
1421 if (atom->memcheck("mask")) bytes += memory->usage(mask,nmax);
1422 if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
1423 if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
1424 if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
1425 if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
1426
1427 if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
1428 if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
1429 if (atom->memcheck("torque"))
1430 bytes += memory->usage(torque,nmax*comm->nthreads,3);
1431 if (atom->memcheck("ellipsoid")) bytes += memory->usage(ellipsoid,nmax);
1432
1433 bytes += nmax_bonus*sizeof(Bonus);
1434
1435 return bytes;
1436 }
1437