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 Contributing authors: Leo Silbert (SNL), Gary Grest (SNL)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_gran_hooke.h"
20 #include <cmath>
21 #include "atom.h"
22 #include "force.h"
23 #include "fix.h"
24 #include "neighbor.h"
25 #include "neigh_list.h"
26 #include "comm.h"
27 #include "memory.h"
28
29 using namespace LAMMPS_NS;
30
31 /* ---------------------------------------------------------------------- */
32
PairGranHooke(LAMMPS * lmp)33 PairGranHooke::PairGranHooke(LAMMPS *lmp) : PairGranHookeHistory(lmp)
34 {
35 no_virial_fdotr_compute = 0;
36 history = 0;
37 }
38
39 /* ---------------------------------------------------------------------- */
40
compute(int eflag,int vflag)41 void PairGranHooke::compute(int eflag, int vflag)
42 {
43 int i,j,ii,jj,inum,jnum;
44 double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz;
45 double radi,radj,radsum,rsq,r,rinv,rsqinv;
46 double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
47 double wr1,wr2,wr3;
48 double vtr1,vtr2,vtr3,vrel;
49 double mi,mj,meff,damp,ccel,tor1,tor2,tor3;
50 double fn,fs,ft,fs1,fs2,fs3;
51 int *ilist,*jlist,*numneigh,**firstneigh;
52
53 ev_init(eflag,vflag);
54
55 // update rigid body info for owned & ghost atoms if using FixRigid masses
56 // body[i] = which body atom I is in, -1 if none
57 // mass_body = mass of each rigid body
58
59 if (fix_rigid && neighbor->ago == 0) {
60 int tmp;
61 int *body = (int *) fix_rigid->extract("body",tmp);
62 double *mass_body = (double *) fix_rigid->extract("masstotal",tmp);
63 if (atom->nmax > nmax) {
64 memory->destroy(mass_rigid);
65 nmax = atom->nmax;
66 memory->create(mass_rigid,nmax,"pair:mass_rigid");
67 }
68 int nlocal = atom->nlocal;
69 for (i = 0; i < nlocal; i++)
70 if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]];
71 else mass_rigid[i] = 0.0;
72 comm->forward_comm_pair(this);
73 }
74
75 double **x = atom->x;
76 double **v = atom->v;
77 double **f = atom->f;
78 double **omega = atom->omega;
79 double **torque = atom->torque;
80 double *radius = atom->radius;
81 double *rmass = atom->rmass;
82 int *mask = atom->mask;
83 int nlocal = atom->nlocal;
84 int newton_pair = force->newton_pair;
85
86 inum = list->inum;
87 ilist = list->ilist;
88 numneigh = list->numneigh;
89 firstneigh = list->firstneigh;
90
91 // loop over neighbors of my atoms
92
93 for (ii = 0; ii < inum; ii++) {
94 i = ilist[ii];
95 xtmp = x[i][0];
96 ytmp = x[i][1];
97 ztmp = x[i][2];
98 radi = radius[i];
99 jlist = firstneigh[i];
100 jnum = numneigh[i];
101
102 for (jj = 0; jj < jnum; jj++) {
103 j = jlist[jj];
104 j &= NEIGHMASK;
105
106 delx = xtmp - x[j][0];
107 dely = ytmp - x[j][1];
108 delz = ztmp - x[j][2];
109 rsq = delx*delx + dely*dely + delz*delz;
110 radj = radius[j];
111 radsum = radi + radj;
112
113 if (rsq < radsum*radsum) {
114 r = sqrt(rsq);
115 rinv = 1.0/r;
116 rsqinv = 1.0/rsq;
117
118 // relative translational velocity
119
120 vr1 = v[i][0] - v[j][0];
121 vr2 = v[i][1] - v[j][1];
122 vr3 = v[i][2] - v[j][2];
123
124 // normal component
125
126 vnnr = vr1*delx + vr2*dely + vr3*delz;
127 vn1 = delx*vnnr * rsqinv;
128 vn2 = dely*vnnr * rsqinv;
129 vn3 = delz*vnnr * rsqinv;
130
131 // tangential component
132
133 vt1 = vr1 - vn1;
134 vt2 = vr2 - vn2;
135 vt3 = vr3 - vn3;
136
137 // relative rotational velocity
138
139 wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;
140 wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
141 wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
142
143 // meff = effective mass of pair of particles
144 // if I or J part of rigid body, use body mass
145 // if I or J is frozen, meff is other particle
146
147 mi = rmass[i];
148 mj = rmass[j];
149 if (fix_rigid) {
150 if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
151 if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
152 }
153
154 meff = mi*mj / (mi+mj);
155 if (mask[i] & freeze_group_bit) meff = mj;
156 if (mask[j] & freeze_group_bit) meff = mi;
157
158 // normal forces = Hookian contact + normal velocity damping
159
160 damp = meff*gamman*vnnr*rsqinv;
161 ccel = kn*(radsum-r)*rinv - damp;
162 if (limit_damping && (ccel < 0.0)) ccel = 0.0;
163
164 // relative velocities
165
166 vtr1 = vt1 - (delz*wr2-dely*wr3);
167 vtr2 = vt2 - (delx*wr3-delz*wr1);
168 vtr3 = vt3 - (dely*wr1-delx*wr2);
169 vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
170 vrel = sqrt(vrel);
171
172 // force normalization
173
174 fn = xmu * fabs(ccel*r);
175 fs = meff*gammat*vrel;
176 if (vrel != 0.0) ft = MIN(fn,fs) / vrel;
177 else ft = 0.0;
178
179 // tangential force due to tangential velocity damping
180
181 fs1 = -ft*vtr1;
182 fs2 = -ft*vtr2;
183 fs3 = -ft*vtr3;
184
185 // forces & torques
186
187 fx = delx*ccel + fs1;
188 fy = dely*ccel + fs2;
189 fz = delz*ccel + fs3;
190 f[i][0] += fx;
191 f[i][1] += fy;
192 f[i][2] += fz;
193
194 tor1 = rinv * (dely*fs3 - delz*fs2);
195 tor2 = rinv * (delz*fs1 - delx*fs3);
196 tor3 = rinv * (delx*fs2 - dely*fs1);
197 torque[i][0] -= radi*tor1;
198 torque[i][1] -= radi*tor2;
199 torque[i][2] -= radi*tor3;
200
201 if (newton_pair || j < nlocal) {
202 f[j][0] -= fx;
203 f[j][1] -= fy;
204 f[j][2] -= fz;
205 torque[j][0] -= radj*tor1;
206 torque[j][1] -= radj*tor2;
207 torque[j][2] -= radj*tor3;
208 }
209
210 if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
211 0.0,0.0,fx,fy,fz,delx,dely,delz);
212 }
213 }
214 }
215
216 if (vflag_fdotr) virial_fdotr_compute();
217 }
218
219 /* ---------------------------------------------------------------------- */
220
single(int i,int j,int,int,double rsq,double,double,double & fforce)221 double PairGranHooke::single(int i, int j, int /*itype*/, int /*jtype*/, double rsq,
222 double /*factor_coul*/, double /*factor_lj*/,
223 double &fforce)
224 {
225 double radi,radj,radsum,r,rinv,rsqinv;
226 double delx,dely,delz;
227 double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3;
228 double vtr1,vtr2,vtr3,vrel;
229 double mi,mj,meff,damp,ccel;
230 double fn,fs,ft;
231
232 double *radius = atom->radius;
233 radi = radius[i];
234 radj = radius[j];
235 radsum = radi + radj;
236
237 // zero out forces if caller requests non-touching pair outside cutoff
238
239 if (rsq >= radsum*radsum) {
240 fforce = 0.0;
241 for (int m = 0; m < single_extra; m++) svector[m] = 0.0;
242 return 0.0;
243 }
244
245 r = sqrt(rsq);
246 rinv = 1.0/r;
247 rsqinv = 1.0/rsq;
248
249 // relative translational velocity
250
251 double **v = atom->v;
252 vr1 = v[i][0] - v[j][0];
253 vr2 = v[i][1] - v[j][1];
254 vr3 = v[i][2] - v[j][2];
255
256 // normal component
257
258 double **x = atom->x;
259 delx = x[i][0] - x[j][0];
260 dely = x[i][1] - x[j][1];
261 delz = x[i][2] - x[j][2];
262
263 vnnr = vr1*delx + vr2*dely + vr3*delz;
264 vn1 = delx*vnnr * rsqinv;
265 vn2 = dely*vnnr * rsqinv;
266 vn3 = delz*vnnr * rsqinv;
267
268 // tangential component
269
270 vt1 = vr1 - vn1;
271 vt2 = vr2 - vn2;
272 vt3 = vr3 - vn3;
273
274 // relative rotational velocity
275
276 double **omega = atom->omega;
277 wr1 = (radi*omega[i][0] + radj*omega[j][0]) * rinv;
278 wr2 = (radi*omega[i][1] + radj*omega[j][1]) * rinv;
279 wr3 = (radi*omega[i][2] + radj*omega[j][2]) * rinv;
280
281 // meff = effective mass of pair of particles
282 // if I or J part of rigid body, use body mass
283 // if I or J is frozen, meff is other particle
284
285 double *rmass = atom->rmass;
286 int *mask = atom->mask;
287
288 mi = rmass[i];
289 mj = rmass[j];
290 if (fix_rigid) {
291 // NOTE: insure mass_rigid is current for owned+ghost atoms?
292 if (mass_rigid[i] > 0.0) mi = mass_rigid[i];
293 if (mass_rigid[j] > 0.0) mj = mass_rigid[j];
294 }
295
296 meff = mi*mj / (mi+mj);
297 if (mask[i] & freeze_group_bit) meff = mj;
298 if (mask[j] & freeze_group_bit) meff = mi;
299
300 // normal forces = Hookian contact + normal velocity damping
301
302 damp = meff*gamman*vnnr*rsqinv;
303 ccel = kn*(radsum-r)*rinv - damp;
304 if (limit_damping && (ccel < 0.0)) ccel = 0.0;
305
306 // relative velocities
307
308 vtr1 = vt1 - (delz*wr2-dely*wr3);
309 vtr2 = vt2 - (delx*wr3-delz*wr1);
310 vtr3 = vt3 - (dely*wr1-delx*wr2);
311 vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3;
312 vrel = sqrt(vrel);
313
314 // force normalization
315
316 fn = xmu * fabs(ccel*r);
317 fs = meff*gammat*vrel;
318 if (vrel != 0.0) ft = MIN(fn,fs) / vrel;
319 else ft = 0.0;
320
321 // set force and return no energy
322
323 fforce = ccel;
324
325
326 // set single_extra quantities
327
328 svector[0] = -ft*vtr1;
329 svector[1] = -ft*vtr2;
330 svector[2] = -ft*vtr3;
331 svector[3] = sqrt(svector[0]*svector[0] +
332 svector[1]*svector[1] +
333 svector[2]*svector[2]);
334 svector[4] = vn1;
335 svector[5] = vn2;
336 svector[6] = vn3;
337 svector[7] = vt1;
338 svector[8] = vt2;
339 svector[9] = vt3;
340
341 return 0.0;
342 }
343