1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/, Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13
14 /* ----------------------------------------------------------------------
15 Contributing authors: Leo Silbert (SNL), Gary Grest (SNL)
16 ------------------------------------------------------------------------- */
17
18 #include "pair_gran_hooke_history.h"
19
20 #include "atom.h"
21 #include "comm.h"
22 #include "error.h"
23 #include "fix.h"
24 #include "fix_dummy.h"
25 #include "fix_neigh_history.h"
26 #include "force.h"
27 #include "memory.h"
28 #include "modify.h"
29 #include "neigh_list.h"
30 #include "neigh_request.h"
31 #include "neighbor.h"
32 #include "update.h"
33
34 #include <cmath>
35 #include <cstring>
36
37 using namespace LAMMPS_NS;
38
39 /* ---------------------------------------------------------------------- */
40
PairGranHookeHistory(LAMMPS * lmp)41 PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp)
42 {
43 single_enable = 1;
44 no_virial_fdotr_compute = 1;
45 centroidstressflag = CENTROID_NOTAVAIL;
46 finitecutflag = 1;
47 history = 1;
48 size_history = 3;
49
50 single_extra = 10;
51 svector = new double[10];
52
53 neighprev = 0;
54
55 nmax = 0;
56 mass_rigid = nullptr;
57
58 // set comm size needed by this Pair if used with fix rigid
59
60 comm_forward = 1;
61
62 // keep default behavior of history[i][j] = -history[j][i]
63
64 nondefault_history_transfer = 0;
65
66 // create dummy fix as placeholder for FixNeighHistory
67 // this is so final order of Modify:fix will conform to input script
68
69 fix_history = nullptr;
70 fix_dummy = (FixDummy *) modify->add_fix("NEIGH_HISTORY_HH_DUMMY" + std::to_string(instance_me) +
71 " all DUMMY");
72 }
73
74 /* ---------------------------------------------------------------------- */
75
~PairGranHookeHistory()76 PairGranHookeHistory::~PairGranHookeHistory()
77 {
78 if (copymode) return;
79
80 delete[] svector;
81
82 if (!fix_history)
83 modify->delete_fix("NEIGH_HISTORY_HH_DUMMY" + std::to_string(instance_me));
84 else
85 modify->delete_fix("NEIGH_HISTORY_HH" + std::to_string(instance_me));
86
87 if (allocated) {
88 memory->destroy(setflag);
89 memory->destroy(cutsq);
90
91 delete[] onerad_dynamic;
92 delete[] onerad_frozen;
93 delete[] maxrad_dynamic;
94 delete[] maxrad_frozen;
95 }
96
97 memory->destroy(mass_rigid);
98 }
99
100 /* ---------------------------------------------------------------------- */
101
compute(int eflag,int vflag)102 void PairGranHookeHistory::compute(int eflag, int vflag)
103 {
104 int i, j, ii, jj, inum, jnum;
105 double xtmp, ytmp, ztmp, delx, dely, delz, fx, fy, fz;
106 double radi, radj, radsum, rsq, r, rinv, rsqinv;
107 double vr1, vr2, vr3, vnnr, vn1, vn2, vn3, vt1, vt2, vt3;
108 double wr1, wr2, wr3;
109 double vtr1, vtr2, vtr3, vrel;
110 double mi, mj, meff, damp, ccel, tor1, tor2, tor3;
111 double fn, fs, fs1, fs2, fs3;
112 double shrmag, rsht;
113 int *ilist, *jlist, *numneigh, **firstneigh;
114 int *touch, **firsttouch;
115 double *shear, *allshear, **firstshear;
116
117 ev_init(eflag, vflag);
118
119 int shearupdate = 1;
120 if (update->setupflag) shearupdate = 0;
121
122 // update rigid body info for owned & ghost atoms if using FixRigid masses
123 // body[i] = which body atom I is in, -1 if none
124 // mass_body = mass of each rigid body
125
126 if (fix_rigid && neighbor->ago == 0) {
127 int tmp;
128 int *body = (int *) fix_rigid->extract("body", tmp);
129 double *mass_body = (double *) fix_rigid->extract("masstotal", tmp);
130 if (atom->nmax > nmax) {
131 memory->destroy(mass_rigid);
132 nmax = atom->nmax;
133 memory->create(mass_rigid, nmax, "pair:mass_rigid");
134 }
135 int nlocal = atom->nlocal;
136 for (i = 0; i < nlocal; i++)
137 if (body[i] >= 0)
138 mass_rigid[i] = mass_body[body[i]];
139 else
140 mass_rigid[i] = 0.0;
141 comm->forward_comm_pair(this);
142 }
143
144 double **x = atom->x;
145 double **v = atom->v;
146 double **f = atom->f;
147 double **omega = atom->omega;
148 double **torque = atom->torque;
149 double *radius = atom->radius;
150 double *rmass = atom->rmass;
151 int *mask = atom->mask;
152 int nlocal = atom->nlocal;
153 int newton_pair = force->newton_pair;
154
155 inum = list->inum;
156 ilist = list->ilist;
157 numneigh = list->numneigh;
158 firstneigh = list->firstneigh;
159 firsttouch = fix_history->firstflag;
160 firstshear = fix_history->firstvalue;
161
162 // loop over neighbors of my atoms
163
164 for (ii = 0; ii < inum; ii++) {
165 i = ilist[ii];
166 xtmp = x[i][0];
167 ytmp = x[i][1];
168 ztmp = x[i][2];
169 radi = radius[i];
170 touch = firsttouch[i];
171 allshear = firstshear[i];
172 jlist = firstneigh[i];
173 jnum = numneigh[i];
174
175 for (jj = 0; jj < jnum; jj++) {
176 j = jlist[jj];
177 j &= NEIGHMASK;
178
179 delx = xtmp - x[j][0];
180 dely = ytmp - x[j][1];
181 delz = ztmp - x[j][2];
182 rsq = delx * delx + dely * dely + delz * delz;
183 radj = radius[j];
184 radsum = radi + radj;
185
186 if (rsq >= radsum * radsum) {
187
188 // unset non-touching neighbors
189
190 touch[jj] = 0;
191 shear = &allshear[3 * jj];
192 shear[0] = 0.0;
193 shear[1] = 0.0;
194 shear[2] = 0.0;
195
196 } else {
197 r = sqrt(rsq);
198 rinv = 1.0 / r;
199 rsqinv = 1.0 / rsq;
200
201 // relative translational velocity
202
203 vr1 = v[i][0] - v[j][0];
204 vr2 = v[i][1] - v[j][1];
205 vr3 = v[i][2] - v[j][2];
206
207 // normal component
208
209 vnnr = vr1 * delx + vr2 * dely + vr3 * delz;
210 vn1 = delx * vnnr * rsqinv;
211 vn2 = dely * vnnr * rsqinv;
212 vn3 = delz * vnnr * rsqinv;
213
214 // tangential component
215
216 vt1 = vr1 - vn1;
217 vt2 = vr2 - vn2;
218 vt3 = vr3 - vn3;
219
220 // relative rotational velocity
221
222 wr1 = (radi * omega[i][0] + radj * omega[j][0]) * rinv;
223 wr2 = (radi * omega[i][1] + radj * omega[j][1]) * rinv;
224 wr3 = (radi * omega[i][2] + radj * omega[j][2]) * rinv;
225
226 // meff = effective mass of pair of particles
227 // if I or J part of rigid body, use body mass
228 // if I or J is frozen, meff is other particle
229
230 mi = rmass[i];
231 mj = rmass[j];
232 if (fix_rigid) {
233 if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
234 if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
235 }
236
237 meff = mi * mj / (mi + mj);
238 if (mask[i] & freeze_group_bit) meff = mj;
239 if (mask[j] & freeze_group_bit) meff = mi;
240
241 // normal forces = Hookian contact + normal velocity damping
242
243 damp = meff * gamman * vnnr * rsqinv;
244 ccel = kn * (radsum - r) * rinv - damp;
245 if (limit_damping && (ccel < 0.0)) ccel = 0.0;
246
247 // relative velocities
248
249 vtr1 = vt1 - (delz * wr2 - dely * wr3);
250 vtr2 = vt2 - (delx * wr3 - delz * wr1);
251 vtr3 = vt3 - (dely * wr1 - delx * wr2);
252 vrel = vtr1 * vtr1 + vtr2 * vtr2 + vtr3 * vtr3;
253 vrel = sqrt(vrel);
254
255 // shear history effects
256
257 touch[jj] = 1;
258 shear = &allshear[3 * jj];
259
260 if (shearupdate) {
261 shear[0] += vtr1 * dt;
262 shear[1] += vtr2 * dt;
263 shear[2] += vtr3 * dt;
264 }
265 shrmag = sqrt(shear[0] * shear[0] + shear[1] * shear[1] + shear[2] * shear[2]);
266
267 // rotate shear displacements
268
269 rsht = shear[0] * delx + shear[1] * dely + shear[2] * delz;
270 rsht *= rsqinv;
271 if (shearupdate) {
272 shear[0] -= rsht * delx;
273 shear[1] -= rsht * dely;
274 shear[2] -= rsht * delz;
275 }
276
277 // tangential forces = shear + tangential velocity damping
278
279 fs1 = -(kt * shear[0] + meff * gammat * vtr1);
280 fs2 = -(kt * shear[1] + meff * gammat * vtr2);
281 fs3 = -(kt * shear[2] + meff * gammat * vtr3);
282
283 // rescale frictional displacements and forces if needed
284
285 fs = sqrt(fs1 * fs1 + fs2 * fs2 + fs3 * fs3);
286 fn = xmu * fabs(ccel * r);
287
288 if (fs > fn) {
289 if (shrmag != 0.0) {
290 shear[0] =
291 (fn / fs) * (shear[0] + meff * gammat * vtr1 / kt) - meff * gammat * vtr1 / kt;
292 shear[1] =
293 (fn / fs) * (shear[1] + meff * gammat * vtr2 / kt) - meff * gammat * vtr2 / kt;
294 shear[2] =
295 (fn / fs) * (shear[2] + meff * gammat * vtr3 / kt) - meff * gammat * vtr3 / kt;
296 fs1 *= fn / fs;
297 fs2 *= fn / fs;
298 fs3 *= fn / fs;
299 } else
300 fs1 = fs2 = fs3 = 0.0;
301 }
302
303 // forces & torques
304
305 fx = delx * ccel + fs1;
306 fy = dely * ccel + fs2;
307 fz = delz * ccel + fs3;
308 f[i][0] += fx;
309 f[i][1] += fy;
310 f[i][2] += fz;
311
312 tor1 = rinv * (dely * fs3 - delz * fs2);
313 tor2 = rinv * (delz * fs1 - delx * fs3);
314 tor3 = rinv * (delx * fs2 - dely * fs1);
315 torque[i][0] -= radi * tor1;
316 torque[i][1] -= radi * tor2;
317 torque[i][2] -= radi * tor3;
318
319 if (newton_pair || j < nlocal) {
320 f[j][0] -= fx;
321 f[j][1] -= fy;
322 f[j][2] -= fz;
323 torque[j][0] -= radj * tor1;
324 torque[j][1] -= radj * tor2;
325 torque[j][2] -= radj * tor3;
326 }
327
328 if (evflag) ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, fx, fy, fz, delx, dely, delz);
329 }
330 }
331 }
332
333 if (vflag_fdotr) virial_fdotr_compute();
334 }
335
336 /* ----------------------------------------------------------------------
337 allocate all arrays
338 ------------------------------------------------------------------------- */
339
allocate()340 void PairGranHookeHistory::allocate()
341 {
342 allocated = 1;
343 int n = atom->ntypes;
344
345 memory->create(setflag, n + 1, n + 1, "pair:setflag");
346 for (int i = 1; i <= n; i++)
347 for (int j = i; j <= n; j++) setflag[i][j] = 0;
348
349 memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
350
351 onerad_dynamic = new double[n + 1];
352 onerad_frozen = new double[n + 1];
353 maxrad_dynamic = new double[n + 1];
354 maxrad_frozen = new double[n + 1];
355 }
356
357 /* ----------------------------------------------------------------------
358 global settings
359 ------------------------------------------------------------------------- */
360
settings(int narg,char ** arg)361 void PairGranHookeHistory::settings(int narg, char **arg)
362 {
363 if (narg != 6 && narg != 7) error->all(FLERR, "Illegal pair_style command");
364
365 kn = utils::numeric(FLERR, arg[0], false, lmp);
366 if (strcmp(arg[1], "NULL") == 0)
367 kt = kn * 2.0 / 7.0;
368 else
369 kt = utils::numeric(FLERR, arg[1], false, lmp);
370
371 gamman = utils::numeric(FLERR, arg[2], false, lmp);
372 if (strcmp(arg[3], "NULL") == 0)
373 gammat = 0.5 * gamman;
374 else
375 gammat = utils::numeric(FLERR, arg[3], false, lmp);
376
377 xmu = utils::numeric(FLERR, arg[4], false, lmp);
378 dampflag = utils::inumeric(FLERR, arg[5], false, lmp);
379 if (dampflag == 0) gammat = 0.0;
380
381 limit_damping = 0;
382 if (narg == 7) {
383 if (strcmp(arg[6], "limit_damping") == 0)
384 limit_damping = 1;
385 else
386 error->all(FLERR, "Illegal pair_style command");
387 }
388
389 if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 || xmu < 0.0 || xmu > 10000.0 ||
390 dampflag < 0 || dampflag > 1)
391 error->all(FLERR, "Illegal pair_style command");
392 }
393
394 /* ----------------------------------------------------------------------
395 set coeffs for one or more type pairs
396 ------------------------------------------------------------------------- */
397
coeff(int narg,char ** arg)398 void PairGranHookeHistory::coeff(int narg, char **arg)
399 {
400 if (narg > 2) error->all(FLERR, "Incorrect args for pair coefficients");
401 if (!allocated) allocate();
402
403 int ilo, ihi, jlo, jhi;
404 utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
405 utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
406
407 int count = 0;
408 for (int i = ilo; i <= ihi; i++) {
409 for (int j = MAX(jlo, i); j <= jhi; j++) {
410 setflag[i][j] = 1;
411 count++;
412 }
413 }
414
415 if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
416 }
417
418 /* ----------------------------------------------------------------------
419 init specific to this pair style
420 ------------------------------------------------------------------------- */
421
init_style()422 void PairGranHookeHistory::init_style()
423 {
424 int i;
425
426 // error and warning checks
427
428 if (!atom->radius_flag || !atom->rmass_flag)
429 error->all(FLERR, "Pair granular requires atom attributes radius, rmass");
430 if (comm->ghost_velocity == 0)
431 error->all(FLERR, "Pair granular requires ghost atoms store velocity");
432
433 // need a granular neigh list
434
435 int irequest = neighbor->request(this, instance_me);
436 neighbor->requests[irequest]->size = 1;
437 if (history) neighbor->requests[irequest]->history = 1;
438
439 dt = update->dt;
440
441 // if history is stored and first init, create Fix to store history
442 // it replaces FixDummy, created in the constructor
443 // this is so its order in the fix list is preserved
444
445 if (history && (fix_history == nullptr)) {
446 auto cmd = fmt::format("NEIGH_HISTORY_HH{} all NEIGH_HISTORY {}", instance_me, size_history);
447 fix_history = (FixNeighHistory *) modify->replace_fix(
448 "NEIGH_HISTORY_HH_DUMMY" + std::to_string(instance_me), cmd, 1);
449 fix_history->pair = this;
450 }
451
452 // check for FixFreeze and set freeze_group_bit
453
454 int ifreeze = modify->find_fix_by_style("^freeze");
455 if (ifreeze < 0)
456 freeze_group_bit = 0;
457 else
458 freeze_group_bit = modify->fix[ifreeze]->groupbit;
459
460 // check for FixRigid so can extract rigid body masses
461 // FIXME: this only catches the first rigid fix, there may be multiple.
462
463 fix_rigid = nullptr;
464 for (i = 0; i < modify->nfix; i++)
465 if (modify->fix[i]->rigid_flag) break;
466 if (i < modify->nfix) fix_rigid = modify->fix[i];
467
468 // check for FixPour and FixDeposit so can extract particle radii
469
470 int ipour = modify->find_fix_by_style("^pour");
471 int idep = modify->find_fix_by_style("^deposit");
472
473 // set maxrad_dynamic and maxrad_frozen for each type
474 // include future FixPour and FixDeposit particles as dynamic
475
476 int itype;
477 for (i = 1; i <= atom->ntypes; i++) {
478 onerad_dynamic[i] = onerad_frozen[i] = 0.0;
479 if (ipour >= 0) {
480 itype = i;
481 onerad_dynamic[i] = *((double *) modify->fix[ipour]->extract("radius", itype));
482 }
483 if (idep >= 0) {
484 itype = i;
485 onerad_dynamic[i] = *((double *) modify->fix[idep]->extract("radius", itype));
486 }
487 }
488
489 double *radius = atom->radius;
490 int *mask = atom->mask;
491 int *type = atom->type;
492 int nlocal = atom->nlocal;
493
494 for (i = 0; i < nlocal; i++)
495 if (mask[i] & freeze_group_bit)
496 onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]], radius[i]);
497 else
498 onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]], radius[i]);
499
500 MPI_Allreduce(&onerad_dynamic[1], &maxrad_dynamic[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
501 MPI_Allreduce(&onerad_frozen[1], &maxrad_frozen[1], atom->ntypes, MPI_DOUBLE, MPI_MAX, world);
502
503 // set fix which stores history info
504
505 if (history) {
506 int ifix = modify->find_fix("NEIGH_HISTORY_HH" + std::to_string(instance_me));
507 if (ifix < 0) error->all(FLERR, "Could not find pair fix neigh history ID");
508 fix_history = (FixNeighHistory *) modify->fix[ifix];
509 }
510 }
511
512 /* ----------------------------------------------------------------------
513 init for one type pair i,j and corresponding j,i
514 ------------------------------------------------------------------------- */
515
init_one(int i,int j)516 double PairGranHookeHistory::init_one(int i, int j)
517 {
518 if (!allocated) allocate();
519
520 // cutoff = sum of max I,J radii for
521 // dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
522
523 double cutoff = maxrad_dynamic[i] + maxrad_dynamic[j];
524 cutoff = MAX(cutoff, maxrad_frozen[i] + maxrad_dynamic[j]);
525 cutoff = MAX(cutoff, maxrad_dynamic[i] + maxrad_frozen[j]);
526 return cutoff;
527 }
528
529 /* ----------------------------------------------------------------------
530 proc 0 writes to restart file
531 ------------------------------------------------------------------------- */
532
write_restart(FILE * fp)533 void PairGranHookeHistory::write_restart(FILE *fp)
534 {
535 write_restart_settings(fp);
536
537 int i, j;
538 for (i = 1; i <= atom->ntypes; i++)
539 for (j = i; j <= atom->ntypes; j++) fwrite(&setflag[i][j], sizeof(int), 1, fp);
540 }
541
542 /* ----------------------------------------------------------------------
543 proc 0 reads from restart file, bcasts
544 ------------------------------------------------------------------------- */
545
read_restart(FILE * fp)546 void PairGranHookeHistory::read_restart(FILE *fp)
547 {
548 read_restart_settings(fp);
549 allocate();
550
551 int i, j;
552 int me = comm->me;
553 for (i = 1; i <= atom->ntypes; i++)
554 for (j = i; j <= atom->ntypes; j++) {
555 if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
556 MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
557 }
558 }
559
560 /* ----------------------------------------------------------------------
561 proc 0 writes to restart file
562 ------------------------------------------------------------------------- */
563
write_restart_settings(FILE * fp)564 void PairGranHookeHistory::write_restart_settings(FILE *fp)
565 {
566 fwrite(&kn, sizeof(double), 1, fp);
567 fwrite(&kt, sizeof(double), 1, fp);
568 fwrite(&gamman, sizeof(double), 1, fp);
569 fwrite(&gammat, sizeof(double), 1, fp);
570 fwrite(&xmu, sizeof(double), 1, fp);
571 fwrite(&dampflag, sizeof(int), 1, fp);
572 }
573
574 /* ----------------------------------------------------------------------
575 proc 0 reads from restart file, bcasts
576 ------------------------------------------------------------------------- */
577
read_restart_settings(FILE * fp)578 void PairGranHookeHistory::read_restart_settings(FILE *fp)
579 {
580 if (comm->me == 0) {
581 utils::sfread(FLERR, &kn, sizeof(double), 1, fp, nullptr, error);
582 utils::sfread(FLERR, &kt, sizeof(double), 1, fp, nullptr, error);
583 utils::sfread(FLERR, &gamman, sizeof(double), 1, fp, nullptr, error);
584 utils::sfread(FLERR, &gammat, sizeof(double), 1, fp, nullptr, error);
585 utils::sfread(FLERR, &xmu, sizeof(double), 1, fp, nullptr, error);
586 utils::sfread(FLERR, &dampflag, sizeof(int), 1, fp, nullptr, error);
587 }
588 MPI_Bcast(&kn, 1, MPI_DOUBLE, 0, world);
589 MPI_Bcast(&kt, 1, MPI_DOUBLE, 0, world);
590 MPI_Bcast(&gamman, 1, MPI_DOUBLE, 0, world);
591 MPI_Bcast(&gammat, 1, MPI_DOUBLE, 0, world);
592 MPI_Bcast(&xmu, 1, MPI_DOUBLE, 0, world);
593 MPI_Bcast(&dampflag, 1, MPI_INT, 0, world);
594 }
595
596 /* ---------------------------------------------------------------------- */
597
reset_dt()598 void PairGranHookeHistory::reset_dt()
599 {
600 dt = update->dt;
601 }
602
603 /* ---------------------------------------------------------------------- */
604
single(int i,int j,int,int,double rsq,double,double,double & fforce)605 double PairGranHookeHistory::single(int i, int j, int /*itype*/, int /*jtype*/, double rsq,
606 double /*factor_coul*/, double /*factor_lj*/, double &fforce)
607 {
608 double radi, radj, radsum;
609 double r, rinv, rsqinv, delx, dely, delz;
610 double vr1, vr2, vr3, vnnr, vn1, vn2, vn3, vt1, vt2, vt3, wr1, wr2, wr3;
611 double mi, mj, meff, damp, ccel;
612 double vtr1, vtr2, vtr3, vrel, shrmag, rsht;
613 double fs1, fs2, fs3, fs, fn;
614
615 double *radius = atom->radius;
616 radi = radius[i];
617 radj = radius[j];
618 radsum = radi + radj;
619
620 if (rsq >= radsum * radsum) {
621 fforce = 0.0;
622 for (int m = 0; m < single_extra; m++) svector[m] = 0.0;
623 return 0.0;
624 }
625
626 r = sqrt(rsq);
627 rinv = 1.0 / r;
628 rsqinv = 1.0 / rsq;
629
630 // relative translational velocity
631
632 double **v = atom->v;
633 vr1 = v[i][0] - v[j][0];
634 vr2 = v[i][1] - v[j][1];
635 vr3 = v[i][2] - v[j][2];
636
637 // normal component
638
639 double **x = atom->x;
640 delx = x[i][0] - x[j][0];
641 dely = x[i][1] - x[j][1];
642 delz = x[i][2] - x[j][2];
643
644 vnnr = vr1 * delx + vr2 * dely + vr3 * delz;
645 vn1 = delx * vnnr * rsqinv;
646 vn2 = dely * vnnr * rsqinv;
647 vn3 = delz * vnnr * rsqinv;
648
649 // tangential component
650
651 vt1 = vr1 - vn1;
652 vt2 = vr2 - vn2;
653 vt3 = vr3 - vn3;
654
655 // relative rotational velocity
656
657 double **omega = atom->omega;
658 wr1 = (radi * omega[i][0] + radj * omega[j][0]) * rinv;
659 wr2 = (radi * omega[i][1] + radj * omega[j][1]) * rinv;
660 wr3 = (radi * omega[i][2] + radj * omega[j][2]) * rinv;
661
662 // meff = effective mass of pair of particles
663 // if I or J part of rigid body, use body mass
664 // if I or J is frozen, meff is other particle
665
666 double *rmass = atom->rmass;
667 int *mask = atom->mask;
668
669 mi = rmass[i];
670 mj = rmass[j];
671 if (fix_rigid) {
672 // NOTE: insure mass_rigid is current for owned+ghost atoms?
673 if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
674 if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
675 }
676
677 meff = mi * mj / (mi + mj);
678 if (mask[i] & freeze_group_bit) meff = mj;
679 if (mask[j] & freeze_group_bit) meff = mi;
680
681 // normal forces = Hookian contact + normal velocity damping
682
683 damp = meff * gamman * vnnr * rsqinv;
684 ccel = kn * (radsum - r) * rinv - damp;
685 if (limit_damping && (ccel < 0.0)) ccel = 0.0;
686
687 // relative velocities
688
689 vtr1 = vt1 - (delz * wr2 - dely * wr3);
690 vtr2 = vt2 - (delx * wr3 - delz * wr1);
691 vtr3 = vt3 - (dely * wr1 - delx * wr2);
692 vrel = vtr1 * vtr1 + vtr2 * vtr2 + vtr3 * vtr3;
693 vrel = sqrt(vrel);
694
695 // shear history effects
696 // neighprev = index of found neigh on previous call
697 // search entire jnum list of neighbors of I for neighbor J
698 // start from neighprev, since will typically be next neighbor
699 // reset neighprev to 0 as necessary
700
701 int jnum = list->numneigh[i];
702 int *jlist = list->firstneigh[i];
703 double *allshear = fix_history->firstvalue[i];
704
705 for (int jj = 0; jj < jnum; jj++) {
706 neighprev++;
707 if (neighprev >= jnum) neighprev = 0;
708 if (jlist[neighprev] == j) break;
709 }
710
711 double *shear = &allshear[3 * neighprev];
712 shrmag = sqrt(shear[0] * shear[0] + shear[1] * shear[1] + shear[2] * shear[2]);
713
714 // rotate shear displacements
715
716 rsht = shear[0] * delx + shear[1] * dely + shear[2] * delz;
717 rsht *= rsqinv;
718
719 // tangential forces = shear + tangential velocity damping
720
721 fs1 = -(kt * shear[0] + meff * gammat * vtr1);
722 fs2 = -(kt * shear[1] + meff * gammat * vtr2);
723 fs3 = -(kt * shear[2] + meff * gammat * vtr3);
724
725 // rescale frictional displacements and forces if needed
726
727 fs = sqrt(fs1 * fs1 + fs2 * fs2 + fs3 * fs3);
728 fn = xmu * fabs(ccel * r);
729
730 if (fs > fn) {
731 if (shrmag != 0.0) {
732 fs1 *= fn / fs;
733 fs2 *= fn / fs;
734 fs3 *= fn / fs;
735 fs *= fn / fs;
736 } else
737 fs1 = fs2 = fs3 = fs = 0.0;
738 }
739
740 // set force and return no energy
741
742 fforce = ccel;
743
744 // set single_extra quantities
745
746 svector[0] = fs1;
747 svector[1] = fs2;
748 svector[2] = fs3;
749 svector[3] = fs;
750 svector[4] = vn1;
751 svector[5] = vn2;
752 svector[6] = vn3;
753 svector[7] = vt1;
754 svector[8] = vt2;
755 svector[9] = vt3;
756
757 return 0.0;
758 }
759
760 /* ---------------------------------------------------------------------- */
761
pack_forward_comm(int n,int * list,double * buf,int,int *)762 int PairGranHookeHistory::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
763 int * /*pbc*/)
764 {
765 int i, j, m;
766
767 m = 0;
768 for (i = 0; i < n; i++) {
769 j = list[i];
770 buf[m++] = mass_rigid[j];
771 }
772 return m;
773 }
774
775 /* ---------------------------------------------------------------------- */
776
unpack_forward_comm(int n,int first,double * buf)777 void PairGranHookeHistory::unpack_forward_comm(int n, int first, double *buf)
778 {
779 int i, m, last;
780
781 m = 0;
782 last = first + n;
783 for (i = first; i < last; i++) mass_rigid[i] = buf[m++];
784 }
785
786 /* ----------------------------------------------------------------------
787 memory usage of local atom-based arrays
788 ------------------------------------------------------------------------- */
789
memory_usage()790 double PairGranHookeHistory::memory_usage()
791 {
792 double bytes = (double) nmax * sizeof(double);
793 return bytes;
794 }
795
796 /* ----------------------------------------------------------------------
797 self-interaction range of particle
798 ------------------------------------------------------------------------- */
799
atom2cut(int i)800 double PairGranHookeHistory::atom2cut(int i)
801 {
802 double cut = atom->radius[i] * 2;
803 return cut;
804 }
805
806 /* ----------------------------------------------------------------------
807 maximum interaction range for two finite particles
808 ------------------------------------------------------------------------- */
809
radii2cut(double r1,double r2)810 double PairGranHookeHistory::radii2cut(double r1, double r2)
811 {
812 double cut = r1 + r2;
813 return cut;
814 }
815