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 "modify.h"
16 #include "style_compute.h" // IWYU pragma: keep
17 #include "style_fix.h" // IWYU pragma: keep
18
19 #include "atom.h"
20 #include "comm.h"
21 #include "compute.h" // IWYU pragma: keep
22 #include "domain.h"
23 #include "error.h"
24 #include "fix.h" // IWYU pragma: keep
25 #include "group.h"
26 #include "input.h"
27 #include "memory.h"
28 #include "region.h"
29 #include "update.h"
30 #include "variable.h"
31
32 #include <cstring>
33
34 using namespace LAMMPS_NS;
35 using namespace FixConst;
36
37 #define DELTA 4
38 #define BIG 1.0e20
39
40 /* ---------------------------------------------------------------------- */
41
Modify(LAMMPS * lmp)42 Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
43 {
44 nfix = maxfix = 0;
45 n_initial_integrate = n_post_integrate = 0;
46 n_pre_exchange = n_pre_neighbor = n_post_neighbor = 0;
47 n_pre_force = n_pre_reverse = n_post_force = 0;
48 n_final_integrate = n_end_of_step = 0;
49 n_energy_couple = n_energy_global = n_energy_atom = 0;
50 n_initial_integrate_respa = n_post_integrate_respa = 0;
51 n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0;
52 n_min_pre_exchange = n_min_pre_force = n_min_pre_reverse = 0;
53 n_min_post_force = n_min_energy = 0;
54 n_timeflag = -1;
55
56 fix = nullptr;
57 fmask = nullptr;
58 list_initial_integrate = list_post_integrate = nullptr;
59 list_pre_exchange = list_pre_neighbor = list_post_neighbor = nullptr;
60 list_pre_force = list_pre_reverse = list_post_force = nullptr;
61 list_final_integrate = list_end_of_step = nullptr;
62 list_energy_couple = list_energy_global = list_energy_atom = nullptr;
63 list_initial_integrate_respa = list_post_integrate_respa = nullptr;
64 list_pre_force_respa = list_post_force_respa = nullptr;
65 list_final_integrate_respa = nullptr;
66 list_min_pre_exchange = list_min_pre_neighbor = list_min_post_neighbor = nullptr;
67 list_min_pre_force = list_min_pre_reverse = list_min_post_force = nullptr;
68 list_min_energy = nullptr;
69
70 end_of_step_every = nullptr;
71
72 list_timeflag = nullptr;
73
74 nfix_restart_global = 0;
75 id_restart_global = style_restart_global = nullptr;
76 state_restart_global = nullptr;
77 used_restart_global = nullptr;
78 nfix_restart_peratom = 0;
79 id_restart_peratom = style_restart_peratom = nullptr;
80 index_restart_peratom = used_restart_peratom = nullptr;
81
82 ncompute = maxcompute = 0;
83 compute = nullptr;
84
85 create_factories();
86 }
87
create_factories()88 void _noopt Modify::create_factories()
89 {
90 // fill map with fixes listed in style_fix.h
91
92 fix_map = new FixCreatorMap();
93
94 #define FIX_CLASS
95 #define FixStyle(key,Class) \
96 (*fix_map)[#key] = &fix_creator<Class>;
97 #include "style_fix.h" // IWYU pragma: keep
98 #undef FixStyle
99 #undef FIX_CLASS
100
101 // fill map with computes listed in style_compute.h
102
103 compute_map = new ComputeCreatorMap();
104
105 #define COMPUTE_CLASS
106 #define ComputeStyle(key,Class) \
107 (*compute_map)[#key] = &compute_creator<Class>;
108 #include "style_compute.h" // IWYU pragma: keep
109 #undef ComputeStyle
110 #undef COMPUTE_CLASS
111 }
112
113 /* ---------------------------------------------------------------------- */
114
~Modify()115 Modify::~Modify()
116 {
117 // delete all fixes
118 // do it via delete_fix() so callbacks in Atom are also updated correctly
119
120 while (nfix) delete_fix(0);
121 memory->sfree(fix);
122 memory->destroy(fmask);
123
124 // delete all computes
125
126 for (int i = 0; i < ncompute; i++) delete compute[i];
127 memory->sfree(compute);
128
129 delete [] list_initial_integrate;
130 delete [] list_post_integrate;
131 delete [] list_pre_exchange;
132 delete [] list_pre_neighbor;
133 delete [] list_post_neighbor;
134 delete [] list_pre_force;
135 delete [] list_pre_reverse;
136 delete [] list_post_force;
137 delete [] list_final_integrate;
138 delete [] list_end_of_step;
139 delete [] list_energy_couple;
140 delete [] list_energy_global;
141 delete [] list_energy_atom;
142 delete [] list_initial_integrate_respa;
143 delete [] list_post_integrate_respa;
144 delete [] list_pre_force_respa;
145 delete [] list_post_force_respa;
146 delete [] list_final_integrate_respa;
147 delete [] list_min_pre_exchange;
148 delete [] list_min_pre_neighbor;
149 delete [] list_min_post_neighbor;
150 delete [] list_min_pre_force;
151 delete [] list_min_pre_reverse;
152 delete [] list_min_post_force;
153 delete [] list_min_energy;
154
155 delete [] end_of_step_every;
156 delete [] list_timeflag;
157
158 restart_deallocate(0);
159
160 delete compute_map;
161 delete fix_map;
162 }
163
164 /* ----------------------------------------------------------------------
165 initialize all fixes and computes
166 ------------------------------------------------------------------------- */
167
init()168 void Modify::init()
169 {
170 int i,j;
171
172 // delete storage of restart info since it is not valid after 1st run
173
174 restart_deallocate(1);
175
176 // init each compute
177 // set invoked_scalar,vector,etc to -1 to force new run to re-compute them
178 // add initial timestep to all computes that store invocation times
179 // since any of them may be invoked by initial thermo
180 // do not clear out invocation times stored within a compute,
181 // b/c some may be holdovers from previous run, like for ave fixes
182
183 for (i = 0; i < ncompute; i++) {
184 compute[i]->init();
185 compute[i]->invoked_scalar = -1;
186 compute[i]->invoked_vector = -1;
187 compute[i]->invoked_array = -1;
188 compute[i]->invoked_peratom = -1;
189 compute[i]->invoked_local = -1;
190 }
191 addstep_compute_all(update->ntimestep);
192
193 // init each fix
194 // should not need to come before compute init
195 // used to b/c temperature computes called fix->dof() in their init,
196 // and fix rigid required its own init before its dof() could be called,
197 // but computes now do their DOF in setup()
198
199 for (i = 0; i < nfix; i++) fix[i]->init();
200
201 // set global flag if any fix has its restart_pbc flag set
202
203 restart_pbc_any = 0;
204 for (i = 0; i < nfix; i++)
205 if (fix[i]->restart_pbc) restart_pbc_any = 1;
206
207 // create lists of fixes to call at each stage of run
208 // needs to happen after init() of computes
209 // b/c a compute::init() can delete a fix, e.g. compute chunk/atom
210
211 list_init(INITIAL_INTEGRATE,n_initial_integrate,list_initial_integrate);
212 list_init(POST_INTEGRATE,n_post_integrate,list_post_integrate);
213 list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange);
214 list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor);
215 list_init(POST_NEIGHBOR,n_post_neighbor,list_post_neighbor);
216 list_init(PRE_FORCE,n_pre_force,list_pre_force);
217 list_init(PRE_REVERSE,n_pre_reverse,list_pre_reverse);
218 list_init(POST_FORCE,n_post_force,list_post_force);
219 list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate);
220 list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step);
221 list_init_energy_couple(n_energy_couple,list_energy_couple);
222 list_init_energy_global(n_energy_global,list_energy_global);
223 list_init_energy_atom(n_energy_atom,list_energy_atom);
224
225 list_init(INITIAL_INTEGRATE_RESPA,n_initial_integrate_respa,list_initial_integrate_respa);
226 list_init(POST_INTEGRATE_RESPA,n_post_integrate_respa,list_post_integrate_respa);
227 list_init(POST_FORCE_RESPA,n_post_force_respa,list_post_force_respa);
228 list_init(PRE_FORCE_RESPA,n_pre_force_respa,list_pre_force_respa);
229 list_init(FINAL_INTEGRATE_RESPA,n_final_integrate_respa,list_final_integrate_respa);
230
231 list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange);
232 list_init(MIN_PRE_NEIGHBOR,n_min_pre_neighbor,list_min_pre_neighbor);
233 list_init(MIN_POST_NEIGHBOR,n_min_post_neighbor,list_min_post_neighbor);
234 list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force);
235 list_init(MIN_PRE_REVERSE,n_min_pre_reverse,list_min_pre_reverse);
236 list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
237 list_init(MIN_ENERGY,n_min_energy,list_min_energy);
238
239 // create list of computes that store invocation times
240
241 list_init_compute();
242
243 // error if any fix or compute is using a dynamic group when not allowed
244
245 for (i = 0; i < nfix; i++)
246 if (!fix[i]->dynamic_group_allow && group->dynamic[fix[i]->igroup])
247 error->all(FLERR,"Fix {} does not allow use with a "
248 "dynamic group",fix[i]->id);
249
250 for (i = 0; i < ncompute; i++)
251 if (!compute[i]->dynamic_group_allow &&
252 group->dynamic[compute[i]->igroup])
253 error->all(FLERR,"Compute {} does not allow use with a "
254 "dynamic group",compute[i]->id);
255
256 // warn if any particle is time integrated more than once
257
258 int nlocal = atom->nlocal;
259 int *mask = atom->mask;
260
261 int *flag = new int[nlocal];
262 for (i = 0; i < nlocal; i++) flag[i] = 0;
263
264 int groupbit;
265 for (i = 0; i < nfix; i++) {
266 if (fix[i]->time_integrate == 0) continue;
267 groupbit = fix[i]->groupbit;
268 for (j = 0; j < nlocal; j++)
269 if (mask[j] & groupbit) flag[j]++;
270 }
271
272 int check = 0;
273 for (i = 0; i < nlocal; i++)
274 if (flag[i] > 1) check = 1;
275
276 delete [] flag;
277
278 int checkall;
279 MPI_Allreduce(&check,&checkall,1,MPI_INT,MPI_SUM,world);
280 if (comm->me == 0 && checkall)
281 error->warning(FLERR,
282 "One or more atoms are time integrated more than once");
283 }
284
285 /* ----------------------------------------------------------------------
286 setup for run, calls setup() of all fixes and computes
287 called from Verlet, RESPA, Min
288 ------------------------------------------------------------------------- */
289
setup(int vflag)290 void Modify::setup(int vflag)
291 {
292 // compute setup needs to come before fix setup
293 // b/c NH fixes need DOF of temperature computes
294 // fix group setup() is special case since populates a dynamic group
295 // needs to be done before temperature compute setup
296
297 for (int i = 0; i < nfix; i++)
298 if (strcmp(fix[i]->style,"GROUP") == 0) fix[i]->setup(vflag);
299
300 for (int i = 0; i < ncompute; i++) compute[i]->setup();
301
302 if (update->whichflag == 1)
303 for (int i = 0; i < nfix; i++) fix[i]->setup(vflag);
304 else if (update->whichflag == 2)
305 for (int i = 0; i < nfix; i++) fix[i]->min_setup(vflag);
306 }
307
308 /* ----------------------------------------------------------------------
309 setup pre_exchange call, only for fixes that define pre_exchange
310 called from Verlet, RESPA, Min, and WriteRestart with whichflag = 0
311 ------------------------------------------------------------------------- */
312
setup_pre_exchange()313 void Modify::setup_pre_exchange()
314 {
315 if (update->whichflag <= 1)
316 for (int i = 0; i < n_pre_exchange; i++)
317 fix[list_pre_exchange[i]]->setup_pre_exchange();
318 else if (update->whichflag == 2)
319 for (int i = 0; i < n_min_pre_exchange; i++)
320 fix[list_min_pre_exchange[i]]->setup_pre_exchange();
321 }
322
323 /* ----------------------------------------------------------------------
324 setup pre_neighbor call, only for fixes that define pre_neighbor
325 called from Verlet, RESPA
326 ------------------------------------------------------------------------- */
327
setup_pre_neighbor()328 void Modify::setup_pre_neighbor()
329 {
330 if (update->whichflag == 1)
331 for (int i = 0; i < n_pre_neighbor; i++)
332 fix[list_pre_neighbor[i]]->setup_pre_neighbor();
333 else if (update->whichflag == 2)
334 for (int i = 0; i < n_min_pre_neighbor; i++)
335 fix[list_min_pre_neighbor[i]]->setup_pre_neighbor();
336 }
337
338 /* ----------------------------------------------------------------------
339 setup post_neighbor call, only for fixes that define post_neighbor
340 called from Verlet, RESPA
341 ------------------------------------------------------------------------- */
342
setup_post_neighbor()343 void Modify::setup_post_neighbor()
344 {
345 if (update->whichflag == 1)
346 for (int i = 0; i < n_post_neighbor; i++)
347 fix[list_post_neighbor[i]]->setup_post_neighbor();
348 else if (update->whichflag == 2)
349 for (int i = 0; i < n_min_post_neighbor; i++)
350 fix[list_min_post_neighbor[i]]->setup_post_neighbor();
351 }
352
353 /* ----------------------------------------------------------------------
354 setup pre_force call, only for fixes that define pre_force
355 called from Verlet, RESPA, Min
356 ------------------------------------------------------------------------- */
357
setup_pre_force(int vflag)358 void Modify::setup_pre_force(int vflag)
359 {
360 if (update->whichflag == 1)
361 for (int i = 0; i < n_pre_force; i++)
362 fix[list_pre_force[i]]->setup_pre_force(vflag);
363 else if (update->whichflag == 2)
364 for (int i = 0; i < n_min_pre_force; i++)
365 fix[list_min_pre_force[i]]->setup_pre_force(vflag);
366 }
367
368 /* ----------------------------------------------------------------------
369 setup pre_reverse call, only for fixes that define pre_reverse
370 called from Verlet, RESPA, Min
371 ------------------------------------------------------------------------- */
372
setup_pre_reverse(int eflag,int vflag)373 void Modify::setup_pre_reverse(int eflag, int vflag)
374 {
375 if (update->whichflag == 1)
376 for (int i = 0; i < n_pre_reverse; i++)
377 fix[list_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
378 else if (update->whichflag == 2)
379 for (int i = 0; i < n_min_pre_reverse; i++)
380 fix[list_min_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
381 }
382
383 /* ----------------------------------------------------------------------
384 1st half of integrate call, only for relevant fixes
385 ------------------------------------------------------------------------- */
386
initial_integrate(int vflag)387 void Modify::initial_integrate(int vflag)
388 {
389 for (int i = 0; i < n_initial_integrate; i++)
390 fix[list_initial_integrate[i]]->initial_integrate(vflag);
391 }
392
393 /* ----------------------------------------------------------------------
394 post_integrate call, only for relevant fixes
395 ------------------------------------------------------------------------- */
396
post_integrate()397 void Modify::post_integrate()
398 {
399 for (int i = 0; i < n_post_integrate; i++)
400 fix[list_post_integrate[i]]->post_integrate();
401 }
402
403 /* ----------------------------------------------------------------------
404 pre_exchange call, only for relevant fixes
405 ------------------------------------------------------------------------- */
406
pre_exchange()407 void Modify::pre_exchange()
408 {
409 for (int i = 0; i < n_pre_exchange; i++)
410 fix[list_pre_exchange[i]]->pre_exchange();
411 }
412
413 /* ----------------------------------------------------------------------
414 pre_neighbor call, only for relevant fixes
415 ------------------------------------------------------------------------- */
416
pre_neighbor()417 void Modify::pre_neighbor()
418 {
419 for (int i = 0; i < n_pre_neighbor; i++)
420 fix[list_pre_neighbor[i]]->pre_neighbor();
421 }
422
423 /* ----------------------------------------------------------------------
424 post_neighbor call, only for relevant fixes
425 ------------------------------------------------------------------------- */
426
post_neighbor()427 void Modify::post_neighbor()
428 {
429 for (int i = 0; i < n_post_neighbor; i++)
430 fix[list_post_neighbor[i]]->post_neighbor();
431 }
432
433 /* ----------------------------------------------------------------------
434 pre_force call, only for relevant fixes
435 ------------------------------------------------------------------------- */
436
pre_force(int vflag)437 void Modify::pre_force(int vflag)
438 {
439 for (int i = 0; i < n_pre_force; i++)
440 fix[list_pre_force[i]]->pre_force(vflag);
441 }
442 /* ----------------------------------------------------------------------
443 pre_reverse call, only for relevant fixes
444 ------------------------------------------------------------------------- */
445
pre_reverse(int eflag,int vflag)446 void Modify::pre_reverse(int eflag, int vflag)
447 {
448 for (int i = 0; i < n_pre_reverse; i++)
449 fix[list_pre_reverse[i]]->pre_reverse(eflag,vflag);
450 }
451
452 /* ----------------------------------------------------------------------
453 post_force call, only for relevant fixes
454 ------------------------------------------------------------------------- */
455
post_force(int vflag)456 void Modify::post_force(int vflag)
457 {
458 for (int i = 0; i < n_post_force; i++)
459 fix[list_post_force[i]]->post_force(vflag);
460 }
461
462 /* ----------------------------------------------------------------------
463 2nd half of integrate call, only for relevant fixes
464 ------------------------------------------------------------------------- */
465
final_integrate()466 void Modify::final_integrate()
467 {
468 for (int i = 0; i < n_final_integrate; i++)
469 fix[list_final_integrate[i]]->final_integrate();
470 }
471
472 /* ----------------------------------------------------------------------
473 end-of-timestep call, only for relevant fixes
474 only call fix->end_of_step() on timesteps that are multiples of nevery
475 ------------------------------------------------------------------------- */
476
end_of_step()477 void Modify::end_of_step()
478 {
479 for (int i = 0; i < n_end_of_step; i++)
480 if (update->ntimestep % end_of_step_every[i] == 0)
481 fix[list_end_of_step[i]]->end_of_step();
482 }
483
484 /* ----------------------------------------------------------------------
485 coupling energy call, only for relevant fixes
486 each thermostsat fix returns this via compute_scalar()
487 ecouple = cumulative energy added to reservoir by thermostatting
488 ------------------------------------------------------------------------- */
489
energy_couple()490 double Modify::energy_couple()
491 {
492 double energy = 0.0;
493 for (int i = 0; i < n_energy_couple; i++)
494 energy += fix[list_energy_couple[i]]->compute_scalar();
495 return energy;
496 }
497
498 /* ----------------------------------------------------------------------
499 global energy call, only for relevant fixes
500 they return energy via compute_scalar()
501 called by compute pe
502 ------------------------------------------------------------------------- */
503
energy_global()504 double Modify::energy_global()
505 {
506 double energy = 0.0;
507 for (int i = 0; i < n_energy_global; i++)
508 energy += fix[list_energy_global[i]]->compute_scalar();
509 return energy;
510 }
511
512 /* ----------------------------------------------------------------------
513 peratom energy call, only for relevant fixes
514 called by compute pe/atom
515 ------------------------------------------------------------------------- */
516
energy_atom(int nlocal,double * energy)517 void Modify::energy_atom(int nlocal, double *energy)
518 {
519 int i,j;
520 double *eatom;
521
522 for (i = 0; i < n_energy_atom; i++) {
523 eatom = fix[list_energy_atom[i]]->eatom;
524 if (!eatom) continue;
525 for (j = 0; j < nlocal; j++) energy[j] += eatom[j];
526 }
527 }
528
529 /* ----------------------------------------------------------------------
530 post_run call
531 ------------------------------------------------------------------------- */
532
post_run()533 void Modify::post_run()
534 {
535 for (int i = 0; i < nfix; i++) fix[i]->post_run();
536
537 // must reset this to its default value, since computes may be added
538 // or removed between runs and with this change we will redirect any
539 // calls to addstep_compute() to addstep_compute_all() instead.
540 n_timeflag = -1;
541 }
542
543 /* ----------------------------------------------------------------------
544 create_attribute call
545 invoked when an atom is added to system during a run
546 necessary so that fixes and computes that store per-atom
547 state can initialize that state for the new atom N
548 computes can store per-atom state via a fix like fix STORE
549 compute has the create_attribute flag, not fix STORE
550 ------------------------------------------------------------------------- */
551
create_attribute(int n)552 void Modify::create_attribute(int n)
553 {
554 for (int i = 0; i < nfix; i++)
555 if (fix[i]->create_attribute) fix[i]->set_arrays(n);
556 for (int i = 0; i < ncompute; i++)
557 if (compute[i]->create_attribute) compute[i]->set_arrays(n);
558 input->variable->set_arrays(n);
559 }
560
561 /* ----------------------------------------------------------------------
562 setup rRESPA pre_force call, only for relevant fixes
563 ------------------------------------------------------------------------- */
564
setup_pre_force_respa(int vflag,int ilevel)565 void Modify::setup_pre_force_respa(int vflag, int ilevel)
566 {
567 for (int i = 0; i < n_pre_force_respa; i++)
568 fix[list_pre_force_respa[i]]->setup_pre_force_respa(vflag,ilevel);
569 }
570
571 /* ----------------------------------------------------------------------
572 1st half of rRESPA integrate call, only for relevant fixes
573 ------------------------------------------------------------------------- */
574
initial_integrate_respa(int vflag,int ilevel,int iloop)575 void Modify::initial_integrate_respa(int vflag, int ilevel, int iloop)
576 {
577 for (int i = 0; i < n_initial_integrate_respa; i++)
578 fix[list_initial_integrate_respa[i]]->
579 initial_integrate_respa(vflag,ilevel,iloop);
580 }
581
582 /* ----------------------------------------------------------------------
583 rRESPA post_integrate call, only for relevant fixes
584 ------------------------------------------------------------------------- */
585
post_integrate_respa(int ilevel,int iloop)586 void Modify::post_integrate_respa(int ilevel, int iloop)
587 {
588 for (int i = 0; i < n_post_integrate_respa; i++)
589 fix[list_post_integrate_respa[i]]->post_integrate_respa(ilevel,iloop);
590 }
591
592 /* ----------------------------------------------------------------------
593 rRESPA pre_force call, only for relevant fixes
594 ------------------------------------------------------------------------- */
595
pre_force_respa(int vflag,int ilevel,int iloop)596 void Modify::pre_force_respa(int vflag, int ilevel, int iloop)
597 {
598 for (int i = 0; i < n_pre_force_respa; i++)
599 fix[list_pre_force_respa[i]]->pre_force_respa(vflag,ilevel,iloop);
600 }
601
602 /* ----------------------------------------------------------------------
603 rRESPA post_force call, only for relevant fixes
604 ------------------------------------------------------------------------- */
605
post_force_respa(int vflag,int ilevel,int iloop)606 void Modify::post_force_respa(int vflag, int ilevel, int iloop)
607 {
608 for (int i = 0; i < n_post_force_respa; i++)
609 fix[list_post_force_respa[i]]->post_force_respa(vflag,ilevel,iloop);
610 }
611
612 /* ----------------------------------------------------------------------
613 2nd half of rRESPA integrate call, only for relevant fixes
614 ------------------------------------------------------------------------- */
615
final_integrate_respa(int ilevel,int iloop)616 void Modify::final_integrate_respa(int ilevel, int iloop)
617 {
618 for (int i = 0; i < n_final_integrate_respa; i++)
619 fix[list_final_integrate_respa[i]]->final_integrate_respa(ilevel,iloop);
620 }
621
622 /* ----------------------------------------------------------------------
623 minimizer pre-exchange call, only for relevant fixes
624 ------------------------------------------------------------------------- */
625
min_pre_exchange()626 void Modify::min_pre_exchange()
627 {
628 for (int i = 0; i < n_min_pre_exchange; i++)
629 fix[list_min_pre_exchange[i]]->min_pre_exchange();
630 }
631
632 /* ----------------------------------------------------------------------
633 minimizer pre-neighbor call, only for relevant fixes
634 ------------------------------------------------------------------------- */
635
min_pre_neighbor()636 void Modify::min_pre_neighbor()
637 {
638 for (int i = 0; i < n_min_pre_neighbor; i++)
639 fix[list_min_pre_neighbor[i]]->min_pre_neighbor();
640 }
641
642 /* ----------------------------------------------------------------------
643 minimizer post-neighbor call, only for relevant fixes
644 ------------------------------------------------------------------------- */
645
min_post_neighbor()646 void Modify::min_post_neighbor()
647 {
648 for (int i = 0; i < n_min_post_neighbor; i++)
649 fix[list_min_post_neighbor[i]]->min_post_neighbor();
650 }
651
652 /* ----------------------------------------------------------------------
653 minimizer pre-force call, only for relevant fixes
654 ------------------------------------------------------------------------- */
655
min_pre_force(int vflag)656 void Modify::min_pre_force(int vflag)
657 {
658 for (int i = 0; i < n_min_pre_force; i++)
659 fix[list_min_pre_force[i]]->min_pre_force(vflag);
660 }
661
662 /* ----------------------------------------------------------------------
663 minimizer pre-reverse call, only for relevant fixes
664 ------------------------------------------------------------------------- */
665
min_pre_reverse(int eflag,int vflag)666 void Modify::min_pre_reverse(int eflag, int vflag)
667 {
668 for (int i = 0; i < n_min_pre_reverse; i++)
669 fix[list_min_pre_reverse[i]]->min_pre_reverse(eflag,vflag);
670 }
671
672 /* ----------------------------------------------------------------------
673 minimizer force adjustment call, only for relevant fixes
674 ------------------------------------------------------------------------- */
675
min_post_force(int vflag)676 void Modify::min_post_force(int vflag)
677 {
678 for (int i = 0; i < n_min_post_force; i++)
679 fix[list_min_post_force[i]]->min_post_force(vflag);
680 }
681
682 /* ----------------------------------------------------------------------
683 minimizer energy/force evaluation, only for relevant fixes
684 return energy and forces on extra degrees of freedom
685 ------------------------------------------------------------------------- */
686
min_energy(double * fextra)687 double Modify::min_energy(double *fextra)
688 {
689 int ifix,index;
690
691 index = 0;
692 double eng = 0.0;
693 for (int i = 0; i < n_min_energy; i++) {
694 ifix = list_min_energy[i];
695 eng += fix[ifix]->min_energy(&fextra[index]);
696 index += fix[ifix]->min_dof();
697 }
698 return eng;
699 }
700
701 /* ----------------------------------------------------------------------
702 store current state of extra minimizer dof, only for relevant fixes
703 ------------------------------------------------------------------------- */
704
min_store()705 void Modify::min_store()
706 {
707 for (int i = 0; i < n_min_energy; i++)
708 fix[list_min_energy[i]]->min_store();
709 }
710
711 /* ----------------------------------------------------------------------
712 manage state of extra minimizer dof on a stack, only for relevant fixes
713 ------------------------------------------------------------------------- */
714
min_clearstore()715 void Modify::min_clearstore()
716 {
717 for (int i = 0; i < n_min_energy; i++)
718 fix[list_min_energy[i]]->min_clearstore();
719 }
720
min_pushstore()721 void Modify::min_pushstore()
722 {
723 for (int i = 0; i < n_min_energy; i++)
724 fix[list_min_energy[i]]->min_pushstore();
725 }
726
min_popstore()727 void Modify::min_popstore()
728 {
729 for (int i = 0; i < n_min_energy; i++)
730 fix[list_min_energy[i]]->min_popstore();
731 }
732
733 /* ----------------------------------------------------------------------
734 displace extra minimizer dof along vector hextra, only for relevant fixes
735 ------------------------------------------------------------------------- */
736
min_step(double alpha,double * hextra)737 void Modify::min_step(double alpha, double *hextra)
738 {
739 int ifix,index;
740
741 index = 0;
742 for (int i = 0; i < n_min_energy; i++) {
743 ifix = list_min_energy[i];
744 fix[ifix]->min_step(alpha,&hextra[index]);
745 index += fix[ifix]->min_dof();
746 }
747 }
748
749 /* ----------------------------------------------------------------------
750 compute max allowed step size along vector hextra, only for relevant fixes
751 ------------------------------------------------------------------------- */
752
max_alpha(double * hextra)753 double Modify::max_alpha(double *hextra)
754 {
755 int ifix,index;
756
757 double alpha = BIG;
758 index = 0;
759 for (int i = 0; i < n_min_energy; i++) {
760 ifix = list_min_energy[i];
761 double alpha_one = fix[ifix]->max_alpha(&hextra[index]);
762 alpha = MIN(alpha,alpha_one);
763 index += fix[ifix]->min_dof();
764 }
765 return alpha;
766 }
767
768 /* ----------------------------------------------------------------------
769 extract extra minimizer dof, only for relevant fixes
770 ------------------------------------------------------------------------- */
771
min_dof()772 int Modify::min_dof()
773 {
774 int ndof = 0;
775 for (int i = 0; i < n_min_energy; i++)
776 ndof += fix[list_min_energy[i]]->min_dof();
777 return ndof;
778 }
779
780 /* ----------------------------------------------------------------------
781 reset minimizer reference state of fix, only for relevant fixes
782 ------------------------------------------------------------------------- */
783
min_reset_ref()784 int Modify::min_reset_ref()
785 {
786 int itmp,itmpall;
787 itmpall = 0;
788 for (int i = 0; i < n_min_energy; i++) {
789 itmp = fix[list_min_energy[i]]->min_reset_ref();
790 if (itmp) itmpall = 1;
791 }
792 return itmpall;
793 }
794
795 /* ----------------------------------------------------------------------
796 add a new fix or replace one with same ID
797 ------------------------------------------------------------------------- */
798
add_fix(int narg,char ** arg,int trysuffix)799 Fix *Modify::add_fix(int narg, char **arg, int trysuffix)
800 {
801 if (narg < 3) error->all(FLERR,"Illegal fix command");
802
803 // cannot define fix before box exists unless style is in exception list
804 // don't like this way of checking for exceptions by adding fixes to list,
805 // but can't think of better way
806 // too late if instantiate fix, then check flag set in fix constructor,
807 // since some fixes access domain settings in their constructor
808 // nullptr must be last entry in this list
809
810 const char *exceptions[] =
811 {"GPU", "OMP", "INTEL", "property/atom", "cmap", "cmap3", "rx",
812 "deprecated", "STORE/KIM", nullptr};
813
814 if (domain->box_exist == 0) {
815 int m;
816 for (m = 0; exceptions[m] != nullptr; m++)
817 if (strcmp(arg[2],exceptions[m]) == 0) break;
818 if (exceptions[m] == nullptr)
819 error->all(FLERR,"Fix command before simulation box is defined");
820 }
821
822 // check group ID
823
824 int igroup = group->find(arg[1]);
825 if (igroup == -1) error->all(FLERR,"Could not find fix group ID");
826
827 // if fix ID exists:
828 // set newflag = 0 so create new fix in same location in fix list
829 // error if new style does not match old style
830 // since can't replace it (all when-to-invoke ptrs would be invalid)
831 // warn if new group != old group
832 // delete old fix, but do not call update_callback(),
833 // since will replace this fix and thus other fix locs will not change
834 // set ptr to a null pointer in case new fix scans list of fixes,
835 // e.g. scan will occur in add_callback() if called by new fix
836 // if fix ID does not exist:
837 // set newflag = 1 so create new fix
838 // extend fix and fmask lists as necessary
839
840 int ifix,newflag;
841 for (ifix = 0; ifix < nfix; ifix++)
842 if (strcmp(arg[0],fix[ifix]->id) == 0) break;
843
844 if (ifix < nfix) {
845 newflag = 0;
846
847 int match = 0;
848 if (strcmp(arg[2],fix[ifix]->style) == 0) match = 1;
849 if (!match && trysuffix && lmp->suffix_enable) {
850 if (lmp->suffix) {
851 std::string estyle = arg[2] + std::string("/") + lmp->suffix;
852 if (estyle == fix[ifix]->style) match = 1;
853 }
854 if (lmp->suffix2) {
855 std::string estyle = arg[2] + std::string("/") + lmp->suffix2;
856 if (estyle == fix[ifix]->style) match = 1;
857 }
858 }
859 if (!match)
860 error->all(FLERR,"Replacing a fix, but new style != old style");
861
862 if (fix[ifix]->igroup != igroup && comm->me == 0)
863 error->warning(FLERR,"Replacing a fix, but new group != old group");
864 delete fix[ifix];
865 fix[ifix] = nullptr;
866
867 } else {
868 newflag = 1;
869 if (nfix == maxfix) {
870 maxfix += DELTA;
871 fix = (Fix **) memory->srealloc(fix,maxfix*sizeof(Fix *),"modify:fix");
872 memory->grow(fmask,maxfix,"modify:fmask");
873 }
874 }
875
876 // create the Fix
877 // try first with suffix appended
878
879 fix[ifix] = nullptr;
880
881 if (trysuffix && lmp->suffix_enable) {
882 if (lmp->suffix) {
883 std::string estyle = arg[2] + std::string("/") + lmp->suffix;
884 if (fix_map->find(estyle) != fix_map->end()) {
885 FixCreator &fix_creator = (*fix_map)[estyle];
886 fix[ifix] = fix_creator(lmp,narg,arg);
887 delete[] fix[ifix]->style;
888 fix[ifix]->style = utils::strdup(estyle);
889 }
890 }
891 if (fix[ifix] == nullptr && lmp->suffix2) {
892 std::string estyle = arg[2] + std::string("/") + lmp->suffix2;
893 if (fix_map->find(estyle) != fix_map->end()) {
894 FixCreator &fix_creator = (*fix_map)[estyle];
895 fix[ifix] = fix_creator(lmp,narg,arg);
896 delete[] fix[ifix]->style;
897 fix[ifix]->style = utils::strdup(estyle);
898 }
899 }
900 }
901
902 if (fix[ifix] == nullptr && fix_map->find(arg[2]) != fix_map->end()) {
903 FixCreator &fix_creator = (*fix_map)[arg[2]];
904 fix[ifix] = fix_creator(lmp,narg,arg);
905 }
906
907 if (fix[ifix] == nullptr)
908 error->all(FLERR,utils::check_packages_for_style("fix",arg[2],lmp));
909
910 // increment nfix (if new)
911
912 if (newflag) nfix++;
913
914 // post_constructor() can call virtual methods in parent or child
915 // which would otherwise not yet be visible in child class
916 // post_constructor() allows new fix to create other fixes
917 // nfix increment must come first so recursive call to add_fix within
918 // post_constructor() will see updated nfix
919
920 fix[ifix]->post_constructor();
921
922 // check if Fix is in restart_global list
923 // if yes, pass state info to the Fix so it can reset itself
924
925 for (int i = 0; i < nfix_restart_global; i++)
926 if (strcmp(id_restart_global[i],fix[ifix]->id) == 0 &&
927 strcmp(style_restart_global[i],fix[ifix]->style) == 0) {
928 fix[ifix]->restart(state_restart_global[i]);
929 used_restart_global[i] = 1;
930 fix[ifix]->restart_reset = 1;
931 if (comm->me == 0)
932 utils::logmesg(lmp,"Resetting global fix info from restart file:\n"
933 " fix style: {}, fix ID: {}\n",
934 fix[ifix]->style,fix[ifix]->id);
935 }
936
937 // check if Fix is in restart_peratom list
938 // if yes, loop over atoms so they can extract info from atom->extra array
939
940 for (int i = 0; i < nfix_restart_peratom; i++)
941 if (strcmp(id_restart_peratom[i],fix[ifix]->id) == 0 &&
942 strcmp(style_restart_peratom[i],fix[ifix]->style) == 0) {
943 used_restart_peratom[i] = 1;
944 for (int j = 0; j < atom->nlocal; j++)
945 fix[ifix]->unpack_restart(j,index_restart_peratom[i]);
946 fix[ifix]->restart_reset = 1;
947 if (comm->me == 0)
948 utils::logmesg(lmp,"Resetting peratom fix info from restart file:\n"
949 " fix style: {}, fix ID: {}\n",
950 fix[ifix]->style,fix[ifix]->id);
951 }
952
953 // set fix mask values
954
955 fmask[ifix] = fix[ifix]->setmask();
956
957 // return pointer to fix
958
959 return fix[ifix];
960 }
961
962 /* ----------------------------------------------------------------------
963 convenience function to allow adding a fix from a single string
964 ------------------------------------------------------------------------- */
965
add_fix(const std::string & fixcmd,int trysuffix)966 Fix *Modify::add_fix(const std::string &fixcmd, int trysuffix)
967 {
968 auto args = utils::split_words(fixcmd);
969 std::vector<char *> newarg(args.size());
970 int i = 0;
971 for (const auto &arg : args) {
972 newarg[i++] = (char *)arg.c_str();
973 }
974 return add_fix(args.size(),newarg.data(),trysuffix);
975 }
976
977 /* ----------------------------------------------------------------------
978 replace replaceID fix with a new fix
979 this is used by callers to preserve ordering of fixes
980 e.g. create replaceID as a FixDummy instance early in the input script
981 replace it later with the desired Fix instance
982 ------------------------------------------------------------------------- */
983
replace_fix(const char * replaceID,int narg,char ** arg,int trysuffix)984 Fix *Modify::replace_fix(const char *replaceID, int narg, char **arg, int trysuffix)
985 {
986 int ifix = find_fix(replaceID);
987 if (ifix < 0) error->all(FLERR,"Modify replace_fix ID {} could not be found", replaceID);
988
989 // change ID, igroup, style of fix being replaced to match new fix
990 // requires some error checking on arguments for new fix
991
992 if (narg < 3) error->all(FLERR,"Illegal replace_fix invocation");
993 int jfix = find_fix(arg[0]);
994 if (jfix >= 0) error->all(FLERR,"Replace_fix ID is already in use");
995
996 delete [] fix[ifix]->id;
997 fix[ifix]->id = utils::strdup(arg[0]);
998
999 int jgroup = group->find(arg[1]);
1000 if (jgroup == -1) error->all(FLERR,"Could not find replace_fix group ID");
1001 fix[ifix]->igroup = jgroup;
1002
1003 delete [] fix[ifix]->style;
1004 fix[ifix]->style = utils::strdup(arg[2]);
1005
1006 // invoke add_fix
1007 // it will find and overwrite the replaceID fix
1008
1009 return add_fix(narg,arg,trysuffix);
1010 }
1011
1012 /* ----------------------------------------------------------------------
1013 convenience function to allow replacing a fix from a single string
1014 ------------------------------------------------------------------------- */
1015
replace_fix(const std::string & oldfix,const std::string & fixcmd,int trysuffix)1016 Fix *Modify::replace_fix(const std::string &oldfix, const std::string &fixcmd, int trysuffix)
1017 {
1018 auto args = utils::split_words(fixcmd);
1019 std::vector<char *> newarg(args.size());
1020 int i = 0;
1021 for (const auto &arg : args) {
1022 newarg[i++] = (char *)arg.c_str();
1023 }
1024 return replace_fix(oldfix.c_str(),args.size(),newarg.data(),trysuffix);
1025 }
1026
1027 /* ----------------------------------------------------------------------
1028 one instance per fix in style_fix.h
1029 ------------------------------------------------------------------------- */
1030
1031 template <typename T>
fix_creator(LAMMPS * lmp,int narg,char ** arg)1032 Fix *Modify::fix_creator(LAMMPS *lmp, int narg, char **arg)
1033 {
1034 return new T(lmp,narg,arg);
1035 }
1036
1037 /* ----------------------------------------------------------------------
1038 modify a Fix's parameters
1039 ------------------------------------------------------------------------- */
1040
modify_fix(int narg,char ** arg)1041 void Modify::modify_fix(int narg, char **arg)
1042 {
1043 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
1044
1045 // lookup Fix ID
1046
1047 int ifix;
1048 for (ifix = 0; ifix < nfix; ifix++)
1049 if (strcmp(arg[0],fix[ifix]->id) == 0) break;
1050 if (ifix == nfix) error->all(FLERR,"Could not find fix_modify ID");
1051
1052 fix[ifix]->modify_params(narg-1,&arg[1]);
1053 }
1054
1055 /* ----------------------------------------------------------------------
1056 delete a Fix from list of Fixes
1057 Atom class must update indices in its list of callbacks to fixes
1058 ------------------------------------------------------------------------- */
1059
delete_fix(const std::string & id)1060 void Modify::delete_fix(const std::string &id)
1061 {
1062 int ifix = find_fix(id);
1063 if (ifix < 0) error->all(FLERR,"Could not find fix ID to delete");
1064 delete_fix(ifix);
1065 }
1066
delete_fix(int ifix)1067 void Modify::delete_fix(int ifix)
1068 {
1069 if ((ifix < 0) || (ifix >= nfix)) return;
1070
1071 // delete instance and move other Fixes and fmask down in list one slot
1072
1073 delete fix[ifix];
1074 atom->update_callback(ifix);
1075
1076 for (int i = ifix+1; i < nfix; i++) fix[i-1] = fix[i];
1077 for (int i = ifix+1; i < nfix; i++) fmask[i-1] = fmask[i];
1078 nfix--;
1079 }
1080
1081 /* ----------------------------------------------------------------------
1082 find a fix by ID
1083 return index of fix or -1 if not found
1084 ------------------------------------------------------------------------- */
1085
find_fix(const std::string & id)1086 int Modify::find_fix(const std::string &id)
1087 {
1088 if (id.empty()) return -1;
1089 for (int ifix = 0; ifix < nfix; ifix++)
1090 if (id == fix[ifix]->id) return ifix;
1091 return -1;
1092 }
1093
1094 /* ----------------------------------------------------------------------
1095 find a fix by style
1096 return index of fix or -1 if not found
1097 ------------------------------------------------------------------------- */
1098
find_fix_by_style(const char * style)1099 int Modify::find_fix_by_style(const char *style)
1100 {
1101 for (int ifix = 0; ifix < nfix; ifix++)
1102 if (utils::strmatch(fix[ifix]->style,style)) return ifix;
1103 return -1;
1104 }
1105
1106 /* ----------------------------------------------------------------------
1107 check for fix associated with package name in compiled list
1108 return 1 if found else 0
1109 used to determine whether LAMMPS was built with
1110 GPU, INTEL, OPENMP packages, which have their own fixes
1111 ------------------------------------------------------------------------- */
1112
check_package(const char * package_fix_name)1113 int Modify::check_package(const char *package_fix_name)
1114 {
1115 if (fix_map->find(package_fix_name) == fix_map->end()) return 0;
1116 return 1;
1117 }
1118
1119 /* ----------------------------------------------------------------------
1120 check if the group indicated by groupbit overlaps with any
1121 currently existing rigid fixes. return 1 in this case otherwise 0
1122 ------------------------------------------------------------------------- */
1123
check_rigid_group_overlap(int groupbit)1124 int Modify::check_rigid_group_overlap(int groupbit)
1125 {
1126 const int * const mask = atom->mask;
1127 const int nlocal = atom->nlocal;
1128 int dim;
1129
1130 int n = 0;
1131 for (int ifix = 0; ifix < nfix; ifix++) {
1132 if (utils::strmatch(fix[ifix]->style,"^rigid")) {
1133 const int * const body = (const int *)fix[ifix]->extract("body",dim);
1134 if ((body == nullptr) || (dim != 1)) break;
1135
1136 for (int i=0; (i < nlocal) && (n == 0); ++i)
1137 if ((mask[i] & groupbit) && (body[i] >= 0)) ++n;
1138 }
1139 }
1140
1141 int n_all = 0;
1142 MPI_Allreduce(&n,&n_all,1,MPI_INT,MPI_SUM,world);
1143
1144 if (n_all > 0) return 1;
1145 return 0;
1146 }
1147
1148 /* ----------------------------------------------------------------------
1149 check if the atoms in the group indicated by groupbit _and_ region
1150 indicated by regionid overlap with any currently existing rigid fixes.
1151 return 1 in this case, otherwise 0
1152 ------------------------------------------------------------------------- */
1153
check_rigid_region_overlap(int groupbit,Region * reg)1154 int Modify::check_rigid_region_overlap(int groupbit, Region *reg)
1155 {
1156 const int * const mask = atom->mask;
1157 const double * const * const x = atom->x;
1158 const int nlocal = atom->nlocal;
1159 int dim;
1160
1161 int n = 0;
1162 reg->prematch();
1163 for (int ifix = 0; ifix < nfix; ifix++) {
1164 if (strncmp("rigid",fix[ifix]->style,5) == 0) {
1165 const int * const body = (const int *)fix[ifix]->extract("body",dim);
1166 if ((body == nullptr) || (dim != 1)) break;
1167
1168 for (int i=0; (i < nlocal) && (n == 0); ++i)
1169 if ((mask[i] & groupbit) && (body[i] >= 0)
1170 && reg->match(x[i][0],x[i][1],x[i][2])) ++n;
1171 }
1172 }
1173
1174 int n_all = 0;
1175 MPI_Allreduce(&n,&n_all,1,MPI_INT,MPI_SUM,world);
1176
1177 if (n_all > 0) return 1;
1178 return 0;
1179 }
1180
1181 /* ----------------------------------------------------------------------
1182 check if the atoms in the selection list (length atom->nlocal,
1183 content: 1 if atom is contained, 0 if not) overlap with currently
1184 existing rigid fixes. return 1 in this case otherwise 0
1185 ------------------------------------------------------------------------- */
1186
check_rigid_list_overlap(int * select)1187 int Modify::check_rigid_list_overlap(int *select)
1188 {
1189 const int nlocal = atom->nlocal;
1190 int dim;
1191
1192 int n = 0;
1193 for (int ifix = 0; ifix < nfix; ifix++) {
1194 if (utils::strmatch(fix[ifix]->style,"^rigid")) {
1195 const int * const body = (const int *)fix[ifix]->extract("body",dim);
1196 if ((body == nullptr) || (dim != 1)) break;
1197
1198 for (int i=0; (i < nlocal) && (n == 0); ++i)
1199 if ((body[i] >= 0) && select[i]) ++n;
1200 }
1201 }
1202
1203 int n_all = 0;
1204 MPI_Allreduce(&n,&n_all,1,MPI_INT,MPI_SUM,world);
1205
1206 if (n_all > 0) return 1;
1207 return 0;
1208 }
1209
1210 /* ----------------------------------------------------------------------
1211 add a new compute
1212 ------------------------------------------------------------------------- */
1213
add_compute(int narg,char ** arg,int trysuffix)1214 Compute *Modify::add_compute(int narg, char **arg, int trysuffix)
1215 {
1216 if (narg < 3) error->all(FLERR,"Illegal compute command");
1217
1218 // error check
1219
1220 for (int icompute = 0; icompute < ncompute; icompute++)
1221 if (strcmp(arg[0],compute[icompute]->id) == 0)
1222 error->all(FLERR,"Reuse of compute ID '{}'",arg[0]);
1223
1224 // extend Compute list if necessary
1225
1226 if (ncompute == maxcompute) {
1227 maxcompute += DELTA;
1228 compute = (Compute **)
1229 memory->srealloc(compute,maxcompute*sizeof(Compute *),"modify:compute");
1230 }
1231
1232 // create the Compute
1233 // try first with suffix appended
1234
1235 compute[ncompute] = nullptr;
1236
1237 if (trysuffix && lmp->suffix_enable) {
1238 if (lmp->suffix) {
1239 std::string estyle = arg[2] + std::string("/") + lmp->suffix;
1240 if (compute_map->find(estyle) != compute_map->end()) {
1241 ComputeCreator &compute_creator = (*compute_map)[estyle];
1242 compute[ncompute] = compute_creator(lmp,narg,arg);
1243 delete[] compute[ncompute]->style;
1244 compute[ncompute]->style = utils::strdup(estyle);
1245 }
1246 }
1247 if (compute[ncompute] == nullptr && lmp->suffix2) {
1248 std::string estyle = arg[2] + std::string("/") + lmp->suffix2;
1249 if (compute_map->find(estyle) != compute_map->end()) {
1250 ComputeCreator &compute_creator = (*compute_map)[estyle];
1251 compute[ncompute] = compute_creator(lmp,narg,arg);
1252 delete[] compute[ncompute]->style;
1253 compute[ncompute]->style = utils::strdup(estyle);
1254 }
1255 }
1256 }
1257
1258 if (compute[ncompute] == nullptr &&
1259 compute_map->find(arg[2]) != compute_map->end()) {
1260 ComputeCreator &compute_creator = (*compute_map)[arg[2]];
1261 compute[ncompute] = compute_creator(lmp,narg,arg);
1262 }
1263
1264 if (compute[ncompute] == nullptr)
1265 error->all(FLERR,utils::check_packages_for_style("compute",arg[2],lmp));
1266
1267 return compute[ncompute++];
1268 }
1269
1270 /* ----------------------------------------------------------------------
1271 convenience function to allow adding a compute from a single string
1272 ------------------------------------------------------------------------- */
1273
add_compute(const std::string & computecmd,int trysuffix)1274 Compute *Modify::add_compute(const std::string &computecmd, int trysuffix)
1275 {
1276 auto args = utils::split_words(computecmd);
1277 std::vector<char *>newarg(args.size());
1278 int i=0;
1279 for (const auto &arg : args) {
1280 newarg[i++] = (char *)arg.c_str();
1281 }
1282 return add_compute(args.size(),newarg.data(),trysuffix);
1283 }
1284
1285
1286 /* ----------------------------------------------------------------------
1287 one instance per compute in style_compute.h
1288 ------------------------------------------------------------------------- */
1289
1290 template <typename T>
compute_creator(LAMMPS * lmp,int narg,char ** arg)1291 Compute *Modify::compute_creator(LAMMPS *lmp, int narg, char **arg)
1292 {
1293 return new T(lmp,narg,arg);
1294 }
1295
1296 /* ----------------------------------------------------------------------
1297 modify a Compute's parameters
1298 ------------------------------------------------------------------------- */
1299
modify_compute(int narg,char ** arg)1300 void Modify::modify_compute(int narg, char **arg)
1301 {
1302 if (narg < 2) error->all(FLERR,"Illegal compute_modify command");
1303
1304 // lookup Compute ID
1305
1306 int icompute;
1307 for (icompute = 0; icompute < ncompute; icompute++)
1308 if (strcmp(arg[0],compute[icompute]->id) == 0) break;
1309 if (icompute == ncompute)
1310 error->all(FLERR,"Could not find compute_modify ID");
1311
1312 compute[icompute]->modify_params(narg-1,&arg[1]);
1313 }
1314
1315 /* ----------------------------------------------------------------------
1316 delete a Compute from list of Computes
1317 ------------------------------------------------------------------------- */
1318
delete_compute(const std::string & id)1319 void Modify::delete_compute(const std::string &id)
1320 {
1321 int icompute = find_compute(id);
1322 if (icompute < 0) error->all(FLERR,"Could not find compute ID to delete");
1323 delete_compute(icompute);
1324 }
1325
delete_compute(int icompute)1326 void Modify::delete_compute(int icompute)
1327 {
1328 if ((icompute < 0) || (icompute >= ncompute)) return;
1329
1330 // delete and move other Computes down in list one slot
1331
1332 delete compute[icompute];
1333 for (int i = icompute+1; i < ncompute; i++) compute[i-1] = compute[i];
1334 ncompute--;
1335 }
1336
1337 /* ----------------------------------------------------------------------
1338 find a compute by ID
1339 return index of compute or -1 if not found
1340 ------------------------------------------------------------------------- */
1341
find_compute(const std::string & id)1342 int Modify::find_compute(const std::string &id)
1343 {
1344 if (id.empty()) return -1;
1345 for (int icompute = 0; icompute < ncompute; icompute++)
1346 if (id == compute[icompute]->id) return icompute;
1347 return -1;
1348 }
1349
1350 /* ----------------------------------------------------------------------
1351 find a compute by style
1352 return index of compute or -1 if not found
1353 ------------------------------------------------------------------------- */
1354
find_compute_by_style(const char * style)1355 int Modify::find_compute_by_style(const char *style)
1356 {
1357 for (int icompute = 0; icompute < ncompute; icompute++)
1358 if (utils::strmatch(compute[icompute]->style,style)) return icompute;
1359 return -1;
1360 }
1361
1362 /* ----------------------------------------------------------------------
1363 clear invoked flag of all computes
1364 called everywhere that computes are used, before computes are invoked
1365 invoked flag used to avoid re-invoking same compute multiple times
1366 and to flag computes that store invocation times as having been invoked
1367 ------------------------------------------------------------------------- */
1368
clearstep_compute()1369 void Modify::clearstep_compute()
1370 {
1371 for (int icompute = 0; icompute < ncompute; icompute++)
1372 compute[icompute]->invoked_flag = Compute::INVOKED_NONE;
1373 }
1374
1375 /* ----------------------------------------------------------------------
1376 loop over computes that store invocation times
1377 if its invoked flag set on this timestep, schedule next invocation
1378 called everywhere that computes are used, after computes are invoked
1379 ------------------------------------------------------------------------- */
1380
addstep_compute(bigint newstep)1381 void Modify::addstep_compute(bigint newstep)
1382 {
1383 // If we are called before the first run init, n_timeflag is not yet
1384 // initialized, thus defer to addstep_compute_all() instead
1385
1386 if (n_timeflag < 0) {
1387 addstep_compute_all(newstep);
1388 return;
1389 }
1390
1391 for (int icompute = 0; icompute < n_timeflag; icompute++)
1392 if (compute[list_timeflag[icompute]]->invoked_flag)
1393 compute[list_timeflag[icompute]]->addstep(newstep);
1394 }
1395
1396 /* ----------------------------------------------------------------------
1397 loop over all computes
1398 schedule next invocation for those that store invocation times
1399 called when not sure what computes will be needed on newstep
1400 do not loop only over n_timeflag, since may not be set yet
1401 ------------------------------------------------------------------------- */
1402
addstep_compute_all(bigint newstep)1403 void Modify::addstep_compute_all(bigint newstep)
1404 {
1405 for (int icompute = 0; icompute < ncompute; icompute++)
1406 if (compute[icompute]->timeflag) compute[icompute]->addstep(newstep);
1407 }
1408
1409 /* ----------------------------------------------------------------------
1410 write to restart file for all Fixes with restart info
1411 (1) fixes that have global state
1412 (2) fixes that store per-atom quantities
1413 ------------------------------------------------------------------------- */
1414
write_restart(FILE * fp)1415 void Modify::write_restart(FILE *fp)
1416 {
1417 int me = comm->me;
1418
1419 int count = 0;
1420 for (int i = 0; i < nfix; i++)
1421 if (fix[i]->restart_global) count++;
1422
1423 if (me == 0) fwrite(&count,sizeof(int),1,fp);
1424
1425 int n;
1426 for (int i = 0; i < nfix; i++)
1427 if (fix[i]->restart_global) {
1428 if (me == 0) {
1429 n = strlen(fix[i]->id) + 1;
1430 fwrite(&n,sizeof(int),1,fp);
1431 fwrite(fix[i]->id,sizeof(char),n,fp);
1432 n = strlen(fix[i]->style) + 1;
1433 fwrite(&n,sizeof(int),1,fp);
1434 fwrite(fix[i]->style,sizeof(char),n,fp);
1435 }
1436 fix[i]->write_restart(fp);
1437 }
1438
1439 count = 0;
1440 for (int i = 0; i < nfix; i++)
1441 if (fix[i]->restart_peratom) count++;
1442
1443 if (me == 0) fwrite(&count,sizeof(int),1,fp);
1444
1445 for (int i = 0; i < nfix; i++)
1446 if (fix[i]->restart_peratom) {
1447 int maxsize_restart = fix[i]->maxsize_restart();
1448 if (me == 0) {
1449 n = strlen(fix[i]->id) + 1;
1450 fwrite(&n,sizeof(int),1,fp);
1451 fwrite(fix[i]->id,sizeof(char),n,fp);
1452 n = strlen(fix[i]->style) + 1;
1453 fwrite(&n,sizeof(int),1,fp);
1454 fwrite(fix[i]->style,sizeof(char),n,fp);
1455 fwrite(&maxsize_restart,sizeof(int),1,fp);
1456 }
1457 }
1458 }
1459
1460 /* ----------------------------------------------------------------------
1461 read in restart file data on all previously defined Fixes with restart info
1462 (1) fixes that have global state
1463 (2) fixes that store per-atom quantities
1464 return maxsize of extra info that will be stored with any atom
1465 ------------------------------------------------------------------------- */
1466
read_restart(FILE * fp)1467 int Modify::read_restart(FILE *fp)
1468 {
1469 // nfix_restart_global = # of restart entries with global state info
1470
1471 int me = comm->me;
1472 if (me == 0) utils::sfread(FLERR,&nfix_restart_global,sizeof(int),1,fp,nullptr,error);
1473 MPI_Bcast(&nfix_restart_global,1,MPI_INT,0,world);
1474
1475 // allocate space for each entry
1476
1477 if (nfix_restart_global) {
1478 id_restart_global = new char*[nfix_restart_global];
1479 style_restart_global = new char*[nfix_restart_global];
1480 state_restart_global = new char*[nfix_restart_global];
1481 used_restart_global = new int[nfix_restart_global];
1482 }
1483
1484 // read each entry and Bcast to all procs
1485 // each entry has id string, style string, chunk of state data
1486
1487 int n;
1488 for (int i = 0; i < nfix_restart_global; i++) {
1489 if (me == 0) utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error);
1490 MPI_Bcast(&n,1,MPI_INT,0,world);
1491 id_restart_global[i] = new char[n];
1492 if (me == 0) utils::sfread(FLERR,id_restart_global[i],sizeof(char),n,fp,nullptr,error);
1493 MPI_Bcast(id_restart_global[i],n,MPI_CHAR,0,world);
1494
1495 if (me == 0) utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error);
1496 MPI_Bcast(&n,1,MPI_INT,0,world);
1497 style_restart_global[i] = new char[n];
1498 if (me == 0) utils::sfread(FLERR,style_restart_global[i],sizeof(char),n,fp,nullptr,error);
1499 MPI_Bcast(style_restart_global[i],n,MPI_CHAR,0,world);
1500
1501 if (me == 0) utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error);
1502 MPI_Bcast(&n,1,MPI_INT,0,world);
1503 state_restart_global[i] = new char[n];
1504 if (me == 0) utils::sfread(FLERR,state_restart_global[i],sizeof(char),n,fp,nullptr,error);
1505 MPI_Bcast(state_restart_global[i],n,MPI_CHAR,0,world);
1506
1507 used_restart_global[i] = 0;
1508 }
1509
1510 // nfix_restart_peratom = # of restart entries with peratom info
1511
1512 int maxsize = 0;
1513
1514 if (me == 0) utils::sfread(FLERR,&nfix_restart_peratom,sizeof(int),1,fp,nullptr,error);
1515 MPI_Bcast(&nfix_restart_peratom,1,MPI_INT,0,world);
1516
1517 // allocate space for each entry
1518
1519 if (nfix_restart_peratom) {
1520 id_restart_peratom = new char*[nfix_restart_peratom];
1521 style_restart_peratom = new char*[nfix_restart_peratom];
1522 index_restart_peratom = new int[nfix_restart_peratom];
1523 used_restart_peratom = new int[nfix_restart_peratom];
1524 }
1525
1526 // read each entry and Bcast to all procs
1527 // each entry has id string, style string, maxsize of one atom's data
1528 // set index = which set of extra data this fix represents
1529
1530 for (int i = 0; i < nfix_restart_peratom; i++) {
1531 if (me == 0) utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error);
1532 MPI_Bcast(&n,1,MPI_INT,0,world);
1533 id_restart_peratom[i] = new char[n];
1534 if (me == 0) utils::sfread(FLERR,id_restart_peratom[i],sizeof(char),n,fp,nullptr,error);
1535 MPI_Bcast(id_restart_peratom[i],n,MPI_CHAR,0,world);
1536
1537 if (me == 0) utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error);
1538 MPI_Bcast(&n,1,MPI_INT,0,world);
1539 style_restart_peratom[i] = new char[n];
1540 if (me == 0) utils::sfread(FLERR,style_restart_peratom[i],sizeof(char),n,fp,nullptr,error);
1541 MPI_Bcast(style_restart_peratom[i],n,MPI_CHAR,0,world);
1542
1543 if (me == 0) utils::sfread(FLERR,&n,sizeof(int),1,fp,nullptr,error);
1544 MPI_Bcast(&n,1,MPI_INT,0,world);
1545 maxsize += n;
1546
1547 index_restart_peratom[i] = i;
1548 used_restart_peratom[i] = 0;
1549 }
1550
1551 return maxsize;
1552 }
1553
1554 /* ----------------------------------------------------------------------
1555 delete all lists of restart file Fix info
1556 if flag set, print list of restart file info not assigned to new fixes
1557 ------------------------------------------------------------------------- */
1558
restart_deallocate(int flag)1559 void Modify::restart_deallocate(int flag)
1560 {
1561 if (nfix_restart_global) {
1562 if (flag && comm->me == 0) {
1563 int i;
1564 for (i = 0; i < nfix_restart_global; i++)
1565 if (used_restart_global[i] == 0) break;
1566 if (i == nfix_restart_global) {
1567 utils::logmesg(lmp,"All restart file global fix info was re-assigned\n");
1568 } else {
1569 utils::logmesg(lmp,"Unused restart file global fix info:\n");
1570 for (i = 0; i < nfix_restart_global; i++) {
1571 if (used_restart_global[i]) continue;
1572 utils::logmesg(lmp," fix style: {}, fix ID: {}\n",
1573 style_restart_global[i],id_restart_global[i]);
1574 }
1575 }
1576 }
1577
1578 for (int i = 0; i < nfix_restart_global; i++) {
1579 delete [] id_restart_global[i];
1580 delete [] style_restart_global[i];
1581 delete [] state_restart_global[i];
1582 }
1583 delete [] id_restart_global;
1584 delete [] style_restart_global;
1585 delete [] state_restart_global;
1586 delete [] used_restart_global;
1587 }
1588
1589 if (nfix_restart_peratom) {
1590 if (flag && comm->me == 0) {
1591 int i;
1592 for (i = 0; i < nfix_restart_peratom; i++)
1593 if (used_restart_peratom[i] == 0) break;
1594 if (i == nfix_restart_peratom) {
1595 utils::logmesg(lmp,"All restart file peratom fix info was re-assigned\n");
1596 } else {
1597 utils::logmesg(lmp,"Unused restart file peratom fix info:\n");
1598 for (i = 0; i < nfix_restart_peratom; i++) {
1599 if (used_restart_peratom[i]) continue;
1600 utils::logmesg(lmp," fix style: {}, fix ID: {}\n",
1601 style_restart_peratom[i],id_restart_peratom[i]);
1602 }
1603 }
1604 }
1605
1606 for (int i = 0; i < nfix_restart_peratom; i++) {
1607 delete [] id_restart_peratom[i];
1608 delete [] style_restart_peratom[i];
1609 }
1610 delete [] id_restart_peratom;
1611 delete [] style_restart_peratom;
1612 delete [] index_restart_peratom;
1613 delete [] used_restart_peratom;
1614 }
1615
1616 nfix_restart_global = nfix_restart_peratom = 0;
1617 }
1618
1619 /* ----------------------------------------------------------------------
1620 create list of fix indices for fixes which match mask
1621 ------------------------------------------------------------------------- */
1622
list_init(int mask,int & n,int * & list)1623 void Modify::list_init(int mask, int &n, int *&list)
1624 {
1625 delete [] list;
1626
1627 n = 0;
1628 for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++;
1629 list = new int[n];
1630
1631 n = 0;
1632 for (int i = 0; i < nfix; i++) if (fmask[i] & mask) list[n++] = i;
1633 }
1634
1635 /* ----------------------------------------------------------------------
1636 create list of fix indices for end_of_step fixes
1637 also create end_of_step_every[]
1638 ------------------------------------------------------------------------- */
1639
list_init_end_of_step(int mask,int & n,int * & list)1640 void Modify::list_init_end_of_step(int mask, int &n, int *&list)
1641 {
1642 delete [] list;
1643 delete [] end_of_step_every;
1644
1645 n = 0;
1646 for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++;
1647 list = new int[n];
1648 end_of_step_every = new int[n];
1649
1650 n = 0;
1651 for (int i = 0; i < nfix; i++)
1652 if (fmask[i] & mask) {
1653 list[n] = i;
1654 end_of_step_every[n++] = fix[i]->nevery;
1655 }
1656 }
1657
1658 /* ----------------------------------------------------------------------
1659 create list of fix indices for fixes that compute reservoir coupling energy
1660 only added to list if fix has ecouple_flag set
1661 ------------------------------------------------------------------------- */
1662
list_init_energy_couple(int & n,int * & list)1663 void Modify::list_init_energy_couple(int &n, int *&list)
1664 {
1665 delete [] list;
1666
1667 n = 0;
1668 for (int i = 0; i < nfix; i++)
1669 if (fix[i]->ecouple_flag) n++;
1670 list = new int[n];
1671
1672 n = 0;
1673 for (int i = 0; i < nfix; i++)
1674 if (fix[i]->ecouple_flag) list[n++] = i;
1675 }
1676
1677 /* ----------------------------------------------------------------------
1678 create list of fix indices for fixes that compute global energy
1679 only added to list if fix has energy_global_flag and thermo_energy set
1680 ------------------------------------------------------------------------- */
1681
list_init_energy_global(int & n,int * & list)1682 void Modify::list_init_energy_global(int &n, int *&list)
1683 {
1684 delete [] list;
1685
1686 n = 0;
1687 for (int i = 0; i < nfix; i++)
1688 if (fix[i]->energy_global_flag && fix[i]->thermo_energy) n++;
1689 list = new int[n];
1690
1691 n = 0;
1692 for (int i = 0; i < nfix; i++)
1693 if (fix[i]->energy_global_flag && fix[i]->thermo_energy) list[n++] = i;
1694 }
1695
1696 /* ----------------------------------------------------------------------
1697 create list of fix indices for fixes that compute peratom energy
1698 only added to list if fix has energy_peratom_flag and thermo_energy set
1699 ------------------------------------------------------------------------- */
1700
list_init_energy_atom(int & n,int * & list)1701 void Modify::list_init_energy_atom(int &n, int *&list)
1702 {
1703 delete [] list;
1704
1705 n = 0;
1706 for (int i = 0; i < nfix; i++)
1707 if (fix[i]->energy_peratom_flag && fix[i]->thermo_energy) n++;
1708 list = new int[n];
1709
1710 n = 0;
1711 for (int i = 0; i < nfix; i++)
1712 if (fix[i]->energy_peratom_flag && fix[i]->thermo_energy) list[n++] = i;
1713 }
1714
1715 /* ----------------------------------------------------------------------
1716 create list of compute indices for computes which store invocation times
1717 ------------------------------------------------------------------------- */
1718
list_init_compute()1719 void Modify::list_init_compute()
1720 {
1721 delete [] list_timeflag;
1722
1723 n_timeflag = 0;
1724 for (int i = 0; i < ncompute; i++)
1725 if (compute[i]->timeflag) n_timeflag++;
1726 list_timeflag = new int[n_timeflag];
1727
1728 n_timeflag = 0;
1729 for (int i = 0; i < ncompute; i++)
1730 if (compute[i]->timeflag) list_timeflag[n_timeflag++] = i;
1731 }
1732
1733 /* ----------------------------------------------------------------------
1734 return # of bytes of allocated memory from all fixes
1735 ------------------------------------------------------------------------- */
1736
memory_usage()1737 double Modify::memory_usage()
1738 {
1739 double bytes = 0;
1740 for (int i = 0; i < nfix; i++)
1741 bytes += fix[i]->memory_usage();
1742 for (int i = 0; i < ncompute; i++)
1743 bytes += compute[i]->memory_usage();
1744 return bytes;
1745 }
1746