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: Trung Nguyen (Northwestern)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_lj_expand_coul_long.h"
20
21 #include <cmath>
22 #include <cstring>
23 #include "atom.h"
24 #include "comm.h"
25 #include "force.h"
26 #include "kspace.h"
27 #include "update.h"
28 #include "respa.h"
29 #include "neighbor.h"
30 #include "neigh_list.h"
31 #include "neigh_request.h"
32 #include "math_const.h"
33 #include "memory.h"
34 #include "error.h"
35
36
37 using namespace LAMMPS_NS;
38 using namespace MathConst;
39
40 #define EWALD_F 1.12837917
41 #define EWALD_P 0.3275911
42 #define A1 0.254829592
43 #define A2 -0.284496736
44 #define A3 1.421413741
45 #define A4 -1.453152027
46 #define A5 1.061405429
47
48 /* ---------------------------------------------------------------------- */
49
PairLJExpandCoulLong(LAMMPS * lmp)50 PairLJExpandCoulLong::PairLJExpandCoulLong(LAMMPS *lmp) : Pair(lmp)
51 {
52 ewaldflag = pppmflag = 1;
53 respa_enable = 1;
54 writedata = 1;
55 ftable = nullptr;
56 qdist = 0.0;
57 cut_respa = nullptr;
58 }
59
60 /* ---------------------------------------------------------------------- */
61
~PairLJExpandCoulLong()62 PairLJExpandCoulLong::~PairLJExpandCoulLong()
63 {
64 if (allocated) {
65 memory->destroy(setflag);
66 memory->destroy(cutsq);
67
68 memory->destroy(cut_lj);
69 memory->destroy(cut_ljsq);
70 memory->destroy(epsilon);
71 memory->destroy(sigma);
72 memory->destroy(lj1);
73 memory->destroy(lj2);
74 memory->destroy(lj3);
75 memory->destroy(lj4);
76 memory->destroy(offset);
77 memory->destroy(shift);
78 }
79 if (ftable) free_tables();
80 }
81
82 /* ---------------------------------------------------------------------- */
83
compute(int eflag,int vflag)84 void PairLJExpandCoulLong::compute(int eflag, int vflag)
85 {
86 int i,ii,j,jj,inum,jnum,itype,jtype,itable;
87 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
88 double fraction,table;
89 double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
90 double grij,expm2,prefactor,t,erfc;
91 int *ilist,*jlist,*numneigh,**firstneigh;
92 double rsq,rshift,rshiftsq,rshift2inv;
93
94 evdwl = ecoul = 0.0;
95 ev_init(eflag,vflag);
96
97 double **x = atom->x;
98 double **f = atom->f;
99 double *q = atom->q;
100 int *type = atom->type;
101 int nlocal = atom->nlocal;
102 double *special_coul = force->special_coul;
103 double *special_lj = force->special_lj;
104 int newton_pair = force->newton_pair;
105 double qqrd2e = force->qqrd2e;
106
107 inum = list->inum;
108 ilist = list->ilist;
109 numneigh = list->numneigh;
110 firstneigh = list->firstneigh;
111
112 // loop over neighbors of my atoms
113
114 for (ii = 0; ii < inum; ii++) {
115 i = ilist[ii];
116 qtmp = q[i];
117 xtmp = x[i][0];
118 ytmp = x[i][1];
119 ztmp = x[i][2];
120 itype = type[i];
121 jlist = firstneigh[i];
122 jnum = numneigh[i];
123
124 for (jj = 0; jj < jnum; jj++) {
125 j = jlist[jj];
126 factor_lj = special_lj[sbmask(j)];
127 factor_coul = special_coul[sbmask(j)];
128 j &= NEIGHMASK;
129
130 delx = xtmp - x[j][0];
131 dely = ytmp - x[j][1];
132 delz = ztmp - x[j][2];
133 rsq = delx*delx + dely*dely + delz*delz;
134 jtype = type[j];
135
136 if (rsq < cutsq[itype][jtype]) {
137 r2inv = 1.0/rsq;
138
139 if (rsq < cut_coulsq) {
140 if (!ncoultablebits || rsq <= tabinnersq) {
141 r = sqrt(rsq);
142 grij = g_ewald * r;
143 expm2 = exp(-grij*grij);
144 t = 1.0 / (1.0 + EWALD_P*grij);
145 erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
146 prefactor = qqrd2e * qtmp*q[j]/r;
147 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
148 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
149 } else {
150 union_int_float_t rsq_lookup;
151 rsq_lookup.f = rsq;
152 itable = rsq_lookup.i & ncoulmask;
153 itable >>= ncoulshiftbits;
154 fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
155 table = ftable[itable] + fraction*dftable[itable];
156 forcecoul = qtmp*q[j] * table;
157 if (factor_coul < 1.0) {
158 table = ctable[itable] + fraction*dctable[itable];
159 prefactor = qtmp*q[j] * table;
160 forcecoul -= (1.0-factor_coul)*prefactor;
161 }
162 }
163 } else forcecoul = 0.0;
164
165 if (rsq < cut_ljsq[itype][jtype]) {
166 r = sqrt(rsq);
167 rshift = r - shift[itype][jtype];
168 rshiftsq = rshift*rshift;
169 rshift2inv = 1.0/rshiftsq;
170 r6inv = rshift2inv*rshift2inv*rshift2inv;
171 forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
172 forcelj *= factor_lj/rshift/r;
173 } else forcelj = 0.0;
174
175 fpair = forcecoul*r2inv + forcelj;
176
177 f[i][0] += delx*fpair;
178 f[i][1] += dely*fpair;
179 f[i][2] += delz*fpair;
180 if (newton_pair || j < nlocal) {
181 f[j][0] -= delx*fpair;
182 f[j][1] -= dely*fpair;
183 f[j][2] -= delz*fpair;
184 }
185
186 if (eflag) {
187 if (rsq < cut_coulsq) {
188 if (!ncoultablebits || rsq <= tabinnersq)
189 ecoul = prefactor*erfc;
190 else {
191 table = etable[itable] + fraction*detable[itable];
192 ecoul = qtmp*q[j] * table;
193 }
194 if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
195 } else ecoul = 0.0;
196
197 if (rsq < cut_ljsq[itype][jtype]) {
198 evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype];
199 evdwl *= factor_lj;
200 } else evdwl = 0.0;
201 }
202
203 if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,ecoul,fpair,delx,dely,delz);
204 }
205 }
206 }
207
208 if (vflag_fdotr) virial_fdotr_compute();
209 }
210
211 /* ---------------------------------------------------------------------- */
212
compute_inner()213 void PairLJExpandCoulLong::compute_inner()
214 {
215 int i,j,ii,jj,inum,jnum,itype,jtype;
216 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
217 double rsq,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
218 double rsw,r,rshift,rshiftsq,rshift2inv;
219 int *ilist,*jlist,*numneigh,**firstneigh;
220
221 double **x = atom->x;
222 double **f = atom->f;
223 double *q = atom->q;
224 int *type = atom->type;
225 int nlocal = atom->nlocal;
226 double *special_coul = force->special_coul;
227 double *special_lj = force->special_lj;
228 int newton_pair = force->newton_pair;
229 double qqrd2e = force->qqrd2e;
230
231 inum = list->inum_inner;
232 ilist = list->ilist_inner;
233 numneigh = list->numneigh_inner;
234 firstneigh = list->firstneigh_inner;
235
236 double cut_out_on = cut_respa[0];
237 double cut_out_off = cut_respa[1];
238
239 double cut_out_diff = cut_out_off - cut_out_on;
240 double cut_out_on_sq = cut_out_on*cut_out_on;
241 double cut_out_off_sq = cut_out_off*cut_out_off;
242
243 // loop over neighbors of my atoms
244
245 for (ii = 0; ii < inum; ii++) {
246 i = ilist[ii];
247 qtmp = q[i];
248 xtmp = x[i][0];
249 ytmp = x[i][1];
250 ztmp = x[i][2];
251 itype = type[i];
252 jlist = firstneigh[i];
253 jnum = numneigh[i];
254
255 for (jj = 0; jj < jnum; jj++) {
256 j = jlist[jj];
257 factor_lj = special_lj[sbmask(j)];
258 factor_coul = special_coul[sbmask(j)];
259 j &= NEIGHMASK;
260
261 delx = xtmp - x[j][0];
262 dely = ytmp - x[j][1];
263 delz = ztmp - x[j][2];
264 rsq = delx*delx + dely*dely + delz*delz;
265
266 if (rsq < cut_out_off_sq) {
267 r2inv = 1.0/rsq;
268 forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
269 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
270
271 jtype = type[j];
272 if (rsq < cut_ljsq[itype][jtype]) {
273 r = sqrt(rsq);
274 rshift = r - shift[itype][jtype];
275 rshiftsq = rshift*rshift;
276 rshift2inv = 1.0/rshiftsq;
277 r6inv = rshift2inv*rshift2inv*rshift2inv;
278 forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
279 forcelj *= factor_lj/rshift/r;
280 } else forcelj = 0.0;
281
282 fpair = forcecoul*r2inv + forcelj;
283 if (rsq > cut_out_on_sq) {
284 rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
285 fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
286 }
287
288 f[i][0] += delx*fpair;
289 f[i][1] += dely*fpair;
290 f[i][2] += delz*fpair;
291 if (newton_pair || j < nlocal) {
292 f[j][0] -= delx*fpair;
293 f[j][1] -= dely*fpair;
294 f[j][2] -= delz*fpair;
295 }
296 }
297 }
298 }
299 }
300
301 /* ---------------------------------------------------------------------- */
302
compute_middle()303 void PairLJExpandCoulLong::compute_middle()
304 {
305 int i,j,ii,jj,inum,jnum,itype,jtype;
306 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
307 double rsq,r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
308 double rsw,rshift,rshiftsq,rshift2inv;
309 int *ilist,*jlist,*numneigh,**firstneigh;
310
311 double **x = atom->x;
312 double **f = atom->f;
313 double *q = atom->q;
314 int *type = atom->type;
315 int nlocal = atom->nlocal;
316 double *special_coul = force->special_coul;
317 double *special_lj = force->special_lj;
318 int newton_pair = force->newton_pair;
319 double qqrd2e = force->qqrd2e;
320
321 inum = list->inum_middle;
322 ilist = list->ilist_middle;
323 numneigh = list->numneigh_middle;
324 firstneigh = list->firstneigh_middle;
325
326 double cut_in_off = cut_respa[0];
327 double cut_in_on = cut_respa[1];
328 double cut_out_on = cut_respa[2];
329 double cut_out_off = cut_respa[3];
330
331 double cut_in_diff = cut_in_on - cut_in_off;
332 double cut_out_diff = cut_out_off - cut_out_on;
333 double cut_in_off_sq = cut_in_off*cut_in_off;
334 double cut_in_on_sq = cut_in_on*cut_in_on;
335 double cut_out_on_sq = cut_out_on*cut_out_on;
336 double cut_out_off_sq = cut_out_off*cut_out_off;
337
338 // loop over neighbors of my atoms
339
340 for (ii = 0; ii < inum; ii++) {
341 i = ilist[ii];
342 qtmp = q[i];
343 xtmp = x[i][0];
344 ytmp = x[i][1];
345 ztmp = x[i][2];
346 itype = type[i];
347 jlist = firstneigh[i];
348 jnum = numneigh[i];
349
350 for (jj = 0; jj < jnum; jj++) {
351 j = jlist[jj];
352 factor_lj = special_lj[sbmask(j)];
353 factor_coul = special_coul[sbmask(j)];
354 j &= NEIGHMASK;
355
356 delx = xtmp - x[j][0];
357 dely = ytmp - x[j][1];
358 delz = ztmp - x[j][2];
359 rsq = delx*delx + dely*dely + delz*delz;
360
361 if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
362 r2inv = 1.0/rsq;
363 forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
364 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
365
366 jtype = type[j];
367 if (rsq < cut_ljsq[itype][jtype]) {
368 r = sqrt(rsq);
369 rshift = r - shift[itype][jtype];
370 rshiftsq = rshift*rshift;
371 rshift2inv = 1.0/rshiftsq;
372 r6inv = rshift2inv*rshift2inv*rshift2inv;
373 forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
374 forcelj *= factor_lj/rshift/r;
375 } else forcelj = 0.0;
376
377 fpair = forcecoul*r2inv + forcelj;
378 if (rsq < cut_in_on_sq) {
379 rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
380 fpair *= rsw*rsw*(3.0 - 2.0*rsw);
381 }
382 if (rsq > cut_out_on_sq) {
383 rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
384 fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
385 }
386
387 f[i][0] += delx*fpair;
388 f[i][1] += dely*fpair;
389 f[i][2] += delz*fpair;
390 if (newton_pair || j < nlocal) {
391 f[j][0] -= delx*fpair;
392 f[j][1] -= dely*fpair;
393 f[j][2] -= delz*fpair;
394 }
395 }
396 }
397 }
398 }
399
400 /* ---------------------------------------------------------------------- */
401
compute_outer(int eflag,int vflag)402 void PairLJExpandCoulLong::compute_outer(int eflag, int vflag)
403 {
404 int i,j,ii,jj,inum,jnum,itype,jtype,itable;
405 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
406 double fraction,table;
407 double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
408 double grij,expm2,prefactor,t,erfc;
409 double rsw,rsq,rshift,rshiftsq,rshift2inv;
410 int *ilist,*jlist,*numneigh,**firstneigh;
411
412 evdwl = ecoul = 0.0;
413 ev_init(eflag,vflag);
414
415 double **x = atom->x;
416 double **f = atom->f;
417 double *q = atom->q;
418 int *type = atom->type;
419 int nlocal = atom->nlocal;
420 double *special_coul = force->special_coul;
421 double *special_lj = force->special_lj;
422 int newton_pair = force->newton_pair;
423 double qqrd2e = force->qqrd2e;
424
425 inum = list->inum;
426 ilist = list->ilist;
427 numneigh = list->numneigh;
428 firstneigh = list->firstneigh;
429
430 double cut_in_off = cut_respa[2];
431 double cut_in_on = cut_respa[3];
432
433 double cut_in_diff = cut_in_on - cut_in_off;
434 double cut_in_off_sq = cut_in_off*cut_in_off;
435 double cut_in_on_sq = cut_in_on*cut_in_on;
436
437 // loop over neighbors of my atoms
438
439 for (ii = 0; ii < inum; ii++) {
440 i = ilist[ii];
441 qtmp = q[i];
442 xtmp = x[i][0];
443 ytmp = x[i][1];
444 ztmp = x[i][2];
445 itype = type[i];
446 jlist = firstneigh[i];
447 jnum = numneigh[i];
448
449 for (jj = 0; jj < jnum; jj++) {
450 j = jlist[jj];
451 factor_lj = special_lj[sbmask(j)];
452 factor_coul = special_coul[sbmask(j)];
453 j &= NEIGHMASK;
454
455 delx = xtmp - x[j][0];
456 dely = ytmp - x[j][1];
457 delz = ztmp - x[j][2];
458 rsq = delx*delx + dely*dely + delz*delz;
459 jtype = type[j];
460
461 if (rsq < cutsq[itype][jtype]) {
462 r2inv = 1.0/rsq;
463
464 if (rsq < cut_coulsq) {
465 if (!ncoultablebits || rsq <= tabinnersq) {
466 r = sqrt(rsq);
467 grij = g_ewald * r;
468 expm2 = exp(-grij*grij);
469 t = 1.0 / (1.0 + EWALD_P*grij);
470 erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
471 prefactor = qqrd2e * qtmp*q[j]/r;
472 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
473 if (rsq > cut_in_off_sq) {
474 if (rsq < cut_in_on_sq) {
475 rsw = (r - cut_in_off)/cut_in_diff;
476 forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw);
477 if (factor_coul < 1.0)
478 forcecoul -=
479 (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw);
480 } else {
481 forcecoul += prefactor;
482 if (factor_coul < 1.0)
483 forcecoul -= (1.0-factor_coul)*prefactor;
484 }
485 }
486 } else {
487 union_int_float_t rsq_lookup;
488 rsq_lookup.f = rsq;
489 itable = rsq_lookup.i & ncoulmask;
490 itable >>= ncoulshiftbits;
491 fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
492 table = ftable[itable] + fraction*dftable[itable];
493 forcecoul = qtmp*q[j] * table;
494 if (factor_coul < 1.0) {
495 table = ctable[itable] + fraction*dctable[itable];
496 prefactor = qtmp*q[j] * table;
497 forcecoul -= (1.0-factor_coul)*prefactor;
498 }
499 }
500 } else forcecoul = 0.0;
501
502 if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) {
503 r = sqrt(rsq);
504 rshift = r - shift[itype][jtype];
505 rshiftsq = rshift*rshift;
506 rshift2inv = 1.0/rshiftsq;
507 r6inv = rshift2inv*rshift2inv*rshift2inv;
508 forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
509 forcelj *= factor_lj/rshift/r;
510 if (rsq < cut_in_on_sq) {
511 rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
512 forcelj *= rsw*rsw*(3.0 - 2.0*rsw);
513 }
514 } else forcelj = 0.0;
515
516 fpair = forcecoul*r2inv + forcelj;
517
518 f[i][0] += delx*fpair;
519 f[i][1] += dely*fpair;
520 f[i][2] += delz*fpair;
521 if (newton_pair || j < nlocal) {
522 f[j][0] -= delx*fpair;
523 f[j][1] -= dely*fpair;
524 f[j][2] -= delz*fpair;
525 }
526
527 if (evflag) ecoul = evdwl = 0.0;
528 if (eflag) {
529 if (rsq < cut_coulsq) {
530 if (!ncoultablebits || rsq <= tabinnersq) {
531 ecoul = prefactor*erfc;
532 if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
533 } else {
534 table = etable[itable] + fraction*detable[itable];
535 ecoul = qtmp*q[j] * table;
536 if (factor_coul < 1.0) {
537 table = ptable[itable] + fraction*dptable[itable];
538 prefactor = qtmp*q[j] * table;
539 ecoul -= (1.0-factor_coul)*prefactor;
540 }
541 }
542 }
543
544 if (rsq < cut_ljsq[itype][jtype]) {
545 r = sqrt(rsq);
546 rshift = r - shift[itype][jtype];
547 rshiftsq = rshift*rshift;
548 rshift2inv = 1.0/rshiftsq;
549 r6inv = rshift2inv*rshift2inv*rshift2inv;
550 evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
551 offset[itype][jtype];
552 evdwl *= factor_lj;
553 }
554 }
555
556 if (vflag) {
557 if (rsq < cut_coulsq) {
558 if (!ncoultablebits || rsq <= tabinnersq) {
559 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
560 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
561 } else {
562 table = vtable[itable] + fraction*dvtable[itable];
563 forcecoul = qtmp*q[j] * table;
564 if (factor_coul < 1.0) {
565 table = ptable[itable] + fraction*dptable[itable];
566 prefactor = qtmp*q[j] * table;
567 forcecoul -= (1.0-factor_coul)*prefactor;
568 }
569 }
570 } else forcecoul = 0.0;
571
572 if (rsq < cut_ljsq[itype][jtype]) {
573 r = sqrt(rsq);
574 rshift = r - shift[itype][jtype];
575 rshiftsq = rshift*rshift;
576 rshift2inv = 1.0/rshiftsq;
577 r6inv = rshift2inv*rshift2inv*rshift2inv;
578 forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
579 forcelj *= factor_lj/rshift/r;
580 } else forcelj = 0.0;
581
582 fpair = forcecoul*r2inv + forcelj;
583 } else fpair = 0.0;
584
585 if (evflag) ev_tally(i,j,nlocal,newton_pair,
586 evdwl,ecoul,fpair,delx,dely,delz);
587 }
588 }
589 }
590 }
591
592 /* ----------------------------------------------------------------------
593 allocate all arrays
594 ------------------------------------------------------------------------- */
595
allocate()596 void PairLJExpandCoulLong::allocate()
597 {
598 allocated = 1;
599 int n = atom->ntypes;
600
601 memory->create(setflag,n+1,n+1,"pair:setflag");
602 for (int i = 1; i <= n; i++)
603 for (int j = i; j <= n; j++)
604 setflag[i][j] = 0;
605
606 memory->create(cutsq,n+1,n+1,"pair:cutsq");
607
608 memory->create(cut_lj,n+1,n+1,"pair:cut_lj");
609 memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq");
610 memory->create(epsilon,n+1,n+1,"pair:epsilon");
611 memory->create(sigma,n+1,n+1,"pair:sigma");
612 memory->create(lj1,n+1,n+1,"pair:lj1");
613 memory->create(lj2,n+1,n+1,"pair:lj2");
614 memory->create(lj3,n+1,n+1,"pair:lj3");
615 memory->create(lj4,n+1,n+1,"pair:lj4");
616 memory->create(offset,n+1,n+1,"pair:offset");
617 memory->create(shift,n+1,n+1,"pair:shift");
618 }
619
620 /* ----------------------------------------------------------------------
621 global settings
622 ------------------------------------------------------------------------- */
623
settings(int narg,char ** arg)624 void PairLJExpandCoulLong::settings(int narg, char **arg)
625 {
626 if (narg < 1 || narg > 2) error->all(FLERR,"Illegal pair_style command");
627
628 cut_lj_global = utils::numeric(FLERR,arg[0],false,lmp);
629 if (narg == 1) cut_coul = cut_lj_global;
630 else cut_coul = utils::numeric(FLERR,arg[1],false,lmp);
631
632 // reset cutoffs that have been explicitly set
633
634 if (allocated) {
635 int i,j;
636 for (i = 1; i <= atom->ntypes; i++)
637 for (j = i; j <= atom->ntypes; j++)
638 if (setflag[i][j]) cut_lj[i][j] = cut_lj_global;
639 }
640 }
641
642 /* ----------------------------------------------------------------------
643 set coeffs for one or more type pairs
644 ------------------------------------------------------------------------- */
645
coeff(int narg,char ** arg)646 void PairLJExpandCoulLong::coeff(int narg, char **arg)
647 {
648 if (narg < 5 || narg > 6)
649 error->all(FLERR,"Incorrect args for pair coefficients");
650 if (!allocated) allocate();
651
652 int ilo,ihi,jlo,jhi;
653 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
654 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
655
656 double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp);
657 double sigma_one = utils::numeric(FLERR,arg[3],false,lmp);
658 double shift_one = utils::numeric(FLERR,arg[4],false,lmp);
659
660 double cut_lj_one = cut_lj_global;
661 if (narg == 6) cut_lj_one = utils::numeric(FLERR,arg[5],false,lmp);
662
663 int count = 0;
664 for (int i = ilo; i <= ihi; i++) {
665 for (int j = MAX(jlo,i); j <= jhi; j++) {
666 epsilon[i][j] = epsilon_one;
667 sigma[i][j] = sigma_one;
668 shift[i][j] = shift_one;
669 cut_lj[i][j] = cut_lj_one;
670 setflag[i][j] = 1;
671 count++;
672 }
673 }
674
675 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
676 }
677
678 /* ----------------------------------------------------------------------
679 init specific to this pair style
680 ------------------------------------------------------------------------- */
681
init_style()682 void PairLJExpandCoulLong::init_style()
683 {
684 if (!atom->q_flag)
685 error->all(FLERR,"Pair style lj/cut/coul/long requires atom attribute q");
686
687 // request regular or rRESPA neighbor list
688
689 int irequest;
690 int respa = 0;
691
692 if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) {
693 if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
694 if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
695 }
696
697 irequest = neighbor->request(this,instance_me);
698
699 if (respa >= 1) {
700 neighbor->requests[irequest]->respaouter = 1;
701 neighbor->requests[irequest]->respainner = 1;
702 }
703 if (respa == 2) neighbor->requests[irequest]->respamiddle = 1;
704
705 cut_coulsq = cut_coul * cut_coul;
706
707 // set rRESPA cutoffs
708
709 if (utils::strmatch(update->integrate_style,"^respa") &&
710 ((Respa *) update->integrate)->level_inner >= 0)
711 cut_respa = ((Respa *) update->integrate)->cutoff;
712 else cut_respa = nullptr;
713
714 // insure use of KSpace long-range solver, set g_ewald
715
716 if (force->kspace == nullptr)
717 error->all(FLERR,"Pair style requires a KSpace style");
718 g_ewald = force->kspace->g_ewald;
719
720 // setup force tables
721
722 if (ncoultablebits) init_tables(cut_coul,cut_respa);
723 }
724
725 /* ----------------------------------------------------------------------
726 init for one type pair i,j and corresponding j,i
727 ------------------------------------------------------------------------- */
728
init_one(int i,int j)729 double PairLJExpandCoulLong::init_one(int i, int j)
730 {
731 if (setflag[i][j] == 0) {
732 epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
733 sigma[i][i],sigma[j][j]);
734 sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
735 cut_lj[i][j] = mix_distance(cut_lj[i][i],cut_lj[j][j]);
736 shift[i][j] = 0.5 * (shift[i][i] + shift[j][j]);
737 }
738
739 // include TIP4P qdist in full cutoff, qdist = 0.0 if not TIP4P
740
741 double cut = MAX(cut_lj[i][j]+shift[i][j],cut_coul+2.0*qdist);
742 cut_ljsq[i][j] = (cut_lj[i][j]+shift[i][j]) *(cut_lj[i][j]+shift[i][j]);
743
744 lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
745 lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
746 lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
747 lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
748
749 if (offset_flag && (cut_lj[i][j] > 0.0)) {
750 double ratio = sigma[i][j] / cut_lj[i][j];
751 offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
752 } else offset[i][j] = 0.0;
753
754 cut_ljsq[j][i] = cut_ljsq[i][j];
755 lj1[j][i] = lj1[i][j];
756 lj2[j][i] = lj2[i][j];
757 lj3[j][i] = lj3[i][j];
758 lj4[j][i] = lj4[i][j];
759 shift[j][i] = shift[i][j];
760 offset[j][i] = offset[i][j];
761
762 // check interior rRESPA cutoff
763
764 if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3])
765 error->all(FLERR,"Pair cutoff < Respa interior cutoff");
766
767 // compute I,J contribution to long-range tail correction
768 // count total # of atoms of type I and J via Allreduce
769
770 if (tail_flag) {
771 int *type = atom->type;
772 int nlocal = atom->nlocal;
773
774 double count[2],all[2];
775 count[0] = count[1] = 0.0;
776 for (int k = 0; k < nlocal; k++) {
777 if (type[k] == i) count[0] += 1.0;
778 if (type[k] == j) count[1] += 1.0;
779 }
780 MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
781
782 double sig2 = sigma[i][j]*sigma[i][j];
783 double sig6 = sig2*sig2*sig2;
784 double rc1 = cut_lj[i][j];
785 double rc2 = rc1*rc1;
786 double rc3 = rc2*rc1;
787 double rc6 = rc3*rc3;
788 double rc9 = rc3*rc6;
789 double shift1 = shift[i][j];
790 double shift2 = shift1*shift1;
791 double shift3 = shift2*shift1;
792 etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * sig6 *
793 ((1.0/9.0 + 2.0*shift1/(10.0*rc1) + shift2/(11.0*rc2))*sig6/rc9 -
794 (1.0/3.0 + 2.0*shift1/(4.0*rc1) + shift2/(5.0*rc2))/rc3);
795 ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * sig6 *
796 ((1.0/9.0 + 3.0*shift1/(10.0*rc1) +
797 3.0*shift2/(11.0*rc2) + shift3/(12.0*rc3))*2.0*sig6/rc9 -
798 (1.0/3.0 + 3.0*shift1/(4.0*rc1) +
799 3.0*shift2/(5.0*rc2) + shift3/(6.0*rc3))/rc3);
800 }
801
802 return cut;
803 }
804
805 /* ----------------------------------------------------------------------
806 proc 0 writes to restart file
807 ------------------------------------------------------------------------- */
808
write_restart(FILE * fp)809 void PairLJExpandCoulLong::write_restart(FILE *fp)
810 {
811 write_restart_settings(fp);
812
813 int i,j;
814 for (i = 1; i <= atom->ntypes; i++)
815 for (j = i; j <= atom->ntypes; j++) {
816 fwrite(&setflag[i][j],sizeof(int),1,fp);
817 if (setflag[i][j]) {
818 fwrite(&epsilon[i][j],sizeof(double),1,fp);
819 fwrite(&sigma[i][j],sizeof(double),1,fp);
820 fwrite(&shift[i][j],sizeof(double),1,fp);
821 fwrite(&cut_lj[i][j],sizeof(double),1,fp);
822 }
823 }
824 }
825
826 /* ----------------------------------------------------------------------
827 proc 0 reads from restart file, bcasts
828 ------------------------------------------------------------------------- */
829
read_restart(FILE * fp)830 void PairLJExpandCoulLong::read_restart(FILE *fp)
831 {
832 read_restart_settings(fp);
833
834 allocate();
835
836 int i,j;
837 int me = comm->me;
838 for (i = 1; i <= atom->ntypes; i++)
839 for (j = i; j <= atom->ntypes; j++) {
840 if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
841 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
842 if (setflag[i][j]) {
843 if (me == 0) {
844 utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error);
845 utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error);
846 utils::sfread(FLERR,&shift[i][j],sizeof(double),1,fp,nullptr,error);
847 utils::sfread(FLERR,&cut_lj[i][j],sizeof(double),1,fp,nullptr,error);
848 }
849 MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
850 MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
851 MPI_Bcast(&shift[i][j],1,MPI_DOUBLE,0,world);
852 MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);
853 }
854 }
855 }
856
857 /* ----------------------------------------------------------------------
858 proc 0 writes to restart file
859 ------------------------------------------------------------------------- */
860
write_restart_settings(FILE * fp)861 void PairLJExpandCoulLong::write_restart_settings(FILE *fp)
862 {
863 fwrite(&cut_lj_global,sizeof(double),1,fp);
864 fwrite(&cut_coul,sizeof(double),1,fp);
865 fwrite(&offset_flag,sizeof(int),1,fp);
866 fwrite(&mix_flag,sizeof(int),1,fp);
867 fwrite(&tail_flag,sizeof(int),1,fp);
868 fwrite(&ncoultablebits,sizeof(int),1,fp);
869 fwrite(&tabinner,sizeof(double),1,fp);
870 }
871
872 /* ----------------------------------------------------------------------
873 proc 0 reads from restart file, bcasts
874 ------------------------------------------------------------------------- */
875
read_restart_settings(FILE * fp)876 void PairLJExpandCoulLong::read_restart_settings(FILE *fp)
877 {
878 if (comm->me == 0) {
879 utils::sfread(FLERR,&cut_lj_global,sizeof(double),1,fp,nullptr,error);
880 utils::sfread(FLERR,&cut_coul,sizeof(double),1,fp,nullptr,error);
881 utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
882 utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
883 utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error);
884 utils::sfread(FLERR,&ncoultablebits,sizeof(int),1,fp,nullptr,error);
885 utils::sfread(FLERR,&tabinner,sizeof(double),1,fp,nullptr,error);
886 }
887 MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world);
888 MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
889 MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
890 MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
891 MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
892 MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
893 MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
894 }
895
896
897 /* ----------------------------------------------------------------------
898 proc 0 writes to data file
899 ------------------------------------------------------------------------- */
900
write_data(FILE * fp)901 void PairLJExpandCoulLong::write_data(FILE *fp)
902 {
903 for (int i = 1; i <= atom->ntypes; i++)
904 fprintf(fp,"%d %g %g %g\n",i,epsilon[i][i],sigma[i][i],shift[i][i]);
905 }
906
907 /* ----------------------------------------------------------------------
908 proc 0 writes all pairs to data file
909 ------------------------------------------------------------------------- */
910
write_data_all(FILE * fp)911 void PairLJExpandCoulLong::write_data_all(FILE *fp)
912 {
913 for (int i = 1; i <= atom->ntypes; i++)
914 for (int j = i; j <= atom->ntypes; j++)
915 fprintf(fp,"%d %d %g %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],shift[i][j],cut_lj[i][j]);
916 }
917
918 /* ---------------------------------------------------------------------- */
919
single(int i,int j,int itype,int jtype,double rsq,double factor_coul,double factor_lj,double & fforce)920 double PairLJExpandCoulLong::single(int i, int j, int itype, int jtype,
921 double rsq,
922 double factor_coul, double factor_lj,
923 double &fforce)
924 {
925 double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor;
926 double fraction,table,forcecoul,forcelj,phicoul,philj;
927 int itable;
928
929 r2inv = 1.0/rsq;
930 if (rsq < cut_coulsq) {
931 if (!ncoultablebits || rsq <= tabinnersq) {
932 r = sqrt(rsq);
933 grij = g_ewald * r;
934 expm2 = exp(-grij*grij);
935 t = 1.0 / (1.0 + EWALD_P*grij);
936 erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
937 prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
938 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
939 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
940 } else {
941 union_int_float_t rsq_lookup_single;
942 rsq_lookup_single.f = rsq;
943 itable = rsq_lookup_single.i & ncoulmask;
944 itable >>= ncoulshiftbits;
945 fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
946 table = ftable[itable] + fraction*dftable[itable];
947 forcecoul = atom->q[i]*atom->q[j] * table;
948 if (factor_coul < 1.0) {
949 table = ctable[itable] + fraction*dctable[itable];
950 prefactor = atom->q[i]*atom->q[j] * table;
951 forcecoul -= (1.0-factor_coul)*prefactor;
952 }
953 }
954 } else forcecoul = 0.0;
955
956 if (rsq < cut_ljsq[itype][jtype]) {
957 r = sqrt(rsq);
958 double rshift = r - shift[itype][jtype];
959 double rshiftsq = rshift*rshift;
960 double rshift2inv = 1.0/rshiftsq;
961 r6inv = rshift2inv*rshift2inv*rshift2inv;
962 forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
963 forcelj *= factor_lj/rshift/r;
964 } else forcelj = 0.0;
965
966 fforce = forcecoul*r2inv + forcelj;
967
968 double eng = 0.0;
969 if (rsq < cut_coulsq) {
970 if (!ncoultablebits || rsq <= tabinnersq)
971 phicoul = prefactor*erfc;
972 else {
973 table = etable[itable] + fraction*detable[itable];
974 phicoul = atom->q[i]*atom->q[j] * table;
975 }
976 if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
977 eng += phicoul;
978 }
979
980 if (rsq < cut_ljsq[itype][jtype]) {
981 philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype];
982 eng += factor_lj*philj;
983 }
984
985 return eng;
986 }
987
988 /* ---------------------------------------------------------------------- */
989
extract(const char * str,int & dim)990 void *PairLJExpandCoulLong::extract(const char *str, int &dim)
991 {
992 dim = 0;
993 if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
994 dim = 2;
995 if (strcmp(str,"epsilon") == 0) return (void *) epsilon;
996 if (strcmp(str,"sigma") == 0) return (void *) sigma;
997 if (strcmp(str,"delta") == 0) return (void *) shift;
998 return nullptr;
999 }
1000