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 (triclinic) : Pieter in 't Veld (SNL)
17 ------------------------------------------------------------------------- */
18
19 #include "domain.h"
20 #include "style_region.h" // IWYU pragma: keep
21
22 #include "atom.h"
23 #include "atom_vec.h"
24 #include "comm.h"
25 #include "error.h"
26 #include "fix.h"
27 #include "fix_deform.h"
28 #include "force.h"
29 #include "kspace.h"
30 #include "lattice.h"
31 #include "memory.h"
32 #include "modify.h"
33 #include "molecule.h"
34 #include "output.h"
35 #include "region.h"
36 #include "thermo.h"
37 #include "universe.h"
38 #include "update.h"
39
40 #include <cstring>
41 #include <cmath>
42
43 using namespace LAMMPS_NS;
44
45 #define BIG 1.0e20
46 #define SMALL 1.0e-4
47 #define DELTAREGION 4
48 #define BONDSTRETCH 1.1
49
50 /* ----------------------------------------------------------------------
51 default is periodic
52 ------------------------------------------------------------------------- */
53
Domain(LAMMPS * lmp)54 Domain::Domain(LAMMPS *lmp) : Pointers(lmp)
55 {
56 box_exist = 0;
57 box_change = 0;
58 deform_flag = deform_vremap = deform_groupbit = 0;
59
60 dimension = 3;
61 nonperiodic = 0;
62 xperiodic = yperiodic = zperiodic = 1;
63 periodicity[0] = xperiodic;
64 periodicity[1] = yperiodic;
65 periodicity[2] = zperiodic;
66
67 boundary[0][0] = boundary[0][1] = 0;
68 boundary[1][0] = boundary[1][1] = 0;
69 boundary[2][0] = boundary[2][1] = 0;
70
71 minxlo = minxhi = 0.0;
72 minylo = minyhi = 0.0;
73 minzlo = minzhi = 0.0;
74
75 triclinic = 0;
76 tiltsmall = 1;
77
78 boxlo[0] = boxlo[1] = boxlo[2] = -0.5;
79 boxhi[0] = boxhi[1] = boxhi[2] = 0.5;
80 xy = xz = yz = 0.0;
81
82 h[3] = h[4] = h[5] = 0.0;
83 h_inv[3] = h_inv[4] = h_inv[5] = 0.0;
84 h_rate[0] = h_rate[1] = h_rate[2] =
85 h_rate[3] = h_rate[4] = h_rate[5] = 0.0;
86 h_ratelo[0] = h_ratelo[1] = h_ratelo[2] = 0.0;
87
88 prd_lamda[0] = prd_lamda[1] = prd_lamda[2] = 1.0;
89 prd_half_lamda[0] = prd_half_lamda[1] = prd_half_lamda[2] = 0.5;
90 boxlo_lamda[0] = boxlo_lamda[1] = boxlo_lamda[2] = 0.0;
91 boxhi_lamda[0] = boxhi_lamda[1] = boxhi_lamda[2] = 1.0;
92
93 lattice = nullptr;
94 char **args = new char*[2];
95 args[0] = (char *) "none";
96 args[1] = (char *) "1.0";
97 set_lattice(2,args);
98 delete [] args;
99
100 nregion = maxregion = 0;
101 regions = nullptr;
102
103 copymode = 0;
104
105 region_map = new RegionCreatorMap();
106
107 #define REGION_CLASS
108 #define RegionStyle(key,Class) \
109 (*region_map)[#key] = ®ion_creator<Class>;
110 #include "style_region.h" // IWYU pragma: keep
111
112 #undef RegionStyle
113 #undef REGION_CLASS
114 }
115
116 /* ---------------------------------------------------------------------- */
117
~Domain()118 Domain::~Domain()
119 {
120 if (copymode) return;
121
122 delete lattice;
123 for (int i = 0; i < nregion; i++) delete regions[i];
124 memory->sfree(regions);
125
126 delete region_map;
127 }
128
129 /* ---------------------------------------------------------------------- */
130
init()131 void Domain::init()
132 {
133 // set box_change flags if box size/shape/sub-domains ever change
134 // due to shrink-wrapping or fixes that change box size/shape/sub-domains
135
136 box_change_size = box_change_shape = box_change_domain = 0;
137
138 // flags for detecting, if multiple fixes try to change the
139 // same box size or shape parameter
140
141 int box_change_x=0, box_change_y=0, box_change_z=0;
142 int box_change_yz=0, box_change_xz=0, box_change_xy=0;
143 Fix **fixes = modify->fix;
144
145 if (nonperiodic == 2) box_change_size = 1;
146 for (int i = 0; i < modify->nfix; i++) {
147 if (fixes[i]->box_change & Fix::BOX_CHANGE_SIZE) box_change_size = 1;
148 if (fixes[i]->box_change & Fix::BOX_CHANGE_SHAPE) box_change_shape = 1;
149 if (fixes[i]->box_change & Fix::BOX_CHANGE_DOMAIN) box_change_domain = 1;
150 if (fixes[i]->box_change & Fix::BOX_CHANGE_X) box_change_x++;
151 if (fixes[i]->box_change & Fix::BOX_CHANGE_Y) box_change_y++;
152 if (fixes[i]->box_change & Fix::BOX_CHANGE_Z) box_change_z++;
153 if (fixes[i]->box_change & Fix::BOX_CHANGE_YZ) box_change_yz++;
154 if (fixes[i]->box_change & Fix::BOX_CHANGE_XZ) box_change_xz++;
155 if (fixes[i]->box_change & Fix::BOX_CHANGE_XY) box_change_xy++;
156 }
157
158 std::string mesg = "Must not have multiple fixes change box parameter ";
159
160 #define CHECK_BOX_FIX_ERROR(par) \
161 if (box_change_ ## par > 1) error->all(FLERR,(mesg + #par))
162
163 CHECK_BOX_FIX_ERROR(x);
164 CHECK_BOX_FIX_ERROR(y);
165 CHECK_BOX_FIX_ERROR(z);
166 CHECK_BOX_FIX_ERROR(yz);
167 CHECK_BOX_FIX_ERROR(xz);
168 CHECK_BOX_FIX_ERROR(xy);
169 #undef CHECK_BOX_FIX_ERROR
170
171 box_change = 0;
172 if (box_change_size || box_change_shape || box_change_domain) box_change = 1;
173
174 // check for fix deform
175
176 deform_flag = deform_vremap = deform_groupbit = 0;
177 for (int i = 0; i < modify->nfix; i++)
178 if (utils::strmatch(modify->fix[i]->style,"^deform")) {
179 deform_flag = 1;
180 if (((FixDeform *) modify->fix[i])->remapflag == Domain::V_REMAP) {
181 deform_vremap = 1;
182 deform_groupbit = modify->fix[i]->groupbit;
183 }
184 }
185
186 // region inits
187
188 for (int i = 0; i < nregion; i++) regions[i]->init();
189 }
190
191 /* ----------------------------------------------------------------------
192 set initial global box
193 assumes boxlo/hi and triclinic tilts are already set
194 expandflag = 1 if need to expand box in shrink-wrapped dims
195 not invoked by read_restart since box is already expanded
196 if don't prevent further expansion, restarted triclinic box
197 with unchanged tilt factors can become a box with atoms outside the box
198 ------------------------------------------------------------------------- */
199
set_initial_box(int expandflag)200 void Domain::set_initial_box(int expandflag)
201 {
202 // error checks for orthogonal and triclinic domains
203
204 if (boxlo[0] >= boxhi[0] || boxlo[1] >= boxhi[1] || boxlo[2] >= boxhi[2])
205 error->one(FLERR,"Box bounds are invalid or missing");
206
207 if (domain->dimension == 2 && (xz != 0.0 || yz != 0.0))
208 error->all(FLERR,"Cannot skew triclinic box in z for 2d simulation");
209
210 // error check or warning on triclinic tilt factors
211
212 if (triclinic) {
213 if ((fabs(xy/(boxhi[0]-boxlo[0])) > 0.5 && xperiodic) ||
214 (fabs(xz/(boxhi[0]-boxlo[0])) > 0.5 && xperiodic) ||
215 (fabs(yz/(boxhi[1]-boxlo[1])) > 0.5 && yperiodic)) {
216 if (tiltsmall)
217 error->all(FLERR,"Triclinic box skew is too large");
218 else if (comm->me == 0)
219 error->warning(FLERR,"Triclinic box skew is large");
220 }
221 }
222
223 // set small based on box size and SMALL
224 // this works for any unit system
225
226 small[0] = SMALL * (boxhi[0] - boxlo[0]);
227 small[1] = SMALL * (boxhi[1] - boxlo[1]);
228 small[2] = SMALL * (boxhi[2] - boxlo[2]);
229
230 // if expandflag, adjust box lo/hi for shrink-wrapped dims
231
232 if (!expandflag) return;
233
234 if (boundary[0][0] == 2) boxlo[0] -= small[0];
235 else if (boundary[0][0] == 3) minxlo = boxlo[0];
236 if (boundary[0][1] == 2) boxhi[0] += small[0];
237 else if (boundary[0][1] == 3) minxhi = boxhi[0];
238
239 if (boundary[1][0] == 2) boxlo[1] -= small[1];
240 else if (boundary[1][0] == 3) minylo = boxlo[1];
241 if (boundary[1][1] == 2) boxhi[1] += small[1];
242 else if (boundary[1][1] == 3) minyhi = boxhi[1];
243
244 if (boundary[2][0] == 2) boxlo[2] -= small[2];
245 else if (boundary[2][0] == 3) minzlo = boxlo[2];
246 if (boundary[2][1] == 2) boxhi[2] += small[2];
247 else if (boundary[2][1] == 3) minzhi = boxhi[2];
248 }
249
250 /* ----------------------------------------------------------------------
251 set global box params
252 assumes boxlo/hi and triclinic tilts are already set
253 ------------------------------------------------------------------------- */
254
set_global_box()255 void Domain::set_global_box()
256 {
257 prd[0] = xprd = boxhi[0] - boxlo[0];
258 prd[1] = yprd = boxhi[1] - boxlo[1];
259 prd[2] = zprd = boxhi[2] - boxlo[2];
260
261 h[0] = xprd;
262 h[1] = yprd;
263 h[2] = zprd;
264 h_inv[0] = 1.0/h[0];
265 h_inv[1] = 1.0/h[1];
266 h_inv[2] = 1.0/h[2];
267
268 prd_half[0] = xprd_half = 0.5*xprd;
269 prd_half[1] = yprd_half = 0.5*yprd;
270 prd_half[2] = zprd_half = 0.5*zprd;
271
272 if (triclinic) {
273 h[3] = yz;
274 h[4] = xz;
275 h[5] = xy;
276 h_inv[3] = -h[3] / (h[1]*h[2]);
277 h_inv[4] = (h[3]*h[5] - h[1]*h[4]) / (h[0]*h[1]*h[2]);
278 h_inv[5] = -h[5] / (h[0]*h[1]);
279
280 boxlo_bound[0] = MIN(boxlo[0],boxlo[0]+xy);
281 boxlo_bound[0] = MIN(boxlo_bound[0],boxlo_bound[0]+xz);
282 boxlo_bound[1] = MIN(boxlo[1],boxlo[1]+yz);
283 boxlo_bound[2] = boxlo[2];
284
285 boxhi_bound[0] = MAX(boxhi[0],boxhi[0]+xy);
286 boxhi_bound[0] = MAX(boxhi_bound[0],boxhi_bound[0]+xz);
287 boxhi_bound[1] = MAX(boxhi[1],boxhi[1]+yz);
288 boxhi_bound[2] = boxhi[2];
289 }
290 }
291
292 /* ----------------------------------------------------------------------
293 set lamda box params
294 assumes global box is defined and proc assignment has been made
295 uses comm->xyz_split or comm->mysplit
296 to define subbox boundaries in consistent manner
297 ------------------------------------------------------------------------- */
298
set_lamda_box()299 void Domain::set_lamda_box()
300 {
301 if (comm->layout != Comm::LAYOUT_TILED) {
302 int *myloc = comm->myloc;
303 double *xsplit = comm->xsplit;
304 double *ysplit = comm->ysplit;
305 double *zsplit = comm->zsplit;
306
307 sublo_lamda[0] = xsplit[myloc[0]];
308 subhi_lamda[0] = xsplit[myloc[0]+1];
309 sublo_lamda[1] = ysplit[myloc[1]];
310 subhi_lamda[1] = ysplit[myloc[1]+1];
311 sublo_lamda[2] = zsplit[myloc[2]];
312 subhi_lamda[2] = zsplit[myloc[2]+1];
313
314 } else {
315 double (*mysplit)[2] = comm->mysplit;
316
317 sublo_lamda[0] = mysplit[0][0];
318 subhi_lamda[0] = mysplit[0][1];
319 sublo_lamda[1] = mysplit[1][0];
320 subhi_lamda[1] = mysplit[1][1];
321 sublo_lamda[2] = mysplit[2][0];
322 subhi_lamda[2] = mysplit[2][1];
323 }
324 }
325
326 /* ----------------------------------------------------------------------
327 set local subbox params for orthogonal boxes
328 assumes global box is defined and proc assignment has been made
329 uses comm->xyz_split or comm->mysplit
330 to define subbox boundaries in consistent manner
331 insure subhi[max] = boxhi
332 ------------------------------------------------------------------------- */
333
set_local_box()334 void Domain::set_local_box()
335 {
336 if (triclinic) return;
337
338 if (comm->layout != Comm::LAYOUT_TILED) {
339 int *myloc = comm->myloc;
340 int *procgrid = comm->procgrid;
341 double *xsplit = comm->xsplit;
342 double *ysplit = comm->ysplit;
343 double *zsplit = comm->zsplit;
344
345 sublo[0] = boxlo[0] + xprd*xsplit[myloc[0]];
346 if (myloc[0] < procgrid[0]-1) subhi[0] = boxlo[0] + xprd*xsplit[myloc[0]+1];
347 else subhi[0] = boxhi[0];
348
349 sublo[1] = boxlo[1] + yprd*ysplit[myloc[1]];
350 if (myloc[1] < procgrid[1]-1) subhi[1] = boxlo[1] + yprd*ysplit[myloc[1]+1];
351 else subhi[1] = boxhi[1];
352
353 sublo[2] = boxlo[2] + zprd*zsplit[myloc[2]];
354 if (myloc[2] < procgrid[2]-1) subhi[2] = boxlo[2] + zprd*zsplit[myloc[2]+1];
355 else subhi[2] = boxhi[2];
356
357 } else {
358 double (*mysplit)[2] = comm->mysplit;
359
360 sublo[0] = boxlo[0] + xprd*mysplit[0][0];
361 if (mysplit[0][1] < 1.0) subhi[0] = boxlo[0] + xprd*mysplit[0][1];
362 else subhi[0] = boxhi[0];
363
364 sublo[1] = boxlo[1] + yprd*mysplit[1][0];
365 if (mysplit[1][1] < 1.0) subhi[1] = boxlo[1] + yprd*mysplit[1][1];
366 else subhi[1] = boxhi[1];
367
368 sublo[2] = boxlo[2] + zprd*mysplit[2][0];
369 if (mysplit[2][1] < 1.0) subhi[2] = boxlo[2] + zprd*mysplit[2][1];
370 else subhi[2] = boxhi[2];
371 }
372 }
373
374 /* ----------------------------------------------------------------------
375 reset global & local boxes due to global box boundary changes
376 if shrink-wrapped, determine atom extent and reset boxlo/hi
377 for triclinic, atoms must be in lamda coords (0-1) before reset_box is called
378 ------------------------------------------------------------------------- */
379
reset_box()380 void Domain::reset_box()
381 {
382 // perform shrink-wrapping
383
384 // nothing to do for empty systems
385
386 if (atom->natoms == 0) return;
387
388 // compute extent of atoms on this proc
389 // for triclinic, this is done in lamda space
390
391 if (nonperiodic == 2) {
392 double extent[3][2],all[3][2];
393
394 extent[2][0] = extent[1][0] = extent[0][0] = BIG;
395 extent[2][1] = extent[1][1] = extent[0][1] = -BIG;
396
397 double **x = atom->x;
398 int nlocal = atom->nlocal;
399
400 for (int i = 0; i < nlocal; i++) {
401 extent[0][0] = MIN(extent[0][0],x[i][0]);
402 extent[0][1] = MAX(extent[0][1],x[i][0]);
403 extent[1][0] = MIN(extent[1][0],x[i][1]);
404 extent[1][1] = MAX(extent[1][1],x[i][1]);
405 extent[2][0] = MIN(extent[2][0],x[i][2]);
406 extent[2][1] = MAX(extent[2][1],x[i][2]);
407 }
408
409 // compute extent across all procs
410 // flip sign of MIN to do it in one Allreduce MAX
411
412 extent[0][0] = -extent[0][0];
413 extent[1][0] = -extent[1][0];
414 extent[2][0] = -extent[2][0];
415
416 MPI_Allreduce(extent,all,6,MPI_DOUBLE,MPI_MAX,world);
417
418 // for triclinic, convert back to box coords before changing box
419
420 if (triclinic) lamda2x(atom->nlocal);
421
422 // in shrink-wrapped dims, set box by atom extent
423 // if minimum set, enforce min box size settings
424 // for triclinic, convert lamda extent to box coords, then set box lo/hi
425 // decided NOT to do the next comment - don't want to sneakily change tilt
426 // for triclinic, adjust tilt factors if 2nd dim is shrink-wrapped,
427 // so that displacement in 1st dim stays the same
428
429 if (triclinic == 0) {
430 if (xperiodic == 0) {
431 if (boundary[0][0] == 2) boxlo[0] = -all[0][0] - small[0];
432 else if (boundary[0][0] == 3)
433 boxlo[0] = MIN(-all[0][0]-small[0],minxlo);
434 if (boundary[0][1] == 2) boxhi[0] = all[0][1] + small[0];
435 else if (boundary[0][1] == 3) boxhi[0] = MAX(all[0][1]+small[0],minxhi);
436 if (boxlo[0] > boxhi[0]) error->all(FLERR,"Illegal simulation box");
437 }
438 if (yperiodic == 0) {
439 if (boundary[1][0] == 2) boxlo[1] = -all[1][0] - small[1];
440 else if (boundary[1][0] == 3)
441 boxlo[1] = MIN(-all[1][0]-small[1],minylo);
442 if (boundary[1][1] == 2) boxhi[1] = all[1][1] + small[1];
443 else if (boundary[1][1] == 3) boxhi[1] = MAX(all[1][1]+small[1],minyhi);
444 if (boxlo[1] > boxhi[1]) error->all(FLERR,"Illegal simulation box");
445 }
446 if (zperiodic == 0) {
447 if (boundary[2][0] == 2) boxlo[2] = -all[2][0] - small[2];
448 else if (boundary[2][0] == 3)
449 boxlo[2] = MIN(-all[2][0]-small[2],minzlo);
450 if (boundary[2][1] == 2) boxhi[2] = all[2][1] + small[2];
451 else if (boundary[2][1] == 3) boxhi[2] = MAX(all[2][1]+small[2],minzhi);
452 if (boxlo[2] > boxhi[2]) error->all(FLERR,"Illegal simulation box");
453 }
454
455 } else {
456 double lo[3],hi[3];
457 if (xperiodic == 0) {
458 lo[0] = -all[0][0]; lo[1] = 0.0; lo[2] = 0.0;
459 lamda2x(lo,lo);
460 hi[0] = all[0][1]; hi[1] = 0.0; hi[2] = 0.0;
461 lamda2x(hi,hi);
462 if (boundary[0][0] == 2) boxlo[0] = lo[0] - small[0];
463 else if (boundary[0][0] == 3) boxlo[0] = MIN(lo[0]-small[0],minxlo);
464 if (boundary[0][1] == 2) boxhi[0] = hi[0] + small[0];
465 else if (boundary[0][1] == 3) boxhi[0] = MAX(hi[0]+small[0],minxhi);
466 if (boxlo[0] > boxhi[0]) error->all(FLERR,"Illegal simulation box");
467 }
468 if (yperiodic == 0) {
469 lo[0] = 0.0; lo[1] = -all[1][0]; lo[2] = 0.0;
470 lamda2x(lo,lo);
471 hi[0] = 0.0; hi[1] = all[1][1]; hi[2] = 0.0;
472 lamda2x(hi,hi);
473 if (boundary[1][0] == 2) boxlo[1] = lo[1] - small[1];
474 else if (boundary[1][0] == 3) boxlo[1] = MIN(lo[1]-small[1],minylo);
475 if (boundary[1][1] == 2) boxhi[1] = hi[1] + small[1];
476 else if (boundary[1][1] == 3) boxhi[1] = MAX(hi[1]+small[1],minyhi);
477 if (boxlo[1] > boxhi[1]) error->all(FLERR,"Illegal simulation box");
478 //xy *= (boxhi[1]-boxlo[1]) / yprd;
479 }
480 if (zperiodic == 0) {
481 lo[0] = 0.0; lo[1] = 0.0; lo[2] = -all[2][0];
482 lamda2x(lo,lo);
483 hi[0] = 0.0; hi[1] = 0.0; hi[2] = all[2][1];
484 lamda2x(hi,hi);
485 if (boundary[2][0] == 2) boxlo[2] = lo[2] - small[2];
486 else if (boundary[2][0] == 3) boxlo[2] = MIN(lo[2]-small[2],minzlo);
487 if (boundary[2][1] == 2) boxhi[2] = hi[2] + small[2];
488 else if (boundary[2][1] == 3) boxhi[2] = MAX(hi[2]+small[2],minzhi);
489 if (boxlo[2] > boxhi[2]) error->all(FLERR,"Illegal simulation box");
490 //xz *= (boxhi[2]-boxlo[2]) / xprd;
491 //yz *= (boxhi[2]-boxlo[2]) / yprd;
492 }
493 }
494 }
495
496 // reset box whether shrink-wrapping or not
497
498 set_global_box();
499 set_local_box();
500
501 // if shrink-wrapped & kspace is defined (i.e. using MSM), call setup()
502 // also call init() (to test for compatibility) ?
503
504 if (nonperiodic == 2 && force->kspace) {
505 //force->kspace->init();
506 force->kspace->setup();
507 }
508
509 // if shrink-wrapped & triclinic, re-convert to lamda coords for new box
510 // re-invoke pbc() b/c x2lamda result can be outside [0,1] due to roundoff
511
512 if (nonperiodic == 2 && triclinic) {
513 x2lamda(atom->nlocal);
514 pbc();
515 }
516 }
517
518 /* ----------------------------------------------------------------------
519 enforce PBC and modify box image flags for each atom
520 called every reneighboring and by other commands that change atoms
521 resulting coord must satisfy lo <= coord < hi
522 MAX is important since coord - prd < lo can happen when coord = hi
523 if fix deform, remap velocity of fix group atoms by box edge velocities
524 for triclinic, atoms must be in lamda coords (0-1) before pbc is called
525 image = 10 or 20 bits for each dimension depending on sizeof(imageint)
526 increment/decrement in wrap-around fashion
527 ------------------------------------------------------------------------- */
528
pbc()529 void Domain::pbc()
530 {
531 int nlocal = atom->nlocal;
532 if (!nlocal) return;
533 int i;
534 imageint idim,otherdims;
535 double *lo,*hi,*period;
536 double **x = atom->x;
537 double **v = atom->v;
538 int *mask = atom->mask;
539 imageint *image = atom->image;
540
541 // verify owned atoms have valid numerical coords
542 // may not if computed pairwise force between 2 atoms at same location
543
544 double *coord;
545 int n3 = 3*nlocal;
546 coord = &x[0][0];
547 int flag = 0;
548 for (i = 0; i < n3; i++)
549 if (!std::isfinite(*coord++)) flag = 1;
550 if (flag) error->one(FLERR,"Non-numeric atom coords - simulation unstable");
551
552 // setup for PBC checks
553
554 if (triclinic == 0) {
555 lo = boxlo;
556 hi = boxhi;
557 period = prd;
558 } else {
559 lo = boxlo_lamda;
560 hi = boxhi_lamda;
561 period = prd_lamda;
562 }
563
564 // apply PBC to each owned atom
565
566 for (i = 0; i < nlocal; i++) {
567 if (xperiodic) {
568 if (x[i][0] < lo[0]) {
569 x[i][0] += period[0];
570 if (deform_vremap && mask[i] & deform_groupbit) v[i][0] += h_rate[0];
571 idim = image[i] & IMGMASK;
572 otherdims = image[i] ^ idim;
573 idim--;
574 idim &= IMGMASK;
575 image[i] = otherdims | idim;
576 }
577 if (x[i][0] >= hi[0]) {
578 x[i][0] -= period[0];
579 x[i][0] = MAX(x[i][0],lo[0]);
580 if (deform_vremap && mask[i] & deform_groupbit) v[i][0] -= h_rate[0];
581 idim = image[i] & IMGMASK;
582 otherdims = image[i] ^ idim;
583 idim++;
584 idim &= IMGMASK;
585 image[i] = otherdims | idim;
586 }
587 }
588
589 if (yperiodic) {
590 if (x[i][1] < lo[1]) {
591 x[i][1] += period[1];
592 if (deform_vremap && mask[i] & deform_groupbit) {
593 v[i][0] += h_rate[5];
594 v[i][1] += h_rate[1];
595 }
596 idim = (image[i] >> IMGBITS) & IMGMASK;
597 otherdims = image[i] ^ (idim << IMGBITS);
598 idim--;
599 idim &= IMGMASK;
600 image[i] = otherdims | (idim << IMGBITS);
601 }
602 if (x[i][1] >= hi[1]) {
603 x[i][1] -= period[1];
604 x[i][1] = MAX(x[i][1],lo[1]);
605 if (deform_vremap && mask[i] & deform_groupbit) {
606 v[i][0] -= h_rate[5];
607 v[i][1] -= h_rate[1];
608 }
609 idim = (image[i] >> IMGBITS) & IMGMASK;
610 otherdims = image[i] ^ (idim << IMGBITS);
611 idim++;
612 idim &= IMGMASK;
613 image[i] = otherdims | (idim << IMGBITS);
614 }
615 }
616
617 if (zperiodic) {
618 if (x[i][2] < lo[2]) {
619 x[i][2] += period[2];
620 if (deform_vremap && mask[i] & deform_groupbit) {
621 v[i][0] += h_rate[4];
622 v[i][1] += h_rate[3];
623 v[i][2] += h_rate[2];
624 }
625 idim = image[i] >> IMG2BITS;
626 otherdims = image[i] ^ (idim << IMG2BITS);
627 idim--;
628 idim &= IMGMASK;
629 image[i] = otherdims | (idim << IMG2BITS);
630 }
631 if (x[i][2] >= hi[2]) {
632 x[i][2] -= period[2];
633 x[i][2] = MAX(x[i][2],lo[2]);
634 if (deform_vremap && mask[i] & deform_groupbit) {
635 v[i][0] -= h_rate[4];
636 v[i][1] -= h_rate[3];
637 v[i][2] -= h_rate[2];
638 }
639 idim = image[i] >> IMG2BITS;
640 otherdims = image[i] ^ (idim << IMG2BITS);
641 idim++;
642 idim &= IMGMASK;
643 image[i] = otherdims | (idim << IMG2BITS);
644 }
645 }
646 }
647 }
648
649 /* ----------------------------------------------------------------------
650 check that point is inside box boundaries, in [lo,hi) sense
651 return 1 if true, 0 if false
652 ------------------------------------------------------------------------- */
653
inside(double * x)654 int Domain::inside(double* x)
655 {
656 double *lo,*hi;
657 double lamda[3];
658
659 if (triclinic == 0) {
660 lo = boxlo;
661 hi = boxhi;
662
663 if (x[0] < lo[0] || x[0] >= hi[0] ||
664 x[1] < lo[1] || x[1] >= hi[1] ||
665 x[2] < lo[2] || x[2] >= hi[2]) return 0;
666 else return 1;
667
668 } else {
669 lo = boxlo_lamda;
670 hi = boxhi_lamda;
671
672 x2lamda(x,lamda);
673
674 if (lamda[0] < lo[0] || lamda[0] >= hi[0] ||
675 lamda[1] < lo[1] || lamda[1] >= hi[1] ||
676 lamda[2] < lo[2] || lamda[2] >= hi[2]) return 0;
677 else return 1;
678
679 }
680
681 }
682
683 /* ----------------------------------------------------------------------
684 check that point is inside nonperiodic boundaries, in [lo,hi) sense
685 return 1 if true, 0 if false
686 ------------------------------------------------------------------------- */
687
inside_nonperiodic(double * x)688 int Domain::inside_nonperiodic(double* x)
689 {
690 double *lo,*hi;
691 double lamda[3];
692
693 if (xperiodic && yperiodic && zperiodic) return 1;
694
695 if (triclinic == 0) {
696 lo = boxlo;
697 hi = boxhi;
698
699 if (!xperiodic && (x[0] < lo[0] || x[0] >= hi[0])) return 0;
700 if (!yperiodic && (x[1] < lo[1] || x[1] >= hi[1])) return 0;
701 if (!zperiodic && (x[2] < lo[2] || x[2] >= hi[2])) return 0;
702 return 1;
703
704 } else {
705 lo = boxlo_lamda;
706 hi = boxhi_lamda;
707
708 x2lamda(x,lamda);
709
710 if (!xperiodic && (lamda[0] < lo[0] || lamda[0] >= hi[0])) return 0;
711 if (!yperiodic && (lamda[1] < lo[1] || lamda[1] >= hi[1])) return 0;
712 if (!zperiodic && (lamda[2] < lo[2] || lamda[2] >= hi[2])) return 0;
713 return 1;
714 }
715
716 }
717
718 /* ----------------------------------------------------------------------
719 warn if image flags of any bonded atoms are inconsistent
720 could be a problem when using replicate or fix rigid
721 ------------------------------------------------------------------------- */
722
image_check()723 void Domain::image_check()
724 {
725 int i,j,k,n,imol,iatom;
726 tagint tagprev;
727
728 // only need to check if system is molecular and some dimension is periodic
729 // if running verlet/split, don't check on KSpace partition since
730 // it has no ghost atoms and thus bond partners won't exist
731
732 if (atom->molecular == Atom::ATOMIC) return;
733 if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return;
734 if (strncmp(update->integrate_style,"verlet/split",12) == 0 &&
735 universe->iworld != 0) return;
736
737 // communicate unwrapped position of owned atoms to ghost atoms
738
739 double **unwrap;
740 memory->create(unwrap,atom->nmax,3,"domain:unwrap");
741
742 double **x = atom->x;
743 imageint *image = atom->image;
744 int nlocal = atom->nlocal;
745
746 for (i = 0; i < nlocal; i++)
747 unmap(x[i],image[i],unwrap[i]);
748
749 comm->forward_comm_array(3,unwrap);
750
751 // compute unwrapped extent of each bond
752 // flag if any bond component is longer than 1/2 of periodic box length
753 // flag if any bond component is longer than non-periodic box length
754 // which means image flags in that dimension were different
755
756 int molecular = atom->molecular;
757
758 int *num_bond = atom->num_bond;
759 tagint **bond_atom = atom->bond_atom;
760 int **bond_type = atom->bond_type;
761 tagint *tag = atom->tag;
762 int *molindex = atom->molindex;
763 int *molatom = atom->molatom;
764 Molecule **onemols = atom->avec->onemols;
765
766 double delx,dely,delz;
767
768 int lostbond = output->thermo->lostbond;
769 int nmissing = 0;
770
771 int flag = 0;
772 for (i = 0; i < nlocal; i++) {
773 if (molecular == Atom::MOLECULAR) n = num_bond[i];
774 else {
775 if (molindex[i] < 0) continue;
776 imol = molindex[i];
777 iatom = molatom[i];
778 n = onemols[imol]->num_bond[iatom];
779 }
780
781 for (j = 0; j < n; j++) {
782 if (molecular == Atom::MOLECULAR) {
783 if (bond_type[i][j] <= 0) continue;
784 k = atom->map(bond_atom[i][j]);
785 } else {
786 if (onemols[imol]->bond_type[iatom][j] < 0) continue;
787 tagprev = tag[i] - iatom - 1;
788 k = atom->map(onemols[imol]->bond_atom[iatom][j]+tagprev);
789 }
790
791 if (k == -1) {
792 nmissing++;
793 if (lostbond == Thermo::ERROR)
794 error->one(FLERR,"Bond atom missing in image check");
795 continue;
796 }
797
798 delx = fabs(unwrap[i][0] - unwrap[k][0]);
799 dely = fabs(unwrap[i][1] - unwrap[k][1]);
800 delz = fabs(unwrap[i][2] - unwrap[k][2]);
801
802 if (xperiodic && delx > xprd_half) flag = 1;
803 if (yperiodic && dely > yprd_half) flag = 1;
804 if (dimension == 3 && zperiodic && delz > zprd_half) flag = 1;
805 if (!xperiodic && delx > xprd) flag = 1;
806 if (!yperiodic && dely > yprd) flag = 1;
807 if (dimension == 3 && !zperiodic && delz > zprd) flag = 1;
808 }
809 }
810
811 int flagall;
812 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
813 if (flagall && comm->me == 0)
814 error->warning(FLERR,"Inconsistent image flags");
815
816 if (lostbond == Thermo::WARN) {
817 int all;
818 MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
819 if (all && comm->me == 0)
820 error->warning(FLERR,"Bond atom missing in image check");
821 }
822
823 memory->destroy(unwrap);
824 }
825
826 /* ----------------------------------------------------------------------
827 warn if end atoms in any bonded interaction
828 are further apart than half a periodic box length
829 could cause problems when bonded neighbor list is built since
830 closest_image() could return wrong image
831 ------------------------------------------------------------------------- */
832
box_too_small_check()833 void Domain::box_too_small_check()
834 {
835 int i,j,k,n,imol,iatom;
836 tagint tagprev;
837
838 // only need to check if system is molecular and some dimension is periodic
839 // if running verlet/split, don't check on KSpace partition since
840 // it has no ghost atoms and thus bond partners won't exist
841
842 if (atom->molecular == Atom::ATOMIC) return;
843 if (!xperiodic && !yperiodic && (dimension == 2 || !zperiodic)) return;
844 if (strncmp(update->integrate_style,"verlet/split",12) == 0 &&
845 universe->iworld != 0) return;
846
847 // maxbondall = longest current bond length
848 // if periodic box dim is tiny (less than 2 * bond-length),
849 // minimum_image() itself may compute bad bond lengths
850 // in this case, image_check() should warn,
851 // assuming 2 atoms have consistent image flags
852
853 int molecular = atom->molecular;
854
855 double **x = atom->x;
856 int *num_bond = atom->num_bond;
857 tagint **bond_atom = atom->bond_atom;
858 int **bond_type = atom->bond_type;
859 tagint *tag = atom->tag;
860 int *molindex = atom->molindex;
861 int *molatom = atom->molatom;
862 Molecule **onemols = atom->avec->onemols;
863 int nlocal = atom->nlocal;
864
865 double delx,dely,delz,rsq;
866 double maxbondme = 0.0;
867
868 int lostbond = output->thermo->lostbond;
869 int nmissing = 0;
870
871 for (i = 0; i < nlocal; i++) {
872 if (molecular == Atom::MOLECULAR) n = num_bond[i];
873 else {
874 if (molindex[i] < 0) continue;
875 imol = molindex[i];
876 iatom = molatom[i];
877 n = onemols[imol]->num_bond[iatom];
878 }
879
880 for (j = 0; j < n; j++) {
881 if (molecular == Atom::MOLECULAR) {
882 if (bond_type[i][j] <= 0) continue;
883 k = atom->map(bond_atom[i][j]);
884 } else {
885 if (onemols[imol]->bond_type[iatom][j] < 0) continue;
886 tagprev = tag[i] - iatom - 1;
887 k = atom->map(onemols[imol]->bond_atom[iatom][j]+tagprev);
888 }
889
890 if (k == -1) {
891 nmissing++;
892 if (lostbond == Thermo::ERROR)
893 error->one(FLERR,"Bond atom missing in box size check");
894 continue;
895 }
896
897 delx = x[i][0] - x[k][0];
898 dely = x[i][1] - x[k][1];
899 delz = x[i][2] - x[k][2];
900 minimum_image(delx,dely,delz);
901 rsq = delx*delx + dely*dely + delz*delz;
902 maxbondme = MAX(maxbondme,rsq);
903 }
904 }
905
906 if (lostbond == Thermo::WARN) {
907 int all;
908 MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
909 if (all && comm->me == 0)
910 error->warning(FLERR,"Bond atom missing in box size check");
911 }
912
913 double maxbondall;
914 MPI_Allreduce(&maxbondme,&maxbondall,1,MPI_DOUBLE,MPI_MAX,world);
915 maxbondall = sqrt(maxbondall);
916
917 // maxdelta = furthest apart 2 atoms in a bonded interaction can be
918 // include BONDSTRETCH factor to account for dynamics
919
920 double maxdelta = maxbondall * BONDSTRETCH;
921 if (atom->nangles) maxdelta = 2.0 * maxbondall * BONDSTRETCH;
922 if (atom->ndihedrals) maxdelta = 3.0 * maxbondall * BONDSTRETCH;
923
924 // warn if maxdelta > than half any periodic box length
925 // since atoms in the interaction could rotate into that dimension
926
927 int flag = 0;
928 if (xperiodic && maxdelta > xprd_half) flag = 1;
929 if (yperiodic && maxdelta > yprd_half) flag = 1;
930 if (dimension == 3 && zperiodic && maxdelta > zprd_half) flag = 1;
931
932 int flagall;
933 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
934 if (flagall && comm->me == 0)
935 error->warning(FLERR,
936 "Bond/angle/dihedral extent > half of periodic box length");
937 }
938
939 /* ----------------------------------------------------------------------
940 check warn if any proc's subbox is smaller than thresh
941 since may lead to lost atoms in comm->exchange()
942 current callers set thresh = neighbor skin
943 ------------------------------------------------------------------------- */
944
subbox_too_small_check(double thresh)945 void Domain::subbox_too_small_check(double thresh)
946 {
947 int flag = 0;
948 if (!triclinic) {
949 if (subhi[0]-sublo[0] < thresh || subhi[1]-sublo[1] < thresh) flag = 1;
950 if (dimension == 3 && subhi[2]-sublo[2] < thresh) flag = 1;
951 } else {
952 double delta = subhi_lamda[0] - sublo_lamda[0];
953 if (delta*prd[0] < thresh) flag = 1;
954 delta = subhi_lamda[1] - sublo_lamda[1];
955 if (delta*prd[1] < thresh) flag = 1;
956 if (dimension == 3) {
957 delta = subhi_lamda[2] - sublo_lamda[2];
958 if (delta*prd[2] < thresh) flag = 1;
959 }
960 }
961
962 int flagall;
963 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
964 if (flagall && comm->me == 0)
965 error->warning(FLERR,"Proc sub-domain size < neighbor skin, "
966 "could lead to lost atoms");
967 }
968
969 /* ----------------------------------------------------------------------
970 minimum image convention in periodic dimensions
971 use 1/2 of box size as test
972 for triclinic, also add/subtract tilt factors in other dims as needed
973 changed "if" to "while" to enable distance to
974 far-away ghost atom returned by atom->map() to be wrapped back into box
975 could be problem for looking up atom IDs when cutoff > boxsize
976 this should not be used if atom has moved infinitely far outside box
977 b/c while could iterate forever
978 e.g. fix shake prediction of new position with highly overlapped atoms
979 use minimum_image_once() instead
980 ------------------------------------------------------------------------- */
981
minimum_image(double & dx,double & dy,double & dz)982 void Domain::minimum_image(double &dx, double &dy, double &dz)
983 {
984 if (triclinic == 0) {
985 if (xperiodic) {
986 while (fabs(dx) > xprd_half) {
987 if (dx < 0.0) dx += xprd;
988 else dx -= xprd;
989 }
990 }
991 if (yperiodic) {
992 while (fabs(dy) > yprd_half) {
993 if (dy < 0.0) dy += yprd;
994 else dy -= yprd;
995 }
996 }
997 if (zperiodic) {
998 while (fabs(dz) > zprd_half) {
999 if (dz < 0.0) dz += zprd;
1000 else dz -= zprd;
1001 }
1002 }
1003
1004 } else {
1005 if (zperiodic) {
1006 while (fabs(dz) > zprd_half) {
1007 if (dz < 0.0) {
1008 dz += zprd;
1009 dy += yz;
1010 dx += xz;
1011 } else {
1012 dz -= zprd;
1013 dy -= yz;
1014 dx -= xz;
1015 }
1016 }
1017 }
1018 if (yperiodic) {
1019 while (fabs(dy) > yprd_half) {
1020 if (dy < 0.0) {
1021 dy += yprd;
1022 dx += xy;
1023 } else {
1024 dy -= yprd;
1025 dx -= xy;
1026 }
1027 }
1028 }
1029 if (xperiodic) {
1030 while (fabs(dx) > xprd_half) {
1031 if (dx < 0.0) dx += xprd;
1032 else dx -= xprd;
1033 }
1034 }
1035 }
1036 }
1037
1038 /* ----------------------------------------------------------------------
1039 minimum image convention in periodic dimensions
1040 use 1/2 of box size as test
1041 for triclinic, also add/subtract tilt factors in other dims as needed
1042 changed "if" to "while" to enable distance to
1043 far-away ghost atom returned by atom->map() to be wrapped back into box
1044 could be problem for looking up atom IDs when cutoff > boxsize
1045 this should not be used if atom has moved infinitely far outside box
1046 b/c while could iterate forever
1047 e.g. fix shake prediction of new position with highly overlapped atoms
1048 use minimum_image_once() instead
1049 ------------------------------------------------------------------------- */
1050
minimum_image(double * delta)1051 void Domain::minimum_image(double *delta)
1052 {
1053 if (triclinic == 0) {
1054 if (xperiodic) {
1055 while (fabs(delta[0]) > xprd_half) {
1056 if (delta[0] < 0.0) delta[0] += xprd;
1057 else delta[0] -= xprd;
1058 }
1059 }
1060 if (yperiodic) {
1061 while (fabs(delta[1]) > yprd_half) {
1062 if (delta[1] < 0.0) delta[1] += yprd;
1063 else delta[1] -= yprd;
1064 }
1065 }
1066 if (zperiodic) {
1067 while (fabs(delta[2]) > zprd_half) {
1068 if (delta[2] < 0.0) delta[2] += zprd;
1069 else delta[2] -= zprd;
1070 }
1071 }
1072
1073 } else {
1074 if (zperiodic) {
1075 while (fabs(delta[2]) > zprd_half) {
1076 if (delta[2] < 0.0) {
1077 delta[2] += zprd;
1078 delta[1] += yz;
1079 delta[0] += xz;
1080 } else {
1081 delta[2] -= zprd;
1082 delta[1] -= yz;
1083 delta[0] -= xz;
1084 }
1085 }
1086 }
1087 if (yperiodic) {
1088 while (fabs(delta[1]) > yprd_half) {
1089 if (delta[1] < 0.0) {
1090 delta[1] += yprd;
1091 delta[0] += xy;
1092 } else {
1093 delta[1] -= yprd;
1094 delta[0] -= xy;
1095 }
1096 }
1097 }
1098 if (xperiodic) {
1099 while (fabs(delta[0]) > xprd_half) {
1100 if (delta[0] < 0.0) delta[0] += xprd;
1101 else delta[0] -= xprd;
1102 }
1103 }
1104 }
1105 }
1106
1107 /* ----------------------------------------------------------------------
1108 minimum image convention in periodic dimensions
1109 use 1/2 of box size as test
1110 for triclinic, also add/subtract tilt factors in other dims as needed
1111 only shift by one box length in each direction
1112 this should not be used if multiple box shifts are required
1113 ------------------------------------------------------------------------- */
1114
minimum_image_once(double * delta)1115 void Domain::minimum_image_once(double *delta)
1116 {
1117 if (triclinic == 0) {
1118 if (xperiodic) {
1119 if (fabs(delta[0]) > xprd_half) {
1120 if (delta[0] < 0.0) delta[0] += xprd;
1121 else delta[0] -= xprd;
1122 }
1123 }
1124 if (yperiodic) {
1125 if (fabs(delta[1]) > yprd_half) {
1126 if (delta[1] < 0.0) delta[1] += yprd;
1127 else delta[1] -= yprd;
1128 }
1129 }
1130 if (zperiodic) {
1131 if (fabs(delta[2]) > zprd_half) {
1132 if (delta[2] < 0.0) delta[2] += zprd;
1133 else delta[2] -= zprd;
1134 }
1135 }
1136
1137 } else {
1138 if (zperiodic) {
1139 if (fabs(delta[2]) > zprd_half) {
1140 if (delta[2] < 0.0) {
1141 delta[2] += zprd;
1142 delta[1] += yz;
1143 delta[0] += xz;
1144 } else {
1145 delta[2] -= zprd;
1146 delta[1] -= yz;
1147 delta[0] -= xz;
1148 }
1149 }
1150 }
1151 if (yperiodic) {
1152 if (fabs(delta[1]) > yprd_half) {
1153 if (delta[1] < 0.0) {
1154 delta[1] += yprd;
1155 delta[0] += xy;
1156 } else {
1157 delta[1] -= yprd;
1158 delta[0] -= xy;
1159 }
1160 }
1161 }
1162 if (xperiodic) {
1163 if (fabs(delta[0]) > xprd_half) {
1164 if (delta[0] < 0.0) delta[0] += xprd;
1165 else delta[0] -= xprd;
1166 }
1167 }
1168 }
1169 }
1170
1171 /* ----------------------------------------------------------------------
1172 return local index of atom J or any of its images that is closest to atom I
1173 if J is not a valid index like -1, just return it
1174 ------------------------------------------------------------------------- */
1175
closest_image(int i,int j)1176 int Domain::closest_image(int i, int j)
1177 {
1178 if (j < 0) return j;
1179
1180 int *sametag = atom->sametag;
1181 double **x = atom->x;
1182 double *xi = x[i];
1183
1184 int closest = j;
1185 double delx = xi[0] - x[j][0];
1186 double dely = xi[1] - x[j][1];
1187 double delz = xi[2] - x[j][2];
1188 double rsqmin = delx*delx + dely*dely + delz*delz;
1189 double rsq;
1190
1191 while (sametag[j] >= 0) {
1192 j = sametag[j];
1193 delx = xi[0] - x[j][0];
1194 dely = xi[1] - x[j][1];
1195 delz = xi[2] - x[j][2];
1196 rsq = delx*delx + dely*dely + delz*delz;
1197 if (rsq < rsqmin) {
1198 rsqmin = rsq;
1199 closest = j;
1200 }
1201 }
1202
1203 return closest;
1204 }
1205
1206 /* ----------------------------------------------------------------------
1207 return local index of atom J or any of its images that is closest to pos
1208 if J is not a valid index like -1, just return it
1209 ------------------------------------------------------------------------- */
1210
closest_image(const double * const pos,int j)1211 int Domain::closest_image(const double * const pos, int j)
1212 {
1213 if (j < 0) return j;
1214
1215 const int * const sametag = atom->sametag;
1216 const double * const * const x = atom->x;
1217
1218 int closest = j;
1219 double delx = pos[0] - x[j][0];
1220 double dely = pos[1] - x[j][1];
1221 double delz = pos[2] - x[j][2];
1222 double rsqmin = delx*delx + dely*dely + delz*delz;
1223 double rsq;
1224
1225 while (sametag[j] >= 0) {
1226 j = sametag[j];
1227 delx = pos[0] - x[j][0];
1228 dely = pos[1] - x[j][1];
1229 delz = pos[2] - x[j][2];
1230 rsq = delx*delx + dely*dely + delz*delz;
1231 if (rsq < rsqmin) {
1232 rsqmin = rsq;
1233 closest = j;
1234 }
1235 }
1236
1237 return closest;
1238 }
1239
1240 /* ----------------------------------------------------------------------
1241 find and return Xj image = periodic image of Xj that is closest to Xi
1242 for triclinic, add/subtract tilt factors in other dims as needed
1243 called by ServerMD class and LammpsInterface in lib/atc.
1244 ------------------------------------------------------------------------- */
1245
closest_image(const double * const xi,const double * const xj,double * const xjimage)1246 void Domain::closest_image(const double * const xi, const double * const xj, double * const xjimage)
1247 {
1248 double dx = xj[0] - xi[0];
1249 double dy = xj[1] - xi[1];
1250 double dz = xj[2] - xi[2];
1251
1252 if (triclinic == 0) {
1253 if (xperiodic) {
1254 if (dx < 0.0) {
1255 while (dx < 0.0) dx += xprd;
1256 if (dx > xprd_half) dx -= xprd;
1257 } else {
1258 while (dx > 0.0) dx -= xprd;
1259 if (dx < -xprd_half) dx += xprd;
1260 }
1261 }
1262 if (yperiodic) {
1263 if (dy < 0.0) {
1264 while (dy < 0.0) dy += yprd;
1265 if (dy > yprd_half) dy -= yprd;
1266 } else {
1267 while (dy > 0.0) dy -= yprd;
1268 if (dy < -yprd_half) dy += yprd;
1269 }
1270 }
1271 if (zperiodic) {
1272 if (dz < 0.0) {
1273 while (dz < 0.0) dz += zprd;
1274 if (dz > zprd_half) dz -= zprd;
1275 } else {
1276 while (dz > 0.0) dz -= zprd;
1277 if (dz < -zprd_half) dz += zprd;
1278 }
1279 }
1280
1281 } else {
1282 if (zperiodic) {
1283 if (dz < 0.0) {
1284 while (dz < 0.0) {
1285 dz += zprd;
1286 dy += yz;
1287 dx += xz;
1288 }
1289 if (dz > zprd_half) {
1290 dz -= zprd;
1291 dy -= yz;
1292 dx -= xz;
1293 }
1294 } else {
1295 while (dz > 0.0) {
1296 dz -= zprd;
1297 dy -= yz;
1298 dx -= xz;
1299 }
1300 if (dz < -zprd_half) {
1301 dz += zprd;
1302 dy += yz;
1303 dx += xz;
1304 }
1305 }
1306 }
1307 if (yperiodic) {
1308 if (dy < 0.0) {
1309 while (dy < 0.0) {
1310 dy += yprd;
1311 dx += xy;
1312 }
1313 if (dy > yprd_half) {
1314 dy -= yprd;
1315 dx -= xy;
1316 }
1317 } else {
1318 while (dy > 0.0) {
1319 dy -= yprd;
1320 dx -= xy;
1321 }
1322 if (dy < -yprd_half) {
1323 dy += yprd;
1324 dx += xy;
1325 }
1326 }
1327 }
1328 if (xperiodic) {
1329 if (dx < 0.0) {
1330 while (dx < 0.0) dx += xprd;
1331 if (dx > xprd_half) dx -= xprd;
1332 } else {
1333 while (dx > 0.0) dx -= xprd;
1334 if (dx < -xprd_half) dx += xprd;
1335 }
1336 }
1337 }
1338
1339 xjimage[0] = xi[0] + dx;
1340 xjimage[1] = xi[1] + dy;
1341 xjimage[2] = xi[2] + dz;
1342 }
1343
1344 /* ----------------------------------------------------------------------
1345 remap the point into the periodic box no matter how far away
1346 adjust 3 image flags encoded in image accordingly
1347 resulting coord must satisfy lo <= coord < hi
1348 MAX is important since coord - prd < lo can happen when coord = hi
1349 for triclinic, point is converted to lamda coords (0-1) before doing remap
1350 image = 10 bits for each dimension
1351 increment/decrement in wrap-around fashion
1352 ------------------------------------------------------------------------- */
1353
remap(double * x,imageint & image)1354 void Domain::remap(double *x, imageint &image)
1355 {
1356 double *lo,*hi,*period,*coord;
1357 double lamda[3];
1358 imageint idim,otherdims;
1359
1360 if (triclinic == 0) {
1361 lo = boxlo;
1362 hi = boxhi;
1363 period = prd;
1364 coord = x;
1365 } else {
1366 lo = boxlo_lamda;
1367 hi = boxhi_lamda;
1368 period = prd_lamda;
1369 x2lamda(x,lamda);
1370 coord = lamda;
1371 }
1372
1373 if (xperiodic) {
1374 while (coord[0] < lo[0]) {
1375 coord[0] += period[0];
1376 idim = image & IMGMASK;
1377 otherdims = image ^ idim;
1378 idim--;
1379 idim &= IMGMASK;
1380 image = otherdims | idim;
1381 }
1382 while (coord[0] >= hi[0]) {
1383 coord[0] -= period[0];
1384 idim = image & IMGMASK;
1385 otherdims = image ^ idim;
1386 idim++;
1387 idim &= IMGMASK;
1388 image = otherdims | idim;
1389 }
1390 coord[0] = MAX(coord[0],lo[0]);
1391 }
1392
1393 if (yperiodic) {
1394 while (coord[1] < lo[1]) {
1395 coord[1] += period[1];
1396 idim = (image >> IMGBITS) & IMGMASK;
1397 otherdims = image ^ (idim << IMGBITS);
1398 idim--;
1399 idim &= IMGMASK;
1400 image = otherdims | (idim << IMGBITS);
1401 }
1402 while (coord[1] >= hi[1]) {
1403 coord[1] -= period[1];
1404 idim = (image >> IMGBITS) & IMGMASK;
1405 otherdims = image ^ (idim << IMGBITS);
1406 idim++;
1407 idim &= IMGMASK;
1408 image = otherdims | (idim << IMGBITS);
1409 }
1410 coord[1] = MAX(coord[1],lo[1]);
1411 }
1412
1413 if (zperiodic) {
1414 while (coord[2] < lo[2]) {
1415 coord[2] += period[2];
1416 idim = image >> IMG2BITS;
1417 otherdims = image ^ (idim << IMG2BITS);
1418 idim--;
1419 idim &= IMGMASK;
1420 image = otherdims | (idim << IMG2BITS);
1421 }
1422 while (coord[2] >= hi[2]) {
1423 coord[2] -= period[2];
1424 idim = image >> IMG2BITS;
1425 otherdims = image ^ (idim << IMG2BITS);
1426 idim++;
1427 idim &= IMGMASK;
1428 image = otherdims | (idim << IMG2BITS);
1429 }
1430 coord[2] = MAX(coord[2],lo[2]);
1431 }
1432
1433 if (triclinic) lamda2x(coord,x);
1434 }
1435
1436 /* ----------------------------------------------------------------------
1437 remap the point into the periodic box no matter how far away
1438 no image flag calculation
1439 resulting coord must satisfy lo <= coord < hi
1440 MAX is important since coord - prd < lo can happen when coord = hi
1441 for triclinic, point is converted to lamda coords (0-1) before remap
1442 ------------------------------------------------------------------------- */
1443
remap(double * x)1444 void Domain::remap(double *x)
1445 {
1446 double *lo,*hi,*period,*coord;
1447 double lamda[3];
1448
1449 if (triclinic == 0) {
1450 lo = boxlo;
1451 hi = boxhi;
1452 period = prd;
1453 coord = x;
1454 } else {
1455 lo = boxlo_lamda;
1456 hi = boxhi_lamda;
1457 period = prd_lamda;
1458 x2lamda(x,lamda);
1459 coord = lamda;
1460 }
1461
1462 if (xperiodic) {
1463 while (coord[0] < lo[0]) coord[0] += period[0];
1464 while (coord[0] >= hi[0]) coord[0] -= period[0];
1465 coord[0] = MAX(coord[0],lo[0]);
1466 }
1467
1468 if (yperiodic) {
1469 while (coord[1] < lo[1]) coord[1] += period[1];
1470 while (coord[1] >= hi[1]) coord[1] -= period[1];
1471 coord[1] = MAX(coord[1],lo[1]);
1472 }
1473
1474 if (zperiodic) {
1475 while (coord[2] < lo[2]) coord[2] += period[2];
1476 while (coord[2] >= hi[2]) coord[2] -= period[2];
1477 coord[2] = MAX(coord[2],lo[2]);
1478 }
1479
1480 if (triclinic) lamda2x(coord,x);
1481 }
1482
1483 /* ----------------------------------------------------------------------
1484 remap xnew to be within half box length of xold
1485 do it directly, not iteratively, in case is far away
1486 for triclinic, both points are converted to lamda coords (0-1) before remap
1487 ------------------------------------------------------------------------- */
1488
remap_near(double * xnew,double * xold)1489 void Domain::remap_near(double *xnew, double *xold)
1490 {
1491 int n;
1492 double *coordnew,*coordold,*period,*half;
1493 double lamdanew[3],lamdaold[3];
1494
1495 if (triclinic == 0) {
1496 period = prd;
1497 half = prd_half;
1498 coordnew = xnew;
1499 coordold = xold;
1500 } else {
1501 period = prd_lamda;
1502 half = prd_half_lamda;
1503 x2lamda(xnew,lamdanew);
1504 coordnew = lamdanew;
1505 x2lamda(xold,lamdaold);
1506 coordold = lamdaold;
1507 }
1508
1509 // iterative form
1510 // if (xperiodic) {
1511 // while (coordnew[0]-coordold[0] > half[0]) coordnew[0] -= period[0];
1512 // while (coordold[0]-coordnew[0] > half[0]) coordnew[0] += period[0];
1513 // }
1514
1515 if (xperiodic) {
1516 if (coordnew[0]-coordold[0] > period[0]) {
1517 n = static_cast<int> ((coordnew[0]-coordold[0])/period[0]);
1518 coordnew[0] -= n*period[0];
1519 }
1520 while (coordnew[0]-coordold[0] > half[0]) coordnew[0] -= period[0];
1521 if (coordold[0]-coordnew[0] > period[0]) {
1522 n = static_cast<int> ((coordold[0]-coordnew[0])/period[0]);
1523 coordnew[0] += n*period[0];
1524 }
1525 while (coordold[0]-coordnew[0] > half[0]) coordnew[0] += period[0];
1526 }
1527
1528 if (yperiodic) {
1529 if (coordnew[1]-coordold[1] > period[1]) {
1530 n = static_cast<int> ((coordnew[1]-coordold[1])/period[1]);
1531 coordnew[1] -= n*period[1];
1532 }
1533 while (coordnew[1]-coordold[1] > half[1]) coordnew[1] -= period[1];
1534 if (coordold[1]-coordnew[1] > period[1]) {
1535 n = static_cast<int> ((coordold[1]-coordnew[1])/period[1]);
1536 coordnew[1] += n*period[1];
1537 }
1538 while (coordold[1]-coordnew[1] > half[1]) coordnew[1] += period[1];
1539 }
1540
1541 if (zperiodic) {
1542 if (coordnew[2]-coordold[2] > period[2]) {
1543 n = static_cast<int> ((coordnew[2]-coordold[2])/period[2]);
1544 coordnew[2] -= n*period[2];
1545 }
1546 while (coordnew[2]-coordold[2] > half[2]) coordnew[2] -= period[2];
1547 if (coordold[2]-coordnew[2] > period[2]) {
1548 n = static_cast<int> ((coordold[2]-coordnew[2])/period[2]);
1549 coordnew[2] += n*period[2];
1550 }
1551 while (coordold[2]-coordnew[2] > half[2]) coordnew[2] += period[2];
1552 }
1553
1554 if (triclinic) lamda2x(coordnew,xnew);
1555 }
1556
1557 /* ----------------------------------------------------------------------
1558 unmap the point via image flags
1559 x overwritten with result, don't reset image flag
1560 for triclinic, use h[] to add in tilt factors in other dims as needed
1561 ------------------------------------------------------------------------- */
1562
unmap(double * x,imageint image)1563 void Domain::unmap(double *x, imageint image)
1564 {
1565 int xbox = (image & IMGMASK) - IMGMAX;
1566 int ybox = (image >> IMGBITS & IMGMASK) - IMGMAX;
1567 int zbox = (image >> IMG2BITS) - IMGMAX;
1568
1569 if (triclinic == 0) {
1570 x[0] += xbox*xprd;
1571 x[1] += ybox*yprd;
1572 x[2] += zbox*zprd;
1573 } else {
1574 x[0] += h[0]*xbox + h[5]*ybox + h[4]*zbox;
1575 x[1] += h[1]*ybox + h[3]*zbox;
1576 x[2] += h[2]*zbox;
1577 }
1578 }
1579
1580 /* ----------------------------------------------------------------------
1581 unmap the point via image flags
1582 result returned in y, don't reset image flag
1583 for triclinic, use h[] to add in tilt factors in other dims as needed
1584 ------------------------------------------------------------------------- */
1585
unmap(const double * x,imageint image,double * y)1586 void Domain::unmap(const double *x, imageint image, double *y)
1587 {
1588 int xbox = (image & IMGMASK) - IMGMAX;
1589 int ybox = (image >> IMGBITS & IMGMASK) - IMGMAX;
1590 int zbox = (image >> IMG2BITS) - IMGMAX;
1591
1592 if (triclinic == 0) {
1593 y[0] = x[0] + xbox*xprd;
1594 y[1] = x[1] + ybox*yprd;
1595 y[2] = x[2] + zbox*zprd;
1596 } else {
1597 y[0] = x[0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
1598 y[1] = x[1] + h[1]*ybox + h[3]*zbox;
1599 y[2] = x[2] + h[2]*zbox;
1600 }
1601 }
1602
1603 /* ----------------------------------------------------------------------
1604 adjust image flags due to triclinic box flip
1605 flip operation is changing box vectors A,B,C to new A',B',C'
1606 A' = A (A does not change)
1607 B' = B + mA (B shifted by A)
1608 C' = C + pB + nA (C shifted by B and/or A)
1609 this requires the image flags change from (a,b,c) to (a',b',c')
1610 so that x_unwrap for each atom is same before/after
1611 x_unwrap_before = xlocal + aA + bB + cC
1612 x_unwrap_after = xlocal + a'A' + b'B' + c'C'
1613 this requires:
1614 c' = c
1615 b' = b - cp
1616 a' = a - (b-cp)m - cn = a - b'm - cn
1617 in other words, for xy flip, change in x flag depends on current y flag
1618 this is b/c the xy flip dramatically changes which tiled image of
1619 simulation box an unwrapped point maps to
1620 ------------------------------------------------------------------------- */
1621
image_flip(int m,int n,int p)1622 void Domain::image_flip(int m, int n, int p)
1623 {
1624 imageint *image = atom->image;
1625 int nlocal = atom->nlocal;
1626
1627 for (int i = 0; i < nlocal; i++) {
1628 int xbox = (image[i] & IMGMASK) - IMGMAX;
1629 int ybox = (image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1630 int zbox = (image[i] >> IMG2BITS) - IMGMAX;
1631
1632 ybox -= p*zbox;
1633 xbox -= m*ybox + n*zbox;
1634
1635 image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
1636 (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
1637 (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
1638 }
1639 }
1640
1641 /* ----------------------------------------------------------------------
1642 return 1 if this proc owns atom with coords x, else return 0
1643 x is returned remapped into periodic box
1644 if image flag is passed, flag is updated via remap(x,image)
1645 if image = nullptr is passed, no update with remap(x)
1646 if shrinkexceed, atom can be outside shrinkwrap boundaries
1647 called from create_atoms() in library.cpp
1648 ------------------------------------------------------------------------- */
1649
ownatom(int,double * x,imageint * image,int shrinkexceed)1650 int Domain::ownatom(int /*id*/, double *x, imageint *image, int shrinkexceed)
1651 {
1652 double lamda[3];
1653 double *coord,*blo,*bhi,*slo,*shi;
1654
1655 if (image) remap(x,*image);
1656 else remap(x);
1657
1658 // if triclinic, convert to lamda coords (0-1)
1659 // for periodic dims, resulting coord must satisfy 0.0 <= coord < 1.0
1660
1661 if (triclinic) {
1662 x2lamda(x,lamda);
1663 if (xperiodic && (lamda[0] < 0.0 || lamda[0] >= 1.0)) lamda[0] = 0.0;
1664 if (yperiodic && (lamda[1] < 0.0 || lamda[1] >= 1.0)) lamda[1] = 0.0;
1665 if (zperiodic && (lamda[2] < 0.0 || lamda[2] >= 1.0)) lamda[2] = 0.0;
1666 coord = lamda;
1667 } else coord = x;
1668
1669 // box and subbox bounds for orthogonal vs triclinic
1670
1671 if (triclinic == 0) {
1672 blo = boxlo;
1673 bhi = boxhi;
1674 slo = sublo;
1675 shi = subhi;
1676 } else {
1677 blo = boxlo_lamda;
1678 bhi = boxhi_lamda;
1679 slo = sublo_lamda;
1680 shi = subhi_lamda;
1681 }
1682
1683 if (coord[0] >= slo[0] && coord[0] < shi[0] &&
1684 coord[1] >= slo[1] && coord[1] < shi[1] &&
1685 coord[2] >= slo[2] && coord[2] < shi[2]) return 1;
1686
1687 // check if atom did not return 1 only b/c it was
1688 // outside a shrink-wrapped boundary
1689
1690 if (shrinkexceed) {
1691 int outside = 0;
1692 if (coord[0] < blo[0] && boundary[0][0] > 1) outside = 1;
1693 if (coord[0] >= bhi[0] && boundary[0][1] > 1) outside = 1;
1694 if (coord[1] < blo[1] && boundary[1][0] > 1) outside = 1;
1695 if (coord[1] >= bhi[1] && boundary[1][1] > 1) outside = 1;
1696 if (coord[2] < blo[2] && boundary[2][0] > 1) outside = 1;
1697 if (coord[2] >= bhi[2] && boundary[2][1] > 1) outside = 1;
1698 if (!outside) return 0;
1699
1700 // newcoord = coords pushed back to be on shrink-wrapped boundary
1701 // newcoord is a copy, so caller's x[] is not affected
1702
1703 double newcoord[3];
1704 if (coord[0] < blo[0] && boundary[0][0] > 1) newcoord[0] = blo[0];
1705 else if (coord[0] >= bhi[0] && boundary[0][1] > 1) newcoord[0] = bhi[0];
1706 else newcoord[0] = coord[0];
1707 if (coord[1] < blo[1] && boundary[1][0] > 1) newcoord[1] = blo[1];
1708 else if (coord[1] >= bhi[1] && boundary[1][1] > 1) newcoord[1] = bhi[1];
1709 else newcoord[1] = coord[1];
1710 if (coord[2] < blo[2] && boundary[2][0] > 1) newcoord[2] = blo[2];
1711 else if (coord[2] >= bhi[2] && boundary[2][1] > 1) newcoord[2] = bhi[2];
1712 else newcoord[2] = coord[2];
1713
1714 // re-test for newcoord inside my sub-domain
1715 // use <= test for upper-boundary since may have just put atom at boxhi
1716
1717 if (newcoord[0] >= slo[0] && newcoord[0] <= shi[0] &&
1718 newcoord[1] >= slo[1] && newcoord[1] <= shi[1] &&
1719 newcoord[2] >= slo[2] && newcoord[2] <= shi[2]) return 1;
1720 }
1721
1722 return 0;
1723 }
1724
1725 /* ----------------------------------------------------------------------
1726 create a lattice
1727 ------------------------------------------------------------------------- */
1728
set_lattice(int narg,char ** arg)1729 void Domain::set_lattice(int narg, char **arg)
1730 {
1731 if (lattice) delete lattice;
1732 lattice = nullptr;
1733 lattice = new Lattice(lmp,narg,arg);
1734 }
1735
1736 /* ----------------------------------------------------------------------
1737 create a new region
1738 ------------------------------------------------------------------------- */
1739
add_region(int narg,char ** arg)1740 void Domain::add_region(int narg, char **arg)
1741 {
1742 if (narg < 2) error->all(FLERR,"Illegal region command");
1743
1744 if (strcmp(arg[1],"delete") == 0) {
1745 delete_region(narg,arg);
1746 return;
1747 }
1748
1749 if (strcmp(arg[1],"none") == 0)
1750 error->all(FLERR,"Unrecognized region style 'none'");
1751
1752 if (find_region(arg[0]) >= 0) error->all(FLERR,"Reuse of region ID");
1753
1754 // extend Region list if necessary
1755
1756 if (nregion == maxregion) {
1757 maxregion += DELTAREGION;
1758 regions = (Region **)
1759 memory->srealloc(regions,maxregion*sizeof(Region *),"domain:regions");
1760 }
1761
1762 // create the Region
1763
1764 if (lmp->suffix_enable) {
1765 if (lmp->suffix) {
1766 std::string estyle = std::string(arg[1]) + "/" + lmp->suffix;
1767 if (region_map->find(estyle) != region_map->end()) {
1768 RegionCreator ®ion_creator = (*region_map)[estyle];
1769 regions[nregion] = region_creator(lmp, narg, arg);
1770 regions[nregion]->init();
1771 nregion++;
1772 return;
1773 }
1774 }
1775
1776 if (lmp->suffix2) {
1777 std::string estyle = std::string(arg[1]) + "/" + lmp->suffix2;
1778 if (region_map->find(estyle) != region_map->end()) {
1779 RegionCreator ®ion_creator = (*region_map)[estyle];
1780 regions[nregion] = region_creator(lmp, narg, arg);
1781 regions[nregion]->init();
1782 nregion++;
1783 return;
1784 }
1785 }
1786 }
1787
1788 if (region_map->find(arg[1]) != region_map->end()) {
1789 RegionCreator ®ion_creator = (*region_map)[arg[1]];
1790 regions[nregion] = region_creator(lmp, narg, arg);
1791 } else error->all(FLERR,utils::check_packages_for_style("region",arg[1],lmp));
1792
1793 // initialize any region variables via init()
1794 // in case region is used between runs, e.g. to print a variable
1795
1796 regions[nregion]->init();
1797 nregion++;
1798 }
1799
1800 /* ----------------------------------------------------------------------
1801 one instance per region style in style_region.h
1802 ------------------------------------------------------------------------- */
1803
1804 template <typename T>
region_creator(LAMMPS * lmp,int narg,char ** arg)1805 Region *Domain::region_creator(LAMMPS *lmp, int narg, char ** arg)
1806 {
1807 return new T(lmp, narg, arg);
1808 }
1809
1810 /* ----------------------------------------------------------------------
1811 delete a region
1812 ------------------------------------------------------------------------- */
1813
delete_region(int narg,char ** arg)1814 void Domain::delete_region(int narg, char **arg)
1815 {
1816 if (narg != 2) error->all(FLERR,"Illegal region command");
1817
1818 int iregion = find_region(arg[0]);
1819 if (iregion == -1) error->all(FLERR,"Delete region ID does not exist");
1820
1821 delete_region(iregion);
1822 }
1823
delete_region(int iregion)1824 void Domain::delete_region(int iregion)
1825 {
1826 if ((iregion < 0) || (iregion >= nregion)) return;
1827
1828 // delete and move other Regions down in list one slot
1829
1830 delete regions[iregion];
1831 for (int i = iregion+1; i < nregion; ++i)
1832 regions[i-1] = regions[i];
1833 nregion--;
1834 }
1835
1836 /* ----------------------------------------------------------------------
1837 return region index if name matches existing region ID
1838 return -1 if no such region
1839 ------------------------------------------------------------------------- */
1840
find_region(const std::string & name)1841 int Domain::find_region(const std::string &name)
1842 {
1843 for (int iregion = 0; iregion < nregion; iregion++)
1844 if (name == regions[iregion]->id) return iregion;
1845 return -1;
1846 }
1847
1848 /* ----------------------------------------------------------------------
1849 return region index if name matches existing region style
1850 return -1 if no such region
1851 ------------------------------------------------------------------------- */
1852
find_region_by_style(const std::string & name)1853 int Domain::find_region_by_style(const std::string &name)
1854 {
1855 for (int iregion = 0; iregion < nregion; iregion++)
1856 if (name == regions[iregion]->style) return iregion;
1857 return -1;
1858 }
1859
1860 /* ----------------------------------------------------------------------
1861 (re)set boundary settings
1862 flag = 0, called from the input script
1863 flag = 1, called from change box command
1864 ------------------------------------------------------------------------- */
1865
set_boundary(int narg,char ** arg,int flag)1866 void Domain::set_boundary(int narg, char **arg, int flag)
1867 {
1868 if (narg != 3) error->all(FLERR,"Illegal boundary command");
1869
1870 char c;
1871 for (int idim = 0; idim < 3; idim++)
1872 for (int iside = 0; iside < 2; iside++) {
1873 if (iside == 0) c = arg[idim][0];
1874 else if (iside == 1 && strlen(arg[idim]) == 1) c = arg[idim][0];
1875 else c = arg[idim][1];
1876
1877 if (c == 'p') boundary[idim][iside] = 0;
1878 else if (c == 'f') boundary[idim][iside] = 1;
1879 else if (c == 's') boundary[idim][iside] = 2;
1880 else if (c == 'm') boundary[idim][iside] = 3;
1881 else {
1882 if (flag == 0) error->all(FLERR,"Illegal boundary command");
1883 if (flag == 1) error->all(FLERR,"Illegal change_box command");
1884 }
1885 }
1886
1887 for (int idim = 0; idim < 3; idim++)
1888 if ((boundary[idim][0] == 0 && boundary[idim][1]) ||
1889 (boundary[idim][0] && boundary[idim][1] == 0))
1890 error->all(FLERR,"Both sides of boundary must be periodic");
1891
1892 if (boundary[0][0] == 0) xperiodic = 1;
1893 else xperiodic = 0;
1894 if (boundary[1][0] == 0) yperiodic = 1;
1895 else yperiodic = 0;
1896 if (boundary[2][0] == 0) zperiodic = 1;
1897 else zperiodic = 0;
1898
1899 // record if we changed a periodic boundary to a non-periodic one
1900
1901 int pflag=0;
1902 if ((periodicity[0] && !xperiodic)
1903 || (periodicity[1] && !yperiodic)
1904 || (periodicity[2] && !zperiodic)) pflag = 1;
1905
1906 periodicity[0] = xperiodic;
1907 periodicity[1] = yperiodic;
1908 periodicity[2] = zperiodic;
1909
1910 nonperiodic = 0;
1911 if (xperiodic == 0 || yperiodic == 0 || zperiodic == 0) {
1912 nonperiodic = 1;
1913 if (boundary[0][0] >= 2 || boundary[0][1] >= 2 ||
1914 boundary[1][0] >= 2 || boundary[1][1] >= 2 ||
1915 boundary[2][0] >= 2 || boundary[2][1] >= 2) nonperiodic = 2;
1916 }
1917
1918 // force non-zero image flags to zero for non-periodic dimensions
1919 // keep track if a change was made, so we can print a warning message
1920
1921 if (pflag) {
1922 pflag = 0;
1923 for (int i=0; i < atom->nlocal; ++i) {
1924 int xbox = (atom->image[i] & IMGMASK) - IMGMAX;
1925 int ybox = (atom->image[i] >> IMGBITS & IMGMASK) - IMGMAX;
1926 int zbox = (atom->image[i] >> IMG2BITS) - IMGMAX;
1927 if ((!xperiodic) && (xbox != 0)) { xbox = 0; pflag = 1; }
1928 if ((!yperiodic) && (ybox != 0)) { ybox = 0; pflag = 1; }
1929 if ((!zperiodic) && (zbox != 0)) { zbox = 0; pflag = 1; }
1930 atom->image[i] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
1931 (((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
1932 (((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
1933 }
1934 int flag_all;
1935 MPI_Allreduce(&pflag,&flag_all, 1, MPI_INT, MPI_SUM, world);
1936 if ((flag_all > 0) && (comm->me == 0))
1937 error->warning(FLERR,"Resetting image flags for non-periodic dimensions");
1938 }
1939 }
1940
1941 /* ----------------------------------------------------------------------
1942 set domain attributes
1943 ------------------------------------------------------------------------- */
1944
set_box(int narg,char ** arg)1945 void Domain::set_box(int narg, char **arg)
1946 {
1947 if (narg < 1) error->all(FLERR,"Illegal box command");
1948
1949 int iarg = 0;
1950 while (iarg < narg) {
1951 if (strcmp(arg[iarg],"tilt") == 0) {
1952 if (iarg+2 > narg) error->all(FLERR,"Illegal box command");
1953 if (strcmp(arg[iarg+1],"small") == 0) tiltsmall = 1;
1954 else if (strcmp(arg[iarg+1],"large") == 0) tiltsmall = 0;
1955 else error->all(FLERR,"Illegal box command");
1956 iarg += 2;
1957 } else error->all(FLERR,"Illegal box command");
1958 }
1959 }
1960
1961 /* ----------------------------------------------------------------------
1962 print box info, orthogonal or triclinic
1963 ------------------------------------------------------------------------- */
1964
print_box(const std::string & prefix)1965 void Domain::print_box(const std::string &prefix)
1966 {
1967 if (comm->me == 0) {
1968 std::string mesg = prefix;
1969 if (triclinic == 0) {
1970 mesg += fmt::format("orthogonal box = ({:.8} {:.8} {:.8}) to "
1971 "({:.8} {:.8} {:.8})\n",boxlo[0],boxlo[1],
1972 boxlo[2],boxhi[0],boxhi[1],boxhi[2]);
1973 } else {
1974 mesg += fmt::format("triclinic box = ({:.8} {:.8} {:.8}) to "
1975 "({:.8} {:.8} {:.8}) with tilt "
1976 "({:.8} {:.8} {:.8})\n",boxlo[0],boxlo[1],
1977 boxlo[2],boxhi[0],boxhi[1],boxhi[2],xy,xz,yz);
1978 }
1979 utils::logmesg(lmp,mesg);
1980 }
1981 }
1982
1983 /* ----------------------------------------------------------------------
1984 format boundary string for output
1985 assume str is 9 chars or more in length
1986 ------------------------------------------------------------------------- */
1987
boundary_string(char * str)1988 void Domain::boundary_string(char *str)
1989 {
1990 int m = 0;
1991 for (int idim = 0; idim < 3; idim++) {
1992 for (int iside = 0; iside < 2; iside++) {
1993 if (boundary[idim][iside] == 0) str[m++] = 'p';
1994 else if (boundary[idim][iside] == 1) str[m++] = 'f';
1995 else if (boundary[idim][iside] == 2) str[m++] = 's';
1996 else if (boundary[idim][iside] == 3) str[m++] = 'm';
1997 }
1998 str[m++] = ' ';
1999 }
2000 str[8] = '\0';
2001 }
2002
2003 /* ----------------------------------------------------------------------
2004 convert triclinic 0-1 lamda coords to box coords for all N atoms
2005 x = H lamda + x0;
2006 ------------------------------------------------------------------------- */
2007
lamda2x(int n)2008 void Domain::lamda2x(int n)
2009 {
2010 double **x = atom->x;
2011
2012 for (int i = 0; i < n; i++) {
2013 x[i][0] = h[0]*x[i][0] + h[5]*x[i][1] + h[4]*x[i][2] + boxlo[0];
2014 x[i][1] = h[1]*x[i][1] + h[3]*x[i][2] + boxlo[1];
2015 x[i][2] = h[2]*x[i][2] + boxlo[2];
2016 }
2017 }
2018
2019 /* ----------------------------------------------------------------------
2020 convert box coords to triclinic 0-1 lamda coords for all N atoms
2021 lamda = H^-1 (x - x0)
2022 ------------------------------------------------------------------------- */
2023
x2lamda(int n)2024 void Domain::x2lamda(int n)
2025 {
2026 double delta[3];
2027 double **x = atom->x;
2028
2029 for (int i = 0; i < n; i++) {
2030 delta[0] = x[i][0] - boxlo[0];
2031 delta[1] = x[i][1] - boxlo[1];
2032 delta[2] = x[i][2] - boxlo[2];
2033
2034 x[i][0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
2035 x[i][1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
2036 x[i][2] = h_inv[2]*delta[2];
2037 }
2038 }
2039
2040 /* ----------------------------------------------------------------------
2041 convert triclinic 0-1 lamda coords to box coords for one atom
2042 x = H lamda + x0;
2043 lamda and x can point to same 3-vector
2044 ------------------------------------------------------------------------- */
2045
lamda2x(double * lamda,double * x)2046 void Domain::lamda2x(double *lamda, double *x)
2047 {
2048 x[0] = h[0]*lamda[0] + h[5]*lamda[1] + h[4]*lamda[2] + boxlo[0];
2049 x[1] = h[1]*lamda[1] + h[3]*lamda[2] + boxlo[1];
2050 x[2] = h[2]*lamda[2] + boxlo[2];
2051 }
2052
2053 /* ----------------------------------------------------------------------
2054 convert box coords to triclinic 0-1 lamda coords for one atom
2055 lamda = H^-1 (x - x0)
2056 x and lamda can point to same 3-vector
2057 ------------------------------------------------------------------------- */
2058
x2lamda(double * x,double * lamda)2059 void Domain::x2lamda(double *x, double *lamda)
2060 {
2061 double delta[3];
2062 delta[0] = x[0] - boxlo[0];
2063 delta[1] = x[1] - boxlo[1];
2064 delta[2] = x[2] - boxlo[2];
2065
2066 lamda[0] = h_inv[0]*delta[0] + h_inv[5]*delta[1] + h_inv[4]*delta[2];
2067 lamda[1] = h_inv[1]*delta[1] + h_inv[3]*delta[2];
2068 lamda[2] = h_inv[2]*delta[2];
2069 }
2070
2071 /* ----------------------------------------------------------------------
2072 convert box coords to triclinic 0-1 lamda coords for one atom
2073 use my_boxlo & my_h_inv stored by caller for previous state of box
2074 lamda = H^-1 (x - x0)
2075 x and lamda can point to same 3-vector
2076 ------------------------------------------------------------------------- */
2077
x2lamda(double * x,double * lamda,double * my_boxlo,double * my_h_inv)2078 void Domain::x2lamda(double *x, double *lamda,
2079 double *my_boxlo, double *my_h_inv)
2080 {
2081 double delta[3];
2082 delta[0] = x[0] - my_boxlo[0];
2083 delta[1] = x[1] - my_boxlo[1];
2084 delta[2] = x[2] - my_boxlo[2];
2085
2086 lamda[0] = my_h_inv[0]*delta[0] + my_h_inv[5]*delta[1] + my_h_inv[4]*delta[2];
2087 lamda[1] = my_h_inv[1]*delta[1] + my_h_inv[3]*delta[2];
2088 lamda[2] = my_h_inv[2]*delta[2];
2089 }
2090
2091 /* ----------------------------------------------------------------------
2092 convert 8 lamda corner pts of lo/hi box to box coords
2093 return bboxlo/hi = bounding box around 8 corner pts in box coords
2094 ------------------------------------------------------------------------- */
2095
bbox(double * lo,double * hi,double * bboxlo,double * bboxhi)2096 void Domain::bbox(double *lo, double *hi, double *bboxlo, double *bboxhi)
2097 {
2098 double x[3];
2099
2100 bboxlo[0] = bboxlo[1] = bboxlo[2] = BIG;
2101 bboxhi[0] = bboxhi[1] = bboxhi[2] = -BIG;
2102
2103 x[0] = lo[0]; x[1] = lo[1]; x[2] = lo[2];
2104 lamda2x(x,x);
2105 bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
2106 bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
2107 bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
2108
2109 x[0] = hi[0]; x[1] = lo[1]; x[2] = lo[2];
2110 lamda2x(x,x);
2111 bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
2112 bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
2113 bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
2114
2115 x[0] = lo[0]; x[1] = hi[1]; x[2] = lo[2];
2116 lamda2x(x,x);
2117 bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
2118 bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
2119 bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
2120
2121 x[0] = hi[0]; x[1] = hi[1]; x[2] = lo[2];
2122 lamda2x(x,x);
2123 bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
2124 bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
2125 bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
2126
2127 x[0] = lo[0]; x[1] = lo[1]; x[2] = hi[2];
2128 lamda2x(x,x);
2129 bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
2130 bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
2131 bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
2132
2133 x[0] = hi[0]; x[1] = lo[1]; x[2] = hi[2];
2134 lamda2x(x,x);
2135 bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
2136 bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
2137 bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
2138
2139 x[0] = lo[0]; x[1] = hi[1]; x[2] = hi[2];
2140 lamda2x(x,x);
2141 bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
2142 bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
2143 bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
2144
2145 x[0] = hi[0]; x[1] = hi[1]; x[2] = hi[2];
2146 lamda2x(x,x);
2147 bboxlo[0] = MIN(bboxlo[0],x[0]); bboxhi[0] = MAX(bboxhi[0],x[0]);
2148 bboxlo[1] = MIN(bboxlo[1],x[1]); bboxhi[1] = MAX(bboxhi[1],x[1]);
2149 bboxlo[2] = MIN(bboxlo[2],x[2]); bboxhi[2] = MAX(bboxhi[2],x[2]);
2150 }
2151
2152 /* ----------------------------------------------------------------------
2153 compute 8 corner pts of my triclinic sub-box
2154 output is in corners, see ordering in lamda_box_corners
2155 ------------------------------------------------------------------------- */
2156
box_corners()2157 void Domain::box_corners()
2158 {
2159 lamda_box_corners(boxlo_lamda,boxhi_lamda);
2160 }
2161
2162 /* ----------------------------------------------------------------------
2163 compute 8 corner pts of my triclinic sub-box
2164 output is in corners, see ordering in lamda_box_corners
2165 ------------------------------------------------------------------------- */
2166
subbox_corners()2167 void Domain::subbox_corners()
2168 {
2169 lamda_box_corners(sublo_lamda,subhi_lamda);
2170 }
2171
2172 /* ----------------------------------------------------------------------
2173 compute 8 corner pts of any triclinic box with lo/hi in lamda coords
2174 8 output corners are ordered with x changing fastest, then y, finally z
2175 could be more efficient if just coded with xy,yz,xz explicitly
2176 ------------------------------------------------------------------------- */
2177
lamda_box_corners(double * lo,double * hi)2178 void Domain::lamda_box_corners(double *lo, double *hi)
2179 {
2180 corners[0][0] = lo[0]; corners[0][1] = lo[1]; corners[0][2] = lo[2];
2181 lamda2x(corners[0],corners[0]);
2182 corners[1][0] = hi[0]; corners[1][1] = lo[1]; corners[1][2] = lo[2];
2183 lamda2x(corners[1],corners[1]);
2184 corners[2][0] = lo[0]; corners[2][1] = hi[1]; corners[2][2] = lo[2];
2185 lamda2x(corners[2],corners[2]);
2186 corners[3][0] = hi[0]; corners[3][1] = hi[1]; corners[3][2] = lo[2];
2187 lamda2x(corners[3],corners[3]);
2188 corners[4][0] = lo[0]; corners[4][1] = lo[1]; corners[4][2] = hi[2];
2189 lamda2x(corners[4],corners[4]);
2190 corners[5][0] = hi[0]; corners[5][1] = lo[1]; corners[5][2] = hi[2];
2191 lamda2x(corners[5],corners[5]);
2192 corners[6][0] = lo[0]; corners[6][1] = hi[1]; corners[6][2] = hi[2];
2193 lamda2x(corners[6],corners[6]);
2194 corners[7][0] = hi[0]; corners[7][1] = hi[1]; corners[7][2] = hi[2];
2195 lamda2x(corners[7],corners[7]);
2196 }
2197