1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 This file is from LAMMPS, but has been modified. Copyright for
36 modification:
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40
41 Copyright of original file:
42 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
43 http://lammps.sandia.gov, Sandia National Laboratories
44 Steve Plimpton, sjplimp@sandia.gov
45
46 Copyright (2003) Sandia Corporation. Under the terms of Contract
47 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
48 certain rights in this software. This software is distributed under
49 the GNU General Public License.
50 ------------------------------------------------------------------------- */
51
52 /* ----------------------------------------------------------------------
53 Contributing authors:
54 Richard Berger (JKU Linz)
55 ------------------------------------------------------------------------- */
56
57 #include <stdio.h>
58 #include <string.h>
59 #include "modify.h"
60 #include "style_compute.h"
61 #include "style_fix.h"
62 #include "atom.h"
63 #include "comm.h"
64 #include "fix.h"
65 #include "compute.h"
66 #include "group.h"
67 #include "update.h"
68 #include "domain.h"
69 #include "memory.h"
70 #include "error.h"
71 #include "force.h"
72 #include <map>
73 #include <string>
74
75 using namespace LAMMPS_NS;
76 using namespace FixConst;
77
78 #define DELTA 4
79
80 #define BIG 1.0e20
81 #define NEXCEPT 4 // change when add to exceptions in add_fix()
82
83 /* ---------------------------------------------------------------------- */
84
Modify(LAMMPS * lmp)85 Modify::Modify(LAMMPS *lmp) :
86 Pointers(lmp),
87 n_timeflag(0)
88 {
89 nfix = maxfix = 0;
90 n_pre_initial_integrate = n_initial_integrate = n_post_integrate = 0;
91 n_pre_exchange = n_pre_neighbor = 0;
92 n_pre_force = n_post_force = 0;
93 n_iterate_implicitly = n_pre_final_integrate = 0;
94 n_final_integrate = n_end_of_step = n_thermo_energy = 0;
95 n_initial_integrate_respa = n_post_integrate_respa = 0;
96 n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0;
97 n_min_pre_exchange = n_min_pre_force = n_min_post_force = n_min_energy = 0;
98 n_min_pre_neighbor = 0;
99
100 fix = NULL;
101 fmask = NULL;
102 list_pre_initial_integrate = list_initial_integrate = list_post_integrate = NULL;
103 list_pre_exchange = list_pre_neighbor = NULL;
104 list_pre_force = list_post_force = NULL;
105 list_iterate_implicitly = list_pre_final_integrate = NULL;
106 list_final_integrate = list_end_of_step = NULL;
107 list_thermo_energy = NULL;
108 list_initial_integrate_respa = list_post_integrate_respa = NULL;
109 list_pre_force_respa = list_post_force_respa = NULL;
110 list_final_integrate_respa = NULL;
111 list_min_pre_exchange = list_min_pre_neighbor = NULL;
112 list_min_pre_force = list_min_post_force = NULL;
113 list_min_energy = NULL;
114
115 end_of_step_every = NULL;
116
117 list_timeflag = NULL;
118
119 nfix_restart_global = 0;
120 id_restart_global = style_restart_global = state_restart_global = NULL;
121 nfix_restart_peratom = 0;
122 id_restart_peratom = style_restart_peratom = NULL;
123 index_restart_peratom = NULL;
124
125 ncompute = maxcompute = 0;
126 compute = NULL;
127
128 timing = 0;
129
130 // fill map with fixes listed in style_fix.h
131
132 fix_map = new std::map<std::string,FixCreator>();
133
134 #define FIX_CLASS
135 #define FixStyle(key,Class) \
136 (*fix_map)[#key] = &fix_creator<Class>;
137 #include "style_fix.h"
138 #undef FixStyle
139 #undef FIX_CLASS
140
141 // fill map with computes listed in style_compute.h
142
143 compute_map = new std::map<std::string,ComputeCreator>();
144
145 #define COMPUTE_CLASS
146 #define ComputeStyle(key,Class) \
147 (*compute_map)[#key] = &compute_creator<Class>;
148 #include "style_compute.h"
149 #undef ComputeStyle
150 #undef COMPUTE_CLASS
151 }
152
153 /* ---------------------------------------------------------------------- */
154
~Modify()155 Modify::~Modify()
156 {
157 // delete all fixes
158 // do it via delete_fix() so callbacks in Atom are also updated correctly
159
160 while (nfix) delete_fix(fix[0]->id);
161 memory->sfree(fix);
162 memory->destroy(fmask);
163
164 // delete all computes
165
166 for (int i = ncompute-1; i >= 0; i--) delete compute[i];
167 memory->sfree(compute);
168
169 delete [] list_pre_initial_integrate;
170 delete [] list_initial_integrate;
171 delete [] list_post_integrate;
172 delete [] list_pre_exchange;
173 delete [] list_pre_neighbor;
174 delete [] list_pre_force;
175 delete [] list_post_force;
176 delete [] list_pre_final_integrate;
177 delete [] list_final_integrate;
178 delete [] list_iterate_implicitly;
179 delete [] list_end_of_step;
180 delete [] list_thermo_energy;
181 delete [] list_initial_integrate_respa;
182 delete [] list_post_integrate_respa;
183 delete [] list_pre_force_respa;
184 delete [] list_post_force_respa;
185 delete [] list_final_integrate_respa;
186 delete [] list_min_pre_exchange;
187 delete [] list_min_pre_neighbor;
188 delete [] list_min_pre_force;
189 delete [] list_min_post_force;
190 delete [] list_min_energy;
191
192 delete [] end_of_step_every;
193 delete [] list_timeflag;
194
195 restart_deallocate();
196 delete compute_map;
197 delete fix_map;
198 }
199
200 /* ----------------------------------------------------------------------
201 initialize all fixes and computes
202 ------------------------------------------------------------------------- */
203
init()204 void Modify::init()
205 {
206 int i,j;
207
208 // delete storage of restart info since it is not valid after 1st run
209
210 restart_deallocate();
211
212 // create lists of fixes to call at each stage of run
213
214 list_init(PRE_INITIAL_INTEGRATE,n_pre_initial_integrate,list_pre_initial_integrate);
215 list_init(INITIAL_INTEGRATE,n_initial_integrate,list_initial_integrate);
216 list_init(POST_INTEGRATE,n_post_integrate,list_post_integrate);
217 list_init_pre_exchange(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange);
218 //list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange);
219 list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor);
220 list_init(PRE_FORCE,n_pre_force,list_pre_force);
221 list_init(POST_FORCE,n_post_force,list_post_force);
222 list_init(PRE_FINAL_INTEGRATE,n_pre_final_integrate,list_pre_final_integrate);
223 list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate);
224 list_init(ITERATE_IMPLICITLY,n_iterate_implicitly,list_iterate_implicitly);
225 list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step);
226 list_init_thermo_energy(THERMO_ENERGY,n_thermo_energy,list_thermo_energy);
227
228 list_init(INITIAL_INTEGRATE_RESPA,
229 n_initial_integrate_respa,list_initial_integrate_respa);
230 list_init(POST_INTEGRATE_RESPA,
231 n_post_integrate_respa,list_post_integrate_respa);
232 list_init(POST_FORCE_RESPA,
233 n_post_force_respa,list_post_force_respa);
234 list_init(PRE_FORCE_RESPA,
235 n_pre_force_respa,list_pre_force_respa);
236 list_init(FINAL_INTEGRATE_RESPA,
237 n_final_integrate_respa,list_final_integrate_respa);
238
239 list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange);
240 list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force);
241 list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
242 list_init(MIN_ENERGY,n_min_energy,list_min_energy);
243
244 // init each fix
245 // not sure if now needs to come before compute init
246 // used to b/c temperature computes called fix->dof() in their init,
247 // and fix rigid required its own init before its dof() could be called,
248 // but computes now do their DOF in setup()
249
250 for (i = 0; i < nfix; i++) fix[i]->init();
251
252 // set global flag if any fix has its restart_pbc flag set
253
254 restart_pbc_any = 0;
255 for (i = 0; i < nfix; i++)
256 if (fix[i]->restart_pbc) restart_pbc_any = 1;
257
258 // create list of computes that store invocation times
259
260 list_init_compute();
261
262 // init each compute
263 // set invoked_scalar,vector,etc to -1 to force new run to re-compute them
264 // add initial timestep to all computes that store invocation times
265 // since any of them may be invoked by initial thermo
266 // do not clear out invocation times stored within a compute,
267 // b/c some may be holdovers from previous run, like for ave fixes
268
269 for (i = 0; i < ncompute; i++) {
270 compute[i]->init();
271 compute[i]->invoked_scalar = -1;
272 compute[i]->invoked_vector = -1;
273 compute[i]->invoked_array = -1;
274 compute[i]->invoked_peratom = -1;
275 compute[i]->invoked_local = -1;
276 }
277 addstep_compute_all(update->ntimestep);
278
279 // warn if any particle is time integrated more than once
280
281 int nlocal = atom->nlocal;
282 int *mask = atom->mask;
283
284 int *flag = new int[nlocal];
285 for (i = 0; i < nlocal; i++) flag[i] = 0;
286
287 int groupbit;
288 for (i = 0; i < nfix; i++) {
289 if (fix[i]->time_integrate == 0) continue;
290 groupbit = fix[i]->groupbit;
291 for (j = 0; j < nlocal; j++)
292 if (mask[j] & groupbit) flag[j]++;
293 }
294
295 int check = 0;
296 for (i = 0; i < nlocal; i++)
297 if (flag[i] > 1) check = 1;
298
299 delete [] flag;
300
301 int checkall;
302 MPI_Allreduce(&check,&checkall,1,MPI_INT,MPI_SUM,world);
303 if (comm->me == 0 && checkall)
304 error->warning(FLERR,
305 "One or more atoms are time integrated more than once");
306 }
307
308 /* ----------------------------------------------------------------------
309 setup for run, calls setup() of all fixes and computes
310 called from Verlet, RESPA, Min
311 ------------------------------------------------------------------------- */
312
setup(int vflag)313 void Modify::setup(int vflag)
314 {
315
316 int nlocal = atom->nlocal;
317 int *mask = atom->mask;
318
319 for (int i = 0; i < nfix; i++)
320 {
321 if (!fix[i]->recent_restart && fix[i]->just_created && fix[i]->create_attribute)
322 {
323 fix[i]->just_created = 0;
324 fix[i]->pre_set_arrays();
325 for(int j = 0; j < nlocal; j++)
326 if(mask[j] & fix[i]->groupbit) fix[i]->set_arrays(j);
327 }
328 else if(fix[i]->just_created) fix[i]->just_created = 0;
329 fix[i]->recent_restart = 0;
330 }
331
332 // compute setup needs to come before fix setup
333 // b/c NH fixes need use DOF of temperature computes
334
335 for (int i = 0; i < ncompute; i++) compute[i]->setup();
336
337 if (update->whichflag == 1)
338 call_method_on_fixes(&Fix::setup, vflag);
339 else if (update->whichflag == 2)
340 call_method_on_fixes(&Fix::min_setup, vflag);
341 }
342
343 /* ----------------------------------------------------------------------
344 setup pre_exchange call, only for fixes that define pre_exchange
345 called from Verlet, RESPA, Min, and WriteRestart with whichflag = 0
346 ------------------------------------------------------------------------- */
347
setup_pre_exchange()348 void Modify::setup_pre_exchange()
349 {
350 if (update->whichflag <= 1)
351 call_method_on_fixes(&Fix::setup_pre_exchange, list_pre_exchange, n_pre_exchange);
352 else if (update->whichflag == 2)
353 call_method_on_fixes(&Fix::min_setup_pre_exchange, list_min_pre_exchange, n_min_pre_exchange);
354 }
355
356 /* ----------------------------------------------------------------------
357 setup pre_neighbor call, only for fixes that define pre_neighbor
358 called from Verlet, RESPA
359 ------------------------------------------------------------------------- */
360
setup_pre_neighbor()361 void Modify::setup_pre_neighbor()
362 {
363 if (update->whichflag == 1)
364 call_method_on_fixes(&Fix::setup_pre_neighbor, list_pre_neighbor, n_pre_neighbor);
365 else if (update->whichflag == 2)
366 call_method_on_fixes(&Fix::min_setup_pre_neighbor, list_min_pre_neighbor, n_min_pre_neighbor);
367 }
368
369 /* ----------------------------------------------------------------------
370 setup pre_force call, only for fixes that define pre_force
371 called from Verlet, RESPA, Min
372 ------------------------------------------------------------------------- */
373
setup_pre_force(int vflag)374 void Modify::setup_pre_force(int vflag)
375 {
376 if (update->whichflag == 1) {
377 call_method_on_fixes(&Fix::setup_pre_force, vflag, list_pre_force, n_pre_force);
378 } else if (update->whichflag == 2) {
379 call_method_on_fixes(&Fix::min_setup_pre_force, vflag, list_min_pre_force, n_min_pre_force);
380 }
381 }
382
383 /* ----------------------------------------------------------------------
384 pre_intial_integrate call, only for relevant fixes
385 ------------------------------------------------------------------------- */
386
pre_initial_integrate()387 void Modify::pre_initial_integrate()
388 {
389 call_method_on_fixes(&Fix::pre_initial_integrate, list_pre_initial_integrate, n_pre_initial_integrate);
390 }
391
392 /* ----------------------------------------------------------------------
393 1st half of integrate call, only for relevant fixes
394 ------------------------------------------------------------------------- */
395
initial_integrate(int vflag)396 void Modify::initial_integrate(int vflag)
397 {
398
399 call_method_on_fixes(&Fix::initial_integrate, vflag, list_initial_integrate, n_initial_integrate);
400 }
401
402 /* ----------------------------------------------------------------------
403 post_integrate call, only for relevant fixes
404 ------------------------------------------------------------------------- */
405
post_integrate()406 void Modify::post_integrate()
407 {
408 call_method_on_fixes(&Fix::post_integrate, list_post_integrate, n_post_integrate);
409 }
410
411 /* ----------------------------------------------------------------------
412 pre_exchange call, only for relevant fixes
413 ------------------------------------------------------------------------- */
414
pre_exchange()415 void Modify::pre_exchange()
416 {
417
418 call_method_on_fixes(&Fix::pre_exchange, list_pre_exchange, n_pre_exchange);
419 }
420
421 /* ----------------------------------------------------------------------
422 pre_neighbor call, only for relevant fixes
423 ------------------------------------------------------------------------- */
424
pre_neighbor()425 void Modify::pre_neighbor()
426 {
427
428 call_method_on_fixes(&Fix::pre_neighbor, list_pre_neighbor, n_pre_neighbor);
429 }
430
431 /* ----------------------------------------------------------------------
432 pre_force call, only for relevant fixes
433 ------------------------------------------------------------------------- */
434
pre_force(int vflag)435 void Modify::pre_force(int vflag)
436 {
437
438 call_method_on_fixes(&Fix::pre_force, vflag, list_pre_force, n_pre_force);
439 }
440
441 /* ----------------------------------------------------------------------
442 post_force call, only for relevant fixes
443 ------------------------------------------------------------------------- */
444
post_force(int vflag)445 void Modify::post_force(int vflag)
446 {
447
448 call_method_on_fixes(&Fix::post_force, vflag, list_post_force, n_post_force);
449 }
450
451 /* ----------------------------------------------------------------------
452 check convergence, only for relevant implicit integration
453 ------------------------------------------------------------------------- */
454
iterate_implicitly()455 bool Modify::iterate_implicitly()
456 {
457 for (int i = 0; i < n_iterate_implicitly; i++)
458 if (fix[list_iterate_implicitly[i]]->iterate_implicitly())
459 return true;
460
461 return false;
462 }
463
464 /* ----------------------------------------------------------------------
465 pre_final_integrate call, only for relevant fixes
466 ------------------------------------------------------------------------- */
467
pre_final_integrate()468 void Modify::pre_final_integrate()
469 {
470 call_method_on_fixes(&Fix::pre_final_integrate, list_pre_final_integrate, n_pre_final_integrate);
471 }
472
473 /* ----------------------------------------------------------------------
474 2nd half of integrate call, only for relevant fixes
475 ------------------------------------------------------------------------- */
476
final_integrate()477 void Modify::final_integrate()
478 {
479 call_method_on_fixes(&Fix::final_integrate, list_final_integrate, n_final_integrate);
480 }
481
482 /* ----------------------------------------------------------------------
483 end-of-timestep call, only for relevant fixes
484 only call fix->end_of_step() on timesteps that are multiples of nevery
485 ------------------------------------------------------------------------- */
486
end_of_step()487 void Modify::end_of_step()
488 {
489 if(timing) {
490 for (int i = 0; i < n_end_of_step; i++) {
491 if (update->ntimestep % end_of_step_every[i] == 0) {
492 const int ifix = list_end_of_step[i];
493 fix[ifix]->begin_time_recording();
494 fix[ifix]->end_of_step();
495 fix[ifix]->end_time_recording();
496 }
497 }
498 }
499 else
500 {
501 for (int i = 0; i < n_end_of_step; i++) {
502 if (update->ntimestep % end_of_step_every[i] == 0) {
503 fix[list_end_of_step[i]]->end_of_step();
504 }
505 }
506 }
507 }
508
509 /* ----------------------------------------------------------------------
510 thermo energy call, only for relevant fixes
511 called by Thermo class
512 compute_scalar() is fix call to return energy
513 ------------------------------------------------------------------------- */
514
thermo_energy()515 double Modify::thermo_energy()
516 {
517 double energy = 0.0;
518 if(timing) {
519 for (int i = 0; i < n_thermo_energy; i++) {
520 const int ifix = list_thermo_energy[i];
521 fix[ifix]->begin_time_recording();
522 energy += fix[ifix]->compute_scalar();
523 fix[ifix]->end_time_recording();
524 }
525 }
526 else
527 {
528 for (int i = 0; i < n_thermo_energy; i++) {
529 energy += fix[list_thermo_energy[i]]->compute_scalar();
530 }
531 }
532 return energy;
533 }
534
535 /* ----------------------------------------------------------------------
536 post_run call
537 ------------------------------------------------------------------------- */
538
post_run()539 void Modify::post_run()
540 {
541 call_method_on_fixes(&Fix::post_run);
542 }
543
544 /* ----------------------------------------------------------------------
545 setup rRESPA pre_force call, only for relevant fixes
546 ------------------------------------------------------------------------- */
547
setup_pre_force_respa(int vflag,int ilevel)548 void Modify::setup_pre_force_respa(int vflag, int ilevel)
549 {
550 call_respa_method_on_fixes(&Fix::setup_pre_force_respa, vflag, ilevel,
551 list_pre_force_respa, n_pre_force_respa);
552 }
553
554 /* ----------------------------------------------------------------------
555 1st half of rRESPA integrate call, only for relevant fixes
556 ------------------------------------------------------------------------- */
557
initial_integrate_respa(int vflag,int ilevel,int iloop)558 void Modify::initial_integrate_respa(int vflag, int ilevel, int iloop)
559 {
560 call_respa_method_on_fixes(&Fix::initial_integrate_respa, vflag, ilevel, iloop,
561 list_initial_integrate_respa, n_initial_integrate_respa);
562 }
563
564 /* ----------------------------------------------------------------------
565 rRESPA post_integrate call, only for relevant fixes
566 ------------------------------------------------------------------------- */
567
post_integrate_respa(int ilevel,int iloop)568 void Modify::post_integrate_respa(int ilevel, int iloop)
569 {
570 call_respa_method_on_fixes(&Fix::post_integrate_respa, ilevel, iloop,
571 list_post_integrate_respa, n_post_integrate_respa);
572 }
573
574 /* ----------------------------------------------------------------------
575 rRESPA pre_force call, only for relevant fixes
576 ------------------------------------------------------------------------- */
577
pre_force_respa(int vflag,int ilevel,int iloop)578 void Modify::pre_force_respa(int vflag, int ilevel, int iloop)
579 {
580 call_respa_method_on_fixes(&Fix::pre_force_respa, vflag, ilevel, iloop,
581 list_pre_force_respa, n_pre_force_respa);
582 }
583
584 /* ----------------------------------------------------------------------
585 rRESPA post_force call, only for relevant fixes
586 ------------------------------------------------------------------------- */
587
post_force_respa(int vflag,int ilevel,int iloop)588 void Modify::post_force_respa(int vflag, int ilevel, int iloop)
589 {
590 call_respa_method_on_fixes(&Fix::post_force_respa, vflag, ilevel, iloop,
591 list_post_force_respa, n_post_force_respa);
592 }
593
594 /* ----------------------------------------------------------------------
595 2nd half of rRESPA integrate call, only for relevant fixes
596 ------------------------------------------------------------------------- */
597
final_integrate_respa(int ilevel,int iloop)598 void Modify::final_integrate_respa(int ilevel, int iloop)
599 {
600 call_respa_method_on_fixes(&Fix::final_integrate_respa, ilevel, iloop,
601 list_final_integrate_respa, n_final_integrate_respa);
602 }
603
604 /* ----------------------------------------------------------------------
605 minimizer pre-exchange call, only for relevant fixes
606 ------------------------------------------------------------------------- */
607
min_pre_exchange()608 void Modify::min_pre_exchange()
609 {
610 call_method_on_fixes(&Fix::min_pre_exchange, list_min_pre_exchange, n_min_pre_exchange);
611 }
612
613 /* ----------------------------------------------------------------------
614 minimizer pre-neighbor call, only for relevant fixes
615 ------------------------------------------------------------------------- */
616
min_pre_neighbor()617 void Modify::min_pre_neighbor()
618 {
619 for (int i = 0; i < n_min_pre_neighbor; i++)
620 fix[list_min_pre_neighbor[i]]->min_pre_neighbor();
621 }
622
623 /* ----------------------------------------------------------------------
624 minimizer pre-force call, only for relevant fixes
625 ------------------------------------------------------------------------- */
626
min_pre_force(int vflag)627 void Modify::min_pre_force(int vflag)
628 {
629 call_method_on_fixes(&Fix::min_pre_force, vflag, list_min_pre_force, n_min_pre_force);
630 }
631
632 /* ----------------------------------------------------------------------
633 minimizer force adjustment call, only for relevant fixes
634 ------------------------------------------------------------------------- */
635
min_post_force(int vflag)636 void Modify::min_post_force(int vflag)
637 {
638 call_method_on_fixes(&Fix::min_post_force, vflag, list_min_post_force, n_min_post_force);
639 }
640
641 /* ----------------------------------------------------------------------
642 minimizer energy/force evaluation, only for relevant fixes
643 return energy and forces on extra degrees of freedom
644 ------------------------------------------------------------------------- */
645
min_energy(double * fextra)646 double Modify::min_energy(double *fextra)
647 {
648 int ifix,index;
649
650 index = 0;
651 double eng = 0.0;
652 for (int i = 0; i < n_min_energy; i++) {
653 ifix = list_min_energy[i];
654 eng += fix[ifix]->min_energy(&fextra[index]);
655 index += fix[ifix]->min_dof();
656 }
657 return eng;
658 }
659
660 /* ----------------------------------------------------------------------
661 store current state of extra dof, only for relevant fixes
662 ------------------------------------------------------------------------- */
663
min_store()664 void Modify::min_store()
665 {
666 call_method_on_fixes(&Fix::min_store, list_min_energy, n_min_energy);
667 }
668
669 /* ----------------------------------------------------------------------
670 mange state of extra dof on a stack, only for relevant fixes
671 ------------------------------------------------------------------------- */
672
min_clearstore()673 void Modify::min_clearstore()
674 {
675 call_method_on_fixes(&Fix::min_clearstore, list_min_energy, n_min_energy);
676 }
677
min_pushstore()678 void Modify::min_pushstore()
679 {
680 call_method_on_fixes(&Fix::min_pushstore, list_min_energy, n_min_energy);
681 }
682
min_popstore()683 void Modify::min_popstore()
684 {
685 call_method_on_fixes(&Fix::min_popstore, list_min_energy, n_min_energy);
686 }
687
688 /* ----------------------------------------------------------------------
689 displace extra dof along vector hextra, only for relevant fixes
690 ------------------------------------------------------------------------- */
691
min_step(double alpha,double * hextra)692 void Modify::min_step(double alpha, double *hextra)
693 {
694 int ifix,index;
695
696 index = 0;
697 for (int i = 0; i < n_min_energy; i++) {
698 ifix = list_min_energy[i];
699 fix[ifix]->min_step(alpha,&hextra[index]);
700 index += fix[ifix]->min_dof();
701 }
702 }
703
704 /* ----------------------------------------------------------------------
705 compute max allowed step size along vector hextra, only for relevant fixes
706 ------------------------------------------------------------------------- */
707
max_alpha(double * hextra)708 double Modify::max_alpha(double *hextra)
709 {
710 int ifix,index;
711
712 double alpha = BIG;
713 index = 0;
714 for (int i = 0; i < n_min_energy; i++) {
715 ifix = list_min_energy[i];
716 double alpha_one = fix[ifix]->max_alpha(&hextra[index]);
717 alpha = MIN(alpha,alpha_one);
718 index += fix[ifix]->min_dof();
719 }
720 return alpha;
721 }
722
723 /* ----------------------------------------------------------------------
724 extract extra dof for minimization, only for relevant fixes
725 ------------------------------------------------------------------------- */
726
min_dof()727 int Modify::min_dof()
728 {
729 int ndof = 0;
730 for (int i = 0; i < n_min_energy; i++)
731 ndof += fix[list_min_energy[i]]->min_dof();
732 return ndof;
733 }
734
735 /* ----------------------------------------------------------------------
736 reset reference state of fix, only for relevant fixes
737 ------------------------------------------------------------------------- */
738
min_reset_ref()739 int Modify::min_reset_ref()
740 {
741 int itmp,itmpall;
742 itmpall = 0;
743 for (int i = 0; i < n_min_energy; i++) {
744 itmp = fix[list_min_energy[i]]->min_reset_ref();
745 if (itmp) itmpall = 1;
746 }
747 return itmpall;
748 }
749
750 /* ----------------------------------------------------------------------
751 add a new fix or replace one with same ID
752 ------------------------------------------------------------------------- */
753
add_fix(int narg,char ** arg,char * suffix)754 void Modify::add_fix(int narg, char **arg, char *suffix)
755 {
756
757 if(update->ntimestep_reset_since_last_run && fix_restart_in_progress())
758 error->all(FLERR,"In case of restart, command 'reset_timestep' must come immediately before 'run'");
759
760 if (narg < 3) error->all(FLERR,"Illegal fix command");
761
762 // cannot define fix before box exists unless style is in exception list
763 // don't like this way of checking for exceptions by adding fixes to list,
764 // but can't think of better way
765 // too late if instantiate fix, then check flag set in fix constructor,
766 // since some fixes access domain settings in their constructor
767 // change NEXCEPT above when add new fix to this list
768
769 const char *exceptions[NEXCEPT] = {"GPU","OMP","property/atom","cmap"};
770
771 if (domain->box_exist == 0) {
772 int m;
773 for (m = 0; m < NEXCEPT; m++)
774 if (strcmp(arg[2],exceptions[m]) == 0) break;
775 if (m == NEXCEPT)
776 error->all(FLERR,"Fix command before simulation box is defined");
777 }
778
779 // check group ID
780
781 int igroup = group->find(arg[1]);
782 if (igroup == -1) error->all(FLERR,"Could not find fix group ID");
783
784 // if fix ID exists:
785 // set newflag = 0 so create new fix in same location in fix list
786 // error if new style does not match old style
787 // since can't replace it (all when-to-invoke ptrs would be invalid)
788 // warn if new group != old group
789 // delete old fix, but do not call update_callback(),
790 // since will replace this fix and thus other fix locs will not change
791 // set ptr to NULL in case new fix scans list of fixes,
792 // e.g. scan will occur in add_callback() if called by new fix
793 // if fix ID does not exist:
794 // set newflag = 1 so create new fix
795 // extend fix and fmask lists as necessary
796
797 int ifix,newflag;
798 for (ifix = 0; ifix < nfix; ifix++)
799 if (strcmp(arg[0],fix[ifix]->id) == 0) break;
800
801 if (ifix < nfix) {
802 newflag = 0;
803
804 if (strncmp(fix[ifix]->style,"insert/",7) == 0)
805 error->all(FLERR,"Using a fix insert/* ID twice, which is not possible. Please use different ID");
806 if (strcmp(arg[2],fix[ifix]->style) != 0)
807 error->all(FLERR,"Replacing a fix, but new style != old style");
808 if (fix[ifix]->igroup != igroup && comm->me == 0)
809 error->warning(FLERR,"Replacing a fix, but new group != old group");
810
811 fix[ifix]->pre_delete(true);
812 delete fix[ifix];
813 fix[ifix] = NULL;
814 } else {
815 newflag = 1;
816 if (nfix == maxfix) {
817 maxfix += DELTA;
818 fix = (Fix **) memory->srealloc(fix,maxfix*sizeof(Fix *),"modify:fix");
819 memory->grow(fmask,maxfix,"modify:fmask");
820 }
821 }
822
823 // create the Fix
824 // try first with suffix appended
825
826 fix[ifix] = NULL;
827
828 if (suffix && lmp->suffix_enable) {
829 char estyle[256];
830 sprintf(estyle,"%s/%s",arg[2],suffix);
831
832 if (fix_map->find(estyle) != fix_map->end()) {
833 FixCreator fix_creator = (*fix_map)[estyle];
834 fix[ifix] = fix_creator(lmp,narg,arg);
835 }
836 }
837
838 if (fix[ifix] == NULL && fix_map->find(arg[2]) != fix_map->end()) {
839 FixCreator fix_creator = (*fix_map)[arg[2]];
840 fix[ifix] = fix_creator(lmp,narg,arg);
841 }
842
843 if (fix[ifix] == NULL && strncmp(arg[2], "mesh/surface", 12) == 0)
844 {
845 FixCreator fix_creator = (*fix_map)["mesh/surface"];
846 fix[ifix] = fix_creator(lmp, narg, arg);
847 }
848
849 if (fix[ifix] == NULL){
850 char * errmsg = new char[30+strlen(arg[2])];
851 sprintf(errmsg,"Invalid fix style: \"%s\"",arg[2]);
852 error->all(FLERR,errmsg);
853 delete [] errmsg;
854 }
855
856 // set fix mask values and increment nfix (if new)
857
858 fmask[ifix] = fix[ifix]->setmask();
859 if (newflag) nfix++;
860
861 fix[ifix]->post_create_pre_restart();
862
863 // check if Fix is in restart_global list
864 // if yes, pass state info to the Fix so it can reset itself
865
866 for (int i = 0; i < nfix_restart_global; i++)
867 {
868
869 if (strcmp(id_restart_global[i],fix[ifix]->id) == 0 &&
870 (strcmp(style_restart_global[i],fix[ifix]->style) == 0 ||
871 (fix[ifix]->accepts_restart_data_from_style && strcmp(style_restart_global[i],fix[ifix]->accepts_restart_data_from_style) == 0)
872 )
873 )
874 {
875 fix[ifix]->restart(state_restart_global[i]);
876 fix[ifix]->recent_restart = 1;
877 if (comm->me == 0) {
878 char *str = (char *) ("Resetting global state of Fix %s Style %s "
879 "from restart file info\n");
880 if (screen) fprintf(screen,str,fix[ifix]->id,fix[ifix]->style);
881 if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style);
882 }
883 }
884 }
885
886 // check if Fix is in restart_peratom list
887 // if yes, loop over atoms so they can extract info from atom->extra array
888
889 for (int i = 0; i < nfix_restart_peratom; i++)
890 {
891 if (strcmp(id_restart_peratom[i],fix[ifix]->id) == 0 &&
892 (strcmp(style_restart_peratom[i],fix[ifix]->style) == 0 ||
893 (fix[ifix]->accepts_restart_data_from_style && strcmp(style_restart_peratom[i],fix[ifix]->accepts_restart_data_from_style) == 0)
894 )
895 )
896 {
897 for (int j = 0; j < atom->nlocal; j++)
898 fix[ifix]->unpack_restart(j,index_restart_peratom[i]);
899 fix[ifix]->recent_restart = 1;
900 fix[ifix]->restart_reset = 1;
901 if (comm->me == 0) {
902 char *str = (char *) ("Resetting per-atom state of Fix %s Style %s "
903 "from restart file info\n");
904 if (screen) fprintf(screen,str,fix[ifix]->id,fix[ifix]->style);
905 if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style);
906
907 }
908 }
909 }
910
911 fix[ifix]->post_create();
912 }
913
914 /* ----------------------------------------------------------------------
915 one instance per fix in style_fix.h
916 ------------------------------------------------------------------------- */
917
918 template <typename T>
fix_creator(LAMMPS * lmp,int narg,char ** arg)919 Fix *Modify::fix_creator(LAMMPS *lmp, int narg, char **arg)
920 {
921 return new T(lmp,narg,arg);
922 }
923
924 /* ----------------------------------------------------------------------
925 modify a Fix's parameters
926 ------------------------------------------------------------------------- */
927
modify_fix(int narg,char ** arg)928 void Modify::modify_fix(int narg, char **arg)
929 {
930 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
931
932 // lookup Fix ID
933
934 int ifix;
935 for (ifix = 0; ifix < nfix; ifix++)
936 if (strcmp(arg[0],fix[ifix]->id) == 0) break;
937 if (ifix == nfix) error->all(FLERR,"Could not find fix_modify ID");
938
939 fix[ifix]->modify_params(narg-1,&arg[1]);
940 }
941
942 /* ----------------------------------------------------------------------
943 delete a Fix from list of Fixes
944 Atom class must update indices in its list of callbacks to fixes
945 ------------------------------------------------------------------------- */
946
delete_fix(const char * id,bool unfixflag)947 void Modify::delete_fix(const char *id, bool unfixflag)
948 {
949 int ifix = find_fix(id);
950 if (ifix < 0) error->all(FLERR,"Could not find fix ID to delete");
951
952 fix[ifix]->pre_delete(unfixflag);
953
954 delete fix[ifix];
955 fix[ifix] = NULL;
956 atom->update_callback(ifix);
957
958 // move other Fixes and fmask down in list one slot
959
960 for (int i = ifix+1; i < nfix; i++) fix[i-1] = fix[i];
961 for (int i = ifix+1; i < nfix; i++) fmask[i-1] = fmask[i];
962 nfix--;
963 }
964
965 /* ----------------------------------------------------------------------
966 find a fix by ID
967 return index of fix or -1 if not found
968 ------------------------------------------------------------------------- */
969
find_fix(const char * id)970 int Modify::find_fix(const char *id)
971 {
972 int ifix;
973 for (ifix = 0; ifix < nfix; ifix++)
974 if (strcmp(id,fix[ifix]->id) == 0) break;
975 if (ifix == nfix) return -1;
976 return ifix;
977 }
978
979 /* ----------------------------------------------------------------------
980 add a new compute
981 ------------------------------------------------------------------------- */
982
add_compute(int narg,char ** arg,char * suffix)983 void Modify::add_compute(int narg, char **arg, char *suffix)
984 {
985 if (narg < 3) error->all(FLERR,"Illegal compute command");
986
987 // error check
988
989 for (int icompute = 0; icompute < ncompute; icompute++)
990 if (strcmp(arg[0],compute[icompute]->id) == 0)
991 error->all(FLERR,"Reuse of compute ID");
992
993 // extend Compute list if necessary
994
995 if (ncompute == maxcompute) {
996 maxcompute += DELTA;
997 compute = (Compute **)
998 memory->srealloc(compute,maxcompute*sizeof(Compute *),"modify:compute");
999 }
1000
1001 // create the Compute
1002 // try first with suffix appended
1003
1004 compute[ncompute] = NULL;
1005
1006 if (suffix && lmp->suffix_enable)
1007 {
1008 char estyle[256];
1009 sprintf(estyle,"%s/%s",arg[2],suffix);
1010 if (compute_map->find(estyle) != compute_map->end())
1011 {
1012 ComputeCreator compute_creator = (*compute_map)[estyle];
1013 int iarg = 0;
1014 compute[ncompute] = compute_creator(lmp, iarg, narg, arg);
1015 }
1016 }
1017
1018 if (compute[ncompute] == NULL &&
1019 compute_map->find(arg[2]) != compute_map->end())
1020 {
1021 int iarg = 0;
1022 ComputeCreator compute_creator = (*compute_map)[arg[2]];
1023 compute[ncompute] = compute_creator(lmp, iarg, narg, arg);
1024 }
1025
1026 if (compute[ncompute] == NULL) error->all(FLERR,"Invalid compute style");
1027
1028 ncompute++;
1029
1030 compute[ncompute-1]->post_create();
1031 }
1032
1033 /* ----------------------------------------------------------------------
1034 one instance per compute in style_compute.h
1035 ------------------------------------------------------------------------- */
1036
1037 template <typename T>
compute_creator(LAMMPS * lmp,int iarg,int narg,char ** arg)1038 Compute *Modify::compute_creator(LAMMPS *lmp, int iarg, int narg, char **arg)
1039 {
1040 return new T(lmp, iarg, narg,arg);
1041 }
1042
1043 /* ----------------------------------------------------------------------
1044 modify a Compute's parameters
1045 ------------------------------------------------------------------------- */
1046
modify_compute(int narg,char ** arg)1047 void Modify::modify_compute(int narg, char **arg)
1048 {
1049 if (narg < 2) error->all(FLERR,"Illegal compute_modify command");
1050
1051 // lookup Compute ID
1052
1053 int icompute;
1054 for (icompute = 0; icompute < ncompute; icompute++)
1055 if (strcmp(arg[0],compute[icompute]->id) == 0) break;
1056 if (icompute == ncompute)
1057 error->all(FLERR,"Could not find compute_modify ID");
1058
1059 compute[icompute]->modify_params(narg-1,&arg[1]);
1060 }
1061
1062 /* ----------------------------------------------------------------------
1063 delete a Compute from list of Computes
1064 ------------------------------------------------------------------------- */
1065
delete_compute(const char * id,bool uncomputeflag)1066 void Modify::delete_compute(const char *id,bool uncomputeflag)
1067 {
1068 int icompute = find_compute(id);
1069 if (icompute < 0) error->all(FLERR,"Could not find compute ID to delete");
1070
1071 compute[icompute]->pre_delete(uncomputeflag);
1072
1073 delete compute[icompute];
1074
1075 // move other Computes down in list one slot
1076
1077 for (int i = icompute+1; i < ncompute; i++) compute[i-1] = compute[i];
1078 ncompute--;
1079 }
1080
1081 /* ----------------------------------------------------------------------
1082 find a compute by ID
1083 return index of compute or -1 if not found
1084 ------------------------------------------------------------------------- */
1085
find_compute(const char * id)1086 int Modify::find_compute(const char *id)
1087 {
1088 int icompute;
1089 for (icompute = 0; icompute < ncompute; icompute++)
1090 if (strcmp(id,compute[icompute]->id) == 0) break;
1091 if (icompute == ncompute) return -1;
1092 return icompute;
1093 }
1094
1095 /* ----------------------------------------------------------------------
1096 clear invoked flag of all computes
1097 called everywhere that computes are used, before computes are invoked
1098 invoked flag used to avoid re-invoking same compute multiple times
1099 and to flag computes that store invocation times as having been invoked
1100 ------------------------------------------------------------------------- */
1101
clearstep_compute()1102 void Modify::clearstep_compute()
1103 {
1104 for (int icompute = 0; icompute < ncompute; icompute++)
1105 compute[icompute]->invoked_flag = 0;
1106 }
1107
1108 /* ----------------------------------------------------------------------
1109 loop over computes that store invocation times
1110 if its invoked flag set on this timestep, schedule next invocation
1111 called everywhere that computes are used, after computes are invoked
1112 ------------------------------------------------------------------------- */
1113
addstep_compute(bigint newstep)1114 void Modify::addstep_compute(bigint newstep)
1115 {
1116 if (list_timeflag == nullptr)
1117 list_init_compute();
1118 for (int icompute = 0; icompute < n_timeflag; icompute++)
1119 if (compute[list_timeflag[icompute]]->invoked_flag)
1120 compute[list_timeflag[icompute]]->addstep(newstep);
1121 }
1122
1123 /* ----------------------------------------------------------------------
1124 loop over all computes
1125 schedule next invocation for those that store invocation times
1126 called when not sure what computes will be needed on newstep
1127 do not loop only over n_timeflag, since may not be set yet
1128 ------------------------------------------------------------------------- */
1129
addstep_compute_all(bigint newstep)1130 void Modify::addstep_compute_all(bigint newstep)
1131 {
1132 for (int icompute = 0; icompute < ncompute; icompute++)
1133 if (compute[icompute]->timeflag) compute[icompute]->addstep(newstep);
1134 }
1135
1136 /* ----------------------------------------------------------------------
1137 write to restart file for all Fixes with restart info
1138 (1) fixes that have global state
1139 (2) fixes that store per-atom quantities
1140 ------------------------------------------------------------------------- */
1141
write_restart(FILE * fp)1142 void Modify::write_restart(FILE *fp)
1143 {
1144 int me = comm->me;
1145
1146 int count = 0;
1147 for (int i = 0; i < nfix; i++)
1148 if (fix[i]->restart_global) count++;
1149
1150 if (me == 0) fwrite(&count,sizeof(int),1,fp);
1151
1152 int n;
1153 for (int i = 0; i < nfix; i++)
1154 if (fix[i]->restart_global) {
1155 if (me == 0) {
1156
1157 n = strlen(fix[i]->id) + 1;
1158 fwrite(&n,sizeof(int),1,fp);
1159 fwrite(fix[i]->id,sizeof(char),n,fp);
1160 n = strlen(fix[i]->style) + 1;
1161 fwrite(&n,sizeof(int),1,fp);
1162 fwrite(fix[i]->style,sizeof(char),n,fp);
1163 }
1164 fix[i]->write_restart(fp);
1165 }
1166
1167 count = 0;
1168 for (int i = 0; i < nfix; i++)
1169 if (fix[i]->restart_peratom) count++;
1170
1171 if (me == 0) fwrite(&count,sizeof(int),1,fp);
1172
1173 for (int i = 0; i < nfix; i++)
1174 if (fix[i]->restart_peratom) {
1175 int maxsize_restart = fix[i]->maxsize_restart();
1176 if (me == 0) {
1177
1178 n = strlen(fix[i]->id) + 1;
1179 fwrite(&n,sizeof(int),1,fp);
1180 fwrite(fix[i]->id,sizeof(char),n,fp);
1181 n = strlen(fix[i]->style) + 1;
1182 fwrite(&n,sizeof(int),1,fp);
1183 fwrite(fix[i]->style,sizeof(char),n,fp);
1184 fwrite(&maxsize_restart,sizeof(int),1,fp);
1185 }
1186 }
1187 }
1188
1189 /* ----------------------------------------------------------------------
1190 read in restart file data on all previously defined Fixes with restart info
1191 (1) fixes that have global state
1192 (2) fixes that store per-atom quantities
1193 return maxsize of extra info that will be stored with any atom
1194 ------------------------------------------------------------------------- */
1195
read_restart(FILE * fp)1196 int Modify::read_restart(FILE *fp)
1197 {
1198 // nfix_restart_global = # of restart entries with global state info
1199
1200 int me = comm->me;
1201 if (me == 0) fread(&nfix_restart_global,sizeof(int),1,fp);
1202 MPI_Bcast(&nfix_restart_global,1,MPI_INT,0,world);
1203
1204 // allocate space for each entry
1205
1206 if (nfix_restart_global) {
1207 id_restart_global = new char*[nfix_restart_global];
1208 style_restart_global = new char*[nfix_restart_global];
1209 state_restart_global = new char*[nfix_restart_global];
1210 }
1211
1212 // read each entry and Bcast to all procs
1213 // each entry has id string, style string, chunk of state data
1214
1215 int n;
1216 for (int i = 0; i < nfix_restart_global; i++) {
1217 if (me == 0) fread(&n,sizeof(int),1,fp);
1218 MPI_Bcast(&n,1,MPI_INT,0,world);
1219 id_restart_global[i] = new char[n];
1220 if (me == 0) fread(id_restart_global[i],sizeof(char),n,fp);
1221 MPI_Bcast(id_restart_global[i],n,MPI_CHAR,0,world);
1222
1223 if (me == 0) fread(&n,sizeof(int),1,fp);
1224 MPI_Bcast(&n,1,MPI_INT,0,world);
1225 style_restart_global[i] = new char[n];
1226 if (me == 0) fread(style_restart_global[i],sizeof(char),n,fp);
1227 MPI_Bcast(style_restart_global[i],n,MPI_CHAR,0,world);
1228
1229 if (me == 0) fread(&n,sizeof(int),1,fp);
1230 MPI_Bcast(&n,1,MPI_INT,0,world);
1231 state_restart_global[i] = new char[n];
1232 if (me == 0) fread(state_restart_global[i],sizeof(char),n,fp);
1233 MPI_Bcast(state_restart_global[i],n,MPI_CHAR,0,world);
1234 }
1235
1236 // nfix_restart_peratom = # of restart entries with peratom info
1237
1238 int maxsize = 0;
1239
1240 if (me == 0) fread(&nfix_restart_peratom,sizeof(int),1,fp);
1241 MPI_Bcast(&nfix_restart_peratom,1,MPI_INT,0,world);
1242
1243 // allocate space for each entry
1244
1245 if (nfix_restart_peratom) {
1246 id_restart_peratom = new char*[nfix_restart_peratom];
1247 style_restart_peratom = new char*[nfix_restart_peratom];
1248 index_restart_peratom = new int[nfix_restart_peratom];
1249 }
1250
1251 // read each entry and Bcast to all procs
1252 // each entry has id string, style string, maxsize of one atom's data
1253 // set index = which set of extra data this fix represents
1254
1255 for (int i = 0; i < nfix_restart_peratom; i++) {
1256 if (me == 0) fread(&n,sizeof(int),1,fp);
1257 MPI_Bcast(&n,1,MPI_INT,0,world);
1258 id_restart_peratom[i] = new char[n];
1259 if (me == 0) fread(id_restart_peratom[i],sizeof(char),n,fp);
1260 MPI_Bcast(id_restart_peratom[i],n,MPI_CHAR,0,world);
1261
1262 if (me == 0) fread(&n,sizeof(int),1,fp);
1263 MPI_Bcast(&n,1,MPI_INT,0,world);
1264 style_restart_peratom[i] = new char[n];
1265 if (me == 0) fread(style_restart_peratom[i],sizeof(char),n,fp);
1266 MPI_Bcast(style_restart_peratom[i],n,MPI_CHAR,0,world);
1267
1268 if (me == 0) fread(&n,sizeof(int),1,fp);
1269 MPI_Bcast(&n,1,MPI_INT,0,world);
1270 maxsize += n;
1271
1272 index_restart_peratom[i] = i;
1273
1274 }
1275
1276 return maxsize;
1277 }
1278
1279 /* ----------------------------------------------------------------------
1280 delete all lists of restart file Fix info
1281 ------------------------------------------------------------------------- */
1282
restart_deallocate()1283 void Modify::restart_deallocate()
1284 {
1285
1286 int n_ms = n_fixes_style("multisphere");
1287 bool have_ms_in_restart = false;
1288
1289 if (nfix_restart_global) {
1290 for (int i = 0; i < nfix_restart_global; i++) {
1291 if(strncmp(style_restart_global[i],"multisphere",11) == 0)
1292 have_ms_in_restart = true;
1293 delete [] id_restart_global[i];
1294 delete [] style_restart_global[i];
1295 delete [] state_restart_global[i];
1296 }
1297 delete [] id_restart_global;
1298 delete [] style_restart_global;
1299 delete [] state_restart_global;
1300 }
1301
1302 if (nfix_restart_peratom) {
1303 for (int i = 0; i < nfix_restart_peratom; i++) {
1304 if(strncmp(style_restart_peratom[i],"multisphere",11) == 0)
1305 have_ms_in_restart = true;
1306 delete [] id_restart_peratom[i];
1307 delete [] style_restart_peratom[i];
1308 }
1309 delete [] id_restart_peratom;
1310 delete [] style_restart_peratom;
1311 delete [] index_restart_peratom;
1312 }
1313
1314 nfix_restart_global = nfix_restart_peratom = 0;
1315
1316 if(0 == n_ms && have_ms_in_restart)
1317 error->all(FLERR,"Restart data contains multi-sphere data, which was not restarted. In order to restart it,\n"
1318 "you have to place a fix multisphere/* command before the first run command in the input script\n");
1319 }
1320
1321 /* ----------------------------------------------------------------------
1322 create list of fix indices for fixes which match mask
1323 ------------------------------------------------------------------------- */
1324
list_init(int mask,int & n,int * & list)1325 void Modify::list_init(int mask, int &n, int *&list)
1326 {
1327 delete [] list;
1328
1329 n = 0;
1330 for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++;
1331 list = new int[n];
1332
1333 n = 0;
1334 for (int i = 0; i < nfix; i++) if (fmask[i] & mask) list[n++] = i;
1335 }
1336
1337 /* ----------------------------------------------------------------------
1338 create list of fix indices for for pre_exchange fixes
1339 have contacthistory fixes always come first so it can copy the data
1340 ------------------------------------------------------------------------- */
1341
list_init_pre_exchange(int mask,int & n,int * & list)1342 void Modify::list_init_pre_exchange(int mask, int &n, int *&list)
1343 {
1344 delete [] list;
1345
1346 n = 0;
1347 for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++;
1348 list = new int[n];
1349
1350 n = 0;
1351
1352 for (int i = 0; i < nfix; i++) if (fmask[i] & mask)
1353 {
1354 if(0 == strncmp(fix[i]->style,"contacthistory",14))
1355 //if(0 == strcmp(fix[i]->style,"contacthistory"))
1356 list[n++] = i;
1357 }
1358
1359 for (int i = 0; i < nfix; i++)
1360 {
1361 if(0 == strncmp(fix[i]->style,"contacthistory",14))
1362 //if(0 == strcmp(fix[i]->style,"contacthistory"))
1363 continue;
1364
1365 if (fmask[i] & mask) list[n++] = i;
1366 }
1367 }
1368
1369 /* ----------------------------------------------------------------------
1370 create list of fix indices for end_of_step fixes
1371 also create end_of_step_every[]
1372 ------------------------------------------------------------------------- */
1373
list_init_end_of_step(int mask,int & n,int * & list)1374 void Modify::list_init_end_of_step(int mask, int &n, int *&list)
1375 {
1376 delete [] list;
1377 delete [] end_of_step_every;
1378
1379 n = 0;
1380 for (int i = 0; i < nfix; i++) if (fmask[i] & mask) n++;
1381 list = new int[n];
1382 end_of_step_every = new int[n];
1383
1384 n = 0;
1385 for (int i = 0; i < nfix; i++)
1386 if (fmask[i] & mask) {
1387 list[n] = i;
1388 end_of_step_every[n++] = fix[i]->nevery;
1389 }
1390 }
1391
1392 /* ----------------------------------------------------------------------
1393 create list of fix indices for thermo energy fixes
1394 only added to list if fix has THERMO_ENERGY mask
1395 and its thermo_energy flag was set via fix_modify
1396 ------------------------------------------------------------------------- */
1397
list_init_thermo_energy(int mask,int & n,int * & list)1398 void Modify::list_init_thermo_energy(int mask, int &n, int *&list)
1399 {
1400 delete [] list;
1401
1402 n = 0;
1403 for (int i = 0; i < nfix; i++)
1404 if (fmask[i] & mask && fix[i]->thermo_energy) n++;
1405 list = new int[n];
1406
1407 n = 0;
1408 for (int i = 0; i < nfix; i++)
1409 if (fmask[i] & mask && fix[i]->thermo_energy) list[n++] = i;
1410 }
1411
1412 /* ----------------------------------------------------------------------
1413 create list of compute indices for computes which store invocation times
1414 ------------------------------------------------------------------------- */
1415
list_init_compute()1416 void Modify::list_init_compute()
1417 {
1418 delete [] list_timeflag;
1419
1420 n_timeflag = 0;
1421 for (int i = 0; i < ncompute; i++)
1422 if (compute[i]->timeflag) n_timeflag++;
1423 list_timeflag = new int[n_timeflag];
1424
1425 n_timeflag = 0;
1426 for (int i = 0; i < ncompute; i++)
1427 if (compute[i]->timeflag) list_timeflag[n_timeflag++] = i;
1428 }
1429
1430 /* ----------------------------------------------------------------------
1431 return # of bytes of allocated memory from all fixes
1432 ------------------------------------------------------------------------- */
1433
memory_usage()1434 bigint Modify::memory_usage()
1435 {
1436 bigint bytes = 0;
1437 for (int i = 0; i < nfix; i++)
1438 bytes += static_cast<bigint> (fix[i]->memory_usage());
1439 for (int i = 0; i < ncompute; i++)
1440 bytes += static_cast<bigint> (compute[i]->memory_usage());
1441 return bytes;
1442 }
1443
1444 /* ======================================================================
1445 helper functions by Richard Berger (JKU)
1446 ========================================================================= */
1447
1448 /* ----------------------------------------------------------------------
1449 calls a member method on all fixes
1450 ------------------------------------------------------------------------- */
1451
call_method_on_fixes(FixMethod method)1452 void Modify::call_method_on_fixes(FixMethod method) {
1453 if(timing) {
1454 for (int i = 0; i < nfix; i++) {
1455 fix[i]->begin_time_recording();
1456 (fix[i]->*method)();
1457 fix[i]->end_time_recording();
1458 }
1459 }
1460 else
1461 {
1462 for (int i = 0; i < nfix; i++) {
1463 (fix[i]->*method)();
1464 }
1465 }
1466 }
1467
1468 /* ----------------------------------------------------------------------
1469 calls a member method on all fixes in the specified list
1470 ------------------------------------------------------------------------- */
1471
call_method_on_fixes(FixMethod method,int * & ilist,int & inum)1472 void Modify::call_method_on_fixes(FixMethod method, int *& ilist, int & inum) {
1473 if(timing) {
1474 for (int i = 0; i < inum; i++) {
1475 const int ifix = ilist[i];
1476 fix[ifix]->begin_time_recording();
1477 (fix[ifix]->*method)();
1478 fix[ifix]->end_time_recording();
1479 }
1480 }
1481 else
1482 {
1483 for (int i = 0; i < inum; i++) {
1484 (fix[ilist[i]]->*method)();
1485 }
1486 }
1487 }
1488
1489 /* ----------------------------------------------------------------------
1490 calls a member method with vflag parameter on all fixes
1491 ------------------------------------------------------------------------- */
1492
call_method_on_fixes(FixMethodWithVFlag method,int vflag)1493 void Modify::call_method_on_fixes(FixMethodWithVFlag method, int vflag) {
1494 if(timing) {
1495 for (int i = 0; i < nfix; i++) {
1496 fix[i]->begin_time_recording();
1497 (fix[i]->*method)(vflag);
1498 fix[i]->end_time_recording();
1499 }
1500 }
1501 else
1502 {
1503 for (int i = 0; i < nfix; i++) {
1504 (fix[i]->*method)(vflag);
1505 }
1506 }
1507 }
1508
1509 /* ----------------------------------------------------------------------
1510 calls a member method with vflag parameter on all fixes in the
1511 specified list
1512 ------------------------------------------------------------------------- */
1513
call_method_on_fixes(FixMethodWithVFlag method,int vflag,int * & ilist,int & inum)1514 void Modify::call_method_on_fixes(FixMethodWithVFlag method, int vflag, int *& ilist, int & inum) {
1515 if(timing) {
1516 for (int i = 0; i < inum; i++) {
1517 const int ifix = ilist[i];
1518 fix[ifix]->begin_time_recording();
1519 (fix[ifix]->*method)(vflag);
1520 fix[ifix]->end_time_recording();
1521 }
1522 }
1523 else
1524 {
1525 for (int i = 0; i < inum; i++) {
1526 (fix[ilist[i]]->*method)(vflag);
1527 }
1528 }
1529 }
1530
1531 /* ----------------------------------------------------------------------
1532 calls a respa member method with 2 int parameters on all fixes in the
1533 specified list
1534 ------------------------------------------------------------------------- */
1535
call_respa_method_on_fixes(FixMethodRESPA2 method,int arg1,int arg2,int * & ilist,int & inum)1536 void Modify::call_respa_method_on_fixes(FixMethodRESPA2 method,
1537 int arg1, int arg2, int *& ilist, int & inum) {
1538 if(timing) {
1539 for (int i = 0; i < inum; i++) {
1540 const int ifix = ilist[i];
1541 fix[ifix]->begin_time_recording();
1542 (fix[ifix]->*method)(arg1, arg2);
1543 fix[ifix]->end_time_recording();
1544 }
1545 }
1546 else
1547 {
1548 for (int i = 0; i < inum; i++) {
1549 (fix[ilist[i]]->*method)(arg1, arg2);
1550 }
1551 }
1552 }
1553
1554 /* ----------------------------------------------------------------------
1555 calls a respa member method with 3 int parameters on all fixes in the
1556 specified list
1557 ------------------------------------------------------------------------- */
1558
call_respa_method_on_fixes(FixMethodRESPA3 method,int arg1,int arg2,int arg3,int * & ilist,int & inum)1559 void Modify::call_respa_method_on_fixes(FixMethodRESPA3 method, int arg1,
1560 int arg2, int arg3, int *& ilist, int & inum) {
1561 if(timing) {
1562 for (int i = 0; i < inum; i++) {
1563 const int ifix = ilist[i];
1564 fix[ifix]->begin_time_recording();
1565 (fix[ifix]->*method)(arg1, arg2, arg3);
1566 fix[ifix]->end_time_recording();
1567 }
1568 }
1569 else
1570 {
1571 for (int i = 0; i < inum; i++) {
1572 (fix[ilist[i]]->*method)(arg1, arg2, arg3);
1573 }
1574 }
1575 }
1576
1577 /* ----------------------------------------------------------------------
1578 Updates all computes that requested it at the end of a run
1579 ------------------------------------------------------------------------- */
1580
update_computes_on_run_end()1581 void Modify::update_computes_on_run_end()
1582 {
1583 for (int i = 0; i < ncompute; i++)
1584 {
1585 if (compute[i]->update_on_run_end())
1586 {
1587 if (compute[i]->scalar_flag)
1588 {
1589 if (!(compute[i]->invoked_flag & INVOKED_SCALAR))
1590 {
1591 compute[i]->compute_scalar();
1592 compute[i]->invoked_flag |= INVOKED_SCALAR;
1593 }
1594 }
1595 if (compute[i]->vector_flag)
1596 {
1597 if (!(compute[i]->invoked_flag & INVOKED_VECTOR))
1598 {
1599 compute[i]->compute_vector();
1600 compute[i]->invoked_flag |= INVOKED_VECTOR;
1601 }
1602 }
1603 if (compute[i]->array_flag)
1604 {
1605 if (!(compute[i]->invoked_flag & INVOKED_ARRAY))
1606 {
1607 compute[i]->compute_array();
1608 compute[i]->invoked_flag |= INVOKED_ARRAY;
1609 }
1610 }
1611 }
1612 }
1613 }
1614