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 "update.h"
16
17 #include "style_integrate.h" // IWYU pragma: keep
18 #include "style_minimize.h" // IWYU pragma: keep
19
20 #include "comm.h"
21 #include "compute.h"
22 #include "integrate.h"
23 #include "error.h"
24 #include "fix.h"
25 #include "force.h"
26 #include "min.h"
27 #include "modify.h"
28 #include "neighbor.h"
29 #include "output.h"
30
31 #include <cstring>
32
33 using namespace LAMMPS_NS;
34
35 /* ---------------------------------------------------------------------- */
36
Update(LAMMPS * lmp)37 Update::Update(LAMMPS *lmp) : Pointers(lmp)
38 {
39 char *str;
40
41 ntimestep = 0;
42 atime = 0.0;
43 atimestep = 0;
44 first_update = 0;
45
46 whichflag = 0;
47 firststep = laststep = 0;
48 beginstep = endstep = 0;
49 restrict_output = 0;
50 setupflag = 0;
51 post_integrate = 0;
52 multireplica = 0;
53
54 eflag_global = vflag_global = -1;
55 eflag_atom = vflag_atom = 0;
56
57 dt_default = 1;
58 dt = 0.0;
59 unit_style = nullptr;
60 set_units("lj");
61
62 integrate_style = nullptr;
63 integrate = nullptr;
64 minimize_style = nullptr;
65 minimize = nullptr;
66
67 integrate_map = new IntegrateCreatorMap();
68
69 #define INTEGRATE_CLASS
70 #define IntegrateStyle(key,Class) \
71 (*integrate_map)[#key] = &integrate_creator<Class>;
72 #include "style_integrate.h" // IWYU pragma: keep
73 #undef IntegrateStyle
74 #undef INTEGRATE_CLASS
75
76 minimize_map = new MinimizeCreatorMap();
77
78 #define MINIMIZE_CLASS
79 #define MinimizeStyle(key,Class) \
80 (*minimize_map)[#key] = &minimize_creator<Class>;
81 #include "style_minimize.h" // IWYU pragma: keep
82 #undef MinimizeStyle
83 #undef MINIMIZE_CLASS
84
85 str = (char *) "verlet";
86 create_integrate(1,&str,1);
87
88 str = (char *) "cg";
89 create_minimize(1,&str,1);
90 }
91
92 /* ---------------------------------------------------------------------- */
93
~Update()94 Update::~Update()
95 {
96 delete [] unit_style;
97
98 delete [] integrate_style;
99 delete integrate;
100
101 delete [] minimize_style;
102 delete minimize;
103
104 delete integrate_map;
105 delete minimize_map;
106 }
107
108 /* ---------------------------------------------------------------------- */
109
init()110 void Update::init()
111 {
112 // init the appropriate integrate and/or minimize class
113 // if neither (e.g. from write_restart) then just return
114
115 if (whichflag == 0) return;
116 if (whichflag == 1) integrate->init();
117 else if (whichflag == 2) minimize->init();
118
119 // only set first_update if a run or minimize is being performed
120
121 first_update = 1;
122 }
123
124 /* ---------------------------------------------------------------------- */
125
set_units(const char * style)126 void Update::set_units(const char *style)
127 {
128 // physical constants from:
129 // https://physics.nist.gov/cuu/Constants/Table/allascii.txt
130 // using thermochemical calorie = 4.184 J
131
132 double dt_old = dt;
133
134 if (strcmp(style,"lj") == 0) {
135 force->boltz = 1.0;
136 force->hplanck = 1.0;
137 force->mvv2e = 1.0;
138 force->ftm2v = 1.0;
139 force->mv2d = 1.0;
140 force->nktv2p = 1.0;
141 force->qqr2e = 1.0;
142 force->qe2f = 1.0;
143 force->vxmu2f = 1.0;
144 force->xxt2kmu = 1.0;
145 force->e_mass = 0.0; // not yet set
146 force->hhmrr2e = 0.0;
147 force->mvh2r = 0.0;
148 force->angstrom = 1.0;
149 force->femtosecond = 1.0;
150 force->qelectron = 1.0;
151
152 dt = 0.005;
153 neighbor->skin = 0.3;
154
155 } else if (strcmp(style,"real") == 0) {
156 force->boltz = 0.0019872067;
157 force->hplanck = 95.306976368;
158 force->mvv2e = 48.88821291 * 48.88821291;
159 force->ftm2v = 1.0 / 48.88821291 / 48.88821291;
160 force->mv2d = 1.0 / 0.602214129;
161 force->nktv2p = 68568.415;
162 force->qqr2e = 332.06371; // see also force->qqr2d_lammps_real
163 force->qe2f = 23.060549;
164 force->vxmu2f = 1.4393264316e4;
165 force->xxt2kmu = 0.1;
166 force->e_mass = 1.0/1836.1527556560675;
167 force->hhmrr2e = 0.0957018663603261;
168 force->mvh2r = 1.5339009481951;
169 force->angstrom = 1.0;
170 force->femtosecond = 1.0;
171 force->qelectron = 1.0;
172
173 dt = 1.0;
174 neighbor->skin = 2.0;
175
176 } else if (strcmp(style,"metal") == 0) {
177 force->boltz = 8.617343e-5;
178 force->hplanck = 4.135667403e-3;
179 force->mvv2e = 1.0364269e-4;
180 force->ftm2v = 1.0 / 1.0364269e-4;
181 force->mv2d = 1.0 / 0.602214129;
182 force->nktv2p = 1.6021765e6;
183 force->qqr2e = 14.399645;
184 force->qe2f = 1.0;
185 force->vxmu2f = 0.6241509647;
186 force->xxt2kmu = 1.0e-4;
187 force->e_mass = 0.0; // not yet set
188 force->hhmrr2e = 0.0;
189 force->mvh2r = 0.0;
190 force->angstrom = 1.0;
191 force->femtosecond = 1.0e-3;
192 force->qelectron = 1.0;
193
194 dt = 0.001;
195 neighbor->skin = 2.0;
196
197 } else if (strcmp(style,"si") == 0) {
198 force->boltz = 1.3806504e-23;
199 force->hplanck = 6.62606896e-34;
200 force->mvv2e = 1.0;
201 force->ftm2v = 1.0;
202 force->mv2d = 1.0;
203 force->nktv2p = 1.0;
204 force->qqr2e = 8.9876e9;
205 force->qe2f = 1.0;
206 force->vxmu2f = 1.0;
207 force->xxt2kmu = 1.0;
208 force->e_mass = 0.0; // not yet set
209 force->hhmrr2e = 0.0;
210 force->mvh2r = 0.0;
211 force->angstrom = 1.0e-10;
212 force->femtosecond = 1.0e-15;
213 force->qelectron = 1.6021765e-19;
214
215 dt = 1.0e-8;
216 neighbor->skin = 0.001;
217
218 } else if (strcmp(style,"cgs") == 0) {
219 force->boltz = 1.3806504e-16;
220 force->hplanck = 6.62606896e-27;
221 force->mvv2e = 1.0;
222 force->ftm2v = 1.0;
223 force->mv2d = 1.0;
224 force->nktv2p = 1.0;
225 force->qqr2e = 1.0;
226 force->qe2f = 1.0;
227 force->vxmu2f = 1.0;
228 force->xxt2kmu = 1.0;
229 force->e_mass = 0.0; // not yet set
230 force->hhmrr2e = 0.0;
231 force->mvh2r = 0.0;
232 force->angstrom = 1.0e-8;
233 force->femtosecond = 1.0e-15;
234 force->qelectron = 4.8032044e-10;
235
236 dt = 1.0e-8;
237 neighbor->skin = 0.1;
238
239 } else if (strcmp(style,"electron") == 0) {
240 force->boltz = 3.16681534e-6;
241 force->hplanck = 0.1519829846;
242 force->mvv2e = 1.06657236;
243 force->ftm2v = 0.937582899;
244 force->mv2d = 1.0;
245 force->nktv2p = 2.94210108e13;
246 force->qqr2e = 1.0;
247 force->qe2f = 1.94469051e-10;
248 force->vxmu2f = 3.39893149e1;
249 force->xxt2kmu = 3.13796367e-2;
250 force->e_mass = 0.0; // not yet set
251 force->hhmrr2e = 0.0;
252 force->mvh2r = 0.0;
253 force->angstrom = 1.88972612;
254 force->femtosecond = 1.0;
255 force->qelectron = 1.0;
256
257 dt = 0.001;
258 neighbor->skin = 2.0;
259
260 } else if (strcmp(style,"micro") == 0) {
261 force->boltz = 1.3806504e-8;
262 force->hplanck = 6.62606896e-13;
263 force->mvv2e = 1.0;
264 force->ftm2v = 1.0;
265 force->mv2d = 1.0;
266 force->nktv2p = 1.0;
267 force->qqr2e = 8.987556e6;
268 force->qe2f = 1.0;
269 force->vxmu2f = 1.0;
270 force->xxt2kmu = 1.0;
271 force->e_mass = 0.0; // not yet set
272 force->hhmrr2e = 0.0;
273 force->mvh2r = 0.0;
274 force->angstrom = 1.0e-4;
275 force->femtosecond = 1.0e-9;
276 force->qelectron = 1.6021765e-7;
277
278 dt = 2.0;
279 neighbor->skin = 0.1;
280
281 } else if (strcmp(style,"nano") == 0) {
282 force->boltz = 0.013806504;
283 force->hplanck = 6.62606896e-4;
284 force->mvv2e = 1.0;
285 force->ftm2v = 1.0;
286 force->mv2d = 1.0;
287 force->nktv2p = 1.0;
288 force->qqr2e = 230.7078669;
289 force->qe2f = 1.0;
290 force->vxmu2f = 1.0;
291 force->xxt2kmu = 1.0;
292 force->e_mass = 0.0; // not yet set
293 force->hhmrr2e = 0.0;
294 force->mvh2r = 0.0;
295 force->angstrom = 1.0e-1;
296 force->femtosecond = 1.0e-6;
297 force->qelectron = 1.0;
298
299 dt = 0.00045;
300 neighbor->skin = 0.1;
301
302 } else error->all(FLERR,"Illegal units command");
303
304 delete [] unit_style;
305 unit_style = utils::strdup(style);
306
307 // check if timestep was changed from default value
308 if (!dt_default && (comm->me == 0)) {
309 error->warning(FLERR,"Changing timestep from {:.6} to {:.6} due to "
310 "changing units to {}", dt_old, dt, unit_style);
311 }
312 dt_default = 1;
313 }
314
315 /* ---------------------------------------------------------------------- */
316
create_integrate(int narg,char ** arg,int trysuffix)317 void Update::create_integrate(int narg, char **arg, int trysuffix)
318 {
319 if (narg < 1) error->all(FLERR,"Illegal run_style command");
320
321 delete [] integrate_style;
322 delete integrate;
323
324 int sflag;
325
326 if (narg-1 > 0) {
327 new_integrate(arg[0],narg-1,&arg[1],trysuffix,sflag);
328 } else {
329 new_integrate(arg[0],0,nullptr,trysuffix,sflag);
330 }
331
332 std::string estyle = arg[0];
333 if (sflag) {
334 estyle += "/";
335 if (sflag == 1) estyle += lmp->suffix;
336 else estyle += lmp->suffix2;
337 }
338 integrate_style = utils::strdup(estyle);
339 }
340
341 /* ----------------------------------------------------------------------
342 create the Integrate style, first with suffix appended
343 ------------------------------------------------------------------------- */
344
new_integrate(char * style,int narg,char ** arg,int trysuffix,int & sflag)345 void Update::new_integrate(char *style, int narg, char **arg,
346 int trysuffix, int &sflag)
347 {
348 if (trysuffix && lmp->suffix_enable) {
349 if (lmp->suffix) {
350 sflag = 1;
351 std::string estyle = style + std::string("/") + lmp->suffix;
352 if (integrate_map->find(estyle) != integrate_map->end()) {
353 IntegrateCreator &integrate_creator = (*integrate_map)[estyle];
354 integrate = integrate_creator(lmp, narg, arg);
355 return;
356 }
357 }
358
359 if (lmp->suffix2) {
360 sflag = 2;
361 std::string estyle = style + std::string("/") + lmp->suffix2;
362 if (integrate_map->find(estyle) != integrate_map->end()) {
363 IntegrateCreator &integrate_creator = (*integrate_map)[estyle];
364 integrate = integrate_creator(lmp, narg, arg);
365 return;
366 }
367 }
368 }
369
370 sflag = 0;
371 if (integrate_map->find(style) != integrate_map->end()) {
372 IntegrateCreator &integrate_creator = (*integrate_map)[style];
373 integrate = integrate_creator(lmp, narg, arg);
374 return;
375 }
376
377 error->all(FLERR,"Illegal integrate style");
378 }
379
380 /* ----------------------------------------------------------------------
381 one instance per integrate style in style_integrate.h
382 ------------------------------------------------------------------------- */
383
384 template <typename T>
integrate_creator(LAMMPS * lmp,int narg,char ** arg)385 Integrate *Update::integrate_creator(LAMMPS *lmp, int narg, char ** arg)
386 {
387 return new T(lmp, narg, arg);
388 }
389
390 /* ---------------------------------------------------------------------- */
391
create_minimize(int narg,char ** arg,int trysuffix)392 void Update::create_minimize(int narg, char **arg, int trysuffix)
393 {
394 if (narg < 1) error->all(FLERR,"Illegal run_style command");
395
396 delete [] minimize_style;
397 delete minimize;
398
399 int sflag;
400 new_minimize(arg[0],narg-1,&arg[1],trysuffix,sflag);
401
402 std::string estyle = arg[0];
403 if (sflag) {
404 estyle += "/";
405 if (sflag == 1) estyle += lmp->suffix;
406 else estyle += lmp->suffix2;
407 }
408 minimize_style = utils::strdup(estyle);
409 }
410
411 /* ----------------------------------------------------------------------
412 create the Minimize style, first with suffix appended
413 ------------------------------------------------------------------------- */
414
new_minimize(char * style,int,char **,int trysuffix,int & sflag)415 void Update::new_minimize(char *style, int /* narg */, char ** /* arg */,
416 int trysuffix, int &sflag)
417 {
418 if (trysuffix && lmp->suffix_enable) {
419 if (lmp->suffix) {
420 sflag = 1;
421 std::string estyle = style + std::string("/") + lmp->suffix;
422 if (minimize_map->find(estyle) != minimize_map->end()) {
423 MinimizeCreator &minimize_creator = (*minimize_map)[estyle];
424 minimize = minimize_creator(lmp);
425 return;
426 }
427 }
428
429 if (lmp->suffix2) {
430 sflag = 2;
431 std::string estyle = style + std::string("/") + lmp->suffix2;
432 if (minimize_map->find(estyle) != minimize_map->end()) {
433 MinimizeCreator &minimize_creator = (*minimize_map)[estyle];
434 minimize = minimize_creator(lmp);
435 return;
436 }
437 }
438 }
439
440 sflag = 0;
441 if (minimize_map->find(style) != minimize_map->end()) {
442 MinimizeCreator &minimize_creator = (*minimize_map)[style];
443 minimize = minimize_creator(lmp);
444 return;
445 }
446
447 error->all(FLERR,"Illegal minimize style");
448 }
449
450 /* ----------------------------------------------------------------------
451 one instance per minimize style in style_minimize.h
452 ------------------------------------------------------------------------- */
453
454 template <typename T>
minimize_creator(LAMMPS * lmp)455 Min *Update::minimize_creator(LAMMPS *lmp)
456 {
457 return new T(lmp);
458 }
459
460 /* ----------------------------------------------------------------------
461 reset timestep as called from input script
462 ------------------------------------------------------------------------- */
463
reset_timestep(int narg,char ** arg)464 void Update::reset_timestep(int narg, char **arg)
465 {
466 if (narg != 1) error->all(FLERR,"Illegal reset_timestep command");
467 bigint newstep = utils::bnumeric(FLERR,arg[0],false,lmp);
468 reset_timestep(newstep);
469 }
470
471 /* ----------------------------------------------------------------------
472 reset timestep
473 called from rerun command and input script (indirectly)
474 ------------------------------------------------------------------------- */
475
reset_timestep(bigint newstep)476 void Update::reset_timestep(bigint newstep)
477 {
478 ntimestep = newstep;
479 if (ntimestep < 0) error->all(FLERR,"Timestep must be >= 0");
480
481 // set atimestep to new timestep
482 // so future update_time() calls will be correct
483
484 atimestep = ntimestep;
485
486 // trigger reset of timestep for output
487 // do not allow any timestep-dependent fixes to be already defined
488
489 output->reset_timestep(ntimestep);
490
491 for (int i = 0; i < modify->nfix; i++) {
492 if (modify->fix[i]->time_depend)
493 error->all(FLERR,
494 "Cannot reset timestep with a time-dependent fix defined");
495 }
496
497 // reset eflag/vflag global so no commands will think eng/virial are current
498
499 eflag_global = vflag_global = -1;
500
501 // reset invoked flags of computes,
502 // so no commands will think they are current between runs
503
504 for (int i = 0; i < modify->ncompute; i++) {
505 modify->compute[i]->invoked_scalar = -1;
506 modify->compute[i]->invoked_vector = -1;
507 modify->compute[i]->invoked_array = -1;
508 modify->compute[i]->invoked_peratom = -1;
509 modify->compute[i]->invoked_local = -1;
510 }
511
512 // clear timestep list of computes that store future invocation times
513
514 for (int i = 0; i < modify->ncompute; i++)
515 if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
516
517 // Neighbor Bin/Stencil/Pair classes store timestamps that need to be cleared
518
519 neighbor->reset_timestep(ntimestep);
520
521 // NOTE: 7Jun12, adding rerun command, don't think this is required
522
523 //for (int i = 0; i < domain->nregion; i++)
524 // if (domain->regions[i]->dynamic_check())
525 // error->all(FLERR,"Cannot reset timestep with a dynamic region defined");
526 }
527
528 /* ----------------------------------------------------------------------
529 update elapsed simulation time
530 called at end of runs or when timestep size changes
531 ------------------------------------------------------------------------- */
532
update_time()533 void Update::update_time()
534 {
535 atime += (ntimestep-atimestep) * dt;
536 atimestep = ntimestep;
537 }
538
539 /* ----------------------------------------------------------------------
540 memory usage of update and integrate/minimize
541 ------------------------------------------------------------------------- */
542
memory_usage()543 double Update::memory_usage()
544 {
545 double bytes = 0;
546 if (whichflag == 1) bytes += integrate->memory_usage();
547 else if (whichflag == 2) bytes += minimize->memory_usage();
548 return bytes;
549 }
550