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: Craig Maloney (UCSB)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_lj_smooth.h"
20
21 #include <cmath>
22 #include "atom.h"
23 #include "comm.h"
24 #include "force.h"
25 #include "neigh_list.h"
26 #include "memory.h"
27 #include "error.h"
28
29
30 using namespace LAMMPS_NS;
31
32 /* ---------------------------------------------------------------------- */
33
PairLJSmooth(LAMMPS * lmp)34 PairLJSmooth::PairLJSmooth(LAMMPS *lmp) : Pair(lmp)
35 {
36 writedata = 1;
37 }
38
39 /* ---------------------------------------------------------------------- */
40
~PairLJSmooth()41 PairLJSmooth::~PairLJSmooth()
42 {
43 if (allocated) {
44 memory->destroy(setflag);
45 memory->destroy(cutsq);
46
47 memory->destroy(cut);
48 memory->destroy(cut_inner);
49 memory->destroy(cut_inner_sq);
50 memory->destroy(epsilon);
51 memory->destroy(sigma);
52 memory->destroy(lj1);
53 memory->destroy(lj2);
54 memory->destroy(lj3);
55 memory->destroy(lj4);
56 memory->destroy(ljsw0);
57 memory->destroy(ljsw1);
58 memory->destroy(ljsw2);
59 memory->destroy(ljsw3);
60 memory->destroy(ljsw4);
61 memory->destroy(offset);
62 }
63 }
64
65 /* ---------------------------------------------------------------------- */
66
compute(int eflag,int vflag)67 void PairLJSmooth::compute(int eflag, int vflag)
68 {
69 int i,j,ii,jj,inum,jnum,itype,jtype;
70 double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
71 double rsq,r2inv,r6inv,forcelj,factor_lj;
72 double r,t,tsq,fskin;
73 int *ilist,*jlist,*numneigh,**firstneigh;
74
75 evdwl = 0.0;
76 ev_init(eflag,vflag);
77
78 double **x = atom->x;
79 double **f = atom->f;
80 int *type = atom->type;
81 int nlocal = atom->nlocal;
82 double *special_lj = force->special_lj;
83 int newton_pair = force->newton_pair;
84
85 inum = list->inum;
86 ilist = list->ilist;
87 numneigh = list->numneigh;
88 firstneigh = list->firstneigh;
89
90 // loop over neighbors of my atoms
91
92 for (ii = 0; ii < inum; ii++) {
93 i = ilist[ii];
94 xtmp = x[i][0];
95 ytmp = x[i][1];
96 ztmp = x[i][2];
97 itype = type[i];
98 jlist = firstneigh[i];
99 jnum = numneigh[i];
100
101 for (jj = 0; jj < jnum; jj++) {
102 j = jlist[jj];
103 factor_lj = special_lj[sbmask(j)];
104 j &= NEIGHMASK;
105
106 delx = xtmp - x[j][0];
107 dely = ytmp - x[j][1];
108 delz = ztmp - x[j][2];
109 rsq = delx*delx + dely*dely + delz*delz;
110 jtype = type[j];
111
112 if (rsq < cutsq[itype][jtype]) {
113 r2inv = 1.0/rsq;
114 if (rsq < cut_inner_sq[itype][jtype]) {
115 r6inv = r2inv*r2inv*r2inv;
116 forcelj = r6inv * (lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
117 } else {
118 r = sqrt(rsq);
119 t = r - cut_inner[itype][jtype];
120 tsq = t*t;
121 fskin = ljsw1[itype][jtype] + ljsw2[itype][jtype]*t +
122 ljsw3[itype][jtype]*tsq + ljsw4[itype][jtype]*tsq*t;
123 forcelj = fskin*r;
124 }
125
126 fpair = factor_lj*forcelj*r2inv;
127
128 f[i][0] += delx*fpair;
129 f[i][1] += dely*fpair;
130 f[i][2] += delz*fpair;
131 if (newton_pair || j < nlocal) {
132 f[j][0] -= delx*fpair;
133 f[j][1] -= dely*fpair;
134 f[j][2] -= delz*fpair;
135 }
136
137 if (eflag) {
138 if (rsq < cut_inner_sq[itype][jtype])
139 evdwl = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]) - offset[itype][jtype];
140 else
141 evdwl = ljsw0[itype][jtype] - ljsw1[itype][jtype]*t -
142 ljsw2[itype][jtype]*tsq/2.0 - ljsw3[itype][jtype]*tsq*t/3.0 -
143 ljsw4[itype][jtype]*tsq*tsq/4.0 - offset[itype][jtype];
144 evdwl *= factor_lj;
145 }
146
147 if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
148 }
149 }
150 }
151
152 if (vflag_fdotr) virial_fdotr_compute();
153 }
154
155 /* ----------------------------------------------------------------------
156 allocate all arrays
157 ------------------------------------------------------------------------- */
158
allocate()159 void PairLJSmooth::allocate()
160 {
161 allocated = 1;
162 int n = atom->ntypes;
163
164 memory->create(setflag,n+1,n+1,"pair:setflag");
165 for (int i = 1; i <= n; i++)
166 for (int j = i; j <= n; j++)
167 setflag[i][j] = 0;
168
169 memory->create(cutsq,n+1,n+1,"pair:cutsq");
170
171 memory->create(cut,n+1,n+1,"pair:cut");
172 memory->create(cut_inner,n+1,n+1,"pair:cut_inner");
173 memory->create(cut_inner_sq,n+1,n+1,"pair:cut_inner_sq");
174 memory->create(epsilon,n+1,n+1,"pair:epsilon");
175 memory->create(sigma,n+1,n+1,"pair:sigma");
176 memory->create(lj1,n+1,n+1,"pair:lj1");
177 memory->create(lj2,n+1,n+1,"pair:lj2");
178 memory->create(lj3,n+1,n+1,"pair:lj3");
179 memory->create(lj4,n+1,n+1,"pair:lj4");
180 memory->create(ljsw0,n+1,n+1,"pair:ljsw0");
181 memory->create(ljsw1,n+1,n+1,"pair:ljsw1");
182 memory->create(ljsw2,n+1,n+1,"pair:ljsw2");
183 memory->create(ljsw3,n+1,n+1,"pair:ljsw3");
184 memory->create(ljsw4,n+1,n+1,"pair:ljsw4");
185 memory->create(offset,n+1,n+1,"pair:offset");
186 }
187
188 /* ----------------------------------------------------------------------
189 global settings
190 ------------------------------------------------------------------------- */
191
settings(int narg,char ** arg)192 void PairLJSmooth::settings(int narg, char **arg)
193 {
194 if (narg != 2) error->all(FLERR,"Illegal pair_style command");
195
196 cut_inner_global = utils::numeric(FLERR,arg[0],false,lmp);
197 cut_global = utils::numeric(FLERR,arg[1],false,lmp);
198
199 if (cut_inner_global <= 0.0 || cut_inner_global > cut_global)
200 error->all(FLERR,"Illegal pair_style command");
201
202 // reset cutoffs that have been explicitly set
203
204 if (allocated) {
205 int i,j;
206 for (i = 1; i <= atom->ntypes; i++)
207 for (j = i; j <= atom->ntypes; j++)
208 if (setflag[i][j]) {
209 cut_inner[i][j] = cut_inner_global;
210 cut[i][j] = cut_global;
211 }
212 }
213 }
214
215 /* ----------------------------------------------------------------------
216 set coeffs for one or more type pairs
217 ------------------------------------------------------------------------- */
218
coeff(int narg,char ** arg)219 void PairLJSmooth::coeff(int narg, char **arg)
220 {
221 if (narg != 4 && narg != 6)
222 error->all(FLERR,"Incorrect args for pair coefficients");
223 if (!allocated) allocate();
224
225 int ilo,ihi,jlo,jhi;
226 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
227 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
228
229 double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp);
230 double sigma_one = utils::numeric(FLERR,arg[3],false,lmp);
231
232 double cut_inner_one = cut_inner_global;
233 double cut_one = cut_global;
234 if (narg == 6) {
235 cut_inner_one = utils::numeric(FLERR,arg[4],false,lmp);
236 cut_one = utils::numeric(FLERR,arg[5],false,lmp);
237 }
238
239 if (cut_inner_one <= 0.0 || cut_inner_one > cut_one)
240 error->all(FLERR,"Incorrect args for pair coefficients");
241
242 int count = 0;
243 for (int i = ilo; i <= ihi; i++) {
244 for (int j = MAX(jlo,i); j <= jhi; j++) {
245 epsilon[i][j] = epsilon_one;
246 sigma[i][j] = sigma_one;
247 cut_inner[i][j] = cut_inner_one;
248 cut[i][j] = cut_one;
249 setflag[i][j] = 1;
250 count++;
251 }
252 }
253
254 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
255 }
256
257 /* ----------------------------------------------------------------------
258 init for one type pair i,j and corresponding j,i
259 ------------------------------------------------------------------------- */
260
init_one(int i,int j)261 double PairLJSmooth::init_one(int i, int j)
262 {
263 if (setflag[i][j] == 0) {
264 epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
265 sigma[i][i],sigma[j][j]);
266 sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]);
267 cut_inner[i][j] = mix_distance(cut_inner[i][i],cut_inner[j][j]);
268 cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
269 }
270
271 cut_inner_sq[i][j] = cut_inner[i][j]*cut_inner[i][j];
272 lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
273 lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
274 lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);
275 lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);
276
277 if (cut_inner[i][j] != cut[i][j]) {
278 double r6inv = 1.0/pow(cut_inner[i][j],6.0);
279 double t = cut[i][j] - cut_inner[i][j];
280 double tsq = t*t;
281 double ratio = sigma[i][j] / cut_inner[i][j];
282 ljsw0[i][j] = 4.0*epsilon[i][j]*(pow(ratio,12.0) - pow(ratio,6.0));
283 ljsw1[i][j] = r6inv*(lj1[i][j]*r6inv-lj2[i][j]) / cut_inner[i][j];
284 ljsw2[i][j] = -r6inv * (13.0*lj1[i][j]*r6inv - 7.0*lj2[i][j]) /
285 cut_inner_sq[i][j];
286 ljsw3[i][j] = -(3.0/tsq) * (ljsw1[i][j] + 2.0/3.0*ljsw2[i][j]*t);
287 ljsw4[i][j] = -1.0/(3.0*tsq) * (ljsw2[i][j] + 2.0*ljsw3[i][j]*t);
288 if (offset_flag)
289 offset[i][j] = ljsw0[i][j] - ljsw1[i][j]*t - ljsw2[i][j]*tsq/2.0 -
290 ljsw3[i][j]*tsq*t/3.0 - ljsw4[i][j]*tsq*tsq/4.0;
291 else offset[i][j] = 0.0;
292 } else {
293 ljsw0[i][j] = 0.0;
294 ljsw1[i][j] = 0.0;
295 ljsw2[i][j] = 0.0;
296 ljsw3[i][j] = 0.0;
297 ljsw4[i][j] = 0.0;
298 double ratio = sigma[i][j] / cut_inner[i][j];
299 if (offset_flag)
300 offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));
301 else offset[i][j] = 0.0;
302 }
303
304 cut_inner[j][i] = cut_inner[i][j];
305 cut_inner_sq[j][i] = cut_inner_sq[i][j];
306 lj1[j][i] = lj1[i][j];
307 lj2[j][i] = lj2[i][j];
308 lj3[j][i] = lj3[i][j];
309 lj4[j][i] = lj4[i][j];
310 ljsw0[j][i] = ljsw0[i][j];
311 ljsw1[j][i] = ljsw1[i][j];
312 ljsw2[j][i] = ljsw2[i][j];
313 ljsw3[j][i] = ljsw3[i][j];
314 ljsw4[j][i] = ljsw4[i][j];
315 offset[j][i] = offset[i][j];
316
317 return cut[i][j];
318 }
319
320 /* ----------------------------------------------------------------------
321 proc 0 writes to restart file
322 ------------------------------------------------------------------------- */
323
write_restart(FILE * fp)324 void PairLJSmooth::write_restart(FILE *fp)
325 {
326 write_restart_settings(fp);
327
328 int i,j;
329 for (i = 1; i <= atom->ntypes; i++)
330 for (j = i; j <= atom->ntypes; j++) {
331 fwrite(&setflag[i][j],sizeof(int),1,fp);
332 if (setflag[i][j]) {
333 fwrite(&epsilon[i][j],sizeof(double),1,fp);
334 fwrite(&sigma[i][j],sizeof(double),1,fp);
335 fwrite(&cut_inner[i][j],sizeof(double),1,fp);
336 fwrite(&cut[i][j],sizeof(double),1,fp);
337 }
338 }
339 }
340
341 /* ----------------------------------------------------------------------
342 proc 0 reads from restart file, bcasts
343 ------------------------------------------------------------------------- */
344
read_restart(FILE * fp)345 void PairLJSmooth::read_restart(FILE *fp)
346 {
347 read_restart_settings(fp);
348 allocate();
349
350 int i,j;
351 int me = comm->me;
352 for (i = 1; i <= atom->ntypes; i++)
353 for (j = i; j <= atom->ntypes; j++) {
354 if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
355 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
356 if (setflag[i][j]) {
357 if (me == 0) {
358 utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error);
359 utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error);
360 utils::sfread(FLERR,&cut_inner[i][j],sizeof(double),1,fp,nullptr,error);
361 utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
362 }
363 MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
364 MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
365 MPI_Bcast(&cut_inner[i][j],1,MPI_DOUBLE,0,world);
366 MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
367 }
368 }
369 }
370
371 /* ----------------------------------------------------------------------
372 proc 0 writes to restart file
373 ------------------------------------------------------------------------- */
374
write_restart_settings(FILE * fp)375 void PairLJSmooth::write_restart_settings(FILE *fp)
376 {
377 fwrite(&cut_inner_global,sizeof(double),1,fp);
378 fwrite(&cut_global,sizeof(double),1,fp);
379 fwrite(&offset_flag,sizeof(int),1,fp);
380 fwrite(&mix_flag,sizeof(int),1,fp);
381 }
382
383 /* ----------------------------------------------------------------------
384 proc 0 reads from restart file, bcasts
385 ------------------------------------------------------------------------- */
386
read_restart_settings(FILE * fp)387 void PairLJSmooth::read_restart_settings(FILE *fp)
388 {
389 int me = comm->me;
390 if (me == 0) {
391 utils::sfread(FLERR,&cut_inner_global,sizeof(double),1,fp,nullptr,error);
392 utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
393 utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
394 utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
395 }
396 MPI_Bcast(&cut_inner_global,1,MPI_DOUBLE,0,world);
397 MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
398 MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
399 MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
400 }
401
402 /* ----------------------------------------------------------------------
403 proc 0 writes to data file
404 ------------------------------------------------------------------------- */
405
write_data(FILE * fp)406 void PairLJSmooth::write_data(FILE *fp)
407 {
408 for (int i = 1; i <= atom->ntypes; i++)
409 fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]);
410 }
411
412 /* ----------------------------------------------------------------------
413 proc 0 writes all pairs to data file
414 ------------------------------------------------------------------------- */
415
write_data_all(FILE * fp)416 void PairLJSmooth::write_data_all(FILE *fp)
417 {
418 for (int i = 1; i <= atom->ntypes; i++)
419 for (int j = i; j <= atom->ntypes; j++)
420 fprintf(fp,"%d %d %g %g %g %g\n",i,j,
421 epsilon[i][j],sigma[i][j],cut_inner[i][j],cut[i][j]);
422 }
423
424 /* ---------------------------------------------------------------------- */
425
single(int,int,int itype,int jtype,double rsq,double,double factor_lj,double & fforce)426 double PairLJSmooth::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
427 double /*factor_coul*/, double factor_lj,
428 double &fforce)
429 {
430 double r2inv,r6inv,forcelj,philj,r,t,tsq,fskin;
431
432 r2inv = 1.0/rsq;
433 if (rsq < cut_inner_sq[itype][jtype]) {
434 r6inv = r2inv*r2inv*r2inv;
435 forcelj = r6inv * (lj1[itype][jtype]*r6inv-lj2[itype][jtype]);
436 } else {
437 r = sqrt(rsq);
438 t = r - cut_inner[itype][jtype];
439 tsq = t*t;
440 fskin = ljsw1[itype][jtype] + ljsw2[itype][jtype]*t +
441 ljsw3[itype][jtype]*tsq + ljsw4[itype][jtype]*tsq*t;
442 forcelj = fskin*r;
443 }
444 fforce = factor_lj*forcelj*r2inv;
445
446 if (rsq < cut_inner_sq[itype][jtype])
447 philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]) - offset[itype][jtype];
448 else
449 philj = ljsw0[itype][jtype] - ljsw1[itype][jtype]*t -
450 ljsw2[itype][jtype]*tsq/2.0 - ljsw3[itype][jtype]*tsq*t/3.0 -
451 ljsw4[itype][jtype]*tsq*tsq/4.0 - offset[itype][jtype];
452 return factor_lj*philj;
453 }
454