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 #include "compute_angle_local.h"
16 #include <cmath>
17 #include <cstring>
18 #include "atom.h"
19 #include "atom_vec.h"
20 #include "molecule.h"
21 #include "update.h"
22 #include "domain.h"
23 #include "force.h"
24 #include "angle.h"
25 #include "input.h"
26 #include "variable.h"
27 #include "math_const.h"
28 #include "memory.h"
29 #include "error.h"
30
31 using namespace LAMMPS_NS;
32 using namespace MathConst;
33
34 #define DELTA 10000
35
36 enum{THETA,ENG,VARIABLE};
37
38 /* ---------------------------------------------------------------------- */
39
ComputeAngleLocal(LAMMPS * lmp,int narg,char ** arg)40 ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
41 Compute(lmp, narg, arg),
42 bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr)
43 {
44 if (narg < 4) error->all(FLERR,"Illegal compute angle/local command");
45
46 if (atom->avec->angles_allow == 0)
47 error->all(FLERR,"Compute angle/local used when angles are not allowed");
48
49 local_flag = 1;
50
51 // style args
52
53 nvalues = narg - 3;
54 bstyle = new int[nvalues];
55 vstr = new char*[nvalues];
56 vvar = new int[nvalues];
57
58 nvalues = 0;
59 tflag = 0;
60 nvar = 0;
61
62 int iarg;
63 for (iarg = 3; iarg < narg; iarg++) {
64 if (strcmp(arg[iarg],"theta") == 0) {
65 bstyle[nvalues++] = THETA;
66 tflag = 1;
67 } else if (strcmp(arg[iarg],"eng") == 0) {
68 bstyle[nvalues++] = ENG;
69 } else if (strncmp(arg[iarg],"v_",2) == 0) {
70 bstyle[nvalues++] = VARIABLE;
71 vstr[nvar] = utils::strdup(&arg[iarg][2]);
72 nvar++;
73 } else break;
74 }
75
76 // optional args
77
78 setflag = 0;
79 tstr = nullptr;
80
81 while (iarg < narg) {
82 if (strcmp(arg[iarg],"set") == 0) {
83 setflag = 1;
84 if (iarg+3 > narg) error->all(FLERR,"Illegal compute angle/local command");
85 if (strcmp(arg[iarg+1],"theta") == 0) {
86 delete [] tstr;
87 tstr = utils::strdup(arg[iarg+2]);
88 tflag = 1;
89 } else error->all(FLERR,"Illegal compute angle/local command");
90 iarg += 3;
91 } else error->all(FLERR,"Illegal compute angle/local command");
92 }
93
94 // error check
95
96 if (nvar) {
97 if (!setflag)
98 error->all(FLERR,"Compute angle/local variable requires a set variable");
99 for (int i = 0; i < nvar; i++) {
100 vvar[i] = input->variable->find(vstr[i]);
101 if (vvar[i] < 0)
102 error->all(FLERR,"Variable name for copute angle/local does not exist");
103 if (!input->variable->equalstyle(vvar[i]))
104 error->all(FLERR,"Variable for compute angle/local is invalid style");
105 }
106
107 if (tstr) {
108 tvar = input->variable->find(tstr);
109 if (tvar < 0)
110 error->all(FLERR,"Variable name for compute angle/local does not exist");
111 if (!input->variable->internalstyle(tvar))
112 error->all(FLERR,"Variable for compute angle/local is invalid style");
113 }
114 } else if (setflag)
115 error->all(FLERR,"Compute angle/local set with no variable");
116
117 // initialize output
118
119 if (nvalues == 1) size_local_cols = 0;
120 else size_local_cols = nvalues;
121
122 nmax = 0;
123 vlocal = nullptr;
124 alocal = nullptr;
125 }
126
127 /* ---------------------------------------------------------------------- */
128
~ComputeAngleLocal()129 ComputeAngleLocal::~ComputeAngleLocal()
130 {
131 delete [] bstyle;
132 for (int i = 0; i < nvar; i++) delete [] vstr[i];
133 delete [] vstr;
134 delete [] vvar;
135
136 delete [] tstr;
137
138 memory->destroy(vlocal);
139 memory->destroy(alocal);
140 }
141
142 /* ---------------------------------------------------------------------- */
143
init()144 void ComputeAngleLocal::init()
145 {
146 if (force->angle == nullptr)
147 error->all(FLERR,"No angle style is defined for compute angle/local");
148
149 if (nvar) {
150 for (int i = 0; i < nvar; i++) {
151 vvar[i] = input->variable->find(vstr[i]);
152 if (vvar[i] < 0)
153 error->all(FLERR,"Variable name for compute angle/local does not exist");
154 }
155
156 if (tstr) {
157 tvar = input->variable->find(tstr);
158 if (tvar < 0)
159 error->all(FLERR,"Variable name for compute angle/local does not exist");
160 }
161 }
162
163 // do initial memory allocation so that memory_usage() is correct
164
165 ncount = compute_angles(0);
166 if (ncount > nmax) reallocate(ncount);
167 size_local_rows = ncount;
168 }
169
170 /* ---------------------------------------------------------------------- */
171
compute_local()172 void ComputeAngleLocal::compute_local()
173 {
174 invoked_local = update->ntimestep;
175
176 // count local entries and compute angle info
177
178 ncount = compute_angles(0);
179 if (ncount > nmax) reallocate(ncount);
180 size_local_rows = ncount;
181 ncount = compute_angles(1);
182 }
183
184 /* ----------------------------------------------------------------------
185 count angles and compute angle info on this proc
186 only count if 2nd atom is the one storing the angle
187 all atoms in interaction must be in group
188 all atoms in interaction must be known to proc
189 if angle is deleted (type = 0), do not count
190 if angle is turned off (type < 0), still count
191 if flag is set, compute requested info about angle
192 if angle is turned off (type < 0), energy = 0.0
193 ------------------------------------------------------------------------- */
194
compute_angles(int flag)195 int ComputeAngleLocal::compute_angles(int flag)
196 {
197 int i,m,na,atom1,atom2,atom3,imol,iatom,atype,ivar;
198 tagint tagprev;
199 double delx1,dely1,delz1,delx2,dely2,delz2;
200 double rsq1,rsq2,r1,r2,c,theta;
201 double *ptr;
202
203 double **x = atom->x;
204 tagint *tag = atom->tag;
205 int *num_angle = atom->num_angle;
206 tagint **angle_atom1 = atom->angle_atom1;
207 tagint **angle_atom2 = atom->angle_atom2;
208 tagint **angle_atom3 = atom->angle_atom3;
209 int **angle_type = atom->angle_type;
210 int *mask = atom->mask;
211
212 int *molindex = atom->molindex;
213 int *molatom = atom->molatom;
214 Molecule **onemols = atom->avec->onemols;
215
216 int nlocal = atom->nlocal;
217 int molecular = atom->molecular;
218
219 // loop over all atoms and their angles
220
221 Angle *angle = force->angle;
222
223 m = 0;
224 for (atom2 = 0; atom2 < nlocal; atom2++) {
225 if (!(mask[atom2] & groupbit)) continue;
226
227 if (molecular == Atom::MOLECULAR) na = num_angle[atom2];
228 else {
229 if (molindex[atom2] < 0) continue;
230 imol = molindex[atom2];
231 iatom = molatom[atom2];
232 na = onemols[imol]->num_angle[iatom];
233 }
234
235 for (i = 0; i < na; i++) {
236 if (molecular == Atom::MOLECULAR) {
237 if (tag[atom2] != angle_atom2[atom2][i]) continue;
238 atype = angle_type[atom2][i];
239 atom1 = atom->map(angle_atom1[atom2][i]);
240 atom3 = atom->map(angle_atom3[atom2][i]);
241 } else {
242 if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue;
243 atype = onemols[imol]->angle_type[atom2][i];
244 tagprev = tag[atom2] - iatom - 1;
245 atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i]+tagprev);
246 atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i]+tagprev);
247 }
248
249 if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
250 if (atom3 < 0 || !(mask[atom3] & groupbit)) continue;
251 if (atype == 0) continue;
252
253 if (!flag) {
254 m++;
255 continue;
256 }
257
258 // theta needed by one or more outputs
259
260 if (tflag) {
261 delx1 = x[atom1][0] - x[atom2][0];
262 dely1 = x[atom1][1] - x[atom2][1];
263 delz1 = x[atom1][2] - x[atom2][2];
264 domain->minimum_image(delx1,dely1,delz1);
265
266 rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
267 r1 = sqrt(rsq1);
268
269 delx2 = x[atom3][0] - x[atom2][0];
270 dely2 = x[atom3][1] - x[atom2][1];
271 delz2 = x[atom3][2] - x[atom2][2];
272 domain->minimum_image(delx2,dely2,delz2);
273
274 rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
275 r2 = sqrt(rsq2);
276
277 // c = cosine of angle
278 // theta = angle in radians
279
280 c = delx1*delx2 + dely1*dely2 + delz1*delz2;
281 c /= r1*r2;
282 if (c > 1.0) c = 1.0;
283 if (c < -1.0) c = -1.0;
284 theta = acos(c);
285 }
286
287 if (nvalues == 1) ptr = &vlocal[m];
288 else ptr = alocal[m];
289
290 if (nvar) {
291 ivar = 0;
292 if (tstr) input->variable->internal_set(tvar,theta);
293 }
294
295 for (int n = 0; n < nvalues; n++) {
296 switch (bstyle[n]) {
297 case THETA:
298 ptr[n] = 180.0*theta/MY_PI;
299 break;
300 case ENG:
301 if (atype > 0) ptr[n] = angle->single(atype,atom1,atom2,atom3);
302 else ptr[n] = 0.0;
303 break;
304 case VARIABLE:
305 ptr[n] = input->variable->compute_equal(vvar[ivar]);
306 ivar++;
307 break;
308 }
309 }
310
311 m++;
312 }
313 }
314
315 return m;
316 }
317
318 /* ---------------------------------------------------------------------- */
319
reallocate(int n)320 void ComputeAngleLocal::reallocate(int n)
321 {
322 // grow vector_local or array_local
323
324 while (nmax < n) nmax += DELTA;
325
326 if (nvalues == 1) {
327 memory->destroy(vlocal);
328 memory->create(vlocal,nmax,"angle/local:vector_local");
329 vector_local = vlocal;
330 } else {
331 memory->destroy(alocal);
332 memory->create(alocal,nmax,nvalues,"angle/local:array_local");
333 array_local = alocal;
334 }
335 }
336
337 /* ----------------------------------------------------------------------
338 memory usage of local data
339 ------------------------------------------------------------------------- */
340
memory_usage()341 double ComputeAngleLocal::memory_usage()
342 {
343 double bytes = (double)nmax*nvalues * sizeof(double);
344 return bytes;
345 }
346