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