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