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 "variable.h"
16
17 #include "atom.h"
18 #include "comm.h"
19 #include "compute.h"
20 #include "domain.h"
21 #include "error.h"
22 #include "fix.h"
23 #include "fix_store.h"
24 #include "group.h"
25 #include "info.h"
26 #include "input.h"
27 #include "lmppython.h"
28 #include "math_const.h"
29 #include "memory.h"
30 #include "modify.h"
31 #include "output.h"
32 #include "random_mars.h"
33 #include "region.h"
34 #include "thermo.h"
35 #include "tokenizer.h"
36 #include "universe.h"
37 #include "update.h"
38
39 #include <cctype>
40 #include <cmath>
41 #include <cstring>
42 #include <unistd.h>
43 #include <unordered_map>
44
45 using namespace LAMMPS_NS;
46 using namespace MathConst;
47
48 #define VARDELTA 4
49 #define MAXLEVEL 4
50 #define MAXLINE 256
51 #define CHUNK 1024
52 #define MAXFUNCARG 6
53
54 #define MYROUND(a) (( a-floor(a) ) >= .5) ? ceil(a) : floor(a)
55
56 enum{ARG,OP};
57
58 // customize by adding a function
59 // if add before XOR:
60 // also set precedence level in constructor and precedence length in *.h
61
62 enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,MODULO,UNARY,
63 NOT,EQ,NE,LT,LE,GT,GE,AND,OR,XOR,
64 SQRT,EXP,LN,LOG,ABS,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
65 RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,LOGFREQ2,
66 LOGFREQ3,STRIDE,STRIDE2,VDISPLACE,SWIGGLE,CWIGGLE,GMASK,RMASK,
67 GRMASK,IS_ACTIVE,IS_DEFINED,IS_AVAILABLE,IS_FILE,
68 VALUE,ATOMARRAY,TYPEARRAY,INTARRAY,BIGINTARRAY,VECTORARRAY};
69
70 // customize by adding a special function
71
72 enum{SUM,XMIN,XMAX,AVE,TRAP,SLOPE};
73
74
75 #define BIG 1.0e20
76
77 // constants for variable expressions. customize by adding new items.
78 // if needed (cf. 'version') initialize in Variable class constructor.
79
80 static std::unordered_map<std::string, double> constants = {
81 {"PI", MY_PI },
82 {"version", -1 },
83 {"yes", 1 },
84 {"no", 0 },
85 {"on", 1 },
86 {"off", 0 },
87 {"true", 1 },
88 {"false", 0 }
89 };
90
91 /* ---------------------------------------------------------------------- */
92
Variable(LAMMPS * lmp)93 Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
94 {
95 MPI_Comm_rank(world,&me);
96
97 nvar = maxvar = 0;
98 names = nullptr;
99 style = nullptr;
100 num = nullptr;
101 which = nullptr;
102 pad = nullptr;
103 reader = nullptr;
104 data = nullptr;
105 dvalue = nullptr;
106 vecs = nullptr;
107
108 eval_in_progress = nullptr;
109
110 randomequal = nullptr;
111 randomatom = nullptr;
112
113 // override initializer since LAMMPS class needs to be instantiated
114
115 constants["version"] = lmp->num_ver;
116
117 // customize by assigning a precedence level
118
119 precedence[DONE] = 0;
120 precedence[OR] = precedence[XOR] = 1;
121 precedence[AND] = 2;
122 precedence[EQ] = precedence[NE] = 3;
123 precedence[LT] = precedence[LE] = precedence[GT] = precedence[GE] = 4;
124 precedence[ADD] = precedence[SUBTRACT] = 5;
125 precedence[MULTIPLY] = precedence[DIVIDE] = precedence[MODULO] = 6;
126 precedence[CARAT] = 7;
127 precedence[UNARY] = precedence[NOT] = 8;
128 }
129
130 /* ---------------------------------------------------------------------- */
131
~Variable()132 Variable::~Variable()
133 {
134 for (int i = 0; i < nvar; i++) {
135 delete [] names[i];
136 delete reader[i];
137 if (style[i] == LOOP || style[i] == ULOOP) delete [] data[i][0];
138 else for (int j = 0; j < num[i]; j++) delete [] data[i][j];
139 delete [] data[i];
140 if (style[i] == VECTOR) memory->destroy(vecs[i].values);
141 }
142 memory->sfree(names);
143 memory->destroy(style);
144 memory->destroy(num);
145 memory->destroy(which);
146 memory->destroy(pad);
147 memory->sfree(reader);
148 memory->sfree(data);
149 memory->sfree(dvalue);
150 memory->sfree(vecs);
151
152 memory->destroy(eval_in_progress);
153
154 delete randomequal;
155 delete randomatom;
156
157 }
158
159 /* ----------------------------------------------------------------------
160 called by variable command in input script
161 ------------------------------------------------------------------------- */
162
set(int narg,char ** arg)163 void Variable::set(int narg, char **arg)
164 {
165 if (narg < 2) error->all(FLERR,"Illegal variable command");
166
167 int replaceflag = 0;
168
169 // DELETE
170 // doesn't matter if variable no longer exists
171
172 if (strcmp(arg[1],"delete") == 0) {
173 if (narg != 2) error->all(FLERR,"Illegal variable command");
174 if (find(arg[0]) >= 0) remove(find(arg[0]));
175 return;
176
177 // INDEX
178 // num = listed args, which = 1st value, data = copied args
179
180 } else if (strcmp(arg[1],"index") == 0) {
181 if (narg < 3) error->all(FLERR,"Illegal variable command");
182 if (find(arg[0]) >= 0) return;
183 if (nvar == maxvar) grow();
184 style[nvar] = INDEX;
185 num[nvar] = narg - 2;
186 which[nvar] = 0;
187 pad[nvar] = 0;
188 data[nvar] = new char*[num[nvar]];
189 copy(num[nvar],&arg[2],data[nvar]);
190
191 // LOOP
192 // 1 arg + pad: num = N, which = 1st value, data = single string
193 // 2 args + pad: num = N2, which = N1, data = single string
194
195 } else if (strcmp(arg[1],"loop") == 0) {
196 if (find(arg[0]) >= 0) return;
197 if (nvar == maxvar) grow();
198 style[nvar] = LOOP;
199 int nfirst = 0,nlast = 0;
200 if (narg == 3 || (narg == 4 && strcmp(arg[3],"pad") == 0)) {
201 nfirst = 1;
202 nlast = utils::inumeric(FLERR,arg[2],false,lmp);
203 if (nlast <= 0) error->all(FLERR,"Illegal variable command");
204 if (narg == 4 && strcmp(arg[3],"pad") == 0) {
205 pad[nvar] = fmt::format("{}",nlast).size();
206 } else pad[nvar] = 0;
207 } else if (narg == 4 || (narg == 5 && strcmp(arg[4],"pad") == 0)) {
208 nfirst = utils::inumeric(FLERR,arg[2],false,lmp);
209 nlast = utils::inumeric(FLERR,arg[3],false,lmp);
210 if (nfirst > nlast || nlast < 0)
211 error->all(FLERR,"Illegal variable command");
212 if (narg == 5 && strcmp(arg[4],"pad") == 0) {
213 pad[nvar] = fmt::format("{}",nlast).size();
214 } else pad[nvar] = 0;
215 } else error->all(FLERR,"Illegal variable command");
216 num[nvar] = nlast;
217 which[nvar] = nfirst-1;
218 data[nvar] = new char*[1];
219 data[nvar][0] = nullptr;
220
221 // WORLD
222 // num = listed args, which = partition this proc is in, data = copied args
223 // error check that num = # of worlds in universe
224
225 } else if (strcmp(arg[1],"world") == 0) {
226 if (narg < 3) error->all(FLERR,"Illegal variable command");
227 if (find(arg[0]) >= 0) return;
228 if (nvar == maxvar) grow();
229 style[nvar] = WORLD;
230 num[nvar] = narg - 2;
231 if (num[nvar] != universe->nworlds)
232 error->all(FLERR,"World variable count doesn't match # of partitions");
233 which[nvar] = universe->iworld;
234 pad[nvar] = 0;
235 data[nvar] = new char*[num[nvar]];
236 copy(num[nvar],&arg[2],data[nvar]);
237
238 // UNIVERSE and ULOOP
239 // for UNIVERSE: num = listed args, data = copied args
240 // for ULOOP: num = N, data = single string
241 // which = partition this proc is in
242 // universe proc 0 creates lock file
243 // error check that all other universe/uloop variables are same length
244
245 } else if (strcmp(arg[1],"universe") == 0 || strcmp(arg[1],"uloop") == 0) {
246 if (strcmp(arg[1],"universe") == 0) {
247 if (narg < 3) error->all(FLERR,"Illegal variable command");
248 if (find(arg[0]) >= 0) return;
249 if (nvar == maxvar) grow();
250 style[nvar] = UNIVERSE;
251 num[nvar] = narg - 2;
252 pad[nvar] = 0;
253 data[nvar] = new char*[num[nvar]];
254 copy(num[nvar],&arg[2],data[nvar]);
255 } else if (strcmp(arg[1],"uloop") == 0) {
256 if (narg < 3 || narg > 4 || (narg == 4 && strcmp(arg[3],"pad") != 0))
257 error->all(FLERR,"Illegal variable command");
258 if (find(arg[0]) >= 0) return;
259 if (nvar == maxvar) grow();
260 style[nvar] = ULOOP;
261 num[nvar] = utils::inumeric(FLERR,arg[2],false,lmp);
262 data[nvar] = new char*[1];
263 data[nvar][0] = nullptr;
264 if (narg == 4) {
265 char digits[12];
266 sprintf(digits,"%d",num[nvar]);
267 pad[nvar] = strlen(digits);
268 } else pad[nvar] = 0;
269 }
270
271 if (num[nvar] < universe->nworlds)
272 error->all(FLERR,"Universe/uloop variable count < # of partitions");
273 which[nvar] = universe->iworld;
274
275 if (universe->me == 0) {
276 FILE *fp = fopen("tmp.lammps.variable","w");
277 if (fp == nullptr)
278 error->one(FLERR,"Cannot open temporary file for world counter: "
279 + utils::getsyserror());
280 fprintf(fp,"%d\n",universe->nworlds);
281 fclose(fp);
282 fp = nullptr;
283 }
284
285 for (int jvar = 0; jvar < nvar; jvar++)
286 if (num[jvar] && (style[jvar] == UNIVERSE || style[jvar] == ULOOP) &&
287 num[nvar] != num[jvar])
288 error->all(FLERR,
289 "All universe/uloop variables must have same # of values");
290
291 // STRING
292 // replace pre-existing var if also style STRING (allows it to be reset)
293 // num = 1, which = 1st value
294 // data = 1 value, string to eval
295
296 } else if (strcmp(arg[1],"string") == 0) {
297 if (narg != 3) error->all(FLERR,"Illegal variable command");
298
299 int maxcopy = strlen(arg[2]) + 1;
300 int maxwork = maxcopy;
301 char *scopy = (char *) memory->smalloc(maxcopy,"var:string/copy");
302 char *work = (char *) memory->smalloc(maxwork,"var:string/work");
303 strcpy(scopy,arg[2]);
304 input->substitute(scopy,work,maxcopy,maxwork,1);
305 memory->sfree(work);
306
307 int ivar = find(arg[0]);
308 if (ivar >= 0) {
309 if (style[ivar] != STRING)
310 error->all(FLERR,"Cannot redefine variable as a different style");
311 delete [] data[ivar][0];
312 copy(1,&scopy,data[ivar]);
313 replaceflag = 1;
314 } else {
315 if (nvar == maxvar) grow();
316 style[nvar] = STRING;
317 num[nvar] = 1;
318 which[nvar] = 0;
319 pad[nvar] = 0;
320 data[nvar] = new char*[num[nvar]];
321 copy(1,&scopy,data[nvar]);
322 }
323 memory->sfree(scopy);
324
325 // GETENV
326 // remove pre-existing var if also style GETENV (allows it to be reset)
327 // num = 1, which = 1st value
328 // data = 1 value, string to eval
329
330 } else if (strcmp(arg[1],"getenv") == 0) {
331 if (narg != 3) error->all(FLERR,"Illegal variable command");
332 if (find(arg[0]) >= 0) {
333 if (style[find(arg[0])] != GETENV)
334 error->all(FLERR,"Cannot redefine variable as a different style");
335 remove(find(arg[0]));
336 }
337 if (nvar == maxvar) grow();
338 style[nvar] = GETENV;
339 num[nvar] = 2;
340 which[nvar] = 0;
341 pad[nvar] = 0;
342 data[nvar] = new char*[num[nvar]];
343 data[nvar][0] = utils::strdup(arg[2]);
344 data[nvar][1] = utils::strdup("(undefined)");
345
346 // SCALARFILE for strings or numbers
347 // which = 1st value
348 // data = 1 value, string to eval
349
350 } else if (strcmp(arg[1],"file") == 0) {
351 if (narg != 3) error->all(FLERR,"Illegal variable command");
352 if (find(arg[0]) >= 0) return;
353 if (nvar == maxvar) grow();
354 style[nvar] = SCALARFILE;
355 num[nvar] = 1;
356 which[nvar] = 0;
357 pad[nvar] = 0;
358 data[nvar] = new char*[num[nvar]];
359 data[nvar][0] = new char[MAXLINE];
360 reader[nvar] = new VarReader(lmp,arg[0],arg[2],SCALARFILE);
361 int flag = reader[nvar]->read_scalar(data[nvar][0]);
362 if (flag) error->all(FLERR,"File variable could not read value");
363
364 // ATOMFILE for numbers
365 // which = 1st value
366 // data = nullptr
367
368 } else if (strcmp(arg[1],"atomfile") == 0) {
369 if (narg != 3) error->all(FLERR,"Illegal variable command");
370 if (find(arg[0]) >= 0) return;
371 if (nvar == maxvar) grow();
372 style[nvar] = ATOMFILE;
373 num[nvar] = 1;
374 which[nvar] = 0;
375 pad[nvar] = 0;
376 data[nvar] = new char*[num[nvar]];
377 data[nvar][0] = nullptr;
378 reader[nvar] = new VarReader(lmp,arg[0],arg[2],ATOMFILE);
379 int flag = reader[nvar]->read_peratom();
380 if (flag) error->all(FLERR,"Atomfile variable could not read values");
381
382 // FORMAT
383 // num = 3, which = 1st value
384 // data = 3 values
385 // 1st is name of variable to eval, 2nd is format string,
386 // 3rd is filled on retrieval
387
388 } else if (strcmp(arg[1],"format") == 0) {
389 if (narg != 4) error->all(FLERR,"Illegal variable command");
390 if (find(arg[0]) >= 0) return;
391 if (nvar == maxvar) grow();
392 style[nvar] = FORMAT;
393 num[nvar] = 3;
394 which[nvar] = 0;
395 pad[nvar] = 0;
396 if (!utils::strmatch(arg[3],"%[0-9 ]*\\.[0-9]+[efgEFG]"))
397 error->all(FLERR,"Incorrect conversion in format string");
398 data[nvar] = new char*[num[nvar]];
399 copy(2,&arg[2],data[nvar]);
400 data[nvar][2] = new char[VALUELENGTH];
401 strcpy(data[nvar][2],"(undefined)");
402
403 // EQUAL
404 // replace pre-existing var if also style EQUAL (allows it to be reset)
405 // num = 2, which = 1st value
406 // data = 2 values, 1st is string to eval, 2nd is filled on retrieval
407
408 } else if (strcmp(arg[1],"equal") == 0) {
409 if (narg != 3) error->all(FLERR,"Illegal variable command");
410 int ivar = find(arg[0]);
411 if (ivar >= 0) {
412 if (style[ivar] != EQUAL)
413 error->all(FLERR,"Cannot redefine variable as a different style");
414 delete [] data[ivar][0];
415 data[ivar][0] = utils::strdup(arg[2]);
416 replaceflag = 1;
417 } else {
418 if (nvar == maxvar) grow();
419 style[nvar] = EQUAL;
420 num[nvar] = 2;
421 which[nvar] = 0;
422 pad[nvar] = 0;
423 data[nvar] = new char*[num[nvar]];
424 data[nvar][0] = utils::strdup(arg[2]);
425 data[nvar][1] = new char[VALUELENGTH];
426 strcpy(data[nvar][1],"(undefined)");
427 }
428
429 // ATOM
430 // replace pre-existing var if also style ATOM (allows it to be reset)
431 // num = 1, which = 1st value
432 // data = 1 value, string to eval
433
434 } else if (strcmp(arg[1],"atom") == 0) {
435 if (narg != 3) error->all(FLERR,"Illegal variable command");
436 int ivar = find(arg[0]);
437 if (ivar >= 0) {
438 if (style[ivar] != ATOM)
439 error->all(FLERR,"Cannot redefine variable as a different style");
440 delete [] data[ivar][0];
441 data[ivar][0] = utils::strdup(arg[2]);
442 replaceflag = 1;
443 } else {
444 if (nvar == maxvar) grow();
445 style[nvar] = ATOM;
446 num[nvar] = 1;
447 which[nvar] = 0;
448 pad[nvar] = 0;
449 data[nvar] = new char*[num[nvar]];
450 data[nvar][0] = utils::strdup(arg[2]);
451 }
452
453 // VECTOR
454 // replace pre-existing var if also style VECTOR (allows it to be reset)
455 // num = 1, which = 1st value
456 // data = 1 value, string to eval
457
458 } else if (strcmp(arg[1],"vector") == 0) {
459 if (narg != 3) error->all(FLERR,"Illegal variable command");
460 int ivar = find(arg[0]);
461 if (ivar >= 0) {
462 if (style[ivar] != VECTOR)
463 error->all(FLERR,"Cannot redefine variable as a different style");
464 delete [] data[ivar][0];
465 data[ivar][0] = utils::strdup(arg[2]);
466 replaceflag = 1;
467 } else {
468 if (nvar == maxvar) grow();
469 style[nvar] = VECTOR;
470 num[nvar] = 1;
471 which[nvar] = 0;
472 pad[nvar] = 0;
473 data[nvar] = new char*[num[nvar]];
474 data[nvar][0] = utils::strdup(arg[2]);
475 }
476
477 // PYTHON
478 // replace pre-existing var if also style PYTHON (allows it to be reset)
479 // num = 2, which = 1st value
480 // data = 2 values, 1st is Python func to invoke, 2nd is filled by invoke
481
482 } else if (strcmp(arg[1],"python") == 0) {
483 if (narg != 3) error->all(FLERR,"Illegal variable command");
484 if (!python->is_enabled())
485 error->all(FLERR,"LAMMPS is not built with Python embedded");
486 int ivar = find(arg[0]);
487 if (ivar >= 0) {
488 if (style[ivar] != PYTHON)
489 error->all(FLERR,"Cannot redefine variable as a different style");
490 delete [] data[ivar][0];
491 data[ivar][0] = utils::strdup(arg[2]);
492 replaceflag = 1;
493 } else {
494 if (nvar == maxvar) grow();
495 style[nvar] = PYTHON;
496 num[nvar] = 2;
497 which[nvar] = 1;
498 pad[nvar] = 0;
499 data[nvar] = new char*[num[nvar]];
500 data[nvar][0] = utils::strdup(arg[2]);
501 data[nvar][1] = new char[VALUELENGTH];
502 strcpy(data[nvar][1],"(undefined)");
503 }
504
505 // INTERNAL
506 // replace pre-existing var if also style INTERNAL (allows it to be reset)
507 // num = 1, for string representation of dvalue, set by retrieve()
508 // dvalue = numeric initialization from 2nd arg, reset by internal_set()
509
510 } else if (strcmp(arg[1],"internal") == 0) {
511 if (narg != 3) error->all(FLERR,"Illegal variable command");
512 int ivar = find(arg[0]);
513 if (ivar >= 0) {
514 if (style[ivar] != INTERNAL)
515 error->all(FLERR,"Cannot redefine variable as a different style");
516 dvalue[nvar] = utils::numeric(FLERR,arg[2],false,lmp);
517 replaceflag = 1;
518 } else {
519 if (nvar == maxvar) grow();
520 style[nvar] = INTERNAL;
521 num[nvar] = 1;
522 which[nvar] = 0;
523 pad[nvar] = 0;
524 data[nvar] = new char*[num[nvar]];
525 data[nvar][0] = new char[VALUELENGTH];
526 dvalue[nvar] = utils::numeric(FLERR,arg[2],false,lmp);
527 }
528
529 } else error->all(FLERR,"Illegal variable command");
530
531 // set name of variable, if not replacing one flagged with replaceflag
532 // name must be all alphanumeric chars or underscores
533
534 if (replaceflag) return;
535
536 if (!utils::is_id(arg[0]))
537 error->all(FLERR,"Variable name '{}' must have only alphanu"
538 "meric characters or underscores",arg[0]);
539 names[nvar] = utils::strdup(arg[0]);
540 nvar++;
541 }
542
543 /* ----------------------------------------------------------------------
544 convenience function to allow defining a variable from a single string
545 ------------------------------------------------------------------------- */
546
set(const std::string & setcmd)547 void Variable::set(const std::string &setcmd)
548 {
549 std::vector<std::string> args = utils::split_words(setcmd);
550 char **newarg = new char*[args.size()];
551 int i=0;
552 for (const auto &arg : args) {
553 newarg[i++] = (char *)arg.c_str();
554 }
555 set(args.size(),newarg);
556 delete[] newarg;
557 }
558
559 /* ----------------------------------------------------------------------
560 INDEX variable created by command-line argument
561 make it INDEX rather than STRING so cannot be re-defined in input script
562 ------------------------------------------------------------------------- */
563
set(char * name,int narg,char ** arg)564 void Variable::set(char *name, int narg, char **arg)
565 {
566 char **newarg = new char*[2+narg];
567 newarg[0] = name;
568 newarg[1] = (char *) "index";
569 for (int i = 0; i < narg; i++) newarg[2+i] = arg[i];
570 set(2+narg,newarg);
571 delete [] newarg;
572 }
573
574 /* ----------------------------------------------------------------------
575 set existing STRING variable to str
576 return 0 if successful
577 return -1 if variable doesn't exist or isn't a STRING variable
578 called via library interface, so external programs can set variables
579 ------------------------------------------------------------------------- */
580
set_string(const char * name,const char * str)581 int Variable::set_string(const char *name, const char *str)
582 {
583 int ivar = find(name);
584 if (ivar < 0) return -1;
585 if (style[ivar] != STRING) return -1;
586 delete [] data[ivar][0];
587 data[ivar][0] = utils::strdup(str);
588 return 0;
589 }
590
591 /* ----------------------------------------------------------------------
592 increment variable(s)
593 return 0 if OK if successfully incremented
594 return 1 if any variable is exhausted, free the variable to allow re-use
595 ------------------------------------------------------------------------- */
596
next(int narg,char ** arg)597 int Variable::next(int narg, char **arg)
598 {
599 int ivar;
600
601 if (narg == 0) error->all(FLERR,"Illegal next command");
602
603 // check that variables exist and are all the same style
604 // exception: UNIVERSE and ULOOP variables can be mixed in same next command
605
606 for (int iarg = 0; iarg < narg; iarg++) {
607 ivar = find(arg[iarg]);
608 if (ivar < 0)
609 error->all(FLERR,"Invalid variable '{}' in next command",
610 arg[iarg]);
611 if (style[ivar] == ULOOP && style[find(arg[0])] == UNIVERSE) continue;
612 else if (style[ivar] == UNIVERSE && style[find(arg[0])] == ULOOP) continue;
613 else if (style[ivar] != style[find(arg[0])])
614 error->all(FLERR,"All variables in next command must have same style");
615 }
616
617 // invalid styles: STRING, EQUAL, WORLD, ATOM, VECTOR, GETENV,
618 // FORMAT, PYTHON, INTERNAL
619
620 int istyle = style[find(arg[0])];
621 if (istyle == STRING || istyle == EQUAL || istyle == WORLD ||
622 istyle == GETENV || istyle == ATOM || istyle == VECTOR ||
623 istyle == FORMAT || istyle == PYTHON || istyle == INTERNAL)
624 error->all(FLERR,"Invalid variable style with next command");
625
626 // if istyle = UNIVERSE or ULOOP, insure all such variables are incremented
627
628 if (istyle == UNIVERSE || istyle == ULOOP)
629 for (int i = 0; i < nvar; i++) {
630 if (style[i] != UNIVERSE && style[i] != ULOOP) continue;
631 int iarg = 0;
632 for (iarg = 0; iarg < narg; iarg++)
633 if (strcmp(arg[iarg],names[i]) == 0) break;
634 if (iarg == narg)
635 error->universe_one(FLERR,"Next command must list all "
636 "universe and uloop variables");
637 }
638
639 // increment all variables in list
640 // if any variable is exhausted, set flag = 1 and remove var to allow re-use
641
642 int flag = 0;
643
644 if (istyle == INDEX || istyle == LOOP) {
645 for (int iarg = 0; iarg < narg; iarg++) {
646 ivar = find(arg[iarg]);
647 which[ivar]++;
648 if (which[ivar] >= num[ivar]) {
649 flag = 1;
650 remove(ivar);
651 }
652 }
653
654 } else if (istyle == SCALARFILE) {
655
656 for (int iarg = 0; iarg < narg; iarg++) {
657 ivar = find(arg[iarg]);
658 int done = reader[ivar]->read_scalar(data[ivar][0]);
659 if (done) {
660 flag = 1;
661 remove(ivar);
662 }
663 }
664
665 } else if (istyle == ATOMFILE) {
666
667 for (int iarg = 0; iarg < narg; iarg++) {
668 ivar = find(arg[iarg]);
669 int done = reader[ivar]->read_peratom();
670 if (done) {
671 flag = 1;
672 remove(ivar);
673 }
674 }
675
676 } else if (istyle == UNIVERSE || istyle == ULOOP) {
677
678 RanMars *random = nullptr;
679
680 uloop_again:
681
682 // wait until lock file can be created and owned by proc 0 of this world
683 // rename() is not atomic in practice, but no known simple fix
684 // means multiple procs can read/write file at the same time (bad!)
685 // random delays help
686 // delay for random fraction of 1 second before first rename() call
687 // delay for random fraction of 1 second before subsequent tries
688 // when successful, read next available index and Bcast it within my world
689
690 int nextindex = -1;
691 if (me == 0) {
692 int seed = 12345 + universe->me + which[find(arg[0])];
693 if (!random) random = new RanMars(lmp,seed);
694 int delay = (int) (1000000*random->uniform());
695 usleep(delay);
696 while (1) {
697 if (!rename("tmp.lammps.variable","tmp.lammps.variable.lock")) break;
698 delay = (int) (1000000*random->uniform());
699 usleep(delay);
700 }
701
702 // if the file cannot be found, we may have a race with some
703 // other MPI rank that has called rename at the same time
704 // and we have to start over.
705 // if the read is short (we need at least one byte) we try reading again.
706
707 FILE *fp;
708 char buf[64];
709 for (int loopmax = 0; loopmax < 100; ++loopmax) {
710 fp = fopen("tmp.lammps.variable.lock","r");
711 if (fp == nullptr) goto uloop_again;
712
713 buf[0] = buf[1] = '\0';
714 fread(buf,1,64,fp);
715 fclose(fp);
716
717 if (strlen(buf) > 0) {
718 nextindex = atoi(buf);
719 break;
720 }
721 delay = (int) (1000000*random->uniform());
722 usleep(delay);
723 }
724 delete random;
725 random = nullptr;
726
727 if (nextindex < 0)
728 error->one(FLERR,"Unexpected error while incrementing uloop "
729 "style variable. Please contact LAMMPS developers.");
730
731 //printf("READ %d %d\n",universe->me,nextindex);
732 fp = fopen("tmp.lammps.variable.lock","w");
733 fprintf(fp,"%d\n",nextindex+1);
734 //printf("WRITE %d %d\n",universe->me,nextindex+1);
735 fclose(fp);
736 fp = nullptr;
737 rename("tmp.lammps.variable.lock","tmp.lammps.variable");
738 if (universe->uscreen)
739 fprintf(universe->uscreen,
740 "Increment via next: value %d on partition %d\n",
741 nextindex+1,universe->iworld);
742 if (universe->ulogfile)
743 fprintf(universe->ulogfile,
744 "Increment via next: value %d on partition %d\n",
745 nextindex+1,universe->iworld);
746 }
747 MPI_Bcast(&nextindex,1,MPI_INT,0,world);
748
749 // set all variables in list to nextindex
750 // must increment all UNIVERSE and ULOOP variables here
751 // error check above tested for this
752
753 for (int iarg = 0; iarg < narg; iarg++) {
754 ivar = find(arg[iarg]);
755 which[ivar] = nextindex;
756 if (which[ivar] >= num[ivar]) {
757 flag = 1;
758 remove(ivar);
759 }
760 }
761 }
762
763 return flag;
764 }
765
766 /* ----------------------------------------------------------------------
767 search for name in list of variables names
768 return index or -1 if not found
769 ------------------------------------------------------------------------- */
770
find(const char * name)771 int Variable::find(const char *name)
772 {
773 if (name == nullptr) return -1;
774 for (int i = 0; i < nvar; i++)
775 if (strcmp(name,names[i]) == 0) return i;
776 return -1;
777 }
778
779 /* ----------------------------------------------------------------------
780 initialize one atom's storage values in all VarReaders via fix STORE
781 called when atom is created
782 ------------------------------------------------------------------------- */
783
set_arrays(int)784 void Variable::set_arrays(int /*i*/)
785 {
786 for (int i = 0; i < nvar; i++)
787 if (reader[i] && style[i] == ATOMFILE)
788 reader[i]->fixstore->vstore[i] = 0.0;
789 }
790
791 /* ----------------------------------------------------------------------
792 called by python command in input script
793 simply pass input script line args to Python class
794 ------------------------------------------------------------------------- */
795
python_command(int narg,char ** arg)796 void Variable::python_command(int narg, char **arg)
797 {
798 if (!python->is_enabled())
799 error->all(FLERR,"LAMMPS is not built with Python embedded");
800 python->command(narg,arg);
801 }
802
803 /* ----------------------------------------------------------------------
804 return 1 if variable is EQUAL or INTERNAL or PYTHON numeric style, 0 if not
805 this is checked before call to compute_equal() to return a double
806 ------------------------------------------------------------------------- */
807
equalstyle(int ivar)808 int Variable::equalstyle(int ivar)
809 {
810 if (style[ivar] == EQUAL || style[ivar] == INTERNAL) return 1;
811 if (style[ivar] == PYTHON) {
812 int ifunc = python->variable_match(data[ivar][0],names[ivar],1);
813 if (ifunc < 0) return 0;
814 else return 1;
815 }
816 return 0;
817 }
818
819 /* ----------------------------------------------------------------------
820 return 1 if variable is ATOM or ATOMFILE style, 0 if not
821 this is checked before call to compute_atom() to return a vector of doubles
822 ------------------------------------------------------------------------- */
823
atomstyle(int ivar)824 int Variable::atomstyle(int ivar)
825 {
826 if (style[ivar] == ATOM || style[ivar] == ATOMFILE) return 1;
827 return 0;
828 }
829
830 /* ----------------------------------------------------------------------
831 return 1 if variable is VECTOR style, 0 if not
832 this is checked before call to compute_vector() to return a vector of doubles
833 ------------------------------------------------------------------------- */
834
vectorstyle(int ivar)835 int Variable::vectorstyle(int ivar)
836 {
837 if (style[ivar] == VECTOR) return 1;
838 return 0;
839 }
840
841 /* ----------------------------------------------------------------------
842 check if variable with name is PYTHON and matches funcname
843 called by Python class before it invokes a Python function
844 return data storage so Python function can return a value for this variable
845 return nullptr if not a match
846 ------------------------------------------------------------------------- */
847
pythonstyle(char * name,char * funcname)848 char *Variable::pythonstyle(char *name, char *funcname)
849 {
850 int ivar = find(name);
851 if (ivar < 0) return nullptr;
852 if (style[ivar] != PYTHON) return nullptr;
853 if (strcmp(data[ivar][0],funcname) != 0) return nullptr;
854 return data[ivar][1];
855 }
856
857 /* ----------------------------------------------------------------------
858 return 1 if variable is INTERNAL style, 0 if not
859 this is checked before call to set_internal() to assure it can be set
860 ------------------------------------------------------------------------- */
861
internalstyle(int ivar)862 int Variable::internalstyle(int ivar)
863 {
864 if (style[ivar] == INTERNAL) return 1;
865 return 0;
866 }
867
868 /* ----------------------------------------------------------------------
869 return ptr to the data text associated with a variable
870 if INDEX or WORLD or UNIVERSE or STRING or SCALARFILE,
871 return ptr to stored string
872 if LOOP or ULOOP, write int to data[0] and return ptr to string
873 if EQUAL, evaluate variable and put result in str
874 if FORMAT, evaluate its variable and put formatted result in str
875 if GETENV, query environment and put result in str
876 if PYTHON, evaluate Python function, it will put result in str
877 if INTERNAL, convert dvalue and put result in str
878 if ATOM or ATOMFILE or VECTOR, return nullptr
879 return nullptr if no variable with name, or which value is bad,
880 caller must respond
881 ------------------------------------------------------------------------- */
882
retrieve(const char * name)883 char *Variable::retrieve(const char *name)
884 {
885 int ivar = find(name);
886 if (ivar < 0) return nullptr;
887 if (which[ivar] >= num[ivar]) return nullptr;
888
889 if (eval_in_progress[ivar])
890 print_var_error(FLERR,"has a circular dependency",ivar);
891
892 eval_in_progress[ivar] = 1;
893
894 char *str = nullptr;
895 if (style[ivar] == INDEX || style[ivar] == WORLD ||
896 style[ivar] == UNIVERSE || style[ivar] == STRING ||
897 style[ivar] == SCALARFILE) {
898 str = data[ivar][which[ivar]];
899 } else if (style[ivar] == LOOP || style[ivar] == ULOOP) {
900 char result[16];
901 if (pad[ivar] == 0) sprintf(result,"%d",which[ivar]+1);
902 else {
903 char padstr[16];
904 sprintf(padstr,"%%0%dd",pad[ivar]);
905 sprintf(result,padstr,which[ivar]+1);
906 }
907 delete [] data[ivar][0];
908 str = data[ivar][0] = utils::strdup(result);
909 } else if (style[ivar] == EQUAL) {
910 double answer = evaluate(data[ivar][0],nullptr,ivar);
911 sprintf(data[ivar][1],"%.15g",answer);
912 str = data[ivar][1];
913 } else if (style[ivar] == FORMAT) {
914 int jvar = find(data[ivar][0]);
915 if (jvar == -1) return nullptr;
916 if (!equalstyle(jvar)) return nullptr;
917 double answer = compute_equal(jvar);
918 sprintf(data[ivar][2],data[ivar][1],answer);
919 str = data[ivar][2];
920 } else if (style[ivar] == GETENV) {
921 const char *result = getenv(data[ivar][0]);
922 if (result == nullptr) result = (const char *) "";
923 delete [] data[ivar][1];
924 str = data[ivar][1] = utils::strdup(result);
925 } else if (style[ivar] == PYTHON) {
926 int ifunc = python->variable_match(data[ivar][0],name,0);
927 if (ifunc < 0)
928 error->all(FLERR,"Python variable {} does not match "
929 "Python function {}", name, data[ivar][0]);
930 python->invoke_function(ifunc,data[ivar][1]);
931 str = data[ivar][1];
932 // if Python func returns a string longer than VALUELENGTH
933 // then the Python class stores the result, query it via long_string()
934 char *strlong = python->long_string(ifunc);
935 if (strlong) str = strlong;
936 } else if (style[ivar] == INTERNAL) {
937 sprintf(data[ivar][0],"%.15g",dvalue[ivar]);
938 str = data[ivar][0];
939 } else if (style[ivar] == ATOM || style[ivar] == ATOMFILE ||
940 style[ivar] == VECTOR) return nullptr;
941
942 eval_in_progress[ivar] = 0;
943
944 return str;
945 }
946
947 /* ----------------------------------------------------------------------
948 return result of equal-style variable evaluation
949 can be EQUAL or INTERNAL style or PYTHON numeric style
950 for PYTHON, don't need to check python->variable_match() error return,
951 since caller will have already checked via equalstyle()
952 ------------------------------------------------------------------------- */
953
compute_equal(int ivar)954 double Variable::compute_equal(int ivar)
955 {
956 if (eval_in_progress[ivar])
957 print_var_error(FLERR,"has a circular dependency",ivar);
958
959 eval_in_progress[ivar] = 1;
960
961 double value = 0.0;
962 if (style[ivar] == EQUAL) value = evaluate(data[ivar][0],nullptr,ivar);
963 else if (style[ivar] == INTERNAL) value = dvalue[ivar];
964 else if (style[ivar] == PYTHON) {
965 int ifunc = python->find(data[ivar][0]);
966 if (ifunc < 0)
967 print_var_error(FLERR,fmt::format("cannot find python function {}",
968 data[ivar][0]),ivar);
969 python->invoke_function(ifunc,data[ivar][1]);
970 value = atof(data[ivar][1]);
971 }
972
973 eval_in_progress[ivar] = 0;
974 return value;
975 }
976
977 /* ----------------------------------------------------------------------
978 return result of immediate equal-style variable evaluation
979 called from Input::substitute()
980 don't need to flag eval_in_progress since is an immediate variable
981 ------------------------------------------------------------------------- */
982
compute_equal(const std::string & str)983 double Variable::compute_equal(const std::string &str)
984 {
985 char *ptr = utils::strdup(str);
986 double val = evaluate(ptr,nullptr,-1);
987 delete[] ptr;
988 return val;
989 }
990
991 /* ----------------------------------------------------------------------
992 compute result of atom-style and atomfile-style variable evaluation
993 only computed for atoms in igroup, else result is 0.0
994 answers are placed every stride locations into result
995 if sumflag, add variable values to existing result
996 ------------------------------------------------------------------------- */
997
compute_atom(int ivar,int igroup,double * result,int stride,int sumflag)998 void Variable::compute_atom(int ivar, int igroup,
999 double *result, int stride, int sumflag)
1000 {
1001 Tree *tree = nullptr;
1002 double *vstore;
1003
1004 if (eval_in_progress[ivar])
1005 print_var_error(FLERR,"has a circular dependency",ivar);
1006
1007 eval_in_progress[ivar] = 1;
1008
1009 if (style[ivar] == ATOM) {
1010 treetype = ATOM;
1011 evaluate(data[ivar][0],&tree,ivar);
1012 collapse_tree(tree);
1013 } else vstore = reader[ivar]->fixstore->vstore;
1014
1015 if (result == nullptr) {
1016 if (style[ivar] == ATOM) free_tree(tree);
1017 eval_in_progress[ivar] = 0;
1018 return;
1019 }
1020
1021 int groupbit = group->bitmask[igroup];
1022 int *mask = atom->mask;
1023 int nlocal = atom->nlocal;
1024
1025 if (style[ivar] == ATOM) {
1026 if (sumflag == 0) {
1027 int m = 0;
1028 for (int i = 0; i < nlocal; i++) {
1029 if (mask[i] & groupbit) result[m] = eval_tree(tree,i);
1030 else result[m] = 0.0;
1031 m += stride;
1032 }
1033
1034 } else {
1035 int m = 0;
1036 for (int i = 0; i < nlocal; i++) {
1037 if (mask[i] & groupbit) result[m] += eval_tree(tree,i);
1038 m += stride;
1039 }
1040 }
1041
1042 } else {
1043 if (sumflag == 0) {
1044 int m = 0;
1045 for (int i = 0; i < nlocal; i++) {
1046 if (mask[i] & groupbit) result[m] = vstore[i];
1047 else result[m] = 0.0;
1048 m += stride;
1049 }
1050
1051 } else {
1052 int m = 0;
1053 for (int i = 0; i < nlocal; i++) {
1054 if (mask[i] & groupbit) result[m] += vstore[i];
1055 m += stride;
1056 }
1057 }
1058 }
1059
1060 if (style[ivar] == ATOM) free_tree(tree);
1061 eval_in_progress[ivar] = 0;
1062 }
1063
1064 /* ----------------------------------------------------------------------
1065 compute result of vector-style variable evaluation
1066 return length of vector and result pointer to vector values
1067 if length == 0 or -1 (mismatch), generate an error
1068 if variable already computed on this timestep, just return
1069 else evaluate the formula and its length, store results in VecVar entry
1070 ------------------------------------------------------------------------- */
1071
compute_vector(int ivar,double ** result)1072 int Variable::compute_vector(int ivar, double **result)
1073 {
1074 Tree *tree = nullptr;
1075 if (vecs[ivar].currentstep == update->ntimestep) {
1076 *result = vecs[ivar].values;
1077 return vecs[ivar].n;
1078 }
1079
1080 if (eval_in_progress[ivar])
1081 print_var_error(FLERR,"has a circular dependency",ivar);
1082
1083 eval_in_progress[ivar] = 1;
1084
1085 treetype = VECTOR;
1086 evaluate(data[ivar][0],&tree,ivar);
1087 collapse_tree(tree);
1088 int nlen = size_tree_vector(tree);
1089 if (nlen == 0)
1090 print_var_error(FLERR,"Vector-style variable has zero length",ivar);
1091
1092 if (nlen < 0)
1093 print_var_error(FLERR,"Inconsistent lengths in vector-style variable",ivar);
1094
1095 // (re)allocate space for results if necessary
1096
1097 if (nlen > vecs[ivar].nmax) {
1098 memory->destroy(vecs[ivar].values);
1099 vecs[ivar].nmax = nlen;
1100 memory->create(vecs[ivar].values,vecs[ivar].nmax,"variable:values");
1101 }
1102
1103 vecs[ivar].n = nlen;
1104 vecs[ivar].currentstep = update->ntimestep;
1105 double *vec = vecs[ivar].values;
1106 for (int i = 0; i < nlen; i++)
1107 vec[i] = eval_tree(tree,i);
1108
1109 free_tree(tree);
1110 eval_in_progress[ivar] = 0;
1111
1112 *result = vec;
1113 return nlen;
1114 }
1115
1116 /* ----------------------------------------------------------------------
1117 set value stored by INTERNAL style ivar
1118 ------------------------------------------------------------------------- */
1119
internal_set(int ivar,double value)1120 void Variable::internal_set(int ivar, double value)
1121 {
1122 dvalue[ivar] = value;
1123 }
1124
1125 /* ----------------------------------------------------------------------
1126 remove Nth variable from list and compact list
1127 delete reader explicitly if it exists
1128 ------------------------------------------------------------------------- */
1129
remove(int n)1130 void Variable::remove(int n)
1131 {
1132 delete [] names[n];
1133 if (style[n] == LOOP || style[n] == ULOOP) delete [] data[n][0];
1134 else for (int i = 0; i < num[n]; i++) delete [] data[n][i];
1135 delete [] data[n];
1136 delete reader[n];
1137
1138 for (int i = n+1; i < nvar; i++) {
1139 names[i-1] = names[i];
1140 style[i-1] = style[i];
1141 num[i-1] = num[i];
1142 which[i-1] = which[i];
1143 pad[i-1] = pad[i];
1144 reader[i-1] = reader[i];
1145 data[i-1] = data[i];
1146 dvalue[i-1] = dvalue[i];
1147 }
1148 nvar--;
1149 }
1150
1151 /* ----------------------------------------------------------------------
1152 make space in arrays for new variable
1153 ------------------------------------------------------------------------- */
1154
grow()1155 void Variable::grow()
1156 {
1157 int old = maxvar;
1158 maxvar += VARDELTA;
1159 names = (char **) memory->srealloc(names,maxvar*sizeof(char *),"var:names");
1160 memory->grow(style,maxvar,"var:style");
1161 memory->grow(num,maxvar,"var:num");
1162 memory->grow(which,maxvar,"var:which");
1163 memory->grow(pad,maxvar,"var:pad");
1164
1165 reader = (VarReader **)
1166 memory->srealloc(reader,maxvar*sizeof(VarReader *),"var:reader");
1167 for (int i = old; i < maxvar; i++) reader[i] = nullptr;
1168
1169 data = (char ***) memory->srealloc(data,maxvar*sizeof(char **),"var:data");
1170 memory->grow(dvalue,maxvar,"var:dvalue");
1171
1172 vecs = (VecVar *) memory->srealloc(vecs,maxvar*sizeof(VecVar),"var:vecvar");
1173 for (int i = old; i < maxvar; i++) {
1174 vecs[i].nmax = 0;
1175 vecs[i].currentstep = -1;
1176 vecs[i].values = nullptr;
1177 }
1178
1179 memory->grow(eval_in_progress,maxvar,"var:eval_in_progress");
1180 for (int i = 0; i < maxvar; i++) eval_in_progress[i] = 0;
1181 }
1182
1183 /* ----------------------------------------------------------------------
1184 copy narg strings from **from to **to, and allocate space for them
1185 ------------------------------------------------------------------------- */
1186
copy(int narg,char ** from,char ** to)1187 void Variable::copy(int narg, char **from, char **to)
1188 {
1189 for (int i = 0; i < narg; i++)
1190 to[i] = utils::strdup(from[i]);
1191 }
1192
1193 /* ----------------------------------------------------------------------
1194 recursive evaluation of a string str
1195 str is an equal-style or atom-style or vector-style formula
1196 containing one or more items:
1197 number = 0.0, -5.45, 2.8e-4, ...
1198 constant = PI, version, yes, no, on, off
1199 thermo keyword = ke, vol, atoms, ...
1200 math operation = (),-x,x+y,x-y,x*y,x/y,x^y,
1201 x==y,x!=y,x<y,x<=y,x>y,x>=y,x&&y,x||y,
1202 sqrt(x),exp(x),ln(x),log(x),abs(x),
1203 sin(x),cos(x),tan(x),asin(x),atan2(y,x),...
1204 group function = count(group), mass(group), xcm(group,x), ...
1205 special function = sum(x),min(x), ...
1206 atom value = x[i], y[i], vx[i], ...
1207 atom vector = x, y, vx, ...
1208 compute = c_ID, c_ID[i], c_ID[i][j]
1209 fix = f_ID, f_ID[i], f_ID[i][j]
1210 variable = v_name, v_name[i]
1211 equal-style variables passes in tree = nullptr:
1212 evaluate the formula, return result as a double
1213 atom-style and vector-style variables pass in tree = non-nullptr:
1214 parse the formula but do not evaluate it
1215 create a parse tree and return it
1216 ------------------------------------------------------------------------- */
1217
evaluate(char * str,Tree ** tree,int ivar)1218 double Variable::evaluate(char *str, Tree **tree, int ivar)
1219 {
1220 int op,opprevious;
1221 double value1,value2;
1222 char onechar;
1223 char *ptr;
1224
1225 double argstack[MAXLEVEL];
1226 Tree *treestack[MAXLEVEL];
1227 int opstack[MAXLEVEL];
1228 int nargstack = 0;
1229 int ntreestack = 0;
1230 int nopstack = 0;
1231
1232 int i = 0;
1233 int expect = ARG;
1234
1235 if (str == nullptr)
1236 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
1237
1238 while (1) {
1239 onechar = str[i];
1240
1241 // whitespace: just skip
1242
1243 if (isspace(onechar)) i++;
1244
1245 // ----------------
1246 // parentheses: recursively evaluate contents of parens
1247 // ----------------
1248
1249 else if (onechar == '(') {
1250 if (expect == OP)
1251 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
1252 expect = OP;
1253
1254 char *contents = nullptr;
1255 i = find_matching_paren(str,i,contents,ivar);
1256 i++;
1257
1258 // evaluate contents and push on stack
1259
1260 if (tree) {
1261 Tree *newtree = nullptr;
1262 evaluate(contents,&newtree,ivar);
1263 treestack[ntreestack++] = newtree;
1264 } else argstack[nargstack++] = evaluate(contents,nullptr,ivar);
1265
1266 delete [] contents;
1267
1268 // ----------------
1269 // number: push value onto stack
1270 // ----------------
1271
1272 } else if (isdigit(onechar) || onechar == '.') {
1273 if (expect == OP)
1274 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
1275 expect = OP;
1276
1277 // istop = end of number, including scientific notation
1278
1279 int istart = i;
1280 while (isdigit(str[i]) || str[i] == '.') i++;
1281 if (str[i] == 'e' || str[i] == 'E') {
1282 i++;
1283 if (str[i] == '+' || str[i] == '-') i++;
1284 while (isdigit(str[i])) i++;
1285 }
1286 int istop = i - 1;
1287
1288 int n = istop - istart + 1;
1289 char *number = new char[n+1];
1290 strncpy(number,&str[istart],n);
1291 number[n] = '\0';
1292
1293 if (tree) {
1294 Tree *newtree = new Tree();
1295 newtree->type = VALUE;
1296 newtree->value = atof(number);
1297 treestack[ntreestack++] = newtree;
1298 } else argstack[nargstack++] = atof(number);
1299
1300 delete [] number;
1301
1302 // ----------------
1303 // letter: c_ID, c_ID[], c_ID[][], f_ID, f_ID[], f_ID[][],
1304 // v_name, v_name[], exp(), xcm(,), x, x[], PI, vol
1305 // ----------------
1306
1307 } else if (isalpha(onechar)) {
1308 if (expect == OP)
1309 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
1310 expect = OP;
1311
1312 // istop = end of word
1313 // word = all alphanumeric or underscore
1314
1315 int istart = i;
1316 while (isalnum(str[i]) || str[i] == '_') i++;
1317 int istop = i-1;
1318
1319 int n = istop - istart + 1;
1320 char *word = new char[n+1];
1321 strncpy(word,&str[istart],n);
1322 word[n] = '\0';
1323
1324 // ----------------
1325 // compute
1326 // ----------------
1327
1328 if (strncmp(word,"c_",2) == 0 || strncmp(word,"C_",2) == 0) {
1329 if (domain->box_exist == 0)
1330 print_var_error(FLERR,"Variable evaluation before "
1331 "simulation box is defined",ivar);
1332
1333 // uppercase used to force access of
1334 // global vector vs global scalar, and global array vs global vector
1335
1336 int lowercase = 1;
1337 if (word[0] == 'C') lowercase = 0;
1338
1339 int icompute = modify->find_compute(word+2);
1340 if (icompute < 0) {
1341 std::string mesg = "Invalid compute ID '";
1342 mesg += (word+2);
1343 mesg += "' in variable formula";
1344 print_var_error(FLERR,mesg,ivar);
1345 }
1346 Compute *compute = modify->compute[icompute];
1347
1348 // parse zero or one or two trailing brackets
1349 // point i beyond last bracket
1350 // nbracket = # of bracket pairs
1351 // index1,index2 = int inside each bracket pair, possibly an atom ID
1352
1353 int nbracket;
1354 tagint index1,index2;
1355 if (str[i] != '[') nbracket = 0;
1356 else {
1357 nbracket = 1;
1358 ptr = &str[i];
1359 index1 = int_between_brackets(ptr,1);
1360 i = ptr-str+1;
1361 if (str[i] == '[') {
1362 nbracket = 2;
1363 ptr = &str[i];
1364 index2 = int_between_brackets(ptr,1);
1365 i = ptr-str+1;
1366 }
1367 }
1368
1369 // c_ID = scalar from global scalar, must be lowercase
1370
1371 if (nbracket == 0 && compute->scalar_flag && lowercase) {
1372
1373 if (update->whichflag == 0) {
1374 if (compute->invoked_scalar != update->ntimestep)
1375 print_var_error(FLERR,"Compute used in variable between "
1376 "runs is not current",ivar);
1377 } else if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) {
1378 compute->compute_scalar();
1379 compute->invoked_flag |= Compute::INVOKED_SCALAR;
1380 }
1381
1382 value1 = compute->scalar;
1383 if (tree) {
1384 Tree *newtree = new Tree();
1385 newtree->type = VALUE;
1386 newtree->value = value1;
1387 treestack[ntreestack++] = newtree;
1388 } else argstack[nargstack++] = value1;
1389
1390 // c_ID[i] = scalar from global vector, must be lowercase
1391
1392 } else if (nbracket == 1 && compute->vector_flag && lowercase) {
1393
1394 if (index1 > compute->size_vector &&
1395 compute->size_vector_variable == 0)
1396 print_var_error(FLERR,"Variable formula compute vector "
1397 "is accessed out-of-range",ivar,0);
1398 if (update->whichflag == 0) {
1399 if (compute->invoked_vector != update->ntimestep)
1400 print_var_error(FLERR,"Compute used in variable between runs "
1401 "is not current",ivar);
1402 } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
1403 compute->compute_vector();
1404 compute->invoked_flag |= Compute::INVOKED_VECTOR;
1405 }
1406
1407 if (compute->size_vector_variable &&
1408 index1 > compute->size_vector) value1 = 0.0;
1409 else value1 = compute->vector[index1-1];
1410 if (tree) {
1411 Tree *newtree = new Tree();
1412 newtree->type = VALUE;
1413 newtree->value = value1;
1414 treestack[ntreestack++] = newtree;
1415 } else argstack[nargstack++] = value1;
1416
1417 // c_ID[i][j] = scalar from global array, must be lowercase
1418
1419 } else if (nbracket == 2 && compute->array_flag && lowercase) {
1420
1421 if (index1 > compute->size_array_rows &&
1422 compute->size_array_rows_variable == 0)
1423 print_var_error(FLERR,"Variable formula compute array "
1424 "is accessed out-of-range",ivar,0);
1425 if (index2 > compute->size_array_cols)
1426 print_var_error(FLERR,"Variable formula compute array "
1427 "is accessed out-of-range",ivar,0);
1428 if (update->whichflag == 0) {
1429 if (compute->invoked_array != update->ntimestep)
1430 print_var_error(FLERR,"Compute used in variable between runs "
1431 "is not current",ivar);
1432 } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) {
1433 compute->compute_array();
1434 compute->invoked_flag |= Compute::INVOKED_ARRAY;
1435 }
1436
1437 if (compute->size_array_rows_variable &&
1438 index1 > compute->size_array_rows) value1 = 0.0;
1439 else value1 = compute->array[index1-1][index2-1];
1440 if (tree) {
1441 Tree *newtree = new Tree();
1442 newtree->type = VALUE;
1443 newtree->value = value1;
1444 treestack[ntreestack++] = newtree;
1445 } else argstack[nargstack++] = value1;
1446
1447 // c_ID = vector from global vector, lowercase or uppercase
1448
1449 } else if (nbracket == 0 && compute->vector_flag) {
1450
1451 if (tree == nullptr)
1452 print_var_error(FLERR,"Compute global vector in "
1453 "equal-style variable formula",ivar);
1454 if (treetype == ATOM)
1455 print_var_error(FLERR,"Compute global vector in "
1456 "atom-style variable formula",ivar);
1457 if (compute->size_vector == 0)
1458 print_var_error(FLERR,"Variable formula compute "
1459 "vector is zero length",ivar);
1460 if (update->whichflag == 0) {
1461 if (compute->invoked_vector != update->ntimestep)
1462 print_var_error(FLERR,"Compute used in variable between "
1463 "runs is not current",ivar);
1464 } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
1465 compute->compute_vector();
1466 compute->invoked_flag |= Compute::INVOKED_VECTOR;
1467 }
1468
1469 Tree *newtree = new Tree();
1470 newtree->type = VECTORARRAY;
1471 newtree->array = compute->vector;
1472 newtree->nvector = compute->size_vector;
1473 newtree->nstride = 1;
1474 treestack[ntreestack++] = newtree;
1475
1476 // c_ID[i] = vector from global array, lowercase or uppercase
1477
1478 } else if (nbracket == 1 && compute->array_flag) {
1479
1480 if (tree == nullptr)
1481 print_var_error(FLERR,"Compute global vector in "
1482 "equal-style variable formula",ivar);
1483 if (treetype == ATOM)
1484 print_var_error(FLERR,"Compute global vector in "
1485 "atom-style variable formula",ivar);
1486 if (compute->size_array_rows == 0)
1487 print_var_error(FLERR,"Variable formula compute array "
1488 "is zero length",ivar);
1489 if (update->whichflag == 0) {
1490 if (compute->invoked_array != update->ntimestep)
1491 print_var_error(FLERR,"Compute used in variable between "
1492 "runs is not current",ivar);
1493 } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) {
1494 compute->compute_array();
1495 compute->invoked_flag |= Compute::INVOKED_ARRAY;
1496 }
1497
1498 Tree *newtree = new Tree();
1499 newtree->type = VECTORARRAY;
1500 newtree->array = &compute->array[0][index1-1];
1501 newtree->nvector = compute->size_array_rows;
1502 newtree->nstride = compute->size_array_cols;
1503 treestack[ntreestack++] = newtree;
1504
1505 // c_ID[i] = scalar from per-atom vector
1506
1507 } else if (nbracket == 1 && compute->peratom_flag &&
1508 compute->size_peratom_cols == 0) {
1509
1510 if (update->whichflag == 0) {
1511 if (compute->invoked_peratom != update->ntimestep)
1512 print_var_error(FLERR,"Compute used in variable "
1513 "between runs is not current",ivar);
1514 } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
1515 compute->compute_peratom();
1516 compute->invoked_flag |= Compute::INVOKED_PERATOM;
1517 }
1518
1519 peratom2global(1,nullptr,compute->vector_atom,1,index1,
1520 tree,treestack,ntreestack,argstack,nargstack);
1521
1522 // c_ID[i][j] = scalar from per-atom array
1523
1524 } else if (nbracket == 2 && compute->peratom_flag &&
1525 compute->size_peratom_cols > 0) {
1526
1527 if (index2 > compute->size_peratom_cols)
1528 print_var_error(FLERR,"Variable formula compute "
1529 "array is accessed out-of-range",ivar,0);
1530 if (update->whichflag == 0) {
1531 if (compute->invoked_peratom != update->ntimestep)
1532 print_var_error(FLERR,"Compute used in variable "
1533 "between runs is not current",ivar);
1534 } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
1535 compute->compute_peratom();
1536 compute->invoked_flag |= Compute::INVOKED_PERATOM;
1537 }
1538
1539 if (compute->array_atom)
1540 peratom2global(1,nullptr,&compute->array_atom[0][index2-1],
1541 compute->size_peratom_cols,index1,
1542 tree,treestack,ntreestack,argstack,nargstack);
1543 else
1544 peratom2global(1,nullptr,nullptr,
1545 compute->size_peratom_cols,index1,
1546 tree,treestack,ntreestack,argstack,nargstack);
1547
1548 // c_ID = vector from per-atom vector
1549
1550 } else if (nbracket == 0 && compute->peratom_flag &&
1551 compute->size_peratom_cols == 0) {
1552
1553 if (tree == nullptr)
1554 print_var_error(FLERR,"Per-atom compute in "
1555 "equal-style variable formula",ivar);
1556 if (treetype == VECTOR)
1557 print_var_error(FLERR,"Per-atom compute in "
1558 "vector-style variable formula",ivar);
1559 if (update->whichflag == 0) {
1560 if (compute->invoked_peratom != update->ntimestep)
1561 print_var_error(FLERR,"Compute used in variable "
1562 "between runs is not current",ivar);
1563 } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
1564 compute->compute_peratom();
1565 compute->invoked_flag |= Compute::INVOKED_PERATOM;
1566 }
1567
1568 Tree *newtree = new Tree();
1569 newtree->type = ATOMARRAY;
1570 newtree->array = compute->vector_atom;
1571 newtree->nstride = 1;
1572 treestack[ntreestack++] = newtree;
1573
1574 // c_ID[i] = vector from per-atom array
1575
1576 } else if (nbracket == 1 && compute->peratom_flag &&
1577 compute->size_peratom_cols > 0) {
1578
1579 if (tree == nullptr)
1580 print_var_error(FLERR,"Per-atom compute in "
1581 "equal-style variable formula",ivar);
1582 if (treetype == VECTOR)
1583 print_var_error(FLERR,"Per-atom compute in "
1584 "vector-style variable formula",ivar);
1585 if (index1 > compute->size_peratom_cols)
1586 print_var_error(FLERR,"Variable formula compute array "
1587 "is accessed out-of-range",ivar,0);
1588 if (update->whichflag == 0) {
1589 if (compute->invoked_peratom != update->ntimestep)
1590 print_var_error(FLERR,"Compute used in variable "
1591 "between runs is not current",ivar);
1592 } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
1593 compute->compute_peratom();
1594 compute->invoked_flag |= Compute::INVOKED_PERATOM;
1595 }
1596
1597 Tree *newtree = new Tree();
1598 newtree->type = ATOMARRAY;
1599 if (compute->array_atom)
1600 newtree->array = &compute->array_atom[0][index1-1];
1601 newtree->nstride = compute->size_peratom_cols;
1602 treestack[ntreestack++] = newtree;
1603
1604 } else if (nbracket == 1 && compute->local_flag) {
1605 print_var_error(FLERR,"Cannot access local data via indexing",ivar);
1606 } else print_var_error(FLERR,
1607 "Mismatched compute in variable formula",ivar);
1608
1609 // ----------------
1610 // fix
1611 // ----------------
1612
1613 } else if (strncmp(word,"f_",2) == 0 || strncmp(word,"F_",2) == 0) {
1614 if (domain->box_exist == 0)
1615 print_var_error(FLERR,"Variable evaluation before "
1616 "simulation box is defined",ivar);
1617
1618 // uppercase used to force access of
1619 // global vector vs global scalar, and global array vs global vector
1620
1621 int lowercase = 1;
1622 if (word[0] == 'F') lowercase = 0;
1623
1624 int ifix = modify->find_fix(word+2);
1625 if (ifix < 0)
1626 print_var_error(FLERR,fmt::format("Invalid fix ID '{}' in variable"
1627 " formula",word+2),ivar);
1628 Fix *fix = modify->fix[ifix];
1629
1630 // parse zero or one or two trailing brackets
1631 // point i beyond last bracket
1632 // nbracket = # of bracket pairs
1633 // index1,index2 = int inside each bracket pair, possibly an atom ID
1634
1635 int nbracket;
1636 tagint index1,index2;
1637 if (str[i] != '[') nbracket = 0;
1638 else {
1639 nbracket = 1;
1640 ptr = &str[i];
1641 index1 = int_between_brackets(ptr,1);
1642 i = ptr-str+1;
1643 if (str[i] == '[') {
1644 nbracket = 2;
1645 ptr = &str[i];
1646 index2 = int_between_brackets(ptr,1);
1647 i = ptr-str+1;
1648 }
1649 }
1650
1651 // f_ID = scalar from global scalar, must be lowercase
1652
1653 if (nbracket == 0 && fix->scalar_flag && lowercase) {
1654
1655 if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
1656 print_var_error(FLERR,"Fix in variable not computed "
1657 "at a compatible time",ivar);
1658
1659 value1 = fix->compute_scalar();
1660 if (tree) {
1661 Tree *newtree = new Tree();
1662 newtree->type = VALUE;
1663 newtree->value = value1;
1664 treestack[ntreestack++] = newtree;
1665 } else argstack[nargstack++] = value1;
1666
1667 // f_ID[i] = scalar from global vector, must be lowercase
1668
1669 } else if (nbracket == 1 && fix->vector_flag && lowercase) {
1670
1671 if (index1 > fix->size_vector &&
1672 fix->size_vector_variable == 0)
1673 print_var_error(FLERR,"Variable formula fix vector is "
1674 "accessed out-of-range",ivar,0);
1675 if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
1676 print_var_error(FLERR,"Fix in variable not computed "
1677 "at a compatible time",ivar);
1678
1679 value1 = fix->compute_vector(index1-1);
1680 if (tree) {
1681 Tree *newtree = new Tree();
1682 newtree->type = VALUE;
1683 newtree->value = value1;
1684 treestack[ntreestack++] = newtree;
1685 } else argstack[nargstack++] = value1;
1686
1687 // f_ID[i][j] = scalar from global array, must be lowercase
1688
1689 } else if (nbracket == 2 && fix->array_flag && lowercase) {
1690
1691 if (index1 > fix->size_array_rows &&
1692 fix->size_array_rows_variable == 0)
1693 print_var_error(FLERR,"Variable formula fix array is "
1694 "accessed out-of-range",ivar,0);
1695 if (index2 > fix->size_array_cols)
1696 print_var_error(FLERR,"Variable formula fix array is "
1697 "accessed out-of-range",ivar,0);
1698 if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
1699 print_var_error(FLERR,"Fix in variable not computed at a "
1700 "compatible time",ivar);
1701
1702 value1 = fix->compute_array(index1-1,index2-1);
1703 if (tree) {
1704 Tree *newtree = new Tree();
1705 newtree->type = VALUE;
1706 newtree->value = value1;
1707 treestack[ntreestack++] = newtree;
1708 } else argstack[nargstack++] = value1;
1709
1710 // f_ID = vector from global vector, lowercase or uppercase
1711
1712 } else if (nbracket == 0 && fix->vector_flag) {
1713
1714 if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
1715 print_var_error(FLERR,"Fix in variable not computed at "
1716 "compatible time",ivar);
1717 if (tree == nullptr)
1718 print_var_error(FLERR,"Fix global vector in "
1719 "equal-style variable formula",ivar);
1720 if (treetype == ATOM)
1721 print_var_error(FLERR,"Fix global vector in "
1722 "atom-style variable formula",ivar);
1723 if (fix->size_vector == 0)
1724 print_var_error(FLERR,"Variable formula "
1725 "fix vector is zero length",ivar);
1726
1727 int nvec = fix->size_vector;
1728 double *vec;
1729 memory->create(vec,nvec,"variable:values");
1730 for (int m = 0; m < nvec; m++)
1731 vec[m] = fix->compute_vector(m);
1732
1733 Tree *newtree = new Tree();
1734 newtree->type = VECTORARRAY;
1735 newtree->array = vec;
1736 newtree->nvector = nvec;
1737 newtree->nstride = 1;
1738 newtree->selfalloc = 1;
1739 treestack[ntreestack++] = newtree;
1740
1741 // f_ID[i] = vector from global array, lowercase or uppercase
1742
1743 } else if (nbracket == 1 && fix->array_flag) {
1744
1745 if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
1746 print_var_error(FLERR,"Fix in variable not computed "
1747 "at a compatible time",ivar);
1748 if (tree == nullptr)
1749 print_var_error(FLERR,"Fix global vector in "
1750 "equal-style variable formula",ivar);
1751 if (treetype == ATOM)
1752 print_var_error(FLERR,"Fix global vector in "
1753 "atom-style variable formula",ivar);
1754 if (fix->size_array_rows == 0)
1755 print_var_error(FLERR,"Variable formula fix array is "
1756 "zero length",ivar);
1757
1758 int nvec = fix->size_array_rows;
1759 double *vec;
1760 memory->create(vec,nvec,"variable:values");
1761 for (int m = 0; m < nvec; m++)
1762 vec[m] = fix->compute_array(m,index1-1);
1763
1764 Tree *newtree = new Tree();
1765 newtree->type = VECTORARRAY;
1766 newtree->array = vec;
1767 newtree->nvector = nvec;
1768 newtree->nstride = 1;
1769 newtree->selfalloc = 1;
1770 treestack[ntreestack++] = newtree;
1771
1772 // f_ID[i] = scalar from per-atom vector
1773
1774 } else if (nbracket == 1 && fix->peratom_flag &&
1775 fix->size_peratom_cols == 0) {
1776
1777 if (update->whichflag > 0 &&
1778 update->ntimestep % fix->peratom_freq)
1779 print_var_error(FLERR,"Fix in variable not computed "
1780 "at a compatible time",ivar);
1781
1782 peratom2global(1,nullptr,fix->vector_atom,1,index1,
1783 tree,treestack,ntreestack,argstack,nargstack);
1784
1785 // f_ID[i][j] = scalar from per-atom array
1786
1787 } else if (nbracket == 2 && fix->peratom_flag &&
1788 fix->size_peratom_cols > 0) {
1789
1790 if (index2 > fix->size_peratom_cols)
1791 print_var_error(FLERR,"Variable formula fix array is "
1792 "accessed out-of-range",ivar,0);
1793 if (update->whichflag > 0 &&
1794 update->ntimestep % fix->peratom_freq)
1795 print_var_error(FLERR,"Fix in variable not computed "
1796 "at a compatible time",ivar);
1797
1798 if (fix->array_atom)
1799 peratom2global(1,nullptr,&fix->array_atom[0][index2-1],
1800 fix->size_peratom_cols,index1,
1801 tree,treestack,ntreestack,argstack,nargstack);
1802 else
1803 peratom2global(1,nullptr,nullptr,
1804 fix->size_peratom_cols,index1,
1805 tree,treestack,ntreestack,argstack,nargstack);
1806
1807 // f_ID = vector from per-atom vector
1808
1809 } else if (nbracket == 0 && fix->peratom_flag &&
1810 fix->size_peratom_cols == 0) {
1811
1812 if (tree == nullptr)
1813 print_var_error(FLERR,"Per-atom fix in "
1814 "equal-style variable formula",ivar);
1815 if (update->whichflag > 0 &&
1816 update->ntimestep % fix->peratom_freq)
1817 print_var_error(FLERR,"Fix in variable not computed at "
1818 "compatible time",ivar);
1819
1820 Tree *newtree = new Tree();
1821 newtree->type = ATOMARRAY;
1822 newtree->array = fix->vector_atom;
1823 newtree->nstride = 1;
1824 treestack[ntreestack++] = newtree;
1825
1826 // f_ID[i] = vector from per-atom array
1827
1828 } else if (nbracket == 1 && fix->peratom_flag &&
1829 fix->size_peratom_cols > 0) {
1830
1831 if (tree == nullptr)
1832 print_var_error(FLERR,"Per-atom fix in "
1833 "equal-style variable formula",ivar);
1834 if (index1 > fix->size_peratom_cols)
1835 print_var_error(FLERR,"Variable formula fix array "
1836 "is accessed out-of-range",ivar,0);
1837 if (update->whichflag > 0 &&
1838 update->ntimestep % fix->peratom_freq)
1839 print_var_error(FLERR,"Fix in variable not computed at "
1840 "compatible time",ivar);
1841
1842 Tree *newtree = new Tree();
1843 newtree->type = ATOMARRAY;
1844 if (fix->array_atom)
1845 newtree->array = &fix->array_atom[0][index1-1];
1846 newtree->nstride = fix->size_peratom_cols;
1847 treestack[ntreestack++] = newtree;
1848
1849 } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar);
1850
1851 // ----------------
1852 // variable
1853 // ----------------
1854
1855 } else if (strncmp(word,"v_",2) == 0) {
1856
1857 int ivar = find(word+2);
1858 if (ivar < 0)
1859 print_var_error(FLERR,fmt::format("Invalid variable reference "
1860 "{} in variable formula",word),ivar);
1861 if (eval_in_progress[ivar])
1862 print_var_error(FLERR,"has a circular dependency",ivar);
1863
1864 // parse zero or one trailing brackets
1865 // point i beyond last bracket
1866 // nbracket = # of bracket pairs
1867 // index = int inside bracket, possibly an atom ID
1868
1869 int nbracket;
1870 tagint index;
1871 if (str[i] != '[') nbracket = 0;
1872 else {
1873 nbracket = 1;
1874 ptr = &str[i];
1875 index = int_between_brackets(ptr,1);
1876 i = ptr-str+1;
1877 }
1878
1879 // v_name = scalar from internal-style variable
1880 // access value directly
1881
1882 if (nbracket == 0 && style[ivar] == INTERNAL) {
1883
1884 value1 = dvalue[ivar];
1885 if (tree) {
1886 Tree *newtree = new Tree();
1887 newtree->type = VALUE;
1888 newtree->value = value1;
1889 treestack[ntreestack++] = newtree;
1890 } else argstack[nargstack++] = value1;
1891
1892 // v_name = scalar from non atom/atomfile & non vector-style variable
1893 // access value via retrieve()
1894
1895 } else if (nbracket == 0 && style[ivar] != ATOM &&
1896 style[ivar] != ATOMFILE && style[ivar] != VECTOR) {
1897
1898 char *var = retrieve(word+2);
1899 if (var == nullptr)
1900 print_var_error(FLERR,"Invalid variable evaluation in "
1901 "variable formula",ivar);
1902 if (tree) {
1903 Tree *newtree = new Tree();
1904 newtree->type = VALUE;
1905 newtree->value = atof(var);
1906 treestack[ntreestack++] = newtree;
1907 } else argstack[nargstack++] = atof(var);
1908
1909 // v_name = per-atom vector from atom-style variable
1910 // evaluate the atom-style variable as newtree
1911
1912 } else if (nbracket == 0 && style[ivar] == ATOM) {
1913
1914 if (tree == nullptr)
1915 print_var_error(FLERR,"Atom-style variable in "
1916 "equal-style variable formula",ivar);
1917 if (treetype == VECTOR)
1918 print_var_error(FLERR,"Atom-style variable in "
1919 "vector-style variable formula",ivar);
1920
1921 Tree *newtree = nullptr;
1922 evaluate(data[ivar][0],&newtree,ivar);
1923 treestack[ntreestack++] = newtree;
1924
1925 // v_name = per-atom vector from atomfile-style variable
1926
1927 } else if (nbracket == 0 && style[ivar] == ATOMFILE) {
1928
1929 if (tree == nullptr)
1930 print_var_error(FLERR,"Atomfile-style variable in "
1931 "equal-style variable formula",ivar);
1932 if (treetype == VECTOR)
1933 print_var_error(FLERR,"Atomfile-style variable in "
1934 "vector-style variable formula",ivar);
1935
1936 Tree *newtree = new Tree();
1937 newtree->type = ATOMARRAY;
1938 newtree->array = reader[ivar]->fixstore->vstore;
1939 newtree->nstride = 1;
1940 treestack[ntreestack++] = newtree;
1941
1942 // v_name = vector from vector-style variable
1943 // evaluate the vector-style variable, put result in newtree
1944
1945 } else if (nbracket == 0 && style[ivar] == VECTOR) {
1946
1947 if (tree == nullptr)
1948 print_var_error(FLERR,"Vector-style variable in "
1949 "equal-style variable formula",ivar);
1950 if (treetype == ATOM)
1951 print_var_error(FLERR,"Vector-style variable in "
1952 "atom-style variable formula",ivar);
1953
1954 double *vec;
1955 int nvec = compute_vector(ivar,&vec);
1956
1957 Tree *newtree = new Tree();
1958 newtree->type = VECTORARRAY;
1959 newtree->array = vec;
1960 newtree->nvector = nvec;
1961 newtree->nstride = 1;
1962 treestack[ntreestack++] = newtree;
1963
1964 // v_name[N] = scalar from atom-style variable
1965 // compute the per-atom variable in result
1966 // use peratom2global to extract single value from result
1967
1968 } else if (nbracket && style[ivar] == ATOM) {
1969
1970 double *result;
1971 memory->create(result,atom->nlocal,"variable:result");
1972 compute_atom(ivar,0,result,1,0);
1973 peratom2global(1,nullptr,result,1,index,
1974 tree,treestack,ntreestack,argstack,nargstack);
1975 memory->destroy(result);
1976
1977 // v_name[N] = scalar from atomfile-style variable
1978
1979 } else if (nbracket && style[ivar] == ATOMFILE) {
1980
1981 peratom2global(1,nullptr,reader[ivar]->fixstore->vstore,1,index,
1982 tree,treestack,ntreestack,argstack,nargstack);
1983
1984 // v_name[N] = scalar from vector-style variable
1985 // compute the vector-style variable, extract single value
1986
1987 } else if (nbracket && style[ivar] == VECTOR) {
1988
1989 double *vec;
1990 int nvec = compute_vector(ivar,&vec);
1991 if (index <= 0 || index > nvec)
1992 print_var_error(FLERR,"Invalid index into "
1993 "vector-style variable",ivar);
1994 int m = index; // convert from tagint to int
1995
1996 if (tree) {
1997 Tree *newtree = new Tree();
1998 newtree->type = VALUE;
1999 newtree->value = vec[m-1];
2000 treestack[ntreestack++] = newtree;
2001 } else argstack[nargstack++] = vec[m-1];
2002
2003 } else print_var_error(FLERR,"Mismatched variable in "
2004 "variable formula",ivar);
2005
2006 // ----------------
2007 // math/group/special function or atom value/vector or
2008 // constant or thermo keyword
2009 // ----------------
2010
2011 } else {
2012
2013 // ----------------
2014 // math or group or special function
2015 // ----------------
2016
2017 if (str[i] == '(') {
2018 char *contents = nullptr;
2019 i = find_matching_paren(str,i,contents,ivar);
2020 i++;
2021
2022 if (math_function(word,contents,tree,treestack,ntreestack,
2023 argstack,nargstack,ivar));
2024 else if (group_function(word,contents,tree,treestack,ntreestack,
2025 argstack,nargstack,ivar));
2026 else if (special_function(word,contents,tree,treestack,ntreestack,
2027 argstack,nargstack,ivar));
2028 else print_var_error(FLERR,fmt::format("Invalid math/group/special "
2029 "function '{}()'in variable "
2030 "formula", word),ivar);
2031 delete [] contents;
2032
2033 // ----------------
2034 // atom value
2035 // ----------------
2036
2037 } else if (str[i] == '[') {
2038 if (domain->box_exist == 0)
2039 print_var_error(FLERR,"Variable evaluation before "
2040 "simulation box is defined",ivar);
2041
2042 ptr = &str[i];
2043 tagint id = int_between_brackets(ptr,1);
2044 i = ptr-str+1;
2045
2046 peratom2global(0,word,nullptr,0,id,
2047 tree,treestack,ntreestack,argstack,nargstack);
2048
2049 // ----------------
2050 // atom vector
2051 // ----------------
2052
2053 } else if (is_atom_vector(word)) {
2054 if (domain->box_exist == 0)
2055 print_var_error(FLERR,"Variable evaluation before "
2056 "simulation box is defined",ivar);
2057
2058 atom_vector(word,tree,treestack,ntreestack);
2059
2060 // ----------------
2061 // constant
2062 // ----------------
2063
2064 } else if (constants.find(word) != constants.end()) {
2065 value1 = constants[word];
2066 if (tree) {
2067 Tree *newtree = new Tree();
2068 newtree->type = VALUE;
2069 newtree->value = value1;
2070 treestack[ntreestack++] = newtree;
2071 } else argstack[nargstack++] = value1;
2072
2073 // ----------------
2074 // thermo keyword
2075 // ----------------
2076
2077 } else {
2078 if (domain->box_exist == 0)
2079 print_var_error(FLERR,"Variable evaluation before "
2080 "simulation box is defined",ivar);
2081
2082 int flag = output->thermo->evaluate_keyword(word,&value1);
2083 if (flag)
2084 print_var_error(FLERR,fmt::format("Invalid thermo keyword '{}' in "
2085 "variable formula",word),ivar);
2086 if (tree) {
2087 Tree *newtree = new Tree();
2088 newtree->type = VALUE;
2089 newtree->value = value1;
2090 treestack[ntreestack++] = newtree;
2091 } else argstack[nargstack++] = value1;
2092 }
2093 }
2094
2095 delete [] word;
2096
2097 // ----------------
2098 // math operator, including end-of-string
2099 // ----------------
2100
2101 } else if (strchr("+-*/^<>=!&|%\0",onechar)) {
2102 if (onechar == '+') op = ADD;
2103 else if (onechar == '-') op = SUBTRACT;
2104 else if (onechar == '*') op = MULTIPLY;
2105 else if (onechar == '/') op = DIVIDE;
2106 else if (onechar == '%') op = MODULO;
2107 else if (onechar == '^') op = CARAT;
2108 else if (onechar == '=') {
2109 if (str[i+1] != '=')
2110 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
2111 op = EQ;
2112 i++;
2113 } else if (onechar == '!') {
2114 if (str[i+1] == '=') {
2115 op = NE;
2116 i++;
2117 } else op = NOT;
2118 } else if (onechar == '<') {
2119 if (str[i+1] != '=') op = LT;
2120 else {
2121 op = LE;
2122 i++;
2123 }
2124 } else if (onechar == '>') {
2125 if (str[i+1] != '=') op = GT;
2126 else {
2127 op = GE;
2128 i++;
2129 }
2130 } else if (onechar == '&') {
2131 if (str[i+1] != '&')
2132 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
2133 op = AND;
2134 i++;
2135 } else if (onechar == '|') {
2136 if (str[i+1] == '|') op = OR;
2137 else if (str[i+1] == '^') op = XOR;
2138 else print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
2139 i++;
2140 } else op = DONE;
2141
2142 i++;
2143
2144 if (op == SUBTRACT && expect == ARG) {
2145 opstack[nopstack++] = UNARY;
2146 continue;
2147 }
2148 if (op == NOT && expect == ARG) {
2149 opstack[nopstack++] = op;
2150 continue;
2151 }
2152
2153 if (expect == ARG)
2154 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
2155 expect = ARG;
2156
2157 // evaluate stack as deep as possible while respecting precedence
2158 // before pushing current op onto stack
2159
2160 while (nopstack && precedence[opstack[nopstack-1]] >= precedence[op]) {
2161 opprevious = opstack[--nopstack];
2162
2163 if (tree) {
2164 Tree *newtree = new Tree();
2165 newtree->type = opprevious;
2166 if ((opprevious == UNARY) || (opprevious == NOT)) {
2167 newtree->first = treestack[--ntreestack];
2168 } else {
2169 newtree->second = treestack[--ntreestack];
2170 newtree->first = treestack[--ntreestack];
2171 }
2172 treestack[ntreestack++] = newtree;
2173
2174 } else {
2175 value2 = argstack[--nargstack];
2176 if (opprevious != UNARY && opprevious != NOT)
2177 value1 = argstack[--nargstack];
2178
2179 if (opprevious == ADD)
2180 argstack[nargstack++] = value1 + value2;
2181 else if (opprevious == SUBTRACT)
2182 argstack[nargstack++] = value1 - value2;
2183 else if (opprevious == MULTIPLY)
2184 argstack[nargstack++] = value1 * value2;
2185 else if (opprevious == DIVIDE) {
2186 if (value2 == 0.0)
2187 print_var_error(FLERR,"Divide by 0 in variable formula",ivar,0);
2188 argstack[nargstack++] = value1 / value2;
2189 } else if (opprevious == MODULO) {
2190 if (value2 == 0.0)
2191 print_var_error(FLERR,"Modulo 0 in variable formula",ivar,0);
2192 argstack[nargstack++] = fmod(value1,value2);
2193 } else if (opprevious == CARAT) {
2194 if (value2 == 0.0)
2195 argstack[nargstack++] = 1.0;
2196 else if ((value1 == 0.0) && (value2 < 0.0))
2197 print_var_error(FLERR,"Invalid power expression in "
2198 "variable formula",ivar,0);
2199 else argstack[nargstack++] = pow(value1,value2);
2200 } else if (opprevious == UNARY) {
2201 argstack[nargstack++] = -value2;
2202 } else if (opprevious == NOT) {
2203 if (value2 == 0.0) argstack[nargstack++] = 1.0;
2204 else argstack[nargstack++] = 0.0;
2205 } else if (opprevious == EQ) {
2206 if (value1 == value2) argstack[nargstack++] = 1.0;
2207 else argstack[nargstack++] = 0.0;
2208 } else if (opprevious == NE) {
2209 if (value1 != value2) argstack[nargstack++] = 1.0;
2210 else argstack[nargstack++] = 0.0;
2211 } else if (opprevious == LT) {
2212 if (value1 < value2) argstack[nargstack++] = 1.0;
2213 else argstack[nargstack++] = 0.0;
2214 } else if (opprevious == LE) {
2215 if (value1 <= value2) argstack[nargstack++] = 1.0;
2216 else argstack[nargstack++] = 0.0;
2217 } else if (opprevious == GT) {
2218 if (value1 > value2) argstack[nargstack++] = 1.0;
2219 else argstack[nargstack++] = 0.0;
2220 } else if (opprevious == GE) {
2221 if (value1 >= value2) argstack[nargstack++] = 1.0;
2222 else argstack[nargstack++] = 0.0;
2223 } else if (opprevious == AND) {
2224 if (value1 != 0.0 && value2 != 0.0) argstack[nargstack++] = 1.0;
2225 else argstack[nargstack++] = 0.0;
2226 } else if (opprevious == OR) {
2227 if (value1 != 0.0 || value2 != 0.0) argstack[nargstack++] = 1.0;
2228 else argstack[nargstack++] = 0.0;
2229 } else if (opprevious == XOR) {
2230 if ((value1 == 0.0 && value2 != 0.0) ||
2231 (value1 != 0.0 && value2 == 0.0)) argstack[nargstack++] = 1.0;
2232 else argstack[nargstack++] = 0.0;
2233 }
2234 }
2235 }
2236
2237 // if end-of-string, break out of entire formula evaluation loop
2238
2239 if (op == DONE) break;
2240
2241 // push current operation onto stack
2242
2243 opstack[nopstack++] = op;
2244
2245 } else print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
2246 }
2247
2248 if (nopstack) print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
2249
2250 // for atom-style variable, return remaining tree
2251 // for equal-style variable, return remaining arg
2252
2253 if (tree) {
2254 if (ntreestack != 1)
2255 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
2256 *tree = treestack[0];
2257 return 0.0;
2258 } else {
2259 if (nargstack != 1)
2260 print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
2261 return argstack[0];
2262 }
2263 }
2264
2265 /* ----------------------------------------------------------------------
2266 one-time collapse of an atom-style variable parse tree
2267 tree was created by one-time parsing of formula string via evaluate()
2268 only keep tree nodes that depend on
2269 ATOMARRAY, TYPEARRAY, INTARRAY, BIGINTARRAY, VECTOR
2270 remainder is converted to single VALUE
2271 this enables optimal eval_tree loop over atoms
2272 customize by adding a function:
2273 sqrt(),exp(),ln(),log(),abs(),sin(),cos(),tan(),asin(),acos(),atan(),
2274 atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
2275 ramp(x,y),stagger(x,y),logfreq(x,y,z),logfreq2(x,y,z),
2276 logfreq3(x,y,z),stride(x,y,z),vdisplace(x,y),swiggle(x,y,z),
2277 cwiggle(x,y,z),gmask(x),rmask(x),grmask(x,y)
2278 ---------------------------------------------------------------------- */
2279
collapse_tree(Tree * tree)2280 double Variable::collapse_tree(Tree *tree)
2281 {
2282 double arg1,arg2,arg3;
2283
2284 if (tree->type == VALUE) return tree->value;
2285 if (tree->type == ATOMARRAY) return 0.0;
2286 if (tree->type == TYPEARRAY) return 0.0;
2287 if (tree->type == INTARRAY) return 0.0;
2288 if (tree->type == BIGINTARRAY) return 0.0;
2289 if (tree->type == VECTORARRAY) return 0.0;
2290
2291 if (tree->type == ADD) {
2292 arg1 = collapse_tree(tree->first);
2293 arg2 = collapse_tree(tree->second);
2294 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2295 tree->type = VALUE;
2296 tree->value = arg1 + arg2;
2297 return tree->value;
2298 }
2299
2300 if (tree->type == SUBTRACT) {
2301 arg1 = collapse_tree(tree->first);
2302 arg2 = collapse_tree(tree->second);
2303 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2304 tree->type = VALUE;
2305 tree->value = arg1 - arg2;
2306 return tree->value;
2307 }
2308
2309 if (tree->type == MULTIPLY) {
2310 arg1 = collapse_tree(tree->first);
2311 arg2 = collapse_tree(tree->second);
2312 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2313 tree->type = VALUE;
2314 tree->value = arg1 * arg2;
2315 return tree->value;
2316 }
2317
2318 if (tree->type == DIVIDE) {
2319 arg1 = collapse_tree(tree->first);
2320 arg2 = collapse_tree(tree->second);
2321 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2322 tree->type = VALUE;
2323 if (arg2 == 0.0) error->one(FLERR,"Divide by 0 in variable formula");
2324 tree->value = arg1 / arg2;
2325 return tree->value;
2326 }
2327
2328 if (tree->type == MODULO) {
2329 arg1 = collapse_tree(tree->first);
2330 arg2 = collapse_tree(tree->second);
2331 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2332 tree->type = VALUE;
2333 if (arg2 == 0.0) error->one(FLERR,"Modulo 0 in variable formula");
2334 tree->value = fmod(arg1,arg2);
2335 return tree->value;
2336 }
2337
2338 if (tree->type == CARAT) {
2339 arg1 = collapse_tree(tree->first);
2340 arg2 = collapse_tree(tree->second);
2341 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2342 tree->type = VALUE;
2343 if (arg2 == 0.0) error->one(FLERR,"Power by 0 in variable formula");
2344 tree->value = pow(arg1,arg2);
2345 return tree->value;
2346 }
2347
2348 if (tree->type == UNARY) {
2349 arg1 = collapse_tree(tree->first);
2350 if (tree->first->type != VALUE) return 0.0;
2351 tree->type = VALUE;
2352 tree->value = -arg1;
2353 return tree->value;
2354 }
2355
2356 if (tree->type == NOT) {
2357 arg1 = collapse_tree(tree->first);
2358 if (tree->first->type != VALUE) return 0.0;
2359 tree->type = VALUE;
2360 if (arg1 == 0.0) tree->value = 1.0;
2361 else tree->value = 0.0;
2362 return tree->value;
2363 }
2364
2365 if (tree->type == EQ) {
2366 arg1 = collapse_tree(tree->first);
2367 arg2 = collapse_tree(tree->second);
2368 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2369 tree->type = VALUE;
2370 if (arg1 == arg2) tree->value = 1.0;
2371 else tree->value = 0.0;
2372 return tree->value;
2373 }
2374
2375 if (tree->type == NE) {
2376 arg1 = collapse_tree(tree->first);
2377 arg2 = collapse_tree(tree->second);
2378 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2379 tree->type = VALUE;
2380 if (arg1 != arg2) tree->value = 1.0;
2381 else tree->value = 0.0;
2382 return tree->value;
2383 }
2384
2385 if (tree->type == LT) {
2386 arg1 = collapse_tree(tree->first);
2387 arg2 = collapse_tree(tree->second);
2388 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2389 tree->type = VALUE;
2390 if (arg1 < arg2) tree->value = 1.0;
2391 else tree->value = 0.0;
2392 return tree->value;
2393 }
2394
2395 if (tree->type == LE) {
2396 arg1 = collapse_tree(tree->first);
2397 arg2 = collapse_tree(tree->second);
2398 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2399 tree->type = VALUE;
2400 if (arg1 <= arg2) tree->value = 1.0;
2401 else tree->value = 0.0;
2402 return tree->value;
2403 }
2404
2405 if (tree->type == GT) {
2406 arg1 = collapse_tree(tree->first);
2407 arg2 = collapse_tree(tree->second);
2408 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2409 tree->type = VALUE;
2410 if (arg1 > arg2) tree->value = 1.0;
2411 else tree->value = 0.0;
2412 return tree->value;
2413 }
2414
2415 if (tree->type == GE) {
2416 arg1 = collapse_tree(tree->first);
2417 arg2 = collapse_tree(tree->second);
2418 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2419 tree->type = VALUE;
2420 if (arg1 >= arg2) tree->value = 1.0;
2421 else tree->value = 0.0;
2422 return tree->value;
2423 }
2424
2425 if (tree->type == AND) {
2426 arg1 = collapse_tree(tree->first);
2427 arg2 = collapse_tree(tree->second);
2428 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2429 tree->type = VALUE;
2430 if (arg1 != 0.0 && arg2 != 0.0) tree->value = 1.0;
2431 else tree->value = 0.0;
2432 return tree->value;
2433 }
2434
2435 if (tree->type == OR) {
2436 arg1 = collapse_tree(tree->first);
2437 arg2 = collapse_tree(tree->second);
2438 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2439 tree->type = VALUE;
2440 if (arg1 != 0.0 || arg2 != 0.0) tree->value = 1.0;
2441 else tree->value = 0.0;
2442 return tree->value;
2443 }
2444
2445 if (tree->type == XOR) {
2446 arg1 = collapse_tree(tree->first);
2447 arg2 = collapse_tree(tree->second);
2448 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2449 tree->type = VALUE;
2450 if ((arg1 == 0.0 && arg2 != 0.0) || (arg1 != 0.0 && arg2 == 0.0))
2451 tree->value = 1.0;
2452 else tree->value = 0.0;
2453 return tree->value;
2454 }
2455
2456 if (tree->type == SQRT) {
2457 arg1 = collapse_tree(tree->first);
2458 if (tree->first->type != VALUE) return 0.0;
2459 tree->type = VALUE;
2460 if (arg1 < 0.0)
2461 error->one(FLERR,"Sqrt of negative value in variable formula");
2462 tree->value = sqrt(arg1);
2463 return tree->value;
2464 }
2465
2466 if (tree->type == EXP) {
2467 arg1 = collapse_tree(tree->first);
2468 if (tree->first->type != VALUE) return 0.0;
2469 tree->type = VALUE;
2470 tree->value = exp(arg1);
2471 return tree->value;
2472 }
2473
2474 if (tree->type == LN) {
2475 arg1 = collapse_tree(tree->first);
2476 if (tree->first->type != VALUE) return 0.0;
2477 tree->type = VALUE;
2478 if (arg1 <= 0.0)
2479 error->one(FLERR,"Log of zero/negative value in variable formula");
2480 tree->value = log(arg1);
2481 return tree->value;
2482 }
2483
2484 if (tree->type == LOG) {
2485 arg1 = collapse_tree(tree->first);
2486 if (tree->first->type != VALUE) return 0.0;
2487 tree->type = VALUE;
2488 if (arg1 <= 0.0)
2489 error->one(FLERR,"Log of zero/negative value in variable formula");
2490 tree->value = log10(arg1);
2491 return tree->value;
2492 }
2493
2494 if (tree->type == ABS) {
2495 arg1 = collapse_tree(tree->first);
2496 if (tree->first->type != VALUE) return 0.0;
2497 tree->type = VALUE;
2498 tree->value = fabs(arg1);
2499 return tree->value;
2500 }
2501
2502 if (tree->type == SIN) {
2503 arg1 = collapse_tree(tree->first);
2504 if (tree->first->type != VALUE) return 0.0;
2505 tree->type = VALUE;
2506 tree->value = sin(arg1);
2507 return tree->value;
2508 }
2509
2510 if (tree->type == COS) {
2511 arg1 = collapse_tree(tree->first);
2512 if (tree->first->type != VALUE) return 0.0;
2513 tree->type = VALUE;
2514 tree->value = cos(arg1);
2515 return tree->value;
2516 }
2517
2518 if (tree->type == TAN) {
2519 arg1 = collapse_tree(tree->first);
2520 if (tree->first->type != VALUE) return 0.0;
2521 tree->type = VALUE;
2522 tree->value = tan(arg1);
2523 return tree->value;
2524 }
2525
2526 if (tree->type == ASIN) {
2527 arg1 = collapse_tree(tree->first);
2528 if (tree->first->type != VALUE) return 0.0;
2529 tree->type = VALUE;
2530 if (arg1 < -1.0 || arg1 > 1.0)
2531 error->one(FLERR,"Arcsin of invalid value in variable formula");
2532 tree->value = asin(arg1);
2533 return tree->value;
2534 }
2535
2536 if (tree->type == ACOS) {
2537 arg1 = collapse_tree(tree->first);
2538 if (tree->first->type != VALUE) return 0.0;
2539 tree->type = VALUE;
2540 if (arg1 < -1.0 || arg1 > 1.0)
2541 error->one(FLERR,"Arccos of invalid value in variable formula");
2542 tree->value = acos(arg1);
2543 return tree->value;
2544 }
2545
2546 if (tree->type == ATAN) {
2547 arg1 = collapse_tree(tree->first);
2548 if (tree->first->type != VALUE) return 0.0;
2549 tree->type = VALUE;
2550 tree->value = atan(arg1);
2551 return tree->value;
2552 }
2553
2554 if (tree->type == ATAN2) {
2555 arg1 = collapse_tree(tree->first);
2556 arg2 = collapse_tree(tree->second);
2557 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2558 tree->type = VALUE;
2559 tree->value = atan2(arg1,arg2);
2560 return tree->value;
2561 }
2562
2563 // random() or normal() do not become a single collapsed value
2564
2565 if (tree->type == RANDOM) {
2566 collapse_tree(tree->first);
2567 collapse_tree(tree->second);
2568 if (randomatom == nullptr) {
2569 int seed = static_cast<int> (collapse_tree(tree->extra[0]));
2570 if (seed <= 0)
2571 error->one(FLERR,"Invalid math function in variable formula");
2572 randomatom = new RanMars(lmp,seed+me);
2573 }
2574 return 0.0;
2575 }
2576
2577 if (tree->type == NORMAL) {
2578 collapse_tree(tree->first);
2579 double sigma = collapse_tree(tree->second);
2580 if (sigma < 0.0)
2581 error->one(FLERR,"Invalid math function in variable formula");
2582 if (randomatom == nullptr) {
2583 int seed = static_cast<int> (collapse_tree(tree->extra[0]));
2584 if (seed <= 0)
2585 error->one(FLERR,"Invalid math function in variable formula");
2586 randomatom = new RanMars(lmp,seed+me);
2587 }
2588 return 0.0;
2589 }
2590
2591 if (tree->type == CEIL) {
2592 arg1 = collapse_tree(tree->first);
2593 if (tree->first->type != VALUE) return 0.0;
2594 tree->type = VALUE;
2595 tree->value = ceil(arg1);
2596 return tree->value;
2597 }
2598
2599 if (tree->type == FLOOR) {
2600 arg1 = collapse_tree(tree->first);
2601 if (tree->first->type != VALUE) return 0.0;
2602 tree->type = VALUE;
2603 tree->value = floor(arg1);
2604 return tree->value;
2605 }
2606
2607 if (tree->type == ROUND) {
2608 arg1 = collapse_tree(tree->first);
2609 if (tree->first->type != VALUE) return 0.0;
2610 tree->type = VALUE;
2611 tree->value = MYROUND(arg1);
2612 return tree->value;
2613 }
2614
2615 if (tree->type == RAMP) {
2616 arg1 = collapse_tree(tree->first);
2617 arg2 = collapse_tree(tree->second);
2618 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2619 tree->type = VALUE;
2620 double delta = update->ntimestep - update->beginstep;
2621 if (delta != 0.0) delta /= update->endstep - update->beginstep;
2622 tree->value = arg1 + delta*(arg2-arg1);
2623 return tree->value;
2624 }
2625
2626 if (tree->type == STAGGER) {
2627 bigint ivalue1 = static_cast<bigint> (collapse_tree(tree->first));
2628 bigint ivalue2 = static_cast<bigint> (collapse_tree(tree->second));
2629 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2630 tree->type = VALUE;
2631 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
2632 error->one(FLERR,"Invalid math function in variable formula");
2633 bigint lower = update->ntimestep/ivalue1 * ivalue1;
2634 bigint delta = update->ntimestep - lower;
2635 if (delta < ivalue2) tree->value = lower+ivalue2;
2636 else tree->value = lower+ivalue1;
2637 return tree->value;
2638 }
2639
2640 if (tree->type == LOGFREQ) {
2641 bigint ivalue1 = static_cast<bigint> (collapse_tree(tree->first));
2642 bigint ivalue2 = static_cast<bigint> (collapse_tree(tree->second));
2643 bigint ivalue3 = static_cast<bigint> (collapse_tree(tree->extra[0]));
2644 if (tree->first->type != VALUE || tree->second->type != VALUE ||
2645 tree->extra[0]->type != VALUE) return 0.0;
2646 tree->type = VALUE;
2647 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
2648 error->one(FLERR,"Invalid math function in variable formula");
2649 if (update->ntimestep < ivalue1) tree->value = ivalue1;
2650 else {
2651 bigint lower = ivalue1;
2652 while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
2653 bigint multiple = update->ntimestep/lower;
2654 if (multiple < ivalue2) tree->value = (multiple+1)*lower;
2655 else tree->value = lower*ivalue3;
2656 }
2657 return tree->value;
2658 }
2659
2660 if (tree->type == LOGFREQ2) {
2661 bigint ivalue1 = static_cast<bigint> (collapse_tree(tree->first));
2662 bigint ivalue2 = static_cast<bigint> (collapse_tree(tree->second));
2663 bigint ivalue3 = static_cast<bigint> (collapse_tree(tree->extra[0]));
2664 if (tree->first->type != VALUE || tree->second->type != VALUE ||
2665 tree->extra[0]->type != VALUE) return 0.0;
2666 tree->type = VALUE;
2667 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 )
2668 error->all(FLERR,"Invalid math function in variable formula");
2669 if (update->ntimestep < ivalue1) tree->value = ivalue1;
2670 else {
2671 tree->value = ivalue1;
2672 double delta = ivalue1*(ivalue3-1.0)/ivalue2;
2673 bigint count = 0;
2674 while (update->ntimestep >= tree->value) {
2675 tree->value += delta;
2676 count++;
2677 if (count % ivalue2 == 0) delta *= ivalue3;
2678 }
2679 }
2680 tree->value = ceil(tree->value);
2681 return tree->value;
2682 }
2683
2684 if (tree->type == LOGFREQ3) {
2685 bigint ivalue1 = static_cast<bigint> (collapse_tree(tree->first));
2686 bigint ivalue2 = static_cast<bigint> (collapse_tree(tree->second));
2687 bigint ivalue3 = static_cast<bigint> (collapse_tree(tree->extra[0]));
2688 if (tree->first->type != VALUE || tree->second->type != VALUE ||
2689 tree->extra[0]->type != VALUE) return 0.0;
2690 tree->type = VALUE;
2691 if (ivalue1 <= 0 || ivalue2 <= 1 || ivalue3 <= 0 ||
2692 ivalue3-ivalue1+1 < ivalue2 )
2693 error->all(FLERR,"Invalid math function in variable formula");
2694 if (update->ntimestep < ivalue1) tree->value = ivalue1;
2695 //else if (update->ntimestep <= ivalue3) {
2696 else {
2697 tree->value = ivalue1;
2698 double logsp = ivalue1;
2699 double factor = pow(((double)ivalue3)/ivalue1, 1.0/(ivalue2-1));
2700 bigint linsp = ivalue1;
2701 while (update->ntimestep >= (tree->value)) {
2702 logsp *= factor;
2703 linsp++;
2704 if (linsp > logsp) tree->value = linsp;
2705 else tree->value = ceil(logsp)-(((int)ceil(logsp)-1)/ivalue3);
2706 }
2707 }
2708 if (update->ntimestep > ivalue3)
2709 error->all(FLERR,"Calls to variable exceeded limit");
2710 return tree->value;
2711 }
2712
2713 if (tree->type == STRIDE) {
2714 bigint ivalue1 = static_cast<bigint> (collapse_tree(tree->first));
2715 bigint ivalue2 = static_cast<bigint> (collapse_tree(tree->second));
2716 bigint ivalue3 = static_cast<bigint> (collapse_tree(tree->extra[0]));
2717 if (tree->first->type != VALUE || tree->second->type != VALUE ||
2718 tree->extra[0]->type != VALUE) return 0.0;
2719 tree->type = VALUE;
2720 if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
2721 error->one(FLERR,"Invalid math function in variable formula");
2722 if (update->ntimestep < ivalue1) tree->value = ivalue1;
2723 else if (update->ntimestep < ivalue2) {
2724 bigint offset = update->ntimestep - ivalue1;
2725 tree->value = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
2726 if (tree->value > ivalue2) tree->value = (double) MAXBIGINT;
2727 } else tree->value = (double) MAXBIGINT;
2728 return tree->value;
2729 }
2730
2731 if (tree->type == STRIDE2) {
2732 bigint ivalue1 = static_cast<bigint> (collapse_tree(tree->first));
2733 bigint ivalue2 = static_cast<bigint> (collapse_tree(tree->second));
2734 bigint ivalue3 = static_cast<bigint> (collapse_tree(tree->extra[0]));
2735 bigint ivalue4 = static_cast<bigint> (collapse_tree(tree->extra[1]));
2736 bigint ivalue5 = static_cast<bigint> (collapse_tree(tree->extra[2]));
2737 bigint ivalue6 = static_cast<bigint> (collapse_tree(tree->extra[3]));
2738 if (tree->first->type != VALUE || tree->second->type != VALUE ||
2739 tree->extra[0]->type != VALUE || tree->extra[1]->type != VALUE ||
2740 tree->extra[2]->type != VALUE || tree->extra[3]->type != VALUE)
2741 return 0.0;
2742 tree->type = VALUE;
2743 if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
2744 error->one(FLERR,"Invalid math function in variable formula");
2745 if (ivalue4 < 0 || ivalue5 < 0 || ivalue6 <= 0 || ivalue4 > ivalue5)
2746 error->one(FLERR,"Invalid math function in variable formula");
2747 if (ivalue4 < ivalue1 || ivalue5 > ivalue2)
2748 error->one(FLERR,"Invalid math function in variable formula");
2749 bigint istep, offset;
2750 if (update->ntimestep < ivalue1) istep = ivalue1;
2751 else if (update->ntimestep < ivalue2) {
2752 if (update->ntimestep < ivalue4 || update->ntimestep > ivalue5) {
2753 offset = update->ntimestep - ivalue1;
2754 istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
2755 if (update->ntimestep < ivalue2 && istep > ivalue4)
2756 tree->value = ivalue4;
2757 } else {
2758 offset = update->ntimestep - ivalue4;
2759 istep = ivalue4 + (offset/ivalue6)*ivalue6 + ivalue6;
2760 if (istep > ivalue5) {
2761 offset = ivalue5 - ivalue1;
2762 istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
2763 if (istep > ivalue2) istep = MAXBIGINT;
2764 }
2765 }
2766 } else istep = MAXBIGINT;
2767 tree->value = istep;
2768 return tree->value;
2769 }
2770
2771 if (tree->type == VDISPLACE) {
2772 arg1 = collapse_tree(tree->first);
2773 arg2 = collapse_tree(tree->second);
2774 if (tree->first->type != VALUE || tree->second->type != VALUE) return 0.0;
2775 tree->type = VALUE;
2776 double delta = update->ntimestep - update->beginstep;
2777 tree->value = arg1 + arg2*delta*update->dt;
2778 return tree->value;
2779 }
2780
2781 if (tree->type == SWIGGLE) {
2782 arg1 = collapse_tree(tree->first);
2783 arg2 = collapse_tree(tree->second);
2784 arg3 = collapse_tree(tree->extra[0]);
2785 if (tree->first->type != VALUE || tree->second->type != VALUE ||
2786 tree->extra[0]->type != VALUE) return 0.0;
2787 tree->type = VALUE;
2788 if (arg3 == 0.0)
2789 error->one(FLERR,"Invalid math function in variable formula");
2790 double delta = update->ntimestep - update->beginstep;
2791 double omega = 2.0*MY_PI/arg3;
2792 tree->value = arg1 + arg2*sin(omega*delta*update->dt);
2793 return tree->value;
2794 }
2795
2796 if (tree->type == CWIGGLE) {
2797 arg1 = collapse_tree(tree->first);
2798 arg2 = collapse_tree(tree->second);
2799 arg3 = collapse_tree(tree->extra[0]);
2800 if (tree->first->type != VALUE || tree->second->type != VALUE ||
2801 tree->extra[0]->type != VALUE) return 0.0;
2802 tree->type = VALUE;
2803 if (arg3 == 0.0)
2804 error->one(FLERR,"Invalid math function in variable formula");
2805 double delta = update->ntimestep - update->beginstep;
2806 double omega = 2.0*MY_PI/arg3;
2807 tree->value = arg1 + arg2*(1.0-cos(omega*delta*update->dt));
2808 return tree->value;
2809 }
2810
2811 // mask functions do not become a single collapsed value
2812
2813 if (tree->type == GMASK) return 0.0;
2814 if (tree->type == RMASK) return 0.0;
2815 if (tree->type == GRMASK) return 0.0;
2816
2817 return 0.0;
2818 }
2819
2820 /* ----------------------------------------------------------------------
2821 evaluate an atom-style or vector-style variable parse tree
2822 index I = atom I or vector index I
2823 tree was created by one-time parsing of formula string via evaluate()
2824 customize by adding a function:
2825 sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan(),
2826 atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
2827 ramp(x,y),stagger(x,y),logfreq(x,y,z),logfreq2(x,y,z),
2828 logfreq3(x,y,z),stride(x,y,z),stride2(x,y,z),vdisplace(x,y),
2829 swiggle(x,y,z),cwiggle(x,y,z),gmask(x),rmask(x),grmask(x,y)
2830 ---------------------------------------------------------------------- */
2831
eval_tree(Tree * tree,int i)2832 double Variable::eval_tree(Tree *tree, int i)
2833 {
2834 double arg,arg1,arg2,arg3;
2835
2836 if (tree->type == VALUE) return tree->value;
2837 if (tree->type == ATOMARRAY) return tree->array[i*tree->nstride];
2838 if (tree->type == TYPEARRAY) return tree->array[atom->type[i]];
2839 if (tree->type == INTARRAY) return (double) tree->iarray[i*tree->nstride];
2840 if (tree->type == BIGINTARRAY) return (double) tree->barray[i*tree->nstride];
2841 if (tree->type == VECTORARRAY) return tree->array[i*tree->nstride];
2842
2843 if (tree->type == ADD)
2844 return eval_tree(tree->first,i) + eval_tree(tree->second,i);
2845 if (tree->type == SUBTRACT)
2846 return eval_tree(tree->first,i) - eval_tree(tree->second,i);
2847 if (tree->type == MULTIPLY)
2848 return eval_tree(tree->first,i) * eval_tree(tree->second,i);
2849 if (tree->type == DIVIDE) {
2850 double denom = eval_tree(tree->second,i);
2851 if (denom == 0.0) error->one(FLERR,"Divide by 0 in variable formula");
2852 return eval_tree(tree->first,i) / denom;
2853 }
2854 if (tree->type == MODULO) {
2855 double denom = eval_tree(tree->second,i);
2856 if (denom == 0.0) error->one(FLERR,"Modulo 0 in variable formula");
2857 return fmod(eval_tree(tree->first,i),denom);
2858 }
2859 if (tree->type == CARAT) {
2860 double exponent = eval_tree(tree->second,i);
2861 if (exponent == 0.0) error->one(FLERR,"Power by 0 in variable formula");
2862 return pow(eval_tree(tree->first,i),exponent);
2863 }
2864 if (tree->type == UNARY) return -eval_tree(tree->first,i);
2865
2866 if (tree->type == NOT) {
2867 if (eval_tree(tree->first,i) == 0.0) return 1.0;
2868 else return 0.0;
2869 }
2870 if (tree->type == EQ) {
2871 if (eval_tree(tree->first,i) == eval_tree(tree->second,i)) return 1.0;
2872 else return 0.0;
2873 }
2874 if (tree->type == NE) {
2875 if (eval_tree(tree->first,i) != eval_tree(tree->second,i)) return 1.0;
2876 else return 0.0;
2877 }
2878 if (tree->type == LT) {
2879 if (eval_tree(tree->first,i) < eval_tree(tree->second,i)) return 1.0;
2880 else return 0.0;
2881 }
2882 if (tree->type == LE) {
2883 if (eval_tree(tree->first,i) <= eval_tree(tree->second,i)) return 1.0;
2884 else return 0.0;
2885 }
2886 if (tree->type == GT) {
2887 if (eval_tree(tree->first,i) > eval_tree(tree->second,i)) return 1.0;
2888 else return 0.0;
2889 }
2890 if (tree->type == GE) {
2891 if (eval_tree(tree->first,i) >= eval_tree(tree->second,i)) return 1.0;
2892 else return 0.0;
2893 }
2894 if (tree->type == AND) {
2895 if (eval_tree(tree->first,i) != 0.0 && eval_tree(tree->second,i) != 0.0)
2896 return 1.0;
2897 else return 0.0;
2898 }
2899 if (tree->type == OR) {
2900 if (eval_tree(tree->first,i) != 0.0 || eval_tree(tree->second,i) != 0.0)
2901 return 1.0;
2902 else return 0.0;
2903 }
2904 if (tree->type == XOR) {
2905 if ((eval_tree(tree->first,i) == 0.0 && eval_tree(tree->second,i) != 0.0)
2906 ||
2907 (eval_tree(tree->first,i) != 0.0 && eval_tree(tree->second,i) == 0.0))
2908 return 1.0;
2909 else return 0.0;
2910 }
2911
2912 if (tree->type == SQRT) {
2913 arg1 = eval_tree(tree->first,i);
2914 if (arg1 < 0.0)
2915 error->one(FLERR,"Sqrt of negative value in variable formula");
2916 return sqrt(arg1);
2917 }
2918 if (tree->type == EXP)
2919 return exp(eval_tree(tree->first,i));
2920 if (tree->type == LN) {
2921 arg1 = eval_tree(tree->first,i);
2922 if (arg1 <= 0.0)
2923 error->one(FLERR,"Log of zero/negative value in variable formula");
2924 return log(arg1);
2925 }
2926 if (tree->type == LOG) {
2927 arg1 = eval_tree(tree->first,i);
2928 if (arg1 <= 0.0)
2929 error->one(FLERR,"Log of zero/negative value in variable formula");
2930 return log10(arg1);
2931 }
2932 if (tree->type == ABS)
2933 return fabs(eval_tree(tree->first,i));
2934
2935 if (tree->type == SIN)
2936 return sin(eval_tree(tree->first,i));
2937 if (tree->type == COS)
2938 return cos(eval_tree(tree->first,i));
2939 if (tree->type == TAN)
2940 return tan(eval_tree(tree->first,i));
2941
2942 if (tree->type == ASIN) {
2943 arg1 = eval_tree(tree->first,i);
2944 if (arg1 < -1.0 || arg1 > 1.0)
2945 error->one(FLERR,"Arcsin of invalid value in variable formula");
2946 return asin(arg1);
2947 }
2948 if (tree->type == ACOS) {
2949 arg1 = eval_tree(tree->first,i);
2950 if (arg1 < -1.0 || arg1 > 1.0)
2951 error->one(FLERR,"Arccos of invalid value in variable formula");
2952 return acos(arg1);
2953 }
2954 if (tree->type == ATAN)
2955 return atan(eval_tree(tree->first,i));
2956 if (tree->type == ATAN2)
2957 return atan2(eval_tree(tree->first,i),eval_tree(tree->second,i));
2958
2959 if (tree->type == RANDOM) {
2960 double lower = eval_tree(tree->first,i);
2961 double upper = eval_tree(tree->second,i);
2962 if (randomatom == nullptr) {
2963 int seed = static_cast<int> (eval_tree(tree->extra[0],i));
2964 if (seed <= 0)
2965 error->one(FLERR,"Invalid math function in variable formula");
2966 randomatom = new RanMars(lmp,seed+me);
2967 }
2968 return randomatom->uniform()*(upper-lower)+lower;
2969 }
2970 if (tree->type == NORMAL) {
2971 double mu = eval_tree(tree->first,i);
2972 double sigma = eval_tree(tree->second,i);
2973 if (sigma < 0.0)
2974 error->one(FLERR,"Invalid math function in variable formula");
2975 if (randomatom == nullptr) {
2976 int seed = static_cast<int> (eval_tree(tree->extra[0],i));
2977 if (seed <= 0)
2978 error->one(FLERR,"Invalid math function in variable formula");
2979 randomatom = new RanMars(lmp,seed+me);
2980 }
2981 return mu + sigma*randomatom->gaussian();
2982 }
2983
2984 if (tree->type == CEIL)
2985 return ceil(eval_tree(tree->first,i));
2986 if (tree->type == FLOOR)
2987 return floor(eval_tree(tree->first,i));
2988 if (tree->type == ROUND)
2989 return MYROUND(eval_tree(tree->first,i));
2990
2991 if (tree->type == RAMP) {
2992 arg1 = eval_tree(tree->first,i);
2993 arg2 = eval_tree(tree->second,i);
2994 double delta = update->ntimestep - update->beginstep;
2995 if (delta != 0.0) delta /= update->endstep - update->beginstep;
2996 arg = arg1 + delta*(arg2-arg1);
2997 return arg;
2998 }
2999
3000 if (tree->type == STAGGER) {
3001 bigint ivalue1 = static_cast<bigint> (eval_tree(tree->first,i));
3002 bigint ivalue2 = static_cast<bigint> (eval_tree(tree->second,i));
3003 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
3004 error->one(FLERR,"Invalid math function in variable formula");
3005 bigint lower = update->ntimestep/ivalue1 * ivalue1;
3006 bigint delta = update->ntimestep - lower;
3007 if (delta < ivalue2) arg = lower+ivalue2;
3008 else arg = lower+ivalue1;
3009 return arg;
3010 }
3011
3012 if (tree->type == LOGFREQ) {
3013 bigint ivalue1 = static_cast<bigint> (eval_tree(tree->first,i));
3014 bigint ivalue2 = static_cast<bigint> (eval_tree(tree->second,i));
3015 bigint ivalue3 = static_cast<bigint> (eval_tree(tree->extra[0],i));
3016 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
3017 error->one(FLERR,"Invalid math function in variable formula");
3018 if (update->ntimestep < ivalue1) arg = ivalue1;
3019 else {
3020 bigint lower = ivalue1;
3021 while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
3022 bigint multiple = update->ntimestep/lower;
3023 if (multiple < ivalue2) arg = (multiple+1)*lower;
3024 else arg = lower*ivalue3;
3025 }
3026 return arg;
3027 }
3028
3029 if (tree->type == LOGFREQ2) {
3030 bigint ivalue1 = static_cast<bigint> (eval_tree(tree->first,i));
3031 bigint ivalue2 = static_cast<bigint> (eval_tree(tree->second,i));
3032 bigint ivalue3 = static_cast<bigint> (eval_tree(tree->extra[0],i));
3033 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 )
3034 error->all(FLERR,"Invalid math function in variable formula");
3035 if (update->ntimestep < ivalue1) arg = ivalue1;
3036 else {
3037 arg = ivalue1;
3038 double delta = ivalue1*(ivalue3-1.0)/ivalue2;
3039 bigint count = 0;
3040 while (update->ntimestep >= arg) {
3041 arg += delta;
3042 count++;
3043 if (count % ivalue2 == 0) delta *= ivalue3;
3044 }
3045 }
3046 arg = ceil(arg);
3047 return arg;
3048 }
3049
3050 if (tree->type == STRIDE) {
3051 bigint ivalue1 = static_cast<bigint> (eval_tree(tree->first,i));
3052 bigint ivalue2 = static_cast<bigint> (eval_tree(tree->second,i));
3053 bigint ivalue3 = static_cast<bigint> (eval_tree(tree->extra[0],i));
3054 if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
3055 error->one(FLERR,"Invalid math function in variable formula");
3056 if (update->ntimestep < ivalue1) arg = ivalue1;
3057 else if (update->ntimestep < ivalue2) {
3058 bigint offset = update->ntimestep - ivalue1;
3059 arg = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
3060 if (arg > ivalue2) arg = (double) MAXBIGINT;
3061 } else arg = (double) MAXBIGINT;
3062 return arg;
3063 }
3064
3065 if (tree->type == STRIDE2) {
3066 bigint ivalue1 = static_cast<bigint> (eval_tree(tree->first,i));
3067 bigint ivalue2 = static_cast<bigint> (eval_tree(tree->second,i));
3068 bigint ivalue3 = static_cast<bigint> (eval_tree(tree->extra[0],i));
3069 bigint ivalue4 = static_cast<bigint> (eval_tree(tree->extra[1],i));
3070 bigint ivalue5 = static_cast<bigint> (eval_tree(tree->extra[2],i));
3071 bigint ivalue6 = static_cast<bigint> (eval_tree(tree->extra[3],i));
3072 if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
3073 error->one(FLERR,"Invalid math function in variable formula");
3074 if (ivalue4 < 0 || ivalue5 < 0 || ivalue6 <= 0 || ivalue4 > ivalue5)
3075 error->one(FLERR,"Invalid math function in variable formula");
3076 if (ivalue4 < ivalue1 || ivalue5 > ivalue2)
3077 error->one(FLERR,"Invalid math function in variable formula");
3078 bigint istep, offset;
3079 if (update->ntimestep < ivalue1) istep = ivalue1;
3080 else if (update->ntimestep < ivalue2) {
3081 if (update->ntimestep < ivalue4 || update->ntimestep > ivalue5) {
3082 offset = update->ntimestep - ivalue1;
3083 istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
3084 if (update->ntimestep < ivalue2 && istep > ivalue4)
3085 tree->value = ivalue4;
3086 } else {
3087 offset = update->ntimestep - ivalue4;
3088 istep = ivalue4 + (offset/ivalue6)*ivalue6 + ivalue6;
3089 if (istep > ivalue5) {
3090 offset = ivalue5 - ivalue1;
3091 istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
3092 if (istep > ivalue2) istep = MAXBIGINT;
3093 }
3094 }
3095 } else istep = MAXBIGINT;
3096 arg = istep;
3097 return arg;
3098 }
3099
3100 if (tree->type == VDISPLACE) {
3101 arg1 = eval_tree(tree->first,i);
3102 arg2 = eval_tree(tree->second,i);
3103 double delta = update->ntimestep - update->beginstep;
3104 arg = arg1 + arg2*delta*update->dt;
3105 return arg;
3106 }
3107
3108 if (tree->type == SWIGGLE) {
3109 arg1 = eval_tree(tree->first,i);
3110 arg2 = eval_tree(tree->second,i);
3111 arg3 = eval_tree(tree->extra[0],i);
3112 if (arg3 == 0.0)
3113 error->one(FLERR,"Invalid math function in variable formula");
3114 double delta = update->ntimestep - update->beginstep;
3115 double omega = 2.0*MY_PI/arg3;
3116 arg = arg1 + arg2*sin(omega*delta*update->dt);
3117 return arg;
3118 }
3119
3120 if (tree->type == CWIGGLE) {
3121 arg1 = eval_tree(tree->first,i);
3122 arg2 = eval_tree(tree->second,i);
3123 arg3 = eval_tree(tree->extra[0],i);
3124 if (arg3 == 0.0)
3125 error->one(FLERR,"Invalid math function in variable formula");
3126 double delta = update->ntimestep - update->beginstep;
3127 double omega = 2.0*MY_PI/arg3;
3128 arg = arg1 + arg2*(1.0-cos(omega*delta*update->dt));
3129 return arg;
3130 }
3131
3132 if (tree->type == GMASK) {
3133 if (atom->mask[i] & tree->ivalue1) return 1.0;
3134 else return 0.0;
3135 }
3136
3137 if (tree->type == RMASK) {
3138 if (domain->regions[tree->ivalue1]->match(atom->x[i][0],
3139 atom->x[i][1],
3140 atom->x[i][2])) return 1.0;
3141 else return 0.0;
3142 }
3143
3144 if (tree->type == GRMASK) {
3145 if ((atom->mask[i] & tree->ivalue1) &&
3146 (domain->regions[tree->ivalue2]->match(atom->x[i][0],
3147 atom->x[i][1],
3148 atom->x[i][2]))) return 1.0;
3149 else return 0.0;
3150 }
3151
3152 return 0.0;
3153 }
3154
3155 /* ----------------------------------------------------------------------
3156 scan entire tree, find size of vectors for vector-style variable
3157 return N for consistent vector size
3158 return 0 for no vector size, caller flags as error
3159 return -1 for inconsistent vector size, caller flags as error
3160 ------------------------------------------------------------------------- */
3161
size_tree_vector(Tree * tree)3162 int Variable::size_tree_vector(Tree *tree)
3163 {
3164 int nsize = 0;
3165 if (tree->type == VECTORARRAY) nsize = tree->nvector;
3166 if (tree->first) nsize = compare_tree_vector(nsize,
3167 size_tree_vector(tree->first));
3168 if (tree->second) nsize = compare_tree_vector(nsize,
3169 size_tree_vector(tree->second));
3170 if (tree->nextra) {
3171 for (int i = 0; i < tree->nextra; i++)
3172 nsize = compare_tree_vector(nsize,size_tree_vector(tree->extra[i]));
3173 }
3174 return nsize;
3175 }
3176
3177 /* ----------------------------------------------------------------------
3178 compare size of two vectors for vector-style variable
3179 return positive size if same or one has no size 0
3180 return -1 error if one is already error or not same positive size
3181 ------------------------------------------------------------------------- */
3182
compare_tree_vector(int i,int j)3183 int Variable::compare_tree_vector(int i, int j)
3184 {
3185 if (i < 0 || j < 0) return -1;
3186 if (i == 0 || j == 0) return MAX(i,j);
3187 if (i != j) return -1;
3188 return i;
3189 }
3190
3191 /* ---------------------------------------------------------------------- */
3192
free_tree(Tree * tree)3193 void Variable::free_tree(Tree *tree)
3194 {
3195 if (tree->first) free_tree(tree->first);
3196 if (tree->second) free_tree(tree->second);
3197 if (tree->nextra) {
3198 for (int i = 0; i < tree->nextra; i++) free_tree(tree->extra[i]);
3199 delete [] tree->extra;
3200 }
3201
3202 if (tree->selfalloc) memory->destroy(tree->array);
3203 delete tree;
3204 }
3205
3206 /* ----------------------------------------------------------------------
3207 find matching parenthesis in str, allocate contents = str between parens
3208 i = left paren
3209 return loc or right paren
3210 ------------------------------------------------------------------------- */
3211
find_matching_paren(char * str,int i,char * & contents,int ivar)3212 int Variable::find_matching_paren(char *str, int i, char *&contents, int ivar)
3213 {
3214 // istop = matching ')' at same level, allowing for nested parens
3215
3216 int istart = i;
3217 int ilevel = 0;
3218 while (1) {
3219 i++;
3220 if (!str[i]) break;
3221 if (str[i] == '(') ilevel++;
3222 else if (str[i] == ')' && ilevel) ilevel--;
3223 else if (str[i] == ')') break;
3224 }
3225 if (!str[i]) print_var_error(FLERR,"Invalid syntax in variable formula",ivar);
3226 int istop = i;
3227
3228 int n = istop - istart - 1;
3229 delete[] contents;
3230 contents = new char[n+1];
3231 strncpy(contents,&str[istart+1],n);
3232 contents[n] = '\0';
3233
3234 return istop;
3235 }
3236
3237 /* ----------------------------------------------------------------------
3238 find int between brackets and return it
3239 return a tagint, since value can be an atom ID
3240 ptr initially points to left bracket
3241 return it pointing to right bracket
3242 error if no right bracket or brackets are empty or index = 0
3243 if varallow = 0: error if any between-bracket chars are non-digits
3244 if varallow = 1: also allow for v_name, where name is variable name
3245 ------------------------------------------------------------------------- */
3246
int_between_brackets(char * & ptr,int varallow)3247 tagint Variable::int_between_brackets(char *&ptr, int varallow)
3248 {
3249 int varflag;
3250 tagint index;
3251
3252 char *start = ++ptr;
3253
3254 if (varallow && utils::strmatch(ptr,"^v_")) {
3255 varflag = 1;
3256 while (*ptr && *ptr != ']') {
3257 if (!isalnum(*ptr) && *ptr != '_')
3258 error->all(FLERR,"Variable name between brackets must be "
3259 "alphanumeric or underscore characters");
3260 ptr++;
3261 }
3262
3263 } else {
3264 varflag = 0;
3265 while (*ptr && *ptr != ']') {
3266 if (!isdigit(*ptr))
3267 error->all(FLERR,"Non digit character between brackets in variable");
3268 ptr++;
3269 }
3270 }
3271
3272 if (*ptr != ']') error->all(FLERR,"Mismatched brackets in variable");
3273 if (ptr == start) error->all(FLERR,"Empty brackets in variable");
3274
3275 *ptr = '\0';
3276
3277 // evaluate index as floating point variable or as tagint via ATOTAGINT()
3278
3279 if (varflag) {
3280 char *id = start+2;
3281 int ivar = find(id);
3282 if (ivar < 0)
3283 error->all(FLERR,"Invalid variable name in variable formula");
3284
3285 char *var = retrieve(id);
3286 if (var == nullptr)
3287 error->all(FLERR,"Invalid variable evaluation in variable formula");
3288 index = static_cast<tagint> (atof(var));
3289
3290 } else index = ATOTAGINT(start);
3291
3292 *ptr = ']';
3293
3294 if (index == 0)
3295 error->all(FLERR,"Index between variable brackets must be positive");
3296 return index;
3297 }
3298
3299 /* ----------------------------------------------------------------------
3300 process a math function in formula
3301 push result onto tree or arg stack
3302 word = math function
3303 contents = str between parentheses with comma-separated args
3304 return 0 if not a match, 1 if successfully processed
3305 customize by adding a math function:
3306 sqrt(),exp(),ln(),log(),abs(),sin(),cos(),tan(),asin(),acos(),atan(),
3307 atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
3308 ramp(x,y),stagger(x,y),logfreq(x,y,z),logfreq2(x,y,z),
3309 logfreq3(x,y,z),stride(x,y,z),stride2(x,y,z,a,b,c),vdisplace(x,y),
3310 swiggle(x,y,z),cwiggle(x,y,z)
3311 ------------------------------------------------------------------------- */
3312
math_function(char * word,char * contents,Tree ** tree,Tree ** treestack,int & ntreestack,double * argstack,int & nargstack,int ivar)3313 int Variable::math_function(char *word, char *contents, Tree **tree,
3314 Tree **treestack, int &ntreestack,
3315 double *argstack, int &nargstack, int ivar)
3316 {
3317 // word not a match to any math function
3318
3319 if (strcmp(word,"sqrt") && strcmp(word,"exp") &&
3320 strcmp(word,"ln") && strcmp(word,"log") &&
3321 strcmp(word,"abs") &&
3322 strcmp(word,"sin") && strcmp(word,"cos") &&
3323 strcmp(word,"tan") && strcmp(word,"asin") &&
3324 strcmp(word,"acos") && strcmp(word,"atan") &&
3325 strcmp(word,"atan2") && strcmp(word,"random") &&
3326 strcmp(word,"normal") && strcmp(word,"ceil") &&
3327 strcmp(word,"floor") && strcmp(word,"round") &&
3328 strcmp(word,"ramp") && strcmp(word,"stagger") &&
3329 strcmp(word,"logfreq") && strcmp(word,"logfreq2") &&
3330 strcmp(word,"logfreq3") && strcmp(word,"stride") &&
3331 strcmp(word,"stride2") && strcmp(word,"vdisplace") &&
3332 strcmp(word,"swiggle") && strcmp(word,"cwiggle"))
3333 return 0;
3334
3335 // parse contents for comma-separated args
3336 // narg = number of args, args = strings between commas
3337
3338 char *args[MAXFUNCARG];
3339 int narg = parse_args(contents,args);
3340
3341 Tree *newtree = nullptr;
3342 double value1,value2;
3343 double values[MAXFUNCARG-2];
3344
3345 if (tree) {
3346 newtree = new Tree();
3347 Tree *argtree = nullptr;
3348 evaluate(args[0],&argtree,ivar);
3349 newtree->first = argtree;
3350 if (narg > 1) {
3351 evaluate(args[1],&argtree,ivar);
3352 newtree->second = argtree;
3353 if (narg > 2) {
3354 newtree->nextra = narg-2;
3355 newtree->extra = new Tree*[narg-2];
3356 for (int i = 2; i < narg; i++) {
3357 evaluate(args[i],&argtree,ivar);
3358 newtree->extra[i-2] = argtree;
3359 }
3360 }
3361 }
3362 treestack[ntreestack++] = newtree;
3363
3364 } else {
3365 value1 = evaluate(args[0],nullptr,ivar);
3366 if (narg > 1) {
3367 value2 = evaluate(args[1],nullptr,ivar);
3368 if (narg > 2) {
3369 for (int i = 2; i < narg; i++)
3370 values[i-2] = evaluate(args[i],nullptr,ivar);
3371 }
3372 }
3373 }
3374
3375 // individual math functions
3376 // customize by adding a function
3377
3378 if (strcmp(word,"sqrt") == 0) {
3379 if (narg != 1)
3380 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3381 if (tree) newtree->type = SQRT;
3382 else {
3383 if (value1 < 0.0)
3384 print_var_error(FLERR,"Sqrt of negative value in "
3385 "variable formula",ivar,0);
3386 argstack[nargstack++] = sqrt(value1);
3387 }
3388
3389 } else if (strcmp(word,"exp") == 0) {
3390 if (narg != 1)
3391 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3392 if (tree) newtree->type = EXP;
3393 else argstack[nargstack++] = exp(value1);
3394 } else if (strcmp(word,"ln") == 0) {
3395 if (narg != 1)
3396 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3397 if (tree) newtree->type = LN;
3398 else {
3399 if (value1 <= 0.0)
3400 print_var_error(FLERR,"Log of zero/negative value in "
3401 "variable formula",ivar,0);
3402 argstack[nargstack++] = log(value1);
3403 }
3404 } else if (strcmp(word,"log") == 0) {
3405 if (narg != 1)
3406 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3407 if (tree) newtree->type = LOG;
3408 else {
3409 if (value1 <= 0.0)
3410 print_var_error(FLERR,"Log of zero/negative value in "
3411 "variable formula",ivar,0);
3412 argstack[nargstack++] = log10(value1);
3413 }
3414 } else if (strcmp(word,"abs") == 0) {
3415 if (narg != 1)
3416 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3417 if (tree) newtree->type = ABS;
3418 else argstack[nargstack++] = fabs(value1);
3419
3420 } else if (strcmp(word,"sin") == 0) {
3421 if (narg != 1)
3422 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3423 if (tree) newtree->type = SIN;
3424 else argstack[nargstack++] = sin(value1);
3425 } else if (strcmp(word,"cos") == 0) {
3426 if (narg != 1)
3427 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3428 if (tree) newtree->type = COS;
3429 else argstack[nargstack++] = cos(value1);
3430 } else if (strcmp(word,"tan") == 0) {
3431 if (narg != 1)
3432 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3433 if (tree) newtree->type = TAN;
3434 else argstack[nargstack++] = tan(value1);
3435
3436 } else if (strcmp(word,"asin") == 0) {
3437 if (narg != 1)
3438 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3439 if (tree) newtree->type = ASIN;
3440 else {
3441 if (value1 < -1.0 || value1 > 1.0)
3442 print_var_error(FLERR,"Arcsin of invalid value in variable formula",ivar,0);
3443 argstack[nargstack++] = asin(value1);
3444 }
3445 } else if (strcmp(word,"acos") == 0) {
3446 if (narg != 1)
3447 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3448 if (tree) newtree->type = ACOS;
3449 else {
3450 if (value1 < -1.0 || value1 > 1.0)
3451 print_var_error(FLERR,"Arccos of invalid value in variable formula",ivar,0);
3452 argstack[nargstack++] = acos(value1);
3453 }
3454 } else if (strcmp(word,"atan") == 0) {
3455 if (narg != 1)
3456 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3457 if (tree) newtree->type = ATAN;
3458 else argstack[nargstack++] = atan(value1);
3459 } else if (strcmp(word,"atan2") == 0) {
3460 if (narg != 2)
3461 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3462 if (tree) newtree->type = ATAN2;
3463 else argstack[nargstack++] = atan2(value1,value2);
3464
3465 } else if (strcmp(word,"random") == 0) {
3466 if (narg != 3)
3467 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3468 if (tree) newtree->type = RANDOM;
3469 else {
3470 if (randomequal == nullptr) {
3471 int seed = static_cast<int> (values[0]);
3472 if (seed <= 0)
3473 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3474 randomequal = new RanMars(lmp,seed);
3475 }
3476 argstack[nargstack++] = randomequal->uniform()*(value2-value1) + value1;
3477 }
3478 } else if (strcmp(word,"normal") == 0) {
3479 if (narg != 3)
3480 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3481 if (tree) newtree->type = NORMAL;
3482 else {
3483 if (value2 < 0.0)
3484 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3485 if (randomequal == nullptr) {
3486 int seed = static_cast<int> (values[0]);
3487 if (seed <= 0)
3488 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3489 randomequal = new RanMars(lmp,seed);
3490 }
3491 argstack[nargstack++] = value1 + value2*randomequal->gaussian();
3492 }
3493
3494 } else if (strcmp(word,"ceil") == 0) {
3495 if (narg != 1)
3496 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3497 if (tree) newtree->type = CEIL;
3498 else argstack[nargstack++] = ceil(value1);
3499
3500 } else if (strcmp(word,"floor") == 0) {
3501 if (narg != 1)
3502 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3503 if (tree) newtree->type = FLOOR;
3504 else argstack[nargstack++] = floor(value1);
3505
3506 } else if (strcmp(word,"round") == 0) {
3507 if (narg != 1)
3508 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3509 if (tree) newtree->type = ROUND;
3510 else argstack[nargstack++] = MYROUND(value1);
3511
3512 } else if (strcmp(word,"ramp") == 0) {
3513 if (narg != 2)
3514 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3515 if (update->whichflag == 0)
3516 print_var_error(FLERR,"Cannot use ramp in "
3517 "variable formula between runs",ivar);
3518 if (tree) newtree->type = RAMP;
3519 else {
3520 double delta = update->ntimestep - update->beginstep;
3521 if (delta != 0.0) delta /= update->endstep - update->beginstep;
3522 double value = value1 + delta*(value2-value1);
3523 argstack[nargstack++] = value;
3524 }
3525
3526 } else if (strcmp(word,"stagger") == 0) {
3527 if (narg != 2)
3528 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3529 if (tree) newtree->type = STAGGER;
3530 else {
3531 bigint ivalue1 = static_cast<bigint> (value1);
3532 bigint ivalue2 = static_cast<bigint> (value2);
3533 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
3534 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3535 bigint lower = update->ntimestep/ivalue1 * ivalue1;
3536 bigint delta = update->ntimestep - lower;
3537 double value;
3538 if (delta < ivalue2) value = lower+ivalue2;
3539 else value = lower+ivalue1;
3540 argstack[nargstack++] = value;
3541 }
3542
3543 } else if (strcmp(word,"logfreq") == 0) {
3544 if (narg != 3)
3545 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3546 if (tree) newtree->type = LOGFREQ;
3547 else {
3548 bigint ivalue1 = static_cast<bigint> (value1);
3549 bigint ivalue2 = static_cast<bigint> (value2);
3550 bigint ivalue3 = static_cast<bigint> (values[0]);
3551 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
3552 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3553 double value;
3554 if (update->ntimestep < ivalue1) value = ivalue1;
3555 else {
3556 bigint lower = ivalue1;
3557 while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
3558 bigint multiple = update->ntimestep/lower;
3559 if (multiple < ivalue2) value = (multiple+1)*lower;
3560 else value = lower*ivalue3;
3561 }
3562 argstack[nargstack++] = value;
3563 }
3564
3565 } else if (strcmp(word,"logfreq2") == 0) {
3566 if (narg != 3)
3567 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3568 if (tree) newtree->type = LOGFREQ2;
3569 else {
3570 bigint ivalue1 = static_cast<bigint> (value1);
3571 bigint ivalue2 = static_cast<bigint> (value2);
3572 bigint ivalue3 = static_cast<bigint> (values[0]);
3573 if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 )
3574 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3575 double value;
3576 if (update->ntimestep < ivalue1) value = ivalue1;
3577 else {
3578 value = ivalue1;
3579 double delta = ivalue1*(ivalue3-1.0)/ivalue2;
3580 bigint count = 0;
3581 while (update->ntimestep >= value) {
3582 value += delta;
3583 count++;
3584 if (count % ivalue2 == 0) delta *= ivalue3;
3585 }
3586 }
3587 argstack[nargstack++] = ceil(value);
3588 }
3589
3590 } else if (strcmp(word,"logfreq3") == 0) {
3591 if (narg != 3)
3592 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3593 if (tree) newtree->type = LOGFREQ3;
3594 else {
3595 bigint ivalue1 = static_cast<bigint> (value1);
3596 bigint ivalue2 = static_cast<bigint> (value2);
3597 bigint ivalue3 = static_cast<bigint> (values[0]);
3598 if (ivalue1 <= 0 || ivalue2 <= 1 || ivalue3 <= 0 ||
3599 ivalue3-ivalue1+1 < ivalue2 )
3600 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3601 double value;
3602 if (update->ntimestep < ivalue1) value = ivalue1;
3603 //else if (update->ntimestep <= ivalue3) {
3604 else {
3605 value = ivalue1;
3606 double logsp = ivalue1;
3607 double factor = pow(((double)ivalue3)/ivalue1, 1.0/(ivalue2-1));
3608 bigint linsp = ivalue1;
3609 while (update->ntimestep >= value) {
3610 logsp *= factor;
3611 linsp++;
3612 if (linsp > logsp) value = linsp;
3613 else value = ceil(logsp)-(((bigint)ceil(logsp)-1)/ivalue3);
3614 }
3615 }
3616 if (update->ntimestep > ivalue3)
3617 error->all(FLERR,"Calls to variable exceeded limit");
3618 argstack[nargstack++] = value;
3619 }
3620
3621 } else if (strcmp(word,"stride") == 0) {
3622 if (narg != 3)
3623 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3624 if (tree) newtree->type = STRIDE;
3625 else {
3626 bigint ivalue1 = static_cast<bigint> (value1);
3627 bigint ivalue2 = static_cast<bigint> (value2);
3628 bigint ivalue3 = static_cast<bigint> (values[0]);
3629 if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
3630 error->one(FLERR,"Invalid math function in variable formula");
3631 double value;
3632 if (update->ntimestep < ivalue1) value = ivalue1;
3633 else if (update->ntimestep < ivalue2) {
3634 bigint offset = update->ntimestep - ivalue1;
3635 value = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
3636 if (value > ivalue2) value = (double) MAXBIGINT;
3637 } else value = (double) MAXBIGINT;
3638 argstack[nargstack++] = value;
3639 }
3640
3641 } else if (strcmp(word,"stride2") == 0) {
3642 if (narg != 6)
3643 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3644 if (tree) newtree->type = STRIDE2;
3645 else {
3646 bigint ivalue1 = static_cast<bigint> (value1);
3647 bigint ivalue2 = static_cast<bigint> (value2);
3648 bigint ivalue3 = static_cast<bigint> (values[0]);
3649 bigint ivalue4 = static_cast<bigint> (values[1]);
3650 bigint ivalue5 = static_cast<bigint> (values[2]);
3651 bigint ivalue6 = static_cast<bigint> (values[3]);
3652 if (ivalue1 < 0 || ivalue2 < 0 || ivalue3 <= 0 || ivalue1 > ivalue2)
3653 error->one(FLERR,"Invalid math function in variable formula");
3654 if (ivalue4 < 0 || ivalue5 < 0 || ivalue6 <= 0 || ivalue4 > ivalue5)
3655 error->one(FLERR,"Invalid math function in variable formula");
3656 if (ivalue4 < ivalue1 || ivalue5 > ivalue2)
3657 error->one(FLERR,"Invalid math function in variable formula");
3658 bigint istep, offset;
3659 if (update->ntimestep < ivalue1) istep = ivalue1;
3660 else if (update->ntimestep < ivalue2) {
3661 if (update->ntimestep < ivalue4 || update->ntimestep > ivalue5) {
3662 offset = update->ntimestep - ivalue1;
3663 istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
3664 if (update->ntimestep < ivalue4 && istep > ivalue4) istep = ivalue4;
3665 } else {
3666 offset = update->ntimestep - ivalue4;
3667 istep = ivalue4 + (offset/ivalue6)*ivalue6 + ivalue6;
3668 if (istep > ivalue5) {
3669 offset = ivalue5 - ivalue1;
3670 istep = ivalue1 + (offset/ivalue3)*ivalue3 + ivalue3;
3671 if (istep > ivalue2) istep = MAXBIGINT;
3672 }
3673 }
3674 } else istep = MAXBIGINT;
3675 double value = istep;
3676 argstack[nargstack++] = value;
3677 }
3678
3679 } else if (strcmp(word,"vdisplace") == 0) {
3680 if (narg != 2)
3681 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3682 if (update->whichflag == 0)
3683 print_var_error(FLERR,"Cannot use vdisplace in "
3684 "variable formula between runs",ivar);
3685 if (tree) newtree->type = VDISPLACE;
3686 else {
3687 double delta = update->ntimestep - update->beginstep;
3688 double value = value1 + value2*delta*update->dt;
3689 argstack[nargstack++] = value;
3690 }
3691
3692 } else if (strcmp(word,"swiggle") == 0) {
3693 if (narg != 3)
3694 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3695 if (update->whichflag == 0)
3696 print_var_error(FLERR,"Cannot use swiggle in "
3697 "variable formula between runs",ivar);
3698 if (tree) newtree->type = CWIGGLE;
3699 else {
3700 if (values[0] == 0.0)
3701 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3702 double delta = update->ntimestep - update->beginstep;
3703 double omega = 2.0*MY_PI/values[0];
3704 double value = value1 + value2*sin(omega*delta*update->dt);
3705 argstack[nargstack++] = value;
3706 }
3707
3708 } else if (strcmp(word,"cwiggle") == 0) {
3709 if (narg != 3)
3710 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3711 if (update->whichflag == 0)
3712 print_var_error(FLERR,"Cannot use cwiggle in "
3713 "variable formula between runs",ivar);
3714 if (tree) newtree->type = CWIGGLE;
3715 else {
3716 if (values[0] == 0.0)
3717 print_var_error(FLERR,"Invalid math function in variable formula",ivar);
3718 double delta = update->ntimestep - update->beginstep;
3719 double omega = 2.0*MY_PI/values[0];
3720 double value = value1 + value2*(1.0-cos(omega*delta*update->dt));
3721 argstack[nargstack++] = value;
3722 }
3723 }
3724
3725 // delete stored args
3726
3727 for (int i = 0; i < narg; i++) delete [] args[i];
3728
3729 return 1;
3730 }
3731
3732 /* ----------------------------------------------------------------------
3733 process a group function in formula with optional region arg
3734 push result onto tree or arg stack
3735 word = group function
3736 contents = str between parentheses with one,two,three args
3737 return 0 if not a match, 1 if successfully processed
3738 customize by adding a group function with optional region arg:
3739 count(group),mass(group),charge(group),
3740 xcm(group,dim),vcm(group,dim),fcm(group,dim),
3741 bound(group,xmin),gyration(group),ke(group),angmom(group,dim),
3742 torque(group,dim),inertia(group,dim),omega(group,dim)
3743 ------------------------------------------------------------------------- */
3744
group_function(char * word,char * contents,Tree ** tree,Tree ** treestack,int & ntreestack,double * argstack,int & nargstack,int ivar)3745 int Variable::group_function(char *word, char *contents, Tree **tree,
3746 Tree **treestack, int &ntreestack,
3747 double *argstack, int &nargstack, int ivar)
3748 {
3749 // word not a match to any group function
3750
3751 if (strcmp(word,"count") && strcmp(word,"mass") &&
3752 strcmp(word,"charge") && strcmp(word,"xcm") &&
3753 strcmp(word,"vcm") && strcmp(word,"fcm") &&
3754 strcmp(word,"bound") && strcmp(word,"gyration") &&
3755 strcmp(word,"ke") && strcmp(word,"angmom") &&
3756 strcmp(word,"torque") && strcmp(word,"inertia") &&
3757 strcmp(word,"omega"))
3758 return 0;
3759
3760 // parse contents for comma-separated args
3761 // narg = number of args, args = strings between commas
3762
3763 char *args[MAXFUNCARG];
3764 int narg = parse_args(contents,args);
3765
3766 // group to operate on
3767
3768 int igroup = group->find(args[0]);
3769 if (igroup == -1) {
3770 std::string mesg = "Group ID '";
3771 mesg += args[0];
3772 mesg += "' in variable formula does not exist";
3773 print_var_error(FLERR,mesg,ivar);
3774 }
3775
3776 // match word to group function
3777
3778 double value = 0.0;
3779
3780 if (strcmp(word,"count") == 0) {
3781 if (narg == 1) value = group->count(igroup);
3782 else if (narg == 2)
3783 value = group->count(igroup,region_function(args[1],ivar));
3784 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3785
3786 } else if (strcmp(word,"mass") == 0) {
3787 if (narg == 1) value = group->mass(igroup);
3788 else if (narg == 2) value = group->mass(igroup,region_function(args[1],ivar));
3789 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3790
3791 } else if (strcmp(word,"charge") == 0) {
3792 if (narg == 1) value = group->charge(igroup);
3793 else if (narg == 2)
3794 value = group->charge(igroup,region_function(args[1],ivar));
3795 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3796
3797 } else if (strcmp(word,"xcm") == 0) {
3798 atom->check_mass(FLERR);
3799 double xcm[3];
3800 if (narg == 2) {
3801 double masstotal = group->mass(igroup);
3802 group->xcm(igroup,masstotal,xcm);
3803 } else if (narg == 3) {
3804 int iregion = region_function(args[2],ivar);
3805 double masstotal = group->mass(igroup,iregion);
3806 group->xcm(igroup,masstotal,xcm,iregion);
3807 } else print_var_error(FLERR,"Invalid group function in "
3808 "variable formula",ivar);
3809 if (strcmp(args[1],"x") == 0) value = xcm[0];
3810 else if (strcmp(args[1],"y") == 0) value = xcm[1];
3811 else if (strcmp(args[1],"z") == 0) value = xcm[2];
3812 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3813
3814 } else if (strcmp(word,"vcm") == 0) {
3815 atom->check_mass(FLERR);
3816 double vcm[3];
3817 if (narg == 2) {
3818 double masstotal = group->mass(igroup);
3819 group->vcm(igroup,masstotal,vcm);
3820 } else if (narg == 3) {
3821 int iregion = region_function(args[2],ivar);
3822 double masstotal = group->mass(igroup,iregion);
3823 group->vcm(igroup,masstotal,vcm,iregion);
3824 } else print_var_error(FLERR,"Invalid group function in "
3825 "variable formula",ivar);
3826 if (strcmp(args[1],"x") == 0) value = vcm[0];
3827 else if (strcmp(args[1],"y") == 0) value = vcm[1];
3828 else if (strcmp(args[1],"z") == 0) value = vcm[2];
3829 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3830
3831 } else if (strcmp(word,"fcm") == 0) {
3832 double fcm[3];
3833 if (narg == 2) group->fcm(igroup,fcm);
3834 else if (narg == 3) group->fcm(igroup,fcm,region_function(args[2],ivar));
3835 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3836 if (strcmp(args[1],"x") == 0) value = fcm[0];
3837 else if (strcmp(args[1],"y") == 0) value = fcm[1];
3838 else if (strcmp(args[1],"z") == 0) value = fcm[2];
3839 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3840
3841 } else if (strcmp(word,"bound") == 0) {
3842 double minmax[6];
3843 if (narg == 2) group->bounds(igroup,minmax);
3844 else if (narg == 3)
3845 group->bounds(igroup,minmax,region_function(args[2],ivar));
3846 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3847 if (strcmp(args[1],"xmin") == 0) value = minmax[0];
3848 else if (strcmp(args[1],"xmax") == 0) value = minmax[1];
3849 else if (strcmp(args[1],"ymin") == 0) value = minmax[2];
3850 else if (strcmp(args[1],"ymax") == 0) value = minmax[3];
3851 else if (strcmp(args[1],"zmin") == 0) value = minmax[4];
3852 else if (strcmp(args[1],"zmax") == 0) value = minmax[5];
3853 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3854
3855 } else if (strcmp(word,"gyration") == 0) {
3856 atom->check_mass(FLERR);
3857 double xcm[3];
3858 if (narg == 1) {
3859 double masstotal = group->mass(igroup);
3860 group->xcm(igroup,masstotal,xcm);
3861 value = group->gyration(igroup,masstotal,xcm);
3862 } else if (narg == 2) {
3863 int iregion = region_function(args[1],ivar);
3864 double masstotal = group->mass(igroup,iregion);
3865 group->xcm(igroup,masstotal,xcm,iregion);
3866 value = group->gyration(igroup,masstotal,xcm,iregion);
3867 } else print_var_error(FLERR,"Invalid group function in "
3868 "variable formula",ivar);
3869
3870 } else if (strcmp(word,"ke") == 0) {
3871 if (narg == 1) value = group->ke(igroup);
3872 else if (narg == 2) value = group->ke(igroup,region_function(args[1],ivar));
3873 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3874
3875 } else if (strcmp(word,"angmom") == 0) {
3876 atom->check_mass(FLERR);
3877 double xcm[3],lmom[3];
3878 if (narg == 2) {
3879 double masstotal = group->mass(igroup);
3880 group->xcm(igroup,masstotal,xcm);
3881 group->angmom(igroup,xcm,lmom);
3882 } else if (narg == 3) {
3883 int iregion = region_function(args[2],ivar);
3884 double masstotal = group->mass(igroup,iregion);
3885 group->xcm(igroup,masstotal,xcm,iregion);
3886 group->angmom(igroup,xcm,lmom,iregion);
3887 } else print_var_error(FLERR,"Invalid group function in "
3888 "variable formula",ivar);
3889 if (strcmp(args[1],"x") == 0) value = lmom[0];
3890 else if (strcmp(args[1],"y") == 0) value = lmom[1];
3891 else if (strcmp(args[1],"z") == 0) value = lmom[2];
3892 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3893
3894 } else if (strcmp(word,"torque") == 0) {
3895 atom->check_mass(FLERR);
3896 double xcm[3],tq[3];
3897 if (narg == 2) {
3898 double masstotal = group->mass(igroup);
3899 group->xcm(igroup,masstotal,xcm);
3900 group->torque(igroup,xcm,tq);
3901 } else if (narg == 3) {
3902 int iregion = region_function(args[2],ivar);
3903 double masstotal = group->mass(igroup,iregion);
3904 group->xcm(igroup,masstotal,xcm,iregion);
3905 group->torque(igroup,xcm,tq,iregion);
3906 } else print_var_error(FLERR,"Invalid group function in "
3907 "variable formula",ivar);
3908 if (strcmp(args[1],"x") == 0) value = tq[0];
3909 else if (strcmp(args[1],"y") == 0) value = tq[1];
3910 else if (strcmp(args[1],"z") == 0) value = tq[2];
3911 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3912
3913 } else if (strcmp(word,"inertia") == 0) {
3914 atom->check_mass(FLERR);
3915 double xcm[3],inertia[3][3];
3916 if (narg == 2) {
3917 double masstotal = group->mass(igroup);
3918 group->xcm(igroup,masstotal,xcm);
3919 group->inertia(igroup,xcm,inertia);
3920 } else if (narg == 3) {
3921 int iregion = region_function(args[2],ivar);
3922 double masstotal = group->mass(igroup,iregion);
3923 group->xcm(igroup,masstotal,xcm,iregion);
3924 group->inertia(igroup,xcm,inertia,iregion);
3925 } else print_var_error(FLERR,"Invalid group function in "
3926 "variable formula",ivar);
3927 if (strcmp(args[1],"xx") == 0) value = inertia[0][0];
3928 else if (strcmp(args[1],"yy") == 0) value = inertia[1][1];
3929 else if (strcmp(args[1],"zz") == 0) value = inertia[2][2];
3930 else if (strcmp(args[1],"xy") == 0) value = inertia[0][1];
3931 else if (strcmp(args[1],"yz") == 0) value = inertia[1][2];
3932 else if (strcmp(args[1],"xz") == 0) value = inertia[0][2];
3933 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3934
3935 } else if (strcmp(word,"omega") == 0) {
3936 atom->check_mass(FLERR);
3937 double xcm[3],angmom[3],inertia[3][3],omega[3];
3938 if (narg == 2) {
3939 double masstotal = group->mass(igroup);
3940 group->xcm(igroup,masstotal,xcm);
3941 group->angmom(igroup,xcm,angmom);
3942 group->inertia(igroup,xcm,inertia);
3943 group->omega(angmom,inertia,omega);
3944 } else if (narg == 3) {
3945 int iregion = region_function(args[2],ivar);
3946 double masstotal = group->mass(igroup,iregion);
3947 group->xcm(igroup,masstotal,xcm,iregion);
3948 group->angmom(igroup,xcm,angmom,iregion);
3949 group->inertia(igroup,xcm,inertia,iregion);
3950 group->omega(angmom,inertia,omega);
3951 } else print_var_error(FLERR,"Invalid group function in "
3952 "variable formula",ivar);
3953 if (strcmp(args[1],"x") == 0) value = omega[0];
3954 else if (strcmp(args[1],"y") == 0) value = omega[1];
3955 else if (strcmp(args[1],"z") == 0) value = omega[2];
3956 else print_var_error(FLERR,"Invalid group function in variable formula",ivar);
3957 }
3958
3959 // delete stored args
3960
3961 for (int i = 0; i < narg; i++) delete [] args[i];
3962
3963 // save value in tree or on argstack
3964
3965 if (tree) {
3966 Tree *newtree = new Tree();
3967 newtree->type = VALUE;
3968 newtree->value = value;
3969 treestack[ntreestack++] = newtree;
3970 } else argstack[nargstack++] = value;
3971
3972 return 1;
3973 }
3974
3975 /* ---------------------------------------------------------------------- */
3976
region_function(char * id,int ivar)3977 int Variable::region_function(char *id, int ivar)
3978 {
3979 int iregion = domain->find_region(id);
3980 if (iregion == -1) {
3981 std::string mesg = "Region ID '";
3982 mesg += id;
3983 mesg += "' in variable formula does not exist";
3984 print_var_error(FLERR,mesg,ivar);
3985 }
3986
3987 // init region in case sub-regions have been deleted
3988
3989 domain->regions[iregion]->init();
3990
3991 return iregion;
3992 }
3993
3994 /* ----------------------------------------------------------------------
3995 process a special function in formula
3996 push result onto tree or arg stack
3997 word = special function
3998 contents = str between parentheses with one,two,three args
3999 return 0 if not a match, 1 if successfully processed
4000 customize by adding a special function:
4001 sum(x),min(x),max(x),ave(x),trap(x),slope(x),
4002 gmask(x),rmask(x),grmask(x,y),next(x)
4003 ------------------------------------------------------------------------- */
4004
special_function(char * word,char * contents,Tree ** tree,Tree ** treestack,int & ntreestack,double * argstack,int & nargstack,int ivar)4005 int Variable::special_function(char *word, char *contents, Tree **tree,
4006 Tree **treestack, int &ntreestack,
4007 double *argstack, int &nargstack, int ivar)
4008 {
4009 double sx,sxx;
4010 double value,sy,sxy;
4011
4012 // word not a match to any special function
4013
4014 if (strcmp(word,"sum") && strcmp(word,"min") && strcmp(word,"max") &&
4015 strcmp(word,"ave") && strcmp(word,"trap") && strcmp(word,"slope") &&
4016 strcmp(word,"gmask") && strcmp(word,"rmask") &&
4017 strcmp(word,"grmask") && strcmp(word,"next") &&
4018 strcmp(word,"is_active") && strcmp(word,"is_defined") &&
4019 strcmp(word,"is_available") && strcmp(word,"is_file"))
4020 return 0;
4021
4022 // parse contents for comma-separated args
4023 // narg = number of args, args = strings between commas
4024
4025 char *args[MAXFUNCARG];
4026 int narg = parse_args(contents,args);
4027
4028 // special functions that operate on global vectors
4029
4030 if (strcmp(word,"sum") == 0 || strcmp(word,"min") == 0 ||
4031 strcmp(word,"max") == 0 || strcmp(word,"ave") == 0 ||
4032 strcmp(word,"trap") == 0 || strcmp(word,"slope") == 0) {
4033
4034 int method = 0;
4035 if (strcmp(word,"sum") == 0) method = SUM;
4036 else if (strcmp(word,"min") == 0) method = XMIN;
4037 else if (strcmp(word,"max") == 0) method = XMAX;
4038 else if (strcmp(word,"ave") == 0) method = AVE;
4039 else if (strcmp(word,"trap") == 0) method = TRAP;
4040 else if (strcmp(word,"slope") == 0) method = SLOPE;
4041
4042 if (narg != 1)
4043 print_var_error(FLERR,"Invalid special function in variable formula",ivar);
4044
4045 Compute *compute = nullptr;
4046 Fix *fix = nullptr;
4047 int index,nvec,nstride;
4048 char *ptr1,*ptr2;
4049 int ivar = -1;
4050
4051 // argument is compute
4052
4053 if (utils::strmatch(args[0],"^c_")) {
4054 ptr1 = strchr(args[0],'[');
4055 if (ptr1) {
4056 ptr2 = ptr1;
4057 index = (int) int_between_brackets(ptr2,0);
4058 *ptr1 = '\0';
4059 } else index = 0;
4060
4061 int icompute = modify->find_compute(&args[0][2]);
4062 if (icompute < 0) {
4063 std::string mesg = "Invalid compute ID '";
4064 mesg += (args[0]+2);
4065 mesg += "' in variable formula";
4066 print_var_error(FLERR,mesg,ivar);
4067 }
4068 compute = modify->compute[icompute];
4069 if (index == 0 && compute->vector_flag) {
4070 if (update->whichflag == 0) {
4071 if (compute->invoked_vector != update->ntimestep)
4072 print_var_error(FLERR,"Compute used in variable between runs "
4073 "is not current",ivar);
4074 } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
4075 compute->compute_vector();
4076 compute->invoked_flag |= Compute::INVOKED_VECTOR;
4077 }
4078 nvec = compute->size_vector;
4079 nstride = 1;
4080 } else if (index && compute->array_flag) {
4081 if (index > compute->size_array_cols)
4082 print_var_error(FLERR,"Variable formula compute array "
4083 "is accessed out-of-range",ivar,0);
4084 if (update->whichflag == 0) {
4085 if (compute->invoked_array != update->ntimestep)
4086 print_var_error(FLERR,"Compute used in variable between runs "
4087 "is not current",ivar);
4088 } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) {
4089 compute->compute_array();
4090 compute->invoked_flag |= Compute::INVOKED_ARRAY;
4091 }
4092 nvec = compute->size_array_rows;
4093 nstride = compute->size_array_cols;
4094 } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar);
4095
4096 // argument is fix
4097
4098 } else if (utils::strmatch(args[0],"^f_")) {
4099 ptr1 = strchr(args[0],'[');
4100 if (ptr1) {
4101 ptr2 = ptr1;
4102 index = (int) int_between_brackets(ptr2,0);
4103 *ptr1 = '\0';
4104 } else index = 0;
4105
4106 int ifix = modify->find_fix(&args[0][2]);
4107 if (ifix < 0) {
4108 std::string mesg = "Invalid fix ID '";
4109 mesg += (args[0]+2);
4110 mesg += "' in variable formula";
4111 print_var_error(FLERR,mesg,ivar);
4112 }
4113 fix = modify->fix[ifix];
4114 if (index == 0 && fix->vector_flag) {
4115 if (update->whichflag > 0 && update->ntimestep % fix->global_freq) {
4116 std::string mesg = "Fix with ID '";
4117 mesg += (args[0]+2);
4118 mesg += "' in variable formula not computed at compatible time";
4119 print_var_error(FLERR,mesg,ivar);
4120 }
4121 nvec = fix->size_vector;
4122 nstride = 1;
4123 } else if (index && fix->array_flag) {
4124 if (index > fix->size_array_cols)
4125 print_var_error(FLERR,"Variable formula fix array "
4126 "is accessed out-of-range",ivar);
4127 if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
4128 print_var_error(FLERR,"Fix in variable not computed at "
4129 "compatible time",ivar);
4130 nvec = fix->size_array_rows;
4131 nstride = fix->size_array_cols;
4132 } else print_var_error(FLERR,"Mismatched fix in variable formula",ivar);
4133
4134 // argument is vector-style variable
4135
4136 } else if (utils::strmatch(args[0],"^v_")) {
4137 ptr1 = strchr(args[0],'[');
4138 if (ptr1) {
4139 ptr2 = ptr1;
4140 index = (int) int_between_brackets(ptr2,0);
4141 *ptr1 = '\0';
4142 } else index = 0;
4143
4144 if (index)
4145 print_var_error(FLERR,"Invalid special function in "
4146 "variable formula",ivar);
4147 ivar = find(&args[0][2]);
4148 if (ivar < 0)
4149 print_var_error(FLERR,"Invalid special function in "
4150 "variable formula",ivar);
4151 if (style[ivar] != VECTOR)
4152 print_var_error(FLERR,"Mis-matched special function variable "
4153 "in variable formula",ivar);
4154 if (eval_in_progress[ivar])
4155 print_var_error(FLERR,"has a circular dependency",ivar);
4156
4157 double *vec;
4158 nvec = compute_vector(ivar,&vec);
4159 nstride = 1;
4160
4161 if ((method == AVE) && (nvec == 0))
4162 print_var_error(FLERR,"Cannot compute average of empty vector",ivar);
4163
4164
4165 } else print_var_error(FLERR,"Invalid special function in "
4166 "variable formula",ivar);
4167
4168 value = 0.0;
4169 if (method == SLOPE) sx = sxx = sy = sxy = 0.0;
4170 else if (method == XMIN) value = BIG;
4171 else if (method == XMAX) value = -BIG;
4172
4173 if (compute) {
4174 double *vec;
4175 if (index) {
4176 if (compute->array) vec = &compute->array[0][index-1];
4177 else vec = nullptr;
4178 } else vec = compute->vector;
4179
4180 int j = 0;
4181 for (int i = 0; i < nvec; i++) {
4182 if (method == SUM) value += vec[j];
4183 else if (method == XMIN) value = MIN(value,vec[j]);
4184 else if (method == XMAX) value = MAX(value,vec[j]);
4185 else if (method == AVE) value += vec[j];
4186 else if (method == TRAP) value += vec[j];
4187 else if (method == SLOPE) {
4188 sx += (double)i;
4189 sy += vec[j];
4190 sxx += (double)i * (double)i;
4191 sxy += (double)i * vec[j];
4192 }
4193 j += nstride;
4194 }
4195 if (method == TRAP) value -= 0.5*vec[0] + 0.5*vec[nvec-1];
4196 }
4197
4198 if (fix) {
4199 double one;
4200 for (int i = 0; i < nvec; i++) {
4201 if (index) one = fix->compute_array(i,index-1);
4202 else one = fix->compute_vector(i);
4203 if (method == SUM) value += one;
4204 else if (method == XMIN) value = MIN(value,one);
4205 else if (method == XMAX) value = MAX(value,one);
4206 else if (method == AVE) value += one;
4207 else if (method == TRAP) value += one;
4208 else if (method == SLOPE) {
4209 sx += (double)i;
4210 sy += one;
4211 sxx += (double)i * (double)i;
4212 sxy += (double)i * one;
4213 }
4214 }
4215 if (method == TRAP) {
4216 if (index) value -= 0.5*fix->compute_array(0,index-1) +
4217 0.5*fix->compute_array(nvec-1,index-1);
4218 else value -= 0.5*fix->compute_vector(0) +
4219 0.5*fix->compute_vector(nvec-1);
4220 }
4221 }
4222
4223 if (ivar >= 0) {
4224 double one;
4225 double *vec = vecs[ivar].values;
4226 for (int i = 0; i < nvec; i++) {
4227 one = vec[i];
4228 if (method == SUM) value += one;
4229 else if (method == XMIN) value = MIN(value,one);
4230 else if (method == XMAX) value = MAX(value,one);
4231 else if (method == AVE) value += one;
4232 else if (method == TRAP) value += one;
4233 else if (method == SLOPE) {
4234 sx += (double) i;
4235 sy += one;
4236 sxx += (double)i * (double)i;
4237 sxy += (double)i * one;
4238 }
4239 }
4240 if (method == TRAP) value -= 0.5*vec[0] + 0.5*vec[nvec-1];
4241 }
4242
4243 if (method == AVE) value /= nvec;
4244
4245 if (method == SLOPE) {
4246 double numerator = nvec*sxy - sx*sy;
4247 double denominator = nvec*sxx - sx*sx;
4248 if (denominator != 0.0) value = numerator/denominator;
4249 else value = BIG;
4250 }
4251
4252 // save value in tree or on argstack
4253
4254 if (tree) {
4255 Tree *newtree = new Tree();
4256 newtree->type = VALUE;
4257 newtree->value = value;
4258 treestack[ntreestack++] = newtree;
4259 } else argstack[nargstack++] = value;
4260
4261 // mask special functions
4262
4263 } else if (strcmp(word,"gmask") == 0) {
4264 if (tree == nullptr)
4265 print_var_error(FLERR,"Gmask function in equal-style "
4266 "variable formula",ivar);
4267 if (narg != 1)
4268 print_var_error(FLERR,"Invalid special function in variable formula",ivar);
4269
4270 int igroup = group->find(args[0]);
4271 if (igroup == -1)
4272 print_var_error(FLERR,"Group ID in variable formula does not exist",ivar);
4273
4274 Tree *newtree = new Tree();
4275 newtree->type = GMASK;
4276 newtree->ivalue1 = group->bitmask[igroup];
4277 treestack[ntreestack++] = newtree;
4278
4279 } else if (strcmp(word,"rmask") == 0) {
4280 if (tree == nullptr)
4281 print_var_error(FLERR,"Rmask function in equal-style "
4282 "variable formula",ivar);
4283 if (narg != 1)
4284 print_var_error(FLERR,"Invalid special function in variable formula",ivar);
4285
4286 int iregion = region_function(args[0],ivar);
4287 domain->regions[iregion]->prematch();
4288
4289 Tree *newtree = new Tree();
4290 newtree->type = RMASK;
4291 newtree->ivalue1 = iregion;
4292 treestack[ntreestack++] = newtree;
4293
4294 } else if (strcmp(word,"grmask") == 0) {
4295 if (tree == nullptr)
4296 print_var_error(FLERR,"Grmask function in equal-style "
4297 "variable formula",ivar);
4298 if (narg != 2)
4299 print_var_error(FLERR,"Invalid special function in variable formula",ivar);
4300
4301 int igroup = group->find(args[0]);
4302 if (igroup == -1)
4303 print_var_error(FLERR,"Group ID in variable formula does not exist",ivar);
4304 int iregion = region_function(args[1],ivar);
4305 domain->regions[iregion]->prematch();
4306
4307 Tree *newtree = new Tree();
4308 newtree->type = GRMASK;
4309 newtree->ivalue1 = group->bitmask[igroup];
4310 newtree->ivalue2 = iregion;
4311 treestack[ntreestack++] = newtree;
4312
4313 // special function for file-style or atomfile-style variables
4314
4315 } else if (strcmp(word,"next") == 0) {
4316 if (narg != 1)
4317 print_var_error(FLERR,"Invalid special function in variable formula",ivar);
4318
4319 int ivar = find(args[0]);
4320 if (ivar < 0) {
4321 std::string mesg = "Variable ID '";
4322 mesg += args[0];
4323 mesg += "' in variable formula does not exist";
4324 print_var_error(FLERR,mesg,ivar);
4325 }
4326
4327 // SCALARFILE has single current value, read next one
4328 // save value in tree or on argstack
4329
4330 if (style[ivar] == SCALARFILE) {
4331 double value = atof(data[ivar][0]);
4332 int done = reader[ivar]->read_scalar(data[ivar][0]);
4333 if (done) remove(ivar);
4334
4335 if (tree) {
4336 Tree *newtree = new Tree();
4337 newtree->type = VALUE;
4338 newtree->value = value;
4339 treestack[ntreestack++] = newtree;
4340 } else argstack[nargstack++] = value;
4341
4342 // ATOMFILE has per-atom values, save values in tree
4343 // copy current per-atom values into result so can read next ones
4344 // set selfalloc = 1 so result will be deleted by free_tree() after eval
4345
4346 } else if (style[ivar] == ATOMFILE) {
4347 if (tree == nullptr)
4348 print_var_error(FLERR,"Atomfile variable in "
4349 "equal-style variable formula",ivar);
4350
4351 double *result;
4352 memory->create(result,atom->nlocal,"variable:result");
4353 memcpy(result,reader[ivar]->fixstore->vstore,atom->nlocal*sizeof(double));
4354
4355 int done = reader[ivar]->read_peratom();
4356 if (done) remove(ivar);
4357
4358 Tree *newtree = new Tree();
4359 newtree->type = ATOMARRAY;
4360 newtree->array = result;
4361 newtree->nstride = 1;
4362 newtree->selfalloc = 1;
4363 treestack[ntreestack++] = newtree;
4364
4365 } else print_var_error(FLERR,"Invalid variable style in "
4366 "special function next",ivar);
4367
4368 } else if (strcmp(word,"is_active") == 0) {
4369 if (narg != 2)
4370 print_var_error(FLERR,"Invalid is_active() function in "
4371 "variable formula",ivar);
4372
4373 Info info(lmp);
4374 value = (info.is_active(args[0],args[1])) ? 1.0 : 0.0;
4375
4376 // save value in tree or on argstack
4377
4378 if (tree) {
4379 Tree *newtree = new Tree();
4380 newtree->type = VALUE;
4381 newtree->value = value;
4382 treestack[ntreestack++] = newtree;
4383 } else argstack[nargstack++] = value;
4384
4385 } else if (strcmp(word,"is_available") == 0) {
4386 if (narg != 2)
4387 print_var_error(FLERR,"Invalid is_available() function in "
4388 "variable formula",ivar);
4389
4390 Info info(lmp);
4391 value = (info.is_available(args[0],args[1])) ? 1.0 : 0.0;
4392
4393 // save value in tree or on argstack
4394
4395 if (tree) {
4396 Tree *newtree = new Tree();
4397 newtree->type = VALUE;
4398 newtree->value = value;
4399 treestack[ntreestack++] = newtree;
4400 } else argstack[nargstack++] = value;
4401
4402 } else if (strcmp(word,"is_defined") == 0) {
4403 if (narg != 2)
4404 print_var_error(FLERR,"Invalid is_defined() function in "
4405 "variable formula",ivar);
4406
4407 Info info(lmp);
4408 value = (info.is_defined(args[0],args[1])) ? 1.0 : 0.0;
4409
4410 // save value in tree or on argstack
4411
4412 if (tree) {
4413 Tree *newtree = new Tree();
4414 newtree->type = VALUE;
4415 newtree->value = value;
4416 treestack[ntreestack++] = newtree;
4417 } else argstack[nargstack++] = value;
4418
4419 } else if (strcmp(word,"is_file") == 0) {
4420 if (narg != 1)
4421 print_var_error(FLERR,"Invalid is_file() function in "
4422 "variable formula",ivar);
4423
4424 FILE *fp = fopen(args[0],"r");
4425 value = (fp == nullptr) ? 0.0 : 1.0;
4426 if (fp) fclose(fp);
4427
4428 // save value in tree or on argstack
4429
4430 if (tree) {
4431 Tree *newtree = new Tree();
4432 newtree->type = VALUE;
4433 newtree->value = value;
4434 treestack[ntreestack++] = newtree;
4435 } else argstack[nargstack++] = value;
4436 }
4437
4438 // delete stored args
4439
4440 for (int i = 0; i < narg; i++) delete [] args[i];
4441
4442 return 1;
4443 }
4444
4445 /* ----------------------------------------------------------------------
4446 extract a global value from a per-atom quantity in a formula
4447 flag = 0 -> word is an atom vector
4448 flag = 1 -> vector is a per-atom compute or fix quantity with nstride
4449 id = global ID of atom, converted to local index
4450 push result onto tree or arg stack
4451 customize by adding an atom vector:
4452 id,mass,type,mol,x,y,z,vx,vy,vz,fx,fy,fz,q
4453 ------------------------------------------------------------------------- */
4454
peratom2global(int flag,char * word,double * vector,int nstride,tagint id,Tree ** tree,Tree ** treestack,int & ntreestack,double * argstack,int & nargstack)4455 void Variable::peratom2global(int flag, char *word,
4456 double *vector, int nstride, tagint id,
4457 Tree **tree, Tree **treestack, int &ntreestack,
4458 double *argstack, int &nargstack)
4459 {
4460 // error check for ID larger than any atom
4461 // int_between_brackets() already checked for ID <= 0
4462
4463 if (atom->map_style == Atom::MAP_NONE)
4464 error->all(FLERR,
4465 "Indexed per-atom vector in variable formula without atom map");
4466
4467 if (id > atom->map_tag_max)
4468 error->all(FLERR,"Variable atom ID is too large");
4469
4470 // if ID does not exist, index will be -1 for all procs,
4471 // and mine will be set to 0.0
4472
4473 int index = atom->map(id);
4474
4475 double mine;
4476 if (index >= 0 && index < atom->nlocal) {
4477
4478 if (flag == 0) {
4479 if (strcmp(word,"id") == 0) mine = atom->tag[index];
4480 else if (strcmp(word,"mass") == 0) {
4481 if (atom->rmass) mine = atom->rmass[index];
4482 else mine = atom->mass[atom->type[index]];
4483 }
4484 else if (strcmp(word,"type") == 0) mine = atom->type[index];
4485 else if (strcmp(word,"mol") == 0) {
4486 if (!atom->molecule_flag)
4487 error->one(FLERR,"Variable uses atom property that isn't allocated");
4488 mine = atom->molecule[index];
4489 }
4490 else if (strcmp(word,"x") == 0) mine = atom->x[index][0];
4491 else if (strcmp(word,"y") == 0) mine = atom->x[index][1];
4492 else if (strcmp(word,"z") == 0) mine = atom->x[index][2];
4493 else if (strcmp(word,"vx") == 0) mine = atom->v[index][0];
4494 else if (strcmp(word,"vy") == 0) mine = atom->v[index][1];
4495 else if (strcmp(word,"vz") == 0) mine = atom->v[index][2];
4496 else if (strcmp(word,"fx") == 0) mine = atom->f[index][0];
4497 else if (strcmp(word,"fy") == 0) mine = atom->f[index][1];
4498 else if (strcmp(word,"fz") == 0) mine = atom->f[index][2];
4499 else if (strcmp(word,"q") == 0) {
4500 if (!atom->q_flag)
4501 error->one(FLERR,"Variable uses atom property that isn't allocated");
4502 mine = atom->q[index];
4503 }
4504 else error->one(FLERR,"Invalid atom vector in variable formula");
4505
4506 } else mine = vector[index*nstride];
4507
4508 } else mine = 0.0;
4509
4510 double value;
4511 MPI_Allreduce(&mine,&value,1,MPI_DOUBLE,MPI_SUM,world);
4512
4513 if (tree) {
4514 Tree *newtree = new Tree();
4515 newtree->type = VALUE;
4516 newtree->value = value;
4517 treestack[ntreestack++] = newtree;
4518 } else argstack[nargstack++] = value;
4519 }
4520
4521 /* ----------------------------------------------------------------------
4522 check if word matches an atom vector
4523 return 1 if yes, else 0
4524 customize by adding an atom vector:
4525 id,mass,type,mol,x,y,z,vx,vy,vz,fx,fy,fz,q
4526 ------------------------------------------------------------------------- */
4527
is_atom_vector(char * word)4528 int Variable::is_atom_vector(char *word)
4529 {
4530 if (strcmp(word,"id") == 0) return 1;
4531 if (strcmp(word,"mass") == 0) return 1;
4532 if (strcmp(word,"type") == 0) return 1;
4533 if (strcmp(word,"mol") == 0) return 1;
4534 if (strcmp(word,"x") == 0) return 1;
4535 if (strcmp(word,"y") == 0) return 1;
4536 if (strcmp(word,"z") == 0) return 1;
4537 if (strcmp(word,"vx") == 0) return 1;
4538 if (strcmp(word,"vy") == 0) return 1;
4539 if (strcmp(word,"vz") == 0) return 1;
4540 if (strcmp(word,"fx") == 0) return 1;
4541 if (strcmp(word,"fy") == 0) return 1;
4542 if (strcmp(word,"fz") == 0) return 1;
4543 if (strcmp(word,"q") == 0) return 1;
4544 return 0;
4545 }
4546
4547 /* ----------------------------------------------------------------------
4548 process an atom vector in formula
4549 push result onto tree
4550 word = atom vector
4551 customize by adding an atom vector:
4552 id,mass,type,mol,x,y,z,vx,vy,vz,fx,fy,fz,q
4553 ------------------------------------------------------------------------- */
4554
atom_vector(char * word,Tree ** tree,Tree ** treestack,int & ntreestack)4555 void Variable::atom_vector(char *word, Tree **tree,
4556 Tree **treestack, int &ntreestack)
4557 {
4558 if (tree == nullptr)
4559 error->all(FLERR,"Atom vector in equal-style variable formula");
4560
4561 Tree *newtree = new Tree();
4562 newtree->type = ATOMARRAY;
4563 newtree->nstride = 3;
4564 treestack[ntreestack++] = newtree;
4565
4566 if (strcmp(word,"id") == 0) {
4567 if (sizeof(tagint) == sizeof(smallint)) {
4568 newtree->type = INTARRAY;
4569 newtree->iarray = (int *) atom->tag;
4570 } else {
4571 newtree->type = BIGINTARRAY;
4572 newtree->barray = (bigint *) atom->tag;
4573 }
4574 newtree->nstride = 1;
4575
4576 } else if (strcmp(word,"mass") == 0) {
4577 if (atom->rmass) {
4578 newtree->nstride = 1;
4579 newtree->array = atom->rmass;
4580 } else {
4581 newtree->type = TYPEARRAY;
4582 newtree->array = atom->mass;
4583 }
4584
4585 } else if (strcmp(word,"type") == 0) {
4586 newtree->type = INTARRAY;
4587 newtree->nstride = 1;
4588 newtree->iarray = atom->type;
4589
4590 } else if (strcmp(word,"mol") == 0) {
4591 if (!atom->molecule_flag)
4592 error->one(FLERR,"Variable uses atom property that isn't allocated");
4593 if (sizeof(tagint) == sizeof(smallint)) {
4594 newtree->type = INTARRAY;
4595 newtree->iarray = (int *) atom->molecule;
4596 } else {
4597 newtree->type = BIGINTARRAY;
4598 newtree->barray = (bigint *) atom->molecule;
4599 }
4600 newtree->nstride = 1;
4601 }
4602
4603 else if (strcmp(word,"x") == 0) newtree->array = &atom->x[0][0];
4604 else if (strcmp(word,"y") == 0) newtree->array = &atom->x[0][1];
4605 else if (strcmp(word,"z") == 0) newtree->array = &atom->x[0][2];
4606 else if (strcmp(word,"vx") == 0) newtree->array = &atom->v[0][0];
4607 else if (strcmp(word,"vy") == 0) newtree->array = &atom->v[0][1];
4608 else if (strcmp(word,"vz") == 0) newtree->array = &atom->v[0][2];
4609 else if (strcmp(word,"fx") == 0) newtree->array = &atom->f[0][0];
4610 else if (strcmp(word,"fy") == 0) newtree->array = &atom->f[0][1];
4611 else if (strcmp(word,"fz") == 0) newtree->array = &atom->f[0][2];
4612
4613 else if (strcmp(word,"q") == 0) {
4614 newtree->nstride = 1;
4615 newtree->array = atom->q;
4616 }
4617 }
4618
4619 /* ----------------------------------------------------------------------
4620 parse string for comma-separated args
4621 store copy of each arg in args array
4622 max allowed # of args = MAXFUNCARG
4623 ------------------------------------------------------------------------- */
4624
parse_args(char * str,char ** args)4625 int Variable::parse_args(char *str, char **args)
4626 {
4627 char *ptrnext;
4628 int narg = 0;
4629 char *ptr = str;
4630
4631 while (ptr && narg < MAXFUNCARG) {
4632 ptrnext = find_next_comma(ptr);
4633 if (ptrnext) *ptrnext = '\0';
4634 args[narg] = utils::strdup(utils::trim(ptr));
4635 narg++;
4636 ptr = ptrnext;
4637 if (ptr) ptr++;
4638 }
4639
4640 if (ptr) error->all(FLERR,"Too many args in variable function");
4641 return narg;
4642 }
4643
4644 /* ----------------------------------------------------------------------
4645 find next comma in str
4646 skip commas inside one or more nested parenthesis
4647 only return ptr to comma at level 0, else nullptr if not found
4648 ------------------------------------------------------------------------- */
4649
find_next_comma(char * str)4650 char *Variable::find_next_comma(char *str)
4651 {
4652 int level = 0;
4653 for (char *p = str; *p; ++p) {
4654 if ('(' == *p) level++;
4655 else if (')' == *p) level--;
4656 else if (',' == *p && !level) return p;
4657 }
4658 return nullptr;
4659 }
4660
4661
4662 /* ----------------------------------------------------------------------
4663 helper routine for printing variable name with error message
4664 ------------------------------------------------------------------------- */
4665
print_var_error(const std::string & srcfile,const int lineno,const std::string & errmsg,int ivar,int global)4666 void Variable::print_var_error(const std::string &srcfile, const int lineno,
4667 const std::string &errmsg, int ivar, int global)
4668 {
4669 if ((ivar >= 0) && (ivar < nvar)) {
4670 std::string msg = fmt::format("Variable {}: ",names[ivar]) + errmsg;
4671 if (global)
4672 error->all(srcfile,lineno,msg);
4673 else
4674 error->one(srcfile,lineno,msg);
4675 } else {
4676 if (global)
4677 error->all(srcfile,lineno,errmsg);
4678 else
4679 error->one(srcfile,lineno,errmsg);
4680 }
4681 }
4682
4683 /* ----------------------------------------------------------------------
4684 debug routine for printing formula tree recursively
4685 ------------------------------------------------------------------------- */
4686
print_tree(Tree * tree,int level)4687 void Variable::print_tree(Tree *tree, int level)
4688 {
4689 printf("TREE %d: %d %g\n",level,tree->type,tree->value);
4690 if (tree->first) print_tree(tree->first,level+1);
4691 if (tree->second) print_tree(tree->second,level+1);
4692 if (tree->nextra)
4693 for (int i = 0; i < tree->nextra; i++) print_tree(tree->extra[i],level+1);
4694 return;
4695 }
4696
4697 /* ----------------------------------------------------------------------
4698 recursive evaluation of string str
4699 called from "if" command in input script
4700 str is a boolean expression containing one or more items:
4701 number = 0.0, -5.45, 2.8e-4, ...
4702 math operation = (),x==y,x!=y,x<y,x<=y,x>y,x>=y,x&&y,x||y
4703 ------------------------------------------------------------------------- */
4704
evaluate_boolean(char * str)4705 double Variable::evaluate_boolean(char *str)
4706 {
4707 int op,opprevious,flag1,flag2;
4708 double value1,value2;
4709 char onechar;
4710 char *str1,*str2;
4711
4712 struct Arg {
4713 int flag; // 0 for numeric value, 1 for string
4714 double value; // stored numeric value
4715 char *str; // stored string
4716 };
4717
4718 Arg argstack[MAXLEVEL];
4719 int opstack[MAXLEVEL];
4720 int nargstack = 0;
4721 int nopstack = 0;
4722
4723 int i = 0;
4724 int expect = ARG;
4725
4726 while (1) {
4727 onechar = str[i];
4728
4729 // whitespace: just skip
4730
4731 if (isspace(onechar)) i++;
4732
4733 // ----------------
4734 // parentheses: recursively evaluate contents of parens
4735 // ----------------
4736
4737 else if (onechar == '(') {
4738 if (expect == OP)
4739 error->all(FLERR,"Invalid Boolean syntax in if command");
4740 expect = OP;
4741
4742 char *contents = nullptr;
4743 i = find_matching_paren(str,i,contents,-1);
4744 i++;
4745
4746 // evaluate contents and push on stack
4747
4748 argstack[nargstack].value = evaluate_boolean(contents);
4749 argstack[nargstack].flag = 0;
4750 nargstack++;
4751
4752 delete [] contents;
4753
4754 // ----------------
4755 // number: push value onto stack
4756 // ----------------
4757
4758 } else if (isdigit(onechar) || onechar == '.' || onechar == '-') {
4759 if (expect == OP)
4760 error->all(FLERR,"Invalid Boolean syntax in if command");
4761 expect = OP;
4762
4763 // set I to end of number, including scientific notation
4764
4765 int istart = i++;
4766 while (isdigit(str[i]) || str[i] == '.') i++;
4767 if (str[i] == 'e' || str[i] == 'E') {
4768 i++;
4769 if (str[i] == '+' || str[i] == '-') i++;
4770 while (isdigit(str[i])) i++;
4771 }
4772
4773 onechar = str[i];
4774 str[i] = '\0';
4775 argstack[nargstack].value = atof(&str[istart]);
4776 str[i] = onechar;
4777
4778 argstack[nargstack++].flag = 0;
4779
4780 // ----------------
4781 // string: push string onto stack
4782 // ----------------
4783
4784 } else if (isalpha(onechar)) {
4785 if (expect == OP)
4786 error->all(FLERR,"Invalid Boolean syntax in if command");
4787 expect = OP;
4788
4789 // set I to end of string
4790
4791 int istart = i++;
4792 while (isalnum(str[i]) || (str[i] == '_') || (str[i] == '/')) i++;
4793
4794 int n = i - istart + 1;
4795 argstack[nargstack].str = new char[n];
4796 onechar = str[i];
4797 str[i] = '\0';
4798 strcpy(argstack[nargstack].str,&str[istart]);
4799 str[i] = onechar;
4800
4801 argstack[nargstack++].flag = 1;
4802
4803 // ----------------
4804 // Boolean operator, including end-of-string
4805 // ----------------
4806
4807 } else if (strchr("<>=!&|\0",onechar)) {
4808 if (onechar == '=') {
4809 if (str[i+1] != '=')
4810 error->all(FLERR,"Invalid Boolean syntax in if command");
4811 op = EQ;
4812 i++;
4813 } else if (onechar == '!') {
4814 if (str[i+1] == '=') {
4815 op = NE;
4816 i++;
4817 } else op = NOT;
4818 } else if (onechar == '<') {
4819 if (str[i+1] != '=') op = LT;
4820 else {
4821 op = LE;
4822 i++;
4823 }
4824 } else if (onechar == '>') {
4825 if (str[i+1] != '=') op = GT;
4826 else {
4827 op = GE;
4828 i++;
4829 }
4830 } else if (onechar == '&') {
4831 if (str[i+1] != '&')
4832 error->all(FLERR,"Invalid Boolean syntax in if command");
4833 op = AND;
4834 i++;
4835 } else if (onechar == '|') {
4836 if (str[i+1] == '|') op = OR;
4837 else if (str[i+1] == '^') op = XOR;
4838 else error->all(FLERR,"Invalid Boolean syntax in if command");
4839 i++;
4840 } else op = DONE;
4841
4842 i++;
4843
4844 if (op == NOT && expect == ARG) {
4845 opstack[nopstack++] = op;
4846 continue;
4847 }
4848
4849 if (expect == ARG)
4850 error->all(FLERR,"Invalid Boolean syntax in if command");
4851 expect = ARG;
4852
4853 // evaluate stack as deep as possible while respecting precedence
4854 // before pushing current op onto stack
4855
4856 while (nopstack && precedence[opstack[nopstack-1]] >= precedence[op]) {
4857 opprevious = opstack[--nopstack];
4858
4859 nargstack--;
4860 flag2 = argstack[nargstack].flag;
4861 value2 = argstack[nargstack].value;
4862 str2 = argstack[nargstack].str;
4863 if (opprevious != NOT) {
4864 nargstack--;
4865 flag1 = argstack[nargstack].flag;
4866 value1 = argstack[nargstack].value;
4867 str1 = argstack[nargstack].str;
4868 }
4869
4870 if (opprevious == NOT) {
4871 if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
4872 if (value2 == 0.0) argstack[nargstack].value = 1.0;
4873 else argstack[nargstack].value = 0.0;
4874 } else if (opprevious == EQ) {
4875 if (flag1 != flag2)
4876 error->all(FLERR,"Invalid Boolean syntax in if command");
4877 if (flag2 == 0) {
4878 if (value1 == value2) argstack[nargstack].value = 1.0;
4879 else argstack[nargstack].value = 0.0;
4880 } else {
4881 if (strcmp(str1,str2) == 0) argstack[nargstack].value = 1.0;
4882 else argstack[nargstack].value = 0.0;
4883 delete [] str1;
4884 delete [] str2;
4885 }
4886 } else if (opprevious == NE) {
4887 if (flag1 != flag2)
4888 error->all(FLERR,"Invalid Boolean syntax in if command");
4889 if (flag2 == 0) {
4890 if (value1 != value2) argstack[nargstack].value = 1.0;
4891 else argstack[nargstack].value = 0.0;
4892 } else {
4893 if (strcmp(str1,str2) != 0) argstack[nargstack].value = 1.0;
4894 else argstack[nargstack].value = 0.0;
4895 delete [] str1;
4896 delete [] str2;
4897 }
4898 } else if (opprevious == LT) {
4899 if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
4900 if (value1 < value2) argstack[nargstack].value = 1.0;
4901 else argstack[nargstack].value = 0.0;
4902 } else if (opprevious == LE) {
4903 if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
4904 if (value1 <= value2) argstack[nargstack].value = 1.0;
4905 else argstack[nargstack].value = 0.0;
4906 } else if (opprevious == GT) {
4907 if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
4908 if (value1 > value2) argstack[nargstack].value = 1.0;
4909 else argstack[nargstack].value = 0.0;
4910 } else if (opprevious == GE) {
4911 if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
4912 if (value1 >= value2) argstack[nargstack].value = 1.0;
4913 else argstack[nargstack].value = 0.0;
4914 } else if (opprevious == AND) {
4915 if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
4916 if (value1 != 0.0 && value2 != 0.0) argstack[nargstack].value = 1.0;
4917 else argstack[nargstack].value = 0.0;
4918 } else if (opprevious == OR) {
4919 if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
4920 if (value1 != 0.0 || value2 != 0.0) argstack[nargstack].value = 1.0;
4921 else argstack[nargstack].value = 0.0;
4922 } else if (opprevious == XOR) {
4923 if (flag2) error->all(FLERR,"Invalid Boolean syntax in if command");
4924 if ((value1 == 0.0 && value2 != 0.0) ||
4925 (value1 != 0.0 && value2 == 0.0))
4926 argstack[nargstack].value = 1.0;
4927 else argstack[nargstack].value = 0.0;
4928 }
4929
4930 argstack[nargstack++].flag = 0;
4931 }
4932
4933 // if end-of-string, break out of entire formula evaluation loop
4934
4935 if (op == DONE) break;
4936
4937 // push current operation onto stack
4938
4939 opstack[nopstack++] = op;
4940
4941 } else error->all(FLERR,"Invalid Boolean syntax in if command");
4942 }
4943
4944 if (nopstack) error->all(FLERR,"Invalid Boolean syntax in if command");
4945 if (nargstack != 1) error->all(FLERR,"Invalid Boolean syntax in if command");
4946 return argstack[0].value;
4947 }
4948
4949 /* ----------------------------------------------------------------------
4950 class to read variable values from a file
4951 for flag = SCALARFILE, reads one value per line
4952 for flag = ATOMFILE, reads set of one value per atom
4953 ------------------------------------------------------------------------- */
4954
VarReader(LAMMPS * lmp,char * name,char * file,int flag)4955 VarReader::VarReader(LAMMPS *lmp, char *name, char *file, int flag) :
4956 Pointers(lmp)
4957 {
4958 me = comm->me;
4959 style = flag;
4960 fp = nullptr;
4961
4962 if (me == 0) {
4963 fp = fopen(file,"r");
4964 if (fp == nullptr)
4965 error->one(FLERR,"Cannot open file variable file {}: {}",
4966 file, utils::getsyserror());
4967 }
4968
4969 // if atomfile-style variable, must store per-atom values read from file
4970 // allocate a new fix STORE, so they persist
4971 // id = variable-ID + VARIABLE_STORE, fix group = all
4972
4973 fixstore = nullptr;
4974 id_fix = nullptr;
4975 buffer = nullptr;
4976
4977 if (style == Variable::ATOMFILE) {
4978 if (atom->map_style == Atom::MAP_NONE)
4979 error->all(FLERR,"Cannot use atomfile-style variable unless an atom map exists");
4980
4981 id_fix = utils::strdup(std::string(name) + "_VARIABLE_STORE");
4982 fixstore = (FixStore *) modify->add_fix(std::string(id_fix) + " all STORE peratom 0 1");
4983 buffer = new char[CHUNK*MAXLINE];
4984 }
4985 }
4986
4987 /* ---------------------------------------------------------------------- */
4988
~VarReader()4989 VarReader::~VarReader()
4990 {
4991 if (me == 0) {
4992 fclose(fp);
4993 fp = nullptr;
4994 }
4995
4996 // check modify in case all fixes have already been deleted
4997
4998 if (fixstore) {
4999 if (modify) modify->delete_fix(id_fix);
5000 delete [] id_fix;
5001 delete [] buffer;
5002 }
5003 }
5004
5005 /* ----------------------------------------------------------------------
5006 read for SCALARFILE style
5007 read next value from file into str for file-style variable
5008 strip comments, skip blank lines
5009 return 0 if successful, 1 if end-of-file
5010 ------------------------------------------------------------------------- */
5011
read_scalar(char * str)5012 int VarReader::read_scalar(char *str)
5013 {
5014 int n;
5015 char *ptr;
5016
5017 // read one string from file
5018
5019 if (me == 0) {
5020 while (1) {
5021 ptr = fgets(str,MAXLINE,fp);
5022 if (!ptr) { n=0; break; } // end of file
5023 ptr[strcspn(ptr,"#")] = '\0'; // strip comment
5024 ptr += strspn(ptr," \t\n\r\f"); // strip leading whitespace
5025 ptr[strcspn(ptr," \t\n\r\f")] = '\0'; // strip trailing whitespace
5026 n = strlen(ptr) + 1;
5027 if (n == 1) continue; // skip if blank line
5028 break;
5029 }
5030 if (n > 0) memmove(str,ptr,n); // move trimmed string back
5031 }
5032 MPI_Bcast(&n,1,MPI_INT,0,world);
5033 if (n == 0) return 1;
5034 MPI_Bcast(str,n,MPI_CHAR,0,world);
5035 return 0;
5036 }
5037
5038 /* ----------------------------------------------------------------------
5039 read snapshot of per-atom values from file
5040 into str for atomfile-style variable
5041 return 0 if successful, 1 if end-of-file
5042 ------------------------------------------------------------------------- */
5043
read_peratom()5044 int VarReader::read_peratom()
5045 {
5046 int i,m,n,nchunk,eof;
5047 tagint tag;
5048 char *ptr,*next;
5049 double value;
5050
5051 // set all per-atom values to 0.0
5052 // values that appear in file will overwrite this
5053
5054 double *vstore = fixstore->vstore;
5055
5056 int nlocal = atom->nlocal;
5057 for (i = 0; i < nlocal; i++) vstore[i] = 0.0;
5058
5059 // read one string from file, convert to Nlines
5060
5061 char str[MAXLINE];
5062 if (me == 0) {
5063 while (1) {
5064 ptr = fgets(str,MAXLINE,fp);
5065 if (!ptr) { n=0; break; } // end of file
5066 ptr[strcspn(ptr,"#")] = '\0'; // strip comment
5067 ptr += strspn(ptr," \t\n\r\f"); // strip leading whitespace
5068 ptr[strcspn(ptr," \t\n\r\f")] = '\0'; // strip trailing whitespace
5069 n = strlen(ptr) + 1;
5070 if (n == 1) continue; // skip if blank line
5071 break;
5072 }
5073 memmove(str,ptr,n); // move trimmed string back
5074 }
5075
5076 MPI_Bcast(&n,1,MPI_INT,0,world);
5077 if (n == 0) return 1;
5078 MPI_Bcast(str,n,MPI_CHAR,0,world);
5079 bigint nlines = utils::bnumeric(FLERR,str,false,lmp);
5080 tagint map_tag_max = atom->map_tag_max;
5081
5082 bigint nread = 0;
5083 while (nread < nlines) {
5084 nchunk = MIN(nlines-nread,CHUNK);
5085 eof = utils::read_lines_from_file(fp,nchunk,MAXLINE,buffer,me,world);
5086 if (eof) return 1;
5087
5088 char *buf = buffer;
5089 for (i = 0; i < nchunk; i++) {
5090 next = strchr(buf,'\n');
5091 *next = '\0';
5092 try {
5093 ValueTokenizer words(buf);
5094 tag = words.next_bigint();
5095 value = words.next_double();
5096 } catch (TokenizerException &e) {
5097 error->all(FLERR,"Invalid atomfile line '{}': {}",
5098 buf,e.what());
5099 }
5100 if ((tag <= 0) || (tag > map_tag_max))
5101 error->all(FLERR,"Invalid atom ID {} in variable "
5102 "file", tag);
5103 if ((m = atom->map(tag)) >= 0) vstore[m] = value;
5104 buf = next + 1;
5105 }
5106
5107 nread += nchunk;
5108 }
5109
5110 return 0;
5111 }
5112