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: Paul Crozier (SNL)
17 ------------------------------------------------------------------------- */
18
19 #include "pair.h"
20
21 #include "atom.h"
22 #include "atom_masks.h"
23 #include "comm.h"
24 #include "compute.h"
25 #include "domain.h"
26 #include "error.h"
27 #include "force.h"
28 #include "kspace.h"
29 #include "math_const.h"
30 #include "memory.h"
31 #include "neighbor.h"
32 #include "suffix.h"
33 #include "update.h"
34
35 #include <cfloat> // IWYU pragma: keep
36 #include <climits> // IWYU pragma: keep
37 #include <cmath>
38 #include <cstring>
39
40 using namespace LAMMPS_NS;
41 using namespace MathConst;
42
43 enum{NONE,RLINEAR,RSQ,BMP};
44
45 // allocate space for static class instance variable and initialize it
46
47 int Pair::instance_total = 0;
48
49 /* ---------------------------------------------------------------------- */
50
Pair(LAMMPS * lmp)51 Pair::Pair(LAMMPS *lmp) : Pointers(lmp)
52 {
53 instance_me = instance_total++;
54
55 eng_vdwl = eng_coul = 0.0;
56
57 comm_forward = comm_reverse = comm_reverse_off = 0;
58
59 single_enable = 1;
60 single_hessian_enable = 0;
61 restartinfo = 1;
62 respa_enable = 0;
63 one_coeff = 0;
64 no_virial_fdotr_compute = 0;
65 writedata = 0;
66 finitecutflag = 0;
67 ghostneigh = 0;
68 unit_convert_flag = utils::NOCONVERT;
69
70 nextra = 0;
71 pvector = nullptr;
72 single_extra = 0;
73 svector = nullptr;
74
75 setflag = nullptr;
76 cutsq = nullptr;
77
78 ewaldflag = pppmflag = msmflag = dispersionflag = tip4pflag = dipoleflag = spinflag = 0;
79 reinitflag = 1;
80 centroidstressflag = CENTROID_SAME;
81
82 // pair_modify settings
83
84 compute_flag = 1;
85 manybody_flag = 0;
86 offset_flag = 0;
87 mix_flag = GEOMETRIC;
88 mixed_flag = 1;
89 tail_flag = 0;
90 etail = ptail = etail_ij = ptail_ij = 0.0;
91 ncoultablebits = 12;
92 ndisptablebits = 12;
93 tabinner = sqrt(2.0);
94 tabinner_disp = sqrt(2.0);
95 ftable = nullptr;
96 fdisptable = nullptr;
97
98 allocated = 0;
99 suffix_flag = Suffix::NONE;
100
101 maxeatom = maxvatom = maxcvatom = 0;
102 eatom = nullptr;
103 vatom = nullptr;
104 cvatom = nullptr;
105
106 num_tally_compute = 0;
107 list_tally_compute = nullptr;
108
109 nelements = nparams = maxparam = 0;
110 elements = nullptr;
111 elem1param = nullptr;
112 elem2param = nullptr;
113 elem3param = nullptr;
114 map = nullptr;
115
116 nondefault_history_transfer = 0;
117 beyond_contact = 0;
118
119 // KOKKOS per-fix data masks
120
121 execution_space = Host;
122 datamask_read = ALL_MASK;
123 datamask_modify = ALL_MASK;
124
125 kokkosable = 0;
126 copymode = 0;
127 }
128
129 /* ---------------------------------------------------------------------- */
130
~Pair()131 Pair::~Pair()
132 {
133 num_tally_compute = 0;
134 memory->sfree((void *) list_tally_compute);
135 list_tally_compute = nullptr;
136
137 if (copymode) return;
138
139 if (elements)
140 for (int i = 0; i < nelements; i++) delete[] elements[i];
141 delete[] elements;
142
143 delete[] map;
144 memory->destroy(eatom);
145 memory->destroy(vatom);
146 memory->destroy(cvatom);
147 }
148
149 /* ----------------------------------------------------------------------
150 modify parameters of the pair style
151 pair_hybrid has its own version of this routine
152 to apply modifications to each of its sub-styles
153 ------------------------------------------------------------------------- */
154
modify_params(int narg,char ** arg)155 void Pair::modify_params(int narg, char **arg)
156 {
157 if (narg == 0) error->all(FLERR,"Illegal pair_modify command");
158
159 int iarg = 0;
160 while (iarg < narg) {
161 if (strcmp(arg[iarg],"mix") == 0) {
162 if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
163 if (strcmp(arg[iarg+1],"geometric") == 0) mix_flag = GEOMETRIC;
164 else if (strcmp(arg[iarg+1],"arithmetic") == 0) mix_flag = ARITHMETIC;
165 else if (strcmp(arg[iarg+1],"sixthpower") == 0) mix_flag = SIXTHPOWER;
166 else error->all(FLERR,"Illegal pair_modify command");
167 iarg += 2;
168 } else if (strcmp(arg[iarg],"shift") == 0) {
169 if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
170 if (strcmp(arg[iarg+1],"yes") == 0) offset_flag = 1;
171 else if (strcmp(arg[iarg+1],"no") == 0) offset_flag = 0;
172 else error->all(FLERR,"Illegal pair_modify command");
173 iarg += 2;
174 } else if (strcmp(arg[iarg],"table") == 0) {
175 if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
176 ncoultablebits = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
177 if (ncoultablebits > (int)sizeof(float)*CHAR_BIT)
178 error->all(FLERR,"Too many total bits for bitmapped lookup table");
179 iarg += 2;
180 } else if (strcmp(arg[iarg],"table/disp") == 0) {
181 if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
182 ndisptablebits = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
183 if (ndisptablebits > (int)sizeof(float)*CHAR_BIT)
184 error->all(FLERR,"Too many total bits for bitmapped lookup table");
185 iarg += 2;
186 } else if (strcmp(arg[iarg],"tabinner") == 0) {
187 if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
188 tabinner = utils::numeric(FLERR,arg[iarg+1],false,lmp);
189 iarg += 2;
190 } else if (strcmp(arg[iarg],"tabinner/disp") == 0) {
191 if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
192 tabinner_disp = utils::numeric(FLERR,arg[iarg+1],false,lmp);
193 iarg += 2;
194 } else if (strcmp(arg[iarg],"tail") == 0) {
195 if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
196 if (strcmp(arg[iarg+1],"yes") == 0) tail_flag = 1;
197 else if (strcmp(arg[iarg+1],"no") == 0) tail_flag = 0;
198 else error->all(FLERR,"Illegal pair_modify command");
199 iarg += 2;
200 } else if (strcmp(arg[iarg],"compute") == 0) {
201 if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command");
202 if (strcmp(arg[iarg+1],"yes") == 0) compute_flag = 1;
203 else if (strcmp(arg[iarg+1],"no") == 0) compute_flag = 0;
204 else error->all(FLERR,"Illegal pair_modify command");
205 iarg += 2;
206 } else if (strcmp(arg[iarg],"nofdotr") == 0) {
207 no_virial_fdotr_compute = 1;
208 ++iarg;
209 } else error->all(FLERR,"Illegal pair_modify command");
210 }
211 }
212
213 /* ---------------------------------------------------------------------- */
214
init()215 void Pair::init()
216 {
217 int i,j;
218
219 if (offset_flag && tail_flag)
220 error->all(FLERR,"Cannot have both pair_modify shift and tail set to yes");
221 if (tail_flag && domain->dimension == 2)
222 error->all(FLERR,"Cannot use pair tail corrections with 2d simulations");
223 if (tail_flag && domain->nonperiodic && comm->me == 0)
224 error->warning(FLERR,"Using pair tail corrections with non-periodic system");
225 if (!compute_flag && tail_flag && comm->me == 0)
226 error->warning(FLERR,"Using pair tail corrections with "
227 "pair_modify compute no");
228 if (!compute_flag && offset_flag && comm->me == 0)
229 error->warning(FLERR,"Using pair potential shift with "
230 "pair_modify compute no");
231
232 // for manybody potentials
233 // check if bonded exclusions could invalidate the neighbor list
234
235 if (manybody_flag && (atom->molecular != Atom::ATOMIC)) {
236 int flag = 0;
237 if (atom->nbonds > 0 && force->special_lj[1] == 0.0 &&
238 force->special_coul[1] == 0.0) flag = 1;
239 if (atom->nangles > 0 && force->special_lj[2] == 0.0 &&
240 force->special_coul[2] == 0.0) flag = 1;
241 if (atom->ndihedrals > 0 && force->special_lj[3] == 0.0 &&
242 force->special_coul[3] == 0.0) flag = 1;
243 if (flag && comm->me == 0)
244 error->warning(FLERR,"Using a manybody potential with "
245 "bonds/angles/dihedrals and special_bond exclusions");
246 }
247
248 // I,I coeffs must be set
249 // init_one() will check if I,J is set explicitly or inferred by mixing
250
251 if (!allocated) error->all(FLERR,"All pair coeffs are not set");
252
253 for (i = 1; i <= atom->ntypes; i++)
254 if (setflag[i][i] == 0) error->all(FLERR,"All pair coeffs are not set");
255
256 // style-specific initialization
257
258 init_style();
259
260 // call init_one() for each I,J
261 // set cutsq for each I,J, used to neighbor
262 // cutforce = max of all I,J cutoffs
263
264 cutforce = 0.0;
265 etail = ptail = 0.0;
266 mixed_flag = 1;
267 double cut;
268
269 for (i = 1; i <= atom->ntypes; i++)
270 for (j = i; j <= atom->ntypes; j++) {
271 if ((i != j) && setflag[i][j]) mixed_flag = 0;
272 cut = init_one(i,j);
273 cutsq[i][j] = cutsq[j][i] = cut*cut;
274 cutforce = MAX(cutforce,cut);
275 if (tail_flag) {
276 etail += etail_ij;
277 ptail += ptail_ij;
278 if (i != j) {
279 etail += etail_ij;
280 ptail += ptail_ij;
281 }
282 }
283 }
284 }
285
286 /* ----------------------------------------------------------------------
287 reset all type-based params by invoking init_one() for each I,J
288 called by fix adapt after it changes one or more params
289 ------------------------------------------------------------------------- */
290
reinit()291 void Pair::reinit()
292 {
293 // generalize this error message if reinit() is used by more than fix adapt
294
295 if (!reinitflag)
296 error->all(FLERR,"Fix adapt interface to this pair style not supported");
297
298 etail = ptail = 0.0;
299
300 for (int i = 1; i <= atom->ntypes; i++)
301 for (int j = i; j <= atom->ntypes; j++) {
302 init_one(i,j);
303 if (tail_flag) {
304 etail += etail_ij;
305 ptail += ptail_ij;
306 if (i != j) {
307 etail += etail_ij;
308 ptail += ptail_ij;
309 }
310 }
311 }
312 }
313
314 /* ----------------------------------------------------------------------
315 init specific to a pair style
316 specific pair style can override this function
317 if needs its own error checks
318 if needs another kind of neighbor list
319 request default neighbor list = half list
320 ------------------------------------------------------------------------- */
321
init_style()322 void Pair::init_style()
323 {
324 neighbor->request(this,instance_me);
325 }
326
327 /* ----------------------------------------------------------------------
328 neighbor callback to inform pair style of neighbor list to use
329 specific pair style can override this function
330 ------------------------------------------------------------------------- */
331
init_list(int,NeighList * ptr)332 void Pair::init_list(int /*which*/, NeighList *ptr)
333 {
334 list = ptr;
335 }
336
337 /* ----------------------------------------------------------------------
338 setup Coulomb force tables used in compute routines
339 ------------------------------------------------------------------------- */
340
init_tables(double cut_coul,double * cut_respa)341 void Pair::init_tables(double cut_coul, double *cut_respa)
342 {
343 int masklo,maskhi;
344 double r,grij,expm2,derfc,egamma,fgamma,rsw;
345 double qqrd2e = force->qqrd2e;
346
347 if (force->kspace == nullptr)
348 error->all(FLERR,"Pair style requires a KSpace style");
349 double g_ewald = force->kspace->g_ewald;
350
351 double cut_coulsq = cut_coul * cut_coul;
352
353 tabinnersq = tabinner*tabinner;
354 init_bitmap(tabinner,cut_coul,ncoultablebits,
355 masklo,maskhi,ncoulmask,ncoulshiftbits);
356
357 int ntable = 1;
358 for (int i = 0; i < ncoultablebits; i++) ntable *= 2;
359
360 // linear lookup tables of length N = 2^ncoultablebits
361 // stored value = value at lower edge of bin
362 // d values = delta from lower edge to upper edge of bin
363
364 if (ftable) free_tables();
365
366 memory->create(rtable,ntable,"pair:rtable");
367 memory->create(ftable,ntable,"pair:ftable");
368 memory->create(ctable,ntable,"pair:ctable");
369 memory->create(etable,ntable,"pair:etable");
370 memory->create(drtable,ntable,"pair:drtable");
371 memory->create(dftable,ntable,"pair:dftable");
372 memory->create(dctable,ntable,"pair:dctable");
373 memory->create(detable,ntable,"pair:detable");
374
375 if (cut_respa == nullptr) {
376 vtable = ptable = dvtable = dptable = nullptr;
377 } else {
378 memory->create(vtable,ntable,"pair:vtable");
379 memory->create(ptable,ntable,"pair:ptable");
380 memory->create(dvtable,ntable,"pair:dvtable");
381 memory->create(dptable,ntable,"pair:dptable");
382 }
383
384 union_int_float_t rsq_lookup;
385 union_int_float_t minrsq_lookup;
386 int itablemin;
387 minrsq_lookup.i = 0 << ncoulshiftbits;
388 minrsq_lookup.i |= maskhi;
389
390 for (int i = 0; i < ntable; i++) {
391 rsq_lookup.i = i << ncoulshiftbits;
392 rsq_lookup.i |= masklo;
393 if (rsq_lookup.f < tabinnersq) {
394 rsq_lookup.i = i << ncoulshiftbits;
395 rsq_lookup.i |= maskhi;
396 }
397 r = sqrtf(rsq_lookup.f);
398 if (msmflag) {
399 egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
400 fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*
401 force->kspace->dgamma(r/cut_coul);
402 } else {
403 grij = g_ewald * r;
404 expm2 = exp(-grij*grij);
405 derfc = erfc(grij);
406 }
407 if (cut_respa == nullptr) {
408 rtable[i] = rsq_lookup.f;
409 ctable[i] = qqrd2e/r;
410 if (msmflag) {
411 ftable[i] = qqrd2e/r * fgamma;
412 etable[i] = qqrd2e/r * egamma;
413 } else {
414 ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
415 etable[i] = qqrd2e/r * derfc;
416 }
417 } else {
418 rtable[i] = rsq_lookup.f;
419 ctable[i] = 0.0;
420 ptable[i] = qqrd2e/r;
421 if (msmflag) {
422 ftable[i] = qqrd2e/r * (fgamma - 1.0);
423 etable[i] = qqrd2e/r * egamma;
424 vtable[i] = qqrd2e/r * fgamma;
425 } else {
426 ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2 - 1.0);
427 etable[i] = qqrd2e/r * derfc;
428 vtable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
429 }
430 if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
431 if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
432 rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
433 ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
434 ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
435 } else {
436 if (msmflag) ftable[i] = qqrd2e/r * fgamma;
437 else ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
438 ctable[i] = qqrd2e/r;
439 }
440 }
441 }
442 minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
443 }
444
445 tabinnersq = minrsq_lookup.f;
446
447 int ntablem1 = ntable - 1;
448
449 for (int i = 0; i < ntablem1; i++) {
450 drtable[i] = 1.0/(rtable[i+1] - rtable[i]);
451 dftable[i] = ftable[i+1] - ftable[i];
452 dctable[i] = ctable[i+1] - ctable[i];
453 detable[i] = etable[i+1] - etable[i];
454 }
455
456 if (cut_respa) {
457 for (int i = 0; i < ntablem1; i++) {
458 dvtable[i] = vtable[i+1] - vtable[i];
459 dptable[i] = ptable[i+1] - ptable[i];
460 }
461 }
462
463 // get the delta values for the last table entries
464 // tables are connected periodically between 0 and ntablem1
465
466 drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]);
467 dftable[ntablem1] = ftable[0] - ftable[ntablem1];
468 dctable[ntablem1] = ctable[0] - ctable[ntablem1];
469 detable[ntablem1] = etable[0] - etable[ntablem1];
470 if (cut_respa) {
471 dvtable[ntablem1] = vtable[0] - vtable[ntablem1];
472 dptable[ntablem1] = ptable[0] - ptable[ntablem1];
473 }
474
475 // get the correct delta values at itablemax
476 // smallest r is in bin itablemin
477 // largest r is in bin itablemax, which is itablemin-1,
478 // or ntablem1 if itablemin=0
479 // deltas at itablemax only needed if corresponding rsq < cut*cut
480 // if so, compute deltas between rsq and cut*cut
481
482 double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp;
483 p_tmp = 0.0;
484 v_tmp = 0.0;
485 itablemin = minrsq_lookup.i & ncoulmask;
486 itablemin >>= ncoulshiftbits;
487 int itablemax = itablemin - 1;
488 if (itablemin == 0) itablemax = ntablem1;
489 rsq_lookup.i = itablemax << ncoulshiftbits;
490 rsq_lookup.i |= maskhi;
491
492 if (rsq_lookup.f < cut_coulsq) {
493 rsq_lookup.f = cut_coulsq;
494 r = sqrtf(rsq_lookup.f);
495 if (msmflag) {
496 egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
497 fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*
498 force->kspace->dgamma(r/cut_coul);
499 } else {
500 grij = g_ewald * r;
501 expm2 = exp(-grij*grij);
502 derfc = erfc(grij);
503 }
504 if (cut_respa == nullptr) {
505 c_tmp = qqrd2e/r;
506 if (msmflag) {
507 f_tmp = qqrd2e/r * fgamma;
508 e_tmp = qqrd2e/r * egamma;
509 } else {
510 f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
511 e_tmp = qqrd2e/r * derfc;
512 }
513 } else {
514 c_tmp = 0.0;
515 p_tmp = qqrd2e/r;
516 if (msmflag) {
517 f_tmp = qqrd2e/r * (fgamma - 1.0);
518 e_tmp = qqrd2e/r * egamma;
519 v_tmp = qqrd2e/r * fgamma;
520 } else {
521 f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2 - 1.0);
522 e_tmp = qqrd2e/r * derfc;
523 v_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
524 }
525 if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
526 if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
527 rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
528 f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
529 c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
530 } else {
531 if (msmflag) f_tmp = qqrd2e/r * fgamma;
532 else f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2);
533 c_tmp = qqrd2e/r;
534 }
535 }
536 }
537
538 drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
539 dftable[itablemax] = f_tmp - ftable[itablemax];
540 dctable[itablemax] = c_tmp - ctable[itablemax];
541 detable[itablemax] = e_tmp - etable[itablemax];
542 if (cut_respa) {
543 dvtable[itablemax] = v_tmp - vtable[itablemax];
544 dptable[itablemax] = p_tmp - ptable[itablemax];
545 }
546 }
547 }
548
549 /* ----------------------------------------------------------------------
550 setup force tables for dispersion used in compute routines
551 ------------------------------------------------------------------------- */
552
init_tables_disp(double cut_lj_global)553 void Pair::init_tables_disp(double cut_lj_global)
554 {
555 int masklo,maskhi;
556 double rsq;
557 double g_ewald_6 = force->kspace->g_ewald_6;
558 double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2;
559
560 tabinnerdispsq = tabinner_disp*tabinner_disp;
561 init_bitmap(tabinner_disp,cut_lj_global,ndisptablebits,
562 masklo,maskhi,ndispmask,ndispshiftbits);
563
564 int ntable = 1;
565 for (int i = 0; i < ndisptablebits; i++) ntable *= 2;
566
567 // linear lookup tables of length N = 2^ndisptablebits
568 // stored value = value at lower edge of bin
569 // d values = delta from lower edge to upper edge of bin
570
571 if (fdisptable) free_disp_tables();
572
573 memory->create(rdisptable,ntable,"pair:rdisptable");
574 memory->create(fdisptable,ntable,"pair:fdisptable");
575 memory->create(edisptable,ntable,"pair:edisptable");
576 memory->create(drdisptable,ntable,"pair:drdisptable");
577 memory->create(dfdisptable,ntable,"pair:dfdisptable");
578 memory->create(dedisptable,ntable,"pair:dedisptable");
579
580 union_int_float_t rsq_lookup;
581 union_int_float_t minrsq_lookup;
582 int itablemin;
583 minrsq_lookup.i = 0 << ndispshiftbits;
584 minrsq_lookup.i |= maskhi;
585
586 for (int i = 0; i < ntable; i++) {
587 rsq_lookup.i = i << ndispshiftbits;
588 rsq_lookup.i |= masklo;
589 if (rsq_lookup.f < tabinnerdispsq) {
590 rsq_lookup.i = i << ndispshiftbits;
591 rsq_lookup.i |= maskhi;
592 }
593 rsq = rsq_lookup.f;
594 double x2 = g2*rsq, a2 = 1.0/x2;
595 x2 = a2*exp(-x2);
596
597 rdisptable[i] = rsq_lookup.f;
598 fdisptable[i] = g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
599 edisptable[i] = g6*((a2+1.0)*a2+0.5)*x2;
600
601 minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
602 }
603
604 tabinnerdispsq = minrsq_lookup.f;
605
606 int ntablem1 = ntable - 1;
607
608 for (int i = 0; i < ntablem1; i++) {
609 drdisptable[i] = 1.0/(rdisptable[i+1] - rdisptable[i]);
610 dfdisptable[i] = fdisptable[i+1] - fdisptable[i];
611 dedisptable[i] = edisptable[i+1] - edisptable[i];
612 }
613
614 // get the delta values for the last table entries
615 // tables are connected periodically between 0 and ntablem1
616
617 drdisptable[ntablem1] = 1.0/(rdisptable[0] - rdisptable[ntablem1]);
618 dfdisptable[ntablem1] = fdisptable[0] - fdisptable[ntablem1];
619 dedisptable[ntablem1] = edisptable[0] - edisptable[ntablem1];
620
621 // get the correct delta values at itablemax
622 // smallest r is in bin itablemin
623 // largest r is in bin itablemax, which is itablemin-1,
624 // or ntablem1 if itablemin=0
625 // deltas at itablemax only needed if corresponding rsq < cut*cut
626 // if so, compute deltas between rsq and cut*cut
627
628 double f_tmp,e_tmp;
629 double cut_lj_globalsq;
630 itablemin = minrsq_lookup.i & ndispmask;
631 itablemin >>= ndispshiftbits;
632 int itablemax = itablemin - 1;
633 if (itablemin == 0) itablemax = ntablem1;
634 rsq_lookup.i = itablemax << ndispshiftbits;
635 rsq_lookup.i |= maskhi;
636
637 if (rsq_lookup.f < (cut_lj_globalsq = cut_lj_global * cut_lj_global)) {
638 rsq_lookup.f = cut_lj_globalsq;
639
640 double x2 = g2*rsq, a2 = 1.0/x2;
641 x2 = a2*exp(-x2);
642 f_tmp = g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
643 e_tmp = g6*((a2+1.0)*a2+0.5)*x2;
644
645 drdisptable[itablemax] = 1.0/(rsq_lookup.f - rdisptable[itablemax]);
646 dfdisptable[itablemax] = f_tmp - fdisptable[itablemax];
647 dedisptable[itablemax] = e_tmp - edisptable[itablemax];
648 }
649 }
650
651 /* ----------------------------------------------------------------------
652 free memory for tables used in Coulombic pair computations
653 ------------------------------------------------------------------------- */
654
free_tables()655 void Pair::free_tables()
656 {
657 memory->destroy(rtable);
658 memory->destroy(drtable);
659 memory->destroy(ftable);
660 memory->destroy(dftable);
661 memory->destroy(ctable);
662 memory->destroy(dctable);
663 memory->destroy(etable);
664 memory->destroy(detable);
665 memory->destroy(vtable);
666 memory->destroy(dvtable);
667 memory->destroy(ptable);
668 memory->destroy(dptable);
669 }
670
671 /* ----------------------------------------------------------------------
672 free memory for tables used in pair computations for dispersion
673 ------------------------------------------------------------------------- */
674
free_disp_tables()675 void Pair::free_disp_tables()
676 {
677 memory->destroy(rdisptable);
678 memory->destroy(drdisptable);
679 memory->destroy(fdisptable);
680 memory->destroy(dfdisptable);
681 memory->destroy(edisptable);
682 memory->destroy(dedisptable);
683 }
684 /* ----------------------------------------------------------------------
685 mixing of pair potential prefactors (epsilon)
686 ------------------------------------------------------------------------- */
687
mix_energy(double eps1,double eps2,double sig1,double sig2)688 double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2)
689 {
690 if (mix_flag == GEOMETRIC)
691 return sqrt(eps1*eps2);
692 else if (mix_flag == ARITHMETIC)
693 return sqrt(eps1*eps2);
694 else if (mix_flag == SIXTHPOWER)
695 return (2.0 * sqrt(eps1*eps2) *
696 pow(sig1,3.0) * pow(sig2,3.0) / (pow(sig1,6.0) + pow(sig2,6.0)));
697 else return 0.0;
698 }
699
700 /* ----------------------------------------------------------------------
701 mixing of pair potential distances (sigma, cutoff)
702 ------------------------------------------------------------------------- */
703
mix_distance(double sig1,double sig2)704 double Pair::mix_distance(double sig1, double sig2)
705 {
706 if (mix_flag == GEOMETRIC)
707 return sqrt(sig1*sig2);
708 else if (mix_flag == ARITHMETIC)
709 return (0.5 * (sig1+sig2));
710 else if (mix_flag == SIXTHPOWER)
711 return pow((0.5 * (pow(sig1,6.0) + pow(sig2,6.0))),1.0/6.0);
712 else return 0.0;
713 }
714
715 /* ---------------------------------------------------------------------- */
716
compute_dummy(int eflag,int vflag)717 void Pair::compute_dummy(int eflag, int vflag)
718 {
719 ev_init(eflag,vflag);
720 }
721
722 /* ---------------------------------------------------------------------- */
723
read_restart(FILE *)724 void Pair::read_restart(FILE *)
725 {
726 if (comm->me == 0)
727 error->warning(FLERR,"Pair style restartinfo set but has no restart support");
728 }
729
730 /* ---------------------------------------------------------------------- */
731
write_restart(FILE *)732 void Pair::write_restart(FILE *)
733 {
734 if (comm->me == 0)
735 error->warning(FLERR,"Pair style restartinfo set but has no restart support");
736 }
737
738 /* -------------------------------------------------------------------
739 register a callback to a compute, so it can compute and accumulate
740 additional properties during the pair computation from within
741 Pair::ev_tally(). ensure each compute instance is registered only once
742 ---------------------------------------------------------------------- */
743
add_tally_callback(Compute * ptr)744 void Pair::add_tally_callback(Compute *ptr)
745 {
746 if (lmp->kokkos)
747 error->all(FLERR,"Cannot yet use compute tally with Kokkos");
748
749 int i,found=-1;
750
751 for (i=0; i < num_tally_compute; ++i) {
752 if (list_tally_compute[i] == ptr)
753 found = i;
754 }
755
756 if (found < 0) {
757 found = num_tally_compute;
758 ++num_tally_compute;
759 void *p = memory->srealloc((void *)list_tally_compute, sizeof(Compute *) * num_tally_compute,
760 "pair:list_tally_compute");
761 list_tally_compute = (Compute **) p;
762 list_tally_compute[num_tally_compute-1] = ptr;
763 }
764 }
765
766 /* -------------------------------------------------------------------
767 unregister a callback to a fix for additional pairwise tallying
768 ---------------------------------------------------------------------- */
769
del_tally_callback(Compute * ptr)770 void Pair::del_tally_callback(Compute *ptr)
771 {
772 int i,found=-1;
773
774 for (i=0; i < num_tally_compute; ++i) {
775 if (list_tally_compute[i] == ptr)
776 found = i;
777 }
778
779 if (found < 0)
780 return;
781
782 // compact the list of active computes
783 --num_tally_compute;
784 for (i=found; i < num_tally_compute; ++i) {
785 list_tally_compute[i] = list_tally_compute[i+1];
786 }
787 }
788
789 /* -------------------------------------------------------------------
790 build element to atom type mapping for manybody potentials
791 also clear and reset setflag[][] array and check missing entries
792 ---------------------------------------------------------------------- */
793
map_element2type(int narg,char ** arg,bool update_setflag)794 void Pair::map_element2type(int narg, char **arg, bool update_setflag)
795 {
796 int i,j;
797 const int ntypes = atom->ntypes;
798
799 // read args that map atom types to elements in potential file
800 // map[i] = which element the Ith atom type is, -1 if "NULL"
801 // nelements = # of unique elements
802 // elements = list of element names
803
804 if (narg != ntypes)
805 error->all(FLERR,"Incorrect args for pair coefficients");
806
807 if (elements) {
808 for (i = 0; i < nelements; i++) delete[] elements[i];
809 delete[] elements;
810 }
811 elements = new char*[ntypes];
812 for (i = 0; i < ntypes; i++) elements[i] = nullptr;
813
814 nelements = 0;
815 map[0] = -1;
816 for (i = 1; i <= narg; i++) {
817 std::string entry = arg[i-1];
818 if (entry == "NULL") {
819 map[i] = -1;
820 continue;
821 }
822 for (j = 0; j < nelements; j++)
823 if (entry == elements[j]) break;
824 map[i] = j;
825 if (j == nelements) {
826 elements[j] = utils::strdup(entry);
827 nelements++;
828 }
829 }
830
831 // if requested, clear setflag[i][j] and set it for type pairs
832 // where both are mapped to elements in map.
833
834 if (update_setflag) {
835
836 int count = 0;
837 for (i = 1; i <= ntypes; i++) {
838 for (j = i; j <= ntypes; j++) {
839 setflag[i][j] = 0;
840 if ((map[i] >= 0) && (map[j] >= 0)) {
841 setflag[i][j] = 1;
842 count++;
843 }
844 }
845 }
846
847 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
848 }
849 }
850
851 /* ----------------------------------------------------------------------
852 setup for energy, virial computation
853 see integrate::ev_set() for bitwise settings of eflag/vflag
854 set the following flags, values are otherwise set to 0:
855 eflag_global != 0 if ENERGY_GLOBAL bit of eflag set
856 eflag_atom != 0 if ENERGY_ATOM bit of eflag set
857 eflag_either != 0 if eflag_global or eflag_atom is set
858 vflag_global != 0 if VIRIAL_PAIR bit of vflag set, OR
859 if VIRIAL_FDOTR bit of vflag is set but no_virial_fdotr = 1
860 vflag_fdotr != 0 if VIRIAL_FDOTR bit of vflag set and no_virial_fdotr = 0
861 vflag_atom != 0 if VIRIAL_ATOM bit of vflag set, OR
862 if VIRIAL_CENTROID bit of vflag set
863 and centroidstressflag != CENTROID_AVAIL
864 cvflag_atom != 0 if VIRIAL_CENTROID bit of vflag set
865 and centroidstressflag = CENTROID_AVAIL
866 vflag_either != 0 if any of vflag_global, vflag_atom, cvflag_atom is set
867 evflag != 0 if eflag_either or vflag_either is set
868 centroidstressflag is set by the pair style to one of these values:
869 CENTROID_SAME = same as two-body stress
870 CENTROID_AVAIL = different and implemented
871 CENTROID_NOTAVAIL = different but not yet implemented
872 ------------------------------------------------------------------------- */
873
ev_setup(int eflag,int vflag,int alloc)874 void Pair::ev_setup(int eflag, int vflag, int alloc)
875 {
876 int i,n;
877
878 eflag_either = eflag;
879 eflag_global = eflag & ENERGY_GLOBAL;
880 eflag_atom = eflag & ENERGY_ATOM;
881
882 vflag_global = vflag & VIRIAL_PAIR;
883 if (vflag & VIRIAL_FDOTR && no_virial_fdotr_compute == 1) vflag_global = 1;
884 vflag_fdotr = 0;
885 if (vflag & VIRIAL_FDOTR && no_virial_fdotr_compute == 0) vflag_fdotr = 1;
886 vflag_atom = vflag & VIRIAL_ATOM;
887 if (vflag & VIRIAL_CENTROID && centroidstressflag != CENTROID_AVAIL) vflag_atom = 1;
888 cvflag_atom = 0;
889 if (vflag & VIRIAL_CENTROID && centroidstressflag == CENTROID_AVAIL) cvflag_atom = 1;
890 vflag_either = vflag_global || vflag_atom || cvflag_atom;
891
892 evflag = eflag_either || vflag_either;
893
894 // reallocate per-atom arrays if necessary
895
896 if (eflag_atom && atom->nmax > maxeatom) {
897 maxeatom = atom->nmax;
898 if (alloc) {
899 memory->destroy(eatom);
900 memory->create(eatom,comm->nthreads*maxeatom,"pair:eatom");
901 }
902 }
903 if (vflag_atom && atom->nmax > maxvatom) {
904 maxvatom = atom->nmax;
905 if (alloc) {
906 memory->destroy(vatom);
907 memory->create(vatom,comm->nthreads*maxvatom,6,"pair:vatom");
908 }
909 }
910 if (cvflag_atom && atom->nmax > maxcvatom) {
911 maxcvatom = atom->nmax;
912 if (alloc) {
913 memory->destroy(cvatom);
914 memory->create(cvatom,comm->nthreads*maxcvatom,9,"pair:cvatom");
915 }
916 }
917
918 // zero accumulators
919 // use force->newton instead of newton_pair
920 // b/c some bonds/dihedrals call pair::ev_tally with pairwise info
921
922 if (eflag_global) eng_vdwl = eng_coul = 0.0;
923 if (vflag_global || vflag_fdotr) for (i = 0; i < 6; i++) virial[i] = 0.0;
924 if (eflag_atom && alloc) {
925 n = atom->nlocal;
926 if (force->newton) n += atom->nghost;
927 for (i = 0; i < n; i++) eatom[i] = 0.0;
928 }
929 if (vflag_atom && alloc) {
930 n = atom->nlocal;
931 if (force->newton) n += atom->nghost;
932 for (i = 0; i < n; i++) {
933 vatom[i][0] = 0.0;
934 vatom[i][1] = 0.0;
935 vatom[i][2] = 0.0;
936 vatom[i][3] = 0.0;
937 vatom[i][4] = 0.0;
938 vatom[i][5] = 0.0;
939 }
940 }
941 if (cvflag_atom && alloc) {
942 n = atom->nlocal;
943 if (force->newton) n += atom->nghost;
944 for (i = 0; i < n; i++) {
945 cvatom[i][0] = 0.0;
946 cvatom[i][1] = 0.0;
947 cvatom[i][2] = 0.0;
948 cvatom[i][3] = 0.0;
949 cvatom[i][4] = 0.0;
950 cvatom[i][5] = 0.0;
951 cvatom[i][6] = 0.0;
952 cvatom[i][7] = 0.0;
953 cvatom[i][8] = 0.0;
954 cvatom[i][9] = 0.0;
955 }
956 }
957
958 // run ev_setup option for TALLY computes
959
960 if (num_tally_compute > 0) {
961 for (int k=0; k < num_tally_compute; ++k) {
962 Compute *c = list_tally_compute[k];
963 c->pair_setup_callback(eflag,vflag);
964 }
965 }
966 }
967
968 /* ----------------------------------------------------------------------
969 set all flags to zero for energy, virial computation
970 called by some complicated many-body potentials that use individual flags
971 to insure no holdover of flags from previous timestep
972 ------------------------------------------------------------------------- */
973
ev_unset()974 void Pair::ev_unset()
975 {
976 evflag = 0;
977
978 eflag_either = 0;
979 eflag_global = 0;
980 eflag_atom = 0;
981
982 vflag_either = 0;
983 vflag_global = 0;
984 vflag_atom = 0;
985 cvflag_atom = 0;
986 vflag_fdotr = 0;
987 }
988
989 /* ----------------------------------------------------------------------
990 tally eng_vdwl and virial into global or per-atom accumulators
991 need i < nlocal test since called by bond_quartic and dihedral_charmm
992 ------------------------------------------------------------------------- */
993
ev_tally(int i,int j,int nlocal,int newton_pair,double evdwl,double ecoul,double fpair,double delx,double dely,double delz)994 void Pair::ev_tally(int i, int j, int nlocal, int newton_pair,
995 double evdwl, double ecoul, double fpair,
996 double delx, double dely, double delz)
997 {
998 double evdwlhalf,ecoulhalf,epairhalf,v[6];
999
1000 if (eflag_either) {
1001 if (eflag_global) {
1002 if (newton_pair) {
1003 eng_vdwl += evdwl;
1004 eng_coul += ecoul;
1005 } else {
1006 evdwlhalf = 0.5*evdwl;
1007 ecoulhalf = 0.5*ecoul;
1008 if (i < nlocal) {
1009 eng_vdwl += evdwlhalf;
1010 eng_coul += ecoulhalf;
1011 }
1012 if (j < nlocal) {
1013 eng_vdwl += evdwlhalf;
1014 eng_coul += ecoulhalf;
1015 }
1016 }
1017 }
1018 if (eflag_atom) {
1019 epairhalf = 0.5 * (evdwl + ecoul);
1020 if (newton_pair || i < nlocal) eatom[i] += epairhalf;
1021 if (newton_pair || j < nlocal) eatom[j] += epairhalf;
1022 }
1023 }
1024
1025 if (vflag_either) {
1026 v[0] = delx*delx*fpair;
1027 v[1] = dely*dely*fpair;
1028 v[2] = delz*delz*fpair;
1029 v[3] = delx*dely*fpair;
1030 v[4] = delx*delz*fpair;
1031 v[5] = dely*delz*fpair;
1032
1033 if (vflag_global) {
1034 if (newton_pair) {
1035 virial[0] += v[0];
1036 virial[1] += v[1];
1037 virial[2] += v[2];
1038 virial[3] += v[3];
1039 virial[4] += v[4];
1040 virial[5] += v[5];
1041 } else {
1042 if (i < nlocal) {
1043 virial[0] += 0.5*v[0];
1044 virial[1] += 0.5*v[1];
1045 virial[2] += 0.5*v[2];
1046 virial[3] += 0.5*v[3];
1047 virial[4] += 0.5*v[4];
1048 virial[5] += 0.5*v[5];
1049 }
1050 if (j < nlocal) {
1051 virial[0] += 0.5*v[0];
1052 virial[1] += 0.5*v[1];
1053 virial[2] += 0.5*v[2];
1054 virial[3] += 0.5*v[3];
1055 virial[4] += 0.5*v[4];
1056 virial[5] += 0.5*v[5];
1057 }
1058 }
1059 }
1060
1061 if (vflag_atom) {
1062 if (newton_pair || i < nlocal) {
1063 vatom[i][0] += 0.5*v[0];
1064 vatom[i][1] += 0.5*v[1];
1065 vatom[i][2] += 0.5*v[2];
1066 vatom[i][3] += 0.5*v[3];
1067 vatom[i][4] += 0.5*v[4];
1068 vatom[i][5] += 0.5*v[5];
1069 }
1070 if (newton_pair || j < nlocal) {
1071 vatom[j][0] += 0.5*v[0];
1072 vatom[j][1] += 0.5*v[1];
1073 vatom[j][2] += 0.5*v[2];
1074 vatom[j][3] += 0.5*v[3];
1075 vatom[j][4] += 0.5*v[4];
1076 vatom[j][5] += 0.5*v[5];
1077 }
1078 }
1079 }
1080
1081 if (num_tally_compute > 0) {
1082 for (int k=0; k < num_tally_compute; ++k) {
1083 Compute *c = list_tally_compute[k];
1084 c->pair_tally_callback(i, j, nlocal, newton_pair,
1085 evdwl, ecoul, fpair, delx, dely, delz);
1086 }
1087 }
1088 }
1089
1090 /* ----------------------------------------------------------------------
1091 tally eng_vdwl and virial into global or per-atom accumulators
1092 can use this version with full neighbor lists
1093 ------------------------------------------------------------------------- */
1094
ev_tally_full(int i,double evdwl,double ecoul,double fpair,double delx,double dely,double delz)1095 void Pair::ev_tally_full(int i, double evdwl, double ecoul, double fpair,
1096 double delx, double dely, double delz)
1097 {
1098 double v[6];
1099
1100 if (eflag_either) {
1101 if (eflag_global) {
1102 eng_vdwl += 0.5*evdwl;
1103 eng_coul += 0.5*ecoul;
1104 }
1105 if (eflag_atom) eatom[i] += 0.5 * (evdwl + ecoul);
1106 }
1107
1108 if (vflag_either) {
1109 v[0] = 0.5*delx*delx*fpair;
1110 v[1] = 0.5*dely*dely*fpair;
1111 v[2] = 0.5*delz*delz*fpair;
1112 v[3] = 0.5*delx*dely*fpair;
1113 v[4] = 0.5*delx*delz*fpair;
1114 v[5] = 0.5*dely*delz*fpair;
1115
1116 if (vflag_global) {
1117 virial[0] += v[0];
1118 virial[1] += v[1];
1119 virial[2] += v[2];
1120 virial[3] += v[3];
1121 virial[4] += v[4];
1122 virial[5] += v[5];
1123 }
1124
1125 if (vflag_atom) {
1126 vatom[i][0] += v[0];
1127 vatom[i][1] += v[1];
1128 vatom[i][2] += v[2];
1129 vatom[i][3] += v[3];
1130 vatom[i][4] += v[4];
1131 vatom[i][5] += v[5];
1132 }
1133 }
1134 }
1135
1136 /* ----------------------------------------------------------------------
1137 tally eng_vdwl and virial into global or per-atom accumulators
1138 for virial, have delx,dely,delz and fx,fy,fz
1139 ------------------------------------------------------------------------- */
1140
ev_tally_xyz(int i,int j,int nlocal,int newton_pair,double evdwl,double ecoul,double fx,double fy,double fz,double delx,double dely,double delz)1141 void Pair::ev_tally_xyz(int i, int j, int nlocal, int newton_pair,
1142 double evdwl, double ecoul,
1143 double fx, double fy, double fz,
1144 double delx, double dely, double delz)
1145 {
1146 double evdwlhalf,ecoulhalf,epairhalf,v[6];
1147
1148 if (eflag_either) {
1149 if (eflag_global) {
1150 if (newton_pair) {
1151 eng_vdwl += evdwl;
1152 eng_coul += ecoul;
1153 } else {
1154 evdwlhalf = 0.5*evdwl;
1155 ecoulhalf = 0.5*ecoul;
1156 if (i < nlocal) {
1157 eng_vdwl += evdwlhalf;
1158 eng_coul += ecoulhalf;
1159 }
1160 if (j < nlocal) {
1161 eng_vdwl += evdwlhalf;
1162 eng_coul += ecoulhalf;
1163 }
1164 }
1165 }
1166 if (eflag_atom) {
1167 epairhalf = 0.5 * (evdwl + ecoul);
1168 if (newton_pair || i < nlocal) eatom[i] += epairhalf;
1169 if (newton_pair || j < nlocal) eatom[j] += epairhalf;
1170 }
1171 }
1172
1173 if (vflag_either) {
1174 v[0] = delx*fx;
1175 v[1] = dely*fy;
1176 v[2] = delz*fz;
1177 v[3] = delx*fy;
1178 v[4] = delx*fz;
1179 v[5] = dely*fz;
1180
1181 if (vflag_global) {
1182 if (newton_pair) {
1183 virial[0] += v[0];
1184 virial[1] += v[1];
1185 virial[2] += v[2];
1186 virial[3] += v[3];
1187 virial[4] += v[4];
1188 virial[5] += v[5];
1189 } else {
1190 if (i < nlocal) {
1191 virial[0] += 0.5*v[0];
1192 virial[1] += 0.5*v[1];
1193 virial[2] += 0.5*v[2];
1194 virial[3] += 0.5*v[3];
1195 virial[4] += 0.5*v[4];
1196 virial[5] += 0.5*v[5];
1197 }
1198 if (j < nlocal) {
1199 virial[0] += 0.5*v[0];
1200 virial[1] += 0.5*v[1];
1201 virial[2] += 0.5*v[2];
1202 virial[3] += 0.5*v[3];
1203 virial[4] += 0.5*v[4];
1204 virial[5] += 0.5*v[5];
1205 }
1206 }
1207 }
1208
1209 if (vflag_atom) {
1210 if (newton_pair || i < nlocal) {
1211 vatom[i][0] += 0.5*v[0];
1212 vatom[i][1] += 0.5*v[1];
1213 vatom[i][2] += 0.5*v[2];
1214 vatom[i][3] += 0.5*v[3];
1215 vatom[i][4] += 0.5*v[4];
1216 vatom[i][5] += 0.5*v[5];
1217 }
1218 if (newton_pair || j < nlocal) {
1219 vatom[j][0] += 0.5*v[0];
1220 vatom[j][1] += 0.5*v[1];
1221 vatom[j][2] += 0.5*v[2];
1222 vatom[j][3] += 0.5*v[3];
1223 vatom[j][4] += 0.5*v[4];
1224 vatom[j][5] += 0.5*v[5];
1225 }
1226 }
1227 }
1228 }
1229
1230 /* ----------------------------------------------------------------------
1231 tally eng_vdwl and virial into global or per-atom accumulators
1232 for virial, have delx,dely,delz and fx,fy,fz
1233 called when using full neighbor lists
1234 ------------------------------------------------------------------------- */
1235
ev_tally_xyz_full(int i,double evdwl,double ecoul,double fx,double fy,double fz,double delx,double dely,double delz)1236 void Pair::ev_tally_xyz_full(int i, double evdwl, double ecoul,
1237 double fx, double fy, double fz,
1238 double delx, double dely, double delz)
1239 {
1240 double evdwlhalf,ecoulhalf,epairhalf,v[6];
1241
1242 if (eflag_either) {
1243 if (eflag_global) {
1244 evdwlhalf = 0.5*evdwl;
1245 ecoulhalf = 0.5*ecoul;
1246 eng_vdwl += evdwlhalf;
1247 eng_coul += ecoulhalf;
1248 }
1249 if (eflag_atom) {
1250 epairhalf = 0.5 * (evdwl + ecoul);
1251 eatom[i] += epairhalf;
1252 }
1253 }
1254
1255 if (vflag_either) {
1256 v[0] = 0.5*delx*fx;
1257 v[1] = 0.5*dely*fy;
1258 v[2] = 0.5*delz*fz;
1259 v[3] = 0.5*delx*fy;
1260 v[4] = 0.5*delx*fz;
1261 v[5] = 0.5*dely*fz;
1262
1263 if (vflag_global) {
1264 virial[0] += v[0];
1265 virial[1] += v[1];
1266 virial[2] += v[2];
1267 virial[3] += v[3];
1268 virial[4] += v[4];
1269 virial[5] += v[5];
1270 }
1271
1272 if (vflag_atom) {
1273 vatom[i][0] += v[0];
1274 vatom[i][1] += v[1];
1275 vatom[i][2] += v[2];
1276 vatom[i][3] += v[3];
1277 vatom[i][4] += v[4];
1278 vatom[i][5] += v[5];
1279 }
1280 }
1281 }
1282
1283 /* ----------------------------------------------------------------------
1284 tally eng_vdwl and virial into global or per-atom accumulators
1285 called by SW and hbond potentials, newton_pair is always on
1286 virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk
1287 ------------------------------------------------------------------------- */
1288
ev_tally3(int i,int j,int k,double evdwl,double ecoul,double * fj,double * fk,double * drji,double * drki)1289 void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul,
1290 double *fj, double *fk, double *drji, double *drki)
1291 {
1292 double epairthird,v[6];
1293
1294 if (eflag_either) {
1295 if (eflag_global) {
1296 eng_vdwl += evdwl;
1297 eng_coul += ecoul;
1298 }
1299 if (eflag_atom) {
1300 epairthird = THIRD * (evdwl + ecoul);
1301 eatom[i] += epairthird;
1302 eatom[j] += epairthird;
1303 eatom[k] += epairthird;
1304 }
1305 }
1306
1307 if (vflag_either) {
1308 v[0] = drji[0]*fj[0] + drki[0]*fk[0];
1309 v[1] = drji[1]*fj[1] + drki[1]*fk[1];
1310 v[2] = drji[2]*fj[2] + drki[2]*fk[2];
1311 v[3] = drji[0]*fj[1] + drki[0]*fk[1];
1312 v[4] = drji[0]*fj[2] + drki[0]*fk[2];
1313 v[5] = drji[1]*fj[2] + drki[1]*fk[2];
1314
1315 if (vflag_global) {
1316 virial[0] += v[0];
1317 virial[1] += v[1];
1318 virial[2] += v[2];
1319 virial[3] += v[3];
1320 virial[4] += v[4];
1321 virial[5] += v[5];
1322 }
1323
1324 if (vflag_atom) {
1325 vatom[i][0] += THIRD*v[0]; vatom[i][1] += THIRD*v[1];
1326 vatom[i][2] += THIRD*v[2]; vatom[i][3] += THIRD*v[3];
1327 vatom[i][4] += THIRD*v[4]; vatom[i][5] += THIRD*v[5];
1328
1329 vatom[j][0] += THIRD*v[0]; vatom[j][1] += THIRD*v[1];
1330 vatom[j][2] += THIRD*v[2]; vatom[j][3] += THIRD*v[3];
1331 vatom[j][4] += THIRD*v[4]; vatom[j][5] += THIRD*v[5];
1332
1333 vatom[k][0] += THIRD*v[0]; vatom[k][1] += THIRD*v[1];
1334 vatom[k][2] += THIRD*v[2]; vatom[k][3] += THIRD*v[3];
1335 vatom[k][4] += THIRD*v[4]; vatom[k][5] += THIRD*v[5];
1336 }
1337 }
1338 }
1339
1340 /* ----------------------------------------------------------------------
1341 tally eng_vdwl and virial into global or per-atom accumulators
1342 called by AIREBO potential, newton_pair is always on
1343 ------------------------------------------------------------------------- */
1344
ev_tally4(int i,int j,int k,int m,double evdwl,double * fi,double * fj,double * fk,double * drim,double * drjm,double * drkm)1345 void Pair::ev_tally4(int i, int j, int k, int m, double evdwl,
1346 double *fi, double *fj, double *fk,
1347 double *drim, double *drjm, double *drkm)
1348 {
1349 double epairfourth,v[6];
1350
1351 if (eflag_either) {
1352 if (eflag_global) eng_vdwl += evdwl;
1353 if (eflag_atom) {
1354 epairfourth = 0.25 * evdwl;
1355 eatom[i] += epairfourth;
1356 eatom[j] += epairfourth;
1357 eatom[k] += epairfourth;
1358 eatom[m] += epairfourth;
1359 }
1360 }
1361
1362 if (vflag_either) {
1363 v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
1364 v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
1365 v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
1366 v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
1367 v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
1368 v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
1369
1370 if (vflag_global) {
1371 virial[0] += v[0];
1372 virial[1] += v[1];
1373 virial[2] += v[2];
1374 virial[3] += v[3];
1375 virial[4] += v[4];
1376 virial[5] += v[5];
1377 }
1378
1379 if (vflag_atom) {
1380 v[0] *= 0.25;
1381 v[1] *= 0.25;
1382 v[2] *= 0.25;
1383 v[3] *= 0.25;
1384 v[4] *= 0.25;
1385 v[5] *= 0.25;
1386
1387 vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
1388 vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
1389 vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
1390 vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
1391 vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2];
1392 vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5];
1393 vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2];
1394 vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5];
1395 }
1396 }
1397 }
1398
1399 /* ----------------------------------------------------------------------
1400 tally ecoul and virial into each of atoms in list
1401 called by TIP4P potential, newton_pair is always on
1402 weight assignments by alpha, so contribution is all to O atom as alpha -> 0.0
1403 key = 0 if neither atom = water O
1404 key = 1 if first atom = water O
1405 key = 2 if second atom = water O
1406 key = 3 if both atoms = water O
1407 ------------------------------------------------------------------------- */
1408
ev_tally_tip4p(int key,int * list,double * v,double ecoul,double alpha)1409 void Pair::ev_tally_tip4p(int key, int *list, double *v,
1410 double ecoul, double alpha)
1411 {
1412 int i;
1413
1414 if (eflag_either) {
1415 if (eflag_global) eng_coul += ecoul;
1416 if (eflag_atom) {
1417 if (key == 0) {
1418 eatom[list[0]] += 0.5*ecoul;
1419 eatom[list[1]] += 0.5*ecoul;
1420 } else if (key == 1) {
1421 eatom[list[0]] += 0.5*ecoul*(1-alpha);
1422 eatom[list[1]] += 0.25*ecoul*alpha;
1423 eatom[list[2]] += 0.25*ecoul*alpha;
1424 eatom[list[3]] += 0.5*ecoul;
1425 } else if (key == 2) {
1426 eatom[list[0]] += 0.5*ecoul;
1427 eatom[list[1]] += 0.5*ecoul*(1-alpha);
1428 eatom[list[2]] += 0.25*ecoul*alpha;
1429 eatom[list[3]] += 0.25*ecoul*alpha;
1430 } else {
1431 eatom[list[0]] += 0.5*ecoul*(1-alpha);
1432 eatom[list[1]] += 0.25*ecoul*alpha;
1433 eatom[list[2]] += 0.25*ecoul*alpha;
1434 eatom[list[3]] += 0.5*ecoul*(1-alpha);
1435 eatom[list[4]] += 0.25*ecoul*alpha;
1436 eatom[list[5]] += 0.25*ecoul*alpha;
1437 }
1438 }
1439 }
1440
1441 if (vflag_either) {
1442 if (vflag_global) {
1443 virial[0] += v[0];
1444 virial[1] += v[1];
1445 virial[2] += v[2];
1446 virial[3] += v[3];
1447 virial[4] += v[4];
1448 virial[5] += v[5];
1449 }
1450
1451 if (vflag_atom) {
1452 if (key == 0) {
1453 for (i = 0; i <= 5; i++) {
1454 vatom[list[0]][i] += 0.5*v[i];
1455 vatom[list[1]][i] += 0.5*v[i];
1456 }
1457 } else if (key == 1) {
1458 for (i = 0; i <= 5; i++) {
1459 vatom[list[0]][i] += 0.5*v[i]*(1-alpha);
1460 vatom[list[1]][i] += 0.25*v[i]*alpha;
1461 vatom[list[2]][i] += 0.25*v[i]*alpha;
1462 vatom[list[3]][i] += 0.5*v[i];
1463 }
1464 } else if (key == 2) {
1465 for (i = 0; i <= 5; i++) {
1466 vatom[list[0]][i] += 0.5*v[i];
1467 vatom[list[1]][i] += 0.5*v[i]*(1-alpha);
1468 vatom[list[2]][i] += 0.25*v[i]*alpha;
1469 vatom[list[3]][i] += 0.25*v[i]*alpha;
1470 }
1471 } else {
1472 for (i = 0; i <= 5; i++) {
1473 vatom[list[0]][i] += 0.5*v[i]*(1-alpha);
1474 vatom[list[1]][i] += 0.25*v[i]*alpha;
1475 vatom[list[2]][i] += 0.25*v[i]*alpha;
1476 vatom[list[3]][i] += 0.5*v[i]*(1-alpha);
1477 vatom[list[4]][i] += 0.25*v[i]*alpha;
1478 vatom[list[5]][i] += 0.25*v[i]*alpha;
1479 }
1480 }
1481 }
1482 }
1483 }
1484
1485 /* ----------------------------------------------------------------------
1486 tally virial into global or per-atom accumulators
1487 called by ReaxFF potential, newton_pair is always on
1488 fi is magnitude of force on atom i, deli is the direction
1489 note that the other atom (j) is not updated, due to newton on
1490 ------------------------------------------------------------------------- */
1491
v_tally2_newton(int i,double * fi,double * deli)1492 void Pair::v_tally2_newton(int i, double *fi, double *deli)
1493 {
1494 double v[6];
1495
1496 v[0] = deli[0]*fi[0];
1497 v[1] = deli[1]*fi[1];
1498 v[2] = deli[2]*fi[2];
1499 v[3] = deli[0]*fi[1];
1500 v[4] = deli[0]*fi[2];
1501 v[5] = deli[1]*fi[2];
1502
1503 if (vflag_global) {
1504 virial[0] += v[0];
1505 virial[1] += v[1];
1506 virial[2] += v[2];
1507 virial[3] += v[3];
1508 virial[4] += v[4];
1509 virial[5] += v[5];
1510 }
1511
1512 if (vflag_atom) {
1513 vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
1514 vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
1515 }
1516 }
1517
1518 /* ----------------------------------------------------------------------
1519 tally virial into global or per-atom accumulators
1520 called by AIREBO potential, newton_pair is always on
1521 fpair is magnitude of force on atom I
1522 ------------------------------------------------------------------------- */
1523
v_tally2(int i,int j,double fpair,double * drij)1524 void Pair::v_tally2(int i, int j, double fpair, double *drij)
1525 {
1526 double v[6];
1527
1528 v[0] = drij[0]*drij[0]*fpair;
1529 v[1] = drij[1]*drij[1]*fpair;
1530 v[2] = drij[2]*drij[2]*fpair;
1531 v[3] = drij[0]*drij[1]*fpair;
1532 v[4] = drij[0]*drij[2]*fpair;
1533 v[5] = drij[1]*drij[2]*fpair;
1534
1535 if (vflag_global) {
1536 virial[0] += v[0];
1537 virial[1] += v[1];
1538 virial[2] += v[2];
1539 virial[3] += v[3];
1540 virial[4] += v[4];
1541 virial[5] += v[5];
1542 }
1543
1544 if (vflag_atom) {
1545 v[0] *= 0.5;
1546 v[1] *= 0.5;
1547 v[2] *= 0.5;
1548 v[3] *= 0.5;
1549 v[4] *= 0.5;
1550 v[5] *= 0.5;
1551 vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
1552 vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
1553 vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
1554 vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
1555 }
1556 }
1557
1558 /* ----------------------------------------------------------------------
1559 tally virial into per-atom accumulators
1560 called by AIREBO and Tersoff potentials, newton_pair is always on
1561 ------------------------------------------------------------------------- */
1562
v_tally3(int i,int j,int k,double * fi,double * fj,double * drik,double * drjk)1563 void Pair::v_tally3(int i, int j, int k, double *fi, double *fj, double *drik, double *drjk)
1564 {
1565 double v[6];
1566
1567 v[0] = (drik[0]*fi[0] + drjk[0]*fj[0]);
1568 v[1] = (drik[1]*fi[1] + drjk[1]*fj[1]);
1569 v[2] = (drik[2]*fi[2] + drjk[2]*fj[2]);
1570 v[3] = (drik[0]*fi[1] + drjk[0]*fj[1]);
1571 v[4] = (drik[0]*fi[2] + drjk[0]*fj[2]);
1572 v[5] = (drik[1]*fi[2] + drjk[1]*fj[2]);
1573
1574 if (vflag_global) {
1575 virial[0] += v[0];
1576 virial[1] += v[1];
1577 virial[2] += v[2];
1578 virial[3] += v[3];
1579 virial[4] += v[4];
1580 virial[5] += v[5];
1581 }
1582
1583 if (vflag_atom) {
1584 v[0] *= THIRD;
1585 v[1] *= THIRD;
1586 v[2] *= THIRD;
1587 v[3] *= THIRD;
1588 v[4] *= THIRD;
1589 v[5] *= THIRD;
1590 vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
1591 vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
1592 vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
1593 vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
1594 vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2];
1595 vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5];
1596 }
1597 }
1598
1599 /* ----------------------------------------------------------------------
1600 tally virial into global or per-atom accumulators
1601 called by AIREBO potential, Tersoff, ReaxFF potentials, newton_pair is always on
1602 ------------------------------------------------------------------------- */
1603
v_tally4(int i,int j,int k,int m,double * fi,double * fj,double * fk,double * drim,double * drjm,double * drkm)1604 void Pair::v_tally4(int i, int j, int k, int m,
1605 double *fi, double *fj, double *fk,
1606 double *drim, double *drjm, double *drkm)
1607 {
1608 double v[6];
1609
1610 v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
1611 v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
1612 v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
1613 v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
1614 v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
1615 v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
1616
1617 if (vflag_global) {
1618 virial[0] += v[0];
1619 virial[1] += v[1];
1620 virial[2] += v[2];
1621 virial[3] += v[3];
1622 virial[4] += v[4];
1623 virial[5] += v[5];
1624 }
1625
1626 if (vflag_atom) {
1627 v[0] *= 0.25;
1628 v[1] *= 0.25;
1629 v[2] *= 0.25;
1630 v[3] *= 0.25;
1631 v[4] *= 0.25;
1632 v[5] *= 0.25;
1633 vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2];
1634 vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5];
1635 vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2];
1636 vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5];
1637 vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2];
1638 vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5];
1639 vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2];
1640 vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5];
1641 }
1642 }
1643
1644 /* ----------------------------------------------------------------------
1645 tally virial into global or per-atom accumulators
1646 called by pair lubricate potential with 6 tensor components
1647 ------------------------------------------------------------------------- */
1648
v_tally_tensor(int i,int j,int nlocal,int newton_pair,double vxx,double vyy,double vzz,double vxy,double vxz,double vyz)1649 void Pair::v_tally_tensor(int i, int j, int nlocal, int newton_pair,
1650 double vxx, double vyy, double vzz,
1651 double vxy, double vxz, double vyz)
1652 {
1653 double v[6];
1654
1655 v[0] = vxx;
1656 v[1] = vyy;
1657 v[2] = vzz;
1658 v[3] = vxy;
1659 v[4] = vxz;
1660 v[5] = vyz;
1661
1662 if (vflag_global) {
1663 if (newton_pair) {
1664 virial[0] += v[0];
1665 virial[1] += v[1];
1666 virial[2] += v[2];
1667 virial[3] += v[3];
1668 virial[4] += v[4];
1669 virial[5] += v[5];
1670 } else {
1671 if (i < nlocal) {
1672 virial[0] += 0.5*v[0];
1673 virial[1] += 0.5*v[1];
1674 virial[2] += 0.5*v[2];
1675 virial[3] += 0.5*v[3];
1676 virial[4] += 0.5*v[4];
1677 virial[5] += 0.5*v[5];
1678 }
1679 if (j < nlocal) {
1680 virial[0] += 0.5*v[0];
1681 virial[1] += 0.5*v[1];
1682 virial[2] += 0.5*v[2];
1683 virial[3] += 0.5*v[3];
1684 virial[4] += 0.5*v[4];
1685 virial[5] += 0.5*v[5];
1686 }
1687 }
1688 }
1689
1690 if (vflag_atom) {
1691 if (newton_pair || i < nlocal) {
1692 vatom[i][0] += 0.5*v[0];
1693 vatom[i][1] += 0.5*v[1];
1694 vatom[i][2] += 0.5*v[2];
1695 vatom[i][3] += 0.5*v[3];
1696 vatom[i][4] += 0.5*v[4];
1697 vatom[i][5] += 0.5*v[5];
1698 }
1699 if (newton_pair || j < nlocal) {
1700 vatom[j][0] += 0.5*v[0];
1701 vatom[j][1] += 0.5*v[1];
1702 vatom[j][2] += 0.5*v[2];
1703 vatom[j][3] += 0.5*v[3];
1704 vatom[j][4] += 0.5*v[4];
1705 vatom[j][5] += 0.5*v[5];
1706 }
1707 }
1708 }
1709
1710 /* ----------------------------------------------------------------------
1711 compute global pair virial via summing F dot r over own & ghost atoms
1712 at this point, only pairwise forces have been accumulated in atom->f
1713 ------------------------------------------------------------------------- */
1714
virial_fdotr_compute()1715 void Pair::virial_fdotr_compute()
1716 {
1717 double **x = atom->x;
1718 double **f = atom->f;
1719
1720 // sum over force on all particles including ghosts
1721
1722 if (neighbor->includegroup == 0) {
1723 int nall = atom->nlocal + atom->nghost;
1724 for (int i = 0; i < nall; i++) {
1725 virial[0] += f[i][0]*x[i][0];
1726 virial[1] += f[i][1]*x[i][1];
1727 virial[2] += f[i][2]*x[i][2];
1728 virial[3] += f[i][1]*x[i][0];
1729 virial[4] += f[i][2]*x[i][0];
1730 virial[5] += f[i][2]*x[i][1];
1731 }
1732
1733 // neighbor includegroup flag is set
1734 // sum over force on initial nfirst particles and ghosts
1735
1736 } else {
1737 int nall = atom->nfirst;
1738 for (int i = 0; i < nall; i++) {
1739 virial[0] += f[i][0]*x[i][0];
1740 virial[1] += f[i][1]*x[i][1];
1741 virial[2] += f[i][2]*x[i][2];
1742 virial[3] += f[i][1]*x[i][0];
1743 virial[4] += f[i][2]*x[i][0];
1744 virial[5] += f[i][2]*x[i][1];
1745 }
1746
1747 nall = atom->nlocal + atom->nghost;
1748 for (int i = atom->nlocal; i < nall; i++) {
1749 virial[0] += f[i][0]*x[i][0];
1750 virial[1] += f[i][1]*x[i][1];
1751 virial[2] += f[i][2]*x[i][2];
1752 virial[3] += f[i][1]*x[i][0];
1753 virial[4] += f[i][2]*x[i][0];
1754 virial[5] += f[i][2]*x[i][1];
1755 }
1756 }
1757
1758 // prevent multiple calls to update the virial
1759 // when a hybrid pair style uses both a gpu and non-gpu pair style
1760 // or when respa is used with gpu pair styles
1761
1762 vflag_fdotr = 0;
1763 }
1764
1765 /* ----------------------------------------------------------------------
1766 write a table of pair potential energy/force vs distance to a file
1767 ------------------------------------------------------------------------- */
1768
write_file(int narg,char ** arg)1769 void Pair::write_file(int narg, char **arg)
1770 {
1771 if (narg != 8 && narg != 10) error->all(FLERR,"Illegal pair_write command");
1772 if (single_enable == 0)
1773 error->all(FLERR,"Pair style does not support pair_write");
1774
1775 // parse arguments
1776
1777 int itype = utils::inumeric(FLERR,arg[0],false,lmp);
1778 int jtype = utils::inumeric(FLERR,arg[1],false,lmp);
1779 if (itype < 1 || itype > atom->ntypes || jtype < 1 || jtype > atom->ntypes)
1780 error->all(FLERR,"Invalid atom types in pair_write command");
1781
1782 int n = utils::inumeric(FLERR,arg[2],false,lmp);
1783
1784 int style = NONE;
1785 if (strcmp(arg[3],"r") == 0) style = RLINEAR;
1786 else if (strcmp(arg[3],"rsq") == 0) style = RSQ;
1787 else if (strcmp(arg[3],"bitmap") == 0) style = BMP;
1788 else error->all(FLERR,"Invalid style in pair_write command");
1789
1790 double inner = utils::numeric(FLERR,arg[4],false,lmp);
1791 double outer = utils::numeric(FLERR,arg[5],false,lmp);
1792 if (inner <= 0.0 || inner >= outer)
1793 error->all(FLERR,"Invalid cutoffs in pair_write command");
1794
1795 // open file in append mode if exists
1796 // add line with DATE: and UNITS: tag when creating new file
1797 // print header in format used by pair_style table
1798
1799 FILE *fp = nullptr;
1800 if (comm->me == 0) {
1801 std::string table_file = arg[6];
1802
1803 // units sanity check:
1804 // - if this is the first time we write to this potential file,
1805 // write out a line with "DATE:" and "UNITS:" tags
1806 // - if the file already exists, print a message about appending
1807 // while printing the date and check that units are consistent.
1808 if (utils::file_is_readable(table_file)) {
1809 std::string units = utils::get_potential_units(table_file,"table");
1810 if (!units.empty() && (units != update->unit_style)) {
1811 error->one(FLERR,"Trying to append to a table file "
1812 "with UNITS: {} while units are {}",
1813 units, update->unit_style);
1814 }
1815 std::string date = utils::get_potential_date(table_file,"table");
1816 utils::logmesg(lmp,"Appending to table file {} with DATE: {}\n", table_file, date);
1817 fp = fopen(table_file.c_str(),"a");
1818 } else {
1819 utils::logmesg(lmp,"Creating table file {} with DATE: {}\n",
1820 table_file, utils::current_date());
1821 fp = fopen(table_file.c_str(),"w");
1822 if (fp) fmt::print(fp,"# DATE: {} UNITS: {} Created by pair_write\n",
1823 utils::current_date(), update->unit_style);
1824 }
1825 if (fp == nullptr)
1826 error->one(FLERR,"Cannot open pair_write file {}: {}",table_file, utils::getsyserror());
1827 fprintf(fp,"# Pair potential %s for atom types %d %d: i,r,energy,force\n",
1828 force->pair_style,itype,jtype);
1829 if (style == RLINEAR)
1830 fprintf(fp,"\n%s\nN %d R %.15g %.15g\n\n",arg[7],n,inner,outer);
1831 if (style == RSQ)
1832 fprintf(fp,"\n%s\nN %d RSQ %.15g %.15g\n\n",arg[7],n,inner,outer);
1833 }
1834
1835 // initialize potentials before evaluating pair potential
1836 // insures all pair coeffs are set and force constants
1837 // also initialize neighbor so that neighbor requests are processed
1838 // NOTE: might be safest to just do lmp->init()
1839
1840 force->init();
1841 neighbor->init();
1842
1843 // if pair style = any of EAM, swap in dummy fp vector
1844
1845 double eamfp[2];
1846 eamfp[0] = eamfp[1] = 0.0;
1847 double *eamfp_hold;
1848
1849 Pair *epair = force->pair_match("^eam",0);
1850 if (epair) epair->swap_eam(eamfp,&eamfp_hold);
1851 if ((comm->me == 0) && (epair))
1852 error->warning(FLERR,"EAM pair style. Table will not include embedding term");
1853
1854 // if atom style defines charge, swap in dummy q vec
1855
1856 double q[2];
1857 q[0] = q[1] = 1.0;
1858 if (narg == 10) {
1859 q[0] = utils::numeric(FLERR,arg[8],false,lmp);
1860 q[1] = utils::numeric(FLERR,arg[9],false,lmp);
1861 }
1862 double *q_hold;
1863
1864 if (atom->q) {
1865 q_hold = atom->q;
1866 atom->q = q;
1867 }
1868
1869 // evaluate energy and force at each of N distances
1870
1871 int masklo,maskhi,nmask,nshiftbits;
1872 if (style == BMP) {
1873 init_bitmap(inner,outer,n,masklo,maskhi,nmask,nshiftbits);
1874 int ntable = 1 << n;
1875 if (comm->me == 0)
1876 fprintf(fp,"\n%s\nN %d BITMAP %.15g %.15g\n\n",arg[7],ntable,inner,outer);
1877 n = ntable;
1878 }
1879
1880 double r,e,f,rsq;
1881 union_int_float_t rsq_lookup;
1882
1883 for (int i = 0; i < n; i++) {
1884 if (style == RLINEAR) {
1885 r = inner + (outer-inner) * i/(n-1);
1886 rsq = r*r;
1887 } else if (style == RSQ) {
1888 rsq = inner*inner + (outer*outer - inner*inner) * i/(n-1);
1889 r = sqrt(rsq);
1890 } else if (style == BMP) {
1891 rsq_lookup.i = i << nshiftbits;
1892 rsq_lookup.i |= masklo;
1893 if (rsq_lookup.f < inner*inner) {
1894 rsq_lookup.i = i << nshiftbits;
1895 rsq_lookup.i |= maskhi;
1896 }
1897 rsq = rsq_lookup.f;
1898 r = sqrt(rsq);
1899 }
1900
1901 if (rsq < cutsq[itype][jtype]) {
1902 e = single(0,1,itype,jtype,rsq,1.0,1.0,f);
1903 f *= r;
1904 } else e = f = 0.0;
1905 if (comm->me == 0) fprintf(fp,"%d %.15g %.15g %.15g\n",i+1,r,e,f);
1906 }
1907
1908 // restore original vecs that were swapped in for
1909
1910 double *tmp;
1911 if (epair) epair->swap_eam(eamfp_hold,&tmp);
1912 if (atom->q) atom->q = q_hold;
1913
1914 if (comm->me == 0) fclose(fp);
1915 }
1916
1917 /* ----------------------------------------------------------------------
1918 define bitmap parameters based on inner and outer cutoffs
1919 ------------------------------------------------------------------------- */
1920
init_bitmap(double inner,double outer,int ntablebits,int & masklo,int & maskhi,int & nmask,int & nshiftbits)1921 void Pair::init_bitmap(double inner, double outer, int ntablebits,
1922 int &masklo, int &maskhi, int &nmask, int &nshiftbits)
1923 {
1924 if (sizeof(int) != sizeof(float))
1925 error->all(FLERR,"Bitmapped lookup tables require int/float be same size");
1926
1927 if (ntablebits > (int)sizeof(float)*CHAR_BIT)
1928 error->all(FLERR,"Too many total bits for bitmapped lookup table");
1929
1930 if (inner >= outer)
1931 error->warning(FLERR,"Table inner cutoff >= outer cutoff");
1932
1933 int nlowermin = 1;
1934 while (!((pow(double(2),(double)nlowermin) <= inner*inner) &&
1935 (pow(double(2),(double)nlowermin+1.0) > inner*inner))) {
1936 if (pow(double(2),(double)nlowermin) <= inner*inner) nlowermin++;
1937 else nlowermin--;
1938 }
1939
1940 int nexpbits = 0;
1941 double required_range = outer*outer / pow(double(2),(double)nlowermin);
1942 double available_range = 2.0;
1943
1944 while (available_range < required_range) {
1945 nexpbits++;
1946 available_range = pow(double(2),pow(double(2),(double)nexpbits));
1947 }
1948
1949 int nmantbits = ntablebits - nexpbits;
1950
1951 if (nexpbits > (int)sizeof(float)*CHAR_BIT - FLT_MANT_DIG)
1952 error->all(FLERR,"Too many exponent bits for lookup table");
1953 if (nmantbits+1 > FLT_MANT_DIG)
1954 error->all(FLERR,"Too many mantissa bits for lookup table");
1955 if (nmantbits < 3) error->all(FLERR,"Too few bits for lookup table");
1956
1957 nshiftbits = FLT_MANT_DIG - (nmantbits+1);
1958
1959 nmask = 1;
1960 for (int j = 0; j < ntablebits+nshiftbits; j++) nmask *= 2;
1961 nmask -= 1;
1962
1963 union_int_float_t rsq_lookup;
1964 rsq_lookup.f = outer*outer;
1965 maskhi = rsq_lookup.i & ~(nmask);
1966 rsq_lookup.f = inner*inner;
1967 masklo = rsq_lookup.i & ~(nmask);
1968 }
1969
1970 /* ---------------------------------------------------------------------- */
1971
hessian_twobody(double fforce,double dfac,double delr[3],double phiTensor[6])1972 void Pair::hessian_twobody(double fforce, double dfac, double delr[3], double phiTensor[6]) {
1973 int m = 0;
1974 for (int k=0; k<3; k++) {
1975 phiTensor[m] = fforce;
1976 for (int l=k; l<3; l++) {
1977 if (l>k) phiTensor[m] = 0;
1978 phiTensor[m++] += delr[k]*delr[l] * dfac;
1979 }
1980 }
1981 }
1982 /* ---------------------------------------------------------------------- */
1983
memory_usage()1984 double Pair::memory_usage()
1985 {
1986 double bytes = (double)comm->nthreads*maxeatom * sizeof(double);
1987 bytes += (double)comm->nthreads*maxvatom*6 * sizeof(double);
1988 bytes += (double)comm->nthreads*maxcvatom*9 * sizeof(double);
1989 return bytes;
1990 }
1991
1992