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 #include "pair_lj_class2_coul_long_cs.h"
16 #include <cmath>
17 #include "atom.h"
18 #include "force.h"
19 #include "neigh_list.h"
20
21 using namespace LAMMPS_NS;
22
23 #define EWALD_F 1.12837917
24 #define EWALD_P 9.95473818e-1
25 #define B0 -0.1335096380159268
26 #define B1 -2.57839507e-1
27 #define B2 -1.37203639e-1
28 #define B3 -8.88822059e-3
29 #define B4 -5.80844129e-3
30 #define B5 1.14652755e-1
31
32 #define EPSILON 1.0e-20
33 #define EPS_EWALD 1.0e-6
34 #define EPS_EWALD_SQR 1.0e-12
35
36 /* ---------------------------------------------------------------------- */
37
PairLJClass2CoulLongCS(LAMMPS * lmp)38 PairLJClass2CoulLongCS::PairLJClass2CoulLongCS(LAMMPS *lmp) : PairLJClass2CoulLong(lmp)
39 {
40 ewaldflag = pppmflag = 1;
41 respa_enable = 0; // TODO: r-RESPA handling is inconsistent and thus disabled until fixed
42 single_enable = 0; // TODO: single function does not match compute
43 writedata = 1;
44 ftable = nullptr;
45 }
46
47 /* ---------------------------------------------------------------------- */
48
compute(int eflag,int vflag)49 void PairLJClass2CoulLongCS::compute(int eflag, int vflag)
50 {
51 int i,j,ii,jj,inum,jnum,itable,itype,jtype;
52 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
53 double fraction,table;
54 double rsq,r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj;
55 double grij,expm2,prefactor,t,erfc,u;
56 double factor_coul,factor_lj;
57 int *ilist,*jlist,*numneigh,**firstneigh;
58
59 evdwl = ecoul = 0.0;
60 ev_init(eflag,vflag);
61
62 double **x = atom->x;
63 double **f = atom->f;
64 double *q = atom->q;
65 int *type = atom->type;
66 int nlocal = atom->nlocal;
67 double *special_coul = force->special_coul;
68 double *special_lj = force->special_lj;
69 int newton_pair = force->newton_pair;
70 double qqrd2e = force->qqrd2e;
71
72 inum = list->inum;
73 ilist = list->ilist;
74 numneigh = list->numneigh;
75 firstneigh = list->firstneigh;
76
77 // loop over neighbors of my atoms
78
79 for (ii = 0; ii < inum; ii++) {
80 i = ilist[ii];
81 qtmp = q[i];
82 xtmp = x[i][0];
83 ytmp = x[i][1];
84 ztmp = x[i][2];
85 itype = type[i];
86 jlist = firstneigh[i];
87 jnum = numneigh[i];
88
89 for (jj = 0; jj < jnum; jj++) {
90 j = jlist[jj];
91 factor_lj = special_lj[sbmask(j)];
92 factor_coul = special_coul[sbmask(j)];
93 j &= NEIGHMASK;
94
95 delx = xtmp - x[j][0];
96 dely = ytmp - x[j][1];
97 delz = ztmp - x[j][2];
98 rsq = delx*delx + dely*dely + delz*delz;
99 jtype = type[j];
100
101 if (rsq < cutsq[itype][jtype]) {
102 rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond;
103 r2inv = 1.0/rsq;
104 if (rsq < cut_coulsq) {
105 if (!ncoultablebits || rsq <= tabinnersq) {
106 r = sqrt(rsq);
107 prefactor = qqrd2e * qtmp*q[j];
108 if (factor_coul < 1.0) {
109 // When bonded parts are being calculated a minimal distance (EPS_EWALD)
110 // has to be added to the prefactor and erfc in order to make the
111 // used approximation functions for the Ewald correction valid
112 grij = g_ewald * (r+EPS_EWALD);
113 expm2 = exp(-grij*grij);
114 t = 1.0 / (1.0 + EWALD_P*grij);
115 u = 1.0 - t;
116 erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
117 prefactor /= (r+EPS_EWALD);
118 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - (1.0-factor_coul));
119 // Additionally r2inv needs to be accordingly modified since the later
120 // scaling of the overall force shall be consistent
121 r2inv = 1.0/(rsq + EPS_EWALD_SQR);
122 } else {
123 grij = g_ewald * r;
124 expm2 = exp(-grij*grij);
125 t = 1.0 / (1.0 + EWALD_P*grij);
126 u = 1.0 - t;
127 erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
128 prefactor /= r;
129 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
130 }
131 } else {
132 union_int_float_t rsq_lookup;
133 rsq_lookup.f = rsq;
134 itable = rsq_lookup.i & ncoulmask;
135 itable >>= ncoulshiftbits;
136 fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
137 table = ftable[itable] + fraction*dftable[itable];
138 forcecoul = qtmp*q[j] * table;
139 if (factor_coul < 1.0) {
140 table = ctable[itable] + fraction*dctable[itable];
141 prefactor = qtmp*q[j] * table;
142 forcecoul -= (1.0-factor_coul)*prefactor;
143 }
144 }
145 } else forcecoul = 0.0;
146
147 if (rsq < cut_ljsq[itype][jtype]) {
148 rinv = sqrt(r2inv);
149 r3inv = r2inv*rinv;
150 r6inv = r3inv*r3inv;
151 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
152 } else forcelj = 0.0;
153
154 fpair = (forcecoul + factor_lj*forcelj) * r2inv;
155
156 f[i][0] += delx*fpair;
157 f[i][1] += dely*fpair;
158 f[i][2] += delz*fpair;
159 if (newton_pair || j < nlocal) {
160 f[j][0] -= delx*fpair;
161 f[j][1] -= dely*fpair;
162 f[j][2] -= delz*fpair;
163 }
164
165 if (eflag) {
166 if (rsq < cut_coulsq) {
167 if (!ncoultablebits || rsq <= tabinnersq)
168 ecoul = prefactor*erfc;
169 else {
170 table = etable[itable] + fraction*detable[itable];
171 ecoul = qtmp*q[j] * table;
172 }
173 if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
174 } else ecoul = 0.0;
175 if (rsq < cut_ljsq[itype][jtype]) {
176 evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
177 offset[itype][jtype];
178 evdwl *= factor_lj;
179 } else evdwl = 0.0;
180 }
181
182 if (evflag) ev_tally(i,j,nlocal,newton_pair,
183 evdwl,ecoul,fpair,delx,dely,delz);
184 }
185 }
186 }
187
188 if (vflag_fdotr) virial_fdotr_compute();
189 }
190
191 /* ---------------------------------------------------------------------- */
192
compute_inner()193 void PairLJClass2CoulLongCS::compute_inner()
194 {
195 int i,j,ii,jj,inum,jnum,itype,jtype;
196 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
197 double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
198 double rsw;
199 int *ilist,*jlist,*numneigh,**firstneigh;
200
201 double **x = atom->x;
202 double **f = atom->f;
203 double *q = atom->q;
204 int *type = atom->type;
205 int nlocal = atom->nlocal;
206 double *special_coul = force->special_coul;
207 double *special_lj = force->special_lj;
208 int newton_pair = force->newton_pair;
209 double qqrd2e = force->qqrd2e;
210
211 inum = list->inum_inner;
212 ilist = list->ilist_inner;
213 numneigh = list->numneigh_inner;
214 firstneigh = list->firstneigh_inner;
215
216 double cut_out_on = cut_respa[0];
217 double cut_out_off = cut_respa[1];
218
219 double cut_out_diff = cut_out_off - cut_out_on;
220 double cut_out_on_sq = cut_out_on*cut_out_on;
221 double cut_out_off_sq = cut_out_off*cut_out_off;
222
223 // loop over neighbors of my atoms
224
225 for (ii = 0; ii < inum; ii++) {
226 i = ilist[ii];
227 qtmp = q[i];
228 xtmp = x[i][0];
229 ytmp = x[i][1];
230 ztmp = x[i][2];
231 itype = type[i];
232 jlist = firstneigh[i];
233 jnum = numneigh[i];
234
235 for (jj = 0; jj < jnum; jj++) {
236 j = jlist[jj];
237 factor_lj = special_lj[sbmask(j)];
238 factor_coul = special_coul[sbmask(j)];
239 j &= NEIGHMASK;
240
241 delx = xtmp - x[j][0];
242 dely = ytmp - x[j][1];
243 delz = ztmp - x[j][2];
244 rsq = delx*delx + dely*dely + delz*delz;
245
246 if (rsq < cut_out_off_sq) {
247 rsq += EPSILON; // Add Epsilon for case: r = 0; Interaction must be removed by special bond;
248 r2inv = 1.0/rsq;
249 forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
250 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
251
252 jtype = type[j];
253 if (rsq < cut_ljsq[itype][jtype]) {
254 rinv = sqrt(r2inv);
255 r3inv = r2inv*rinv;
256 r6inv = r3inv*r3inv;
257 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
258 } else forcelj = 0.0;
259
260 fpair = (forcecoul + factor_lj*forcelj) * r2inv;
261 if (rsq > cut_out_on_sq) {
262 rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
263 fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
264 }
265
266 f[i][0] += delx*fpair;
267 f[i][1] += dely*fpair;
268 f[i][2] += delz*fpair;
269 if (newton_pair || j < nlocal) {
270 f[j][0] -= delx*fpair;
271 f[j][1] -= dely*fpair;
272 f[j][2] -= delz*fpair;
273 }
274 }
275 }
276 }
277 }
278
279 /* ---------------------------------------------------------------------- */
280
compute_middle()281 void PairLJClass2CoulLongCS::compute_middle()
282 {
283 int i,j,ii,jj,inum,jnum,itype,jtype;
284 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fpair;
285 double rsq,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
286 double rsw;
287 int *ilist,*jlist,*numneigh,**firstneigh;
288
289 double **x = atom->x;
290 double **f = atom->f;
291 double *q = atom->q;
292 int *type = atom->type;
293 int nlocal = atom->nlocal;
294 double *special_coul = force->special_coul;
295 double *special_lj = force->special_lj;
296 int newton_pair = force->newton_pair;
297 double qqrd2e = force->qqrd2e;
298
299 inum = list->inum_middle;
300 ilist = list->ilist_middle;
301 numneigh = list->numneigh_middle;
302 firstneigh = list->firstneigh_middle;
303
304 double cut_in_off = cut_respa[0];
305 double cut_in_on = cut_respa[1];
306 double cut_out_on = cut_respa[2];
307 double cut_out_off = cut_respa[3];
308
309 double cut_in_diff = cut_in_on - cut_in_off;
310 double cut_out_diff = cut_out_off - cut_out_on;
311 double cut_in_off_sq = cut_in_off*cut_in_off;
312 double cut_in_on_sq = cut_in_on*cut_in_on;
313 double cut_out_on_sq = cut_out_on*cut_out_on;
314 double cut_out_off_sq = cut_out_off*cut_out_off;
315
316 // loop over neighbors of my atoms
317
318 for (ii = 0; ii < inum; ii++) {
319 i = ilist[ii];
320 qtmp = q[i];
321 xtmp = x[i][0];
322 ytmp = x[i][1];
323 ztmp = x[i][2];
324 itype = type[i];
325 jlist = firstneigh[i];
326 jnum = numneigh[i];
327
328 for (jj = 0; jj < jnum; jj++) {
329 j = jlist[jj];
330 factor_lj = special_lj[sbmask(j)];
331 factor_coul = special_coul[sbmask(j)];
332 j &= NEIGHMASK;
333
334 delx = xtmp - x[j][0];
335 dely = ytmp - x[j][1];
336 delz = ztmp - x[j][2];
337 rsq = delx*delx + dely*dely + delz*delz;
338
339 if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
340 r2inv = 1.0/rsq;
341 forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
342 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*forcecoul;
343
344 jtype = type[j];
345 if (rsq < cut_ljsq[itype][jtype]) {
346 rinv = sqrt(r2inv);
347 r3inv = r2inv*rinv;
348 r6inv = r3inv*r3inv;
349 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
350 } else forcelj = 0.0;
351
352 fpair = (forcecoul + factor_lj*forcelj) * r2inv;
353 if (rsq < cut_in_on_sq) {
354 rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
355 fpair *= rsw*rsw*(3.0 - 2.0*rsw);
356 }
357 if (rsq > cut_out_on_sq) {
358 rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
359 fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
360 }
361
362 f[i][0] += delx*fpair;
363 f[i][1] += dely*fpair;
364 f[i][2] += delz*fpair;
365 if (newton_pair || j < nlocal) {
366 f[j][0] -= delx*fpair;
367 f[j][1] -= dely*fpair;
368 f[j][2] -= delz*fpair;
369 }
370 }
371 }
372 }
373 }
374
375 /* ---------------------------------------------------------------------- */
376
compute_outer(int eflag,int vflag)377 void PairLJClass2CoulLongCS::compute_outer(int eflag, int vflag)
378 {
379 int i,j,ii,jj,inum,jnum,itype,jtype,itable;
380 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
381 double fraction,table;
382 double r,rinv,r2inv,r3inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
383 double grij,expm2,prefactor,t,erfc,u;
384 double rsw;
385 int *ilist,*jlist,*numneigh,**firstneigh;
386 double rsq;
387
388 evdwl = ecoul = 0.0;
389 ev_init(eflag,vflag);
390
391 double **x = atom->x;
392 double **f = atom->f;
393 double *q = atom->q;
394 int *type = atom->type;
395 int nlocal = atom->nlocal;
396 double *special_coul = force->special_coul;
397 double *special_lj = force->special_lj;
398 int newton_pair = force->newton_pair;
399 double qqrd2e = force->qqrd2e;
400
401 inum = list->inum;
402 ilist = list->ilist;
403 numneigh = list->numneigh;
404 firstneigh = list->firstneigh;
405
406 double cut_in_off = cut_respa[2];
407 double cut_in_on = cut_respa[3];
408
409 double cut_in_diff = cut_in_on - cut_in_off;
410 double cut_in_off_sq = cut_in_off*cut_in_off;
411 double cut_in_on_sq = cut_in_on*cut_in_on;
412
413 // loop over neighbors of my atoms
414
415 for (ii = 0; ii < inum; ii++) {
416 i = ilist[ii];
417 qtmp = q[i];
418 xtmp = x[i][0];
419 ytmp = x[i][1];
420 ztmp = x[i][2];
421 itype = type[i];
422 jlist = firstneigh[i];
423 jnum = numneigh[i];
424
425 for (jj = 0; jj < jnum; jj++) {
426 j = jlist[jj];
427 factor_lj = special_lj[sbmask(j)];
428 factor_coul = special_coul[sbmask(j)];
429 j &= NEIGHMASK;
430
431 delx = xtmp - x[j][0];
432 dely = ytmp - x[j][1];
433 delz = ztmp - x[j][2];
434 rsq = delx*delx + dely*dely + delz*delz;
435 jtype = type[j];
436
437 if (rsq < cutsq[itype][jtype]) {
438 r2inv = 1.0/rsq;
439
440 if (rsq < cut_coulsq) {
441 if (!ncoultablebits || rsq <= tabinnersq) {
442 r = sqrt(rsq);
443 grij = g_ewald * r;
444 expm2 = exp(-grij*grij);
445 t = 1.0 / (1.0 + EWALD_P*grij);
446 u = 1. - t;
447 erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;
448 prefactor = qqrd2e * qtmp*q[j]/r;
449 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - 1.0);
450 if (rsq > cut_in_off_sq) {
451 if (rsq < cut_in_on_sq) {
452 rsw = (r - cut_in_off)/cut_in_diff;
453 forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw);
454 if (factor_coul < 1.0)
455 forcecoul -=
456 (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw);
457 } else {
458 forcecoul += prefactor;
459 if (factor_coul < 1.0)
460 forcecoul -= (1.0-factor_coul)*prefactor;
461 }
462 }
463 } else {
464 union_int_float_t rsq_lookup;
465 rsq_lookup.f = rsq;
466 itable = rsq_lookup.i & ncoulmask;
467 itable >>= ncoulshiftbits;
468 fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
469 table = ftable[itable] + fraction*dftable[itable];
470 forcecoul = qtmp*q[j] * table;
471 if (factor_coul < 1.0) {
472 table = ctable[itable] + fraction*dctable[itable];
473 prefactor = qtmp*q[j] * table;
474 forcecoul -= (1.0-factor_coul)*prefactor;
475 }
476 }
477 } else forcecoul = 0.0;
478
479 if (rsq < cut_ljsq[itype][jtype] && rsq > cut_in_off_sq) {
480 rinv = sqrt(r2inv);
481 r3inv = r2inv*rinv;
482 r6inv = r3inv*r3inv;
483 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
484 if (rsq < cut_in_on_sq) {
485 rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
486 forcelj *= rsw*rsw*(3.0 - 2.0*rsw);
487 }
488 } else forcelj = 0.0;
489
490 fpair = (forcecoul + forcelj) * r2inv;
491
492 f[i][0] += delx*fpair;
493 f[i][1] += dely*fpair;
494 f[i][2] += delz*fpair;
495 if (newton_pair || j < nlocal) {
496 f[j][0] -= delx*fpair;
497 f[j][1] -= dely*fpair;
498 f[j][2] -= delz*fpair;
499 }
500
501 if (eflag) {
502 if (rsq < cut_coulsq) {
503 if (!ncoultablebits || rsq <= tabinnersq) {
504 ecoul = prefactor*erfc;
505 if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
506 } else {
507 table = etable[itable] + fraction*detable[itable];
508 ecoul = qtmp*q[j] * table;
509 if (factor_coul < 1.0) {
510 table = ptable[itable] + fraction*dptable[itable];
511 prefactor = qtmp*q[j] * table;
512 ecoul -= (1.0-factor_coul)*prefactor;
513 }
514 }
515 } else ecoul = 0.0;
516
517 if (rsq < cut_ljsq[itype][jtype]) {
518 rinv = sqrt(r2inv);
519 r3inv = r2inv*rinv;
520 r6inv = r3inv*r3inv;
521 evdwl = r6inv*(lj3[itype][jtype]*r3inv-lj4[itype][jtype]) -
522 offset[itype][jtype];
523 evdwl *= factor_lj;
524 } else evdwl = 0.0;
525 }
526
527 if (vflag) {
528 if (rsq < cut_coulsq) {
529 if (!ncoultablebits || rsq <= tabinnersq) {
530 forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
531 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
532 } else {
533 table = vtable[itable] + fraction*dvtable[itable];
534 forcecoul = qtmp*q[j] * table;
535 if (factor_coul < 1.0) {
536 table = ptable[itable] + fraction*dptable[itable];
537 prefactor = qtmp*q[j] * table;
538 forcecoul -= (1.0-factor_coul)*prefactor;
539 }
540 }
541 } else forcecoul = 0.0;
542
543 if (rsq <= cut_in_off_sq) {
544 rinv = sqrt(r2inv);
545 r3inv = r2inv*rinv;
546 r6inv = r3inv*r3inv;
547 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
548 } else if (rsq <= cut_in_on_sq) {
549 rinv = sqrt(r2inv);
550 r3inv = r2inv*rinv;
551 r6inv = r3inv*r3inv;
552 forcelj = r6inv * (lj1[itype][jtype]*r3inv - lj2[itype][jtype]);
553 }
554 fpair = (forcecoul + factor_lj*forcelj) * r2inv;
555 }
556
557 if (evflag) ev_tally(i,j,nlocal,newton_pair,
558 evdwl,ecoul,fpair,delx,dely,delz);
559 }
560 }
561 }
562 }
563
564