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 #include <cmath>
47 #include <stdlib.h>
48 #include <string.h>
49 #include "ctype.h"
50 #include "pair_hybrid.h"
51 #include "atom.h"
52 #include "force.h"
53 #include "pair.h"
54 #include "neighbor.h"
55 #include "neigh_request.h"
56 #include "update.h"
57 #include "comm.h"
58 #include "memory.h"
59 #include "error.h"
60
61 using namespace LAMMPS_NS;
62
63 /* ---------------------------------------------------------------------- */
64
PairHybrid(LAMMPS * lmp)65 PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp)
66 {
67 nstyles = 0;
68 styles = NULL;
69 keywords = NULL;
70 multiple = NULL;
71
72 outerflag = 0;
73 }
74
75 /* ---------------------------------------------------------------------- */
76
~PairHybrid()77 PairHybrid::~PairHybrid()
78 {
79 if (nstyles) {
80 for (int m = 0; m < nstyles; m++) delete styles[m];
81 for (int m = 0; m < nstyles; m++) delete [] keywords[m];
82 }
83 delete [] styles;
84 delete [] keywords;
85 delete [] multiple;
86
87 delete [] svector;
88
89 if (allocated) {
90 memory->destroy(setflag);
91 memory->destroy(cutsq);
92 memory->destroy(cutghost);
93 memory->destroy(nmap);
94 memory->destroy(map);
95 }
96 }
97
98 /* ----------------------------------------------------------------------
99 call each sub-style's compute() or compute_outer() function
100 accumulate sub-style global/peratom energy/virial in hybrid
101 for global vflag = 1:
102 each sub-style computes own virial[6]
103 sum sub-style virial[6] to hybrid's virial[6]
104 for global vflag = 2:
105 call sub-style with adjusted vflag to prevent it calling
106 virial_fdotr_compute()
107 hybrid calls virial_fdotr_compute() on final accumulated f
108 ------------------------------------------------------------------------- */
109
compute(int eflag,int vflag)110 void PairHybrid::compute(int eflag, int vflag)
111 {
112 int i,j,m,n;
113
114 // if no_virial_fdotr_compute is set and global component of
115 // incoming vflag = 2, then
116 // reset vflag as if global component were 1
117 // necessary since one or more sub-styles cannot compute virial as F dot r
118
119 if (no_virial_fdotr_compute && vflag % 4 == 2) vflag = 1 + vflag/4 * 4;
120
121 if (eflag || vflag) ev_setup(eflag,vflag);
122 else evflag = vflag_fdotr = eflag_global = vflag_global =
123 eflag_atom = vflag_atom = 0;
124
125 // check if global component of incoming vflag = 2
126 // if so, reset vflag passed to substyle as if it were 0
127 // necessary so substyle will not invoke virial_fdotr_compute()
128
129 int vflag_substyle;
130 if (vflag % 4 == 2) vflag_substyle = vflag/4 * 4;
131 else vflag_substyle = vflag;
132
133 for (m = 0; m < nstyles; m++) {
134
135 // invoke compute() unless
136 // outerflag is set and sub-style has a compute_outer() method
137
138 if (outerflag && styles[m]->respa_enable)
139 styles[m]->compute_outer(eflag,vflag_substyle);
140 else styles[m]->compute(eflag,vflag_substyle);
141
142 if (eflag_global) {
143 eng_vdwl += styles[m]->eng_vdwl;
144 eng_coul += styles[m]->eng_coul;
145 }
146 if (vflag_global) {
147 for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n];
148 }
149 if (eflag_atom) {
150 n = atom->nlocal;
151 if (force->newton_pair) n += atom->nghost;
152 double *eatom_substyle = styles[m]->eatom;
153 for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i];
154 }
155 if (vflag_atom) {
156 n = atom->nlocal;
157 if (force->newton_pair) n += atom->nghost;
158 double **vatom_substyle = styles[m]->vatom;
159 for (i = 0; i < n; i++)
160 for (j = 0; j < 6; j++)
161 vatom[i][j] += vatom_substyle[i][j];
162 }
163 }
164
165 if (vflag_fdotr) virial_fdotr_compute();
166 }
167
168 /* ---------------------------------------------------------------------- */
169
compute_inner()170 void PairHybrid::compute_inner()
171 {
172 for (int m = 0; m < nstyles; m++)
173 if (styles[m]->respa_enable) styles[m]->compute_inner();
174 }
175
176 /* ---------------------------------------------------------------------- */
177
compute_middle()178 void PairHybrid::compute_middle()
179 {
180 for (int m = 0; m < nstyles; m++)
181 if (styles[m]->respa_enable) styles[m]->compute_middle();
182 }
183
184 /* ---------------------------------------------------------------------- */
185
compute_outer(int eflag,int vflag)186 void PairHybrid::compute_outer(int eflag, int vflag)
187 {
188 outerflag = 1;
189 compute(eflag,vflag);
190 outerflag = 0;
191 }
192
193 /* ----------------------------------------------------------------------
194 allocate all arrays
195 ------------------------------------------------------------------------- */
196
allocate()197 void PairHybrid::allocate()
198 {
199 allocated = 1;
200 int n = atom->ntypes;
201
202 memory->create(setflag,n+1,n+1,"pair:setflag");
203 for (int i = 1; i <= n; i++)
204 for (int j = i; j <= n; j++)
205 setflag[i][j] = 0;
206
207 memory->create(cutsq,n+1,n+1,"pair:cutsq");
208 memory->create(cutghost,n+1,n+1,"pair:cutghost");
209
210 memory->create(nmap,n+1,n+1,"pair:nmap");
211 memory->create(map,n+1,n+1,nstyles,"pair:map");
212 for (int i = 1; i <= n; i++)
213 for (int j = i; j <= n; j++)
214 nmap[i][j] = 0;
215 }
216
217 /* ----------------------------------------------------------------------
218 create one pair style for each arg in list
219 ------------------------------------------------------------------------- */
220
settings(int narg,char ** arg)221 void PairHybrid::settings(int narg, char **arg)
222 {
223
224 if (narg < 1) error->all(FLERR,"Illegal pair_style command");
225
226 // delete old lists, since cannot just change settings
227
228 if (nstyles) {
229 for (int m = 0; m < nstyles; m++) delete styles[m];
230 delete [] styles;
231 for (int m = 0; m < nstyles; m++) delete [] keywords[m];
232 delete [] keywords;
233 }
234
235 if (allocated) {
236 memory->destroy(setflag);
237 memory->destroy(cutsq);
238 memory->destroy(cutghost);
239 memory->destroy(nmap);
240 memory->destroy(map);
241 }
242 allocated = 0;
243
244 // allocate list of sub-styles as big as possibly needed if no extra args
245
246 styles = new Pair*[narg];
247 keywords = new char*[narg];
248 multiple = new int[narg];
249
250 memset(styles, 0, sizeof(Pair*)*narg);
251
252 // allocate each sub-style
253 // call settings() with set of args that are not pair style names
254 // use force->pair_map to determine which args these are
255
256 int iarg,jarg,dummy;
257
258 iarg = 0;
259 nstyles = 0;
260 while (iarg < narg) {
261 if (strcmp(arg[iarg],"hybrid") == 0)
262 error->all(FLERR,"Pair style hybrid cannot have hybrid as an argument");
263 if (strcmp(arg[iarg],"none") == 0)
264 error->all(FLERR,"Pair style hybrid cannot have none as an argument");
265
266 // find next pair style
267 int n = strlen(arg[iarg]) + 1;
268 keywords[nstyles] = new char[n];
269 strcpy(keywords[nstyles],arg[iarg]);
270 jarg = iarg + 1;
271 while (jarg < narg && !force->pair_map->count(arg[jarg])) jarg++;
272
273 int nremaining = jarg - iarg - 1;
274 char ** remaining_args = &arg[iarg+1];
275 styles[nstyles] = force->new_pair(arg[iarg],lmp->suffix,dummy);
276 styles[nstyles]->settings(nremaining,remaining_args);
277 iarg = jarg;
278 nstyles++;
279 }
280
281 // multiple[i] = 1 to M if sub-style used multiple times, else 0
282
283 for (int i = 0; i < nstyles; i++) {
284 int count = 0;
285 for (int j = 0; j < nstyles; j++) {
286 if (strcmp(keywords[j],keywords[i]) == 0) count++;
287 if (j == i) multiple[i] = count;
288 }
289 if (count == 1) multiple[i] = 0;
290 }
291
292 // set pair flags from sub-style flags
293
294 flags();
295 }
296
297 /* ----------------------------------------------------------------------
298 set top-level pair flags from sub-style flags
299 ------------------------------------------------------------------------- */
300
flags()301 void PairHybrid::flags()
302 {
303 int m;
304
305 // set comm_forward, comm_reverse, comm_reverse_off to max of any sub-style
306
307 for (m = 0; m < nstyles; m++) {
308 if (styles[m]) comm_forward = MAX(comm_forward,styles[m]->comm_forward);
309 if (styles[m]) comm_reverse = MAX(comm_reverse,styles[m]->comm_reverse);
310 if (styles[m]) comm_reverse_off = MAX(comm_reverse_off,
311 styles[m]->comm_reverse_off);
312 }
313
314 // single_enable = 1 if any sub-style is set
315 // respa_enable = 1 if any sub-style is set
316 // manybody_flag = 1 if any sub-style is set
317 // no_virial_fdotr_compute = 1 if any sub-style is set
318 // ghostneigh = 1 if any sub-style is set
319 // ewaldflag, pppmflag, msmflag, dispersionflag, tip4pflag = 1
320 // if any sub-style is set
321
322 single_enable = 0;
323 for (m = 0; m < nstyles; m++) {
324 if (styles[m]->single_enable) single_enable = 1;
325 if (styles[m]->respa_enable) respa_enable = 1;
326 if (styles[m]->manybody_flag) manybody_flag = 1;
327 if (styles[m]->no_virial_fdotr_compute) no_virial_fdotr_compute = 1;
328 if (styles[m]->ghostneigh) ghostneigh = 1;
329 if (styles[m]->ewaldflag) ewaldflag = 1;
330 if (styles[m]->pppmflag) pppmflag = 1;
331 if (styles[m]->msmflag) msmflag = 1;
332 if (styles[m]->dispersionflag) dispersionflag = 1;
333 if (styles[m]->tip4pflag) tip4pflag = 1;
334 }
335
336 // single_extra = min of all sub-style single_extra
337 // allocate svector
338
339 single_extra = styles[0]->single_extra;
340 for (m = 1; m < nstyles; m++)
341 single_extra = MIN(single_extra,styles[m]->single_extra);
342
343 if (single_extra) {
344 delete [] svector;
345 svector = new double[single_extra];
346 }
347 }
348
349 /* ----------------------------------------------------------------------
350 set coeffs for one or more type pairs
351 ------------------------------------------------------------------------- */
352
coeff(int narg,char ** arg)353 void PairHybrid::coeff(int narg, char **arg)
354 {
355 if (narg < 3) error->all(FLERR,"Incorrect args for pair coefficients");
356 if (!allocated) allocate();
357
358 int ilo,ihi,jlo,jhi;
359 force->bounds(arg[0],atom->ntypes,ilo,ihi);
360 force->bounds(arg[1],atom->ntypes,jlo,jhi);
361
362 // 3rd arg = pair sub-style name
363 // 4th arg = pair sub-style index if name used multiple times
364 // allow for "none" as valid sub-style name
365
366 int multflag = 0;
367 int m;
368
369 for (m = 0; m < nstyles; m++) {
370 multflag = 0;
371 if (strcmp(arg[2],keywords[m]) == 0) {
372 if (multiple[m]) {
373 multflag = 1;
374 if (narg < 4) error->all(FLERR,"Incorrect args for pair coefficients");
375 if (!isdigit(arg[3][0]))
376 {
377
378 error->all(FLERR,"Incorrect args for pair coefficients");
379 }
380 int index = force->inumeric(FLERR,arg[3]);
381 if (index == multiple[m]) break;
382 else continue;
383 } else break;
384 }
385 }
386
387 int none = 0;
388 if (m == nstyles) {
389 if (strcmp(arg[2],"none") == 0) none = 1;
390 else error->all(FLERR,"Pair coeff for hybrid has invalid style");
391 }
392
393 // move 1st/2nd args to 2nd/3rd args
394 // if multflag: move 1st/2nd args to 3rd/4th args
395 // just copy ptrs, since arg[] points into original input line
396
397 arg[2+multflag] = arg[1];
398 arg[1+multflag] = arg[0];
399
400 // invoke sub-style coeff() starting with 1st remaining arg
401
402 if (!none) styles[m]->coeff(narg-1-multflag,&arg[1+multflag]);
403
404 // if sub-style only allows one pair coeff call (with * * and type mapping)
405 // then unset setflag/map assigned to that style before setting it below
406 // in case pair coeff for this sub-style is being called for 2nd time
407
408 if (!none && styles[m]->one_coeff)
409 for (int i = 1; i <= atom->ntypes; i++)
410 for (int j = i; j <= atom->ntypes; j++)
411 if (nmap[i][j] && map[i][j][0] == m) {
412 setflag[i][j] = 0;
413 nmap[i][j] = 0;
414 }
415
416 // set setflag and which type pairs map to which sub-style
417 // if sub-style is none: set hybrid setflag, wipe out map
418 // else: set hybrid setflag & map only if substyle setflag is set
419 // previous mappings are wiped out
420
421 int count = 0;
422 for (int i = ilo; i <= ihi; i++) {
423 for (int j = MAX(jlo,i); j <= jhi; j++) {
424 if (none) {
425 setflag[i][j] = 1;
426 nmap[i][j] = 0;
427 count++;
428 } else if (styles[m]->setflag[i][j]) {
429 setflag[i][j] = 1;
430 nmap[i][j] = 1;
431 map[i][j][0] = m;
432 count++;
433 }
434 }
435 }
436
437 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
438 }
439
440 /* ----------------------------------------------------------------------
441 init specific to this pair style
442 ------------------------------------------------------------------------- */
443
init_style()444 void PairHybrid::init_style()
445 {
446 int i,m,itype,jtype,used,istyle,skip;
447
448 // error if a sub-style is not used
449
450 int ntypes = atom->ntypes;
451
452 for (istyle = 0; istyle < nstyles; istyle++) {
453 used = 0;
454 for (itype = 1; itype <= ntypes; itype++)
455 for (jtype = itype; jtype <= ntypes; jtype++)
456 for (m = 0; m < nmap[itype][jtype]; m++)
457 if (map[itype][jtype][m] == istyle) used = 1;
458 if (used == 0) error->all(FLERR,"Pair hybrid sub-style is not used");
459 }
460
461 // each sub-style makes its neighbor list request(s)
462
463 for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style();
464
465 // create skip lists for each pair neigh request
466 // any kind of list can have its skip flag set at this stage
467
468 for (i = 0; i < neighbor->nrequest; i++) {
469 if (!neighbor->requests[i]->pair) continue;
470
471 // istyle = associated sub-style
472
473 for (istyle = 0; istyle < nstyles; istyle++)
474 if (styles[istyle] == neighbor->requests[i]->requestor) break;
475
476 // allocate iskip and ijskip
477 // initialize so as to skip all pair types
478 // set ijskip = 0 if type pair matches any entry in sub-style map
479 // set ijskip = 0 if mixing will assign type pair to this sub-style
480 // will occur if type pair is currently unassigned
481 // and both I,I and J,J are assigned to single sub-style
482 // and sub-style for both I,I and J,J match istyle
483 // set iskip = 1 only if all ijskip for itype are 1
484
485 int *iskip = new int[ntypes+1];
486 int **ijskip;
487 memory->create(ijskip,ntypes+1,ntypes+1,"pair_hybrid:ijskip");
488
489 for (itype = 1; itype <= ntypes; itype++)
490 for (jtype = 1; jtype <= ntypes; jtype++)
491 ijskip[itype][jtype] = 1;
492
493 for (itype = 1; itype <= ntypes; itype++)
494 for (jtype = itype; jtype <= ntypes; jtype++) {
495 for (m = 0; m < nmap[itype][jtype]; m++)
496 if (map[itype][jtype][m] == istyle)
497 ijskip[itype][jtype] = ijskip[jtype][itype] = 0;
498 if (nmap[itype][jtype] == 0 &&
499 nmap[itype][itype] == 1 && map[itype][itype][0] == istyle &&
500 nmap[jtype][jtype] == 1 && map[jtype][jtype][0] == istyle)
501 ijskip[itype][jtype] = ijskip[jtype][itype] = 0;
502 }
503
504 for (itype = 1; itype <= ntypes; itype++) {
505 iskip[itype] = 1;
506 for (jtype = 1; jtype <= ntypes; jtype++)
507 if (ijskip[itype][jtype] == 0) iskip[itype] = 0;
508 }
509
510 // if any skipping occurs
511 // set request->skip and copy iskip and ijskip into request
512 // else delete iskip and ijskip
513
514 skip = 0;
515 for (itype = 1; itype <= ntypes; itype++)
516 for (jtype = 1; jtype <= ntypes; jtype++)
517 if (ijskip[itype][jtype] == 1) skip = 1;
518
519 if (skip) {
520 neighbor->requests[i]->skip = 1;
521 neighbor->requests[i]->iskip = iskip;
522 neighbor->requests[i]->ijskip = ijskip;
523 } else {
524 delete [] iskip;
525 memory->destroy(ijskip);
526 }
527 }
528
529 // combine sub-style neigh list requests and create new ones if needed
530
531 modify_requests();
532 }
533
534 /* ----------------------------------------------------------------------
535 init for one type pair i,j and corresponding j,i
536 ------------------------------------------------------------------------- */
537
init_one(int i,int j)538 double PairHybrid::init_one(int i, int j)
539 {
540 // if I,J is not set explicitly:
541 // perform mixing only if I,I sub-style = J,J sub-style
542 // also require I,I and J,J are both assigned to single sub-style
543
544 if (setflag[i][j] == 0) {
545 if (nmap[i][i] != 1 || nmap[j][j] != 1 || map[i][i][0] != map[j][j][0])
546 error->one(FLERR,"All pair coeffs are not set");
547 nmap[i][j] = 1;
548 map[i][j][0] = map[i][i][0];
549 }
550
551 // call init/mixing for all sub-styles of I,J
552 // set cutsq in sub-style just as Pair::init() does via call to init_one()
553 // set cutghost for I,J and J,I just as sub-style does
554 // sum tail corrections for I,J
555 // return max cutoff of all sub-styles assigned to I,J
556 // if no sub-styles assigned to I,J (pair_coeff none), cutmax = 0.0 returned
557
558 double cutmax = 0.0;
559 cutghost[i][j] = cutghost[j][i] = 0.0;
560 if (tail_flag) etail_ij = ptail_ij = 0.0;
561
562 nmap[j][i] = nmap[i][j];
563
564 for (int k = 0; k < nmap[i][j]; k++) {
565 map[j][i][k] = map[i][j][k];
566 double cut = styles[map[i][j][k]]->init_one(i,j);
567 styles[map[i][j][k]]->cutsq[i][j] =
568 styles[map[i][j][k]]->cutsq[j][i] = cut*cut;
569 if (styles[map[i][j][k]]->ghostneigh)
570 cutghost[i][j] = cutghost[j][i] =
571 MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]);
572 if (tail_flag) {
573 etail_ij += styles[map[i][j][k]]->etail_ij;
574 ptail_ij += styles[map[i][j][k]]->ptail_ij;
575 }
576 cutmax = MAX(cutmax,cut);
577 }
578
579 return cutmax;
580 }
581
582 /* ----------------------------------------------------------------------
583 combine sub-style neigh list requests and create new ones if needed
584 ------------------------------------------------------------------------- */
585
modify_requests()586 void PairHybrid::modify_requests()
587 {
588 int i,j;
589 NeighRequest *irq,*jrq;
590
591 // loop over pair requests only
592 // if list is skip list and not copy, look for non-skip list of same kind
593 // if one exists, point at that one via otherlist
594 // else make new non-skip request of same kind and point at that one
595 // don't bother to set ID for new request, since pair hybrid ignores list
596 // only exception is half_from_full:
597 // ignore it, turn off skip, since it will derive from its skip parent
598 // after possible new request creation, unset skip flag and otherlist
599 // for these derived lists: granhistory, rRESPA inner/middle
600 // this prevents neighbor from treating them as skip lists
601 // copy list check is for pair style = hybrid/overlay
602 // which invokes this routine
603
604 for (i = 0; i < neighbor->nrequest; i++) {
605 if (!neighbor->requests[i]->pair) continue;
606
607 irq = neighbor->requests[i];
608 if (irq->skip == 0 || irq->copy) continue;
609 if (irq->half_from_full) {
610 irq->skip = 0;
611 continue;
612 }
613
614 for (j = 0; j < neighbor->nrequest; j++) {
615 if (!neighbor->requests[j]->pair) continue;
616 jrq = neighbor->requests[j];
617 if (irq->same_kind(jrq) && jrq->skip == 0) break;
618 }
619
620 if (j < neighbor->nrequest) irq->otherlist = j;
621 else {
622 int newrequest = neighbor->request(this);
623 neighbor->requests[newrequest]->copy_request(irq);
624 irq->otherlist = newrequest;
625 }
626
627 if (irq->granhistory || irq->respainner || irq->respamiddle) {
628 irq->skip = 0;
629 irq->otherlist = -1;
630 }
631 }
632 }
633
634 /* ----------------------------------------------------------------------
635 proc 0 writes to restart file
636 ------------------------------------------------------------------------- */
637
write_restart(FILE * fp)638 void PairHybrid::write_restart(FILE *fp)
639 {
640 fwrite(&nstyles,sizeof(int),1,fp);
641
642 // each sub-style writes its settings, but no coeff info
643
644 int n;
645 for (int m = 0; m < nstyles; m++) {
646 n = strlen(keywords[m]) + 1;
647 fwrite(&n,sizeof(int),1,fp);
648 fwrite(keywords[m],sizeof(char),n,fp);
649 styles[m]->write_restart_settings(fp);
650 }
651 }
652
653 /* ----------------------------------------------------------------------
654 proc 0 reads from restart file, bcasts
655 ------------------------------------------------------------------------- */
656
read_restart(FILE * fp,const int major,const int minor)657 void PairHybrid::read_restart(FILE *fp, const int major, const int minor)
658 {
659 int me = comm->me;
660 if (me == 0) fread(&nstyles,sizeof(int),1,fp);
661 MPI_Bcast(&nstyles,1,MPI_INT,0,world);
662
663 // allocate list of sub-styles
664
665 styles = new Pair*[nstyles];
666 keywords = new char*[nstyles];
667 multiple = new int[nstyles];
668
669 // each sub-style is created via new_pair()
670 // each reads its settings, but no coeff info
671
672 int n,dummy;
673 for (int m = 0; m < nstyles; m++) {
674 if (me == 0) fread(&n,sizeof(int),1,fp);
675 MPI_Bcast(&n,1,MPI_INT,0,world);
676 keywords[m] = new char[n];
677 if (me == 0) fread(keywords[m],sizeof(char),n,fp);
678 MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
679 styles[m] = force->new_pair_from_restart(fp, keywords[m],lmp->suffix,dummy);
680 styles[m]->read_restart_settings(fp, major, minor);
681 }
682
683 // multiple[i] = 1 to M if sub-style used multiple times, else 0
684
685 for (int i = 0; i < nstyles; i++) {
686 int count = 0;
687 for (int j = 0; j < nstyles; j++) {
688 if (strcmp(keywords[j],keywords[i]) == 0) count++;
689 if (j == i) multiple[i] = count;
690 }
691 if (count == 1) multiple[i] = 0;
692 }
693
694 // set pair flags from sub-style flags
695
696 flags();
697 }
698
699 /* ----------------------------------------------------------------------
700 call sub-style to compute single interaction
701 error if sub-style does not support single() call
702 since overlay could have multiple sub-styles, sum results explicitly
703 ------------------------------------------------------------------------- */
704
single(int i,int j,int itype,int jtype,double rsq,double factor_coul,double factor_lj,double & fforce)705 double PairHybrid::single(int i, int j, int itype, int jtype,
706 double rsq, double factor_coul, double factor_lj,
707 double &fforce)
708 {
709 if (nmap[itype][jtype] == 0)
710 error->one(FLERR,"Invoked pair single on pair style none");
711
712 double fone;
713 fforce = 0.0;
714 double esum = 0.0;
715
716 for (int m = 0; m < nmap[itype][jtype]; m++) {
717 if (rsq < styles[map[itype][jtype][m]]->cutsq[itype][jtype]) {
718 if (styles[map[itype][jtype][m]]->single_enable == 0)
719 error->one(FLERR,"Pair hybrid sub-style does not support single call");
720
721 esum += styles[map[itype][jtype][m]]->
722 single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fone);
723 fforce += fone;
724
725 // copy substyle extra values into hybrid's svector
726
727 if (single_extra && styles[map[itype][jtype][m]]->single_extra)
728 for (m = 0; m < single_extra; m++)
729 svector[m] = styles[map[itype][jtype][m]]->svector[m];
730 }
731 }
732
733 return esum;
734 }
735
736 /* ----------------------------------------------------------------------
737 modify parameters of the pair style
738 call modify_params of PairHybrid
739 also pass command args to each sub-style of hybrid
740 ------------------------------------------------------------------------- */
741
modify_params(int narg,char ** arg)742 void PairHybrid::modify_params(int narg, char **arg)
743 {
744 Pair::modify_params(narg,arg);
745 for (int m = 0; m < nstyles; m++) styles[m]->modify_params(narg,arg);
746 }
747
748 /* ----------------------------------------------------------------------
749 extract a ptr to a particular quantity stored by pair
750 pass request thru to sub-styles
751 return first non-NULL result except for cut_coul request
752 for cut_coul, insure all non-NULL results are equal since required by Kspace
753 ------------------------------------------------------------------------- */
754
extract(const char * str,int & dim)755 void *PairHybrid::extract(const char *str, int &dim)
756 {
757 void *cutptr = NULL;
758 void *ptr;
759 double cutvalue = 0.0;
760
761 for (int m = 0; m < nstyles; m++) {
762 ptr = styles[m]->extract(str,dim);
763 if (ptr && strcmp(str,"cut_coul") == 0) {
764 double *p_newvalue = (double *) ptr;
765 double newvalue = *p_newvalue;
766 if (cutptr && newvalue != cutvalue)
767 error->all(FLERR,
768 "Coulomb cutoffs of pair hybrid sub-styles do not match");
769 cutptr = ptr;
770 cutvalue = newvalue;
771 } else if (ptr) return ptr;
772 }
773
774 if (strcmp(str,"cut_coul") == 0) return cutptr;
775 return NULL;
776 }
777
778 /* ---------------------------------------------------------------------- */
779
reset_dt()780 void PairHybrid::reset_dt()
781 {
782 for (int m = 0; m < nstyles; m++) styles[m]->reset_dt();
783 }
784
785 /* ----------------------------------------------------------------------
786 check if itype,jtype maps to sub-style
787 ------------------------------------------------------------------------- */
788
check_ijtype(int itype,int jtype,char * substyle)789 int PairHybrid::check_ijtype(int itype, int jtype, char *substyle)
790 {
791 for (int m = 0; m < nmap[itype][jtype]; m++)
792 if (strcmp(keywords[map[itype][jtype][m]],substyle) == 0) return 1;
793 return 0;
794 }
795
796 /* ----------------------------------------------------------------------
797 memory usage of each sub-style
798 ------------------------------------------------------------------------- */
799
memory_usage()800 double PairHybrid::memory_usage()
801 {
802 double bytes = maxeatom * sizeof(double);
803 bytes += maxvatom*6 * sizeof(double);
804 for (int m = 0; m < nstyles; m++) bytes += styles[m]->memory_usage();
805 return bytes;
806 }
807