1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing author: Greg Wagner (SNL)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_meam.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "error.h"
24 #include "force.h"
25 #include "meam.h"
26 #include "memory.h"
27 #include "neigh_list.h"
28 #include "neigh_request.h"
29 #include "neighbor.h"
30 #include "potential_file_reader.h"
31 #include "tokenizer.h"
32
33 #include <cstring>
34 #include <memory>
35
36 using namespace LAMMPS_NS;
37
38 #define MAXLINE 1024
39
40 static const int nkeywords = 22;
41 static const char *keywords[] = {
42 "Ec","alpha","rho0","delta","lattce",
43 "attrac","repuls","nn2","Cmin","Cmax","rc","delr",
44 "augt1","gsmooth_factor","re","ialloy",
45 "mixture_ref_t","erose_form","zbl",
46 "emb_lin_neg","bkgd_dyn", "theta"};
47
48 /* ---------------------------------------------------------------------- */
49
PairMEAM(LAMMPS * lmp)50 PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp)
51 {
52 single_enable = 0;
53 restartinfo = 0;
54 one_coeff = 1;
55 manybody_flag = 1;
56 centroidstressflag = CENTROID_NOTAVAIL;
57
58 allocated = 0;
59
60 nlibelements = 0;
61 meam_inst = new MEAM(memory);
62 scale = nullptr;
63
64 // set comm size needed by this Pair
65
66 comm_forward = 38;
67 comm_reverse = 30;
68 }
69
70 /* ----------------------------------------------------------------------
71 free all arrays
72 check if allocated, since class can be destructed when incomplete
73 ------------------------------------------------------------------------- */
74
~PairMEAM()75 PairMEAM::~PairMEAM()
76 {
77 delete meam_inst;
78
79 if (allocated) {
80 memory->destroy(setflag);
81 memory->destroy(cutsq);
82 memory->destroy(scale);
83 }
84 }
85
86 /* ---------------------------------------------------------------------- */
87
compute(int eflag,int vflag)88 void PairMEAM::compute(int eflag, int vflag)
89 {
90 int i,ii,n,inum_half,errorflag;
91 int *ilist_half,*numneigh_half,**firstneigh_half;
92 int *numneigh_full,**firstneigh_full;
93
94 ev_init(eflag,vflag);
95
96 // neighbor list info
97
98 inum_half = listhalf->inum;
99 ilist_half = listhalf->ilist;
100 numneigh_half = listhalf->numneigh;
101 firstneigh_half = listhalf->firstneigh;
102 numneigh_full = listfull->numneigh;
103 firstneigh_full = listfull->firstneigh;
104
105 // strip neighbor lists of any special bond flags before using with MEAM
106 // necessary before doing neigh_f2c and neigh_c2f conversions each step
107
108 if (neighbor->ago == 0) {
109 neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half);
110 neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full);
111 }
112
113 // check size of scrfcn based on half neighbor list
114
115 int nlocal = atom->nlocal;
116 int nall = nlocal + atom->nghost;
117
118 n = 0;
119 for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]];
120
121 meam_inst->meam_dens_setup(atom->nmax, nall, n);
122
123 double **x = atom->x;
124 double **f = atom->f;
125 int *type = atom->type;
126 int ntype = atom->ntypes;
127
128 // 3 stages of MEAM calculation
129 // loop over my atoms followed by communication
130
131 int offset = 0;
132 errorflag = 0;
133
134 for (ii = 0; ii < inum_half; ii++) {
135 i = ilist_half[ii];
136 meam_inst->meam_dens_init(i,ntype,type,map,x,
137 numneigh_half[i],firstneigh_half[i],
138 numneigh_full[i],firstneigh_full[i],
139 offset);
140 offset += numneigh_half[i];
141 }
142
143 comm->reverse_comm_pair(this);
144
145 meam_inst->meam_dens_final(nlocal,eflag_either,eflag_global,eflag_atom,
146 &eng_vdwl,eatom,ntype,type,map,scale,errorflag);
147 if (errorflag)
148 error->one(FLERR,"MEAM library error {}",errorflag);
149
150 comm->forward_comm_pair(this);
151
152 offset = 0;
153
154 // vptr is first value in vatom if it will be used by meam_force()
155 // else vatom may not exist, so pass dummy ptr
156
157 double **vptr = nullptr;
158 if (vflag_atom) vptr = vatom;
159
160 for (ii = 0; ii < inum_half; ii++) {
161 i = ilist_half[ii];
162 meam_inst->meam_force(i,eflag_global,eflag_atom,vflag_global,
163 vflag_atom,&eng_vdwl,eatom,ntype,type,map,scale,x,
164 numneigh_half[i],firstneigh_half[i],
165 numneigh_full[i],firstneigh_full[i],
166 offset,f,vptr,virial);
167 offset += numneigh_half[i];
168 }
169
170 if (vflag_fdotr) virial_fdotr_compute();
171 }
172
173 /* ---------------------------------------------------------------------- */
174
allocate()175 void PairMEAM::allocate()
176 {
177 allocated = 1;
178 int n = atom->ntypes;
179
180 memory->create(setflag,n+1,n+1,"pair:setflag");
181 memory->create(cutsq,n+1,n+1,"pair:cutsq");
182 memory->create(scale,n+1,n+1,"pair:scale");
183
184 map = new int[n+1];
185 }
186
187 /* ----------------------------------------------------------------------
188 global settings
189 ------------------------------------------------------------------------- */
190
settings(int narg,char **)191 void PairMEAM::settings(int narg, char ** /*arg*/)
192 {
193 if (narg != 0) error->all(FLERR,"Illegal pair_style command");
194 }
195
196 /* ----------------------------------------------------------------------
197 set coeffs for one or more type pairs
198 ------------------------------------------------------------------------- */
199
coeff(int narg,char ** arg)200 void PairMEAM::coeff(int narg, char **arg)
201 {
202 int m,n;
203
204 if (!allocated) allocate();
205
206 if (narg < 6) error->all(FLERR,"Incorrect args for pair coefficients");
207
208 // insure I,J args are * *
209
210 if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
211 error->all(FLERR,"Incorrect args for pair coefficients");
212
213 // check for presence of first meam file
214
215 std::string lib_file = utils::get_potential_file_path(arg[2]);
216 if (lib_file.empty())
217 error->all(FLERR,"Cannot open MEAM library file {}",lib_file);
218
219 // find meam parameter file in arguments:
220 // first word that is a file or "NULL" after the MEAM library file
221 // we need to extract at least one element, so start from index 4
222
223 int paridx=-1;
224 std::string par_file;
225 for (int i = 4; i < narg; ++i) {
226 if (strcmp(arg[i],"NULL") == 0) {
227 par_file = "NULL";
228 paridx = i;
229 break;
230 }
231 par_file = utils::get_potential_file_path(arg[i]);
232 if (!par_file.empty()) {
233 paridx=i;
234 break;
235 }
236 }
237 if (paridx < 0) error->all(FLERR,"No MEAM parameter file in pair coefficients");
238 if ((narg - paridx - 1) != atom->ntypes)
239 error->all(FLERR,"Incorrect args for pair coefficients");
240
241 // MEAM element names between 2 filenames
242 // nlibelements = # of MEAM elements
243 // elements = list of unique element names
244
245 if (nlibelements) {
246 libelements.clear();
247 mass.clear();
248 }
249
250 nlibelements = paridx - 3;
251 if (nlibelements < 1) error->all(FLERR,"Incorrect args for pair coefficients");
252 if (nlibelements > maxelt)
253 error->all(FLERR,"Too many elements extracted from MEAM "
254 "library (current limit: {}). Increase "
255 "'maxelt' in meam.h and recompile.", maxelt);
256
257 for (int i = 0; i < nlibelements; i++) {
258 libelements.push_back(arg[i+3]);
259 mass.push_back(0.0);
260 }
261
262 // read MEAM library and parameter files
263 // pass all parameters to MEAM package
264 // tell MEAM package that setup is done
265
266 read_files(lib_file,par_file);
267 meam_inst->meam_setup_done(&cutmax);
268
269 // read args that map atom types to MEAM elements
270 // map[i] = which element the Ith atom type is, -1 if not mapped
271
272 for (int i = 4 + nlibelements; i < narg; i++) {
273 m = i - (4+nlibelements) + 1;
274 int j;
275 for (j = 0; j < nlibelements; j++)
276 if (libelements[j] == arg[i]) break;
277 if (j < nlibelements) map[m] = j;
278 else if (strcmp(arg[i],"NULL") == 0) map[m] = -1;
279 else error->all(FLERR,"Incorrect args for pair coefficients");
280 }
281
282 // clear setflag since coeff() called once with I,J = * *
283
284 n = atom->ntypes;
285 for (int i = 1; i <= n; i++)
286 for (int j = i; j <= n; j++)
287 setflag[i][j] = 0;
288
289 // set setflag i,j for type pairs where both are mapped to elements
290 // set mass for i,i in atom class
291
292 int count = 0;
293 for (int i = 1; i <= n; i++) {
294 for (int j = i; j <= n; j++) {
295 if (map[i] >= 0 && map[j] >= 0) {
296 setflag[i][j] = 1;
297 if (i == j) atom->set_mass(FLERR,i,mass[map[i]]);
298 count++;
299 }
300 scale[i][j] = 1.0;
301 }
302 }
303
304 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
305 }
306
307 /* ----------------------------------------------------------------------
308 init specific to this pair style
309 ------------------------------------------------------------------------- */
310
init_style()311 void PairMEAM::init_style()
312 {
313 if (force->newton_pair == 0)
314 error->all(FLERR,"Pair style MEAM requires newton pair on");
315
316 // need full and half neighbor list
317
318 int irequest_full = neighbor->request(this,instance_me);
319 neighbor->requests[irequest_full]->id = 1;
320 neighbor->requests[irequest_full]->half = 0;
321 neighbor->requests[irequest_full]->full = 1;
322 int irequest_half = neighbor->request(this,instance_me);
323 neighbor->requests[irequest_half]->id = 2;
324 }
325
326 /* ----------------------------------------------------------------------
327 neighbor callback to inform pair style of neighbor list to use
328 half or full
329 ------------------------------------------------------------------------- */
330
init_list(int id,NeighList * ptr)331 void PairMEAM::init_list(int id, NeighList *ptr)
332 {
333 if (id == 1) listfull = ptr;
334 else if (id == 2) listhalf = ptr;
335 }
336
337 /* ----------------------------------------------------------------------
338 init for one type pair i,j and corresponding j,i
339 ------------------------------------------------------------------------- */
340
init_one(int i,int j)341 double PairMEAM::init_one(int i, int j)
342 {
343 if (setflag[i][j] == 0) scale[i][j] = 1.0;
344 scale[j][i] = scale[i][j];
345 return cutmax;
346 }
347
348 /* ---------------------------------------------------------------------- */
349
read_files(const std::string & globalfile,const std::string & userfile)350 void PairMEAM::read_files(const std::string &globalfile,
351 const std::string &userfile)
352 {
353 read_global_meam_file(globalfile);
354 read_user_meam_file(userfile);
355 }
356
357 /* ---------------------------------------------------------------------- */
358
read_global_meam_file(const std::string & globalfile)359 void PairMEAM::read_global_meam_file(const std::string &globalfile)
360 {
361 // allocate parameter arrays
362 std::vector<lattice_t> lat(nlibelements);
363 std::vector<int> ielement(nlibelements);
364 std::vector<int> ibar(nlibelements);
365 std::vector<double> z(nlibelements);
366 std::vector<double> atwt(nlibelements);
367 std::vector<double> alpha(nlibelements);
368 std::vector<double> b0(nlibelements);
369 std::vector<double> b1(nlibelements);
370 std::vector<double> b2(nlibelements);
371 std::vector<double> b3(nlibelements);
372 std::vector<double> alat(nlibelements);
373 std::vector<double> esub(nlibelements);
374 std::vector<double> asub(nlibelements);
375 std::vector<double> t0(nlibelements);
376 std::vector<double> t1(nlibelements);
377 std::vector<double> t2(nlibelements);
378 std::vector<double> t3(nlibelements);
379 std::vector<double> rozero(nlibelements);
380 std::vector<bool> found(nlibelements, false);
381
382 // open global meamf file on proc 0
383
384 if (comm->me == 0) {
385 PotentialFileReader reader(lmp, globalfile, "MEAM", " library");
386 char * line;
387
388 const int params_per_line = 19;
389 int nset = 0;
390
391 while ((line = reader.next_line(params_per_line))) {
392 try {
393 ValueTokenizer values(line, "' \t\n\r\f");
394
395 // read each set of params from global MEAM file
396 // one set of params can span multiple lines
397 // store params if element name is in element list
398 // if element name appears multiple times, only store 1st entry
399 std::string element = values.next_string();
400
401 // skip if element name isn't in element list
402
403 int index;
404 for (index = 0; index < nlibelements; index++)
405 if (libelements[index] == element) break;
406 if (index == nlibelements) continue;
407
408 // skip if element already appeared (technically error in library file, but always ignored)
409
410 if (found[index]) continue;
411 found[index] = true;
412
413 // map lat string to an integer
414 std::string lattice_type = values.next_string();
415
416 if (!MEAM::str_to_lat(lattice_type.c_str(), true, lat[index]))
417 error->one(FLERR,"Unrecognized lattice type in MEAM "
418 "library file: {}", lattice_type);
419
420 // store parameters
421
422 z[index] = values.next_double();
423 ielement[index] = values.next_int();
424 atwt[index] = values.next_double();
425 alpha[index] = values.next_double();
426 b0[index] = values.next_double();
427 b1[index] = values.next_double();
428 b2[index] = values.next_double();
429 b3[index] = values.next_double();
430 alat[index] = values.next_double();
431 esub[index] = values.next_double();
432 asub[index] = values.next_double();
433 t0[index] = values.next_double();
434 t1[index] = values.next_double();
435 t2[index] = values.next_double();
436 t3[index] = values.next_double();
437 rozero[index] = values.next_double();
438 ibar[index] = values.next_int();
439
440 if (!isone(t0[index]))
441 error->one(FLERR,"Unsupported parameter in MEAM library file: t0!=1");
442
443 // z given is ignored: if this is mismatched, we definitely won't do what the user said -> fatal error
444 if (z[index] != MEAM::get_Zij(lat[index]))
445 error->one(FLERR,"Mismatched parameter in MEAM library file: z!=lat");
446
447 nset++;
448 } catch (TokenizerException &e) {
449 error->one(FLERR, e.what());
450 }
451 }
452
453 // error if didn't find all elements in file
454
455 if (nset != nlibelements) {
456 std::string msg = "Did not find all elements in MEAM library file, missing:";
457 for (int i = 0; i < nlibelements; i++)
458 if (!found[i]) {
459 msg += " ";
460 msg += libelements[i];
461 }
462 error->one(FLERR,msg);
463 }
464 }
465
466 // distribute complete parameter sets
467 MPI_Bcast(lat.data(), nlibelements, MPI_INT, 0, world);
468 MPI_Bcast(ielement.data(), nlibelements, MPI_INT, 0, world);
469 MPI_Bcast(ibar.data(), nlibelements, MPI_INT, 0, world);
470 MPI_Bcast(z.data(), nlibelements, MPI_DOUBLE, 0, world);
471 MPI_Bcast(atwt.data(), nlibelements, MPI_DOUBLE, 0, world);
472 MPI_Bcast(alpha.data(), nlibelements, MPI_DOUBLE, 0, world);
473 MPI_Bcast(b0.data(), nlibelements, MPI_DOUBLE, 0, world);
474 MPI_Bcast(b1.data(), nlibelements, MPI_DOUBLE, 0, world);
475 MPI_Bcast(b2.data(), nlibelements, MPI_DOUBLE, 0, world);
476 MPI_Bcast(b3.data(), nlibelements, MPI_DOUBLE, 0, world);
477 MPI_Bcast(alat.data(), nlibelements, MPI_DOUBLE, 0, world);
478 MPI_Bcast(esub.data(), nlibelements, MPI_DOUBLE, 0, world);
479 MPI_Bcast(asub.data(), nlibelements, MPI_DOUBLE, 0, world);
480 MPI_Bcast(t0.data(), nlibelements, MPI_DOUBLE, 0, world);
481 MPI_Bcast(t1.data(), nlibelements, MPI_DOUBLE, 0, world);
482 MPI_Bcast(t2.data(), nlibelements, MPI_DOUBLE, 0, world);
483 MPI_Bcast(t3.data(), nlibelements, MPI_DOUBLE, 0, world);
484 MPI_Bcast(rozero.data(), nlibelements, MPI_DOUBLE, 0, world);
485
486 // pass element parameters to MEAM package
487
488 meam_inst->meam_setup_global(nlibelements, lat.data(), ielement.data(), atwt.data(),
489 alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(),
490 alat.data(), esub.data(), asub.data(), t0.data(), t1.data(),
491 t2.data(), t3.data(), rozero.data(), ibar.data());
492
493 // set element masses
494
495 for (int i = 0; i < nlibelements; i++) mass[i] = atwt[i];
496 }
497
498 /* ---------------------------------------------------------------------- */
499
read_user_meam_file(const std::string & userfile)500 void PairMEAM::read_user_meam_file(const std::string &userfile)
501 {
502 // done if user param file is "NULL"
503
504 if (userfile == "NULL") return;
505
506 // open user param file on proc 0
507
508 std::shared_ptr<PotentialFileReader> reader;
509
510 if (comm->me == 0) {
511 reader = std::make_shared<PotentialFileReader>(lmp, userfile, "MEAM");
512 }
513
514 // read settings
515 // pass them one at a time to MEAM package
516 // match strings to list of corresponding ints
517 char * line = nullptr;
518 char buffer[MAXLINE];
519
520 while (1) {
521 int which;
522 int nindex, index[3];
523 double value;
524 int nline;
525 if (comm->me == 0) {
526 line = reader->next_line();
527 if (line == nullptr) {
528 nline = -1;
529 } else nline = strlen(line) + 1;
530 } else {
531 line = buffer;
532 }
533
534 MPI_Bcast(&nline,1,MPI_INT,0,world);
535 if (nline<0) break;
536 MPI_Bcast(line,nline,MPI_CHAR,0,world);
537
538 ValueTokenizer values(line, "=(), '\t\n\r\f");
539 int nparams = values.count();
540 std::string keyword = values.next_string();
541
542 for (which = 0; which < nkeywords; which++)
543 if (keyword == keywords[which]) break;
544 if (which == nkeywords)
545 error->all(FLERR,"Keyword {} in MEAM parameter file not "
546 "recognized", keyword);
547
548 nindex = nparams - 2;
549 for (int i = 0; i < nindex; i++) index[i] = values.next_int() - 1;
550
551 // map lattce_meam value to an integer
552 if (which == 4) {
553 std::string lattice_type = values.next_string();
554 lattice_t latt;
555 if (!MEAM::str_to_lat(lattice_type, false, latt))
556 error->all(FLERR, "Unrecognized lattice type in MEAM "
557 "parameter file: {}", lattice_type);
558 value = latt;
559 }
560 else value = values.next_double();
561
562 // pass single setting to MEAM package
563
564 int errorflag = 0;
565 meam_inst->meam_setup_param(which,value,nindex,index,&errorflag);
566 if (errorflag) {
567 const char *descr[] = { "has an unknown error",
568 "is out of range (please report a bug)",
569 "expected more indices",
570 "has out of range element index"};
571 if ((errorflag < 0) || (errorflag > 3)) errorflag = 0;
572 error->all(FLERR,"Error in MEAM parameter file: keyword {} {}", keyword, descr[errorflag]);
573 }
574 }
575 }
576
577 /* ---------------------------------------------------------------------- */
578
pack_forward_comm(int n,int * list,double * buf,int,int *)579 int PairMEAM::pack_forward_comm(int n, int *list, double *buf,
580 int /*pbc_flag*/, int * /*pbc*/)
581 {
582 int i,j,k,m;
583
584 m = 0;
585 for (i = 0; i < n; i++) {
586 j = list[i];
587 buf[m++] = meam_inst->rho0[j];
588 buf[m++] = meam_inst->rho1[j];
589 buf[m++] = meam_inst->rho2[j];
590 buf[m++] = meam_inst->rho3[j];
591 buf[m++] = meam_inst->frhop[j];
592 buf[m++] = meam_inst->gamma[j];
593 buf[m++] = meam_inst->dgamma1[j];
594 buf[m++] = meam_inst->dgamma2[j];
595 buf[m++] = meam_inst->dgamma3[j];
596 buf[m++] = meam_inst->arho2b[j];
597 buf[m++] = meam_inst->arho1[j][0];
598 buf[m++] = meam_inst->arho1[j][1];
599 buf[m++] = meam_inst->arho1[j][2];
600 buf[m++] = meam_inst->arho2[j][0];
601 buf[m++] = meam_inst->arho2[j][1];
602 buf[m++] = meam_inst->arho2[j][2];
603 buf[m++] = meam_inst->arho2[j][3];
604 buf[m++] = meam_inst->arho2[j][4];
605 buf[m++] = meam_inst->arho2[j][5];
606 for (k = 0; k < 10; k++) buf[m++] = meam_inst->arho3[j][k];
607 buf[m++] = meam_inst->arho3b[j][0];
608 buf[m++] = meam_inst->arho3b[j][1];
609 buf[m++] = meam_inst->arho3b[j][2];
610 buf[m++] = meam_inst->t_ave[j][0];
611 buf[m++] = meam_inst->t_ave[j][1];
612 buf[m++] = meam_inst->t_ave[j][2];
613 buf[m++] = meam_inst->tsq_ave[j][0];
614 buf[m++] = meam_inst->tsq_ave[j][1];
615 buf[m++] = meam_inst->tsq_ave[j][2];
616 }
617
618 return m;
619 }
620
621 /* ---------------------------------------------------------------------- */
622
unpack_forward_comm(int n,int first,double * buf)623 void PairMEAM::unpack_forward_comm(int n, int first, double *buf)
624 {
625 int i,k,m,last;
626
627 m = 0;
628 last = first + n;
629 for (i = first; i < last; i++) {
630 meam_inst->rho0[i] = buf[m++];
631 meam_inst->rho1[i] = buf[m++];
632 meam_inst->rho2[i] = buf[m++];
633 meam_inst->rho3[i] = buf[m++];
634 meam_inst->frhop[i] = buf[m++];
635 meam_inst->gamma[i] = buf[m++];
636 meam_inst->dgamma1[i] = buf[m++];
637 meam_inst->dgamma2[i] = buf[m++];
638 meam_inst->dgamma3[i] = buf[m++];
639 meam_inst->arho2b[i] = buf[m++];
640 meam_inst->arho1[i][0] = buf[m++];
641 meam_inst->arho1[i][1] = buf[m++];
642 meam_inst->arho1[i][2] = buf[m++];
643 meam_inst->arho2[i][0] = buf[m++];
644 meam_inst->arho2[i][1] = buf[m++];
645 meam_inst->arho2[i][2] = buf[m++];
646 meam_inst->arho2[i][3] = buf[m++];
647 meam_inst->arho2[i][4] = buf[m++];
648 meam_inst->arho2[i][5] = buf[m++];
649 for (k = 0; k < 10; k++) meam_inst->arho3[i][k] = buf[m++];
650 meam_inst->arho3b[i][0] = buf[m++];
651 meam_inst->arho3b[i][1] = buf[m++];
652 meam_inst->arho3b[i][2] = buf[m++];
653 meam_inst->t_ave[i][0] = buf[m++];
654 meam_inst->t_ave[i][1] = buf[m++];
655 meam_inst->t_ave[i][2] = buf[m++];
656 meam_inst->tsq_ave[i][0] = buf[m++];
657 meam_inst->tsq_ave[i][1] = buf[m++];
658 meam_inst->tsq_ave[i][2] = buf[m++];
659 }
660 }
661
662 /* ---------------------------------------------------------------------- */
663
pack_reverse_comm(int n,int first,double * buf)664 int PairMEAM::pack_reverse_comm(int n, int first, double *buf)
665 {
666 int i,k,m,last;
667
668 m = 0;
669 last = first + n;
670 for (i = first; i < last; i++) {
671 buf[m++] = meam_inst->rho0[i];
672 buf[m++] = meam_inst->arho2b[i];
673 buf[m++] = meam_inst->arho1[i][0];
674 buf[m++] = meam_inst->arho1[i][1];
675 buf[m++] = meam_inst->arho1[i][2];
676 buf[m++] = meam_inst->arho2[i][0];
677 buf[m++] = meam_inst->arho2[i][1];
678 buf[m++] = meam_inst->arho2[i][2];
679 buf[m++] = meam_inst->arho2[i][3];
680 buf[m++] = meam_inst->arho2[i][4];
681 buf[m++] = meam_inst->arho2[i][5];
682 for (k = 0; k < 10; k++) buf[m++] = meam_inst->arho3[i][k];
683 buf[m++] = meam_inst->arho3b[i][0];
684 buf[m++] = meam_inst->arho3b[i][1];
685 buf[m++] = meam_inst->arho3b[i][2];
686 buf[m++] = meam_inst->t_ave[i][0];
687 buf[m++] = meam_inst->t_ave[i][1];
688 buf[m++] = meam_inst->t_ave[i][2];
689 buf[m++] = meam_inst->tsq_ave[i][0];
690 buf[m++] = meam_inst->tsq_ave[i][1];
691 buf[m++] = meam_inst->tsq_ave[i][2];
692 }
693
694 return m;
695 }
696
697 /* ---------------------------------------------------------------------- */
698
unpack_reverse_comm(int n,int * list,double * buf)699 void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf)
700 {
701 int i,j,k,m;
702
703 m = 0;
704 for (i = 0; i < n; i++) {
705 j = list[i];
706 meam_inst->rho0[j] += buf[m++];
707 meam_inst->arho2b[j] += buf[m++];
708 meam_inst->arho1[j][0] += buf[m++];
709 meam_inst->arho1[j][1] += buf[m++];
710 meam_inst->arho1[j][2] += buf[m++];
711 meam_inst->arho2[j][0] += buf[m++];
712 meam_inst->arho2[j][1] += buf[m++];
713 meam_inst->arho2[j][2] += buf[m++];
714 meam_inst->arho2[j][3] += buf[m++];
715 meam_inst->arho2[j][4] += buf[m++];
716 meam_inst->arho2[j][5] += buf[m++];
717 for (k = 0; k < 10; k++) meam_inst->arho3[j][k] += buf[m++];
718 meam_inst->arho3b[j][0] += buf[m++];
719 meam_inst->arho3b[j][1] += buf[m++];
720 meam_inst->arho3b[j][2] += buf[m++];
721 meam_inst->t_ave[j][0] += buf[m++];
722 meam_inst->t_ave[j][1] += buf[m++];
723 meam_inst->t_ave[j][2] += buf[m++];
724 meam_inst->tsq_ave[j][0] += buf[m++];
725 meam_inst->tsq_ave[j][1] += buf[m++];
726 meam_inst->tsq_ave[j][2] += buf[m++];
727 }
728 }
729
730 /* ----------------------------------------------------------------------
731 memory usage of local atom-based arrays
732 ------------------------------------------------------------------------- */
733
memory_usage()734 double PairMEAM::memory_usage()
735 {
736 double bytes = 11 * meam_inst->nmax * sizeof(double);
737 bytes += (double)(3 + 6 + 10 + 3 + 3 + 3) * meam_inst->nmax * sizeof(double);
738 bytes += (double)3 * meam_inst->maxneigh * sizeof(double);
739 return bytes;
740 }
741
742 /* ----------------------------------------------------------------------
743 strip special bond flags from neighbor list entries
744 are not used with MEAM
745 need to do here so Fortran lib doesn't see them
746 done once per reneighbor so that neigh_f2c and neigh_c2f don't see them
747 ------------------------------------------------------------------------- */
748
neigh_strip(int inum,int * ilist,int * numneigh,int ** firstneigh)749 void PairMEAM::neigh_strip(int inum, int *ilist,
750 int *numneigh, int **firstneigh)
751 {
752 int i,j,ii,jnum;
753 int *jlist;
754
755 for (ii = 0; ii < inum; ii++) {
756 i = ilist[ii];
757 jlist = firstneigh[i];
758 jnum = numneigh[i];
759 for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK;
760 }
761 }
762
763 /* ---------------------------------------------------------------------- */
764
extract(const char * str,int & dim)765 void *PairMEAM::extract(const char *str, int &dim)
766 {
767 dim = 2;
768 if (strcmp(str,"scale") == 0) return (void *) scale;
769 return nullptr;
770 }
771