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:  Quy-Dong To, Universite Gustave Eiffel, FRANCE
17    Email: toquydong at gmail.com
18 ------------------------------------------------------------------------- */
19 
20 #include "fix_wall_reflect_stochastic.h"
21 
22 #include "atom.h"
23 #include "comm.h"
24 #include "domain.h"
25 #include "error.h"
26 #include "force.h"
27 #include "lattice.h"
28 #include "random_mars.h"
29 
30 #include <cmath>
31 #include <cstring>
32 
33 using namespace LAMMPS_NS;
34 using namespace FixConst;
35 
36 enum{NONE,DIFFUSIVE,MAXWELL,CCL};
37 
38 /* ---------------------------------------------------------------------- */
39 
40 FixWallReflectStochastic::
FixWallReflectStochastic(LAMMPS * lmp,int narg,char ** arg)41 FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) :
42   FixWallReflect(lmp, narg, arg), random(nullptr)
43 {
44   if (narg < 8) error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
45 
46   if (domain->triclinic != 0)
47     error->all(FLERR, "Fix wall/reflect/stochastic cannot be used with "
48                "triclinic simulation box");
49 
50   dynamic_group_allow = 1;
51 
52   // parse args
53 
54   int arginc;
55 
56   nwall = 0;
57   int scaleflag = 1;
58   rstyle = NONE;
59 
60   if (strcmp(arg[3],"diffusive") == 0) {
61     rstyle = DIFFUSIVE;
62     arginc = 6;
63   } else if (strcmp(arg[3],"maxwell") == 0) {
64     rstyle = MAXWELL;
65     arginc = 7;
66   } else if (strcmp(arg[3],"ccl") == 0) {
67     rstyle = CCL;
68     arginc = 9;
69   } else error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
70 
71 
72   seedfix = utils::inumeric(FLERR,arg[4],false,lmp);
73   if (seedfix <= 0) error->all(FLERR,"Random seed must be a positive number");
74 
75   int iarg = 5;
76   while (iarg < narg) {
77     if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) ||
78         (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) ||
79         (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) {
80       if (iarg+arginc > narg)
81         error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
82 
83       int newwall;
84       if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO;
85       else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI;
86       else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO;
87       else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI;
88       else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO;
89       else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI;
90 
91       for (int m = 0; (m < nwall) && (m < 6); m++)
92         if (newwall == wallwhich[m])
93           error->all(FLERR,
94                      "Fix wall/reflect/stochastic command wall defined twice");
95 
96       wallwhich[nwall] = newwall;
97       if (strcmp(arg[iarg+1],"EDGE") == 0) {
98         wallstyle[nwall] = EDGE;
99         int dim = wallwhich[nwall] / 2;
100         int side = wallwhich[nwall] % 2;
101         if (side == 0) coord0[nwall] = domain->boxlo[dim];
102         else coord0[nwall] = domain->boxhi[dim];
103       } else {
104         wallstyle[nwall] = CONSTANT;
105         coord0[nwall] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
106       }
107 
108       walltemp[nwall]= utils::numeric(FLERR,arg[iarg+2],false,lmp);
109 
110       for (int dir = 0; dir < 3; dir++) {
111         wallvel[nwall][dir]= utils::numeric(FLERR,arg[iarg+dir+3],false,lmp);
112         int dim = wallwhich[nwall] / 2;
113         if ((wallvel[nwall][dir] !=0) & (dir == dim))
114           error->all(FLERR,"The wall velocity must be tangential");
115 
116         // DIFFUSIVE = no accommodation coeffs
117         // MAXWELL = one for all dimensions
118         // CCL = one for each dimension
119 
120         if (rstyle == CCL)
121           wallaccom[nwall][dir]= utils::numeric(FLERR,arg[iarg+dir+6],false,lmp);
122         else if (rstyle == MAXWELL)
123           wallaccom[nwall][dir]= utils::numeric(FLERR,arg[iarg+6],false,lmp);
124       }
125 
126       nwall++;
127       iarg += arginc;
128 
129     } else if (strcmp(arg[iarg],"units") == 0) {
130       if (iarg+2 > narg)
131         error->all(FLERR,"Illegal wall/reflect/stochastic command");
132       if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
133       else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
134       else error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
135       iarg += 2;
136     } else error->all(FLERR,"Illegal fix wall/reflect/stochastic command");
137   }
138 
139   // error check
140 
141   if (nwall == 0) error->all(FLERR,"Illegal fix wall command");
142 
143   for (int m = 0; m < nwall; m++) {
144     if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
145       error->all(FLERR,"Cannot use fix wall/reflect/stochastic "
146                  "in periodic dimension");
147     if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
148       error->all(FLERR,"Cannot use fix wall/reflect/stochastic "
149                  "in periodic dimension");
150     if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
151       error->all(FLERR,"Cannot use fix wall/reflect/stochastic "
152                  "in periodic dimension");
153   }
154 
155   for (int m = 0; m < nwall; m++)
156     if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2)
157       error->all(FLERR,
158                  "Cannot use fix wall/reflect/stochastic zlo/zhi "
159                  "for a 2d simulation");
160 
161   // scale factors for CONSTANT walls, VARIABLE wall is not allowed
162 
163   int flag = 0;
164   for (int m = 0; m < nwall; m++)
165     if (wallstyle[m] != EDGE) flag = 1;
166 
167   if (flag) {
168     if (scaleflag) {
169       xscale = domain->lattice->xlattice;
170       yscale = domain->lattice->ylattice;
171       zscale = domain->lattice->zlattice;
172     }
173     else xscale = yscale = zscale = 1.0;
174 
175     for (int m = 0; m < nwall; m++) {
176       if (wallstyle[m] != CONSTANT) continue;
177       if (wallwhich[m] < YLO) coord0[m] *= xscale;
178       else if (wallwhich[m] < ZLO) coord0[m] *= yscale;
179       else coord0[m] *= zscale;
180     }
181   }
182 
183   // random number generator
184 
185   random = new RanMars(lmp,seedfix + comm->me);
186 }
187 
188 /* ---------------------------------------------------------------------- */
189 
~FixWallReflectStochastic()190 FixWallReflectStochastic::~FixWallReflectStochastic()
191 {
192   delete random;
193 }
194 
195 /* ---------------------------------------------------------------------- */
196 
wall_particle(int m,int which,double coord)197 void FixWallReflectStochastic::wall_particle(int m, int which, double coord)
198 {
199   int i, dir, dim, side, sign;
200   double factor,timecol,difftest,theta;
201 
202   double *rmass = atom->rmass;
203   double *mass = atom->mass;
204 
205   // coord = current position of wall
206 
207   double **x = atom->x;
208   double **v = atom->v;
209   int *mask = atom->mask;
210   int *type = atom->type;
211   int nlocal = atom->nlocal;
212 
213   dim = which / 2;
214   side = which % 2;
215 
216   for (i = 0; i < nlocal; i++) {
217     if (!(mask[i] & groupbit)) continue;
218 
219     sign = 0;
220     if ((side == 0) & (x[i][dim] < coord)) sign = 1;
221     else if ((side == 1) & (x[i][dim] > coord)) sign = -1;
222     if (sign == 0) continue;
223 
224     // theta = kT/m
225 
226     if (rmass) theta = force->boltz*walltemp[m]/(rmass[i]*force->mvv2e);
227     else theta = force->boltz*walltemp[m]/(mass[type[i]]*force->mvv2e);
228     factor = sqrt(theta);
229 
230     // time travelling after collision
231 
232     timecol = (x[i][dim]-coord) / v[i][dim];
233 
234     // only needed for Maxwell model
235 
236     if (rstyle == MAXWELL) difftest = random->uniform();
237 
238     for (dir = 0; dir < 3; dir++) {
239 
240       // pushing back atoms to wall position, assign new reflected velocity
241 
242       x[i][dir] -= v[i][dir]*timecol;
243 
244       // diffusive reflection
245 
246       if (rstyle  == DIFFUSIVE) {
247         if (dir != dim)
248           v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor);
249         else v[i][dir] =  sign*random->rayleigh(factor);
250 
251       // Maxwell reflection
252 
253       } else if (rstyle  == MAXWELL) {
254         if (difftest < wallaccom[m][dir]) {
255           if (dir != dim)
256             v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor);
257           else v[i][dir] =  sign*random->rayleigh(factor);
258         } else {
259           if (dir == dim) v[i][dir] = -v[i][dir];
260         }
261 
262       // Cercignani Lampis reflection
263 
264       } else if (rstyle  == CCL) {
265         if (dir != dim)
266           v[i][dir] = wallvel[m][dir] +
267             random->gaussian((1-wallaccom[m][dir]) *
268                              (v[i][dir] - wallvel[m][dir]),
269                              sqrt((2-wallaccom[m][dir]) *
270                                   wallaccom[m][dir]*theta));
271         else v[i][dir] = random->besselexp(theta,wallaccom[m][dir],v[i][dir]);
272       }
273 
274       // update new position due to the new velocity
275 
276       if (dir != dim) x[i][dir] += v[i][dir]*timecol;
277       else x[i][dir] = coord + v[i][dir]*timecol;
278     }
279   }
280 }
281