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
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 /* ----------------------------------------------------------------------
47 Contributing authors: Mark Stevens (SNL), Paul Crozier (SNL)
48 ------------------------------------------------------------------------- */
49
50 #include <stdlib.h>
51 #include <string.h>
52 #include "respa.h"
53 #include "neighbor.h"
54 #include "atom.h"
55 #include "domain.h"
56 #include "comm.h"
57 #include "force.h"
58 #include "pair.h"
59 #include "bond.h"
60 #include "angle.h"
61 #include "dihedral.h"
62 #include "improper.h"
63 #include "kspace.h"
64 #include "output.h"
65 #include "update.h"
66 #include "modify.h"
67 #include "compute.h"
68 #include "fix_respa.h"
69 #include "timer.h"
70 #include "memory.h"
71 #include "error.h"
72 #include "signal_handling.h"
73
74 using namespace LAMMPS_NS;
75
76 /* ---------------------------------------------------------------------- */
77
Respa(LAMMPS * lmp,int narg,char ** arg)78 Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
79 {
80 if (narg < 1) error->all(FLERR,"Illegal run_style respa command");
81
82 nlevels = force->inumeric(FLERR,arg[0]);
83 if (nlevels < 1) error->all(FLERR,"Respa levels must be >= 1");
84
85 if (narg < nlevels) error->all(FLERR,"Illegal run_style respa command");
86 loop = new int[nlevels];
87 for (int iarg = 1; iarg < nlevels; iarg++) {
88 loop[iarg-1] = force->inumeric(FLERR,arg[iarg]);
89 if (loop[iarg-1] <= 0) error->all(FLERR,"Illegal run_style respa command");
90 }
91 loop[nlevels-1] = 1;
92
93 // set level at which each force is computed
94 // argument settings override defaults
95
96 level_bond = level_angle = level_dihedral = level_improper = -1;
97 level_pair = level_kspace = -1;
98 level_inner = level_middle = level_outer = -1;
99
100 int iarg = nlevels;
101 while (iarg < narg) {
102 if (strcmp(arg[iarg],"bond") == 0) {
103 if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
104 level_bond = force->inumeric(FLERR,arg[iarg+1]) - 1;
105 iarg += 2;
106 } else if (strcmp(arg[iarg],"angle") == 0) {
107 if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
108 level_angle = force->inumeric(FLERR,arg[iarg+1]) - 1;
109 iarg += 2;
110 } else if (strcmp(arg[iarg],"dihedral") == 0) {
111 if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
112 level_dihedral = force->inumeric(FLERR,arg[iarg+1]) - 1;
113 iarg += 2;
114 } else if (strcmp(arg[iarg],"improper") == 0) {
115 if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
116 level_improper = force->inumeric(FLERR,arg[iarg+1]) - 1;
117 iarg += 2;
118 } else if (strcmp(arg[iarg],"pair") == 0) {
119 if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
120 level_pair = force->inumeric(FLERR,arg[iarg+1]) - 1;
121 iarg += 2;
122 } else if (strcmp(arg[iarg],"inner") == 0) {
123 if (iarg+4 > narg) error->all(FLERR,"Illegal run_style respa command");
124 level_inner = force->inumeric(FLERR,arg[iarg+1]) - 1;
125 cutoff[0] = force->numeric(FLERR,arg[iarg+2]);
126 cutoff[1] = force->numeric(FLERR,arg[iarg+3]);
127 iarg += 4;
128 } else if (strcmp(arg[iarg],"middle") == 0) {
129 if (iarg+4 > narg) error->all(FLERR,"Illegal run_style respa command");
130 level_middle = force->inumeric(FLERR,arg[iarg+1]) - 1;
131 cutoff[2] = force->numeric(FLERR,arg[iarg+2]);
132 cutoff[3] = force->numeric(FLERR,arg[iarg+3]);
133 iarg += 4;
134 } else if (strcmp(arg[iarg],"outer") == 0) {
135 if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
136 level_outer = force->inumeric(FLERR,arg[iarg+1]) - 1;
137 iarg += 2;
138 } else if (strcmp(arg[iarg],"kspace") == 0) {
139 if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
140 level_kspace = force->inumeric(FLERR,arg[iarg+1]) - 1;
141 iarg += 2;
142 } else error->all(FLERR,"Illegal run_style respa command");
143 }
144
145 // cannot specify both pair and inner/middle/outer
146
147 if (level_pair >= 0 &&
148 (level_inner >= 0 || level_middle >= 0 || level_outer >= 0))
149 error->all(FLERR,"Cannot set both respa pair and inner/middle/outer");
150
151 // if either inner and outer is specified, then both must be
152
153 if ((level_inner >= 0 && level_outer == -1) ||
154 (level_outer >= 0 && level_inner == -1))
155 error->all(FLERR,"Must set both respa inner and outer");
156
157 // middle cannot be set without inner/outer
158
159 if (level_middle >= 0 && level_inner == -1)
160 error->all(FLERR,"Cannot set respa middle without inner/outer");
161
162 // set defaults if user did not specify level
163 // bond to innermost level
164 // angle same as bond, dihedral same as angle, improper same as dihedral
165 // pair to outermost level if no inner/middle/outer
166 // inner/middle/outer have no defaults
167 // kspace same as pair or outer
168
169 if (level_bond == -1) level_bond = 0;
170 if (level_angle == -1) level_angle = level_bond;
171 if (level_dihedral == -1) level_dihedral = level_angle;
172 if (level_improper == -1) level_improper = level_dihedral;
173 if (level_pair == -1 && level_inner == -1) level_pair = nlevels-1;
174 if (level_kspace == -1 && level_pair >= 0) level_kspace = level_pair;
175 if (level_kspace == -1 && level_pair == -1) level_kspace = level_outer;
176
177 // print respa levels
178
179 if (comm->me == 0) {
180 if (screen) {
181 fprintf(screen,"Respa levels:\n");
182 for (int i = 0; i < nlevels; i++) {
183 fprintf(screen," %d =",i+1);
184 if (level_bond == i) fprintf(screen," bond");
185 if (level_angle == i) fprintf(screen," angle");
186 if (level_dihedral == i) fprintf(screen," dihedral");
187 if (level_improper == i) fprintf(screen," improper");
188 if (level_pair == i) fprintf(screen," pair");
189 if (level_inner == i) fprintf(screen," pair-inner");
190 if (level_middle == i) fprintf(screen," pair-middle");
191 if (level_outer == i) fprintf(screen," pair-outer");
192 if (level_kspace == i) fprintf(screen," kspace");
193 fprintf(screen,"\n");
194 }
195 }
196 if (logfile) {
197 fprintf(logfile,"Respa levels:\n");
198 for (int i = 0; i < nlevels; i++) {
199 fprintf(logfile," %d =",i+1);
200 if (level_bond == i) fprintf(logfile," bond");
201 if (level_angle == i) fprintf(logfile," angle");
202 if (level_dihedral == i) fprintf(logfile," dihedral");
203 if (level_improper == i) fprintf(logfile," improper");
204 if (level_pair == i) fprintf(logfile," pair");
205 if (level_inner == i) fprintf(logfile," pair-inner");
206 if (level_middle == i) fprintf(logfile," pair-middle");
207 if (level_outer == i) fprintf(logfile," pair-outer");
208 if (level_kspace == i) fprintf(logfile," kspace");
209 fprintf(logfile,"\n");
210 }
211 }
212 }
213
214 // check that levels are in correct order
215
216 if (level_angle < level_bond || level_dihedral < level_angle ||
217 level_improper < level_dihedral)
218 error->all(FLERR,"Invalid order of forces within respa levels");
219 if (level_pair >= 0) {
220 if (level_pair < level_improper || level_kspace < level_pair)
221 error->all(FLERR,"Invalid order of forces within respa levels");
222 }
223 if (level_pair == -1 && level_middle == -1) {
224 if (level_inner < level_improper || level_outer < level_inner ||
225 level_kspace < level_outer)
226 error->all(FLERR,"Invalid order of forces within respa levels");
227 }
228 if (level_pair == -1 && level_middle >= 0) {
229 if (level_inner < level_improper || level_middle < level_inner ||
230 level_outer < level_inner || level_kspace < level_outer)
231 error->all(FLERR,"Invalid order of forces within respa levels");
232 }
233
234 // warn if any levels are devoid of forces
235
236 int flag = 0;
237 for (int i = 0; i < nlevels; i++)
238 if (level_bond != i && level_angle != i && level_dihedral != i &&
239 level_improper != i && level_pair != i && level_inner != i &&
240 level_middle != i && level_outer != i && level_kspace != i) flag = 1;
241 if (flag && comm->me == 0)
242 error->warning(FLERR,"One or more respa levels compute no forces");
243
244 // check cutoff consistency if inner/middle/outer are enabled
245
246 if (level_inner >= 0 && cutoff[1] < cutoff[0])
247 error->all(FLERR,"Respa inner cutoffs are invalid");
248 if (level_middle >= 0 && (cutoff[3] < cutoff[2] || cutoff[2] < cutoff[1]))
249 error->all(FLERR,"Respa middle cutoffs are invalid");
250
251 // set outer pair of cutoffs to inner pair if middle is not enabled
252
253 if (level_inner >= 0 && level_middle < 0) {
254 cutoff[2] = cutoff[0];
255 cutoff[3] = cutoff[1];
256 }
257
258 // allocate other needed arrays
259
260 newton = new int[nlevels];
261 step = new double[nlevels];
262 }
263
264 /* ---------------------------------------------------------------------- */
265
~Respa()266 Respa::~Respa()
267 {
268 delete [] loop;
269 delete [] newton;
270 delete [] step;
271 }
272
273 /* ----------------------------------------------------------------------
274 initialization before run
275 ------------------------------------------------------------------------- */
276
init()277 void Respa::init()
278 {
279 Integrate::init();
280
281 // warn if no fixes
282
283 if (modify->nfix == 0 && comm->me == 0)
284 error->warning(FLERR,"No fixes defined, atoms won't move");
285
286 // create fix needed for storing atom-based respa level forces
287 // will delete it at end of run
288
289 char **fixarg = new char*[4];
290 fixarg[0] = (char *) "RESPA";
291 fixarg[1] = (char *) "all";
292 fixarg[2] = (char *) "RESPA";
293 fixarg[3] = new char[8];
294 sprintf(fixarg[3],"%d",nlevels);
295 modify->add_fix(4,fixarg);
296 delete [] fixarg[3];
297 delete [] fixarg;
298 fix_respa = (FixRespa *) modify->fix[modify->nfix-1];
299
300 // insure respa inner/middle/outer is using Pair class that supports it
301
302 if (level_inner >= 0)
303 if (force->pair && force->pair->respa_enable == 0)
304 error->all(FLERR,"Pair style does not support rRESPA inner/middle/outer");
305
306 // virial_style = 1 (explicit) since never computed implicitly like Verlet
307
308 virial_style = 1;
309
310 // setup lists of computes for global and per-atom PE and pressure
311
312 ev_setup();
313
314 // detect if fix omp is present and will clear force arrays
315
316 int ifix = modify->find_fix("package_omp");
317 if (ifix >= 0) external_force_clear = 1;
318
319 // set flags for what arrays to clear in force_clear()
320 // need to clear additionals arrays if they exist
321
322 torqueflag = 0;
323 if (atom->torque_flag) torqueflag = 1;
324 erforceflag = 0;
325 if (atom->erforce_flag) erforceflag = 1;
326 e_flag = 0;
327 if (atom->e_flag) e_flag = 1;
328 rho_flag = 0;
329 if (atom->rho_flag) rho_flag = 1;
330
331 // step[] = timestep for each level
332
333 step[nlevels-1] = update->dt;
334 for (int ilevel = nlevels-2; ilevel >= 0; ilevel--)
335 step[ilevel] = step[ilevel+1]/loop[ilevel];
336
337 // set newton flag for each level
338
339 for (int ilevel = 0; ilevel < nlevels; ilevel++) {
340 newton[ilevel] = 0;
341 if (force->newton_bond) {
342 if (level_bond == ilevel || level_angle == ilevel ||
343 level_dihedral == ilevel || level_improper == ilevel)
344 newton[ilevel] = 1;
345 }
346 if (force->newton_pair) {
347 if (level_pair == ilevel || level_inner == ilevel ||
348 level_middle == ilevel || level_outer == ilevel)
349 newton[ilevel] = 1;
350 }
351 }
352
353 // orthogonal vs triclinic simulation box
354
355 triclinic = domain->triclinic;
356 }
357
358 /* ----------------------------------------------------------------------
359 setup before run
360 ------------------------------------------------------------------------- */
361
setup()362 void Respa::setup()
363 {
364 if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n");
365
366 update->setupflag = 1;
367
368 // setup domain, communication and neighboring
369 // acquire ghosts
370 // build neighbor lists
371
372 atom->setup();
373 modify->setup_pre_exchange();
374 if (triclinic) domain->x2lamda(atom->nlocal);
375 domain->pbc();
376 domain->reset_box();
377 comm->setup();
378 if (neighbor->style) neighbor->setup_bins();
379 comm->exchange();
380 if (atom->sortfreq > 0) atom->sort();
381 comm->borders();
382 if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
383 domain->image_check();
384 domain->box_too_small_check();
385 modify->setup_pre_neighbor();
386 neighbor->build();
387 neighbor->ncalls = 0;
388
389 // compute all forces
390
391 ev_set(update->ntimestep);
392
393 for (int ilevel = 0; ilevel < nlevels; ilevel++) {
394 force_clear(newton[ilevel]);
395 modify->setup_pre_force_respa(vflag,ilevel);
396 if (level_pair == ilevel && pair_compute_flag)
397 force->pair->compute(eflag,vflag);
398 if (level_inner == ilevel && pair_compute_flag)
399 force->pair->compute_inner();
400 if (level_middle == ilevel && pair_compute_flag)
401 force->pair->compute_middle();
402 if (level_outer == ilevel && pair_compute_flag)
403 force->pair->compute_outer(eflag,vflag);
404 if (level_bond == ilevel && force->bond)
405 force->bond->compute(eflag,vflag);
406 if (level_angle == ilevel && force->angle)
407 force->angle->compute(eflag,vflag);
408 if (level_dihedral == ilevel && force->dihedral)
409 force->dihedral->compute(eflag,vflag);
410 if (level_improper == ilevel && force->improper)
411 force->improper->compute(eflag,vflag);
412 if (level_kspace == ilevel && force->kspace) {
413 force->kspace->setup();
414 if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
415 }
416 if (newton[ilevel]) comm->reverse_comm();
417 copy_f_flevel(ilevel);
418 }
419
420 modify->setup(vflag);
421 sum_flevel_f();
422 output->setup();
423 update->setupflag = 0;
424 }
425
426 /* ----------------------------------------------------------------------
427 setup without output
428 flag = 0 = just force calculation
429 flag = 1 = reneighbor and force calculation
430 ------------------------------------------------------------------------- */
431
setup_minimal(int flag)432 void Respa::setup_minimal(int flag)
433 {
434 update->setupflag = 1;
435
436 // setup domain, communication and neighboring
437 // acquire ghosts
438 // build neighbor lists
439
440 if (flag) {
441 modify->setup_pre_exchange();
442 if (triclinic) domain->x2lamda(atom->nlocal);
443 domain->pbc();
444 domain->reset_box();
445 comm->setup();
446 if (neighbor->style) neighbor->setup_bins();
447 comm->exchange();
448 comm->borders();
449 if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
450 domain->image_check();
451 domain->box_too_small_check();
452 modify->setup_pre_neighbor();
453 neighbor->build();
454 neighbor->ncalls = 0;
455 }
456
457 // compute all forces
458
459 ev_set(update->ntimestep);
460
461 for (int ilevel = 0; ilevel < nlevels; ilevel++) {
462 force_clear(newton[ilevel]);
463 modify->setup_pre_force_respa(vflag,ilevel);
464 if (level_pair == ilevel && pair_compute_flag)
465 force->pair->compute(eflag,vflag);
466 if (level_inner == ilevel && pair_compute_flag)
467 force->pair->compute_inner();
468 if (level_middle == ilevel && pair_compute_flag)
469 force->pair->compute_middle();
470 if (level_outer == ilevel && pair_compute_flag)
471 force->pair->compute_outer(eflag,vflag);
472 if (level_bond == ilevel && force->bond)
473 force->bond->compute(eflag,vflag);
474 if (level_angle == ilevel && force->angle)
475 force->angle->compute(eflag,vflag);
476 if (level_dihedral == ilevel && force->dihedral)
477 force->dihedral->compute(eflag,vflag);
478 if (level_improper == ilevel && force->improper)
479 force->improper->compute(eflag,vflag);
480 if (level_kspace == ilevel && force->kspace) {
481 force->kspace->setup();
482 if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
483 }
484 if (newton[ilevel]) comm->reverse_comm();
485 copy_f_flevel(ilevel);
486 }
487
488 modify->setup(vflag);
489 sum_flevel_f();
490 update->setupflag = 0;
491 }
492
493 /* ----------------------------------------------------------------------
494 run for N steps
495 ------------------------------------------------------------------------- */
496
run(int n)497 void Respa::run(int n)
498 {
499 bigint ntimestep;
500
501 for (int i = 0; i < n; i++) {
502
503 ntimestep = ++update->ntimestep;
504 ev_set(ntimestep);
505
506 recurse(nlevels-1);
507
508 if (modify->n_end_of_step) modify->end_of_step();
509
510 if (ntimestep == output->next) {
511 timer->stamp();
512 sum_flevel_f();
513 output->write(update->ntimestep);
514 timer->stamp(TIME_OUTPUT);
515 }
516
517 if (SignalHandler::request_quit && !SignalHandler::request_write_restart)
518 break;
519 }
520 }
521
522 /* ----------------------------------------------------------------------
523 delete rRESPA fix at end of run, so its atom arrays won't persist
524 ------------------------------------------------------------------------- */
525
cleanup()526 void Respa::cleanup()
527 {
528 modify->post_run();
529 modify->delete_fix("RESPA");
530 domain->box_too_small_check();
531 update->update_time();
532 }
533
534 /* ---------------------------------------------------------------------- */
535
reset_dt()536 void Respa::reset_dt()
537 {
538 step[nlevels-1] = update->dt;
539 for (int ilevel = nlevels-2; ilevel >= 0; ilevel--)
540 step[ilevel] = step[ilevel+1]/loop[ilevel];
541 }
542
543 /* ---------------------------------------------------------------------- */
544
recurse(int ilevel)545 void Respa::recurse(int ilevel)
546 {
547 copy_flevel_f(ilevel);
548
549 for (int iloop = 0; iloop < loop[ilevel]; iloop++) {
550
551 modify->initial_integrate_respa(vflag,ilevel,iloop);
552 if (modify->n_post_integrate_respa)
553 modify->post_integrate_respa(ilevel,iloop);
554
555 if (ilevel) recurse(ilevel-1);
556
557 // at outermost level, check on rebuilding neighbor list
558 // at innermost level, communicate
559 // at middle levels, do nothing
560
561 if (ilevel == nlevels-1) {
562 int nflag = neighbor->decide();
563 if (nflag) {
564 if (modify->n_pre_exchange) modify->pre_exchange();
565 if (triclinic) domain->x2lamda(atom->nlocal);
566 domain->pbc();
567 if (domain->box_change) {
568 domain->reset_box();
569 comm->setup();
570 if (neighbor->style) neighbor->setup_bins();
571 }
572 timer->stamp();
573 comm->exchange();
574 if (atom->sortfreq > 0 &&
575 update->ntimestep >= atom->nextsort) atom->sort();
576 comm->borders();
577 if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
578 timer->stamp(TIME_COMM);
579 if (modify->n_pre_neighbor) modify->pre_neighbor();
580 neighbor->build();
581 timer->stamp(TIME_NEIGHBOR);
582 } else if (ilevel == 0) {
583 timer->stamp();
584 comm->forward_comm();
585 timer->stamp(TIME_COMM);
586 }
587
588 } else if (ilevel == 0) {
589 timer->stamp();
590 comm->forward_comm();
591 timer->stamp(TIME_COMM);
592 }
593
594 // force computations
595 // important that ordering is same as Verlet
596 // so that any order dependencies are the same
597 // when potentials are invoked at same level
598
599 force_clear(newton[ilevel]);
600 if (modify->n_pre_force_respa)
601 modify->pre_force_respa(vflag,ilevel,iloop);
602
603 timer->stamp();
604 if (level_pair == ilevel && pair_compute_flag) {
605 force->pair->compute(eflag,vflag);
606 timer->stamp(TIME_PAIR);
607 }
608 if (level_inner == ilevel && pair_compute_flag) {
609 force->pair->compute_inner();
610 timer->stamp(TIME_PAIR);
611 }
612 if (level_middle == ilevel && pair_compute_flag) {
613 force->pair->compute_middle();
614 timer->stamp(TIME_PAIR);
615 }
616 if (level_outer == ilevel && pair_compute_flag) {
617 force->pair->compute_outer(eflag,vflag);
618 timer->stamp(TIME_PAIR);
619 }
620 if (level_bond == ilevel && force->bond) {
621 force->bond->compute(eflag,vflag);
622 timer->stamp(TIME_BOND);
623 }
624 if (level_angle == ilevel && force->angle) {
625 force->angle->compute(eflag,vflag);
626 timer->stamp(TIME_BOND);
627 }
628 if (level_dihedral == ilevel && force->dihedral) {
629 force->dihedral->compute(eflag,vflag);
630 timer->stamp(TIME_BOND);
631 }
632 if (level_improper == ilevel && force->improper) {
633 force->improper->compute(eflag,vflag);
634 timer->stamp(TIME_BOND);
635 }
636 if (level_kspace == ilevel && kspace_compute_flag) {
637 force->kspace->compute(eflag,vflag);
638 timer->stamp(TIME_KSPACE);
639 }
640
641 if (newton[ilevel]) {
642 comm->reverse_comm();
643 timer->stamp(TIME_COMM);
644 }
645
646 if (modify->n_post_force_respa)
647 modify->post_force_respa(vflag,ilevel,iloop);
648 modify->final_integrate_respa(ilevel,iloop);
649 }
650
651 copy_f_flevel(ilevel);
652 }
653
654 /* ----------------------------------------------------------------------
655 clear force on own & ghost atoms
656 ------------------------------------------------------------------------- */
657
force_clear(int newtonflag)658 void Respa::force_clear(int newtonflag)
659 {
660 if (external_force_clear) return;
661
662 // clear global force array
663 // nall includes ghosts only if newton flag is set
664
665 int nall;
666 if (newtonflag) nall = atom->nlocal + atom->nghost;
667 else nall = atom->nlocal;
668
669 size_t nbytes = sizeof(double) * nall;
670
671 if (nbytes > 0 ) {
672 memset(&(atom->f[0][0]),0,3*nbytes);
673 if (torqueflag) memset(&(atom->torque[0][0]),0,3*nbytes);
674 if (erforceflag) memset(&(atom->erforce[0]), 0, nbytes);
675 if (e_flag) memset(&(atom->de[0]), 0, nbytes);
676 if (rho_flag) memset(&(atom->drho[0]), 0, nbytes);
677 }
678 }
679
680 /* ----------------------------------------------------------------------
681 copy force components from atom->f to FixRespa->f_level
682 ------------------------------------------------------------------------- */
683
copy_f_flevel(int ilevel)684 void Respa::copy_f_flevel(int ilevel)
685 {
686 double ***f_level = fix_respa->f_level;
687 double **f = atom->f;
688 int n = atom->nlocal;
689
690 for (int i = 0; i < n; i++) {
691 f_level[i][ilevel][0] = f[i][0];
692 f_level[i][ilevel][1] = f[i][1];
693 f_level[i][ilevel][2] = f[i][2];
694 }
695 }
696
697 /* ----------------------------------------------------------------------
698 copy force components from FixRespa->f_level to atom->f
699 ------------------------------------------------------------------------- */
700
copy_flevel_f(int ilevel)701 void Respa::copy_flevel_f(int ilevel)
702 {
703 double ***f_level = fix_respa->f_level;
704 double **f = atom->f;
705 int n = atom->nlocal;
706
707 for (int i = 0; i < n; i++) {
708 f[i][0] = f_level[i][ilevel][0];
709 f[i][1] = f_level[i][ilevel][1];
710 f[i][2] = f_level[i][ilevel][2];
711 }
712 }
713
714 /* ----------------------------------------------------------------------
715 sum all force components from FixRespa->f_level to create full atom->f
716 ------------------------------------------------------------------------- */
717
sum_flevel_f()718 void Respa::sum_flevel_f()
719 {
720 copy_flevel_f(0);
721
722 double ***f_level = fix_respa->f_level;
723 double **f = atom->f;
724 int n = atom->nlocal;
725
726 for (int ilevel = 1; ilevel < nlevels; ilevel++) {
727 for (int i = 0; i < n; i++) {
728 f[i][0] += f_level[i][ilevel][0];
729 f[i][1] += f_level[i][ilevel][1];
730 f[i][2] += f_level[i][ilevel][2];
731 }
732 }
733 }
734