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