1 /* ----------------------------------------------------------------------
2     This is the
3 
4     ██╗     ██╗ ██████╗  ██████╗  ██████╗ ██╗  ██╗████████╗███████╗
5     ██║     ██║██╔════╝ ██╔════╝ ██╔════╝ ██║  ██║╚══██╔══╝██╔════╝
6     ██║     ██║██║  ███╗██║  ███╗██║  ███╗███████║   ██║   ███████╗
7     ██║     ██║██║   ██║██║   ██║██║   ██║██╔══██║   ██║   ╚════██║
8     ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║  ██║   ██║   ███████║
9     ╚══════╝╚═╝ ╚═════╝  ╚═════╝  ╚═════╝ ╚═╝  ╚═╝   ╚═╝   ╚══════╝®
10 
11     DEM simulation engine, released by
12     DCS Computing Gmbh, Linz, Austria
13     http://www.dcs-computing.com, office@dcs-computing.com
14 
15     LIGGGHTS® is part of CFDEM®project:
16     http://www.liggghts.com | http://www.cfdem.com
17 
18     Core developer and main author:
19     Christoph Kloss, christoph.kloss@dcs-computing.com
20 
21     LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22     License, version 2 or later. It is distributed in the hope that it will
23     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25     received a copy of the GNU General Public License along with LIGGGHTS®.
26     If not, see http://www.gnu.org/licenses . See also top-level README
27     and LICENSE files.
28 
29     LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30     the producer of the LIGGGHTS® software and the CFDEM®coupling software
31     See http://www.cfdem.com/terms-trademark-policy for details.
32 
33 -------------------------------------------------------------------------
34     Contributing author and copyright for this file:
35     This file is from LAMMPS, but has been modified. Copyright for
36     modification:
37 
38     Copyright 2012-     DCS Computing GmbH, Linz
39     Copyright 2009-2012 JKU Linz
40 
41     Copyright of original file:
42     LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
43     http://lammps.sandia.gov, Sandia National Laboratories
44     Steve Plimpton, sjplimp@sandia.gov
45 
46     Copyright (2003) Sandia Corporation.  Under the terms of Contract
47     DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
48     certain rights in this software.  This software is distributed under
49     the GNU General Public License.
50 ------------------------------------------------------------------------- */
51 
52 #ifndef LMP_DOMAIN_H
53 #define LMP_DOMAIN_H
54 
55 #include <cmath>
56 #include "pointers.h"
57 #include "error.h"
58 #include "comm.h"
59 #include "vector_liggghts.h"
60 #include "neighbor.h"
61 #include "atom.h"
62 #include "math_extra_liggghts.h"
63 
64 #define SMALL_DMBRDR 1.0e-8
65 
66 namespace LAMMPS_NS {
67 
68 class Domain : protected Pointers {
69  public:
70   int box_exist;                         // 0 = not yet created, 1 = exists
71   int dimension;                         // 2 = 2d, 3 = 3d
72   int nonperiodic;                       // 0 = periodic in all 3 dims
73                                          // 1 = periodic or fixed in all 6
74                                          // 2 = shrink-wrap in any of 6
75   int xperiodic,yperiodic,zperiodic;     // 0 = non-periodic, 1 = periodic
76   int periodicity[3];                    // xyz periodicity as array
77 
78   int boundary[3][2];                    // settings for 6 boundaries
79                                          // 0 = periodic
80                                          // 1 = fixed non-periodic
81                                          // 2 = shrink-wrap non-periodic
82                                          // 3 = shrink-wrap non-per w/ min
83 
84   int triclinic;                         // 0 = orthog box, 1 = triclinic
85   int tiltsmall;                         // 1 if limit tilt, else 0
86 
87                                          // orthogonal box
88   double xprd,yprd,zprd;                 // global box dimensions
89   double xprd_half,yprd_half,zprd_half;  // half dimensions
90   double prd[3];                         // array form of dimensions
91   double prd_half[3];                    // array form of half dimensions
92 
93                                          // triclinic box
94                                          // xprd,xprd_half,prd,prd_half =
95                                          // same as if untilted
96   double prd_lamda[3];                   // lamda box = (1,1,1)
97   double prd_half_lamda[3];              // lamda half box = (0.5,0.5,0.5)
98 
99   double boxlo[3],boxhi[3];              // orthogonal box global bounds
100 
101                                          // triclinic box
102                                          // boxlo/hi = same as if untilted
103   double boxlo_lamda[3],boxhi_lamda[3];  // lamda box = (0,1)
104   double boxlo_bound[3],boxhi_bound[3];  // bounding box of tilted domain
105   double corners[8][3];                  // 8 corner points
106 
107                                          // orthogonal box & triclinic box
108   double minxlo,minxhi;                  // minimum size of global box
109   double minylo,minyhi;                  // when shrink-wrapping
110   double minzlo,minzhi;                  // tri only possible for non-skew dims
111 
112                                          // orthogonal box
113   double sublo[3],subhi[3];              // sub-box bounds on this proc
114 
115                                          // triclinic box
116                                          // sublo/hi = undefined
117   double sublo_lamda[3],subhi_lamda[3];  // bounds of subbox in lamda
118 
119                                          // triclinic box
120   double xy,xz,yz;                       // 3 tilt factors
121   double h[6],h_inv[6];                           // shape matrix in Voigt notation
122   double h_rate[6],h_ratelo[3];          // rate of box size/shape change
123 
124   int box_change;                // 1 if any of next 3 flags are set, else 0
125   int box_change_size;           // 1 if box size changes, 0 if not
126   int box_change_shape;          // 1 if box shape changes, 0 if not
127 
128   int box_change_domain;         // 1 if proc sub-domains change, 0 if not
129 
130   int deform_flag;                // 1 if fix deform exist, else 0
131   int deform_vremap;              // 1 if fix deform remaps v, else 0
132   int deform_groupbit;            // atom group to perform v remap for
133 
134   class Lattice *lattice;                  // user-defined lattice
135 
136   int nregion;                             // # of defined Regions
137   int maxregion;                           // max # list can hold
138   class Region **regions;                  // list of defined Regions
139 
140   bool is_wedge;
141 
142   Domain(class LAMMPS *);
143   virtual ~Domain();
144   virtual void init();
145   void set_initial_box();
146   virtual void set_global_box();
147   virtual void set_lamda_box();
148   virtual void set_local_box();
149   virtual void reset_box();
150   virtual void pbc();
151   void image_check();
152   void box_too_small_check();
153   void minimum_image(double &, double &, double &);
154   void minimum_image(double *);
155   int closest_image(int, int);
156   void closest_image(const double * const, const double * const,
157                      double * const);
158   void remap(double *, tagint &);
159   void remap(double *);
160   void remap_near(double *, double *);
161   void unmap(double *, tagint);
162   void unmap(double *, tagint, double *);
163   void image_flip(int, int, int);
164   void set_lattice(int, char **);
165   void add_region(int, char **);
166   void delete_region(int, char **);
167   int find_region(const char *);
168   virtual void set_boundary(int, char **, int);
169   void set_box(int, char **);
170   virtual void print_box(const char *);
171   void boundary_string(char *);
172 
173   virtual void lamda2x(int);
174   virtual void x2lamda(int);
175   virtual void lamda2x(double *, double *);
176   virtual void x2lamda(double *, double *);
177   void x2lamda(double *, double *, double *, double *);
178   void bbox(double *, double *, double *, double *);
179   void box_corners();
180 
181   // minimum image convention check
182   // return 1 if any distance > 1/2 of box size
183   // inline since called from neighbor build inner loop
184 
minimum_image_check(double dx,double dy,double dz)185   inline int minimum_image_check(double dx, double dy, double dz) {
186     if (xperiodic && fabs(dx) > xprd_half) return 1;
187     if (yperiodic && fabs(dy) > yprd_half) return 1;
188     if (zperiodic && fabs(dz) > zprd_half) return 1;
189     return 0;
190   }
191 
192   int is_in_domain(double* pos);
193   int is_in_subdomain(double* pos);
194   int is_in_extended_subdomain(double* pos);
195   double dist_subbox_borders(double* pos);
196   void min_subbox_extent(double &min_extent,int &dim);
197   int is_periodic_ghost(int i);
198   bool is_owned_or_first_ghost(int i);
199 
is_in_domain_wedge(double * pos)200   virtual int is_in_domain_wedge(double* pos) { UNUSED(pos); return 0; }
is_in_subdomain_wedge(double * pos)201   virtual int is_in_subdomain_wedge(double* pos) { UNUSED(pos); return 0; }
is_in_extended_subdomain_wedge(double * pos)202   virtual int is_in_extended_subdomain_wedge(double* pos) { UNUSED(pos); return 0; }
dist_subbox_borders_wedge(double * pos)203   virtual double dist_subbox_borders_wedge(double* pos) { UNUSED(pos); return 0.; }
is_periodic_ghost_wedge(int i)204   virtual int is_periodic_ghost_wedge(int i) { UNUSED(i); return 0;}
205 
206  private:
207   double small[3];                  // fractions of box lengths
208 };
209 
210 #include "domain_I.h"
211 
212 }
213 
214 #endif
215 
216 /* ERROR/WARNING messages:
217 
218 E: Box bounds are invalid
219 
220 The box boundaries specified in the read_data file are invalid.  The
221 lo value must be less than the hi value for all 3 dimensions.
222 
223 E: Cannot skew triclinic box in z for 2d simulation
224 
225 Self-explanatory.
226 
227 E: Triclinic box skew is too large
228 
229 The displacement in a skewed direction must be less than half the box
230 length in that dimension.  E.g. the xy tilt must be between -half and
231 +half of the x box length.  This constraint can be relaxed by using
232 the box tilt command.
233 
234 W: Triclinic box skew is large
235 
236 The displacement in a skewed direction is normally required to be less
237 than half the box length in that dimension.  E.g. the xy tilt must be
238 between -half and +half of the x box length.  You have relaxed the
239 constraint using the box tilt command, but the warning means that a
240 LAMMPS simulation may be inefficient as a result.
241 
242 E: Illegal simulation box
243 
244 The lower bound of the simulation box is greater than the upper bound.
245 
246 E: Bond atom missing in image check
247 
248 The 2nd atom in a particular bond is missing on this processor.
249 Typically this is because the pairwise cutoff is set too short or the
250 bond has blown apart and an atom is too far away.
251 
252 W: Inconsistent image flags
253 
254 The image flags for a pair on bonded atoms appear to be inconsistent.
255 Inconsistent means that when the coordinates of the two atoms are
256 unwrapped using the image flags, the two atoms are far apart.
257 Specifically they are further apart than half a periodic box length.
258 Or they are more than a box length apart in a non-periodic dimension.
259 This is usually due to the initial data file not having correct image
260 flags for the 2 atoms in a bond that straddles a periodic boundary.
261 They should be different by 1 in that case.  This is a warning because
262 inconsistent image flags will not cause problems for dynamics or most
263 LAMMPS simulations.  However they can cause problems when such atoms
264 are used with the fix rigid or replicate commands.
265 
266 E: Bond atom missing in box size check
267 
268 The 2nd atoms needed to compute a particular bond is missing on this
269 processor.  Typically this is because the pairwise cutoff is set too
270 short or the bond has blown apart and an atom is too far away.
271 
272 W: Bond/angle/dihedral extent > half of periodic box length
273 
274 This is a restriction because LAMMPS can be confused about which image
275 of an atom in the bonded interaction is the correct one to use.
276 "Extent" in this context means the maximum end-to-end length of the
277 bond/angle/dihedral.  LAMMPS computes this by taking the maximum bond
278 length, multiplying by the number of bonds in the interaction (e.g. 3
279 for a dihedral) and adding a small amount of stretch.
280 
281 E: Illegal ... command
282 
283 Self-explanatory.  Check the input script syntax and compare to the
284 documentation for the command.  You can use -echo screen as a
285 command-line option when running LAMMPS to see the offending line.
286 
287 E: Reuse of region ID
288 
289 A region ID cannot be used twice.
290 
291 E: Invalid region style
292 
293 The choice of region style is unknown.
294 
295 E: Delete region ID does not exist
296 
297 Self-explanatory.
298 
299 E: Both sides of boundary must be periodic
300 
301 Cannot specify a boundary as periodic only on the lo or hi side.  Must
302 be periodic on both sides.
303 
304 */
305