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 Contributing authors: Eugen Rozic (University College London)
16 ------------------------------------------------------------------------- */
17
18 #include "pair_cosine_squared.h"
19
20 #include "atom.h"
21 #include "comm.h"
22 #include "error.h"
23 #include "force.h"
24 #include "math_const.h"
25 #include "memory.h"
26 #include "neigh_list.h"
27
28 #include <cmath>
29 #include <cstring>
30
31 using namespace LAMMPS_NS;
32 using namespace MathConst;
33
34 /* ---------------------------------------------------------------------- */
35
PairCosineSquared(LAMMPS * lmp)36 PairCosineSquared::PairCosineSquared(LAMMPS *lmp) : Pair(lmp)
37 {
38 writedata = 1;
39 }
40
41 /* ---------------------------------------------------------------------- */
42
~PairCosineSquared()43 PairCosineSquared::~PairCosineSquared()
44 {
45 if (allocated) {
46 memory->destroy(setflag);
47 memory->destroy(cutsq);
48
49 memory->destroy(epsilon);
50 memory->destroy(sigma);
51 memory->destroy(w);
52 memory->destroy(cut);
53 memory->destroy(wcaflag);
54
55 memory->destroy(lj12_e);
56 memory->destroy(lj6_e);
57 memory->destroy(lj12_f);
58 memory->destroy(lj6_f);
59 }
60 }
61
62 /* ----------------------------------------------------------------------
63 allocate all arrays
64 ------------------------------------------------------------------------- */
65
allocate()66 void PairCosineSquared::allocate()
67 {
68 allocated = 1;
69 int n = atom->ntypes;
70 memory->create(setflag, n+1, n+1, "pair:setflag");
71 memory->create(cutsq, n+1, n+1, "pair:cutsq");
72
73 memory->create(cut, n+1, n+1, "pair:cut");
74 memory->create(epsilon, n+1, n+1, "pair:epsilon");
75 memory->create(sigma, n+1, n+1, "pair:sigma");
76 memory->create(w, n+1, n+1, "pair:w");
77 memory->create(wcaflag, n+1, n+1, "pair:wcaflag");
78
79 memory->create(lj12_e, n+1, n+1, "pair:lj12_e");
80 memory->create(lj6_e, n+1, n+1, "pair:lj6_e");
81 memory->create(lj12_f, n+1, n+1, "pair:lj12_f");
82 memory->create(lj6_f, n+1, n+1, "pair:lj6_f");
83
84 for (int i = 1; i <= n; i++) {
85 for (int j = i; j <= n; j++) {
86 setflag[i][j] = 0;
87 wcaflag[i][j] = 0;
88 }
89 }
90 }
91
92 /* ----------------------------------------------------------------------
93 global settings
94 ------------------------------------------------------------------------- */
95
settings(int narg,char ** arg)96 void PairCosineSquared::settings(int narg, char **arg)
97 {
98 if (narg != 1) {
99 error->all(FLERR, "Illegal pair_style command (wrong number of params)");
100 }
101
102 cut_global = utils::numeric(FLERR, arg[0],false,lmp);
103
104 // reset cutoffs that have been explicitly set
105
106 if (allocated) {
107 int i, j;
108 for (i = 1; i <= atom->ntypes; i++)
109 for (j = i+1; j <= atom->ntypes; j++)
110 if (setflag[i][j])
111 cut[i][j] = cut_global;
112 }
113 }
114
115
116 /* ----------------------------------------------------------------------
117 set coeffs for one or more type pairs
118 ------------------------------------------------------------------------- */
119
coeff(int narg,char ** arg)120 void PairCosineSquared::coeff(int narg, char **arg)
121 {
122 if (narg < 4 || narg > 6)
123 error->all(FLERR, "Incorrect args for pair coefficients (too few or too many)");
124
125 if (!allocated)
126 allocate();
127
128 int ilo, ihi, jlo, jhi;
129 utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
130 utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
131
132 double epsilon_one = utils::numeric(FLERR, arg[2],false,lmp);
133 double sigma_one = utils::numeric(FLERR, arg[3],false,lmp);
134
135 double cut_one = cut_global;
136 double wca_one = 0;
137 if (narg == 6) {
138 cut_one = utils::numeric(FLERR, arg[4],false,lmp);
139 if (strcmp(arg[5], "wca") == 0) {
140 wca_one = 1;
141 } else {
142 error->all(FLERR, "Incorrect args for pair coefficients (unknown option)");
143 }
144 } else if (narg == 5) {
145 if (strcmp(arg[4], "wca") == 0) {
146 wca_one = 1;
147 } else {
148 cut_one = utils::numeric(FLERR, arg[4],false,lmp);
149 }
150 }
151
152 if (cut_one < sigma_one) {
153 error->all(FLERR, "Incorrect args for pair coefficients (cutoff < sigma)");
154 } else if (cut_one == sigma_one) {
155 if (wca_one == 0) {
156 error->all(FLERR, "Incorrect args for pair coefficients (cutoff = sigma w/o wca)");
157 } else {
158 error->warning(FLERR, "Cosine/squared set to WCA only (cutoff = sigma)");
159 }
160 }
161
162 int count = 0;
163 for (int i = ilo; i <= ihi; i++) {
164 for (int j = MAX(jlo,i); j <= jhi; j++) {
165 epsilon[i][j] = epsilon_one;
166 sigma[i][j] = sigma_one;
167 cut[i][j] = cut_one;
168 wcaflag[i][j] = wca_one;
169 setflag[i][j] = 1;
170 count++;
171 }
172 }
173
174 if (count == 0)
175 error->all(FLERR, "Incorrect args for pair coefficients (none set)");
176 }
177
178 /* ----------------------------------------------------------------------
179 init specific to this pair style (unnecessary)
180 ------------------------------------------------------------------------- */
181
182 /*
183 void PairCosineSquared::init_style()
184 {
185 neighbor->request(this,instance_me);
186 }
187 */
188
189 /* ----------------------------------------------------------------------
190 init for one type pair i,j and corresponding j,i
191 ------------------------------------------------------------------------- */
192
init_one(int i,int j)193 double PairCosineSquared::init_one(int i, int j)
194 {
195 if (setflag[i][j] == 0)
196 error->all(FLERR, "Mixing not supported in pair_style cosine/squared");
197
198 epsilon[j][i] = epsilon[i][j];
199 sigma[j][i] = sigma[i][j];
200 cut[j][i] = cut[i][j];
201 wcaflag[j][i] = wcaflag[i][j];
202
203 w[j][i] = w[i][j] = cut[i][j] - sigma[i][j];
204
205 if (wcaflag[i][j]) {
206 lj12_e[j][i] = lj12_e[i][j] = epsilon[i][j] * pow(sigma[i][j], 12.0);
207 lj6_e[j][i] = lj6_e[i][j] = 2.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
208 lj12_f[j][i] = lj12_f[i][j] = 12.0 * epsilon[i][j] * pow(sigma[i][j], 12.0);
209 lj6_f[j][i] = lj6_f[i][j] = 12.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
210 }
211
212 // Note: cutsq is set in pair.cpp
213
214 return cut[i][j];
215 }
216
217 /* ----------------------------------------------------------------------
218 this is here to throw errors & warnings for given options
219 ------------------------------------------------------------------------- */
220
modify_params(int narg,char ** arg)221 void PairCosineSquared::modify_params(int narg, char **arg)
222 {
223 Pair::modify_params(narg, arg);
224
225 int iarg = 0;
226 while (iarg < narg) {
227 if (strcmp(arg[iarg], "mix") == 0) {
228 error->all(FLERR, "pair_modify mix not supported for pair_style cosine/squared");
229 } else if (strcmp(arg[iarg], "shift") == 0) {
230 error->warning(FLERR, "pair_modify shift has no effect on pair_style cosine/squared");
231 offset_flag = 0;
232 } else if (strcmp(arg[iarg], "tail") == 0) {
233 error->warning(FLERR, "pair_modify tail has no effect on pair_style cosine/squared");
234 tail_flag = 0;
235 }
236 iarg++;
237 }
238 }
239
240 /* ----------------------------------------------------------------------
241 proc 0 writes to restart file
242 ------------------------------------------------------------------------- */
243
write_restart(FILE * fp)244 void PairCosineSquared::write_restart(FILE *fp)
245 {
246 write_restart_settings(fp);
247
248 int i, j;
249 for (i = 1; i <= atom->ntypes; i++)
250 for (j = i; j <= atom->ntypes; j++) {
251 fwrite(&setflag[i][j], sizeof(int), 1, fp);
252 if (setflag[i][j]) {
253 fwrite(&epsilon[i][j], sizeof(double), 1, fp);
254 fwrite(&sigma[i][j], sizeof(double), 1, fp);
255 fwrite(&cut[i][j], sizeof(double), 1, fp);
256 fwrite(&wcaflag[i][j], sizeof(int), 1, fp);
257 }
258 }
259 }
260
261 /* ----------------------------------------------------------------------
262 proc 0 reads from restart file, bcasts
263 ------------------------------------------------------------------------- */
264
read_restart(FILE * fp)265 void PairCosineSquared::read_restart(FILE *fp)
266 {
267 read_restart_settings(fp);
268 allocate();
269
270 int i,j;
271 int me = comm->me;
272 for (i = 1; i <= atom->ntypes; i++) {
273 for (j = i; j <= atom->ntypes; j++) {
274 if (me == 0)
275 utils::sfread(FLERR,&setflag[i][j], sizeof(int), 1, fp,nullptr,error);
276 MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
277 if (setflag[i][j]) {
278 if (me == 0) {
279 utils::sfread(FLERR,&epsilon[i][j], sizeof(double), 1, fp,nullptr,error);
280 utils::sfread(FLERR,&sigma[i][j], sizeof(double), 1, fp,nullptr,error);
281 utils::sfread(FLERR,&cut[i][j], sizeof(double), 1, fp,nullptr,error);
282 utils::sfread(FLERR,&wcaflag[i][j], sizeof(int), 1, fp,nullptr,error);
283 }
284 MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
285 MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world);
286 MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
287 MPI_Bcast(&wcaflag[i][j], 1, MPI_INT, 0, world);
288 }
289 }
290 }
291 }
292
293 /* ----------------------------------------------------------------------
294 proc 0 writes to restart file
295 ------------------------------------------------------------------------- */
296
write_restart_settings(FILE * fp)297 void PairCosineSquared::write_restart_settings(FILE *fp)
298 {
299 fwrite(&cut_global, sizeof(double), 1, fp);
300 }
301
302 /* ----------------------------------------------------------------------
303 proc 0 reads from restart file, bcasts
304 ------------------------------------------------------------------------- */
305
read_restart_settings(FILE * fp)306 void PairCosineSquared::read_restart_settings(FILE *fp)
307 {
308 int me = comm->me;
309 if (me == 0) {
310 utils::sfread(FLERR,&cut_global, sizeof(double), 1, fp,nullptr,error);
311 }
312 MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
313 }
314
315 /* ----------------------------------------------------------------------
316 proc 0 writes to data file
317 ------------------------------------------------------------------------- */
318
write_data(FILE * fp)319 void PairCosineSquared::write_data(FILE *fp)
320 {
321 for (int i = 1; i <= atom->ntypes; i++)
322 fprintf(fp, "%d %g %g %g %s\n", i, epsilon[i][i], sigma[i][i], cut[i][i],
323 (wcaflag[i][i] ? "wca" : ""));
324 }
325
326 /* ----------------------------------------------------------------------
327 proc 0 writes all pairs to data file
328 ------------------------------------------------------------------------- */
329
write_data_all(FILE * fp)330 void PairCosineSquared::write_data_all(FILE *fp)
331 {
332 for (int i = 1; i <= atom->ntypes; i++)
333 for (int j = i; j <= atom->ntypes; j++)
334 fprintf(fp, "%d %d %g %g %g %s\n", i, j, epsilon[i][j], sigma[i][j],
335 cut[i][j], (wcaflag[i][j] ? "wca" : ""));
336 }
337
338 /* ---------------------------------------------------------------------- */
339
compute(int eflag,int vflag)340 void PairCosineSquared::compute(int eflag, int vflag)
341 {
342 int i, j, ii, jj, inum, jnum, itype, jtype;
343 int *ilist, *jlist, *numneigh, **firstneigh;
344 double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
345 double r, rsq, r2inv, r6inv;
346 double factor_lj, force_lj, force_cos, cosone;
347
348 ev_init(eflag, vflag);
349 evdwl = 0.0;
350
351 double **x = atom->x;
352 double **f = atom->f;
353 int *type = atom->type;
354 int nlocal = atom->nlocal;
355 double *special_lj = force->special_lj;
356 int newton_pair = force->newton_pair;
357
358 inum = list->inum;
359 ilist = list->ilist;
360 numneigh = list->numneigh;
361 firstneigh = list->firstneigh;
362
363 // loop over neighbors of my atoms
364
365 for (ii = 0; ii < inum; ii++) {
366 i = ilist[ii];
367 xtmp = x[i][0];
368 ytmp = x[i][1];
369 ztmp = x[i][2];
370 itype = type[i];
371 jlist = firstneigh[i];
372 jnum = numneigh[i];
373
374 for (jj = 0; jj < jnum; jj++) {
375 j = jlist[jj];
376 factor_lj = special_lj[sbmask(j)];
377 j &= NEIGHMASK;
378
379 delx = xtmp - x[j][0];
380 dely = ytmp - x[j][1];
381 delz = ztmp - x[j][2];
382 rsq = delx*delx + dely*dely + delz*delz;
383 jtype = type[j];
384
385 if (rsq < cutsq[itype][jtype]) {
386
387 /*
388 This is exactly what the "single" method does, in fact it could be called
389 here instead of repeating the code but here energy calculation is optional
390 so a little bit of calculation is possibly saved
391 */
392
393 r = sqrt(rsq);
394
395 if (r <= sigma[itype][jtype]) {
396 if (wcaflag[itype][jtype]) {
397 r2inv = 1.0/rsq;
398 r6inv = r2inv*r2inv*r2inv;
399 force_lj = r6inv*(lj12_f[itype][jtype]*r6inv - lj6_f[itype][jtype]);
400 fpair = factor_lj*force_lj*r2inv;
401 if (eflag) {
402 evdwl = factor_lj*r6inv*
403 (lj12_e[itype][jtype]*r6inv - lj6_e[itype][jtype]);
404 if (sigma[itype][jtype] == cut[itype][jtype]) {
405 // this is the WCA-only case (it requires this shift by definition)
406 evdwl += factor_lj*epsilon[itype][jtype];
407 }
408 }
409 } else {
410 fpair = 0.0;
411 if (eflag) {
412 evdwl = -factor_lj*epsilon[itype][jtype];
413 }
414 }
415 } else {
416 force_cos = -(MY_PI*epsilon[itype][jtype] / (2.0*w[itype][jtype])) *
417 sin(MY_PI*(r-sigma[itype][jtype]) / w[itype][jtype]);
418 fpair = factor_lj*force_cos / r;
419 if (eflag) {
420 cosone = cos(MY_PI*(r-sigma[itype][jtype]) / (2.0*w[itype][jtype]));
421 evdwl = -factor_lj*epsilon[itype][jtype]*cosone*cosone;
422 }
423 }
424
425 f[i][0] += delx*fpair;
426 f[i][1] += dely*fpair;
427 f[i][2] += delz*fpair;
428
429 if (newton_pair || j < nlocal) {
430 f[j][0] -= delx*fpair;
431 f[j][1] -= dely*fpair;
432 f[j][2] -= delz*fpair;
433 }
434
435 if (evflag)
436 ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
437 }
438 }
439 }
440
441 if (vflag_fdotr)
442 virial_fdotr_compute();
443 }
444
445 /* ----------------------------------------------------------------------
446 This is used be pair_write;
447 it is called only if rsq < cutsq[itype][jtype], no need to check that
448 ------------------------------------------------------------------------- */
449
single(int,int,int itype,int jtype,double rsq,double,double factor_lj,double & fforce)450 double PairCosineSquared::single(int /* i */, int /* j */, int itype, int jtype, double rsq,
451 double /* factor_coul */, double factor_lj,
452 double &fforce)
453 {
454 double r, r2inv, r6inv, cosone, force, energy;
455
456 r = sqrt(rsq);
457
458 if (r <= sigma[itype][jtype]) {
459 if (wcaflag[itype][jtype]) {
460 r2inv = 1.0/rsq;
461 r6inv = r2inv*r2inv*r2inv;
462 force = r6inv*(lj12_f[itype][jtype]*r6inv - lj6_f[itype][jtype])*r2inv;
463 energy = r6inv*(lj12_e[itype][jtype]*r6inv - lj6_e[itype][jtype]);
464 if (sigma[itype][jtype] == cut[itype][jtype]) {
465 // this is the WCA-only case (it requires this shift by definition)
466 energy += epsilon[itype][jtype];
467 }
468 } else {
469 force = 0.0;
470 energy = -epsilon[itype][jtype];
471 }
472 } else {
473 cosone = cos(MY_PI*(r-sigma[itype][jtype]) / (2.0*w[itype][jtype]));
474 force = -(MY_PI*epsilon[itype][jtype] / (2.0*w[itype][jtype])) *
475 sin(MY_PI*(r-sigma[itype][jtype]) / w[itype][jtype]) / r;
476 energy = -epsilon[itype][jtype]*cosone*cosone;
477 }
478 fforce = factor_lj*force;
479 return factor_lj*energy;
480 }
481
482