1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/ Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13
14 /* ----------------------------------------------------------------------
15 Contributing author: Trung Nguyen (Northwestern)
16 ------------------------------------------------------------------------- */
17
18 #include "pair_lj_cut_coul_msm_dielectric.h"
19
20 #include "atom.h"
21 #include "atom_vec_dielectric.h"
22 #include "error.h"
23 #include "force.h"
24 #include "kspace.h"
25 #include "math_const.h"
26 #include "memory.h"
27 #include "neigh_list.h"
28 #include "neigh_request.h"
29 #include "neighbor.h"
30
31 #include <cmath>
32 #include <cstring>
33
34 using namespace LAMMPS_NS;
35 using namespace MathConst;
36
37 #define EPSILON 1e-6
38
39 /* ---------------------------------------------------------------------- */
40
PairLJCutCoulMSMDielectric(LAMMPS * lmp)41 PairLJCutCoulMSMDielectric::PairLJCutCoulMSMDielectric(LAMMPS *lmp) : PairLJCutCoulLong(lmp)
42 {
43 ewaldflag = pppmflag = 0;
44 msmflag = 1;
45 respa_enable = 0;
46 cut_respa = nullptr;
47
48 nmax = 0;
49 ftmp = nullptr;
50 efield = nullptr;
51 }
52
53 /* ---------------------------------------------------------------------- */
54
~PairLJCutCoulMSMDielectric()55 PairLJCutCoulMSMDielectric::~PairLJCutCoulMSMDielectric()
56 {
57 if (ftmp) memory->destroy(ftmp);
58 memory->destroy(efield);
59 }
60
61 /* ---------------------------------------------------------------------- */
62
compute(int eflag,int vflag)63 void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
64 {
65 int i, ii, j, jj, inum, jnum, itype, jtype, itable;
66 double qtmp, etmp, xtmp, ytmp, ztmp, delx, dely, delz, evdwl, ecoul, fpair;
67 double fpair_i, fpair_j;
68 double fraction, table;
69 double r, r2inv, r6inv, forcecoul, forcelj, factor_coul, factor_lj;
70 double egamma, fgamma, prefactor, prefactorE, efield_i;
71 int *ilist, *jlist, *numneigh, **firstneigh;
72 double rsq;
73 int eflag_old = eflag;
74
75 if (force->kspace->scalar_pressure_flag && vflag) {
76 if (vflag > 2)
77 error->all(FLERR,
78 "Must use 'kspace_modify pressure/scalar no' "
79 "to obtain per-atom virial with kspace_style MSM");
80
81 if (atom->nmax > nmax) {
82 if (ftmp) memory->destroy(ftmp);
83 nmax = atom->nmax;
84 memory->create(ftmp, nmax, 3, "pair:ftmp");
85 }
86 memset(&ftmp[0][0], 0, nmax * 3 * sizeof(double));
87
88 // must switch on global energy computation if not already on
89
90 if (eflag == 0 || eflag == 2) { eflag++; }
91 }
92
93 if (!efield || atom->nmax > nmax) {
94 if (efield) memory->destroy(efield);
95 nmax = atom->nmax;
96 memory->create(efield, nmax, 3, "pair:efield");
97 }
98
99 evdwl = ecoul = 0.0;
100 ev_init(eflag, vflag);
101
102 double **x = atom->x;
103 double **f = atom->f;
104 double *q = atom->q;
105 double *eps = atom->epsilon;
106 double **norm = atom->mu;
107 double *curvature = atom->curvature;
108 double *area = atom->area;
109 int *type = atom->type;
110 int nlocal = atom->nlocal;
111 double *special_coul = force->special_coul;
112 double *special_lj = force->special_lj;
113 int newton_pair = force->newton_pair;
114 double qqrd2e = force->qqrd2e;
115
116 inum = list->inum;
117 ilist = list->ilist;
118 numneigh = list->numneigh;
119 firstneigh = list->firstneigh;
120
121 // loop over neighbors of my atoms
122
123 for (ii = 0; ii < inum; ii++) {
124 i = ilist[ii];
125 qtmp = q[i];
126 xtmp = x[i][0];
127 ytmp = x[i][1];
128 ztmp = x[i][2];
129 etmp = eps[i];
130 itype = type[i];
131 jlist = firstneigh[i];
132 jnum = numneigh[i];
133
134 // self term Eq. (55) for I_{ii} and Eq. (52) and in Barros et al
135 double curvature_threshold = sqrt(area[i]);
136 if (curvature[i] < curvature_threshold) {
137 double sf = curvature[i] / (4.0 * MY_PIS * curvature_threshold) * area[i] * q[i];
138 efield[i][0] = sf * norm[i][0];
139 efield[i][1] = sf * norm[i][1];
140 efield[i][2] = sf * norm[i][2];
141 } else {
142 efield[i][0] = efield[i][1] = efield[i][2] = 0;
143 }
144
145 for (jj = 0; jj < jnum; jj++) {
146 j = jlist[jj];
147 factor_lj = special_lj[sbmask(j)];
148 factor_coul = special_coul[sbmask(j)];
149 j &= NEIGHMASK;
150
151 delx = xtmp - x[j][0];
152 dely = ytmp - x[j][1];
153 delz = ztmp - x[j][2];
154 rsq = delx * delx + dely * dely + delz * delz;
155 jtype = type[j];
156
157 if (rsq < cutsq[itype][jtype]) {
158 r2inv = 1.0 / rsq;
159
160 if (rsq < cut_coulsq && rsq > EPSILON) {
161 if (!ncoultablebits || rsq <= tabinnersq) {
162 r = sqrt(rsq);
163 prefactor = qqrd2e * qtmp * q[j] / r;
164 egamma = 1.0 - (r / cut_coul) * force->kspace->gamma(r / cut_coul);
165 fgamma = 1.0 + (rsq / cut_coulsq) * force->kspace->dgamma(r / cut_coul);
166 forcecoul = prefactor * fgamma;
167 if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
168
169 prefactorE = q[j] / r;
170 efield_i = prefactorE * fgamma;
171 if (factor_coul < 1.0) efield_i -= (1.0 - factor_coul) * prefactorE;
172
173 } else {
174 union_int_float_t rsq_lookup;
175 rsq_lookup.f = rsq;
176 itable = rsq_lookup.i & ncoulmask;
177 itable >>= ncoulshiftbits;
178 fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
179 table = ftable[itable] + fraction * dftable[itable];
180 forcecoul = qtmp * q[j] * table;
181 efield_i = q[j] * table / qqrd2e;
182 if (factor_coul < 1.0) {
183 table = ctable[itable] + fraction * dctable[itable];
184 prefactor = qtmp * q[j] * table;
185 forcecoul -= (1.0 - factor_coul) * prefactor;
186
187 prefactorE = q[j] * table / qqrd2e;
188 efield_i -= (1.0 - factor_coul) * prefactorE;
189 }
190 }
191 } else
192 forcecoul = 0.0;
193
194 if (rsq < cut_ljsq[itype][jtype]) {
195 r6inv = r2inv * r2inv * r2inv;
196 forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
197 } else
198 forcelj = 0.0;
199
200 if (!(force->kspace->scalar_pressure_flag && vflag)) {
201
202 fpair_i = (forcecoul * etmp + factor_lj * forcelj) * r2inv;
203 f[i][0] += delx * fpair_i;
204 f[i][1] += dely * fpair_i;
205 f[i][2] += delz * fpair_i;
206
207 efield_i *= (etmp * r2inv);
208 efield[i][0] += delx * efield_i;
209 efield[i][1] += dely * efield_i;
210 efield[i][2] += delz * efield_i;
211
212 if (newton_pair && j >= nlocal) {
213 fpair_j = (forcecoul * eps[j] + factor_lj * forcelj) * r2inv;
214 f[j][0] -= delx * fpair_j;
215 f[j][1] -= dely * fpair_j;
216 f[j][2] -= delz * fpair_j;
217 }
218 } else {
219
220 // separate LJ and Coulombic forces
221
222 fpair = (factor_lj * forcelj) * r2inv;
223
224 f[i][0] += delx * fpair;
225 f[i][1] += dely * fpair;
226 f[i][2] += delz * fpair;
227 if (newton_pair) {
228 f[j][0] -= delx * fpair;
229 f[j][1] -= dely * fpair;
230 f[j][2] -= delz * fpair;
231 }
232
233 fpair_i = (forcecoul * etmp) * r2inv;
234 ftmp[i][0] += delx * fpair_i;
235 ftmp[i][1] += dely * fpair_i;
236 ftmp[i][2] += delz * fpair_i;
237
238 efield_i *= (etmp * r2inv);
239 efield[i][0] += delx * efield_i;
240 efield[i][1] += dely * efield_i;
241 efield[i][2] += delz * efield_i;
242
243 if (newton_pair && j >= nlocal) {
244 fpair_j = (forcecoul * eps[j]) * r2inv;
245 ftmp[j][0] -= delx * fpair_j;
246 ftmp[j][1] -= dely * fpair_j;
247 ftmp[j][2] -= delz * fpair_j;
248 }
249 }
250
251 if (eflag) {
252 if (rsq < cut_coulsq) {
253 if (!ncoultablebits || rsq <= tabinnersq)
254 ecoul = prefactor * (etmp + eps[j]) * egamma;
255 else {
256 table = etable[itable] + fraction * detable[itable];
257 ecoul = qtmp * q[j] * (etmp + eps[j]) * table;
258 }
259 if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
260 } else
261 ecoul = 0.0;
262
263 if (eflag_old && rsq < cut_ljsq[itype][jtype]) {
264 evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
265 evdwl *= factor_lj;
266 } else
267 evdwl = 0.0;
268 }
269
270 if (evflag) ev_tally_full(i, evdwl, ecoul, fpair_i, delx, dely, delz);
271 }
272 }
273 }
274
275 if (vflag_fdotr) virial_fdotr_compute();
276
277 if (force->kspace->scalar_pressure_flag && vflag) {
278 for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul / 3.0;
279 for (int i = 0; i < nmax; i++) {
280 f[i][0] += ftmp[i][0];
281 f[i][1] += ftmp[i][1];
282 f[i][2] += ftmp[i][2];
283 }
284 }
285 }
286
287 /* ---------------------------------------------------------------------- */
288
single(int i,int j,int itype,int jtype,double rsq,double factor_coul,double factor_lj,double & fforce)289 double PairLJCutCoulMSMDielectric::single(int i, int j, int itype, int jtype, double rsq,
290 double factor_coul, double factor_lj, double &fforce)
291 {
292 double r2inv, r6inv, r, egamma, fgamma, prefactor;
293 double fraction, table, forcecoul, forcelj, phicoul, philj;
294 int itable;
295
296 r2inv = 1.0 / rsq;
297 if (rsq < cut_coulsq) {
298 if (!ncoultablebits || rsq <= tabinnersq) {
299 r = sqrt(rsq);
300 prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
301 egamma = 1.0 - (r / cut_coul) * force->kspace->gamma(r / cut_coul);
302 fgamma = 1.0 + (rsq / cut_coulsq) * force->kspace->dgamma(r / cut_coul);
303 forcecoul = prefactor * fgamma;
304 if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
305 } else {
306 union_int_float_t rsq_lookup_single;
307 rsq_lookup_single.f = rsq;
308 itable = rsq_lookup_single.i & ncoulmask;
309 itable >>= ncoulshiftbits;
310 fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
311 table = ftable[itable] + fraction * dftable[itable];
312 forcecoul = atom->q[i] * atom->q[j] * table;
313 if (factor_coul < 1.0) {
314 table = ctable[itable] + fraction * dctable[itable];
315 prefactor = atom->q[i] * atom->q[j] * table;
316 forcecoul -= (1.0 - factor_coul) * prefactor;
317 }
318 }
319 } else
320 forcecoul = 0.0;
321
322 if (rsq < cut_ljsq[itype][jtype]) {
323 r6inv = r2inv * r2inv * r2inv;
324 forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
325 } else
326 forcelj = 0.0;
327
328 fforce = (forcecoul + factor_lj * forcelj) * r2inv;
329
330 double eng = 0.0;
331 if (rsq < cut_coulsq) {
332 if (!ncoultablebits || rsq <= tabinnersq)
333 phicoul = prefactor * egamma;
334 else {
335 table = etable[itable] + fraction * detable[itable];
336 phicoul = atom->q[i] * atom->q[j] * table;
337 }
338 if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor;
339 eng += phicoul;
340 }
341
342 if (rsq < cut_ljsq[itype][jtype]) {
343 philj = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
344 eng += factor_lj * philj;
345 }
346
347 return eng;
348 }
349
350 /* ----------------------------------------------------------------------
351 init specific to this pair style
352 ------------------------------------------------------------------------- */
353
init_style()354 void PairLJCutCoulMSMDielectric::init_style()
355 {
356 avec = (AtomVecDielectric *) atom->style_match("dielectric");
357 if (!avec) error->all(FLERR, "Pair lj/cut/coul/msm/dielectric requires atom style dielectric");
358
359 int irequest = neighbor->request(this, instance_me);
360 neighbor->requests[irequest]->half = 0;
361 neighbor->requests[irequest]->full = 1;
362
363 cut_coulsq = cut_coul * cut_coul;
364
365 // insure use of KSpace long-range solver, set g_ewald
366
367 if (force->kspace == nullptr) error->all(FLERR, "Pair style requires a KSpace style");
368 g_ewald = force->kspace->g_ewald;
369
370 // setup force tables
371
372 if (ncoultablebits) init_tables(cut_coul, cut_respa);
373 }
374
375 /* ---------------------------------------------------------------------- */
376
extract(const char * str,int & dim)377 void *PairLJCutCoulMSMDielectric::extract(const char *str, int &dim)
378 {
379 dim = 0;
380 if (strcmp(str, "cut_coul") == 0) return (void *) &cut_coul;
381 dim = 2;
382 if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
383 if (strcmp(str, "sigma") == 0) return (void *) sigma;
384 return nullptr;
385 }
386