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 /* ----------------------------------------------------------------------
16 FixPOEMS is a LAMMPS interface to the POEMS coupled multi-body simulator
17 POEMS authors: Rudranarayan Mukherjee (mukher@rpi.edu)
18 Kurt Anderson (anderk5@rpi.edu)
19 ------------------------------------------------------------------------- */
20
21 #include "fix_poems.h"
22
23 #include <cmath>
24 #include <cstring>
25
26 #include "workspace.h"
27 #include "atom.h"
28 #include "domain.h"
29 #include "update.h"
30 #include "respa.h"
31 #include "modify.h"
32 #include "force.h"
33 #include "group.h"
34 #include "comm.h"
35 #include "citeme.h"
36 #include "memory.h"
37 #include "error.h"
38 #include "math_eigen.h"
39
40 using namespace LAMMPS_NS;
41 using namespace FixConst;
42
43 #define MAXBODY 2 // currently 2 since only linear chains allowed
44 #define DELTA 128
45 #define TOLERANCE 1.0e-6
46 #define EPSILON 1.0e-7
47
48 static const char cite_fix_poems[] =
49 "fix poems command:\n\n"
50 "@Article{Mukherjee08,\n"
51 " author = {R. M. Mukherjee, P. S. Crozier, S. J. Plimpton, K. S. Anderson},\n"
52 " title = {Substructured molecular dynamics using multibody dynamics algorithms},\n"
53 " journal = {Intl.~J.~Non-linear Mechanics},\n"
54 " year = 2008,\n"
55 " volume = 43,\n"
56 " pages = {1045--1055}\n"
57 "}\n\n";
58
59 /* ----------------------------------------------------------------------
60 define rigid bodies and joints, initiate POEMS
61 ------------------------------------------------------------------------- */
62
FixPOEMS(LAMMPS * lmp,int narg,char ** arg)63 FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) :
64 Fix(lmp, narg, arg), step_respa(nullptr), natom2body(nullptr),
65 atom2body(nullptr), displace(nullptr), nrigid(nullptr), masstotal(nullptr),
66 xcm(nullptr), vcm(nullptr), fcm(nullptr), inertia(nullptr), ex_space(nullptr),
67 ey_space(nullptr), ez_space(nullptr), angmom(nullptr), omega(nullptr),
68 torque(nullptr), sum(nullptr), all(nullptr), jointbody(nullptr),
69 xjoint(nullptr), freelist(nullptr), poems(nullptr)
70 {
71 if (lmp->citeme) lmp->citeme->add(cite_fix_poems);
72
73 int i,j,ibody;
74
75 time_integrate = 1;
76 rigid_flag = 1;
77 virial_global_flag = virial_peratom_flag = 1;
78 centroidstressflag = CENTROID_NOTAVAIL;
79 thermo_virial = 1;
80 dof_flag = 1;
81
82 MPI_Comm_rank(world,&me);
83
84 // perform initial allocation of atom-based arrays
85 // register with atom class
86
87 natom2body = nullptr;
88 atom2body = nullptr;
89 displace = nullptr;
90 grow_arrays(atom->nmax);
91 atom->add_callback(Atom::GROW);
92
93 // initialize each atom to belong to no rigid bodies
94
95 int nlocal = atom->nlocal;
96 for (int i = 0; i < nlocal; i++) natom2body[i] = 0;
97
98 // create an atom map if one doesn't exist already
99 // readfile() and jointbuild() use global atom IDs
100
101 int mapflag = 0;
102 if (atom->map_style == Atom::MAP_NONE) {
103 mapflag = 1;
104 atom->map_init();
105 atom->map_set();
106 }
107
108 // parse command-line args
109 // set natom2body, atom2body for all atoms and nbody = # of rigid bodies
110 // atoms must also be in fix group to be in a body
111
112 if (narg < 4) error->all(FLERR,"Illegal fix poems command");
113
114 // group = arg has list of groups
115
116 if (strcmp(arg[3],"group") == 0) {
117 nbody = narg-4;
118 if (nbody <= 0) error->all(FLERR,"Illegal fix poems command");
119
120 int *igroups = new int[nbody];
121 for (ibody = 0; ibody < nbody; ibody++) {
122 igroups[ibody] = group->find(arg[ibody+4]);
123 if (igroups[ibody] == -1)
124 error->all(FLERR,"Could not find fix poems group ID");
125 }
126
127 int *mask = atom->mask;
128
129 for (int i = 0; i < nlocal; i++) {
130 if (mask[i] & groupbit)
131 for (ibody = 0; ibody < nbody; ibody++)
132 if (mask[i] & group->bitmask[igroups[ibody]]) {
133 if (natom2body[i] < MAXBODY) atom2body[i][natom2body[i]] = ibody;
134 natom2body[i]++;
135 }
136 }
137
138 delete [] igroups;
139
140 // file = read bodies from file
141 // file read doesn't pay attention to fix group,
142 // so after read, reset natom2body = 0 if atom is not in fix group
143
144 } else if (strcmp(arg[3],"file") == 0) {
145
146 readfile(arg[4]);
147
148 int *mask = atom->mask;
149 for (int i = 0; i < nlocal; i++)
150 if (!(mask[i] & groupbit)) natom2body[i] = 0;
151
152 // each molecule in fix group is a rigid body
153 // maxmol = largest molecule ID
154 // ncount = # of atoms in each molecule (have to sum across procs)
155 // nbody = # of non-zero ncount values
156 // use nall as incremented ptr to set atom2body[] values for each atom
157
158 } else if (strcmp(arg[3],"molecule") == 0) {
159 if (narg != 4) error->all(FLERR,"Illegal fix poems command");
160 if (atom->molecular == Atom::ATOMIC)
161 error->all(FLERR,
162 "Must use a molecular atom style with fix poems molecule");
163
164 int *mask = atom->mask;
165 tagint *molecule = atom->molecule;
166 int nlocal = atom->nlocal;
167
168 tagint maxmol_tag = -1;
169 for (i = 0; i < nlocal; i++)
170 if (mask[i] & groupbit) maxmol_tag = MAX(maxmol_tag,molecule[i]);
171
172 tagint itmp;
173 MPI_Allreduce(&maxmol_tag,&itmp,1,MPI_LMP_TAGINT,MPI_MAX,world);
174 if (itmp+1 > MAXSMALLINT)
175 error->all(FLERR,"Too many molecules for fix poems");
176 int maxmol = (int) itmp;
177
178 int *ncount;
179 memory->create(ncount,maxmol+1,"rigid:ncount");
180 for (i = 0; i <= maxmol; i++) ncount[i] = 0;
181
182 for (i = 0; i < nlocal; i++)
183 if (mask[i] & groupbit) ncount[molecule[i]]++;
184
185 int *nall;
186 memory->create(nall,maxmol+1,"rigid:ncount");
187 MPI_Allreduce(ncount,nall,maxmol+1,MPI_INT,MPI_SUM,world);
188
189 nbody = 0;
190 for (i = 0; i <= maxmol; i++)
191 if (nall[i]) nall[i] = nbody++;
192 else nall[i] = -1;
193
194 for (i = 0; i < nlocal; i++) {
195 natom2body[i] = 0;
196 if (mask[i] & groupbit) {
197 natom2body[i] = 1;
198 atom2body[i][0] = nall[molecule[i]];
199 }
200 }
201
202 memory->destroy(ncount);
203 memory->destroy(nall);
204
205 } else error->all(FLERR,"Illegal fix poems command");
206
207 // error if no bodies
208 // error if any atom in too many bodies
209
210 if (nbody == 0) error->all(FLERR,"No rigid bodies defined");
211
212 int flag = 0;
213 for (int i = 0; i < nlocal; i++)
214 if (natom2body[i] > MAXBODY) flag = 1;
215 int flagall;
216 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
217 if (flagall)
218 error->all(FLERR,"Atom in too many rigid bodies - boost MAXBODY");
219
220 // create all nbody-length arrays
221
222 nrigid = new int[nbody];
223 masstotal = new double[nbody];
224 memory->create(xcm,nbody,3,"poems:xcm");
225 memory->create(vcm,nbody,3,"poems:vcm");
226 memory->create(fcm,nbody,3,"poems:fcm");
227 memory->create(inertia,nbody,3,"poems:inertia");
228 memory->create(ex_space,nbody,3,"poems:ex_space");
229 memory->create(ey_space,nbody,3,"poems:ey_space");
230 memory->create(ez_space,nbody,3,"poems:ez_space");
231 memory->create(angmom,nbody,3,"poems:angmom");
232 memory->create(omega,nbody,3,"poems:omega");
233 memory->create(torque,nbody,3,"poems:torque");
234
235 memory->create(sum,nbody,6,"poems:sum");
236 memory->create(all,nbody,6,"poems:all");
237
238 // nrigid[n] = # of atoms in Nth rigid body
239 // double count joint atoms as being in multiple bodies
240 // error if one or zero atoms
241
242 int *ncount = new int[nbody];
243 for (ibody = 0; ibody < nbody; ibody++) ncount[ibody] = 0;
244
245 for (i = 0; i < nlocal; i++)
246 for (j = 0; j < natom2body[i]; j++)
247 ncount[atom2body[i][j]]++;
248
249 MPI_Allreduce(ncount,nrigid,nbody,MPI_INT,MPI_SUM,world);
250 delete [] ncount;
251
252 for (ibody = 0; ibody < nbody; ibody++)
253 if (nrigid[ibody] <= 1) error->all(FLERR,"One or zero atoms in rigid body");
254
255 // build list of joint connections and check for cycles and trees
256
257 jointbuild();
258
259 // delete temporary atom map
260
261 if (mapflag) {
262 atom->map_delete();
263 atom->map_style = Atom::MAP_NONE;
264 }
265
266 // create POEMS instance
267
268 poems = new Workspace;
269
270 // compute per body forces and torques inside final_integrate() by default
271
272 earlyflag = 0;
273
274 // print statistics
275
276 int nsum = 0;
277 for (ibody = 0; ibody < nbody; ibody++) nsum += nrigid[ibody];
278 nsum -= njoint;
279
280 if (me == 0)
281 utils::logmesg(lmp,"{} clusters, {} bodies, {} joints, {} atoms\n",
282 ncluster,nbody,njoint,nsum);
283 }
284
285 /* ----------------------------------------------------------------------
286 free all memory for rigid bodies, joints, and POEMS
287 ------------------------------------------------------------------------- */
288
~FixPOEMS()289 FixPOEMS::~FixPOEMS()
290 {
291 // if atom class still exists:
292 // unregister this fix so atom class doesn't invoke it any more
293
294 if (atom) atom->delete_callback(id,Atom::GROW);
295
296 // delete locally stored arrays
297
298 memory->destroy(natom2body);
299 memory->destroy(atom2body);
300 memory->destroy(displace);
301
302 // delete nbody-length arrays
303
304 delete [] nrigid;
305 delete [] masstotal;
306 memory->destroy(xcm);
307 memory->destroy(vcm);
308 memory->destroy(fcm);
309 memory->destroy(inertia);
310 memory->destroy(ex_space);
311 memory->destroy(ey_space);
312 memory->destroy(ez_space);
313 memory->destroy(angmom);
314 memory->destroy(omega);
315 memory->destroy(torque);
316
317 memory->destroy(sum);
318 memory->destroy(all);
319
320 // delete joint arrays
321
322 memory->destroy(jointbody);
323 memory->destroy(xjoint);
324 delete [] freelist;
325
326 // delete POEMS object
327
328 delete poems;
329 }
330
331 /* ---------------------------------------------------------------------- */
332
setmask()333 int FixPOEMS::setmask()
334 {
335 int mask = 0;
336 mask |= INITIAL_INTEGRATE;
337 mask |= FINAL_INTEGRATE;
338 mask |= PRE_NEIGHBOR;
339 mask |= POST_FORCE;
340 mask |= INITIAL_INTEGRATE_RESPA;
341 mask |= FINAL_INTEGRATE_RESPA;
342 mask |= POST_FORCE_RESPA;
343 return mask;
344 }
345
346 /* ---------------------------------------------------------------------- */
347
init()348 void FixPOEMS::init()
349 {
350 int i,ibody;
351
352 // warn if more than one POEMS fix
353 // if earlyflag, warn if any post-force fixes come after POEMS fix
354
355 int count = 0;
356 for (int i = 0; i < modify->nfix; i++)
357 if (strcmp(modify->fix[i]->style,"poems") == 0) count++;
358 if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one fix poems");
359
360 if (earlyflag) {
361 int pflag = 0;
362 for (i = 0; i < modify->nfix; i++) {
363 if (strcmp(modify->fix[i]->style,"poems") == 0) pflag = 1;
364 if (pflag && (modify->fmask[i] & POST_FORCE) &&
365 !modify->fix[i]->rigid_flag) {
366 if (comm->me == 0)
367 error->warning(FLERR,std::string("Fix ") + modify->fix[i]->id
368 + std::string(" alters forces after fix poems"));
369 }
370 }
371 }
372
373 // error if npt,nph fix comes before rigid fix
374
375 for (i = 0; i < modify->nfix; i++) {
376 if (strcmp(modify->fix[i]->style,"npt") == 0) break;
377 if (strcmp(modify->fix[i]->style,"nph") == 0) break;
378 }
379 if (i < modify->nfix) {
380 for (int j = i; j < modify->nfix; j++)
381 if (strcmp(modify->fix[j]->style,"poems") == 0)
382 error->all(FLERR,"POEMS fix must come before NPT/NPH fix");
383 }
384
385 // timestep info
386
387 dtv = update->dt;
388 dtf = 0.5 * update->dt * force->ftm2v;
389 dthalf = 0.5 * update->dt;
390
391 // rRESPA info
392
393 if (utils::strmatch(update->integrate_style,"^respa")) {
394 step_respa = ((Respa *) update->integrate)->step;
395 nlevels_respa = ((Respa *) update->integrate)->nlevels;
396 }
397
398 // compute masstotal & center-of-mass xcm of each rigid body
399 // only count joint atoms in 1st body
400
401 int *type = atom->type;
402 imageint *image = atom->image;
403 double *mass = atom->mass;
404 double **x = atom->x;
405 double **v = atom->v;
406 int nlocal = atom->nlocal;
407
408 double xprd = domain->xprd;
409 double yprd = domain->yprd;
410 double zprd = domain->zprd;
411
412 int xbox,ybox,zbox;
413 double massone;
414
415 for (ibody = 0; ibody < nbody; ibody++)
416 for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
417
418 for (i = 0; i < nlocal; i++) {
419 if (natom2body[i]) {
420 ibody = atom2body[i][0];
421 xbox = (image[i] & IMGMASK) - IMGMAX;
422 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
423 zbox = (image[i] >> IMG2BITS) - IMGMAX;
424 massone = mass[type[i]];
425 sum[ibody][0] += (x[i][0] + xbox*xprd) * massone;
426 sum[ibody][1] += (x[i][1] + ybox*yprd) * massone;
427 sum[ibody][2] += (x[i][2] + zbox*zprd) * massone;
428 sum[ibody][3] += massone;
429 sum[ibody][4] += massone *
430 (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
431 }
432 }
433
434 MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
435
436 total_ke = 0.0;
437 for (ibody = 0; ibody < nbody; ibody++) {
438 masstotal[ibody] = all[ibody][3];
439 xcm[ibody][0] = all[ibody][0]/masstotal[ibody];
440 xcm[ibody][1] = all[ibody][1]/masstotal[ibody];
441 xcm[ibody][2] = all[ibody][2]/masstotal[ibody];
442 total_ke += 0.5 * all[ibody][4];
443 }
444
445 // compute 6 moments of inertia of each body
446 // only count joint atoms in 1st body
447 // dx,dy,dz = coords relative to center-of-mass
448
449 double dx,dy,dz;
450
451 for (ibody = 0; ibody < nbody; ibody++)
452 for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
453
454 for (i = 0; i < nlocal; i++) {
455 if (natom2body[i]) {
456 ibody = atom2body[i][0];
457
458 xbox = (image[i] & IMGMASK) - IMGMAX;
459 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
460 zbox = (image[i] >> IMG2BITS) - IMGMAX;
461 dx = x[i][0] + xbox*xprd - xcm[ibody][0];
462 dy = x[i][1] + ybox*yprd - xcm[ibody][1];
463 dz = x[i][2] + zbox*zprd - xcm[ibody][2];
464 massone = mass[type[i]];
465
466 sum[ibody][0] += massone * (dy*dy + dz*dz);
467 sum[ibody][1] += massone * (dx*dx + dz*dz);
468 sum[ibody][2] += massone * (dx*dx + dy*dy);
469 sum[ibody][3] -= massone * dx*dy;
470 sum[ibody][4] -= massone * dy*dz;
471 sum[ibody][5] -= massone * dx*dz;
472 }
473 }
474
475 MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
476
477 // inertia = 3 eigenvalues = principal moments of inertia
478 // ex_space,ey_space,ez_space = 3 eigenvectors = principal axes of rigid body
479
480 double **tensor,**evectors;
481 memory->create(tensor,3,3,"fix_rigid:tensor");
482 memory->create(evectors,3,3,"fix_rigid:evectors");
483
484 int ierror;
485 double ez0,ez1,ez2;
486
487 for (ibody = 0; ibody < nbody; ibody++) {
488 tensor[0][0] = all[ibody][0];
489 tensor[1][1] = all[ibody][1];
490 tensor[2][2] = all[ibody][2];
491 tensor[0][1] = tensor[1][0] = all[ibody][3];
492 tensor[1][2] = tensor[2][1] = all[ibody][4];
493 tensor[0][2] = tensor[2][0] = all[ibody][5];
494
495 ierror = MathEigen::jacobi3(tensor,inertia[ibody],evectors);
496 if (ierror) error->all(FLERR,"Insufficient Jacobi rotations for POEMS body");
497
498 ex_space[ibody][0] = evectors[0][0];
499 ex_space[ibody][1] = evectors[1][0];
500 ex_space[ibody][2] = evectors[2][0];
501
502 ey_space[ibody][0] = evectors[0][1];
503 ey_space[ibody][1] = evectors[1][1];
504 ey_space[ibody][2] = evectors[2][1];
505
506 ez_space[ibody][0] = evectors[0][2];
507 ez_space[ibody][1] = evectors[1][2];
508 ez_space[ibody][2] = evectors[2][2];
509
510 // if any principal moment < scaled EPSILON, error
511 // this is b/c POEMS cannot yet handle degenerate bodies
512
513 double max;
514 max = MAX(inertia[ibody][0],inertia[ibody][1]);
515 max = MAX(max,inertia[ibody][2]);
516
517 if (inertia[ibody][0] < EPSILON*max ||
518 inertia[ibody][1] < EPSILON*max ||
519 inertia[ibody][2] < EPSILON*max)
520 error->all(FLERR,"Rigid body has degenerate moment of inertia");
521
522 // enforce 3 evectors as a right-handed coordinate system
523 // flip 3rd evector if needed
524
525 ez0 = ex_space[ibody][1]*ey_space[ibody][2] -
526 ex_space[ibody][2]*ey_space[ibody][1];
527 ez1 = ex_space[ibody][2]*ey_space[ibody][0] -
528 ex_space[ibody][0]*ey_space[ibody][2];
529 ez2 = ex_space[ibody][0]*ey_space[ibody][1] -
530 ex_space[ibody][1]*ey_space[ibody][0];
531
532 if (ez0*ez_space[ibody][0] + ez1*ez_space[ibody][1] +
533 ez2*ez_space[ibody][2] < 0.0) {
534 ez_space[ibody][0] = -ez_space[ibody][0];
535 ez_space[ibody][1] = -ez_space[ibody][1];
536 ez_space[ibody][2] = -ez_space[ibody][2];
537 }
538 }
539
540 // free temporary memory
541
542 memory->destroy(tensor);
543 memory->destroy(evectors);
544
545 // displace = initial atom coords in basis of principal axes
546 // only set joint atoms relative to 1st body
547 // set displace = 0.0 for atoms not in any rigid body
548
549 for (i = 0; i < nlocal; i++) {
550 if (natom2body[i]) {
551 ibody = atom2body[i][0];
552
553 xbox = (image[i] & IMGMASK) - IMGMAX;
554 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
555 zbox = (image[i] >> IMG2BITS) - IMGMAX;
556 dx = x[i][0] + xbox*xprd - xcm[ibody][0];
557 dy = x[i][1] + ybox*yprd - xcm[ibody][1];
558 dz = x[i][2] + zbox*zprd - xcm[ibody][2];
559
560 displace[i][0] = dx*ex_space[ibody][0] + dy*ex_space[ibody][1] +
561 dz*ex_space[ibody][2];
562 displace[i][1] = dx*ey_space[ibody][0] + dy*ey_space[ibody][1] +
563 dz*ey_space[ibody][2];
564 displace[i][2] = dx*ez_space[ibody][0] + dy*ez_space[ibody][1] +
565 dz*ez_space[ibody][2];
566 } else displace[i][0] = displace[i][1] = displace[i][2] = 0.0;
567 }
568
569 // test for valid principal moments & axes
570 // recompute moments of inertia around new axes
571 // only count joint atoms in 1st body
572 // 3 diagonal moments should equal principal moments
573 // 3 off-diagonal moments should be 0.0
574 // (ddx,ddy,ddz) is projection of atom within rigid body onto principal axes
575 // 6 moments use (ddx,ddy,ddz) displacements from principal axes
576
577 for (ibody = 0; ibody < nbody; ibody++)
578 for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
579
580 double ddx,ddy,ddz;
581
582 for (i = 0; i < nlocal; i++) {
583 if (natom2body[i]) {
584 ibody = atom2body[i][0];
585
586 xbox = (image[i] & IMGMASK) - IMGMAX;
587 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
588 zbox = (image[i] >> IMG2BITS) - IMGMAX;
589 dx = x[i][0] + xbox*xprd - xcm[ibody][0];
590 dy = x[i][1] + ybox*yprd - xcm[ibody][1];
591 dz = x[i][2] + zbox*zprd - xcm[ibody][2];
592 massone = mass[type[i]];
593
594 ddx = dx*ex_space[ibody][0] + dy*ex_space[ibody][1] +
595 dz*ex_space[ibody][2];
596 ddy = dx*ey_space[ibody][0] + dy*ey_space[ibody][1] +
597 dz*ey_space[ibody][2];
598 ddz = dx*ez_space[ibody][0] + dy*ez_space[ibody][1] +
599 dz*ez_space[ibody][2];
600
601 sum[ibody][0] += massone * (ddy*ddy + ddz*ddz);
602 sum[ibody][1] += massone * (ddx*ddx + ddz*ddz);
603 sum[ibody][2] += massone * (ddx*ddx + ddy*ddy);
604 sum[ibody][3] -= massone * ddx*ddy;
605 sum[ibody][4] -= massone * ddy*ddz;
606 sum[ibody][5] -= massone * ddx*ddz;
607 }
608 }
609
610 MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
611
612 for (ibody = 0; ibody < nbody; ibody++) {
613 if (fabs(all[ibody][0]-inertia[ibody][0]) > TOLERANCE ||
614 fabs(all[ibody][1]-inertia[ibody][1]) > TOLERANCE ||
615 fabs(all[ibody][2]-inertia[ibody][2]) > TOLERANCE)
616 error->all(FLERR,"Bad principal moments");
617 if (fabs(all[ibody][3]) > TOLERANCE ||
618 fabs(all[ibody][4]) > TOLERANCE ||
619 fabs(all[ibody][5]) > TOLERANCE)
620 error->all(FLERR,"Bad principal moments");
621 }
622 }
623
624 /* ----------------------------------------------------------------------
625 compute initial rigid body info
626 make setup call to POEMS
627 ------------------------------------------------------------------------- */
628
setup(int vflag)629 void FixPOEMS::setup(int vflag)
630 {
631 int i,n,ibody;
632
633 // vcm = velocity of center-of-mass of each rigid body
634 // angmom = angular momentum of each rigid body
635 // only count joint atoms in 1st body
636
637 int *type = atom->type;
638 imageint *image = atom->image;
639 double *mass = atom->mass;
640 double **x = atom->x;
641 double **v = atom->v;
642 int nlocal = atom->nlocal;
643
644 double xprd = domain->xprd;
645 double yprd = domain->yprd;
646 double zprd = domain->zprd;
647
648 int xbox,ybox,zbox;
649 double massone,dx,dy,dz;
650
651 for (ibody = 0; ibody < nbody; ibody++)
652 for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
653
654 for (i = 0; i < nlocal; i++) {
655 if (natom2body[i]) {
656 ibody = atom2body[i][0];
657 massone = mass[type[i]];
658
659 xbox = (image[i] & IMGMASK) - IMGMAX;
660 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
661 zbox = (image[i] >> IMG2BITS) - IMGMAX;
662 dx = x[i][0] + xbox*xprd - xcm[ibody][0];
663 dy = x[i][1] + ybox*yprd - xcm[ibody][1];
664 dz = x[i][2] + zbox*zprd - xcm[ibody][2];
665
666 sum[ibody][0] += v[i][0] * massone;
667 sum[ibody][1] += v[i][1] * massone;
668 sum[ibody][2] += v[i][2] * massone;
669 sum[ibody][3] += dy * massone*v[i][2] - dz * massone*v[i][1];
670 sum[ibody][4] += dz * massone*v[i][0] - dx * massone*v[i][2];
671 sum[ibody][5] += dx * massone*v[i][1] - dy * massone*v[i][0];
672 }
673 }
674
675 MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
676
677 for (ibody = 0; ibody < nbody; ibody++) {
678 vcm[ibody][0] = all[ibody][0]/masstotal[ibody];
679 vcm[ibody][1] = all[ibody][1]/masstotal[ibody];
680 vcm[ibody][2] = all[ibody][2]/masstotal[ibody];
681 angmom[ibody][0] = all[ibody][3];
682 angmom[ibody][1] = all[ibody][4];
683 angmom[ibody][2] = all[ibody][5];
684 }
685
686 // virial setup before call to set_v
687
688 v_init(vflag);
689
690 // set velocities from angmom & omega
691
692 for (ibody = 0; ibody < nbody; ibody++)
693 omega_from_mq(angmom[ibody],ex_space[ibody],ey_space[ibody],
694 ez_space[ibody],inertia[ibody],omega[ibody]);
695 set_v();
696
697 // guestimate virial as 2x the set_v contribution
698
699 if (evflag) {
700 if (vflag_global)
701 for (n = 0; n < 6; n++) virial[n] *= 2.0;
702 if (vflag_atom) {
703 for (i = 0; i < nlocal; i++)
704 for (n = 0; n < 6; n++)
705 vatom[i][n] *= 2.0;
706 }
707 }
708
709 // use post_force() to compute initial fcm & torque
710
711 compute_forces_and_torques();
712 //post_force(vflag);
713
714 // setup for POEMS
715
716 poems->MakeSystem(nbody,masstotal,inertia,xcm,vcm,omega,
717 ex_space,ey_space,ez_space,
718 njoint,jointbody,xjoint,nfree,freelist,
719 dthalf,dtv,force->ftm2v,total_ke);
720 }
721
722 /* ----------------------------------------------------------------------
723 update vcm,omega by 1/2 step and xcm,orientation by full step
724 set x,v of body atoms accordingly
725 ---------------------------------------------------------------------- */
726
initial_integrate(int vflag)727 void FixPOEMS::initial_integrate(int vflag)
728 {
729 // perform POEMS integration
730
731 poems->LobattoOne(xcm,vcm,omega,torque,fcm,ex_space,ey_space,ez_space);
732
733 // virial setup before call to set_xv
734
735 v_init(vflag);
736
737 // set coords and velocities of atoms in rigid bodies
738
739 set_xv();
740 }
741
742 /* ---------------------------------------------------------------------- */
743
post_force(int)744 void FixPOEMS::post_force(int /* vflag */)
745 {
746 if (earlyflag) compute_forces_and_torques();
747 }
748
749 /* ----------------------------------------------------------------------
750 compute fcm,torque on each rigid body
751 only count joint atoms in 1st body
752 ------------------------------------------------------------------------- */
753
compute_forces_and_torques()754 void FixPOEMS::compute_forces_and_torques()
755 {
756 int i,ibody;
757 int xbox,ybox,zbox;
758 double dx,dy,dz;
759
760 imageint *image = atom->image;
761 double **x = atom->x;
762 double **f = atom->f;
763 int nlocal = atom->nlocal;
764
765 double xprd = domain->xprd;
766 double yprd = domain->yprd;
767 double zprd = domain->zprd;
768
769 for (ibody = 0; ibody < nbody; ibody++)
770 for (i = 0; i < 6; i++) sum[ibody][i] = 0.0;
771
772 for (i = 0; i < nlocal; i++) {
773 if (natom2body[i]) {
774 ibody = atom2body[i][0];
775
776 sum[ibody][0] += f[i][0];
777 sum[ibody][1] += f[i][1];
778 sum[ibody][2] += f[i][2];
779
780 xbox = (image[i] & IMGMASK) - IMGMAX;
781 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
782 zbox = (image[i] >> IMG2BITS) - IMGMAX;
783 dx = x[i][0] + xbox*xprd - xcm[ibody][0];
784 dy = x[i][1] + ybox*yprd - xcm[ibody][1];
785 dz = x[i][2] + zbox*zprd - xcm[ibody][2];
786
787 sum[ibody][3] += dy*f[i][2] - dz*f[i][1];
788 sum[ibody][4] += dz*f[i][0] - dx*f[i][2];
789 sum[ibody][5] += dx*f[i][1] - dy*f[i][0];
790 }
791 }
792
793 MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world);
794
795 for (ibody = 0; ibody < nbody; ibody++) {
796 fcm[ibody][0] = all[ibody][0];
797 fcm[ibody][1] = all[ibody][1];
798 fcm[ibody][2] = all[ibody][2];
799 torque[ibody][0] = all[ibody][3];
800 torque[ibody][1] = all[ibody][4];
801 torque[ibody][2] = all[ibody][5];
802 }
803 }
804
805 /* ----------------------------------------------------------------------
806 update vcm,omega by last 1/2 step
807 set v of body atoms accordingly
808 ------------------------------------------------------------------------- */
809
final_integrate()810 void FixPOEMS::final_integrate()
811 {
812 if (!earlyflag) compute_forces_and_torques();
813
814 /*
815 for (int ibody = 0; ibody < nbody; ibody++) {
816 if (ibody == 0) {
817 printf("FI %d %g %g %g\n",ibody,fcm[ibody][0],fcm[ibody][1],fcm[ibody][2]);
818 printf("TQ %d %g %g %g\n",ibody,
819 torque[ibody][0],torque[ibody][1],torque[ibody][2]);
820 }
821 }
822 */
823
824 // perform POEMS integration
825
826 poems->LobattoTwo(vcm,omega,torque,fcm);
827
828 // set velocities of atoms in rigid bodies
829 // virial is already setup from initial_integrate
830
831 set_v();
832 }
833
834 /* ---------------------------------------------------------------------- */
835
initial_integrate_respa(int vflag,int ilevel,int)836 void FixPOEMS::initial_integrate_respa(int vflag, int ilevel, int /* iloop */)
837 {
838 dtv = step_respa[ilevel];
839 dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
840 dthalf = 0.5 * step_respa[ilevel];
841
842 if (ilevel == 0) initial_integrate(vflag);
843 else final_integrate();
844 }
845
846 /* ---------------------------------------------------------------------- */
847
post_force_respa(int vflag,int ilevel,int)848 void FixPOEMS::post_force_respa(int vflag, int ilevel, int /* iloop */)
849 {
850 if (ilevel == nlevels_respa-1) post_force(vflag);
851 }
852
853 /* ---------------------------------------------------------------------- */
854
final_integrate_respa(int ilevel,int)855 void FixPOEMS::final_integrate_respa(int ilevel, int /* iloop */)
856 {
857 dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
858 final_integrate();
859 }
860
861 /* ----------------------------------------------------------------------
862 remap xcm of each rigid body back into periodic simulation box
863 done during pre_neighbor so will be after call to pbc()
864 and after fix_deform::pre_exchange() may have flipped box
865 if don't do this, then atoms of a body which drifts far away
866 from a triclinic box will be remapped back into box
867 with huge displacements when the box tilt changes via set_x()
868 NOTE: cannot do this by changing xcm of each body in cluster
869 or even 1st body in cluster
870 b/c POEMS library does not see xcm but only sets xcm
871 so remap needs to be coordinated with POEMS library
872 thus this routine does nothing for now
873 ------------------------------------------------------------------------- */
874
pre_neighbor()875 void FixPOEMS::pre_neighbor() {}
876
877 /* ----------------------------------------------------------------------
878 count # of degrees-of-freedom removed by fix_poems for atoms in igroup
879 ------------------------------------------------------------------------- */
880
dof(int igroup)881 int FixPOEMS::dof(int igroup)
882 {
883 int groupbit = group->bitmask[igroup];
884
885 // ncount = # of atoms in each rigid body that are also in group
886 // only count joint atoms as part of first body
887
888 int *mask = atom->mask;
889 int nlocal = atom->nlocal;
890
891 int *ncount = new int[nbody];
892 for (int ibody = 0; ibody < nbody; ibody++) ncount[ibody] = 0;
893
894 for (int i = 0; i < nlocal; i++)
895 if (mask[i] & groupbit)
896 if (natom2body[i]) ncount[atom2body[i][0]]++;
897
898 int *nall = new int[nbody];
899 MPI_Allreduce(ncount,nall,nbody,MPI_INT,MPI_SUM,world);
900
901 // remove 3N - 6 dof for each rigid body if at least 2 atoms are in igroup
902
903 int n = 0;
904 for (int ibody = 0; ibody < nbody; ibody++)
905 if (nall[ibody] > 2) n += 3*nall[ibody] - 6;
906
907 // subtract 3 additional dof for each joint if atom is also in igroup
908
909 int m = 0;
910 for (int i = 0; i < nlocal; i++)
911 if (natom2body[i] > 1 && (mask[i] & groupbit)) m += 3*(natom2body[i]-1);
912 int mall;
913 MPI_Allreduce(&m,&mall,1,MPI_INT,MPI_SUM,world);
914 n += mall;
915
916 // delete local memory
917
918 delete [] ncount;
919 delete [] nall;
920
921 return n;
922 }
923
924 /* ----------------------------------------------------------------------
925 adjust xcm of each cluster due to box deformation
926 called by various fixes that change box size/shape
927 flag = 0/1 means map from box to lamda coords or vice versa
928 NOTE: cannot do this by changing xcm of each body in cluster
929 or even 1st body in cluster
930 b/c POEMS library does not see xcm but only sets xcm
931 so deform needs to be coordinated with POEMS library
932 thus this routine does nothing for now
933 ------------------------------------------------------------------------- */
934
deform(int)935 void FixPOEMS::deform(int /* flag */) {}
936
937 /* ---------------------------------------------------------------------- */
938
readfile(char * file)939 void FixPOEMS::readfile(char *file)
940 {
941 FILE *fp;
942
943 if (me == 0) {
944 fp = fopen(file,"r");
945 if (fp == nullptr)
946 error->one(FLERR,"Cannot open fix poems file {}: {}",
947 file, utils::getsyserror());
948 }
949
950 nbody = 0;
951 char *line = nullptr;
952 int maxline = 0;
953 char *ptr;
954 int nlocal = atom->nlocal;
955 int i,id,nlen;
956
957 while (1) {
958 if (me == 0) nlen = readline(fp,&line,&maxline);
959 MPI_Bcast(&nlen,1,MPI_INT,0,world);
960 if (nlen == 0) break;
961 MPI_Bcast(line,nlen,MPI_CHAR,0,world);
962
963 ptr = strtok(line," ,\t\n\0");
964 if (ptr == nullptr || ptr[0] == '#') continue;
965 ptr = strtok(nullptr," ,\t\n\0");
966
967 while ((ptr = strtok(nullptr," ,\t\n\0"))) {
968 id = atoi(ptr);
969 i = atom->map(id);
970 if (i < 0 || i >= nlocal) continue;
971 if (natom2body[i] < MAXBODY) atom2body[i][natom2body[i]] = nbody;
972 natom2body[i]++;
973 }
974 nbody++;
975 }
976
977 memory->destroy(line);
978 if (me == 0) fclose(fp);
979 }
980
981 /* ---------------------------------------------------------------------- */
982
readline(FILE * fp,char ** pline,int * pmaxline)983 int FixPOEMS::readline(FILE *fp, char **pline, int *pmaxline)
984 {
985 int n = 0;
986 char *line = *pline;
987 int maxline = *pmaxline;
988
989 while (1) {
990 if (n+1 >= maxline) {
991 maxline += DELTA;
992 memory->grow(line,maxline,"fix_poems:line");
993 }
994 if (fgets(&line[n],maxline-n,fp) == nullptr) {
995 n = 0;
996 break;
997 }
998 n = strlen(line);
999 if (n < maxline-1 || line[n-1] == '\n') break;
1000 }
1001
1002 *pmaxline = maxline;
1003 *pline = line;
1004 return n;
1005 }
1006
1007 /* ----------------------------------------------------------------------
1008 build list of joints and error check for cycles and trees
1009 ------------------------------------------------------------------------- */
1010
jointbuild()1011 void FixPOEMS::jointbuild()
1012 {
1013 int i,j;
1014
1015 // convert atom2body into list of joint atoms on this proc
1016 // mjoint = # of joint atoms in this proc
1017 // an atom in N rigid bodies, infers N-1 joints between 1st body and others
1018 // mylist = [0],[1] = 2 body indices, [2] = global ID of joint atom
1019
1020 tagint *tag = atom->tag;
1021 int nlocal = atom->nlocal;
1022
1023 int mjoint = 0;
1024 for (i = 0; i < nlocal; i++) {
1025 if (natom2body[i] <= 1) continue;
1026 mjoint += natom2body[i]-1;
1027 }
1028
1029 tagint **mylist = nullptr;
1030 if (mjoint) memory->create(mylist,mjoint,3,"poems:mylist");
1031
1032 mjoint = 0;
1033 for (i = 0; i < nlocal; i++) {
1034 if (natom2body[i] <= 1) continue;
1035 for (j = 1; j < natom2body[i]; j++) {
1036 mylist[mjoint][0] = atom2body[i][0];
1037 mylist[mjoint][1] = atom2body[i][j];
1038 mylist[mjoint][2] = tag[i];
1039 mjoint++;
1040 }
1041 }
1042
1043 // jlist = mylist concatenated across all procs via MPI_Allgatherv
1044
1045 MPI_Allreduce(&mjoint,&njoint,1,MPI_INT,MPI_SUM,world);
1046 tagint **jlist = nullptr;
1047 if (njoint) memory->create(jlist,njoint,3,"poems:jlist");
1048
1049 int nprocs;
1050 MPI_Comm_size(world,&nprocs);
1051
1052 int *recvcounts = new int[nprocs];
1053 int tmp = 3*mjoint;
1054 MPI_Allgather(&tmp,1,MPI_INT,recvcounts,1,MPI_INT,world);
1055
1056 int *displs = new int[nprocs];
1057 displs[0] = 0;
1058 for (i = 1; i < nprocs; i++) displs[i] = displs[i-1] + recvcounts[i-1];
1059
1060 // allgather the local joint lists
1061 // 2 versions in case mjoint is 0 on this proc
1062
1063 if (njoint) {
1064 if (mjoint)
1065 MPI_Allgatherv(mylist[0],3*mjoint,MPI_LMP_TAGINT,jlist[0],
1066 recvcounts,displs,MPI_LMP_TAGINT,world);
1067 else
1068 MPI_Allgatherv(nullptr,3*mjoint,MPI_LMP_TAGINT,jlist[0],
1069 recvcounts,displs,MPI_LMP_TAGINT,world);
1070 }
1071
1072 delete [] recvcounts;
1073 delete [] displs;
1074
1075 // warning if no joints
1076
1077 if (njoint == 0 && me == 0)
1078 error->warning(FLERR,
1079 "No joints between rigid bodies, use fix rigid instead");
1080
1081 // sort joint list in ascending order by body indices
1082 // check for loops in joint connections between rigid bodies
1083 // check for trees = same body in more than 2 joints
1084
1085 sortlist(njoint,jlist);
1086
1087 if (loopcheck(nbody,njoint,jlist))
1088 error->all(FLERR,"Cyclic loop in joint connections");
1089
1090 int *bodyflag = new int[nbody];
1091 for (i = 0; i < nbody; i++) bodyflag[i] = 0;
1092 for (i = 0; i < njoint; i++) {
1093 bodyflag[jlist[i][0]]++;
1094 bodyflag[jlist[i][1]]++;
1095 }
1096 for (i = 0; i < nbody; i++)
1097 if (bodyflag[i] > 2)
1098 error->all(FLERR,"Tree structure in joint connections");
1099 delete [] bodyflag;
1100
1101 // allocate and setup joint arrays
1102 // jointbody stores body indices from 1 to Nbody to pass to POEMS
1103 // each proc sets myjoint if it owns joint atom
1104 // MPI_Allreduce gives all procs the xjoint coords
1105
1106 jointbody = nullptr;
1107 xjoint = nullptr;
1108 double **myjoint = nullptr;
1109 if (njoint) {
1110 memory->create(jointbody,njoint,2,"poems:jointbody");
1111 memory->create(xjoint,njoint,3,"poems:xjoint");
1112 memory->create(myjoint,njoint,3,"poems:myjoint");
1113 }
1114
1115 double **x = atom->x;
1116
1117 for (i = 0; i < njoint; i++) {
1118 jointbody[i][0] = jlist[i][0] + 1;
1119 jointbody[i][1] = jlist[i][1] + 1;
1120 j = atom->map(jlist[i][2]);
1121 if (j >= 0 && j < nlocal) {
1122 myjoint[i][0] = x[j][0];
1123 myjoint[i][1] = x[j][1];
1124 myjoint[i][2] = x[j][2];
1125 } else myjoint[i][0] = myjoint[i][1] = myjoint[i][2] = 0.0;
1126 }
1127
1128 if (njoint)
1129 MPI_Allreduce(myjoint[0],xjoint[0],3*njoint,MPI_DOUBLE,MPI_SUM,world);
1130
1131 // compute freelist of nfree single unconnected bodies
1132 // POEMS could do this itself
1133
1134 int *mark = new int[nbody];
1135 for (i = 0; i < nbody; i++) mark[i] = 1;
1136 for (i = 0; i < njoint; i++) {
1137 mark[jointbody[i][0]-1] = 0;
1138 mark[jointbody[i][1]-1] = 0;
1139 }
1140
1141 nfree = 0;
1142 for (i = 0; i < nbody; i++)
1143 if (mark[i]) nfree++;
1144 if (nfree) freelist = new int[nfree];
1145 else freelist = nullptr;
1146 nfree = 0;
1147 for (i = 0; i < nbody; i++)
1148 if (mark[i]) freelist[nfree++] = i + 1;
1149 delete [] mark;
1150
1151 // free memory local to this routine
1152
1153 memory->destroy(mylist);
1154 memory->destroy(jlist);
1155 memory->destroy(myjoint);
1156 }
1157
1158 /* ----------------------------------------------------------------------
1159 sort joint list (Numerical Recipes shell sort)
1160 sort criterion: sort on 1st body, if equal sort on 2nd body
1161 ------------------------------------------------------------------------- */
1162
sortlist(int n,tagint ** list)1163 void FixPOEMS::sortlist(int n, tagint **list)
1164 {
1165 int i,j,flag;
1166 tagint v0,v1,v2;
1167
1168 int inc = 1;
1169 while (inc <= n) inc = 3*inc + 1;
1170
1171 do {
1172 inc /= 3;
1173 for (i = inc+1; i <= n; i++) {
1174 v0 = list[i-1][0];
1175 v1 = list[i-1][1];
1176 v2 = list[i-1][2];
1177 j = i;
1178 flag = 0;
1179 if (list[j-inc-1][0] > v0 ||
1180 (list[j-inc-1][0] == v0 && list[j-inc-1][1] > v1)) flag = 1;
1181 while (flag) {
1182 list[j-1][0] = list[j-inc-1][0];
1183 list[j-1][1] = list[j-inc-1][1];
1184 list[j-1][2] = list[j-inc-1][2];
1185 j -= inc;
1186 if (j <= inc) break;
1187 flag = 0;
1188 if (list[j-inc-1][0] > v0 ||
1189 (list[j-inc-1][0] == v0 && list[j-inc-1][1] > v1)) flag = 1;
1190 }
1191 list[j-1][0] = v0;
1192 list[j-1][1] = v1;
1193 list[j-1][2] = v2;
1194 }
1195 } while (inc > 1);
1196 }
1197
1198 /* ----------------------------------------------------------------------
1199 check for cycles in list of joint connections between rigid bodies
1200 treat as graph: vertex = body, edge = joint between 2 bodies
1201 ------------------------------------------------------------------------- */
1202
loopcheck(int nvert,int nedge,tagint ** elist)1203 int FixPOEMS::loopcheck(int nvert, int nedge, tagint **elist)
1204 {
1205 int i,j,k;
1206
1207 // ecount[i] = # of vertices connected to vertex i via edge
1208 // elistfull[i][*] = list of vertices connected to vertex i
1209
1210 int *ecount = new int[nvert];
1211 for (i = 0; i < nvert; i++) ecount[i] = 0;
1212 for (i = 0; i < nedge; i++) {
1213 ecount[elist[i][0]]++;
1214 ecount[elist[i][1]]++;
1215 }
1216
1217 int emax = 0;
1218 for (i = 0; i < nvert; i++) emax = MAX(emax,ecount[i]);
1219
1220 int **elistfull;
1221 memory->create(elistfull,nvert,emax,"poems:elistfull");
1222 for (i = 0; i < nvert; i++) ecount[i] = 0;
1223 for (i = 0; i < nedge; i++) {
1224 elistfull[elist[i][0]][ecount[elist[i][0]]++] = elist[i][1];
1225 elistfull[elist[i][1]][ecount[elist[i][1]]++] = elist[i][0];
1226 }
1227
1228 // cycle detection algorithm
1229 // mark = 0/1 marking of each vertex, all initially unmarked
1230 // outer while loop:
1231 // if all vertices are marked, no cycles, exit loop
1232 // push an unmarked vertex on stack and mark it, parent is -1
1233 // while stack is not empty:
1234 // pop vertex I from stack
1235 // loop over vertices J connected to I via edge
1236 // if J is parent (vertex that pushed I on stack), skip it
1237 // else if J is marked, a cycle is found, return 1
1238 // else push J on stack and mark it, parent is I
1239 // increment ncluster each time stack empties since that is new cluster
1240
1241 int *parent = new int[nvert];
1242 int *mark = new int[nvert];
1243 for (i = 0; i < nvert; i++) mark[i] = 0;
1244
1245 int nstack = 0;
1246 int *stack = new int[nvert];
1247 ncluster = 0;
1248
1249 while (1) {
1250 for (i = 0; i < nvert; i++)
1251 if (mark[i] == 0) break;
1252 if (i == nvert) break;
1253 stack[nstack++] = i;
1254 mark[i] = 1;
1255 parent[i] = -1;
1256
1257 while (nstack) {
1258 i = stack[--nstack];
1259 for (k = 0; k < ecount[i]; k++) {
1260 j = elistfull[i][k];
1261 if (j == parent[i]) continue;
1262 if (mark[j]) return 1;
1263 stack[nstack++] = j;
1264 mark[j] = 1;
1265 parent[j] = i;
1266 }
1267 }
1268 ncluster++;
1269 }
1270
1271 // free memory local to this routine
1272
1273 delete [] ecount;
1274 memory->destroy(elistfull);
1275 delete [] parent;
1276 delete [] mark;
1277 delete [] stack;
1278
1279 return 0;
1280 }
1281
1282 /* ----------------------------------------------------------------------
1283 compute omega from angular momentum
1284 w = omega = angular velocity in space frame
1285 wbody = angular velocity in body frame
1286 set wbody component to 0.0 if inertia component is 0.0
1287 otherwise body can spin easily around that axis
1288 project space-frame angular momentum onto body axes
1289 and divide by principal moments
1290 ------------------------------------------------------------------------- */
1291
omega_from_mq(double * m,double * ex,double * ey,double * ez,double * inertia,double * w)1292 void FixPOEMS::omega_from_mq(double *m, double *ex, double *ey, double *ez,
1293 double *inertia, double *w)
1294 {
1295 double wbody[3];
1296
1297 if (inertia[0] == 0.0) wbody[0] = 0.0;
1298 else wbody[0] = (m[0]*ex[0] + m[1]*ex[1] + m[2]*ex[2]) / inertia[0];
1299 if (inertia[1] == 0.0) wbody[1] = 0.0;
1300 else wbody[1] = (m[0]*ey[0] + m[1]*ey[1] + m[2]*ey[2]) / inertia[1];
1301 if (inertia[2] == 0.0) wbody[2] = 0.0;
1302 else wbody[2] = (m[0]*ez[0] + m[1]*ez[1] + m[2]*ez[2]) / inertia[2];
1303
1304 w[0] = wbody[0]*ex[0] + wbody[1]*ey[0] + wbody[2]*ez[0];
1305 w[1] = wbody[0]*ex[1] + wbody[1]*ey[1] + wbody[2]*ez[1];
1306 w[2] = wbody[0]*ex[2] + wbody[1]*ey[2] + wbody[2]*ez[2];
1307 }
1308
1309 /* ----------------------------------------------------------------------
1310 set space-frame coords and velocity of each atom in each rigid body
1311 x = Q displace + Xcm, mapped back to periodic box
1312 v = Vcm + (W cross (x - Xcm))
1313 ------------------------------------------------------------------------- */
1314
set_xv()1315 void FixPOEMS::set_xv()
1316 {
1317 int ibody;
1318 int xbox,ybox,zbox;
1319 double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
1320 double vr[6];
1321
1322 imageint *image = atom->image;
1323 double **x = atom->x;
1324 double **v = atom->v;
1325 double **f = atom->f;
1326 double *mass = atom->mass;
1327 int *type = atom->type;
1328 int nlocal = atom->nlocal;
1329
1330 double xprd = domain->xprd;
1331 double yprd = domain->yprd;
1332 double zprd = domain->zprd;
1333
1334 // set x and v of each atom
1335 // only set joint atoms for 1st rigid body they belong to
1336
1337 for (int i = 0; i < nlocal; i++) {
1338 if (natom2body[i] == 0) continue;
1339 ibody = atom2body[i][0];
1340
1341 xbox = (image[i] & IMGMASK) - IMGMAX;
1342 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1343 zbox = (image[i] >> IMG2BITS) - IMGMAX;
1344
1345 // save old positions and velocities for virial
1346
1347 if (evflag) {
1348 x0 = x[i][0] + xbox*xprd;
1349 x1 = x[i][1] + ybox*yprd;
1350 x2 = x[i][2] + zbox*zprd;
1351
1352 v0 = v[i][0];
1353 v1 = v[i][1];
1354 v2 = v[i][2];
1355 }
1356
1357 // x = displacement from center-of-mass, based on body orientation
1358 // v = vcm + omega around center-of-mass
1359
1360 x[i][0] = ex_space[ibody][0]*displace[i][0] +
1361 ey_space[ibody][0]*displace[i][1] +
1362 ez_space[ibody][0]*displace[i][2];
1363 x[i][1] = ex_space[ibody][1]*displace[i][0] +
1364 ey_space[ibody][1]*displace[i][1] +
1365 ez_space[ibody][1]*displace[i][2];
1366 x[i][2] = ex_space[ibody][2]*displace[i][0] +
1367 ey_space[ibody][2]*displace[i][1] +
1368 ez_space[ibody][2]*displace[i][2];
1369
1370 v[i][0] = omega[ibody][1]*x[i][2] - omega[ibody][2]*x[i][1] +
1371 vcm[ibody][0];
1372 v[i][1] = omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2] +
1373 vcm[ibody][1];
1374 v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] +
1375 vcm[ibody][2];
1376
1377 // add center of mass to displacement
1378 // map back into periodic box via xbox,ybox,zbox
1379
1380 x[i][0] += xcm[ibody][0] - xbox*xprd;
1381 x[i][1] += xcm[ibody][1] - ybox*yprd;
1382 x[i][2] += xcm[ibody][2] - zbox*zprd;
1383
1384 // virial = unwrapped coords dotted into body constraint force
1385 // body constraint force = implied force due to v change minus f external
1386 // assume f does not include forces internal to body
1387 // 1/2 factor b/c final_integrate contributes other half
1388 // assume per-atom contribution is due to constraint force on that atom
1389
1390 if (evflag) {
1391 massone = mass[type[i]];
1392 fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
1393 fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
1394 fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
1395
1396 vr[0] = 0.5*fc0*x0;
1397 vr[1] = 0.5*fc1*x1;
1398 vr[2] = 0.5*fc2*x2;
1399 vr[3] = 0.5*fc1*x0;
1400 vr[4] = 0.5*fc2*x0;
1401 vr[5] = 0.5*fc2*x1;
1402
1403 v_tally(1,&i,1.0,vr);
1404 }
1405 }
1406 }
1407
1408 /* ----------------------------------------------------------------------
1409 set space-frame velocity of each atom in a rigid body
1410 v = Vcm + (W cross (x - Xcm))
1411 ------------------------------------------------------------------------- */
1412
set_v()1413 void FixPOEMS::set_v()
1414 {
1415 int ibody;
1416 int xbox,ybox,zbox;
1417 double dx,dy,dz;
1418 double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
1419 double vr[6];
1420
1421 double *mass = atom->mass;
1422 double **f = atom->f;
1423 double **x = atom->x;
1424 double **v = atom->v;
1425 int *type = atom->type;
1426 imageint *image = atom->image;
1427 int nlocal = atom->nlocal;
1428
1429 double xprd = domain->xprd;
1430 double yprd = domain->yprd;
1431 double zprd = domain->zprd;
1432
1433 // set v of each atom
1434 // only set joint atoms for 1st rigid body they belong to
1435
1436 for (int i = 0; i < nlocal; i++) {
1437 if (natom2body[i] == 0) continue;
1438 ibody = atom2body[i][0];
1439
1440 dx = ex_space[ibody][0]*displace[i][0] +
1441 ey_space[ibody][0]*displace[i][1] +
1442 ez_space[ibody][0]*displace[i][2];
1443 dy = ex_space[ibody][1]*displace[i][0] +
1444 ey_space[ibody][1]*displace[i][1] +
1445 ez_space[ibody][1]*displace[i][2];
1446 dz = ex_space[ibody][2]*displace[i][0] +
1447 ey_space[ibody][2]*displace[i][1] +
1448 ez_space[ibody][2]*displace[i][2];
1449
1450 // save old velocities for virial
1451
1452 if (evflag) {
1453 v0 = v[i][0];
1454 v1 = v[i][1];
1455 v2 = v[i][2];
1456 }
1457
1458 v[i][0] = omega[ibody][1]*dz - omega[ibody][2]*dy + vcm[ibody][0];
1459 v[i][1] = omega[ibody][2]*dx - omega[ibody][0]*dz + vcm[ibody][1];
1460 v[i][2] = omega[ibody][0]*dy - omega[ibody][1]*dx + vcm[ibody][2];
1461
1462 // virial = unwrapped coords dotted into body constraint force
1463 // body constraint force = implied force due to v change minus f external
1464 // assume f does not include forces internal to body
1465 // 1/2 factor b/c initial_integrate contributes other half
1466 // assume per-atom contribution is due to constraint force on that atom
1467
1468 if (evflag) {
1469 massone = mass[type[i]];
1470 fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
1471 fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
1472 fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
1473
1474 xbox = (image[i] & IMGMASK) - IMGMAX;
1475 ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1476 zbox = (image[i] >> IMG2BITS) - IMGMAX;
1477
1478 x0 = x[i][0] + xbox*xprd;
1479 x1 = x[i][1] + ybox*yprd;
1480 x2 = x[i][2] + zbox*zprd;
1481
1482 vr[0] = 0.5*fc0*x0;
1483 vr[1] = 0.5*fc1*x1;
1484 vr[2] = 0.5*fc2*x2;
1485 vr[3] = 0.5*fc1*x0;
1486 vr[4] = 0.5*fc2*x0;
1487 vr[5] = 0.5*fc2*x1;
1488
1489 v_tally(1,&i,1.0,vr);
1490 }
1491 }
1492 }
1493
1494 /* ----------------------------------------------------------------------
1495 allocate local atom-based arrays
1496 ------------------------------------------------------------------------- */
1497
grow_arrays(int nmax)1498 void FixPOEMS::grow_arrays(int nmax)
1499 {
1500 memory->grow(natom2body,nmax,"fix_poems:natom2body");
1501 memory->grow(atom2body,nmax,MAXBODY,"fix_poems:atom2body");
1502 memory->grow(displace,nmax,3,"fix_poems:displace");
1503 }
1504
1505 /* ----------------------------------------------------------------------
1506 copy values within local atom-based arrays
1507 ------------------------------------------------------------------------- */
1508
copy_arrays(int i,int j,int)1509 void FixPOEMS::copy_arrays(int i, int j, int /* delflag */)
1510 {
1511 natom2body[j] = natom2body[i];
1512 for (int k = 0; k < natom2body[j]; k++) atom2body[j][k] = atom2body[i][k];
1513 displace[j][0] = displace[i][0];
1514 displace[j][1] = displace[i][1];
1515 displace[j][2] = displace[i][2];
1516 }
1517
1518 /* ----------------------------------------------------------------------
1519 memory usage of local atom-based arrays
1520 ------------------------------------------------------------------------- */
1521
memory_usage()1522 double FixPOEMS::memory_usage()
1523 {
1524 int nmax = atom->nmax;
1525 double bytes = (double)nmax * sizeof(int);
1526 bytes += (double)nmax*MAXBODY * sizeof(int);
1527 bytes += (double)nmax*3 * sizeof(double);
1528 return bytes;
1529 }
1530
1531 /* ----------------------------------------------------------------------
1532 pack values in local atom-based arrays for exchange with another proc
1533 ------------------------------------------------------------------------- */
1534
pack_exchange(int i,double * buf)1535 int FixPOEMS::pack_exchange(int i, double *buf)
1536 {
1537 int m = 0;
1538 buf[m++] = static_cast<double> (natom2body[i]);
1539 for (int j = 0; j < natom2body[i]; j++)
1540 buf[m++] = static_cast<double> (atom2body[i][j]);
1541 buf[m++] = displace[i][0];
1542 buf[m++] = displace[i][1];
1543 buf[m++] = displace[i][2];
1544 return m;
1545 }
1546
1547 /* ----------------------------------------------------------------------
1548 unpack values in local atom-based arrays from exchange with another proc
1549 ------------------------------------------------------------------------- */
1550
unpack_exchange(int nlocal,double * buf)1551 int FixPOEMS::unpack_exchange(int nlocal, double *buf)
1552 {
1553 int m = 0;
1554 natom2body[nlocal] = static_cast<int> (buf[m++]);
1555 for (int i = 0; i < natom2body[nlocal]; i++)
1556 atom2body[nlocal][i] = static_cast<int> (buf[m++]);
1557 displace[nlocal][0] = buf[m++];
1558 displace[nlocal][1] = buf[m++];
1559 displace[nlocal][2] = buf[m++];
1560 return m;
1561 }
1562
1563 /* ---------------------------------------------------------------------- */
1564
modify_param(int narg,char ** arg)1565 int FixPOEMS::modify_param(int narg, char **arg)
1566 {
1567 if (strcmp(arg[0],"bodyforces") == 0) {
1568 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
1569 if (strcmp(arg[1],"early") == 0) earlyflag = 1;
1570 else if (strcmp(arg[1],"late") == 0) earlyflag = 0;
1571 else error->all(FLERR,"Illegal fix_modify command");
1572 return 2;
1573 }
1574
1575 return 0;
1576 }
1577
1578 /* ---------------------------------------------------------------------- */
1579
reset_dt()1580 void FixPOEMS::reset_dt()
1581 {
1582 dtv = update->dt;
1583 dtf = 0.5 * update->dt * force->ftm2v;
1584 dthalf = 0.5 * update->dt;
1585 }
1586