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 #include "displace_atoms.h"
16 
17 #include "atom.h"
18 #include "atom_vec_body.h"
19 #include "atom_vec_ellipsoid.h"
20 #include "atom_vec_line.h"
21 #include "atom_vec_tri.h"
22 #include "comm.h"
23 #include "domain.h"
24 #include "error.h"
25 #include "group.h"
26 #include "input.h"
27 #include "irregular.h"
28 #include "lattice.h"
29 #include "math_const.h"
30 #include "math_extra.h"
31 #include "memory.h"
32 #include "modify.h"
33 #include "random_park.h"
34 #include "variable.h"
35 
36 #include <cmath>
37 #include <cstring>
38 
39 using namespace LAMMPS_NS;
40 using namespace MathConst;
41 
42 enum{MOVE,RAMP,RANDOM,ROTATE};
43 
44 /* ---------------------------------------------------------------------- */
45 
DisplaceAtoms(LAMMPS * lmp)46 DisplaceAtoms::DisplaceAtoms(LAMMPS *lmp) : Command(lmp)
47 {
48   mvec = nullptr;
49 }
50 
51 /* ---------------------------------------------------------------------- */
52 
~DisplaceAtoms()53 DisplaceAtoms::~DisplaceAtoms()
54 {
55   memory->destroy(mvec);
56 }
57 
58 /* ---------------------------------------------------------------------- */
59 
command(int narg,char ** arg)60 void DisplaceAtoms::command(int narg, char **arg)
61 {
62   int i;
63 
64   if (domain->box_exist == 0)
65     error->all(FLERR,"Displace_atoms command before simulation box is defined");
66   if (narg < 2) error->all(FLERR,"Illegal displace_atoms command");
67   if (modify->nfix_restart_peratom)
68     error->all(FLERR,"Cannot displace_atoms after "
69                "reading restart file with per-atom info");
70 
71   if (comm->me == 0) utils::logmesg(lmp,"Displacing atoms ...\n");
72 
73   // group and style
74 
75   igroup = group->find(arg[0]);
76   if (igroup == -1) error->all(FLERR,"Could not find displace_atoms group ID");
77   groupbit = group->bitmask[igroup];
78 
79   if (modify->check_rigid_group_overlap(groupbit))
80     error->warning(FLERR,"Attempting to displace atoms in rigid bodies");
81 
82   int style = -1;
83   if (strcmp(arg[1],"move") == 0) style = MOVE;
84   else if (strcmp(arg[1],"ramp") == 0) style = RAMP;
85   else if (strcmp(arg[1],"random") == 0) style = RANDOM;
86   else if (strcmp(arg[1],"rotate") == 0) style = ROTATE;
87   else error->all(FLERR,"Illegal displace_atoms command");
88 
89   // set option defaults
90 
91   scaleflag = 1;
92 
93   // read options from end of input line
94 
95   if (style == MOVE) options(narg-5,&arg[5]);
96   else if (style == RAMP) options(narg-8,&arg[8]);
97   else if (style == RANDOM) options(narg-6,&arg[6]);
98   else if (style == ROTATE) options(narg-9,&arg[9]);
99 
100   // setup scaling
101 
102   double xscale,yscale,zscale;
103   if (scaleflag) {
104     xscale = domain->lattice->xlattice;
105     yscale = domain->lattice->ylattice;
106     zscale = domain->lattice->zlattice;
107   }
108   else xscale = yscale = zscale = 1.0;
109 
110   // move atoms by 3-vector or specified variable(s)
111 
112   if (style == MOVE) {
113     move(0,arg[2],xscale);
114     move(1,arg[3],yscale);
115     move(2,arg[4],zscale);
116   }
117 
118   // move atoms in ramped fashion
119 
120   if (style == RAMP) {
121 
122     int d_dim = 0;
123     if (strcmp(arg[2],"x") == 0) d_dim = 0;
124     else if (strcmp(arg[2],"y") == 0) d_dim = 1;
125     else if (strcmp(arg[2],"z") == 0) d_dim = 2;
126     else error->all(FLERR,"Illegal displace_atoms ramp command");
127 
128     double d_lo,d_hi;
129     if (d_dim == 0) {
130       d_lo = xscale*utils::numeric(FLERR,arg[3],false,lmp);
131       d_hi = xscale*utils::numeric(FLERR,arg[4],false,lmp);
132     } else if (d_dim == 1) {
133       d_lo = yscale*utils::numeric(FLERR,arg[3],false,lmp);
134       d_hi = yscale*utils::numeric(FLERR,arg[4],false,lmp);
135     } else if (d_dim == 2) {
136       d_lo = zscale*utils::numeric(FLERR,arg[3],false,lmp);
137       d_hi = zscale*utils::numeric(FLERR,arg[4],false,lmp);
138     }
139 
140     int coord_dim = 0;
141     if (strcmp(arg[5],"x") == 0) coord_dim = 0;
142     else if (strcmp(arg[5],"y") == 0) coord_dim = 1;
143     else if (strcmp(arg[5],"z") == 0) coord_dim = 2;
144     else error->all(FLERR,"Illegal displace_atoms ramp command");
145 
146     double coord_lo,coord_hi;
147     if (coord_dim == 0) {
148       coord_lo = xscale*utils::numeric(FLERR,arg[6],false,lmp);
149       coord_hi = xscale*utils::numeric(FLERR,arg[7],false,lmp);
150     } else if (coord_dim == 1) {
151       coord_lo = yscale*utils::numeric(FLERR,arg[6],false,lmp);
152       coord_hi = yscale*utils::numeric(FLERR,arg[7],false,lmp);
153     } else if (coord_dim == 2) {
154       coord_lo = zscale*utils::numeric(FLERR,arg[6],false,lmp);
155       coord_hi = zscale*utils::numeric(FLERR,arg[7],false,lmp);
156     }
157 
158     double **x = atom->x;
159     int *mask = atom->mask;
160     int nlocal = atom->nlocal;
161 
162     double fraction,dramp;
163 
164     for (i = 0; i < nlocal; i++) {
165       if (mask[i] & groupbit) {
166         fraction = (x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo);
167         fraction = MAX(fraction,0.0);
168         fraction = MIN(fraction,1.0);
169         dramp = d_lo + fraction*(d_hi - d_lo);
170         x[i][d_dim] += dramp;
171       }
172     }
173   }
174 
175   // move atoms randomly
176   // makes atom result independent of what proc owns it via random->reset()
177 
178   if (style == RANDOM) {
179     RanPark *random = new RanPark(lmp,1);
180 
181     double dx = xscale*utils::numeric(FLERR,arg[2],false,lmp);
182     double dy = yscale*utils::numeric(FLERR,arg[3],false,lmp);
183     double dz = zscale*utils::numeric(FLERR,arg[4],false,lmp);
184     int seed = utils::inumeric(FLERR,arg[5],false,lmp);
185     if (seed <= 0) error->all(FLERR,"Illegal displace_atoms random command");
186 
187     double **x = atom->x;
188     int *mask = atom->mask;
189     int nlocal = atom->nlocal;
190 
191     for (i = 0; i < nlocal; i++) {
192       if (mask[i] & groupbit) {
193         random->reset(seed,x[i]);
194         x[i][0] += dx * 2.0*(random->uniform()-0.5);
195         x[i][1] += dy * 2.0*(random->uniform()-0.5);
196         x[i][2] += dz * 2.0*(random->uniform()-0.5);
197       }
198     }
199 
200     delete random;
201   }
202 
203   // rotate atoms by right-hand rule by theta around R
204   // P = point = vector = point of rotation
205   // R = vector = axis of rotation
206   // R0 = runit = unit vector for R
207   // D = X - P = vector from P to X
208   // C = (D dot R0) R0 = projection of atom coord onto R line
209   // A = D - C = vector from R line to X
210   // B = R0 cross A = vector perp to A in plane of rotation
211   // A,B define plane of circular rotation around R line
212   // X = P + C + A cos(theta) + B sin(theta)
213 
214   if (style == ROTATE) {
215     double theta_new;
216     double axis[3],point[3],qrotate[4],qnew[4];
217     double a[3],b[3],c[3],d[3],disp[3],runit[3];
218     double *quat;
219 
220     int dim = domain->dimension;
221     point[0] = xscale*utils::numeric(FLERR,arg[2],false,lmp);
222     point[1] = yscale*utils::numeric(FLERR,arg[3],false,lmp);
223     point[2] = zscale*utils::numeric(FLERR,arg[4],false,lmp);
224     axis[0] = utils::numeric(FLERR,arg[5],false,lmp);
225     axis[1] = utils::numeric(FLERR,arg[6],false,lmp);
226     axis[2] = utils::numeric(FLERR,arg[7],false,lmp);
227     double theta = utils::numeric(FLERR,arg[8],false,lmp);
228     if (dim == 2 && (axis[0] != 0.0 || axis[1] != 0.0))
229       error->all(FLERR,"Invalid displace_atoms rotate axis for 2d");
230 
231     double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
232     if (len == 0.0)
233       error->all(FLERR,"Zero length rotation vector with displace_atoms");
234     runit[0] = axis[0]/len;
235     runit[1] = axis[1]/len;
236     runit[2] = axis[2]/len;
237 
238     double angle = MY_PI*theta/180.0;
239     double cosine = cos(angle);
240     double sine = sin(angle);
241 
242     double qcosine = cos(0.5*angle);
243     double qsine = sin(0.5*angle);
244     qrotate[0] = qcosine;
245     qrotate[1] = runit[0]*qsine;
246     qrotate[2] = runit[1]*qsine;
247     qrotate[3] = runit[2]*qsine;
248 
249     double ddotr;
250 
251     // flags for additional orientation info stored by some atom styles
252 
253     int ellipsoid_flag = atom->ellipsoid_flag;
254     int line_flag = atom->line_flag;
255     int tri_flag = atom->tri_flag;
256     int body_flag = atom->body_flag;
257 
258     int theta_flag = 0;
259     int quat_flag = 0;
260     if (line_flag) theta_flag = 1;
261     if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1;
262 
263     // AtomVec pointers to retrieve per-atom storage of extra quantities
264 
265     AtomVecEllipsoid *avec_ellipsoid =
266       (AtomVecEllipsoid *) atom->style_match("ellipsoid");
267     AtomVecLine *avec_line = (AtomVecLine *) atom->style_match("line");
268     AtomVecTri *avec_tri = (AtomVecTri *) atom->style_match("tri");
269     AtomVecBody *avec_body = (AtomVecBody *) atom->style_match("body");
270 
271     double **x = atom->x;
272     int *ellipsoid = atom->ellipsoid;
273     int *line = atom->line;
274     int *tri = atom->tri;
275     int *body = atom->body;
276     int *mask = atom->mask;
277     int nlocal = atom->nlocal;
278     imageint *image = atom->image;
279 
280     for (i = 0; i < nlocal; i++) {
281       if (mask[i] & groupbit) {
282         // unwrap coordinate and reset image flags accordingly
283         domain->unmap(x[i],image[i]);
284         image[i] = ((imageint) IMGMAX << IMG2BITS) |
285           ((imageint) IMGMAX << IMGBITS) | IMGMAX;
286 
287         d[0] = x[i][0] - point[0];
288         d[1] = x[i][1] - point[1];
289         d[2] = x[i][2] - point[2];
290         ddotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2];
291         c[0] = ddotr*runit[0];
292         c[1] = ddotr*runit[1];
293         c[2] = ddotr*runit[2];
294         a[0] = d[0] - c[0];
295         a[1] = d[1] - c[1];
296         a[2] = d[2] - c[2];
297         b[0] = runit[1]*a[2] - runit[2]*a[1];
298         b[1] = runit[2]*a[0] - runit[0]*a[2];
299         b[2] = runit[0]*a[1] - runit[1]*a[0];
300         disp[0] = a[0]*cosine  + b[0]*sine;
301         disp[1] = a[1]*cosine  + b[1]*sine;
302         disp[2] = a[2]*cosine  + b[2]*sine;
303         x[i][0] = point[0] + c[0] + disp[0];
304         x[i][1] = point[1] + c[1] + disp[1];
305         if (dim == 3) x[i][2] = point[2] + c[2] + disp[2];
306 
307         // theta for lines
308 
309         if (theta_flag && line[i] >= 0.0) {
310           theta_new = fmod(avec_line->bonus[line[i]].theta+angle,MY_2PI);
311           avec_line->bonus[atom->line[i]].theta = theta_new;
312         }
313 
314         // quats for ellipsoids, tris, and bodies
315 
316         if (quat_flag) {
317           quat = nullptr;
318           if (ellipsoid_flag && ellipsoid[i] >= 0)
319             quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
320           else if (tri_flag && tri[i] >= 0)
321             quat = avec_tri->bonus[tri[i]].quat;
322           else if (body_flag && body[i] >= 0)
323             quat = avec_body->bonus[body[i]].quat;
324           if (quat) {
325             MathExtra::quatquat(qrotate,quat,qnew);
326             quat[0] = qnew[0];
327             quat[1] = qnew[1];
328             quat[2] = qnew[2];
329             quat[3] = qnew[3];
330           }
331         }
332       }
333     }
334   }
335 
336   // move atoms back inside simulation box and to new processors
337   // use remap() instead of pbc() in case atoms moved a long distance
338   // use irregular() in case atoms moved a long distance
339 
340   double **x = atom->x;
341   imageint *image = atom->image;
342   int nlocal = atom->nlocal;
343   for (i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
344 
345   if (domain->triclinic) domain->x2lamda(atom->nlocal);
346   domain->reset_box();
347   Irregular *irregular = new Irregular(lmp);
348   irregular->migrate_atoms(1);
349   delete irregular;
350   if (domain->triclinic) domain->lamda2x(atom->nlocal);
351 
352   // check if any atoms were lost
353 
354   bigint natoms;
355   bigint nblocal = atom->nlocal;
356   MPI_Allreduce(&nblocal,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
357   if (natoms != atom->natoms && comm->me == 0)
358     error->warning(FLERR,"Lost atoms via displace_atoms: original {} "
359                    "current {}",atom->natoms,natoms);
360 }
361 
362 /* ----------------------------------------------------------------------
363    move atoms either by specified numeric displacement or variable evaluation
364 ------------------------------------------------------------------------- */
365 
move(int idim,char * arg,double scale)366 void DisplaceAtoms::move(int idim, char *arg, double scale)
367 {
368   double **x = atom->x;
369   int *mask = atom->mask;
370   int nlocal = atom->nlocal;
371 
372   if (strstr(arg,"v_") != arg) {
373     double delta = scale*utils::numeric(FLERR,arg,false,lmp);
374     for (int i = 0; i < nlocal; i++)
375       if (mask[i] & groupbit) x[i][idim] += delta;
376 
377   } else {
378     int ivar = input->variable->find(arg+2);
379     if (ivar < 0)
380       error->all(FLERR,"Variable name for displace_atoms does not exist");
381 
382     if (input->variable->equalstyle(ivar)) {
383       double delta = scale * input->variable->compute_equal(ivar);
384       for (int i = 0; i < nlocal; i++)
385         if (mask[i] & groupbit) x[i][idim] += delta;
386     } else if (input->variable->atomstyle(ivar)) {
387       if (mvec == nullptr) memory->create(mvec,nlocal,"displace_atoms:mvec");
388       input->variable->compute_atom(ivar,igroup,mvec,1,0);
389       for (int i = 0; i < nlocal; i++)
390         if (mask[i] & groupbit) x[i][idim] += scale*mvec[i];
391     } else error->all(FLERR,"Variable for displace_atoms is invalid style");
392   }
393 }
394 
395 /* ----------------------------------------------------------------------
396    parse optional parameters at end of displace_atoms input line
397 ------------------------------------------------------------------------- */
398 
options(int narg,char ** arg)399 void DisplaceAtoms::options(int narg, char **arg)
400 {
401   if (narg < 0) error->all(FLERR,"Illegal displace_atoms command");
402 
403   int iarg = 0;
404   while (iarg < narg) {
405     if (strcmp(arg[iarg],"units") == 0) {
406       if (iarg+2 > narg) error->all(FLERR,"Illegal displace_atoms command");
407       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
408       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
409       else error->all(FLERR,"Illegal displace_atoms command");
410       iarg += 2;
411     } else error->all(FLERR,"Illegal displace_atoms command");
412   }
413 }
414