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: Paul Crozier (SNL)
16 ------------------------------------------------------------------------- */
17
18 #include "pair_coul_long.h"
19
20 #include "atom.h"
21 #include "comm.h"
22 #include "error.h"
23 #include "force.h"
24 #include "kspace.h"
25 #include "memory.h"
26 #include "neigh_list.h"
27 #include "neighbor.h"
28
29 #include <cmath>
30 #include <cstring>
31
32 using namespace LAMMPS_NS;
33
34 #define EWALD_F 1.12837917
35 #define EWALD_P 0.3275911
36 #define A1 0.254829592
37 #define A2 -0.284496736
38 #define A3 1.421413741
39 #define A4 -1.453152027
40 #define A5 1.061405429
41
42 /* ---------------------------------------------------------------------- */
43
PairCoulLong(LAMMPS * lmp)44 PairCoulLong::PairCoulLong(LAMMPS *lmp) : Pair(lmp)
45 {
46 ewaldflag = pppmflag = 1;
47 ftable = nullptr;
48 qdist = 0.0;
49 cut_respa = nullptr;
50 }
51
52 /* ---------------------------------------------------------------------- */
53
~PairCoulLong()54 PairCoulLong::~PairCoulLong()
55 {
56 if (copymode) return;
57
58 if (allocated) {
59 memory->destroy(setflag);
60 memory->destroy(cutsq);
61
62 memory->destroy(scale);
63 }
64 if (ftable) free_tables();
65 }
66
67 /* ---------------------------------------------------------------------- */
68
compute(int eflag,int vflag)69 void PairCoulLong::compute(int eflag, int vflag)
70 {
71 int i, j, ii, jj, inum, jnum, itable, itype, jtype;
72 double qtmp, xtmp, ytmp, ztmp, delx, dely, delz, ecoul, fpair;
73 double fraction, table;
74 double r, r2inv, forcecoul, factor_coul;
75 double grij, expm2, prefactor, t, erfc;
76 int *ilist, *jlist, *numneigh, **firstneigh;
77 double rsq;
78
79 ecoul = 0.0;
80 ev_init(eflag, vflag);
81
82 double **x = atom->x;
83 double **f = atom->f;
84 double *q = atom->q;
85 int *type = atom->type;
86 int nlocal = atom->nlocal;
87 double *special_coul = force->special_coul;
88 int newton_pair = force->newton_pair;
89 double qqrd2e = force->qqrd2e;
90
91 inum = list->inum;
92 ilist = list->ilist;
93 numneigh = list->numneigh;
94 firstneigh = list->firstneigh;
95
96 // loop over neighbors of my atoms
97
98 for (ii = 0; ii < inum; ii++) {
99 i = ilist[ii];
100 qtmp = q[i];
101 xtmp = x[i][0];
102 ytmp = x[i][1];
103 ztmp = x[i][2];
104 itype = type[i];
105 jlist = firstneigh[i];
106 jnum = numneigh[i];
107
108 for (jj = 0; jj < jnum; jj++) {
109 j = jlist[jj];
110 factor_coul = special_coul[sbmask(j)];
111 j &= NEIGHMASK;
112
113 delx = xtmp - x[j][0];
114 dely = ytmp - x[j][1];
115 delz = ztmp - x[j][2];
116 rsq = delx * delx + dely * dely + delz * delz;
117 jtype = type[j];
118
119 if (rsq < cut_coulsq) {
120 r2inv = 1.0 / rsq;
121 if (!ncoultablebits || rsq <= tabinnersq) {
122 r = sqrt(rsq);
123 grij = g_ewald * r;
124 expm2 = exp(-grij * grij);
125 t = 1.0 / (1.0 + EWALD_P * grij);
126 erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
127 prefactor = qqrd2e * scale[itype][jtype] * qtmp * q[j] / r;
128 forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
129 if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
130 } else {
131 union_int_float_t rsq_lookup;
132 rsq_lookup.f = rsq;
133 itable = rsq_lookup.i & ncoulmask;
134 itable >>= ncoulshiftbits;
135 fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
136 table = ftable[itable] + fraction * dftable[itable];
137 forcecoul = scale[itype][jtype] * qtmp * q[j] * table;
138 if (factor_coul < 1.0) {
139 table = ctable[itable] + fraction * dctable[itable];
140 prefactor = scale[itype][jtype] * qtmp * q[j] * table;
141 forcecoul -= (1.0 - factor_coul) * prefactor;
142 }
143 }
144
145 fpair = forcecoul * r2inv;
146
147 f[i][0] += delx * fpair;
148 f[i][1] += dely * fpair;
149 f[i][2] += delz * fpair;
150 if (newton_pair || j < nlocal) {
151 f[j][0] -= delx * fpair;
152 f[j][1] -= dely * fpair;
153 f[j][2] -= delz * fpair;
154 }
155
156 if (eflag) {
157 if (!ncoultablebits || rsq <= tabinnersq)
158 ecoul = prefactor * erfc;
159 else {
160 table = etable[itable] + fraction * detable[itable];
161 ecoul = scale[itype][jtype] * qtmp * q[j] * table;
162 }
163 if (factor_coul < 1.0) ecoul -= (1.0 - factor_coul) * prefactor;
164 }
165
166 if (evflag) ev_tally(i, j, nlocal, newton_pair, 0.0, ecoul, fpair, delx, dely, delz);
167 }
168 }
169 }
170
171 if (vflag_fdotr) virial_fdotr_compute();
172 }
173
174 /* ----------------------------------------------------------------------
175 allocate all arrays
176 ------------------------------------------------------------------------- */
177
allocate()178 void PairCoulLong::allocate()
179 {
180 allocated = 1;
181 int n = atom->ntypes;
182
183 memory->create(setflag, n + 1, n + 1, "pair:setflag");
184 for (int i = 1; i <= n; i++)
185 for (int j = i; j <= n; j++) setflag[i][j] = 0;
186
187 memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
188
189 memory->create(scale, n + 1, n + 1, "pair:scale");
190 }
191
192 /* ----------------------------------------------------------------------
193 global settings
194 ------------------------------------------------------------------------- */
195
settings(int narg,char ** arg)196 void PairCoulLong::settings(int narg, char **arg)
197 {
198 if (narg != 1) error->all(FLERR, "Illegal pair_style command");
199
200 cut_coul = utils::numeric(FLERR, arg[0], false, lmp);
201 }
202
203 /* ----------------------------------------------------------------------
204 set coeffs for one or more type pairs
205 ------------------------------------------------------------------------- */
206
coeff(int narg,char ** arg)207 void PairCoulLong::coeff(int narg, char **arg)
208 {
209 if (narg != 2) error->all(FLERR, "Incorrect args for pair coefficients");
210 if (!allocated) allocate();
211
212 int ilo, ihi, jlo, jhi;
213 utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
214 utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
215
216 int count = 0;
217 for (int i = ilo; i <= ihi; i++) {
218 for (int j = MAX(jlo, i); j <= jhi; j++) {
219 scale[i][j] = 1.0;
220 setflag[i][j] = 1;
221 count++;
222 }
223 }
224
225 if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
226 }
227
228 /* ----------------------------------------------------------------------
229 init specific to this pair style
230 ------------------------------------------------------------------------- */
231
init_style()232 void PairCoulLong::init_style()
233 {
234 if (!atom->q_flag) error->all(FLERR, "Pair style lj/cut/coul/long requires atom attribute q");
235
236 neighbor->request(this, instance_me);
237
238 cut_coulsq = cut_coul * cut_coul;
239
240 // insure use of KSpace long-range solver, set g_ewald
241
242 if (force->kspace == nullptr) error->all(FLERR, "Pair style requires a KSpace style");
243 g_ewald = force->kspace->g_ewald;
244
245 // setup force tables
246
247 if (ncoultablebits) init_tables(cut_coul, nullptr);
248 }
249
250 /* ----------------------------------------------------------------------
251 init for one type pair i,j and corresponding j,i
252 ------------------------------------------------------------------------- */
253
init_one(int i,int j)254 double PairCoulLong::init_one(int i, int j)
255 {
256 scale[j][i] = scale[i][j];
257 return cut_coul + 2.0 * qdist;
258 }
259
260 /* ----------------------------------------------------------------------
261 proc 0 writes to restart file
262 ------------------------------------------------------------------------- */
263
write_restart(FILE * fp)264 void PairCoulLong::write_restart(FILE *fp)
265 {
266 write_restart_settings(fp);
267
268 for (int i = 1; i <= atom->ntypes; i++)
269 for (int j = i; j <= atom->ntypes; j++) {
270 fwrite(&setflag[i][j], sizeof(int), 1, fp);
271 if (setflag[i][j]) fwrite(&scale[i][j], sizeof(double), 1, fp);
272 }
273 }
274
275 /* ----------------------------------------------------------------------
276 proc 0 reads from restart file, bcasts
277 ------------------------------------------------------------------------- */
278
read_restart(FILE * fp)279 void PairCoulLong::read_restart(FILE *fp)
280 {
281 read_restart_settings(fp);
282
283 allocate();
284
285 int i, j;
286 int me = comm->me;
287 for (i = 1; i <= atom->ntypes; i++)
288 for (j = i; j <= atom->ntypes; j++) {
289 if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
290 MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
291 if (setflag[i][j]) {
292 if (me == 0) utils::sfread(FLERR, &scale[i][j], sizeof(double), 1, fp, nullptr, error);
293 MPI_Bcast(&scale[i][j], 1, MPI_DOUBLE, 0, world);
294 }
295 }
296 }
297
298 /* ----------------------------------------------------------------------
299 proc 0 writes to restart file
300 ------------------------------------------------------------------------- */
301
write_restart_settings(FILE * fp)302 void PairCoulLong::write_restart_settings(FILE *fp)
303 {
304 fwrite(&cut_coul, sizeof(double), 1, fp);
305 fwrite(&offset_flag, sizeof(int), 1, fp);
306 fwrite(&mix_flag, sizeof(int), 1, fp);
307 fwrite(&ncoultablebits, sizeof(int), 1, fp);
308 fwrite(&tabinner, sizeof(double), 1, fp);
309 }
310
311 /* ----------------------------------------------------------------------
312 proc 0 reads from restart file, bcasts
313 ------------------------------------------------------------------------- */
314
read_restart_settings(FILE * fp)315 void PairCoulLong::read_restart_settings(FILE *fp)
316 {
317 if (comm->me == 0) {
318 utils::sfread(FLERR, &cut_coul, sizeof(double), 1, fp, nullptr, error);
319 utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
320 utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
321 utils::sfread(FLERR, &ncoultablebits, sizeof(int), 1, fp, nullptr, error);
322 utils::sfread(FLERR, &tabinner, sizeof(double), 1, fp, nullptr, error);
323 }
324 MPI_Bcast(&cut_coul, 1, MPI_DOUBLE, 0, world);
325 MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
326 MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
327 MPI_Bcast(&ncoultablebits, 1, MPI_INT, 0, world);
328 MPI_Bcast(&tabinner, 1, MPI_DOUBLE, 0, world);
329 }
330
331 /* ---------------------------------------------------------------------- */
332
single(int i,int j,int,int,double rsq,double factor_coul,double,double & fforce)333 double PairCoulLong::single(int i, int j, int /*itype*/, int /*jtype*/, double rsq,
334 double factor_coul, double /*factor_lj*/, double &fforce)
335 {
336 double r2inv, r, grij, expm2, t, erfc, prefactor;
337 double fraction, table, forcecoul, phicoul;
338 int itable;
339
340 r2inv = 1.0 / rsq;
341 if (!ncoultablebits || rsq <= tabinnersq) {
342 r = sqrt(rsq);
343 grij = g_ewald * r;
344 expm2 = exp(-grij * grij);
345 t = 1.0 / (1.0 + EWALD_P * grij);
346 erfc = t * (A1 + t * (A2 + t * (A3 + t * (A4 + t * A5)))) * expm2;
347 prefactor = force->qqrd2e * atom->q[i] * atom->q[j] / r;
348 forcecoul = prefactor * (erfc + EWALD_F * grij * expm2);
349 if (factor_coul < 1.0) forcecoul -= (1.0 - factor_coul) * prefactor;
350 } else {
351 union_int_float_t rsq_lookup;
352 rsq_lookup.f = rsq;
353 itable = rsq_lookup.i & ncoulmask;
354 itable >>= ncoulshiftbits;
355 fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
356 table = ftable[itable] + fraction * dftable[itable];
357 forcecoul = atom->q[i] * atom->q[j] * table;
358 if (factor_coul < 1.0) {
359 table = ctable[itable] + fraction * dctable[itable];
360 prefactor = atom->q[i] * atom->q[j] * table;
361 forcecoul -= (1.0 - factor_coul) * prefactor;
362 }
363 }
364 fforce = forcecoul * r2inv;
365
366 if (!ncoultablebits || rsq <= tabinnersq)
367 phicoul = prefactor * erfc;
368 else {
369 table = etable[itable] + fraction * detable[itable];
370 phicoul = atom->q[i] * atom->q[j] * table;
371 }
372 if (factor_coul < 1.0) phicoul -= (1.0 - factor_coul) * prefactor;
373
374 return phicoul;
375 }
376
377 /* ---------------------------------------------------------------------- */
378
extract(const char * str,int & dim)379 void *PairCoulLong::extract(const char *str, int &dim)
380 {
381 if (strcmp(str, "cut_coul") == 0) {
382 dim = 0;
383 return (void *) &cut_coul;
384 }
385 if (strcmp(str, "scale") == 0) {
386 dim = 2;
387 return (void *) scale;
388 }
389 return nullptr;
390 }
391