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: Abdoreza Ershadinia, a.ershadinia at gmail.com
17 ------------------------------------------------------------------------- */
18
19 #include "fix_wall_ees.h"
20
21 #include "math_extra.h"
22 #include "atom.h"
23 #include "atom_vec_ellipsoid.h"
24 #include "error.h"
25
26 #include <cmath>
27
28 using namespace LAMMPS_NS;
29 using namespace FixConst;
30
31 /* ---------------------------------------------------------------------- */
32
FixWallEES(LAMMPS * lmp,int narg,char ** arg)33 FixWallEES::FixWallEES(LAMMPS *lmp, int narg, char **arg) :
34 FixWall(lmp, narg, arg) {}
35
36 /* ---------------------------------------------------------------------- */
37
precompute(int m)38 void FixWallEES::precompute(int m)
39 {
40 coeff1[m] = ( 2. / 4725. ) * epsilon[m] * pow(sigma[m],12.0);
41 coeff2[m] = ( 1. / 24. ) * epsilon[m] * pow(sigma[m],6.0);
42
43 coeff3[m] = ( 2. / 315. ) * epsilon[m] * pow(sigma[m],12.0);
44 coeff4[m] = ( 1. / 3. ) * epsilon[m] * pow(sigma[m],6.0);
45
46 coeff5[m] = ( 4. / 315. ) * epsilon[m] * pow(sigma[m],12.0);
47 coeff6[m] = ( 1. / 12. ) * epsilon[m] * pow(sigma[m],6.0);
48 }
49
50 /* ---------------------------------------------------------------------- */
init()51 void FixWallEES::init()
52 {
53 avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
54 if (!avec)
55 error->all(FLERR,"Fix wall/ees requires atom style ellipsoid");
56
57 // check that all particles are finite-size ellipsoids
58 // no point particles allowed, spherical is OK
59
60 int *ellipsoid = atom->ellipsoid;
61 int *mask = atom->mask;
62 int nlocal = atom->nlocal;
63
64 for (int i = 0; i < nlocal; i++)
65 if (mask[i] & groupbit)
66 if (ellipsoid[i] < 0)
67 error->one(FLERR,"Fix wall/ees requires extended particles");
68
69 FixWall::init();
70 }
71
72
73 /* ----------------------------------------------------------------------
74 interaction of all particles in group with a wall
75 m = index of wall coeffs
76 which = xlo,xhi,ylo,yhi,zlo,zhi
77 error if any particle is on or behind wall
78 ------------------------------------------------------------------------- */
79
wall_particle(int m,int which,double coord)80 void FixWallEES::wall_particle(int m, int which, double coord)
81 {
82 double delta;
83
84 double **x = atom->x;
85 double **f = atom->f;
86 double **tor = atom->torque;
87
88 avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
89 AtomVecEllipsoid::Bonus *bonus = avec->bonus;
90 int *ellipsoid = atom->ellipsoid;
91 int *mask = atom->mask;
92 int nlocal = atom->nlocal;
93
94 int dim = which / 2;
95 int side = which % 2;
96 if (side == 0) side = -1;
97
98 int onflag = 0;
99
100 for (int i = 0; i < nlocal; i++)
101 if (mask[i] & groupbit) {
102
103 if (side < 0)
104 delta = x[i][dim] - coord;
105 else
106 delta = coord - x[i][dim];
107
108 if (delta >= cutoff[m])
109 continue;
110
111 double A[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
112 double tempvec[3]= {0,0,0};
113 double sigman = 0.0, sigman2 = 0.0;
114 double nhat[3] = {0,0,0};
115
116 nhat[dim]=-1*side;
117 nhat[(dim+1)%3] = 0 ;
118 nhat[(dim+2)%3] = 0 ;
119
120
121 double* shape = bonus[ellipsoid[i]].shape;;
122 MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A);
123 MathExtra::transpose_matvec(A,nhat,tempvec);
124 for (int k = 0; k<3; k++) tempvec[k] *= shape[k];
125 for (int k = 0; k<3 ; k++) sigman2 += tempvec[k]*tempvec[k];
126 sigman = sqrt(sigman2);
127
128 if (delta <= sigman) {
129 onflag = 1;
130 continue;
131 }
132
133
134 double fwall = 0.0, twall = 0.0;
135 double delta2 = 0.0, delta3 = 0.0, delta4 = 0.0, delta5 = 0.0, delta6 = 0.0;
136 double sigman3 = 0.0, sigman4 = 0.0, sigman5 = 0.0, sigman6 = 0.0;
137 double hhss = 0.0, hhss2 = 0.0, hhss4 = 0.0, hhss7 = 0.0, hhss8 = 0.0;
138 double hps = 0.0;
139 double hms = 0.0;
140
141 double tempvec2[3]= {0,0,0};
142
143 double SAn[3] = {0,0,0};
144 double that[3] = {0,0,0};
145
146 double Lx[3][3] = {{0,0,0},{0,0,-1},{0,1,0}};
147 double Ly[3][3] = {{0,0,1},{0,0,0},{-1,0,0}};
148 double Lz[3][3] = {{0,-1,0},{1,0,0},{0,0,0}};
149
150
151 for (int k = 0; k<3; k++) SAn[k] = tempvec[k];
152
153 sigman3 = sigman2 * sigman;
154 sigman4 = sigman2 * sigman2;
155 sigman5 = sigman4 * sigman;
156 sigman6 = sigman3 * sigman3;
157
158
159 delta2 = delta * delta;
160 delta3 = delta2 * delta;
161 delta4 = delta2 * delta2;
162 delta5 = delta3 * delta2;
163 delta6 = delta3 * delta3;
164
165 hhss = delta2 - sigman2;
166 hhss2 = hhss * hhss;
167 hhss4 = hhss2 * hhss2;
168 hhss8 = hhss4 * hhss4;
169 hhss7 = hhss4 * hhss2 * hhss;
170
171 hps = delta + sigman;
172 hms = delta - sigman;
173
174 fwall = side*(
175 -1*coeff4[m]/hhss2 +
176 coeff3[m] * ( 21*delta6 + 63*delta4*sigman2 + 27*delta2*sigman4 + sigman6 ) / hhss8
177 );
178 f[i][dim] -= fwall;
179
180 ewall[0] += -1*coeff2[m] * ( 4*delta/sigman2/hhss + 2*log(hms/hps)/sigman3 ) +
181 coeff1[m] * ( 35*delta5 + 70*delta3*sigman2 + 15*delta*sigman4 ) / hhss7;
182
183 ewall[m+1] += fwall;
184
185 twall = coeff6[m] * ( 6.*delta3/sigman4/hhss2 - 10.*delta/sigman2/hhss2 + 3.*log(hms/hps)/sigman5 ) +
186 coeff5[m] * ( 21.*delta5 + 30.*delta3*sigman2 + 5.*delta*sigman4 ) / hhss8 ;
187
188
189 MathExtra::matvec(Lx,nhat,tempvec);
190 MathExtra::transpose_matvec(A,tempvec,tempvec2);
191 for (int k = 0; k<3; k++) tempvec2[k] *= shape[k];
192 that[0] = MathExtra::dot3(SAn,tempvec2);
193
194 MathExtra::matvec(Ly,nhat,tempvec);
195 MathExtra::transpose_matvec(A,tempvec,tempvec2);
196 for (int k = 0; k<3; k++) tempvec2[k] *= shape[k];
197 that[1] = MathExtra::dot3(SAn,tempvec2);
198
199 MathExtra::matvec(Lz,nhat,tempvec);
200 MathExtra::transpose_matvec(A,tempvec,tempvec2);
201 for (int k = 0; k < 3; k++) tempvec2[k] *= shape[k];
202 that[2] = MathExtra::dot3(SAn,tempvec2);
203
204
205 for (int j = 0; j<3 ; j++)
206 tor[i][j] += twall * that[j];
207
208 }
209
210 if (onflag) error->one(FLERR,"Particle on or inside fix wall surface");
211 }
212