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 author: Tony Sheh (U Michigan), Trung Dac Nguyen (U Michigan)
17 references: Kamberaj et al., J. Chem. Phys. 122, 224114 (2005)
18 Miller et al., J Chem Phys. 116, 8649-8659 (2002)
19 ------------------------------------------------------------------------- */
20
21 /* ----------------------------------------------------------------------
22 Contributing author:
23 Morteza Jalalvand (IASBS) jalalvand.m AT gmail.com
24
25 This is an extension of fix/rigid/nve to SPH/SDPD particles
26 You can see the original copyright notice of fix/rigid authors above
27 Note that the Kamberaj paper was related to the nvt variant
28 and all codes relevant to that has been removed
29 ------------------------------------------------------------------------- */
30
31 #include "fix_rigid_meso.h"
32 #include "math_extra.h"
33 #include "atom.h"
34 #include "domain.h"
35 #include "memory.h"
36 #include "error.h"
37
38 using namespace LAMMPS_NS;
39 using namespace FixConst;
40
41 /* ---------------------------------------------------------------------- */
42
FixRigidMeso(LAMMPS * lmp,int narg,char ** arg)43 FixRigidMeso::FixRigidMeso (LAMMPS *lmp, int narg, char **arg) :
44 FixRigid (lmp, narg, arg)
45 {
46 scalar_flag = 0;
47 size_array_cols = 28;
48
49 if ((atom->esph_flag != 1) || (atom->rho_flag != 1))
50 error->all (FLERR, "fix rigid/meso command requires atom_style with"
51 " both energy and density");
52
53 if (langflag || tstat_flag)
54 error->all (FLERR,"Can not use thermostat with fix rigid/meso");
55
56 if (pstat_flag)
57 error->all (FLERR,"Can not use barostat with fix rigid/meso");
58
59 // memory allocation and initialization
60
61 memory->create(conjqm,nbody,4,"rigid_nh:conjqm");
62 }
63
64 /* ---------------------------------------------------------------------- */
65
~FixRigidMeso()66 FixRigidMeso::~FixRigidMeso () {
67 memory->destroy(conjqm);
68 }
69
70 /* ---------------------------------------------------------------------- */
71
setmask()72 int FixRigidMeso::setmask () {
73 int mask = 0;
74 mask |= INITIAL_INTEGRATE;
75 mask |= FINAL_INTEGRATE;
76 mask |= PRE_NEIGHBOR;
77 return mask;
78 }
79
80 /* ---------------------------------------------------------------------- */
81
setup(int vflag)82 void FixRigidMeso::setup (int vflag) {
83 FixRigid::setup(vflag);
84
85 double mbody[3];
86 for (int ibody = 0; ibody < nbody; ibody++) {
87 MathExtra::transpose_matvec (ex_space[ibody],ey_space[ibody],ez_space[ibody],
88 angmom[ibody],mbody);
89 MathExtra::quatvec (quat[ibody],mbody,conjqm[ibody]);
90 conjqm[ibody][0] *= 2.0;
91 conjqm[ibody][1] *= 2.0;
92 conjqm[ibody][2] *= 2.0;
93 conjqm[ibody][3] *= 2.0;
94 }
95 }
96
97 /* ----------------------------------------------------------------------
98 perform preforce velocity Verlet integration
99 see Kamberaj paper for step references
100 ------------------------------------------------------------------------- */
101
initial_integrate(int vflag)102 void FixRigidMeso::initial_integrate (int vflag) {
103 double dtfm,mbody[3],tbody[3],fquat[4];
104 double dtf2 = dtf * 2.0;
105
106 // update xcm, vcm, quat, conjqm and angmom
107
108 for (int ibody = 0; ibody < nbody; ibody++) {
109
110 // step 1.1 - update vcm by 1/2 step
111
112 dtfm = dtf / masstotal[ibody];
113 vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
114 vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
115 vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
116
117 // step 1.2 - update xcm by full step
118
119 xcm[ibody][0] += dtv * vcm[ibody][0];
120 xcm[ibody][1] += dtv * vcm[ibody][1];
121 xcm[ibody][2] += dtv * vcm[ibody][2];
122
123 // step 1.3 - apply torque (body coords) to quaternion momentum
124
125 torque[ibody][0] *= tflag[ibody][0];
126 torque[ibody][1] *= tflag[ibody][1];
127 torque[ibody][2] *= tflag[ibody][2];
128
129 MathExtra::transpose_matvec (ex_space[ibody],ey_space[ibody],ez_space[ibody],
130 torque[ibody],tbody);
131 MathExtra::quatvec (quat[ibody],tbody,fquat);
132
133 conjqm[ibody][0] += dtf2 * fquat[0];
134 conjqm[ibody][1] += dtf2 * fquat[1];
135 conjqm[ibody][2] += dtf2 * fquat[2];
136 conjqm[ibody][3] += dtf2 * fquat[3];
137
138 // step 1.4 to 1.13 - use no_squish rotate to update p and q
139
140 MathExtra::no_squish_rotate (3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
141 MathExtra::no_squish_rotate (2,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
142 MathExtra::no_squish_rotate (1,conjqm[ibody],quat[ibody],inertia[ibody],dtv);
143 MathExtra::no_squish_rotate (2,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
144 MathExtra::no_squish_rotate (3,conjqm[ibody],quat[ibody],inertia[ibody],dtq);
145
146 // update exyz_space
147 // transform p back to angmom
148 // update angular velocity
149
150 MathExtra::q_to_exyz (quat[ibody],ex_space[ibody],ey_space[ibody],
151 ez_space[ibody]);
152 MathExtra::invquatvec (quat[ibody],conjqm[ibody],mbody);
153 MathExtra::matvec (ex_space[ibody],ey_space[ibody],ez_space[ibody],
154 mbody,angmom[ibody]);
155
156 angmom[ibody][0] *= 0.5;
157 angmom[ibody][1] *= 0.5;
158 angmom[ibody][2] *= 0.5;
159
160 MathExtra::angmom_to_omega (angmom[ibody],ex_space[ibody],ey_space[ibody],
161 ez_space[ibody],inertia[ibody],omega[ibody]);
162 }
163
164 // virial setup before call to set_xv
165
166 v_init(vflag);
167
168 // set coords/orient and velocity/rotation of atoms in rigid bodies
169 // from quarternion and omega
170
171 set_xv();
172 }
173
174 /* ---------------------------------------------------------------------- */
175
final_integrate()176 void FixRigidMeso::final_integrate () {
177 int ibody;
178 double dtfm;
179 double mbody[3],tbody[3],fquat[4];
180
181 double dtf2 = dtf * 2.0;
182
183 // late calculation of forces and torques (if requested)
184
185 if (!earlyflag) compute_forces_and_torques();
186
187 // update vcm and angmom
188 // fflag,tflag = 0 for some dimensions in 2d
189
190 for (ibody = 0; ibody < nbody; ibody++) {
191
192 // update vcm by 1/2 step
193
194 dtfm = dtf / masstotal[ibody];
195
196 vcm[ibody][0] += dtfm * fcm[ibody][0] * fflag[ibody][0];
197 vcm[ibody][1] += dtfm * fcm[ibody][1] * fflag[ibody][1];
198 vcm[ibody][2] += dtfm * fcm[ibody][2] * fflag[ibody][2];
199
200 // update conjqm, then transform to angmom, set velocity again
201 // virial is already setup from initial_integrate
202
203 torque[ibody][0] *= tflag[ibody][0];
204 torque[ibody][1] *= tflag[ibody][1];
205 torque[ibody][2] *= tflag[ibody][2];
206
207 MathExtra::transpose_matvec (ex_space[ibody],ey_space[ibody],
208 ez_space[ibody],torque[ibody],tbody);
209 MathExtra::quatvec (quat[ibody],tbody,fquat);
210
211 conjqm[ibody][0] += dtf2 * fquat[0];
212 conjqm[ibody][1] += dtf2 * fquat[1];
213 conjqm[ibody][2] += dtf2 * fquat[2];
214 conjqm[ibody][3] += dtf2 * fquat[3];
215
216 MathExtra::invquatvec (quat[ibody],conjqm[ibody],mbody);
217 MathExtra::matvec (ex_space[ibody],ey_space[ibody],ez_space[ibody],
218 mbody,angmom[ibody]);
219
220 angmom[ibody][0] *= 0.5;
221 angmom[ibody][1] *= 0.5;
222 angmom[ibody][2] *= 0.5;
223
224 MathExtra::angmom_to_omega (angmom[ibody],ex_space[ibody],ey_space[ibody],
225 ez_space[ibody],inertia[ibody],omega[ibody]);
226 }
227
228 // set velocity/rotation of atoms in rigid bodies
229 // virial is already setup from initial_integrate
230
231 set_v();
232 }
233
234 /* ----------------------------------------------------------------------
235 set space-frame coords and velocity of each atom in each rigid body
236 set orientation and rotation of extended particles
237 x = Q displace + Xcm, mapped back to periodic box
238 v = Vcm + (W cross (x - Xcm))
239 ------------------------------------------------------------------------- */
240
set_xv()241 void FixRigidMeso::set_xv () {
242 int ibody;
243 int xbox,ybox,zbox;
244 double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
245 double xy,xz,yz;
246 double vr[6];
247
248 double **x = atom->x;
249 double **v = atom->v;
250 double **vest = atom->vest;
251 double **f = atom->f;
252 double *esph = atom->esph;
253 double *desph = atom->desph;
254 double *rho = atom->rho;
255 double *drho = atom->drho;
256 double *rmass = atom->rmass;
257 double *mass = atom->mass;
258 int *type = atom->type;
259 int nlocal = atom->nlocal;
260
261 double xprd = domain->xprd;
262 double yprd = domain->yprd;
263 double zprd = domain->zprd;
264
265 if (triclinic) {
266 xy = domain->xy;
267 xz = domain->xz;
268 yz = domain->yz;
269 }
270
271 // set x and v of each atom
272
273 for (int i = 0; i < nlocal; i++) {
274 if (body[i] < 0) continue;
275
276 // half-step update of particle internal energy and density
277 esph[i] += dtf * desph[i];
278 rho[i] += dtf * drho[i];
279
280 ibody = body[i];
281
282 xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
283 ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
284 zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
285
286 // save old positions and velocities for virial
287
288 if (evflag) {
289 if (triclinic == 0) {
290 x0 = x[i][0] + xbox*xprd;
291 x1 = x[i][1] + ybox*yprd;
292 x2 = x[i][2] + zbox*zprd;
293 } else {
294 x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
295 x1 = x[i][1] + ybox*yprd + zbox*yz;
296 x2 = x[i][2] + zbox*zprd;
297 }
298 }
299
300 v0 = v[i][0];
301 v1 = v[i][1];
302 v2 = v[i][2];
303
304 // x = displacement from center-of-mass, based on body orientation
305 // v = vcm + omega around center-of-mass
306 // vest = 2*v - v_old
307
308 MathExtra::matvec (ex_space[ibody],ey_space[ibody],
309 ez_space[ibody],displace[i],x[i]);
310
311 v[i][0] = omega[ibody][1]*x[i][2] - omega[ibody][2]*x[i][1] +
312 vcm[ibody][0];
313 v[i][1] = omega[ibody][2]*x[i][0] - omega[ibody][0]*x[i][2] +
314 vcm[ibody][1];
315 v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] +
316 vcm[ibody][2];
317
318 vest[i][0] = 2*v[i][0] - v0;
319 vest[i][1] = 2*v[i][1] - v1;
320 vest[i][2] = 2*v[i][2] - v2;
321
322 // add center of mass to displacement
323 // map back into periodic box via xbox,ybox,zbox
324 // for triclinic, add in box tilt factors as well
325
326 if (triclinic == 0) {
327 x[i][0] += xcm[ibody][0] - xbox*xprd;
328 x[i][1] += xcm[ibody][1] - ybox*yprd;
329 x[i][2] += xcm[ibody][2] - zbox*zprd;
330 } else {
331 x[i][0] += xcm[ibody][0] - xbox*xprd - ybox*xy - zbox*xz;
332 x[i][1] += xcm[ibody][1] - ybox*yprd - zbox*yz;
333 x[i][2] += xcm[ibody][2] - zbox*zprd;
334 }
335
336 // virial = unwrapped coords dotted into body constraint force
337 // body constraint force = implied force due to v change minus f external
338 // assume f does not include forces internal to body
339 // 1/2 factor b/c final_integrate contributes other half
340 // assume per-atom contribution is due to constraint force on that atom
341
342 if (evflag) {
343 if (rmass) massone = rmass[i];
344 else massone = mass[type[i]];
345 fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
346 fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
347 fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
348
349 vr[0] = 0.5*x0*fc0;
350 vr[1] = 0.5*x1*fc1;
351 vr[2] = 0.5*x2*fc2;
352 vr[3] = 0.5*x0*fc1;
353 vr[4] = 0.5*x0*fc2;
354 vr[5] = 0.5*x1*fc2;
355
356 v_tally(1,&i,1.0,vr);
357 }
358 }
359
360 // set orientation, omega, angmom of each extended particle
361
362 if (extended) {
363 // TBD
364 }
365 }
366
367 /* ----------------------------------------------------------------------
368 set space-frame velocity of each atom in a rigid body
369 set omega and angmom of extended particles
370 v = Vcm + (W cross (x - Xcm))
371 ------------------------------------------------------------------------- */
372
set_v()373 void FixRigidMeso::set_v () {
374 int xbox,ybox,zbox;
375 double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone;
376 double xy,xz,yz;
377 double delta[3],vr[6];
378
379 double **x = atom->x;
380 double **v = atom->v;
381 double **f = atom->f;
382 double *esph = atom->esph;
383 double *desph = atom->desph;
384 double *rho = atom->rho;
385 double *drho = atom->drho;
386 double *rmass = atom->rmass;
387 double *mass = atom->mass;
388 int *type = atom->type;
389 int nlocal = atom->nlocal;
390
391 double xprd = domain->xprd;
392 double yprd = domain->yprd;
393 double zprd = domain->zprd;
394 if (triclinic) {
395 xy = domain->xy;
396 xz = domain->xz;
397 yz = domain->yz;
398 }
399
400 // set v of each atom
401
402 for (int i = 0; i < nlocal; i++) {
403 if (body[i] < 0) continue;
404
405 // half-step update of particle internal energy and density
406 esph[i] += dtf * desph[i];
407 rho[i] += dtf * drho[i];
408
409 const int ibody = body[i];
410
411 MathExtra::matvec (ex_space[ibody],ey_space[ibody],
412 ez_space[ibody],displace[i],delta);
413
414 // save old velocities for virial
415
416 if (evflag) {
417 v0 = v[i][0];
418 v1 = v[i][1];
419 v2 = v[i][2];
420 }
421
422 v[i][0] = omega[ibody][1]*delta[2] - omega[ibody][2]*delta[1] +
423 vcm[ibody][0];
424 v[i][1] = omega[ibody][2]*delta[0] - omega[ibody][0]*delta[2] +
425 vcm[ibody][1];
426 v[i][2] = omega[ibody][0]*delta[1] - omega[ibody][1]*delta[0] +
427 vcm[ibody][2];
428
429 // virial = unwrapped coords dotted into body constraint force
430 // body constraint force = implied force due to v change minus f external
431 // assume f does not include forces internal to body
432 // 1/2 factor b/c initial_integrate contributes other half
433 // assume per-atom contribution is due to constraint force on that atom
434
435 if (evflag) {
436 if (rmass) massone = rmass[i];
437 else massone = mass[type[i]];
438 fc0 = massone*(v[i][0] - v0)/dtf - f[i][0];
439 fc1 = massone*(v[i][1] - v1)/dtf - f[i][1];
440 fc2 = massone*(v[i][2] - v2)/dtf - f[i][2];
441
442 xbox = (xcmimage[i] & IMGMASK) - IMGMAX;
443 ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX;
444 zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX;
445
446 if (triclinic == 0) {
447 x0 = x[i][0] + xbox*xprd;
448 x1 = x[i][1] + ybox*yprd;
449 x2 = x[i][2] + zbox*zprd;
450 } else {
451 x0 = x[i][0] + xbox*xprd + ybox*xy + zbox*xz;
452 x1 = x[i][1] + ybox*yprd + zbox*yz;
453 x2 = x[i][2] + zbox*zprd;
454 }
455
456 vr[0] = 0.5*x0*fc0;
457 vr[1] = 0.5*x1*fc1;
458 vr[2] = 0.5*x2*fc2;
459 vr[3] = 0.5*x0*fc1;
460 vr[4] = 0.5*x0*fc2;
461 vr[5] = 0.5*x1*fc2;
462
463 v_tally(1,&i,1.0,vr);
464 }
465 }
466
467 // set omega, angmom of each extended particle
468
469 if (extended) {
470 // TBD
471 }
472 }
473
474 /* ----------------------------------------------------------------------
475 return attributes of a rigid body
476 19 values per body
477 xcm = 0,1,2; vcm = 3,4,5; fcm = 6,7,8;
478 quat = 9,10,11,12; omega = 13,14,15; torque = 16,17,18;
479 inertia = 19,20,21; angmom = 22,23,24;
480 image = 25,26,27
481 ------------------------------------------------------------------------- */
482
compute_array(int i,int j)483 double FixRigidMeso::compute_array (int i, int j) {
484 if (j < 3) return xcm[i][j];
485 if (j < 6) return vcm[i][j-3];
486 if (j < 9) return fcm[i][j-6];
487 if (j < 13) return quat[i][j-9];
488 if (j < 16) return omega[i][j-13];
489 if (j < 19) return torque[i][j-16];
490 if (j < 22) return inertia[i][j-19];
491 if (j < 25) return angmom[i][j-22];
492 if (j == 25) return (imagebody[i] & IMGMASK) - IMGMAX;
493 if (j == 26) return (imagebody[i] >> IMGBITS & IMGMASK) - IMGMAX;
494 return (imagebody[i] >> IMG2BITS) - IMGMAX;
495 }
496